计算智能-粒子群搜索JAVA实现-学习笔记

工资十号结

一、简介

起源于对鸟群寻找食物行为的研究。
通俗可以理解为，如果要求一二维函数f(x, y)的最大/小值，先随机找一批点，让他们每个自己按一定规则移动，来试图求得目标值。
这里注意，经过证明，粒子群搜索算法是收敛的，但不一定收敛于最值或者极值。说白了就是简单有效，但不完全有效。

二、实现要素



1. 粒子群和粒子

经过通俗的讲解，我们可以了解每个粒子应该有以下变量和功能：
（1）[变量]当前粒子的坐标位置
（2）[变量]此粒子的速度（用于粒子的移动，下一个位置=当前位置+速度，一开始初始化也是随机的）
（3）[变量]此粒子历史最优值的位置（用于“按一定规则移动”，下面会讲）
（4）[变量]此粒子历史最优值（用于“按一定规则移动”，下面会讲）
（3）[功能]更新粒子的位置
（4）[功能]评估此粒子当前位置的适应值（适应值的计算是自定义的，比如求函数最大值那么直接用函数值即可，此时适应值越大越好）

2. 按一定规则移动

这里采用的规则是，对于一个粒子考虑自己走过得历史最优位置，和所有粒子的历史最优位置，进行移动。当然还有别的移动规则，就不赘述了。
这里的移动规则公式是：
v = w*v+c1*r1*(pBest-x)+c2*r2*(gBest-x)  \\ .\\ w: 惯量权重，一般取[0,1]间的数\\ c1、c2: 加速系数，通常取固定值2.0\\ r1、r2: [0,1]区间内随机数\\ pBest: 一个粒子历史最优位置\\ gBest: 所有粒子历史最优位置
这里的x、v、pBest、gBest都是向量形式。

3. 流程

要素就这么一点，整体流程就不用流程图叙述了，还不如文字看的简洁明白。
STEP1. 先随机初始化一批粒子的位置和速度
STEP2. 评估所有粒子当前位置的适应值，过程中求得每个粒子自己的历史最优值、位置，所有粒子/全局的历史最优值、位置（比如函数求最大值，评估函数是函数本身，那么值越大越好）
STEP3. 按规则更新所有粒子的位置
STEP4. 不断重复“评估记录”、“更新位置”两个功能
STEP5. 达到最高迭代次数退出

二、code（JAVA）

1. 评估函数接口、类

评估函数接口
public interface IFunction {
    public abstract double eval(double[] x);
}评估函数类
这里是求次函数的最小值。


import static java.lang.Math.*;

public class EvalFunction1 implements IFunction {

    @Override
    public double eval(double[] x) {
        double sumx = 0;
        for (double v : x) {
            sumx += pow(v, 2);
        }
        double up = pow(sin(sqrt(sumx)), 2) - 0.5;
        double down = pow(1 + 0.001*sumx, 2);
        return 0.5 + up / down;
    }
}

2. 粒子类Particle

import java.util.Random;

public class Particle {
    public static double w;
    public static double c1;
    public static double c2;
    public static int dim;
    public static int f;
    public static IFunction evalFunc;
    public static double[] range;
    public static Random random;
    public double[] x;
    public double[] v;
    public double[] pBestx;
    public double pBestFitness;

    /**构造、初始化粒子的位置和速度
     */
    Particle(){
        pBestFitness = f == 1? -1e9 : 1e9;  // 历史最优值初始化
        random = new Random();
        // 随机初始化粒子的速度和位置
        x = new double[dim];
        v = new double[dim];
        pBestx = new double[dim];
        for (int i = 0; i < dim; ++i) {
            x = random.nextDouble()*(range[1] - range[0]) + range[0];
            v = random.nextDouble()*(range[1] - range[0]) + range[0];
        }
    }

    /**评估，同时记录历史最优位置
     */
    public double evaluate(){
        double fitness = evalFunc.eval(this.x);
        if ((f==1&&fitness>pBestFitness)||(f!=1&&fitness<pBestFitness)){
            pBestFitness = fitness;
            System.arraycopy(x, 0, pBestx, 0, x.length);
        }
        return fitness;
    }

    /**更新粒子的速度和位置
     * @param gBest 全局最优位置
     */
    public void update(double[] gBest){
        double r1 = random.nextDouble();
        double r2 = random.nextDouble();
        for (int i = 0; i < v.length; i++) {
            v = w*v + c1*r1*(pBestx-x) + c2*r2*(gBest-x);
            x += v;

            // 如果出界: 速度跨越范围，位置出界
            if(v>range[1]-range[0]) v = random.nextDouble()*(range[1]-range[0]) + range[0];
            if(x>range[1]) x = range[1];
            else if(x<range[0]) x = range[0];
        }
    }
}

3. 算法类PSOAlgo

public class PSOAlgo {
    public int f;
    public int dim;
    public int step;
    public Particle[] pars;
    public double[] gBestx;
    public double gBestFitness;

    /**
     * @param parCount 粒子群数
     * @param dim x维度
     * @param ran 取值范围 [小, 大]
     * @param w 权重
     * @param c1 加速参数c1
     * @param c2 加速参数c2
     * @param flag 1 代表求最大值，-1代表求最小值
     * @param step 迭代步数
     * @param evalFunc 评估函数
     */
    PSOAlgo(int parCount, int dim, double[] ran, double w, double c1, double c2, int flag, int step, IFunction evalFunc){
        Particle.evalFunc = evalFunc;
        Particle.f = flag;
        Particle.w = w;
        Particle.c1 = c1;
        Particle.c2 = c2;
        Particle.dim = dim;
        Particle.range = ran;
        this.step = step;
        this.f = flag;
        this.dim = dim;
        this.gBestx = new double[dim];

        // 初始化粒子群
        pars = new Particle[parCount];
        for (int i = 0; i < parCount; i++) {
            pars = new Particle();
        }

        // 初次评估
        gBestFitness = f == 1? -1e9 : 1e9;  // 历史全局最优值初始化
        evalParticles();
    }

    public void evalParticles(){
        double fitness;
        for (Particle par : pars) {
            fitness = par.evaluate();
            if ((f == 1 && fitness > gBestFitness) || (f != 1 && fitness < gBestFitness)) {
                gBestFitness = fitness;
                System.arraycopy(par.x, 0, gBestx, 0, dim);
            }
        }
    }

    public void flyParticles(){
        for (Particle par : pars) {
            par.update(gBestx);
        }
    }

    public void run(){
        for (int i = 0; i < step; ++i) {
            flyParticles();
            evalParticles();
            System.out.println("第" + i+1 + "次迭代，全局最优适应值为" + gBestFitness + "解为: \n");
            for (int j = 0; j < dim; j++) {
                System.out.printf("x%d=%f,", j+1, gBestx[j]);
            }
            System.out.println();
        }
    }
}

4. 测试类TestMyPSO

public class TestMyPSO {
    public static void main(String[] args) {
        EvalFunction1 evalFunc = new EvalFunction1();
        double[] ran = new double[]{-10, 10};
        PSOAlgo mypso = new PSOAlgo(
                100, 10, ran, 0.8, 2, 2, -1, 20, evalFunc
        );
        mypso.run();
    }
}运行结果如下：

迪士尼在逃驸马

大佬的肯定是更新的动力哈哈[调皮][耶]