一篇文章搞懂PSO（粒子群算法）理论讲解＋Python代码 ...

kirin77

写在前面，我在自学PSO的过程中，发现很多知乎文章都是片段化的，没有统一起来一起讲，讲理论的讲理论，讲实例的讲实例。而且很多都是讲单目标无约束问题的优化，对于有约束问题的最优化问题却没有找到一篇怎张能够讲的很清楚的，本文将在粒子群算法理论结合带有约束问题的最优化问题利用Python代码进行讲解。其中有约束问题可以通过
阅读上面的罚函数法，写的非常好。进而转换成无约束问题，所以我们的实例中举一个分别用连续性PSO和离散型PSO算法解无约束最优化的问题以及用离散型算法求解0-1背包问题。欢迎大家评论区交流。当然如果有具体的优化问题需要让我来解决，欢迎大家付费咨询
一、粒子群算法理论

1.首先了解粒子群算法能够解决什么样的问题？

粒子群算法是一种解决最优化问题的通用方法，目前在数学上还没有给出严格的证明，但是它的实际效果具有求解速度快，数值相对稳定，算法简单等特点，是在数学建模中常用的一种算法。它主要分为离散型粒子群算法和连续性粒子群算法，其中连续性粒子群算法可以解决我们通常遇到的函数最优化问题。离散型粒子群算法可以解决像0-1背包问题这种动态规划问题。
2.粒子群算法仿生理论

鸟类在捕食过程中，鸟群成员可以通过个体之间的信息交流与共享获得其他成员的发现与飞行经历。在食物源零星分布并且不可预测的条件下，这种协作机制所带来的优势是决定性的，远远大于对食物的竞争所引起的劣势。粒子群算法受鸟类捕食行为的启发并对这种行为进行模仿，将优化问题的搜索空间类比于鸟类的飞行空间，将每只鸟抽象为一个粒子，粒子无质量、无体积，用以表征问题的一个可行解，优化问题所要搜索到的最优解则等同于鸟类寻找的食物源。粒子群算法为每个粒子制定了与鸟类运动类似的简单行为规则，使整个粒子群的运动表现出与鸟类捕食相似的特性，从而可以求解复杂的优化问题。
3.粒子群算法建模

粒子群优化算法源自对鸟群捕食行为的研究：一群鸟在区域中随机搜索食物，所有鸟知道自己当前位置离食物多远（这个距离就是适应度函数值），那么搜索的最简单有效的策略就是搜寻目前离食物最近的鸟的周围区域。粒子群算法利用这种模型得到启示并应用于解决优化问题。在粒子群算法中，每个优化问题的潜在解都是搜索空间中的一只鸟，称之为粒子（在连续性例子群算法中，每一个粒子就代表一个可行解）。所有的粒子都有一个由被优化的函数决定的适应度值，每个粒子还有一个速度决定它们飞翔的方向和距离。然后，粒子们就追随当前的最优粒子在解空间中搜索。
粒子群算法首先在给定的解空间中随机初始化粒子群，待优化问题的变量数决定了解空间的维数（这里表示每个粒子的维度，后面我们会用D表示维度）。每个粒子有了初始位置与初始速度（x表示当前位置，v表示速度，粒子群算法最重要的是对速度公式的更新，因为速度决定了更新的方向，进而决定解的值），然后通过迭代寻优。在每一次迭代中，每个粒子通过跟踪两个“极值”来更新自己在解空间中的空间位置与飞行速度：一个极值就是单个粒子本身在迭代过程中找到的最优解粒子，这个粒子叫作个体极值；（这里的意思是，每个粒子具有记忆性，每个粒子子记得两个值，一个是全部的粒子距离目标的最优值，另一个值是这个粒子曾经历史的最优值）另一个极值是种群所有粒子在迭代过程中所找到的最优解粒子，这个粒子是全局极值。
4.粒子群算法的特点

粒子群算法本质是一种随机搜索算法，它是一种新兴的智能优化技术。该算法能以较大概率收敛于全局最优解。实践证明，它适合在动态、多目标优化环境中寻优，与传统优化算法相比，具有较快的计算速度和更好的全局搜索能力
另外粒子群算法具有并行性，而且由于每个粒子具有记忆性，因此在算法迭代结束时，不仅可以得到最优值，还可以得到次优值。因此对于调度和决策问题可以给出多种有意义的方案。另外，粒子群算法对种群大小不敏感，性能下降也不是很大。
二、基本粒子群算法

假设在一个D维的目标搜索空间中，有N个粒子组成一个群落，其中第i个粒子表示为一个D维的向量：


第i个粒子的“飞行”速度也是一个D维的向量，记为


第i个粒子迄今为止搜索到的最优位置称为个体极值，记为


整个粒子群迄今为止搜索到的最优位置为全局极值，记为


在找到这两个最优值时，粒子根据如下的式（6.5）和式（6.6）来更新自己的速度和位置



这个公式相当重要，为PSO核心部分

其中：c1和c2为学习因子，也称加速常数；r1和r2为［0，1］范围内的均匀随机数；j=1，2，…，D；vij是粒子的速度，vij∈［-vmax，vmax］，vmax是常数，由用户设定来限制粒子的速度。r1和r2是介于0和1之间的随机数，增加了粒子飞行的随机性。式（6.5）右边由三部分组成：第一部分为“惯性”或“动量”部分，反映了粒子的运动“习惯”，代表粒子有维持自己先前速度的趋势；第二部分为“认知”部分，反映了粒子对自身历史经验的记忆或回忆，代表粒子有向自身历史最佳位置逼近的趋势；第三部分为“社会”部分，反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或邻域历史最佳位置逼近的趋势。
上面的是基本粒子群算法，下面的标准粒子群算法的更新是这样的，就是在速度的前面加一个权重系数。使得权重系数是一个动态调整的。保证在迭代前期保证足够的全局搜索能力，迭代后期能够专注于局部最优解的搜索能力。





还有更多的w更新公式，其他的就自己查吧，我就是这一种方法

三、离散粒子群算法

先一句话解释一下离散粒子群算法：就是将实数进行二进制编码，然后依然保持了位置速度更新公式，但是此时速度的每一维表示对应位置相应维度的取1的概率。那么速度怎么转换成概率呢？其实就是像logistic回归一样，通过logistics函数将实数映射乘0-1区间的概率值。
更专业的术语表述如下：
粒子在状态空间的取值和变化只限于0和1两个值，而速度的每一维vij代表位置每一位xij取值为1的可能性。因此，在连续粒子群中的vij更新公式依然保持不变，但是pbest和gbest只在［0，1］内取值。其位置更新等式表示如下：


其中，r是从均匀分布rand（0，1）中产生的随机数。
四、粒子群算法流程

1）初始化粒子群，包括群体规模N，每个粒子的位置xi和速度vi。
2）计算每个粒子的适应度值fit［i］。
3）对每个粒子，用它的适应度值fit［i］和个体极值pbest（i）比较。如果fit［i］＜pbest（i），则用fit［i］替换掉pbest（i）。
4）对每个粒子，用它的适应度值fit［i］和全局极值gbest比较。如果fit［i］＜gbest，则用fit［i］替换gbest。（5）迭代更新粒子的速度vi和位置xi。
6）进行边界条件处理。
7）判断算法终止条件是否满足：若是，则结束算法并输出优化结果；否则返回步骤（2）。
流程图如下所示：


五、代码示例讲解

1.使用连续PSO求解：

例6.2　求函数f（x，y）=3cos（xy）+x+y2的最小值，其中x的取值范围为［-4，4］，y的取值范围为［-4，4］
首先用Python绘制一下函数图像如下：
# 求函数f(x,y) = 3*cos(x * y) + x + y**2的最小值，其中-4 <= x <= 4, -4 <= y <= 4

# 首先绘制这个函数的三维图像
import numpy as np
import matplotlib.pyplot as plt

X = np.arange(-4 ,4 ,0.01)
Y = np.arange(-4 ,4 ,0.01)
x, y = np.meshgrid(X ,Y)
Z = 3*np.cos(x * y) + x + y**2

# 作图
fig = plt.figure(figsize=(10,15))
ax3 = plt.axes(projection = "3d")
ax3.plot_surface(x,y,Z ,cmap = "rainbow")
# ax3.contour(x ,y ,Z ,zdim = "z" ,offset=-2 ,cmap = "rainbow")
plt.show()输出结果如下：


可以看到局部极值点还是非常多的，下面是粒子群算法代码：每行代码基本上都有注释，应该算很良心了，如果不了解的欢迎评论区交流
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
# 设置字体和设置负号
matplotlib.rc("font", family="KaiTi")
matplotlib.rcParams["axes.unicode_minus"] = False
# 初始化种群，群体规模，每个粒子的速度和规模
N = 100 # 种群数目
D = 2 # 维度
T = 200 # 最大迭代次数
c1 = c2 = 1.5 # 个体学习因子与群体学习因子
w_max = 0.8 # 权重系数最大值
w_min = 0.4 # 权重系数最小值
x_max = 4 # 每个维度最大取值范围，如果每个维度不一样，那么可以写一个数组，下面代码依次需要改变
x_min = -4 # 同上
v_max = 1 # 每个维度粒子的最大速度
v_min = -1 # 每个维度粒子的最小速度


# 定义适应度函数
def func(x):
    return 3 * np.cos(x[0] * x[1]) + x[0] + x[1] ** 2


# 初始化种群个体
x = np.random.rand(N, D) * (x_max - x_min) + x_min # 初始化每个粒子的位置
v = np.random.rand(N, D) * (v_max - v_min) + v_min # 初始化每个粒子的速度

# 初始化个体最优位置和最优值
p = x # 用来存储每一个粒子的历史最优位置
p_best = np.ones((N, 1))  # 每行存储的是最优值
for i in range(N): # 初始化每个粒子的最优值，此时就是把位置带进去，把适应度值计算出来
    p_best = func(x[i, :])

# 初始化全局最优位置和全局最优值
g_best = 100 #设置真的全局最优值
gb = np.ones(T) # 用于记录每一次迭代的全局最优值
x_best = np.ones(D) # 用于存储最优粒子的取值

# 按照公式依次迭代直到满足精度或者迭代次数
for i in range(T):
    for j in range(N):
        # 个更新个体最优值和全局最优值
        if p_best[j] > func(x[j,:]):
            p_best[j] = func(x[j,:])
            p[j,:] = x[j,:].copy()
        # p_best[j] = func(x[j,:]) if func(x[j,:]) < p_best[j] else p_best[j]
        # 更新全局最优值
        if g_best > p_best[j]:
            g_best = p_best[j]
            x_best = x[j,:].copy()   # 一定要加copy，否则后面x[j,:]更新也会将x_best更新
        # 计算动态惯性权重
        w = w_max - (w_max - w_min) * i / T
        # 更新位置和速度
        v[j, :] = w * v[j, :] + c1 * np.random.rand(1) * (p[j, :] - x[j, :]) + c2 * np.random.rand(1) * (x_best - x[j, :])
        x[j, :] = x[j, :] + v[j, :]
        # 边界条件处理
        for ii in range(D):
            if (v[j, ii] > v_max) or (v[j, ii] < v_min):
                v[j, ii] = v_min + np.random.rand(1) * (v_max - v_min)
            if (x[j, ii] > x_max) or (x[j, ii] < x_min):
                x[j, ii] = x_min + np.random.rand(1) * (x_max - x_min)
    # 记录历代全局最优值
    gb = g_best
print("最优值为", gb[T - 1],"最优位置为",x_best)
plt.plot(range(T),gb)
plt.xlabel("迭代次数")
plt.ylabel("适应度值")
plt.title("适应度进化曲线")
plt.show()运行结果如下：


最优值为 -6.407604473773    最优位置为 [-3.99989325  0.75605927]
2.使用离散PSO求解：

用离散粒子群算法求函数f（x）=x+6sin（4x）+9cos（5x）的最小值，其中x的取值范围为［0，9］
函数作图代码如下：
# 用离散粒子群算法求解函数的最小值 f(x) = x + 6*np.sin(4*x) + 9*np.cos(5*x) ,0<= x <= 9
# 这是一个多个局部极值的函数
# 下面进行简单绘图

matplotlib.rc("font", family="KaiTi")
matplotlib.rcParams["axes.unicode_minus"] = False

x = np.arange(0 ,9 ,0.01)
y = x + 6*np.sin(4*x) + 9*np.cos(5*x)
plt.figure(figsize=(15,15))
plt.plot(x, y)
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("f(x) = x + 6sin(4x) + 9cos(5x)")
plt.show()运行结果如下：


可以看到，这是一个多极值的一元函数。下面用离散型粒子群算法进行求解。
# 离散粒子群算法

# 定义适应度函数
def func2(x):
    x = int("x",2)
    return x + 6*np.sin(4*x) + 9*np.cos(5*x)

# 初始化粒子群相关参数
N = 100
D = 20
T = 200
c1 = 1.5
c2 = 1.5
w_max = 0.8
w_min = 0.8
x_max = 9
x_min = 0
v_max = 10
v_min = -10
x = np.random.randint(0,2,[N, D])
v = (v_max - v_min) * np.random.rand(N,D) + v_min
vx = np.zeros_like(v)

# 初始化每个粒子的适应度值
p = x # 用来存储每个粒子的最佳位置
p_best = np.ones(N) # 用来存储每个粒子的适应度值
for i in range(N):
    p_best = func2(x[i,:])
#     p[i,:] = x[j,:]
   

# 初始化全局最优位置与最优值
x_best = np.ones(D)
for i in range(N):
    if p_best < g_best:
        g_best = p_best
        x_best = x[i,:].copy()
g_best = 100

gb = np.ones(T) # 用来存储每依次迭代的最优值
for i in range(T):
    for j in range(N):
        # 更新每个个体最优值和最优位置
        if p_best[j] > func2(x[j,;]):
            p_best[j] = func2(x[j,:])
            p[j,:] = x[j,:].copy()
        # 更新全局最优位置和最优值
        if p_best[j] > g_best:
            g_best = p_best[j]
            x_best = x[j,:].copy()
        # 计算动态惯性权重
        w = w_max - (w_max - w_min) * i / T
        # 更新速度, 因为位置需要后面进行概率判断更新
        v[j,:] = w * v[j,:] + c1 * np.random.rand(1)*(p_best[j] - x[j,:]) + c2 * np.random.rand(1)*(x_best - x[j,:])
        # 边界条件处理
        for jj in range(D):
            if (v[j,jj] > v_max) or (v[j,jj] < v_min):
                v[j,jj] = v_min + np.random.rand(1)*(v_max - v_min)
        # 进行概率计算并且更新位置
        vx[j,:] = 1/(1 + np.exp(-v[j,:]))
        for ii in range(D):
            x[j ,ii] = 1 if vx[j, ii] > np.random.rand(1) else 0
    gb = g_best
   
print("最优值为", gb[T - 1],"最优位置为",int("{}".format(x_best)))
plt.plot(range(T),gb)
plt.xlabel("迭代次数")
plt.ylabel("适应度值")
plt.title("适应度进化曲线")
plt.show()
实验结果如下：


最优值为 -10.419856432242188 最优位置为 4.371515628352764
3.离散PSO求解0-1背包问题

0-1背包问题。有N件物品和一个容量为V的背包。第i件物品的体积是c（i），价值是w（i）。求解将哪些物品放入背包可使物品的体积总和不超过背包的容量，且价值总和最大。假设物品数量为10，背包的容量为300。每件物品的体积为［95，75，23，73，50，22，6，57，89，98］，每件物品的价值为［89，59，19，43，100，72，44，16，7，64］。
我们的主要想法是初始化速度v和二进制编码的种群粒子位置x，其中1表示选择该物品，0表示不选择该物品。取适应度值为选择物品的价值总和，计算个体适应度值，当物品体积总和大于背包容量时，对适应度值进行惩罚计算（这就是罚函数的一种，下面的适应度函数回具体展示我是怎么处理当不满足条件时适应度函数的编写）
直接上代码：

# 离散粒子群算法，本质就是将用二进制编码的方式进行解决，将离散问题空间映射到连续粒子运动空间中，其他速度位置更新原则还是和之前连续性粒子群
# 算法保持一致，但是不同的是粒子在状态空间的取值和变化只限于0和1两个值，而速度表示位置取值为 1 的可能性，其表达形式和logistics回归相似
# 连续问题也可以离散化，离散化后收敛速度变快，但是运行时间变长

# 粒子群算法流程
# 1.初始化粒子群，包括群体规模N，每个粒子的位置xi 和速度vi
# 2.计算每个粒子的适应度值fit
# 3.通过适应度值取计算每个个体的适应度值以及全局适应度值
# 4.迭代更新每个粒子的速度和位置
# 5.进行边界条件处理
# 6、判断算法条件是否终止，否则返回步骤2

# 离散粒子群算法
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
# 设置字体和设置负号
matplotlib.rc("font", family="KaiTi")
matplotlib.rcParams["axes.unicode_minus"] = False

# 定义适应度函数
def func2(x):
    count = 0
    vv = 0
    for i in range(len(x)):
        if x == 1:
            vv += value
            count += volumn
    if count <= 300:
        return vv
    else:
        return  vv-1000


# 初始化粒子群相关参数
N = 100
D = 10
T = 200
c1 = 1.5
c2 = 1.5
w_max = 0.8
w_min = 0.8
# x_max = 9
# x_min = 0
v_max = 10
v_min = -10
x = np.random.randint(0, 2, [N, D])
v = (v_max - v_min) * np.random.rand(N, D) + v_min
vx = np.random.rand(N,D) # 这个是将速度转换成概率的矩阵
value = [89 ,59, 19 ,43 ,100 ,72 ,44 ,16 ,7 ,64]
volumn = [95 ,75 ,23 ,73 ,50 ,22 ,6 ,57 ,89 ,98]

# 初始化每个粒子的适应度值
p = x  # 用来存储每个粒子的最佳位置
p_best = np.ones(N)  # 用来存储每个粒子的适应度值
for i in range(N):
    p_best = func2(x[i, :])
#     p[i,:] = x[j,:]


# 初始化全局最优位置与最优值
g_best = 100
x_best = np.ones(D)
for i in range(N):
    if p_best > g_best:
        g_best = p_best
        x_best = x[i, :].copy()


gb = np.ones(T)  # 用来存储每依次迭代的最优值
for i in range(T):
    for j in range(N):
        # 更新每个个体最优值和最优位置
        if p_best[j] < func2(x[j,:]):
            p_best[j] = func2(x[j, :])
            p[j, :] = x[j, :].copy()
        # 更新全局最优位置和最优值
        if p_best[j] > g_best:
            g_best = p_best[j]
            x_best = x[j, :].copy()
        # 计算动态惯性权重
        w = w_max - (w_max - w_min) * i / T
        # 更新速度, 因为位置需要后面进行概率判断更新
        v[j, :] = w * v[j, :] + c1 * np.random.rand(1) * (p_best[j] - x[j, :]) + c2 * np.random.rand(1) * (x_best - x[j, :])
        # 边界条件处理
        for jj in range(D):
            if (v[j, jj] > v_max) or (v[j, jj] < v_min):
                v[j, jj] = v_min + np.random.rand(1) * (v_max - v_min)
        # 进行概率计算并且更新位置
        vx[j, :] = 1 / (1 + np.exp(-v[j, :]))
        for ii in range(D):
            r = np.random.rand(1)
            x[j, ii] = 1 if vx[j, ii] > r else 0
    gb = g_best

print("最优值为", gb[T - 1], "最优位置为", x_best)
plt.plot(range(T), gb)
plt.xlabel("迭代次数")
plt.ylabel("适应度值")
plt.title("适应度进化曲线")
plt.show()实验结果：


最优值为 388.0 最优位置为 [1 0 1 0 1 1 1 0 0 1]
其中的1表示选择该物品，0表示不选择该物品。
至此我们的讲解就结束了，但是我们还是总结一下
六、总结

PSO算法设计要点：
1.主要是适应度函数的改写。
（1）像0-1背包问题，当不满足约束条件时。给与适应度函数一个较大的惩罚。
（2）约束问题改写为无约束问题，主要使用罚函数法，罚函数法介绍说明主要参考知乎链接。简单说一下，主要是适应度函数的构造，
如果是不等式，那么用到的是max或者min函数。如果是等式则用到绝对值函数。具体使用什么函数则要考问题是取最大值还是最小值。
总之一句话就是在可行域之内，罚函数的值与目标函数的值一致。在可行域之外，罚函数的值根据问题的优化方向进行适当惩罚。
（3）默认算法都是往最小值方向进行优化，如果往最大值方向进行优化，只需要在算法中更改大于号小于号的方向。

Ylisar

666[机智]

DomDomm

[赞同]