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

Posted 2023-03-10 IE06

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=1
    m2_ind=1
    for i in [1:root-1;root+1:n]
        if arr[i]<m1
            m2=m1
            m1=arr[i]
            m2_ind=m1_ind
            m1_ind=i
        elseif arr[i]<m2
            m2=arr[i]
            m2_ind=i
        end
    end
    return 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:n
        if x[i]==1
            append!(leaves,i)
        end
    end
    max_w = 0
    max_m1 = 1
    max_m2 = 1
    for leaf_ind in 1:length(leaves)
        leaf = leaves[leaf_ind]
        _,m2_ind,_,m2 = minimum_two(distmat[leaf,:],leaf)
        w = c+m2
        if w > max_w
            max_w = w
            max_m1 = leaf
            max_m2 = m2_ind
        end
    end
    max_v = [x[i] for i in 1:size(distmat)[1]].-2
    max_v[max_m1]+=1
    max_v[max_m2]+=1
    return 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]=2
    x[m1_ind]+=1
    x[m2_ind]+=1
    w = c+m1+m2-2*sum(pi)
    v = [x[i] for i in 1:n].-2
    distmat = distmat.-pi'.-pi
    return w,v
end

2.3 梯度下降代码

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

function ascent(c)
    n = size(c)[1]
    pi = zeros(n) # 初始化优化参数pi
    best_pi = pi
    t = 1 # 优化步长
    best_w,v,m1,m2 = minimum1tree(c, pi) # 初始化w和v
    best_deg = sum(v.*v) # 初始化目标函数
    last_v = v 
    period = max(floor(Int,n/2), 100)
    initial_period = period
    initialPhase = true
    while t > 0.01
        p = 1
        while p<=period
            for i in 1:n;if v[i] == 0;last_v[i] = 0; end;end
            pi = pi+t*(0.7*v+0.3*last_v)
            last_v = v
            w, 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, pi
            elseif (w > best_w && deg <= best_deg)
                @info("** T = $t, Period = $period, BestW = $w, Norm = $deg, m1 = $m1, m2 = $m2")
                best_w = w
                best_pi = pi
                best_deg = deg
                if 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 = false
                p = 1
                t = t*3/4 # 第一阶段过后，开始逐步收缩步长
            end
            p+=1
        end
        t = t/2 # 每个阶段迭代完成后，都收缩步长和迭代次数进行下轮迭代
        period = floor(Int,period/2)
    end
    return 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

此时结果如下：

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