粒子群算法求解旅行商问题

客官°小女子只卖身不卖艺 2022-11-15 06:25 260阅读 0赞

## 算法原理 ##

旅行商问题是一个经典的NP问题，假设有N个城市，需要确定一个访问顺序，使得每个城市都访问一面，最后回到起点城市，且保证行走的总距离最短。  
       假设随机生成10个城市坐标，城市之间的距离用欧式距离表示。

在以往例子中，粒子xi的值是具有意义的值，而这里需要优化序列的顺序，因此粒子xi如何表示序列以及序列的改进是粒子群算法的一个关键所在。同时也应该认识到，虽然智能算法的原理比较容易理解，但在不同应用里，巧妙的构造解xi的表示形式，以及解xi的进化方式，才是关键。

在求解TSP问题时，可以用一个序列s表示城市的访问顺序，因此xi可以表示成如下形式：  
                            xi=(1,2,3,4,5)  
       式中表示有5个城市的TSP问题，其访问顺序为1-3-2-5-4-1.因此全部城市的所有可能序列构成了问题的搜索空间，例子位置的更新意味着粒子xi从所有城市的一种序列si变成另一种序列sj。  
       那么序列是如何变化得，这里引入交换子 和交换序列的概念。  
      交换子定义为s=Swap x(i,j)表示交换子s作用在序列x上使位置i的元素和位置j的元素相互交换，如x=(1,3,2,5,4)，对其添加一个交换操作s=Swap x(1,3)，其新的序列为x=(2,3,1,5,4)，在粒子群算法中，交换子可以看作是粒子的速度，能改变x的位置。  
      交换序列 定义为一组有前后顺序的交换子的集合，即ss=\[Swap1,Swap2…\]，表示对序列x进行多次交换子交换操作，如x=(1,3,2,5,4)，交换序列ss=\[(3,2),(1,5)\]，当序列进行交换序列后，新的序列为x=(4,2,3,5,1)  
      解释了交换子和交换序列后，重新定义粒子的速度和更新公式：  
![20210331142440391.png_pic_center][]  
      式中Pbesti-xi表示有一个交换序列ss，它使xi经过ss变换得到Pbest，同理Gbesti-xi也是一个交换序列ss，符号⨁表示两个交换序列合并一个交换序列，如\[(3,2),(1,5)\]⨁\[(4,2)\]=\[(3,2),(1,5),(4,2)\]；符号+表示对一个序列按照交换序列ss执行交换操作。因此位置更新公式表示一个粒子经过交换序列Vi(t+1) 作用得到新的序列。r1,r2表示0-1之间的随机数，表示一定概率进行交换子操作。

粒子群算法原理路见往期：[粒子群算法求解无约束优化问题 源码实现][Link 1]

## 算例 ##

假设随机生成10个城市坐标

import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    import seaborn as sns
    
    sns.set_style("whitegrid")
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
    mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
    
    # 固定随机数种子
    np.random.seed(1234)
    
    # 一些参数
    city_num = 10  # 城市的数量
    iter_num = 1000  # 算法最大迭代次数
    
    # 随机生成city_num个城市的坐标，注意是不放回抽样
    X = np.random.choice(list(range(1, 100)), size=city_num, replace=False)
    Y = np.random.choice(list(range(1, 100)), size=city_num, replace=False)
    
    plt.scatter(X, Y, color='r')
    plt.title('城市坐标图')
    plt.show()

![20210331134548488.png_pic_center][]

步骤1：初始化参数

import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    import seaborn as sns
    
    sns.set_style("whitegrid")
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
    mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
    
    # 固定随机数种子
    np.random.seed(1234)
    
    # 一些参数
    city_num = 10  # 城市的数量
    size = 50  # 种群大小，即粒子的个数
    r1 = 0.7  # pbest-xi 的保留概率
    r2 = 0.8  # gbest-xi 的保留概率
    iter_num = 1000  # 算法最大迭代次数
    fitneess_value_list = []  # 每一步迭代的最优解
    
    # 随机生成city_num个城市的坐标，注意是不放回抽样
    X = np.random.choice(list(range(1, 100)), size=city_num, replace=False)
    Y = np.random.choice(list(range(1, 100)), size=city_num, replace=False)

步骤2：距离函数 和适应度函数

# 计算城市之间的距离
    def calculate_distance(X, Y):
        """ 计算城市两辆之间的欧式距离，结果用numpy矩阵存储 :param X: 城市的X坐标，np.array数组 :param Y: 城市的Y坐标，np.array数组 """
        distance_matrix = np.zeros((city_num, city_num))
        for i in range(city_num):
            for j in range(city_num):
                if i == j:
                    continue
                dis = np.sqrt((X[i] - X[j]) ** 2 + (Y[i] - Y[j]) ** 2)  # 欧式距离计算
                distance_matrix[i][j] = dis
        return distance_matrix
    
    def fitness_func(distance_matrix, xi):
        """ 适应度函数，计算目标函数值. :param distance: 城市的距离矩阵 :param xi: PSO的一个解 :return: 目标函数值，即总距离 """
        total_distance = 0
        for i in range(1, city_num):
            start = xi[i - 1]
            end = xi[i]
            total_distance += distance_matrix[start][end]
        total_distance += distance_matrix[end][xi[0]]  # 从最后一个城市回到出发城市
        return total_distance

步骤3：绘制最优解的图像函数

def plot_tsp(gbest):
        """绘制最优解的图形"""
        plt.scatter(X, Y, color='r')
        for i in range(1, city_num):
            start_x, start_y = X[gbest[i - 1]], Y[gbest[i - 1]]
            end_x, end_y = X[gbest[i]], Y[gbest[i]]
            plt.plot([start_x, end_x], [start_y, end_y], color='b', alpha=0.8)
        start_x, start_y = X[gbest[0]], Y[gbest[0]]
        plt.plot([start_x, end_x], [start_y, end_y], color='b', alpha=0.8)

步骤4：交换序列函数 和交换操作

> import numpy as np  
> xi=np.array(\[1,2,3,4,5,6,7,8,9,10\])  
> xbest=np.array(\[2,1,4,3,6,5,7,8,9,10\])  
> velocity\_ss = \[\]  
> for i in range(len(xi)):  
>       if xi\[i\] != xbest\[i\]:  
>       j = np.where(xi == xbest\[i\])\[0\]\[0\]  
>        print(j)  
>       xi\[i\], xi\[j\] = xi\[j\], xi\[i\] \#表示去掉重复计算，不同处只计算一处  
> \# j=1,3,5

#交换序列函数
    def get_ss(xbest, xi, r):#r=r1或r2
        """ 计算交换序列，即 x2 经过交换序列ss得到x1，对应PSO速度更新公式的： r1(pbest-xi) 和 r2(gbest-xi) :param xbest: pbest 或者 gbest :param xi: 例子当前解 :return: xbest#最优解 """
        velocity_ss = []
        for i in range(len(xi)):
            if xi[i] != xbest[i]:#如果当前索引粒子元素不等于最优粒子元素
                j = np.where(xi == xbest[i])[0][0]#找到当前粒子和最优粒子不同的元素所在索引
                so = (i, j, r)  # 得到交换子,表示以r的概率对i,j进行操作
                velocity_ss.append(so)
                xi[i], xi[j] = xi[j], xi[i]  # 执行交换操作
        return velocity_ss#返回交换子序列
     #交换操作，返回序列xi解
    def do_ss(xi, ss):
        """ 执行交换操作 :param xi: :param ss: 由交换子组成的交换序列 :return: """
        for i, j, r in ss:
            rand = np.random.random()
            if rand <= r:
                xi[i], xi[j] = xi[j], xi[i]
        return xi

步骤5：PSO  
      第一步：计算城市之间的距离矩阵

# 计算城市之间的距离矩阵
    distance_matrix = calculate_distance(X, Y)

第二步：初始化种群的各个粒子的位置，作为个体的历史最优解pbest。用一个 50\*10 的矩阵表示种群，每行表示一个粒子, 每一行是0-9不重复随机数，表示城市访问的顺序

XX = np.zeros((size, city_num), dtype=np.int)
    for i in range(size):
        XX[i] = np.random.choice(list(range(city_num)), size=city_num, replace=False)

第三步：计算每个粒子对应适应度

pbest = XX#初始化粒子最优解
    pbest_fitness = np.zeros((size, 1))
    for i in range(size):
        pbest_fitness[i] = fitness_func(distance_matrix, xi=XX[i])

第四步：计算全局适应度和对应的解gbest

gbest = XX[pbest_fitness.argmin()]
    gbest_fitness = pbest_fitness.min()

第五步：迭代寻优

# 记录算法迭代效果
    fitneess_value_list.append(gbest_fitness)
    
    # 下面开始迭代
    for i in range(iter_num):#遍历迭代次数
        for j in range(size):  # 对每个粒子迭代
            pbesti = pbest[j].copy()  # 此处要用copy，不然出现浅拷贝问题
            xi = XX[j].copy()
            # 计算交换序列，即 v = r1(pbest-xi) + r2(gbest-xi)
            ss1 = get_ss(pbesti, xi, r1)
    
            ss2 = get_ss(gbest, xi, r2)
            ss = ss1 + ss2
            # 执行交换操作，即 x = x + v
            xi = do_ss(xi, ss)
            # 判断是否更优
            fitness_new = fitness_func(distance_matrix, xi)
            fitness_old = pbest_fitness[j]
            if fitness_new < fitness_old:
                pbest_fitness[j] = fitness_new
                pbest[j] = xi
        # 判断全局是否更优# 记录算法迭代效果
        # fitneess_value_list.append(gbest_fitness)
        gbest_fitness_new = pbest_fitness.min()
        gbest_new = pbest[pbest_fitness.argmin()]
        if gbest_fitness_new < gbest_fitness:
            gbest_fitness = gbest_fitness_new
            gbest = gbest_new
        # 加入到列表
        fitneess_value_list.append(gbest_fitness)

第六步：绘制结果

# 打印迭代的结果
    print('迭代最优结果是：', gbest_fitness)
    print('迭代最优变量是：', gbest)
    # 迭代最优结果是： 230.344
    # 迭代最优变量是： [5 8 2 3 6 1 7 0 4 9]
    
    # 绘制TSP路径图
    plot_tsp(gbest)
    plt.title('TSP路径规划结果')
    plt.show()

![20210331154637247.png_pic_center][] ![20210331154823160.png_pic_center][]

往期TSP回顾  
[基于组合遗传粒子群算法的旅行商问题求解][Link 2]  
[遗传算法求最短路径（旅行商问题）python实现][python]

![在这里插入图片描述][20210331155102422.jpg_pic_center]  
作者：电气余登武  
![2021033115513752.jpg_pic_center][]

[20210331142440391.png_pic_center]: /images/20221022/2ac250cb29844d1880dda3144262aff1.png
[Link 1]: https://blog.csdn.net/KOBEYU652453/article/details/115298847
[20210331134548488.png_pic_center]: /images/20221022/021726fd2019490d82351eec070fda10.png
[20210331154637247.png_pic_center]: https://img-blog.csdnimg.cn/20210331154637247.png#pic_center
[20210331154823160.png_pic_center]: /images/20221022/334f9ac72b5644ddbecad9d4a1a03077.png
[Link 2]: https://blog.csdn.net/KOBEYU652453/article/details/114830751
[python]: https://blog.csdn.net/KOBEYU652453/article/details/114787101
[20210331155102422.jpg_pic_center]: /images/20221022/98d9cbc60d194df48ac46757fa502642.png
[2021033115513752.jpg_pic_center]: https://img-blog.csdnimg.cn/2021033115513752.jpg#pic_center