使用蚁群算法求解 TSP 问题

发表于 29 Jul 2020 总计 6010字阅读时间 21分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

导言

问题背景

旅行商问题，即 TSP 问题（Traveling Salesman Problem）又译为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访 n 个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。等价于求图的最短哈密尔顿回路问题。

求解思路

蚁群算法

蚁群算法（Ant Colony Algorithm）最初于 1992 年由意大利学者 M.Dorigo 等人提出，它是一种模拟自然界中真实蚁群觅食行为的仿生优化算法。研究发现：每只蚂蚁觅食时在走过的路线上会留下一种称为信息素的物质，蚂蚁之间靠感知这种物质的浓度进行信息传递。蚂蚁在选择路径时总是倾向于朝信息索浓度高的方向移动，而距离短的路径上走过的蚂蚁多，留下的信息素也多，后续蚂蚁选择它的概率也会越大；其他路径上的信息素会随着时间的推移不断挥发，这样就形成了一种正反馈机制，最后整个蚁群聚集到最短路径上。

人工蚁群算法模拟了这一过程。每只蚂蚁在解空间独立地搜索可行解，解越好留下的信息素越多，随着算法推进，较优解路径上的信息素增多，选择它的蚂蚁也随之增多，最终收敛到最优或近似最优的解上。求解 TSP 问题的蚁群算法中的人工蚂蚁具有以下特点：

概率性地选择下一条路径，该概率与路径长度和路径上的信息素浓度有关；
为了保证解的逻辑可行，蚂蚁不允许选择已经走过的路径；
蚂蚁走过一条路径时会在该路径上面分泌一种叫做信息素的物质。

蚂蚁系统采用 $\tau_{i,j}(t)$ 来模仿 $t$ 时刻路径 $i\to j$ 上面的信息残留量，即信息素浓度。类似于蚂蚁觅食过程，每条路径上面的信息素会挥发，如果有蚂蚁经过的时候，信息素的浓度会相应增加。因此，蚂蚁系统中的信息素浓度的更新公式为：

\[\tau_{i,j}(t+n)=\rho\cdot\tau_{i,j}(t)+\Delta\tau_{i,j}\]

上式中 $\rho$ 是一轮遍历后信息素的保留量，$1-\rho$ 即为挥发因子；$\Delta\tau_{i,j}$ 表示表示一次旅行（遍历完所有城市）后，所有路径 i 到 j 的蚂蚁留下的信息素总量，即

\[\Delta\tau_{i,j}=\sum_k\tau_{i,j}^k\]

上式中 $\tau_{i,j}^k$ 表示第 k 只蚂蚁在路径 $i\to j$ 上面留下的信息素量。有如下三种生成模型：

蚁量模型： $\tau_{i,j}^k=\begin{cases} \frac{Q}{d_{i,j}},\,\text{蚂蚁}\,k\,\text{经过路径}\,i\to j \\ 0,\,\text{else} \end{cases}$
蚁周模型： $\tau_{i,j}^k=\begin{cases} \frac{Q}{L_k},\, \text{蚂蚁}\,k\,\text{经过路径}\,i\to j \\ 0 ,\,\text{else} \end{cases}$
蚁密模型： $\tau_{i,j}^k=\begin{cases} Q ,\,\text{蚂蚁}\,k\,\text{经过路径}\,i\to j \\ 0 ,\,\text{else} \end{cases}$

式中，$Q$ 代表蚂蚁单次分泌的信息素数量，是常数；$L_k$ 为蚂蚁已经走过路径的总长度。可以看出，蚁量模型使用到了当前路径的信息，蚂周模型使用了整体路径的信息，蚁密模型没有用到路径长度的信息。本文使用蚁量模型。

一般来说有了信息素浓度的更新公式，就可以直接给出蚂蚁对每条路径的选择概率了。然而，为了更好的利用 TSP 问题自身的性质，M.Dorigo 等引入了一个启发项：$\eta_{i,j}=\frac{1}{d_{i,j}}$ 。通过结合信息素浓度和启发因子，可以得到蚂蚁选择路径 i 到 j 的概率为：

\[p_{i,j}^k(t)= \begin{cases} \frac{\left[\tau_{i,j}(t)\right]^\alpha\left[\eta_{i,j}\right]^\beta} {\sum_{k\in \text{allowed}}\left[\tau_{i,k}(t)\right]^\alpha\left[\eta_{i,k}\right]^\beta} ,\,j\in \text{allowed} \\ 0,\,\text{else} \end{cases}\]

实验过程

实验环境

所用机器型号为 VAIO Z Flip 2016。

Intel(R) Core(TM) i7-6567U CPU @3.30GHZ 3.31GHz
8.00GB RAM
Windows 10 2004 19041.264, 64-bit
- Visual Studio Code 1.45.1
  - Python 2020.5.80290：去年九月底发布的 VSCode Python 插件支持在编辑器窗口内原生运行 juyter nootbook 了，非常赞！
  - Remote - WSL 0.44.2：配合 WSL，在 Windows 上获得 Linux 接近原生环境的体验。
- Windows Subsystem for Linux [Ubuntu 20.04 LTS]：WSL 是以软件的形式运行在 Windows 下的 Linux 子系统，是近些年微软推出来的新工具，可以在 Windows 系统上原生运行 Linux。
  - Python 3.8.2 64-bit：安装在 WSL 中。
    - jupyter==1.0.0
    - numpy==1.18.4
    - matplotlib==3.2.1

import time
import io
import numpy
import random
from matplotlib import pyplot

下载并获得数据集`eil51.tsp`

选择的问题和数据集是 eil51。

curl -O http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/tsp/eil51.tsp.gz
gzip -d eil51.tsp.gz

根据 http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/STSP.html 这个问题的最优解是 426。

预处理工作

此处完成了如下的工作：

从文件中读取各顶点的二维坐标。
预处理原图的邻接矩阵。这是由于，求二维点对的欧式距离需要引入比较耗时的平方运算，对其预处理之后可以大大减少程序运行时间。
对于给定的求解方法，多次求解并输出若干次运行后所得的最好解、平均值、标准差等等，并对求解过程和结果进行可视化。

def get_coords(s):
    with io.open(s, 'r') as f:
        coords = []
        coords_flag = 0
        for s in f.readlines():
            if s.strip() == 'EOF':
                coords_flag = 0
            if coords_flag:
                t = s.split()
                coords.append([int(t[1]), int(t[2])])
            if s.strip() == 'NODE_COORD_SECTION':
                coords_flag = 1
        return coords


def cal_length(solution, adjmap):
    length = 0
    for i in range(len(solution)):
        length += int(adjmap[solution[i-1]][solution[i]])
    return length


def tsp(method, test_time=10, coords=get_coords('eil101.tsp')):
    adjmap = [[numpy.linalg.norm(numpy.array(ci)-numpy.array(cj))
               for cj in coords]for ci in coords]
    length = []
    solution = []
    run_time = []
    ans = []
    for i in range(test_time):
        start = time.perf_counter()
        solution1, length1 = method(adjmap)
        run_time.append(time.perf_counter()-start)
        ans.append(length1[-1])
        if i == 0 or length[-1] > length1[-1]:
            length = length1
            solution = solution1

    print('以下是各次运行时间（s）及求解时间的均值、方差、标准差：')
    print(run_time)
    print(numpy.mean(run_time))
    print(numpy.var(run_time))
    print(numpy.std(run_time))

    print('以下是各次运行求得解的大小及其均值、方差、标准差：')
    print(ans)
    print(numpy.mean(ans))
    print(numpy.var(ans))
    print(numpy.std(ans))
    print(solution)

    print('以下是所求最优解及其可视化：')
    print(length[-1])
    x = [coords[i][0]for i in solution]
    y = [coords[i][1]for i in solution]
    x.append(x[0])
    y.append(y[0])
    pyplot.plot(x, y)
    pyplot.show()

    print('以下是最优解的收敛过程：')
    pyplot.plot(length)
    pyplot.show()

蚁群算法

# 蚂蚁个数，信息素重要程度因子，启发函数重要程度因子，信息素的挥发速度
def ag(adjmap, itermax=300, numant=10, alpha=1, beta=3, rho=0.1, q=1):
    numcity = len(adjmap)  # 城市个数
    etatable = numpy.array(
        [[1/(adjmap[i][j]+1) if i != j else 0 for j in range(numcity)]for i in range(numcity)])
    # 启发函数矩阵，表示蚂蚁从城市i转移到矩阵j的期望程度
    pheromonetable = numpy.ones((numcity, numcity))  # 信息素矩阵

    pathbest = []
    lengthbest = []

    for _iter in range(itermax):
        changepheromonetable = numpy.zeros((numcity, numcity))
        for i in range(numant):
            solution = [random.choice(range(numcity))]
            unvisited = set(range(numcity))  # 未访问标记
            length = 0
            for j in range(numcity - 1):
                unvisited.remove(solution[j])
                listunvisited = list(unvisited)
                solution.append(random.choices(listunvisited, k=1, weights=[numpy.power(
                    pheromonetable[solution[j]][k], alpha) * numpy.power(etatable[solution[j]][k], beta) for k in listunvisited])[0])
                # 每次用轮盘法选择下一个要访问的城市
            length = cal_length(solution, adjmap)
            if len(lengthbest) == 0 or lengthbest[-1] > length:
                lengthbest.append(length)
                pathbest = solution
            else:
                lengthbest.append(lengthbest[-1])
            for j in range(numcity - 1):
                changepheromonetable[solution[j-1]][solution[j]] += q / length
        pheromonetable = (1 - rho) * pheromonetable + \
            changepheromonetable  # 更新信息素
    return pathbest, lengthbest


tsp(ag, coords=get_coords('eil51.tsp'))

结果分析

蚁群算法

以下是各次运行时间（s）及求解时间的均值、方差、标准差：
[25.543905400000085, 26.20327139999995, 26.227915800000005, 25.777282900000046, 25.608007600000064, 25.970474299999978, 25.994625400000018, 24.738855700000045, 25.376835500000084, 31.898372999999992]
26.333954700000028
3.6159814329762767
1.9015734098309949
以下是各次运行求得解的大小及其均值、方差、标准差：
[446, 427, 447, 426, 439, 443, 446, 444, 436, 441]
439.5
52.65
7.256031973468694
[3, 46, 11, 45, 50, 26, 0, 31, 10, 37, 4, 48, 8, 49, 15, 1, 21, 7, 25, 30, 27, 2, 19, 34, 35, 28, 20, 33, 29, 9, 38, 32, 44, 14, 43, 36, 16, 41, 39, 18, 40, 12, 24, 13, 23, 42, 6, 22, 47, 5, 17]
以下是所求最优解及其可视化：
426

可视化

以下是最优解的收敛过程：

收敛过程

从时间上来看，由于蚁群算法生成一个新解的时间更慢（要通过禁忌表判断点是否在路径上；然后基于轮盘赌的方式选择下一个要访问的点），因此运行时间比之前实验中做的模拟退火和遗传算法都要慢了很多。蚂蚁容易聚集在信息素浓度较高的地方，这导致可以看到蚁群算法很容易收敛到一个局部最优解上，运行十次只有一次得到了最优解。

使用蚁群算法求解 TSP 问题

导言

问题背景

求解思路

蚁群算法

实验过程

实验环境

下载并获得数据集eil51.tsp

预处理工作

蚁群算法

结果分析

蚁群算法

主要参考文献

下载并获得数据集`eil51.tsp`