Geatpy遺傳演算法在曲線尋優上的初步探究

摘要： 園子裡關於遺傳演算法的教案不少，但基於geatpy框架的並未多見，故分享此文以作參考，還望廣大園友多多指教！ Geatpy出自三所名校聯合團隊之手，是遺傳演算法領域的權威框架（python），其效率之高、應用領域之廣遠勝諸多第三方工具，此處不作贅述，直接上鍊接： 官網：http:/...

園子裡關於遺傳演算法的教案不少，但基於geatpy框架的並未多見，故分享此文以作參考，還望廣大園友多多指教！

Geatpy出自三所名校聯合團隊之手，是遺傳演算法領域的權威框架（python），其效率之高、應用領域之廣遠勝諸多第三方工具，此處不作贅述，直接上鍊接：

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

原始碼：https://github.com/geatpy-dev/geatpy/tree/master/geatpy

使用 Geatpy 需要安裝geatpy模組（pip install geatpy），linux下如果裝完後import時出現報錯，可以下載我帖尾連結裡的wheel檔案進行安裝。

言歸正傳，根據經典的遺傳演算法流程，無外乎這幾個步驟： 種群初始化 ->（ 適應度評價 -> 遴選 -> 交叉 -> 變異 ）<- 迴圈進化直至終止條件達標。

當然，有關遺傳演算法的原理和過程不做深討，本文旨在剖析遺傳演算法在陣曲線尋優上的高效應用，我寫了一個簡單示例來幫助大家更好的理解，程式碼如下：

# -*- coding: utf-8 -*-
"""punishing.py - 罰函式demo"""

import numpy as np

def punishing(LegV, FitnV):
FitnV[np.where(LegV == 0)[0]] = 0
return FitnV

# -*- coding: utf-8 -*-
""" aimfc.py即目標函式，本例通過輸入每一代的染色體，由自定義評價函式計算與目標曲線的面積差，作為目標函式值ObjV來輸出 """

import numpy as np

def MakeObjCurve(width):
''' 建立目標曲線，此處定義為一組正弦波拼接序列 '''
n1 = width//3
n2 = width - 2*n1
x1=np.cos(np.arange(0,n1)) * 1
x2=np.cos(np.arange(0,n1)) * 4
x3=np.cos(np.arange(0,n2)) * 2
ObjCurve=np.hstack((x1,x2,x3))
return ObjCurve

def CalScore(chrom):
''' 返回染色體與目標曲線之間的面積的倒數作為評分值 '''
objCurve = MakeObjCurve(len(chrom))
area = chrom - objCurve
area *= 10**5#調整係數確保分值不受小數項干擾
score = 1 / np.dot(area, area)#計算差值的平方和以簡化求面積過程
return score

def myEvaFunc(chroms):
''' 自定義評價函式，以評分值作為目標函式值 '''
scores = []
for chrom in chroms:
score = CalScore(chrom)
scores.append(score)
scores = np.array([scores]).T
return scores

def aimfuc(Phen, LegV):

ObjV = myEvaFunc(Phen)
exIdx = np.argmin(ObjV[:, 0])

# 懲罰方法2： 標記非可行解在可行性列向量中對應的值為0，並編寫punishing罰函式來修改非可行解的適應度。
# 也可以不寫punishing，因為Geatpy內建的演算法模板及核心已經對LegV標記為0的個體的適應度作出了修改。
# 使用punishing罰函式實質上是對非可行解個體的適應度作進一步的修改
LegV[exIdx] = 0 # 對非可行解作出標記，使其在可行性列向量中對應的值為0，此處標記的是得分最小項

return [ObjV, LegV]

# -*- coding: utf-8 -*-
""" main.py即主函式，本例僅用於演示“已知曲線尋優”的過程 """

import numpy as np
import geatpy as ga
import time
import matplotlib.pyplot as plt

def search_objects(directory):
directory=os.path.normpath(directory)#規格化，防止分隔符造成的差異
if not os.path.isdir(directory):
raise IOError("The directory '"+"' doesn't exist!")
objects={}
for curdir,substrs,files in os.walk(directory):
for jpeg in (file for file in files if file.endswith('.csv')):
path=os.path.join(curdir,jpeg)
label=path.split(os.path.sep)[-2]
if label not in objects:
objects[label]=[]
objects[label].append(path)
return objects

def sga_mps_real_templet(AIM_M, AIM_F, PUN_M, PUN_F, FieldDRs, problem, maxormin, MAXGEN, NIND, SUBPOP, GGAP, selectStyle, recombinStyle, recopt, pm, distribute, drawing = 1):
    """ 基於多種群獨立進化單目標程式設計模板(實值編碼)，各種群獨立將父子兩代合併進行選擇，採取精英保留機制 """

#==========================初始化配置===========================
GGAP = 0.5 # 因為父子合併後選擇，因此要將代溝設為0.5以維持種群規模
# 獲取目標函式和罰函式
aimfuc = getattr(AIM_M, AIM_F) # 獲得目標函式
if PUN_F is not None:
punishing = getattr(PUN_M, PUN_F) # 獲得罰函式
NVAR = FieldDRs[0].shape[1] # 得到控制變數的個數
# 定義全域性進化記錄器，初始值為nan
pop_trace = (np.zeros((MAXGEN ,2)) * np.nan)
pop_trace[:, 0] = 0
# 定義變數記錄器，記錄控制變數值，初始值為nan
var_trace = (np.zeros((MAXGEN ,NVAR)) * np.nan)
"""=========================開始遺傳演算法進化======================="""
start_time = time.time() # 開始計時
# 對於各個網格分別進行進化，採用全域性進化記錄器記錄最優值
for index in range(len(FieldDRs)): # 遍歷各個子種群，各子種群獨立進化，互相不競爭
FieldDR = FieldDRs[index]
if problem == 'R':
Chrom = ga.crtrp(NIND, FieldDR) # 生成初始種群
elif problem == 'I':
Chrom = ga.crtip(NIND, FieldDR)
LegV = np.ones((NIND, 1)) # 初始化種群的可行性列向量
[ObjV, LegV] = aimfuc(Chrom, LegV) # 求初始種群的目標函式值
repnum = 0 # 初始化重複個體數為0
ax = None # 儲存上一幀圖形
gen = 0
badCounter = 0 # 用於記錄在“遺忘策略下”被忽略的代數
# 開始進化！！
while gen < MAXGEN:
if badCounter >= 10 * MAXGEN: # 若多花了10倍的迭代次數仍沒有可行解出現，則跳出
break
# 進行遺傳運算元，生成子代
SelCh = ga.recombin(recombinStyle, Chrom, recopt, SUBPOP) # 重組
if problem == 'R':
SelCh = ga.mutbga(SelCh,FieldDR, pm) # 變異
if repnum > Chrom.shape[0] * 0.01: # 當最優個體重複率高達1%時，進行一次高斯變異
SelCh = ga.mutgau(SelCh, FieldDR, pm) # 高斯變異
elif problem == 'I':
SelCh = ga.mutint(SelCh, FieldDR, pm)
LegVSel = np.ones((SelCh.shape[0], 1)) # 初始化育種種群的可行性列向量
[ObjVSel, LegVSel] = aimfuc(SelCh, LegVSel) # 求育種種群的目標函式值
# 父子合併
Chrom = np.vstack([Chrom, SelCh])
ObjV = np.vstack([ObjV, ObjVSel])
LegV = np.vstack([LegV, LegVSel])
FitnV = ga.ranking(maxormin * ObjV, LegV, None, SUBPOP) # 適應度評價
if PUN_F is not None:
FitnV = punishing(LegV, FitnV) # 呼叫懲罰函式
repnum = len(np.where(ObjV[np.argmax(FitnV)] == ObjV)[0]) # 計算最優個體重複數
# 記錄進化過程
bestIdx = np.argmax(FitnV)
if (LegV[bestIdx] != 0) and ((np.isnan(pop_trace[gen,1])) or ((maxormin == 1) & (pop_trace[gen,1] >= ObjV[bestIdx])) or ((maxormin == -1) & (pop_trace[gen,1] <= ObjV[bestIdx]))):
feasible = np.where(LegV != 0)[0] # 排除非可行解
pop_trace[gen,0] += np.sum(ObjV[feasible]) / ObjV[feasible].shape[0] / len(FieldDRs) # 記錄種群個體平均目標函式值
pop_trace[gen,1] = ObjV[bestIdx] # 記錄當代目標函式的最優值
var_trace[gen,:] = Chrom[bestIdx, :] # 記錄當代最優的控制變數值
# 繪製動態圖
if drawing == 2:
ax = ga.sgaplot(pop_trace[:,[1]],'子種群'+str(index+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) # 根據相鄰距離修改適應度，突出相鄰距離大的個體，以增加種群的多樣性
[Chrom, ObjV, LegV] = ga.selecting(selectStyle, Chrom, FitnV, GGAP, SUBPOP, ObjV, LegV) # 選擇
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 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))
# 繪圖
if drawing != 0:
ga.trcplot(pop_trace, [['種群個體平均目標函式值', '種群最優個體目標函式值']])
# 返回進化記錄器、變數記錄器以及執行時間
return [pop_trace, var_trace, times, best_gen]

# 獲取函式介面地址
AIM_M = __import__('aimfuc')
PUN_M = __import__('punishing')
POP_SIZE = 300# 種群高度
CHROM_LENGTH = 20# 染色體寬度
max_generation = 150# 進化代數
chrom_bottom = -4#染色體數值下限
chrom_top = 4#染色體數值上限

# 變數設定
x = []; b = []
for i in range(CHROM_LENGTH):
x.append([chrom_bottom, chrom_top]) # 自變數的範圍
b.append([0, 0]) # 自變數是否包含下界
ranges=np.vstack(x).T # 生成自變數的範圍矩陣
borders = np.vstack(b).T # 生成自變數的邊界矩陣
precisions = [1]*CHROM_LENGTH # 在二進位制/格雷碼編碼中代表自變數的編碼精度，當控制變數是連續型時，根據crtfld參考資料，該變數只表示邊界精度，故設定為一定的正數即可
# 生成網格化後的區域描述器集合
FieldDRs = []
for i in range(1):
FieldDRs.append(ga.crtfld(ranges, borders, precisions))

# 呼叫程式設計模板(設定problem = 'R'處理實數型變數問題，詳見該演算法模板的原始碼)
[pop_trace, var_trace, times, best_gen] = sga_mps_real_templet(AIM_M, 'aimfuc', PUN_M, 'punishing',
 FieldDRs, problem = 'R', maxormin = -1, MAXGEN = max_generation, NIND = POP_SIZE, SUBPOP = 1, GGAP = 0.9, \
 selectStyle = 'tour', recombinStyle = 'xovdprs', recopt = 0.9, pm = 0.3, distribute = True, drawing = 1)

bstChrom = var_trace[best_gen]
objCurve = AIM_M.MakeObjCurve(CHROM_LENGTH)

plt.ion()
fig = plt.figure('曲線尋優演示',facecolor='lightgray')
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
ax1.set_title("Evaluation Map")
ax1.grid(axis='y', linestyle=':')
for i in range(max_generation):
if i%5==0:
ax1.plot(var_trace[i], 'o-')
ax2.cla()
ax2.set_title("最優染色體[gen:%i]"%(i+1))
ax2.plot(var_trace[i], 'o-', c='dodgerblue')
plt.pause(0.001)
ax2.cla()
ax2.grid(axis='y', linestyle=':')
ax2.plot(objCurve, 'o-', c='orangered', label='目標曲線')
ax2.plot(bstChrom, 'o-', c='dodgerblue', label='最優染色體[gen:%i]'%(best_gen+1))
plt.legend()
plt.ioff()
plt.show()

請注意，此處我已 將模板函式單獨放到主函式中 以便大家更好的理解，返回值中增加了最優代數以便後續圖例的顯示。

本例採用的進化模板是 sga_mps_real_templet ，基於多種群進化單目標（實數值），用於實現 尋找目標曲線 。

目標曲線的定義函式在aimfc.py 的 MakeObjCurve函式中 ，本例為3段振幅不同的cos函式拼接而成的模擬曲線，寬度20。

種群初始值設定：種群高度300、染色體寬度20（與目標曲線寬度保持一致）、進化代數150、染色體數值上下限[-4,4]（與 目標 曲線的上下限保持一致）。

接著我們開始尋優， 通過 CalScore 函式計算每代種群的每條染色體與目標函式之間的差值，經過一定的係數轉換得到評分值（差值越大，評分越低），以單目標（1列）形式輸出，由geatpy的ranking函式來決定適應度評價，然後繼續遴選、交叉、變異，如此迴圈往復，直至達到近似目標值時終止（本例設定為150代時終止），尋優過程如下圖：

最優的目標函式值為：7.287567405118076e-09
有效進化代數：150
最優的一代是第 148 代
時間已過 1.3259999752044678 秒

可以看到從50代左右優化曲線開始顯著上揚，直到130代左右逐漸平緩，並且耗時非常少，來看尋優結果圖：

可以看到除了第6個點的數值有細微差異之外，其他點幾乎都是吻合的，基本實現了目標曲線的尋求。

之所以丟擲本例，最重要的一點在於geatpy遺傳演算法 不僅能尋求已知的目標函式 ， 還可以通過自定義的評分體系或第三方介面來參與實現尋優 過程，只需將 CalScore 函式稍作改動即可，以上。

【wheel檔案】: https://pan.baidu.com/s/1BwLq_m3Dd5RMqatvTYXrAw 提取碼: vgkz