一、人工蜂群算法简介
1 蜜蜂采蜜
自然界中的蜜蜂总能在任何环境下以极高的效率找到优质蜜源，且能适应环境的改变。蜜蜂群的采蜜系统由蜜源、雇佣蜂、非雇佣蜂三部分组成，其中一个蜜源的优劣有很多要素，如蜜源花蜜量的大小、离蜂巢距离的远近、提取的难易程度等；雇佣蜂和特定的蜜源联系并将蜜源信息以一定概率形式告诉同伴；非雇佣蜂的职责是寻找待开采的蜜源，分为跟随蜂和侦查蜂两类，跟随峰是在蜂巢等待而侦查蜂是探测蜂巢周围的新蜜源。蜜蜂采蜜时，蜂巢中的一部分蜜蜂作为侦查蜂，不断并随机地在蜂巢附近寻找蜜源，如果发现了花蜜量超过某个阈值的蜜源，则此侦査蜂变为雇佣蜂开始釆蜜，采蜜完成后飞回蜂巢跳摇摆舞告知跟随峰。摇摆舞是蜜蜂之间交流信息的一种基本形式，它传达了有关蜂巢周围蜜源的重要信息如蜜源方向及离巢距离等，跟随峰利用这些信息准确评价蜂巢周围的蜜源质量。当雇佣蜂跳完摇摆舞之后，就与蜂巢中的一些跟随蜂一起返回原蜜源采蜜，跟随蜂数量取决于蜜源质量。以这种方式，蜂群能快速且有效地找到花蜜量最高的蜜源。
蜜蜂采蜜的群体智能就是通过不同角色之间的交流转换及协作来实现的。具体采蜜过程如图所示。在最初阶段，蜜蜂是以侦查蜂的形式出现，且对蜂巢周闱的蜜源没有任何了解，由于蜜蜂内在动机和外在的条件不同侦查蜂有两种选择：①成为雇佣蜂，开始在蜂巢周围随机搜索蜜源，如图中路线②成为跟随峰，在观察完摇摆舞后开始搜索蜜源，如图中路线。假设发现两个蜜源和，在发现蜜源后，该侦查蜂变成一只雇佣蜂，雇佣蜂利用其自身属性记住蜜源的位置，并立刻投入到采蜜中。采蜜完成后蜜蜂带着满载花蜜返回蜂巢，将花蜜卸载到卸蜜房，卸载完成后雇佣蜂有三种可能的行为①放弃自己发现的花蜜量不高的蜜源，变成一个不受约束的非雇佣蜂，如图中的路线；②在招募区跳摇摆舞，招募一些待在蜂巢中跟随峰，带领其再次返回所发现的蜜源如图中的路线；③不招募其他蜜蜂，继续回到原来的蜜源采蜜如图中的路线。在现实生活中并不是所有的蜜蜂一开始就立刻采蜜，另外大多数蜜蜂在一次采蜜完成后都会选择回到招募区跳摇摆舞来招募更多的蜜蜂去采蜜。
2 算法模型
人工蜂群算法就是模拟蜜蜂的采蜜过程而提出的一种新型智能优化算法，它也是由食物源、雇佣蜂和非雇佣蜂三部分组成。
食物源：食物源即为蜜源。在任何一个优化问题中，问题的可行解都是以一定形式给出的。在人工蜂群算法中，食物源就是待求优化问题的可行解，是人工蜂群算法中所要处理的基本对象。食物源的优劣即可行解的好坏是用蜜源花蜜量的大小即适应度来评价的。
雇佣蜂：雇佣蜂即为引领蜂与食物源的位置相对应，一个食物源对应一个引领蜂。在人工蜂群算法中，食物源的个数与引领蜂的个数相等；引领蜂的任务是发现食物源信息并以一定的概率与跟随蜂分享；概率的计算即为人工蜂群算法中的选择策略，一般是根据适应度值以轮盘赌的方法计算。
非雇佣蜂：非雇佣蜂包括跟随蜂和侦査蜂跟随蜂在蜂巢的招募区内根据引领蜂提供的蜜源信息来选择食物源，而侦查蜂是在蜂巢附近寻找新的食物源。在人工蜂群算法中，跟随蜂依据引领蜂传递的信息，在食物源附近搜索新食物源，并进行贪婪选择。若一个食物源在经过次后仍未被更新，则此引领蜂变成侦査蜂，侦查蜂寻找新的食物源代替原来的食物源。
3 算法搜索过程
4 分类
人工蜂群算法中将人工蜂群分为引领蜂、跟随蜂和侦查蜂三类，每一次搜索过程中，引领蜂和跟随蜂是先后开采食物源，即寻找最优解，而侦查蜂是观察是否陷入局部最优，若陷入局部最优则随机地搜索其它可能的食物源。每个食物源代表问题一个可能解，食物源的花蜜量对应相应解的质量（适应度值fiti）。
（1）人工蜂群算法搜索过程中，首先需要初始化，其中包括确定种群数、最大迭代次数MCN、、控制参数limit和确定搜索空间即解的范围，在搜索空间中随机生成初始解xi(i=1,2,3,……,SN)，SN为食物源个数，每个解xi是一个D维的向量，D是问题的维数。初始化之后，整个种群将进行引领蜂、跟随蜂和侦查蜂搜寻过程的重复循环，直到达到最大迭代次数MCN或误差允许值 ε。
（2）在搜索过程开始阶段，每个引领蜂由式(2-3)产生一个新解即新食物源，
vij=xij+Φij(xij-xkj) (2-3)
式中，k∈﹛1,2,...,SN﹜,j∈{1,2,...,D},且k ≠i;Φij为[-1,1]之间的随机数。计算新解的fiti并评价它，若新解的fiti优于旧解，则引领蜂记住新解忘记旧解。反之，它将保留旧解。
（3）在所有引领蜂完成搜寻过程之后，引领蜂会在招募区跳摇摆舞把解的信息及信息与跟随蜂分享。跟随蜂根据式计算每个解的选择概率，
pi=fiti/∑k=1SNfitk。 (2-4)
然后在区间[-1,1]内随机产生一个数，如果解的概率值大于该随机数，则跟随蜂由式(2-3)产生一个新解，并检验新解的fiti，若新解的fiti比之前好，则跟随蜂将记住新解忘掉旧解；反之，它将保留旧解。
四、在所有跟随蜂完成搜寻过程之后，如果一个解经过limit次循环仍然没有被进一步更新，那么就认为此解陷入局部最优，该食物源就会被舍弃。设食物源xi被舍弃，则此食物源对应的引领蜂转成一个侦查蜂。侦察蜂由(2-5)式产生一个新的食物源代替它。
xij=xminj+rand(0,1)(xmaxj-xminj) (2-5)
其中j∈{1,2....,D}。然后返回引领蜂搜索过程，开始重复循环。
五、人工蜂群算法的食物源质量一般是越大越好，即适应度值越大越好，而对应于要优化的问题，需要分两种情况考虑：即最小值问题、最大值问题。设fi是优化问题的目标函数，所以若优化最小值问题时，适应度函数为fi的变形，一般用式(2-6)表示；若优化最大值问题，适应度函数即目标函数。
fiti={1+abs(fi) fi>=01/1+fi fi>0 (2-6)
人工蜂群算法在评价食物源时一般进行贪婪选择按式(2-7)进行。
vi={xi fit(xi)<=fit(vi)vi fit(vi)>fit(xi) (2-7)
人工蜂群算法就是通过循环搜索，最终找到最优食物源或最优解。
5 算法步骤
人工蜂群算法具体实现步骤：
步骤1：初始化种群：初始化各个参数，蜂群总数SN、食物源被采集次数即最大迭代次数MCN及控制参数limit，确定问题搜索范围，并且在搜索范围内随机产生初始解xi(i=1,2,...SN) 。
步骤2：计算并评估每个初始解的适应度。
步骤3：设定循环条件并开始循环
步骤4：引领蜂对解xi按照式（2-3）进行邻域搜索产生新解（食物源）vi，并计算其适应度值；
步骤5：按照式(2-7)进行贪婪选择：如果vi的适应度值优于xi，则利用vi替换xi，将vi作为当前最好的解，否则保留xi不变；
步骤6：根据式(2-4)计算食物源的概率pi；
步骤7：跟随蜂依照概率pi选择解或食物源，按照式(2-3)搜索产生新解（食物源）vi，并计算其适应度。
步骤8：按式(2-7)进行贪婪选择；如果vi的适应度优于xi，则用vi代替xi，将vi作为当前最好解，否则保留xi不变；
步骤9：判断是否有要放弃的解。若有，则侦查蜂按式(2-5)随机产生新解将其替换；
步骤10：记录到目前为止的最优解；
步骤11：判断是否满足循环终止条件，若满足，循环结束，输出最优解，否则返回步骤4继续搜索。
二、部分源代码
%%%%%%%%%%%%
%主程序开始
clc;
clear;
TestProblem='UF7';
nVar=10;
fobj = cec09(TestProblem);
xrange = xboundary(TestProblem, nVar);
% Lower bound and upper bound
lb=xrange(:,1)';
ub=xrange(:,2)';
FoodNumber = 100;
Limit = 100;
maxCycle = 200;
runtimes=1;
GD=zeros(runtimes,1);
IGD=zeros(runtimes,1);
Spread=zeros(runtimes,1);
TPF=xlsread('Case3.xls',TestProblem);
%种群初始化
M=2; %目标函数值个数 %函数维度
%蜂群初始化
Range = repmat((ub - lb),[FoodNumber 1]);
Lower = repmat(lb,[FoodNumber 1]);
for r=1:runtimes
Foods = rand(FoodNumber,nVar).*Range+Lower;
Funcvalue = zeros(FoodNumber,M); %记录种群中所有个体的目标函数值
trial = zeros(1,FoodNumber);
%%%%%计算所有个体的所有目标函数值
for i=1:FoodNumber
Funcvalue(i,:)=fobj(Foods(i,:)')'; %计算个体所有目标函数的函数值 %Funcvalue 中每一行代表一个个体的所有目标函数值
end
%%%%%%%%%%%%%外部档案集初始化 %%%%%%%%%%%%%%
pop=[Foods,Funcvalue]; %将种群个体和其对应的目标函数值组合成一个矩阵 %每一行前6列为每维自变量值，7和8列分别为对应的两个目标函数的值
Rpop=quick_sort(pop);%% 外部档案Rpop
%快速排序，传入参数 %全部种群，种群中所有个体的目标函数值，以及当前优化的问题序号 %非支配集初始化
Npop=Rpop; %外部档案和非支配解集中为所有非支配解及其对应的目标函数值
iter=1;
j=1;
while((iter<=maxCycle))
%%%%%%%%%%%%%采蜜蜂----邻域搜索，寻找局部最优解
for i=1:(FoodNumber)
%在外部档案中随机选取一个解，对当前的邻域搜索进行指导
[P,Q]=size(Rpop);
s=1+fix(rand*P);
Param2Change=fix(rand*nVar)+1; %随机选取一维
sol=Foods(i,:);
sol(Param2Change)=Foods(i,Param2Change)+(Rpop(s,Param2Change)-Foods(i,Param2Change))*(rand-0.5)*2; %对当前解的随机一维进行领域搜索指导
%使用外部档案中的随机解对当前解进行更新
ind=find(sol<lb);
sol(ind)=lb(ind);
ind=find(sol>ub);
sol(ind)=ub(ind);
FuncSol=fobj(sol')'; %计算邻域解的目标函数值
%%%%%%%%%以pareto支配为标准的新贪婪准则 %%%%%%%%%%
if Dominates(FuncSol',Funcvalue(i,:)')
Foods(i,:)=sol;
Funcvalue(i,:)=FuncSol;
trial(i)=0;
else
trial(i)=trial(i)+1;
end
end
prob=calprob(Foods,Funcvalue); %计算跟随概率
%%%%%%%%%%%%%%% 观察蜂模式，根据概率随机选择蜜源 %%%%%%%%%%%
i=1;
t=0;
while(t<FoodNumber) %循环完后用完FoodNumber个观察蜂
if(rand<prob(i)) %%%如果成立，才往下执行，用去一只跟随蜂
t=t+1;
sol(Param2Change)=Foods(i,Param2Change)+(Rpop(s,Param2Change)-Foods(i,Param2Change))*(rand-0.5)*2;%使用外部档案中的随机解对当前解进行更新
ind=find(sol<lb);
sol(ind)=lb(ind);
ind=find(sol>ub);
sol(ind)=ub(ind);
FuncSol=fobj(sol')';
%%%%%%%%% 以pareto支配为标准的新贪婪准则 %%%%%%%%
if Dominates(FuncSol',Funcvalue(i,:)')
Foods(i,:)=sol;
Funcvalue(i,:)=FuncSol;
trial(i)=0;
else
trial(i)=trial(i)+1;
end
end
i=i+1;
if(i==(FoodNumber)+1)
i=1;
end
end
%%%%%%% 更新外部档案 %%%%%%%%
pop=[Foods,Funcvalue]; %将种群和目标函数值写入同一矩阵
Npop=quick_sort(pop); %将当前种群的所有非支配解写入非支配解集
temp_set=[Rpop;Npop]; %将Npop加入外部档案，构成新的临时外部档案temp_set
Rpop=quick_sort(temp_set); %新的外部档案集
while(size(Rpop,1))>100 %根据拥挤距离对外部档案进行裁剪，使其中的解不多于100个
crowd_dictance=[];
crowd_dictance=crowd_distance_measure(Rpop);
min_distance_num=1; %找外部档案中拥挤距离最小的个体
for i=2:size(crowd_dictance,1)
if crowd_dictance(i)<crowd_dictance(min_distance_num)
min_distance_num=i;
end
end
for i=min_distance_num:(size(Rpop,1)-1)
Rpop(i,:)=Rpop(i+1,:); %把拥挤距离最小的支配解删除，同时把最后一个支配解移动一次。删除一个拥挤距离最小的支配解且新生成的支配解集仍然连续。
end %从拥挤距离最小的非支配解的位置开始 左移，后一个非支配解覆盖前一个非支配解。
Rpop=Rpop(1:size(Rpop,1)-1,:); %size(Rpop,1)-1减一是为了把最后一个非支配解给去除，因为它已经被移动到倒数第二个位置上！
end
%%%%%%%%%%%%%绘制优化结果 %%%%%%%%%%%%%%
plot(Rpop(:,nVar+1),Rpop(:,nVar+2),'ro');
GD(r)=GD_matlab(Rpop(:,nVar+1:nVar+2),TPF);
IGD(r)=IGD_matlab(Rpop(:,nVar+1:nVar+2),TPF);
Spread(r)=Spread_matlab(Rpop(:,nVar+1:nVar+2),TPF);
end
%%%%%%%%%主程序结束 %%%%%%%%%%
%%%%%%%%%%%寻找支配解%%%%%%%%
function non_dom=quick_sort(pop)
npop=size(pop,1);
K=10;
pop(:,K+3)=0;
for i=1:npop
pop(i,K+3)=0;
for j=1:i-1
if pop(j,K+3)==0
if Dominates(pop(i,K+1:K+2)',pop(j,K+1:K+2)')
pop(j,K+3)=1;
elseif Dominates(pop(j,K+1:K+2)',pop(i,K+1:K+2)')
pop(i,K+3)=1;
break;
end
end
end
end
ind=find(pop(:,K+3)==0);
non_dom=pop(ind,1:K+2);
% [N,P]=size(pop);
% M=2; %目标函数值的个数
% V=10; %每个目标函数值的维数
% non_dom=[];
% individual=pop; %pop：种群（自变量+目标函数值），行数为种群数量
% while ~isempty(individual) %individual为空时停止循环
% temp1=[];
% temp2=[];
% flag=0;
% %Number of individuals that dominate this individual 支配i的解的个数
% %individual(i).n=0;
% %Individuals which this individual dominate 被i支配的解
% %individual(i).p=[]
% for i=2:size(individual,1) %size(individual,1)为individual的行数
% dom_less=0;
% dom_equal=0;
% dom_more=0;
% for j=1:M %M=2，两个目标函数
% if(individual(1,V+j)<individual(i,V+j)) %V=6，每个目标函数均为六维 %%pareto支配判断
% dom_less=dom_less+1; %解1有可能支配解i
% elseif (individual(1,V+j)==individual(i,V+j))
% dom_equal=dom_equal+1; %互不支配
% else
% dom_more=dom_more+1; %解1有可能被解i支配
% end
% end
% if dom_less==0 && dom_equal~=M %该行为真时，表示解i支配解1
% flag=flag+1; %flag中存放解i支配解1的个数
% temp1=[temp1;individual(i,:)]; %temp1中存放支配解i
% elseif dom_more==0 && dom_equal~=M %该行为真时，表示解1支配解i
% temp2=[temp2;individual(1,:)]; %temp2中存放支配解1
% else %其余情况：两个解互不支配
% temp1=[temp1;individual(i,:)]; %temp1中也存放有可能支配的解i
% end %则称j优于i，则把此时的j放入individual(i)_p中
%此时temp1中存放的就是通过上次比较后剩下的非支配解（非支配解即不受其他解支配）
end
function out=calvalue(x,pro) %计算目标函数值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if pro==1
f=[];
f(1)=x(1); %目标函数1
sum=0;
for i= 2:6
sum=sum+x(i)/5;
end
g_x=1+9*sum; %中间函数
f(2)=g_x*(1-(f(1)/(g_x))^0.5);%目标函数值2
out=f;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif pro==2
f=[];
f(1)=x(1); %目标函数1
sum=0;
for i= 2:6
sum=sum+x(i)/5;
end
g_x=1+9*sum; %中间函数
f(2)=g_x*(1-(f(1)/(g_x))^2);%目标函数值2
out=f;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif pro==3
f=[];
f(1)=x(1); %目标函数1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif pro==6
f=[];
f(1)=1-exp((-4)*x(1))*(sin(6*pi*x(1)))^6; %目标函数1
sum=0;
for i= 2:6
sum=sum+x(i)/5;
end
g_x=1+9*sum^0.25; %中间函数
f(2)=g_x*(1-(f(1)/(g_x))^2);%目标函数值2
out=f;
end
end
三、运行结果
四、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
