遺傳演算法框架Geatpy學習之——基於網格化處理的多種群進化優化及其在含等式約束的優化問題中的應用
Geatpy是由華南理工大學、華南農業大學、德州奧斯汀公立大學學生聯合團隊開發的一款Python上的遺傳和進化演算法高效能權威框架。其效能遠高於matlab遺傳演算法工具箱及類似的諸如gatbx、GEATbx、gaot等第三方工具箱。尤其適合需要應用遺傳或其他進化演算法求解建模問題、與當前的其他專案結合、抑或是需要一個平臺來從事進化演算法研究與實驗的使用者。
在前面幾篇文章中，我們已經介紹了高效能Python遺傳和進化演算法框架——Geatpy的使用及一些案例。
（注意：Geatpy已於2018.10.03更新至1.0.7版本，這是一個重要更新，新版本完善了對約束優化問題的求解，重構了多目標優化非支配排序以及帕累託前沿搜尋演算法的核心，並提供網格化的單目標優化求解、多種群競爭的單目標優化演算法模板）。
詳細的Geatpy官方教程文件在github可以檢視：
================分割線==================
下面就等式約束問題展開論述：
含等式約束的目標優化問題是常見的比較難用遺傳演算法求解的一類問題。經典的處理方法有以下三類：
一. 直接採取罰函式來處理等式約束：
這種方法相對簡單，直接把遺傳演算法中不滿足等式約束條件的個體賦予最低的適應度，使其在種群進化中逐漸被淘汰。
然而，罰函式在處理不等式約束較為成功，但面對等式約束這類可行解空間遠比非可行解空間小的問題時，無論罰函式對處於非可行解的個體加以多麼嚴厲的懲罰，種群進化過程中會很難跳出非可行解區域，這就極有可能導致遺傳演算法出現連可行解都找不到
二. 將等式約束化成一個新的目標，用多目標優化演算法去處理等式約束問題：
例如一個等式約束為 x + y + z = 1，那麼，可以設計一個新的目標函式：
該目標表示 x + y + z - 1的絕對值越小越好。該目標越接近0，說明越接近滿足約束條件。當該目標的值為0時，對應的解就能滿足約束條件了。
這種方法可以極大地拓寬所研究問題的可行解的區域，從而使得遺傳演算法更方便地進行最優解的搜尋。
但是：這種方法也存在弊端。用這種方法對可行解的拓寬程度相對較大，並且，由於使用多目標優化演算法求解，它得到的結果是一個帕累托最優解集。而問題最容易出在：求解得到的帕累托最優解集裡面有時候會不包含
因此方法二也會導致演算法的不穩定。
三. 採用數學方法降維消除等式約束：
該方法是利用數學推導，利用等式約束，將其中一個或多個變數用其餘變數線性表出或非線性表出，從而使模型降維，同時也成功消除了等式約束，或是成功地將等式約束化為另一類不等式約束。
這種方法在效果上是最好的。但其侷限性是：當對於複雜的等式約束、無法將其中一個或多個變數用其餘變數線性表出或非線性表出時，該方法無法使用。
因此，應優先使用方法三，當不可行時再考慮方法一和二。
================分割線==================
採用Geatpy來求解這類含等式約束優化問題：
Geatpy解決這類約束優化問題極為簡單。首先建立目標函式，寫在“aimfuc.py”檔案中：
# -*- coding: utf-8 -*-
"""
aimfc.py - 目標函式demo
描述:
Geatpy的目標函式遵循本案例的定義方法，
若要改變目標函式的輸入引數、輸出引數的格式，則需要修改或自定義演算法模板
"""
def aimfuc(Phen):
x1 = Phen[:, [0]]
x2 = Phen[:, [1]]
x3 = 1 - x1 - x2 # 將x1 + x2 + x3 = 1的等式約束降維化處理
f = 4 * x1 + 2 * x2 + x3
return f
Geatpy v1.0.7版本之後，對罰函式作出了重要改進：罰函式不僅應該返回修改後的適應度矩陣，還應該返回不滿足約束條件的個體所在的下標。因此，編寫罰函式如下，寫在"punishing.py"檔案中：
# -*- coding: utf-8 -*-
"""
punishing.py - 罰函式demo
描述:
Geatpy的罰函式遵循本案例的定義方法，
若要改變罰函式的輸入引數、輸出引數的格式，則需要修改或自定義演算法模板
"""
import numpy as np
def punishing(Phen, FitnV):
x1 = Phen[:, [0]]
x2 = Phen[:, [1]]
x3 = 1 - x1 - x2
# 約束條件
idx1 = np.where(2 * x1 + x2 > 1)[0]
idx2 = np.where(x1 + 2 * x3 > 2)[0]
idx3 = np.where(x1 + x2 > 1)[0]
# 懲罰
FitnV[idx1] = 0
FitnV[idx2] = 0
FitnV[idx3] = 0
exIdx = np.unique(np.hstack([idx1, idx2, idx3])) # 得到非可行解在種群中的下標
return [FitnV, exIdx]
最後編寫執行指令碼"main.py"：
# -*- coding: utf-8 -*-
"""
執行指令碼main.py
描述：
該demo是展示如何計算帶等式約束的單目標優化問題:
max 4 * x1 + 2 * x2 + x3
s.t.
2 * x1 + x2 <= 1
x1 + 2 * x3 > 2
x1 + x2 + x3 == 1
0 <= x1 <= 1
0 <= x2 <= 1
0 <= x3 <= 2
其中目標函式寫在aimfuc.py檔案中，約束條件寫在罰函式檔案punishing.py中
本案例通過降維的方法，將等式約束化成了不等式約束，大大拓寬了可行解的空間，方便遺傳演算法求解
此外，本案例展示了利用多種群競爭的進化演算法模板sga_mpc_real_templet瞭解決該問題。
"""
import numpy as np
import geatpy as ga
# 獲取函式介面地址
AIM_M = __import__('aimfuc')
PUN_M = __import__('punishing')
# 變數設定
x1 = [0, 1] # 自變數1的範圍
x2 = [0, 1] # 自變數2的範圍
b1 = [1, 1] # 自變數1是否包含下界
b2 = [1, 1] # 自變數2是否包含上界
ranges=np.vstack([x1, x2]).T # 生成自變數的範圍矩陣
borders = np.vstack([b1, b2]).T # 生成自變數的邊界矩陣
precisions = [2] * 2 # 自變數的編碼精度，由於控制變數是連續型的，並且演算法模板中使用實值編碼，因此這裡設定成大於0的任意值即可
newRanges = ga.meshrng(ranges, gridnum = 2) # 對控制變數範圍進行網格化，網格邊長為2
# 生成網格化後的區域描述器集合
FieldDRs = []
for i in range(len(newRanges)):
FieldDRs.append(ga.crtfld(newRanges[i], borders, precisions))
# 呼叫程式設計模板
[pop_trace, var_trace, times] = ga.sga_mpc_real_templet(AIM_M, 'aimfuc', PUN_M,\
'punishing', FieldDRs, problem = 'R', maxormin = -1, MAXGEN = 50, NIND = 100,\
SUBPOP = 1, GGAP = 0.9, selectStyle = 'tour', recombinStyle = 'xovdprs',\
recopt = 0.9, pm = 0.3)
執行"main.py",得到結果如下：
最優的目標函式值為： 2.5 最優的控制變數值為： 0.5 0.0 最優的一代是第 5 代 時間已過 0.5979366302490234 秒
分析：sga_mpc_real_templet是Geatpy v1.0.7版本新增的一個內建演算法模板，其原始碼見：
該模板實現了基於多種群競爭進化的單目標優化演算法。由於多種群是在相互競爭中進化的，因此收斂速度極快。並且，每個種群的最大進化代數也不需要設定太大，本例中設定50代就夠了。
Geatpy v1.0.7另外還提供了多種群獨立進化的單目標優化演算法模板:sga_mps_real_templet。與上面的不同的是該演算法模板中種群是互相獨立進化，最後再挑選最優個體的，而不是在競爭中進化。
在呼叫多種群進化演算法模板sga_mpc_real_templet和sga_mps_real_templet前，常根據各控制變數的範圍先對其進行網格化，使其網格化成若干個小的範圍，然後不同範圍用不同的種群去進化搜尋。變數範圍網格化的庫函式為：meshrng，用法詳見：https://github.com/geatpy-dev/geatpy/blob/master/geatpy/doc/API/meshrng/meshrng.pdf
更多的demo詳見：
歡迎繼續跟進，感謝！