运筹系列68：julia启发式求解tsp问题

1. 基础函数

1.1 读取数据和绘图

将读取数据/绘图等库安装好：

using Gadfly
using LinearAlgebra
using TravelingSalesmanHeuristics
using TSPLIB
tsp = readTSPLIB(:a280)
n = tsp.dimension
plot_instance(tsp) = plot(x = tsp.nodes[:,1], y = tsp.nodes[:,2], Geom.point, Guide.xlabel(nothing), Guide.ylabel(nothing))
function plot_solution(tsp, path, extras = [])
    ptspath = tsp.nodes[path,:]
    plot(x = ptspath[:,1], y = ptspath[:,2], Geom.point, Geom.path, Guide.xlabel(nothing), Guide.ylabel(nothing), extras...)
end
distmat = [norm(tsp.nodes[i,:] - tsp.nodes[j,:]) for i in 1:n, j in 1:n]

还有一个是树的绘制，代码如下：

for i in 1:length(mst)
    l = tsp.nodes[collect(mst[i]),:]
    PyPlot.plot(l[:,1], l[:,2], color="r",linewidth = 0.5, linestyle="--")
end

1.2 求解调用

第一种调用方法是使用quality_factor：

@time path, cost = solve_tsp(distmat; quality_factor = 85)

第二种是指定构造型启发式算法：

path_nn, cost_nn = nearest_neighbor(distmat; firstcity = 1, do2opt = false)
path_nn2opt, cost_nn2opt = nearest_neighbor(distmat; firstcity = 1, do2opt = true)
path_fi, cost_fi = farthest_insertion(distmat; firstcity = 1, do2opt = false)
path_fi2opt, cost_fi2opt = farthest_insertion(distmat; firstcity = 1, do2opt = true) 

第三种是指定改进型启发式算法：

path_sa, cost_sa = simulated_annealing(distmat; init_path = path_nn, num_starts = 10)

2. 获取下界

可以调用方法如下：

TravelingSalesmanHeuristics.vertwise_bound(distmat)
TravelingSalesmanHeuristics.hkinspired_bound(distmat)

2.1 vertwise下界

第一种下界非常简单，每个点取最短的两条边加起来即可。

# the cost of a tour must be >= the sum over all vertices of
# the cost of the cheapest edge leaving that vertex
# likewise for the cheapest edge entering that vertex
# since we must go to and leave each vertex
function vertwise_bound(distmat::AbstractMatrix{T}) where {T<:Real}
    # the simple code below would tend to pick out the 0 costs on the diagonal
    # so make a doctored copy of the distance matrix with high costs on the diagonal
    m = maximum(distmat)
    distmat_nodiag = distmat + m * I
    leaving = sum(minimum(distmat_nodiag, dims = 2))
    entering = sum(minimum(distmat_nodiag, dims = 1))
    return maximum([leaving, entering])
end

2.2 1-tree下界

这里要用到最小生成树：

# 最小生成树
# returns a (n-1) long Vector of Tuple{Int, Int} where each tuple is an edge in the MST
# and the total weight of the tree
# the matrix passed in must be symmetric or you won't get out the minimum spanning tree
function minspantree(dm::AbstractMatrix{T}) where {T<:Real} # accepts views
    mst_edges = Vector{Tuple{Int, Int}}()
    mst_cost = zero(T)
    n = size(dm, 1)

    # we keep a running list of the distance from each vertex to the partly formed tree
    # rather than using 0 for vertices already in the tree, we use a large value so that we
    # can find the closest non-tree vertex via call to Julia's `findmin`.
    bigval = maximum(dm) + one(T)
    tree_dists = dm[1,:] # distance to tree
    closest_tree_verts = ones(Int, n