遗传算法介绍并附上Python代码

Posted 2020-11-26 alittlecomputer

以前在学校的时候，写了一个博客介绍遗传算法，并通过Matlab实现了该算法。但是很多读者反馈说代码运行不起来。

我因为现在没有Matlab工具了，而且Matlab还是挺贵的，所以还是用Python实现以下遗传算法，供大家学习。

写这个博客主要是因为我在学校学习遗传算法的时候，没学会，只是大概知道，但是完全用不起来。但是这个算法还是很有用的，所以还是想着记录一下，让在课堂上没学会的同学，可以在这里学会，并且可以着手。

1、遗传算法介绍

遗传算法，模拟达尔文进化论的自然选择和遗传学机理的生物进化构成的计算模型，一种不断选择优良个体的算法。谈到遗传，想想自然界动物遗传是怎么来的，自然主要过程包括染色体的选择，交叉，变异（不明白这个的可以去看看生物学），这些操作后，保证了以后的个体基本上是最优的，那么以后再继续这样下去，就可以一直最优了。

2、解决的问题

先说说自己要解决的问题吧，遗传算法很有名，自然能解决的问题很多了，在原理上不变的情况下，只要改变模型的应用环境和形式，基本上都可以。但是遗传算法主要还是解决优化类问题，尤其是那种不能直接解出来的很复杂的问题，而实际情况通常也是这样的。

本部分主要为了了解遗传算法的应用，选择一个复杂的二维函数来进行遗传算法优化，函数显示为y=10*sin(5*x)+7*abs(x-5)+10,这个函数图像为：

 

 Python实现代码：

import numpy as np
import matplotlib.pyplot as plt
x=np.arange(0,10,0.001)
y=10*np.sin(5*x) + 7*np.abs(x-5) + 10
plt.plot(x,y)
plt.show()

View Code

 怎么样，还是有一点复杂的吧，当然你还可以任意假设和编写，只要符合就可以。那么现在问你，要你一下求出最大值你能求出来吗？（这个貌似可以，很容易写出来，求个导数然后，----如果再复杂一点估计就不行了）这类问题如果用遗传算法或者其他优化方法就很简单了，为什么呢？说白了，其实就是计算机太笨了，同时计算速度又超快，举个例子吧，我把x等分成100万份，再一下子都带值进去算，求出对应的100万个y的值，再比较他们的大小，找到最大值不就可以了吗，很笨吧，人算是不可能的，但是计算机可以。而遗传算法也是很笨的一个个搜索，只不过加了一点什么东西，就是人为的给它算的方向和策略，让它有目的的算，这也就是算法了。

3、如何开始

不明白遗传算法的会问怎么开始呢？恩，其实每个算法都有自己的开始方式，遗传算法也是，首先是选择个体了。我们知道一个种群中可能只有一个个体吗？不可能吧，肯定很多才对，这样相互结合的机会才多，产生的后代才会多种多样，才会有更好的优良基因，有利于种群的发展。那么算法也是如此，当然个体多少是个问题，一般来说20-100之间我觉得差不多了。那么个体究竟是什么呢？在我们这个问题中自然就是x值了。其他情况下，个体就是所求问题的变量，这里我们假设个体数选100个，也就是开始选100个不同的x值，不明白的话就假设是100个猴子吧。好了，现在有了100个猴子组成的一个种群，那么这个种群应该怎么发展才能越来越好？说到这，我们想想，如何定义这个越来越好呢？这个应该有一个评价指标吧。在我们的这个问题中，好像是对应的y值越大越好是吧。我们甚至可以给他们排个名来决定哪些好哪些不好。我们把这个叫做对于个体的适应度，这应该算是算法的后半部分才对。

4、编码

好，什么是编码？其实编码就是把自变量（x）换一下形式而已，在这个形式下，更容易操作其他过程（比如交叉，变异什么的）而已。举个例子吧，加入我们取x=1, 2, 3，我们可以把x编码成x=a, b, c,我就说123对应就是abc，为什么要这样做呢？比如问题里面你能够获取的都是abc的组合之类的，那么这样编码以后，你就可以再返回去成123来操作了。一般的编码都是些什么呢？二进制编码，自然数编码，矩阵编码。。等等，很多，不详细写了。而用的最多的可以说是二进制编码，感觉这和人体DNA，基因的排列很相似。想想DNA怎么排的？不就是在两条长链上一对一排的吗？那么什么是二进制编码？很简单，就是1,0,1,0对应的来组合排列而已。比如：1100100010, 0011001001等等，这些都是位数长度为10的二进制编码。再想想1在计算机的二进制形式是什么？如果八位来表示的话，是不是就是0000 0001；8是不是就是0000 1000；以此类推，那么我们这里也是这样，把对应的x值换算成这种编码形式，我们这里可以看到x的范围是0-10吧，如果按照计算机这样的方式是不是到0000 1010这里就完事了？想想这样多短，前面四位都没有用上多浪费呀，那么要想都用上怎么办呢？也很简单，我们把0000 0001不认为是1不就可以了吗？因为1111 1111是255，那么如果说每一份为1/255的话，那么0000 0001不就是1/255了吗？这个时候1怎样表示了？不就是1111 1111吗？好了我们把范围扩大一些吧，每一份不是1/255,而是1/255*10,那么这个时候最大值是多少？是不是10，恩，这样x编码的范围不就在0-10之间了吗？这里就又有问题了，想想这样做的话，x的最小精度是多少？就是1/255*10.虽然很小，但是在0-1/255*10之间的数你能不能取到？无论如何都取不到吧。那么就又来一个问题，怎样去扩大这个精度呢？如果要保持0-10不变的话，只能增加位数了，把9位编码，10位，20位，100位，哇，够大了吧，变成100个0,1组合，很恐怖吧，事实上，究竟是多少要视情况而定，一般20位左右感觉就可以了，虽然说越大越好，但是太大了消耗内存，速度慢了，不值得。本题中，我们设置它为一个变量，先暂时取为10来实验。好了，如果还不明白为什么要编码看下面的吧。知道了交叉与变异，你就知道了。

5、关于交叉与变异

先说变异，什么是变异？简单，基因发生突变就叫变异，有了编码的概念，那就在编码基础上来说变异。首先就讲最简单的变异，就是个体的变异。现在以10位长的编码来说，比如把x=3编码一下，随便假设为11000 10010吧，好了，在变异操作时，假设第5位变异了（说一下变异就是一位或者多位1或0变成0或1，也只能0,1之间变，没办法啊），那么这个时候变成什么了？是不是11001 10010再反编码回去成x是多少呢？那肯定不是3了，变了呀，是多少是肯定可以反算回去的，这里懒得算了，就假设为3.213吧，发没发现，这样一来，x是不是变了？既然变了就好啊，带到原函数（适应度函数）里面比较这两个x值对应的哪个y值大一些，如果后面变异后的大一些是不是就是说产生了好的变异啊，就可以在下一次个体选择的时候选择它了。那么想想很多x来一起变异会怎么样呢？肯定会生成很多很多的解吧，反复这么做会怎么样呢？只要每次都保留最优解的话，我来个循环100万次，也总能找到解吧，当然这么多次得花多久，也不太合适，这还只是一个点位在进行变异，如果每次我让多个点位变异呢？哇，又不可思议了，变化更大了吧。当然，变异不止如此，更多的去看专业的论文吧。知道了变异是干什么的，剩下的都好说了，好了，这还是变异，想想自然界遗传中除了变异还有什么？交叉吧，那么交叉又是什么？

学过生物的都知道，动物交配时，部分染色体干什么了？是不是交叉了？就是把相应部分的基因交换了，你的给我，我的给你，很有爱吧。再以编码为例，比如现在随便从100个x值中选取两个吧，假设正好选中了x=3和4，对应的编码假设是11001 10101 和00101 01011，那么怎么交叉呢？我们知道每次交叉的染色体通常是一块一块的，恩，这里在算法设计上也来一块一块的吧。比如说就把位置在2,3,4号的编码给整体交叉了，那么x=3对应的位置是100吧，x=4对应的位置是010吧，好，交换以后x=3对应的位置就变成了010，x=4对应的位置就变成100，加回去就变成什么了？x=3是不是就是10101 10101，x=4是不是就是01001 01011了。而现在，把他们再反编码回去还是x=3和x=4吗？显然又不是了吧（当然也有小概率是一样的吧，很小）。那是什么？不想算，还是假设吧，假设为3.234和4.358吧，好了新的个体是不是又来了？恩，同理，带到适应度函数里面去吧，再取优秀个体，完事。同样，有些专门研究这种算法的开发出来各种各样的交叉方式，什么一个个体的前3个与后一个个体的后三个交叉，中间几位来交叉等等，总之就是生产新个体。而这样做的目的在哪呢？无非是三个字，随机性，充分保证生产新个体具有随机性，你说你的x=3变异后为3.2,3.2距离3那么近，在一些存在局部最优解的问题上就永远跳不出局部最优解，相反，你的x=1一下子变异成了x=5，哇，好大的变化，一下从这头到了那头，这对于算法的广阔搜索能力来说是非常好的。

交叉前：   交叉后：

讲完了这部分，现在知道了为什么要编码了吧？如果你不编码，你说你想要你的x=3怎么去变异？怎么去交叉？当然也不是没有方法，比如你生成一个小的随机数加到x=3上，但是你想想这两种方法哪一个更具有随机性、普遍性？显然的。而更多的时候交叉与变异是在一起操作的，先交叉，再变异是普遍遗传算法的操作步骤。

 

6、关于选择的问题

说完了上面的部分，再说说选择吧，选择是什么？就是优胜劣汰。好的留下来，差的走人，在自然界中直接gg了是吧。不停地选择使的种群一直朝着较好的方向行走。

对应到本问题来说，遗传算法的选择是什么样子？在前面说到，每次交叉或者变异是不是产生了新的个体？如果这些个体都保留下来的话，种群是要爆炸的，第一次循环可能有100个x，第二次循环就200个个体，再来那么10万次，哇哦，多少了，好多。显然不可能吧。而且在算法里面，我们还规定的是每次循环都必须保证都是100个个体。那么必须在200个个体中剔除100个吧。好了，问题来了，如何剔除呢？有人说很简单，排名吧，取前100号x不就可以了吗？排名这个东西真的好吗？我就不信，凭什么差一点的不能选上，搞不好在下一次变异中一下子冲到了第一呢？这个问题在选择上也有一些对应的规则，最通用的就是轮盘赌法，简单来说就是一种概率选择法（当然还有许多其他的方法，感兴趣的自己搜相关的文献吧，我也没用过）。什么是轮盘赌法呢？就是把对应所有y值（适应度函数值）加起来，再用各自的y值去除以这个sum值，这样是不是谁的概率大谁的概率小就很清楚了？然后再随机生成一个0-1的概率值p，谁在p的范围里面是不是就选择谁，比如说x=3时在100个x中y的值最大，那么选择它的概率是不是就最大，比如说是0.1（0.1小吗？不小了好吧，想想其他的会是什么，都比0.1小，那么从概率上讲，选100次的话，是不是就有10次选择x=3，其他的都不足10次是吧，那么在下一次100个种群个体中就有10个x=3了，再来一回可能就有20个x=3的个体了。再就是30个，最后就只剩下100个x=3的个体了，它自己在哪里交叉变异已经没有意义了，如果到了这个时候就意味着这个算法可以结束了）。再详细点，如下图所示吧：现在要在下面三个大类中选取100个x个体，轮盘赌转100次以后，是不是个体数落在s3中的个体多一些，选择的原理就是这样，再不明白直接后面的程序吧。

 7、还差点什么呢？

至此，感觉也差不多了吧，选择完后再重复上述步骤交叉，变异等等，那么什么时候是个头了？很简单，办法就是迭代次数，迭代10次看一下结果，20次，看一下结果，30次，40次，100次，当次数达到一定程度以后，优秀的个体越来越多，大都集中在最优解附近，即使变异或者交叉了也是在这个最优解附近，没有影响的。在下一次选择就又变回来了，那么至此就真的结束了。比如说先来结果吧，该问题按我的思路做完后，迭代100次变成什么样子了？看图：

 

 

代码：

import numpy as np
import random
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
import heapq
import time

# 求染色体长度
def getEncodeLength(decisionvariables, delta):
    # 将每个变量的编码长度放入数组
    # 这个要看有几个自变量，假如自变量的个数是3，那么返回的lengths的元素就有3个
    lengths = []
    for decisionvar in decisionvariables:
        uper = decisionvar[1]
        low = decisionvar[0]
        # res返回一个数组
        res = fsolve(lambda x:((uper-low)/delta-2**x+1), 30)
        # ceil()向上取整
        length = int(np.ceil(res[0]))
        lengths.append(length)
    return lengths


def getinitialPopulation(length, populationSize):
    # 根据DNA长度和种群大小初始化种群
    chromsomes = np.zeros((populationSize, length), dtype=np.int)
    for popusize in range(populationSize):
        # np.random.randint产生[0,2)之间的随机整数，第三个数表示随机数的数量
        chromsomes[popusize, :] = np.random.randint(0, 2, length)
    return chromsomes

# 染色体解码得到表现形的解
def getDecode(population,encodeLength,decisionVariables, delta):
    # 把染色体计算成数字
    populationsize = population.shape[0]
    length = len(encodeLength)
    decodeVariables = np.zeros((populationsize, length), dtype=np.float)
    for i,populationchild in enumerate(population):
        start = 0
        for j, lengthchild in enumerate(encodeLength):
            power = lengthchild - 1
            decimal = 0
            for k in range(start, start+lengthchild):
                decimal += populationchild[k] * (2**power)
                power = power - 1
            start = lengthchild
            lower = decisionVariables[j][0]
            uper = decisionVariables[j][1]
            decodevalue = lower+decimal*(uper-lower)/(2**lengthchild-1)
            decodeVariables[i][j] = decodevalue
    return decodeVariables

# 得到每个个体的适应度值及累计概率
def getFitnessValue(func, decode):
    popusize, decisionvar = decode.shape
    fitnessValue = np.zeros((popusize, 1))
    for popunum in range(popusize):
        fitnessValue[popunum][0] = func(decode[popunum][0])
    probability = fitnessValue/np.sum(fitnessValue)
    cum_probability = np.cumsum(probability)
    return fitnessValue, cum_probability

def selectNewPopulation(decodepopu,cum_probability):
    m,n = decodepopu.shape
    newPopulation = np.zeros((m,n))
    for i in range(m):
        randomnum = np.random.random()
        for j in range(m):
            if(randomnum < cum_probability[j]):
                newPopulation[i] = decodepopu[j]
                break
    return newPopulation


def crossNewPopulation(newpopu, prob):
    m, n = newpopu.shape
    numbers = np.uint8(m*prob)
    if numbers %2 != 0:
        numbers = numbers + 1
    updatepopulation = np.zeros((m,n),dtype=np.uint8)
    index = random.sample(range(m),numbers)
    for i in range(m):
        if not index.__contains__(i):
            updatepopulation[i] = newpopu[i]
    j = 0
    while j < numbers:
        crosspoint = np.random.randint(0,n,1)
        crosspoint = crosspoint[0]
        updatepopulation[index[j]][0:crosspoint] = newpopu[index[j]][0:crosspoint]
        updatepopulation[index[j]][crosspoint:] = newpopu[index[j+1]][crosspoint:]
        updatepopulation[index[j+1]][0:crosspoint] = newpopu[index[j+1]][0:crosspoint]
        updatepopulation[index[j+1]][crosspoint:] = newpopu[index[j]][crosspoint:]
        j = j + 2
    return updatepopulation


def mutation(crosspopulation,mutaprob):
    mutationpopu = np.copy(crosspopulation)
    m,n = crosspopulation.shape
    mutationnums = np.uint8(m*n*mutaprob)
    mutationindex = random.sample(range(m*n),mutationnums)
    for geneindex in mutationindex:
        row = np.uint8(np.floor(geneindex/n))
        colume = geneindex % n
        if mutationpopu[row][colume] == 0:
            mutationpopu[row][colume] =1
        else:
            mutationpopu[row][colume] = 0
    return mutationpopu


def findMaxPopulation(population, maxevaluation, maxSize):
    maxevalue = maxevaluation.flatten()
    maxevaluelist = maxevalue.tolist()
    maxIndex = map(maxevaluelist.index, heapq.nlargest(100, maxevaluelist))
    index = list(maxIndex)
    colume = population.shape[1]
    maxPopulation = np.zeros((maxSize, colume))
    i = 0
    for ind in index:
        maxPopulation[i] = population[ind]
        i = i + 1
    return maxPopulation


def fitnessFunction():
    return lambda x: 10*np.sin(5*x) + 7*np.abs(x-5) + 10





def main():
    optimalvalue = []
    optimalvariables = []

    decisionVariables = [[0.0, 10.0]]
    delta = 0.0001
    EncodeLength = getEncodeLength(decisionVariables, delta)
    print(EncodeLength)

    initialPopuSize = 100
    population = getinitialPopulation(sum(EncodeLength), initialPopuSize)
    print(population)
    # 最大进化代数
    maxgeneration = 30

    # 交叉概率
    prob = 0.8

    # 变异概率
    mutationprob = 0.1

    # 新生成的种群数量
    maxPopuSize = 100

    for generation in range(maxgeneration):
        decode = getDecode(population, EncodeLength, decisionVariables, delta)

        evaluation, cum_proba = getFitnessValue(fitnessFunction(), decode)

        newpopulations = selectNewPopulation(population, cum_proba)

        crossPopulation = crossNewPopulation(newpopulations, prob)

        mutationpopulation = mutation(crossPopulation, mutationprob)

        totalpopalation = np.vstack((population, mutationpopulation))

        final_decode = getDecode(totalpopalation, EncodeLength, decisionVariables,delta)

        final_evaluation, final_cumprob = getFitnessValue(fitnessFunction(), final_decode)

        population = findMaxPopulation(totalpopalation, final_evaluation, maxPopuSize)

        optimalvalue.append(np.max(final_evaluation))

        index = np.where(final_evaluation == max(final_evaluation))

        optimalvariables.append(list(final_decode[index[0][0]]))
        decodeTemp = getDecode(population, EncodeLength, decisionVariables, delta)
        evaluationTemp, cum_probaTemp = getFitnessValue(fitnessFunction(), decodeTemp)
        plt.plot(decodeTemp,evaluationTemp,"*")
        x = np.arange(0.0, 10, 0.001)
        y = 10 * np.sin(5 * x) + 7 * np.abs(x - 5) + 10
        plt.plot(x, y)
        plt.show()
        time.sleep(1)


    x = [i for i in range(maxgeneration)]
    y = [optimalvalue[i] for i in range(maxgeneration)]
    plt.plot(x,y)
    plt.show()
    optimalval = np.max(optimalvalue)
    index = np.where(optimalvalue == max(optimalvalue))

    optimavar = optimalvariables[index[0][0]]
    return optimalval, optimavar


if __name__ == \'__main__\':
    optval, optvar = main()
    print("x1:",optvar[0])
    print("maxValue:",optval)

View Code

 

 

知乎上看到一个，https://zhuanlan.zhihu.com/p/35986593