遺傳算法求數值函數的最值

0. 引言

設有函數：

f(x) = x + 10*sin(5*x) + 7*cos(4*x);

其圖像容易畫出，如下所示：

先要求求該函數的最大值，讀者可能已經有了很多種思路，本文介紹遺傳算法是如何解決此類問題的。

1. 遺傳算法簡介

如果不關心算法的實現細節的話，遺傳算法可以使用如下的流程描述。

這基本是借鑒生物種群的自然演化規律而抽象得到的流程圖。下面分別解釋流程圖中的各個步驟。

編碼

眾所周知的，生物中都有保存其遺傳信息的物質——染色體，染色體中的記錄的遺傳信息——基因決定了一個生物個體的表現性狀，可以簡單地認為一個染色體唯一確定一個生物個體。

染色體的本質是一種蛋白質，其結構非常復雜，用計算機模擬是不現實的。因此遺傳算法簡單地把染色體視為一串二進制位組成的序列（即該序列中的所有元素取自集合{0,1}）。

假設有命名為A的生物，規定其染色體序列長度為8，這時可以取得一個個體：

[0,1,0,0,1,0,1,1]

初始化種群

顯然一個生物種群不可能只包含一個個體，因此可以規定種群的大小，為方便說明原理，這是規定其大小為4。初始種群可以隨機生成。如下：

[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]
[1,0,1,0,0,0,0,1]
[0,0,0,1,0,1,0,1]

評估種群中的個體適應度

生物種群的存續與興衰決定於其對於環境的適應能力。對於使用二進制位序列表示的個體，總可以找出一個函數來評估這個個體對環境的適應能力。可以規定函數：

f(x) = int(x)
# x是輸入的一個個體，如[0,1,0,0,1,0,1,1]
# 返回值是該二進制序列表示的整數

此外還可以規定上述函數的返回值（一個整數）越大，個體的適應度越高。顯然適應度最高的個體是：

[1,1,1,1,1,1,1,1]

選擇

選擇操作用於選出一個種群中存活概率更高的個體，這些個體組成一個新的種群。顯然個體適應度越高的個體存活概率更高。因此可以使用賭輪盤算法選出新種群。賭輪盤算法的原理不是本文的重點，因此略去，其實現代碼見后文。

交叉

交叉操作也是對生物遺傳過程的抽象。設有兩個生物個體，這兩個個體的生殖繁衍過程中，其染色體會進行交叉交換，生成新的染色體，繁衍得來的新的個體攜帶新的染色體。

設有兩個生物個體如下：

[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]

設對這兩個生物個體有一定的概率發生交叉操作，而交叉操作必然發生在二進制位序列中的某一位置。這一位置的選擇也是隨機的。假設從第5位（從1開始計數）開始交叉，其步驟如下：

[1,1,1,1,1,1,0,0]
[1,1,1,1,0,1,1,0]
->
[1,1,1,1,1],[1,0,0]
[1,1,1,1,0],[1,1,0]
->
[1,1,1,1,1],[1,1,0]
[1,1,1,1,0],[1,0,0]
->
[1,1,1,1,1,1,1,0]
[1,1,1,1,0,1,0,0]

變異

生物在自然環境中受環境的影響，有可能發生遺傳信息的永久改變，也就是變異。遺傳算法規定，二進制位序列的變異指的是在一定的概率下，隨機選擇該序列中的一位取反。假設在第2位發生了變異，步驟如下：

[0,1,0,0,1,0,1,1]
->
[0],1,[0,0,1,0,1,1]
->
[0],0,[0,0,1,0,1,1]
->
[0,0,0,0,1,0,1,1]

小結

上述的各個步驟按照流程圖所示依次執行，直到一個種群中所有的個體都相同（即收斂），如：

[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]
[1,0,1,0,0,0,0,1]

或者主循環的迭代次數超過規定次數，算法停止，輸出結果。

Matlab入門

本文中的算法使用Matlab描述，因此需要對Matlab有一定的了解。但是幫讀者入門Matlab不是我們的任務，因此推薦以下手冊，可以幫助讀者快速掌握Matlab：
https://ww2.mathworks.cn/help/pdf_doc/matlab/getstart_zh_CN.pdf

Matlab是收費軟件，而且對於運行該軟件的性能有一定要求，因此可以使用開源免費且占用資源較少的Octave作為替代品，其官網如下：
https://www.gnu.org/software/octave/

Octave使用的語法規則與Matlab基本相同，不兼容之處見下：
https://en.wikibooks.org/wiki/MATLAB_Programming/Differences_between_Octave_and_MATLAB

本文所附的代碼在Octave上運行通過，並未在Matlab上進行測試，如果有需要使用Matlab運行代碼，需要自己檢查兼容性。

算法實現

function retval = main (input1, input2)
  % 主函數，input1,input2,以及retval都沒有實際作用
  clear all;
  close all;
  clc;
  
  NP = 50; % 種群大小
  L = 16; % 染色體序列長度
  Pc = 0.8; % 交叉概率
  Pm = 0.1; % 變異概率
  G = 50; % 主循環迭代次數
  Xs = 10; % 所求函數的定義域的最大值
  Xx = 0; % 所求函數的定義域的最小值
  S = 10; % 什么意思自己悟吧，我只是不想硬編碼成10
  
  pop = initPop(NP,L);
  
  for i = 1:G
    dpop = decode(pop);
    xpop = getX(dpop,L,Xx,Xs);
    fit = fitValue(xpop);
    if all(fit == fit(1)) % 如果收斂，停止主循環
      break;
    endif
    [maxFit,maxIndex] = max(fit);
    maxX(i) = getX(decode(pop(maxIndex,:)),L,Xx,Xs);
    maxY(i) = maxFit;
    pop = select(pop,fit,xpop);    
    pop = crossOver(pop,Pc);    
    pop = mutation(pop,Pm);
  endfor
  
  
  x = Xx:(Xs-Xx)/(2^L-1)/S:Xs;
  y = targetFunc(x);  
  subplot(2,1,1);
  plot(x,y); % 繪制函數圖像
  hold on;
  plot(maxX,maxY,'r*'); % 繪制歷次迭代生成的最大值點
  
  [px,py] = size(maxY);
  x = 1:py;
  y = maxY;  
  subplot(2,1,2);
  plot(x,y); % 繪制歷次生成的最大值點關於迭代次數的折線圖
  hold off;
  
endfunction


function pop = initPop(NP,L)
  pop = round(rand(NP,L));
endfunction

function res = targetFunc(x)
  res = x + 10*sin(5*x) + 7*cos(4*x);
endfunction

function res = decode(pop)
  % 用於將二進制位序列轉化成整數
  [px,py] = size(pop);
  res = ones(px,1);
  for j = 1:py
    res(:,j) = 2.^(py-j)*pop(:,j);
  endfor
  res = sum(res,2);
endfunction

function res = getX(dpop,L,Xx,Xs)
  % 縮放二進制位序列表示成的整數到
  % 定義域中
  res = dpop./(2^L-1).*(Xs-Xx);
endfunction

function res = fitValue(xpop)
  res = targetFunc(xpop);
endfunction

function newpop = select(pop,fit,xpop)
  % 賭輪盤算法的實現
  [maxFit,maxIndex] = max(fit);
  [minFit,minIndex] = min(fit);
  [px,py] = size(pop);
  newpop = ones(px,py);
  fit = (fit-minFit)./(maxFit-minFit);
  fit = fit./sum(fit);
  cumfit = cumsum(fit);
  ms = sort(rand(px,1));
  fiti = 1;
  newi = 1;
  while newi <= px
    if ms(newi) < cumfit(fiti)
      newpop(newi,:) = pop(fiti,:);
      newi = newi + 1;
    else
      fiti = fiti + 1;
    endif
  endwhile
endfunction

function newpop = crossOver(pop,Pc)
  % 交叉操作
  [px,py] = size(pop);
  newpop = ones(size(pop));
  for i = 1:2:px-1
    if rand<Pc
      cpoint = ceil(rand*py);
      newpop(i,:) = [pop(i,1:cpoint),pop(i+1,cpoint+1:py)];
      newpop(i+1,:) = [pop(i+1,1:cpoint),pop(i,cpoint+1:py)];
     else
      newpop(i,:) = pop(i,:);
      newpop(i+1,:) = pop(i+1,:);
    endif
  endfor
endfunction

function newpop = mutation(pop,Pm)
  % 變異操作
  [px,py] = size(pop);
  newpop = ones(px,py);
  for i = 1:px
    if rand<Pm
      mpoint=ceil(rand*py);
      newpop(i,mpoint) = 1 - pop(i,mpoint);
      newpop(i,:) = pop(i,:);
    else
      newpop(i,:) = pop(i,:);
    endif
  endfor
endfunction

繪圖結果如下：