【优化算法】蝴蝶优化算法（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.

tic

if ~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 : Keep
        chromKeep(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.popsize
        if popindex <= numButterfly1
            Population1(popindex).chrom = Population(popindex).chrom;
        else
            Population2(popindex-numButterfly1).chrom = Population(popindex).chrom;
        end
    end
    % % % % % % % % % % %    End of Divide the whole population into two subpopulations  % % %%%
    %%
    
    % % % % % % % % % % % %%            Migration operator          % % % % % % % % % % % %%%%
    %% Migration operator
    for k1 = 1 : numButterfly1
        for parnum1 = 1 : OPTIONS.numVar
            r1 = rand*period;
            if r1 <= partition
                r2 = round(numButterfly1 * rand + 0.5);
                Land1(k1,parnum1) = Population1(r2).chrom(parnum1);
            else
                r3 = round(numButterfly2 * rand + 0.5);
                Land1(k1,parnum1) = Population2(r3).chrom(parnum1);
            end
        end %% for parnum1
        NewPopulation1(k1).chrom =  Land1(k1,:);
    end  %% for k1
    % % % % % % % % % % % %%%       End of Migration operator      % % % % % % % % % % % %%%
    %%
    
    % % % % % % % % % % % %             Evaluate NewPopulation1          % % % % % % % % % % % %%
    %% Evaluate NewPopulation1
    SavePopSize = OPTIONS.popsize;
    OPTIONS.popsize = numButterfly1;
    % Make sure each individual is legal.
    NewPopulation1 = FeasibleFunction(OPTIONS, NewPopulation1);
    % Calculate cost
    NewPopulation1 = CostFunction(OPTIONS, NewPopulation1);
    OPTIONS.popsize = SavePopSize;
    % % % % % % % % % % % %       End of Evaluate NewPopulation1      % % % % % % % % % % % %%
    %%
    
    % % % % % % % % % % % %             Butterfly adjusting operator          % % % % % % % % % % % %%
    %% Butterfly adjusting operator
    for k2 = 1 : numButterfly2
        scale = maxStepSize/(GenIndex^2); %Smaller step for local walk
        StepSzie = 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);
            else
                r4 = round(numButterfly2*rand + 0.5);
                Land2(k2,parnum2) =  Population2(r4).chrom(1);
                if (rand > BAR) % Butterfly-Adjusting rate
                    Land2(k2,parnum2) =  Land2(k2,parnum2) +  scale*(delataX(parnum2)-0.5);
                end
            end
        end  %% for parnum2
        NewPopulation2(k2).chrom =  Land2(k2,:);
    end %% for k2
    % % % % % % % % % % % %       End of Butterfly adjusting operator      % % % % % % % % % % % %
    %%
    
    % % % % % % % % % % % %             Evaluate NewPopulation2          % % % % % % % % % % % %%
    %% Evaluate NewPopulation2
    SavePopSize = OPTIONS.popsize;
    OPTIONS.popsize = numButterfly2;
    % Make sure each individual is legal.
    NewPopulation2 = FeasibleFunction(OPTIONS, NewPopulation2);
    % Calculate cost
    NewPopulation2 = 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 Population
    Population = CombinePopulation(OPTIONS, NewPopulation1, NewPopulation2);
    % Sort from best to worst
    Population = 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 : Keep
        Population(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 worst
    Population = PopSort(Population);
    % Compute the average cost
    [AverageCost, nLegal] = ComputeAveCost(Population);
    % Display info to screen
    MinCost = [MinCost Population(1).cost];
    AvgCost = [AvgCost AverageCost];
    if DisplayFlag
        disp(['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 distribution
    fx = tan(pi * rand(1,StepSize));
    delataX(i) = sum(fx);
end


  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
  
24
  
25
  
26
  
27
  
28
  
29
  
30
  
31
  
32
  
33
  
34
  
35
  
36
  
37
  
38
  
39
  
40
  
41
  
42
  
43
  
44
  
45
  
46
  
47
  
48
  
49
  
50
  
51
  
52
  
53
  
54
  
55
  
56
  
57
  
58
  
59
  
60
  
61
  
62
  
63
  
64
  
65
  
66
  
67
  
68
  
69
  
70
  
71
  
72
  
73
  
74
  
75
  
76
  
77
  
78
  
79
  
80
  
81
  
82
  
83
  
84
  
85
  
86
  
87
  
88
  
89
  
90
  
91
  
92
  
93
  
94
  
95
  
96
  
97
  
98
  
99
  
100
  
101
  
102
  
103
  
104
  
105
  
106
  
107
  
108
  
109
  
110
  
111
  
112
  
113
  
114
  
115
  
116
  
117
  
118
  
119
  
120
  
121
  
122
  
123
  
124
  
125
  
126
  
127
  
128
  
129
  
130
  
131
  
132
  
133
  
134
  
135
  
136
  
137
  
138
  
139
  
140
  
141
  
142
  
143
  
144
  
145
  
146
  
147
  
148
  
149
  
150
  
151
  
152
  
153
  
154
  
155
  
156
  
157
  
158
  
159
  
160
  
161
  
162
  
163
  
164
  
165
  
166
  
167
  
168
  
169
  
170
  
171
  
172
  
173
  
174
  
175
  
176
  
177
  
178
  
179
  
180
  
181
  
182
  
183
  
184
  
185
  
186
  
187
  
188
  
189
  
190
  
191
  
192
  
193
 

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

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

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/117442075