粒子群算法实例_Python简易流程版

本篇文章以实现如下需求为例，用Python实现粒子群算法：求解 y=sin(10πx）/x x在[1,2] 之间的最大值。所展示代码无脑复制粘贴即可运行。

一、导入第三方库

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

上面都是常用的第三方库，你既然都开始看算法了，这几个库应该非常熟悉了。

二、初始化粒子群算法的相关参数

c1 = 1.41
c2 = 1.41
epochs = 50   # 迭代次数
pop_size = 10  # 个体数

Vmax = 0.5  # 速度的最大值
Vmin = -0.5  # 速度的最小值
pop_max = 2  # x的最大值
pop_min = 1  # x的最小值

上述参数中，c1，c2的值可以随意设置。

三、定义目标函数

# 定义函数，x为数组
def fun(x):
    if type(x) == list:
        a = np.dot(10*np.pi,x)
        b = np.sin(a) / x
        return list(b)
    else:
        return np.sin(10 * np.pi * x) / x

这里定义这个函数的主要作用是求解方程$$y=sin(10\ast \pi \ast x）/x$$的值。这里做了一个判断，如果输入的是列表，用np计算，如果是数值，则直接运算。

四、初始化粒子数和速度

pop,V = [],[]
for i in range(pop_size):
    pop.append((random() + 1) / 2 +1)
    V.append(0.5 * random())
    fitness = fun(pop)  # 计算每个初始化个体所对应的函数值

这里是对粒子数和运动速度进行初始化话。
比如，pop（粒子数）初始化如下：

V 初始化如下：

fitness 是初始化的 pop 值根据 fun（x）运算后的结果：

五、挑选个体最优解和全局最优解

bestfitness = max(fitness)   # 10 个随机的个体中，最大的函数值
bestindex = fitness.index(bestfitness)  #个体中，最大函数值对应的索引
zbest = pop[bestindex]  # 得到最大函数值，所对应的个体，即哪个个体的函数值最大
gbest = pop    # 存储所有的个体值
fitnessgbest = fitness  # 存储 10 个个体所对应的函数值
fitnesszbest = bestfitness # 存储最大的函数值

• bestfitness 是保存初始化中 函数值最大的数：
  
• bestindex 是保存函数值最大数的索引值：
  
• zbest 是得到最大函数值的 x 的值。
  
• gbest 是存储所有的 pop 值
  
• fitnessgbest 表示的是10个粒子数所求得的函数值
  
• fitnesszbest 表示的是这10个粒子群中对应最大的函数值。
  

六、迭代优化

yy = []
for i in range(epochs):
    temp = []
    for j in range(pop_size):
        # 速度更新
        V[j] = V[j] + c1*random()*(gbest[j]-pop[j]) + c2*random()*(zbest - pop[j])
        if V[j] > Vmax: V[j] = Vmax
        if V[j] < Vmin: V[j] = Vmin

        # 个体更新
        pop[j] = pop[j] + V[j]
        if pop[j] > pop_max: pop[j] = pop_max
        if pop[j] < pop_min: pop[j] = pop_min

        # 全局最优值更新
        temp.append(fun(pop[j]))
        # fitness[j] = fun(pop[j])

    for j in range(pop_size):
        # 个体最优值
        if temp[j] > fitnessgbest[j]:
            gbest[j] = pop[j]
            fitnessgbest[j] = temp[j]

        # 全局最优值
        if temp[j] > fitnesszbest:
            zbest = pop[j]
            fitnesszbest = temp[j]
    yy.append(fitnesszbest)

这里受篇幅影响就不在具体展示，大家可以自己debug一下，其思想就是：

• 第一个小循环主要是将每个粒子进行更新后，再次计算函数值，并进行更新；
• 第二个小循环主要是用更新后的函数值与之前的函数值进行比较，如果大于粒子数之前的函数值，就替换所对应的值，如果大于全局最优值，则将全局最优值进行替换。

这里需要注意的是，在第一个小循环上要设置 temp 来接受新的函数值，不能拿之前的fitness 列表接受，因为 “=” 是浅拷贝。

七、可视化图像

x = np.arange(1,2,0.01)
y = np.sin(10 * np.pi * x) / x
plt.subplot(1,2,1)
plt.scatter(x=zbest,y=fitnesszbest,c='b',marker='*')
plt.plot(x,y,'r-')

plt.subplot(1,2,2)
iters = np.arange(1,epochs+1)
plt.plot(iters,yy,'r-')
plt.show()

具体运行结果如下图：