【建模算法】基于粒子群算法求解TSP问题（Python实现）

【建模算法】基于粒子群算法（PSO）求解TSP问题（Python实现）

TSP (traveling salesman problem，旅行商问题)是典型的NP完全问题，即其最坏情况下的时间复杂度随着问题规模的增大按指数方式增长，到目前为止还未找到一个多项式时间的有效算法。本文探讨了基于粒子群算法求解TSP问题的Python实现。

一、问题描述

​ 本案例以31个城市为例，假定31个城市的位置坐标如表1所列。寻找出一条最短的遍历31个城市的路径。

城市编号	X坐标	Y坐标	城市编号	X坐标	Y坐标
1	1.304	2.312	17	3.918	2.179
2	3.639	1.315	18	4.061	2.37
3	4.177	2.244	19	3.78	2.212
4	3.712	1.399	20	3.676	2.578
5	3.488	1.535	21	4.029	2.838
6	3.326	1.556	22	4.263	2.931
7	3.238	1.229	23	3.429	1.908
8	4.196	1.044	24	3.507	2.376
9	4.312	0.79	25	3.394	2.643
10	4.386	0.57	26	3.439	3.201
11	3.007	1.97	27	2.935	3.24
12	2.562	1.756	28	3.14	3.55
13	2.788	1.491	29	2.545	2.357
14	2.381	1.676	30	2.778	2.826
15	1.332	0.695	31	2.37	2.975
16	3.715	1.678

二、粒子群算法主要步骤

1. 初始化粒子群。包括粒子的初始位置及速度，惯性因子等参数值

粒子数M一般选取20~40个，不过对于一些特殊的难题需要更多的粒子，粒子数量越多，搜索范围就越广，越容易找到全局最优解，但是也会带来更大的计算消耗。

2. 评价各个粒子的初始适应值。

3. 将初始的适应值作为各个粒子的局部最优解，保存各粒子的最优位置。并找到其中的最优值，作为全局最优解的初值，并记录其位置

4. 更新粒子速度及位置

5. 计算更新后粒子的适应值，更新每个粒子的局部最优值以及整个粒子群的全局最优值。

6. 重复4~5直至满足迭代结束条件

   流程图如下：
   

三、粒子群算法求解TSP问题

1.TSP问题

**旅行商问题（Traveling Salesman Problem，TSP）**是一个典型的NP问题。

问题描述：一个销售员需要访问31个城市，恰好访问每个城市一次，并最终回到出发城市，销售员希望整个旅行的路程最短（或者成本定义为旅行所需的全部费用是他旅行经过的的各边费用之和，而销售员希望使整个旅行费用最低）。从图论的角度来看，该问题实质是在一个带权完全无向图中，找一个权值最小的Hamilton回路。

2.位置更新（position update）

在TSP问题中，我们所求的是一个序列，需要优化的是这个序列的顺序，粒子$x_i$ 如何定义是关键，我们将其定义为当前粒子 i所走的顺序序列，如：

$x_i=(1,3,5,2,4)$代表着该TSP问题一共有五个城市需要访问，当前粒子i的访问顺序为1→3→5→2→4，

因此全部城市的访问顺序的所有可能序列就构成了问题的可行空间，**粒子位置的更新就意味着粒子$x_i$ 从一种序列$s_i$ 变化成另外一种序列 $s_j$。

3.速度更新（velocity update）

四、求解结果

最优路线与最优值：

最优轨迹图：

优化过程：

当然，PSO是一个基于概率的随机自搜索算法，每次的搜索结果都不尽相同，但是若算法的收敛性表现的较好，也可以认为对算法的设计是合理的。

五、Python源代码

#粒子群算法求解31座城市TSP问题完整代码：
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import time

#计算距离矩阵
def clac_distance(X, Y):
    """
    计算两个城市之间的欧氏距离，二范数
    :param X: 城市X的坐标.np.array数组
    :param Y: 城市Y的坐标.np.array数组
    :return:
    """
    distance_matrix = np.zeros((city_num, city_num))
    for i in range(city_num):
        for j in range(city_num):
            if i == j:
                continue

            distance = np.sqrt((X[i] - X[j]) ** 2 + (Y[i] - Y[j]) ** 2)
            distance_matrix[i][j] = distance

    return distance_matrix

#定义总距离(路程即适应度值)
def fitness_func(distance_matrix, x_i):
    """
    适应度函数
    :param distance_matrix: 城市距离矩阵
    :param x_i: PSO的一个解（路径序列）
    :return:
    """
    total_distance = 0
    for i in range(1, city_num):
        start_city = x_i[i - 1]
        end_city = x_i[i]
        total_distance += distance_matrix[start_city][end_city]
    total_distance += distance_matrix[x_i[-1]][x_i[0]]  # 从最后的城市返回出发的城市

    return total_distance

#定义速度更新函数
def get_ss(x_best, x_i, r):
    """
    计算交换序列，即x2结果交换序列ss得到x1，对应PSO速度更新公式中的 r1(pbest-xi) 和 r2(gbest-xi)
    :param x_best: pbest or gbest
    :param x_i: 粒子当前的解
    :param r: 随机因子
    :return:
    """
    velocity_ss = []
    for i in range(len(x_i)):
        if x_i[i] != x_best[i]:
            j = np.where(x_i == x_best[i])[0][0]
            so = (i, j, r)  # 得到交换子
            velocity_ss.append(so)
            x_i[i], x_i[j] = x_i[j], x_i[i]  # 执行交换操作

    return velocity_ss

# 定义位置更新函数
def do_ss(x_i, ss):
    """
    执行交换操作
    :param x_i:
    :param ss: 由交换子组成的交换序列
    :return:
    """
    for i, j, r in ss:
        rand = np.random.random()
        if rand <= r:
            x_i[i], x_i[j] = x_i[j], x_i[i]
    return x_i

def draw(best):
    result_x = [0 for col in range(city_num+1)]
    result_y = [0 for col in range(city_num+1)]
    
    for i in range(city_num):
        result_x[i] = city_x[best[i]]
        result_y[i] = city_y[best[i]]
    result_x[city_num] = result_x[0]
    result_y[city_num] = result_y[0]
    plt.xlim(0, 5)  # 限定横轴的范围
    plt.ylim(0, 4)  # 限定纵轴的范围
    plt.plot(result_x, result_y, marker='>', mec='r', mfc='w',label=u'路线')
    plt.legend()  # 让图例生效
    plt.margins(0)
    plt.subplots_adjust(bottom=0.15)
    for i in range(len(best)):
        plt.text(result_x[i] + 0.05, result_y[i] + 0.05, str(best[i]+1), color='red')
    plt.xlabel('横坐标')
    plt.ylabel('纵坐标')
    plt.title('轨迹图')
    plt.show()
     
def print_route(route):
    result_cur_best=[]
    for i in route:
        result_cur_best+=[i]
    for i in range(len(result_cur_best)):
        result_cur_best[i] += 1
    result_path = result_cur_best
    result_path.append(result_path[0])
    return result_path

if __name__=="__main__":
    
    #读取31座城市坐标
    coord = []
    with open("data.txt", "r") as lines:
        lines = lines.readlines()
    for line in lines:
        xy = line.split()
        coord.append(xy)
    coord = np.array(coord)
    w, h = coord.shape
    coordinates = np.zeros((w, h), float)
    for i in range(w):
        for j in range(h):
            coordinates[i, j] = float(coord[i, j])
    city_x=coordinates[:,0]
    city_y=coordinates[:,1]
    #城市数量
    city_num = coordinates.shape[0]
    
    # 参数设置
    size = 500       #粒子数量
    r1 = 0.7         #个体学习因子
    r2 = 0.8         #社会学习因子
    iter_max_num = 1000    #迭代次数
    fitness_value_lst = []

    distance_matrix = clac_distance(city_x, city_y)

    # 初始化种群各个粒子的位置，作为个体的历史最优pbest
    pbest_init = np.zeros((size, city_num), dtype=np.int64)
    for i in range(size):
        pbest_init[i] = np.random.choice(list(range(city_num)), size=city_num, replace=False)

    # 计算每个粒子对应的适应度
    pbest = pbest_init
    pbest_fitness = np.zeros((size, 1))
    for i in range(size):
        pbest_fitness[i] = fitness_func(distance_matrix, x_i=pbest_init[i])

    # 计算全局适应度和对应的gbest
    gbest = pbest_init[pbest_fitness.argmin()]
    gbest_fitness = pbest_fitness.min()

    # 记录算法迭代效果
    fitness_value_lst.append(gbest_fitness)

    # 迭代过程
    for i in range(iter_max_num):
        # 控制迭代次数
        for j in range(size):
            # 遍历每个粒子
            pbest_i = pbest[j].copy()
            x_i = pbest_init[j].copy()

            # 计算交换序列，即 v = r1(pbest-xi) + r2(gbest-xi)
            ss1 = get_ss(pbest_i, x_i, r1)
            ss2 = get_ss(gbest, x_i, r2)
            ss = ss1 + ss2
            x_i = do_ss(x_i, ss)

            fitness_new = fitness_func(distance_matrix, x_i)
            fitness_old = pbest_fitness[j]
            if fitness_new < fitness_old:
                pbest_fitness[j] = fitness_new
                pbest[j] = x_i

            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           
        fitness_value_lst.append(gbest_fitness)

    # 输出迭代结果
    print("最优路线：", print_route(gbest))
    print("最优值：", gbest_fitness)
    
    #绘图    
    sns.set_style('whitegrid')
    plt.rcParams['font.sans-serif'] = 'SimHei'  # 设置中文显示
    plt.rcParams['axes.unicode_minus'] = False
    draw(gbest)
    plt.figure(2)
    plt.plot(fitness_value_lst)
    plt.title('优化过程')
    plt.ylabel('最优值')
    plt.xlabel('迭代次数({}->{})'.format(0, iter_max_num))
    plt.show()