Python遺傳演算法工具箱Geatpy的使用（一）求解帶約束的單目標優化
知乎：Python遺傳演算法工具箱Geatpy的使用（一）求解帶約束的單目標優化
前言
https://github.com/geatpy-dev/geatpy/tree/master/geatpy/doc，
geatpy官網：http://www.geatpy.com
若有錯誤之處，懇請同志們指正和諒解，謝謝！
因為是基於geatpy遺傳和進化演算法工具箱，所以下文的代碼部分在執行前，需要安裝geatpy：
pip install geatpy
下面切入主題：
?
本文只研究經典的遺傳演算法，力求最後能夠解決各種帶約束單目標優化問題，並能夠很輕鬆地進行擴展，讓大家不僅學到演算法理論，還能輕鬆地通過「複製粘貼」就能夠將相關遺傳演算法代碼結合到各類不同的現實問題的求解當中。
?
下面就來詳細講一下相關的理論和代碼實現：
正文
一. 基礎術語：
先介紹一下遺傳演算法的幾個基礎的術語，分別是：」個體「、」種群「、」編碼與解碼「、」目標函數值「、」適應度值「。
1.個體：「個體」其實是一個抽象的概念，與之有關的術語有：
（1）個體染色體：即對決策變數編碼後得到的行向量。
比如：有兩個決策變數x1=1,x2=2，各自用3位的二進位串去表示的話，寫成染色體就是：
（2）個體表現型：即對個體染色體進行解碼後，得到的直接指代各個控制變數的值的行向量。
比如對上面的染色體「0 0 1 0 1 0」進行解碼，得到 「1 2」，它就是個體的表現型，可看出該個體存儲著兩個變數，值分別是1和2。
注意：遺傳演算法中可以進行「實值編碼」，即可以不用二進位編碼，直接用變數的實際值來作為染色體。這個時候，個體的染色體數值上是等於個體的表現型的。
（3）個體可行性：即標識該個體是否是可行解。在遺傳演算法中，搜索空間內是允許出現非可行解的，此時為了標識哪些是可行解，哪些是非可行解，通常給可行解的個體標記1，非可行解的個體標記0。
2. 種群：「種群」也是一個抽象的概念，與之有關的術語有：
（1）種群染色體矩陣(Chrom)：它每一行對應一個個體的染色體。此時會發出疑問：一個個體可以有多條染色體嗎？答：可以有，但一般沒有必要，一條染色體就可以存儲很多決策變數的信息了，如果要用到多條染色體，可以用兩個種群來表示。
例如：（這裡由於矩陣的括弧比較難輸入，故省略）
?
它表示有3個個體（因為有3行），因此有3條染色體。此時需要注意：這些染色體代表決策變數的什麼值，我們是不知道的，我們也不知道該種群的染色體採用的是什麼編碼。染色體具體代表了什麼，取決於我們採用什麼方式去解碼。假如我們採用的是二進位的解碼方式，並約定上述的種群染色體矩陣中前3列代表第一個決策變數，後3列代表第二個決策變數，那麼，該種群就可以解碼成：（這裡由於矩陣的括弧比較難輸入，故省略）
（2）種群表現型矩陣(Phen)：它每一行對應一個個體的表現型。比如上圖就是根據Chrom種群染色體矩陣解碼得到的種群表現型矩陣。同樣地，當種群染色體採用的是「實值編碼」時，種群染色體矩陣與表現型矩陣實際上是一樣的。
（3）種群可行性列向量(LegV)：它每一行對應一個個體的可行性。比如：（這裡由於矩陣的括弧比較難輸入，故省略）
它表示上面種群中，第一個個體是非可行解（即個體表現型不符合約束條件），而第2、第3個個體都是可行解（標記為1）。
下面是代碼實現：
代碼1. 實數值種群的創建：
import numpy as np
from geatpy import crtrp
help(crtrp) # 查看幫助
# 定義種群規模（個體數目）
Nind = 4
# 創建「區域描述器」，表明有4個決策變數，範圍分別是[-3.1, 4.2], [-2, 2],[0, 1],[3, 3]
FieldDR=np.array([[-3.1, -2, 0, 3],
[4.2, 2, 1, 3]])
# 調用crtrp函數創建實數值種群
Chrom=crtrp(Nind, FieldDR)
print(Chrom)
代碼1的運行結果：
代碼2. 二進位種群的創建：
上面說過，二進位種群的染色體具體代表決策變數的什麼含義是不由染色體本身決定的，而是由解碼方式決定的。因此，創建二進位種群，實際上只需要創建只有0和1元素組成的整數值種群即可。因此crtip和crtbp兩個函數都可以完成該創建過程。
import numpy as np
from geatpy import crtip
help(crtip) # 查看幫助
# 定義種群規模（個體數目）
Nind = 4
# 創建「區域描述器」，表明有4個決策變數，範圍分別是[-3.1, 4.2], [-2, 2],[0, 1],[3, 3]
FieldDR=np.array([[0, 0, 0, 0, 0],
[1, 1, 1, 1, 1]])
# 調用crtip函數創建整數值種群
Chrom=crtip(Nind, FieldDR)
print(Chrom)
# 方法二：
#from geatpy import crtbp
#help(crtbp) # 查看幫助
## 定義種群規模（個體數目）
#Nind = 4
## 定義染色體長度
#Lind = 5
## 調用crtbp函數創建二進位值種群
#Chrom=crtbp(Nind, Lind)
#print(Chrom)
代碼2運行結果：
?
3. 編碼與解碼上面講過，實值編碼後染色體不需要解碼，但對於二進位編碼，則需要解碼才能得到染色體的表現型。有了表現型，才可以代入目標函數計算出種群各個個體的目標函數值。
上面也講到，二進位種群的染色體具體是代表多少個決策變數，各個決策變數具體是什麼值，是不確定的，取決於解碼的方式。因此解碼前我們要對解碼的方式作出聲明。如何做呢？我們可以採取一個「解碼矩陣」（或叫「區域描述器」）來指明具體如何解碼：
我們把這個解碼矩陣命名為FieldD，它是具有以下結構的列向量：
?其中，lens 包含染色體的每個子染色體的長度。sum(lens) 等於染色體長度。
lb 和ub 分別代表每個變數的上界和下界。
codes 指明染色體子串用的是標準二進位編碼還是格雷編碼。codes[i] = 0 表示第i個變數使用的是標準二進位編碼；codes[i] = 1 表示使用格雷編碼。
scales 指明每個子串用的是算術刻度還是對數刻度。scales[i] = 0 為算術刻度，scales[i] = 1 為對數刻度。對數刻度可以用於變數的範圍較大而且不確定的情況，對於大範圍的參數邊界，對數刻度讓搜索可用較少的位數，從而減少了遺傳演算法的計算量。
lbin 和ubin 指明了變數是否包含其範圍的邊界。0 表示不包含邊界；1 表示包含邊界。如對於範圍為(0, 1]的某個決策變數，其lbin值就為0，ubin值就為1。
下面是解碼的代碼實現，本例我們想讓二進位染色體解碼後表示整數型的決策變數，即決策變數是整數：
代碼3：
import numpy as np
from geatpy import crtbp
from geatpy import bs2int
help(crtbp) # 查看幫助
help(bs2int)
# 定義種群規模（個體數目）
Nind = 4
# 定義染色體長度
Lind = 5
# 調用crtbp函數創建二進位值種群
Chrom = crtbp(Nind, Lind)
print(染色體矩陣 =
, Chrom)
# 創建區域描述器
FieldD = np.array([[3, 2], # 各決策變數編碼後所佔位數
[0, 0], # 各決策變數的範圍下界
[7, 3], # 各決策變數的範圍上界
[0, 0], # 各決策變數採用什麼編碼方式(0為二進位編碼)
[0, 0], # 各決策變數是否採用對數刻度(0為採用算術刻度)
[1, 1], # 各決策變數的範圍是否包含下界(對bs2int實際無效，詳見help(bs2int))
[1, 1]])# 各決策變數的範圍是否包含上界(對bs2int實際無效)
Phen = bs2int(Chrom, FieldD)
print(表現型矩陣 =
, Phen)
運行效果如下：
?
4.目標函數值：種群的目標函數值存在一個矩陣裡面(一般命名為ObjV)，它每一行對應一個個體的目標函數值。對於單目標而言，這個目標函數值矩陣只有1列，而對於多目標而言，就有多列了，比如下面就是一個含兩個目標的種群目標函數值矩陣：
?
（這裡Nind表示種群的規模，即種群含多少個個體；Nvar表示決策變數的個數）下面來跑一下代碼，接著代碼3，在對二進位染色體解碼成整數值種群後，我們希望計算出f(x,y)=x+y這個目標函數值。此處不得不說一下，上文也說過，我們可以通過一個「種群可行性列向量」來標記哪些個體不符合約束條件。比如現在我想讓x為7的個體標記為不合法，完整代碼如下：
代碼4：
import numpy as np
from geatpy import crtbp
from geatpy import bs2int
def aimfuc(Phen, LegV):
x = Phen[:, [0]] # 取出種群表現型矩陣的第一列，得到所有個體的決策變數x
y = Phen[:, [1]] # 取出種群表現型矩陣的第二列，得到所有個體的決策變數y
f = x + y # 計算目標函數值
exIdx = np.where(x == 7)[0] # 找到x為7的個體在種群中的下標，
LegV[exIdx] = 0 # 標記其不合法
return f, LegV # 返回目標函數值和更新後的可行性列向量
help(crtbp) # 查看幫助
help(bs2int)
# 定義種群規模（個體數目）
Nind = 4
# 定義染色體長度
Lind = 5
# 調用crtbp函數創建二進位值種群
Chrom = crtbp(Nind, Lind)
print(染色體矩陣 =
, Chrom)
# 創建區域描述器
FieldD = np.array([[3, 2], # 各決策變數編碼後所佔位數
[0, 0], # 各決策變數的範圍下界
[7, 3], # 各決策變數的範圍上界
[0, 0], # 各決策變數採用什麼編碼方式(0為二進位編碼)
[0, 0], # 各決策變數是否採用對數刻度(0為採用算術刻度)
[1, 1], # 各決策變數的範圍是否包含下界(對bs2int實際無效，詳見help(bs2int))
[1, 1]])# 各決策變數的範圍是否包含上界(對bs2int實際無效)
Phen = bs2int(Chrom, FieldD)
print(表現型矩陣 =
, Phen)
LegV = np.ones((Chrom.shape[0], 1)) # 初始化可行性列向量
[ObjV, LegV] = aimfuc(Phen, LegV) # 計算目標函數值，同時更新可行性列向量
print(目標函數值矩陣 =
, ObjV)
print(種群可行性列向量 =
, LegV)
運行結果如下：
觀察該結果，我們發現種群第2個個體的x為7，因此種群可行性列向量的第2行為0。
疑問：對非可行解進行標記有什麼用呢？
答：標記非可行解，在含約束條件的優化問題中有用。對於含約束條件的優化問題，我們可以採用罰函數的方式來」懲罰「不滿足約束條件的非可行解。懲罰的方式有2種，一種是直接修改非可行解對應個體的目標函數值，使得它的目標值比其他可行解要」差「。用這種懲罰方法，標不標記非可行解已經不重要了，但是這種方法會有潛在的風險：一旦懲罰力度不夠，在進化過程中出現某一代裡面種群的所有個體都不滿足約束條件，此時就會對尋優的結果造成」干擾「：當需要記錄每一代的最優個體，進化完成後統計哪一代是最優時，上面的情況會有可能導致把最優結果判斷成正好是那個所有個體都不滿足約束條件的那一代，從而導致遺傳演算法尋優結果變成是一個非可行解，這就不符合我們尋優的初衷了。
因此，這種情況下，我們可以採用另一種懲罰方法：把種群中的非可行解個體的可行性標記為0，從而在尋找種群中的最優個體時，排除這些標記了0的」非法個體「。
此時要注意：給非可行解的可行性標記為0，並不意味著需要在進化過程中不給這些個體保留到下一代。因為有些情況下，最優解會出現在非可行域的附近！因此，我們一般只需降低這些非可行解個體的適應度，使之有更小的概率保留到下一代即可，而不是嚴格排除它們。上一段所說的」排除「，僅僅是指對種群的最優個體做記錄時不考慮非可行解。
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, LegV)
print(種群適應度矩陣 =
, FitnV)
運行結果：
?分析這個結果我們發現，由於第2個個體的可行性標記了0，因此得到適應度為0。而目標函數值最小的是第1和第3個個體（目標函數值為4），對應的適應度結果為1.66666...，可見遵循「最小化目標」的約定，即目標函數值越小，適應度越大。
好了，基本的術語和用法講完後，下面講一下遺傳演算法的基本運算元：
二. 遺傳演算法基本運算元：
我們不用破費時間去寫底層的遺傳運算元，因為geatpy工具箱提供豐富的遺傳運算元，例如：
可以用help(運算元名)來查看相關的用法，也可以參考官方API文檔，查看更詳細的用法和例子：
https://github.com/geatpy-dev/geatpy/tree/master/geatpy/doc/Geatpy-API
下面講一下理論：
1.選擇：
選擇是指根據個體適應度從種群中挑選個體保留到下一代的過程。其對應的是生物學中的」自然選擇」。
選擇運算元有：「輪盤賭選擇」、「隨機抽樣選擇」、「錦標賽選擇」、「本地選擇」、「截斷選擇」等等，這裡不展開敘述了，可以參考：
https://github.com/geatpy-dev/geatpy/blob/master/geatpy/doc/EA-introduction/%E7%AC%AC%E5%9B%9B%E7%AB%A0%EF%BC%9A%E9%80%89%E6%8B%A9.pdf
這裡要注意：遺傳演算法選擇出的後代是可以有重複的。
下面以低級選擇函數：錦標賽選擇運算元(tour)為例，使用help(tour)查看其API，或參見：
https://github.com/geatpy-dev/geatpy/blob/master/geatpy/doc/Geatpy-API/tour/tour.pdf
代碼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(個體的適應度為：
, FitnV)
print(選擇出的個體的下標為：
, chooseIdx)
運行結果：
光這樣還不夠，這裡只是得出了選擇個體的下標，我們還需要得到選擇後的種群染色體矩陣，這樣才能進行其他遺傳操作。於是，把代碼6修改一下，使用「高級選擇函數」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 = selecting(tour,Chrom, FitnV) # 使用tour錦標賽選擇運算元
print(個體的適應度為：
, FitnV)
print(選擇後得到的種群染色矩陣為：
, SelCh)
代碼7運行結果如下：
?將代碼7中的tour換成工具箱中的其他選擇運算元的名稱（如etour, rws, sus），就可以使用相應的選擇運算元進行選擇。
2.重組：
在很多的國內教材、博客文章、論文中，只提到「交叉」的概念。事實上，遺傳演算法的「交叉」是屬於「重組」運算元裡面的。因為交叉運算元經常使用，而對於「離散重組」、「中間重組」、「線性重組」等等，因為用的很少，所以我們常常只談「交叉」運算元了。交叉運算元實際上是「值互換重組」(Values exchanged recombination)。這裡就不展開敘述了，可以參考：
https://github.com/geatpy-dev/geatpy/blob/master/geatpy/doc/EA-introduction/%E7%AC%AC%E4%BA%94%E7%AB%A0%EF%BC%9A%E9%87%8D%E7%BB%84.pdf
與重組有關的遺傳演算法參數是「重組概率」（對於交叉而言就是「交叉概率」）(Pc)，它是指兩個個體的染色體之間發生重組的概率。
下面以兩點交叉(xovdp)為例，使用help(xovdp)查看其API，或參見：
https://github.com/geatpy-dev/geatpy/blob/master/geatpy/doc/Geatpy-API/xovdp/xovdp.pdf
代碼8：
from geatpy import xovdp
from geatpy import crtbp
help(xovdp)
OldChrom=crtbp(5,6) #調用crtbp創建一個5行6列的二進位種群矩陣
print(交叉前種群染色矩陣為：
, OldChrom)
NewChrom = xovdp(OldChrom, 1) # 設交叉概率為1
print(交叉後種群染色矩陣為：
, NewChrom)
代碼8運行結果如下：
?由此可見，xovdp是將相鄰的兩個個體進行交叉的，有人認為應該隨機選擇交叉個體。事實上，在遺傳演算法進化過程中，各個位置的個體是什麼，本身是隨機的，不必要在交叉這裡再增添一個隨機（當然，可以在執行xovdp兩點交叉之前，將種群染色體矩陣按行打亂順序然後再交叉，從而達到隨機選擇交叉個體的目的）。
3.變異：
變異是指通過改變父代染色體中的一部分基因來形成新的子代染色體的過程。它能夠保持種群的多樣性，降低遺傳演算法陷入局部最優解風險。最經典的兩種突變方法是實數值突變和離散突變。另外還有「互換突變」和其他有「特殊功能」的突變演算法。這裡就不展開敘述了，可以參考：
https://github.com/geatpy-dev/geatpy/blob/master/geatpy/doc/EA-introduction/%E7%AC%AC%E5%85%AD%E7%AB%A0%EF%BC%9A%E5%8F%98%E5%BC%82.pdf
與變異有關的遺傳演算法參數是「變異概率」(Pm)，它是指種群個體發生變異的概率。
下面以實值突變(mutbga)為例，編寫代碼實現實數值編碼的種群染色體的實值突變，使用help(mutbga)查看API，或參見：
https://github.com/geatpy-dev/geatpy/blob/master/geatpy/doc/Geatpy-API/mutbga/mutbga.pdf
代碼9：
import numpy as np
from geatpy import crtrp
from geatpy import mutbga
help(mutbga)
# 創建區域描述器FieldDR
FieldDR = np.array([
[8,7], # 變數下界
[10,10]]) # 變數上界
OldChrom = crtrp(3, FieldDR) # 調用crtrp函數創建實數值種群
print(變異前種群染色矩陣為：
, OldChrom)
NewChrom = mutbga(OldChrom, FieldDR, 0.5, 1) # 設置變異概率為0.5，
print(變異後種群染色矩陣為：
, NewChrom)
碼9的運行結果如下：
?4.重插入：
經典遺傳演算法通過選擇、重組和變異後，我們得到的是育種後代，此時育種後代的個體數有可能會跟父代種群的個體數不相同。這時，為了保持種群的規模，這些育種後代可以重新插入到父代中，替換父代種群的一部分個體，或者丟棄一部分育種個體，最終形成子代種群。
在geatpy工具箱中，採用"reins"函數來進行重插入，相關詳細用法可以用"help"命令查看，也可以參見：
https://github.com/geatpy-dev/geatpy/tree/master/geatpy/doc/Geatpy-API/reins
例如：現有四個決策變數，範圍分別是[-10,10]、[-5,5]、[-3,3]、[-1,1]。創建一個含有這4 個變數的6 個個體的實數值種群Chrom，同時再創建一個含有2 個個體的實數值種群SelCh（假定它就是「交叉變異後的育種種群」）來重插入到Chrom 中。
代碼10：
import numpy as np
from geatpy import crtrp
from geatpy import reins
help(reins)
FieldDR = np.array([[-10,-5,-3,-1],[10, 5, 3, 1]]) # 創建區域描述器
Chrom = crtrp(6, FieldDR) # 創建含有6個個體的種群，把它看作父代種群
# 創建列向量來存儲父代種群個體的目標函數值
FitnVCh = np.array([[21,22,23,16,15,24]]).T
SelCh=crtrp(2, FieldDR) # 創建含有2個個體的種群，看成是待重插入的育種種群
# 把育種個體重插入到父代種群中
print(重插入前種群染色矩陣為：
, Chrom)
NewChrom = reins(Chrom, SelCh, 1, 1, 1, FitnVCh) # 注意，傳入的Chrom是傳引用，重插入后里外都會變
print(育種種群染色矩陣為：
, SelCh)
print(重插入後種群染色矩陣為：
, NewChrom)
代碼10的運行結果如下：
?由此可見，原種群中適應度較低的2個個體被育種個體給替換掉了。
在一些改進的遺傳演算法中，可以不使用重插入。比如多目標優化NSGA2演算法中，通過對父代交叉、變異，將得到的子代與父代進行合併，再從合併的種群中選擇保留到下一代的若干個體。此時就不需要進行重插入了。用這種父子合併選擇的方法，是一種很好的精英保留策略。
講完上面的基本術語以及遺傳演算法基本運算元後，我們就可以來利用遺傳演算法的「套路」編寫一個「遺傳演算法模板」了：
三.完整實現遺傳演算法：
上文提到遺傳演算法程序可以寫成這個樣子：
a
a
其核心是演算法模板。在遺傳演算法模板里，我們根據遺傳演算法的「套路」，進行：初始化種群、目標函數值計算、適應度評價、選擇、重組、變異、記錄各代最優個體等操作。geatpy工具箱內置有開源的「套路模板」，源代碼參見：
https://github.com/geatpy-dev/geatpy/tree/master/geatpy/source-code/templets
這些內置演算法模板有詳細的輸入輸出參數說明，以及遺傳演算法整個流程的完整注釋，它們可以應對簡單或複雜的、單目標優化的、多目標優化的、約束優化的、組合優化的等等的問題。
但為了學習，我這裡自定義一個「演算法模板」，以解決一個帶約束的單目標優化問題：
編寫代碼 11、12、13，分別放在同一個文件夾下：
代碼11（目標函數aimfuc.py）（這裡要回顧一下前面，Phen是種群表現型矩陣，存儲的是種群所有個體的表現型，而不是單個個體。因而計算得到的目標函數值矩陣也是包含所有個體的目標函數值）：
# -*- coding: utf-8 -*-
"""
aimfc.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 aimfuc(Phen, LegV):
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]
LegV[exIdx1] = 0
LegV[exIdx2] = 0
return [f, LegV]
代碼12(自定義演算法模板templet.py)（由於偷懶，直接改geatpy工具箱的內置演算法模板實現）：
# -*- coding: utf-8 -*-
import numpy as np
import geatpy as ga
import time
import sys
def templet(AIM_M, AIM_F, PUN_M, PUN_F, FieldD, problem, maxormin, MAXGEN, NIND, SUBPOP, GGAP, selectStyle, recombinStyle, recopt, pm, distribute, drawing = 1):
"""
templet.py - 自定義單目標編程模板(二進位/格雷編碼種群)
語法：
該函數除了參數drawing外，不設置可預設參數。當某個參數需要預設時，在調用函數時傳入None即可。
比如當沒有罰函數時，則在調用編程模板時將第3、4個參數設置為None即可，如：
templet(AIM_M, aimfuc, None, None, ..., maxormin)
輸入參數：
AIM_M - 目標函數的地址，由AIM_M = __import__(目標函數所在文件名)語句得到
目標函數規範定義：[f,LegV] = aimfuc(Phen,LegV)
其中Phen是種群的表現型矩陣, LegV為種群的可行性列向量,f為種群的目標函數值矩陣
AIM_F : str - 目標函數名
PUN_M - 罰函數的地址，由PUN_M = __import__(罰函數所在文件名)語句得到
罰函數規範定義： newFitnV = punishing(LegV, FitnV)
其中LegV為種群的可行性列向量, FitnV為種群個體適應度列向量
一般在罰函數中對LegV為0的個體進行適應度懲罰，返回修改後的適應度列向量newFitnV
PUN_F : str - 罰函數名
FieldD : array - 二進位/格雷碼種群區域描述器，
描述種群每個個體的染色體長度和如何解碼的矩陣，它有以下結構:
[lens; (int) 每個控制變數編碼後在染色體中所佔的長度
lb; (float) 指明每個變數使用的下界
ub; (float) 指明每個變數使用的上界
codes; (0:binary | 1:gray) 指明子串是怎麼編碼的,
0為標準二進位編碼,1為各類編碼
scales; (0: rithmetic | 1:logarithmic) 指明每個子串是否使用對數或算術刻度,
1為使用對數刻度，2為使用算術刻度
lbin; (0:excluded | 1:included)
ubin] (0:excluded | 1:included)
lbin和ubin指明範圍中是否包含每個邊界。
選擇lbin=0或ubin=0,表示範圍中不包含相應邊界。
選擇lbin=1或ubin=1,表示範圍中包含相應邊界。
problem : str - 表明是整數問題還是實數問題，I表示是整數問題，R表示是實數問題
maxormin int - 最小最大化標記，1表示目標函數最小化；-1表示目標函數最大化
MAXGEN : int - 最大遺傳代數
NIND : int - 種群規模，即種群中包含多少個個體
SUBPOP : int - 子種群數量，即對一個種群劃分多少個子種群
GGAP : float - 代溝，表示子代與父代染色體及性狀不相同的概率
selectStyle : str - 指代所採用的低級選擇運算元的名稱，如rws(輪盤賭選擇運算元)
recombinStyle: str - 指代所採用的低級重組運算元的名稱，如xovsp(單點交叉)
recopt : float - 交叉概率
pm : float - 重組概率
distribute : bool - 是否增強種群的分布性（可能會造成收斂慢）
drawing : int - (可選參數)，0表示不繪圖，1表示繪製最終結果圖。默認drawing為1
輸出參數：
pop_trace : array - 種群進化記錄器(進化追蹤器),
第0列記錄著各代種群最優個體的目標函數值
第1列記錄著各代種群的適應度均值
第2列記錄著各代種群最優個體的適應度值
var_trace : array - 變數記錄器，記錄著各代種群最優個體的變數值，每一列對應一個控制變數
times : float - 進化所用時間
模板使用注意：
1.本模板調用的目標函數形如：[ObjV,LegV] = aimfuc(Phen,LegV),
其中Phen表示種群的表現型矩陣, LegV為種群的可行性列向量(詳見Geatpy數據結構)
2.本模板調用的罰函數形如: newFitnV = punishing(LegV, FitnV),
其中FitnV為用其他演算法求得的適應度
若不符合上述規範，則請修改演算法模板或自定義新演算法模板
3.關於maxormin: geatpy的內核函數全是遵循「最小化目標」的約定的，即目標函數值越小越好。
當需要優化最大化的目標時，需要設置maxormin為-1。
本演算法模板是正確使用maxormin的典型範例，其具體用法如下：
當調用的函數傳入參數包含與「目標函數值矩陣」有關的參數(如ObjV,ObjVSel,NDSetObjV等)時，
查看該函數的參考資料(可用help命令查看，也可到官網上查看相應的教程)，
裡面若要求傳入前對參數乘上maxormin,則需要乘上。
裡面若要求對返回參數乘上maxormin進行還原，
則調用函數返回得到的相應參數需要乘上maxormin進行還原，否則其正負號就會被改變。
"""
#==========================初始化配置==========================="""
# 獲取目標函數和罰函數
aimfuc = getattr(AIM_M, AIM_F) # 獲得目標函數
if PUN_F is not None:
punishing = getattr(PUN_M, PUN_F) # 獲得罰函數
NVAR = FieldD.shape[1] # 得到控制變數的個數
# 定義進化記錄器，初始值為nan
pop_trace = (np.zeros((MAXGEN ,2)) * np.nan)
# 定義變數記錄器，記錄控制變數值，初始值為nan
var_trace = (np.zeros((MAXGEN ,NVAR)) * np.nan)
ax = None # 存儲上一幀圖形
"""=========================開始遺傳演算法進化======================="""
Lind = np.sum(FieldD[0, :]) # 種群染色體長度
Chrom = ga.crtbp(NIND, Lind) # 生成初始種群
if problem == R:
variable = ga.bs2rv(Chrom, FieldD) # 解碼
elif problem == I:
if np.any(FieldD >= sys.maxsize):
variable = ga.bs2int(Chrom, FieldD).astype(object) # 解碼
else:
variable = ga.bs2int(Chrom, FieldD).astype(int64) # 解碼
LegV = np.ones((NIND, 1)) # 生成可行性列向量，元素為1表示對應個體是可行解，0表示非可行解
[ObjV, LegV] = aimfuc(variable, LegV) # 求種群的目標函數值
gen = 0
badCounter = 0 # 用於記錄在「遺忘策略下」被忽略的代數
# 開始進化！！
start_time = time.time() # 開始計時
while gen < MAXGEN:
if badCounter >= 10 * MAXGEN: # 若多花了10倍的迭代次數仍沒有可行解出現，則跳出
break
FitnV = ga.ranking(maxormin * ObjV, LegV, None, SUBPOP)
if PUN_F is not None:
FitnV = punishing(LegV, FitnV) # 調用罰函數
# 記錄進化過程
bestIdx = np.argmax(FitnV) # 獲取最優個體的下標
if LegV[bestIdx] != 0:
feasible = np.where(LegV != 0)[0] # 排除非可行解
pop_trace[gen,0] = np.sum(ObjV[feasible]) / ObjV[feasible].shape[0] # 記錄種群個體平均目標函數值
pop_trace[gen,1] = ObjV[bestIdx] # 記錄當代目標函數的最優值
var_trace[gen,:] = variable[bestIdx, :] # 記錄當代最優的控制變數值
# 繪製動態圖
if drawing == 2:
ax = ga.sgaplot(pop_trace[:,[1]],種群最優個體目標函數值, False, ax, gen)
badCounter = 0 # badCounter計數器清零
else:
gen -= 1 # 忽略這一代（遺忘策略）
badCounter += 1
if distribute == True: # 若要增強種群的分布性（可能會造成收斂慢）
idx = np.argsort(ObjV[:, 0], 0)
dis = np.diff(ObjV[idx,0]) / (np.max(ObjV[idx,0]) - np.min(ObjV[idx,0]) + 1)# 差分計算距離的修正偏移量
dis = np.hstack([dis, dis[-1]])
dis = dis + np.min(dis) # 修正偏移量+最小量=修正絕對量
FitnV[idx, 0] *= np.exp(dis) # 根據相鄰距離修改適應度，突出相鄰距離大的個體，以增加種群的多樣性
# 進行遺傳運算元
SelCh=ga.selecting(selectStyle, Chrom, FitnV, GGAP, SUBPOP) # 選擇
SelCh=ga.recombin(recombinStyle, SelCh, recopt, SUBPOP) # 對所選個體進行重組
SelCh=ga.mutbin(SelCh,pm) # 變異
# 計算種群適應度
if problem == R:
variable = ga.bs2rv(SelCh, FieldD) # 解碼
elif problem == I:
if np.any(FieldD >= sys.maxsize):
variable = ga.bs2int(SelCh, FieldD).astype(object) # 解碼
else:
variable = ga.bs2int(SelCh, FieldD).astype(int64)
LegVSel = np.ones((SelCh.shape[0], 1)) # 初始化育種種群的可行性列向量
[ObjVSel, LegVSel] = aimfuc(variable, LegVSel) # 求後代的目標函數值
FitnVSel = ga.ranking(maxormin * ObjVSel, LegVSel, None, SUBPOP) # 計算育種種群的適應度
if PUN_F is not None:
FitnVSel = punishing(LegVSel, FitnVSel) # 調用罰函數
# 重插入
[Chrom, ObjV, LegV]=ga.reins(Chrom, SelCh, SUBPOP, 1, 1, FitnV, FitnVSel, ObjV, ObjVSel, LegV, LegVSel)
# 計算新一代種群的控制變數解碼值
if problem == R:
variable = ga.bs2rv(Chrom, FieldD) # 解碼
elif problem == I:
if np.any(FieldD >= sys.maxsize):
variable = ga.bs2int(Chrom, FieldD).astype(object) # 解碼
else:
variable = ga.bs2int(SelCh, FieldD).astype(int64)
gen += 1
end_time = time.time() # 結束計時
times = end_time - start_time
# 後處理進化記錄器
delIdx = np.where(np.isnan(pop_trace))[0]
pop_trace = np.delete(pop_trace, delIdx, 0)
var_trace = np.delete(var_trace, delIdx, 0)
if pop_trace.shape[0] == 0:
raise RuntimeError(error: no feasible solution. (有效進化代數為0，沒找到可行解。))
# 繪圖
if drawing != 0:
ga.trcplot(pop_trace, [[種群個體平均目標函數值, 種群最優個體目標函數值]])
# 輸出結果
if maxormin == 1:
best_gen = np.argmin(pop_trace[:, 1]) # 記錄最優種群是在哪一代
best_ObjV = np.min(pop_trace[:, 1])
elif maxormin == -1:
best_gen = np.argmax(pop_trace[:, 1]) # 記錄最優種群是在哪一代
best_ObjV = np.max(pop_trace[:, 1])
print(最優的目標函數值為：%s%(best_ObjV))
print(最優的控制變數值為：)
for i in range(NVAR):
print(var_trace[best_gen, i])
print(有效進化代數：%s%(pop_trace.shape[0]))
print(最優的一代是第 %s 代%(best_gen + 1))
print(時間已過 %s 秒%(times))
# 返回進化記錄器、變數記錄器以及執行時間
return [pop_trace, var_trace, times]
代碼13(執行腳本main.py)：
運行結果如下：
?
四.後記：最後十分感謝SCUT、SCAU和Austin碩博聯合團隊提供的高性能實用型遺傳和進化演算法工具箱，它提供開源的進化演算法框架為遺傳演算法求解單目標/多目標優化、約束優化、組合優化等等給出了相當準確和快捷的解決方案。據網上其他相關博客的分析，在某些問題的求解上（如多目標優化），geatpy的運行效率相當的高，比matlab遺傳演算法工具箱、DEAP框架等快一個甚至數個數量級（視問題複雜情況而言）。下面放一張geatpy的層次結構和庫函數調用關係，以此記錄方便查看：
?
下面附一張一位在職的朋友使用犀牛軟體(Rhinoceros)結合geatpy工具箱進行產品優化設計的截圖：
?
很多工程軟體都提供Python介面，當需要用到進化優化時，就可以編寫Python代碼進行優化了。後續我將繼續學習和挖掘該工具箱的更多深入的用法。希望這篇文章在幫助自己記錄學習點滴之餘，也能幫助大家！
推薦閱讀：
