%使用遗传算法求解
% f(x1,x2)=x1^2+x2^2-0.3*cos(3*pi*x1)-0.4*cos(4*pi*x2)+0.7的极小值,并画出每一代中个体适配函
% 数的平均值和最小值随迭代次数的变化关系

%程序总体说明：
% 编码方式说明： 由于期望所求的函数极小值包含了平方项，而正弦和余弦项的范围只能是-1---1，因此优化范围的x1,x2只能在-1--1之间，编码时
% 首先将x1,x2加1变到0--2之间，之后为了编码需要，设采样精度是0.01 ,将0--2之间的数乘以100,得到0--200之间的整数，另外为了编码的需要，
%设编码为8位，编码成0--255的二进制码。,每个个体用16位的二进制串表示，的前8位是x1,后8位是x2
%适配制计算说明：%根据编码方式，将每一个个体的编码变化到-1--1之间，之后按照函数f计算
%可以选择  1--直接利用函数的倒数作为适配值计算； 2--使用排序法计算， 3--使用1-J(i)/sum(J)方法计算
%复制操所机制：1--轮盘赌方法  ； 2--保留最优淘汰法
%交叉机制: 使用掩码法进行交叉操作.

clear all;
%使用
disp \'f(x1,x2)=x1^2+x2^2-0.3*cos(3*pi*x1)-0.4*cos(4*pi*x2)+0.7\'
disp \'输入初始化参数\'
inum=input(\'请输入种群大小(回车为缺省值30个个体):    \');
if isempty(inum)
    inum=50;   %种群中的个体数
end
count=input(\'请输入迭代次数(回车为缺省值100代):     \');
if isempty(count)
    count=100;    %迭代次数
end
pcross=input(\'请输入交叉概率(回车为缺省值0.6):      \');
if isempty(pcross)
    pcross=0.6;   %交叉概率
end
pchange=input(\'请输入变异概率(回车为缺省值0.001):      \');
if isempty(pchange)
    pchange=0.001;     %变异概率
end
disp \'程序实现机制选择\'
%程序实现机制选择
Fitflag=input(\'请输入适配度的计算方法(1为1/f(x1,x2), 2为排序法,3为1-J(i)/sum(J),缺省为2 ):       \'  );
if (isempty(Fitflag))||(Fitflag>3)
    Fitflag=2;
end
Copyflag=input(\'请输入适配度的计算方法(1为轮盘赌法,2为最优保存淘汰法,缺省为2):          \');
if (isempty(Copyflag))||(Copyflag>2)
    Copyflag=2;
end

%程序计算开始
constk=5; %用来控制适配制之间的常数
%初始种群
Individual=rand(inum,16);       % 表示的个体的矩阵
Individual=round(Individual);   %变成01 串
J=zeros(inum,1);                %存储各个个体的函数值
Jorder=zeros(inum,1);        %函数值的大小顺序
Adap=zeros(inum,1);          %适配度
MeanAdap=zeros(count,1);        %存储每一代的适配度均值
MinAdap=zeros(count,1);         %存储每一代的适配度最小值
%初始化随便找一个Xg;

%迭代
for k=1:count
    [inum1 col]=size(Individual);
    for i=1:inum1   %每一个个体计算         
    %适配度计算
        %转化成十进制数
        x1=0;
        x2=0;
        for j=1:8
            ex=8-j;
            x1=x1+Individual(i,j)*2^ex;   % 计算x1  0--256
            x2=x2+Individual(i,j+8)*2^ex; % 计算x2  0--256           
        end  %end of j
        x1=x1/100-1.28;   % 变换 x1到 -1.28:1.28
        x2=x2/100-1.28;   % 变换 x2到 -1.28:1.28
        %适配度计算
        J(i,k)=(x1^2+x2^2-0.3*cos(3*pi*x1)-0.4*cos(4*pi*x2)+0.7);
    end   %end of i
    switch Fitflag
        case 1        %1为1/f(x1,x2), 2为排序法,3为1-J(i)/sum(J)
            Adap=1/J(:,k);            
        case 2
            %排序,按照顺序计算适配值
            [B,Jorder] = sort(J(:,k));  %次序是从小到大排列
            for i=1:inum1
                tempi=Jorder(i);
                Adap(tempi)=constk*(inum1-i)/inum1;
            end            
        case 3
            tempj=ones(inum1,1);
            Adap=tempj-J(:,k)/sum(J(:,k));       
    end       %end of switch            
    MeanAdap(k)=mean(Adap);
    MinAdap(k)=min(Adap);
    MeanJ(k)=mean(J(:,k));
    [MinJ(k),Indexmin]=min(J(:,k));
    %得到最优解，保存在Xg中
    if (k==1)
        Xg(1)=MinJ(k);
        MinIndividual=Individual(Indexmin,:);
    else
        if (MinJ(k)<Xg(k-1))   
            Xg(k)=MinJ(k);
            MinIndividual=Individual(Indexmin,:);
        else
            Xg(k)=Xg(k-1);
        end
    end     % end of if
    Adap1=Adap/sum(Adap);   %适配度的比例
    TempIndividual=Individual;
    Individual=zeros(inum,16);    
    %复制
    switch Copyflag
        case 1     %1为轮盘赌法,2为最优保存淘汰法
            %随机产生一个数
            for i=1:inum1
                if (i==1)
                    route(i)=Adap1(i);        %产生轮盘
                else
                    route(i)=route(i-1)+Adap1(i);
                end        
            end
            for i=1:inum         %复制产生inum个个体
                prob=rand;
                for j=1:inum1    % 根据轮盘赌概率决定哪个个体被复制
                    if (j==1)            %第一个个体的概率
                        if (prob<route(j))
                            Individual(i,:)=TempIndividual(j,:);
                            break;
                        else
                            continue;
                        end    %第一个节点判断完成
                    else          %不是第一个
                        if (route(j-1)<prob)&&(prob<route(j))
                            Individual(i,:)=TempIndividual(j,:);
                            break;
                        else
                            continue;
                        end                   
                    end      % end of  j==1
                end      %end of j
            end         %end of i   
            % end of case 1  轮盘赌法复制结束          
        case 2
             [B,Jorder] = sort(J(:,k));  %次序是从小到大排列
              for i=1:inum
                  Individual(i,:)=TempIndividual(Jorder(i),:);        
              end         %end of i    
              % end of case 2 最优个体保存法复制 结束
    end   % end of switch  复制结束
    %交叉
    ptemp=rand;
    if (ptemp<pcross)            %产生交叉
        %掩码交叉
        rand1=randperm(inum);   %产生随机数,然后每两个两两交叉.
        for i=1:inum/2
            p1=rand1(2*(i-1)+1);
            p2=rand1(2*(i-1)+2);
            %产生父母
            parent1=Individual(p1,:);
            parent2=Individual(p2,:);
            %产生个体
              %产生掩码
            maskcode=rand(inum,16);    %产生掩码
            maskcode=round(maskcode);   %变成01 串        
            for j=1:col
                if (maskcode(j)==1)   %是1则父母的位交换
                    Child1(j)=parent2(j);     %子1
                    Child2(j)=parent1(j);         %子2
                else   %是0不交换
                    Child1(j)=parent1(j);     %子1
                    Child2(j)=parent2(j);         %子2
                end    % end if
            end   % 按位交叉变异   end of j       
        end  %end of i        
        Individual(inum+1,:)=Child1;
        Individual(inum+2,:)=Child2;
    else         %不产生交叉,就随机生成2个子代供选择
        Temp=rand(1,16);       % 表示的个体
        Individual(inum+1,:)=round(Temp);   %变成01 串
        Temp1=rand(1,16);       % 表示的个体
        Individual(inum+2,:)=round(Temp1);   %变成01 串
    end   %end of if 交叉
    %变异
    ptempchange=rand;
    if (ptempchange<pchange)    %发生变异
       rand2=randperm(inum);
       p1=rand2(1);         %取出随机产生的一个个体变异
       randbyte=randperm(16);
       Individual(p1,randbyte)=abs(Individual(p1,randbyte)-1);    % 减一取绝对值,变化0,1
    end  % end of if 变异
end  %迭代结束
%显示结果
figure(\'color\',\'white\');
k=1:count;
plot(k,MeanJ);
title(\'函数值均值随迭代次数的变化\');
figure(\'color\',\'white\');
plot(k,MinJ);
title(\'函数值最小值随迭代次数的变化\');
disp \'最小的函数值\'
min(Xg)
x1=0;
x2=0;
for j=1:8
    ex=8-j;
    x1=x1+MinIndividual(j)*2^ex;   % 计算x1  0--256
    x2=x2+MinIndividual(j+8)*2^ex; % 计算x2  0--256           
end  %end of j
disp \'取得最小值时的变量x1,x2的值\'
 x1=x1/100-1.28   % 变换 x1到 -1.28:1.28
 x2=x2/100-1.28   % 变换 x2到 -1.28:1.28
 %适配度计算
Result=(x1^2+x2^2-0.3*cos(3*pi*x1)-0.4*cos(4*pi*x2)+0.7)
figure(\'color\',\'white\');
plot(k,Xg);
t=strcat(\'最优解随迭代次数的变化,最优解= \',num2str(Xg(count)));
title(t);

figure(\'color\',\'white\');
plot(k,MeanAdap);
title(\'适配值的均值\');

figure(\'color\',\'white\');
plot(k,MinAdap);
title(\'适配值的最小值\');