浅谈PSO算法（粒子群优化算法）

本文最初写于 2021年11月29日。

一、算法流程

粒子群优化算法（Particle Swarm Optimization, PSO）是一种启发式进化计算技术，来源于对简化社会群体智能行为模型的模拟，是由Kennedy和Eberhart提出的一种进化计算方法。PSO算法属于万能启发式搜索算法，能够在没有得知太多信息的情况下，有效的搜索具有庞大解空间的问题并找到候选解，但同时不保证其找到的最佳解为真实的最佳解。

粒子群优化算法模仿昆虫、兽群、鸟群和鱼群等的群集行为这些群体按照一种合作的方式寻找食物群体中的每个成员通过学习它自身的经验和其他成员的经验来不断改变其搜索模式。

具体地说，PSO算法是对生物种群觅食行为的仿真算法描述了数量为 $N$ 的生物在空间中以不同的速度运动，例如每只鸟的飞行速度和运动方向都依赖于自身的经验和群体的经验，对在每一维空间中的速度和运动方向调整。

粒子在多维空间中改变速度以及下一次运动方向的公式如下： \[ v_{\mathrm{id}}^{k}=w\, v_{\mathrm{id}}^{k-1}+c_1 r_1\left(\mathrm{pbest}_{\mathrm{id}}-x_{\mathrm{id}}^{k-1}\right)+c_2 r_2\left(\mathrm{gbest}_{\mathrm{id}}-x_{\mathrm{id}}^{k-1}\right) \] 位置更新公式为 \[ x_{i d}^k=x_{i d}^{k-1}+v_{i d}^{k-1} \] 其中 $v_{\mathrm{id}}^k$ 为第 $k$ 次迭代粒子 $i$ 飞行速度矢量的第 $d$ 维分量，$x_{i d}^k$ 为第 $k$ 次迭代粒子 $i$ 位置矢量的第 $d$ 维分量，$c_1, c_2$ 为学习因子，$r_1, r_2$ 为取值 $[0,1]$ 间的随机函数，$w$ 为惯性权重，$\mathrm{pbest}$ 为粒子 $i$ 个体经历过的最好位置，$\mathrm{gbest}$ 为种群所经历过的最好位置。

通常情况下，第 $d$ 维的位置变化范围限定在 $\left[\min\{x_{\mathrm{id}}\}, \max\{x_{\mathrm{id}}\}\right]$ 内；

同时，第 $d$ 维的速度变化范围限定在 $\left[-\max\{v_{\mathrm{id}}\}, \max\{v_{\mathrm{id}}\}\right]$ 内。

算法流程：

初始化一群微粒，包括随机的位置和速度
评价每个微粒的适应度
根据适应度更新 pbest 和 gbest，并更新粒子的位置和速度
未达到结束条件（通常为足够好的适应值）则回到2 。

流程图（有点糊）

二、例题

模板题：洛谷P3382 【模板】三分法

题意：

如题，给出一个 $N$ 次函数，保证在范围 $[l, r]$ 内存在一点 $x$，使得 $[l, x]$ 上单调增，$[x, r]$ 上单调减。试求出 $x$ 的值。

输入格式：

第一行一次包含一个正整数 $N$ 和两个实数 $l, r$，含义如题目描述所示。

第二行包含 $N + 1$ 个实数，从高到低依次表示该 $N$ 次函数各项的系数。

输出格式：

输出为一行，包含一个实数，即为 $x$ 的值。若你的答案满足以下二者之一，则算正确： - 你的答案 $x'$ 与标准答案 $x$ 的相对或绝对误差不超过 $10^{-5}$。 - 你的答案 $x'$ 与标准答案 $x$ 对应的函数值，即 $f(x') $ 和 $f(x)$ 的相对或绝对误差不超过 $10^{-5}$。

数据范围：

$6 \le N \le 13$，函数系数均在 $[-100,100]$ 内且至多 $15$ 位小数，$|l|,|r|\leq 10$ 且至多 $15$ 位小数。$l\leq r$。

考虑用PSO算法求解该问题。没什么特别的技巧，直接放代码了。

#include <bits/stdc++.h>
using namespace std;
const int cnt=100;
const double c1=2,c2=1;
double w=0.7;
int n;
double l,r,a[15];
double cal(double x)
{
	double ans=0;
	for(int i=0; i<=n; i++)
		ans=ans*x+a[i];
	return ans;
}
double Rand(){return (double)rand()/RAND_MAX;}
struct node
{
	double v,x,y,pval,pbest;
}b[105];
double gval=-1e233,gbest;
void update(int i)
{
	b[i].v=w*b[i].v+c1*Rand()*(gbest-b[i].x)+c2*Rand()*(b[i].pbest-b[i].x);
	b[i].x+=b[i].v;
	if(b[i].x<l)b[i].v*=-1,b[i].x=l;
	if(b[i].x>r)b[i].v*=-1,b[i].x=r;
	if(abs(b[i].v)>r)b[i].v=0;
	b[i].y=cal(b[i].x);
	if(b[i].y>b[i].pval)
	{
		b[i].pval=b[i].y;
		b[i].pbest=b[i].x;
	}
}
signed main()
{
	scanf("%d%lf%lf",&n,&l,&r);
	for(int i=0; i<=n; i++)
		scanf("%lf",&a[i]);
	srand(time(0));
	for(int i=1; i<=cnt; i++)
	{
		b[i].x=b[i].pbest=l+Rand()*(r-l);
		b[i].v=0;
		b[i].y=b[i].pval=cal(b[i].x);
		if(gval<b[i].y)
		{
			gval=b[i].y;
			gbest=b[i].x;
		}
	}
	for(;w>=0.3;w-=0.004)
		for(int i=1; i<=cnt; i++)
		{
			update(i);
			if(gval<b[i].pval)
			{
				gval=b[i].pval;
				gbest=b[i].pbest;
			}
		}
	printf("%.5lf\n",gbest);
	return 0;
}