PSO算法概述I

来源于文献[1]，对其进行整理

起源与背景

PSO算法最简单的模型是于1987年提出的名为Boid（Bird-oid）的模型，此模型被设计用于模拟鸟的行为，也是PSO算法的雏形。Boid（Bird-oid）模型的描述如下：

1. 鸟类的每个个体都由笛卡尔坐标系中的一个点表示，该点随机分配了初始速度和位置。
2. 根据“最接近的速度匹配规则”运行程序，以使一个个体的速度与其最接近的邻居的速度相同。
3. 以相同的方式进行迭代，所有点将快速地具有相同的速度。

由于这个初代模型过于简单并偏离实际，引入了一个随机变量将其添加至速度项，使整个模拟接近实际情况。Heppner设计了一个“cornfield model”来模拟一群鸟类的觅食行为。假设平面上有一个“cornfield model”，即食物的位置，开始时鸟类随机散布在平面上。 为了找到食物的位置，它们按照以下规则移动。

1. 假设食物（cornfield）的位置坐标为$(x_0,y_0)$，单个鸟类个体的位置坐标和速度坐标分别为$(x,y)$和$(v_x,v_y)$；当前个体的位置与食物之间的距离用于衡量当前个体的性能（距食物的距离越近，性能越好，相反，性能更差。）。
2. 假设每个个体都有记忆能力，并可以记住它曾经到达的最佳位置，记为$p_{best}$。$a$是速度调节常数。$rand$表示$U[0,1]$中的随机数，可以根据以下规则设置速度项的变化：
   如果$x>p_{best}^x$，则$v_x=v_x-rand\times a$，否则，$v_x=v_x+rand\times a$。
   如果$y>p_{best}^y$，则$v_y=v_y-rand\times a$，否则，$v_y=v_y+rand\times a$。
   $p_{best}^x$表示的是个体在$x$维度上曾经到达的最佳位置。
3. 然后假设鸟群模型可以以某种方式进行通信，并且每个个体都能够知道并记住迄今为止整个群体的最佳位置，记为$g_{best}$。$b$是速度调节常数，速度更新与上述类似：
   如果$x>g_{best}^x$，则$v_x=v_x-rand\times b$，否则，$v_x=v_x+rand\times b$。
   如果$y>g_{best}^y$，则$v_y=v_y-rand\times b$，否则，$v_y=v_y+rand\times b$。
   $g_{best}^x$表示的是整个群体在$x$维度上的最佳位置。
4. 计算机仿真结果表明，当$a/b$较大时，所有个体都会迅速聚集到“食物”。 相反，如果$a/b$小，粒子将不稳定且缓慢地聚集在“食物”周围。 通过这种简单的模拟，可以发现群体可以快速找到最佳点。 受此模型的启发，Kennedy和Eberhart设计了一种进化优化算法，经过一连串的试验和错误，他们最终将基本算法固定如下：
   $$v_x=v_x+2\ast rand\ast (p_{best}^x-x)+2\ast rand\ast (g_{best}^x-x)\text{\hspace{1em}(1.1)}$$
   $$x=x+v_x\text{\hspace{1em}(1.2)}$$
   他们将每个个体抽象为一个没有质量和体积，只有速度和位置的粒子，因此他们将此算法称为“粒子群优化算法”。
   PSO算法的流程图如下所示

在连续空间坐标系下，粒子群算法用数学方法描述如下：
假设群大小为$N$，在$D$维空间中每个粒子的位置向量为$X_i=(x_{i1},x_{i2},\cdot \cdot \cdot ,x_{id},\cdot \cdot \cdot ,x_{iD})$，速度向量为$V_i=(v_{i1},v_{i2},\cdot \cdot \cdot ,v_{id},\cdot \cdot \cdot ,v_{iD})$，个体的最佳位置（即粒子历史的最佳位置）是$P_i=(p_{i1},p_{i2},\cdot \cdot \cdot ,p_{id},\cdot \cdot \cdot ,p_{iD})$，种群的最佳位置表示为$P_g=(p_{g1},p_{g2},\cdot \cdot \cdot ,p_{gd},\cdot \cdot \cdot ,p_{gD})$。 以最小化问题为例，在PSO算法的初始版本中，个体最佳位置的更新公式为：
$$p_{i,t+1}^d=\{\begin{array}{l}{x}_{i,t+1}^d,iff(X_{i,t+1}^{})<f(P_{i,t})\\ p_{i,t}^d,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\text{\hspace{1em}(1.1)}$$
而群体的最佳位置是所有个体的最佳位置。

个体的速度和位置的更新公式分别表示如下：

$$\begin{array}{rl}{v}_{i,t+1}^d=& v_{i,t}^d+c_1\ast \mathrm{rand}\ast \left(p_{i,t}^d-x_{i,t}^d\right)+c_2\ast \mathrm{rand}\ast \left(p_{g,t}^d-x_{i,t}^d\right)\end{array}\text{\hspace{1em}(1.3)}$$

$$x_{i,t+1}^d=x_{i,t}^d+v_{i,t+1}^d\text{\hspace{1em}(1.4)}$$

由于上述初始版本的PSO在优化问题上不是很有效，因此在提出初始算法不久后，1998年Shi和Eberhart提出了一种改进的PSO算法，他们将惯性权重$\omega$引入进速度更新公式，新速度更新公式变为
$$\begin{array}{rl}{v}_{i,t+1}^d=& \omega \ast v_{i,t}^d+c_1\ast \mathrm{rand}\ast \left(p_{i,t}^d-x_{i,t}^d\right)+c_2\ast \mathrm{rand}\ast \left(p_{g,t}^d-x_{i,t}^d\right)\end{array}\text{\hspace{1em}(1.5)}$$

尽管此改进算法的复杂度与初始版本几乎相同，但它极大地提高了算法性能。 因此，它获得了广泛的应用。 通常，将修改后的算法称为规范PSO算法（canonical PSO algorithm），而初始版本称为原始PSO算法（original PSO algorithm）。

通过分析PSO算法的收敛行为，在2002年Clerc和Kennedy引入了PSO算法的一种变体，其压缩因子为$\chi$，从而确保了收敛性并提高了收敛速度。 然后，速度更新公式变为
$$\begin{array}{r}{v}_{i,t+1}^d=\chi (v_{i,t}^d+{\varphi }_1\ast \mathrm{rand}\ast \left(p_{i,t}^d-x_{i,t}^d\right)+{\varphi }_2\ast \mathrm{rand}\ast \left(p_{g,t}^d-x_{i,t}^d\right))\end{array}\text{\hspace{1em}(1.6)}$$
显然，公式（1.5）和公式（1.6）之间没有本质区别。 如果选择了适当的参数，则两个公式相同。

PSO算法有两个版本，分别称为全局版本和局部版本。在全局模型中，粒子跟踪的两个极端是其自身的最优位置$p_{best}$和群的最优位置$g_{best}$。
因此，在局部版本中，除了跟踪其自身的最优位置$p_{best}$外，粒子不跟踪种群的最优位置$g_{best}$，而是跟踪其拓扑邻域中所有粒子的最优位置$n_{best}$。对于本地版本，速度更新方程（1.5）变为：
$$\begin{array}{r}{v}_{i,t+1}^d=\omega \ast v_{i,t}^d+c_1\ast \mathrm{rand}\ast \left(p_{i,t}^d-x_{i,t}^d\right)+c_2\ast \mathrm{rand}\ast \left(p_{l,t}^d-x_{i,t}^d\right)\end{array}\text{\hspace{1em}(1.7)}$$
其中$p_l$是局部邻域中的最佳位置。

粒子的迭代方案

在每一代中，任何粒子的迭代过程如上图所示。从社会学的角度分析速度更新公式，可以看出在这个更新公式中，第一部分是粒子先前速度的影响。这意味着粒子对其当前的运动状态充满信心，并根据其自身的速度进行惯性运动，因此参数$\omega$称为惯性权重。 第二部分取决于粒子当前位置与其自身最佳位置之间的距离，称为“认知”项。它是指粒子自身的思维，即粒子自身的经验所产生的运动。因此，参数$c1$被称为认知学习因子（也称为认知加速因子）。第三部分依赖于粒子当前位置与群中全局（或局部）最优位置之间的距离，称为“社会”因素。它是指粒子之间的信息共享与协作，即粒子的运动来源于粒子群中其他粒子的经验。它通过认知来模拟好粒子的运动，因此参数$c2$被称为社会学习因子（也称为社会加速因子）。

PSO是一种随机并行优化算法。其优点可以概括如下：它不需要优化函数微分、导数和连续；其收敛速度快；算法简单，易于通过编程实现。不幸的是，它也有一些缺点${}^{[2]}$：

1. 对于具有多个局部极值的函数，它可能落入局部极值并且无法获得正确的结果。导致这种现象的原因有两个：一个是优化函数的特性，另一个是粒子的多样性迅速消失，导致过早收敛。这两个因素通常是密不可分的。
2. 由于缺乏良好的搜索方法的配合，PSO算法无法获得满意的结果。原因是PSO算法没有充分利用在计算过程中获得的信息。相反，在每次迭代期间，它仅使用群体最优值和个体最优值的信息。
3. 尽管PSO算法提供了全局搜索的可能性，但它不能保证收敛到全局最优值。
4. PSO算法是一种元启发式仿生优化算法，到目前为止还没有严格的理论基础。它仅通过简化和模拟某些群体的搜索现象进行设计，但是它既没有从原理上解释为什么该算法有效，也没有指定其适用范围。

因此，PSO算法通常适用于一类高维且不需要获得非常准确的解决方案的优化问题。

目前对PSO算法的研究有很多种，可分为以下八类：

1. 对PSO算法进行理论分析，试图了解其工作机理。
2. 改变它的结构，努力获得更好的性能。
3. 研究了不同参数配置对PSO算法的影响。
4. 研究了不同拓扑结构对PSO算法的影响。
5. 研究了并行PSO算法。
6. 研究了离散PSO算法。
7. 研究了PSO算法在多目标优化中的应用。
8. 将粒子群算法应用于各个工程领域。

以下是利用PSO解决一个连续函数的Java代码，目标函数是最小化$f(x,y)=x^3+2x{y}^2-15x-12y$ ${}^{[3]}$

import java.util.ArrayList;

public class PSOInContinuousFunction {
    int particleNumber = 50;						//粒子数量
    public ArrayList<ParticleProperties> positionArray = new ArrayList<>(particleNumber);			//位置数组
    public ArrayList<ParticleProperties>velocityArray = new ArrayList<>(particleNumber);			//速度数组
    public ArrayList<ParticleProperties>pBest = new ArrayList<>(particleNumber);		//粒子最优解
    ParticleProperties gBest = new ParticleProperties(0.0, 0.0,Double.MAX_VALUE);		//全局最优解
    double velocityMax = 0.1;				//最大速度
    double velocityMin = -0.1;
    int c1=3;                   //认知学习因子
    int c2=3;					//社会学习因子

    public static void main(String[] args) {
        PSOInContinuousFunction e = new PSOInContinuousFunction();
        e.ParticleSwarmOptimization(200);
    }

    //适应度函数
    public void calculationFitnessFunction() {
        for(int i = 0; i< particleNumber; i++) {
            positionArray.get(i).functionValue = (positionArray.get(i).x* positionArray.get(i).x* positionArray.get(i).x + 2*positionArray.get(i).x* positionArray.get(i).y* positionArray.get(i).y - 15* positionArray.get(i).x -12* positionArray.get(i).y);
        }
    }
    public PSOInContinuousFunction()
    {
        for(int i = 0; i< particleNumber; i++)
        {
            positionArray.add(i,new ParticleProperties(Math.random(),Math.random(),Double.MAX_VALUE));
            velocityArray.add(i,new ParticleProperties(Math.random()* velocityMax,Math.random()* velocityMax,Double.MAX_VALUE));
        }
        calculationFitnessFunction();
        //分别初始化粒子与粒子群体的最优值
        gBest.functionValue = Integer.MAX_VALUE;
        for(int i = 0; i< particleNumber; i++)
        {
            pBest.add(i, positionArray.get(i));
            if(positionArray.get(i).functionValue < gBest.functionValue)
            {
                gBest = new ParticleProperties(positionArray.get(i));
            }
        }
        System.out.println("初始解 -- gBest.x = "+ gBest.x+" gBest.y = "+ gBest.y+" "+" gBest.f = "+ gBest.functionValue);
    }
    //粒子群主体算法	iterationNumber-迭代次数
    public void ParticleSwarmOptimization(int iterationNumber) {
        for(int i=0;i<iterationNumber;i++) {
            double w = 0.2;//惯性权重
            for(int j = 0; j< particleNumber; j++) {
                //更新粒子中速度向量中的每一个元素
                double elementFirst =w * velocityArray.get(j).x+c1*Math.random()*(pBest.get(j).x- positionArray.get(j).x)+c2*Math.random()*(gBest.x- positionArray.get(j).x);
                double elementSecond =w * velocityArray.get(j).y+c1*Math.random()*(pBest.get(j).y- positionArray.get(j).y)+c2*Math.random()*(gBest.y- positionArray.get(j).y);

                if (elementFirst > velocityMax)
                {
                    elementFirst = velocityMax;
                }else if(elementFirst < velocityMin)
                {
                    elementFirst = velocityMin;
                }
                if (elementSecond > velocityMax)
                {
                    elementSecond = velocityMax;
                }else if(elementSecond < velocityMin)
                {
                    elementSecond = velocityMin;
                }

                velocityArray.set(j,new ParticleProperties(elementFirst, elementSecond, Double.MAX_VALUE));

                //更新粒子的位置
                positionArray.get(j).x += velocityArray.get(j).x;
                positionArray.get(j).y += velocityArray.get(j).y;
            }

            //计算每个粒子的适应度
            calculationFitnessFunction();

            //更新个体极值和群体极值
            for(int j = 0; j< particleNumber; j++) {
                if(positionArray.get(j).functionValue < pBest.get(j).functionValue) {
                    pBest.set(j, positionArray.get(j));
                }
                if(positionArray.get(j).functionValue < gBest.functionValue) {
                    gBest = new ParticleProperties(positionArray.get(j));
                }
            }
            System.out.println("第"+(i+1)+"代 -- gBest.x = "+ gBest.x+" gBest.y = "+ gBest.y+" "+" gBest.f = "+ gBest.functionValue);
        }
    }
}

//粒子的状态类
class ParticleProperties {
    public double x;//x变量
    public double y;//y变量
    public double functionValue = 0.0;//适应度，即求解函数值

    ParticleProperties(double x, double y, double functionValue) {
        this.x = x;
        this.y = y;
        this.functionValue = functionValue;
    }

    public ParticleProperties(ParticleProperties particleProperties) {
        this.x = particleProperties.x;
        this.y = particleProperties.y;
        this.functionValue = particleProperties.functionValue;
    }
}

参考文献

[1] Wang D, Tan D, Liu L. Particle swarm optimization algorithm: an overview[J]. Soft Computing, 2018, 22(2): 387-408.

[2] Wang W. Research on particle swarm optimization algorithm and its application[J]. Southwest Jiaotong University, Doctor Degree Dissertation, 2012: 36-37.

[3] https://blog.csdn.net/qq_27124771/article/details/80945337