GGAP = .8;          %代間隔，即每代產生多少個新個體
INSR = .9;           %插入率，即插入多少個子代
XOVR =  1;         %交叉率
SP = 2;              %選擇壓力
MUTR = 1;          %變異率
MIGR = 0.2;          %子代之間的遷移率
MIGGEN = 20;         %用于遷移的子代數目

TERMEXACT = 1e-4;    %當達到最小值時的終止值

SEL_F = 'sus';       %選擇函數名
XOV_F = 'recdis';    %個體組合函數名
MUT_F = 'mutbga';    %變異函數名

OBJ_F = 'objharv';   %目標函數名

%得到目標函數的邊界
   FieldDR = feval(OBJ_F,[],1);

%根據在目標函數中定義的變量的數目計算SUBPOP和NIND
   NVAR = size(FieldDR,2);        %變量的數目
   SUBPOP = 2 * floor(sqrt(NVAR));   %子代的數目
   NIND = 20 + 5 * floor(NVAR/50);   %每個子代中個體的數目
   MAXGEN = 300 * floor(sqrt(NVAR)); %最大代數
   MUTR = MUTR / NVAR;               %依賴于NVAR的變異率

%得到目標函數的最小值
   GlobalMin = feval(OBJ_F,[],3);

%目標函數圖形輸出的標題
   FigTitle = [feval(OBJ_F,[],2) '   (' int2str(SUBPOP) ':' int2str(MAXGEN) ') '];

%清除Best和存儲矩陣內容
   %初始化存儲最優結果的矩陣
Best = NaN * ones(MAXGEN,3);
Best(:,3) = zeros(size(Best,1),1);
%存儲最優個體的矩陣
IndAll = [];

%生成實值初始種群
   Chrom = crtrp(SUBPOP*NIND,FieldDR);

%計數變量歸零
   gen = 0;
   termopt = 0;

%計算給定種群的目標函數值
   ObjV = feval(OBJ_F,Chrom);
   % 記錄目標函數估計的次數
   Best(gen+1,3) = Best(gen+1,3) + NIND;

%迭代運算直至符合終止條件或達到了最大代數
   while ((gen < MAXGEN) & (termopt == 0)),

   %保存最優和平均的目標函數值以及最優的個體
      [Best(gen+1,1),ix] = min(ObjV);
      Best(gen+1,2) = mean(ObjV);
      IndAll = [IndAll; Chrom(ix,:)];

   %適應度分配
      FitnV = ranking(ObjV,[2 0],SUBPOP);
            
   %從種群中選擇個體
      SelCh = select(SEL_F, Chrom, FitnV, GGAP, SUBPOP);
      
   %重新組合選擇的個體
      SelCh=recombin(XOV_F, SelCh, XOVR, SUBPOP);

   %變異
      SelCh=mutate(MUT_F, SelCh, FieldDR, [MUTR], SUBPOP);

   %計算對應子代的目標函數值
      ObjVOff = feval(OBJ_F,SelCh);
      Best(gen+1,3) = Best(gen+1,3) + size(SelCh,1);

   %在種群仲插入最優的子代以替換最劣的父代
      [Chrom, ObjV] = reins(Chrom, SelCh, SUBPOP, [1 INSR], ObjV, ObjVOff);

      gen=gen+1;
   %檢驗, 如果最優的目標函數值接近于終止值，

%計算兩者的差
ActualMin = abs(min(ObjV) - GlobalMin);

%如果最優的目標函數值小于終止值，
if ((ActualMin < (TERMEXACT * abs(GlobalMin))) | (ActualMin < TERMEXACT))
         termopt = 1;
      end   

   %子代遷移
      if ((termopt ~= 1) & (rem(gen,MIGGEN) == 0))
         [Chrom, ObjV] = migrate(Chrom, SUBPOP, [MIGR, 1, 0], ObjV);
      end
   end

%結果
   %疊加目標函數估計的次數
   Results = cumsum(Best(1:gen,3));
   %函數估計的次數、均值和最優結果
   Results = [Results Best(1:gen,2) Best(1:gen,1)];
   
%繪圖
   plot(IndAll(gen-4:gen,:)');
   xlabel('u(k)');
   ylabel('f(u)');