【数模/启发式算法】遗传算法

简介

 遗传算法（Genetic Algorithm，GA）最早是由美国的 John holland于20世纪70年代提出，该算法是根据大自然中生物体进化规律而设计提出的。是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。

符号说明

符号	含义
$n$	种群个体个数
$a$	求解区间左端点
$b$	求解区间右端点
$length$	编码长度
$X$	编码串的二进制转十进制值
$val$	编码串的映射值
$pc$	‘基因’交叉的概率
$pm$	‘基因’变异的概率

核心思想

 遗传算法借鉴了达尔文的进化论和孟德尔的遗传学说。其本质是一种并行、高效、全局搜索的方法。其能够在搜索过程中自动获取和积累搜索空间的知识，并自适应的控制搜索过程以获得最优解。

 遗传算法使用“适者生存”的原则，在潜在的解决方案种群中逐次产生一个近似最优的方案。在遗传算法的每一代中，根据个体在问题域中的适应度和从自然学说中借鉴来的再造方法进行个体选择，产生一个新的近似解。在这个过程中导致种群的进化，得到的新个体比原个体更能适应环境，就像自然界中的改造一样。

流程图

文章使用到的测试函数

一元函数：$$y=-abs(xsin(x)cos(2x)-2xsin(3x)+3xsin(4x)),x\in [0,50]$$参考最小值为：$min(y)=-219.501$

二元函数：$$y=x_1^2+x_2^2-x_1{x}_2-10x_1-4x_2+60$$参考最小值为：$min(y)=8.0$

遗传算法基本原理

编码解码

常见的编码方式为使用二进制进行编码，并利用解码函数将二进制码映射到可行域中。对于不同的问题编码方式会作出相应调整。

当串长$length=n$，可行域为$[a,b]$时。我们将可行域划分为$2^n-1$份，其中一个二进制串s（s的10进制值$\in [0,2^{length}-1]$）所映射值的计算公式为：$$val=a+X\ast \frac{b-a}{2^{length}-1}$$其中，$X$表示s的十进制值。能够直观理解，将原区间划分为$2^{length}-1$份，每一份的值为$\frac{b-a}{2^{length}-1}$。

因此，串长决定了每个二进制串映射值的精度，串长越长则划分越细精度越高。（注意：每个串表示的精度越高，也间接影响到遗传算法计算结果的精度。通常串长取20-30较为合适。）

代码中的二进制映射函数：

# 二进制 ”转“（映射）十进制
# 公式：val = 区间begin + 串十进制值 * 区间长度 / (2^length - 1)
# 注意：为适应多元函数，返回值为一个包含var个实数值的向量
def m_decode(self, indiv):
    xpos = [0] * self.var
    # var个自变量
    for i in range(self.var):
        xpos[i] = self.lb[i] + int(indiv[i], 2) * (self.ub[i] - self.lb[i]) / (2 ** self.length - 1)

    return xpos

“基因”复制

为了得到最优解，算法模拟基因编码的复制。并控制较为优秀的基因要多复制（表现为复制概率更大），而较差一点的基因要少复制（表现为复制概率较小）。这里的优秀与否取决于对应编码的适应度。

那么，如何根据适应度来赋予复制概率，并按照概率来进行复制呢？

假设，我们在解决一个函数最大值的问题。此时适应度即为编码对应的函数值，假设fit(a) = 1，fit(b) = 2，fit(c ) = 3，fit(d) = 4。对每个编码适应度进行归一化，即可得到对应的概率P为：P(a) = 1/10, P(b) = 2/10, P(c ) = 3/10, P(d) = 4/10。

有了编码对应的概率，我们要如何依概率进行复制呢？
我们将概率进行一次累加得到[1/10, 3/10, 6/10, 10/10]，此时我们取一随机数，若随机数落入区间为$r\le 1/10$则编码1需要复制、随机数落入区间为$1/10<r\le 3/10$则编码2需要复制…

我们只需要得到一组随机数，并将其排序依次对累加概率进行比较，即可按概率进行复制。此方法相当于将区间问题转化为右端点问题，简化了原问题。

基因复制代码（预处理部分为了保证函数值越小，概率越大）：

# 基因（编码）复制
def m_copy(self):
    new_popu = [[''] * self.var for _ in range(self.n)]
    # 计算复制每个个体复制基因的概率， 概率越大复制概率越大
    ## 按照取值赋概率， 取值越小概率越大。 因此对所有取值先预处理
    ### 预处理：首先要非负，因为概率取值原理是归一化

    minval = min(self.fit)
    addval = 0
    if minval <= 0:
        # 原值要加上addval， 保证取倒数后分母不能为0。 取倒数原因：值越小需要概率越大
        addval = -minval + 1

    tem = list(map(lambda x : 1 / (x + addval), self.fit))
    s_tem = sum(tem)
    # 归一化计算概率
    p = list(map(lambda x : x / s_tem, tem))

    # 概率求累加, 为了方便让概率大的多复制
    cum_p = [0] * self.n
    cum_p[0] = p[0]
    for i in range(1, self.n):
        cum_p[i] = cum_p[i - 1] + p[i]

    # 制造一组有序的随机数， 随机数落入累加区间即需要复制
    choice = [random.random() for _ in range(self.n)]
    choice.sort()

    i = 0
    j = 0
    while j < self.n and i < self.n:
        # 只要当前随机数小于右端点就说明落入该区间，需要复制
        if choice[j] <= cum_p[i]:
            new_popu[j] = deepcopy(self.popu[i])
            j = j + 1
        else: # 随机数大于右端点，需要和之后的区间进行比较
            i = i + 1

    return new_popu

“基因”交叉

使用三段式交叉方式，即掐头，去尾，交换中间（可以参考上图）。
代码中对相邻两编码进行“基因”交叉操作，所交换部分由随机数控制。

基因交叉代码：

# 基因（编码）交叉
def m_cross(self, new_popu):
    # 相邻两个个体同种基因（同一自变量编码）进行交叉
    i = 1
    while i < self.n:
        if random.random() < self.pc:
            # 取出要进行交叉的两个个体
            indiv1 = deepcopy(new_popu[i - 1])
            indiv2 = deepcopy(new_popu[i])

            # 逐变量进行编码交叉
            # hmax 划分值：self.length - 1， 因为下边从1开始
            hmax = self.length - 1
            for j in range(self.var):
                h1 = random.randint(0, hmax)
                h2 = random.randint(0, hmax)

                if h1 > h2:
                    h1, h2 = h2, h1

                # 交叉划分构成的区间
                new_popu[i - 1][j] = indiv1[j][0 : h1] + indiv2[j][h1 : h2 + 1] + indiv1[j][h2 + 1:]
                new_popu[i][j] = indiv2[j][0 : h1] + indiv1[j][h1 : h2 + 1] + indiv2[j][h2 + 1:]

        i = i + 2

    return new_popu

“基因”变异

由于采用二进制编码，因此“基因”变异十分简单，可以描述为：随机选择位点进行取反操作即可，因为二进制串只有0和1。

“基因”变异代码：

# 基因（编码）变异
def m_mutation(self, new_popu):
    # n个个体
    for i in range(self.n): # var个自变量
        for j in range(self.var):
            if random.random() < self.pm:
                # 随机挑选2个数变异
                for k in range(2):
                    # 要变异的下标
                    idx = random.randint(0, self.length - 1)
                    if new_popu[i][j][idx] == '0':
                        new_popu[i][j] = new_popu[i][j][0:idx] + '1' + new_popu[i][j][idx + 1:]
                    else:
                        new_popu[i][j] = new_popu[i][j][0:idx] + '0' + new_popu[i][j][idx + 1:]

    return new_popu

遗传算法代码

Python版本：

import random
from copy import deepcopy
import math

def func(x):
    # 一元函数测试
    return -abs(x[0] * math.sin(x[0]) * math.cos(2 * x[0]) - 2 * x[0] * math.sin(3 * x[0]) + 3 * x[0] * math.sin(4 * x[0]))

    # # 二元函数测试
    # return x[0]**2 + x[1]**2 - x[0]*x[1] - 10 * x[0] - 4 * x[1] + 60

class GA:

    def __init__(self, func, n, var = 1, length = 20, iter = 50, lb = None, ub = None):
        """ 默认寻找最小值，以及对应自变量取值。 若需要寻找最大值，对目标函数乘-1，并对最终结果乘-1即可
        Args:
            :param func: type: 函数， des: 所要求解优化问题的目标函数
            :param n: type: int， des: 粒子群粒子的个数
            :param var: type: int， des: 函数中自变量的个数，即：几元函数
            :param length: type: int， des: 用于编码的串长了。
            :param iter: type: int， des: 迭代次数
            :param lb: type: list(double)， des: 每一种自变量的下界，注意应和自变量一一对应
            :param ub: type: list(double)， des: 每一种自变量的上界，注意应和自变量一一对应
        """
        if lb is None:
            lb = []
        if ub is None:
            ub = []

        # 目标函数
        self.func = func
        # 初始化种群个数
        self.n = n
        # 初始化编码串长
        self.length = length
        # 初始化变量种类数
        self.var = var
        # 初始化迭代次数
        self.iter = iter
        # 初始化自变量范围 len(lb) = len(ub) = var
        self.lb = lb
        self.ub = ub

        # 交叉概率和变异概率 0.6 - 0.1
        self.pc = 0.6
        self.pm = 0.1

        # 初始化适应度
        self.fit = [0] * n

        # 产生初始群体 每一行为一个种群个体的var个字符串, 因为每一个个体需要var个变量来描述
        ## var等于1时即为一元函数
        self.popu = [[''] * var for _ in range(n)]
        # 使用随机数初始个体
        for i in range(n): # n个个体
            for j in range(var): # 每个个体的自变量个数
                self.popu[i][j] = str(bin(random.randint(0, int(2 ** length))))[2:]
                if len(self.popu[i][j]) < length:
                    self.popu[i][j] = '0' * (length - len(self.popu[i][j])) + self.popu[i][j]

        # 最优解对应的对应自变量取值
        self.x = self.m_decode(self.popu[0])
        self.fit[0] = self.func(self.x)

        # 用初代值计算一组解
        for i in range(1, n):
            now_x = self.m_decode(self.popu[i])
            self.fit[i] = self.func(now_x)
            if self.fit[i] < self.func(self.x): # 以求解最小值的方式更新当前最优解
                self.x = deepcopy(now_x)



    # 二进制 ”转“（映射）十进制
    # 公式：val = 区间begin + 串十进制值 * 区间长度 / (2^length - 1)
    # 注意：为适应多元函数，返回值为一个包含var个实数值的向量
    def m_decode(self, indiv):
        xpos = [0] * self.var
        # var个自变量
        for i in range(self.var):
            xpos[i] = self.lb[i] + int(indiv[i], 2) * (self.ub[i] - self.lb[i]) / (2 ** self.length - 1)

        return xpos


    # 基因（编码）复制
    def m_copy(self):
        new_popu = [[''] * self.var for _ in range(self.n)]
        # 计算复制每个个体复制基因的概率， 概率越大复制概率越大
        ## 按照取值赋概率， 取值越小概率越大。 因此对所有取值先预处理
        ### 预处理：首先要非负，因为概率取值原理是归一化

        minval = min(self.fit)
        addval = 0
        if minval <= 0:
            # 原值要加上addval， 保证取倒数后分母不能为0。 取倒数原因：值越小需要概率越大
            addval = -minval + 1

        tem = list(map(lambda x : 1 / (x + addval), self.fit))
        s_tem = sum(tem)
        # 归一化计算概率
        p = list(map(lambda x : x / s_tem, tem))

        # 概率求累加, 为了方便让概率大的多复制
        cum_p = [0] * self.n
        cum_p[0] = p[0]
        for i in range(1, self.n):
            cum_p[i] = cum_p[i - 1] + p[i]

        # 制造一组有序的随机数， 随机数落入累加区间即需要复制
        choice = [random.random() for _ in range(self.n)]
        choice.sort()

        i = 0
        j = 0
        while j < self.n and i < self.n:
            # 只要当前随机数小于右端点就说明落入该区间，需要复制
            if choice[j] <= cum_p[i]:
                new_popu[j] = deepcopy(self.popu[i])
                j = j + 1
            else: # 随机数大于右端点，需要和之后的区间进行比较
                i = i + 1

        return new_popu


    # 基因（编码）交叉
    def m_cross(self, new_popu):
        # 相邻两个个体同种基因（同一自变量编码）进行交叉
        i = 1
        while i < self.n:
            if random.random() < self.pc:
                # 取出要进行交叉的两个个体
                indiv1 = deepcopy(new_popu[i - 1])
                indiv2 = deepcopy(new_popu[i])

                # 逐变量进行编码交叉
                # hmax 划分值：self.length - 1， 因为下边从1开始
                hmax = self.length - 1
                for j in range(self.var):
                    h1 = random.randint(0, hmax)
                    h2 = random.randint(0, hmax)

                    if h1 > h2:
                        h1, h2 = h2, h1

                    # 交叉划分构成的区间
                    new_popu[i - 1][j] = indiv1[j][0 : h1] + indiv2[j][h1 : h2 + 1] + indiv1[j][h2 + 1:]
                    new_popu[i][j] = indiv2[j][0 : h1] + indiv1[j][h1 : h2 + 1] + indiv2[j][h2 + 1:]

            i = i + 2

        return new_popu


    # 基因（编码）变异
    def m_mutation(self, new_popu):
        # n个个体
        for i in range(self.n): # var个自变量
            for j in range(self.var):
                if random.random() < self.pm:
                    # 随机挑选2个数变异
                    for k in range(2):
                        # 要变异的下标
                        idx = random.randint(0, self.length - 1)
                        if new_popu[i][j][idx] == '0':
                            new_popu[i][j] = new_popu[i][j][0:idx] + '1' + new_popu[i][j][idx + 1:]
                        else:
                            new_popu[i][j] = new_popu[i][j][0:idx] + '0' + new_popu[i][j][idx + 1:]

        return new_popu


    def run(self):
        for k in range(1, self.iter + 1):

            # 控制交叉概率线性递减 [0.6 - 0.2]
            self.pc = 0.6 - 0.4 * k / self.iter
            # 控制变异概率线性递减 [0.3 - 0.1]
            self.pm = 0.3 - 0.2 * k / self.iter

            # "基因"复制
            new_popu = self.m_copy()
            # "基因"交叉
            new_popu = self.m_cross(new_popu)
            # "基因"变异
            new_popu = self.m_mutation(new_popu)

            for i in range(self.n):

                new_fit = self.func(self.m_decode(new_popu[i]))
                if new_fit < self.fit[i]:
                    # 概率保留优秀基因，因为较差的基因也有可能更新出最优解
                    if random.random() < 0.5:
                        self.fit[i] = new_fit
                        self.popu[i] = deepcopy(new_popu[i])

                    self.x = self.m_decode(new_popu[i])

        print(f'最优解为：{self.func(self.x)}')
        print(f'最优解对应自变量取值为：{self.x}')



"""一元函数测试"""
ga = GA(func, 50, 1, 20, 100, [0], [50])


"""二元函数测试"""
# ga = GA(func, 50, 2, 20, 100, [-15, -15], [15, 15])

ga.run()

以一元函数测试为例，对应输出为（结果较为不错）：

最优解为：-219.4965723300476
最优解对应自变量取值为：[47.55863910545264]