%% 算法说明：荧火虫算法
clc; %清屏
clear all; %清除变量
format long; %确定精度

%% 各参数初始化开始
domx=[-5.12,5.12;-5.12,5.12];
% domx=[-2.048,2.048;-2.048,2.048]; %解空间
rho=0.4; %荧光素挥发因子
gamma=0.6; %适应度提取比例
beta=0.08; %邻域变化率
nt=5; %邻域阀值(邻域荧火虫数)
s=0.03; %步长
lo=5; %荧光素浓度
rs=5.12; %感知半径
r0=5.12;
% rs=2.048; %感知半径
% r0=2.048; %决策半径
%各参数初始化结束
iter_max=200; %最大迭代次数
%% 分配空间开始

m=size(domx,1); %函数空间维数
global n;
n=50; %种群规模
address=zeros(n,m); %分配荧火虫地址空间
value=zeros(n,1); %分配适应度存放空间
li=zeros(n,1); %分配荧光素存放空间
rdi=zeros(n,1); %分配荧火虫决策半径存放空间
uu=zeros(1,iter_max);
%% 荧火虫常量初始化开始

for i=1:m
address(:,i)=(domx(i,1)+(domx(i,2)-domx(i,1))*rand(n,1)); %初始化地址
% address(i,:)=5*rands(n,1); %随机产生初萤火虫所在位置
end

f=fun(address);
x=-5.12:.1:5.12;
% x=-2.048:.05:2.048;
[x,y] = meshgrid(x);
figure(1);
plot3(address(:,1),address(:,2),f,'k*')
hold on;
grid on;
z=-(x.^2-10*cos(2*pi.*x)+10+y.^2-10*cos(2*pi.*y)+10); %Rastrigin'函数
% z=-(x.^2+y.^2); %目标函数 J2
% z=-(20+x.^2-10*cos(2*pi.*x)+y.^2-10*cos(2*pi.*y)); %目标函数 J1
% z=-(-20*exp(-0.2*sqrt((x.^2+y.^2)/2))-exp((cos(2*pi*x)+cos(2*pi*y))/2)+20+exp(1)); %%目
标函数
% z=-(0.5+sin(sqrt(x.^2 + y.^2).^2-0.5)./(1+0.001*(x.^2 + y.^2)).^2);%目标函数 J3
mesh(x,y,z)
xlabel('x轴');
ylabel('y轴');
zlabel('z轴');
title('萤火虫初始分布图');

li(:,1)=lo; %荧光素初始值
rdi(:,1)=r0; %决策半径初始值
t=1; %迭代累计量
%荧火虫常量初始化结束

%% 迭代开始
while(t<=iter_max)

li=(1-rho)*li+ gamma * fun(address); %.更新荧光素值li
%各荧火虫移动过程开始
for i=1:n
Ord_number=[]; %存放荧火虫序号
for j=1:n
if (norm(address(j,:)-address(i,:))Ord_number(numel(Ord_number)+1)=j;
end
end

%计算Ord_number各元素被选择概率，确定j位置
if ~isempty(Ord_number) %先判断Ord_number个数不为空
Ord_num_li=li(Ord_number,1); %选出Ord_number荧光素
Sum_Ord_li=sum(Ord_num_li); %Ord_number荧光素和
Mol=Ord_num_li-li(i,1); %分子(lj-li)
Den=Sum_Ord_li-li(i,1); %分母
Pij=Mol./Den; %计算各元素被选择概率
Pij=cumsum(Pij); %累计
Pij=Pij./Pij(end); %归一化
Pos=find(randj=Ord_number(Pos(1)); %确定 j 的位置
%荧火虫i向位置j移动
address(i,:)=address(i,:)+ s*(address(j,:)-address(i,:))/norm(address(j,:)-address(i,:));
address(i,:)=range1(address(i,:),domx);%限制范围
end
%更新决策半径
rdi(i)=rdi(i)+beta*(nt-length(Ord_number));
rdi(i)=min(rs,max(0,rdi(i)));
uu(t)=-max(fun(address));
end
t=t+1;
% iter_max迭代结束

g=[];
g=fun(address);
x=-5.12:.1:5.12;
% x=-2.048:.05:2.048;
[x,y] = meshgrid(x);
figure(2);
plot3(address(:,1),address(:,2),g,'k*');
hold on;
z=-(x.^2-10*cos(2*pi.*x)+10+y.^2-10*cos(2*pi.*y)+10); %Rastrigin'函数
% z=-(x.^2+y.^2);%目标函数 J2
% z=-(20+x.^2-10*cos(2*pi*x)+y.^2-10*cos(2*pi*y)); %目标函数 J1
% z=-(-20*exp(-0.2*sqrt((x.^2+y.^2)/2))-exp((cos(2*pi*x)+cos(2*pi*y))/2)+20+exp(1)); %%目
标函数
% z=-(0.5+sin(sqrt(x.^2 + y.^2).^2-0.5)./(1+0.001*(x.^2 + y.^2)).^2);%目标函数 J3
grid on;
mesh(x,y,z);
title('萤火虫运动轨迹');
xlabel('x轴');
ylabel('y轴');
zlabel('z轴');
hold off;
end

figure('name','收敛曲线');
plot(uu)
title(['萤火虫最优值曲线 ' '终止次数＝' num2str(iter_max)]);
xlabel('迭代次数');
ylabel('最优值');

%% 输出最优结果
value=fun(address);
disp('最优值为：')
num=find(value==max(value));%最大值序号
MaxValue=-max(value)
disp('最优解为：')
BestAddress = address(num,:)

%目标函数
function y = fun(x)
% y=x(:,1).^2+x(:,2).^2;
% y=0.5+sin(sqrt(x(:,1).^2 + x(:,2).^2).^2-0.5)./(1+0.001*(x(:,1).^2 + x(:,2).^2)).^2;%目标函数 J3
% y=20+x(:,1).^2-10*cos(2*pi*x(:,1))+x(:,2).^2-10*cos(2*pi*x(:,2)); %目标函数 J1
% y=-20*exp(-0.2*sqrt((x(:,1).^2+x(:,2).^2)/2))...
% -exp((cos(2*pi*x(:,1))+cos(2*pi*x(:,2)))/2)+20+exp(1); %%目标函数
y=x(:,1).^2-10*cos(2*pi.*x(:,1))+10+x(:,2).^2-10*cos(2*pi.*x(:,2))+10; %Rastrigin'函数
y=-y;