这一节我们想要探讨如何在遗传算法中处理约束,这部分内容主要是对Coello Coello大神的经典文章《Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art》的再加工以及代码实现。想要更加深入了解的同学,强烈推荐阅读一下这篇文章。
用遗传算法求解约束优化问题时,罚函数法是最常见的一种有效手段。罚函数法的基本思想是通过在适应度函数中对违反约束的个体施加惩罚,将约束优化问题转化为求解无约束优化问题。
罚函数法在使用中面临的最大问题是如何选取合适的惩罚函数:
在通常的罚函数设计中,有三种惩罚的思路:
无论违反约束的程度如何,只要违反了约束,就施加以同样力度的惩罚
根据违反约束程度的不同加以惩罚;
根据不可行解“修复”的难度或者距离可行域的距离施加惩罚。
罚函数法主要有以下几种:
死亡惩罚就是对所有违反约束的解个体分配一个极差的适应度,例如归零或者在最小化问题中设为一个很大的值。通过这样操作,在后续的选择和变异操作中,不可行解就会被筛选掉。
死亡惩罚是在进化算法中非常流行的惩罚策略,其优点在于实施简单,但是缺点也不少:首先,它只能适用于可行域占解空间比例较大的情况,如果可行域很小,那么初始生成个体很可能没有落于其中的,那么就会全部死亡,使算法陷入停滞;其次,它没有从不可行解中利用任何信息,不能为下一步的进化提供指导,也就是前人的牺牲对于后人完全没有意义。各类研究都表明死亡惩罚是一种比较弱的方法,基本所有其他惩罚方法都能在结果或搜索效率上吊锤它。
利用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。
静态惩罚会将违反约束的程度作为一种考量,在加权后纳入到适应度函数当中。权重在迭代过程当中保持不变,这就是“静态”名字的由来。相比于死亡惩罚绝不允许个体踏出可行域,静态惩罚允许一定程度的由外至内的搜索,因此其搜索效果会相对较好 – 尤其考虑到很多约束优化的最优解会落在可行域边界和顶点上。
动态惩罚是对于静态惩罚的一种改进,它的思想是在迭代过程中动态改变惩罚系数,在扩大搜索范围和保证收敛性之间取得动态平衡。
退火惩罚实际上也是一种动态惩罚,只是它的惩罚系数调整的方式来自于模拟退火算法,以一种类似模拟退火中的temperature schedule的方式来动态调整惩罚系数。
自适应惩罚是对动态惩罚的一种演进,它的主要思路在于从迭代过程中取得反馈,用这个反馈来指导迭代系数的调整。
用简单的函数来测试各类约束下的GA算法是否能够如预期般工作
目标函数:
准备测试的惩罚方式有:
#!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(