数学建模之【遗传算法】

遗传算法

一、介绍

​ 遗传算法是用于解决最优化问题的一种搜索算法。遗传算法借用了生物学里达尔文的进化理论:“适者生存,优胜劣汰”,将该理论以算法的形式表现出来就是遗传算法的过程。

二、基本原理

1)程序流程

​ 遗传算法是通过大量备选解的变换、迭代和变异,在解空间中并行动态地进行全局搜索的最优化方法,是模拟生物基因遗传的做法,算法的基本原理通过编码组成的群体组成初始群体后,遗传操作的任务就算对群体的个体按照他们对环境适应度(适应度)评估,施加了一定的操作之后从而实现优胜劣汰的进化过程,算法的基本程序框图:

enter

2)物竞天择

​ 我们可以假设有一个种群,这个种群存在的目标是帮我寻找一个函数F(x)的极值(可能是max or min),首先我们随机生成一组种群,假设我对变量x的取值范围为【0,8】,这是我们的变量域,因为我们要玩一些不一样的操作,先跟着我,后面你会恍然大悟!我们会有一个二进制域【0,2^8】,我们的二进制域对应的就是我们的变量域,我们将数字8看成8位的二进制数,可以得到2的八次方为:256,但因为我们在二进制当中我们是从0开始的,0-255,有256个数,我们的八位二进制数最多只能取到255,取不到256。

"二进制" 0 0 0 0   0 0 0 0
"二进制数 x-> " 1 1 1 1  1 1 1 1
"x转化为十进制" 1*2^0 + 1*2^1 + 1*2^2 + 1*2^3 + 1*2^4 + 1*2^5 + 1*2^6 + 1*2^7 = 255

​ 因为我们实数x的范围只有【0,8】,我们将二进制转化为十进制之后,还有求解二进制域到变量域值,二进制域变量域之间的映射关系为:

"变量域x_"  x_ = x*8/2^8

​ 如此,我们将我们的8看作是染色体的长度,这样我们就可以做一些基因重组和基因变异的操作了!!!,我们先来一组基因重组图:

在这里插入图片描述

​ 是不是很熟悉,这个不就是高中生物的杂交实验吗?如果我们把DDdd看成我们染色体的01呢?我们再形象一点

# 定义我们的染色体
x1 = [0 1 0 0 1 0 0 1]
x2 = [1 1 0 1 0 0 0 1]
#OK 我们现在有两个种群了(抽象一点)
# 第一个种群 x1 的染色体为 0 1 0 0  1 0 0 1
# 第一个种群 x2 的染色体为 1 1 0 1  0 0 0 1

#我们尝试做一下杂交实验?学一学孟德尔的八年抗战
"第一代:" 0 1 0 0    1 0 0 1      *      1 1 0 1    0 0 0 1
                   |                              |
                   v                              v
        
"配子:" [0 1 0 0]    [1 0 0 1 ]        [1 1 0 1 ]   [0 0 0 1]  
														 #用x1的第一个配子和x2的第二个配子进行基因重组
														 #用x1的第二个配子和x2的第一个配子进行基因重组

"F1代:" [0 1 0 0][0 0 0 1]              [1 1 0 1 ][1 0 0 1] 

new_x1 = [0 1 0 0 0 0 0 1]
new_x2 = [1 1 0 1 1 0 0 1]

# 不知道你有没有恍然大悟

​ 有了二进制的操作,我们做变异的时候也是同样的道理,我们只需要把0变异为1,1变异为0,就能完成我们的变异算法。

3)算法步骤

​ 求解函数:f(x) = 9sin(5x) + 7cos(4*x)的最大值,x的取值范围为【0,8】

3.1初始化种群

​ 我们先初始化我们的种群,先初始化一些随机解,因为我们要在搜索空间内寻找我们的最优解,可参考粒子群算法的思想。

def initpop(popsize, chromlength):
    """
    Parameters
    ----------
    popsize :  TYPE
        种群的数目.
    chromlength : TYPE
        表示染色体的长度(二值数的长度), 长度的大小取决于变量的二进制编码的长度     

    Returns 
    -------
    pop : TYPE
        随机种群
    """
    pop=np.round(np.random.rand(popsize,chromlength))
    return pop

3.2二进制编码转化为十进制数

​ 设计我们的算法,将二进制的矩阵转化为十进制数,为了得到我们的二进制域。

def decodebinary(pop):
    """
    二进制数转十进制数函数
    Parameters
    ----------
    pop : TYPE
          初始化种群.

    Returns
    -------
    pop2 种群的染色体二进制数转十进制数

    """
    px,py = pop.shape
    pop1 = np.zeros((px,py))
    for i in range(py):
        pop1[:,i] = (2**(py-i-1))*pop[:,i]
    pop2=np.sum(pop1, 1)#对pop1向量的每行求和 标识位0代表每列求和,标识位1代表每行求和
    
    return pop2

def encode(bestpop):
    """
    解码x的值

    Parameters
    ----------
    bestpop : TYPE
        最好的个体.

    Returns
    -------
    x : TYPE
        最优X的取值.

    """
    
    temp = decodebinary(pop)
    x=temp*15/32767
    return x

​ 我们的pop2是一个转化为十进制数之后的种群,可以理解为种群的表现型,我们的0和1的操作反应的是种群的基因型,因为我们的变量不可能只有一种,扩展我们的思维,我们先设置的染色体长度,因为我们的变量只有一个,它的值域在【0,8】,如果我们有两个变量,不同变量对应不同值域,我们的染色体的长度是否也要切片计算?仔细思考。

def decodechrom(pop,spoint,length):
    """
    将二进制编码转化为十进制数
    Parameters
    ----------
    pop : TYPE
        DESCRIPTION.
    spoint : TYPE
        染色体的起始位.
    length : TYPE
        DESCRIPTION.
    Returns
    -------
    pop2 : TYPE
        DESCRIPTION.

    """
    #对于多个变量而言,如有两个变量,采用20为表示,每个变量10位,则第一个变量从1开始,另一个变量从11开始。本例为一个变量
    #这句话的意思就是加入目标函数需要两个变量,则我可将染色体的数量拆为两个,然后遗传迭代
    #值得注意的是,我的染色体的长度也要跟着变量的数量改变,呈倍数关系
    pop1=pop[:,spoint:spoint+length]
    pop2=decodebinary(pop1)
    
    return pop2

​ 这里的spoint代表的是我们的一个终止位,可以对我们的变量进行一个切片,自行理解。

3.3实现目标函数的计算

​ 实现目标函数的计算,就是把x的值带进去,很简单。

def calobjvalue(pop):
    """
    实现目标函数的计算
    Parameters
    ----------
    pop : TYPE
        种群.

    Returns
    -------
    objvalue : TYPE
        返回目标函数值.

    """
    temp = decodechrom(pop, 0, 15)
    x=temp*15/32767
    objvalue = 9*np.sin(5*x)+7*np.cos(4*x)
    objvalue = np.reshape(objvalue,(-1,1))#以行的形式输出
    return objvalue

3.4剔除0以外的值

​ 计算个体的适应值,在calobjvalue已经计算好,需要将小于0的个体删除,方便后续的概率计算。

def calfitvalue(objvalue):
    """
    计算个体的适应值,在calobjvalue已经计算好,需要将小于0的个体删除,方便后续的概率计算

    Parameters
    ----------
    objvalue : TYPE
        目标函数值.
    Returns
    -------
    fitvalue : TYPE
        个体适应值.

    """
    global Cmin
    Cmin = 0
    fitvalue = np.zeros((objvalue.shape[0],objvalue.shape[1]))
    
    px, py = objvalue.shape
    for i in range(px):
        if objvalue[i,0] + Cmin > 0:
            temp = objvalue[i,0] + Cmin
        else:
            temp = 0
        
        fitvalue[i,0]=temp
    return fitvalue

3.5选择算法

​ 选择算法我们采用轮盘赌,决定哪些个体可以进入下一代,用轮盘赌选择复制

def selection(pop, fitvalue):
    """
    选择函数 选择复制,决定哪些个体可以进入下一代
    采用轮盘赌选择

    Parameters
    ----------
    pop : TYPE
        种群的个体.
    fitvalue : TYPE
        个体适应值.

    Returns
    -------
    newpop : dict
        新的种群.

    """
    totalfit = np.sum(fitvalue)
    fitvalue_pro = fitvalue / (totalfit + 0.000001)
    fitvalue_pro_cumnsum = np.cumsum(fitvalue_pro)
    fitvalue_pro_cumnsum = np.reshape(fitvalue_pro_cumnsum,(-1,1))#以行的形式
    
    px, py = pop.shape 
    #轮盘随机概率 从小到大排序
    ms = np.sort(np.random.rand(px,1),0)
    
    fitin = 1 - 1   #pop种群 第几代个体 因为python的下标从0开始。所以是第一代的索引是0 故用1-1
    newin = 1 - 1  #pop种群 第几代个体
    newpop_dict = {}
    #我愿称为[适者生存,优胜劣汰] while循环
    while newin <= px-1:
        if ms[newin, 0] < fitvalue_pro_cumnsum[fitin, 0]:
            newpop_dict[newin] = pop[fitin,:]
            newin += 1
        else:
            fitin += 1
    
    newpop = np.zeros((newin,py))
    for i in range(newin):
        newpop[i,:] = newpop_dict[i]
    
    return newpop

3.6交叉算法(基因重组)

​ 物竞天择是遗传学的核心,相应的,选择和变异也是遗传算法的核心,我们对染色体进行一个基因重组的操作,具体代码如下:

def crossover(pop, pc):
    """
    交叉算法 实现基因重组
    Parameters
    ----------
    pop : TYPE
        种群.
    pc : TYPE
        交叉概率.
    Returns
    -------
    newpop : TYPE
        DESCRIPTION.

    """
    px, py = pop.shape 
    newpop = np.zeros((px, py))
    # seletion_litst = [i for i in range(px)]
    # sl = seletion_litst[0:px:2]
    for i in range(px-1):
        #是否能够进行基因重组
        if pc > np.random.rand(1)[0]:
            cpoint = int(np.round(np.random.rand(1)[0] * py) - 1)
            if cpoint == 0:
                cpoint = 1
                
            newpop[i, :][0:cpoint] = pop[i, 0:cpoint]
            newpop[i, :][cpoint+1:py] = pop[i+1, cpoint+1:py]
            
            newpop[i+1, :][0:cpoint] = pop[i+1, 0:cpoint]
            newpop[i+1, :][cpoint+1:py] = pop[i, cpoint+1:py]
            
        else:
            newpop[i,:] = pop[i,:]
            newpop[i+1,:] = pop[i+1,:]
    
    return newpop

3.7变异算法(基因突变)

​ 0变1,1变0 你上你也行

def mutation(pop, pm):
    """
    变异算法 实现基因突变

    Parameters
    ----------
    pop : TYPE
        种群.
    pm : TYPE
        变异概率.

    Returns
    -------
    newpop : TYPE
        新的变异种群.

    """
    px, py = pop.shape 
    newpop = np.zeros((px, py))
    
    for i in range(px):
        if pm > np.random.rand(1)[0]:
            mpoint = int(np.round(np.random.rand(1)[0] * py) - 1)
            if mpoint == 0:
                mpoint = 1
            newpop[i,:] = pop[i,:]
            if newpop[i, mpoint] == 0:
                newpop[i, mpoint] = 1

        newpop[i,:] = pop[i,:]
    return newpop

3.7最优计算(基因突变)

​ 简单的排序取最大(不要忘记我们粒子群算法里面的全局极值和个体极值)

def best(pop, fitvalue):
    """
    最优的个体及其适应值
    Parameters
    ----------
    pop : TYPE
        种群.
    fitvalue : TYPE
        适应值.
    Returns
    -------
    bestindividual : 最大适应值
        DESCRIPTION.
    bestfit : TYPE
        最大适应值的个体.

    """
    px, py = pop.shape 
    bestindividual = pop[0,:]
    bestfit = fitvalue[0]
    
    for i in range(1,px):
        if fitvalue[i] > bestfit:
            bestindividual = pop[i,:]
            bestfit = fitvalue[i]
    
    return bestindividual,bestfit

3.完整代码

​ 不要忘记导入我们的numpy 和matplotlib库,将上面的函数和导入的库,放在主函数的上面,然后运行代码,不过我这里的取值是【0,15】,上面说【0,8】便于读者理解。

import numpy as np
import matplotlib.pyplot as plt

# 遗传算法
def initpop(popsize, chromlength):
    """
    Parameters
    ----------
    popsize :  TYPE
        种群的数目.
    chromlength : TYPE
        表示染色体的长度(二值数的长度), 长度的大小取决于变量的二进制编码的长度     

    Returns 
    -------
    pop : TYPE
        随机种群

    """
    pop=np.round(np.random.rand(popsize,chromlength))
    return pop
    
def decodebinary(pop):
    """
    二进制数转十进制数函数
    Parameters
    ----------
    pop : TYPE
          初始化种群.

    Returns
    -------
    pop2 种群的染色体二进制数转十进制数

    """
    px,py = pop.shape
    pop1 = np.zeros((px,py))
    for i in range(py):
        pop1[:,i] = (2**(py-i-1))*pop[:,i]
    pop2=np.sum(pop1, 1)#对pop1向量的每行求和 标识位0代表每列求和,标识位1代表每行求和
    
    return pop2

def decodechrom(pop,spoint,length):
    """
    将二进制编码转化为十进制数
    Parameters
    ----------
    pop : TYPE
        DESCRIPTION.
    spoint : TYPE
        染色体的起始位.
    length : TYPE
        DESCRIPTION.

    Returns
    -------
    pop2 : TYPE
        DESCRIPTION.

    """
    #对于多个变量而言,如有两个变量,采用20为表示,每个变量10位,则第一个变量从1开始,另一个变量从11开始。本例为一个变量
    #这句话的意思就是加入目标函数需要两个变量,则我可将染色体的数量拆为两个,然后遗传迭代
    #值得注意的是,我的染色体的长度也要跟着变量的数量改变,呈倍数关系
    pop1=pop[:,spoint:spoint+length]
    pop2=decodebinary(pop1)
    
    return pop2

def encode(bestpop):
    """
    解码x的值

    Parameters
    ----------
    bestpop : TYPE
        最好的个体.

    Returns
    -------
    x : TYPE
        最优X的取值.

    """
    
    temp = decodebinary(pop)
    x=temp*15/32767
    return x

def calobjvalue(pop):
    """
    实现目标函数的计算
    Parameters
    ----------
    pop : TYPE
        种群.

    Returns
    -------
    objvalue : TYPE
        返回目标函数值.

    """
    temp = decodechrom(pop, 0, 15)
    x=temp*15/32767
    objvalue = 9*np.sin(5*x)+7*np.cos(4*x)
    objvalue = np.reshape(objvalue,(-1,1))#以行的形式输出
    return objvalue

def calfitvalue(objvalue):
    """
    计算个体的适应值,在calobjvalue已经计算好,需要将小于0的个体删除,方便后续的概率计算

    Parameters
    ----------
    objvalue : TYPE
        目标函数值.

    Returns
    -------
    fitvalue : TYPE
        个体适应值.

    """
    global Cmin
    Cmin = 0
    fitvalue = np.zeros((objvalue.shape[0],objvalue.shape[1]))
    px, py = objvalue.shape
    for i in range(px):
        if objvalue[i,0] + Cmin > 0:
            temp = objvalue[i,0] + Cmin
        else:
            temp = 0
        
        fitvalue[i,0]=temp
    return fitvalue

def selection(pop, fitvalue):
    """
    选择函数 选择复制,决定哪些个体可以进入下一代
    采用轮盘赌选择

    Parameters
    ----------
    pop : TYPE
        种群的个体.
    fitvalue : TYPE
        个体适应值.

    Returns
    -------
    newpop : dict
        新的种群.

    """
    totalfit = np.sum(fitvalue)
    fitvalue_pro = fitvalue / (totalfit + 0.000001)
    fitvalue_pro_cumnsum = np.cumsum(fitvalue_pro)
    fitvalue_pro_cumnsum = np.reshape(fitvalue_pro_cumnsum,(-1,1))#以行的形式
    
    px, py = pop.shape 
    #轮盘随机概率 从小到大排序
    ms = np.sort(np.random.rand(px,1),0)
    
    fitin = 1 - 1   #pop种群 第几代个体 因为python的下标从0开始。所以是第一代的索引是0 故用1-1
    newin = 1 - 1  #pop种群 第几代个体
    newpop_dict = {}
    #我愿称为[适者生存,优胜劣汰] while循环
    while newin <= px-1:
        if ms[newin, 0] < fitvalue_pro_cumnsum[fitin, 0]:
            newpop_dict[newin] = pop[fitin,:]
            newin += 1
        else:
            fitin += 1
    
    newpop = np.zeros((newin,py))
    for i in range(newin):
        newpop[i,:] = newpop_dict[i]
    
    return newpop
    


def crossover(pop, pc):
    """
    交叉算法 实现基因重组
    Parameters
    ----------
    pop : TYPE
        种群.
    pc : TYPE
        交叉概率.

    Returns
    -------
    newpop : TYPE
        DESCRIPTION.

    """
    px, py = pop.shape 
    newpop = np.zeros((px, py))
    # seletion_litst = [i for i in range(px)]
    # sl = seletion_litst[0:px:2]
    for i in range(px-1):
        #是否能够进行基因重组
        if pc > np.random.rand(1)[0]:
            cpoint = int(np.round(np.random.rand(1)[0] * py) - 1)
            if cpoint == 0:
                cpoint = 1
                
            newpop[i, :][0:cpoint] = pop[i, 0:cpoint]
            newpop[i, :][cpoint+1:py] = pop[i+1, cpoint+1:py]
            
            newpop[i+1, :][0:cpoint] = pop[i+1, 0:cpoint]
            newpop[i+1, :][cpoint+1:py] = pop[i, cpoint+1:py]
            
        else:
            newpop[i,:] = pop[i,:]
            newpop[i+1,:] = pop[i+1,:]
    
    return newpop

def mutation(pop, pm):
    """
    变异算法 实现基因突变

    Parameters
    ----------
    pop : TYPE
        种群.
    pm : TYPE
        变异概率.

    Returns
    -------
    newpop : TYPE
        新的变异种群.

    """
    px, py = pop.shape 
    newpop = np.zeros((px, py))
    
    for i in range(px):
        if pm > np.random.rand(1)[0]:
            mpoint = int(np.round(np.random.rand(1)[0] * py) - 1)
            if mpoint == 0:
                mpoint = 1
            newpop[i,:] = pop[i,:]
            if newpop[i, mpoint] == 0:
                newpop[i, mpoint] = 1

        newpop[i,:] = pop[i,:]
    return newpop

def best(pop, fitvalue):
    """
    最优的个体及其适应值
    Parameters
    ----------
    pop : TYPE
        种群.
    fitvalue : TYPE
        适应值.

    Returns
    -------
    bestindividual : 最大适应值
        DESCRIPTION.
    bestfit : TYPE
        最大适应值的个体.

    """
    px, py = pop.shape 
    bestindividual = pop[0,:]
    bestfit = fitvalue[0]
    
    for i in range(1,px):
        if fitvalue[i] > bestfit:
            bestindividual = pop[i,:]
            bestfit = fitvalue[i]
    
    return bestindividual,bestfit



    
if __name__ == "__main__":
    popsize=30; #群体大小
    chromlength=15; #字符串长度(染色体的长度)
    pc=0.7; #交叉概率 
    pm=0.005 #变异概率
    
    #初始化种群
    pop=initpop(popsize, chromlength)
    poptest=pop
    #开始迭代 
    epoch = 200
    
    x_ = []
    y_ = []
    
    for i in range(epoch):
        objvalue = calobjvalue(pop)     #计算目标函数值
        fitvalue = calfitvalue(objvalue)    #计算群体中每个个体的适应度
        
        newpop = selection(pop, fitvalue)
        newpop1 = crossover(newpop, pc)
        newpop2 = mutation(newpop1, pm)
        
        objvalue = calobjvalue(newpop2)     #计算目标函数值
        fitvalue = calfitvalue(objvalue)    #计算群体中每个个体的适应度
        
        bestindividual,bestfit = best(newpop2, fitvalue) #求出群体中适应值最大的个体及其适应值
        x_.append(bestfit)
        pop = newpop2
        
        plt.plot(x_, lw=3, label='funcotio_max_value_x')
        plt.show()
        
    reshape_best_infrivedual = np.reshape(bestindividual,(1,-1))#以行的形式
    best_x_value = encode(reshape_best_infrivedual)
    plt.plot(x_, lw=3, label='funcotio_max_value_x')
    plt.legend()
    plt.show()
    
    print("最好的个体是 ",bestindividual)
    print("函数:f(x) = 9sin(5x)+7cos(4x)的极值为 ",bestfit[0]) #最后拟合之后随便取一个值即可
    print("x的取值应为: ",best_x_value[0]) #最后拟合之后随便取一个值即可

4)运行结果

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-gbMtHJ8N-1650549372734)(genetic.assets/image-20220421215451632.png)]

最好的个体是 [1. 0. 0. 0. 0. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0.]

函数:f(x) = 9sin(5x)+7cos(4x)的极值为 15.999199583116948

x的取值应为: 7.851802117984557

评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值