%在比較好的結果上繼續改動080319
%防止近親繁殖采取交叉控制
%080310采用自適應調整算子(pc=pc1;pm=pm1;pm_plus1=pm1;)pc\pm all adapted
%080224浮點數編碼遺傳算法
%080217made by wjh
%Generic Algorithm for function f(x1,x2) optimum
%f(x1,x2) optimum is a example from a GA book writen by zhou ming
clc
clear all;
close all;
for t_times=1:20
%Parameters
CS_Size=240;   
%zq_k=1;
G=100;     
CodeL=5;
b_zhishu=3;

G_z=100;
if t_times<=3
   b_zhishu=1;
   bianyi_z=10;
elseif 3<t_times<7
         b_zhishu=3;
          bianyi_z=10;
elseif 7<=t_times<10
     b_zhishu=10;
     bianyi_z=10
else
    b_zhishu=10;
     bianyi_z=100
end


pc1=0.9;pc2=0.5;pm1=0.2;pm2=0.001;
%pm=0.05;pc=0.8;
pc=pc1;pm=pm1;pm_plus1=pm1;

    E=rand(CS_Size,CodeL);    %Initial Code產生初始種群
    edge_CanshuDelta=[6 4 0.2 0.01 4];  %delta_thetab[H r1 b thetab  fai]
        for k=1:CS_Size
        E_begin(k,:)=edge_CanshuDelta*2.*E(k,:)-edge_CanshuDelta;
        end
 
        f=zeros(1,CS_Size);
        for kf1=1:(CS_Size)
        f(kf1)=F_err(E_begin(kf1,:));
        end


%Main Program
for k=1:1:G
   
f_central(k)=sum(f)/CS_Size;
f_min(k)=min(f);
time(k)=k;
n_add=0;
%選擇

  %Arranging fi small to bigger
[Oderfi,Indexfi]=sort(f);
E_begin=[ E_begin(Indexfi(1:60),:);E_begin(Indexfi(1:160),:); E_begin(Indexfi(61:80),:);];
f=[f(1:60),f(1:160),f(61:80)];
    
%--------------------------------------------------------------------------

peidui_r=rand(1,CS_Size);
peidui_hou=rand(1,CS_Size);
[Oderpeidui_qian,Indexpeidui_qian]=sort(peidui_r);          %產生隨機配對表
[Oderpeidui_hou,Indexpeidui_hou]=sort(peidui_hou);   
  for s=1:1:CS_Size/2
      
        f_pd_temp_min=min(f(Indexpeidui_qian(s)),f(Indexpeidui_hou(s)));
        if f_pd_temp_min>f_central(k)
            pc=pc1-(pc1-pc2)*(f_central(k)-f_pd_temp_min)/(f_central(k));%-f_min(k)
        else
            pc=pc1;
        end
        
	if rand(1)<pc
     
            n_add=n_add+1;
            r=rand(1);
           E_begin(CS_Size+n_add,:)= r*E_begin(Indexpeidui_qian(s),:)+(1-r)*E_begin(Indexpeidui_hou(s),:);
           if  E_begin(CS_Size+n_add,:)<-edge_CanshuDelta
               temp=E_begin(CS_Size+n_add,:)<-edge_CanshuDelta;
               E_begin(CS_Size+n_add,:)= E_begin(CS_Size+n_add,:).*(temp)+edge_CanshuDelta.*(1-temp);
           elseif E_begin(CS_Size+n_add,:)>edge_CanshuDelta
               
               temp=E_begin(CS_Size+n_add,:)>edge_CanshuDelta;
               E_begin(CS_Size+n_add,:)= E_begin(CS_Size+n_add,:).*(temp)+edge_CanshuDelta.*(1-temp);
           else
               
           end
            
         
              n_add=n_add+1;
           
              E_begin( CS_Size+n_add,:)=  r*E_begin(Indexpeidui_hou(s),:)+(1-r)*E_begin(Indexpeidui_qian(s),:) ;
              if  E_begin(CS_Size+n_add,:)<-edge_CanshuDelta
               temp=E_begin(CS_Size+n_add,:)<-edge_CanshuDelta;
               E_begin(CS_Size+n_add,:)= E_begin(CS_Size+n_add,:).*(temp)+edge_CanshuDelta.*(1-temp);
           elseif E_begin(CS_Size+n_add,:)>edge_CanshuDelta
               
               temp=E_begin(CS_Size+n_add,:)>edge_CanshuDelta;
               E_begin(CS_Size+n_add,:)= E_begin(CS_Size+n_add,:).*(temp)+edge_CanshuDelta.*(1-temp);
           else
               
              end
    end

    %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
           if f(Indexpeidui_qian(s))>f_central(k)
            pm=pm1-(pm1-pm2)*(f_central(k)-f(Indexpeidui_qian(s)))/(f_central(k));%-f_min(k)
           else
               pm=pm1;
            end
            if f(Indexpeidui_hou(s))>f_central(k)
           pm_plus1=pm1-(pm1-pm2)*(f_central(k)-f(Indexpeidui_hou(s)))/(f_central(k));%-f_min(k)
            else
                 pm_plus1=pm1;
                
            end
     for km=1:CodeL
                        if rand(1)<pm 
                          n_add=n_add+1;   
                if rand(1)<=0.5
                E_begin( CS_Size+n_add,km)=E_begin(s,km)+(edge_CanshuDelta(km)-E_begin(s,km))*rand(1)*(1-k/G_z)^b_zhishu/bianyi_z;% 
                % E_begin(s,km)=E_begin(s,km)+(edge_CanshuDelta(km)-E_begin(s,km))*rand(1);
                else
                 E_begin( CS_Size+n_add,km)=E_begin(s,km)-(E_begin(s,km)-edge_CanshuDelta(km))*rand(1)*(1-k/G_z)^b_zhishu/bianyi_z;% 
                 % E_begin(s,km)=E_begin(s,km)-(E_begin(s,km)-edge_CanshuDelta(km))*rand(1); 
                end
                        end
            
                        if rand(1)<pm_plus1
                           n_add=n_add+1;
                if rand(1)<=0.5
                E_begin( CS_Size+n_add,km)=E_begin(s+CS_Size/2,km)+(edge_CanshuDelta(km)-E_begin(s+CS_Size/2,km))*rand(1)*(1-k/G_z)^b_zhishu/bianyi_z; 
                % E_begin(s,km)=E_begin(s,km)+(edge_CanshuDelta(km)-E_begin(s,km))*rand(1);
                else
                 E_begin( CS_Size+n_add,km)=E_begin(s+CS_Size/2,km)-(E_begin(s+CS_Size/2,km)-edge_CanshuDelta(km))*rand(1)*(1-k/G_z)^b_zhishu/bianyi_z; 
                 % E_begin(s,km)=E_begin(s,km)-(E_begin(s,km)-edge_CanshuDelta(km))*rand(1); 
                end
                        end
     end
              
     
     %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%變異操作

  end
    
              for kf2=1:(CS_Size+n_add)
                f(kf2)=F_err(E_begin(kf2,:));
              end
               
                if (min(f))<1e-7
                        break
                end


  
%****** Step 1 : Evaluate BestJ ******


                         %Fitness Function
[Oderfi,Indexfi]=sort(f);     %Arranging fi small to bigger
Bestfi(k)=Oderfi(1);           %Let Bestfi=max(fi)
E_G(k,:)=E_begin(Indexfi(1),:); 
E_begin=E_begin(Indexfi(1:CS_Size),:);      %Let BestS=E(m), m is the Indexfi belong to max(fi)
f=f(Indexfi(1:CS_Size));
n_add=0;

if mod(k,50)==0
    k
    Bestfi
    E_G
    
end

end
Bestfi_times(t_times,:)=Bestfi;           %Let Bestfi=max(fi)
E_G_times(:,:,t_times)=E_G;
f_central_times(t_times,:)=f_central;

end

save E_G_times E_G_times
save Bestfi_times Bestfi_times   
save  f_central_times  f_central_times