clear;clc;close all;

%自变量取值范围

range_x=[ones(1,1),-ones(1,1)]*500;

%维度

n=size(range_x,1);

%种群数量

gn=400;

%遗传迭代时进入下一代的数量

m=50 ;

%从记忆细胞种拿出的个数

l=200;

%迭代次数

times=200;

%免疫作用次数

t=10;

%记忆细胞

remember=zeros(n,gn);

for k=1:n

remember(k,:)=(rand(1,gn))*(range_x(k,2)-range_x(k,1))+range_x(k,1);

end

tic

for p=1:t

tem=zeros(n,gn-l);

for k=1:n

tem(k,:)=(rand(1,gn-l))*(range_x(k,2)-range_x(k,1))+range_x(k,1);

end

group=[remember(:,randperm(gn,l)),tem];

%设置当前最优解

best_value=zeros(1,times);

for k=1:times

y=f(group);

%全部变成正值

if min(y)<0

tem=y-min(y)*1.0001;

else

tem=y+0.1;

end

%值越小适应越好

tem=1./tem;

child=zeros(n,gn);

%挑选m个种群进入下一代

for i=1:m

%轮盘赌选择，适应大的选择概率大

temp=zeros(1,gn-i+1);

for j=1:gn-i+1

temp(j)=sum(tem(1:j));

end

temp=temp/temp(gn-i+1);

%保留最合适的物种

choose=find(temp>rand(1),1);

child(:,i)=group(:,choose);

group=[group(:,1:choose-1),group(:,choose+1:end)];

tem=[tem(1:choose-1),tem(choose+1:end)];

end

%染色体交换,保留的物种产生后代时发生基因重组

for i=1:floor((gn-m)/2)

exchange=randperm(m,2);

a=rand(n,1);

child(:,i*2-1+m)=a.*child(:,exchange(1))+(1-a).*child(:,exchange(2));

child(:,i*2+m)=(1-a).*child(:,exchange(1))+a.*child(:,exchange(2));

end

if mod(gn-m,2)==1

exchange=randperm(m,2);

child(:,gn)=(child(:,exchange(1))+child(:,exchange(2)))/2;

end

%基因重组的过程中可能发生染色体变异

if rand(1)<0.1

exchange=randperm(gn-m,1);

a=rand(1);

for j=1:n

child(j,exchange+m)=a.*child(j,exchange+m)+(1-a).*(rand(1)*(range(j,2)-range(j,1))+range(j,1));

end

end

%重组之后后代变成当前种群

group=child;

best_value(k)=min(f(group));

if k>5&&abs(best_value(k)-best_value(k-5))<1e-5

break;

end

end

%将本次免疫过程中最好的和记忆细胞比较，选取最好的作为记忆细胞

if min(f(group))<=min(f(remember))

[~,index]=min(f(group));

remember=ones(n,gn).*repmat(group(:,index),1,gn);

else

[~,index]=min(f(remember));

remember=ones(n,gn).*repmat(group(:,index),1,gn);

end

end

time=toc;

disp(['用时：',num2str(time),'秒'])

[mini,index]=min(f(remember));

disp(['fmin=',num2str(mini)]);

for k=1:n

disp(['x',num2str(k),'=',num2str(remember(k,index))]);

end

if n==1

hold on;

plot(remember(index),mini,'ro');

plot_x=range_x(1):(range_x(2)-range_x(1))/1000:range_x(2);

plot_y=f(plot_x);

plot(plot_x,plot_y);

text((range_x(1)+range_x(2))/2,max(plot_y)+0.1*(max(plot_y)-min(plot_y)),['用时：',num2str(time),'秒']);

hold off;

end

if n==2

%所求最小值的函数

func=@(x1,x2)x1.*sin(sqrt(abs(x1)))+x2.*sin(sqrt(abs(x2)));

plot_x=range_x(1,1):(range_x(1,2)-range_x(1,1))/1000:range_x(1,2);

plot_y=range_x(2,1):(range_x(2,2)-range_x(2,1))/1000:range_x(2,2);

[plot_x,plot_y] =meshgrid(plot_x,plot_y);

plot_z=func(plot_x,plot_y);

surf(plot_x,plot_y,plot_z);

xlabel('x1');

ylabel('x2');

zlabel('y');

hold on;

plot3(remember(1,index),remember(2,index),mini,'ko')

text((range_x(1,1)+range_x(1,2))/2,(range_x(2,1)+range_x(2,2))/2,max(max(plot_z))+0.5*(max(max(plot_z))-min(min(plot_z))),['用时：',num2str(time),'秒']);

hold off;

end