%使用遺傳算法求解
% f(x1,x2)=x1^2+x2^2-0.3*cos(3*pi*x1)-0.4*cos(4*pi*x2)+0.7的極小值,并畫出每一代中個體適配函
% 數(shù)的平均值和最小值隨迭代次數(shù)的變化關(guān)系

%程序總體說明：
% 編碼方式說明： 由于期望所求的函數(shù)極小值包含了平方項，而正弦和余弦項的范圍只能是-1---1，因此優(yōu)化范圍的x1,x2只能在-1--1之間，編碼時
% 首先將x1,x2加1變到0--2之間，之后為了編碼需要，設(shè)采樣精度是0.01 ,將0--2之間的數(shù)乘以100,得到0--200之間的整數(shù)，另外為了編碼的需要，
%設(shè)編碼為8位，編碼成0--255的二進(jìn)制碼。,每個個體用16位的二進(jìn)制串表示，的前8位是x1,后8位是x2
%適配制計算說明：%根據(jù)編碼方式，將每一個個體的編碼變化到-1--1之間，之后按照函數(shù)f計算
%可以選擇  1--直接利用函數(shù)的倒數(shù)作為適配值計算； 2--使用排序法計算， 3--使用1-J(i)/sum(J)方法計算
%復(fù)制操所機制：1--輪盤賭方法  ； 2--保留最優(yōu)淘汰法
%交叉機制: 使用掩碼法進(jìn)行交叉操作.

clear all;
%使用
disp '輸入初始化參數(shù)'
inum=input('請輸入種群大小(回車為缺省值30個個體):    ');
if isempty(inum)
    inum=50;   %種群中的個體數(shù)
end
count=input('請輸入迭代次數(shù)(回車為缺省值100代):     ');
if isempty(count)
    count=100;    %迭代次數(shù)
end
pcross=input('請輸入交叉概率(回車為缺省值0.6):      ');
if isempty(pcross)
    pcross=0.6;   %交叉概率 
end
pchange=input('請輸入變異概率(回車為缺省值0.001):      ');
if isempty(pchange)
    pchange=0.001;     %變異概率
end
disp '程序?qū)崿F(xiàn)機制選擇'
%程序?qū)崿F(xiàn)機制選擇
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為最優(yōu)保存淘汰法,缺省為2):          ');
if (isempty(Copyflag))||(Copyflag>2)
    Copyflag=2;
end

%程序計算開始
constk=5; %用來控制適配制之間的常數(shù)
%初始種群
Individual=rand(inum,16);       % 表示的個體的矩陣
Individual=round(Individual);   %變成01 串
J=zeros(inum,1);                %存儲各個個體的函數(shù)值
Jorder=zeros(inum,1);        %函數(shù)值的大小順序 
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   %每一個個體計算         
    %適配度計算
        %轉(zhuǎn)化成十進(jìn)制數(shù)
        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));
    %得到最優(yōu)解，保存在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);    
    %復(fù)制
    switch Copyflag
        case 1     %1為輪盤賭法,2為最優(yōu)保存淘汰法
            %隨機產(chǎn)生一個數(shù)
            for i=1:inum1
                if (i==1)
                    route(i)=Adap1(i);        %產(chǎn)生輪盤
                else
                    route(i)=route(i-1)+Adap1(i);
                end        
            end
            for i=1:inum         %復(fù)制產(chǎn)生inum個個體 
                prob=rand;
                for j=1:inum1    % 根據(jù)輪盤賭概率決定哪個個體被復(fù)制
                    if (j==1)            %第一個個體的概率
                        if (prob<route(j))
                            Individual(i,:)=TempIndividual(j,:);
                            break;
                        else
                            continue;
                        end    %第一個節(jié)點判斷完成
                    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  輪盤賭法復(fù)制結(jié)束          
        case 2
             [B,Jorder] = sort(J(:,k));  %次序是從小到大排列
              for i=1:inum
                  Individual(i,:)=TempIndividual(Jorder(i),:);        
              end         %end of i    
              % end of case 2 最優(yōu)個體保存法復(fù)制 結(jié)束
    end   % end of switch  復(fù)制結(jié)束 
    %交叉
    ptemp=rand;
    if (ptemp<pcross)            %產(chǎn)生交叉
        %掩碼交叉
        rand1=randperm(inum);   %產(chǎn)生隨機數(shù),然后每兩個兩兩交叉.
        for i=1:inum/2
            p1=rand1(2*(i-1)+1);
            p2=rand1(2*(i-1)+2);
            %產(chǎn)生父母
            parent1=Individual(p1,:);
            parent2=Individual(p2,:);
            %產(chǎn)生個體
              %產(chǎn)生掩碼
            maskcode=rand(inum,16);    %產(chǎn)生掩碼
            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         %不產(chǎn)生交叉,就隨機生成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)    %發(fā)生變異
       rand2=randperm(inum);
       p1=rand2(1);         %取出隨機產(chǎn)生的一個個體變異
       randbyte=randperm(16);
       Individual(p1,randbyte)=abs(Individual(p1,randbyte)-1);    % 減一取絕對值,變化0,1
    end  % end of if 變異
end  %迭代結(jié)束
%顯示結(jié)果
figure('color','white');
k=1:count;
plot(k,MeanJ);
title('函數(shù)值均值隨迭代次數(shù)的變化');
figure('color','white');
plot(k,MinJ);
title('函數(shù)值最小值隨迭代次數(shù)的變化');
disp '最小的函數(shù)值'
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('最優(yōu)解隨迭代次數(shù)的變化,最優(yōu)解= ',num2str(Xg(count)));
title(t);

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

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