Matlab與遺傳算法

文章目錄

• 引入
• 1 步驟描述
• 1.1 編碼
• • 1.2 解碼
  • 1.3 交配
  • 1.4 突變
  • 1.5 倒位
  • 1.6 個體適應度評估
  • 1.7 複製
• 2 示例
• • 2.1 問題1
  • • 2.1.1 代碼
    • 2.1.2 求解步驟

引入

  J. Holland教授於1975年提出Simple Genetic Algorithm (SGA) ¹算法，其由編解碼、個體適應度評估和遺傳運算三大模塊構成；遺傳運算包括染色體複製、交叉、變異及倒位等。
  參考書目：卓金武等，MATLAB在數學建模中的應用。

1 步驟描述

1.1 編碼

  遺傳算法的編碼 (基因型)有浮點編碼和二進制編碼兩種，這裏只對後者進行介紹。設某一參數的取值範圍爲$(L,U)$，使用長度爲$k$的二進制編碼表示該參數，則它共有$2^k$種不一樣的編碼。該參數編碼時的對應關係爲：

$$\begin{array}{l}000000000000000000=0\to L\\ 000000000000000001=1\to L+\delta \\ 000000000000000010=2\to L+2\delta \\ 000000000000000011=3\to L+3\delta \\ \cdots \cdots \\ 1111111111111111111=2^k-1\to U\end{array}$$其中

$$\delta =\frac{U-L}{2^k-1}.$$

1.2 解碼

  解碼 (表現型)的目的是爲了將不直觀的二進制數據還原成十進制。設某一個體的二進制編碼爲$b_k{b}_{k-1}\cdots b_2{b}_1$，則相應解碼公式爲：

$$x=L+(\sum _{i=1}^k{b}_i{2}^{i-1})\frac{U-L}{2^k-1}.$$  例如，設有參數$x\in [2,4]$，現用$5$位二進制數對$x$進行編碼，可得$2^5=32$條染色體：

$$00000,00001,\cdots ,11111$$  對於任意二進制數據，例如$x_{22}=10101$，解碼後對應的$x$值爲：

$$x=2+21\times \frac{4-2}{2^5-1}=3.3548.$$

1.3 交配

  交配運算是使用單點或多點進行交叉的算子。首先用隨機數產生一個或多個交配點位置，而後兩個個體在交配點位置互換部分基因，造成兩個子個體。例以下圖中，父染色體$S_1=01001011$、$S_2=10010101$交換其後四位基因，獲得$S_1^{\prime }=01000101$、$S_2^{\prime }=10011011$兩條子染色體。

1.4 突變

  突變運算是使用基本位進行基因突變。爲了不在算法迭代後期出現種羣過早收斂，對於二進制基因碼組成的個體種羣，實行基因碼的小概率翻轉。例如，$S=11001101$第$3$位翻轉，即$S=11001101\to 11101101=S^{\prime }$。$S^{\prime }$可看做是$S$的子染色體。

1.5 倒位

  倒位運算是指一個染色體某區段正常排列順序發生$180^{\circ }$的顛倒，形成染色體內的DNA從新排列。例如對$S$下劃線的部分進行倒位：$S=11\underline{0011}01\Rightarrow S^{\prime }=11\underline{1100}01$。

1.6 個體適應度評估

  **遺傳算法依照與個體適應度成正比的概率決定種羣中各個個體遺傳到下一代的機會。**個體適應度大的個體更容易被遺傳到下一代。一般，求目標函數最大值的問題能夠直接把目標函數做爲檢測個體適應度大小的函數。

1.7 複製

  複製運算是根據個體適應度大小決定其下代遺傳的可能性。設種羣中個體總數爲$N$，個體$i$的適應度爲$f_i$則個體$i$被選取的概率爲：

$$P_i=\frac{f_i}{{\sum }_{k=1}^N{f}_k}.$$  當個體複製概率決定後，再產生$[0,1]$區間的均勻隨機數來決定哪一個個體參加複製。

2 示例

2.1 問題1

$$\begin{gathered} maxf(x_1,x_2)=21.5+x_1\sin (4\pi x_1)+x_2\sin (20\pi x_2) \\ \text{s.t.}\{\begin{array}{l}-3.0\le x_1\le 12.1;\\ 4.1\le x_2\le 5.8.\end{array} \end{gathered}$$其三維圖以下：

2.1.1 代碼

  這裏給出完整求解代碼，具體的每一步將在代碼後進行介紹。

function f = test()
    %% 參數設置
    para.num = 10; % 初始個體數量，建議<100
    para.num_loop = 1000; % 迭代次數，建議100~500
    para.len_x1 = 18; % x1編碼長度
    para.len_x2 = 15; % x2編碼長度
    para.l_x1 = -3.0; % x1下界
    para.l_x2 = 4.1; % x2下界
    para.u_x1 = 12.1; % x1上界
    para.u_x2 = 5.8; % x2上界
    para.pc = 0.25; % 交配機率，建議0.4~0.99
    para.pm = 0.01; % 基因突變機率，建議0.0001~0.2
    para.num_pc = floor(0.25 * para.num); % 染色體交配數量
    
    %% 參數計算
    para.len = para.len_x1 + para.len_x2;
    para.num_pc = para.num_pc - mod(para.num_pc, 2); % 染色體交配數量默認設置爲偶數
    f.obj = 0;
    f.para = [0, 0];

    %% 步驟1.1：生成隨機個體
    ini.geni = get_init_geni(para.num, para.len);
    
    for loop = 1: para.num_loop
        fprintf('%d-th looping...\n', loop);
       %% 步驟1.2：解碼
        ini.decode = get_decode(ini.geni, para);

        %% 步驟2.1：計算適應度值
        ini.eval = zeros(1, para.num);
        for i = 1: para.num
           ini.eval(i) = get_eval(ini.decode(i, 1), ini.decode(i, 2));
        end
        ini.sum_eval = sum(ini.eval);
        temp_max = max(ini.eval);
		
		%% 保存當前最優解
        f.obj = max(temp_max, f.obj);
        f.para = ini.decode(ini.eval == temp_max, :);

        %% 步驟2.2：計算選取機率
        ini.prob = ini.eval / ini.sum_eval;
        ini.sum_prob = cumsum(ini.prob);
        
        %% 步驟3：新種羣複製
        temp_rand = rand(1, para.num);
        temp_choose = zeros(1, para.num);
        for i = 1: para.num
            temp_idx = find(ini.sum_prob > temp_rand(i));
            temp_choose(i) = temp_idx(1);
        end
        ini.geni = ini.geni(temp_choose, :);
        
        %% 步驟4.1：肯定交配個體
        temp_rand = rand(1, para.num);
        [~, temp_min_idx] = sort(temp_rand);
        temp_min_idx = temp_min_idx(1 : para.num_pc);
        
        %% 步驟4.2：開始交配
        for i = 1: floor(para.num_pc / 2)
           % 隨機產生交配點
           temp_idx = round(rand(1,1) * (para.len - 1)) + 1;
           temp_geni = ini.geni(temp_min_idx(2 * i - 1), temp_idx: para.len);
           ini.geni(temp_min_idx(2 * i - 1), temp_idx: para.len) = ini.geni(temp_min_idx(2 * i), temp_idx: para.len);
           ini.geni(temp_min_idx(2 * i), temp_idx: para.len) = temp_geni;
        end
        
        %% 步驟5：基因突變
        temp_rand = rand(1, para.num * para.len);
        temp_rand = find(temp_rand <= para.pm);
        temp_idx = floor(temp_rand / para.len) + 1;
        temp_idx(temp_idx == para.num + 1) = para.num;
        temp_geni_idx = mod(temp_rand, para.len);
        temp_geni_idx(temp_geni_idx == 0) = para.len;
        
        if size(temp_idx, 2) == 0
            disp(loop);
            continue;
        end
        
        for i = 1: size(temp_idx, 2)
           if ini.geni(temp_idx(i), temp_geni_idx(i)) == 0
              ini.geni(temp_idx(i), temp_geni_idx(i)) = 1;
           else
               ini.geni(temp_idx(i), temp_geni_idx(i)) = 0;
           end
        end
        
    end
end

function f = get_init_geni(num, len)
    %% 生成初始隨機個體
    %   num：初始個體數
    %   len：基因長度
    f = zeros(num, len);
    for i = 1: num
       f(i, :) = rand(1, len) < 0.5; 
    end
end

function f = get_decode(geni, para)
    %% 獲取解碼值
    f = zeros(para.num, 2);
    for i = 1: para.num
        f(i, 1) = decode(geni(i, 1: para.len_x1), para.l_x1, para.u_x1);
        f(i, 2) = decode(geni(i, para.len_x1 + 1: para.len), para.l_x2, para.u_x2);
    end
end

function f = decode(geni, l, u)
    %% 解碼每個個體
    m = size(geni, 2);
    sum_bit = 0;
    for i = 1: m
       sum_bit = sum_bit + geni(i) * 2^(i - 1); 
    end
    f = l + sum_bit * (u - l) / (2^m - 1);
end

function f = get_eval(x1, x2)
    %% 適應度函數，對於極大值問題，適應度函數即優化目標函數
    f = 21.5 + x1 * sin(4 * pi * x1) + x2 * sin(20 * pi * x2);
end

  這裏給出迭代一千次的結果：

obj: 38.7462
    para: [10.6994 4.8154]

  適應度的變化以下圖：

2.1.2 求解步驟

  這裏對詳細的步驟進行介紹，整體的步驟與第一節中步驟、代碼步驟一致 (假設初始個體數爲$10$)。

  1）編 碼 解 碼
  編碼的目的是將變量轉換爲二進制數串。數串的長度$m_j$取決於所要求的精度$p$，即小數點保留位數：

$$2^{m_j-1}<(U-L)\times 10^p\le 2^{m_j}-1$$所以，對於上題的約束條件，取精度爲$4$，有：

$$\{\begin{array}{l}-3.0\le x_1\le 12.1\\ [12.1-(-3.0)]\times 10^4=151000\\ 2^{m_1-1}<151000\le 2^{m_1}-1\\ m_1=18\\ 4.1\le x_1\le 5.8\\ [5.8-4.1]\times 10^4=17000\\ 2^{m_2-1}<17000\le 2^{m_2}-1\\ m_1=15\\ m=m_1+m_2=33\end{array}$$即本例中任意染色體數串的長度均爲$33$，其中前$18$位表示$x_1$，後$15$位表示$x_2$。解碼時，只需明確當前變量的上下界、數串、數串及起始位置。

  2）評 價 個 體 適 應 度
  對於極大值函數，適應度函數爲該函數便可；極小值時，只需目標函數總體取負值。所以，每一個個體 (染色體)的適應度，即將解碼後的變量帶入適應度函數便可。
  2.1）計算每一個個體的適應度$f_1$，須要保留當前最佳適應度即相應參數；
  2.2）計算適應度總和$F$；
  2.3）計算每一個個體被複制的機率：

$$P_k=\frac{f_k}{F}.$$  2.4）計算每一個個體被複制的累加機率：

$$Q_k=\sum _{j=1}^k{P}_k.$$該函數對應matlab中的cumsum函數。

  3）新 種 羣 復 制
  3.1）爲每一個個體生成屬於$0\sim 1$的隨機數，假設十個個體的相關係數以下 (圖片源自引入中所提書籍)：

  每次選擇時，根據當前個體的隨機數進行選擇，例如隨機數$0.300000$，其大於$Q_3$小於$Q_4$，所以$U_4$被選擇。顯然，複製的樣本源自原樣本，且新樣本中可能有重複個體。這樣選擇的依據在於適應度大的個體，其機率區間$[Q_k,Q_{k+1}]$也更大，即更易被選中。

  4）新 種 羣 交 配
  4.1）肯定交配染色體數量，例如交配機率$P_c=0.25$，則交配數量爲2，代碼中爲了簡便，最終的交配數量將強制爲偶數；
  4.2）肯定交配對象，首先爲每一個個體生成[0, 1]區間隨機數，選取符合交配數量的最小數的個體做爲交配對象；
  4.3）交配，首先生成一個介於$0\sim (m-1)$的整數，例如產生的數爲$1$，則將選中個體的$[1:]$區間的數值進行交互，示例以下 (圖片源自引入中所提書籍)：

  5）基 因 突 變
  本例中，基因的總長度爲$10\times 33=330$，假設突變機率爲$P_m=0.01$：首先爲生成330個介於$0\sim 1$的隨機數，並進行編號；獲取小於$P_m$的隨機數的編號，並經過如下方式得到突變基因所在的位置 (圖片源自引入中所提書籍)：

  至此，全新的子代產生，並回退至步驟2，直至達到給定迭代次數。

---

1. J. Holland, Adaptation in Natural and Artificial Systems. ↩︎

本文同步分享在 博客「因吉」（CSDN）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。