布谷鸟算法浅谈与简单应用

简介

布谷鸟算法是由剑桥大学Xin-She Yang教授和S.Deb于2009年提出的一种新兴的启发算法，是一种通过模拟自然界当中布谷鸟（也就是杜鹃，故该算法也称为杜鹃算法）在繁育后代的行为而提出的一种搜索算法。

本文章将以在工程实践当中的生产车间调度问题为例，对该算法进行一个简单介绍，并给出基于JAVA语言的实现。

原型及流程

布谷鸟算法实际上来源于自然界当中布谷鸟繁育后代的行为。

自然界当中布谷鸟布谷鸟在繁育下一代当中无固定配偶固定配偶，也不自己营巢和孵卵，而是将卵产于大苇莺、麻雀、灰喜鹊、伯劳、棕头鸦雀、北红尾鸲、棕扇尾莺等各类雀形目鸟类巢中，由这些鸟替它带孵带育。

自然，在这些寄生的鸟巢中，有的后代长势较好，有的后代长势较差，也存在部分后代被宿主发现而清除的情况。

后代变为本代后，会在出生点周围随机寻找新的巢穴继续重复繁衍工作。

而布谷鸟算法的核心实际上就是对上述流程的模拟，对其中的一些关键问题进行了一定的约束：

一只布谷鸟只产生一个后代。
后代清除：后代有一定几率会被清除，需要在被清除的后代周围的巢穴重新繁育保持种群数量稳定。
本代寻找新的巢穴的策略：按照莱维飞行规律，进行移动
每一代当中长势最好的后代保留，及仍在上一代繁衍的巢穴中继续繁衍后代。
总种群数量固定。

简单的抽象为带有流程的语言，可以描述为一下步骤：

有一群初始布谷鸟在宿主鸟的鸟巢下了一组蛋
计算这些蛋的孵化和生长情况，记录下一代中适应度最好的后代，同时此时下一代适应度最好的个体也是历代中适应度最好的个体。
下一代生长为本代，本代依据一定的规律（莱维飞行）寻找一组新鸟巢下蛋。注意，本代中在生长期最好的那一只布谷鸟仍在上一代下蛋的地点下蛋。
按照c的概率，一些蛋被宿主发现并被清除。对于这些蛋，本代需要在这些蛋的周围一些的鸟巢下蛋从而保持总种群数量稳定。
记录这些蛋的孵化和生长情况，记录本代中适应度最好的后代，若适应度最好的后代比历代中所有的个体适应度都好，则更新历代最好适应度个体。
若迭代次数超过规定上线，输出历代中适应度最好的个体并结束，否则返回3进入下一代。

在上面的描述当中，实际的下蛋可以被我们理解成为当前解，而“下一代中适应度最好的后代”，则可以被我们理解为当前迭代周期中的最优解，“历代中最好适应度个体”则可以被抽象为全局最优解，即我们想要最终求得的答案。

上述步骤种第四步骤随机抛弃部分解与第五步骤计算适应度的过程，在网上资料中有两种观点。第一种是先计算适应度后抛弃，此种需要在抛弃时避免将本代中适应度最好的解抛弃，第二种是先随机抛弃后计算适应度。实际验证下，从结果来看作者暂未发现两种流程的差异，故两种均应该可行。

莱维飞行

在每轮的后代中，除了最好的后代会继续在当前“巢穴”（当前最优解）中继续繁衍后代外，其余的后代都需要寻找新的“巢穴”，具体在布谷鸟算法中的实际方法是以莱维飞行的路径搜索方案对每组解进行更新。

什么是莱维飞行呢？莱维飞行实际上就是是自然界中鸟类飞行路径，其特点简单来说就是，在一定的移动次数下长移动少，短移动次数多。

若不能理解，想象一下自然界当中猛禽捕食搜索猎物的时候的场景：在一个小区域内短距离频繁的移动搜索猎物，之后突然会移动到离当前距离很远的另一个小区域，之后继续在这个区域进行短距离频繁移动搜索猎物。莱维飞行时的移动与之类似，就是短距离多次移动，长距离少次移动。

莱维飞行的稳定分布积分形式在1937年由P. Levy确定：
$\frac{1}{\pi} \int_0^{+\infty}exp(-\beta |k|^\lambda)cos(ks) dk$

但对该积分并未有一个明确的解析，不过当 $s\gg s_0 >0，即s \rightarrow \infty$ 时:

$\approx \frac{\lambda \beta \Gamma(\lambda)sin(\frac{\pi \lambda}{2})}{\pi}.\frac{1}{s^{1+\lambda}}$

针对上述的近似解，历史上有许多学者根据此提出用于服从莱维分布的随机数实现方法，其中在布谷鸟算法种用于实现的方法主要是采用Mantegna在1994年提出的一种用正态分布求解随机数的方法，有时候也叫Mantegna方法，生成满足莱维分布的随机步长的方法如下：
$s=\frac{u}{|v|^{\frac{1}{\beta}}}$
其中 $u～N(0,\sigma^2) ,v～N(0,1)$ , $\sigma = \{ \frac{\Gamma(1+\beta)sin(\frac{\pi \beta}{2})}{ \Gamma(\frac{1+\beta}{2})2^{\frac{\beta-1}{2}}} \}^{\frac{1}{\beta}}$ ，其中通常 $\beta$ 取值在1到2之间，为1.5时算法效率最好。

下面给出一个在Matlab种用Mantegna方法模拟二维平面莱维飞行的代码：

% Mantegna方法模拟萊维飞行
% author zhaoyuqiang 
% 参考：https://blog.csdn.net/zyqblog/article/details/80905019
x = [0,0];
y = [0,0];
beta = 1.5;
sigma_u = (gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
sigma_v = 1;
for i=1:1000u = normrnd(0,sigma_u);v = normrnd(0,sigma_v);s = u/(abs(v))^(1/beta);x(:,1) = x(:,2);x(:,2) = x(:,1)+1*s;u = normrnd(0,sigma_u);v = normrnd(0,sigma_v);s = u/(abs(v))^(1/beta);y(:,1) = y(:,2);y(:,2) = y(:,1)+1*s;plot(x,y);hold on;
end
axis square;

结果如下图：

不难发现，其与我们之前描述的路径一致，即为短距离移动多，长距离移动少。

使用布谷鸟算法解决作业车间调度问题（JSP）

下面介绍布谷鸟算法在作业车间调度问题上的实际应用。

问题描述

作业车间调度问题（Job Shop Scheduling, JSP）是最经典的几个NP-hard问题之一。其应用领域极其广泛，涉及航母调度，机场飞机调度，港口码头货船调度，汽车加工流水线等。

JSP问题描述：一个加工系统有M台机器，要求加工N个作业，其中，作业i包含工序数为Li。令 $L=\begin{matrix} \sum_{i=1}^N k^2 \end{matrix}$ ，则L为任务集的总工序数。其中，各工序的加工时间已确定，并且每个作业必须按照工序的先后顺序加工。调度的任务是安排所有作业的加工调度排序，约束条件被满足的同时，使性能指标得到优化。

作业车间调度在本系统中需要考虑如下约束:

$Cons_1$ ：每道工序在指定的机器上加工，且必须在其前一道工序加工完成后才能开始加工；
$Cons_2$ ：某一时刻1台机器只能加工1个作业；
$Cons_3$ ：每个作业只能在1台机器上加工1次；
$Cons_4$ ：各作业的工序顺序和加工时间已知，不随加工排序的改变而改变。

本问题要求算法目标为所有作业加工完成总时间最短

问题的编码与解码

若由n个工件，所有工件最多需要m道工序完成，则JSP的一个可行序列可以被描述成一个 $n * m$ 的整数串。本文先通过一定的操作生成一系列初始编码值，后通过对于编码值进行解码来获取满足约束的可行解序列，通过对于可行解序列的好坏适应度情况评估来对编码值进行算法操作。

下面使用一个 $n = 2$ , $m = 3$ 的示例来对编码和解码进行说明。

首先在定义域内随机生成一个 $2 * 3$ 的实数组如 $（ 0.3 ， 0.2 ， 0.8 ， 0.1 ， 0.6 ， 0.4 ）$ 并对此这些元素按大小排序。规定数值小的顺序靠前且最小的数编号为0，则此时排序结果为：

[2，1，5，0，4，3]

将此数组每个元素分别除以最大工序数 m ，并取下整。则可得到一组可行序列 A ：

A =[0，0，1，0，1，1]

显然，可行序列 A 的加工路线为工件0的第1个工序，工件0的第2个工序，工件1的第1个工序，工件0的第3个工序，工件1的第2个工序，工件1的第3个工序。

如此，我就就通过一个可操作的实数序列转换成一个满足约束条件的工序。

但对于所有工件，其的操作序列可能最大工序次数小于m,我们只需将其正常加入，对于缺少的工序数量我们添加相同数量的且持续时间为0的无关工序即可，在进行计算工序时间时忽略这些持续时间为0的工序即可。

程序流程

算法流程基本同上述部分中描述基本一致，现描述如下：

读入数据，初始化程序，设定好工件数n，最大工件加工流程m，全局最优序列bestans,全局最优解ansCost，迭代次数step，种群数量Bugu.n。
随机Bugu.n组初始解，每组解当中存在 $n * m$ 个随机实数。
遍历初始解，得到当前最少工时数bestCost，与当前最优解下标bestIndex，更新bestans与ansCost。
若超过迭代次数，进入步骤10，否则继续向下执行。
对当前所有种群进行莱维飞行，所有解的解空间中的实数。注意对于bestIndex位置的种群不进行更新。
将25%数量的种群视为无效解，再其周围长度为0.4，方向为0-1均匀分布的正负方向上重新更新该点的值。
更新当前最优解下标bestIndex与当前最少工时数bestCost。
更新全局最优解ansCost与全局最优解。
回到步骤4。
输出结果。结束。

代码实现（by JAVA）

代码质量并不是特别好，请各位看官见谅。

package jspNew;import java.util.*;import org.apache.commons.math3.special.Gamma;/** 重点关注test类中的myRead方法中对于数据读入的部分。* 同时注意在Data类中reVec方法对于结果的返回。* 全局最优解 在主流程类Run中的run方法的ansCost中 全局最优机器任务序列，由reVec方法返回*/class test{public static void myRead() {Scanner cin = new Scanner(System.in);//全局定义Data.n = cin.nextInt(); //TODO 此处需要读取n个工件Data.m = cin.nextInt(); //TODO 此处需要读取最大工序数 每个工件最多有m个工序Pair[][] t = new Pair[Data.n][Data.m];			//暂存的任务列表Bugu.n = cin.nextInt();	//TODO 此处需要读取种群数量Bugu.step = cin.nextInt();	//TODO 此处需要读取迭代次数//t[0][0] = Pair.make_pair(1, 3);	//读取数据for(int i = 0 ; i < Data.n ; i++) {int tm = cin.nextInt(); //TODO 此处需要读取第i个工件的工序数for(int j = 0 ; j < Data.m ; j++) {if(j >= tm)t[i][j] = Pair.make_pair(0, 0);else {int a = cin.nextInt(),b = cin.nextInt(); //TODO a是第i个工件 在第j道工序 需要在 a机器上执行 b时间t[i][j] = Pair.make_pair(a, b);}}}//创建全局工件，并将数组赋值（从暂存的任务列表中获取）Data.init(Data.n);for(int i = 0 ; i < Data.n ; i++) {Data.myData[i] = Data.get_data(t[i]);}
//		int[] tans = new int[]{1,2,2,2,0,0,1,0,1};
//		System.out.println(Data.getCost(tans));}
}//实际上是任务类
class Pair{public int machine,time;Pair(int machine,int time){this.machine = machine;this.time = time;}//工厂模式static Pair make_pair(int machine,int time) {return new Pair(machine,time);}
}//2020年7月11日
class three extends Pair{int third;three(int a,int b,int c){super(a, b);third = c;}static three make_three(int a,int b,int c) {return new three(a,b,c);}
}class Data{public static Data[] myData;			//全局工件public static int n,m;					//全局工件数和最大任务数                                                                                                                                                                                                                                                            任务数量public static double Myrandom() {return Math.random();}public static void init(int nn) {n = nn;myData = new Data[nn];}public Pair[] task;						//每个工件的任务序列public static Data get_data(Pair[] t) {	//返回了一个工件Data tt= new Data();tt.task = new Pair[m];for(int i = 0 ; i < t.length ; i++) {tt.task[i] = t[i]; }return tt;}public static int getCost(int[] t) {Map<Integer,Integer> map = new HashMap<Integer,Integer>();	//记录每个工件做到了第几个步骤。int[] lastRun = new int[n];			//每个工件上个任务的完成时间。int[] lastMachine = new int[m];		//每个机器上个工件的完成时间。for(int i = 0; i < n ;i++) {map.put(i, 0);					//一开始都是第0个步骤。}int ans = -1;for(int i = 0 ; i < t.length ; i++) {int now = t[i];					//获取当前的工件编号int pairNum = map.get(now);		//获取工件编号的步骤。map.remove(now);				//移除工件上一个工件编号。map.put(now, pairNum+1);		//更新工件做到的步骤Pair nowPair = myData[now].task[pairNum];		//获取当前工件的步骤的任务。int nowMachine = nowPair.machine;				//获取任务机器。int nowCost = nowPair.time;						//获取任务时间。if(nowCost == 0)continue;int newTime = Math.max(lastMachine[nowMachine], lastRun[now]) + nowCost;	//从当前工件上一个任务结束时间//和待放入机器的最后一个工件完成时间中挑个大的。lastMachine[nowMachine] = lastRun[now] = newTime;
//			System.out.println("newTime:" + newTime
//					+" nowMachine:" + nowMachine + 
//					" nowCost:" + nowCost);ans = Math.max(ans, newTime);	//获取最长时间}return ans;}//2020年7月11日static Vector<three>[] reVec(int[] in){Vector<three>[] t = new Vector[n*m];int[] vis = new int[n];									//工序数int[] machineLastTime = new int[m];int[] taskLastTime = new int[n];for(int i = 0 ; i < t.length ; i++) {t[i] = new Vector<three>();}for(int i = 0 ; i < in.length ; i++) {int aim = in[i];int index = vis[aim]++;int machine = myData[aim].task[index].machine;//System.out.println(machine);int time = myData[aim].task[index].time;//System.out.println(time);if(time == 0)continue;int maxStartTime = taskLastTime[aim] > machineLastTime[machine] ? taskLastTime[aim] : machineLastTime[machine];/* System.out.println("aim:"+aim + "\n" + "index:" + index + "\n" + "machine:" + machine + "\n" +"time:" + time + "\n"+"taskLastTime:" +taskLastTime[aim] + "\n"+"machineLastTime" +machineLastTime[machine]+ "\n-------------------------"); */taskLastTime[aim] = machineLastTime[machine] = maxStartTime + time;t[machine].add(three.make_three(maxStartTime, taskLastTime[aim],aim));}//test//System.out.println(m);for(int i = 0 ; i < m ; i++) {System.out.print("machine "+i + ":" + '\n');for(int j = 0 ; j < t[i].size() ; j++)System.out.println("task :" + t[i].get(j).third + " start:" + t[i].get(j).machine + " end:" + t[i].get(j).time);}//TODO 返回了一个Vector数组，以 机器编号为 vector数组下标索引//对于每个Vector，中存放了若干个三元组，表示一个机器的 任务序列//对于每个任务序列中的每个任务 以 <工件名称,开始任务时间,结束任务时间>储存，//三元组的数据结构为three，可见代码中的three部分。return t;}
}class see{//解码，将编码解码成可行任务序列。public static int[] recode(double[] tt) {double[]t = new double[tt.length];	//临时数组，防止原数组被改。for(int i = 0 ; i < t.length ; i++)t[i] = tt[i];Map<Double,Integer> map = new HashMap<Double,Integer>();	//记录以下数组编号。for(int i = 0 ; i < t.length ; i++)map.put(t[i], i);Arrays.sort(t);					//排个序int[] ans = new int[t.length];	for(int i = 0 ; i < t.length ; i++) {ans[map.get(t[i])] = i;		//将排序后的顺序填入原来的位置}for(int i = 0 ; i < t.length ; i++) {int t1 = (int)((double)ans[i] / (double)Data.m);ans[i] = t1;}
//		for(int i = 0 ; i < ans.length ; i++)
//			System.out.println(ans[i]);
//		System.out.println();return ans;}//得到一个从u到v的正态分布public static double NormalDistribution(double u,double v){java.util.Random random = new java.util.Random();return Math.sqrt(v)*random.nextGaussian()+u;}//得到伽马函数public static double getGamma(double t) {return Gamma.gamma(t);}//获取u分布static double u,beta = 1.5;;private static void getU(){u = (getGamma(1.0 + beta)*Math.sin((Math.PI * beta) / 2.0) )/ (getGamma(((1.0 + beta) / 2.0) * beta * Math.pow(2.0, (beta - 1.0) / 2.0)));}//获取步长public static double getS() {getU();double newu = NormalDistribution(0,u);double newv = Math.pow(Math.abs(NormalDistribution(0,1)),1.0/beta);return newu/newv;}public static double randomWalk() {return  1 * 1 * getD();}//获取方向 https://blog.csdn.net/sj2050/article/details/98496868#_43public static double getD() {double r = (int)Data.Myrandom() * (100) + 1;if(r < 50) return -1.0;return 1.0;}static public void seeBest(double[] t,int bestCost) {System.out.println("本次布谷鸟找到的最优解是：" + bestCost);System.out.println("对应的序列是：");int[] ansShow = see.recode(t);for(int i= 0 ; i < ansShow.length ; i++) {System.out.printf("%3d ", ansShow[i]);}System.out.println();}static public void seeIndex(double[] t) {int[] ansShow = see.recode(t);for(int i= 0 ; i < ansShow.length ; i++) {System.out.printf("%3d ", ansShow[i]);}System.out.println();}
}class Bugu{static public int n;public double []ans;static public int step;Bugu(int n){ans = new double[n];for(int i = 0 ; i < n ; i++) {ans[i] = Data.Myrandom();}}static public void seeBugu(Bugu[] a) {for(int i = 0 ; i < a.length ; i++) {System.out.print("bugu["+i+"]:");for(int j = 0 ; j < a[i].ans.length ; j++) {System.out.print(a[i].ans[j] + " ");}System.out.println();}System.out.println();}}class run {public static void Run() {test.myRead();//Bugu.n = 4;int[] cost = new int[Bugu.n];Bugu[] bugu = new Bugu[Bugu.n];double[] bestans = new double[Data.m * Data.n];int ansCost = 0x3f3f3f;for(int i = 0 ; i < Bugu.n ; i++) {bugu[i] = new Bugu(Data.n * Data.m);//System.out.println();cost[i] = Data.getCost(see.recode(bugu[i].ans));//System.out.println(cost[i]);}int bestIndex = -1,bestCost = 0x3f3f3f;for(int i = 0 ; i < cost.length ; i++)if(cost[i] < bestCost) {bestCost = cost[i];bestIndex = i;for(int k = 0 ; k < bugu[i].ans.length; k++)bestans[k] = bugu[i].ans[k];}ansCost = bestCost;//see.seeBest(bestans, ansCost);//if(ansCost > -1)//	return;while(Bugu.step > 0) {//Bugu.seeBugu(bugu);Bugu.step--;//找个窝下蛋，更新窝的位置。即莱维飞行。for(int i = 0 ; i < Bugu.n ; i++) {//System.out.println(updata);for(int j = 0 ; j < bugu[i].ans.length ; j++) {double s = see.getS(),d = see.getD();	//s是步长，d是方向double updata = s * d;bugu[i].ans[j] += Math.abs(bestIndex - i)*(updata);//}	}//在差的解当中，25%的蛋被宿主发现for(int i = 0 ; i < Bugu.n ; i++ ) {if(Math.random() < 0.25) {for(int j = 0 ; j < bugu[i].ans.length ; j++) {double newUpdata = see.randomWalk();bugu[i].ans[j] += 0.4 * newUpdata;}}}//寻找本轮最佳bestCost = 0x3f3f3f;for(int i = 0 ; i < cost.length ; i++) {cost[i] = Data.getCost(see.recode(bugu[i].ans));//System.out.println("cost[i]:" + cost[i]);//see.seeIndex(bugu[i].ans);if(cost[i] < bestCost) {bestCost = cost[i];bestIndex = i;}}//本轮最佳比全局最佳好，则更新。if(bestCost < ansCost) {ansCost = bestCost;for(int i = 0 ; i < bestans.length ; i++)bestans[i] = bugu[bestIndex].ans[i];}//see.seeBest(bestans, ansCost);//Bugu.seeBugu(bugu);}//see.seeBest(bestans, ansCost);	//这一行只是用来打印结果的，按需删除//ansCost就是全局最优解的值，bestans就是最优解的序列。System.out.println("本次布谷鸟找到的最优解是：" + ansCost);Data.reVec(see.recode(bestans));}
}public class Jsp {public static void main(String[] argv) {//test.myRead();run.Run();//test.print();}
}
/*
3
3
4
1000
3
0 3
1 2
2 2
3
0 2
2 1
1 4
2
1 4
2 3
*/

算例结果

现以一个 $10 * 5$ 的算例进行测试：

实际上，在种群数量为5，迭代次数为20000次时，已经可以得到总时间为711的最优解，而在上述例子将迭代次数与种群数量大大提高的情况下，结果如下：

本次布谷鸟找到的最优解是：690
machine 0:
task :4 start:0 end:83
task :0 start:83 end:136
task :7 start:156 end:216
task :1 start:216 end:237
task :3 start:237 end:292
task :5 start:292 end:384
task :6 start:390 end:483
task :8 start:483 end:527
task :9 start:527 end:623
task :2 start:656 end:668
machine 1:
task :0 start:0 end:21
task :5 start:21 end:75
task :3 start:75 end:152
task :6 start:154 end:241
task :7 start:241 end:282
task :8 start:282 end:331
task :9 start:331 end:406
task :4 start:406 end:425
task :2 start:546 end:588
task :1 start:588 end:659
machine 2:
task :5 start:75 end:118
task :7 start:118 end:156
task :9 start:196 end:239
task :4 start:239 end:303
task :6 start:303 end:390
task :1 start:423 end:449
task :3 start:449 end:515
task :8 start:527 end:625
task :2 start:625 end:656
task :0 start:656 end:690
machine 3:
task :6 start:0 end:69
task :4 start:83 end:117
task :9 start:117 end:196
task :8 start:196 end:213
task :1 start:237 end:289
task :2 start:289 end:328
task :7 start:328 end:352
task :0 start:352 end:407
task :3 start:515 end:592
task :5 start:592 end:654
machine 4:
task :9 start:0 end:77
task :6 start:77 end:154
task :5 start:154 end:233
task :0 start:233 end:328
task :3 start:328 end:407
task :1 start:407 end:423
task :8 start:423 end:448
task :2 start:448 end:546
task :4 start:546 end:583
task :7 start:583 end:666

不难看出，大大提高了迭代次数确实使结果更加精确，但并没有特别大的提升。

实际上布谷鸟算法对算法中实际参数并不是特别敏感。