运筹系列68：TSP问题Held-Karp下界的julia实现

1. 介绍

Held-Karp下界基于1tree下界，但是增加了点权重，如下图

通过梯度下降的方法找到最优的$\pi$。
这里用到的1tree有下面几种：

1. 全部点用来生成最小生成树，再加上所有叶子结点第二短的边中数值最大的那个
2. 任意选一个点，选取它最短的两边边；然后剩下的点生成最小生成树
3. 和2类似，但是枚举所有的点。

2. 代码分析

首先是根据距离列表arr，获得当前节点root最近的两条边

function minimum_two(arr,root)n = length(arr)m1=arr[1]m2=arr[1]m1_ind=1m2_ind=1for i in [1:root-1;root+1:n]if arr[i]<m1m2=m1m1=arr[i]m2_ind=m1_indm1_ind=ielseif arr[i]<m2m2=arr[i]m2_ind=iendendreturn m1_ind,m2_ind,m1,m2
end

2.1 第一种1tree

function minimum1tree(distmat,pi)distmat = distmat.+pi.+pi'mst, c = TravelingSalesmanHeuristics.minspantree(distmat)x = counter(cat([m[1] for m in mst],[m[2] for m in mst],dims=1))n = size(distmat)[1]leaves = []for i in 1:nif x[i]==1append!(leaves,i)endendmax_w = 0max_m1 = 1max_m2 = 1for leaf_ind in 1:length(leaves)leaf = leaves[leaf_ind]_,m2_ind,_,m2 = minimum_two(distmat[leaf,:],leaf)w = c+m2if w > max_wmax_w = wmax_m1 = leafmax_m2 = m2_indendendmax_v = [x[i] for i in 1:size(distmat)[1]].-2max_v[max_m1]+=1max_v[max_m2]+=1return max_w-2*sum(pi),max_v,max_m1,max_m2
end

2.2 第二种1tree

function minimum1tree(distmat,pi,first_node)distmat = distmat.+pi.+pi'm1_ind,m2_ind,m1,m2 = minimum_two(distmat[first_node,:],first_node)n = size(distmat)[1]left_nodes = [1:first_node-1;first_node+1:n]mst, c = TravelingSalesmanHeuristics.minspantree(distmat[left_nodes,left_nodes])x = counter(cat([left_nodes[m[1]] for m in mst],[left_nodes[m[2]] for m in mst],dims=1))x[first_node]=2x[m1_ind]+=1x[m2_ind]+=1w = c+m1+m2-2*sum(pi)v = [x[i] for i in 1:n].-2distmat = distmat.-pi'.-pireturn w,v
end

2.3 梯度下降代码

这里使用的是第一种minimum1tree，第二种类似。

function ascent(c)n = size(c)[1]pi = zeros(n) # 初始化优化参数pibest_pi = pit = 1 # 优化步长best_w,v,m1,m2 = minimum1tree(c, pi) # 初始化w和vbest_deg = sum(v.*v) # 初始化目标函数last_v = v period = max(floor(Int,n/2), 100)initial_period = periodinitialPhase = truewhile t > 0.01p = 1while p<=periodfor i in 1:n;if v[i] == 0;last_v[i] = 0; end;endpi = pi+t*(0.7*v+0.3*last_v)last_v = vw, v, m1, m2 = minimum1tree(c, pi) deg = sum(v.*v)if deg == 0@info("* T = $t, Period = $period, BestW = $w, Norm = $deg, m1 = $m1, m2 = $m2")return w, pielseif (w > best_w && deg <= best_deg)@info("** T = $t, Period = $period, BestW = $w, Norm = $deg, m1 = $m1, m2 = $m2")best_w = wbest_pi = pibest_deg = degif initialPhase;t *= 2;end # 增加步长if p == period;period*=2;end # 增加迭代次数elseif initialPhase && p > initial_period /2@info("* T = $t, Period = $period, BestW = $w, Norm = $deg")initialPhase = falsep = 1t = t*3/4 # 第一阶段过后，开始逐步收缩步长endp+=1endt = t/2 # 每个阶段迭代完成后，都收缩步长和迭代次数进行下轮迭代period = floor(Int,period/2)endreturn best_w, best_pi
end

3. 实测结果

我们使用TSPLIB的att48进行观测，最优解为33522.0。梯度下降打印信息如下：

[ Info: ** T = 1, Period = 100, BestW = 29266.0, Norm = 34, m1 = 2, m2 = 42
[ Info: ** T = 2, Period = 100, BestW = 29333.000000000004, Norm = 34, m1 = 2, m2 = 42
[ Info: ** T = 4, Period = 100, BestW = 29463.4, Norm = 32, m1 = 2, m2 = 42
[ Info: ** T = 8, Period = 100, BestW = 29954.6, Norm = 32, m1 = 2, m2 = 42
[ Info: ** T = 16, Period = 100, BestW = 30398.2, Norm = 26, m1 = 2, m2 = 42
[ Info: ** T = 32, Period = 100, BestW = 31464.399999999998, Norm = 18, m1 = 2, m2 = 42
[ Info: ** T = 64, Period = 100, BestW = 31744.999999999993, Norm = 18, m1 = 2, m2 = 42
[ Info: ** T = 128, Period = 100, BestW = 32487.200000000004, Norm = 18, m1 = 29, m2 = 5
[ Info: * T = 256, Period = 100, BestW = 28475.4, Norm = 62
[ Info: ** T = 96.0, Period = 50, BestW = 32514.400000000012, Norm = 18, m1 = 5, m2 = 29
[ Info: ** T = 96.0, Period = 50, BestW = 32873.40000000001, Norm = 10, m1 = 5, m2 = 29
[ Info: ** T = 96.0, Period = 50, BestW = 33086.20000000001, Norm = 10, m1 = 5, m2 = 29
[ Info: ** T = 96.0, Period = 50, BestW = 33115.8, Norm = 8, m1 = 29, m2 = 5
[ Info: ** T = 48.0, Period = 25, BestW = 33168.00000000001, Norm = 8, m1 = 29, m2 = 34
[ Info: ** T = 48.0, Period = 25, BestW = 33292.399999999994, Norm = 6, m1 = 29, m2 = 34
[ Info: ** T = 24.0, Period = 12, BestW = 33391.399999999994, Norm = 6, m1 = 29, m2 = 34
[ Info: ** T = 24.0, Period = 12, BestW = 33410.6, Norm = 6, m1 = 29, m2 = 34
[ Info: ** T = 12.0, Period = 6, BestW = 33411.200000000004, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 12.0, Period = 12, BestW = 33417.19999999999, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 12.0, Period = 12, BestW = 33421.200000000004, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 12.0, Period = 24, BestW = 33423.6, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 12.0, Period = 24, BestW = 33424.799999999996, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 12.0, Period = 24, BestW = 33424.80000000001, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 12.0, Period = 24, BestW = 33435.200000000004, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 6.0, Period = 12, BestW = 33441.00000000001, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 3.0, Period = 6, BestW = 33442.399999999994, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 3.0, Period = 6, BestW = 33443.4, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 3.0, Period = 12, BestW = 33444.299999999996, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 1.5, Period = 12, BestW = 33444.9, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 1.5, Period = 12, BestW = 33445.350000000006, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 1.5, Period = 12, BestW = 33446.59999999999, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 0.375, Period = 3, BestW = 33447.31249999999, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 0.1875, Period = 3, BestW = 33447.55625, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 0.1875, Period = 6, BestW = 33447.706249999996, Norm = 2, m1 = 29, m2 = 34
[ Info: ** T = 0.09375, Period = 3, BestW = 33447.825, Norm = 2, m1 = 29, m2 = 34

此时结果如下：

最优结果如下，已非常接近。