【python】利用python实现简单粒子群优化算法实例

前言

最近在准备复（bai）习（lan）智能信息处理ing……想着通过复现一些PPT上简单算法例子的方式加深自己对算法的理解，也作为大家使用粒子群这一算法的思路与代码参考hhh

算法简述

粒子群算法（PSO）是1995年由美国社会心理学家J.Kennedy和电气工程师R.Eberhart共同提出，其基本思想是受到鸟类群体行为研究结果的启发，由于计算的快速性和易于实现从而引起重视。

PSO算法是通过随机化一群粒子的位置和运动速度，通过群体之间的联系和不同微粒之间个体的速率差异，通过一定的迭代规则，使其不断趋近于优化目标的真实解或解集，即求取最（极）大值或是最（极）小值的过程。

在这里，我们可以将每一个微粒个体看作是在n维搜索空间中的一个没有质量和体积的微粒，并在搜索空间中以一定的速度飞行，而飞行的速度由个体的飞行经验和群体的飞行经验进行动态调整。

我们以求函数最小值的优化问题帮助理解粒子群优化算法并尝试通过代码实现。

实例

今有如下Rosenbrock函数，求解其最小值：

$$f(x,y)=100(y-x^2{)}^2+(1-x{)}^2$$

其中，限定定义域
$$x,y\in [-30,30]$$

粒子最大速度$V_{max}=60$

函数大致图示如下：

由于问题的维数为2（x和y），所以我们每个粒子的速度和位置都应该是2维的向量，例如：

对于第i个粒子，有：

$$locatio{n}_i=(x_i,y_i)$$
$$v=(v_{xi},v_{yi})$$

紧接着对位置和速度进行随机数初始化，限定范围为[-30,30]即可。

接下来就是计算每个粒子的适应度（fitness）。

适应度这一概念应用比较广泛，在遗传算法等一些其它优化算法中，这一名词也经常出现。在这里它表征的是每个粒子对于优化问题结果的评价指标，比如，以本例来说，我们希望找到函数的最小值以及取得最小值的点，那么我们就可以把本次实例的粒子适应度确定为粒子当前位置的函数值，函数值越小，意味着我们可能越接近最小值，即离我们的优化目标越近。

综上，我们可以通过直接求每个粒子的当前函数值来作为粒子当前的适应度。

在求取适应度值后，选取对于当前问题的局部最优解（个体历史最优解）和全局最优解（群体历史最优解）。

其中局部最优解指的是每一轮迭代后，所有粒子当前状态的适应度最大（小）值；全局最优解指的是所有粒子从开始迭代优化一直到现在的这一大段时间里出现过的适应度最大（小）值；当算法完成迭代轮次后，应该输出的是最终状态的全局最优解以及取得最优解的点。

接下来是更新粒子的速度和位置。

我们采用如下速度和位置的更新函数：

$$v_{k+1}=w\ast v_k+c_1\ast \xi (p_k-x_k)+c_2\ast \eta (p_g-x_k)$$
$$x_{k+1}=x_k+v_{k+1}$$

其中$w$为惯性权重（调节对解空间搜索范围），本例取1，$c_1$和$c_2$系数为加速度常数，本例均取2，$\xi$和$\eta$为（0，1）之间的随机数。

通过上述更新式完成对粒子位置和速度的更新（要注意一下，本例限定速度在任何一个维度不能超过60，但对于位置，一般不强行拉回解空间）。

重复上述步骤，直至达到迭代次数或迭代结果与真实值之间差值小于某一给定误差。

代码部分

这里直接给出全部过程代码，关键位置给出注释，欢迎大家有问题随时来讨论

# -*- coding:utf-8 -*-
# 2022/11/29
# Jonathan.K.Wolf
import numpy as np

"""
每一个粒子的形式：
z:
[x, y]
"""
# 种群规模
m = 5
# 迭代参数
c0 = 1
c1 = 2
c2 = 2


# 定义函数
def Rosenbrock(z):
    sum = 100.0 * np.power(z[1] - np.power(z[0], 2), 2) + np.power((1 - z[0]), 2)
    return sum


# 初始化种群
def init_group(size):
    group = []
    for i in range(size):
        group.append([])
    for i in range(size):
        for j in range(2):
            group[i].append(np.random.uniform(-30.0, 30.0))
    return group


# 初始化速度
def init_speed(size):
    speed = []
    for i in range(size):
        speed.append([])
    for i in range(size):
        for j in range(2):
            speed[i].append(np.random.uniform(-60.0, 60.0))
    return speed


# 计算适应度
def calculate_fitness(group):
    fitness = []
    for item in range(len(group)):
        fitness.append(Rosenbrock(group[item]))
    return fitness, fitness.index(min(fitness))


# 更新位置和速度
def update_loc_speed(group, speed, pk, pg):
    new_group = []
    new_speed = []
    for item in range(len(group)):
        new_group.append([])
        new_speed.append([])
    for item in range(len(group)):
        for index in range(len(group[0])):
            new_speed[item].append(
                speed[item][index] + 2 * np.random.rand() * (pk[index] - group[item][index]) + 2 * np.random.rand() * (
                        pg[index] - group[item][index]))
            new_group[item].append(group[item][index] + new_speed[item][index])
    # 限定速度范围
    for item in range(len(speed)):
        for index in range(len(speed[0])):
            if new_speed[item][index] > 60.0:
                new_speed[item][index] = 60.0
            if new_speed[item][index] < -60.0:
                new_speed[item][index] = -60.0
    return new_group, new_speed


# 主要过程迭代
def process():
    epoch = 100000
    group_real = init_group(5)
    speed_real = init_speed(5)
    # 历史最优解
    history_best = group_real[0]
    # 当前最优解
    p = []

    # print(group_real)
    # print(speed_real)
    for i in range(epoch):
        fitness_real, pg_ = calculate_fitness(group_real)
        p.append(group_real[pg_])
        if Rosenbrock(p[i]) <= Rosenbrock(history_best):
            history_best = p[i]
        group_real, speed_real = update_loc_speed(group_real, speed_real, p[i], history_best)
        if i % 100 == 0:
            # print("current epoch:%d, current solution:%e" % (i, Rosenbrock(history_best)))
            print("current epoch:{}, current solution:{}".format(i, history_best))
    print(
        "the answer after {} epoch is ({},{}), and current value is {}".format(epoch, history_best[0], history_best[1],
                                                                               Rosenbrock(history_best)))


if __name__ == '__main__':
    process()

运行结果

100000轮次下得到最小值点及结果如图，与真实点（1，1）和最小值0相差无几，成功印证了算法的可用性

当然本篇只是实现最基本的粒子群优化算法，针对于优化问题的更优算法和改进方式还有很多，篇幅有限，这里不再详细解读了hhh

最后也欢迎各位大佬有问题随时批评指正~