粒子群算法

粒子群算法（Particle Swarm Optimization, PSO）是一种模拟鸟群觅食行为的群体智能优化算法，1995 年由 Kennedy 和 Eberhart 提出。它实现简单、参数少、收敛快、不依赖梯度信息，是数学建模竞赛中求解连续变量全局优化与模型参数寻优问题的"性价比之王"——当你面对一个没有显式导数、又可能有多个局部极值的黑箱目标函数时，PSO 往往是最值得先尝试的算法之一。本文从原理、适用场景、评价指标、可视化诊断到可运行代码（手写 PSO + scipy 库对照），完整梳理 PSO 的竞赛实战用法。

一、算法含义

1.1 通俗理解：鸟群觅食的隐喻

设想一群鸟在一片陌生的田野上找食物。田野里食物分布未知（可能有好几处，但只有一处最丰盛）。每只鸟都不知道食物在哪，但都知道两件事：

1. 自己的经验：自己飞过的所有位置中，哪一处离食物最近（个体历史最优，pbest）；
2. 群体的经验：整个鸟群到目前为止发现的最好位置在哪（群体历史最优，gbest）。

于是每只鸟下一步往哪飞，由三股力量合成：惯性（继续沿原来方向飞）＋自我认知（向自己历史最好的位置靠拢）＋社会认知（向全群最好的位置靠拢）。鸟群就这样在"各自探索 + 相互通报"中，逐渐聚集到食物最丰盛的地方。这就是粒子群算法的全部思想。

对应到优化问题：

• 每只鸟 = 一个粒子（particle），代表一个候选解；
• 鸟的位置 = 决策变量向量 $x_i$（即候选解本身）；
• 鸟的速度 = 速度向量 $v_i$（下一步移动的方向与步长）；
• 位置处食物的丰盛程度 = 适应度 $f(x_i)$（目标函数值，最小化问题中越小越好）；
• 每只鸟记忆中的"自己最好的位置" = 个体最优 $p_{best,i}$；
• 全群已知的"最好的位置" = 全局最优 $g_{best}$。

粒子没有质量、没有体积，可以瞬间转向；算法不需要目标函数的导数，只需要能计算函数值，因此特别适合"黑箱"目标函数。

1.2 核心公式：速度-位置更新

PSO 每代对每个粒子做两件事：更新速度，再更新位置：

$v_i^{t+1} = w\, v_i^t + c_1 r_1 \left( p_{best,i}^t - x_i^t \right) + c_2 r_2 \left( g_{best}^t - x_i^t \right)$

$x_i^{t+1} = x_i^t + v_i^{t+1}$

逐项解释速度更新公式（右端三项，对应鸟群的三股力量）：

项	名称	含义	作用
$w v_i^t$	惯性项	保留上一代速度的 $w$ 倍	维持原有飞行方向；$w$ 越大粒子越"固执"，全局搜索（探索）能力强；$w$ 越小越容易转向，局部精细搜索（开发）能力强
$c_1 r_1 (p_{best,i}^t - x_i^t)$	个体认知项（自我学习）	指向"自己历史最佳位置"的拉力	把粒子拉回自己去过的最好的地方，保持个体多样性
$c_2 r_2 (g_{best}^t - x_i^t)$	社会认知项（群体学习）	指向"全群历史最佳位置"的拉力	把粒子拉向群体最好的地方，使全群向最优解聚集

其中：

• $w$ 是惯性权重，控制上一代速度的影响；
• $c_1, c_2$ 是学习因子（加速系数），控制个体经验与群体经验的影响强度；
• $r_1, r_2$ 是 $[0,1]$ 上均匀分布的随机数，每个粒子每一维每一代都重新抽取——随机性让粒子走不同的路线，避免全体"排队走一条道"。

几个有意思的特殊情形（帮助理解各项作用）：

• 若 $c_1 = c_2 = 0$：粒子只按惯性匀速飞行，直到撞上搜索域边界，算法退化为"随机漫步"；
• 若 $c_1 = 0$（没有个体经验，只有社会学习）：全体粒子只向 gbest 冲，收敛很快但极易早熟陷入局部最优——被称为"只有社会模型"；
• 若 $c_2 = 0$（只有个体经验，没有社会学习）：各粒子互不交流、各找各的，等价于 $N$ 个独立的随机爬山，收敛很慢；
• 若 $w = 0$：粒子完全"没有记忆"，速度只由当前位置与两个最优位置的差决定。

1.3 pbest 与 gbest：算法中的"记忆"

PSO 与遗传算法最本质的区别在于记忆：每个粒子记住自己到过的最好位置。

• 个体最优 pbest（personal best）：粒子 $i$ 从第 0 代到第 $t$ 代所经历的最优位置。更新规则（最小化问题）：

$p_{best,i}^{t+1} = \begin{cases} p_{best,i}^{t}, & f(x_i^{t+1}) \ge f(p_{best,i}^{t}) \\[4pt] x_i^{t+1}, & f(x_i^{t+1}) < f(p_{best,i}^{t}) \end{cases}$

• 全局最优 gbest（global best）：全体粒子的 pbest 中最优的那个：

$g_{best}^{t+1} = \arg\min_{p \in \{p_{best,1}^{t+1}, \dots, p_{best,N}^{t+1}\}} f(p)$

拓扑结构还分两种：全局版 PSO（星型拓扑，每个粒子都能看到全群的 gbest，收敛快、易早熟）与局部版 PSO（环形拓扑，粒子只与左右邻居交流，收敛慢、不易早熟）。竞赛中默认使用全局版即可，遇到早熟问题再考虑局部版。

1.4 算法流程（伪代码）

输入：目标函数 f(x)（最小化），粒子数 N，维数 D，搜索域 [LB, UB]^D，
      惯性权重起止值 w_start、w_end，学习因子 c1、c2，最大迭代次数 T，速度上限 vmax
1:  随机初始化：x_i ~ U(LB, UB)，v_i ~ U(-vmax, vmax)，i = 1, ..., N
2:  计算适应度 f(x_i)；置 pbest_i = x_i；gbest = argmin_i f(pbest_i)
3:  for t = 0 to T-1:
4:      w(t) = w_start - (w_start - w_end) * t / (T-1)      # 惯性权重线性递减
5:      for i = 1 to N:
6:          r1, r2 ~ U(0, 1)                                 # 每个维度独立抽样
7:          v_i = w(t)*v_i + c1*r1*(pbest_i - x_i) + c2*r2*(gbest - x_i)
8:          v_i = clip(v_i, -vmax, vmax)                     # 速度限幅
9:          x_i = clip(x_i + v_i, LB, UB)                    # 位置更新并夹回搜索域
10:         if f(x_i) < f(pbest_i): pbest_i = x_i            # 更新个体最优
11:     gbest = argmin_i f(pbest_i)                          # 更新全局最优
输出：全局最优位置 gbest 与最优适应度 f(gbest)

复杂度：每次迭代评估 $N$ 次函数值，总评估次数 $N \times T$。算法存储量仅 $O(ND)$，实现起来不超过 100 行代码。

1.5 参数设置建议

参数	推荐取值	说明
粒子数 $N$	低维（$D \le 5$）：2050；中高维：50100	经验公式 $N \approx 10D$ 可作起点；太少多样性不足易早熟，太多浪费计算
惯性权重 $w$	0.9 → 0.4 线性递减	Shi & Eberhart (1998) 的经典方案：前期大 $w$ 全局探索，后期小 $w$ 精细开发；线性递减公式 $w(t) = w_{start} - (w_{start} - w_{end}) \cdot \frac{t}{T-1}$
学习因子 $c_1, c_2$	$c_1 = c_2 = 2$	Kennedy & Eberhart 的原始建议，适用面最广；也可取 $c_1 = c_2 = 1.49445$ 配合收缩因子（见第八节）
速度上限 $v_{max}$	搜索域宽度的 10%~20%	本例 $v_{max} = 0.2 \times (5.12 - (-5.12)) = 2.048$；太大会"飞过"最优解来回震荡，太小则搜索缓慢
迭代次数 $T$	竞赛常用 100~1000	可配合停止准则：连续若干代 gbest 变化小于阈值 $\varepsilon$（如 $10^{-8}$）即提前停止
随机种子	固定 seed	保证结果可复现；论文中必须报告多次独立运行的统计量（见 3.4）

1.6 优缺点

优点：

• 实现极简单：核心就两条更新公式，几十行代码即可完成，出错率低；
• 不依赖梯度：目标函数不可导、不连续、甚至是黑箱（只能求值）都能用；
• 参数少：核心参数只有 $N, w, c_1, c_2$，调参负担小；
• 收敛快：前期通常几十代就能进入最优解附近，比遗传算法、模拟退火通常更快；
• 天然并行：粒子之间相互独立，容易向量化（numpy 矩阵运算）或并行加速；
• 连续优化效果好：在参数标定、拟合、路径规划等连续问题上久经考验。

缺点：

• 易早熟：群体可能过早聚集到某个局部最优（尤其是多峰函数），且一旦聚集很难再散开；
• 随机性强：单次运行结果有波动，必须多次运行取均值/标准差才有说服力；
• 理论收敛性弱：不像遗传算法有较完善的收敛性分析，PSO 的全局收敛保证较弱（原始版本甚至不能保证收敛）；
• 不直接适用于离散问题：位置/速度天然是连续量，离散问题需改造为二进制 PSO 等变体；
• 精度有限：收敛到最优解"附近"容易，收敛到高精度（如 $10^{-10}$ 级）慢，需要与局部搜索混合。

二、何时使用（适用场景与条件）

2.1 适用场景

1. 连续变量全局优化：决策变量是实数，目标函数多峰（有多个局部极值），传统梯度法容易陷入局部最优时；
2. 模型参数寻优/标定：模型里有几个无法解析求解的参数（如 SEIR 传染病模型的传播率、Logistic 增长模型的参数、微分方程参数反演），以拟合误差最小为目标搜索参数；
3. 黑箱/无梯度目标：目标函数不可导（含绝对值、max、分段）、仿真软件给出的输出（"仿真优化"）、或无法写出解析式时；
4. 与机器学习结合：如 SVM 的 $C, \gamma$ 超参数寻优、随机森林参数、神经网络结构参数（此时 PSO 与网格搜索、贝叶斯优化是同类工具）；
5. 其他连续子问题：多模型融合的权重优化、资源连续分配、无人设备路径的连续参数化（如轨迹控制点坐标）。

2.2 竞赛典型题目

• 参数拟合类：给出观测数据，模型形式已知但参数未知（如"建立病毒传播模型并估计传播参数"），把参数当作粒子、把拟合误差当作适应度；
• 调度与分配类：加工顺序确定后对加工时间/投放量的连续优化部分；
• 路径规划类：把轨迹离散成控制点，用 PSO 优化控制点坐标，适应度为路径长度 + 避障惩罚；
• PID 控制器参数整定：控制类赛题中优化 $K_p, K_i, K_d$；
• 灰色模型/预测模型的参数修正：如 GM(1,1) 背景值系数、ARIMA 参数等；
• 作为基准算法之一出现在对比实验中（与 GA、SA 对比，体现工作量与严谨性）。

2.3 使用前提

• 目标函数可求值（每次给一组决策变量能算出适应度），不要求可导；
• 决策变量连续且能给出有意义的搜索域上下界（$[LB, UB]$ 的选取直接影响算法效果，界太宽浪费搜索）；
• 目标是最小化或最大化（最大化问题加个负号即可转换）；
• 能接受随机近似解：PSO 给出的是"很好的解"，不是"证明最优的解"，出题要求严格最优时不适合。

2.4 不适用情形

情形	原因与对策
纯离散组合问题（如 0-1 背包、TSP 的访问顺序）	位置/速度是连续量，需改用二进制 PSO（BPSO）或改用 GA/蚁群
高精度要求（如误差必须 < $10^{-10}$）	PSO 后期收敛慢，应"PSO 粗搜 + 局部搜索精化"混合
变量维数极高（成千上万维）	粒子群在超高维空间效率骤降，需降维或换算法
有解析解或梯度可用的光滑凸问题	直接用梯度法/牛顿法更准更快，PSO 是"杀鸡用牛刀"
严格全局最优证明类问题	PSO 无法证明最优性，需用分支定界、动态规划等精确算法

2.5 与遗传算法 / 模拟退火 / 蚁群算法的对比选择

方法	核心机制	记忆机制	参数数量	收敛速度	全局搜索能力	主要适用
PSO 粒子群	速度-位置更新，向 pbest/gbest 飞行	有（pbest/gbest）	少（$w, c_1, c_2, N$）	快	中（易早熟）	连续全局优化、参数寻优
GA 遗传算法	选择、交叉、变异	无（种群整体进化）	中（交叉率、变异率、编码等）	中	强（多样性机制完备）	离散/组合优化、编码灵活的问题
SA 模拟退火	以概率接受劣解、温度缓慢下降	有（当前最优解）	少（初温、降温系数）	慢	强（理论上可全局收敛）	小规模问题、可接受较长计算时间
ACO 蚁群算法	信息素积累与挥发、路径概率选择	有（信息素矩阵）	中（$\alpha, \beta, \rho$ 等）	中	中	离散路径问题（TSP、VRP、路径规划）

竞赛选择口诀：变量连续、想快出结果 → PSO；变量离散、问题有编码结构 → GA；规模小、要理论保证 → SA；本质是找路径 → ACO。也可以两两组合成混合算法（如 GA 粗搜 + PSO 精搜），增加论文亮点。

三、算法指标

PSO 是随机算法，评价它必须"既看最好、又看平均、还要看稳定"。下面每个指标给出：中文名、公式、方向、如何解读。符号定义见第五节。

3.1 历史最优适应度（全局最优适应度）

$f_{gbest}(t) = \min_{1 \le i \le N} f\left(p_{best,i}(t)\right)$

• 方向：最小化问题中越小越好（最大化问题取相反方向）；
• 解读：第 $t$ 代结束时算法"目前找到的最好解"的质量。它是每代都记录的序列，最终报告值取 $f_{gbest}^* = f_{gbest}(T)$。收敛曲线图（第四节图①）画的就是这条序列。注意 $f_{gbest}(t)$ 关于 $t$ 单调不增（gbest 只会被更优的解替换），所以曲线必然单调下降——这是检查实现是否有 bug 的一个快速办法。

3.2 每代平均适应度

$\bar{f}(t) = \frac{1}{N}\sum_{i=1}^{N} f\left(x_i^t\right)$

• 方向：越小越好；
• 解读：全体粒子当前适应度的平均水平，是种群"整体质量"的体温计。它与 $f_{gbest}(t)$ 的差距反映种群多样性：差距大说明粒子散布广、还在探索；两者几乎重合说明粒子已聚集到一处（对多峰问题要警惕"全体挤在局部最优"）。它不单调（粒子可能飞到更差的位置），与 gbest 曲线的单调性形成对照。

3.3 收敛代数

$T_{conv} = \min\left\{ t \in \{0, 1, \dots, T\} : f_{gbest}(t) \le f_{gbest}^* + \varepsilon \right\}$

• 方向：越小越好（本文取 $\varepsilon = 10^{-8}$）；
• 解读：gbest 首次达到"最终精度"所需的迭代代数。它回答"算法到底要跑多久"：如果 $T_{conv}$ 远小于 $T$，说明预算的迭代次数绰绰有余，可以适当减少 $T$ 或粒子数；如果 $T_{conv} = T$（直到最后一轮才不再进步），说明还没收敛透，应加大 $T$ 或调参。论文中常写"算法在约 $T_{conv}$ 代内收敛"。

3.4 多次运行最优值的均值与标准差（随机算法必须报告）

$\bar{f}^* = \frac{1}{M}\sum_{k=1}^{M} f_k^*, \qquad s_{f^*} = \sqrt{\frac{1}{M-1}\sum_{k=1}^{M}\left(f_k^* - \bar{f}^*\right)^2}$

其中 $f_k^*$ 是第 $k$ 次独立运行得到的最终最优适应度，$M$ 为独立运行次数（本文 $M = 10$）。

• 方向：均值越小越好；标准差越小越好；
• 解读：均值度量算法的平均求解质量，标准差度量稳定性。只报一次运行结果对随机算法是没有说服力的——同一份代码换个随机种子结果可能完全不同。竞赛论文的标准做法：独立运行 10 次，报告"最优值均值 $\bar{f}^*$、标准差 $s_{f^*}$（以及最好值/最差值）"。标准差相对均值很小（如 $s_{f^*}/\bar{f}^* < 1\%$）说明算法稳定可复现。

3.5 与已知全局最优的差距（误差）

$gap = f_{gbest}^* - f^*$

• 方向：越小越好（最小化问题中 gap $\ge 0$）；
• 解读：当测试函数（如 Rastrigin 函数的全局最优 $f^* = 0$）或模拟数据的真值已知时，gap 是绝对精度的直接度量。在仿真验证中常用它证明"算法找到了真解"；在真实赛题中真值未知，gap 退化为与"已知最好已知解"的比较（如与文献结果、与 scipy 库结果对比）。

3.6 惯性权重 w 的敏感性

$w$ 是 PSO 中最影响"探索 vs 开发"平衡的参数，通常通过控制变量实验考察：固定 $N, c_1, c_2, T, v_{max}$ 与随机种子不变，仅改变 $w$ 的取值方式（如 $w = 0.4$ 固定、$w = 0.7$ 固定、$w = 0.9$ 固定、$w: 0.9 \to 0.4$ 线性递减），比较各自的最终最优值与收敛曲线。若不同 $w$ 下结果相差不大，说明算法对该参数不敏感（稳健），论文中可一笔带过；若差异明显，则应在论文中说明取值的依据并给出对比图（第四节图④）。本例中 $w$ 的影响见 7.1 节的实验结果。

3.7 其他可选指标

• 成功率：$P_{succ} = \dfrac{1}{M}\sum_{k=1}^{M} \mathbb{1}\left[ f_k^* - f^* \le \varepsilon_{acc} \right]$，即 $M$ 次运行中达到精度 $\varepsilon_{acc}$（如 $10^{-6}$）的比例，度量算法可靠性；
• 函数评估次数：$N_{fev} = N \times T$（PSO 每代对全体粒子评估一次），度量计算成本，与收敛代数配合可计算"达到收敛实际花了多少次评估"；
• 运行时间：秒级单次耗时，用于算法间横向比较。

3.8 指标汇总表

指标	中文名	公式	方向/范围	一句话解读
$f_{gbest}(t)$	历史最优适应度	$\min_i f(p_{best,i}(t))$	越小越好	算法当前找到的最好解，收敛曲线即它的轨迹，随 $t$ 单调不增
$\bar{f}(t)$	每代平均适应度	$\frac{1}{N}\sum_i f(x_i^t)$	越小越好	种群整体质量与多样性；与 gbest 的差距反映探索程度
$T_{conv}$	收敛代数	$\min\{t : f_{gbest}(t) \le f_{gbest}^* + \varepsilon\}$	越小越好	首次达到最终精度所需代数，衡量"要跑多久"
$\bar{f}^*$	多次运行最优值均值	$\frac{1}{M}\sum_k f_k^*$	越小越好	平均求解质量，随机算法必报
$s_{f^*}$	多次运行最优值标准差	$\sqrt{\frac{1}{M-1}\sum_k (f_k^* - \bar{f}^*)^2}$	越小越好	稳定性/可复现性，随机算法必报
gap	与已知全局最优的差距	$f_{gbest}^* - f^*$	越小越好，$\ge 0$	绝对精度，仅真值已知时可用
$P_{succ}$	成功率	$\frac{1}{M}\sum_k \mathbb{1}[f_k^* - f^* \le \varepsilon_{acc}]$	$[0,1]$ 越大越好	达到给定精度的运行比例
$N_{fev}$	函数评估次数	$N \times T$	越小越省时	计算成本

四、可视化图表

PSO 的随机与迭代特性决定了"光报一个最优值"远远不够，画图看过程才是诊断算法健康度的关键。以下 4 张图覆盖"收敛稳定性、粒子运动、种群分布演化、参数敏感性"，全部代码见第六节，图片自动保存到 figures/ 目录。

4.1 四张图速查表

图名（输出文件）	用途	关键解读点
① 最优适应度收敛曲线（pso_convergence.png）	展示 10 次独立运行 gbest 随迭代的变化，检验收敛性与稳定性	曲线应呈阶梯式下降（快速下降段与局部极小平台段交替）；10 条线（对数坐标下）最终汇聚到同一水平 = 算法稳定（本例全部收敛到全局最优 0）；若某条线停在明显更高的平台 = 该次运行陷入局部最优；均值曲线（黑粗线）代表平均水平
② 粒子运动轨迹图（pso_trajectory.png）	在 Rastrigin 等高线图上画出若干代表粒子从初始位置到最终位置的运动路径	轨迹初期大幅跳跃（全局搜索），后期在原点附近小步密集（局部精细搜索）；轨迹会穿越多个"盆地"（局部极小）最终落入中心盆地；圆点为初始位置、叉号为最终位置、五角星为全局最优点，叉号应聚集在五角星附近
③ 迭代三阶段粒子分布图（pso_distribution.png）	对比迭代初期/中期/后期全体粒子的空间分布，展示"聚集"过程	初期：粒子均匀散布全搜索域（距原点距离中位数约 4.6）；中期：粒子明显向中心收缩，最优粒子已落入局部极小盆地（gbest 在 $f \approx 1$ 处停滞）；后期：50 个粒子中 48 个进入原点附近的 $10^{-4}$ 级小盆地，仅个别滞留局部极小——这就是"鸟群聚拢"的可视化
④ 惯性权重 w 对比图（pso_w_sensitivity.png）	展示不同 $w$ 设置下收敛曲线的差异，说明参数影响	$w=0.9$ 固定：曲线下降最慢且最终停在高位（本例 $6.3\times10^{-3}$），惯性太大"刹不住车"、无法精细收敛；$w=0.4$、$w=0.7$ 固定与线性递减：本例均收敛到 0；结合文献，固定小 $w$ 在多峰问题上早熟风险更高，故推荐 0.9→0.4 线性递减

4.2 每张图"好"与"异常"的特征

图① 收敛曲线（多次运行）

• 好图特征：对数坐标下曲线呈阶梯式下降——快速下降段（跨越多个盆地）与平台段（在某个局部极小处短暂停滞）交替出现，最终降至全局最优附近；多条曲线最终高度几乎一致；均值曲线平滑下降。
• 异常特征一（陷入局部最优）：个别曲线停在明显高于其他曲线的平台（如 Rastrigin 的 $f \approx 1$ 级），且后期纹丝不动；
• 异常特征二（未收敛）：所有曲线到最后一轮仍在缓慢下降，说明 $T$ 不够；
• 异常特征三（震荡）：曲线上下抖动（gbest 理论单调，出现抖动说明实现有 bug——gbest 更新条件写反或把平均适应度画成了 gbest）。

图② 粒子运动轨迹图

• 好图特征：轨迹呈"先长途奔袭、后原地打转"的形态——初期每一步跨度大（速度大、惯性大），后期轨迹在最优解附近缠绕成小团；
• 好图特征：轨迹有绕开"盆地"壁面的趋势，最终端（叉号）落在同一深色（低适应度）盆地内；
• 异常特征：轨迹像布朗运动一样无方向、始终在全域乱窜 → 说明粒子没有"向最优飞行"的行为（$c_1, c_2$ 太小或 $w$ 恒为 0）；轨迹全部直线冲边界 → 速度爆炸或限幅失效。

图③ 三阶段分布图

• 好图特征：三面板从左到右呈现"均匀散布 → 多簇并立 → 单簇聚拢"的清晰演化，最后一簇的中心与五角星（全局最优）重合；
• 异常特征一（早熟）：后期粒子聚成单簇，但簇中心停在某个局部极小盆地而非全局最优盆地；
• 异常特征二（未聚集）：后期粒子仍散布全域，说明还没收敛（$T$ 不够、$w$ 太大或 $c$ 太小）。

图④ w 敏感性对比图

• 好图特征：能看出 $w$ 对最终精度的实际影响——本例中 $w=0.9$ 固定的曲线下降最慢且最终停在高位（$6.3\times10^{-3}$），说明大惯性下粒子"刹不住车"、难以在最优解附近精细搜索；$w=0.4$ 固定、$w=0.7$ 固定与线性递减三条曲线均收敛到 0 附近；
• 异常特征：四条曲线几乎完全重合 → 说明问题对 $w$ 不敏感（这也是有效信息，论文里写"参数稳健"）；某条曲线发散或停在极高值 → 该参数设置失效（如 $w$ 过大导致震荡）。

五、符号说明

符号	含义	示例/单位
$N$	粒子数（种群规模）	$N = 50$
$D$（或 $d$）	决策变量维数	2 维、5 维
$x_i^t$	第 $i$ 个粒子在第 $t$ 代的位置向量（候选解）	$(3.2, -1.1)$
$v_i^t$	第 $i$ 个粒子在第 $t$ 代的速度向量	$(0.8, -0.3)$
$w$	惯性权重	$0.9 \to 0.4$ 线性递减
$c_1$	个体认知（学习）因子	$c_1 = 2$
$c_2$	社会认知（学习）因子	$c_2 = 2$
$r_1, r_2$	$[0,1]$ 均匀随机数，每维每代重新抽样	$r \sim U(0,1)$
$p_{best,i}$	粒子 $i$ 的个体历史最优位置	更新规则见 1.3
$g_{best}$	全体粒子的历史最优位置（全局最优）	本例收敛于 $(0, 0)$
$f(\cdot)$	目标函数（适应度函数），最小化	Rastrigin 函数
$f_{gbest}(t)$	第 $t$ 代的历史最优适应度	见 3.1
$\bar{f}(t)$	第 $t$ 代种群平均适应度	见 3.2
$v_{max}$	速度上限（限幅值）	$0.2 \times (UB - LB) = 2.048$
$[LB, UB]$	搜索域上下界	$[-5.12, 5.12]$
$T$	最大迭代次数	$T = 200$
$t$	当前迭代代数（从 1 起编号）	$t = 1, 2, \dots, T$
$T_{conv}$	收敛代数	见 3.3
$M$	独立运行次数	$M = 10$
$f^*$	已知全局最优值	Rastrigin：$f^* = 0$
$f_k^*$	第 $k$ 次运行的最终最优适应度	见 3.4
$\bar{f}^*$	$M$ 次运行最优值的均值	见 3.4
$s_{f^*}$	$M$ 次运行最优值的样本标准差	见 3.4
$\varepsilon$	收敛判定阈值	$10^{-8}$
gap	与已知全局最优的差距	gap $= f_{gbest}^* - f^*$
$\chi$	收缩因子（带收缩因子的 PSO）	$\chi \approx 0.72984$（见第八节）

六、可运行程序（完整代码）

环境要求：Python 3.12，依赖 numpy、scipy、matplotlib（均已在项目 .venv 中安装，无需 pyswarm 等额外库）。以下所有代码块按顺序拼接保存为 pso_demo.py 后运行即可：控制台打印第三节的全部指标，并在 figures/ 子目录生成 4 张诊断图（pso_convergence.png、pso_trajectory.png、pso_distribution.png、pso_w_sensitivity.png）。无图形界面环境可设 MPLBACKEND=Agg 跳过 plt.show() 的窗口弹出，图片照常保存。

# -*- coding: utf-8 -*-
"""
================================================================
粒子群算法（PSO）完整示例：手写实现 + scipy 库对照
----------------------------------------------------------------
主实例：2 维 Rastrigin 函数
    f(x) = 20 + Σ (x_i² − 10·cos(2πx_i))，x_i ∈ [−5.12, 5.12]
    已知全局最优：f = 0，位于原点 x* = (0, 0)
附加实例：5 维 Rastrigin 函数（展示维数升高后求解难度上升）
库实现对照：scipy.optimize.differential_evolution（作为成熟进化算法的参照）
输出：控制台打印全部指标（历史最优、平均适应度、收敛代数、多次运行
      均值/标准差、与全局最优的差距、w 敏感性），并生成 4 张图
依赖：numpy、scipy、matplotlib（Python 3.12，无其他第三方库）
================================================================
"""

# ========== 0. 导入库与全局设置 ==========
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution   # 库对照：差分进化（全局优化）

np.random.seed(42)   # 固定全局随机种子；PSO 内部另行传入独立种子，保证完全可复现

# ---- matplotlib 中文显示设置（防止图内中文乱码）----
plt.rcParams["font.sans-serif"] = ["PingFang SC", "Arial Unicode MS", "SimHei"]
plt.rcParams["axes.unicode_minus"] = False   # 让负号"-"正常显示

# ---- 图片输出目录（相对当前工作目录的 figures/ 子目录）----
FIG_DIR = "figures"
os.makedirs(FIG_DIR, exist_ok=True)

# ========== 1. 目标函数：Rastrigin 函数（多峰、含大量局部极小） ==========
def rastrigin(x):
    """Rastrigin 函数：f(x) = 10d + Σ(x_i² − 10·cos(2πx_i))，d 维。
    输入 x 形状 (d,) 或 (n, d)，沿最后一维求和，返回形状 () 或 (n,)。
    全局最优 f = 0 位于原点，是检验全局优化算法"抗早熟"能力的经典函数。"""
    x = np.asarray(x, dtype=float)
    d = x.shape[-1]
    A = 10.0
    return A * d + np.sum(x ** 2 - A * np.cos(2.0 * np.pi * x), axis=-1)

LB, UB = -5.12, 5.12            # 搜索域上下界（每个维度相同）
TRUE_OPT = 0.0                  # 已知全局最优值 f* = 0
TRUE_OPT_X = np.zeros(2)        # 已知全局最优点（2 维时为原点）
print(f"目标函数：Rastrigin，搜索域 [{LB}, {UB}]²，已知全局最优 f* = {TRUE_OPT}")

# ========== 2. 手写 PSO 实现（完整类） ==========
class ParticleSwarm:
    """标准粒子群算法（惯性权重线性递减 + 速度限幅）

    参数
    ----
    n_particles : 粒子数（种群规模），默认 50
    dim         : 决策变量维数，默认 2
    bounds      : 每个维度的搜索区间 (下界, 上界)，默认 (-5.12, 5.12)
    max_iter    : 最大迭代次数，默认 200
    w_start, w_end : 惯性权重 w 的起止值，随迭代线性递减（0.9 -> 0.4）
    c1, c2      : 个体认知与社会认知学习因子，默认均为 2
    vmax        : 速度上限；默认取搜索域宽度的 20%
    seed        : 随机种子（保证结果可复现）
    """

    def __init__(self, n_particles=50, dim=2, bounds=(-5.12, 5.12),
                 max_iter=200, w_start=0.9, w_end=0.4,
                 c1=2.0, c2=2.0, vmax=None, seed=None):
        self.n_particles = n_particles
        self.dim = dim
        self.lb = np.full(dim, bounds[0])
        self.ub = np.full(dim, bounds[1])
        self.max_iter = max_iter
        self.w_start, self.w_end = w_start, w_end
        self.c1, self.c2 = c1, c2
        # 速度上限：默认取搜索域宽度的 20%（本例 = 0.2 * 10.24 = 2.048）
        self.vmax = vmax if vmax is not None else 0.2 * (self.ub - self.lb)
        self.rng = np.random.RandomState(seed)

    def _init_swarm(self):
        """随机初始化：位置在搜索域内均匀分布，速度在 [-vmax, vmax] 内均匀分布"""
        X = self.rng.uniform(self.lb, self.ub, (self.n_particles, self.dim))
        V = self.rng.uniform(-self.vmax, self.vmax, (self.n_particles, self.dim))
        return X, V

    def optimize(self, func, track_idx=(0, 7, 19, 42), verbose=False):
        """运行 PSO，返回结果字典（含历史信息，供画图与指标计算使用）

        track_idx : 需要记录完整运动轨迹的代表粒子编号（用于画轨迹图）
        """
        X, V = self._init_swarm()

        # ---- 初始化：评估适应度，初始化 pbest 与 gbest ----
        fit = np.asarray(func(X), dtype=float)          # 当前适应度 f(x)
        pbest = X.copy()                                 # 个体历史最优位置
        pbest_fit = fit.copy()                           # 个体历史最优适应度
        g = int(np.argmin(pbest_fit))                    # 全局最优粒子下标
        gbest = pbest[g].copy()                          # 全局最优位置
        gbest_fit = pbest_fit[g]                         # 全局最优适应度

        # ---- 历史记录（用于指标统计与画图）----
        gbest_hist = [gbest_fit]                         # 每代全局最优适应度
        mean_hist = [float(np.mean(fit))]                # 每代种群平均适应度
        gbest_pos_hist = [gbest.copy()]                  # 每代全局最优位置
        X_hist = [X.copy()]                              # 每代全部粒子位置
        track = {i: [X[i].copy()] for i in track_idx}    # 代表粒子的位置轨迹

        # ---- 主迭代循环 ----
        for t in range(self.max_iter):
            # (1) 惯性权重线性递减：w_start -> w_end（Shi & Eberhart, 1998）
            w = self.w_start - (self.w_start - self.w_end) * t / max(1, self.max_iter - 1)

            # (2) 速度更新：v = w*v + c1*r1*(pbest - x) + c2*r2*(gbest - x)
            r1 = self.rng.rand(self.n_particles, self.dim)
            r2 = self.rng.rand(self.n_particles, self.dim)
            V = (w * V
                 + self.c1 * r1 * (pbest - X)
                 + self.c2 * r2 * (gbest - X))

            # (3) 速度限幅：防止速度爆炸，保证算法数值稳定
            V = np.clip(V, -self.vmax, self.vmax)

            # (4) 位置更新：x = x + v，越界粒子夹回搜索域边界
            X = np.clip(X + V, self.lb, self.ub)

            # (5) 评估新位置，更新 pbest 与 gbest
            fit = np.asarray(func(X), dtype=float)
            better = fit < pbest_fit                  # 比自身历史最优更优的粒子
            pbest[better] = X[better]
            pbest_fit[better] = fit[better]
            g = int(np.argmin(pbest_fit))
            if pbest_fit[g] < gbest_fit:              # 全局最优被刷新
                gbest = pbest[g].copy()
                gbest_fit = pbest_fit[g]

            # (6) 记录历史（每代都存，供画图使用）
            gbest_hist.append(gbest_fit)
            mean_hist.append(float(np.mean(fit)))
            gbest_pos_hist.append(gbest.copy())
            X_hist.append(X.copy())
            for i in track_idx:
                track[i].append(X[i].copy())

            # (7) 每 20 代打印一次进度（verbose=True 时）
            if verbose and ((t + 1) == 1 or (t + 1) % 20 == 0 or t == self.max_iter - 1):
                print(f"  第 {t + 1:3d} 代 | w = {w:.3f} | "
                      f"gbest = {gbest_fit:.6e} | 平均 = {mean_hist[-1]:.6e}")

        return {
            "gbest": gbest, "gbest_fit": gbest_fit,
            "gbest_hist": np.array(gbest_hist),
            "mean_hist": np.array(mean_hist),
            "gbest_pos_hist": np.array(gbest_pos_hist),
            "X_hist": np.array(X_hist),
            "track": {i: np.array(v) for i, v in track.items()},
        }

# ========== 3. 单次运行演示（观察逐代变化过程） ==========
print("\n" + "=" * 72)
print("实例一：2 维 Rastrigin 函数 —— 手写 PSO（粒子数 50，迭代 200 次，seed=0）")
print("=" * 72)
pso1 = ParticleSwarm(n_particles=50, dim=2, max_iter=200, seed=0)
res1 = pso1.optimize(rastrigin, verbose=True)
print(f"\n单次运行结果：最优适应度 f* = {res1['gbest_fit']:.6e}")
print(f"最优位置 x* = ({res1['gbest'][0]:.6f}, {res1['gbest'][1]:.6f})")
print(f"与全局最优 f=0 的差距：{res1['gbest_fit'] - TRUE_OPT:.6e}")

# ========== 4. 多次独立运行：随机算法必须报告的统计指标 ==========
N_RUNS = 10                                   # 独立运行次数 M
final_bests, conv_gens, all_hist, all_res = [], [], [], []
for run_id in range(N_RUNS):
    pso = ParticleSwarm(n_particles=50, dim=2, max_iter=200, seed=run_id)
    res = pso.optimize(rastrigin)
    fb = res["gbest_fit"]
    final_bests.append(fb)
    # 收敛代数：gbest 首次进入 [最终最优值 + 1e-8] 以内的代数（见 3.3 节定义）
    cg = int(np.argmax(res["gbest_hist"] <= fb + 1e-8))
    conv_gens.append(cg)
    all_hist.append(res["gbest_hist"])
    all_res.append(res)

final_bests = np.array(final_bests)
conv_gens = np.array(conv_gens)
print("\n" + "=" * 72)
print(f"实例一：{N_RUNS} 次独立运行统计（随机算法必须报告均值与标准差）")
print("=" * 72)
print(f"每次运行最优适应度：{[f'{v:.6e}' for v in final_bests]}")
print(f"最优适应度均值 = {final_bests.mean():.6e}，标准差 = {final_bests.std(ddof=1):.6e}")
print(f"最好的一次 = {final_bests.min():.6e}，最差的一次 = {final_bests.max():.6e}")
print(f"各次运行收敛代数：{conv_gens.tolist()}，中位数 = {np.median(conv_gens):.0f} 代")
print(f"与全局最优 f=0 的差距（均值意义上）= {final_bests.mean() - TRUE_OPT:.6e}")
# 选出"最终适应度取中位数"的那次运行作为代表性运行，供逐代指标与画图使用
rep_id = int(np.argsort(final_bests)[N_RUNS // 2])
print(f"代表性运行（最终适应度排序取中位）：第 {rep_id + 1} 次（seed={rep_id}）")

# ========== 5. 代表性运行的逐代指标与代表粒子轨迹端点 ==========
rep_res = all_res[rep_id]
print("\n代表性运行逐代指标（历史最优 gbest 与种群平均适应度）：")
print(f"{'代数':>5s} | {'gbest 适应度':>14s} | {'平均适应度':>14s}")
for t in [1, 20, 50, 100, 150, 200]:
    print(f"{t:5d} | {rep_res['gbest_hist'][t - 1]:14.6e} | {rep_res['mean_hist'][t - 1]:14.6e}")

print(f"\n代表性运行最优位置 x* = ({rep_res['gbest'][0]:.6f}, {rep_res['gbest'][1]:.6f})，"
      f"适应度 = {rep_res['gbest_fit']:.6e}")
print("代表粒子的初始位置与最终位置（轨迹图图2 中的 4 条轨迹）：")
for pid, traj in rep_res["track"].items():
    print(f"  粒子 {pid + 1}：初始 ({traj[0, 0]:+.3f}, {traj[0, 1]:+.3f})"
          f" -> 最终 ({traj[-1, 0]:+.3f}, {traj[-1, 1]:+.3f})")

# ========== 6. 附加实例：5 维 Rastrigin 函数（维数升高后难度上升） ==========
print("\n" + "=" * 72)
print("实例二：5 维 Rastrigin 函数 —— 手写 PSO（粒子数 50，迭代 200 次，seed=42）")
print("=" * 72)
pso5 = ParticleSwarm(n_particles=50, dim=5, max_iter=200, seed=42)
res5 = pso5.optimize(rastrigin)
print(f"最优适应度 f* = {res5['gbest_fit']:.6e}，与全局最优 0 的差距 = {res5['gbest_fit'] - TRUE_OPT:.6e}")
print(f"最优位置 x* = {np.round(res5['gbest'], 4)}")

# ========== 7. 库实现对照：scipy.optimize.differential_evolution ==========
# 说明：差分进化（differential evolution）是 scipy 内置的另一种成熟全局优化
# 算法（进化算法类），此处用于同一函数上作参照，验证手写 PSO 结果的正确性。
print("\n" + "=" * 72)
print("对照实验：scipy.optimize.differential_evolution（差分进化，作为参照）")
print("=" * 72)
res_de = differential_evolution(
    rastrigin, bounds=[(LB, UB)] * 2, seed=42,
    tol=1e-10, polish=True, disp=False)
print(f"最优适应度 f* = {res_de.fun:.6e}，最优位置 x* = ({res_de.x[0]:.6f}, {res_de.x[1]:.6f})")
print(f"函数评估次数 = {res_de.nfev}（PSO 为 50 粒子 * 200 代 = 10000 次）")

# ========== 8. 参数敏感性：惯性权重 w 的四种取值对比 ==========
print("\n" + "=" * 72)
print("参数敏感性实验：惯性权重 w 的不同设置（粒子数 50，迭代 200，seed=0）")
print("=" * 72)
w_settings = [
    (0.9, 0.9, "w = 0.9 固定（大惯性，探索强）"),
    (0.4, 0.4, "w = 0.4 固定（小惯性，开发强）"),
    (0.7, 0.7, "w = 0.7 固定（折中）"),
    (0.9, 0.4, "w：0.9 → 0.4 线性递减（推荐）"),
]
w_curves, w_labels, w_finals = [], [], []
for ws, we, label in w_settings:
    pso = ParticleSwarm(n_particles=50, dim=2, max_iter=200,
                        w_start=ws, w_end=we, seed=0)
    res = pso.optimize(rastrigin)
    w_curves.append(res["gbest_hist"])
    w_labels.append(label)
    w_finals.append(res["gbest_fit"])
    print(f"{label:<32s} -> 最终最优适应度 = {res['gbest_fit']:.6e}")
print("结论：w 影响探索与开发的平衡，解读见第七节 7.1。")

# ========== 9. 可视化：4 张图（保存至 figures/ 目录并显示） ==========
# ---- 等高线底图（Rastrigin 在 [-5.12, 5.12]² 上的适应度曲面）----
X2D, Y2D = np.meshgrid(np.linspace(LB, UB, 400), np.linspace(LB, UB, 400))
Z2D = rastrigin(np.dstack([X2D, Y2D]))         # 形状 (400, 400) 的适应度矩阵
levels = np.arange(0, 90, 10)                  # 等高线层次

# ---- 图 1：最优适应度收敛曲线（10 次独立运行 + 均值曲线）----
fig1, ax1 = plt.subplots(figsize=(9, 5.5))
for k, h in enumerate(all_hist):
    ax1.plot(h, lw=1.2, alpha=0.55, label=f"第 {k + 1} 次运行")
mean_curve = np.mean(np.array(all_hist), axis=0)
ax1.plot(mean_curve, "k-", lw=2.5, label="10 次运行均值")
ax1.axhline(TRUE_OPT, color="r", ls="--", lw=1.2, label=f"全局最优 f = {TRUE_OPT}")
ax1.set_yscale("log")
ax1.set_xlabel("迭代次数")
ax1.set_ylabel("全局最优适应度 gbest（对数坐标）")
ax1.set_title("图1  最优适应度收敛曲线：10 次独立运行")
ax1.legend(fontsize=8, ncol=2)
ax1.grid(alpha=0.3)
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, "pso_convergence.png"), dpi=150)
print("已保存图1：figures/pso_convergence.png")

# ---- 图 2：粒子运动轨迹图（等高线 + 代表粒子轨迹 + 全局最优点）----
fig2, ax2 = plt.subplots(figsize=(7, 6))
cs = ax2.contourf(X2D, Y2D, Z2D, levels=levels, cmap="viridis", alpha=0.85)
fig2.colorbar(cs, ax=ax2, label="适应度 f(x)")
track_colors = ["tab:red", "tab:blue", "tab:orange", "tab:purple"]
for j, (pidx, traj) in enumerate(rep_res["track"].items()):
    ax2.plot(traj[:, 0], traj[:, 1], color=track_colors[j], lw=1.2, alpha=0.9,
             label=f"粒子 {pidx + 1}")
    ax2.scatter(traj[0, 0], traj[0, 1], marker="o", s=45,
                color=track_colors[j], zorder=5)          # 初始位置：圆点
    ax2.scatter(traj[-1, 0], traj[-1, 1], marker="x", s=60,
                color=track_colors[j], zorder=5)          # 最终位置：叉号
ax2.scatter(*TRUE_OPT_X, marker="*", s=320, color="gold", edgecolors="k", zorder=6,
            label=f"全局最优点 (0, 0)，f = {TRUE_OPT}")
ax2.scatter(*rep_res["gbest"], marker="^", s=90, color="white", edgecolors="k", zorder=6,
            label=f"PSO 找到的最优点 ({rep_res['gbest'][0]:.2f}, {rep_res['gbest'][1]:.2f})")
ax2.set_xlim(LB, UB)
ax2.set_ylim(LB, UB)
ax2.set_xlabel("$x_1$")
ax2.set_ylabel("$x_2$")
ax2.set_title("图2  粒子运动轨迹（圆点=初始，叉号=最终，五角星=全局最优）")
ax2.legend(fontsize=8, loc="lower right")
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, "pso_trajectory.png"), dpi=150)
print("已保存图2：figures/pso_trajectory.png")

# ---- 图 3：迭代初期 / 中期 / 后期粒子分布三面板对比 ----
fig3, axes = plt.subplots(1, 3, figsize=(15, 4.6))
snap_gens = [1, 50, 200]                    # 取第 1、50、200 代（初期/中期/后期）
snap_titles = ["迭代初期（第 1 代）", "迭代中期（第 50 代）", "迭代后期（第 200 代）"]
for ax, tg, title in zip(axes, snap_gens, snap_titles):
    ax.contourf(X2D, Y2D, Z2D, levels=levels, cmap="viridis", alpha=0.85)
    X_all = rep_res["X_hist"][tg - 1]       # 该代全体粒子的位置
    ax.scatter(X_all[:, 0], X_all[:, 1], s=10, color="white",
               edgecolors="k", linewidths=0.4, alpha=0.85, label="粒子")
    ax.scatter(*rep_res["gbest_pos_hist"][tg - 1], marker="*", s=260,
               color="gold", edgecolors="k", zorder=6, label="全局最优 gbest")
    ax.set_xlim(LB, UB)
    ax.set_ylim(LB, UB)
    ax.set_xlabel("$x_1$")
    ax.set_ylabel("$x_2$")
    ax.set_title(title)
    ax.legend(fontsize=8, loc="upper right")
fig3.suptitle("图3  粒子空间分布演化：从均匀散布到聚集于全局最优", fontsize=13)
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, "pso_distribution.png"), dpi=150)
print("已保存图3：figures/pso_distribution.png")

# ---- 图 4：惯性权重 w 不同取值下的收敛曲线对比 ----
fig4, ax4 = plt.subplots(figsize=(9, 5.5))
for curve, label in zip(w_curves, w_labels):
    ax4.plot(curve, lw=1.8, label=label)
ax4.set_yscale("log")
ax4.set_xlabel("迭代次数")
ax4.set_ylabel("全局最优适应度 gbest（对数坐标）")
ax4.set_title("图4  惯性权重 w 不同取值下的收敛曲线对比")
ax4.legend(fontsize=9)
ax4.grid(alpha=0.3)
fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, "pso_w_sensitivity.png"), dpi=150)
print("已保存图4：figures/pso_w_sensitivity.png")

plt.show()   # 显示全部图片（无图形界面环境可设 MPLBACKEND=Agg 跳过）

七、结果解读与注意事项

7.1 实例运行结果解读（与第六节程序输出对应）

（1）单次运行（seed=0）：第 1 代 gbest 适应度即达 $4.46$（50 个随机初始点中的最好者），随后阶梯式下降：第 20 代 $8.53 \times 10^{-1}$；第 60~80 代在 $3.44 \times 10^{-3}$ 附近短暂停滞（粒子群在局部极小盆地"集结"）；第 100 代 $1.83 \times 10^{-4}$；第 140 代 $6.38 \times 10^{-9}$；第 180 代 $1.28 \times 10^{-13}$；最终在第 200 代达到 $0$（双精度意义下的精确 0）。最优位置收敛到 $(-0.000000, -0.000000)$，与 Rastrigin 函数的全局最优点原点完全一致，与全局最优 $f=0$ 的差距为 0——这是"算法找到了真解"的最直接证据。

（2）10 次独立运行统计：10 次运行的最终最优适应度全部为 0，最优值均值 = 0、标准差 = 0，最好值与最差值均为 0；各次收敛代数依次为 139、116、134、141、131、142、141、141、129、139，中位数 139 代（范围 116~142）。这说明"50 粒子 × 200 代"的配置对 2 维 Rastrigin 非常充裕，算法在该参数下高度稳定、可复现；同时收敛代数约 140 代，说明 200 代的迭代预算基本"物尽其用"。注意：若粒子数或迭代次数不足（如迭代减半），就会出现部分运行陷入局部极小、均值与标准差大于 0 的情形（见 7.2 第 1 条）——随机算法必须报告均值/标准差的原因正在于此。

（3）代表性运行的逐代指标（第 6 次运行，seed=5）：gbest 从第 1 代的 $7.81$ 下降到第 200 代的 0。其中第 50 代 gbest $= 1.00$——正好停在 Rastrigin 大量局部极小（$f \approx 1$ 的盆地）上，直到第 100 代才突破到 $4.56 \times 10^{-5}$，第 150 代 $6.93 \times 10^{-11}$：这个"停滞—突破"过程正是社会认知项持续拉动群体脱离局部极小的体现。平均适应度从第 1 代的 $37.8$ 一路降到第 200 代的 $0.70$：最终仍有 1~2 个粒子停留在 $f \approx 1$ 的局部极小盆地，说明种群保留着残余多样性，这是 PSO 的正常形态——gbest 只代表"最好的那个粒子"，不等于全体粒子都已收敛。

（4）轨迹图（图②）的解读：4 条代表粒子轨迹清晰展示"先探索、后开发"的两阶段：粒子 1 从 $(-2.85, +3.80)$ 出发，先以大幅步长跨越多个等高线"盆地"，随后绕行几圈后在原点附近以小步幅缠绕收敛；粒子 8（$+3.89, -2.31$）、粒子 20（$+1.43, +4.97$）、粒子 43（$-4.09, -1.19$）同样从搜索域各处出发，最终全部落到原点。轨迹初期跨度大（惯性权重 $w$ 大、速度大），后期步幅越来越小（$w$ 递减 + 粒子相互靠拢），最终叉号与五角星（全局最优点）重合。

（5）三阶段分布图（图③）的解读：第 1 代粒子均匀散布整个搜索域（距原点距离中位数约 4.6）；第 50 代粒子已明显向中心收缩（距离中位数 1.3），最优粒子落在 $f \approx 1$ 的局部极小盆地；第 200 代 50 个粒子中 48 个进入原点附近的 $10^{-4}$ 级小盆地、距离中位数恰好为 0，仅个别粒子滞留局部极小——鸟群"聚拢到食物最丰盛处"的过程一目了然。

（6）5 维实例：同样的粒子数与迭代次数下，5 维 Rastrigin 的最优适应度为 $6.36 \times 10^{-10}$（与全局最优 0 的差距即 $6.36 \times 10^{-10}$），最优位置各维均约等于 0——算法仍然找到了全局最优。但 gbest 曲线从第 1 代的 $31.2$、第 50 代的 $4.69$、第 100 代的 $0.682$ 缓慢下降到第 200 代的 $6.72 \times 10^{-10}$，收敛代数 190，几乎用满全部迭代预算。与 2 维（精确 0）相比精度低约 10 个数量级：这正是"维数灾难"在全局优化中的体现，论文中遇到高维问题应加大粒子数与迭代次数，并如实报告差距。

（7）库对照：scipy 的差分进化在同一函数上求得 $f^* = 0$（精确 0），最优位置 $(0.000000, -0.000000)$，仅用 1743 次函数评估（PSO 为 $50 \times 200 = 10000$ 次）。两种独立算法互相印证了"Rastrigin 全局最优为 0"以及手写 PSO 结果的正确性；评估次数差异仅作参考（两者搜索机制与停止准则不同）。

（8）w 敏感性（图④）：$w = 0.9$ 固定 → $6.27 \times 10^{-3}$（最差），最优位置停在 $(0.006, 0)$——惯性太大，粒子"刹不住车"，只能在原点附近来回震荡而无法精细收敛；$w = 0.4$ 固定 → 0；$w = 0.7$ 固定 → $7.11 \times 10^{-14}$；$w$ 线性递减（0.9→0.4）→ 0。本例中除"恒为大惯性"外其余设置均命中全局最优；结合文献经验，固定小 $w$ 在多峰问题上早熟风险更高，因此惯性权重推荐采用 0.9 → 0.4 线性递减，并在论文中用图④这类敏感性实验佐证参数选择的合理性。

7.2 常见坑与对策

1. 过早收敛（早熟）：多峰函数上粒子群挤进局部最优盆地后，三个力互相平衡，很难再跳出来。对策：增大粒子数；采用 lBest 环形拓扑；$w$ 取递减方案而不是固定小值；与局部搜索/变异算子混合；多次运行取最好。
2. 速度爆炸：$w$、$c_1$、$c_2$ 取值过大时速度可能指数增长，粒子来回震荡甚至越过最优解。对策：必须加速度限幅 $v_{max}$（取搜索域宽的 10%~20%），实现中漏掉 np.clip 是新手最常见的 bug。
3. 参数不当：$w$ 恒为 0.9 收敛慢、$w$ 恒为 0.4 易早熟、$c_1 \ne c_2$ 严重失衡会破坏搜索行为。对策：先用推荐值（$w: 0.9 \to 0.4$，$c_1 = c_2 = 2$），再按图④的方法做敏感性实验微调。
4. 粒子数过少：如只用 10 个粒子搜 5 维空间，种群多样性不足，几乎必然早熟。对策：低维 2050 个，高维 50100 个；不放心就做"粒子数-最优值"敏感性实验。
5. 边界处理不当：越界粒子直接丢弃会造成有效搜索域缩小。本文采用"夹回边界"（clip）；也可用"反射"或"随机重置"，但要在论文中说明。
6. 只跑一次就下结论：换一个随机种子结果可能完全不同。对策：固定 seed 保证可复现 + 独立运行 10 次报告均值/标准差（见 3.4）。
7. 适应度量纲问题：目标函数值跨越多个数量级（如从 $10^4$ 到 $10^{-10}$）时，直接看线性坐标的收敛曲线会"看不到后期进展"，收敛曲线应使用对数坐标（本文图①④即如此）。
8. 把 PSO 当精确算法用：PSO 擅长快速定位最优解"附近"，把最后几位的精度磨出来很慢。高精度需求应接"PSO 粗搜 + 局部搜索（如 scipy 的 BFGS）精化"。

7.3 竞赛论文写作建议（话术模板）

模板 1：模型参数寻优的标准写法

为确定模型（式 (12)）中无法解析求解的参数 $\boldsymbol{\theta}$，本文采用粒子群算法进行参数寻优：以观测值与模型输出的均方误差最小为目标函数，设置粒子数 $N = 50$，最大迭代次数 $T = 200$，惯性权重 $w$ 由 0.9 线性递减至 0.4，学习因子 $c_1 = c_2 = 2$，速度上限取搜索域宽度的 20%。考虑到算法的随机性，独立运行 10 次：最优适应度均值 $\bar{f}^* = 0$，标准差 $s_{f^*} = 0$，收敛代数中位数 139（算法在约 140 代内收敛），结果表明该参数配置下算法收敛稳定、结果可复现。最终取 10 次运行中的最优解 $\boldsymbol{\theta}^*$ 作为模型参数估计。

模板 2：结果汇报表（论文中可直接使用，数值对应本文实例）

指标	数值
粒子数 $N$ / 最大迭代 $T$	50 / 200
独立运行次数 $M$	10
最优适应度均值 $\bar{f}^*$	0（10 次全部命中全局最优）
最优适应度标准差 $s_{f^*}$	0
最好值 / 最差值	0 / 0
收敛代数（各次 / 中位数）	116~142 / 139
与已知全局最优的差距	0
对照算法（差分进化）结果	$f^* = 0$（1743 次评估）

模板 3：收敛性/稳定性说明

图 × 给出了 10 次独立运行的收敛曲线（对数坐标）。可见各次运行的目标函数值均随迭代次数单调下降，并在约 140 代内先后收敛至全局最优，最终汇聚于同一水平，说明算法在给定参数下收敛稳定；图 × 的粒子空间分布演化进一步表明，粒子群在迭代后期高度聚集于全局最优解附近（48/50 个粒子进入最优解的 $10^{-4}$ 邻域），验证了算法的有效性。

写作提醒：论文中务必交代"随机种子固定 + 多次运行统计"两个细节；把 PSO 与至少一种其他算法（如 GA、差分进化）做对照，能显著提升评审对算法工作的认可度；所有图表（本文 4 张图即标准配置）应放入正文或附录并逐图引用。

八、延伸阅读

8.1 二进制粒子群算法（BPSO）

Kennedy & Eberhart (1997) 提出：把位置每一维映射为取 0/1 的概率，速度经过 sigmoid 函数 $\sigma(v) = 1/(1+e^{-v})$ 变换为"该维取 1 的概率"，再按概率随机置位。BPSO 使 PSO 能处理 0-1 背包、特征选择等离散组合问题，代价是损失了原始 PSO 的部分几何直观，且速度含义变为"概率倾向"。

8.2 带收缩因子的 PSO（Clerc PSO）

Clerc & Kennedy (2002) 对原始更新公式引入收缩因子 $\chi$：

$v_i^{t+1} = \chi\left[ v_i^t + c_1 r_1 (p_{best,i} - x_i^t) + c_2 r_2 (g_{best} - x_i^t) \right], \quad \chi = \frac{2}{\left| 2 - \varphi - \sqrt{\varphi^2 - 4\varphi} \right|} \approx 0.72984$

其中 $\varphi = c_1 + c_2 > 4$（常取 $c_1 = c_2 = 2.05$，$\varphi = 4.1$）。收缩因子在数学上等价于"惯性权重 + 速度限幅"的另一种写法，可以省略速度限幅步骤且收敛行为有理论保证，被证明在不少基准问题上优于线性递减 $w$ 的标准版。

8.3 量子粒子群算法（QPSO）

Sun 等人 (2004) 受量子力学启发提出：粒子的状态不用"位置 + 速度"描述，而是用一个吸引子（pbest 与 gbest 的加权平均）周围的概率分布（$\delta$ 势阱模型）描述，通过蒙特卡洛采样更新位置。QPSO 去掉了速度概念，只有一个可调参数，全局搜索能力通常强于标准 PSO，是论文中常用的改进点。

8.4 多目标粒子群算法（MOPSO）

标准 PSO 只能优化单目标。MOPSO（Coello Coello et al., 2004）针对多目标问题做了三处改造：用外部档案（external archive）保存搜索过程中发现的非支配解；用自适应网格维护解的分布均匀性；gbest 改为从档案中按密度（稀疏区域优先）随机选取，引导粒子均匀覆盖 Pareto 前沿。MOPSO 与 NSGA-II 类似，是竞赛中处理多目标优化问题的备选方案。

8.5 其他方向

• 自适应 PSO：让 $w$ 或 $c_1, c_2$ 随迭代进度自适应调整（如模糊 PSO、APSO）；
• 混合算法：PSO 与局部搜索（BFGS/爬山）、模拟退火、遗传算法杂交，取长补短——"PSO 粗定位 + 局部精化"是工程上最实用的组合；
• 参考来源：Kennedy & Eberhart, Particle Swarm Optimization, 1995；Shi & Eberhart, A modified particle swarm optimizer, 1998（惯性权重递减）；Clerc & Kennedy, 2002（收缩因子）。