烟花算法实现求解函数极值问题（一维）

本文的算法参考北京大学的Ying Tan and Yuanchun Zhu的《Fireworks Algorithm for Optimization》，实现使用java求解xsin(x)cos(2x)-2xsin(3x)+3xsin(4*x)在[0,50]求解最小值的问题。

概要:受烟花爆炸观测的启发，提出了一种新的群体智能算法，称为烟花算法（FA），用于复杂函数的全局优化。在提出的FA中，采用了两种类型的爆炸（搜索）过程，并设计了保持火花多样性的机制。为了证明FA的有效性，在九个基准测试函数上进行了大量实验，以将FA与两种粒子群优化（PSO）算法（即标准PSO和克隆PSO）进行比较。结果表明，所提出的FA在收敛速度和全局解精度方面明显优于PSO的两种变体。

Inspired by observing fireworks explosion, a novel swarm intelligence algorithm, called Fireworks Algorithm (FA), is proposed for global optimization of complex functions. In the proposed FA, two types of explosion (search) processes are employed, and the mechanisms for keeping diversity of sparks are also well designed. In order to demonstrate the validation of the FA, a number of experiments were conducted on nine benchmark test functions to compare the FA with two variants of particle swarm optimization (PSO) algorithms, namely Standard PSO and Clonal PSO. It turns out from the results that the proposed FA clearly outperforms the two variants of the PSOs in both convergence speed and global solution accuracy.

import java.util.*;

import static java.lang.Math.*;

public class FWAlg {
    static List<Firework> fireworks = null;
    static int initFireworkNums = 50;
    public static double upLimit = 50;
    public static double downLimit = 0;
    public static double ymax = 0;
    public static double ymin = Double.MAX_VALUE;
    public static double xmax = -1;
    public static double xmin = -1;
    final public static double kesi = Double.MIN_VALUE;

    final public static double a = 1;
    final public static double b = 3;
    final public static double m = 3;
    final public static double AHat = 5;
    public static int bestIndex = 0;
    public static double f(double x){
        return -(x*sin(x)*cos(2*x)-2*x*sin(3*x)+3*x*sin(4*x)-1000);//修正函数，因为是求最小值，所以要先往下拉下去1000，再沿着x轴翻转，保证适应度越大越好
    }
    static double sHat[] = null;
    static double A[] = null;
    static int i = 0;
    Random random = new Random();
    static FWAlg fwAlg = new FWAlg();
    List<Spark> sparks = new ArrayList<>();
    public static void main(String[] args) {
        initRandomFireWork();
        while(true){
            initS_A();
            fwAlg.getSparkLoc();
            fwAlg.getGaussianSpark();
            fwAlg.seleteNewSpark();
            PrintData();
            RefreshData();
        }
    }
    public static void initRandomFireWork(){
        fireworks = new ArrayList<>();
        ymax = 0;
        ymin = 0;
        for (int i = 0; i < initFireworkNums; i++) {
            double loc = random()*(upLimit-downLimit)+downLimit;//随机生成位置
            double y = f(loc);
            fireworks.add(fwAlg.new Firework(loc,y));//内部类构造方法
            if(ymax < y){//寻找最优适应度解
                bestIndex = i;
            }
            xmax = ymax < y ? loc : xmax;
            ymax = ymax < y ? y : ymax;
            xmin = ymin > y ? loc : xmin;
            ymin = ymin > y ? y : ymin;
        }
    }
    public static void RefreshData(){
        ymax = 0;
        ymin = 0;
        for (int i = 0; i < fireworks.size(); i++) {
            double loc = fireworks.get(i).getLoc();//随机生成位置
            double y = fireworks.get(i).getY();
            if(ymax < y){//寻找最优适应度解
                bestIndex = i;
            }
            xmax = ymax < y ? loc : xmax;
            ymax = ymax < y ? y : ymax;
            xmin = ymin > y ? loc : xmin;
            ymin = ymin > y ? y : ymin;
        }
    }
    public static void PrintData(){
        System.out.println("第"+i+++"轮"+"最好的xmax:"+xmax+"ymax:"+(-ymax+1000));
    }
    public static void initS_A(){
        double siSum = 0;//si分母累加的部分
        double s[] = new double[fireworks.size()];
        sHat = new double[fireworks.size()];
        A = new double[fireworks.size()];
        double ASum = 0;//a分母累加的部分
        for (int i = 0; i < fireworks.size(); i++) {//对应eq(2)
            Firework firework = fireworks.get(i);
            s[i] = ymax - firework.getY();
            siSum += s[i];
            s[i] += kesi;
            //计算增幅向量
            A[i] = firework.getY() - ymin;
            ASum += A[i];
            A[i] += kesi;
        }
        for(int i = 0; i < fireworks.size();i++){//对应eq(3) s[i]是没有m的
            s[i] = s[i]/(siSum+kesi);
            if(s[i] < a) sHat[i] = round(a * m);
            else if (s[i] > b) sHat[i] = round(b * m);
            else sHat[i] = round(s[i] * m);
            A[i] = AHat*A[i]/(ASum+kesi);
        }
    }

	//得到火花
    public void getSparkLoc(){
        sparks = new ArrayList<>();//
        for(int index = 0;index < fireworks.size();index++){
            Firework firework = fireworks.get(index);
            for (int i = 0; i < sHat[index]; i++) {
                //由于是单维度所以没有z,多维度将在后面的博客给出
                double loc = firework.getLoc();
                double h = A[index]*(random()*2+-1);
                loc+=h;
                if(loc>upLimit||loc<downLimit){//越界判断
                    loc = downLimit+abs(loc)%(upLimit-downLimit);
                }
                Spark spark = this.new Spark(loc,f(loc));
                sparks.add(spark);
            }
        }
    }
    //得到高斯变异火花
    public void getGaussianSpark(){
        for (int i = 0; i < m; i++) {
            int index = (int)round(random()*fireworks.size())%fireworks.size();//随机挑出一个烟花
            Firework firework = fireworks.get(index);
            double gaussian = random.nextGaussian();
            double loc = firework.getLoc()*gaussian;
            if(loc>upLimit||loc<downLimit){//越界判断
                loc = downLimit+abs(loc)%(upLimit-downLimit);
            }
            Spark spark = this.new Spark(loc,f(loc));
            sparks.add(spark);
        }
    }
    public void seleteNewSpark(){
        double sum = 0;
        List<Firework> newFireWorks = new LinkedList<>();
        Firework firework = fireworks.remove(bestIndex);//将最好的加进来
        newFireWorks.add(firework);
        fireworks.addAll(sparks);//将所有烟花，候选的加入候选集
        boolean isVisit[] = new boolean[fireworks.size()];
        double R[] = new double[fireworks.size()];
        for (int i = 0; i < fireworks.size(); i++) {
            R[i] = 0;
            Firework fireI = fireworks.get(i);
            for (int j = 0; j < fireworks.size(); j++) {
                Firework fireJ = fireworks.get(j);
                R[i]+=abs(fireI.getLoc()-fireJ.getLoc());//算出距离
            }
        }
        while(newFireWorks.size()<initFireworkNums){
            sum = 0;
            for (int i = 0; i < fireworks.size(); i++) {
                if(!isVisit[i]){
                    sum+=R[i];
                }
            }
            double prob = random()*sum;
            double curProb = 0;
            for(int i = 0;i < fireworks.size();i++){
                if(!isVisit[i]){
                    curProb+=R[i];
                    if(prob<curProb){//轮盘命中
                        newFireWorks.add(fireworks.get(i));
                        isVisit[i] = true;
                        break;
                    }
                }
            }
        }
        fireworks = newFireWorks;//构建完新的
    }
    public class Firework{
        double loc;
        double y;

        public Firework(double loc, double y) {
            this.loc = loc;
            this.y = y;
        }

        public double getLoc() {
            return loc;
        }

        public void setLoc(double loc) {
            this.loc = loc;
        }

        public double getY() {
            return y;
        }

        public void setY(double y) {
            this.y = y;
        }

        @Override
        public String toString() {
            return "Firework{" +
                    "loc=" + loc +
                    ", y=" + y +
                    '}';
        }
    }
    public class Spark extends Firework{
        public Spark(double loc, double y) {
            super(loc, y);
        }
    }
}

结果
经过 6 轮迭代以后得到很好的解

第0轮最好的xmax:46.56815244206834ymax:-192.84594041260607
第1轮最好的xmax:46.620654460215036ymax:-207.45448898181667
第2轮最好的xmax:46.65741883422353ymax:-213.37479319040176
第3轮最好的xmax:46.70717915472167ymax:-215.41771493754823
第4轮最好的xmax:46.70717915472167ymax:-215.41771493754823
第5轮最好的xmax:46.69039581563346ymax:-215.50884808933165

对比了一下pso和这个算法
这个算法一轮大概会有150个左右的花火，6轮下来有6*150个判断。
而pso要得到这个程度的收敛，在种群大小为50的情况下，需要50个轮次以上才能得到这种程度的解。也就是进行了50 * 50个判断，确实是有点快。