遗传算法（GA）的python实现

遗传算法（Genetic Algorithm，简称GA），属于启发式算法的一种，利用该算法我们能够在一定的变量空间中搜索我们需要的最优解。

遗传算法核心就是模拟了自然界中的“物竞天择、自然选择”的生存法则，去解决我们日常生活中的问题。遗传算法中的主要步骤和我们所学习的生物学知识一样，包括了自然选择、染色体交叉、碱基对变异。那么我们将会以这样的方式进行描述。

一、种群（Pop）、DNA和碱基对

（1）种群Pop

首先，我们明白了遗传算法就是模拟了自然界的适者生存的法则。那么我们就得明白，首先肯定有一个概念叫做种群，或者叫做物种，我们将这个种群起一个名字叫做Pop。那么这个种群内一共有若干个个体。如果打个比方的话，就像是人类种群一样，一个人类种群会有几十亿个个体。

并且每一个人都是独特的独一无二的存在的，那么是什么决定了个体的独特性呢？

（2）DNA和染色体

没错就是染色体，每个人拥有的46条染色体决定了每个人的独特性。假设我们需要解决一个函数，那么在我们想要解决的问题中，什么能够代表一个输入呢？对了！那就是这四个自变量。换句话来说，这就是这个问题的特征变量或者是决策变量。那么在这个问题中，每个“个体”的DNA就是。

（3）碱基对

这时，我们知道了什么是遗传算法中的DNA，那么又怎么样去表示碱基对呢？很简单，在计算机中二进制能够完美的代表一个数字，那么我们将代表DNA数值的这样一串二进制数组描述为碱基对，那么这样一个种群就建立啦！

但是现在还存在一个问题，我们的每一个变量都有不同的定义域，那么怎么去映射成二进制呢？其实我们可以自定义碱基对的长度，例如碱基对为10位，那么10位二进制可以表示的数值范围是[0, 2^10)即[0, 1023]，我们将的定义域映射到二进制范围中，最终转换为二进制数组。这就是DNA编码。并且我们不仅要能够将DNA进行编码，也要能根据二进制进行DNA解码。编码过程，如下图所示：

例如，我们要表示10这个数字，我们首先使用

我们向下取整，我们就能够用546的二进制表示10这个数字。

那么546.821又怎么能够转化为10呢？

这样我们又能够将其转化回去，那么我们的具体的代码是如何实现的呢？其具体代码如下：

#将编码后的DNA翻译回来（解码）
    def translateDNA(self, index=None):
        W_vector = np.array([2**i for i in range(self.dna_size)]).reshape((self.dna_size, 1)) # shape = (dna_size, 1)
        if index != None:
            binary_vector = self.pop[index].dot(W_vector).reshape(-1) # shape = (pop_size, nums)
            for i in range(self.nums):
                binary_vector[i] = binary_vector[i] / ((2**self.dna_size)/self.var_len[0])
                binary_vector[i] = binary_vector[i] + self.bound[0][0]
            return binary_vector
            # print(binary_vector)
            # exit()
        else:
            # W_vector = [[0],[2],[4],[8],[16],[32] ..., [2^dna_size]]
            # (100, 4, 32) · (32, 1) = (100, 4, 1) # 先乘后加，然后就将dna转换成了一个十进制的数字
            binary_vector = self.pop.dot(W_vector).reshape(self.pop.shape[0:2]) # shape = (pop_size, nums)
            for i in range(self.pop_size):
                for j in range(self.nums):
                    binary_vector[i, j] = binary_vector[i, j] / ((2**self.dna_size)/self.var_len[j])
                    binary_vector[i, j] = binary_vector[i, j] + self.bound[j][0]
            return binary_vector

二、自然选择

我们如何去筛选掉“不适应环境”的个体呢？然而我们在代码中需要有一个“评判标准”去看哪个个体需要被淘汰。

在遗传算法中，这个评判标准就是适应度函数。我们给每一个个体计算出一个适应度值，适应度值越大就说明该个体更应该被选择。我们将适应度值转换为概率，按照这个概率进行采样，这就模拟出了自然界中的自然选择啦！

那么在具体问题中，适应度函数就是我们要求界的目标函数啦。其具体代码如下：

#得到适应度
    def get_fitness(self):
        x = self.translateDNA()
        res = []
        for xx in x:
            temp = self.F(xx)[0]
            temp = temp if temp > 0 else 0
            res.append(temp)
        self.fitness = res
        return res
    
    #得到个体适应度
    def get_fitness_by_index(self, index=None):
        x = self.translateDNA(index)
        res = self.F(x)[0]
        return res

    #自然选择,选择较好的个体
    def select(self):
        fitness = self.get_fitness() # 获取适应度值
        # print(fitness)
        # exit()
        # np.random.choice
        # 第一个参数 a：随机数的来源，可以是数组、元组，但一定是一维的。
        # 第二个参数 size：随机数的个数，默认值为1
        # replace：元素是否可以重复采样
        # p：一个数组，大小与a相同，代表着每个元素选取的概率。默认每个概率一样
        # 假设种群大小为n，此代码为在0-n中选取n个元素， 可重复取值，根据适应度的概率进行随机抽样。
        choice = np.random.choice(np.arange(self.pop_size), size=self.pop_size, replace=True, p=fitness/np.sum(fitness)) # 筛选出优质种群
        self.pop = self.pop[choice] # 进行自然选择

三、染色体交叉

什么是染色体交叉，就是在种群中，有两个个体发生了交配，那么他们的DNA中碱基对就会按照一定概率发生交换，这个过程我们就将其称之为染色体交叉。那么其具体代码实现如下：

这个代码中有一个优化，就是说，如果我新交配产生的子代还不如原来的个体，那么我们就将其返回原来的个体。这样我们能够减少模型的震荡，能够加速收敛。

#染色体交叉，选择个体进行
    def crossover(self):
        # arr_max = heapq.nlargest(int(self.pop_size*0.1),self.fitness)#获取前五大的值并排序
        # best_ten_percent_people_index = list(map(self.fitness.index,arr_max))#获取前五大的值下标 
        for i, people in enumerate(self.pop):
            # if i in best_ten_percent_people_index: 
            #     continue 
            temp = people.copy()
            temp_fit = self.get_fitness_by_index(i)
            if np.random.rand() < self.cross_rate: # 染色体交叉概率
                mate = np.random.randint(0, self.pop.shape[0], size=1) # 选择一个个体进行染色体交叉
                cross_points = np.random.randint(0, 2, size=(self.nums ,self.dna_size)).astype(bool) # 随机生成要交叉的位，n条染色体，dna_size个碱基对
                people[cross_points] = self.pop[mate, cross_points] # 布尔索引然后进行染色体交叉
            cur_fit = self.get_fitness_by_index(i)
            if cur_fit < temp_fit: self.pop[i] = temp

四、碱基对变异

我们的DNA在个体的生存过程中，我们会受到环境的影响，那么我们的碱基对偶尔也换发生变异，这就称之为基因变异。

在计算机中我们的碱基对并没有那么复杂，没有A、T、C、G，我们的碱基对就是二进制数组，我们仅需要将0变成1，将1变成0，那就完成了基因变异的过程了。其具体代码如下：

这里也和上面的交叉一样，做了一个同样的小技巧。

#基因变异,dna中的某个碱基发生变异
    def mutate(self):
        # arr_max = heapq.nlargest(int(self.pop_size*0.1),self.fitness)#获取前五大的值并排序
        # best_ten_percent_people_index = list(map(self.fitness.index,arr_max))#获取前五大的值下标 
        for i, people in enumerate(self.pop):
            # if i in best_ten_percent_people_index: 
            #     continue 
            temp = people.copy()
            temp_fit = self.get_fitness_by_index(i)
            for dna in people: # 多条DNA
                for point in range(self.dna_size): # 遍历碱基对
                    if np.random.rand() < self.mutation: # 变异概率
                        dna[point] = 1 if dna[point] == 0 else 0 # 基因变异
            cur_fit = self.get_fitness_by_index(i)
            if cur_fit < temp_fit: self.pop[i] = temp

五、种群进化

这就相对简单，就是自然选择、染色体交叉和变异全流程。其具体代码如下：

#进化，单次进化包括自然选择，染色体交叉和基因变异
    def evolution(self):
        for e in range(self.epochs):
            self.gen = e
            self.select()
            self.crossover()
            self.mutate()
            self.log()

六、python代码实现（以）为例）

每行代码都有详细的注释！！！
针对不同的目标函数需要进行修改。

实验结果

"""
测试函数 y = -(x-1.2)**2 + 52.0
"""
# 先导入所需要的库
#可以用来计算算法运行时间
import datetime
#数据运算，操作库
import numpy as np
import heapq

# 遗传算法类
class GA():
    def __init__(self, nums, bound, func, dna_size=None, cross_rate=0.8, mutation=0.003, pop_size=300, epochs=500):
        # 构造函数
        # 参数列表：
        # - nums：需要编码的数，也就是需要nums个染色体
        # - bound：每一个编码数的边界，输入为数组，n * 2  例如：[(min, nax), (min, max), (min, max),...]
        # - func: 需要优化的函数(实际上就是适应度函数，即get_fitness函数)
        # - DNA_SIZE：二进制位数，表示一个编码
        # - cross_rate：交叉概率
        # - mutation：变异概率
        # - pop_size：种群大小
        # - epochs：迭代次数
        self.nums = nums
        self.bound = np.array(bound)
        self.F = func 
        self.dna_size = dna_size
        self.cross_rate = cross_rate 
        self.mutation = mutation 
        self.pop_size = pop_size
        self.epochs = epochs 
        self.startTime = datetime.datetime.now() # 运行开始时间
        self.gen = 0 # 当前迭代次数
        self.fitness = [0]

        if self.bound.shape[0] != self.nums or self.bound.shape[1] != 2: 
            raise Exception(f'Please check the length of bound is the same with nums,  bound:{self.bound.shape[0]} != nums:{self.nums}')

        for i, (min_bound, max_bound) in enumerate(self.bound):
            if max_bound < min_bound:
                raise Exception(f'Sorry, bound[{i}]: (min:{min_bound}, max:{max_bound})is not allowed!')

        min_nums, max_nums = np.array(list(zip(*self.bound))) # 返回值是两个数组，将最大值和最小值分离
        self.var_len = var_len = max_nums-min_nums # 最大值和最小值之间的差距
        bits = np.ceil(np.log2(var_len+1)) # 用多少位数可以表示所有的编码变量

        if self.dna_size == None: self.dna_size = int(np.max(bits))

        # 初始化种群
        self.init_pop()
        self.copy_pop = self.pop.copy() # 种群备份，用于重置种群
        self.k = 0

    # 初始化种群
    def init_pop(self, pop=None):
        if pop != None:
            pop = np.array(pop)
            if len(pop.shape) != 3:
                raise Exception(f"Except get 3D vector, but get {len(pop.shape)} dim, please check it! ")
            self.pop = pop
        else:
            # self.pop = np.random.randint((0, 21), size=(self.pop_size, self.nums), dtype=np.float32) 
            # 随机抽样
            random_sample = np.random.uniform(bound[0][0], bound[0][1]-0.1, size=(self.pop_size, self.nums))
            self.pop = np.zeros((self.pop_size, self.nums, self.dna_size))
            for i in range(self.pop_size):
                for j in range(self.nums):
                    # 编码方式
                    # 将[0，21] -> [0, 2^dna_size]进行映射
                    # 所以单位长度是：((2**self.dna_size) / self.var_len[j])
                    # 长度是：(random_sample[i,j] - bound[j][0])
                    # 然后将映射后的数字转换成二进制，并转换成列表
                    num = int(round((random_sample[i,j] - bound[j][0]) * ((2**self.dna_size) / self.var_len[j])))
                    # if num >= 2**self.dna_size: num = 2**self.dna_size-1
                    #用python自带的格式化转化为前面空0的二进制字符串，然后拆分成列表
                    bin_num = ('{:0'+str(self.dna_size)+'b}').format(num) # 将num转换成dna_size位的二进制
                    self.pop[i,j] = np.array(list(bin_num)) # 将二进制转换成数组

    #将编码后的DNA翻译回来（解码）
    def translateDNA(self, index=None):
        W_vector = np.array([2**i for i in range(self.dna_size)]).reshape((self.dna_size, 1)) # shape = (dna_size, 1)
        if index != None:
            binary_vector = self.pop[index].dot(W_vector).reshape(-1) # shape = (pop_size, nums)
            for i in range(self.nums):
                binary_vector[i] = binary_vector[i] / ((2**self.dna_size)/self.var_len[0])
                binary_vector[i] = binary_vector[i] + self.bound[0][0]
            return binary_vector
            # print(binary_vector)
            # exit()
        else:
            # W_vector = [[0],[2],[4],[8],[16],[32] ..., [2^dna_size]]
            # (100, 4, 32) · (32, 1) = (100, 4, 1) # 先乘后加，然后就将dna转换成了一个十进制的数字
            binary_vector = self.pop.dot(W_vector).reshape(self.pop.shape[0:2]) # shape = (pop_size, nums)
            for i in range(self.pop_size):
                for j in range(self.nums):
                    binary_vector[i, j] = binary_vector[i, j] / ((2**self.dna_size)/self.var_len[j])
                    binary_vector[i, j] = binary_vector[i, j] + self.bound[j][0]
            return binary_vector

    #得到适应度
    def get_fitness(self):
        x = self.translateDNA()
        res = []
        for xx in x:
            temp = self.F(xx)[0]
            temp = temp if temp > 0 else 0
            res.append(temp)
        self.fitness = res
        return res
    
    #得到个体适应度
    def get_fitness_by_index(self, index=None):
        x = self.translateDNA(index)
        res = self.F(x)[0]
        return res

    #自然选择,选择较好的个体
    def select(self):
        fitness = self.get_fitness() # 获取适应度值
        # print(fitness)
        # exit()
        # np.random.choice
        # 第一个参数 a：随机数的来源，可以是数组、元组，但一定是一维的。
        # 第二个参数 size：随机数的个数，默认值为1
        # replace：元素是否可以重复采样
        # p：一个数组，大小与a相同，代表着每个元素选取的概率。默认每个概率一样
        # 假设种群大小为n，此代码为在0-n中选取n个元素， 可重复取值，根据适应度的概率进行随机抽样。
        choice = np.random.choice(np.arange(self.pop_size), size=self.pop_size, replace=True, p=fitness/np.sum(fitness)) # 筛选出优质种群
        self.pop = self.pop[choice] # 进行自然选择

    #染色体交叉，选择个体进行
    def crossover(self):
        # arr_max = heapq.nlargest(int(self.pop_size*0.1),self.fitness)#获取前五大的值并排序
        # best_ten_percent_people_index = list(map(self.fitness.index,arr_max))#获取前五大的值下标 
        for i, people in enumerate(self.pop):
            # if i in best_ten_percent_people_index: 
            #     continue 
            temp = people.copy()
            temp_fit = self.get_fitness_by_index(i)
            if np.random.rand() < self.cross_rate: # 染色体交叉概率
                mate = np.random.randint(0, self.pop.shape[0], size=1) # 选择一个个体进行染色体交叉
                cross_points = np.random.randint(0, 2, size=(self.nums ,self.dna_size)).astype(bool) # 随机生成要交叉的位，n条染色体，dna_size个碱基对
                people[cross_points] = self.pop[mate, cross_points] # 布尔索引然后进行染色体交叉
            cur_fit = self.get_fitness_by_index(i)
            if cur_fit < temp_fit: self.pop[i] = temp

    #基因变异,dna中的某个碱基发生变异
    def mutate(self):
        # arr_max = heapq.nlargest(int(self.pop_size*0.1),self.fitness)#获取前五大的值并排序
        # best_ten_percent_people_index = list(map(self.fitness.index,arr_max))#获取前五大的值下标 
        for i, people in enumerate(self.pop):
            # if i in best_ten_percent_people_index: 
            #     continue 
            temp = people.copy()
            temp_fit = self.get_fitness_by_index(i)
            for dna in people: # 多条DNA
                for point in range(self.dna_size): # 遍历碱基对
                    if np.random.rand() < self.mutation: # 变异概率
                        dna[point] = 1 if dna[point] == 0 else 0 # 基因变异
            cur_fit = self.get_fitness_by_index(i)
            if cur_fit < temp_fit: self.pop[i] = temp

    #打印当前状态日志
    def log(self):
        res = self.translateDNA()
        self.fitness = self.get_fitness()
        max_index = np.argmax(self.fitness)
        print(f"epoch:[{self.gen}/{self.epochs}], pop_size:{self.pop_size}, mean fit_val:{np.mean(self.fitness)}, max fit_val:{np.max(self.fitness)}, best x:{res[max_index]}")

    #进化，单次进化包括自然选择，染色体交叉和基因变异
    def evolution(self):
        for e in range(self.epochs):
            self.gen = e
            self.select()
            self.crossover()
            self.mutate()
            self.log()

    #重置
    def reset(self):
        self.pop = self.copy_pop.copy()

    # 设置k
    def set_k(self, k):
        self.k = k 

    # 记录最佳切分点
    def best_cp(self):
        res = self.translateDNA()
        self.fitness = self.get_fitness()
        max_index = np.argmax(self.fitness)
        return (self.fitness[max_index], res[max_index])

# 适应度函数
def Fun(x):
    # return -(x-1.2)**2 + 52.000000001
    return -(x-1.2)**2 + 52.0

if __name__ == "__main__":
    # 染色体个数，或者是求解问题的变量数
    nums = 1 
    # 每个变量的边界
    bound = [(-3,3)]*nums
    # 适应度值计算函数
    func = Fun 
    # DNA编码数量，就是一条DNA中碱基对的数量, 如果设置为None，就是自动适应
    dna_size = 128 
    # 交叉概率一般为0.4-0.99
    # 变异概率一般为0.001%-0.1%           
    # 增大变异概率可以防止早熟
    # 交叉运算决定了遗传算法的全局搜索能力
    # 变异运算决定了遗传算法的局部搜索能力
    cross_rate=0.8 # 交叉概率 0.8
    mutation=0.1 # 变异概率 0.003
    # 种群大小
    pop_size=500
    # 每一个k值遗传算法的迭代次数
    epochs=100

    # 创建遗传算法实例
    g = GA(nums, bound, func, dna_size, cross_rate, mutation, pop_size, epochs)
    g.evolution() # 进化算法

    print()
    print(f"Total cost time: {datetime.datetime.now() - g.startTime}.")