人工蚂蚁算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%function [x,y, minvalue] = AA(func)% Example [x, y,minvalue] = AA('Foxhole')clc;tic;subplot(2,2,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot 1 draw(func);title([func, ' Function']);%初始化各参数Ant=100;%蚂蚁规模ECHO=200;%迭代次数step=0.01*rand(1);%局部搜索时的步长temp=[0,0];%各子区间长度start1=-100;end1=100;start2=-100;end2=100;Len1=(end1-start1)/Ant;Len2=(end2-start2)/Ant;%P = 0.2;%初始化蚂蚁位置for i=1:AntX(i,1)=(start1+(end1-start1)*rand(1));X(i,2)=(start2+(end2-start2)*rand(1));%func=AA_Foxhole_Func(X(i,1),X(i,2));val=feval(func,[X(i,1),X(i,2)]);T0(i)=exp(-val);%初始信息素,随函数值大,信息素浓度小,反之亦然 %%%%%***************************************** ****************************end;%至此初始化完成for Echo=1:ECHO %开始寻优%P0函数定义,P0为全局转移选择因子a1=0.9;b1=(1/ECHO)*2*log(1/2);f1=a1*exp(b1*Echo);a2=0.225;b2=(1/ECHO)*2*log(2);f2=a2*exp(b2*Echo);if Echo<=(ECHO/2)P0=f1;elseP0=f2;end;%P函数定义，P为信息素蒸发系数a3=0.1; b3=(1/ECHO).*log(9);P=a3*exp(b3*Echo);lamda=0.10+(0.14-0.1)*rand(1);%全局转移步长参数Wmax=1.0+(1.4-1.0)*rand(1);%步长更新参数上限Wmin=0.2+(0.8-0.2)*rand(1);%步长更新参数下限%寻找初始最优值T_Best=T0(1);for j=1:Antif T0(j)>=T_BestT_Best=T0(j);BestIndex=j;end;end;W=Wmax-(Wmax-Wmin)*(Echo/ECHO); %局部搜索步长更新参数for j_g=1:Ant %全局转移概率求取,当该蚂蚁随在位置不是bestindex时if j_g~=BestIndexr=T0(BestIndex)-T0(j_g);Prob(j_g)=exp(r)/exp(T0(BestIndex));else%当j_g=BestIndex的时候进行局部搜索if rand(1)<0.5temp(1,1)=X(BestIndex,1)+W*step; temp(1,2)=X(BestIndex,2)+W*step;elsetemp(1,1)=X(BestIndex,1)-W*step; temp(1,2)=X(BestIndex,2)-W*step;end;Prob(j_g)=0;%bestindex的蚂蚁不进行全局转移end;X1_T=temp(1,1);X2_T=temp(1,2);X1_B=X(BestIndex,1);X2_B=X(BestIndex,2);%func1 =AA_Foxhole_Func(X1_T,X2_T); %%%%%%%%%%%********* ******************************************%F1_T=func1;F1_T=feval(func,[X(i,1),X(i,2)]);F1_B=feval(func,[X1_B,X2_B]);%F1_T=(X1_T-1).^2+(X2_T-2.2).^2+1;%func2 =AA_Foxhole_Func(X1_B,X2_B); %%%%%%%%%%%%%******** *******************************************%F1_B=func2;%F1_B=(X1_B-1).^2+(X2_B-2.2).^2+1;if exp(-F1_T)>exp(-F1_B)X(BestIndex,1)=temp(1,1);X(BestIndex,2)=temp(1,2);end;end;for j_g_tr=1:Antif Prob(j_g_tr)<P0X(j_g_tr,1)=X(j_g_tr,1)+lamda*(X(BestIndex,1)-X(j_g_tr,1));%Xi=Xi+lamda*(Xbest-Xi)21X(j_g_tr,2)=X(j_g_tr,2)+lamda*(X(BestIndex,2)-X(j_g_tr,2));%Xi=Xi+lamda*(Xbest-Xi)X(j_g_tr,1)=bound(X(j_g_tr,1),start1,end1);X(j_g_tr,2)=bound(X(j_g_tr,2),start2,end2);elseX(j_g_tr,1)=X(j_g_tr,1)+((-1)+2*rand(1))*Len1;%Xi=Xi+rand(-1,1)*Len1X(j_g_tr,2)=X(j_g_tr,2)+((-1)+2*rand(1))*Len2;%Xi=Xi+rand(-1,1)*Len2X(j_g_tr,1)=bound(X(j_g_tr,1),start1,end1);X(j_g_tr,2)=bound(X(j_g_tr,2),start2,end2);end;end;%信息素更新subplot(2,2,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Plot 1bar([X(BestIndex,1) X(BestIndex,2)],0.25);%colormap (cool);axis([0 3 -40 40 ]) ;title ({date;['Iteration ',num2str(Echo)]});xlabel(['Min_x = ',num2str(X(BestIndex,1)),' ', 'Min_y = ', num2str(X(BestIndex,2))]);for t_t=1:Ant%func=AA_Foxhole_Func(X(t_t,1),X(t_t,2)); val1=feval(func,[X(t_t,1),X(t_t,2)]);T0(t_t)=(1-P)*T0(t_t)+(exp(-val1)); %**************************************** *********************************end;[c_iter,i_iter]=max(T0); %求取每代全局最优解 minpoint_iter=[X(i_iter,1),X(i_iter,2)];%func3 =AA_Foxhole_Func(X(i_iter,1),X(i_iter,2)); %%%%%%% %%*********************************************** ****************************val2=feval(func,[X(i_iter,1),X(i_iter,2)]); minvalue_iter= val2;%minvalue_iter=(X(i_iter,1)-1).^2+(X(i_iter,2)-2.2).^2+1;min_local(Echo)=minvalue_iter;%保存每代局部最优解subplot(2,2,3);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Plot 2plot(X(BestIndex,1),X(BestIndex,2),'rs','MarkerFa ceColor','r', 'MarkerSize',8), grid on;title (['Global Min Value = ',num2str(minvalue_iter)]);hold on;plot(X(:,1),X(:,2),'g.'), pause (0.02); hold off ;axis([-100 100 -100 100]);grid on;%将每代全局最优解存到min_global矩阵中if Echo >= 2if min_local(Echo)<min_global(Echo-1)min_global(Echo)=min_local(Echo);elsemin_global(Echo)=min_global(Echo-1);end;elsemin_global(Echo)=minvalue_iter;end;subplot(2,2,4);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 3min_global=min_global';index(:,1)=1:ECHO;plot(Echo, min_global(Echo),'y*')%axis([0 ECHO 0 10]);hold on;title ([func,' (X) = ',num2str(minvalue_iter)],'Color','r');xlabel('iteration');ylabel('f(x)');grid on;end;%ECHO循环结束[c_max,i_max]=max(T0);minpoint=[X(i_max,1),X(i_max,2)];%func3 =AA_Foxhole_Func(X(i_max,1),X(i_max,2)); %%%****** ************************************************* ******************%minvalue = func3;minvalue=feval(func,[X(i_max,1),X(i_max,2)]);x=X(BestIndex,1);y=X(BestIndex,2);runtime=toc21人工免疫算法function [x,y,fx,vfx,vmfit,P,vpm] =AI(func,gen,n,pm,per);% Example [x,y,fx] = AI('Foxhole')subplot(2,2,1);draw(func);title( [func, ' Function']);if nargin == 1,% gen = 200; n = round(size(P,1)/2); pm =0.0005; per = 0.0; fat = 10;%gen = 250; n = size(P,1); pm = 0.01; per = 0.0; fat = .1;P = cadeia(200,44,0,0,0);gen = 40; n = size(P,1); pm = 0.2; per = 0.0; fat = 0.1;end;while n <= 0,n = input('n has to be at least one. Type a new value for n: ');end;xmin=-100;xmax=100;ymin=-100;ymax=100;x = decode(P(:,1:22),xmin,xmax); y =decode(P(:,23:end),ymin,ymax);%fit = eval(f);%fit=AI_Foxhole_Func(x,y);%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%fit=feval(func,[x' y']);%imprime(1,vxp,vyp,vzp,x,y,fit,1,1);% Hypermutation controlling parameterspma = pm; itpm = gen; pmr = 0.8;% General defintionsvpm = []; vfx = []; vmfit = []; valfx = 1; [N,L] = size(P); it = 0; PRINT = 1;% Generationswhile it <= gen & valfx <= 100,x = decode(P(:,1:22),xmin,xmax); y =decode(P(:,23:end),ymin,ymax); T = []; cs = [];%fit = eval(f);%fit=AI_Foxhole_Func(x,y);%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%fit=feval(func,[x' y']);[a,ind] = sort(fit);valx = x(ind(end-n+1:end)); valy = y(ind(end-n+1:end));fx = a(end-n+1:end); % n best individuals (maximization)% Reproduction [T,pcs] = reprod(n,fat,N,ind,P,T);% HypermutationM = rand(size(T,1),L) <= pm;T = T - 2 .* (T.*M) + M;T(pcs,:) = P(fliplr(ind(end-n+1:end)),:);% New Re-Selection (Multi-peak solution)x = decode(T(:,1:22),xmin,xmax); y =decode(T(:,23:end),ymin,ymax);%fit=AI_Foxhole_Func(x,y);%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%fit=feval(func,[x' y']);%fit = eval(f);pcs = [0 pcs];for i=1:n,[out(i),bcs(i)] =min(fit(pcs(i)+1:pcs(i+1))); % Mimimazion problem %%%*************************bcs(i) = bcs(i) + pcs(i);end;P(fliplr(ind(end-n+1:end)),:) = T(bcs,:);% Editing (Repertoire shift)nedit = round(per*N); it = it + 1;P(ind(1:nedit),:) = cadeia(nedit,L,0,0,0);pm = pmcont(pm,pma,pmr,it,itpm); valfx =min(fx); %*************************************** **********************vpm = [vpm pm]; vfx = [vfx valfx]; vmfit = [vmfit mean(fit)];disp(sprintf('It.: %d pm: %.4f x: %2.2fy: %2.2f Av.: %2.2ff(x,y): %2.3f',it,pm,valx(1),valy(1),vmfit(1),val fx));subplot(2,2,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 1bar([valx(1) valy(1)],0.25);axis([0 3 -40 40 ]) ;title (['Iteration ', num2str(it)]); pause (0.1);subplot(2,2,3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 2plot(valx(1),valy(1),'rs','MarkerFaceColor','r', 'MarkerSize',8)hold on;%plot(x(:,1),x(:,2),'k.');set(gca,'Color','g')hold off;grid on;axis([-100 100 -100 100 ]) ;21title(['Global Min =',num2str(valfx)]);xlabel(['Min_x= ',num2str(valx(1)),'Min_y= ',num2str(valy(1))]);subplot(2,2,4); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 3 plot(it,valfx ,'k.')axis([0 gen 0 10]);hold on;title ([func ,'(X) = ', num2str(valfx)]);xlabel('iteration');ylabel('f(x)');grid on;end; % end while%imprime(PRINT,vxp,vyp,vzp,x,y,fit,it,1);x = valx(1); y = valy(1); fx =min(fx); %***********************************************************************% x = P(ind(end),1:22); y = P(ind(end),23:44); fx= max(fx);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% Plot 4% --------------------- %% INTERNAL SUBFUNCTIONS% --------------------- %% Printfunction [] =imprime(PRINT,vx,vy,vz,x,y,fx,it,mit);% x,fx -> actual values% vxplot, vplot -> original (base) functionif PRINT == 1,if rem(it,mit) == 0,mesh(vx,vy,vz); hold on; axis([-100 100 -100 100 0 500]);xlabel('x'); ylabel('y'); zlabel('f(x,y)');plot3(x,y,fx,'k*'); drawnow; hold off;end;end;% Reproductionfunction [T,pcs] = reprod(n,fat,N,ind,P,T);% n -> number of clones% fat -> multiplying factor% ind -> best individuals% T -> temporary population% pcs -> final position of each cloneif n == 1,cs = N;T = ones(N,1) * P(ind(1),:);else,for i=1:n,% cs(i) = round(fat*N/i);cs(i) = round(fat*N); pcs(i) = sum(cs);T = [T; ones(cs(i),1) * P(ind(end-i+1),:)];end;end;% Control of pmfunction [pm] = pmcont(pm,pma,pmr,it,itpm);% pma -> initial value% pmr -> control rate% itpm -> iterations for restoringif rem(it,itpm) == 0,pm = pm * pmr;if rem(it,10*itpm) == 0,pm = pma;end;end;% Decodify bitstringsfunction x = decode(v,min,max);% x -> real value (precision: 6)% v -> binary string (length: 22)v = fliplr(v); s = size(v);aux = 0:1:21; aux = ones(s(1),1)*aux;x1 = sum((v.*2.^aux)');x = min + (max-min)*x1 ./ 4194303;function [ab,ag] = cadeia(n1,s1,n2,s2,bip)%default parameter value seetingif nargin == 2,n2 = n1; s2 = s1; bip = 1;elseif nargin == 4,bip = 1;end;% Antibody (Ab) chainsab = 2 .* rand(n1,s1) - 1;%create n1 row s1 column array, its value range is between -1 or 1if bip == 1,ab = hardlims(ab);else,ab = hardlim(ab);end;% Antigen (Ag) chainsag = 2 .* rand(n2,s2) - 1;if bip == 1,ag = hardlims(ag);else,ag = hardlim(ag);end;% End Function CADEIA21%------免疫粒子群优化算法（Artificial Immune - Particle Swarm Optimization）function [x,y,Result]=PSO_AI(func)% Example [x, y,minvalue] = PSO_AI('Foxhole') clc;subplot(2,2,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot 1 draw(func);title([func, ' Function']);ticformat long;%------给定初始化条件----------------------------------------------c1=1.4962; %学习因子1c2=1.4962; %学习因子2w=0.7298; %惯性权重MaxDT=200; %最大迭代次数D=2; %搜索空间维数（未知数个数）N=100; %初始化群体个体数目eps=10^(-20); %设置精度(在已知最小值时候用)DS=10; %每隔DS次循环就检查最优个体是否变优replaceP=0.6; %粒子的概率大于replaceP将被免疫替换minD=1e-015; %粒子间的最小距离Psum=0; %个体最佳的和range=100;count = 0;%------初始化种群的个体------------for i=1:Nfor j=1:Dx(i,j)=-range+2*range*rand; %随机初始化位置v(i,j)=randn; %随机初始化速度endend%------先计算各个粒子的适应度，并初始化Pi和Pg----------------------for i=1:N%p(i)=Foxhole(x(i,:),D); %fitness是计算各个粒子适应度的函数，见文件fitness.m %%%%%%%%****************************** ***********************p(i)=feval(func,x(i,:));y(i,:)=x(i,:);endpg=x(1,:); %Pg为全局最优for i=2:Niffeval(func,x(i,:))<feval(func,pg) %%%%******** ************************************************* ************************************* pg=x(i,:);endend%------进入主要循环，按照公式依次迭代，直到满足精度要求------------for t=1:MaxDTfor i=1:Nv(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));x(i,:)=x(i,:)+v(i,:);iffeval(func,x(i,:))<p(i) %%%%%%%%%%%%%%*********** ************************************************* **********************************p(i)=feval(func,x(i,:)); %%%%%%%%%%%%%%%******** ************************************************* *****************************y(i,:)=x(i,:);endifp(i)<feval(func,pg) %%%%%%%%%%%%%%%%%%%%******** ************************************************* *********************************pg=y(i,:);subplot(2,2,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 1bar(pg,0.25);axis([0 3 -40 40 ]) ;title (['Iteration ', num2str(t)]); pause (0.1);subplot(2,2,3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot 2plot(pg(1,1),pg(1,2),'rs','MarkerFaceColor','r', 'MarkerSize',8)hold on;plot(x(:,1),x(:,2),'k.');set(gca,'Color','g')hold off;grid on;axis([-100 100 -100 100 ]) ;title(['Global Min =',num2str(p(i))]);xlabel(['Min_x= ',num2str(pg(1,1)),' Min_y= ',num2str(pg(1,2))]);endendPbest(t)=feval(func,pg) ; %%%%%%%%************* ************************************************* **********************************************21% if Foxhole(pg,D)<eps %如果结果满足精度要求则跳出循环% break;% end%-----------开始进行免疫－－－－－－－－－－－－－－－－if t>DSif mod(t,DS)==0 && (Pbest(t-DS+1)-Pbest(t))<1e-020 %如果连续DS代数，群体中的最优没有明显变优，则进行免疫.%在函数测试的过程中发现，经过一定代数的更新，个体最优不完全相等，但变化非常非常小，%我认为这个时候也应用免疫了，所以我没有用“Pbest(t-DS+1)＝Pbest(t)”作为判断条件，%不过“(Pbest(t-DS+1)-Pbest(t))<1e-020”是否合理也值得探讨。