【优化算法】蝴蝶优化算法（MBO）【含Matlab源码 952期】

一、获取代码方式

获取代码方式1：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

获取代码方式2：
通过紫极神光博客主页开通CSDN会员，凭支付凭证，私信博主，可获得此代码。

获取代码方式3：
完整代码已上传我的资源：【优化算法】蝴蝶优化算法（MBO）【含Matlab源码 952期】

备注：开通CSDN会员，仅只能免费获得1份代码（有效期为开通日起，三天内有效）；
订阅紫极神光博客付费专栏，可免费获得2份代码（有效期为订阅日起，三天内有效）；

二、蝴蝶优化算法（MBO）简介

1 介绍
蝴蝶优化算法(butterfly optimization algorithm, BOA)是Arora 等人于2019年提出的一种元启发式智能算法。该算法受到了蝴蝶觅食和交配行为的启发，蝴蝶接收/感知并分析空气中的气味，以确定食物来源/交配伙伴的潜在方向。
蝴蝶利用它们的嗅觉、视觉、味觉、触觉和听觉来寻找食物和伴侣，这些感觉也有助于它们从一个地方迁徙到另一个地方，逃离捕食者并在合适的地方产卵。在所有感觉中，嗅觉是最重要的，它帮助蝴蝶寻找食物（通常是花蜜）。蝴蝶的嗅觉感受器分散在蝴蝶的身体部位，如触角、腿、触须等。这些感受器实际上是蝴蝶体表的神经细胞，被称为化学感受器。它引导蝴蝶寻找最佳的交配对象，以延续强大的遗传基因。雄性蝴蝶能够通过信息素识别雌性蝴蝶，信息素是雌性蝴蝶发出的气味分泌物，会引起特定的反应。
通过观察，发现蝴蝶对这些来源的位置有非常准确的判断。此外，它们可以辨识出不同的香味，并感知它们的强度。蝴蝶会产生与其适应度相关的某种强度的香味，即当蝴蝶从一个位置移动到另一个位置时，它的适应度会相应地变化。当蝴蝶感觉到另一只蝴蝶在这个区域散发出更多的香味时，就会去靠近，这个阶段被称为全局搜索。另外一种情况，当蝴蝶不能感知大于它自己的香味时，它会随机移动，这个阶段称为局部搜索。

2 香味
为了理解BOA中的香味是如何计算的，首先需要理解，像气味、声音、光、温度等这样的模态是如何计算的。感知、处理这些模态需要知道三个重要的术语：感觉模态C、刺激强度I和幂指数a。在感觉模态中，感觉意味着测量能量的形式并以类似方式对其进行处理，而模态是指传感器使用的原始输入。不同的形态可以是气味，声音，光线，温度，在BOA中，模态是香味。I是物理刺激的大小。在BOA中，I与蝴蝶/解决方案的适应度相关。这意味着，当一只蝴蝶散发出更多的香味时，周围的其他蝴蝶可以感知到并被吸引。幂是强度增加的指数。参数a允许正则表达式、线性响应和响应压缩。响应扩展是当I增加时，香味（f）比I增长更快。响应压缩是当I增加时，f比I增长慢。线性响应是当I增加时，f成比例地增加。经实验证明，有时随着刺激的增强，昆虫对刺激变化的敏感性变得越来越低。因此在BOA中，为了估计I的大小，使用了响应压缩。
蝴蝶的自然现象基于两个重要问题：I的变化和f的表示。简单地说，蝴蝶的I与编码后的目标函数相关联。但是，f是相对的，即应该由其他蝴蝶来感知。史蒂文斯幂定律中，为了将气味与其他形式区别开来，使用了C。现在，当I较少的蝴蝶向I较多的蝴蝶移动时，f比I增加得更快。因此，我们应该允许f随幂指数参数a实现的吸收程度而变化。在BOA中，香味被表示为刺激物的物理强度的函数，如下所示：

3 具体算法
为了用搜索算法演示上述讨论，将蝴蝶的上述特征理想化如下：
（1）所有的蝴蝶都可以发出气味，这使蝴蝶间相互吸引。
（2）每只蝴蝶都会随机移动或朝最好的蝴蝶移动，散发出更多的芳香。
（3）蝴蝶的刺激强度受目标函数的景观影响或决定。
该算法分为三个阶段：(1)初始化阶段、(2)迭代阶段和(3)结束阶段。
在BOA的每次运行中，首先执行初始化阶段，然后进行迭代搜索，最后在找到最优解时终止算法。BOA中使用的参数值也会被分配，设置这些值后，算法将继续创建初始蝴蝶种群以进行优化。由于在BOA的模拟过程中蝴蝶总数保持不变，分配了一个固定大小的内存来存储信息。蝴蝶的位置是在搜索空间中随机生成的，并计算和存储它们的香味和适应值。这样就完成了初始化阶段，算法开始了迭代阶段，该阶段使用创建的人工蝶形执行搜索。算法的第二阶段，即迭代阶段，由算法执行多次迭代。在每次迭代中，解空间中的所有蝶形都移到新位置，然后重新评估其适应性值。算法首先计算解空间中不同位置的所有蝴蝶的适应度值。那么这些蝴蝶就会利用式1在自己的位置产生香味。该算法有两个关键步骤，即全局搜索阶段和局部搜索阶段。在全局搜索阶段，蝴蝶向最合适的蝴蝶/解g∗迈出一步，该蝴蝶/解g可以用公式(2)来表示。

这里，g∗表示在当前迭代的所有解中找到的当前最佳解；fi表示第i只蝴蝶的香味，r是[0，1]中的随机数。局部搜索阶段可以表示为

其中，xjt和xkt是解空间中的第j个蝴蝶和第k个蝴蝶。
蝴蝶寻找食物、交配伙伴可以在局部和全局范围内发生。考虑到地理上的接近和各种其他因素，如雨、风等，在整个交配伙伴或蝴蝶的觅食活动中，寻找食物可能占很大比例。因此，在BOA中使用切换概率p来在普通全局搜索和密集局部搜索之间切换。
在未达到停止标准之前，一直进行迭代。迭代结束的标准可以有多个，如使用的最大CPU时间、达到的最大迭代次数、没有改进的最大迭代次数、达到错误率的特定值或任何其他适当的标准。当迭代阶段结束时，算法输出具有最佳适应度的最优解。

三、部分源代码

function [MinCost] = MBO(ProblemFunction, DisplayFlag, RandSeed)% Monarch Butterfly Optimization (MBO) software for minimizing a general function
% The fixed generation is considered as termination condition.% INPUTS: ProblemFunction is the handle of the function that returns
%         the handles of the initialization, cost, and feasibility functions.
%         DisplayFlag = true or false, whether or not to display and plot results.
%         ProbFlag = true or false, whether or not to use probabilities to update emigration rates.
%         RandSeed = random number seed
% OUTPUTS: MinCost = array of best solution, one element for each generation
%          Hamming = final Hamming distance between solutions
% CAVEAT: The "ClearDups" function that is called below replaces duplicates with randomly-generated
%         individuals, but it does not then recalculate the cost of the replaced individuals.ticif ~exist('ProblemFunction', 'var')ProblemFunction = @Ackley;
end
if ~exist('DisplayFlag', 'var')DisplayFlag = true;
end
if ~exist('RandSeed', 'var')RandSeed = round(sum(100*clock));
end[OPTIONS, MinCost, AvgCost, InitFunction, CostFunction, FeasibleFunction, ...MaxParValue, MinParValue, Population] = Init(DisplayFlag, ProblemFunction, RandSeed);% % % % % % % % % % % %             Initial parameter setting          % % % % % % % % % % % %%%%
%% Initial parameter setting
Keep = 2; % elitism parameter: how many of the best habitats to keep from one generation to the next
maxStepSize = 1.0;        %Max Step size
partition = OPTIONS.partition;
numButterfly1 = ceil(partition*OPTIONS.popsize);  % NP1 in paper
numButterfly2 = OPTIONS.popsize - numButterfly1; % NP2 in paper
period = 1.2; % 12 months in a year
Land1 = zeros(numButterfly1, OPTIONS.numVar);
Land2 = zeros(numButterfly2, OPTIONS.numVar);
BAR = partition; % you can change the BAR value in order to get much better performance
% % % % % % % % % % % %       End of Initial parameter setting       % % % % % % % % % % % %%
%%% % % % % % % % % % % %             Begin the optimization loop        % % % % % % % % % %%%%
% Begin the optimization loop
for GenIndex = 1 : OPTIONS.Maxgen% % % % % % % % % % % %            Elitism Strategy           % % % % % % % % % % % %%%%%%% Save the best monarch butterflis in a temporary array.for j = 1 : KeepchromKeep(j,:) = Population(j).chrom;costKeep(j) = Population(j).cost;end% % % % % % % % % % % %       End of  Elitism Strategy      % % % % % % % % % % % %%%%%%% % % % % % % % % % % %    Divide the whole population into two subpopulations % % % %%%%% Divide the whole population into Population1 (Land1) and Population2 (Land2)% according to their fitness.% The monarch butterflies in Population1 are better than or equal to Population2.% Of course, we can randomly divide the whole population into Population1 and Population2.% We do not test the different performance between two ways.for popindex = 1 : OPTIONS.popsizeif popindex <= numButterfly1Population1(popindex).chrom = Population(popindex).chrom;elsePopulation2(popindex-numButterfly1).chrom = Population(popindex).chrom;endend% % % % % % % % % % %    End of Divide the whole population into two subpopulations  % % %%%%%% % % % % % % % % % % %%            Migration operator          % % % % % % % % % % % %%%%%% Migration operatorfor k1 = 1 : numButterfly1for parnum1 = 1 : OPTIONS.numVarr1 = rand*period;if r1 <= partitionr2 = round(numButterfly1 * rand + 0.5);Land1(k1,parnum1) = Population1(r2).chrom(parnum1);elser3 = round(numButterfly2 * rand + 0.5);Land1(k1,parnum1) = Population2(r3).chrom(parnum1);endend %% for parnum1NewPopulation1(k1).chrom =  Land1(k1,:);end  %% for k1% % % % % % % % % % % %%%       End of Migration operator      % % % % % % % % % % % %%%%%% % % % % % % % % % % %             Evaluate NewPopulation1          % % % % % % % % % % % %%%% Evaluate NewPopulation1SavePopSize = OPTIONS.popsize;OPTIONS.popsize = numButterfly1;% Make sure each individual is legal.NewPopulation1 = FeasibleFunction(OPTIONS, NewPopulation1);% Calculate costNewPopulation1 = CostFunction(OPTIONS, NewPopulation1);OPTIONS.popsize = SavePopSize;% % % % % % % % % % % %       End of Evaluate NewPopulation1      % % % % % % % % % % % %%%%% % % % % % % % % % % %             Butterfly adjusting operator          % % % % % % % % % % % %%%% Butterfly adjusting operatorfor k2 = 1 : numButterfly2scale = maxStepSize/(GenIndex^2); %Smaller step for local walkStepSzie = ceil(exprnd(2*OPTIONS.Maxgen,1,1));delataX = LevyFlight(StepSzie,OPTIONS.numVar);for parnum2 = 1:OPTIONS.numVar,if (rand >= partition)Land2(k2,parnum2) = Population(1).chrom(parnum2);elser4 = round(numButterfly2*rand + 0.5);Land2(k2,parnum2) =  Population2(r4).chrom(1);if (rand > BAR) % Butterfly-Adjusting rateLand2(k2,parnum2) =  Land2(k2,parnum2) +  scale*(delataX(parnum2)-0.5);endendend  %% for parnum2NewPopulation2(k2).chrom =  Land2(k2,:);end %% for k2% % % % % % % % % % % %       End of Butterfly adjusting operator      % % % % % % % % % % % %%%% % % % % % % % % % % %             Evaluate NewPopulation2          % % % % % % % % % % % %%%% Evaluate NewPopulation2SavePopSize = OPTIONS.popsize;OPTIONS.popsize = numButterfly2;% Make sure each individual is legal.NewPopulation2 = FeasibleFunction(OPTIONS, NewPopulation2);% Calculate costNewPopulation2 = CostFunction(OPTIONS, NewPopulation2);OPTIONS.popsize = SavePopSize;% % % % % % % % % % % %       End of Evaluate NewPopulation2      % % % % % % % % % % % %%%%% % % % % % %  Combine two subpopulations into one and rank monarch butterflis       % % % % % %%% Combine Population1 with Population2 to generate a new PopulationPopulation = CombinePopulation(OPTIONS, NewPopulation1, NewPopulation2);% Sort from best to worstPopulation = PopSort(Population);% % % % % %     End of Combine two subpopulations into one and rank monarch butterflis  % %% % %%%% % % % % % % % % % % %            Elitism Strategy          % % % % % % % % % % % %%% %% %%% Replace the worst with the previous generation's elites.n = length(Population);for k3 = 1 : KeepPopulation(n-k3+1).chrom = chromKeep(k3,:);Population(n-k3+1).cost = costKeep(k3);end % end for k3% % % % % % % % % % % %     End of  Elitism Strategy      % % % % % % % % % % % %%% %% %%%% % % % % % % % % %           Precess and output the results          % % % % % % % % % % % %%%% Sort from best to worstPopulation = PopSort(Population);% Compute the average cost[AverageCost, nLegal] = ComputeAveCost(Population);% Display info to screenMinCost = [MinCost Population(1).cost];AvgCost = [AvgCost AverageCost];if DisplayFlagdisp(['The best and mean of Generation # ', num2str(GenIndex), ' are ',...num2str(MinCost(end)), ' and ', num2str(AvgCost(end))]);end% % % % % % % % % % %    End of Precess and output the results     %%%%%%%%%% %% %%%end % end for GenIndex
Conclude1(DisplayFlag, OPTIONS, Population, nLegal, MinCost, AvgCost);toc% % % % % % % % % %     End of Monarch Butterfly Optimization implementation     %%%% %% %
%%function [delataX] = LevyFlight(StepSize, Dim)%Allocate matrix for solutions
delataX = zeros(1,Dim);%Loop over each dimension
for i=1:Dim% Cauchy distributionfx = tan(pi * rand(1,StepSize));delataX(i) = sum(fx);
end

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.
[3]蝴蝶优化算法