人工蜂群演算法

小z同學發表於2021-06-06

演算法背景

人工蜂群演算法 (Artificial Bee Colony, ABC) 是由 Karaboga 於 2005 年提出的一種新穎的基於叢集智慧的全域性優化演算法,其直觀背景來源於蜂群的採蜜行為。它的主要特點是不需要了解問題的特殊資訊,只需要對問題進行優劣的比較,通過各人工蜂個體的區域性尋優行為,最終在群體中使全域性最優值突現出來,有著較快的收斂速度。

蜜蜂是一種群居昆蟲,雖然單個昆蟲的行為極其簡單,但是由單個簡單的個體所組成的群體卻表現出極其複雜的行為。真實的蜜蜂種群能夠在任何環境下,以極高的效率從食物源(花朵)中採集花蜜;同時,它們能適應環境的改變。

搜尋流程

演算法的呼叫過程如下:

初始化所有蜜源

記錄最優蜜源

while:

	僱傭蜂對所有蜜源進行鄰域搜尋(避免飢餓效應)

	計算輪盤度,判斷蜜源質量

	觀察蜂對優質蜜源進行鄰域搜尋(加速演算法收斂)

	記錄最優蜜源

	偵查蜂放棄枯竭蜜源進行全域性搜尋(跳出區域性最優)

	記錄最優蜜源

end

其中僱傭蜂和觀察蜂有著相似的邏輯,特別在對指定蜜源進行鄰域搜尋時,兩者的邏輯是完全的一樣的:

  1. 基於原有蜜源進行鄰域突變
  2. 保證鄰域突變的有效性
  3. 若為優質突變,則進行蜜源替換
  4. 若為劣質突變,則進行蜜源開採

但是演算法的設計者們卻特意區分出兩種不同的邏輯,其原因可以從實現程式碼中看出。在進行領域搜尋時,對指定蜜源的選擇和限定是關鍵所在,它暗示了僱傭蜂和觀察蜂的區別以及承擔的不同角色。

首先對於僱傭蜂的角色,其指定蜜源的方式簡單粗暴,對每一個蜜源進行遍歷指定。通過這種方式進行鄰域搜尋,是建立整個演算法的基礎核心。

而對於觀察蜂角色,它是根據輪盤賭策略進行蜜源的指定,也就是說,蜜源越是優質,其被指定的、被進行領域搜尋的概率就越高。通過這種正向反饋的過程,加速了整個演算法的收斂性,可以幫助我們在多個區域性中快速找到最優解。

如此看來觀察蜂似乎是僱傭蜂的進化版,觀察蜂似乎可以完全替代僱傭蜂?其實不然。觀察蜂角色在進行快速收斂、對優質蜜源進行了較多照顧的同時,劣質的蜜源可能會被忽略,從而產生飢餓效應。而僱傭蜂的角色讓所有的蜜源,不論優劣,都有能夠持續搜尋的機會,僱傭蜂和觀察蜂相互配合,使得演算法能夠均衡高效的執行。

最後的角色是偵查蜂。當持續的搜尋使得蜜源枯竭的時候,偵查蜂對枯竭的蜜源進行放棄,跳出原有的區域性空間,在全域性空間中隨機探索新的蜜源,之後轉變為僱傭蜂重複上述過程。

Java 實現

主函式

public class test {
	// 建立蜂群
    static beeColony bee = new beeColony();

    public static void main(String[] args) {
        int iter = 0;
        int run = 0;
        int j = 0;
        double mean = 0;
        //srand(time(NULL));

		// 重複計算,獲得平均情況下的全域性最優解
        for (run = 0; run < bee.runtime; run++) {
          	// 初始化蜜源
            bee.initial();
            bee.MemorizeBestSource();
            for (iter = 0; iter < bee.maxCycle; iter++) {
                // 僱傭蜂的邏輯: 更新所有蜜源(確保能跳出區域性蜜源)
                bee.SendEmployedBees();
                bee.CalculateProbabilities();
                // 觀察蜂的邏輯: 好蜜源正向反饋(加速收斂的作用)
                bee.SendOnlookerBees();
                bee.MemorizeBestSource();
                // 僅僅重製枯竭的蜜源 (區域性優質蜜源枯竭的快)
                bee.SendScoutBees();
            }
            // 列印最優的座標,以及最終的函式值
            for (j = 0; j < bee.D; j++) {
                //System.out.println("GlobalParam[%d]: %f\n",j+1,GlobalParams[j]);
                System.out.println("GlobalParam[" + (j + 1) + "]:" + bee.GlobalParams[j]);
            }
            //System.out.println("%d. run: %e \n",run+1,GlobalMin);
            System.out.println((run + 1) + ".run:" + bee.GlobalMin);
            bee.GlobalMins[run] = bee.GlobalMin;
            mean = mean + bee.GlobalMin;
        }
        mean = mean / bee.runtime;
        //System.out.println("Means of %d runs: %e\n",runtime,mean);
        System.out.println("Means  of " + bee.runtime + "runs: " + mean);
    }
}

演算法函式

import java.lang.Math;

public class beeColony {


    /* Control Parameters of ABC algorithm*/
    int NP = 20; /* The number of colony size (employed bees+onlooker bees)*/
    int FoodNumber = NP / 2; /*The number of food sources equals the half of the colony size*/
    int limit = 100;  /*A food source which could not be improved through "limit" trials is abandoned by its employed bee*/
    int maxCycle = 2500; /*The number of cycles for foraging {a stopping criteria}*/

    /* Problem specific variables*/
    int D = 100; /*The number of parameters of the problem to be optimized*/
    double lb = -5.12; /*lower bound of the parameters. */
    double ub = 5.12; /*upper bound of the parameters. lb and ub can be defined as arrays for the problems of which parameters have different bounds*/


    int runtime = 30;  /*Algorithm can be run many times in order to see its robustness*/

    int dizi1[] = new int[10];
    double Foods[][] = new double[FoodNumber][D];        /*Foods is the population of food sources. Each row of Foods matrix is a vector holding D parameters to be optimized. The number of rows of Foods matrix equals to the FoodNumber*/
    double f[] = new double[FoodNumber];        /*f is a vector holding objective function values associated with food sources */
    double fitness[] = new double[FoodNumber];      /*fitness is a vector holding fitness (quality) values associated with food sources*/
    double trial[] = new double[FoodNumber];         /*trial is a vector holding trial numbers through which solutions can not be improved*/
    double prob[] = new double[FoodNumber];          /*prob is a vector holding probabilities of food sources (solutions) to be chosen*/
    double solution[] = new double[D];            /*New solution (neighbour) produced by v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij}) j is a randomly chosen parameter and k is a randomlu chosen solution different from i*/


    double ObjValSol;              /*Objective function value of new solution*/
    double FitnessSol;              /*Fitness value of new solution*/
    int neighbour, param2change;                   /*param2change corrresponds to j, neighbour corresponds to k in equation v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij})*/

    double GlobalMin;                       /*Optimum solution obtained by ABC algorithm*/
    double GlobalParams[] = new double[D];                   /*Parameters of the optimum solution*/
    double GlobalMins[] = new double[runtime];
    /*GlobalMins holds the GlobalMin of each run in multiple runs*/
    double r; /*a random number in the range [0,1)*/

    /*a function pointer returning double and taking a D-dimensional array as argument */
    /*If your function takes additional arguments then change function pointer definition and lines calling "...=function(solution);" in the code*/


//	typedef double (*FunctionCallback)(double sol[D]);  

    /*benchmark functions */

//	double sphere(double sol[D]);
//	double Rosenbrock(double sol[D]);
//	double Griewank(double sol[D]);
//	double Rastrigin(double sol[D]);

    /*Write your own objective function name instead of sphere*/
//	FunctionCallback function = &sphere;

    /*Fitness function*/
    double CalculateFitness(double fun) {
        double result = 0;
        if (fun >= 0) {
            result = 1 / (fun + 1);
        } else {

            result = 1 + Math.abs(fun);
        }
        return result;
    }

    /*The best food source is memorized*/
    void MemorizeBestSource() {
        int i, j;

        for (i = 0; i < FoodNumber; i++) {
            if (f[i] < GlobalMin) {
                GlobalMin = f[i];
                for (j = 0; j < D; j++)
                    GlobalParams[j] = Foods[i][j];
            }
        }
    }

    /*Variables are initialized in the range [lb,ub]. If each parameter has different range, use arrays lb[j], ub[j] instead of lb and ub */
    /* Counters of food sources are also initialized in this function*/


    void init(int index) {
        int j;
        for (j = 0; j < D; j++) {
            r = ((double) Math.random() * 32767 / ((double) 32767 + (double) (1)));
            Foods[index][j] = r * (ub - lb) + lb;
            solution[j] = Foods[index][j];
        }
        f[index] = calculateFunction(solution);
        fitness[index] = CalculateFitness(f[index]);
        trial[index] = 0;
    }


    /*All food sources are initialized */
    void initial() {
        int i;
        for (i = 0; i < FoodNumber; i++) {
            init(i);
        }
        GlobalMin = f[0];
        for (i = 0; i < D; i++)
            GlobalParams[i] = Foods[0][i];


    }

    void SendEmployedBees() {
        int i, j;
        /*Employed Bee Phase*/
        for (i = 0; i < FoodNumber; i++) {
            /*The parameter to be changed is determined randomly*/
            r = ((double) Math.random() * 32767 / ((double) (32767) + (double) (1)));
            param2change = (int) (r * D);

            /*A randomly chosen solution is used in producing a mutant solution of the solution i*/
            r = ((double) Math.random() * 32767 / ((double) (32767) + (double) (1)));
            neighbour = (int) (r * FoodNumber);

            /*Randomly selected solution must be different from the solution i*/
            // while(neighbour==i)
            // {
            // r = (   (double)Math.random()*32767 / ((double)(32767)+(double)(1)) );
            // neighbour=(int)(r*FoodNumber);
            // }
            for (j = 0; j < D; j++)
                solution[j] = Foods[i][j];

            /*v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij}) */
            r = ((double) Math.random() * 32767 / ((double) (32767) + (double) (1)));
            solution[param2change] = Foods[i][param2change] + (Foods[i][param2change] - Foods[neighbour][param2change]) * (r - 0.5) * 2;

            /*if generated parameter value is out of boundaries, it is shifted onto the boundaries*/
            if (solution[param2change] < lb)
                solution[param2change] = lb;
            if (solution[param2change] > ub)
                solution[param2change] = ub;
            ObjValSol = calculateFunction(solution);
            FitnessSol = CalculateFitness(ObjValSol);

            /*a greedy selection is applied between the current solution i and its mutant*/
            if (FitnessSol > fitness[i]) {

                /*If the mutant solution is better than the current solution i, replace the solution with the mutant and reset the trial counter of solution i*/
                trial[i] = 0;
                for (j = 0; j < D; j++)
                    Foods[i][j] = solution[j];
                f[i] = ObjValSol;
                fitness[i] = FitnessSol;
            } else {   /*if the solution i can not be improved, increase its trial counter*/
                trial[i] = trial[i] + 1;
            }


        }

        /*end of employed bee phase*/

    }

    /* A food source is chosen with the probability which is proportioal to its quality*/
    /*Different schemes can be used to calculate the probability values*/
    /*For example prob(i)=fitness(i)/sum(fitness)*/
    /*or in a way used in the metot below prob(i)=a*fitness(i)/max(fitness)+b*/
    /*probability values are calculated by using fitness values and normalized by dividing maximum fitness value*/
    void CalculateProbabilities() {
        int i;
        double maxfit;
        maxfit = fitness[0];
        for (i = 1; i < FoodNumber; i++) {
            if (fitness[i] > maxfit)
                maxfit = fitness[i];
        }

        for (i = 0; i < FoodNumber; i++) {
            prob[i] = (0.9 * (fitness[i] / maxfit)) + 0.1;
        }

    }

    // just add choose tactic, the center logic is as same as sendEmployedBees
    void SendOnlookerBees() {

        int i, j, t;
        i = 0;
        t = 0;
        /*onlooker Bee Phase*/
        while (t < FoodNumber) {

            r = ((double) Math.random() * 32767 / ((double) (32767) + (double) (1)));
            if (r < prob[i]) /*choose a food source depending on its probability to be chosen*/ {
                t++;

                /*The parameter to be changed is determined randomly*/
                r = ((double) Math.random() * 32767 / ((double) (32767) + (double) (1)));
                param2change = (int) (r * D);

                /*A randomly chosen solution is used in producing a mutant solution of the solution i*/
                r = ((double) Math.random() * 32767 / ((double) (32767) + (double) (1)));
                neighbour = (int) (r * FoodNumber);

                /*Randomly selected solution must be different from the solution i*/
                while (neighbour == i) {
                    //System.out.println(Math.random()*32767+"  "+32767);
                    r = ((double) Math.random() * 32767 / ((double) (32767) + (double) (1)));
                    neighbour = (int) (r * FoodNumber);
                }
                for (j = 0; j < D; j++)
                    solution[j] = Foods[i][j];

                /*v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij}) */
                r = ((double) Math.random() * 32767 / ((double) (32767) + (double) (1)));
                solution[param2change] = Foods[i][param2change] + (Foods[i][param2change] - Foods[neighbour][param2change]) * (r - 0.5) * 2;

                /*if generated parameter value is out of boundaries, it is shifted onto the boundaries*/
                if (solution[param2change] < lb)
                    solution[param2change] = lb;
                if (solution[param2change] > ub)
                    solution[param2change] = ub;
                ObjValSol = calculateFunction(solution);
                FitnessSol = CalculateFitness(ObjValSol);

                /*a greedy selection is applied between the current solution i and its mutant*/
                if (FitnessSol > fitness[i]) {
                    /*If the mutant solution is better than the current solution i, replace the solution with the mutant and reset the trial counter of solution i*/
                    trial[i] = 0;
                    for (j = 0; j < D; j++)
                        Foods[i][j] = solution[j];
                    f[i] = ObjValSol;
                    fitness[i] = FitnessSol;
                } else {   /*if the solution i can not be improved, increase its trial counter*/
                    trial[i] = trial[i] + 1;
                }
            } /*if */
            i++;
            if (i == FoodNumber)
                i = 0;
        }/*while*/

        /*end of onlooker bee phase     */
    }

    /*determine the food sources whose trial counter exceeds the "limit" value. In Basic ABC, only one scout is allowed to occur in each cycle*/
    void SendScoutBees() {
        int maxtrialindex, i;
        maxtrialindex = 0;
        // 遍歷當前所有食物源,尋找最大的trail
        for (i = 1; i < FoodNumber; i++) {
            if (trial[i] > trial[maxtrialindex])
                maxtrialindex = i;
        }
        if (trial[maxtrialindex] >= limit) {
            init(maxtrialindex);
        }
    }

		// target function for calculate
    double calculateFunction(double sol[]) {
        return Rastrigin(sol);
    }

    double sphere(double sol[]) {
        int j;
        double top = 0;
        for (j = 0; j < D; j++) {
            top = top + sol[j] * sol[j];
        }
        return top;
    }

    double Rosenbrock(double sol[]) {
        int j;
        double top = 0;
        for (j = 0; j < D - 1; j++) {
            top = top + 100 * Math.pow((sol[j + 1] - Math.pow((sol[j]), (double) 2)), (double) 2) + Math.pow((sol[j] - 1), (double) 2);
        }
        return top;
    }

    double Griewank(double sol[]) {
        int j;
        double top1, top2, top;
        top = 0;
        top1 = 0;
        top2 = 1;
        for (j = 0; j < D; j++) {
            top1 = top1 + Math.pow((sol[j]), (double) 2);
            top2 = top2 * Math.cos((((sol[j]) / Math.sqrt((double) (j + 1))) * Math.PI) / 180);

        }
        top = (1 / (double) 4000) * top1 - top2 + 1;
        return top;
    }

    double Rastrigin(double sol[]) {
        int j;
        double top = 0;

        for (j = 0; j < D; j++) {
            top = top + (Math.pow(sol[j], (double) 2) - 10 * Math.cos(2 * Math.PI * sol[j]) + 10);
        }
        return top;
    }
}


參考連結

[1] https://abc.erciyes.edu.tr/
[2] https://en.wikipedia.org/wiki/Test_functions_for_optimization
[3] https://www.pianshen.com/article/729179041/
[4] https://www.jianshu.com/p/ebd436d27cf8

相關文章