使用模拟退火和遗传算法求解 TSP 问题

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

---

摘要： 本文分别使用模拟退火算法和遗传算法解决 TSP 问题，TSP 问题的规模大小为 101 个城市。实验中采用多种邻域操作的局部搜索 local search 策略求解；并在在局部搜索策略的基础上，加入模拟退火 simulated annealing 策略，并比较两者的效果；求得的解不超过最优值的 10％，并提供可视化，便于观察路径的变化和交叉程度。随后使用遗传算法求解，设计较好的交叉操作，并且引入相同的局部搜索操作；和之前的模拟退火算法进行比较，得出设计高效遗传算法的一些经验，并比较单点搜索和多点搜索的优缺点。最终通过该实验得出结论，局部搜索容易陷入局部最优解，而模拟退火相较于局部搜索较少依赖于初始解的选择而具有更好的鲁棒性；遗传算法在本实验中的小算例上表现不如前两者，但是在更大规模的数据集上值得期待。

导言

问题背景

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

求解思路

局部搜索

局部搜索是解决最优化问题的一种启发式算法。对于某些计算起来非常复杂的最优化问题，比如各种 NP 完全问题，要找到最优解需要的时间随问题规模呈指数增长，因此诞生了各种启发式算法来退而求其次寻找次优解，属于近似算法（Approximate algorithms），以精度换时间。局部搜索就是其中的一种方法。

局部搜索算法从一个初始解开始，通过邻域动作，产生其邻居解，判断邻居解的质量，根据某种策略，来选择邻居解，重复上述过程，至到达终止条件。

模拟退火

模拟退火法（Simulated Annealing）是克服爬山法容易陷入局部最优解的缺点的有效方法。所谓退火是冶金专家为了达到某些特种晶体结构重复将金属加热或冷却的过程，该过程由系统的温度进行控制。

模拟退火法的基本思想是，在系统朝着能量减小的趋势这样一个变化过程中，偶尔允许系统跳到能量较高的状态，以避开局部极小点，最终稳定到全局最小点。

遗传算法

遗传算法是解决搜索问题的一种通用算法，对于各种通用问题都可以使用。在遗传算法中，这几个特征以一种特殊的方式组合在一起：基于染色体群的并行搜索，带有猜测性质的选择操作、交换操作和突变操作。这种特殊的组合方式将遗传算法与其它搜索算法区别开来。

遗传算法通常首先组成一组候选解，然后依据某些适应性条件测算这些候选解的适应度，并根据适应度保留某些候选解、放弃其他候选解，最后对保留的候选解进行某些操作，生成新的候选解。

实验过程

实验环境

所用机器型号为 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

下载并获得数据集eil101.tsp

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

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

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

预处理工作

此处完成了如下的工作：

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

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()

局部搜索

在这里我设计了四种策略用于产生邻居解：

1. 两点交换 swap_pos：随机找到路径上的两个点并将其交换
2. 区间反转 reverse_pos：随机选中路径上的一个区间将其反转
3. 首尾对调 flip_pos：将路径随机划分成连续的三段，将首尾两段对调
4. 螺旋对调 rotate_pos：将路径随机划分成连续的三段，将前两段对调

原始序列	两点交换	区间反转	首尾对调	螺旋对调
1	1	1	7	3
2	2	2	8	4
3	6	6	9	5
4	4	5	3	6
5	5	4	4	1
6	3	3	5	2
7	7	7	6	7
8	8	8	1	8
9	9	9	2	9

def variation(solution):
    def swap_pos(solution):
        new_s = solution.copy()
        pos = random.sample(range(len(new_s)), 2)
        pos.sort()
        new_s[pos[0]], new_s[pos[1]] = new_s[pos[1]], new_s[pos[0]]
        return new_s

    def reverse_pos(solution):
        new_s = solution.copy()
        pos = random.sample(range(len(new_s)), 2)
        pos.sort()
        new_s1 = new_s[0:pos[0]:1]
        new_s2 = new_s[pos[0]:pos[1]:1]
        new_s3 = new_s[pos[1]:len(new_s):1]
        return new_s1+list(reversed(new_s2))+new_s3

    def flip_pos(solution):
        new_s = solution.copy()
        pos = random.sample(range(len(new_s)), 2)
        pos.sort()
        new_s = new_s[pos[1]:len(new_s):1] + \
            new_s[pos[0]:pos[1]:1]+new_s[0:pos[0]:1]
        return new_s

    def rotate_pos(solution):
        new_s = solution.copy()
        pos = random.sample(range(len(new_s)), 2)
        pos.sort()
        new_s = new_s[pos[0]:pos[1]:1] + \
            new_s[0:pos[0]:1]+new_s[pos[1]:len(new_s):1]
        return new_s
    shuffer = random.choice([swap_pos, reverse_pos, flip_pos, rotate_pos])
    return shuffer(solution)


def local_search(adjmap, num_epoch=100000):
    # 构造一组初始解
    solution = [i for i in range(len(adjmap))]
    length = [cal_length(solution, adjmap)]
    for i in range(1, num_epoch+1):
        new_s = variation(solution)
        new_l = cal_length(new_s, adjmap)
        if new_l < length[-1]:
            solution = new_s
            length.append(new_l)
        else:
            length.append(length[-1])
    return solution, length


tsp(local_search)

模拟退火

在局部搜索的基础上，不改变局部搜索的策略，仅仅添加上温度的控制，以此来让算法在一定的概率下接受差解。

在这里我设定起始温度 t_beg 为 100 度，接受温度为 0.01，退火系数为 0.99，每个温度内循环 100 次。这样一共产生了 $100\times \log_{0.99}\frac{0.01}{100}\approx 91642$ 个新解，接近于局部搜索中设定的十万次，便于对比。

def simulated_annealing(adjmap, t_beg=100, t_end=0.01, q=0.99, l=100):
    # 构造一组初始解
    solution = [i for i in range(len(adjmap))]
    length = [cal_length(solution, adjmap)]
    t = t_beg
    while t > t_end:
        for i in range(l):
            new_s = variation(solution)
            new_l = cal_length(new_s, adjmap)
            pre_l = cal_length(solution, adjmap)
            df = new_l-pre_l
            if df < 0 or numpy.exp(-df / t) > random.random():  # 退火，一定概率下接受差解
                solution = new_s
        length.append(cal_length(solution, adjmap))
        t *= q
    return solution, length


tsp(simulated_annealing)

遗传算法

算法的实现思路如下：

1. 计算开始时，随机初始化一定数目的个体，并计算每个个体的适应度值，产生第一代（初始种群），随后反复执行第二步到第五步；
2. 按照适应度值选择若干（chose_num）个体，作为亲代（适者生存繁衍）；
3. 亲代产生子代，并按一定概率 p_c 进行交叉操作；
4. 产生的子代按一定概率 p_m 变异；
5. 当种群数量超过环境容纳上限 up_bound 时淘汰适应度差（路径长）的个体。

在第四步中我复用了之前局部搜索的变异操作，但是仍然需要实现第三部中「染色体交叉」的操作。

举例：例如，如果双亲的路径编码分别如下：

双亲 1	双亲 2
1	5
2	4
3	6
4	9
5	2
6	1
7	7
8	8
9	3

选中斜体部分序列进行交换，得到原始子代的序列：

原始子代 1	原始子代 2
1	5
2	4
6	3
9	4
2	5
1	6
7	7
8	8
9	3

然后我们需要确定交换序列部分的映射关系：

子代 1	子代 2
6	3
9	4
2	5
1	6

将其应用到未变换部分从而得到合法的子代序列：原始子代 1 的未交叉部分发生了 $1\to 6\to 3,2\to 5,9\to 4$ 的变化；原始子代 2 的未交叉部分发生了 $5\to 2,4\to 9,3\to 6\to 1$ 的变化。注意到当发生映射后仍然和交叉部分冲突时需要多次应用映射规则直至不冲突。最终得到合法的子代个体。

合法子代 1	合法子代 2
3	2
5	9
6	3
9	4
2	5
1	6
7	7
8	8
4	1

我设定的繁衍轮数 num_epoch 为 25000，每次选择的亲代数量 chose_num 为 4，种群上限 up_bound 值为 10，产生染色体交叉的概率 p_c 是 0.9，产生染色体变异的概率 p_m 为 0.4。这样产生的新个体数量为 $25000\times 4=100000$，与之前局部搜索的数量相同，便于对比。

def cross(p1, p2):
    a = p1.copy()
    b = p2.copy()

    pos = random.sample(range(len(p1)), 2)
    pos.sort()

    # print pos[0],end
    # 建立映射关系
    cross_map = {}
    is_exist = False
    # 初步映射
    for i in range(pos[0], pos[1]):
        if a[i] not in cross_map.keys():
            cross_map[a[i]] = []
        if b[i] not in cross_map.keys():
            cross_map[b[i]] = []

        cross_map[a[i]].append(b[i])
        cross_map[b[i]].append(a[i])

    # 处理传递映射 如1:[6],6:[3,1]-->1:[6,3,1],6:[3,1]
    # 计算子串中元素出现的个数，个数为2，则该元素为传递的中间结点，如如1:[6],6:[3,1],‘6’出现的次数为2
    appear_times = {}
    for i in range(pos[0], pos[1]):
        if a[i] not in appear_times.keys():
            appear_times[a[i]] = 0
        if b[i] not in appear_times.keys():
            appear_times[b[i]] = 0

        appear_times[a[i]] += 1
        appear_times[b[i]] += 1

        if a[i] == b[i]:
            appear_times[a[i]] -= 1

    for k, v in appear_times.items():
        if v == 2:
            values = cross_map[k]
            for key in values:
                cross_map[key].extend(values)
                cross_map[key].append(k)
                cross_map[key].remove(key)
                cross_map[key] = list(set(cross_map[key]))

    # 使用映射关系交叉
    # 先映射选中的子串
    temp = a[pos[0]:pos[1]].copy()
    a[pos[0]:pos[1]] = b[pos[0]:pos[1]]
    b[pos[0]:pos[1]] = temp

    # 根据映射规则映射剩下的子串
    seg_a = a[pos[0]:pos[1]]
    seg_b = b[pos[0]:pos[1]]

    remain = list(range(pos[0]))
    remain.extend(range(pos[1], len(a)))

    for i in remain:
        keys = cross_map.keys()
        if a[i] in keys:
            for fi in cross_map[a[i]]:
                if fi not in seg_a:
                    a[i] = fi
                    break
        if b[i] in keys:
            for fi in cross_map[b[i]]:
                if fi not in seg_b:
                    b[i] = fi
                    break
    return [a, b]


def genetic_algorithm(adjmap, num_epoch=25000, chose_num=4, up_bound=10, p_c=0.9, p_m=0.4):
    # 构造初始种群
    w = [i for i in range(len(adjmap))]
    population = []
    for i in range(up_bound):
        population.append(w.copy())
        w = variation(w)
    length = []
    for i in range(num_epoch):
        # 适应者繁殖
        w = [cal_length(p, adjmap) for p in population]
        sum = 0
        for ww in w:
            sum += ww
        chose = random.choices(population, weights=[
                               numpy.exp(-ww/sum) for ww in w], k=chose_num)
        for j in range(1, chose_num, 2):
            parent = [chose[j-1], chose[j]]
            if random.random() < p_c:
                parent = cross(parent[0], parent[1])
            for ch in parent:
                if random.random() < p_m:
                    ch = variation(ch)
                population.append(ch)
        # 不适者淘汰
        population.sort(key=lambda sol: cal_length(sol, adjmap))
        while len(population) > up_bound:
            population.pop()
        length.append(cal_length(population[0], adjmap))
    return population[0], length


tsp(genetic_algorithm)

结果分析

局部搜索

以下是各次运行时间（s）及求解时间的均值、方差、标准差：
[2.6684435999995912, 2.604407300001185, 2.5684375000018917, 2.378306499998871, 2.226840600000287, 2.380468999999721, 2.2538437000002887, 2.300632200000109, 2.3810795999997936, 2.278879900000902]
2.404133990000264
0.021914453292613316
0.14803531096536837
以下是各次运行求得解的大小及其均值、方差、标准差：
[696.1930716255447, 680.5388755363587, 691.153152755653, 680.3436846688317, 710.2411460614444, 695.9271346529877, 709.0199083731379, 696.2208983797532, 681.3700501082578, 700.982603576283]
694.199052573825
108.65058008479028
10.423558897266819
[37, 13, 43, 90, 99, 36, 97, 91, 96, 86, 41, 42, 14, 56, 1, 12, 94, 93, 5, 100, 26, 27, 52, 57, 39, 25, 20, 72, 71, 3, 55, 74, 73, 21, 40, 22, 66, 38, 24, 54, 53, 23, 28, 77, 33, 34, 70, 64, 65, 19, 50, 8, 80, 32, 78, 2, 76, 67, 79, 11, 75, 49, 0, 68, 51, 88, 95, 98, 58, 92, 84, 4, 83, 82, 59, 17, 81, 6, 87, 30, 69, 29, 31, 89, 62, 9, 61, 10, 63, 48, 35, 46, 18, 47, 45, 7, 44, 16, 60, 15, 85]
以下是所求最优解及其可视化：
680.3436846688317

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

可以看出，在最初的 20000 次迭代中可以快速收敛，收敛的速度不断变慢；40000 次之后已经陷入一个局部最优解所在的区域，并且逐渐趋于稳定。

从最优解的可视化上可以很直观看到，得到的解没有出现交叉的情况，属于比较稳定的结构。

模拟退火

以下是各次运行时间（s）及求解时间的均值、方差、标准差：
[4.028491099998064, 4.096584099999745, 4.395690099998319, 4.359389299999748, 4.236625900000945, 4.337160200000653, 4.295187499999884, 4.29896159999771, 4.315229199997702, 4.238674100000935]
4.26019930999937
0.012135203608892125
0.1101599001855581
以下是各次运行求得解的大小及其均值、方差、标准差：
[689.7491434950963, 677.6316535178441, 710.2786571395188, 695.0104107819034, 698.7898420494348, 682.7677502467288, 682.8278839175407, 687.828082173087, 692.5436782781678, 684.0435322027494]
690.1470633802071
81.68238425936455
9.037830727523312
[16, 85, 37, 13, 43, 99, 97, 36, 91, 96, 86, 1, 56, 41, 42, 14, 40, 21, 73, 72, 20, 71, 74, 55, 38, 66, 22, 3, 24, 54, 53, 23, 28, 77, 33, 8, 34, 70, 64, 65, 19, 29, 69, 50, 80, 32, 78, 2, 76, 67, 79, 11, 75, 49, 0, 68, 30, 87, 6, 81, 47, 18, 10, 61, 9, 31, 89, 62, 63, 48, 35, 46, 45, 7, 44, 82, 59, 17, 88, 51, 26, 27, 100, 52, 25, 39, 57, 12, 5, 93, 94, 58, 95, 98, 92, 84, 90, 15, 60, 4, 83]
以下是所求最优解及其可视化：
677.6316535178441

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

从最优解的收敛过程中，我们可以很显著的看出模拟退火算法的一些特征：一开始系统处于「高能态」的时候，产生解的解非常的不稳定，甚至在最开始的时候有一个非常陡峭的增幅；在震荡了一段时间之后，系统开始趋于平和，并在降温 600 次之后稳定下来。

从结果的可视化上来看，路径上出现了一个交叉，不过其结果仍然比局部裸搜索的结果要好一些。这说明模拟退火算法具有更佳的鲁棒性，不依赖于初始解，并且可以通过调整模拟退火的环境变量来探索更大的解空间；而局部搜索则不然，如果初始解选择的不好或者数据比较特殊的话就非常尴尬了。

此外，由于计算过程和逻辑比局部搜索更加复杂，模拟退火的用时也多于局部搜索，是后者的两倍。

遗传算法

以下是各次运行时间（s）及求解时间的均值、方差、标准差：
[14.54444439999861, 14.35590890000094, 15.490778699997463, 13.975171499998396, 13.794036399998731, 14.393401699999231, 14.615606200000911, 14.529433499999868, 14.908401900000172, 14.58785749999879]
14.51950406999931
0.1974510827827878
0.4443546812882563
以下是各次运行求得解的大小及其均值、方差、标准差：
[716.6000944269756, 689.9624497677642, 751.6800288206394, 722.4561170763457, 694.3045029478357, 703.2937838680601, 728.1576655214293, 717.2876165930192, 716.2627410735281, 710.8916651889701]
715.0896665284566
279.2546395035423
16.710913784217258
[26, 68, 27, 25, 39, 20, 72, 71, 73, 21, 40, 74, 55, 22, 66, 38, 3, 24, 54, 53, 11, 79, 67, 23, 28, 33, 77, 78, 2, 76, 75, 49, 0, 50, 32, 80, 8, 34, 70, 64, 65, 19, 29, 69, 31, 89, 62, 63, 10, 18, 48, 35, 45, 46, 47, 81, 6, 87, 61, 9, 30, 51, 17, 59, 82, 7, 44, 16, 83, 4, 60, 15, 85, 37, 42, 14, 56, 86, 41, 13, 43, 90, 84, 99, 36, 97, 92, 98, 95, 58, 91, 96, 94, 93, 5, 88, 12, 1, 57, 52, 100]
以下是所求最优解及其可视化：
689.9624497677642

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

在繁衍达到 15000 代的时候已经找到了比较稳定的优解。最终得到的结果如图所示，没有交叉部分，比之前模拟退火的结果略差一些但相差极少可以接受。

可以看到，和局部搜索或是模拟退火相比，遗传算法的用时多了很多，新增的交叉操作要取计算映射关系，是比较耗时的。

遗传算法的模型在单次迭代中繁衍的过程没有循环依赖，天然的适合用并行的算法去优化，这既是优势也是劣势：我这里选择的问题规模不够大，并且我自己的笔记本内存仅有 8GB，很难同时计算并存储大规模的种群，限制了遗传算法的用武之地。可以预见，当问题规模进一步扩大，并且有充足计算资源的时候，种群算法一定能取得更优的表现。

结论

本次实验虽然不算很难，但是还是遇到了一些坑点和感悟。这里我整理一下：

1. 变异操作。最开始我还设置了「轮换」的操作：将路径上所有点顺时针旋转若干单位。但是这么做反而让收敛的速度降低了。后来，我突然想到，在旅行商问题中，起点是可以任意选择的，这就是说轮换操作在本问题中是画蛇添足了。
2. 种群「劣汰」。一开始我按照这篇博客里的思路实现种群算法，然而实现代码运行极慢。经过思考，我发现他并没有实现遗传算法中的「劣汰」过程，导致种群无节制的扩张，这当然是不符合实际情况的。这不由得让我由衷佩服起造物主的精巧规则，缺少哪一部分都不足以形成合理的结构。
3. 「基因库的多样性」。最开始我在遗传算法初始化种群的时候，使用了若干和局部搜索中相同的初始解作为起始种群。然而这样做的话求解的收敛速度非常慢，几万轮子代仍然和亲代区别不大。因此我在生成起始种群时强制增加了多次变异操作，从而增加了这个种群基因库的多样性。
4. 尽量多使用库提供的接口，可以极大简化代码实现，从而把注意力放到算法本身上面去。比如按照比重选取亲代（适者生存繁衍）的过程，网上的很多代码都是手动实现，而我使用了 random.choice 函数就可以很方便的解决问题。最终我的代码行数不到参考资料中的三分之一，可以说是非常让人满意了。

主要参考文献

• TSPLIB
• 现代优化算法三部曲 | 模拟退火算法 - nitric acid 的文章 - 知乎
• 遗传算法解决 TSP 问题(python)_NLPTutorial-CSDN 博客_遗传算法 tsppython