function [BestPop,Trace]=fmaxga(FUN,LB,UB,eranum,popsize,pcross,pmutation,options)
% [BestPop,Trace]=fmaxga(FUN,LB,UB,eranum,popsize,pcross,pmutation) 
% Finds a  maximum of a function of several variables.
% fmaxga solves problems of the form:  
%      max F(X)  subject to:  LB <= X <= UB                            
%  BestPop--------最優(yōu)的群體即為最優(yōu)的染色體群
%  Trace----------最佳染色體所對應(yīng)的目標(biāo)函數(shù)值
%  FUN------------目標(biāo)函數(shù)
%  LB-------------自變量下限
%  UB-------------自變量上限
%  eranum---------種群的代數(shù),取100--1000(默認1000)
%  popsize--------每一代種群的規(guī)模；此可取50--100(默認50)
%  pcross---------交叉的概率,此概率一般取0.5--0.85之間較好(默認0.8)
%  pmutation------變異的概率,該概率一般取0.05-0.2左右較好(默認0.1)
%  options--------1×2矩陣,options(1)=0二進制編碼(默認0),option(1)~=0十進制編碼,option(2)設(shè)定求解精度(默認1e-4)
%                 
% 例如測試Shaffer's F6函數(shù),自變量下限[-100,-100],上限[100,100],當(dāng)x=[0 0]時,MaxF6=1
% 運行得到相當(dāng)好的結(jié)果:自變量為 0.00033379-4.7684e-005 時,最優(yōu)值 1.000000
% 對應(yīng)染色體是：1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1

T1=clock;
if nargin<3, error('FMAXGA requires at least three input arguments'); end
if nargin==3, eranum=1000;popsize=50;pcross=0.8;pmutation=0.1;options=[0 1e-4];end
if nargin==4, popsize=50;pcross=0.8;pmutation=0.1;options=[0 1e-4];end
if nargin==5, pcross=0.8;pmutation=0.1;options=[0 1e-4];end
if nargin==6, pmutation=0.1;options=[0 1e-4];end
if nargin==7, options=[0 1e-4];end
if find((LB-UB)>0)
   error('數(shù)據(jù)輸入錯誤,請重新輸入(LB<UB):');
end
s=sprintf('程序運行需要約%.4f 秒鐘時間,請稍等......',(eranum*popsize*40/(1000*50)));
disp(s);
bounds=[LB;UB]';bits=[];
precision=options(2);%由求解精度確定二進制編碼長度
bits=ceil(log2((bounds(:,2)-bounds(:,1))' ./ precision));%由設(shè)定精度劃分區(qū)間
[Pop]=initpop(popsize,bits);%初始化種群
[m,n]=size(Pop);
pm0=pmutation;
BestPop=zeros(eranum,n);Trace=zeros(eranum,length(bits)+1);%分配初始解空間
i=1;
while i<=eranum
    for j=1:m
        value(j)=feval(FUN(1,:),(b2f(Pop(j,:),bounds,bits)));%計算適應(yīng)度
    end
    [MaxValue,Index]=max(value);
    BestPop(i,:)=Pop(Index,:);
    Trace(i,1)=MaxValue;
    Trace(i,(2:length(bits)+1))=b2f(BestPop(i,:),bounds,bits);
    [selectpop]=SelectChrom(FUN,Pop,bounds,bits);%選擇
    [CrossOverPop]=CrossOver(selectpop,pcross);%交叉
    [NewPop]=Mutation(CrossOverPop,pmutation);%變異
    Pop=NewPop;%更新
    pmutation=pm0+(i^4)*(pcross/2-pm0)/(eranum^4); %隨著種群向前進化，逐步增大變異率
    p(i)=pmutation;
    i=i+1;
end
t=1:eranum;
plot(t,Trace(:,1)');
title('函數(shù)優(yōu)化的遺傳算法');xlabel('進化世代數(shù)(eranum)');ylabel('每一代最優(yōu)適應(yīng)度(maxfitness)');
[MaxFval,I]=max(Trace(:,1));
X=Trace(I,(2:length(bits)+1));
hold on;  plot(I,MaxFval,'*');
text(I+5,MaxFval,['FMAX=' num2str(MaxFval)]);
str1=sprintf('進化到 %d 代 ,自變量為 %s 時,得本次求解的最優(yōu)值 %f\n對應(yīng)染色體是：%s',...
             I,num2str(X),MaxFval,num2str(BestPop(I,:)));
disp(str1);
%figure(2);plot(t,p);%繪制變異值增大過程
T2=clock;
CostTime=T2-T1;
if CostTime(6)<0
    CostTime(6)=CostTime(6)+60; CostTime(5)=CostTime(5)-1;
end
if CostTime(5)<0
    CostTime(5)=CostTime(5)+60;CostTime(4)=CostTime(4)-1;
end  %像這種程序當(dāng)然不考慮運行上小時啦
str2=sprintf('程序運行耗時 %d 小時 %d 分鐘 %.4f 秒',CostTime(4),CostTime(5),CostTime(6));
disp(str2);

%初始化種群，采用二進制編碼
function [pop]=initpop(popsize,bits)
len=sum(bits);
pop(1,:)=zeros(1,len);%The whole zero encoding
for i=2:popsize-1
    pop(i,:)=round(rand(1,len));
end
pop(popsize,:)=ones(1,len);%The whole one encoding

%解碼
function [fval] = b2f(bval,bounds,bits)
% fval   - 表征各變量的十進制數(shù)
% bval   - 表征各變量的二進制編碼串
% bounds - 各變量的取值范圍
% bits   - 各變量的二進制編碼長度
scale=(bounds(:,2)-bounds(:,1))'./(2.^bits-1); %The range of the variables
numV=size(bounds,1);
cs=[0 cumsum(bits)]; 
for i=1:numV
  a=bval((cs(i)+1):cs(i+1));
  fval(i)=sum(2.^(size(a,2)-1:-1:0).*a)*scale(i)+bounds(i,1);
end

%選擇操作
function [selectpop]=SelectChrom(FUN,pop,bounds,bits)%計算各個體的適應(yīng)度并采用輪盤賭進行選擇
[m,n]=size(pop);
for i=1:m
    fit(i)=feval(FUN(1,:),(b2f(pop(i,:),bounds,bits)));%以函數(shù)值為適應(yīng)度
end
selectprob=fit/sum(fit);%選擇概率
prob=cumsum(selectprob);%累計選擇概率
sumprob=[0 prob];
for i=1:m
    selectpop(i,:)=pop(length(find(rand>=sumprob)),:);
end    
    
%交叉操作
function [NewPop]=CrossOver(OldPop,pcross)%OldPop為父代種群，pcross為交叉概率
[m,n]=size(OldPop);
r=rand(1,m);
y1=find(r<pcross);
y2=find(r>=pcross);
len=length(y1);
if len>2&mod(len,2)==1%如果用來進行交叉的染色體的條數(shù)為奇數(shù)，將其調(diào)整為偶數(shù)
    y2(length(y2)+1)=y1(len);
    y1(len)=[];
end
if length(y1)>=2
   for i=0:2:length(y1)-2
       [NewPop(y1(i+1),:),NewPop(y1(i+2),:)]=EqualCrossOver(OldPop(y1(i+1),:),OldPop(y1(i+2),:));
   end     
end
NewPop(y2,:)=OldPop(y2,:);

function [children1,children2]=EqualCrossOver(parent1,parent2)
%采用均勻交叉 例：
%父1：0 1 1 1 0 0 1 1 0 1 0
%父2：1 0 1 0 1 1 0 0 1 0 1
%掩碼：0 1 1 0 0 0 1 1 0 1 0
%交叉后新個體：
%子1：1 1 1 0 1 1 1 1 1 1 1 
%子2：0 0 1 1 0 0 0 0 0 0 0
L=length(parent1);
hidecode=round(rand(1,L));%隨機生成掩碼，如hidecode=[0 1 1 0 0 0 1 1 0 1 0];
children1=zeros(1,L);
children2=zeros(1,L);
children1(find(hidecode==1))=parent1(find(hidecode==1));%掩碼為1，父1為子1提供基因
children1(find(hidecode==0))=parent2(find(hidecode==0));%掩碼為0，父2為子1提供基因
children2(find(hidecode==1))=parent2(find(hidecode==1));%掩碼為1，父2為子2提供基因
children2(find(hidecode==0))=parent1(find(hidecode==0));%掩碼為0，父1為子2提供基因

%變異操作
function [NewPop]=Mutation(OldPop,pmutation)
[m,n]=size(OldPop);
r=rand(1,m);
position=find(r<=pmutation);
len=length(position);
if len>=1
   for i=1:len
       k=unidrnd(n,1,1); %設(shè)置變異點數(shù)，一般設(shè)置1點
       for j=1:length(k)
           if OldPop(position(i),k(j))==1
              OldPop(position(i),k(j))=0;
           else
              OldPop(position(i),k(j))=1;
           end
       end
   end
end
NewPop=OldPop;