蜜蜂算法（BA）

        中文互联网上面好像没有关于蜜蜂算法的详细介绍，正好最近在使用蜜蜂算法，在次介绍一下。蜜蜂算法（Bee Algorithm, BA）是由Pham教授于2005年提出的一种优化算法。该算法也属于自然启发式算法的一种，这里要注意不能将其与人工蜂群算法（Artificial Bee Colony Algorithm, ABC）弄混淆了。

蜜蜂算法的工作流如图所示，

1、首先初始化侦察蜂群（n个）

2、计算该蜂群中蜜蜂的fitness适应值再排序

3、按照fitness值的大小选出最优（最大或最小取决于具体问题）蜜蜂（e个），将此蜜蜂的采集点标定为精英采集点

4、按照fitness值的大小选出次采集点（m-e个），所以总共的站点有m个

5、决定领域搜索范围，这是组合优化中的概念。分为离散问题和连续问题，要具体问题具体分析。

6、派出雇佣蜂到选出的站点去进行领域搜索，精英采集点的雇佣蜂可以比次优采集点的雇佣蜂多一些。

7、如果雇佣蜂在领域搜索的时候找到了更好的采集点（拥有更好的fitness值），那么就取代掉原来的蜜蜂，将自己的位置作为新的采集点。

8、将派出除了采集点之外的n-m个蜜蜂去重新搜索，并且计算fitness值，生成新的蜂群（n个），回到第二步。跳出的循环的条件可以是迭代次数上限也可以是最优fitness值达到要求。

        下图是连续问题下的实例图。在该例子中，生成蜜蜂的办法是生成区间内随机大小的x，得到的fitness值就是他在函数上的映射y，取y最大的为采集点，派出侦察蜂。。。

        离散问题就拿旅行商（TSP）问题举例，生成蜜蜂就是代表途径城市顺序的一维数组，fitness值就是经过所有城市所需要的总时间。在离散问题中，领域搜索一般是用交换或者插入算子来决定的，当然还有其他的算子，但是用的最多的是这两种。

        总之就是要跟要解决的实际问题相结合，但是大致连续和离散的思路都和上面差不多。

下面贴出用matlab实现蜜蜂算法的代码，这是我在GitHub上面找到的。我的蜜蜂算法就是按照他的结构写出来的。

% Bees Algorithm (BA) Economic Dispatching (ED)
% P = power plants
% CTotal = the cost
% q= violation
% PL = power loss
% PTotal = all powers
% PD = power demand
clc;
clear;
close all;

%% Problem 
model=MakeModel();
CostFunction=@(x) Cost(x,model);        % Cost Function
nVar=model.nPlant;                      % Number of Decision Variables
VarSize=[1 nVar];                       % Size of Decision Variables Matrix
VarMin=0;                               % Lower Bound of Variables
VarMax=1;                               % Upper Bound of Variables

%% Bees Algorithm Parameters
MaxIt = 50;                               % Maximum Number of Iterations
nScoutBee = 5;                            % Number of Scout Bees
nSelectedSite = round(0.5*nScoutBee);     % Number of Selected Sites
nEliteSite = round(0.4*nSelectedSite);    % Number of Selected Elite Sites
nSelectedSiteBee = round(0.5*nScoutBee);  % Number of Recruited Bees for Selected Sites
nEliteSiteBee = 2*nSelectedSiteBee;       % Number of Recruited Bees for Elite Sites
r = 0.1*(VarMax-VarMin);	              % Neighborhood Radius
rdamp = 0.95;                             % Neighborhood Radius Damp Rate

%% Start
% Empty Bee Structure
empty_bee.Position = [];
empty_bee.Cost = [];
empty_bee.Sol = [];
% Initialize Bees Array
bee = repmat(empty_bee, nScoutBee, 1);
% Create New Solutions
for i = 1:nScoutBee
bee(i).Position = unifrnd(VarMin, VarMax, VarSize);
[bee(i).Cost bee(i).Sol] = CostFunction(bee(i).Position);
end
% Sort
[~, SortOrder] = sort([bee.Cost]);
bee = bee(SortOrder);
% Update Best Solution Ever Found
BestSol = bee(1);
% Array to Hold Best Cost Values
BestCost = zeros(MaxIt, 1);

%% Bees Algorithm Body
for it = 1:MaxIt
% Elite Sites
for i = 1:nEliteSite
bestnewbee.Cost = inf;
for j = 1:nEliteSiteBee
newbee.Position = BeeDance(bee(i).Position, r);
[newbee.Cost newbee.Sol] = CostFunction(newbee.Position);
if newbee.Cost<bestnewbee.Cost
bestnewbee = newbee;
end
end
if bestnewbee.Cost<bee(i).Cost
bee(i) = bestnewbee;
end
end
% Selected Non-Elite Sites
for i = nEliteSite+1:nSelectedSite
bestnewbee.Cost = inf;
for j = 1:nSelectedSiteBee
newbee.Position = BeeDance(bee(i).Position, r);
[newbee.Cost newbee.Sol] = CostFunction(newbee.Position);
if newbee.Cost<bestnewbee.Cost
bestnewbee = newbee;
end
end
if bestnewbee.Cost<bee(i).Cost
bee(i) = bestnewbee;
end
end
% Non-Selected Sites
for i = nSelectedSite+1:nScoutBee
bee(i).Position = unifrnd(VarMin, VarMax, VarSize);
[bee(i).Cost bee(i).Sol] = CostFunction(bee(i).Position);
end
% Sort
[~, SortOrder] = sort([bee.Cost]);
bee = bee(SortOrder);
% Update Best Solution Ever Found
BestSol = bee(1);
% Store Best Cost Ever Found
BestCost(it) = BestSol.Cost;
% Display Iteration Information
disp(['In Iteration No ' num2str(it) ': Bees Cost is = ' num2str(BestCost(it))]);
% Damp Neighborhood Radius
r = r*rdamp;
end

% Plot ITR
plot(BestCost,'-or', 'LineWidth', 2);
xlabel('Bees Iteration');
ylabel('Best Cost Value');
ax = gca; 
ax.FontSize = 12; 
set(gca,'Color','y')
legend({'Bees ED'},'FontSize',12,'FontWeight','bold','TextColor','r');
%% Results
BestSol
%
BestSol.Sol
%
Res=BestSol.Sol.PTotal-BestSol.Sol.PL-model.PD

function P=Parse(x,model)
PminActual=model.Plants.PminActual;
PmaxActual=model.Plants.PmaxActual;
P=PminActual+(PmaxActual-PminActual).*x;
PZ=model.Plants.PZ;
nPlant=model.nPlant;
for i=1:nPlant
for j=1:numel(PZ{i})
if P(i)>PZ{i}{j}(1) && P(i)<PZ{i}{j}(2)
% Correction
if P(i)<(PZ{i}{j}(1)+PZ{i}{j}(2))/2
P(i)=PZ{i}{j}(1);
else
P(i)=PZ{i}{j}(2);
end
end
end
end
end

function out=MCalc(P,model)
alpha=model.Plants.alpha;
beta=model.Plants.beta;
gamma=model.Plants.gamma;
C=alpha+beta.*P+gamma.*P.*P;
% The Cost
CTotal=sum(C);
B=model.B;
B0=model.B0;
B00=model.B00;
% Power Loss
PL=P*B*P'+B0*P'+B00;
% All Powers
PTotal=sum(P);
% Power Demand
PD=model.PD;
PowerBalanceViolation=max(1-(PTotal-PL)/PD,0);

%% Violation
q=5; %Violation (more the better) 
%
z=CTotal*(1+q*PowerBalanceViolation);
out.P=P;
out.PTotal=PTotal;
out.C=C;
out.CTotoal=CTotal;
out.PL=PL;
out.PowerBalanceViolation=PowerBalanceViolation;
out.z=z;
end

function model=MakeModel()
model.PD=1263;

model.Plants.Pmin=[100 50 80 50 50 50];
model.Plants.Pmax=[500 200 300 150 200 120];
model.Plants.alpha=[240 200 220 200 220 190];
model.Plants.beta=[7 10 8.5 11 10.5 12];
model.Plants.gamma=[0.007 0.0095 0.009 0.009 0.008 0.0075];
model.Plants.P0=[440 170 200 150 190 110];
model.Plants.UR=[80 50 65 50 50 50];
model.Plants.DR=[120 90 100 90 90 90];

model.Plants.PminActual = max(model.Plants.Pmin,model.Plants.P0-model.Plants.DR);
model.Plants.PmaxActual = min(model.Plants.Pmax,model.Plants.P0+model.Plants.UR);

model.Plants.PZ{1}={[210 240],[350 380]};
model.Plants.PZ{2}={[90 110],[140 160]};
model.Plants.PZ{3}={[150 170],[210 240]};
model.Plants.PZ{4}={[80 90],[110 120]};
model.Plants.PZ{5}={[90 110],[140 150]};
model.Plants.PZ{6}={[75 85],[100 105]};

model.nPlant=numel(model.Plants.alpha);

model.B=[ 0.0017  0.0012  0.0007 -0.0001 -0.0005 -0.0002
0.0012  0.0014  0.0009  0.0001 -0.0006 -0.0001
0.0007  0.0009  0.0031  0.0000 -0.0010 -0.0006
-0.0001  0.0001  0.0000  0.0024 -0.0006 -0.0008
-0.0005 -0.0006 -0.0010 -0.0006  0.0129 -0.0002
-0.0002 -0.0001 -0.0006 -0.0008 -0.0002  0.0150]/40;

model.B0=1e-3*[-0.3908 -0.1279 0.7047 0.0591 0.2161 -0.6635];

model.B00=0.056;

end

function [z, out]=Cost(x,model)
P=Parse(x,model);
out=MCalc(P,model);
z=out.z;
end

function y = BeeDance(x, r)
nVar = numel(x);
k = randi([1 nVar]);
y = x;
y(k) = x(k)+unifrnd(-r, r);
end