上一篇讲了Geatpy的库函数和数据结构https://blog.csdn.net/qq_33353186/article/details/82020507
熟悉了之后，我尝试编写一个进化算法模板解决一个约束优化问题：
在这个问题里，我们很容易知道目标函数的最小值为0，而这是个目标函数最大化问题，因此，我们可以对不符合约束条件的解的目标函数值设为0。而且不需要写罚函数了，直接把约束条件写在目标函数接口文件里即可。
编写目标函数接口文件：
""" aimfuc.py """
import numpy as np
# aim function
def aimfuc(variables, LegV): # LegV为种群的可行性列向量，详见"Geatpy数据结构"官方教程
x1 = variables[:, [0]] # get x1
x2 = variables[:, [1]] # get x2
f = 5 * x1 + 8 * x2
idx1 = np.where(x1 + x2 > 6)[0]
idx2 = np.where(5 * x1 + 9 * x2 > 45)[0]
# 采用惩罚方法1对非可行解进行惩罚，而不标记LegV为0（更多惩罚方法详见Geatpy官方教程的"函数接口及进化算法模板"章节
f[idx1] = 0
f[idx2] = 0
return [f, LegV]
下面一步步来编写进化算法模板，以解决这个约束优化问题。
编写算法模板时一个要注意的地方是：调用目标函数时，传入的参数必须和目标函数的参数列表一致。
我们试一下一种改进的遗传算法，即不用遗传算法的经典流程，而是在进化过程中，先对父代进行交叉和变异产生子代，然后父代和子代直接合并成一个大的种群，再从这个大的种群中挑选50%的优秀个体存活下来。在下一次进化迭代中，这50%的个体继续沿用上面的流程进行交叉、变异、选择。算法流程如下：
这种改进遗传算法的好处是具有很强的精英保留能力。可以使算法快速收敛，得到最优结果。
为了使模板函数能够用于更多的场景，能同时解决整数型和实数型控制变量的问题。于是编写进化算法模板如下：
# -*- coding: utf-8 -*-
import numpy as np
import geatpy as ga
import time
def templet(AIM_M, AIM_F, PUN_M, PUN_F, FieldDR, problem, maxormin, MAXGEN, NIND, SUBPOP, GGAP, selectStyle, recombinStyle, recopt, pm, distribute, drawing = 1):
"""
templet.py - 自定义的单目标编程模板(实值编码)
本模板实现改进单目标编程模板(实值编码)，将父子两代合并进行选择，增加了精英保留机制
语法：
该函数除了drawing外，不设置可缺省参数。当某个参数需要缺省时，在调用函数时传入None即可。
比如当没有罚函数时，则在调用编程模板时将第3、4个参数设置为None即可，如：
sga_new_real_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 - 罚函数名
FieldDR : array - 实际值种群区域描述器
[lb; (float) 指明每个变量使用的下界
ub] (float) 指明每个变量使用的上界
注：不需要考虑是否包含变量的边界值。在crtfld中已经将是否包含边界值进行了处理
本函数生成的矩阵的元素值在FieldDR的[下界, 上界)之间
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表示绘制最终结果图，2表示绘制动画。默认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'进行还原，否则其正负号就会被改变。
"""
"""==========================初始化配置==========================="""
GGAP = 0.5 # 因为父子合并后选择，因此要将代沟设为0.5以维持种群规模
# 获取目标函数和罚函数
aimfuc = getattr(AIM_M, AIM_F) # 获得目标函数
if PUN_F is not None:
punishing = getattr(PUN_M, PUN_F) # 获得罚函数
NVAR = FieldDR.shape[1] # 得到控制变量的个数
# 定义进化记录器，初始值为nan
pop_trace = (np.zeros((MAXGEN ,2)) * np.nan)
# 定义变量记录器，记录控制变量值，初始值为nan
var_trace = (np.zeros((MAXGEN ,NVAR)) * np.nan)
repnum = 0 # 初始化重复个体数为0
ax = None # 存储上一帧图形
"""=========================开始遗传算法进化======================="""
if problem == 'R':
Chrom = ga.crtrp(NIND, FieldDR) # 生成初始种群
elif problem == 'I':
Chrom = ga.crtip(NIND, FieldDR)
LegV = np.ones((NIND, 1)) # 生成可行性列向量，元素为1表示对应个体是可行解，0表示非可行解
[ObjV, LegV] = aimfuc(Chrom, LegV) # 求初代的目标函数值
gen = 0
badCounter = 0 # 用于记录在“遗忘策略下”被忽略的代数
# 开始进化！！
start_time = time.time() # 开始计时
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 distribute == True and 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) # 调用罚函数
# 记录进化过程
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,:] = Chrom[bestIdx, :] # 记录当代最优的控制变量值
repnum = len(np.where(ObjV[bestIdx] == ObjV)[0]) # 计算最优个体重复数
# 绘制动态图
if drawing == 2:
ax = ga.sgaplot(pop_trace[:,[1]],'种群最优个体目标函数值', False, ax, gen)
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
# 绘图
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])
if np.isnan(best_ObjV):
raise RuntimeError('error: no feasible solution. (没找到可行解。)')
print('最优的目标函数值为：', best_ObjV)
print('最优的控制变量值为：')
for i in range(NVAR):
print(var_trace[best_gen, i])
print('最优的一代是第', best_gen + 1, '代')
print('时间已过', times, '秒')
# 返回进化记录器、变量记录器以及执行时间
return [pop_trace, var_trace, times]
好了，现在该进化算法模板可以用于解决这类约束优化问题了。控制变量既可以是整数，也可以是实数。同时，通过设置maxormin的值来控制是最小化模板还是最大化目标。
下面编写执行脚本：
""" main.py """
import numpy as np
import geatpy as ga
from templet import templet # 导入自定义的进化算法模板templet
# 获取目标函数地址
AIM_M = __import__('aimfuc') # aimfuc.py要和该执行脚本放在同一目录下
# 变量设置
x1 = [0, 100] # 自变量1的范围
x2 = [0, 100] # 自变量2的范围
b1 = [1, 1] # 自变量1是否包含下界
b2 = [1, 1] # 自变量2是否包含上界
ranges = np.vstack([x1, x2]).T # 生成自变量的范围矩阵
borders = np.vstack([b1, b2]).T # 生成自变量的边界矩阵
FieldDR = ga.crtfld(ranges, borders) # 生成区域描述器
# 调用自定义的进化算法模板templet
[pop_trace, var_trace, times] = templet(AIM_M, 'aimfuc', None, None, FieldDR, problem = 'I', maxormin = -1, MAXGEN = 1000, NIND = 100, SUBPOP = 1, GGAP = 0.9, selectStyle = 'tour', recombinStyle = 'xovdp', recopt = 0.9, pm = 0.3, distribute = True, drawing = 1)因为这道题的控制变量x1和x2的区间都是单侧的，因此我们设个大一点的区间以便搜索最优解。我们采用锦标赛选择算子，并设置最大进化代数为1000，种群规模为100，使用两点交叉，交叉和变异概率分别是0.9和0.1。
运行结果如下：
最优的目标函数值为： 40.0
最优的控制变量值为：
0.0
5.0
最优的一代是第 85 代
时间已过 1.387918472290039 秒
至此，我们知道了如何自定义进化算法模板实现遗传算法了。在Geatpy中，内置了好几种进化算法模板。详细源码可以参见https://github.com/geatpy-dev/geatpy/tree/master/geatpy/source-code/templets
使用Geatpy内置模板会使上面的例子变得更加简单，只需编写目标函数接口和执行脚本就可以了。比如上面的例子，使用Geatpy框架，真正只需几行代码便可以解决：
""" main.py """
import numpy as np
import geatpy as ga
# define aim function
def aimfuc(variables, LegV):
x1 = variables[:, [0]] # get x1
x2 = variables[:, [1]] # get x2
f = 5 * x1 + 8 * x2
idx1 = np.where(x1 + x2 > 6)[0]
idx2 = np.where(5 * x1 + 9 * x2 > 45)[0]
f[idx1] = 0; f[idx2] = 0
return [f, LegV] # 对结果进行转置，使目标函数值矩阵符合Geatpy数据结构
# 获取目标函数地址
AIM_M = __import__('main')
# 变量设置
x1 = [0, 100] # 自变量1的范围
x2 = [0, 100] # 自变量2的范围
b1 = [1, 1] # 自变量1是否包含下界
b2 = [1, 1] # 自变量2是否包含上界
ranges = np.vstack([x1, x2]).T # 生成自变量的范围矩阵
borders = np.vstack([b1, b2]).T # 生成自变量的边界矩阵
FieldDR = ga.crtfld(ranges, borders) # 生成区域描述器
# 调用Geatpy内置进化算法模板
[pop_trace, var_trace, times] = ga.sga_new_real_templet(AIM_M, 'aimfuc', None, None, FieldDR, problem = 'I', maxormin = -1, MAXGEN = 1000, NIND = 100, SUBPOP = 1, GGAP = 0.9, selectStyle = 'tour', recombinStyle = 'xovdp', recopt = 0.9, pm = 0.3, distribute = True, drawing = 1)
关于罚函数的写法：“罚函数”原本是指数学上的对非可行解进行惩罚的“数学公式”，而在Geatpy编程中，罚函数一般指独立出来的对种群非可行解的适应度进行惩罚的函数”punishing”。
如何正确对非可行解进行惩罚？Geatpy教程的“函数接口及进化算法模板”章节有详细的叙述。下面概况如下：
惩罚方法1：在目标函数aimfuc 中，当找到非可行解个体时，当即对该个体的目标函数值进行惩罚。若为最小化目标，则惩罚时设置一个很大的目标函数值；反之设置一个很小的目标函数值。
惩罚方法2：在目标函数aimfuc 中，当找到非可行解个体时，并不当即对该个体的目标函数值进行惩罚，而是修改其在LegV 上对应位置的值为0，同时编写罚函数punishing，对LegV 为0 的个体加以惩罚——降低其适应度。事实上，在Geatpy 内置的算法模板中，已经对LegV 为0 的个体的适应度加以一定的惩罚，因此若是使用内置模板，则不需要编写罚函数punishing。此时若仍要编写罚函数punishing 的话，起到的是辅助性的适应度惩罚。
特别注意：如果采取上面的方法1 对非可行解进行惩罚，在修改非可行解的目标函数值时，必须设置一个“极大或极小”的值，即当为最小化目标时，要给非可行解设置一个绝对地比所有可行解还要大的值；反之要设置一个绝对比所有可行解小的值。否则会容易出现“被欺骗”的现象：即某一代的所有个体全是非可行解，而此时因为惩罚力度不足，在后续的进化中，再也没有可行解比该非可行解修改后的目标函数值要优秀(即更大或更小)。此时就会让进化算法“被欺骗”，得出一个非可行解的“最优”搜索结果。因此，在使用方法1 的同时，可以同时对LegV 加以标记为0，两种方法结合着对非可行解进行惩罚。
下面试一下采用惩罚方法2进行惩罚，加入一个罚函数punishing，代码修改如下：
""" main.py """
import numpy as np
import geatpy as ga
# 定义目标函数
def aimfuc(variables, LegV):
x1 = variables[:, [0]] # get x1
x2 = variables[:, [1]] # get x2
f = 5 * x1 + 8 * x2
# 约束条件
idx1 = np.where(x1 + x2 > 6)[0]
idx2 = np.where(5 * x1 + 9 * x2 > 45)[0]
exIdx = np.unique(np.hstack([idx1, idx2])) # 得到非可行解个体的下标
LegV[exIdx] = 0 # 标记非可行解在种群可行性列向量中对应的值为0(0表示非可行解，1表示可行解)
return [f, LegV]
def punishing(LegV, FitnV):
FitnV[np.where(LegV == 0)[0]] = 0 # 惩罚非可行解个体的适应度
return FitnV
if __name__ == "__main__":
AIM_M = __import__('main') # 获取目标函数地址
PUN_M = __import__('main') # 获取罚函数地址
# 变量设置
x1 = [0, 100] # 自变量1的范围
x2 = [0, 100] # 自变量2的范围
b1 = [1, 1] # 自变量1是否包含下界
b2 = [1, 1] # 自变量2是否包含上界
ranges = np.vstack([x1, x2]).T # 生成自变量的范围矩阵
borders = np.vstack([b1, b2]).T # 生成自变量的边界矩阵
FieldDR = ga.crtfld(ranges, borders) # 生成区域描述器
# 调用Geatpy内置进化算法模板
[pop_trace, var_trace, times] = ga.sga_new_real_templet(AIM_M, 'aimfuc', PUN_M, 'punishing', FieldDR, problem = 'I', maxormin = -1, MAXGEN = 1000, NIND = 100, SUBPOP = 1, GGAP = 0.9, selectStyle = 'tour', recombinStyle = 'xovdp', recopt = 0.9, pm = 0.3, distribute = True, drawing = 1)运行结果如下：
最优的目标函数值为： 40.0
最优的控制变量值为：
0.0
5.0
最优的一代是第 76 代
时间已过 2.6823768615722656 秒
在后面的版本中，Geatpy将不断丰富和改进内置的进化算法模板，以Geatpy的算法性能。比如添加MOEA/D多目标优化算法模板等。有相关想法的也可以email：jazzbin@geatpy.com。
更多教程可以参见官网http://www.geatpy.com ，以及参见Geatpy的源码和demo示例：https://github.com/geatpy-dev/geatpy/tree/master/geatpy/demo 你将有更多的收获。
在下一篇，我们将讲述如何使用Geatpy解决一些实际的有趣问题。
https://blog.csdn.net/qq_33353186/article/details/82047692 欢迎继续跟进，谢谢大家！
