用python3实现粒子群优化算法（PSO）

梦里梦外; 2022-05-24 02:45 650阅读 0赞

**【注】：完整源码在我的[github][]上，找不找得到就看自己咯 ^ \_ ^**

粒子群优化算法（Particle Swarm Optimization，PSO）属于进化算法的一种，是通过模拟鸟群捕食行为设计的。从随机解出发，通过迭代寻找最优解，通过适应度来评价解的品质。设想这样一个场景：一群鸟在随机搜索食物。在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻目前离食物最近的鸟的周围区域。  
所有的粒子具有以下两个属性：速度、位置。  
PSO 初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个"极值"来更新自己。第一个就是粒子本身所找到的最优解pbest。另一个极值是整个种群目前找到的最优解，即全局极值gbest。粒子通过下面的公式来更新自己的速度和位置：

速度：  
 v i + 1 = w ∗ v i + c 1 ∗ r a n d 1 ∗ ( p b e s t i − x i ) + c 2 ∗ r a n d 2 ∗ ( g b e s t i − x i ) v\_\{i+1\}=w\*v\_i + c\_1\*rand\_1\*(pbest\_i-x\_i) + c\_2\*rand\_2\*(gbest\_i-x\_i) vi\+1=w∗vi\+c1∗rand1∗(pbesti−xi)\+c2∗rand2∗(gbesti−xi)

位置：  
 x i = x i + v i + 1 x\_i=x\_i+v\_\{i+1\} xi=xi\+vi\+1

以上两式中：  
 w w w为惯性因子，一般取1  
 c 1 、 c 2 c\_1、 c\_2 c1、c2为学习因子，一般取2  
 r a n d 1 、 r a n d 2 rand\_1、rand\_2 rand1、rand2为两个（0,1）之间的随机数  
 v i 和 x i v\_i和x\_i vi和xi分别表示粒子第i维的速度和位置  
 p b e s t i 、 g b e s t i pbest\_i、gbest\_i pbesti、gbesti分别表示某个粒子最好位置第i维的值、整个种群最好位置第i维的值

注意：以上两式是针对粒子的某一维进行跟新的。对粒子的每一维，都要用上述的式子进行更新

说了这么多，下面把主要部分代码上上来（不上全是因为伸手党太多，哈哈）：

import numpy as np
    import random
    
    class PSO:
        def __init__(self, dim, size, iter_num, x_max, max_vel, best_fitness_value=float('Inf'), C1 = 2, C2 = 2, W = 1):
            self.C1 = C1
            self.C2 = C2
            self.W = W
            self.dim = dim  # 粒子的维度
            self.size = size  # 粒子个数
            self.iter_num = iter_num  # 迭代次数
            self.x_max = x_max
            self.max_vel = max_vel  # 粒子最大速度
            self.best_fitness_value = best_fitness_value
            self.best_position = [0.0 for i in range(dim)]  # 种群最优位置
            self.fitness_val_list = []  # 每次迭代最优适应值
    
            # 对种群进行初始化
            self.Particle_list = [Particle(self.x_max, self.max_vel, self.dim) for i in range(self.size)]
    
        # 更新速度
        def update_vel(self, part):
            for i in range(self.dim):
                vel_value = self.W * part.get_vel()[i] + self.C1 * random.random() * (part.get_best_pos()[i] - part.get_pos()[i]) \
                            + self.C2 * random.random() * (self.get_bestPosition()[i] - part.get_pos()[i])
                if vel_value > self.max_vel:
                    vel_value = self.max_vel
                elif vel_value < -self.max_vel:
                    vel_value = -self.max_vel
                part.set_vel(i, vel_value)
    
        # 更新位置
        def update_pos(self, part):
            for i in range(self.dim):
                pos_value = part.get_pos()[i] + part.get_vel()[i]
                part.set_pos(i, pos_value)
            value = fit_fun(part.get_pos())
            if value < part.get_fitness_value():
                part.set_fitness_value(value)
                for i in range(self.dim):
                    part.set_best_pos(i, part.get_pos()[i])
            if value < self.get_bestFitnessValue():
                self.set_bestFitnessValue(value)
                for i in range(self.dim):
                    self.set_bestPosition(i, part.get_pos()[i])
    
        def update(self):
            for i in range(self.iter_num):
                for part in self.Particle_list:
                    self.update_vel(part)  # 更新速度
                    self.update_pos(part)  # 更新位置
                self.fitness_val_list.append(self.get_bestFitnessValue())  # 每次迭代完把当前的最优适应度存到列表
            return self.fitness_val_list, self.get_bestPosition()

在上述代码中，跟粒子的速度设置了最大值，确保粒子速度的变化在某一范围内。接下来，我们可以用以上代码来进行测试了，这里测试用的是Hölder table function测试函数，在我的另一篇博客[用python绘制评估优化算法性能的测试函数][python]里有提到，其图像如下：

![这里写图片描述][70]

它的最优解有如下4个：  
![这里写图片描述][70 1]

接下了我们通过以下代码进行测试：

from OptAlgorithm.PSO import PSO
    import matplotlib.pyplot as plt
    import numpy as np
    
    
    dim = 2
    size = 20
    iter_num = 1000
    x_max = 10
    max_vel = 0.5
    
    pso = PSO(dim, size, iter_num, x_max, max_vel)
    fit_var_list, best_pos = pso.update()
    print("最优位置:" + str(best_pos))
    print("最优解:" + str(fit_var_list[-1]))
    plt.plot(np.linspace(0, iter_num, iter_num), fit_var_list, c="R", alpha=0.5)
    plt.show()

以上代码中，设置粒子维度为2维，种群大小即粒子个数为20，迭代次数为1000，粒子初始化的位置在（-10,10）之间，粒子的最大速度为0.5

以下是一次算法的收敛图像结果，：  
![这里写图片描述][70 2]

打印的最优解和适应度如下：  
![这里写图片描述][70 3]

可以看出其最优位置是函数四个最优解里面的一个，运行多次会出现不同的结果。甚至可能会得到如下图所示的函数收敛图像：  
![这里写图片描述][70 4]  
注意看此时函数的适应度在-9.5左右，继续迭代也没有发生变化，此时即陷入了局部最优。这篇博客主要介绍粒子群算法，这个情况这里就不继续讨论了。

[github]: https://github.com/TimePickerWang
[python]: https://blog.csdn.net/wang454592297/article/details/80336753
[70]: /images/20220524/3557889f251a4954a00504e26dda02aa.png
[70 1]: /images/20220524/d6271df475984f7682f0f9ec4f67f730.png
[70 2]: /images/20220524/a6d5a740694845a898d1bd07ec131a5b.png
[70 3]: /images/20220524/10f8924e37e04eafaf0c75eb3b9a9682.png
[70 4]: /images/20220524/04dcc0b1841347628f05109299c01ab5.png