粒子群算法（PSO）：从鸟群觅食到Rastrigin函数优化，手写一个完整的粒子群优化器粒子群算法（Particle S

摘要

粒子群算法（Particle Swarm Optimization）PSO ，参考鸟类觅食行为结合人类社会，核心假设是通过社会信息分享，个体可以从整个群体中获益，从而在寻找最优解中获取优势，常用于求解非线性、多峰函数的全局最优解，在机器学习超参数调优、工程优化、路径规划等领域有广泛应用。

引言

在计算机上存在许多优化问题，比如高数中的函数求极值问题，机器学习的超参数调优问题，还有无人机的路径规划问题。

计算爆炸：随着问题规模增大，精确算法（如穷举法、分支定界法）的计算时间会呈指数级增长。例如，为20个城市找最短路径，穷举所有可能性就需要检查约 2.4×10¹⁸ 条路径，远超任何计算机的处理能力。
现实约束苛刻：许多实际问题的目标函数（如成本、能耗）非线性、不连续、不可微，或者变量混杂了整数和实数。传统基于梯度或线性假设的数学规划方法根本无法处理。

假设你和同伴们在外面寻宝，分别在不同的位置进行寻宝，此时有一位同伴挖到宝了并通知了所有人，你就会倾向于往他寻到宝的方向进行寻宝，同时你又有自己的主见和经验，结果你会综合所有信息更新自己的寻宝方向。

在这个场景中，每个寻宝者就是一个“粒子”，它当前的位置对应解空间中的一个候选解，它移动的方向和快慢对应“速度”。粒子根据自身找到过的最好位置（个体最优）和整个队伍共享的最好位置（全局最优），动态调整自己的搜索方向——这一过程不断重复，直到找到“宝藏”（全局最优解）。

核心原理

PSO 算法起源对鸟群觅食、鱼群巡游等生物群体社会行为的模拟，其核心假设是：通过社会信息的共享，个体可以从整个群体的经验中获益，从而在寻找食物（最优解）时获得进化优势。

粒子表示

一个粒子 $i$ 在 $D$ 维空间中的位置： $x_i = (x_{i1},x_{i2}, \dots ,x_{iD})$
粒子速度： $v_i = (v_{i1},v_{i2}, \dots ,v_{iD})$
个体历史最优位置： $pbest$
全局历史最优位置： $gbest$
粒子适应度：位置越好，适应度越高

更新方程

粒子会结合个人经验（个体历史最优位置）和社会经验（全局历史最优位置）去更新自己的位置和速度从而不断靠近全局最优解。

速度更新：

1.首次提出：在 1995 年 PSO 算法首次提出时的速度更新公式：

\pmb{v_i^{t+1}} = \pmb{v_i^{t}} + 2*rand()*(pbest-\pmb{x_i}) + 2*rand()*(gbest-\pmb{x_i})

$pbest-x_i ：$ 粒子向自身最优位置 (pbest) 加速。
$gbest-x_i ：$ 粒子向全局最优位置 (gbest) 加速。
$rand() ：$ 随机因子区间[0,1]，增加搜索的随机性。
$2*rand() ：$ 让随机因子的期望值为1，使得粒子在每次迭代中平均有一半的概率“飞过”目标区域，从而增强探索能力。

补充：论文中实验移除原始速度发现算法在寻找全局最优时非常低效。

\pmb{v_i^{t+1}} = 2*rand()*(pbest-\pmb{x_i}) + 2*rand()*(gbest-\pmb{x_i})

2.惯性权重添加：会发现这个公式和现代使用的公式有一定的差别，原始公式中没有现在的 $w$ 、 $c$ 和 $r$ ，原始固定了 $w = 1，c = 2 ，r = rand()$ ，现代公式是在 1998 年《A Modified Particle Swarm Optimizer》提出：

\pmb{v_i^{t+1}} = w * \pmb{v_i^{t}} + c_1*r_1*(pbest-\pmb{x_i}) + c_2*r_2*(gbest-\pmb{x_i})

在 1995 年的原始版本中，粒子的速度更新非常剧烈，容易导致“飞过头”或在最优解附近震荡。

平衡搜索能力： $w$ 负责平衡全局探索和局部开发。
动态调整策略： 论文建议 $w$ 不应该是固定不变的。最好的做法是让 $w$ 随迭代次数 t 线性递减（例如从 0.9 降到 0.4）。
- 前期（高 $w$ ）： 粒子飞得快，主要进行大范围的粗略搜索。
- 后期（低 $w$ ）： 粒子飞得稳，主要在当前发现的最优区域进行细致寻找。

位置更新：

计算得到下一轮 (t+1) 的移动速度，粒子就会朝着新一轮的移动速度 $v_i^{t+1}$ 前进：

\pmb{x_i^{t+1}}=\pmb{x_i^t}+\pmb{v_i^{t+1}}

案例：粒子 A 受原本移动的方向 u ，个体历史最优位置 C ，全局历史最优位置 B 影响，最终通过公式计算速度方向是 v 。

算法流程

通过讲解 Pyhton 实现 PSO 进一步了解 PSO 算法。

定义适应度函数

PSO 的作用就是优化和求解全局最优解，所以需要有一个适应度函数作为度量标准。

class ObjectiveFunction:
    def __call__(self, x): 
        return 函数计算结果

初始化全局参数

粒子群数量 N ：简单函数粒子数量不需要太多，往往设置为 N = 20
维度 dim ：适应度函数自变量的数量，f(x) 或 f(x,y)
迭代次数 ger ：粒子群移动的次数
位置区间 limit : 适应度函数的定义域
速度区间 v_limit : 粒子单次移动的距离，要远小于位置区间limit否则容易飞出边界
惯性系数 w ：常用 w=0.9
学习因子 c1，c2 : 常用 c1,c2=2

N = 100  # 种群粒子数量
dim = 2  # 维度
ger = 100  # 迭代次数
limit = [-5, 5]  # 范围区间
v_limit = [-1, 1]  # 速度区间
w = 0.9  # 惯性系数
c1 = 2.0  # 自我学习因子
c2 = 2.0  # 群体学习因子

初始化粒子参数

粒子初始位置 x ：随机分布粒子的初始位置在位置区间内，形状为（N，dim）
粒子初始速度 v ：随机设置粒子的初始速度在速度区间内，形状为（N，dim）
个体历史最优位置 p_best : 刚开始没有个体经验，设置为初始位置，形状为（N，dim）
个体历史最优适应度 p_best_evaluate ：个体最优位置的适应度，形状为（N，1）
全局历史最优位置 g_best ：取 p_best_evaluate 最值的位置，形状为（1，dim）
全局历史最优适应度 g_best_evaluate ：全局最优位置的适应度，取 p_best_evaluate 最值，形状为（1，）

x = np.random.uniform(limit[0], limit[1], (N, dim))
v = np.random.uniform(v_limit[0], v_limit[1], (N, dim))
p_best = np.copy(x)
p_best_evaluate = f(x[:, 0], x[:, 1]) # 二维举例
g_best_idx = np.argmin(p_best_evaluate)
g_best = np.copy(p_best[g_best_idx])
g_best_evaluate = p_best_evaluate[g_best_idx]

循环迭代

当前轮次 t : 当前粒子执行的次数，初始化 t = 0
初始化随机变量 r1 r2 : [0,1]之间的随机数，增加随机性
计算新速度 v ：根据公式计算新的速度，保证速度在速度区间内
计算新位置 x ：根据公式计算新的位置，保证位置在位置区间内
计算新适应值 e ：将新位置代入适应度函数得到新适应度 e
更新个体最优 ：对比新适应度和个体历史最优适应度，条件更新其位置和适应度
更新全局最优 ：对比新适应度和全局历史最优适应度，条件更新其位置和适应度
轮次加 1 ：本轮结束 t++
动态更新 w （可选） ：动态更新惯性权重随轮次上升而下降（如0.9->0.4）

t = 0
while t < ger:
    # 更新惯性权重
    w = 0.9 - (0.5 * t / ger)
    # 计算新速度
    r1, r2 = np.random.rand(N, dim), np.random.rand(N, dim)
    v = w * v + c1 * r1 * (p_best - x) + c2 * r2 * (g_best - x)
    v = np.clip(v, v_limit[0], v_limit[1])
    # 计算新位置
    x = np.clip(x + v, limit[0], limit[1])
    # 计算新适应值
    evaluate_new = f(x)
    # 更新个体历史最优（使用NumPy的掩码Masking机制进行向量化操作）
    mask = evaluate_new < p_best_evaluate
    p_best[mask] = x[mask]
    p_best_evaluate[mask] = evaluate_new[mask]
    # 更新全局历史最优
    min_idx = np.argmin(evaluate_new)
    if np.min(p_best_evaluate) < g_best_evaluate:
        best_idx = np.argmin(p_best_evaluate)
        g_best = np.copy(p_best[best_idx])
        g_best_evaluate = p_best_evaluate[best_idx]
    # 轮次加1
    t += 1

手写 PSO

实现手写 PSO 算法解决函数求最值问题，这是 PSO 最常用的场景之一，优化求取函数的最值，如损失函数最小值。

实现 PSO 计算 Rastrigin 二维多峰函数最小值，位置（0，0），最小值 0 ，算法整体实现按照上面介绍的算法流程：

"""
    摘要：使用粒子群算法 PSO 寻找维度为 2 且区间有限的函数 f(x，y) 最小值。
    引言：常规找区间内函数最小值，通过求导或通过遍历区间得到。遍历区间的时间复杂度极高， PSO 参考鸟群觅食行为，核心假设为通过
        社会信息分享，个体可以从整个群体中获益，从而在寻找最优解中获取优势。
    核心原理：
        使用场景：优化目标达到最优值（最小值/最大值）；
        适应度：粒子位置计算的结果越小/大，适应度越高；
        个人位置 site ：粒子当前的位置（坐标）；
        个人历史最优位置 p_best ：粒子在移动过程中适应度最高的位置；
        全局历史最优位置 g_best : 所有粒子在移动过程中适应度最高的位置。
    数学推导：
        粒子通过 t 轮移动接近最优值
        粒子种群 x ，粒子位置 x_i(t) ，粒子区间x_max
        粒子速度 v_i(t) ，最大速度 v_max （v<=v_max）
        v_i(t+1) = 惯性系数 * v_i(t) + c_1 * r1 * (p_best - x_i(t)) + c_2 * r2 * (g_best - x_i(t))
        x_i(t+1) = x_i(t) + v_i(t+1)
"""

import numpy as np


# Rastrigin 函数（多峰）
class ObjectiveFunction:

    def __call__(self, x):
        n=x.shape[-1]
        return 10*n + np.sum(x**2-10*np.cos(2*np.pi*x),axis=1)


if __name__ == '__main__':
    # 1.定义函数
    f = ObjectiveFunction()

    # 2.初始化全局参数
    N = 100  # 种群粒子数量
    dim = 2  # 维度
    ger = 100  # 迭代次数
    limit = [-5, 5]  # 范围区间
    v_limit = [-1, 1]  # 速度区间
    w = 0.9  # 惯性系数
    c1 = 2.0  # 自我学习因子
    c2 = 2.0  # 群体学习因子

    # 3.初始化粒子参数
    x = np.random.uniform(limit[0], limit[1], (N, dim))
    v = np.random.uniform(v_limit[0], v_limit[1], (N, dim))
    p_best = np.copy(x)
    p_best_evaluate = f(x)
    g_best_idx = np.argmin(p_best_evaluate)
    g_best = np.copy(p_best[g_best_idx])
    g_best_evaluate = p_best_evaluate[g_best_idx]

    # 4.种群移动
    t = 0
    while t < ger:
        # 更新惯性权重
        w = 0.9 - (0.5 * t / ger)
        # 4.1 粒子移动
        r1, r2 = np.random.rand(N, dim), np.random.rand(N, dim)
        v = w * v + c1 * r1 * (p_best - x) + c2 * r2 * (g_best - x)
        # 计算新位置
        v = np.clip(v, v_limit[0], v_limit[1])
        x = np.clip(x + v, limit[0], limit[1])
        # 计算适应值
        evaluate_new = f(x)
        # 4.2更新粒子历史最优
        mask = evaluate_new < p_best_evaluate
        p_best[mask] = x[mask]
        p_best_evaluate[mask] = evaluate_new[mask]
        # 4.3更新全局历史最优
        min_idx = np.argmin(evaluate_new)
        if np.min(p_best_evaluate) < g_best_evaluate:
            best_idx = np.argmin(p_best_evaluate)
            g_best = np.copy(p_best[best_idx])
            g_best_evaluate = p_best_evaluate[best_idx]
        # 轮次加1
        t += 1

    print(f'粒子历史最优位置：{p_best}\n最优适应度：{p_best_evaluate}')
    print(f'种群历史最优位置：{g_best}\n最优适应度：{g_best_evaluate}')

绘图函数

初始化粒子随机位置展示，可以看到部分粒子随机的位置已经非常接近最小值点。

粒子会不断往最小值点靠近

最终绝大部分的粒子都会收缩到最小值点附近，图像历史最优适应度小于零是因为浮点数精度问题

绘图函数很长很复杂，我是使用AI辅助生成的，想要了解的可以到Study944-GitHub获取。

数学推导

Clerc & Kennedy (2002) 通过数学分析证明，如果粒子要在无限次的迭代中最终停下来（收敛），而不是漫无目的地在空间中无限飞行（发散），那么惯性权重 $w$ 、个体学习因子 $c_1$ 和社会学习因子 $c_2$ 之间必须满足特定的数值比例关系。

w>\frac{1}{2}(c_1+c_2)-1

并且隐含要求 $w < 1$ （以保证能量衰减）。

$\frac{1}{2}(c_1+c_2)-1$ ：粒子受到的“向心力”（被 $pbest$ 和 $gbest$ 拉的力量）对例子移动增加贡献。
$w$ ：代表了粒子的“惯性保持能力”，w越大，粒子越倾向于保持原有运动方向。
结论：如果惯性权重（保持速度的能力）相对于加速度系数设置得太小，粒子可能无法有效平衡拉力，导致轨迹出现不稳定的振荡甚至发散。

手写 PSO 时使用 $w=0.9$ ， $c_1 ,c_2=2$ ，遵循公式我的 $w>1$ 才合理，那为什么我上面的手写实现 PSO 没有遵循？

因为这个公式的核心是用在参数固定且搜索空间无限的场景下，用于限制粒子的移动速度避免发散到处飞，我的实现中参数 $w$ 是会收缩的，还设置了速度和位置区间限制，超过区间会截断。

pyswarms

pyswarms 是 Python 提供实现优化算法的第三方库，其中包含大部分优化算法的实现，如粒子群算法，遗传算法和蚁群算法，其中与本文手写 PSO 对应的是 pyswarms.single.GlobalBestPSO。

下载：pip install pyswarms

核心对象类 Swarm

Swarm 是 pyswarms 库中的核心数据结构，它封装了整个粒子群的所有状态信息。

position ：所有粒子的位置矩阵，形状（n,dim）
velocity ：所有粒子的速度矩阵，形状（n,dim）
n_particles : 粒子数量，整数
dimensions : 问题维度，整数
options ：配置参数字典，如 {'c1': 2.0, 'c2': 2.0, 'w': 0.9}
pbest_pos ：所有粒子个体历史最优位置，矩阵形状（n,dim）
pbest_cos ：所有粒子个体历史最优适应度，矩阵形状（n,dim）
best_pos ：全局历史最优位置，形状（dim,1）
best_cos ：全局历史最优适应度，浮点数
current_cos ：当前所有粒子的适应度，矩阵形状（n,）

@attrs
class Swarm(object):
    # Required attributes
    position = attrib(type=np.ndarray, validator=instance_of(np.ndarray))
    velocity = attrib(type=np.ndarray, validator=instance_of(np.ndarray))
    # With defaults
    n_particles = attrib(type=int, validator=instance_of(int))
    dimensions = attrib(type=int, validator=instance_of(int))
    options = attrib(type=dict, default={}, validator=instance_of(dict))
    pbest_pos = attrib(type=np.ndarray, validator=instance_of(np.ndarray))
    best_pos = attrib(
        type=np.ndarray,
        default=np.array([]),
        validator=instance_of(np.ndarray),
    )
    pbest_cost = attrib(
        type=np.ndarray,
        default=np.array([]),
        validator=instance_of(np.ndarray),
    )
    best_cost = attrib(
        type=float, default=np.inf, validator=instance_of((int, float))
    )
    current_cost = attrib(
        type=np.ndarray,
        default=np.array([]),
        validator=instance_of(np.ndarray),
    )

    @n_particles.default
    def n_particles_default(self):
        return self.position.shape[0]

    @dimensions.default
    def dimensions_default(self):
        return self.position.shape[1]

    @pbest_pos.default
    def pbest_pos_default(self):
        return self.position

主类 GlobalBestPSO

最外层直接使用的类是 GlobalBestPSO ，对其进行初始化后调用 optimize 函数传入参数即可获取粒子群中全局历史最优位置和历史最优适应度。

class GlobalBestPSO(SwarmOptimizer):
#...省略一大段代码
    def optimize(
        self, objective_func, iters, n_processes=None, verbose=True, **kwargs
    ):
        # ...省略部分代码
        # 初始化粒子群个体历史最优适应度为最大值
        self.swarm.pbest_cost = np.full(self.swarm_size[0], np.inf)
        ftol_history = deque(maxlen=self.ftol_iter)
        for i in self.rep.pbar(iters, self.name) if verbose else range(iters):
            # 计算更新粒子群当前的适应度
            self.swarm.current_cost = compute_objective_function(self.swarm, objective_func, pool=pool, **kwargs)
            # 计算更新粒子群个体历史最优位置和适应度
            self.swarm.pbest_pos, self.swarm.pbest_cost = compute_pbest(self.swarm)
            # 计算更新粒子群全局历史最优位置和适应度
            best_cost_yet_found = self.swarm.best_cost # 保存历史用于展示
            self.swarm.best_pos, self.swarm.best_cost = self.top.compute_gbest(self.swarm)
            # fmt: on
            if verbose:
                self.rep.hook(best_cost=self.swarm.best_cost)
            # 保存历史用于展示
            hist = self.ToHistory(
                best_cost=self.swarm.best_cost,
                mean_pbest_cost=np.mean(self.swarm.pbest_cost),
                mean_neighbor_cost=self.swarm.best_cost,
                position=self.swarm.position,
                velocity=self.swarm.velocity,
            )
            self._populate_history(hist)
            # 收敛判断，会提前判断是否已收敛，收敛则提前终止
            relative_measure = self.ftol * (1 + np.abs(best_cost_yet_found))
            delta = (
                np.abs(self.swarm.best_cost - best_cost_yet_found)
                < relative_measure
            )
            if i < self.ftol_iter:
                ftol_history.append(delta)
            else:
                ftol_history.append(delta)
                if all(ftol_history):
                    break
            # 更新配置参数，如惯性权重w随轮次iternow下降0.9->0.4
            self.swarm.options = self.oh(
                self.options, iternow=i, itermax=iters
            )
            # 更新速度
            self.swarm.velocity = self.top.compute_velocity(
                self.swarm, self.velocity_clamp, self.vh, self.bounds
            )
            # 更新位置
            self.swarm.position = self.top.compute_position(
                self.swarm, self.bounds, self.bh
            )
        # 获取全局历史最优位置和适应度
        final_best_cost = self.swarm.best_cost.copy()
        final_best_pos = self.swarm.pbest_pos[
            self.swarm.pbest_cost.argmin()
        ].copy()
        # 打日志
        self.rep.log(
            "Optimization finished | best cost: {}, best pos: {}".format(
                final_best_cost, final_best_pos
            ),
            lvl=log_level,
        )
        # 关闭线程池
        if n_processes is not None:
            pool.close()
        return (final_best_cost, final_best_pos)

核心实现 operators

核心逻辑实现在 GlobalBestPSO ，但只是通过逻辑串联实现了 PSO ，具体的实现在 operators 工具中，包括计算粒子适应度，计算更新个体最优位置和适应度和全局最优位置和适应度。

""" 计算更新个体历史最优位置和适应度 """
def compute_pbest(swarm):
    try:
        dimensions = swarm.dimensions
        # mask_cost 是一个布尔数组，形状为 (N,)，满足条件为True说明找到更好的，用于更新 pbest_cost
        mask_cost = swarm.current_cost < swarm.pbest_cost
        # mask_cost 是一个布尔矩阵，形状为 (N,dim)，用于更新 pbest_pos
        mask_pos = np.repeat(mask_cost[:, np.newaxis], dimensions, axis=1)
        # np.where(条件,结果1(True),结果2(False))
        new_pbest_pos = np.where(~mask_pos, swarm.pbest_pos, swarm.position)
        new_pbest_cost = np.where(
            ~mask_cost, swarm.pbest_cost, swarm.current_cost
        )
    # 异常处理
    except AttributeError:
        rep.logger.exception(
            "Please pass a Swarm class. You passed {}".format(type(swarm))
        )
        raise
    else:
        return (new_pbest_pos, new_pbest_cost)

""" 计算更新粒子速度 """
def compute_velocity(swarm, clamp, vh, bounds=None):
    try:
        swarm_size = swarm.position.shape
        c1 = swarm.options["c1"]
        c2 = swarm.options["c2"]
        w = swarm.options["w"]
        # 分别计算朝个体历史最优方向速度和全局历史最优方向速度
        cognitive = (
            c1
            * np.random.uniform(0, 1, swarm_size)
            * (swarm.pbest_pos - swarm.position)
        )
        social = (
            c2
            * np.random.uniform(0, 1, swarm_size)
            * (swarm.best_pos - swarm.position)
        )
        # 计算新速度
        temp_velocity = (w * swarm.velocity) + cognitive + social
        updated_velocity = vh(
            temp_velocity, clamp, position=swarm.position, bounds=bounds
        )
    # 异常处理
    except AttributeError:
        rep.logger.exception(
            "Please pass a Swarm class. You passed {}".format(type(swarm))
        )
        raise
    except KeyError:
        rep.logger.exception("Missing keyword in swarm.options")
        raise
    else:
        return updated_velocity

""" 计算更新粒子位置 """
def compute_position(swarm, bounds, bh):
    try:
        temp_position = swarm.position.copy()
        temp_position += swarm.velocity

        if bounds is not None:
            temp_position = bh(temp_position, bounds)

        position = temp_position
    except AttributeError:
        rep.logger.exception(
            "Please pass a Swarm class. You passed {}".format(type(swarm))
        )
        raise
    else:
        return position

""" 计算更新粒子适应度 """
def compute_objective_function(swarm, objective_func, pool=None, **kwargs):
    # 非线程池做法
    if pool is None:
        return objective_func(swarm.position, **kwargs)
    # 线程池做法
    else:
        results = pool.map(
            partial(objective_func, **kwargs),
            np.array_split(swarm.position, pool._processes),
        )
        return np.concatenate(results)

对比pyswarms和手写效果

对比在同等情况下 pyswarms 的效果和手写的效果，通过比较二者的历史最佳适应度。

补充：GlobalBestPSO 配置动态更新惯性权重 $w$ 的设置非常麻烦远不及自己手动更新方便。

"""
    对比手写 PSO 与 Python 自带 pyswarms
"""

from pyswarms.single import GlobalBestPSO
import numpy as np


# 多峰函数类
class ObjectiveFunction:

    def __call__(self, x):
        n = x.shape[1]
        return 10 * n + np.sum(x ** 2 - 10 * np.cos(2 * np.pi * x), axis=1)


# 多峰函数
def rastrigin_func(x):
    n = x.shape[1]
    return 10 * n + np.sum(x ** 2 - 10 * np.cos(2 * np.pi * x), axis=1)


if __name__ == '__main__':
    # 1.参数定义
    N = 100  # 种群粒子数量
    dim = 2  # 维度
    ger = 100  # 迭代次数
    limit = [-5, 5]  # 范围区间
    v_limit = [-1, 1]  # 速度区间
    w = 0.9  # 惯性系数
    c1 = 2.0  # 自我学习因子
    c2 = 2.0  # 群体学习因子

    # 2.GlobalBestPSO
    optimizer = GlobalBestPSO(
        n_particles=N,
        dimensions=dim,
        options={'c1': c1, 'c2': c2, 'w': w},
        bounds=(np.array([limit[0]] * dim), np.array([limit[1]] * dim)),
        velocity_clamp=(v_limit[0], v_limit[1]),
        vh_strategy="unmodified",
        bh_strategy="intermediate",
    )
    # cost, pos = optimizer.optimize(
    #     ,rastrigin_func
    #     iters=100,
    #     verbose=False
    # )
    # 3. 手动迭代循环
    # 由于GlobalBestPSO配置动态变化w非常麻烦，使用手动更新实现
    cost, pos = 0.0, np.zeros(dim)
    for i in range(ger):
        # 计算当前迭代的 w 
        current_w = 0.9 - (0.9 - 0.4) * (i / ger)
        # 更新惯性权重
        optimizer.options['w'] = current_w
        cost, pos = optimizer.optimize(
            rastrigin_func,
            iters=1,
            verbose=False
        )

    print(f'pyswarms最优适应度: {cost:.11e}')
    print(f'pyswarms最优位置: [{pos[0]:.11e}, {pos[1]:.11e}]')

    # 3.手写
    f = ObjectiveFunction()
    x = np.random.uniform(limit[0], limit[1], (N, dim))
    v = np.random.uniform(v_limit[0], v_limit[1], (N, dim))
    p_best = np.copy(x)
    p_best_evaluate = f(x)
    g_best_idx = np.argmin(p_best_evaluate)
    g_best = np.copy(p_best[g_best_idx])
    g_best_evaluate = p_best_evaluate[g_best_idx]

    # 种群移动
    t = 0
    while t < ger:
        # 更新惯性权重
        w = 0.9 - (0.5 * t / ger)
        # 粒子移动
        r1, r2 = np.random.rand(N, dim), np.random.rand(N, dim)
        v = w * v + c1 * r1 * (p_best - x) + c2 * r2 * (g_best - x)
        # 计算新位置
        v = np.clip(v, v_limit[0], v_limit[1])
        x = np.clip(x + v, limit[0], limit[1])
        # 计算适应值
        evaluate_new = f(x)
        # 更新粒子历史最优
        mask = evaluate_new < p_best_evaluate
        p_best[mask] = x[mask]
        p_best_evaluate[mask] = evaluate_new[mask]
        # 更新全局历史最优
        min_idx = np.argmin(evaluate_new)
        if np.min(p_best_evaluate) < g_best_evaluate:
            best_idx = np.argmin(p_best_evaluate)
            g_best = np.copy(p_best[best_idx])
            g_best_evaluate = p_best_evaluate[best_idx]
        # 轮次加1
        t += 1
    print(f'手写最优适应度: {g_best_evaluate}')
    print(f'手写最优位置: {g_best}')

    print('手写效果好' if g_best_evaluate < cost else 'pyswarms效果好')

参数分析实验

综合上面的介绍，可以得到 PSO 的核心参数有 惯性权重 $w$ 、个体学习因子 $c_1$ 和社会学习因子 $c_2$ ，还有种群规模 N 和迭代次数 ger。

惯性权重 $w$

固定 c1 = c2 = 2.0，对比ω = 0.9、ω = 0.7、ω = 0.4以及线性衰减ω = 0.9→0.4四条曲线，横轴迭代次数，纵轴gBest适应度。

"""
    固定 c1 = c2 = 2.0，对比ω = 0.9、ω = 0.7、ω = 0.4以及线性衰减ω = 0.9 → 0.4四条曲线，横轴迭代次数，纵轴gBest适应度。
"""
import numpy as np
import matplotlib.pyplot as plt


# Rastrigin 函数（多峰）
class ObjectiveFunction:

    def __call__(self, x):
        n = x.shape[1]
        return 10 * n + np.sum(x ** 2 - 10 * np.cos(2 * np.pi * x), axis=1)


def run_pso(w_strategy, N=50, dim=2, ger=100, c1=2.0, c2=2.0):
    # 1.定义函数
    f = ObjectiveFunction()

    # 2.定义参数
    N = N  # 种群粒子数量
    dim = dim  # 维度
    ger = ger  # 迭代次数
    limit = [-5, 5]  # 范围区间
    v_limit = [-1, 1]  # 速度区间
    # v_limit = [-10, 10]  # 速度区间
    c1 = c1  # 自我学习因子
    c2 = c2  # 群体学习因子

    # 3.初始化种群粒子
    np.random.seed(42)
    x = np.random.uniform(limit[0], limit[1], (N, dim))
    v = np.random.uniform(v_limit[0], v_limit[1], (N, dim))
    p_best = np.copy(x)
    p_best_evaluate = f(x)
    g_best_idx = np.argmin(p_best_evaluate)
    g_best = np.copy(p_best[g_best_idx])
    g_best_evaluate = p_best_evaluate[g_best_idx]
    g_best_evaluate_his = []

    # 4.种群移动
    t = 0
    while t < ger:
        # 更新惯性权重
        if w_strategy == 'liner':
            w = 0.9 - (0.5 * t / ger)
        else:
            w = w_strategy
        # 4.1 粒子移动
        r1, r2 = np.random.rand(N, dim), np.random.rand(N, dim)
        v = w * v + c1 * r1 * (p_best - x) + c2 * r2 * (g_best - x)
        # 计算新位置
        v = np.clip(v, v_limit[0], v_limit[1])
        x = np.clip(x + v, limit[0], limit[1])
        # 计算适应值
        evaluate_new = f(x)
        # 4.2更新粒子历史最优
        mask = evaluate_new < p_best_evaluate
        p_best[mask] = x[mask]
        p_best_evaluate[mask] = evaluate_new[mask]
        # 4.3更新全局历史最优
        if np.min(p_best_evaluate) < g_best_evaluate:
            best_idx = np.argmin(p_best_evaluate)
            g_best = np.copy(p_best[best_idx])
            g_best_evaluate = p_best_evaluate[best_idx]
        g_best_evaluate_his.append(g_best_evaluate)
        # 轮次加1
        t += 1

    print(f'粒子历史最优位置：{p_best}\n最优适应度：{p_best_evaluate}')
    print(f'种群历史最优位置：{g_best}\n最优适应度：{g_best_evaluate}')
    return g_best_evaluate_his


if __name__ == '__main__':
    # 分别计算ω = 0.9、ω = 0.7、ω = 0.4以及线性衰减ω = 0.9 → 0.4
    ger = 100 # [100,1000]
    dim = 2 # [2,100]
    his_04 = run_pso(w_strategy=0.4, ger=ger,dim=dim)
    his_07 = run_pso(w_strategy=0.7, ger=ger,dim=dim)
    his_09 = run_pso(w_strategy=0.9, ger=ger,dim=dim)
    his_liner = run_pso(w_strategy='liner', ger=ger,dim=dim)
    # 5.绘图
    plt.figure(figsize=(10, 6))
    plt.plot(range(ger), his_04, label='$\omega = 0.4$')
    plt.plot(range(ger), his_07, label='$\omega = 0.7$')
    plt.plot(range(ger), his_09, label='$\omega = 0.9$')
    plt.plot(range(ger), his_liner, label='$\omega = 0.9 -> 0.4$')
    plt.title(f'Comparison of Different Inertia Weights $\omega$ on Rastrigin Function dim = {dim}')
    plt.xlabel('Iterations')
    plt.ylabel('gBest Fitness Score')
    plt.legend()
    plt.yscale('log')
    plt.grid(True, which="both", ls="-", alpha=0.5)
    plt.show()

结果

进行了一组对照试验，分别在 dim = 2 / 100 , ger = 100 / 1000四种情况的实现，从理论上来说惯性权重随迭代次数线性递减的效果应该最优，因为其结合了前期的探索和后期的开发，理应获得最低的适应度（即搜索效果最好）。

dim = 2 , ger = 100

结果表明在该条件下 $w=0.4$ 效果最佳，但可以认为 $线性w$ 的缩减速度 0.9 -> 0.4 收缩过慢，导致后期的细致搜索时间不够，可以尝试增加迭代次数 ger 。

dim = 2 , ger = 1000

增加迭代次数 ger = 1000，结果还是没有很大变化， $w=0.4$ 效果最佳，但是可以清晰看到 $w=0.9$ 时收缩效果一般，是因为惯性系数过大导致粒子速度过大没办法精确移动到最小值点。

dim = 100 , ger = 100

添加怀疑 $w=0.4$ 效果好是因为维度问题，维度过低问题简单，即使搜索的速度低也能逐步接近最小值点，修改 dim = 100，复杂化适应度函数。

结果依旧表明 $w=0.4$ 效果最好，再次测试轮次因素。

dim = 100 , ger = 1000

结果表明在大部分条件下综合下来 $w=0.4$ 效果最佳，其次才是 $线性递减的w$ ，还是怀疑是不是迭代次数过多了的问题，下面再次测试。

dim = 200 , ger = 10

在同等的迭代次数下，不同惯性系数的效果相差不大，因为开始的搜索很大程度受限于运气成分，随机生成的粒子位置可能就在最小值点附近，所以参考价值不大。

最终结果表明 $w=0.4$ 效果最佳，其次才是 $线性递减的w$ ，合理提出结论：

解除速度限制：尝试不限制速度增长时 $w=0.4$ 效果是否还是最优，发现实验结果变化不大。

“快速收敛优于过度探索”：对于特定具有明显全局凹面特征的目标函数（如 Rastrigin），复杂的权重递减策略可能引入不必要的搜索延迟，而较小的固定惯性权重配合严格的速度限制，能提供更高效的寻优路径。

个体学习因子 c_1 和社会学习因子 c_2

设置 c1=1,c2=3 | c1=2,c2=2 | c1=3,c2=1三组参数进行对比，固定 w=0.4，低维（dim=2）时搜索空间简单，三种参数组合都能较快找到全局最优，因此曲线差异不大。高维（dim=200）时差异才被显著放大。

dim=2

在低维情况下，三组学习因子的适应度大小还是收缩情况都差不多。

dim=200

在高维情况下，个体学习因子 c1 的重要性就会提高，社会学习因子 c2 的重要性就会下降。

当 $c_1=3, c_2=1$ 时（绿线），收敛曲线不仅更低，而且下降更平滑，没有出现蓝线（ $c_1=1, c_2=3$ ）那种明显的长时间停滞。

搜索步长的差异：
- $c_1 > c_2$ ：粒子运动轨迹更偏向于在自身周围徘徊探路，运动具有更强的“自我纠偏”能力。
- $c_2 > c_1$ ：粒子表现出更强的“趋群性”，一旦 $gbest$ 停滞，整个种群就会陷入集体瘫痪。

前沿进展

1. 综合学习粒子群（CLPSO - Comprehensive Learning PSO）

解决痛点：标准 PSO 极易陷入局部最优（早熟收敛），因为每个粒子只向自己的 $pbest$ 和全局的 $gbest$ 学习。

核心机制：粒子不再只盯着全局最优，而是根据一个学习概率 $Pc$ ，从种群中所有粒子的 $pbest$ 中随机选择维度进行学习。
多样性来源：每个维度可能学习自不同的粒子。这种“杂交”式的学习方式极大增强了种群的探索能力（Exploration），防止因跟随某个错误的 $gbest$ 而全军覆没。
适用场景：维数极高、局部极值点极多的复杂多峰函数。

2. 多目标粒子群（MOPSO - Multi-Objective PSO）

解决痛点：现实中的问题往往有多个相互冲突的目标（例如：要求 KNN 准确率最高的同时，计算耗时最少）。

核心挑战：不存在单一的最优解，而是一个帕累托最优集（Pareto Front） 。
存档机制（External Archive） ：MOPSO 会维护一个外部存档，专门记录迭代过程中发现的所有非劣解（即没有被其他解全面打败的解）。
引导者选择：在更新速度时，粒子不再选择一个 $gbest$ ，而是通过“拥挤度距离”或“轮盘赌”从存档中选出一个作为引导者，以保证最终解在帕累托前沿分布均匀。

3. 并行 PSO（Parallel PSO）

解决痛点：当适应度函数计算极其耗时，单机串行速度无法忍受。

A. 基于 MapReduce 的并行（适用于超大规模数据）

Map 阶段：将种群划分为多个子群，分配到不同的计算节点执行位置更新和适应度计算。
Reduce 阶段：收集各节点的 $pbest$ ，对比产生全局 $gbest$ ，再分发回各节点。
特点：适合分布式系统，容错性强，但通信开销较大。

B. 基于 GPU 的并行（适用于计算密集型）

细粒度并行：利用 CUDA 将数以万计的粒子计算分配到 GPU 核心。每个线程处理一个粒子的位置/速度更新。
速度优势：在适应度函数本身涉及大量矩阵运算（如 KNN 的距离矩阵计算）时，GPU 方案通常能比 CPU 快 10-100 倍。

总结：

PSO以其“无梯度、群体协作”的简洁范式，在非线性、多峰函数优化中展现了强大的全局搜索能力。其核心思想——个体认知与社会学习的平衡——已在CLPSO、MOPSO等变体中不断演化。理解PSO的参数敏感性（ω、c1、c2）和收敛条件，是将其成功应用于机器学习超参数调优等实际问题的关键。

参考文献

Kennedy J, Eberhart R. Particle swarm optimization[C]. Proceedings of ICNN'95, 1995.
Shi Y, Eberhart R. A modified particle swarm optimizer[C]. IEEE World Congress on Computational Intelligence, 1998.
Clerc M, Kennedy J. The particle swarm - explosion, stability, and convergence in a multidimensional complex space[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(1): 58-73.
pyswarms documentation: pyswarms.readthedocs.io/

粒子群算法（PSO）：从鸟群觅食到Rastrigin函数优化，手写一个完整的粒子群优化器

摘要

引言