Python遺傳演算法工具箱的使用（一）求解帶約束的單目標優化

Strong_wind發表於2019-07-31

Python演算法優化

加了個小目錄~方便定位檢視~

前言

前言

網上有很多部落格講解遺傳演算法，但是大都只是“點到即止”，雖然給了一些程式碼實現，但也是“淺嘗輒止”，沒能很好地幫助大家進行擴充套件應用，抑或是進行深入的研究。

這是我的開篇之作~之前沒有寫部落格的習慣，一般是將筆記存本地，但久而久之發現回看不便，而且無法與大家交流和學習。現特此寫下開篇之作，若有疏漏之處，敬請指正，謝謝！

本文對遺傳演算法的原理進行梳理，相關程式碼是基於國內高校學生聯合團隊開源的高效能實用型進化演算法工具箱geatpy來實現，部分教程引用了geatpy的官方文件：

http://geatpy.com/index.php/details/，

geatpy官網：http://www.geatpy.com

若有錯誤之處，懇請同志們指正和諒解，謝謝！

因為是基於geatpy遺傳和進化演算法工具箱，所以下文的程式碼部分在執行前，需要安裝geatpy：

pip install geatpy

安裝時會自動根據系統版本匹配下載安裝對應的版本。這裡就有個小坑：如果最新版Geatpy沒有與當前版本相匹配的包的話，會自動下載舊版的包。而舊版的包在Linux和Mac下均不可用。

安裝好後，在Python中輸出版本檢查是否是最新版（version 2.5.0）：

import geatpy as ea
print(ea.__version__)

下面切入主題：

自然界生物在周而復始的繁衍中，基因的重組、變異等，使其不斷具有新的性狀，以適應複雜多變的環境，從而實現進化。遺傳演算法精簡了這種複雜的遺傳過程而抽象出一套數學模型，用較為簡單的編碼方式來表現複雜的現象，並通過簡化的遺傳過程來實現對複雜搜尋空間的啟發式搜尋，最終能夠在較大的概率下找到全域性最優解，同時與生俱來地支援平行計算。

下圖展示了常規遺傳演算法(左側) 和某種在平行計算下的遺傳演算法(右側) 的流程。

本文只研究經典的遺傳演算法，力求最後能夠解決各種帶約束單目標優化問題，並能夠很輕鬆地進行擴充套件，讓大家不僅學到演算法理論，還能輕鬆地通過“複製貼上”就能夠將相關遺傳演算法程式碼結合到各類不同的現實問題的求解當中。

從上面的遺傳演算法流程圖可以直觀地看出，遺傳演算法是有一套完整的“固定套路”的，我們可以把這個“套路”寫成一個“演算法模板”，即把：種群初始化、計算適應度、選擇、重組、變異、生成新一代、記錄並輸出等等這些基本不需要怎麼變的“套路”寫在一個函式裡面，而經常要變的部分：變數範圍、遺傳演算法引數等寫在這個函式外面，對於要求解的目標函式，由於在遺傳進化的過程中需要進行呼叫目標函式進行計算，因此可以把目標函式、約束條件寫在另一個函式裡面。

另外我們還可以發現，在遺傳演算法的“套路”裡面，執行的“初始化種群”、“選擇”、“重組”、“變異”等等，其實是一個一個的“運算元”，在geatpy工具箱裡，已經提供現行的多種多樣的進化運算元了，因此直接呼叫即可。

Geatpy工具箱提供一個物件導向的進化演算法框架，因此一個完整的遺傳演算法程式就可以寫成這個樣子：

關於演算法框架的詳細介紹可參見：http://geatpy.com/index.php/geatpy%E6%95%99%E7%A8%8B/

下面就來詳細講一下相關的理論和程式碼實現：

正文

一. 基礎術語：

先介紹一下遺傳演算法的幾個基礎的術語，分別是：”個體“、”種群“、”編碼與解碼“、”目標函式值“、”適應度值“。

1.個體：“個體”其實是一個抽象的概念，與之有關的術語有：

（1）個體染色體：即對決策變數編碼後得到的行向量。

比如：有兩個決策變數x1=1,x2=2，各自用3位的二進位制串去表示的話，寫成染色體就是：

（2）個體表現型：即對個體染色體進行解碼後，得到的直接指代各個控制變數的值的行向量。

比如對上面的染色體“0 0 1 0 1 0”進行解碼，得到 “1 2”，它就是個體的表現型，可看出該個體儲存著兩個變數，值分別是1和2。

注意：遺傳演算法中可以進行“實值編碼”，即可以不用二進位制編碼，直接用變數的實際值來作為染色體。這個時候，個體的染色體數值上是等於個體的表現型的。

（3）染色體區域描述器：用於規定染色體每一位元素範圍，詳細描述見下文。

2. 種群：“種群”也是一個抽象的概念，與之有關的術語有：

（1）種群染色體矩陣(Chrom)：它每一行對應一個個體的染色體。此時會發出疑問：一個個體可以有多條染色體嗎？答：可以有，但一般沒有必要，一條染色體就可以儲存很多決策變數的資訊了，如果要用到多條染色體，可以用兩個種群來表示。

例如：

它表示有3個個體（因為有3行），因此有3條染色體。此時需要注意：這些染色體代表決策變數的什麼值，我們是不知道的，我們也不知道該種群的染色體採用的是什麼編碼。染色體具體代表了什麼，取決於我們採用什麼方式去解碼。假如我們採用的是二進位制的解碼方式，並約定上述的種群染色體矩陣中前3列代表第一個決策變數，後3列代表第二個決策變數，那麼，該種群染色體就可以解碼成：

（2）種群表現型矩陣(Phen)：它每一行對應一個個體的表現型。比如上圖就是根據Chrom種群染色體矩陣解碼得到的種群表現型矩陣。同樣地，當種群染色體採用的是“實值編碼”時，種群染色體矩陣與表現型矩陣實際上是一樣的。

（3）種群個體違反約束程度矩陣(CV)：它每一行對應一個個體，每一列對應一種約束條件（可以是等式約束或不等式約束）。CV矩陣中元素小於或等於0表示對應個體滿足對應的約束條件，大於0則表示不滿足，且越大表示違反約束條件的程度越高。比如有兩個約束條件：

如何計算CV矩陣？可以建立兩個列向量CV1和CV2，然後把它們左右拼合而成一個CV矩陣。

假設x1、x2、x3均為儲存著種群所有個體的決策變數值的列向量（這裡可以利用種群表現型矩陣Phen得到，比如x1=Phen[:, [0]]；x2=Phen[:, [1]])；x3=Phen[:, [2]])，這樣就可以得到種群所有個體對應的x1、x2和x3）。

那麼：

比如在某一代中，種群表現型矩陣Phen為：

則有：

此時CV矩陣的值為：

由此可見，第一個個體滿足兩個約束條件；第二個個體違反了2個約束條件；第三和第四個個體滿足第一個約束條件但違反了第二個約束條件。

下面看下如何用程式碼來生成一個種群染色體矩陣：

程式碼1. 實整數值種群染色體矩陣的建立：

import numpy as np
from geatpy import crtpc
help(crtpc) # 檢視幫助
# 定義種群規模（個體數目）
Nind = 4
Encoding = 'RI' # 表示採用“實整數編碼”，即變數可以是連續的也可以是離散的
# 建立“區域描述器”，表明有4個決策變數，範圍分別是[-3.1, 4.2], [-2, 2],[0, 1],[3, 3]，
# FieldDR第三行[0,0,1,1]表示前兩個決策變數是連續型的，後兩個變數是離散型的
FieldDR=np.array([[-3.1, -2, 0, 3],
                  [ 4.2,  2, 1, 5],
                  [ 0,    0, 1, 1]])
# 呼叫crtri函式建立實數值種群
Chrom=crtpc(Encoding, Nind, FieldDR)
print(Chrom)

程式碼1的執行結果：

這裡要插入講一下“區域描述器”（見程式碼1中的FieldDR），它是用於描述種群染色體所表示的決策變數的一些資訊，如變數範圍、連續/離散性。另外還有一種區域描述器(FIeldD)，用於描述二進位制/格雷碼的種群染色體。FieldDR和FieldD兩個合稱“Field”，又可以認為它們是“譯碼矩陣”。FieldDR具有以下的資料結構：

程式碼1中的FieldDR矩陣的第三行即為這裡的varTypes。它如果為0，表示對應的決策變數是連續型的變數；為1表示對應的是離散型變數。

另一種則是用於二進位制/格雷編碼種群染色體解碼的譯碼矩陣FieldD，它是具有以下結構的矩陣：

其中，lens, lb, ub, codes, scales, lbin, ubin, varTypes均為長度等於決策變數個數的行向量。

lens 包含染色體的每個子染色體的長度。sum(lens) 等於染色體長度。

lb 和ub 分別代表每個決策變數的上界和下界。

codes 指明染色體子串用的是二進位制編碼還是格雷編碼。codes[i] = 0 表示第i 個變數使用的是標準二進位制編碼；codes[i] = 1 表示使用格雷編碼。

scales 指明每個子串用的是算術刻度還是對數刻度。scales[i] = 0 為算術刻度，scales[i] = 1 為對數刻度。對數刻度可以用於變數的範圍較大而且不確定的情況，對於大範圍的引數邊界，對數刻度讓搜尋可用較少的位數，從而減少了遺傳演算法的計算量。（注意：當某個變數是對數刻度時，其取值範圍中不能有0，即要麼上下界都大於0，要麼上下界都小於0。）

從2.5.0版本開始，取消了對對數刻度的支援，該引數暫時保留，但不在起作用。

lbin 和ubin 指明瞭變數是否包含其範圍的邊界。0 表示不包含邊界；1 表示包含邊界。

varTypes 指明瞭決策變數的型別，元素為0 表示對應位置的決策變數是連續型變數；1 表示對應的是離散型變數。

對於二進位制編碼，二進位制種群的染色體具體代表決策變數的什麼含義是不由染色體本身決定的，而是由解碼方式決定的。因此在建立二進位制種群染色體之初就要設定好譯碼矩陣（又稱“區域描述器”）。

因此，可以通過以下程式碼構建一個二進位制種群染色體矩陣：

程式碼2. 二進位制種群染色體矩陣的建立：

import numpy as np
from geatpy import crtpc
help(crtpc) # 檢視幫助
# 定義種群規模（個體數目）
Nind = 4
Encoding = 'BG' # 表示採用“實整數編碼”，即變數可以是連續的也可以是離散的
# 建立“譯碼矩陣”
FieldD = np.array([[3, 2], # 各決策變數編碼後所佔二進位制位數，此時染色體長度為3+2=5
                   [0, 0], # 各決策變數的範圍下界
                   [7, 3], # 各決策變數的範圍上界
                   [0, 0], # 各決策變數採用什麼編碼方式(0為二進位制編碼，1為格雷編碼)
                   [0, 0], # 各決策變數是否採用對數刻度(0為採用算術刻度)
                   [1, 1], # 各決策變數的範圍是否包含下界(對bs2int實際無效，詳見help(bs2int))
                   [1, 1], # 各決策變數的範圍是否包含上界(對bs2int實際無效)
                   [0, 0]])# 表示兩個決策變數都是連續型變數（0為連續1為離散）
# 呼叫crtpc函式來根據編碼方式和譯碼矩陣來建立種群染色體矩陣
Chrom=crtpc(Encoding, Nind, FieldD)
print(Chrom)

程式碼2執行結果：

3. 編碼與解碼

對於實整數編碼（即上面程式碼1所建立的實整數種群染色體），它是不需要解碼，染色體直接就對應著它所代表的決策變數值。而對於程式碼2生成的二進位制種群染色體矩陣，它需要根據譯碼矩陣FieldD來進行解碼。在程式碼2後面新增以下語句即可解碼：

程式碼3（上接程式碼2）：

from geatpy import bs2ri
help(bs2ri)
Phen = bs2ri(Chrom, FieldD)
print('表現型矩陣 = \n', Phen)

執行效果如下：

4.目標函式值：

種群的目標函式值存在一個矩陣裡面(一般命名為ObjV)，它每一行對應一個個體的目標函式值。對於單目標而言，這個目標函式值矩陣只有1列，而對於多目標而言，就有多列了，比如下面就是一個含兩個目標的種群目標函式值矩陣：

（這裡Nind表示種群的規模，即種群含多少個個體；Nvar表示決策變數的個數）

下面來跑一下程式碼，接著程式碼3，在對二進位制染色體解碼成整數值種群後，我們希望計算出f(x,y)=x+y這個目標函式值。同時設定一個等式約束：要求x + y = 3。於是完整程式碼如下：

程式碼4：

import numpy as np
from geatpy import crtpc
from geatpy import bs2ri

def aim(Phen):
    x = Phen[:, [0]] # 取出種群表現型矩陣的第一列，得到所有個體的決策變數x
    y = Phen[:, [1]] # 取出種群表現型矩陣的第二列，得到所有個體的決策變數y
    CV = np.abs(x + y - 3) # 生成種群個體違反約束程度矩陣CV，以處理等式約束：x + y == 3
    f = x + y # 計算目標函式值
    return f, CV # 返回目標函式值矩陣

# 定義種群規模（個體數目）
Nind = 4
Encoding = 'BG' # 表示採用“實整數編碼”，即變數可以是連續的也可以是離散的
# 建立“譯碼矩陣”
FieldD = np.array([[3, 2], # 各決策變數編碼後所佔二進位制位數，此時染色體長度為3+2=5
                   [0, 0], # 各決策變數的範圍下界
                   [7, 3], # 各決策變數的範圍上界
                   [0, 0], # 各決策變數採用什麼編碼方式(0為二進位制編碼，1為格雷編碼)
                   [0, 0], # 各決策變數是否採用對數刻度(0為採用算術刻度)
                   [1, 1], # 各決策變數的範圍是否包含下界(對bs2int實際無效，詳見help(bs2int))
                   [1, 1], # 各決策變數的範圍是否包含上界(對bs2int實際無效)
                   [0, 0]])# 表示兩個決策變數都是連續型變數（0為連續1為離散）
# 呼叫crtpc函式來根據編碼方式和譯碼矩陣來建立種群染色體矩陣
Chrom=crtpc(Encoding, Nind, FieldD)
print('二進位制染色體矩陣 = \n', Chrom)

# 解碼
Phen = bs2ri(Chrom, FieldD)
print('表現型矩陣 = \n', Phen)

# 計算目標函式值矩陣
ObjV, CV = aim(Phen)
print('目標函式值矩陣 = \n', ObjV)
print('CV矩陣 = \n', CV)

執行結果如下：

由上面對CV矩陣的描述可知，第三個個體的CV值為0，表示第三個個體滿足x+y=3這個等式約束。其他都大於0，表示不滿足該約束。

疑問：CV矩陣有什麼用呢？

答：CV矩陣既可用於標記非可行解，在含約束條件的優化問題中有用，又可用於度量種群個體違反各個約束條件的程度的高低。對於含約束條件的優化問題，我們可以採用罰函式或者是可行性法則來進行處理。罰函式法這裡就不展開贅述了，最簡單的罰函式可以是直接找到非可行解個體的索引，然後修改其對應的ObjV的目標函式值即可。

對於可行性法則，它需要計算每個個體違反約束的程度，並把結果儲存在種群類的CV矩陣中。CV矩陣的每一行對應一個個體、每一列對應一個約束條件（可以是等式約束也可以是不等式約束），CV矩陣中元素小於或等於0表示對應個體滿足對應的約束條件，否則是違反對應的約束條件，大於0的值越大，表示違反約束的程度越高。生成CV標記之後，在後面呼叫適應度函式計算適應度時，只要把CV矩陣作為函式傳入引數傳進函式體，就會自動根據CV矩陣所描述的種群個體違反約束程度來計算出合適的種群個體適應度。

5.適應度值：

適應度值通俗來說就是對種群個體的”生存能力的評價“。對於簡單的單目標優化，我們可以簡單地把目標函式值直接當作是適應度值（注意：當用geatpy遺傳和進化演算法工具箱時，則需要對目標函式值加個負號才能簡單地把它當作適應度值，因為geatpy遵循的是”目標函式值越小，適應度值越大“的約定。）

對於多目標優化，則需要根據“非支配排序”或是其他方法來確定種群個體的適應度值，本文不對其展開敘述。

種群適應度(FitnV)：它是一個列向量，每一行代表一個個體的適應度值：

（這裡Nind表示種群的規模，即種群含多少個個體）

geatpy遺傳和進化演算法工具箱裡面有幾個函式可以計算種群個體的適應度，分別是：

ranking、indexing、powing、scaling。具體的用法，可以用help命令檢視，如help(ranking)。

下面接著程式碼4，利用ranking（基於目標函式值排序的適應度分配）計算種群的適應度：

程式碼5（接著程式碼4）：

from geatpy import ranking
help(ranking)
FitnV = ranking(ObjV, CV)
print('種群適應度矩陣 = \n', FitnV)

執行結果：

分析這個結果我們發現，由於第1、2、4個體違反約束條件，而第三個個體滿足約束條件，因此第3個個體的適應度最高。而在第1、2、4個體中，個體1的目標函式值最大，因此適應度最低。可見遵循“最小化目標”的約定，即目標函式值越小，適應度越大。

好了，基本的術語和用法講完後，下面講一下遺傳演算法的基本運算元：

二. 遺傳演算法基本運算元：

我們不用破費時間去寫底層的遺傳運算元，因為geatpy工具箱提供豐富的進化運算元，以下所列運算元不僅限於遺傳運算元：

(如果是在iPython 控制檯中呼叫視覺化繪圖函式（例如使用Spyder 開發工具），一般影像會預設顯示在控制檯或者是開發工具中。此時可以在iPython控制檯下執行%matplotlib 來設定把影像顯示在一個獨立視窗中。)

對於多目標優化，Geatpy中內建以下運算元：

可以用help(運算元名)來檢視對應的API文件，檢視更詳細的用法和例子。

下面講一下理論：

1.選擇：

在進化演算法中存在兩個階段的選擇。第一次是參與進化操作的個體的選擇。這個階段的選擇可以是基於個體適應度的、也可以是完全隨機地選擇交配個體。一旦個體被選中，那麼它們就會參與交叉、變異等進化操作。未被選中的個體不會參與到進化操作中。

第二次是常被稱為“重插入”或“環境選擇”的選擇，它是指在個體經過交叉、變異等進化操作所形成的子代（或稱“育種個體”）後用某種方法來保留到下一代從而形成新一代種群的過程。這個選擇過程對應的是生物學中的” 自然選擇”。它可以是顯性地根據適應度（再次注意：適應度並不等價於目標函式值）來進行選擇的，也可以是隱性地根據適應度（即不刻意去計算個體適應度）來選擇。例如在多目標優化的NSGA-II 演算法中，父代與子代合併後，處於帕累託分層中第一層級的個體以及處於臨界層中的
且擁擠距離最大的若干個個體被保留到下一代。這個過程就沒有顯性地去計算每個個體的適應度。

經典的選擇運算元有：“輪盤賭選擇”、“隨機抽樣選擇”、“錦標賽選擇”、“本地選擇”、“截斷選擇”、“一對一生存者競爭選擇”等等，這裡不展開敘述了，可以參考：

http://geatpy.com/index.php/2019/07/28/%E7%AC%AC%E5%9B%9B%E7%AB%A0%EF%BC%9A%E9%80%89%E6%8B%A9/

這裡要注意：遺傳演算法選擇出的後代是可以有重複的。

下面以低階選擇函式：錦標賽選擇運算元(tour)為例，使用help(tour)檢視其API，得到：

實戰演練如下：

程式碼6：

import numpy as np
from geatpy import tour

help(tour)
FitnV = np.array([[1.2],[0.8],[2.1], [3.2],[0.6],[2.2],[1.7],[0.2]])
chooseIdx = tour(FitnV, 6)
print('個體的適應度為：\n', FitnV)
print('選擇出的個體的下標為：\n', chooseIdx)

執行結果：

光這樣還不夠，這裡只是得出了選擇個體的下標，如果我們需要得到被選中個體的染色體，同時嘗試改用高階選擇函式“selecting”來呼叫低階選擇運算元“tour”來進行選擇，則可以如下操作：

程式碼7：

import numpy as np
from geatpy import selecting

help(selecting)
Chrom=np.array([[1,11,21],
[2,12,22],
[3,13,23],
[4,14,24],
[5,15,25],
[6,16,26],
[7,17,27],
[8,18,28]])

FitnV = np.array([[1.2],[0.8],[2.1], [3.2],[0.6],[2.2],[1.7],[0.2]])
SelCh = Chrom[selecting('tour', FitnV, 6), :] # 使用'tour'錦標賽選擇運算元，同時片取Chrom得到所選擇個體的染色體
print('個體的適應度為：\n', FitnV)
print('選擇後得到的種群染色矩陣為：\n', SelCh)

程式碼7執行結果如下：

將程式碼7中的'tour'換成工具箱中的其他選擇運算元的名稱（如etour, rws, sus），就可以使用相應的選擇運算元進行選擇。

2.重組：

在很多的國內教材、部落格文章、論文中，只提到“交叉”的概念。事實上，遺傳演算法的“交叉”是屬於“重組”運算元裡面的。因為交叉運算元經常使用，而對於“離散重組”、“中間重組”、“線性重組”等等，因為用的很少，所以我們常常只談“交叉”運算元了。交叉運算元實際上是“值互換重組”(Values exchanged recombination)。這裡就不展開敘述了，可以參考：

http://geatpy.com/index.php/2019/07/28/%E7%AC%AC%E4%BA%94%E7%AB%A0%EF%BC%9A%E9%87%8D%E7%BB%84/

與重組有關的遺傳演算法引數是“重組概率”（對於交叉而言就是“交叉概率”）(Pc)，它是指兩個個體的染色體之間發生重組的概率。

下面以兩點交叉(xovdp)為例，類似上面的“tour”那樣我們使用help(xovdp)檢視其API：

程式碼實戰演練如下：

程式碼8：

from geatpy import xovdp
import numpy as np

help(xovdp)
OldChrom=np.array([[1,1,1,1,1],[1,1,1,1,1],[0,0,0,0,0],[0,0,0,0,0]]) #建立一個二進位制種群染色體矩陣
print('交叉前種群染色矩陣為：\n', OldChrom)
NewChrom = xovdp(OldChrom, 1) # 設交叉概率為1
print('交叉後種群染色矩陣為：\n', NewChrom)

程式碼8執行結果如下：

由此可見，xovdp是將前一半個體和後一半個體進行配對交叉的，有人認為應該隨機選擇交叉個體。事實上，在遺傳演算法進化過程中，各個位置的個體是什麼，本身是隨機的，不必要在交叉這裡再增添一個隨機（當然，可以在執行xovdp兩點交叉之前，將種群染色體矩陣按行打亂順序然後再交叉，從而滿足自身需求）。

3.變異：

下面以均勻突變(mutuni)為例，編寫程式碼實現實數值編碼的種群染色體的均勻突變，使用help(mutuni)檢視API文件。

程式碼9：

from mutuni import mutuni
import numpy as np
# 自定義種群染色體矩陣，表示有3個個體，且染色體元素直接表示變數的值（即實值編碼）
OldChrom = np.array([[9,10],
                     [10,10],
                     [10,10]])
# 建立區域描述器（又稱譯碼矩陣）
FieldDR = np.array([[7,8],
                    [10,10],
                    [1, 1]])
# 此處設編碼方式為實值編碼中的“實整數編碼”RI，表示染色體可代表實數和整數
NewChrom = mutuni('RI', OldChrom, FieldDR, 1)
print(NewChrom)

程式碼9的執行結果如下：

4.重插入：

經典遺傳演算法通過選擇、重組和變異後，我們得到的是育種後代，此時育種後代的個體數有可能會跟父代種群的個體數不相同。這時，為了保持種群的規模，這些育種後代可以重新插入到父代中，替換父代種群的一部分個體，或者丟棄一部分育種個體，最終形成新一代種群。前面曾提到過“重插入”也是一種選擇，但它是環境選擇，是用於生成新一代種群的；而前面在交叉變異之前的選擇是用於選擇個體參與交叉變異，那個選擇常被稱作“抽樣”。

現考慮使用精英個體保留的遺傳演算法“EGA”，則重插入操作如下：

程式碼10：

from mutuni import mutuni
import numpy as np
# 自定義父代種群染色體（僅作為例子）：
Chrom = np.array([[1.1, 1.3],
                  [2.4, 1.2],
                  [3,   2.1],
                  [4,   3.1]])
# 若父代個體的適應度為：
FitnV = np.array([[1],
                 [2],
                 [3],
                 [4]])
# 考慮採用“精英保留策略”的遺傳演算法，此時從父代選擇出4-1=3個個體，經過交叉變異後假設子代的染色體為：
offChrom = np.array([[2.1, 2.3],
                     [2.3, 2.2],
                     [3.4, 1.1]])
# 假設直接把目標函式值當作適應度，且認為適應度越大越好。則通過以下程式碼重插入生成新一代種群：
bestIdx = np.argmax(FitnV) # 得到父代精英個體的索引
NewChrom = np.vstack([Chrom[bestIdx, :], offChrom]) # 得到新一代種群的染色體矩陣
print('新一代種群的染色體矩陣為：\n', NewChrom)

在“EGA”中，假設父代種群規模為N，則選擇出(N-1)個個體進行交叉變異，然後找出父代的精英個體，用著個精英個體和交叉變異得到的子代個體進行“拼合”，得到新一代種群。

除了這種重插入生成新一代種群的方法外，還有“父子兩代個體合併選擇”等更優秀的生成新一代種群的方法，這裡就不一一贅述了。

講完上面的基本術語以及遺傳演算法基本運算元後，我們就可以來利用遺傳演算法的“套路”編寫一個遺傳演算法求解問題的程式了：

三.完整實現遺傳演算法：

上文提到遺傳演算法程式可以寫成這個樣子：

其核心是演算法模板類。在遺傳演算法模板裡，我們根據遺傳演算法的“套路”，進行：初始化種群、目標函式值計算、適應度評價、選擇、重組、變異、記錄各代最優個體等操作。geatpy工具箱內建有開源的“套路模板”，原始碼參見：

https://github.com/geatpy-dev/geatpy/tree/master/geatpy/source-code/templets

這些內建演算法模板有詳細的輸入輸出引數說明，以及遺傳演算法整個流程的完整註釋，它們可以應對簡單或複雜的、單目標優化的、多目標優化的、約束優化的、組合優化的等等的問題。

但為了學習，我這裡先不採用框架，直接利用工具箱提供的庫函式來寫一個帶精英個體保留的遺傳演算法。這樣程式碼量比較大，但有利於入門。

編寫程式碼 11、12，分別放在同一個資料夾下：

程式碼11（目標函式aimfuc.py）（這裡要回顧一下前面，Phen是種群表現型矩陣，儲存的是種群所有個體的表現型，而不是單個個體。因而計算得到的目標函式值矩陣也是包含所有個體的目標函式值）：

# -*- coding: utf-8 -*-
"""
aimfunc.py - 目標函式檔案
描述:
    目標：max f = 21.5 + x1 * np.sin(4 * np.pi * x1) + x2 * np.sin(20 * np.pi * x2)
    約束條件：
        x1 != 10
        x2 != 5
        x1 ∈ [-3, 12.1] # 變數範圍是寫在遺傳演算法的引數設定裡面
        x2 ∈ [4.1, 5.8]
"""

import numpy as np

def aimfunc(Phen, CV):
    x1 = Phen[:, [0]] # 獲取表現型矩陣的第一列，得到所有個體的x1的值
    x2 = Phen[:, [1]]
    f = 21.5 + x1 * np.sin(4 * np.pi * x1) + x2 * np.sin(20 * np.pi * x2)
    exIdx1 = np.where(x1 == 10)[0] # 因為約束條件之一是x1不能為10，這裡把x1等於10的個體找到
    exIdx2 = np.where(x2 == 5)[0]
    CV[exIdx1] = 1
    CV[exIdx2] = 1
    return [f, CV]

然後是編寫演算法：

程式碼12：

# -*- coding: utf-8 -*-
"""main.py"""
import numpy as np
import geatpy as ea # 匯入geatpy庫
from aimfunc import aimfunc # 匯入自定義的目標函式
import time

"""============================變數設定============================"""
x1 = [-3, 12.1] # 第一個決策變數範圍
x2 = [4.1, 5.8] # 第二個決策變數範圍
b1 = [1, 1] # 第一個決策變數邊界，1表示包含範圍的邊界，0表示不包含
b2 = [1, 1] # 第二個決策變數邊界，1表示包含範圍的邊界，0表示不包含
ranges=np.vstack([x1, x2]).T # 生成自變數的範圍矩陣，使得第一行為所有決策變數的下界，第二行為上界
borders=np.vstack([b1, b2]).T # 生成自變數的邊界矩陣
varTypes = np.array([0, 0]) # 決策變數的型別，0表示連續，1表示離散
"""==========================染色體編碼設定========================="""
Encoding = 'BG' # 'BG'表示採用二進位制/格雷編碼
codes = [0, 0] # 決策變數的編碼方式，設定兩個0表示兩個決策變數均使用二進位制編碼
precisions =[4, 4] # 決策變數的編碼精度，表示二進位制編碼串解碼後能表示的決策變數的精度可達到小數點後6位
scales = [0, 0] # 0表示採用算術刻度，1表示採用對數刻度
FieldD = ea.crtfld(Encoding,varTypes,ranges,borders,precisions,codes,scales) # 呼叫函式建立譯碼矩陣
"""=========================遺傳演算法引數設定========================"""
NIND      = 100; # 種群個體數目
MAXGEN    = 200; # 最大遺傳代數
maxormins = [-1] # 列表元素為1則表示對應的目標函式是最小化，元素為-1則表示對應的目標函式是最大化
maxormins = np.array(maxormins) # 轉化為Numpy array行向量
selectStyle = 'rws' # 採用輪盤賭選擇
recStyle  = 'xovdp' # 採用兩點交叉
mutStyle  = 'mutbin' # 採用二進位制染色體的變異運算元
Lind = int(np.sum(FieldD[0, :])) # 計算染色體長度
pc        = 0.7 # 交叉概率
pm        = 1/Lind  # 變異概率
obj_trace = np.zeros((MAXGEN, 2)) # 定義目標函式值記錄器
var_trace = np.zeros((MAXGEN, Lind)) # 染色體記錄器，記錄歷代最優個體的染色體
"""=========================開始遺傳演算法進化========================"""
start_time = time.time() # 開始計時
Chrom = ea.crtpc(Encoding, NIND, FieldD) # 生成種群染色體矩陣
variable = ea.bs2real(Chrom, FieldD) # 對初始種群進行解碼
CV = np.zeros((NIND, 1)) # 初始化一個CV矩陣（此時因為未確定個體是否滿足約束條件，因此初始化元素為0，暫認為所有個體是可行解個體）
ObjV, CV = aimfunc(variable, CV) # 計算初始種群個體的目標函式值
FitnV = ea.ranking(ObjV, CV, maxormins) # 根據目標函式大小分配適應度值
best_ind = np.argmax(FitnV) # 計算當代最優個體的序號
# 開始進化
for gen in range(MAXGEN):
    SelCh = Chrom[ea.selecting(selectStyle,FitnV,NIND-1),:] # 選擇
    SelCh = ea.recombin(recStyle, SelCh, pc) # 重組
    SelCh = ea.mutate(mutStyle, Encoding, SelCh, pm) # 變異
    # 把父代精英個體與子代的染色體進行合併，得到新一代種群
    Chrom = np.vstack([Chrom[best_ind, :], SelCh])
    Phen = ea.bs2real(Chrom, FieldD) # 對種群進行解碼(二進位制轉十進位制)
    ObjV, CV = aimfunc(Phen, CV) # 求種群個體的目標函式值
    FitnV = ea.ranking(ObjV, CV, maxormins) # 根據目標函式大小分配適應度值
    # 記錄
    best_ind = np.argmax(FitnV) # 計算當代最優個體的序號
    obj_trace[gen,0]=np.sum(ObjV)/ObjV.shape[0] #記錄當代種群的目標函式均值
    obj_trace[gen,1]=ObjV[best_ind] #記錄當代種群最優個體目標函式值
    var_trace[gen,:]=Chrom[best_ind,:] #記錄當代種群最優個體的染色體
# 進化完成
end_time = time.time() # 結束計時
ea.trcplot(obj_trace, [['種群個體平均目標函式值', '種群最優個體目標函式值']]) # 繪製影像
"""============================輸出結果============================"""
best_gen = np.argmax(obj_trace[:, [1]])
print('最優解的目標函式值：', obj_trace[best_gen, 1])
variable = ea.bs2real(var_trace[[best_gen], :], FieldD) # 解碼得到表現型（即對應的決策變數值）
print('最優解的決策變數值為：')
for i in range(variable.shape[1]):
    print('x'+str(i)+'=',variable[0, i])
print('用時：', end_time - start_time, '秒')

執行程式碼12得到如下結果：

終於，我們把遺傳演算法完整地實現了，但擴充套件性還不夠高。下面學習下如何使用Geatpy提供的進化演算法框架來求解上述問題：（關於使用框架來優化的介紹可詳見http://geatpy.com/index.php/geatpy%E6%95%99%E7%A8%8B/）

在這裡我們可以回顧以下在本文開頭提到的採用遺傳演算法的“套路”來程式設計求解問題的基本流程：

其中執行指令碼和問題類是需要編寫的，演算法模板類我直接呼叫Geatpy內建的"soea_EGA_templet"（帶精英個體保留的單目標遺傳演算法模板）。下面開始編寫程式碼：

程式碼13：問題類"MyProblem"的編寫：

# -*- coding: utf-8 -*-
"""
MyProblem.py
該案例展示了一個簡單的連續型決策變數最大化目標的單目標優化問題。
max f = x * np.sin(10 * np.pi * x) + 2.0
s.t.
-1 <= x <= 2
"""

import numpy as np
import geatpy as ea

class MyProblem(ea.Problem): # 繼承Problem父類
    def __init__(self):
        name = 'MyProblem' # 初始化name（函式名稱，可以隨意設定）
        M = 1 # 初始化M（目標維數）
        maxormins = [-1] # 初始化maxormins（目標最小最大化標記列表，1：最小化該目標；-1：最大化該目標）
        Dim = 2 # 初始化Dim（決策變數維數）
        varTypes = [0] * Dim # 初始化varTypes（決策變數的型別，元素為0表示對應的變數是連續的；1表示是離散的）
        lb = [-3, 4.1] # 決策變數下界
        ub = [12.1, 5.8] # 決策變數上界
        lbin = [1] * Dim # 決策變數下邊界
        ubin = [1] * Dim # 決策變數上邊界
        # 呼叫父類構造方法完成例項化
        ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
    
    def aimFunc(self, pop): # 目標函式
        x1 = pop.Phen[:, [0]] # 獲取表現型矩陣的第一列，得到所有個體的x1的值
        x2 = pop.Phen[:, [1]]
        f = 21.5 + x1 * np.sin(4 * np.pi * x1) + x2 * np.sin(20 * np.pi * x2)
        exIdx1 = np.where(x1 == 10)[0] # 因為約束條件之一是x1不能為10，這裡把x1等於10的個體找到
        exIdx2 = np.where(x2 == 5)[0]
        pop.CV = np.zeros((pop.sizes, 2))
        pop.CV[exIdx1, 0] = 1
        pop.CV[exIdx2, 1] = 1
        pop.ObjV = f # 計算目標函式值，賦值給pop種群物件的ObjV屬性

第二步：編寫執行指令碼“main.py”

# -*- coding: utf-8 -*-
"""main.py"""
import numpy as np
import geatpy as ea # import geatpy
from MyProblem import MyProblem # 匯入自定義問題介面

"""==================================例項化問題物件================================"""
problem = MyProblem() # 生成問題物件
"""==================================種群設定================================"""
Encoding = 'BG'       # 編碼方式
NIND = 100            # 種群規模
Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders) # 建立區域描述器
population = ea.Population(Encoding, Field, NIND) # 例項化種群物件（此時種群還沒被初始化，僅僅是完成種群物件的例項化）
"""==================================演算法引數設定================================"""
myAlgorithm = ea.soea_EGA_templet(problem, population) # 例項化一個演算法模板物件
myAlgorithm.MAXGEN = 200 # 最大進化代數
"""=======================呼叫演算法模板進行種群進化=============================="""
[population, obj_trace, var_trace] = myAlgorithm.run() # 執行演算法模板
population.save() # 把最後一代種群的資訊儲存到檔案中
# 輸出結果
best_gen = np.argmax(obj_trace[:, 1]) # 記錄最優種群是在哪一代
best_ObjV = obj_trace[best_gen, 1]
print('最優的目標函式值為：%s'%(best_ObjV))
print('最優的控制變數值為：')
for i in range(var_trace.shape[1]):
    print(var_trace[best_gen, i])
print('有效進化代數：%s'%(obj_trace.shape[0]))
print('最優的一代是第 %s 代'%(best_gen + 1))
print('評價次數：%s'%(myAlgorithm.evalsNum))
print('時間已過 %s 秒'%(myAlgorithm.passTime))

執行"main.py"執行指令碼即可得到以下結果：

程式碼解析：在“main.py”執行指令碼中，一開始需要例項化一個問題物件。然後是種群物件的例項化。在例項化種群物件前，需要設定種群的編碼方式Encoding、種群規模NIND，並且生成區域描述器Field（或稱譯碼矩陣），因為種群類的構造方法中需要至少用到這三個引數（詳見“Population.py”中種群類的構造方法）。

在完成了問題類物件和種群物件的例項化後，將其傳入演算法模板類的構造方法來例項化一個演算法模板物件。這裡我例項化的是單目標優化的EGA演算法（即帶精英個體保留的遺傳演算法）的模板類物件，即程式碼中的"soea_EGA_templet"。裡面的進化演算法具體是如何操作的，可詳見https://github.com/geatpy-dev/geatpy/blob/master/geatpy/templates/soeas/GA/EGA/soea_EGA_templet.py

採用Geatpy提供的進化演算法框架可以既能最大程度地描述清楚所要求解的問題，而且與進化演算法是高度脫耦的，即上面在編寫問題類的時候完全不需要管後面採用什麼演算法、採用什麼樣編碼的種群，只需把問題描述清楚即可。

而且，遺傳演算法有個好處是：目標函式可以寫得相當複雜，可以解決各種複雜的問題，比如神經網路。以BP神經網路為例，可以把神經網路的引數作為決策變數，神經網路的訓練誤差作為目標函式值，只需把上面的例子修改一下就行了。

而且，一般而言我們不需要像我剛剛最開始那樣刻意去手寫進化演算法，可以直接呼叫geatpy內建的演算法模板就可以解決問題了。geatpy工具箱提供這些內建的演算法模板：

應用Geatpy求解數學建模、工業設計、資源排程等實際優化問題的的朋友們可以直接使用這些演算法模板快速解決各種靈活的優化問題。

四.後記：

最後十分感謝由Geatpy團隊提供的高效能實用型遺傳和進化演算法工具箱，它提供開源的進化演算法框架為遺傳演算法求解單目標/多目標優化、約束優化、組合優化等等給出了相當準確和快捷的解決方案。據稱，geatpy的執行效能相當的高，遠高於matlab的遺傳演算法工具箱、以及採用JAVA、matlab或者Python編寫的一些進化優化平臺或框架，比如jMetal、platemo、pymoo、deap等，後面有時間我將進行詳細的效能對比實驗分析，有相關經驗的讀者也可以自行對比效能。而且依我的體驗來看，這是我網上到處找程式碼、找資料學習、碰了無數次壁後所看到的比較易學易用的工具箱了。

最後值得注意的是Geatpy的頂層是物件導向的進化演算法框架，底層是程式導向的進化運算元。下面放一張geatpy的UML類圖、演算法呼叫的層次結構和庫函式呼叫關係圖，以此記錄方便檢視：