函数接口及进化算法模板
之前的章节中，展示了利用Geatpy 简单进化算法模板解决多元单峰函数最值的搜索问题。在那时，我们把所有的代码都放在一个脚本文件里。对于解决简单的问题这种方式能够轻松胜任，但非常不利于重构和把Geatpy 与其他算法或项目进行融合。本章我们将介绍如何通过编写函数接口以及精简的进化算法模板来实现低耦合的编程。
1.函数接口
函数接口是目标函数和罚函数的统称，一般写在独立的Python 源文件里，你也可以直接把它们定义在脚本中。目标函数传入种群表现型矩阵Phen 以及种群可行性列向量LegV ，计算各个个体的目标函数值，若有约束条件，此时将对不满足约束条件的个体在LegV 上对应的值设置为0。最后返回种群的目标函数值列向量ObjV 以及修改后的种群可行性列向量LegV ；罚函数传入种群可行性列向量LegV 和适应度值列向量FitnV ，此时对于LegV 中值为0 的个体(即非可行解个体)“加以惩罚”，降低其适应度，最后返回新的适应度列向量NewFitnV 。
在Geatpy 有两种方式来对非可行解的个体进行惩罚：
1) 在目标函数aimfuc 中，当找到非可行解个体时，当即对该个体的目标函数值进行惩罚。若为最小化目标，则惩罚时设置一个很大的目标函数值；反之设置一个很小的目标函数值。
2) 在目标函数aimfuc 中，当找到非可行解个体时，并不当即对该个体的目标函数值进行惩罚，而是修改其在LegV 上对应位置的值为0，同时编写罚函数punishing，对LegV 为0 的个体加以惩罚——降低其适应度。事实上，在Geatpy 内置的算法模板中，已经对LegV 为0 的个体的适应度加以一定的惩罚，因此若是使用内置模板，则不需要 编写罚函数punishing。此时若仍要编写罚函数punishing 的话，起到的是辅助性的适应度惩罚。
(这里的“罚函数punishing”跟数学上的“罚函数”含义是不一样的。后者是纯粹的数学公式，而前者是值使用Geatpy 时自定义的一个名为’punishing’ 的根据可行性列向量LegV 来对非可行解的适应度FitnV 进行惩罚的一个函数。)
特别注意：如果采取上面的方法1 对非可行解进行惩罚，在修改非可行解的目标函数值时，必须设置一个“极大或极小”的值，即当为最小化目标时，要给非可行解设置一个绝对地比所有可行解还要大的值；反之要设置一个绝对比所有可行解小的值。否则会容易出现“被欺骗”的现象：即某一代的所有个体全是非可行解，而此时因为惩罚力 度不足，在后续的进化中，再也没有可行解比该非可行解修改后的目标函数值要优秀(即更大或更小)。此时就会让进化算法“被欺骗”，得出一个非可行解的“最优”搜索结果。因此，在使用方法1 的同时，可以同时对LegV 加以标记为0，两种方法结合着对非可行解进行惩罚。
例：设计目标函数和罚函数，分别命名为aimfuc 和punishing。其中目标函数实现的是2 个变量的平方和，罚函数实现的是惩罚值为0 的变量。
注意：若使用Geatpy 的内置算法模板，目标函数和罚函数必须按照下述的输入输出格式进行定义。
"""目标函数aimfuc.py"""
import numpy as np
def aimfuc(Phen, LegV): # 传入种群染色体矩阵解码后的基因表现型矩阵
x1 = Phen[:, [0]] # 从Phen中片取得到x1变量
x2 = Phen[:, [1]]
ObjV = x1*x1 + x2*x2
exIdx = np.where(Phen == 0)[0] # 得到非可行解在种群中的位置
# 采用方法2对非可行解进行标记，在punishing中对其进行惩罚
LegV[exIdx] = 0 # 标记非可行解为0
return [ObjV,LegV]
传入的Phen 代表种群染色体的表现型，它是一个矩阵，LegV 代表种群个体的可行性列向量(其中元素值为0 表示对应个体是非可行解，为1 表示对应的是可行解)。注意返回的是ObjV 和LegV 两个参数，前者是一个numpy 的array 类型列向量，每一列对应一个个体的目标函数值。
假设传入aimfuc 函数的Phen 的值如下： LegV的值如下： 那么调用aimfuc(Phen, LegV) 后，将会得到种群个体对应的目标函数值为： 新的LegV为： 进一步地理解numpy 的array 类型变量的维度，我们执行print(Phen.shape) 和print(ObjV.shape)输出它们的维度信息，可看到结果分别(3, 2) 和(3, 1)。说明Phen 是3 行2 列的矩阵(实际上是数组类型)；ObjV 是3 行1 列的列向量。
下面继续编写罚函数punishing，对LegV 标记了0 的变量进行惩罚，使其适应度比当前种群最小适应度值还要小至少50%：
"""罚函数punishing.py"""
import numpy as np
def punishing(LegV, FitnV):
FitnV[np.where(LegV == 0)[0]] = np.min(FitnV) // 2 # 取整除法
return FitnV
下面是完整执行上述例子的代码：
"""test.py"""
import numpy as np
from punishing import punishing
from aimfuc import aimfuc
import geatpy as ga
# 创建Phen代表种群的基因表现型矩阵，注意array的维度
Phen = np.array([
[1, 2],
[4, 0],
[3, 4]])
LegV = np.ones((3, 1)) # 初始化种群可行性列向量
[ObjV, LegV] = aimfuc(Phen, LegV)
# 调用ranking计算适应度(因为若给ranking传入LegV参数时，ranking函数会自动
# 对LegV标记为0的非可行解进行惩罚，这里为了展示punishing的作用，因此调用
# ranking时不传入LegV)
FitnV = ga.ranking(ObjV)
FitnV=punishing(LegV,FitnV) # 得到新的适应度以及非可行解的下标
print(FitnV) # 输出结果
输出结果为： 若不调用punishing，将会得到： 可见种群第二个个体成功地被“惩罚”，其适应度变成了0。 注意事项： 在采用方法1 来惩罚非可行解时，要注意一个很容易出错的地方，请看下面的例子： 满足如下约束 这是一个双目标优化问题，其中约束条件不再是简单的范围区间，而是多个变量的多个约束不等式。此时我们想让不符合约束条件的解对应的目标函数值f1变大，同时f2 变小，于是可以设计以下罚函数： 其中，若满足约束条件gi(x)(i = 1; 2; :::; 6)，取p1 = p2 = 1，反之，取p1 < 0 和p2 > 0 的随机数。
这里看上去并无问题，但实际上，我们不建议取p1 < 0 的随机数，因为它会改变数值的符号。因为目标函数f1 是恒不大于0 的，如果罚函数遍历所有约束条件来依次让xi 不满足约束条件的目标函数值乘上p1，那么此时稍有不慎就会极有可能出现目标函数值多次被修改的情况，因为不满足各个约束条件的非可行解可能会有交集，而p1 < 0，所以此时着意味着这些目标函数值会被反复变换正负号，这将导致罚函数反而将目标函数值变得更小的错误（“惩罚出错”)。对此，一定不能依次地对不满足各个约束条件的非可行解进行惩罚，而是要取不满足各个约束条件的非可行解的全集，对这个全集中的非可行解进行惩罚，这样能避免上述问题的“惩罚出错”的情况出现。另一种解决办法是像上文所说的，让非可行解设置一个绝对比所有可行解要大的值。这个需要数学上证明才可去像这样设值。这样也能避免“惩罚出错”的情况出现。
因此，更建议利用惩罚方法2 来对可行解进行惩罚。
2.进化算法模板
前面我们已经介绍过进化算法模板的概念和重要性。下面我们从头开始自定义一个遗传算法模板，命名为mintemp1，并编写测试脚本，调用该遗传算法模板来解决搜索
的最小值，其中x1 与x2 分别是[-5, 0)U(0, 5] 以及(2, 10] 的整数。
模板中调用的Geatpy 内置函数的相关用法可以在”Geatpy 函数” 章节中找到详细的讲解。
## tyest
# -*- coding: utf-8 -*-
"""自定义进化算法模板mintemp1.py"""
import numpy as np
import geatpy as ga # 导入geatpy库
import time
def mintemp1(AIM_M, AIM_F, PUN_M, PUN_F, ranges, borders, MAXGEN,
NIND, SUBPOP, GGAP, selectStyle, recombinStyle, recopt, pm,
maxormin):
"""==========================初始化配置==========================="""
# 获取目标函数和罚函数
aimfuc = getattr(AIM_M, AIM_F) # 获得目标函数
punishing = getattr(PUN_M, PUN_F) # 获得罚函数
FieldDR = ga.crtfld(ranges, borders) # 初始化区域描述器
NVAR = ranges.shape[1] # 得到控制变量的个数
# 定义进化记录器，初始值为nan
pop_trace = (np.zeros((MAXGEN ,3)) * np.nan).astype('int64')
var_trace = (np.zeros((MAXGEN ,NVAR)) * np.nan).astype('int64')
"""=========================开始遗传算法进化======================="""
Chrom = ga.crtip(NIND, FieldDR) #根据区域描述器FieldDR生成整数型初始种群
LegV = np.ones((NIND, 1)) #生成可行性列向量，元素为1表示对应个体是可行解，0表示非可行解
[ObjV, LegV] = aimfuc(Chrom, LegV) #计算种群目标函数值，同时更新LegV
start_time = time.time() # 开始计时
# 开始进化！！
for gen in range(MAXGEN):
FitnV = ga.ranking(maxormin * ObjV, LegV) # 计算种群适应度
FitnV = punishing(LegV, FitnV) # 调用罚函数
# 记录进化过程
bestIdx = np.argmax(FitnV)
if LegV[bestIdx] != 0:
feasible = np.where(LegV != 0)[0] # 排除非可行解
# 记录当代种群的适应度均值
pop_trace[gen, 1] = np.sum(FitnV[feasible]) / FitnV[feasible].shape[0]
# 记录当代种群最优个体的目标函数值
pop_trace[gen, 0] = ObjV[bestIdx]
# 记录当代种群的最优个体的适应度值
pop_trace[gen, 2] = FitnV[bestIdx]
# 记录当代种群最优个体的变量值
var_trace[gen, :] = Chrom[bestIdx, :]# 进行遗传操作！！
SelCh=ga.selecting(selectStyle, Chrom, FitnV, GGAP, SUBPOP) #选择
SelCh=ga.recombin(recombinStyle, SelCh, recopt, SUBPOP) #交叉
SelCh=ga.mutint(SelCh, FieldDR, pm) # 实值变异
LegVSel = np.ones((SelCh.shape[0], 1)) #创建育种个体的可行性列向量
[ObjVSel, LegVSel] = aimfuc(SelCh, LegVSel)#求育种个体的目标函数值
FitnVSel = punishing(LegVSel, FitnV) # 调用罚函数
[Chrom,ObjV,LegV] =ga.reins(Chrom,SelCh,SUBPOP,1,1,FitnV,FitnVSel,ObjV,ObjVSel,LegV,LegVSel) #重插入
end_time = time.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)
# 返回进化记录器、变量记录器以及执行时间
return [pop_trace, var_trace, end_time - start_time]
main py文件
# -*- coding: utf-8 -*-
"""
执行脚本main.py
描述：
该demo是展示如何计算带约束的单目标优化问题
本案例通过自定义算法模板"mintemp1.py"来解决该问题
其中目标函数写在aimfuc.py文件中，约束条件写在罚函数文件punishing.py中
"""
import numpy as np
import geatpy as gea # 导入geatpy库
from mintemp1 import mintemp1 # 导入自定义的编程模板
%matplotlib inline
# 获取函数接口地址
AIM_M = __import__('aimfuc') # 目标函数
PUN_M = __import__('punishing') # 罚函数
"""============================变量设置============================"""
x1 = [-5, 5]; x2 = [2, 10] # 自变量的范围
b1 = [1, 1] ; b2 = [0, 1] # 自变量的边界
ranges=np.vstack([x1, x2]).T # 生成自变量的范围矩阵
borders=np.vstack([b1, b2]).T # 生成自变量的边界矩阵
"""========================遗传算法参数设置========================="""
NIND = 50 # 种群规模
MAXGEN = 500 # 最大遗传代数
GGAP = 0.8 # 代沟：子代与父代的重复率为(1-GGAP)
selectStyle = 'rws' # 遗传算法的选择方式设为"rws"——轮盘赌选择
recombinStyle = 'xovdp' # 遗传算法的重组方式，设为两点交叉
recopt = 0.9 # 交叉概率
pm = 0.01 # 变异概率
SUBPOP = 1 # 设置种群数为1
maxormin = 1 # 设置标记表明这是最小化目标
"""=======================调用编程模板进行种群进化==================="""
# 调用编程模板进行种群进化，得到种群进化和变量的追踪器以及运行时间
[pop_trace, var_trace, times] = mintemp1(AIM_M, 'aimfuc', PUN_M,
'punishing', ranges, borders, MAXGEN, NIND, SUBPOP, GGAP,
selectStyle, recombinStyle, recopt, pm, maxormin)
"""=========================绘图及输出结果========================="""
# 传入pop_trace进行绘图
gea.trcplot(pop_trace, [['各代种群最优目标函数值'],
['各代种群个体平均适应度值', '各代种群最优个体适应度值']],
['demo_result1', 'demo_result2'])
# 输出结果
best_gen = np.argmin(pop_trace[:, 0]) # 记录最优种群是在哪一代
print('最优的目标函数值为：', np.min(pop_trace[:, 0]))
print('最优的控制变量值为：')
for i in range(var_trace.shape[1]):
print(var_trace[best_gen, i])
print('最优的一代是第',best_gen + 1,'代')
print('用时：', times, '秒')
运行结果有问题