文章介绍了旅行商问题(TSP)的解决方法，包括Christofides算法构建的上界和使用Blossom算法的示例，以及通过松弛问题和子回路约束寻找下界。此外，讨论了如何通过添加约束获取整数解，如polyhedralcombinatorics和分支定界法，并提到了PQ树和shrinking在寻找subtourconstraints中的应用。
1. 问题概述
tsp问题的数学模型如下：
1.1 给定upbound的Christofides方法
这是可以给出上界的一个方法，可以证明构造出的路线不超过最优路线的1.5倍。步骤为：
1）构造MST（最小生成树）
2）将里面的奇点连接起来构成欧拉回路称为完美匹配。Edmonds给出了多项式时间内构造最小代价完美匹配的方法，其长度不超过最优解的1.5倍。
证明方法也很直观，奇点最短路径可以拆分成两条完美匹配，其中总有一条的长度≤\le≤最优奇点路径长度/2≤\le≤最优完整路径长度/2
3）按欧拉回路顺序逐个扫描点，跳过重复经过的点，即可构造一条完整的路径。
下面是blossom算法的使用示例：
using BlossomV,TravelingSalesmanHeuristics,TSPLIB,Distances,DataStructures
tsp = readTSP("vlsi/bcl380.tsp"); # bays29
N = tsp.dimension
distmat = [euclidean(tsp.nodes[i,:],tsp.nodes[j,:]) for i in 1:N, j in 1:N]
mst = TravelingSalesmanHeuristics.minspantree(distmat)[1]
x = counter(cat([m[1] for m in mst],[m[2] for m in mst],dims=1))
odds = []
for xi in x
if xi[2]%2==1
append!(odds,xi[1])
end
end
mat = Matching(Float64, length(odds))
for i in 1:length(odds)
for j in 1:i-1
add_edge(mat,i-1,j-1,distmat[odds[i],odds[j]])
end
end
solve(mat)
using PyPlot
tree = mst
PyPlot.scatter(tsp.nodes[:,1],tsp.nodes[:,2],color="black",s=10)
for i in 1:length(tree)
l = tsp.nodes[[tree[i][1],tree[i][2]],:]
PyPlot.plot(l[:,1], l[:,2], color="black",linewidth = 1)
end
for i in 1:length(odds)
j = get_match(mat,i-1)+1
if j>i
l = tsp.nodes[[odds[i],odds[j]],:]
PyPlot.plot(l[:,1],l[:,2],color="r",linewidth = 0.5,linestyle="--")
end
end
1.2. 给定lowerbound的松弛问题
首先看一个例子：
我们定义xijx_{ij}xij为边ij是否在路径上。根据TSP问题的要求，每个点只能路过一次，因此每个点连着两条边，有：
1.3 MIP求解方法
可以参考TravelingSalesmanExact 包：
using TravelingSalesmanExact, GLPK
tsp = readTSPLIB(:att48)
cities = [tsp.nodes[i,:] for i in 1:tsp.dimension];
@time tour, cost = get_optimal_tour(cities, GLPK.Optimizer)
我们自己实现的代码如下。如果我们能够调用MIP求解器，那么我们可以使用combination包把所有可能的子集罗列出来求解。当然也可以松弛subtour约束后使用行生成求解：
using TSPLIB,JuMP, GLPK, Distances
tsp = readTSP("25.tsp")
N = tsp.dimension
distmat = [euclidean(tsp.nodes[i,:],tsp.nodes[j,:]) for i in 1:N, j in 1:N]
m = Model(GLPK.Optimizer)
@variable(m, x[1:N,1:N], Bin)
@objective(m, Min, sum(x[i,j]*distmat[i,j] for i=1:N,j=1:N))
for i=1:N
@constraint(m, x[i,i] == 0)
@constraint(m, sum(x[i,1:N]) == 1)
end
for j=1:N
@constraint(m, sum(x[1:N,j]) == 1)
end
for f=1:N, t=1:N
@constraint(m, x[f,t]+x[t,f] <= 1)
end
optimize!(m)
结果如下：
接下来不断增加subtour约束重新求解
function is_tsp_solved(m)
g = JuMP.value.(x)
N = size(g,1)
path = Int[]
push!(path,1)
while true
v, idx = findmax(g[path[end],:])
if idx == path[1]
break
else
push!(path,idx)
end
end
n = length(path)
if n < N
@constraint(m,sum(x[path,path])<=n-1)
return false
end
return true
end
while !is_tsp_solved(m)
optimize!(m)
end
最终结果如下
整数规划还是比较费时间的，att48用了177s才求出最优解。下面进行优化，每次同时把所有的subtour都添加进来。我们修改后的代码如下：
using JuMP, GLPK, Distances
function get_cycles(perm_matrix)
N = size(perm_matrix, 1)
remaining_inds = Set(1:N)
cycles = Vector{Int}[]
while length(remaining_inds) > 0
cycle = find_cycle(perm_matrix, first(remaining_inds))
push!(cycles, cycle)
setdiff!(remaining_inds, cycle)
end
return cycles
end
function find_cycle(perm_matrix, starting_ind = 1)
cycle = [starting_ind]
prev_ind = ind = starting_ind
while true
next_ind = findfirst(>(0.5), @views(perm_matrix[ind, 1:prev_ind-1]))
if isnothing(next_ind)
next_ind = findfirst(>(0.5), @views(perm_matrix[ind, prev_ind+1:end])) +
prev_ind
end
next_ind == starting_ind && break
push!(cycle, next_ind)
prev_ind, ind = ind, next_ind
end
return cycle
end
function is_tsp_solved(m,x)
cycles = get_cycles(JuMP.value.(x))
if size(cycles,1)>1
for cycle in cycles
@constraint(m,sum(x[cycle,cycle])<=2*size(cycle,1)-2)
end
return false
end
return true
end
function solve_tsp_with_mip(tsp)
N = tsp.dimension
distmat = [euclidean(tsp.nodes[i,:],tsp.nodes[j,:]) for i in 1:N, j in 1:N]
m = Model(GLPK.Optimizer)
@variable(m, x[1:N,1:N], Symmetric, Bin)
@objective(m, Min, sum(x[i,j]*distmat[i,j] for i=1:N,j=1:i))
for i=1:N
@constraint(m, x[i,i] == 0)
@constraint(m, sum(x[i,:]) == 2)
end
optimize!(m)
while !is_tsp_solved(m,x)
optimize!(m)
end
return JuMP.value.(x), objective_value(m)
end
@time solve_tsp_with_mip(tsp)
1.4 BFF建模方法
给路径上的边定义一个数，第一条边为1，第二条边为2，……，并且定义方向，则除了depot外每个点的所有边加起来为1，设计模型如下：
using DataFrames, Plots, LinearAlgebra, Random
using TSPLIB,LKH,Distances,PyPlot
tsp=readTSPLIB(:att48)
c = tsp.dimension
distance = 1000;
xpos = tsp.nodes[:,1]; ypos = tsp.nodes[:,2];
Nodes = DataFrames.DataFrame(;k = 1:c,x = xpos,y = ypos,);
Distance = zeros(Int(c*(c-1)/2),3)
global cont = 1
for i = 1:size(Nodes,1)
for j = i+1:size(Nodes,1)
global Distance[cont,:] =
[i,j,floor(0.5+norm([(Nodes.x[i,1]-Nodes.x[j,1]),(Nodes.y[i,1]-Nodes.y[j,1])]))]
global cont = cont+1;
end
end
L = size(Distance,1)
branches = DataFrames.DataFrame(; k = Distance[:,1],m = Distance[:,2], lkm = Distance[:,3],)
# 添加model求解代码
using JuMP, HiGHS
b = size(branches, 1) # 边的个数
A = zeros(c, b) # 点是否为变的起点或者终点
for l = 1:b
k = Int(branches.k[l, 1]);
m = Int(branches.m[l, 1])
A[k, l] = 1; A[m, l] = -1;
end
Branch = 1:b; Cities = 1:c;
model = Model(HiGHS.Optimizer)
set_silent(model)
sd = ones(c, 1); sd[1, 1] = 0; smax = zeros(c, 1); smax[1, 1] = c; flmax = c;
@variable(model, x[l in Branch], Bin) # 是否包含边
@variable(model, f[l in Branch]) # 边上的流量
@variable(model, s[k in Cities] >= 0)
@objective(model, Min, sum(branches.lkm[l, 1] * x[l] for l in Branch)) # 总距离最短
for k in Cities
@constraint(model, s[k] <= smax[k])
@constraint(model, sum(abs(A[k, l]) * x[l] for l in Branch) == 2) # 每个点包含2个边
@constraint(model, s[k] - sd[k] == sum(-A[k, l] * f[l] for l in Branch)) # depot的流为47（因为有一个回来的48），其他点的流为1
end
for l in Branch
@constraint(model, -flmax * x[l] <= f[l]) # 只有有边才能有流
@constraint(model, f[l] <= flmax * x[l])
end
@constraint(model, sum(x[i] for i in Branch) == c) # 总共有c条边
@time JuMP.optimize!(model)
@show objective_value(model)
# 添加可视化代码
global sol = [value(x[l]) for l in Branch];
global sou = [value(s[k]) for k in Cities];
global flow = [value(f[l]) for l in Branch];
for l = 1:size(branches, 1)
if sol[l, 1] >= 1/2
local index1 = findall(x -> x == branches.k[l, 1], Nodes.k)
local index2 = findall(x -> x == branches.m[l, 1], Nodes.k)
local x = [Nodes.x[index1, 1]; Nodes.x[index2, 1]];
local y = [Nodes.y[index1, 1]; Nodes.y[index2, 1]];
plot!(x, y, linewidth=3, label="", linecolor=:black)
end
end
scatter!(Nodes.x, Nodes.y, label="", mc=:blue, ms=5, ma=2)
1.5 BFF加速
这里做一些模型优化，首先松弛调depot的约束，其次每个点只考虑附近的K个点，模型如下：
using DataFrames, Random, TSPLIB,Distances, PyPlot, ProgressBars
plt.figure(figsize=(20, 10))
function plot_tree(tsp,mst)
s = tsp.nodes;scatter(s[:,1],s[:,2]);for m in mst;plot(s[[m...],1],s[[m...],2],c="b");end #for p in 1:tsp.dimension;text(s[p,1],s[p,2],p);end;
end
function greedyConstruct(dist)
path = [1];for _ in 1:size(dist)[1];c_dist = dist[path[end],:];c_dist[path].=MaxNum;next_id = findmin(c_dist)[2];push!(path,next_id);end
return path
end
neighbour_num = 15
MaxNum = 1000000
tsp=readTSPLIB(:d493)
println("Generate branches")
c = tsp.dimension
xpos = tsp.nodes[:,1]; ypos = tsp.nodes[:,2];
Nodes = DataFrames.DataFrame(;k = 1:c,x = xpos,y = ypos,);
branch_list = [Set() for i in 1:c]
for i in tqdm(1:c)
dist = [floor.(0.5+euclidean(tsp.nodes[i,:],tsp.nodes[j,:])) for j in 1:c];dist[i] = MaxNum
for j in partialsortperm(dist, 1:neighbour_num)
if i<j
push!(branch_list[i],(i,j,dist[j]))
else
push!(branch_list[j],(j,i,dist[j]))
end
end
end
println("Generate greedy route")
dist_matrix = [floor.(0.5+euclidean(tsp.nodes[i,:],tsp.nodes[j,:])) for i in 1:tsp.dimension,j in 1:tsp.dimension];
path = greedyConstruct(dist_matrix);
for idx in tqdm(1:length(path))
i = path[idx]
j = path[idx==length(path) ? 1 : idx+1]
if i<j
push!(branch_list[i],(i,j,dist_matrix[i,j]))
else
push!(branch_list[j],(j,i,dist_matrix[i,j]))
end
end
branches = DataFrames.DataFrame(union(branch_list...),[:k,:m,:lkm]);
println("Finish branches")
# 添加model求解代码
using JuMP, Gurobi
b = size(branches, 1) # 边的个数
A = zeros(c, b) # 点是否为变的起点或者终点
for l = 1:b
k = Int(branches.k[l, 1]);
m = Int(branches.m[l, 1])
A[k, l] = 1; A[m, l] = -1;
end
Branch = 1:b; Cities = 2:c;
model = Model(Gurobi.Optimizer)
println("start")
#set_silent(model)
set_time_limit_sec(model, 150.0)
@variable(model, x[l in Branch], Bin) # 是否包含边
@variable(model, f[l in Branch]) # 边上的流量
@objective(model, Min, sum(branches.lkm[l, 1] * x[l] for l in Branch)) # 总距离最短
for k in Cities
@constraint(model, sum(abs(A[k, l]) * x[l] for l in Branch) == 2) # 每个点包含2个边
@constraint(model, 1 == sum(A[k, l] * f[l] for l in Branch))
end
for l in Branch
@constraint(model, -c * x[l] <= f[l]) # 只有有边才能有流
@constraint(model, f[l] <= c * x[l])
end
@constraint(model, sum(x[i] for i in Branch) == c) # 总共有c条边
@time JuMP.optimize!(model)
@show objective_value(model)
println("finish")
xv = value.(x);mst = []
for i in 1:size(branches, 1);if xv[i]>0.99;push!(mst,(branches[i,"k"],branches[i,"m"]));end;end
plot_tree(tsp,mst)
1.6 一个简单的优化算法
我们可以拆除一些边，然后重新优化。这里选取距离最长的10%的边，取消后重新组合优化，代码如下：
# 选取10%的距离最远的点进行重新优化
xv = value.(x);mst = []
for i in 1:size(branches, 1);if xv[i]>0.99;push!(mst,(branches[i,"lkm"],branches[i,"k"],branches[i,"m"]));end;end
sort!(mst);
# 选取10%距离最远的边进行拆解
percent_to_remove = 0.1
num_edges_to_remove = round(Int, length(mst) * percent_to_remove)
# 函数：查找连通分量
function find_connected_components(edges, total_nodes)
# 构建邻接表
visit_num = zeros(total_nodes)
adj_list = [[] for _ in 1:total_nodes]
for (_, u, v) in edges
push!(adj_list[u], v)
push!(adj_list[v], u)
visit_num[u] += 1
visit_num[v] += 1
end
visited = falses(total_nodes)
components = []
function dfs(node, component)
visited[node] = true
push!(component, node)
for neighbor in adj_list[node]
if !visited[neighbor]
dfs(neighbor, component)
end
end
end
for node in 1:total_nodes
if !visited[node]
component = []
dfs(node, component)
sort!(component,by = x->visit_num[x])
push!(components, component)
end
end
return components
end
# 按照距离从大到小排序，选取最远的10%的边
farthest_edges = sort(mst, by=x->x[1], rev=true)[1:num_edges_to_remove]
println("\n移除的边数量: ", num_edges_to_remove)
# for edge in farthest_edges;println("距离: ", edge[1], " 点对: (", edge[2], ", ", edge[3], ")");end
# 创建被移除边的集合
removed_edges_set = Set()
for edge in farthest_edges;push!(removed_edges_set, edge[2] < edge[3] ? (edge[2], edge[3]) : (edge[3], edge[2]));end
# 创建剩余的边集合
remaining_edges = []
for edge in mst
edge_key = edge[2] < edge[3] ? (edge[2], edge[3]) : (edge[3], edge[2])
if !(edge_key in removed_edges_set)
push!(remaining_edges, edge)
end
end
# 查找连通分量
components = find_connected_components(remaining_edges, c)
println("\n连通分量数量: ", length(components))
# 为每个连通分量创建虚拟边（选择第一个点和最后一个点作为起点和终点）
virtual_edges = []
for (i, component) in enumerate(components)
if length(component) >= 1
# 为简化，我们选择分量中的第一个点和最后一个点作为虚拟边的端点
# 在实际应用中，可能需要更复杂的方法来选择端点
start_node = component[1]
end_node = component[min(2,length(component))]
# 计算虚拟边的距离
virtual_distance = round(Int, euclidean(tsp.nodes[start_node, :], tsp.nodes[end_node, :]))
push!(virtual_edges, (virtual_distance, start_node, end_node, i))
# println("连通分量 ", i, ": 包含 ", length(component), " 个点, 虚拟边: (", start_node, ", ", end_node, "), 距离: ", virtual_distance)
end
end
# 创建新的节点集：所有连通分量的起点和终点
new_nodes = Set()
for (_, start, stop, _) in virtual_edges
push!(new_nodes, start)
push!(new_nodes, stop)
end
new_nodes = collect(new_nodes)
sort!(new_nodes)
println("\n新的节点集大小: ", length(new_nodes))
# 创建节点映射：原始节点ID -> 新节点索引
node_map = Dict{Int, Int}()
for (i, node) in enumerate(new_nodes)
node_map[node] = i
end
# 为新的节点集创建完全连接的边
new_branches = []
new_c = length(new_nodes)
for i in 1:new_c-1
for j in i+1:new_c
u = new_nodes[i]
v = new_nodes[j]
distance = round(Int, euclidean(tsp.nodes[u, :], tsp.nodes[v, :]))
push!(new_branches, (u, v, distance))
end
end
# 转换为DataFrame
new_branches_df = DataFrame(k=[b[1] for b in new_branches],
m=[b[2] for b in new_branches],
lkm=[b[3] for b in new_branches])
println("\n新的边数量: ", size(new_branches_df, 1))
# 使用新的节点集和边重新建立模型并求解
println("\n重新优化路径...")
new_b = size(new_branches_df, 1)
new_A = zeros(new_c, new_b)
for l = 1:new_b
k = findfirst(==(new_branches_df.k[l]), new_nodes)
m = findfirst(==(new_branches_df.m[l]), new_nodes)
new_A[k, l] = 1; new_A[m, l] = -1
end
new_Branch = 1:new_b
new_Cities = 2:new_c
new_model = Model(Gurobi.Optimizer)
set_silent(new_model)
set_time_limit_sec(new_model, 150.0)
@variable(new_model, new_x[l in new_Branch], Bin)
@variable(new_model, new_f[l in new_Branch])
@objective(new_model, Min, sum(new_branches_df.lkm[l] * new_x[l] for l in new_Branch))
for k in new_Cities
@constraint(new_model, sum(abs(new_A[k, l]) * new_x[l] for l in new_Branch) == 2)
@constraint(new_model, 1 == sum(new_A[k, l] * new_f[l] for l in new_Branch))
end
ve = [(v[2],v[3]) for v in virtual_edges]
for l in new_Branch
# 虚拟边的x一定为1
if (new_branches_df[l, "k"], new_branches_df[l, "m"]) in ve|| (new_branches_df[l, "m"], new_branches_df[l, "k"]) in ve
@constraint(new_model, new_x[l] == 1)
end
@constraint(new_model, -new_c * new_x[l] <= new_f[l])
@constraint(new_model, new_f[l] <= new_c * new_x[l])
end
@constraint(new_model, sum(new_x[i] for i in new_Branch) == new_c)
@time JuMP.optimize!(new_model)
new_optimal_value = objective_value(new_model)
println("\n新节点集的最优路径总距离: ", new_optimal_value)
# 收集新的解
new_xv = value.(new_x)
new_optimal_edges = []
for i in 1:size(new_branches_df, 1)
if new_xv[i] > 0.99
push!(new_optimal_edges, (new_branches_df[i, "lkm"],
new_branches_df[i, "k"],
new_branches_df[i, "m"]))
end
end
# 重建完整路径：结合连通分量内部路径和新优化的虚拟边连接
function reconstruct_full_path(components, component_edges, virtual_edges_optimal)
full_path = []
# 构建连通分量的映射：节点 -> 所属分量
node_to_component = Dict()
for (i, component) in enumerate(components)
for node in component
node_to_component[node] = i
end
end
# 构建连通分量的边映射
component_edges_map = Dict()
virtual_edges = []
for (i, c) in enumerate(components)
component_edges_map[i] = []
push!(virtual_edges,(c[1],c[min(2,length(c))]))
push!(virtual_edges,(c[min(2,length(c))],c[1]))
end
for edge in component_edges
u, v = edge[2], edge[3]
comp_id = node_to_component[u] # u和v属于同一个分量
push!(component_edges_map[comp_id], edge)
end
# 将虚拟边的最优连接加入完整路径
for (dist, u, v) in virtual_edges_optimal
if !((u,v) in virtual_edges)
push!(full_path, (dist, u, v))
end
end
# 将各分量内部的边加入完整路径
for (_, edges) in component_edges_map
append!(full_path, edges)
end
return full_path
end
# 重建完整路径
full_path = reconstruct_full_path(components, remaining_edges, new_optimal_edges)
# 计算完整路径的总距离（包括虚拟边和分量内部边）
full_path_distance = sum(edge[1] for edge in full_path)
improvement = (objective_value(model) - full_path_distance) / objective_value(model) * 100
println("原路径总距离: ", objective_value(model))
println("完整路径总距离: ", full_path_distance)
# 绘制原图
figure(1)
plt.figure(figsize=(20, 10))
s = tsp.nodes;for m in mst;plot(s[[m[2], m[3]],1],s[[m[2], m[3]],2],c="r",ls = "dotted", linewidth=2);end
plot_tree(tsp, [(m[2], m[3]) for m in remaining_edges])
colors = ["r", "g", "b", "c", "m", "y", "k"]
for (i, component) in enumerate(components)
color_idx = (i-1) % length(colors) + 1
# 高亮显示该分量中的点
scatter(tsp.nodes[component, 1], tsp.nodes[component, 2], c=colors[color_idx], s=100, alpha=0.5)
end
title("original TSP (distance: " * string(round(objective_value(model))) * ")")
# 绘制移除边后的图
figure(2)
plt.figure(figsize=(20, 10))
plot_tree(tsp, [(m[2], m[3]) for m in remaining_edges])
# 高亮显示连通分量
colors = ["r", "g", "b", "c", "m", "y", "k"]
for (i, component) in enumerate(components)
color_idx = (i-1) % length(colors) + 1
# 高亮显示该分量中的点
scatter(tsp.nodes[component, 1], tsp.nodes[component, 2], c=colors[color_idx], s=100, alpha=0.5)
end
for (_, u, v) in new_optimal_edges
if !((u,v) in ve || (v,u) in ve)
plot(tsp.nodes[[u, v], 1], tsp.nodes[[u, v], 2], c="r", ls = "dotted",linewidth=2)
end
end
title("reconstructed TSP (distance: " * string(round(full_path_distance)) * ")")
show()
2. 行生成算法
2.1 子回路约束
如下图，上面的模型无法避免subtour
为了消除两个子路径的情形，再增添约束条件：
要想完全枚举出所有这样的约束条件很难，在城市数目为10时，不等式数目就达到51043900866个。因此我们常采用行生成（也叫割平面）的方法，从松弛问题开始，每次求解完后，若发现有子路径，则不断增添拆散子路径的约束条件即可。示意如下，其中红色表示0.5，黑色表示1：
也可以使用最大流最小割来确定subtour
如上图，我们随意选定两个点，得到它们之间的最大流，这个值等于最小割。因此，我们计算n-1个最大流，如果得到小于2的值，那么找到对应的最小割加入约束集。
实现代码如下：
using TravelingSalesmanExact, GLPK,TSPLIB,Distances,JuMP
function get_tours(perm_matrix, thresh = 0.1)
N = size(perm_matrix, 1)
remaining_inds = Set(1:N)
tours = Vector{Int}[]
while length(remaining_inds) > 0
tour = find_tour(perm_matrix, first(remaining_inds), thresh)
push!(tours, tour)
setdiff!(remaining_inds, tour)
end
return tours
end
function find_tour(perm_matrix, starting_ind, thresh)
tour = []
search = [starting_ind]
n = size(perm_matrix,1)
while length(search)>0
v=pop!(search)
if !(v in tour)
push!(tour,v)
for i in 1:n
if perm_matrix[v,i]>=thresh && !(i in tour)
push!(search, i)
end
end
end
end
return tour
end
tsp = readTSPLIB(:bays29); # bays29
N = tsp.dimension
distmat = [euclidean(tsp.nodes[i,:],tsp.nodes[j,:]) for i in 1:N, j in 1:N]
m = Model(GLPK.Optimizer)
@variable(m, 0<=x[1:N,1:N]<=1, Symmetric)
@objective(m, Min, sum(x[i,j]*distmat[i,j] for i=1:N,j=1:i))
for i=1:N
@constraint(m, x[i,i] == 0)
@constraint(m, sum(x[i,:]) == 2)
end
@time optimize!(m)
thresh = 0.1
tours = get_tours(JuMP.value.(x), thresh)
while length(tours) > 1
for tour in tours
@constraint(m,sum(x[tour,tour])<=2*size(tour,1)-2)
end
@time optimize!(m)
tours = get_tours(JuMP.value.(x), thresh)
end
# println(JuMP.objective_value(m))
plot_matrix(tsp,round.(JuMP.value.(x); digits=1))
最终结果如下：
2.2 梳子约束
接下来，为了取整，我们再添加4集合约束：
后面逐步增加个数，称为梳子不等式，每条回路穿过边界的数目至少是3k+1次，其中k是梳齿的数目。
实战中可以使用启发式方法，如下图，对每个红色子图构造梳子约束：
生成梳子的研究可以参考：Sylvia Boyd/Sally Cockburn/ Danielle Vella；Daniel Espinoza/Marcos Goycoolea。
实现代码如下：
function find_comb(sol, odds, starting_ind = 1)
comb = Vector{Vector{Int}}()
teeths = Vector{Vector{Int}}()
tour = Vector{Int}()
search = [starting_ind]
n = size(odds,1)
N = size(sol,1)
while length(search)>0
v = odds[pop!(search)]
if !(v in tour)
push!(tour,v)
for i in 1:n
if sol[v,odds[i]]>0.01 && sol[v,odds[i]]<0.99 && !(odds[i] in tour)
push!(search, i)
end
end
teeth = [v]
for i in 1:N
if sol[v,i]==1
push!(teeth, i)
end
end
if length(teeth)>1
push!(teeths,teeth)
end
end
end
for teeth in teeths
if !(teeth[2] in tour)
push!(comb,teeth)
end
end
push!(comb,tour)
return comb
end
function get_combs(sol)
N = size(sol, 1)
odds = Vector{Int}()
for i in 1:N
for j in 1:N
if 0< sol[i,j] && sol[i,j]<1
if !(i in odds);push!(odds,i);end
if !(j in odds);push!(odds,j);end
end
end
end
combs = Vector{Vector{Vector{Int}}}()
while length(odds) > 0
comb = find_comb(sol, odds)
push!(combs, comb)
setdiff!(odds, comb[end])
end
return combs
end
combs = get_combs(JuMP.value.(x))
while length(combs)>0
for comb in combs
k = length(comb)-1
@constraint(m,sum([sum(x[teeth,setdiff(1:N,teeth)]) for teeth in comb])>=3*k+1)
end
optimize!(m)
combs = get_combs(JuMP.value.(x))
end
plot_matrix(tsp,round.(JuMP.value.(x); digits=1))
2.3 blossom不等式
可以看作subtour约束的加强版。选取顶点数目为奇数的簇，完美匹配中至少有一条边穿过簇的边界（相对应的，subtour约束的右边是2）。
3. 分支定界法
3.1 polyhedral combinatorics
我们以一个42城市为例，使用LP和subtour约束，得到结果如下：
虚线部分表示数值为0.5，其他部分数值为1，其中long path是城市30-42。
添加如下两个约束:
可以获得最终的整数解。
上面这种添加约束的方法，称为polyhedral combinatorics，即寻找所有整数解满足但是LP最优解不满足的约束。
当添加约束的方法对lb的改善低于一定阈值时，我们开始使用branch-and-bound方法来求解问题。
3.2 Tight sets和PQ树
tight set指的是满足x(S,V−S)=2x(S,V-S)=2x(S,V−S)=2的SSS。
我们首先选择一个城市exterior，然后定义W=V−exteriorW=V-{exterior}W=V−exterior。
定义T为一颗PQ树，D(u)为u的所有叶子节点，B(T)为所有D(u)的集合。如果B(T)的所有元素都是tight set，那么我们称T和x相容。
例如下面的这棵PQ树就和3.1节的x相容（选择exterior为42）
3.3 使用shrinking寻找subtour constraints
shrinking指的是将集合SSS替换为点σ\sigmaσ，相对应的，边集合做如下替换：
如果x满足subtour constraints，那么x相容的PQ树的所有u节点构成的x[u]也满足subtour constraints。
3.4 分枝定界法概述
分枝定界法是指将问题拆分成多个子问题来求解的一种方法。
如果只用subtour约束，我们可以结合分枝定界算法来处理整数约束条件：
分枝定界的关键是如何快速确定下界，否则膨胀速度会非常快。因此将其与分割法结合起来，得到新的分枝切割法。
3.5 代码实现
Node算子用来存储搜索树的状态。其中level等于path的长度，path是当前节点已经访问过的vertex清单，bound则是当前的lb。
这里的bound函数是一种启发式方法，等于当前路径的总长度，再加上往后走两步的最小值。
struct Node
level::Int
path::Vector{Int64}
bound::Int
end
function totaldist(adj_mat::Array{Int64,2},t::Vector{Int64} )
n = length(t)
sum([adj_mat[t[i],t[i+1]] for i in 1:n-1])+adj_mat[t[n],t[1]]
end
function bound(adj_mat::Array{Int64,2}, path::Vector{Int64} )
_bound = 0
n = size(adj_mat)[1]
determined, last = path[1:end-1], path[end]
remain = setdiff(1:n,path)
for i in 1:length(path)-1;_bound += adj_mat[path[i],path[i + 1]];end
_bound += minimum([adj_mat[last,i] for i in remain])
p = [path[1];remain]
for r in remain
_bound+=minimum([adj_mat[r,i] for i in setdiff(p,r)])
end
return _bound
end;
这里用priorityQueue存储节点，用Queue也是一样的。
分枝条件为bound<ub，往下搜索所有没有探访过的节点，使用函数setdiff(1:n,v.path)。当然这里可以尝试将搜索范围缩小，比如仅搜索最近的一些节点，不过就不保证最优性了。
当搜索到level==n-1时，获得一个可行解，并且停止往下探索。此时如果路径长度比ub还短，则更新ub。
function solve(adj_mat::Array{Int64,2},ub::Int64 = 10^9)
optimal_tour = Vector{Int64}()
optimal_length = 0
n = size(adj_mat)[1]
PQ = PriorityQueue{Node,Int}()
path = Vector{Int64}([1])
v = Node(1,path,bound(adj_mat,path))
enqueue!(PQ,v,v.bound)
while length(PQ)>0
v = dequeue!(PQ)
if v.bound<ub
level = v.level+1
b = 0
for i in setdiff(1:n,v.path)
path = [v.path;i]
if level==n-1 #终止条件
push!(path,setdiff(1:n,path)[1])
_len = totaldist(adj_mat,path)
if _len < ub
ub = _len
optimal_length = _len
optimal_tour = path
end
else # 进行分叉
b = bound(adj_mat,path)
if b < ub # 分枝条件
enqueue!(PQ,Node(level,path,b),b)
end
end
end
end
end
optimal_tour,optimal_length
end
solve([0 14 4 10 20;14 0 7 8 7;4 5 0 7 16;11 7 9 0 2;18 7 17 4 0])
输出([1, 4, 5, 2, 3], 30)。
TSP时一个NPhard问题，当点数增多时，使用b&b的算法性能会急速下降。
3.6 另一种b&b算法
参考1980年Ton VOLGENANT and Roy JONKER的《A branch and bound algorithm for the symmetric traveling salesman probli m based on the 1-tree relaxation》，基于minimum-1-tree进行分枝定界。
首先是基于mst的ascent算法：
其中z(T)是mst的长度。我们的目标是迭代π\piπ得到最优的lowerbound，即这里的w(π)w(\pi)w(π)。
Held and Karp首先使用了ascent method，并且在Held, Wolfe and Crowder论文中进行了优化。Sraith and
Thompson开发了变种算法。令L为lower bound，U为upper bound，迭代算法为：
Held, Wolfe and Crowder证明了，当U等于optimal长度，并且p足够小时，上面的计算是收敛的。
由于上面的收敛速度比较慢，因此常用的迭代公式为：
令R为必须在结果中的边，F为不可以出现在结果中的边，则分枝策略为：
找到degree大于2的点p，然后选择其中的两条边e1和e2进行分枝：
至于upper bound，这里使用Christofides的简化版本：
紧接着使用Lin algorithm，等价于insertion+2opt算法。
测试结果如下：
4. 动态规划法
4.1 概述
定义c(s,k)c(s,k)c(s,k)为当前在kkk，待访问点的集合sss，最后返回城市0的最短路径，那么Bellman方程为：
c(s,k)=mini∈s{c(s−{i},i)+di,k}c(s,k)=\min_{i \in s}\{c(s-\{i\},i)+d_{i,k}\}c(s,k)=mini∈s{c(s−{i},i)+di,k}
为了使用方便，这里使用一个trick，即使用二进制来表达，用位运算符来计算，称作set bits：
- 左移和右移运算符可以快速计算2的幂：每左移一位，相当于该数乘以2；每右移一位，相当于该数除以2。因此，1 << k等价于2k2^k2k。假设S中包含k1,...,knk_1,...,k_nk1,...,kn，则我们可以将s等价替换位S=2k1+...+2knS=2^{k_1}+...+2^{k_n}S=2k1+...+2kn
- 按位或运算符|：可以用来计算集合的并集
- 按位与运算符&和取反运算符 ~：可以用来计算集合的差集
我们有初始状态c({0},k)=(d0,k,0)c(\{0\},k)=(d_{0,k},0)c({0},k)=(d0,k,0)。
逐步扩大set的尺寸，遍历所有可能的subset，使用bellman方程迭代计算。我们以n=5为例进行计算，初始状态为：
第一轮结束后，状态清单为：
第三轮结束后，状态清单为：
第四轮结束后，状态清单为：
4.2 Julia代码示例
using Combinatorics
function held_karp(dists::Array{Int64, 2})
n = size(dists,1)
C = Dict{Tuple{Int64, Int64},Tuple{Int64, Int64}}()
for k in 1:n-1
C[(1 << k, k)] = (dists[n,k], n)
end
for subset_size in 2:n-1
for subset in combinations(1:n-1,subset_size)
bits = 0
for bit in subset;bits |= 1 << bit;end
for k in subset
prev = bits & ~(1 << k)
res = Array{Tuple{Int64, Int64},1}()
for m in subset
if m == k;continue;end
push!(res,(C[(prev,m)][1]+dists[m,k],m))
end
C[(bits,k)] = minimum(res)
end
end
end
bits = 1<<n - 2
res = Array{Tuple{Int64, Int64},1}()
for k in 1:n-1
push!(res, (C[(bits, k)][1] + dists[k,n], k))
end
opt, parent = minimum(res)
path = []
for i in 1:n-1
push!(path,parent)
new_bits = bits & ~(1 << parent)
_, parent = C[(bits, parent)]
bits = new_bits
end
path
end
4.3 python代码示例
import itertools
import random
import sys
def generate_distances(n):
dists = [[0] * n for i in range(n)]
for i in range(n):
for j in range(i+1, n):
dists[i][j] = dists[j][i] = random.randint(1, 99)
return dists
def held_karp(dists):
"""
Implementation of Held-Karp, an algorithm that solves the Traveling
Salesman Problem using dynamic programming with memoization.
Parameters:
dists: distance matrix
Returns:
A tuple, (cost, path).
"""
n = len(dists)
# Maps each subset of the nodes to the cost to reach that subset, as well
# as what node it passed before reaching this subset.
# Node subsets are represented as set bits.
C = {}
# Set transition cost from initial state
for k in range(1, n):
C[(1 << k, k)] = (dists[0][k], 0)
# Iterate subsets of increasing length and store intermediate results
# in classic dynamic programming manner
for subset_size in range(2, n):
for subset in itertools.combinations(range(1, n), subset_size):
# Set bits for all nodes in this subset
bits = 0
for bit in subset:
bits |= 1 << bit
# Find the lowest cost to get to this subset
for k in subset:
prev = bits & ~(1 << k)
res = []
for m in subset:
if m == 0 or m == k: # 不允许回到0，不允许为当前点
continue
res.append((C[(prev, m)][0] + dists[m][k], m))
C[(bits, k)] = min(res)
# We're interested in all bits but the least significant (the start state)
bits = (2**n - 1) - 1
# Calculate optimal cost
res = []
for k in range(1, n):
res.append((C[(bits, k)][0] + dists[k][0], k))
opt, parent = min(res)
# Backtrack to find full path
path = []
for i in range(n - 1):
path.append(parent)
new_bits = bits & ~(1 << parent)
_, parent = C[(bits, parent)]
bits = new_bits
# Add implicit start state
path.append(0)
return opt, list(reversed(path))
def held_karp2(dists):
n = len(dists)
C = {}
for k in range(1, n):
C[(1, k)] = (dists[0][k], 0)
for subset_size in range(1, n):
for m in range(1, n):
for subset in itertools.combinations(list(range(1,m))+list(range(m+1,n)), subset_size):
bits = 1
for bit in subset:
bits |= 1 << bit
res = []
for k in subset: # 决策变量
next = bits & ~(1 << k)
res.append((C[(next, k)][0] + dists[k][m], k))
C[(bits, m)] = min(res)
# We're interested in all bits but the least significant (the start state)
bits = 2**n - 1
# Calculate optimal cost
res = []
for k in range(1, n):
res.append((C[(bits & ~(1 << k), k)][0] + dists[k][0], k))
opt, parent = min(res)
# Backtrack to find full path
path = []
for i in range(n-1):
path.append(parent)
bits = bits & ~(1 << parent)
_, parent = C[(bits, parent)]
# Add implicit start state
path.append(0)
return opt, list(reversed(path))
def held_karp3(dists):
n = len(dists)
C = {}
for k in range(1, n):
C[(frozenset([k]), k)] = (dists[0][k], 0)
for subset_size in range(2, n):
for subset in itertools.combinations(range(1, n), subset_size):
for k in subset:
prev = frozenset(set(subset) - {k})
res = []
for m in subset:
if m == 0 or m == k:
continue
res.append((C[(prev, m)][0] + dists[m][k], m))
C[(frozenset(subset), k)] = min(res)
# We're interested in all bits but the least significant (the start state)
bits = set(list(range(1,n)))
# Calculate optimal cost
res = []
for k in range(1, n):
res.append((C[(frozenset(bits), k)][0] + dists[k][0], k))
opt, parent = min(res)
# Backtrack to find full path
path = []
for i in range(n - 1):
path.append(parent)
new_bits = bits - {parent}
_, parent = C[(frozenset(bits), parent)]
bits = new_bits
# Add implicit start state
path.append(0)
return opt, list(reversed(path))
简化版如下：
# memoization for top down recursion
memo = [[-1]*(1 << (n+1)) for _ in range(n+1)]
def fun(i, mask):
# base case
# if only ith bit and 1st bit is set in our mask,
# it implies we have visited all other nodes already
if mask == ((1 << i) | 3):
return dist[1][i]
# memoization
if memo[i][mask] != -1:
return memo[i][mask]
res = 10**9 # result of this sub-problem
# we have to travel all nodes j in mask and end the path at ith node
# so for every node j in mask, recursively calculate cost of
# travelling all nodes in mask
# except i and then travel back from node j to node i taking
# the shortest path take the minimum of all possible j nodes
for j in range(1, n+1):
if (mask & (1 << j)) != 0 and j != i and j != 1:
res = min(res, fun(j, mask & (~(1 << i))) + dist[j][i])
memo[i][mask] = res # storing the minimum value
return res
# Driver program to test above logic
ans = 10**9
for i in range(1, n+1):
# try to go from node 1 visiting all nodes in between to i
# then return from i taking the shortest route to 1
ans = min(ans, fun(i, (1 << (n+1))-1) + dist[i][1])
print("The cost of most efficient tour = " + str(ans))
5. reduce matrix法
5.1. 概述
reduced matrix本质是一个branch&bound方法，其基本步骤为：
- 矩阵reduce的长度就是tsp问题的最优解
- 每一步，我们都要决定如果路径u到v采用的话，最小的总长度是多少。做法是：将第u排和第v列全部设置为infinity，然后reduce矩阵，将reduce出来的值和uv的长度驾到已经计算出来的最小路径长度上
- 如果路径上已访问完所有的点，则更新ub
5.2 例子说明
我们用一个例子来说明：
距离矩阵为：
我们用行表示从每个点出去的距离，用列表示从每个点进来的距离。
首先计算从每个点出去的最小值，如下：
将最小值从相应的边上减去后，我们再来看每个点进来的最小值：
我们得到reduce cost，可以作为tsp的lb = (10 + 10 + 15 + 20 + 5 + 10) = 70
接下来我们对1->2进行分枝，以上面的矩阵为基础，删除第1行和第2列：
继续分枝下去，分枝图为：
扫一扫
5626
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
