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代码进行优化了。后续我将继续学习和挖掘该工具箱的更多深入的用法。希望这篇文章在帮助自己记录学习点滴之余，也能帮助大家！
推荐阅读：
