TSP、MTSP问题遗传算法详细解读及python实现

写在前面

遗传算法是一种求解NPC问题的启发式算法，属于仿生进化算法族的一员。仿生进化算法是受生物行为启发而发明的智能优化算法，往往是人们发现某种生物的个体虽然行为较为简单，但生物集群通过某种原理却能表现出智能行为。于是不同的人研究不同的生物行为原理，受到启发而发明出新的仿生进化算法。比如免疫优化算法，蚁群算法，模拟退火算法等，这些算法以后也会简单介绍。
本文的主题是遗传算法，该算法也是受到生物行为启发。物竞天择，适者生存，优胜劣汰，是该优化算法的核心思想。笔者在业务中需要用到遗传算法求解TSP问题，但是网上能查找到的资料对遗传算法的讲解不够通俗易懂，往往上来就是遗传变异交叉，对于我这样的初学者来说有点不知所云，于是不得不直接看源码，一行一行地理解代码的意思，才弄懂了原理。这种方法对于初学者和编程基础薄弱者颇为困难，而且费时费力，苦不堪言。同时，由于读者可能熟练掌握的是不同的语言，因此若代码是某一种语言编写的，那么掌握其他语言的读者很可能难以吸收，浪费了资源。此外，网上关于TSP问题的资料很多，但是关于MTSP问题的资料却凤毛麟角。因此有了创作本文的意图，旨在用最通俗详尽的语言深入浅出地解释遗传算法解TSP、MTSP问题的原理及应用

遗传算法解TSP问题原理

一、TSP问题

旅行商问题，即TSP问题（Traveling Salesman Problem）又译为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。
想要求解出TSP问题的最优解，目前唯一的方法是穷举出所有的路径。然而，路径的数量级是n!，也就是目标点数量的阶乘。当n为14时，n!已经大于800亿。当n更大，为30，40 时，更是天文数字，即使计算机一秒钟计算一亿次，其求解时间也远大于我们的寿命。因此，穷举法解TSP问题根本不可行，需要别的启发式算法来代替。

二、遗传算法解TSP原理

2.1

笔者认为，大多数的遗传算法讲解让人难以理解的原因在于其术语名称的意义不明确，往往与其实际意义相去甚远。譬如，什么是适应度函数，什么是遗传，变异，交叉。乍一听总有不知所云的感觉，这跟TSP问题求最短路径有什么关系。因此本文将从TSP问题出发，讲解遗传算法的思想和原理。

2.2

对于TSP问题，我们拿到的是n个目标点的坐标，作为例子，我们视为面对的是n=10个城市。有了坐标，我们便可以计算出n个城市之间的距离，得到一个n*n的矩阵，这个矩阵表示的是城市之间两两连接形成的无向图，边的权重是城市之间的距离，而TSP问题则是从图中找出一个包含所有城市的最小的环。

2.3

问题明确了，接下来就是遗传算法登场了。

2.3.1 种群初始化

首先是初始化种群，TSP里的初始化种群其实就是一个长度为n=10的包含1~n(这个例子里为10)且不含重复元素的序列，其意义就是一个人从某个点出发，随机访问下一个未访问过的城市，直到所有的城市都访问完毕，他再从最后一个城市返回出发城市。他的轨迹就是一个包含了所有城市且没有重复访问的环。种群数量设为M，那么该种群初始化的意义就是M个人独立、随机、不重复地访问一遍各个城市，单个个体轨迹构成一个环。

2.3.2 适应度计算

M个个体的轨迹得到了，这些轨迹是随机游走得到的，当然很可能远远不包含属于最优解，而是比最优解坏得多。但是，随机游走的路径也有大有小，有好有坏。因此需要有一个函数衡量个体的好坏，也就是环路径的长短。
那么这个函数就被称为适应度函数，其功能是衡量个体的好坏。对于TSP问题，其适应度函数当然是距离相关的了。个体的好坏，衡量标准就是该个体序列的路径长度。路径越长，个体越“坏”，路径越短，个体越“好”。因此，设序列为x,那么该序列的路径长度便是d(x)，而适应度函数则应该取为1/d(x)，适应度越大，个体越优。

2.3.3 选择

有了每个个体的适应度，就能评价每个个体的好坏。物竞天择，优胜劣汰，将优秀的个体选择出来进行交配，以期得到更好的个体，并由此不断进化，一代代传承，后代不断比前一代变得更好，最终收敛，种群中的某个个体达到了优秀的极限，便是最优解。
对于TSP问题，选择的具体操作是，计算出所有个体的适应度，也就是路径距离的倒数。然后将所有个体的适应度归一化，得到概率。然后从数量为n的个体中以轮盘赌的形式选择出若干个个体，视为优秀个体。
举个例子，为了简单起见，设种群大小m=4。并计算出了四个个体的适应度分别为1，2，3，4。对适应度进行归一化得到概率：0.1，0.2，0.3，0.4。那么这四个概率就可以构成一个轮盘，每个个体对应的被选择的概率分别为0.1，0.2，0.3，0.4。随后，在这个转盘上转动m次，此处为4。将选中的个体放入一个集合，未选中的个体抛弃。
注意，这是一种重复抽样的选择方式，并且重复抽到的个体不去重。每次抽样，轮盘都是相同的轮盘，并没有改变。比如，在这个例子中，种群数量为4，那么就要抽取四次，那么0.4对应的个体很可能被多次抽样，而0.1对应的个体很可能未被抽到，直接淘汰。那么被选择的个体集合就很可能含有多个适应度为4的个体，并且被选择的集合的种群数量仍旧为m，只不过抛弃了较坏的个体，而可能保留了多个优秀的个体。个体越优秀，被重复选择的概率就越大。

2.3.4 交叉（遗传）

遗传算法中的交叉属于遗传算法的核心以及关键步骤，其思想受启发于生物遗传中的染色体交叉。
对于TSP问题，有如下的交叉方式进行交叉。

1. Partial-Mapped Crossover（部分映射交叉）
   

2. Order Crossover（顺序交叉）

3. Position-based Crossover（基于位置的交叉）

那么，在这个被选择的大小为种群数量m的新的种群中，如何通过上述方式交叉呢？这里有一个超参数pc(probability of crossover)。将种群打乱之后，随机两两组合。每对组合以pc的概率采用上述方式之一进行交叉，以（1-pc）的概率不进行交叉，而是二者直接进入下一代。pc的值一般在0.5到0.9之间。
通过上述方式，便能产生一个同样大小为m的新的种群。此刻，选择、交叉的过程已经完成，还剩最后一步操作——变异。

2.3.5 变异

对于上一步交叉得到的大小为m的种群，以pm的概率选择出部分个体进行变异操作。选择出来的个体序列随机选择两个位置上的数进行交换。如，其中一个被选择的个体的序列为196278354，那么在这个序列中随机选择两个位置上的数，比如第二个数9和第四个数2，进行交换。那么得到的变异后的序列就是：126978354，变异操作完成。

2.3.6 传承

至此，种群的迭代更新方式基本阐述完毕，还剩下最后一步操作，传承。这是说，将上一轮中最优的个体（对于TSP问题，也就是路径距离最短的序列），保留下来，并用他替换掉新产生的种群中最差的个体。这一步的意义在于，上一轮中最优的个体被选择的几率也最大，但是它与其他个体交叉之后得到的新个体不一定优于它自己。如果这样，那么这个最优的个体便被覆盖、迭代掉了，这样很可能造成下一代中的所有个体都比上一代的最优个体差，优秀的基因被淘汰。为了避免这样的情况，需要这一步传承，保证每次迭代之后的最优个体不坏于之前出现的所有个体。

三、遗传算法解TSP问题python实现

对于单纯的TSP问题，用python求解非常简便，原因在于python已有十分强大的遗传算法工具包，用户只需要将需要求解的目标点位置转化成n*n的距离矩阵distance_matrix，然后直接调用python里的库sko.GA即可，下面是示例代码。读者不用担心直接调包无法了解算法实现的细节，在下一节的MTSP问题本文将基于纯python手写实现遗传算法。

// An highlighted block
def cal_total_distance(routine):
    '''The objective function. input routine, return total distance.
    cal_total_distance(np.arange(num_points))
    '''
    num_points, = routine.shape
    return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
from sko.GA import GA_TSP
#num_points为城市数量，size_pop为种群数量，max_iter为迭代次数
#prob_mut为变异率
ga_tsp = GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=2000, prob_mut=1)
best_points, best_distance = ga_tsp.run()

其中best_points返回的是最优路径的顺序，best_distance返回的是最优路径的距离。

用python工具包实现TSP问题就是这样简洁容易。在下一节，要讨论的是更为复杂的MTSP问题。

四、遗传算法解MTSP问题

4.1问题描述

问题描述：m个旅行商去旅游 n个城市，一种规定是都必须从同一个出发点出发，而且返回原出发点，需要将所有的城市遍历完毕，每个城市只能游历一次，但是为了路径最短可以路过这个城市多次。这个就是多旅行商问题。是在TSP问题的基础上进行了扩展。

4.2 问题分析

关于MTSP问题，可以找到的资料比TSP要少很多。对于该问题的求解，可以基于TSP问题的思路进行修改。首先需要明确具体的业务需求，是多个旅行商从不同的城市出发，遍历所有的目标点并回到自己的原点，还是都从同一个点出发回到所有起点，还是回到同一个点但该点不是起点。另外，每个旅行商访问城市的数量是任意的，还是两个旅行商需要访问的城市数目是大致相等的，还是规定了每个旅行商访问的数量具体是多少？
不论问题是哪种，其思想是差不多的。本文以多个旅行商从同一点出发，回到同一个起点，且每人访问的城市数量相同的情形。首先，跟TSP问题一样，初始化大小为M的种群，每个个体是一个不含重复点的长度为n的序列。然后重点是设置断点，将该序列拆分成m段，也就是旅行商的数量。此后，将m段序列的首尾都添上起点，形成m个旅行商各自的环路径，然后计算各个环的距离的总和的倒数，作为适应度函数的值，后续操作按照TSP问题的遗传算法求解即可。在本例中，考虑的m取2，两位旅行商。断点取序列的中间，也就是序列平均分为两半。若读者想不限制各位旅行商的访问城市的数目相等，只需要设置断点在不同的位置即可。

4.2 MTSP问题python实现

对于该MTSP问题，笔者没有找到相关的python实现包，因此只能手撸代码如下：

import numpy as np
from matplotlib import pyplot as plt
import math
import copy
import random
def init_population(length, num):
    li = list(range(length))
    print(li)
    population = []
    for i in range(num):
        random.shuffle(li)
        population.append(copy.deepcopy(li))
    return population
def aimFunction(entity, DMAT, break_points):
    """
    目标函数
    :param entity: 个体
    :param DMAT: 距离矩阵
    :param break_points: 切断点
    :return:适应度函数值
    """
    distance = 0
    break_points.insert(0, 0)
    break_points.append(len(entity)) 
    routes = []
    for i in range(len(break_points) - 1):
        routes.append(entity[break_points[i]:break_points[i + 1]])
    # print(routes)
    #上面代码的作用是将路径拆成m段。break_points的长度设置为k时，则将路径拆分成了k+1段，此处的k 取的1。
    for route in routes:
        if 0 in route:
            route.remove(0)
        route.insert(0,0)
        route.append(0)#此处给每段路径添上首尾点0，因为本例所有旅行商都从0这个城市出发，具体从哪里出发根据实际问题修改该设定。
        for i in range(len(route)-1):
            distance += DMAT[route[i],route[i+1]]

    return 1.0/distance
#返回种群所有个体的适应度列表
def fitness(population, DMAT, break_points, aimFunction):
    """
    适应度
    :param population: 种群
    :param DMAT: 距离矩阵
    :param break_points:切断点
    :param aimFunction: 目标函数
    :return:
    """

    value = []
    for i in range(len(population)):
        value.append(aimFunction(population[i], DMAT, copy.deepcopy(break_points)))
        # weed out negative value
        if value[i] < 0:
            value[i] = 0
    return value
#这里传入的是种群和每个个体的适应度
#这里的轮盘赌与蚁群算法的有一定区别。这里对适应度归一化得到概率之后，每个个体被选中的概率就是这个概率
#对每次被选中的个体的数目没有限制，完全随机，限制的是选择次数n与种群数目相同，使得新的种群数量与旧的种群一致
def selection(population, value):
    # 轮盘赌选择
    fitness_sum = []
    for i in range(len(value)):
        if i == 0:
            fitness_sum.append(value[i])
        else:
            fitness_sum.append(fitness_sum[i - 1] + value[i])

    for i in range(len(fitness_sum)):
        fitness_sum[i] /= sum(value)

    # select new population
    population_new = []
    for i in range(len(value)):
        rand = np.random.uniform(0, 1)
        for j in range(len(value)):
            if j == 0:
                if 0 < rand and rand <= fitness_sum[j]:
                    population_new.append(population[j])

            else:
                if fitness_sum[j - 1] < rand and rand <= fitness_sum[j]:
                    population_new.append(population[j])
    return population_new
#对于被选中的双亲，随机两两组合。并以pc的概率交配
#若没有被选择交配，则直接双亲进入下一代。如果被选择，则交换中间片段。
def amend(entity, low, high):
    """
    修正个体
    :param entity: 个体
    :param low: 交叉点最低处
    :param high: 交叉点最高处
    :return:
    """
    length = len(entity)
    cross_gene = entity[low:high]  # 交叉基因
    not_in_cross = []  # 应交叉基因
    raw = entity[0:low] + entity[high:]  # 非交叉基因


    for i in range(length):
        if not i in cross_gene:
            not_in_cross.append(i)

    error_index = []
    for i in range(len(raw)):
        if raw[i] in not_in_cross:
            not_in_cross.remove(raw[i])
        else:
            error_index.append(i)
    for i in range(len(error_index)):
        raw[error_index[i]] = not_in_cross[i]

    entity = raw[0:low] + cross_gene + raw[low:]

    return entity


def crossover(population_new, pc):
    """
    交叉
    :param population_new: 种群
    :param pc: 交叉概率
    :return:
    """
    half = int(len(population_new) / 2)
    father = population_new[:half]
    mother = population_new[half:]
    np.random.shuffle(father)
    np.random.shuffle(mother)
    offspring = []
    for i in range(half):
        if np.random.uniform(0, 1) <= pc:
            # cut1 = np.random.randint(0, len(population_new[0]))
            # if cut1 >len(father[i]) -5:
            #     cut2 = cut1-5
            # else:
            #     cut2 = cut1+5
            cut1 = 0
            cut2 = np.random.randint(0, len(population_new[0]))
            if cut1 > cut2:
                cut1, cut2 = cut2, cut1
            if cut1 == cut2:
                son = father[i]
                daughter = mother[i]
            else:
                son = father[i][0:cut1] + mother[i][cut1:cut2] + father[i][cut2:]
                son = amend(son, cut1, cut2)
                daughter = mother[i][0:cut1] + father[i][cut1:cut2] + mother[i][cut2:]
                daughter = amend(daughter, cut1, cut2)

        else:
            son = father[i]
            daughter = mother[i]
        offspring.append(son)
        offspring.append(daughter)

    return offspring
#这里的变异是最简单的变异法则，直接随机选取两个位置上的数进行交换
def mutation(offspring, pm):
    for i in range(len(offspring)):
        if np.random.uniform(0, 1) <= pm:
            position1 = np.random.randint(0, len(offspring[i]))
            position2 = np.random.randint(0, len(offspring[i]))
            # print(offspring[i])
            offspring[i][position1],offspring[i][position2] = offspring[i][position2],offspring[i][position1]
            # print(offspring[i])
    return offspring
#主函数
if __name__ == '__main__':

#这里的graph为距离矩阵，需要自己定义后传入

    DMAT=graph
    break_points = [len(graph)//2]#这里是将所有城市从中间断开成两段。读者可以根据需求切断成任意段（每段代表一个旅行商），并且可以在任意位置切断。
    population = init_population(len(graph), 90)

    t = []
    dic_result={}
    for i in range(20000):
        # selection
        value = fitness(population, DMAT, break_points, aimFunction)
        population_new = selection(population, value)
        # crossover
        offspring = crossover(population_new, 0.65)
        # mutation
        population = mutation(offspring, 0.02)
        # if i % 1 == 0:
        #     show_population(population)
        result = []
        for j in range(len(population)):
            result.append(1.0 / aimFunction(population[j], DMAT, copy.deepcopy(break_points)))

        t.append(min(result))
        min_entity = population[result.index(min(result))]
        routes = []
        break_points_plt = copy.deepcopy(break_points)
        break_points_plt.insert(0, 0)
        break_points_plt.append(len(min_entity))
        for i in range(len(break_points_plt) - 1):
            routes.append(min_entity[break_points_plt[i]:break_points_plt[i + 1]])
        for route in routes:
            if 0 in route:
                route.remove(0)
            route.insert(0,0)
            route.append(0)
        dic_result[min(result)]=routes
        if i % 400 == 0:
            min_entity = population[result.index(min(result))]
            routes = []
            break_points_plt = copy.deepcopy(break_points)
            break_points_plt.insert(0, 0)
            break_points_plt.append(len(min_entity))
            for i in range(len(break_points_plt) - 1):
                routes.append(min_entity[break_points_plt[i]:break_points_plt[i + 1]])
            for route in routes:
                if 0 in route:
                    route.remove(0)
                route.insert(0,0)
                route.append(0)
            for route in routes:
                print(route)
    print(min(t))#每次迭代的最优路径
    print('最短路径距离之和为',min(t),'米')
    plt.plot(t)#画出每次迭代的最优个体对应的环路程
  

这里我传入的是一个22*22的距离矩阵，读者在使用时需要自定义矩阵和断点。该矩阵第一行元素如下：
[0.00000e+00, 1.00000e+00, 6.83910e+04, 6.81450e+04, 4.29630e+04,
4.13760e+04, 7.52490e+04, 7.09670e+04, 7.82330e+04, 7.34120e+04,
4.08460e+04, 4.09140e+04, 8.49820e+04, 8.49350e+04, 8.18330e+04,
5.04900e+04, 4.71150e+04, 3.36550e+04, 3.35980e+04, 3.16970e+04,
6.96790e+04, 7.21160e+04, 6.93410e+04, 5.62990e+04, 5.16090e+04]
断点设置在中间，也就是22整除2，将22个点均分成两份，每份11个点，然后两位旅行商分别从第一个点出发，巡游各自需要访问的点，再回到起点。具体迭代过程在上文中已经阐述。最终得到的结果plot如下：

五、总结

遗传算法是一种仿生进化启发式算法，该算法不保证搜索到最优解，但是借鉴生物进化中优胜劣汰，物竞天择的思想，随机初始化个体之后使用适应度函数对个体的优劣进行评测，然后选择较为优秀的个体进行交配，较劣的个体淘汰。人们期望优秀的个体之间交配能产生更优的个体，如此代代传承进化下去，得到越来越优秀的个体。当然，优秀的个体之间交叉之后，并不保证其后代一定优于双亲，因此遗传算法也可能收敛到局部最小值或者近优值。当然，这也是所有仿生进化算法的特性。目前还没有能完全解决NPC问题的有效算法。