粒子群优化算法Python实现

粒子群优化算法原理参考：
如何直观形象地理解粒子群算法？
最优化算法之粒子群算法（PSO）
参数设置
关于粒子群算法的两种惯性权重

粒子群优化实现思路：
1、知道粒子群优化原理后，首先需要设置初始参数：粒子个数n、社会学习因子c1、个体学习因子c2、收敛条件loss、惯性因子w
2、确定需要优化函数的维度，给出每个维度下参数的取值范围
例如：

对于该函数有三个维度x,y,z，在进行粒子群优化时需要给出每个维度的取值区间x=[-10, 10], y=[-80, 20], z=[-9, 2]
3、在步骤2中给出取值范围后，先随机初始化每个粒子的初始位置Position（每个粒子初始化的位置都在步骤2给的范围内），将每个粒子的初始速度V0设为0，将每个粒子的个体最优解Pbest设为与初始化位置Position相等，Pbest=Position。
4、计算每一个粒子当前位置下目标函数f的值Value，记录最小值所对应的粒子index
5、更新个体最优解Pbest。将4中算出的当前时刻的Value，与上一时刻计算的Value记为lastValue进行比较，看哪一个粒子的lastValue大于Value，那么就将这个粒子的当前时刻坐标设为这个粒子的Pbest，否则不变。即：（如果当上一步值比当前值大，那么就将Pbest的坐标值改变，也就是每一个粒子Pbest的值等于该粒子当前时刻历史数据中所处位置对应目标函数最小值的位置）
6、更新全局最优解Pz。全局最优解为4中index对应的粒子的位置
7、根据下面两个公式更新粒子的速度与位置

8、重复4-8直到loss收敛为止。

粒子群优化Python实现：

import numpy as np
import random


class ParticleSwarmOptimization:
    """粒子群优化算法"""
    def __init__(self, n=20, c1=2, c2=2, loss=0.001, w=0.6):
        """
        :param n: 粒子个数 一般为[10, 40]个粒子。根据实际情况而定可以取更多
        :param c1: 社会学习因子 通常在[0, 4]区间范围，一般取2
        :param c2: 个体学习因子 通常在[0, 4]区间范围，一般取2
        :param loss: 判断是否收敛条件
        :param w:惯性因子 较大时全局搜索能力强，较小时局部搜索能力强（本文为定值） 参考：https://blog.csdn.net/qq_38883271/article/details/101538704
        """
        self.n = n
        self.c1 = c1
        self.c2 = c2
        self.loss = loss
        self.w = w

    def fit(self, region):
        """
        初始化变量
        :param region: 每个维度的可行域，该变量的格式为字典
        例如：region = {0:[20, 50], 1:[30, 100],....}
            其中0为第一个维度，其可行域为[20， 50]
            其中1为第二个维度，其可行域为[30, 100]
        """
        self.partical = {}  # 存储每个粒子初始化的位置结果
        # 为每个粒子进行初始化位置
        for num in range(self.n):  # 初始化self.n个粒子
            init = []  # 用于存储每个粒子的初始化参数
            for ndim in range(len(region)):  # 每个粒子初始化len(region)个维度
                init.append(random.uniform(region[ndim][0], region[ndim][1]))
            self.partical[num] = init

        # 初始化每个粒子的初始速度V为0。
        self.V = {}
        for temp in range(self.n):
            self.V[temp] = 0

        # 设每个粒子的个体最优解Pb为当前初始化的位置（也可以设为随机位置）
        self.Pb = self.partical.copy()

    def predict(self, callback):
        """
        计算最终最优结果
        :param callback: 回调函数，用于优化的目标函数
        :return:
        """
        # 存储上一步的位置信息
        last_position = self.partical.copy()
        last_value = []  # 存储上一步目标函数值，value的上一步值
        # 存储位置变化绝对值
        loss = 100
        # 循环更新每个粒子的位置，直到位置变化小于某一值
        loop = 0  # 循环次数
        while True:
            # 计算每一个粒子当前位置下目标函数的值,并将最小值所对应的粒子记录下来index
            min = callback(self.partical[0])  # 设最小值为min
            index = 0  # 最小值所对应的粒子索引
            value = []  # 存储每个粒子所处位置对应目标函数的值
            for temp in range(self.n):  # 对每一个粒子进行遍历
                value.append(callback(self.partical[temp]))  # 计算每一个粒子当前时刻位置所对应目标函数的值，存储到value中
                if min > value[temp]:
                    min = value[temp]
                    index = temp
            if loop == 0:
                last_value = value.copy()
            else:
                # 更新个体最优解Pb
                for temp in range(self.n):
                    if last_value[temp] > value[temp]:  # 如果当上一步值比当前值大，那么就将self.Pb的坐标值改变，也就是每一个粒子self.Pb的值等于该粒子当前时刻历史数据中所处位置对应目标函数最小值的位置
                        last_value[temp] = value[temp]
                        self.Pb[temp] = self.partical[temp]
            # 将全局最优解设为index粒子对应的坐标
            Pz = self.partical[index]
            # 更新每个粒子的坐标
            for temp in range(self.n):
                # 更新速度
                if temp != index:
                    self.V[temp] = self.w*np.array(self.V[temp]) + self.c1*random.random()*(np.array(Pz) - np.array(self.partical[temp]))\
                                   + self.c2*random.random()*(np.array(self.Pb[temp]) - np.array(self.partical[temp]))  # 这个地方谁减谁不能搞反， 不然出错
                    # 更新位置
                    self.partical[temp] = np.array(self.partical[temp]) + np.array(self.V[temp])
                    # 计算更新前后位置变化
                    if loss > np.sqrt(np.mean(np.square(np.array(self.partical[temp]) - np.array(last_position[temp])))):
                        loss = np.sqrt(np.mean(np.square(np.array(self.partical[temp]) - np.array(last_position[temp]))))

            # 判断位置变化，是否跳出循环，是否收敛
            print(f"loop:{loop}  loss:{loss}")
            if loss < self.loss:
                break
            last_position = self.partical.copy()
            loop += 1

        return Pz

粒子群优化算法实例分析：

对于上图目标函数进行参数寻优：

def callback(x):
    """回调函数"""
    sum = 0
    for i in range(len(x)):
        sum += (x[i]-i)*(x[i]-i)
    sum /= 100

    return sum


def run():
    # 实例化粒子群优化算法
    pso = ParticleSwarmOptimization(n=10, c1=2, c2=2, loss=0.001, w=0.6)
    # 定义变量取值范围
    region = {0:[-50, 50], 1:[-50, 50], 2:[-40, 100]}
    pso.fit(region)
    result = pso.predict(callback)
    print(result)


if __name__ == "__main__":
    run()

运行结果：

[2.00528254e-03 9.99545977e-01 2.00752164e+00]

由上面运行结果可以看出：实际结果应该为：[0, 1, 2],运行结果为[0.002005, 0.9995, 2.0075],结果相差不大