基于遗传算法（deap）的非线性函数寻优与编码方式浅谈

前言

  借非线性函数寻优问题总结遗传算法中个体编码的方式与在DEAP库中的写法。

问题

  问题为一个非线性函数在[-4.5,4.5]内求最小值。
  该问题中有两个自变量，即一个个体（染色体）中有2个基因，这样可以推导至N基因组成的个体编码方式。
  遗传算法中常见的编码方式为二进制编码法与实数编码法，先贴出两种方式对应的 完整代码。

二进制编码

import numpy as np
from deap import base, tools, creator, algorithms
import random

#定义问题
creator.create('FitnessMin',base.Fitness,weights=(-1.0,))#单变量，求最小值
creator.create('Individual',list,fitness = creator.FitnessMin)#创建individual类

GENE = 48#基因长度
toolbox  = base.Toolbox()
toolbox.register('Binary',random.randint,0,1)
toolbox.register('individual',tools.initRepeat,creator.Individual,toolbox.Binary,GENE)
#print(toolbox.individual())#打印一个个体检查

low = [-4.5,-4.5]
up = [4.5,4.5]
def decode(ind,low=low,up=up,genLen=GENE):
    #染色体里包含两个解，所以将染色体一分为二再解码
    subGeneLen = int(GENE/2)
    x1 = int(''.join(str(_) for _ in ind[:subGeneLen]),2)
    x2 = int(''.join(str(_) for _ in ind[subGeneLen:]),2)
    x1Rescaled = low[0] + x1 * (up[0] - low[0])/(2**subGeneLen - 1)
    x2Rescaled = low[1] + x2 * (up[1] - low[0])/(2**subGeneLen - 1)
    return x1Rescaled,x2Rescaled

def evaluate(ind):
    x1,x2 = decode(ind)
    f = (1.5-x1+x1*x2)*(1.5-x1+x1*x2) + \
        (2.25-x1+x1*x2*x2)*(2.25-x1+x1*x2*x2) +\
        (2.625-x1+x1*x2*x2*x2)*(2.625-x1+x1*x2*x2*x2)
    return f,
toolbox.register('evaluate',evaluate)

N_POP = 200#种群内个体数量，参数过小，搜索速度过慢
toolbox.register('population',tools.initRepeat,list,toolbox.individual)
pop = toolbox.population(n=N_POP)
#for ind in pop:#打印一个种群检查
#   print(ind)

toolbox.register('select',tools.selTournament,tournsize=3)
toolbox.register('mate',tools.cxTwoPoint)
toolbox.register('mutate',tools.mutFlipBit,indpb=0.5)

NGEN = 200#迭代步数，参数过小，在收敛之前就结束搜索
CXPB = 0.8#交叉概率，参数过小，族群不能有效更新
MUTPB = 0.2#突变概率，参数过小，容易陷入局部最优

invalid_ind = [ind for ind in pop if not ind.fitness.valid]
fitnesses = toolbox.map(toolbox.evaluate,invalid_ind)
for ind ,fitness in zip(invalid_ind,fitnesses):
    ind.fitness.values = fitness

for gen in range(NGEN):
    offspring = toolbox.select(pop,N_POP)#先选择一次
    offspring = [toolbox.clone(_) for _ in offspring]#再克隆一次

    for ind1,ind2 in zip(offspring[::2],offspring[1::2]):#交叉
        if random.random() < CXPB:
            toolbox.mate(ind1,ind2)
            del ind1.fitness.values
            del ind2.fitness.values
    for ind in offspring:#变异
        if random.random() <MUTPB:
            toolbox.mutate(ind)
            del ind.fitness.values

    invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
    fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
    for fitness, ind in zip(fitnesses, invalid_ind):
        ind.fitness.values = fitness

    combinedPop = pop + offspring#将子代与父代结合起来
    pop = tools.selBest(combinedPop,N_POP)#再从子代与父代的结合中选择出适应度最高的一批作为新的种群


bestInd = tools.selBest(pop,1)[0]#选择出最好的个体编号
bestFit = bestInd.fitness.values[0]#最好的个体适应度
print('best solution: '+str(decode(bestInd)))
print('best fitness: '+str(bestFit))

实数编码

import numpy as np
from deap import base, tools, creator, algorithms
import random

#定义问题
creator.create('FitnessMin',base.Fitness,weights=(-1.0,))#单变量，求最小值
creator.create('Individual',list,fitness = creator.FitnessMin)#创建individual类

low = [-4.5,-4.5]
up = [4.5,4.5]
def genInd(low,up):
    return[random.uniform(low[0],up[0]),random.uniform(low[1],up[1])]
toolbox  = base.Toolbox()
toolbox.register('genInd',genInd,low,up)
toolbox.register('individual',tools.initIterate,creator.Individual,toolbox.genInd)
#print(toolbox.individual())#打印一个个体检查


def evaluate(ind):
    f = (1.5-ind[0]+ind[0]*ind[1])*(1.5-ind[0]+ind[0]*ind[1]) + \
        (2.25-ind[0]+ind[0]*ind[1]*ind[1])*(2.25-ind[0]+ind[0]*ind[1]*ind[1]) +\
        (2.625-ind[0]+ind[0]*ind[1]*ind[1]*ind[1])*(2.625-ind[0]+ind[0]*ind[1]*ind[1]*ind[1])
    return f,
toolbox.register('evaluate',evaluate)

N_POP = 200#种群内个体数量，参数过小，搜索速度过慢
toolbox.register('population',tools.initRepeat,list,toolbox.individual)
pop = toolbox.population(n=N_POP)
#for ind in pop:#打印一个种群检查
#   print(ind)

toolbox.register('select',tools.selTournament,tournsize=3)
toolbox.register('mate',tools.cxSimulatedBinaryBounded, eta=20, low=low, up=up)
toolbox.register('mutate',tools.mutPolynomialBounded, eta=20, low=low, up=up, indpb=0.2)

NGEN = 200#迭代步数，参数过小，在收敛之前就结束搜索
CXPB = 0.8#交叉概率，参数过小，族群不能有效更新
MUTPB = 0.2#突变概率，参数过小，容易陷入局部最优

invalid_ind = [ind for ind in pop if not ind.fitness.valid]
fitnesses = toolbox.map(toolbox.evaluate,invalid_ind)
for ind ,fitness in zip(invalid_ind,fitnesses):
    ind.fitness.values = fitness

for gen in range(NGEN):
    offspring = toolbox.select(pop,N_POP)#先选择一次
    offspring = [toolbox.clone(_) for _ in offspring]#再克隆一次

    for ind1,ind2 in zip(offspring[::2],offspring[1::2]):#交叉
        if random.random() < CXPB:
            toolbox.mate(ind1,ind2)
            del ind1.fitness.values
            del ind2.fitness.values
    for ind in offspring:#变异
        if random.random() <MUTPB:
            toolbox.mutate(ind)
            del ind.fitness.values
     
    invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
    fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
    for fitness, ind in zip(fitnesses, invalid_ind):
        ind.fitness.values = fitness
    combinedPop = pop + offspring#将子代与父代结合起来
    pop = tools.selBest(combinedPop,N_POP)#再从子代与父代的结合中选择出适应度最高的一批作为新的种群


bestInd = tools.selBest(pop,1)[0]#选择出最好的个体编号
bestFit = bestInd.fitness.values[0]#最好的个体适应度
print('best solution: '+str(bestInd))
print('best fitness: '+str(bestFit))

二进制编码在DEAP程序中写法

  先从复杂一点的二进制编码记录。二进制编码法的显性特点有：需要自定义基因长度，需要为其单独编码解码函数。
  首先确定基因长度，先从单个解开始分析。自变量的取值范围为[-4.5,4.5]，总范围长度就是4.5-(-4.5)=9；希望答案的精确度为6位，则所有的可取值为9*1e6=9000000，那么基因的位数n一定符合2的n次方-1>9000000,2的23次方为8388608<9000000，2的24次方为16777216>9000000，因此单个基因长度应该为24位，一个个体有2个基因，因此一个个体的长度应该是48位。

GENE = 48#基因长度
toolbox  = base.Toolbox()
toolbox.register('Binary',random.randint,0,1)
toolbox.register('individual',tools.initRepeat,creator.Individual,toolbox.Binary,n=GENE)

  二进制编码，每一位非0即1，因此使用random.randint()函数，将0和1当做边界，就可以随机生成每一位的码。

toolbox.register('Binary',random.randint,0,1)

个体生成函数tools.initRepeat用法

  在个体生成中二进制编码使用的方法为tools.initRepeat()，和实数编码方式使用的方法不同（很重要）。
  复习一下DEAP库中工具箱注册方式：下句代码中，individual()是被注册成的函数名称，今后的个体生成都靠调用toolbox.individual()来实现；而individual()方法的本质是tools.initRepeat()方法；creator.Individual、toolbox.Binary、n=GENE则是tools.initRepeat()方法的三个传参。
  第一个传参creator.Individual的作用为令该个体确定为用于求解最大值问题还是最小值问题，即确定其适应度fitness.value是正值还是负值。
  第二个传参toolbox.Binary的作用为确定该个体的每位生成方法。
  第三个传参n=GENE即为该个体的位数。函数名repeat顾名思义，是将某种方法形式重复进行一定次数，也就是将每个位生成的方法重复使用n次。

toolbox.register('individual',tools.initRepeat,creator.Individual,toolbox.Binary,n=GENE)

解码函数写法

  解码函数的写法与前文一元函数寻优文章里的原理相同，主要区别为个体（染色体）中有多个解，需要先将个体拆分为单个基因，再解码。
  subGeneLen为单个基因的长度（24），第一个基因为ind[:subGeneLen]，第二个基因为ind[subGeneLen:]。同理，如果有n个基因，第一个基因为ind[:subGeneLen]，第二个基因为ind[subGeneLen:subGeneLen*2]，第三个基因为ind[subGeneLen*2:subGeneLen*3]…第n个基因为ind[subGeneLen*(n-1):]。

low = [-4.5,-4.5]
up = [4.5,4.5]
def decode(ind,low=low,up=up,genLen=GENE):
    #染色体里包含两个解，所以将染色体一分为二再解码
    subGeneLen = int(GENE/2)
    x1 = int(''.join(str(_) for _ in ind[:subGeneLen]),2)
    x2 = int(''.join(str(_) for _ in ind[subGeneLen:]),2)
    x1Rescaled = low[0] + x1 * (up[0] - low[0])/(2**subGeneLen - 1)
    x2Rescaled = low[1] + x2 * (up[1] - low[0])/(2**subGeneLen - 1)
    return x1Rescaled,x2Rescaled

评价函数写法

  评价函数较为简单，将解码函数返回的两个值赋值给两个变量代入原方程即可。

def evaluate(ind):
    x1,x2 = decode(ind)
    f = (1.5-x1+x1*x2)*(1.5-x1+x1*x2) + \
        (2.25-x1+x1*x2*x2)*(2.25-x1+x1*x2*x2) +\
        (2.625-x1+x1*x2*x2*x2)*(2.625-x1+x1*x2*x2*x2)
    return f,

实数编码在DEAP程序中写法

  实数编码法相较于二进制编码法更简洁，实数编码法通常直接构造一个列表（list）作为个体使用。
  该问题的个体由两个基因组成，因此每个列表中应当有两个元素，其形式应为list[x1,x2]。
  首先用两个列表记录x1和x2的上下界，方便之后为随机函数赋值。

low = [-4.5,-4.5]
up = [4.5,4.5]

  然后定义生成个体的方法。两个实数，也就是x1，x2,都是在各自的取值范围内随机生成的，因此使用均匀分布随机函数random.uniform()来生成。random.uniform()方法有两个传参，分别对应变量下界与上界，返回值为边界内等概率选取的实数值。

def genInd(low,up):
    return[random.uniform(low[0],up[0]),random.uniform(low[1],up[1])]

  自己定义的方法不能直接用于deap的个体生成程序中，因此还需要对其在工具箱中进行一次注册：

toolbox.register('genInd',genInd,low,up)

  完整代码如下，可以看到个体生成函数不是initRepeat而是initIterate。

low = [-4.5,-4.5]
up = [4.5,4.5]
def genInd(low,up):
    return[random.uniform(low[0],up[0]),random.uniform(low[1],up[1])]
toolbox  = base.Toolbox()
toolbox.register('genInd',genInd,low,up)
toolbox.register('individual',tools.initIterate,creator.Individual,toolbox.genInd)

个体生成函数 tools.initIterate用法

  工具箱注册方式不再重复说明。主要分析tools.initIterate()与tools.initRepeat()的不同。
  tools.initIterate()可以看出其参数中没有n，也就是说该方法无法确定个体的具体长度（不做重复操作），其余的地方与tools.initRepeat()无区别，个人认为tools.initIterate()从用法上讲就是不做基因位数重复操作的tools.initRepeat()。

toolbox.register('individual',tools.initIterate,creator.Individual,toolbox.genInd)

评价函数写法

  对于每个个体ind而言，ind[0]就是x1，ind[1]就是x2，代入到方程中即可。

def evaluate(ind):
    f = (1.5-ind[0]+ind[0]*ind[1])*(1.5-ind[0]+ind[0]*ind[1]) + \
        (2.25-ind[0]+ind[0]*ind[1]*ind[1])*(2.25-ind[0]+ind[0]*ind[1]*ind[1]) +\
        (2.625-ind[0]+ind[0]*ind[1]*ind[1]*ind[1])*(2.625-ind[0]+ind[0]*ind[1]*ind[1]*ind[1])
    return f,
toolbox.register('evaluate',evaluate)

  当然也可以先存进临时变量中，这样代码与二进制编码法的代码相似度更高些。

def evaluate(ind):
    x1 = ind[0]
    x2 = ind[1]
    f = (1.5-x1+x1*x2)*(1.5-x1+x1*x2) + \
        (2.25-x1+x1*x2*x2)*(2.25-x1+x1*x2*x2) +\
        (2.625-x1+x1*x2*x2*x2)*(2.625-x1+x1*x2*x2*x2)
    return f,

不同编码方式下变异方式的选择

  从代码中可以看到不同编码方式的遗传算法代码使用了不同的交叉变异方式，经过一段时间项目工程测试与加深理解后，得到以下结论：
  二进制编码方式的变异方式选择较为自由，而带边界的实数编码方法对变异方式选择有特定要求。
  二进制编码法中，当基因长度确定了之后，变量的边界和精度就已经确定完毕，无论在基因序列中做怎样的切片与位变换操作，都不会对变量的边界产生影响。
  实数编码中，基因确实是在边界范围内随机生成的，原生基因本身不可能出现越界的情况。
  但是在实数变异后，基因有相当概率跨越边界，使得其值超出预定范围。此种错误的直观结果是：能够得到运算结果，但是参与运算的数据会出现明显的数值膨胀情况。
  因此，综合个人的学习情况认为，在带边界的实数编码法中，如果实数是浮点型，应当使用DEAP库中的有界多项式突变；如果实数是整形，应当使用DEAP库中的均匀整数突变。其余在传参中与边界无关的变异函数，经过一些测试后都发现不适合带边界的实数编码法。

总结

  单纯从数值计算的效率上看，二进制编码法与实数编码法相比多了解码的步骤，运算复杂度是相对增加的。
  但是二进制编码可以人为设定计算的精确程度，相比之下灵活程度更高。
  对程序员来讲，实数编码法较为简易，使用起来方便。
  此外，对于自变量的取值范围，可以并不当做“约束条件”来看待，因为无论有没有其他约束条件，自变量的取值范围都是在个体生成时需要确定的。如果自变量之间除了组成函数求值外，还存在其他方程式关系，这些方程式才是真正的“约束条件”，是需要写在惩罚函数中的。