运筹系列65:使用Julia精确求解tsp问题

文章介绍了旅行商问题(TSP)的解决方法,包括Christofides算法构建的上界和使用Blossom算法的示例,以及通过松弛问题和子回路约束寻找下界。此外,讨论了如何通过添加约束获取整数解,如polyhedralcombinatorics和分支定界法,并提到了PQ树和shrinking在寻找subtourconstraints中的应用。
摘要由CSDN通过智能技术生成

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的松弛问题

首先看一个例子:
在这里插入图片描述
我们定义 x i j x_{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)

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 ) = 2 x(S,V-S)=2 x(S,VS)=2 S S S
我们首先选择一个城市exterior,然后定义 W = V − e x t e r i o r W=V-{exterior} W=Vexterior
定义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指的是将集合 S S S替换为点 σ \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)为当前在 k k k,待访问点的集合 s s s,最后返回城市0的最短路径,那么Bellman方程为:
c ( s , k ) = min ⁡ i ∈ s { c ( s − { i } , i ) + d i , k } c(s,k)=\min_{i \in s}\{c(s-\{i\},i)+d_{i,k}\} c(s,k)=minis{c(s{i},i)+di,k}
为了使用方便,这里使用一个trick,即使用二进制来表达,用位运算符来计算,称作set bits:

  1. 左移和右移运算符可以快速计算2的幂:每左移一位,相当于该数乘以2;每右移一位,相当于该数除以2。因此,1 << k等价于 2 k 2^k 2k。假设S中包含 k 1 , . . . , k n k_1,...,k_n k1,...,kn,则我们可以将s等价替换位 S = 2 k 1 + . . . + 2 k n S=2^{k_1}+...+2^{k_n} S=2k1+...+2kn
  2. 按位或运算符|:可以用来计算集合的并集
  3. 按位与运算符&和取反运算符 ~:可以用来计算集合的差集

我们有初始状态 c ( { 0 } , k ) = ( d 0 , 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方法,其基本步骤为:

  1. 矩阵reduce的长度就是tsp问题的最优解
  2. 每一步,我们都要决定如果路径u到v采用的话,最小的总长度是多少。做法是:将第u排和第v列全部设置为infinity,然后reduce矩阵,将reduce出来的值和uv的长度驾到已经计算出来的最小路径长度上
  3. 如果路径上已访问完所有的点,则更新ub

5.2 例子说明

我们用一个例子来说明:
在这里插入图片描述
距离矩阵为:
在这里插入图片描述
我们用行表示从每个点出去的距离,用列表示从每个点进来的距离。
首先计算从每个点出去的最小值,如下:
在这里插入图片描述
将最小值从相应的边上减去后,我们再来看每个点进来的最小值:
在这里插入图片描述
我们得到reduce cost,可以作为tsp的lb = (10 + 10 + 15 + 20 + 5 + 10) = 70
接下来我们对1->2进行分枝,以上面的矩阵为基础,删除第1行和第2列:
在这里插入图片描述
在这里插入图片描述

继续分枝下去,分枝图为:
在这里插入图片描述

  • 2
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值