算法原理

      旅行商问题是一个经典的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)
      解释了交换子和交换序列后，重新定义粒子的速度和更新公式：

      式中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之间的随机数，表示一定概率进行交换子操作。

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

算例

假设随机生成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()

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

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

作者：电气余登武