中文互联网上面好像没有关于蜜蜂算法的详细介绍,正好最近在使用蜜蜂算法,在次介绍一下。蜜蜂算法(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