基于DEAP库的python进化算法-6.遗传算法中的约束处理

前言

这一节我们想要探讨如何在遗传算法中处理约束，这部分内容主要是对Coello Coello大神的经典文章《Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art》的再加工以及代码实现。想要更加深入了解的同学，强烈推荐阅读一下这篇文章。

一.罚函数法简介和实现思路

1.罚函数法的基本思想

用遗传算法求解约束优化问题时，罚函数法是最常见的一种有效手段。罚函数法的基本思想是通过在适应度函数中对违反约束的个体施加惩罚，将约束优化问题转化为求解无约束优化问题。

罚函数法在使用中面临的最大问题是如何选取合适的惩罚函数：

• 惩罚系数较大，族群会更加集中在可行域中，而不鼓励向不可行域探索。当惩罚系数过大，容易使算法收敛于局部最优；
• 惩罚系数较小，族群会更积极在不可行域中进行大量探索，一定程度上能帮助寻找全局最优，但也有浪费算力的风险。当惩罚系数过小，算法可能会收敛于不可行解。

在通常的罚函数设计中，有三种惩罚的思路：

• 无论违反约束的程度如何，只要违反了约束，就施加以同样力度的惩罚

• 根据违反约束程度的不同加以惩罚；

• 根据不可行解“修复”的难度或者距离可行域的距离施加惩罚。

2.实现思路

罚函数法主要有以下几种：

• 死亡惩罚（Death Penalty）
  死亡惩罚简单粗暴，对于非可行解，直接将其适应度函数调整到极容易被淘汰的大小（常用归零），使其“死亡”。
• 静态惩罚（Static Penalty）
  顾名思义，在静态惩罚中，通过在适应度函数中整合惩罚项，降低违反约束的个体的适应度，在计算中惩罚系数为常数，不会随着算法迭代变化。
• 动态惩罚（Dynamic Penalty）
  在动态惩罚中，惩罚系数会随着代数变化。通常为保证收敛性，惩罚系数会随着迭代进行而逐渐增大。
• 退火惩罚（Annealing Penalty）
  结合模拟退火的思想，惩罚系数在一定代数内保持稳定，呈阶梯式调整。
• 自适应惩罚（Adaptive Penalty）
  自适应惩罚的思想是从遗传算法迭代的过程中获得反馈，由得到的解质量来判定是要加大还是较小惩罚系数。

二.各类罚函数详解与代码实现

1.死亡惩罚

死亡惩罚就是对所有违反约束的解个体分配一个极差的适应度，例如归零或者在最小化问题中设为一个很大的值。通过这样操作，在后续的选择和变异操作中，不可行解就会被筛选掉。

死亡惩罚是在进化算法中非常流行的惩罚策略，其优点在于实施简单，但是缺点也不少：首先，它只能适用于可行域占解空间比例较大的情况，如果可行域很小，那么初始生成个体很可能没有落于其中的，那么就会全部死亡，使算法陷入停滞；其次，它没有从不可行解中利用任何信息，不能为下一步的进化提供指导，也就是前人的牺牲对于后人完全没有意义。各类研究都表明死亡惩罚是一种比较弱的方法，基本所有其他惩罚方法都能在结果或搜索效率上吊锤它。
利用DEAP提供的装饰器实现死亡惩罚-代码实现

deap.tools.DeltaPenalty(feasibility, delta[, distance])
        这个装饰器为无效的个体返回惩罚适应度，并为有效的个体返回原始适应度值。
        惩罚的适应度由一个常数因子delta加上一个(可选的)距离惩罚组成。如果提供了距离函数，则返回的值将随着个体离开有效区域而增长     
        参数：
            feasibility：返回任何个体的有效性状态函数
            delta:为无效个体返回的常量或常量数组
            distance:返回个体与给定有效点之间距离的函数。距离函数也可以返回一个长度序列，等于目标的数量，
                   	以不同的方式影响多目标适配(可选，默认为0)    
            return 返回求值函数的装饰器
            
decorate(alias, decorator[, decorator[, ...]])
            用指定的装饰器装饰别名，别名必须是当前工具箱中的注册函数
            参数：
                alias：别名
                decorator：一个或多个函数装饰器。
                            如果提供了多个装饰器，则将按顺序应用它们，最后一个装饰器将装饰所有其他装饰器

def feasible(ind):
    # 判定解是否满足约束条件
    # 如果满足约束条件，返回True，否则返回False
    if feasible:
        return True
    return False

## 用装饰器修饰评价函数
toolbox.register('evaluate', evalFit) # 不加约束的评价函数
# 假设时最小化问题，并且已知的最小值远小于1e3
toolbox.decorate('evaluate', tools.DeltaPenalty(feasible, 1e3)) # death penalty
# 这样添加装饰器之后，在feasible返回True的时候，评价函数会返回evalFit的返回值；否则会返回1e3。

2.静态惩罚(Static Penalty)

静态惩罚会将违反约束的程度作为一种考量，在加权后纳入到适应度函数当中。权重在迭代过程当中保持不变，这就是“静态”名字的由来。相比于死亡惩罚绝不允许个体踏出可行域，静态惩罚允许一定程度的由外至内的搜索，因此其搜索效果会相对较好 – 尤其考虑到很多约束优化的最优解会落在可行域边界和顶点上。

3.动态惩罚（Dynamic Penalty）

动态惩罚是对于静态惩罚的一种改进，它的思想是在迭代过程中动态改变惩罚系数，在扩大搜索范围和保证收敛性之间取得动态平衡。

4.退火惩罚（Annealing Penalty）

退火惩罚实际上也是一种动态惩罚，只是它的惩罚系数调整的方式来自于模拟退火算法，以一种类似模拟退火中的temperature schedule的方式来动态调整惩罚系数。

5.自适应惩罚（Adaptive Penalty）

自适应惩罚是对动态惩罚的一种演进，它的主要思路在于从迭代过程中取得反馈，用这个反馈来指导迭代系数的调整。

三.简单测试

用简单的函数来测试各类约束下的GA算法是否能够如预期般工作
目标函数：

准备测试的惩罚方式有:

• 死亡惩罚
• 静态惩罚
• 动态惩罚

1.施加死亡惩罚

#!usr/bin/env python
# -*- coding:utf-8 _*-
"""
@author: liujie
@software: PyCharm
@file: 死亡惩罚.py
@time: 2020/11/26 20:46
"""
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import bernoulli
from deap import creator,tools,base

# random.seed(42)

# 定义问题
creator.create('FitnessMin',base.Fitness,weights=(-1.0,))       # 单变量优化最小值
creator.create('Individual',list,fitness = creator.FitnessMin)

# 个体编码
def uniform(low,up):
    # 均匀分布生成个体
    return [random.uniform(low[0],up[0]),random.uniform(low[1],up[1])]

gen_size = 2
# 两个变量下界
low = [-10] * gen_size
# 两个变量上界
up = [10] * gen_size
# 生成个体
toolbox = base.Toolbox()
toolbox.register('Attr_float',uniform,low,up)
toolbox.register('Individual',tools.initIterate,creator.Individual,toolbox.Attr_float)
toolbox.register('Population',tools.initRepeat,list,toolbox.Individual)
# 生成初始种群
pop_size = 100
pop = toolbox.Population(n=pop_size)
# print(pop)

# 死亡惩罚实现约束
# 评价函数
def eval(ind):
    return (ind[0]**2 + ind[1]**2),

# 死亡惩罚
def feasible(ind):
    # 判定解是否可行
    if (ind[0]**2 + ind[1]**2) < 25:
        return True
    return False

toolbox.register('evaluate',eval)
toolbox.decorate('evaluate',tools.DeltaPenality(feasible,100))
# 在feasible函数满足True时，评价函数会返回eval值，否则返回100

# 注册进化过程中需要的工具：配种选择、交叉、变异
toolbox.register('select',tools.selTournament,tournsize=2)  # 锦标赛选择缺k
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)
# 用数据记录算法迭代过程
stats = tools.Statistics(key= lambda ind : ind.fitness.values)
stats.register('avg',np.mean)
stats.register('std',np.std)
stats.register('min',np.min)
stats.register('max',np.max)

# 创建日志对象
logbook = tools.Logbook()

# 进化迭代，不妨设置如下参数
N_Gen = 50
CXPB = 0.8
MUTPB = 0.2

# 评价初代种群
fitnesses = map(toolbox.evaluate,pop)
for ind,fit in zip(pop,fitnesses):
    ind.fitness.values = fit
record = stats.compile(pop)
logbook.record(gen=0,**record)

# 族群迭代-从第二代开始
for gen in range(1,N_Gen+1):
    # 配种选择
    selectTour = toolbox.select(pop,pop_size*2)
    # 复制个体，供交叉使用
    selectInd = list(map(toolbox.clone,selectTour))
    # print(selectInd)
    # 对选出的个体两两进行交叉，对于被改变的个体，删除其适应度
    for child1,child2 in zip(selectInd[::2],selectInd[1::2]):
        if random.random() < CXPB:
            toolbox.mate(child1,child2)
            del child1.fitness.values
            del child2.fitness.values

    # 变异
    for mutate in selectInd:
        if random.random() < MUTPB:
            toolbox.mutate(mutate)
            del mutate.fitness.values

    # 对于被改变的个体，重新计算其适应度
    invalid_ind = [ind for ind in selectInd if not ind.fitness.valid