一种旅行商问题（TSP）具有高近似比的启发式解法的python实现

简介

文本实现了一个自研的TSP问题近似解法。在规模5～100的问题上近似比在1.1左右。另外在规模为1000的一个测试用例上测得近似比为：1.0742784198679693。时间复杂度O(n^3)。

基本思想：基于最近邻的思想，每次选择一个最近且加入后不会产生边交叉的顶点。
但是有可能不存在这样的顶点。所以当不存在时，选择最近的顶点，且把该顶点加入后，去除边交叉。
直到所有顶点都被选择为止。

注意：这里的边交叉是指：路径中存在一条子路径，可以通过翻转该子路径后使得原路径变短。如下图。

实验结果

本文对5～100规模的问题随机产生100个测试用例测得其近似比统计情况如图。横坐标表示问题的规模，即TSP中顶点的个数。纵坐标表示近似比。其中橙色线条表示的是中位数。可以发现，在规模20以下几乎能得到最优解。

代码

以下我代码复制到您的项目中可直接使用。

首先是定义一个工具类 SimpleGraph。命名为graph.py

from typing import *
import numpy as np

class SimpleGraph:
    """
    简单图结构
    """
    def __init__(self, size:int, bidirectional:bool = False):
        """
        构造函数
        :param size: 图的大小
        :param bidirectional: 是否是无向图，默认是有向图
        """
        assert size > 0
        self.size = size
        self.bidirectional = bidirectional
        # 邻接矩阵
        self.neighbourMat = np.zeros((size, size), dtype=np.int)
        # 权值矩阵
        self.weightMat = np.zeros((size, size), dtype=np.float)

    def getEdgeWeight(self, s:int, t:int)->float:
        """
        获取边的权重
        :param s: 源点
        :param t: 目标点
        :return: 权重值
        """
        return self.weightMat[s, t]

    def setEdgeWeight(self, s:int, t:int, weight:float):
        """
        设置边的权重
        :param s: 源点
        :param t: 目标点
        :param weight: 权重值
        :param bidirectional: 是否双向设置权重（仅仅适用于无向图）
        :return: 无
        """
        self.weightMat[s, t] = weight
        if(self.bidirectional):
            self.weightMat[t, s] = weight

    def getVerticeWeight(self, index:int)->float:
        """
        获取顶点权重
        :param index: 顶点编号
        :return: 权值
        """
        return self.weightMat[index, index]

    def setVerticeWeight(self, index:int, weight:float):
        """
        设置顶点权重
        :param index: 顶点编号
        :return: 无
        """
        self.weightMat[index, index] = weight

    def connect(self, s:int, t:int):
        """
        连接两个顶点
        :param s: 源点
        :param t: 目标点
        :return: 无
        """
        assert s != t
        self.neighbourMat[s, t] = 1
        if (self.bidirectional):
            self.neighbourMat[t, s] = 1

    def disconnect(self, s:int, t:int):
        """
        取消连接两个顶点
        :param s: 源点
        :param t: 目标点
        :return: 无
        """
        assert s != t
        self.neighbourMat[s, t] = 0
        if (self.bidirectional):
            self.neighbourMat[t, s] = 0

    def connected(self,  s:int, t:int)->bool:
        """
        判断两个顶点是否连接
        :param s: 源点
        :param t: 目标点
        :return: 是否连接
        """
        if not self.bidirectional:
            return s != t and self.neighbourMat[s, t] == 1
        else:
            return s != t and self.neighbourMat[s, t] == 1 and self.neighbourMat[t, s] == 1

    def getSubGraph(self, vertices:List[int]):
        """
        指定顶点获取其子图
        :param vertices: 顶点的索引
        :return: 子图
        """
        subgraph:SimpleGraph = SimpleGraph(len(vertices), self.bidirectional)
        for i in range(len(subgraph)):
            for j in range(len(subgraph)):
                subgraph.neighbourMat[i, j] = self.neighbourMat[vertices[i], vertices[j]]
                subgraph.weightMat[i, j] = self.weightMat[vertices[i], vertices[j]]
        return subgraph

    def __len__(self)->int:
        """
        图的大小
        :return: 大小
        """
        return self.size

然后是本文算法的代码：

from typing import *
from graph import SimpleGraph
import numpy as np

def swap(ls:List[int], fromIndex:int, toIndex:int):
    """
    列表交换函数
    :param ls: 要交换的列表
    :param fromIndex: 要交换的开始位置
    :param toIndex: 要交换的结束位置
    :return: 无
    """
    while fromIndex < toIndex:
        t = ls[fromIndex]
        ls[fromIndex] = ls[toIndex]
        ls[toIndex] = t
        fromIndex += 1
        toIndex -= 1
        
def dot2graph(dots:List[Tuple[float, float]])->SimpleGraph:
    """
    把二维平面的点转换成无向简单图
    :param dots: 二维平面的点
    :return: 简单图
    """
    graph: SimpleGraph = SimpleGraph(len(dots), bidirectional=True)
    for i in range(len(dots)):
        for j in range(i + 1, len(dots)):
            graph.connect(i, j)
            dx = dots[i][0] - dots[j][0]
            dy = dots[i][1] - dots[j][1]
            weight = math.sqrt(dx ** 2 + dy ** 2)
            graph.setEdgeWeight(i, j, weight)
    return graph

def optimize(graph:SimpleGraph, circle:List[int]):
    """
    对已有回圈的优化, 根据已有的回路去除回路中存在的边交叉
    :param graph: 简单图
    :param circle: 回圈
    :return: 无
    """
    size = len(graph)
    count = 0
    while True:
        count += 1
        existsCross = False
        for i in range(size):
            for j in range(size):
                # i , j指的都是顶点的编号，同时指代它跟着的那条边
                # 😊 ---i--> 😊 ---....--- 😊 --j--> 😊
                if i != j:
                    dot1 = i
                    dot2 = (i + 1) % size
                    dot3 = j
                    dot4 = (j + 1) % size
                    curCost = graph.getEdgeWeight(circle[dot1], circle[dot2])\
                              + graph.getEdgeWeight(circle[dot3], circle[dot4])
                    newCost = graph.getEdgeWeight(circle[dot1], circle[dot3])\
                              + graph.getEdgeWeight(circle[dot2], circle[dot4])
                    if newCost < curCost:
                        existsCross = True
                        if dot2 < dot3:
                            swap(circle, dot2, dot3)
                        else:
                            swap(circle, dot2 - size, dot3)
                        break # 如果加了这个break，效果会差些，时间复杂度降到O(n^3)左右，实际测试为O(n^2)。不加两者都会增加一个量级
        if not existsCross:
            break

    print(count)


def tsp(dots:List[Tuple[float, float]])->List[int]:
    """
    启发式方法，选择没有交叉的结点
    :param dots: 一系列的点的坐标，点之间的距离表示代价
    :return: 一系列点的编号，代表得到的哈密顿环
    """
    #构造简单图
    graph:SimpleGraph = dot2graph(dots)

    # 开始贪心地搜索
    path =  [0]
    remained = [i for i in range(1, len(graph))]
    while len(remained) > 0:
        minCost = float("inf")
        minIndex = -1
        minNoCrossCost = float("inf")
        minNoCrossIndex = -1
        for j in range(len(remained)) :
            cost = graph.getEdgeWeight(path[-1], remained[j])
            if cost < minCost:
                minCost = cost
                minIndex = j
            # 检测cross
            detectedCross = False
            for k in range(0, len(path) - 2):
                currentTwoEdgeCost = graph.getEdgeWeight(path[k], path[k+1]) + graph.getEdgeWeight(path[-1], remained[j])
                newTwoEdgeCost = graph.getEdgeWeight(path[k], path[-1]) + graph.getEdgeWeight(path[k+1], remained[j])
                if newTwoEdgeCost < currentTwoEdgeCost:
                    detectedCross = True
                    break
            if not detectedCross and cost < minNoCrossCost:
                minNoCrossCost = cost
                minNoCrossIndex = j

        if minNoCrossIndex != -1:
            path.append(remained[minNoCrossIndex])
            del remained[minNoCrossIndex]
        else:
            path.append(remained[minIndex])
            del remained[minIndex]

            # 清理交叉(理论复杂度O(n^3)， 实际运行复杂度O(n^2))
            while True:
                existsCross = False
                for i in range(0, len(path) - 3):
                    j = len(path) - 1
                    while i <= j - 2:
                        curCost = graph.getEdgeWeight(path[i], path[i+1]) + graph.getEdgeWeight(path[j-1], path[j])
                        newCost = graph.getEdgeWeight(path[i], path[j-1]) + graph.getEdgeWeight(path[i+1], path[j])
                        if newCost < curCost:
                            existsCross = True
                            swap(path, i+1, j-1)
                            break
                        else:
                            j -= 1
                if not existsCross:
                    break
    # 最后再优化一下
    optimize(graph, path)
    return path

使用示例

下面是一个规模为20的使用示例，采用随机产生的点。算法实际可以适用于对称TSP。

    def test_tsp3(self):
        dots = []
        xs = []
        ys = []
        for i in range(20):
            x = random()
            y = random()
            dots.append((x, y))
            xs.append(x)
            ys.append(y)
        path = tsp(dots)  #在这里调用一下
        path.append(path[0])
        plt.scatter(np.array(xs)[path], np.array(ys)[path])
        for i in range(len(xs)):
            plt.annotate(str(path[i]),
                         xy=(xs[path[i]], ys[path[i]]),
                         xytext=(xs[path[i]] + 0.01, ys[path[i]] + 0.01))
            # 这里xy是需要标记的坐标
        plt.plot(np.array(xs)[path], np.array(ys)[path])
        plt.show()

运行结果：