TSP问题——PSO粒子群算法（代码+个人理解）

 大致的原理可以参考这篇 文章 简单了解一下，我就直接进入代码环节了

首先来看一个简单的应用：求x,y,z的最小平方和

import numpy as np


def fitness_value(bird_site, dimension):
    value = 0
    for i in range(dimension):
        value += bird_site[i]**2
    return value


def create_birds(pop_size, dimension, x_max, v_max):
    ibest_site = np.zeros((pop_size, dimension))
    ibest_value = [0] * pop_size
    best_value = 1000000
    birds_site = np.zeros((pop_size, dimension))
    birds_speed = np.zeros((pop_size, dimension))
    birds_fitness = [0] * pop_size
    for i in range(pop_size):
        birds_site[i] = np.random.uniform(-x_max, x_max, (1, dimension))
        birds_speed[i] = np.random.uniform(-v_max, v_max, (1, dimension))
        birds_fitness[i] = fitness_value(birds_site[i], dimension)
        ibest_site[i] = birds_site[i]
        ibest_value[i] = birds_fitness[i]
        if birds_fitness[i] < best_value:
            best_value = birds_fitness[i]
            best_site = birds_site[i]

    return birds_site, birds_speed, birds_fitness, ibest_site, ibest_value, best_value, best_site


def update(birds_site, birds_speed, birds_fitness, ibest_site, ibest_value,
           best_value, best_site, w, c1, c2, x_max, v_max, update_number, hope_value):
    for count in range(update_number):
        for i in range(len(ibest_value)):
            for j in range(len(best_site)):
                birds_site[i][j] += birds_speed[i][j]
                if birds_site[i][j] > x_max:
                    birds_site[i][j] = x_max
                if birds_site[i][j] < -x_max:
                    birds_site[i][j] = -x_max
            birds_fitness[i] = fitness_value(birds_site[i], len(best_site))
            if birds_fitness[i] < ibest_value[i]:
                ibest_value[i] = birds_fitness[i]
                ibest_site[i] = birds_site[i]
                if ibest_value[i] < best_value:
                    best_value = ibest_value[i]
                    best_site = ibest_site[i]
                    if best_value < hope_value:
                        return best_value, best_site

        """*********更新速度********"""
        for i in range(len(ibest_value)):
            for j in range(len(best_site)):
                # 可使用Clerc收缩因子
                birds_speed[i][j] = w * birds_speed[i][j] + c1 * np.random.rand() * (ibest_site[i][j] - birds_site[i][j]) \
                                    + c2 * np.random.rand() * (best_site[j] - birds_site[i][j])
                if birds_speed[i][j] > v_max:
                    birds_speed[i][j] = v_max
                if birds_speed[i][j] < -v_max:
                    birds_speed[i][j] = -v_max

    return best_value, best_site




def main():
    dimension = 3
    """********设置超参数********"""
    pop_size = 10
    w = 1    # 速度记忆
    c1 = 2   # 自我程度
    c2 = 2   # 社会程度

    """********设置实际范围*********"""
    x_max = 100    # 某一维度最大范围（也可分开设置界限）
    v_max = 50
    update_number = 1000
    hope_value = 0.0001  # 期望停止值
    birds_site, birds_speed, birds_fitness, ibest_site, ibest_value, best_value, best_site = create_birds(pop_size, dimension, x_max, v_max)
    best_value, best_site = update(birds_site, birds_speed, birds_fitness, ibest_site, ibest_value, best_value,
                                   best_site, w, c1, c2, x_max, v_max, update_number, hope_value)
    print(best_value, "\n",  best_site)





if __name__=="__main__":
    main()

运行结果： 

对于实数范围，粒子群算法拥有不错的效果，因此能够更加精确的向全局最优解收敛。

但是在TSP问题中，每一个个体的变量要求是整型数，把每个粒子的位置和速度直接换成整型肯定是不行的，你就想想一张开放世界的大地图只能固定传送到里面几个点，这样能出最优解就见了鬼了。所以粒子的位置依然是用浮点数表示，然后我们给每一个粒子标记上自己对应的city的编号，再根据粒子的位置的绝对值大小来排序，这样同时也实现了对TSP路径的规划。

import random
import numpy as np
import matplotlib.pyplot as plt

# randomly generate the map with constraint of [-100, 100]
def gen_cities(city_num, random_state=True):
    if random_state:
        cities = (np.random.uniform(0, 2, (city_num, 2))-1)*200
        np.savetxt('cities.txt', cities)
    else:
        cities = np.loadtxt('cities.txt')
    return cities


def fitness_value(bird_site, dimension, distances):
    route = bird_site.copy()
    value = 0
    index = [0] * dimension
    for i in range(dimension):
        minimum = 10000
        for j in range(dimension):
            if route[j] < minimum:
                minimum = route[j]
                index[i] = j
        route[index[i]] = 10000

    for i in range(dimension-1):
        value += distances[index[i]][index[i+1]]
    value += distances[0][index[dimension-1]]
    return value


def get_distance_mateix(cities):
    distances = np.zeros((len(cities), len(cities)))
    for i in range(len(cities)):
        for j in range(len(cities)):
            distances[i][j] = np.sqrt((cities[i][0]-cities[j][0])**2 + (cities[i][1]-cities[j][1])**2)

    return distances


def create_birds(pop_size, dimension, x_max, v_max, distances):
    ibest_site = np.zeros((pop_size, dimension))
    ibest_value = [0] * pop_size
    best_value = 1000000
    birds_site = np.zeros((pop_size, dimension))
    birds_speed = np.zeros((pop_size, dimension))
    birds_fitness = [0] * pop_size
    for i in range(pop_size):
        birds_site[i] = np.random.uniform(-x_max, x_max, (1, dimension))
        birds_speed[i] = np.random.uniform(-v_max, v_max, (1, dimension))
        birds_fitness[i] = fitness_value(birds_site[i], dimension, distances)
        ibest_site[i] = birds_site[i]
        ibest_value[i] = birds_fitness[i]
        if birds_fitness[i] < best_value:
            best_value = birds_fitness[i]
            best_site = birds_site[i]

    return birds_site, birds_speed, birds_fitness, ibest_site, ibest_value, best_value, best_site

# get distance matrix



# plot cities on a graph
def plot_map(cities, path):
    plt.figure(figsize=(10,10), dpi=100)
    plt.scatter(cities[:,0], cities[:,1], s=150)
    for i in range(len(cities)-1):
        plt.plot([cities[path[i]][0], cities[path[i+1]][0]], [cities[path[i]][1], cities[path[i+1]][1]], c='r')
    plt.plot([cities[path[0]][0], cities[path[len(cities)-1]][0]], [cities[path[0]][1], cities[path[len(cities)-1]][1]], c='r')
    plt.savefig('cities.png')

# find fittest path called every generation
def find_fittest(routes_length, pop_size):
    key = 10000000000
    fittest = 0
    for i in range(pop_size):
        if routes_length[i] < key:
            key = routes_length[i]
            fittest = i
    return fittest

def update(birds_site, birds_speed, birds_fitness, ibest_site, ibest_value,
           best_value, best_site, w, c1, c2, x_max, v_max, update_number, hope_value, distances):
    for count in range(update_number):
        for i in range(len(ibest_value)):
            for j in range(len(best_site)):
                birds_site[i][j] += birds_speed[i][j]
                if birds_site[i][j] > x_max:
                    birds_site[i][j] = x_max
                if birds_site[i][j] < -x_max:
                    birds_site[i][j] = -x_max
            birds_fitness[i] = fitness_value(birds_site[i], len(best_site), distances)
            if birds_fitness[i] < ibest_value[i]:
                ibest_value[i] = birds_fitness[i]
                ibest_site[i] = birds_site[i]
                if ibest_value[i] < best_value:
                    best_value = ibest_value[i]
                    best_site = ibest_site[i]
                    if best_value < hope_value:
                        return best_value, best_site

        """*********更新速度********"""
        for i in range(len(ibest_value)):
            for j in range(len(best_site)):
                # 可使用Clerc收缩因子
                birds_speed[i][j] = w * birds_speed[i][j] + c1 * np.random.rand() * (ibest_site[i][j] - birds_site[i][j]) \
                                    + c2 * np.random.rand() * (best_site[j] - birds_site[i][j])
                if birds_speed[i][j] > v_max:
                    birds_speed[i][j] = v_max
                if birds_speed[i][j] < -v_max:
                    birds_speed[i][j] = -v_max

    return best_value, best_site




def PSO(cities, distances, city_num):
    dimension = city_num
    """********设置超参数********"""
    pop_size = 30
    w = 1    # 速度记忆
    c1 = 2   # 自我程度
    c2 = 2   # 社会程度

    """********设置实际范围*********"""
    x_max = 100    # 某一维度最大范围（也可分开设置界限）# 地图大小
    v_max = 30
    update_number = 1000
    hope_value = 100  # 期望停止值
    birds_site, birds_speed, birds_fitness, ibest_site, ibest_value, best_value, best_site = create_birds(pop_size, dimension, x_max, v_max, distances)
    best_value, best_site = update(birds_site, birds_speed, birds_fitness, ibest_site, ibest_value, best_value,
                                   best_site, w, c1, c2, x_max, v_max, update_number, hope_value, distances)
    print(best_value)
    index = [0] * dimension
    for i in range(dimension):
        minimum = 10000
        for j in range(dimension):
            if best_site[j] < minimum:
                minimum = best_site[j]
                index[i] = j
        best_site[index[i]] = 10000

    print(index)
    plot_map(cities, index)


def main():
    city_num = 8
    cities = gen_cities(city_num, random_state=True)
    distances = get_distance_mateix(cities)
    PSO(cities, distances, city_num)

if __name__=="__main__":
    main()
    plt.show()

代码直接用就行，成功率极高

运行结果（8点与15点，粒子数量30，其他超参数默认）：

只需要代码内容的直接复制上面就可以了，以下是我对粒子群算法在PSO问题的缺陷的探讨

首先，每一个粒子（城市）对应的适应度是随机的，也就是说适应度对于这个粒子本身来说是无意义的，但是对与其他粒子其实是有联系的。现在我创造一个定义：

粒子b对粒子a的适应度为 |f(a)-f(b)| / f(a)

（f(x)就是每个粒子自己的位置x，这个“粒子b对粒子a的适应度为”的定义形式可以是任意的，只要能表示粒子a，b之间在地图上的距离就行，这个地图就是你自己设的x_max）

把每一个粒子看成一只鸟，TSP问题在PSO中相当于每一只鸟在飞行中寻找最适合自己的位置排成一列.（如下图）假设在一个8只鸟组成的鸟群中a,b,c是一个最好的顺序，对于a鸟来说，它找到了最适合自己的一个位置，如果这个时候飞过b,c,d三只鸟，a鸟利用了|f(a)-f(b)| / f(a)

这个式子判断出了b,c鸟应该排在后面两个位置，这样d鸟就不能把b,c的位置占掉了，然而b,c鸟排在a鸟后面之后，其他的鸟还没有排好，因此b,c鸟只能在a鸟后面两个位置乱飞。有时候是a,b,c，有时候是a,c,b，但是在|f(a)-f(b)| / f(a)里面，a鸟眼里b和c长得都很像，所以a鸟不能百分百确定到底是a,b,c正确还是a,c,b正确，最终变成了a,c,b。于是出现了下面的一个交叉三角形

而下面这个图是两个交叉出来的三角形，其实就是两个abc的序列像上面的方式排错位置了（为方便观看就分两次看）

  

第一个序列abc排成了acb

第二个序列abc排成了acb

有人可能会发现第一个序列中b既没有和a相连，又没有和c相连，那为什么叫acb呢？因为a后面第一个是c。这个时候就需要第二个序列帮忙了，也就是说纠正这个错误至少需要两个acb序列改正成为abc序列：

      首先a1接上b1（1号线），这个时候a1接了三条线多了一条，3号线断开，为什么不断d点那一条呢，因为d点不在序列a1b1c1也不在a2b2c2，所以d点是正确点不用理会。这个时候a2也接了三条，a1点因为已经重新接过了所以成为了修改点，所以断4号线。然后a2接上b2（2号线），a2又多了一条，同理只能断5号线了。至此，两个acb序列均已被改正成abc序列。这个时候再观察到c2点和e点都只连了一条线，它们自然就连上了。（断的优先级：修改点<正确点<错误点）

于是就有了全局最优解。

在鸟的数量多的时候还会出现交叉的四边形：

把中间两个点看成一个点其实就是交叉三角形了。

上面就是 我对我的猜想的自证，这就是我想说的PSO在解决TSP问题中的弊端：

       由于个体速度具有很强的随机性，因此导致个体差距变得非常模糊，让鸟群在每一次飞行的位置更换的操作并不直接，会浪费掉很多步操作导致收敛的速度太慢。遗传算法因为每一步的操作都是直接有意义的，所以才要控制它的收敛速度。

这算是我脑洞大开自己胡思乱想的吧，有理解不当的欢迎大佬在评论区指出。