粒子群优化算法及其java实现

Posted 2021-12-02 numen-fan

憋了两周终于把开题报告憋出来了，再一次证明自己不适合搞学术，哎……，花了点时间把报告中提到的粒子群算法看了看，看了些资料，用java跑起来。

算法简介

粒子群算法最先由Barnhart博士和Kennedy博士于1995 年提出，是一种源于对鸟群捕食行为的研究而发明的进化计算技术，原理是模仿鸟群寻觅食物的搜索过程，设想鸟群在一定区域搜寻食物，在不知道食物确切位置的情况下，鸟群依靠群体中个体判断距离食物的远近程度来调节飞行方向和飞行速度，最终通过群体的经验和自身记忆的智慧找到食物。

算法原理

算法描述

算法流程图

算法的实现(java)

• Particle.java文件

package com.jiajia.pso;

import java.util.Random;

/**
 * @ClassName: Particle
 * @Author: fanjiajia
 * @Date: 2019/5/13 上午11:01
 * @Version: 1.0
 * @Description:
 */
public class Particle {

    //维数
    public  int dimension = 2;

    //粒子的位置
    public double[] X = new double[dimension];

    //局部最好位置
    public double[] pbest = new double[dimension];

    //粒子的速度
    public double[] V = new double[dimension];

    //最大速度
    public double Vmax = 2;

    //适应值
    public double fitness;

    /**
     * 根据当前位置计算适应值
     * @return newFitness
     */
    public double calculateFitness() {
        //1.Ackley's function:
        //double newFitness = -20*Math.pow(Math.E,(-0.2*Math.sqrt(0.5*(X[0]*X[0]+X[1]*X[1]))))-Math.pow(Math.E,(0.5*(Math.cos(2*Math.PI*X[0])+Math.cos(2*Math.PI*X[1]))))+Math.E+20;

        //2.Sphere function
        //double newFitness = X[0]*X[0]+X[1]*X[1];

        //3.Rosenbrock function
        double newFitness = 100*(Math.pow((X[1]-X[0]*X[0]),2))+Math.pow((X[0]-1), 2);

        return newFitness;
    }


    /**
     * 初始化自己的位置和pbest
     */
    public void initialX() {
        for(int i=0;i<dimension;i++) {
            X[i] = new Random().nextInt(50);
            pbest[i] = X[i];
        }
    }
    /**
     * 初始化自己的速度
     */
    public void initialV() {
        for(int i=0;i<dimension;i++) {
            double tmp = new Random().nextDouble();//随机产生一个0~1的随机小数
            V[i] = tmp*4+(-2);
        }
    }
}

• PSO.java

package com.jiajia.pso;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;


import java.io.File;
import java.io.FileWriter;
import java.io.IOException;

/**
 * @ClassName: PSO
 * @Author: fanjiajia
 * @Date: 2019/5/13 上午11:02
 * @Version: 1.0
 * @Description:
 */
public class PSO {

    private static double[] gbest;//全局最优位置

    private static double gbest_fitness = Double.MAX_VALUE;//全局最优位置对应的fitness

    private static int particle_num = 20;//粒子数

    private static int N = 500;//迭代次数

    private static int c1,c2 = 2;

    private static double w = 1.4;//惯性因子

    private static List<Particle> particles = new ArrayList<Particle>();//粒子群

    private static List<Double> fittessList = new ArrayList<>(N);

    /**
     * 主程序入口
     * @param args
     */
    public static void main(String[] args) {
        process();
    }

    /**
     * 初始化所有粒子
     */
    public static void initialParticles() {
        for(int i=0;i<particle_num;i++) {
            Particle particle = new Particle();
            particle.initialX();
            particle.initialV();
            particle.fitness = particle.calculateFitness();
            particles.add(particle);
        }
    }

    /**
     * update gbest
     */
    public static void updateGbest() {
        double fitness = Double.MAX_VALUE;
        int index = 0;
        for(int i=0;i<particle_num;i++) { // 找到群体中适应值最小的粒子
            if(particles.get(i).fitness<fitness) {
                index = i;
                fitness = particles.get(i).fitness;
            }
        }
        if(fitness<gbest_fitness) { // 如果个体适应值小于全局适应值，更新全局的最优值为个体最优值
            gbest = particles.get(index).pbest.clone();
            gbest_fitness = fitness;
        }
    }

    /**
     * 跟新每个粒子的速度
     */
    public static void updateV(int n) {
        for(Particle particle:particles) {
            for(int i=0;i<particle.dimension;i++) {
                double v =(0.9 - n*(0.9-0.4)/N) * particle.V[i]+c1*rand()*(particle.pbest[i]-particle.X[i])+c2*rand()*(gbest[i]-particle.X[i]);
                if(v>particle.Vmax) // 判断速度是否超过最大的速度
                    v = particle.Vmax;
                else if(v<-particle.Vmax) // 比最大速度的相反数小
                    v = -particle.Vmax;
                particle.V[i] = v;//更新Vi
            }
        }
    }

    /**
     * 更新每个粒子的位置和pbest
     */
    public static void updateX() {
        for(Particle particle:particles) {
            for(int i=0;i<particle.dimension;i++) {
                particle.X[i] = particle.X[i] + particle.V[i];
            }
            double newFitness = particle.calculateFitness();//新的适应值
            //如果新的适应值比原来的小则跟新fitness和pbest
            if(newFitness<particle.fitness) {
                particle.pbest = particle.X.clone();
                particle.fitness = newFitness;
            }
        }
    }

    /**
     * 算法主要流程
     */
    public static void process() {
        int n = 0;
        initialParticles();
        updateGbest();
        while(n++<N) {
            updateV(n);
            updateX();
            updateGbest();
            fittessList.add(gbest_fitness);
            System.out.println(n+".当前gbest:("+gbest[0]+","+gbest[1]+")  fitness="+gbest_fitness);
        }
        write2File();
    }
    
    /**
     * 返回一个0~1的随机数
     * @return
     */
    public static double rand() {
        return new Random().nextDouble();
    }
}

代码参考了其他资料，后面有说明，但是对其中部分进行了改进。

在Particle（粒子类）中设定了三个适应函数，Ackley,Sphere,Rosenbrock,关于这三个函数的介绍可以参考https://en.wikipedia.org/wiki/Test_functions_for_optimization，这里面列出来了很多优化测试函数，很多的paper在设计了优化策略或者改进相应的优化策略之后，都会利用其中的函数进行测试。
这里用到的函数是Rosenbrock:

可以看出这里最小值在（1，1，，，1）处取的。
为了看到相应的效果，这里将全局适应值写到txt文件中，并利用python绘制出来（莫名的感觉繁琐，要是python，哪有这么麻烦，只可惜最后的实验都是java写）。

上面是$w$为1时，即惯性系数为1时的收敛结果，可以看出，算法前期搜索很快，后期较慢，且偶尔会陷入局部最优解里面。

惯性因子w的优化

惯性因子$w$代表受上一次粒子速度的影响程度，$w$越大，收敛越快，但容易错过最优解。$w$越小，收敛较慢，容易陷入局部最优解，出不来。因此改进$w$成为很多改进的焦点，其中采用较多的是Shi建议的线性递减权值策略，通常将$w$设定在[0.4，0.9]之间：

采用线性递减权值策略后得到的收敛效果：

可以看出前期收敛直线下降，且不容易陷入局部最优，最后达到全局收敛。
除了上面的线性递减权值策略，还有自适应权值策略，随机权重策略。详见参考资料[4];

参考资料

• 粒子群算法原理及Matlab实现
• 基本PSO算法实现(Java)
• Y.Shi. A Modified Particle Swarm Optimizer. 1998
• Test functions for optimization

• 武装. 几种改进的智能优化算法及其应用[M].科学技术文献出版社, 2018.

  最后

  生命不息，生活好难！！！