nga.m

用小生境遗传算法实现求解rosenbrock函数的最优值

clear all;
clc;
tic
bound1=[-10 10];%x1可行解空间
bound2=[-10 10];%x2可行解空间
bounds=[bound1;bound2];%可行解期间
numbits=10;%变量用20位2进制编码
chrobits=2*numbits;%染色体基因为数
M=80;%群体规模
N=25;%精英保留的个数
T=1000;%生长代数
pc=0.8;%交叉概率
pm=0.1;%变异概率
L=1.5;%小生镜之间的参数
penalty=10^(-30);%惩罚数
startpop=zeros(M,chrobits+2);%初始群体
Nichepop=zeros(N+M,chrobits+2);%小生镜
scale=(bounds(:,2)-bounds(:,1))./(2.^numbits-1);%可行解的精度
cs=0:numbits:chrobits;
endpop = zeros(M,chrobits+2); %下一代群体
c1 = zeros(1,chrobits+2); 	%An individual
c2 = zeros(1,chrobits+2); 	%An individual
%初始化群体
for i=1:50
    for j=1:chrobits
        if rand>0.5
            startpop(i,j)=1;
        end
    end
end
%Evaluating the fitniss of the individuals
for i=1:M
    fitness=fit(startpop(i,1:chrobits));
    startpop(i,chrobits+1)=fitness;
end
%降序排列
[result,IX]=sort(startpop(:,chrobits+1),'descend');
startpop=startpop(IX,:);
totalFitness=sum(startpop(:,chrobits+1));
startpop(:,chrobits+2)=startpop(:,chrobits+1)/totalFitness;
Nichepop(1:N,:)=startpop(IX(1:N),:);
t=0;%当前代数
while(t<T)
    %选择算子
    selectpro=cumsum(startpop(:,chrobits+2));
    rands=sort(rand(M,1),'ascend');
    index1=1;
    index2=1;
    while index1<=M
        if rands(index1)<selectpro(index2)
            endpop(index1,:)=startpop(index2,:);
            index1=index1+1;
        else
            index2=index2+1;
        end
    end
    %交叉算子
    %随机配对
    index=1:M;
    for i=1:M
        point=round(M*rand+0.5);
        temp=index(i);
        index(i)=index(point);
        index(point)=temp;
    end
    %交叉运算
    for i=1:2:M
        if rand<pc
            point=round(chrobits*rand+0.5);%随机产生交叉点
            b1=endpop(index(i),:);
            b2=endpop(index(i+1),:);
            c1=[b1(1:point),b2(point+1:chrobits+2)];
            c2=[b2(1:point),b1(point+1:chrobits+2)];
            endpop(index(i),:)=c1;
            endpop(index(i+1),:)=c2;
            endpop(index(i),chrobits+1)=fit(c1(1:chrobits));
            endpop(index(i+1),chrobits+1)=fit(c2(1:chrobits));
        end
    end
    %变异算子
    for i=1:M
        if rand<pm
            point=round(rand*chrobits+0.5);
            endpop(i,point)=1-endpop(i,point);
            endpop(i,chrobits+1)=fit(endpop(i,1:chrobits));
        end
    end
    Nichepop(N+1:N+M,:)=endpop;
    [result,IX]=sort(Nichepop(:,chrobits+1),'descend');
   Nichepop=Nichepop(IX,:);
   %小生镜淘汰运算
   %二进制转换为格雷码
   for i=1:M+N
       Nichepopgl(i,chrobits)=Nichepop(i,chrobits);
       for j=chrobits-1:-1:1
           if Nichepop(i,j+1)==Nichepop(i,j)
               Nichepopgl(i,j)=0;
           else 
               Nichepopgl(i,j)=1;
           end
       end
   end 
   %小生镜
   for i=1:M+N-1
       for j=(i+1):M+N
           hamdistance=0;
           for k=1:chrobits
               hamdistance=hamdistance+(Nichepop(i,k)-Nichepop(j,k))^2;
           end
           if sqrt(hamdistance)<L
               Nichepop(j,chrobits+1)=penalty;
           end
       end
   end   
   %格雷码转换为二进制   
   for i=1:M+N
       Nichepop(i,chrobits)=Nichepopgl(i,chrobits);
       for j=chrobits-1:-1:1
           if Nichepop(i,j+1)==Nichepopgl(i,j)
               Nichepop(i,j)=0;
           else 
               Nichepop(i,j)=1;
           end
       end
   end      
   [result,IX]=sort(Nichepop(:,chrobits+1),'descend');
   Nichepop=Nichepop(IX,:);
   startpop=Nichepop(1:M,:);%适应度高的
   Nichepop(1:N,:)=startpop(1:N,:);
   totalfitness=sum(startpop(:,chrobits+1));
   startpop(:,chrobits+2)=startpop(:,chrobits+1)/totalfitness;
   t=t+1;
   chroinfo(t,1)=startpop(1,chrobits+1);
   chroinfo(t,2)=mean(startpop(:,chrobits+1));
   chroinfo(t,3)=std(startpop(:,chrobits+1));
end
toc;
plot(chroinfo(:,1),'r+');
text(100,chroinfo(100,1)+0.5,['Red + point represent the max fitness']);  
hold on;
plot(chroinfo(:,2),'gp');
text(150,chroinfo(150,2)+1,['Green pentagram point represent the mean fitness']); 
hold on;
plot(chroinfo(:,3),'b*');
text(250,chroinfo(250,3)+1,['Blue * point represent the std fitness']); 
title('NicheGA');
xlabel('Generation');
ylabel('Statistical values of Fitness');
%下面的程序主要是用于求函数的全局最优解和对应的x1和x2
%funx矩阵的第一列、第二列和第三列表示x1、x2和对应的函数值
for i=1:M
    for j=1:2
        varbits=startpop(i,(cs(j)+1):cs(j+1));
        x(j)=sum(2.^(size(varbits,2)-1:-1:0).*varbits)*scale(j)+bounds(j,1);
        funx(i,j)=x(j);
    end
    val1=0;
    val2=0;
    for k=1:5
        val1=val1+k*cos((k+1)*x(1)+k);
        val2=val2+k*cos((k+1)*x(2)+k);
    end
    funx(i,3)=val1*val2;
end