【优化算法】Python实现面向对象的遗传算法-编程知识

遗传算法

遗传算法(Genetic Algorithm)属于智能优化算法的一种，本质上是模拟自然界中种群的演化来寻求问题的最优解。与之相似的还有模拟退火、粒子群、蚁群等算法。

在具体介绍遗传算法之前，我们先来了解一些知识🧀

DNA：携带有合成RNA和蛋白质所必需的遗传信息的生物大分子，是生物体发育和正常运作必不可少的生物大分子。一般情况下，是以双螺旋结构存在。

现实中的DNA由碱基、脱氧核糖和磷酸双分子层组成，两条脱氧核苷酸链通过碱基间的氢键形成的碱基对相连，形成稳定的结构。碱基由四种，腺嘌呤，鸟嘌呤，胸腺嘧啶，胞嘧啶。

那么如何在计算机中对DNA进行编码？也不需要想的过于复杂，我们只需要表达每个位置上的碱基就行啦！譬如，一条DNA链可以写作：012332100123；每个数字对应一个碱基的映射。为了提高运行速度，我们将其以二进制的形式进行简化表达，即一条DNA链可以看做一串二进制文本：01011101010

个体与种群

我们将每条DNA链看做一个个体，实际上，也就是问题的一个解。譬如，我们寻找映射 $F (X, Y)$ 的最优解，其中一个可能的值 $X = 6$ ，便是一个个体。

种群就是个体的集合。注意，同一种群内部发生信息交换，不同种群之间不会发生信息交换。例如， $X$ 的全集不会与 $Y$ 的全集发生交换，DNA交换只会发生在 $X$ 集合或 $Y$ 集合内部。

遗传：生物的DNA来自于父母，一般情况下由父亲提供X或Y染色体，母亲提供X染色体

假设现在有两个个体：

$1001\ 0001$ 和 $00011001$ 分别作为父亲和母亲，生下了一个新的个体，那么该个体的DNA将由父亲和母亲来决定，例如前四位由父亲提供，后四位由母亲提供，那么该子代个体就是：

$1001\ 1001$

变异：在DNA的某个位置发生了变化

因为是二进制表达的DNA，那么所谓的变化就是某一位由0到1，或是由1到0

自然选择：优胜劣汰，将会选择更加适宜的个体

个体的适宜度，实际上也就是满足函数最优解的值。适宜度越高，该个体在下一轮的自然选择中越容易存活，从而保留自身的DNA。

下面我们将来说一下如何实现一个GA算法~

一、创建属于我们的种群

在一切的开始，我们为种群制定一个规则，比方说这个种群的大小，DNA链的长度~

Population_Nums=200 # 种群有200个个体
DNA_Size=16 # DNA链的长度
Range=[[-3,3],[-3,3]] # 自变量的值域范围# 初始化
import numpy as np 
# pop维度为(n,Population_Nums,DNA_Size)
# 其中n表示有几个种群，也就是自变量
pop=np.random.randint(0,2,size=(len(Range),Population_Nums,DNA_Size))

种群的数量呢决定了收敛的速度，但是也有可能因此陷入局部最优解，并降低运行速度(注意跟收敛速度的区别)

DNA链的长度实际上决定了精度。这句话如何理解呢？我们要来看如何将一串DNA转译成我们需要的信息~

假设有一串DNA链： $10101$ ，我们的映射函数为 $F (X, Y)$ ，其中 $X$ 的值域为 $[- 5, 5]$ ，想一想如何进行转译呢？

为了方便计算和模拟遗传变异，我们采用了二进制作为个体DNA。而我们想要的结果是十进制，那就需要先将二进制转为十进制！ $1*2^{4}+0*2^3+1*2^2+0*2^1+1*2^0=21$ ，这样就来到了十进制。对于种群的基因，只需要除以一个最大值，即 $11111$ ，或者说 $2^n$ ，就可以压缩到区间 $[0, 1]$ ，然后再通过区间匹配到实际值域区间中。

这段写成代码的话，可以是这样：

def decoding(pop):deList=[]for idx,i in enumerate(pop):deList.append(i.dot(2**np.arange(DNA_Size))/float(2**DNA_Size)*(Range[idx][1]-Range[idx][0])+Range[idx][0])return deList

好啦，那么现在我们就需要评估一个个体的适宜度，这也是自然选择中最重要的部分。适宜度越大的个体，越容易在下一轮的选择中存活。

假设函数为：

def F(x,y):return 3*(1-x)**2*np.exp(-(x**2)-(y+1)**2)-10*(x/5-x**3-y**5)*np.exp(-x**2-y**2)-1/3**np.exp(-(x+1)**2-y**2)

假设优化目标是求这个函数的极大值，那么我们的适宜度就应该是个体DNA转为十进制编码后，带入函数的结果，这个结果越大，说明适宜度越高。

def fitnetss(pop):deList=decoding(pop)pred=F(*deList)return (pred-pred.min())+1e-3

后面这个1e-3的实际含义是，让每一个个体都有机会，而不是绝对肯定或绝对否定哪个个体。

二、遗传和变异

在遗传部分，我们设置了一个参数，用来控制遗传发生的比例。毕竟有些个体并没有后代~

在变异部分，我们同样也有一个较小的参数，用来控制变异发生的可能性。

def mutation(pop,rate=1e-3):# 变异将随机发生for i in pop:if np.random.rand()<rate:# 随机一个个体发生变异i[np.random.randint(0,DNA_Size)]^=1

变异的作用是跳出局部最优解，下面是进行变异的三次结果：

(x,y):  -0.03668099860435703 1.499994903802568
(x,y):  -0.013274610833800438 1.6933678801875045
(x,y):  0.05119961805341333 1.4999723732455

而下面是不进行变异的三次结果：

(x,y):  -0.027307929236169315 1.5981612562037268
(x,y):  0.18948037561657305 1.4062327388663727
(x,y):  -0.0962074456338553 1.4998157322296937

可以发现，发生了变异后，结果稳定在[0,1.5]，而不是陷入部分最优解。

遗传过程的算法可以描述如下：

遍历种群中的每个个体，并将该个体A作为父母个体
有一定概率该个体可以随机跟种群中的其他个体B发生基因交换（甚至包括它自己，但这对结果并没有影响，只是降低了遗传概率）
发生基因交换时，随机选择DNA的断点，断点前半部分由个体A提供，后半段由个体B提供

def crossover(pop,rate=0.7):# 注意这里只与自身种群发生变化new_pop=[]for idx in pop.shape[0]:children=[]for father in pop[idx]:if np.random.rand()<rate:child=fathermother=pop[idx][np.random.randn(0,Population_Nums)]# 随机选择发生互换的碱基对choicePoint=np.random.randn(0,DNA_Size)child[choicePoint:]=mother[choicePoint:]# 发生变异children.append(mutation(child))return chidren

三、自然选择

这部分将会根据个体的适宜度分配权值，决定该个体基因出现在下一轮概率。

def select(pop,fitness):pop_s=[]for i in pop:pop_s.append(i[np.random.choice(np.array(Population_Nums),size=Population_Nums,replace=True,p=(fitness)/(fitness.sum()))])return pop_s

四、基于面向对象的遗传算法

现在，我们就要将这些东西封装成一个类啦，提高复用性和稳定性。

首先是构造函数，就是先写入一些常量。

class GA(object):def __init__(self,popN=2000,DNA_Size=16,Epochs=500,crossRate=0.8,mutationRate=0.005):self.popN=popN # 种群数量self.DNA_Size=DNA_Size # DNA长度self.Epochs=Epochs # 迭代次数self.crossRate=crossRate # 交叉遗传概率self.mutationRate=mutationRate  # 变异概率self.Range=None # 输入数据的值域# 例如：[[-3,3],[2,5],[1,9]] 这表示第一个变量的值域是[-3,3]，第二个是[2,5]self.plot_=[] # 保留每轮的最优值self.bestScore=None # 最佳得分self.best=None # 最佳个体

然后，需要提供一个输入函数的接口：

   def fit_function(self,f,range):self.f=fself.Range=range# 初始化种群self.pop=np.random.randint(0,2,size=(len(range),self.popN,self.DNA_Size))self.plot_=[]

解码方法：

   def decoding(self):deList = []for idx, i in enumerate(self.pop):deList.append(i.dot(2 ** np.arange(self.DNA_Size)) / float(2 ** self.DNA_Size) * (self.Range[idx][1] - self.Range[idx][0]) + self.Range[idx][0])return deList

适应值：

    def fitness(self):deList=self.decoding()pred=self.f(*deList)return (pred-pred.min())+1e-3

变异：

    def mutation(self,pop):if np.random.rand()<self.mutationRate:pop[np.random.randint(0,self.DNA_Size)]^=1return pop

交叉遗传：

    def crossover(self):for idx in range(self.pop.shape[0]):for _,father in enumerate(self.pop[idx]):child = fatherif np.random.rand()<self.crossRate:mother=self.pop[idx][np.random.randint(0,self.popN)]crossPoint=np.random.randint(0,self.DNA_Size)child[crossPoint:]=mother[crossPoint:]self.pop[idx][_]=self.mutation(child)

自然选择：

    def select(self,fitness):pops=[]for i in self.pop:pops.append(i[np.random.choice(self.popN,size=self.popN,replace=True,p=fitness/(fitness.sum()))])return pops

打印信息：

    def getInfo(self):print('最优参数为： ',[i for i in self.best])print("最优结果为: ",self.bestScore)

提供一个训练接口和绘图接口供使用者调用：

    def train(self,plot=False):for _ in range(self.Epochs):self.crossover() # 交叉变异f=self.fitness() # 计算适宜度max_fit = np.argmax(f) # 获取最大适宜度下标k=[(i[max_fit].dot(2**np.arange(self.DNA_Size))/float(2**self.DNA_Size))*(self.Range[idx][1]-self.Range[idx][0])+self.Range[idx][0] for idx,i in enumerate(self.pop)] # 获取最佳个体的十进制值bs=self.f(*k) # 计算该值的适宜度# 请注意，适宜度并不代表函数结果，适宜度是一个相对的值# 记录全局最优结果if self.bestScore==None or bs>self.bestScore:self.bestScore=bsself.best=kself.pop = np.array(self.select(f)) # 自然选择if plot:self.plot_.append(bs)self.getInfo()def plot(self):if self.plot_==[]:passplt.plot([i for i in range(self.Epochs)],self.plot_)plt.xlabel("Epochs")plt.ylabel("BestValue")plt.show()

最终的结果如下：

在这里插入图片描述

可以看到在25个Epoch左右就开始收敛了。

完整代码

import numpy as np
import matplotlib.pyplot as pltclass GA(object):def __init__(self,popN=2000,DNA_Size=16,Epochs=500,crossRate=0.8,mutationRate=0.005):self.popN=popNself.DNA_Size=DNA_Sizeself.Epochs=Epochsself.crossRate=crossRateself.mutationRate=mutationRateself.Range=Noneself.plot_=[]self.bestScore=Noneself.best=Nonedef fit_function(self,f,range):self.f=fself.Range=rangeself.pop=np.random.randint(0,2,size=(len(range),self.popN,self.DNA_Size))self.plot_=[]def decoding(self):deList = []for idx, i in enumerate(self.pop):deList.append(i.dot(2 ** np.arange(self.DNA_Size)) / float(2 ** self.DNA_Size) * (self.Range[idx][1] - self.Range[idx][0]) + self.Range[idx][0])return deListdef fitness(self):deList=self.decoding()pred=self.f(*deList)return (pred-pred.min())+1e-3def mutation(self,pop):if np.random.rand()<self.mutationRate:pop[np.random.randint(0,self.DNA_Size)]^=1return popdef crossover(self):for idx in range(self.pop.shape[0]):for _,father in enumerate(self.pop[idx]):child = fatherif np.random.rand()<self.crossRate:mother=self.pop[idx][np.random.randint(0,self.popN)]crossPoint=np.random.randint(0,self.DNA_Size)child[crossPoint:]=mother[crossPoint:]self.pop[idx][_]=self.mutation(child)def select(self,fitness):pops=[]for i in self.pop:pops.append(i[np.random.choice(self.popN,size=self.popN,replace=True,p=fitness/(fitness.sum()))])return popsdef getInfo(self):print('最优参数为： ',[i for i in self.best])print("最优结果为: ",self.bestScore)def train(self,plot=False):for _ in range(self.Epochs):self.crossover()f=self.fitness()max_fit = np.argmax(f)k=[(i[max_fit].dot(2**np.arange(self.DNA_Size))/float(2**self.DNA_Size))*(self.Range[idx][1]-self.Range[idx][0])+self.Range[idx][0] for idx,i in enumerate(self.pop)]bs=self.f(*k)if self.bestScore==None or bs>self.bestScore:self.bestScore=bsself.best=kself.pop = np.array(self.select(f))if plot:self.plot_.append(bs)self.getInfo()return self.bestdef plot(self):if self.plot_==[]:passplt.plot([i for i in range(self.Epochs)],self.plot_)plt.xlabel("Epochs")plt.ylabel("BestValue")plt.show()if __name__ == '__main__':ga=GA(popN=200,DNA_Size=16,Epochs=200)def F(x, y):return 3 * (1 - x) ** 2 * np.exp(-(x ** 2) - (y + 1) ** 2) - 10 * (x / 5 - x ** 3 - y ** 5) * np.exp(-x ** 2 - y ** 2) - 1 / 3 ** np.exp(-(x + 1) ** 2 - y ** 2)Range = [[-3, 3], [-3, 3]]ga.fit_function(F,Range)ga.train(True)ga.plot()