粒子群算法及其Python实现

文章目录

算法原理

粒子群算法，缩写为PSO(Particle Swarm Optimization)，是一种非线性寻优算法，其特点是实现简单、收敛速度快，对多元函数的局部最优有较好的克服能力。

所谓粒子群，就是一群粒子，每个粒子都有自己的位置和速度，记第 i i i个粒子的位置为 x ⃗ i \vec x_i x i​，速度为 v ⃗ i \vec v_i v i​。如果没有任何外加条件，这群粒子的轨迹，将完全由某一时刻的位置和速度决定。

而想要通PSO进行优化的函数，则可理解为对粒子的某种约束。例如，现有一个二元函数 y = x 1 2 + x 2 2 y=x_1^2+x_2^2 y=x12​+x22​，要找到这个函数的最小值，根据当前每个粒子的位置都可以得到一个 y y y，记每个粒子的 y y y值为 y i = x i 1 2 + x i 2 2 y_i=x_{i1}^2+x_{i2}^2 yi​=xi12​+xi22​，那么就当前来说，必定有一个 y y y值最小的粒子，这个粒子会成为其他粒子的风向标，大家都纷纷向这颗粒子靠拢，以便让自己也达到最优。

另一方面，每当经历一轮迭代，第 i i i个粒子都会经历一个地方，如果粒子有记忆的话，它会记住自己的 y y y值最小时的参数，从而不必完全地随波逐流，向当前最优粒子靠拢。

也就是说，粒子的速度，将由三个因素决定，一是粒子当前的速度；二是粒子群中最优粒子的位置；三则是这个粒子曾经去过的最优的位置。

记 X ⃗ i \vec X_i X i​为第 i i i个粒子历史上最好的位置； b ⃗ j \vec b_j b j​为第 j j j次迭代得到的最好结果，则第 i i i个粒子的速度表达式可以记作

v ⃗ i , j + 1 = v ⃗ i j + c 0 r 0 ( b ⃗ j − x ⃗ i j ) + c 1 r 1 ( X ⃗ i − x ⃗ i j ) \vec v_{i,j+1} = \vec v_{ij} + c_0r_0(\vec b_{j}-\vec x_{ij})+c_1r_1(\vec X_i-\vec x_{ij}) v i,j+1​=v ij​+c0​r0​(b j​−x ij​)+c1​r1​(X i​−x ij​)

其中 c 0 , c 1 c_0,c_1 c0​,c1​分别是全局最优和历史最优系数， r 0 , r 1 r_0,r_1 r0​,r1​则是随机数。

Python实现

首先，需要创建一个类或者字典来表示粒子，这里用类来实现

import numpy as np
uniRand = np.random.uniform

class Particle:
    def __init__(self, N, xRange, vRange):
        # 生成(N)维参数
        self.x = uniRand(*xRange, (N,))
        self.v = uniRand(*vRange, (N,))
        self.best = np.inf
        self.xBest = np.zeros((N,))

接下来实现粒子群，粒子群无非是多个粒子，所以要有一个存放粒子的列表。故其初始化函数如下，其中best和xBest分别存放全局最优解和全局最优参数。

from copy import deepcopy
rand = np.random.rand
class Swarm:
    # pNum 粒子个数，N粒子维度，即参数个数
    # wRange 最大最小权重
    def __init__(self, pNum, N, xRange, vRange, wRange, c):
        self.ps = [Particle(N, xRange, vRange) for _ in range(pNum)]
        self.bestPs = deepcopy(self.ps)
        self.best = np.inf  #全局最优值
        self.xBest = np.zeros((N,)) #全局最优参数
        self.wRange = wRange
        self.c0, self.c1 = c
        self.N = N

然后将参数更新的规则封入一个函数中，这里包括三个内容的更新，首先是从当前计算结果中，挑选出最优结果，作为新的全局最优值；然后针对单个粒子，更新粒子的历史最优值；最后，也是最复杂的一步，即更新每个粒子的速度和位置。

在写好单步更新函数之后，将其写入循环中，就是粒子群算法的主干了，代码如下。

# 下面的函数写在class Swarm中
# 为了书写方便，这里顶个写，注意在Swarm中的缩进
    # 参数更新 func为待优化函数
def oneStep(self, func, w=1):
    for p in self.ps:
        y = func(p.x)
        # 更新历史最优值
        if y < p.best:
            p.best = y
            p.xBest = deepcopy(p.x)

        # 更新粒子群最优情况
        if y < self.best:
            self.best = y
            self.xBest = deepcopy(p.x)
            self.bestPs = deepcopy(self.ps)

    for p in self.ps:
        iw = uniRand(*self.wRange, 1)[0]
        us = self.c0 * rand(self.N) * (p.xBest - p.x)
        vs = self.c1 * rand(self.N) * (self.xBest - p.x)
        p.v = iw * p.v + us + vs
        p.x = p.x + p.v
    return self.best

# nMsg为输出提示的周期
# iter为迭代次数
def optimize(self, func, nMsg, iter):
    for i in range(iter):
        best = self.oneStep()
        if i % nMsg == 0:
            print(f"第{i}次迭代最小值为{best}")

得益于pyhton的函数式特性，可以方便地将函数作为参数传入，从而在优化函数optimize中，直接把被优化的函数func用作输入参数。

算法测试

最后，代码写好之后，可以测试一下，假设现有一个比较复杂的多元函数，其中 x i x_i xi​为不同维度的参数，令 i = 1..10 i=1..10 i=1..10，即下面是一个10元函数。

y = ∑ i = 1 10 i cos ⁡ ( i x i 5 ) y = \sum_{i=1}^{10} i\cos(\frac{ix_i}{5}) y=i=1∑10​icos(5ixi​​)

测试代码写为

def test(xs):
    s = 0.0
    for i in range(len(xs)):
        s += np.cos(i*xs[i]/5)*(i+1)
    return s

if __name__ == "__main__":
    xRange = (-3,3)
    vRange = (-1,1)
    wRange = (0.5,1)
    C = (1.5, 1.5)
    s = Swarm(20, 10, xRange, vRange, wRange, C)
    s.optimize(test, 20, 200)
    print("最佳位置在：\n", s.xBest)

测试结果为

第1次迭代最小值为-29.665606552469008
第21次迭代最小值为-36.351923008184464
第41次迭代最小值为-47.14461157300818
第61次迭代最小值为-52.25291216470125
第81次迭代最小值为-52.964363058825406
第101次迭代最小值为-52.99124052250346
第121次迭代最小值为-52.99966161164592
第141次迭代最小值为-52.99990804268494
第161次迭代最小值为-52.99997494166753
第181次迭代最小值为-52.9999871143357
最佳位置在：
[ 28.01165182 -15.70858655   7.85375734   5.23563916  -3.92708994
  -3.14138732  -2.61760754   2.24427198   5.89048506  -1.74515961]

可见其收敛速度还是很快的。