PSO粒子群优化算法

算法介绍

粒子群优化（Particle Swarm Optimization, PSO），又称粒子群算法、微粒群算法，是由 J. Kennedy 和 R. C. Eberhart 等[1]于1995年开发的一种演化计算技术，来源于对一个简化社会模型的模拟。其中“群（swarm）”来源于微粒群符合 M. M. Millonas 在开发应用于人工生命（artificial life）的模型时所提出的群体智能的5个基本原则。“粒子（particle）”是一个折衷的选择，因为既需要将群体中的成员描述为没有质量、没有体积的，同时也需要描述它的速度和加速状态。 ——维基百科

我用一个小故事来简明扼要的给大家讲讲粒子群优化算法的基本原理和实现过程，仅代表个人理解，并不一定严谨，如有错误还请指出~

上课了，老师问小明，在一个鸟群中，鸟类是如何寻找到最优的觅食地点的呢？

小明说，是因为鸟儿会自己寻找最优的觅食地点！

老师回答道，不错，但是单个鸟类个体可以搜寻的空间是有限的，但是鸟类总是喜欢群体觅食，与个体想比，群体有什么优势呢？

小明恍然大悟，：我懂了！老师，一群鸟可以去各自想去的地方寻找食物，并互相交流告诉其它鸟自己找到的最好的地点在哪，这样鸟群就可以有很大的搜寻面积的同时找到整个鸟群搜寻过的面积中的最优解了！

其实这就是粒子群优化算法的基本思想了，下面我们来介绍一下这个搜索方法是如何被定义成为一个完整的算法的~

算法原理

首先我们知道，想要实现上面的觅食模拟，对于一个鸟群，我们需要知道的基本因素有以下三个鸟群中鸟类的数量，鸟群中个体找到的最优位置，鸟群整体的最优位置。

我们来依次定义每一个内容。

1. 粒子的数量（鸟的数量）

粒子的数量是比较好定义的，这里我们假设一共有m个粒子，这样我们就完成了算法的第一步，定义了粒子的数量为n。

2.粒子的最优位置（鸟群中个体搜寻到的最优位置）

如何定义一个粒子的位置呢？

我们知道，在地球表面，我们可以用经纬度来定位，这是因为地球表面可以假想为一个二维曲面，只需要两个坐标（经纬度）即可确定曲面上的唯一点。

学过高中物理的可能知道，在一维空间，也就是直线中，直线上一个点的位置就是这个点的坐标，用一个数就可以表示，

而在三维空间中，我们就需要三维坐标系，也就是三个坐标来加以表示了。

那么到D维的空间究竟需要多少个坐标来加以确定呢？

这里简单给出结论，容易推知，确定D维空间中的一个点，需要D个坐标。

因此，对于任意一个粒子$X_i$，他在D维空间的坐标应该由D个坐标表示，在算法中定义为：
$$
Xi =(x_{i1}, x_{i2}, …, x_{iD})
$$
这样我们就定义好了一个粒子所在的位置，我们也就可以记录一个粒子i经历过的最好的位置$P_i$（pbest）为：
$$
P_i = (p_{i1}, p_{i2}, …, p_{iD})
$$
而整体粒子搜寻到的最优位置，也可以用同样的方法记录$P_g$(gbest):
$$
P_g=best(P_1,P_2,…,P_m)
$$
定义好这些，我们就可以进行最优解的搜寻了吗？

我们遗憾的发现，似乎并不行，尽管我们定义好了鸟类个体的数量、个体的位置、最优位置和群体最优位置，这些鸟是如何按照自己各自的喜好去搜寻食物的？这依然是一个需要思考的问题，那么该如何定义鸟类的搜寻过程和动作呢？

我们很自然的就会想到，高中物理是如何描述一只鸟的运动的？

最简单的方法就是，你告诉我一只鸟的平均速度和运动的时间，我就可以得到这只鸟下一时刻的位置。

3.粒子的速度（鸟的速度）

为了得到每一只鸟的运动，以及下一时刻的位置，我们给出鸟每一时刻间平均速度的描述，学过高中物理运动的合成与分解的我们也可以很容易的知道，给出粒子每个时刻每个维度的分速度，我们就可以得到粒子的合速度，因此粒子的速度可以被描述为：
$$
V_i = (v_{i1},v_{i2},…,v_{iD})
$$
每一个时刻，粒子位置的改变于是就可以描述为：
$$
X_{id+1}=X_{id}+v_{id+1}\tag{1a}
$$
因为$X_{1d}$的位置即为初始位置，此时可以理解为粒子并没有运动，因此还没有平均速度，因此此刻的位置加上下一时刻的平均速度就可以得到粒子下一时刻的位置。

这样我们就得到了一个最基础的，每个粒子朝着固定方向搜寻并返回最优解的搜索算法，这样的算法显然是机械且搜索空间不足的，因此我们需要让粒子的运动更多样化，因此在每一次运动之间，许多PSO算法的变体给粒子引入了加速常数$c_1c_2$和惯性常数$w$:
$$
v_{id+1} = w∙v_{id}+c_1∙rand()∙(p_{id}-x_{id})+c_2∙Rand()∙(p_{gd}-x_{id})\tag{1b}
$$
$rand()$和$Rand()$是两个在[0,1]范围里变化的随机值，也就是说，鸟类首先会根据两个固定的加速常数，随机的朝着个体认为的最优解位置和群体认为的最优解位置进行加速调整，并在原本的惯性常数的定义下，朝着原本的方向运动一段距离，这里用高中的运动的合成和分解来理解会相当合适。

举例来说，单个粒子（鸟）的新速度会由三个速度合成，原本方向的惯性速度，朝个体最优解调整的速度，朝全体最优解方向调整的速度。

如图所示，黄色为原本的惯性速度，A、B分别为群体和个体最优解的位置，最后粒子的运动方向就是综合考量多个因素合成的速度方向，也就是最后的C‘。

到这，我们已经完成了鸟类进行搜寻空间的所有准备，可以开始正式的介绍算法了。

算法流程

标准的PSO算法流程如下：

1. 初始化一群微粒（群体规模为 m ），包括随机的位置和速度；
2. 评价每个微粒的适应度；
3. 对每个微粒，将它的适应值和它经历过的最好位置 pbest 的作比较，如果较好，则将其作为当前的最好位置pbest；
4. 对每个微粒，将它的适应值和全局所经历最好位置 gbest 的作比较，如果较好，则重新设置gbest的索引号；
5. 根据方程（1）变化微粒的速度和位置；
6. 如未达到结束条件（通常为足够好的适应值或达到一个预设最大代数 Gmax ），回到2.。

算法参数

PSO 参数包括：群体规模 m ，惯性权重 w ，加速常数 $c_1$ 和$ c_2$ ，最大速度 $V_{max}$，最大代数 $G_{max}$ 。

$V_{max}$决定在当前位置与最好位置之间的区域的分辨率（或精度）。如果 $V_{max}$太高，微粒可能会飞过好解，如果 $V_{max}$太小，微粒不能进行足够的探索，导致陷入局部优值。该限制有三个目的：防止计算溢出；实现人工学习和态度转变；决定问题空间搜索的粒度。

惯性权重w使微粒保持运动的惯性，使其有扩展搜索空间的趋势，有能力探索新的区域。

加速常数$c_1$ 和$ c_2$代表将每个微粒推向 pbest 和 gbest 位置的统计加速项的权重。低的值允许微粒在被拉回来之前可以在目标区域外徘徊，而高的值导致微粒突然的冲向或者越过目标区域。

如果没有后两部分，即 c1 = c2 = 0，微粒将一直以当前的速度飞行，直到到达边界。由于它只能搜索有限的区域，将很难找到好的解。

如果没有第一部分，即 w = 0，则速度只取决于微粒当前的位置和它们历史最好位置 pbest 和 gbest ，速度本身没有记忆性。假设一个微粒位于全局最好位置，它将保持静止。而其它微粒则飞向它本身最好位置 pbest 和全局最好位置 gbest 的加权中心。在这种条件下，微粒群将统计的收缩到当前的全局最好位置，更像一个局部算法。

在加上第一部分后，微粒有扩展搜索空间的趋势，即第一部分有全局搜索的能力。这也使得w的作用为针对不同的搜索问题，调整算法全局和局部搜索能力的平衡。

如果没有第二部分，即 c1 = 0，则微粒没有认知能力，也就是“只有社会（social-only）”的模型。在微粒的相互作用下，有能力到达新的搜索空间。它的收敛速度比标准版本更快，但是对复杂问题，比标准版本更容易陷入局部优值点。

如果没有第三部分，即 c2 = 0，则微粒之间没有社会信息共享，也就是“只有认知（cognition-only）”的模型。因为个体间没有交互，一个规模为m的群体等价于m个单个微粒的运行。因而得到解的几率非常小。

代码实现

这里以pyswarm库实现的PSO算法为例介绍PSO的算法代码是如何实现的。这个库里PSO代码的函数功能介绍还是比较详细的，有需要的朋友可以直接去看源码。

整体代码如下：

import numpy as np

def pso(func, lb, ub, ieqcons=[], f_ieqcons=None, args=(), kwargs={}, 
        swarmsize=100, omega=0.5, phip=0.5, phig=0.5, maxiter=100, 
        minstep=1e-8, minfunc=1e-8, debug=False):
    
    assert len(lb) == len(ub)
    assert hasattr(func, '__call__')
    lb = np.array(lb)
    ub = np.array(ub)
    assert np.all(ub > lb)
   
    vhigh = np.abs(ub - lb)
    vlow = -vhigh
    
    obj = lambda x: func(x, *args, **kwargs)
    if f_ieqcons is None:
        if not len(ieqcons):
            if debug:
                print('No constraints given.')
            cons = lambda x: np.array([0])
        else:
            if debug:
                print('Converting ieqcons to a single constraint function')
            cons = lambda x: np.array([y(x, *args, **kwargs) for y in ieqcons])
    else:
        if debug:
            print('Single constraint function given in f_ieqcons')
        cons = lambda x: np.array(f_ieqcons(x, *args, **kwargs))
        
    def is_feasible(x):
        check = np.all(cons(x) >= 0)
        return check
        
    S = swarmsize
    D = len(lb)
    x = np.random.rand(S, D)
    v = np.zeros_like(x)
    p = np.zeros_like(x)
    fp = np.zeros(S)
    g = []
    fg = 1e100
    
    for i in range(S):
        x[i, :] = lb + x[i, :] * (ub - lb)
        p[i, :] = x[i, :]
        fp[i] = obj(p[i, :])
        if i == 0:
            g = p[0, :].copy()
        if fp[i] < fg and is_feasible(p[i, :]):
            fg = fp[i]
            g = p[i, :].copy()
        v[i, :] = vlow + np.random.rand(D) * (vhigh - vlow)
       
    it = 1
    while it <= maxiter:
        rp = np.random.uniform(size=(S, D))
        rg = np.random.uniform(size=(S, D))
        for i in range(S):
            v[i, :] = omega * v[i, :] + phip * rp[i, :] * (p[i, :] - x[i, :]) + \
                      phig * rg[i, :] * (g - x[i, :])
            x[i, :] = x[i, :] + v[i, :]
            mark1 = x[i, :] < lb
            mark2 = x[i, :] > ub
            x[i, mark1] = lb[mark1]
            x[i, mark2] = ub[mark2]
            fx = obj(x[i, :])
            if fx < fp[i] and is_feasible(x[i, :]):
                p[i, :] = x[i, :].copy()
                fp[i] = fx
                if fx < fg:
                    if debug:
                        print('New best for swarm at iteration {:}: {:} {:}'.format(it, x[i, :], fx))
                    tmp = x[i, :].copy()
                    stepsize = np.sqrt(np.sum((g - tmp) ** 2))
                    if np.abs(fg - fx) <= minfunc:
                        print('Stopping search: Swarm best objective change less than {:}'.format(minfunc))
                        return tmp, fx
                    elif stepsize <= minstep:
                        print('Stopping search: Swarm best position change less than {:}'.format(minstep))
                        return tmp, fx
                    else:
                        g = tmp.copy()
                        fg = fx
        if debug:
            print('Best after iteration {:}: {:} {:}'.format(it, g, fg))
        it += 1

    print('Stopping search: maximum iterations reached --> {:}'.format(maxiter))
    
    if not is_feasible(g):
        print("However, the optimization couldn't find a feasible design. Sorry")
    return g, fg

算法有三个必须参数和12个可选\默认参数。

参数	类型	说明	默认值
func	function	要最小化的函数	-
lb	array	设计变量的下界	-
ub	array	设计变量的上界	-
ieqcons	list	函数列表，长度为n，满足ieqcons[j](x,*args) >= 0.0 表示优化成功	[]
f_ieqcons	function	返回一个1-D数组，其中每个元素在成功优化时必须大于或等于0.0。如果指定了f_ieqcons，则忽略ieqcons	None
args	tuple	传递给目标和约束函数的附加参数	()
kwargs	dict	传递给目标和约束函数的附加关键字参数	{}
swarmsize	int	群体中的粒子数量	100
omega	scalar	粒子速度缩放因子	0.5
phip	scalar	从粒子已知最佳位置搜索的缩放因子	0.5
phig	scalar	从群体已知最佳位置搜索的缩放因子	0.5
maxiter	int	群体搜索的最大迭代次数	100
minstep	scalar	群体最佳位置的最小步长，搜索终止	1e-8
minfunc	scalar	群体最佳目标值的最小变化，搜索终止	1e-8
debug	boolean	如果为True，每次迭代将显示进度语句	False

返回值

返回值	类型	说明
g	array	群体已知的最佳位置（最优设计）
f	scalar	在g处的函数目标值

算法使用

那么如何使用这个算法呢？

很简单，利用这个已经写好的算法库，我们定义出你需要进行优化的目标，然后输入进算法计算得到输出就可以啦~

我们来给出一个简单的示例：

假设我们的鸟群在一个二维空间搜寻食物，这个空间长宽都是20，空间中的食物和所在的位置有特定的函数关系，即不同位置会对应一个食物的数量，假设这个食物的数量和位置的关系是$F(x)=-x_0^2-x_1^2$我们可以定义为$x_0,x_1∈[-10,10]$,优化函数就是对应的这个关系函数，因为这个库是求函数的最小值，我们需要对我们的求食物最大的目标做一个小小的调整，也就是加上负号，具体的简单代码实现就如下所示：

import numpy as np
from pyswarm import pso

# 定义目标函数
def objective_function(x):
    return x[0]**2 + x[1]**2

# 定义变量的上下界
lower_bounds = [-10, -10]
upper_bounds = [10, 10]

# 调用PSO函数
best_position, best_value = pso(
    func=objective_function,
    lb=lower_bounds,
    ub=upper_bounds,
    swarmsize=30,  # 粒子数量
    omega=0.5,     # 惯性权重
    phip=0.5,      # 个体学习因子
    phig=0.5,      # 群体学习因子
    maxiter=100,   # 最大迭代次数
    debug=True     # 打开调试信息
)

print("最佳位置:", best_position)
print("最佳值:", best_value)

尝试运行几次，结果如下所示：

最佳位置 (x0, x1)	最佳值
[-1.2223 × 10⁻⁴, 8.967 × 10⁻⁶]	1.5021 × 10⁻⁸
[-7.6302 × 10⁻⁵, 1.0833 × 10⁻⁴]	1.7558 × 10⁻⁸
[-5.1334 × 10⁻⁴, -4.8416 × 10⁻⁴]	4.9793 × 10⁻⁷
[7.9078 × 10⁻⁵, 3.1941 × 10⁻⁵]	7.2735 × 10⁻⁹

我们知道，这个函数的最小值位置其实是[0,0],最小值其实是0，算法每次的结果都并不相同，但是都十分接近最优的结果。

模型不是直接得到最小值的原因是因为PSO优化算法并不是直接搜寻到最优解，PSO算法是一种元启发算法（英文：metaheuristic）, 又称 万能启发式算法、万用启发式算法。在计算机科学和数学优化中，元启发是一种高级的程序或启发式算法，专门用于搜索、生成或选取一个启发式结果（局部搜索算法），该结果可以为一个最优化问题提供足够好的求解，尤其适用于信息不完备或者计算能力受限时的最优化问题。

到这里，我们已经理解粒子群优化算法的基本思想以及该算法的简单使用了，接下来我要介绍一下我们研究生数学建模比赛（2024年C题）中使用到的PSO算法用到的小tips了，后面是拓展知识，但是实现起来也并不困难，仅供各位参考~

研究生数学建模大赛2024C题第五问

问题五 磁性元件的最优化条件

在磁性元件的设计与优化领域内，磁芯损耗固然是一个不容忽视的核心评价指标，但在工程实践中，为了实现磁性元件整体性能的卓越与最优化，需要综合考虑多个评价指标，其中，传输磁能就是重要的评价指标之一，因此，同时考虑磁芯损耗与传输磁能这二个评价指标，对于指导磁性元件的设计方向、优化其性能表现，具有重要的理论及实践意义。

请以问题四构建的磁芯损耗预测模型为目标函数，同时考虑传输磁能这个重要指标（由于传输磁能概念的复杂性，我们仅以频率与磁通密度峰值的乘积来衡量传输磁能大小），利用附件一中的实验数据，建立优化模型，分析在什么条件下（温度、频率、波形、磁通密度峰值及磁芯材料），能达到最小的磁芯损耗以及具有最大的传输磁能？

温度，取4个值：25、50、70、90，单位：摄氏度；

频率，取值范围：50000—500000，单位：赫兹；

磁芯损耗，单位：每立方米瓦特；

励磁波形类型：正弦波、三角波和梯形波；

峰值，从磁通密度采样点计算得到，共1024个采样点（一个周期时间内，相同间距采样），单位：特斯拉；

第四问中我们已经设计了一个机器学习模型，可以接收这五个值并计算得到对应的磁芯损耗

显然，我们的目标是最小化磁芯损耗，机器学习模型也被我们建模出来了，因此很自然的，我们将机器学习模型作为我们的目标函数，直接采用粒子群优化算法，问题不就迎刃而解了吗？

不过不知道你发现了没有，机器学习的模型输入有很多的分类指标，这些指标并不能输入为连续变量，如果你输入一个连续变量，这就意味着你的机器学习模型训练出来的预测值的准确度不一定可靠了，因为你的输入并不符合模型训练的过程。

因此，引入了一个新的问题——PSO优化算法如何引入分类变量进行优化空间搜寻？

下面是我们的解决办法：

PSO优化算法进行分类变量优化空间搜寻

针对这个问题，我定义了一个函数，用来将PSO优化算法的连续值优化为分类变量作为机器学习的输入以得到有效的输入和可靠的预测值。

首先我们定义好上下界，也就是粒子在空间中可以探索的范围：

这里温度是有序变量0、1、2、3代表逐渐升高的四种温度，频率和峰值是连续变量，波形和材料是独热编码的0/1变量。

参数	下界 (lb)	上界 (ub)	说明
温度 (°C)	-1	3	温度
频率 (Hz)	-3	3	频率
峰值	-3	3	峰值
励磁波形_正弦波	0	1	励磁波形_正弦波
励磁波形_三角波	0	1	励磁波形_三角波
励磁波形_梯形波	0	1	励磁波形_梯形波
磁芯材料_材料1	0	1	磁芯材料_材料1
磁芯材料_材料2	0	1	磁芯材料_材料2
磁芯材料_材料3	0	1	磁芯材料_材料3
磁芯材料_材料4	0	1	磁芯材料_材料4

在函数中，温度向上取整，将-1—3的连续变量转化为了0、1、2、3的分类变量，对独热编码，取四个材料中的最大值为1，其它为0，原则上，对PSO的四个取值应当有限制，但对于简单实现来说，这样的实现从输入和结果上并不会表现出明显的错误，因为时间有限，没有进一步进行限制。

def extract_and_normalize_parameters(best_solution, stats):
    """
    提取最佳解中的各个参数，并进行归一化处理

    参数：
    - best_solution: 最佳解数组，形状为(10,)
    - stats: 归一化参数字典

    返回：
    - 提取的参数和归一化后的频率和峰值
    """
    # 提取最佳解中的各个参数
   T = math.ceil(best_solution[0])  # 温度，向上取整
    F = best_solution[1]  # 频率
    waveforms = best_solution[2:5]
    max_waveform_index = np.argmax(waveforms)
    w_sine = 1 if max_waveform_index == 0 else 0
    w_triangle = 1 if max_waveform_index == 1 else 0
    w_trapezoid = 1 if max_waveform_index == 2 else 0
    Pk = best_solution[5]  # 峰值
    materials = best_solution[6:10]
    max_material_index = np.argmax(materials)
    β1 = 1 if max_material_index == 0 else 0
    β2 = 1 if max_material_index == 1 else 0
    β3 = 1 if max_material_index == 2 else 0
    β4 = 1 if max_material_index == 3 else 0

    # 获取频率和峰值的均值和标准差
    freq_mean = stats["频率，Hz"]["mean"]
    freq_std = stats["频率，Hz"]["std"]
    pk_mean = stats["峰值"]["mean"]
    pk_std = stats["峰值"]["std"]
    # 对频率F和峰值Pk进行z归一化
    f = (F - freq_mean) / freq_std
    pk = (Pk - pk_mean) / pk_std

    return T, F, w_sine, w_triangle, w_trapezoid, Pk, β1, β2, β3, β4, f, pk

经过函数处理后，我们就得到了可以正常输入进模型的不同变量，把这些变量输入进模型得到对应的预测值，并将预测值和频率和峰值的乘积一同放入优化函数中，我们就完成了一个可以正确处理机器学习函数以及分类变量的多目标优化的PSO算法实现~

我们的目标函数（适应度函数）的代码实现如下：

def fitness_function(x):
    
    

    # 检查x的长度是否足够
    if len(x) < 10:
        raise ValueError("x必须包含至少10个元素")

    # 处理温度，生成x最接近的整数作为T的输入
    T = math.ceil(x[0])  # x[0]对应温度，向上取整

    # 处理频率
    F = x[1]  # x[1]对应频率

    # 处理波形类型，x[2]到x[4]为波形类型
    waveforms = x[2:5]
    max_waveform_index = np.argmax(waveforms)
    w_sine = 1 if max_waveform_index == 0 else 0
    w_triangle = 1 if max_waveform_index == 1 else 0
    w_trapezoid = 1 if max_waveform_index == 2 else 0

    # 处理峰值
    Pk = x[5]  # x[5]对应峰值

    # 处理材料类型，x[6]到x[9]为材料类型
    materials = x[6:10]
    max_material_index = np.argmax(materials)
    β1 = 1 if max_material_index == 0 else 0
    β2 = 1 if max_material_index == 1 else 0
    β3 = 1 if max_material_index == 2 else 0
    β4 = 1 if max_material_index == 3 else 0
    

    # 获取频率和峰值的均值和标准差
    freq_mean = stats["频率，Hz"]["mean"]
    freq_std = stats["频率，Hz"]["std"]
    pk_mean = stats["峰值"]["mean"]
    pk_std = stats["峰值"]["std"]
    # 对频率F和峰值Pk进行z归一化
    f = (F - freq_mean) / np.sqrt(freq_std)
    pk = (Pk - pk_mean) / np.sqrt(pk_std)


    data = pd.DataFrame({
        '温度，oC': [T],
        '频率，Hz': [f],
        '励磁波形_正弦波': [w_sine],
        '励磁波形_三角波': [w_triangle],
        '励磁波形_梯形波': [w_trapezoid],
        '峰值': [pk],
        '磁芯材料_材料1': [β1],
        '磁芯材料_材料2': [β2],
        '磁芯材料_材料3': [β3],
        '磁芯材料_材料4': [β4]
    })

    select_features = ['温度，oC', '频率，Hz', '励磁波形_正弦波', '励磁波形_三角波', '励磁波形_梯形波', '峰值', '磁芯材料_材料1', '磁芯材料_材料2', '磁芯材料_材料3', '磁芯材料_材料4']

    # 模拟机器学习模型 Y 的输出
    Y = predict_with_xgboost(data, select_features)

    # 计算 pk * f
    product = Pk * F

    # 这里假设我们要最小化 Y 和最大化 pk * f
    # 适应度函数：最小化 Y - 权重 * 最大化 product
    return w_1*Y - w_2 * product  # 权重可调整

理解浅薄，希望可以为你提供一点思路~~~