python学智能算法（五）|差分进化算法：原理认识和极小值分析

【1】引言

前序已经学习了模拟退火算法和遗传算法，相关文章链接为：

python学智能算法（一）|模拟退火算法：原理解释和最小值求解_模拟退火算法python-CSDN博客

python学智能算法（二）|模拟退火算法：进阶分析-CSDN博客

python学智能算法（三）|模拟退火算法：深层分析-CSDN博客

python学智能算法（四）|遗传算法：原理认识和极大值分析-CSDN博客

在此基础上，今天开始学习差分进化算法。

差分进化算法在遗传算法的基础上进行了新的变化，所以学习起来顺理成章。

【2】原理解释

在遗传算法中，数据的操作顺序为：选择，交叉，变异。这个顺序和生物的进化类似，实际操作也获得了很好的效果。

• 但差分进化算法不同，它的操作顺序几乎是完全反过来：
• 第一步，选定一个种群，对应一个多维数组；
• 第二步，在种群中挑选出三个成员，其中两个成员进行相减作差计算，这也是差分的本质。获得的差乘以一个系数以后和剩余成员叠加在一起。这一步的本质是“变异”。由于每个成员本质上代表的是一维数组，所以变异获得的依然是一维数组；
• 第三步，变异后的成员会和预定的目标作比较。成员和预定目标都是数组，比较的过程中成员的部分会留下，另一部分会被预定目标中的元素所取代，这一步的本质是“交叉”。交叉获得的也是一维数组。
• 第四步，交叉后的数组会被代入目标函数进行计算，如果获得的计算值相比之前变小，那这个数组和代入目标函数的计算值都会留下；如果计算值比最佳值还小，那计算值会直接赋给最佳值，数组也会赋给最佳解。这一步的本质是“选择”。选择获得更好的计算值和更好的解。
• 第五步，输出最佳值和最佳解。

【3】代码实现

为在代码上实现这个过程，需要进一步细化。

【3.1】准备工作

首先引入必要的模块：

import numpy as np #引入必要模块

【3.2】目标函数

# 定义目标函数，这里是求向量元素平方和
# 这个函数可以求一系列的x，所有自变量x会组成数组
def objective_function(x):
    return np.sum(x ** 2)

定义目标函数为一系列x的平方和。

【3.3】差分进化算法

# 差分进化算法实现
def differential_evolution(objective_function, bounds, pop_size=20, max_iter=100, F=0.5, CR=0.7):
    # bounds就是初始化的自变量取值范围，dim获取bounds的行数
    dim = len(bounds)
    # 初始化种群，随机生成在边界范围内的个体
    # 取bounds数组的第一列bounds[:, 0]为最小值界限
    # 取bounds数组的第二列bounds[:, 1]为最大值界限
    # 随机数的生成规则是从均匀分布中取值，先取列，population第i列的取值上下限，是由bounds的第i行确定
    # 按照均匀分布的原则，逐列随机生成位于左闭右开区间[bounds[:, 0],bounds[:, 1])上的随机数
    # dim本来是bounds的行数，现在转化为population的列数
    # 每一列的随机数重复取pop_size次，就组成了行列为(pop_size, dim)的新数组
    population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))
    print('population=',population)
    # 计算初始种群每个个体的适应度
    # 虽然population是一个(pop_size, dim)形式的数组，但经过计算之后，会获得pop_size个fitness
    # fitness在这被定义成一个数组，是一个一维的pop_size列数组
    fitness = np.array([objective_function(ind) for ind in population])
    # 找出初始种群中适应度最优（值最小）的个体索引
    # 这个索引只是一个整数，也就是fitness中最小值对应population中行的排序
    best_index = np.argmin(fitness)
    # 记录初始最优解
    # 按照索引直接获取最佳的解，这个解是一个行向量，一共有dim列
    best_solution = population[best_index]
    # 记录初始最优适应度
    # 按照索引直接获取最佳的适应度，这个适应度只有一个值
    best_fitness = fitness[best_index]

    # 以上部分为前置工作，相当于给定一些初始值
    # 以下部分为循环计算，会和前置工作获得的参数进行比较
    # 迭代优化
    for _ in range(max_iter):
        for i in range(pop_size):
            # 随机选择三个不同个体
            # 在0到pop_size-1的范围内，逐个取值
            # choices取值的组合形式为，不包含i的任何组合
            choices = [j for j in range(pop_size) if j != i]
            # 在上一步获得的choices中，随机且不重复地挑选3个数，命名为r1, r2, r3
            r1, r2, r3 = np.random.choice(choices, 3, replace=False)
            # 变异操作
            # 把上一步取得的三个数叠加在一起，后两个数求差再乘以一个系数，这就是差分
            # 实际操作中，是把第r2行的数据和第r3行的数据作差，然后乘以系数，再和第r1行的数据叠加
            # 所以新的mutant也是一个数组，population有dim列，mutant也有dim列
            mutant = population[r1] + F * (population[r2] - population[r3])
            # 边界处理，确保变异后的个体在边界内,最小值不会小于bounds[:, 0]，最大值不会超过bounds[:, 1]
            mutant = np.clip(mutant, bounds[:, 0], bounds[:, 1])
            # 交叉操作
            # 使用np.random.rand取dim个介于[0,1)之间的随机数
            # 当任何一个位置的随机数<CR时，这个位置会被赋值True，否则赋值False
            cross_points = np.random.rand(dim) < CR
            # 这是一个条件比较函数，对应cross_points上取True的位置，会给trial赋值mutant, 否则赋值population[i]
            # 这个population[i]会回到循环顶部的i，population在循环之前已经定义好了
            # trial也是一个一行dim列的数据
            trial = np.where(cross_points, mutant, population[i])
            # 计算试验个体的适应度
            # 上一步计算出来的trial是一个一行dim列的值，它们会被代入目标函数进行计算
            trial_fitness = objective_function(trial)
            # 选择操作
            if trial_fitness < fitness[i]:
                # 如果计算出来的数值比之前预先获得的适应度更小
                # 就将trail行全部赋值给population[i],计算出来的数值赋值给第i个适应度
                population[i] = trial
                fitness[i] = trial_fitness
                if trial_fitness < best_fitness:
                    # 如果计算所得的数值还比之前获得的最佳适应度小
                    # 就将计算所得值直接赋值给最佳适应度
                    # 将trial赋值给最佳解
                    best_fitness = trial_fitness
                    best_solution = trial
                # 由于上述的计算过程会循环pop_size次，随意获得的trial也会有(pop_size,dim)个
    #输出最佳解和最佳适应度
    return best_solution, best_fitness

代码可以分三部分：预设参数，代际循环和代内循环。

【3.3.1】预设参数

首先是预设参数：

def differential_evolution(objective_function, bounds, pop_size=20, max_iter=100, F=0.5, CR=0.7):
    # bounds就是初始化的自变量取值范围，dim获取bounds的行数
    dim = len(bounds)
    # 初始化种群，随机生成在边界范围内的个体
    # 取bounds数组的第一列bounds[:, 0]为最小值界限
    # 取bounds数组的第二列bounds[:, 1]为最大值界限
    # 随机数的生成规则是从均匀分布中取值，先取列，population第i列的取值上下限，是由bounds的第i行确定
    # 按照均匀分布的原则，逐列随机生成位于左闭右开区间[bounds[:, 0],bounds[:, 1])上的随机数
    # dim本来是bounds的行数，现在转化为population的列数
    # 每一列的随机数重复取pop_size次，就组成了行列为(pop_size, dim)的新数组
    population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))
    print('population=',population)
    # 计算初始种群每个个体的适应度
    # 虽然population是一个(pop_size, dim)形式的数组，但经过计算之后，会获得pop_size个fitness
    # fitness在这被定义成一个数组，是一个一维的pop_size列数组
    fitness = np.array([objective_function(ind) for ind in population])
    # 找出初始种群中适应度最优（值最小）的个体索引
    # 这个索引只是一个整数，也就是fitness中最小值对应population中行的排序
    best_index = np.argmin(fitness)
    # 记录初始最优解
    # 按照索引直接获取最佳的解，这个解是一个行向量，一共有dim列
    best_solution = population[best_index]
    # 记录初始最优适应度
    # 按照索引直接获取最佳的适应度，这个适应度只有一个值
    best_fitness = fitness[best_index]

    # 以上部分为前置工作，相当于给定一些初始值

预设参数部分，首先会读取一个预设数组的行数，预设数组类似于自变量。

求多个自变量的平方和，把这些自变量放在一起，就可以组成一个自变量数组。

# 初始化种群，随机生成在边界范围内的个体
# 取bounds数组的第一列bounds[:, 0]为最小值界限
# 取bounds数组的第二列bounds[:, 1]为最大值界限
# 随机数的生成规则是从均匀分布中取值，先取列，population第i列的取值上下限，是由bounds的第i行确定
# 按照均匀分布的原则，逐列随机生成位于左闭右开区间[bounds[:, 0],bounds[:, 1])上的随机数
# dim本来是bounds的行数，现在转化为population的列数
# 每一列的随机数重复取pop_size次，就组成了行列为(pop_size, dim)的新数组
population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))

然后用这个数组来生成新的随机数：

population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))

这行代码的pop_size是提前预设好的，相当于种群的大小；

随机数的生成是在自变量数组的范围内选定的，bounds[:, 0]代表的是第一列数据，作为取值的左区间，bounds[:, 1]代表的是第二列数据，作为取值的右区间。取到的随机数组成pop_size行dim列的population数组。

# 计算初始种群每个个体的适应度
# 虽然population是一个(pop_size, dim)形式的数组，但经过计算之后，会获得pop_size个fitness
# fitness在这被定义成一个数组，是一个一维的pop_size列数组
fitness = np.array([objective_function(ind) for ind in population])

之后定义了一个适应度数组fitness，它存储了将population数组中获得的数据逐行代入目标函数获得的计算值。 fitness是一个一维的pop_size列数组。

# 找出初始种群中适应度最优（值最小）的个体索引
# 这个索引只是一个整数，也就是fitness中最小值对应population中行的排序
best_index = np.argmin(fitness)

之后用np.argmin()函数，取到fitness中最小值对应的位置索引，比如最小值是第3个数，索引就是2（因为索引从0开始，第n个值对应的位置索引为n-1）。

# 记录初始最优解
# 按照索引直接获取最佳的解，这个解是一个行向量，一共有dim列
best_solution = population[best_index]

之后依据索引，直接找回population中计算获得这个最小fitness所对应的行。

# 记录初始最优适应度
# 按照索引直接获取最佳的适应度，这个适应度只有一个值
best_fitness = fitness[best_index]

同时将最佳的fitness通过索引记录下来。

索引是为了找出population中获得最小fitness所对应的行，所以为了获得索引，需要使用best_index = np.argmin(fitness)函数。这个过程也没有直接一步计算获得best_fitness，而是再次通过索引找出最小的fitness。

【3.3.2】循环

然后是循环，包括代际循环和代内循环，由于这两个循环嵌套在一起，所以要一起解读：

for _ in range(max_iter):
    for i in range(pop_size):
        # 随机选择三个不同个体
        # 在0到pop_size-1的范围内，逐个取值
        # choices取值的组合形式为，不包含i的任何组合
        choices = [j for j in range(pop_size) if j != i]
        # 在上一步获得的choices中，随机且不重复地挑选3个数，命名为r1, r2, r3
        r1, r2, r3 = np.random.choice(choices, 3, replace=False)
        # 变异操作
        # 把上一步取得的三个数叠加在一起，后两个数求差再乘以一个系数，这就是差分
        # 实际操作中，是把第r2行的数据和第r3行的数据作差，然后乘以系数，再和第r1行的数据叠加
        # 所以新的mutant也是一个数组，population有dim列，mutant也有dim列
        mutant = population[r1] + F * (population[r2] - population[r3])

首先看in range(pop_size)的代内循环，这里先输出0到pop_size-1的整数中随机除i之外的数字：

choices = [j for j in range(pop_size) if j != i]

需要注意的是，挑选出来的数不包括i。

然后再任选三个数：

            # 在上一步获得的choices中，随机且不重复地挑选3个数，命名为r1, r2, r3
            r1, r2, r3 = np.random.choice(choices, 3, replace=False)

之后进行变异操作：

            # 变异操作
            # 把上一步取得的三个数叠加在一起，后两个数求差再乘以一个系数，这就是差分
            # 实际操作中，是把第r2行的数据和第r3行的数据作差，然后乘以系数，再和第r1行的数据叠加
            # 所以新的mutant也是一个数组，population有dim列，mutant也有dim列
            mutant = population[r1] + F * (population[r2] - population[r3])

变异就是进行差分操作，上一步挑选出来的三个数，会分别对应population中的三行，其中两行进行差分后乘以系数叠加到剩余行。

# 边界处理，确保变异后的个体在边界内,最小值不会小于bounds[:, 0]，最大值不会超过bounds[:, 1]
            mutant = np.clip(mutant, bounds[:, 0], bounds[:, 1])

变异后可能会超出边界限制，所以还要进行边界约束。

            # 交叉操作
            # 使用np.random.rand取dim个介于[0,1)之间的随机数
            # 当任何一个位置的随机数<CR时，这个位置会被赋值True，否则赋值False
            cross_points = np.random.rand(dim) < CR
            # 这是一个条件比较函数，对应cross_points上取True的位置，会给trial赋值mutant, 否则赋值population[i]
            # 这个population[i]会回到循环顶部的i，population在循环之前已经定义好了
            # trial也是一个一行dim列的数据
            trial = np.where(cross_points, mutant, population[i])
            # 计算试验个体的适应度
            # 上一步计算出来的trial是一个一行dim列的值，它们会被代入目标函数进行计算
            trial_fitness = objective_function(trial)

之后进行的交叉操作需要放在一起看：

先用np.random.rand()函数生成一个基于[0,1)之间的dim列的随机数据数组，约束为dim列是因为population本身也是dim列。

然后将生成的数组和预先设定好的参数CR作比较，如果小于CR，在那个位置赋值true，反之赋值Flase，这也会组成一个dim列的一维数组。

之后基于上一步获得的True和False数组，调用np.where()函数，根据True和False所在的位置，对变异获得的数组mutant和预设获得的第i行数组population[i]进行交叉。

因为变异数组mutant要和第i行数组population[i]进行交叉，所以进行变异的时候，没有选择第i行数据参加变异。

之后将交叉获得的数组trial代入目标函数，会获得此时的适应度trial_fitness。

然后就是选择操作：

    # 选择操作
            if trial_fitness < fitness[i]:
                # 如果计算出来的数值比之前预先获得的适应度更小
                # 就将trail行全部赋值给population[i],计算出来的数值赋值给第i个适应度
                population[i] = trial
                fitness[i] = trial_fitness
                if trial_fitness < best_fitness:
                    # 如果计算所得的数值还比之前获得的最佳适应度小
                    # 就将计算所得值直接赋值给最佳适应度
                    # 将trial赋值给最佳解
                    best_fitness = trial_fitness
                    best_solution = trial
                # 由于上述的计算过程会循环pop_size次，所以获得的trial也会有(pop_size,dim)个

选择的原则是，新计算的适应度更小，则留下交叉数组和对应计算获得的适应度；如果新计算的适应度比最佳适应度还小，则继续将交叉数组定义为最佳解对应的数组，将对应计算获得的适应度定义为最佳适应度。

    return best_solution, best_fitness

之后将最佳值一起返回即可。

【3.4】算法调用

算法调用只需要给定自变量数组范围，直接调用函数即可：

# 定义变量边界，这里是二维问题，每个变量范围是 [-5, 5]，bounds是一个二行二列的数组
bounds = np.array([[-5, 5], [-5, 5]])
# 运行差分进化算法
best_sol, best_fit = differential_evolution(objective_function, bounds)
print("最优解:", best_sol)
print("最优适应度:", best_fit)

【3.5】运行效果

代码运行的效果为：

图1 运行效果

此时的完整代码为：

import numpy as np #引入必要模块

# 定义目标函数，这里是求向量元素平方和
# 这个函数可以求一系列的x，所有自变量x会组成数组
def objective_function(x):
    return np.sum(x ** 2)

# 差分进化算法实现
def differential_evolution(objective_function, bounds, pop_size=20, max_iter=100, F=0.5, CR=0.7):
    # bounds就是初始化的自变量取值范围，dim获取bounds的行数
    dim = len(bounds)
    # 初始化种群，随机生成在边界范围内的个体
    # 取bounds数组的第一列bounds[:, 0]为最小值界限
    # 取bounds数组的第二列bounds[:, 1]为最大值界限
    # 随机数的生成规则是从均匀分布中取值，先取列，population第i列的取值上下限，是由bounds的第i行确定
    # 按照均匀分布的原则，逐列随机生成位于左闭右开区间[bounds[:, 0],bounds[:, 1])上的随机数
    # dim本来是bounds的行数，现在转化为population的列数
    # 每一列的随机数重复取pop_size次，就组成了行列为(pop_size, dim)的新数组
    population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))
    print('population=',population)
    # 计算初始种群每个个体的适应度
    # 虽然population是一个(pop_size, dim)形式的数组，但经过计算之后，会获得pop_size个fitness
    # fitness在这被定义成一个数组，是一个一维的pop_size列数组
    fitness = np.array([objective_function(ind) for ind in population])
    # 找出初始种群中适应度最优（值最小）的个体索引
    # 这个索引只是一个整数，也就是fitness中最小值对应population中行的排序
    best_index = np.argmin(fitness)
    # 记录初始最优解
    # 按照索引直接获取最佳的解，这个解是一个行向量，一共有dim列
    best_solution = population[best_index]
    # 记录初始最优适应度
    # 按照索引直接获取最佳的适应度，这个适应度只有一个值
    best_fitness = fitness[best_index]

    # 以上部分为前置工作，相当于给定一些初始值
    # 以下部分为循环计算，会和前置工作获得的参数进行比较
    # 迭代优化
    for _ in range(max_iter):
        for i in range(pop_size):
            # 随机选择三个不同个体
            # 在0到pop_size-1的范围内，逐个取值
            # choices取值的组合形式为，不包含i的顺序组合
            choices = [j for j in range(pop_size) if j != i]
            # 在上一步获得的choices中，随机且不重复地挑选3个数，命名为r1, r2, r3
            r1, r2, r3 = np.random.choice(choices, 3, replace=False)
            # 变异操作
            # 把上一步取得的三个数叠加在一起，后两个数求差再乘以一个系数，这就是差分
            # 实际操作中，是把第r2行的数据和第r3行的数据作差，然后乘以系数，再和第r1行的数据叠加
            # 所以新的mutant也是一个数组，population有dim列，mutant也有dim列
            mutant = population[r1] + F * (population[r2] - population[r3])
            # 边界处理，确保变异后的个体在边界内,最小值不会小于bounds[:, 0]，最大值不会超过bounds[:, 1]
            mutant = np.clip(mutant, bounds[:, 0], bounds[:, 1])
            # 交叉操作
            # 使用np.random.rand取dim个介于[0,1)之间的随机数
            # 当任何一个位置的随机数<CR时，这个位置会被赋值True，否则赋值False
            cross_points = np.random.rand(dim) < CR
            # 这是一个条件比较函数，对应cross_points上取True的位置，会给trial赋值mutant, 否则赋值population[i]
            # 这个population[i]会回到循环顶部的i，population在循环之前已经定义好了
            # trial也是一个一行dim列的数据
            trial = np.where(cross_points, mutant, population[i])
            # 计算试验个体的适应度
            # 上一步计算出来的trial是一个一行dim列的值，它们会被代入目标函数进行计算
            trial_fitness = objective_function(trial)
            # 选择操作
            if trial_fitness < fitness[i]:
                # 如果计算出来的数值比之前预先获得的适应度更小
                # 就将trail行全部赋值给population[i],计算出来的数值赋值给第i个适应度
                population[i] = trial
                fitness[i] = trial_fitness
                if trial_fitness < best_fitness:
                    # 如果计算所得的数值还比之前获得的最佳适应度小
                    # 就将计算所得值直接赋值给最佳适应度
                    # 将trial赋值给最佳解
                    best_fitness = trial_fitness
                    best_solution = trial
                # 由于上述的计算过程会循环pop_size次，所以获得的trial也会有(pop_size,dim)个
    #输出最佳解和最佳适应度
    return best_solution, best_fitness

# 定义变量边界，这里是二维问题，每个变量范围是 [-5, 5]，bounds是一个二行二列的数组
bounds = np.array([[-5, 5], [-5, 5]])
# 运行差分进化算法
best_sol, best_fit = differential_evolution(objective_function, bounds)
print("最优解:", best_sol)
print("最优适应度:", best_fit)

【4】原理分析

在完成第一轮代码解析之后，可以再回头看差分进化算法代码：

# 差分进化算法实现
def differential_evolution(objective_function, bounds, pop_size=20, max_iter=100, F=0.5, CR=0.7):
    # bounds就是初始化的自变量取值范围，dim获取bounds的行数
    dim = len(bounds)
    # 初始化种群，随机生成在边界范围内的个体
    # 取bounds数组的第一列bounds[:, 0]为最小值界限
    # 取bounds数组的第二列bounds[:, 1]为最大值界限
    # 随机数的生成规则是从均匀分布中取值，先取列，population第i列的取值上下限，是由bounds的第i行确定
    # 按照均匀分布的原则，逐列随机生成位于左闭右开区间[bounds[:, 0],bounds[:, 1])上的随机数
    # dim本来是bounds的行数，现在转化为population的列数
    # 每一列的随机数重复取pop_size次，就组成了行列为(pop_size, dim)的新数组
    population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))
    print('population=',population)
    # 计算初始种群每个个体的适应度
    # 虽然population是一个(pop_size, dim)形式的数组，但经过计算之后，会获得pop_size个fitness
    # fitness在这被定义成一个数组，是一个一维的pop_size列数组
    fitness = np.array([objective_function(ind) for ind in population])
    # 找出初始种群中适应度最优（值最小）的个体索引
    # 这个索引只是一个整数，也就是fitness中最小值对应population中行的排序
    best_index = np.argmin(fitness)
    # 记录初始最优解
    # 按照索引直接获取最佳的解，这个解是一个行向量，一共有dim列
    best_solution = population[best_index]
    # 记录初始最优适应度
    # 按照索引直接获取最佳的适应度，这个适应度只有一个值
    best_fitness = fitness[best_index]

    # 以上部分为前置工作，相当于给定一些初始值
    # 以下部分为循环计算，会和前置工作获得的参数进行比较
    # 迭代优化
    for _ in range(max_iter):
        for i in range(pop_size):
            # 随机选择三个不同个体
            # 在0到pop_size-1的范围内，逐个取值
            # choices取值的组合形式为，不包含i的顺序组合
            choices = [j for j in range(pop_size) if j != i]
            # 在上一步获得的choices中，随机且不重复地挑选3个数，命名为r1, r2, r3
            r1, r2, r3 = np.random.choice(choices, 3, replace=False)
            # 变异操作
            # 把上一步取得的三个数叠加在一起，后两个数求差再乘以一个系数，这就是差分
            # 实际操作中，是把第r2行的数据和第r3行的数据作差，然后乘以系数，再和第r1行的数据叠加
            # 所以新的mutant也是一个数组，population有dim列，mutant也有dim列
            mutant = population[r1] + F * (population[r2] - population[r3])
            # 边界处理，确保变异后的个体在边界内,最小值不会小于bounds[:, 0]，最大值不会超过bounds[:, 1]
            mutant = np.clip(mutant, bounds[:, 0], bounds[:, 1])
            # 交叉操作
            # 使用np.random.rand取dim个介于[0,1)之间的随机数
            # 当任何一个位置的随机数<CR时，这个位置会被赋值True，否则赋值False
            cross_points = np.random.rand(dim) < CR
            # 这是一个条件比较函数，对应cross_points上取True的位置，会给trial赋值mutant, 否则赋值population[i]
            # 这个population[i]会回到循环顶部的i，population在循环之前已经定义好了
            # trial也是一个一行dim列的数据
            trial = np.where(cross_points, mutant, population[i])
            # 计算试验个体的适应度
            # 上一步计算出来的trial是一个一行dim列的值，它们会被代入目标函数进行计算
            trial_fitness = objective_function(trial)
            # 选择操作
            if trial_fitness < fitness[i]:
                # 如果计算出来的数值比之前预先获得的适应度更小
                # 就将trail行全部赋值给population[i],计算出来的数值赋值给第i个适应度
                population[i] = trial
                fitness[i] = trial_fitness
                if trial_fitness < best_fitness:
                    # 如果计算所得的数值还比之前获得的最佳适应度小
                    # 就将计算所得值直接赋值给最佳适应度
                    # 将trial赋值给最佳解
                    best_fitness = trial_fitness
                    best_solution = trial
                # 由于上述的计算过程会循环pop_size次，所以获得的trial也会有(pop_size,dim)个
    #输出最佳解和最佳适应度
    return best_solution, best_fitness

代码里面增加了大量的注释带辅助理解，总体上分三个部分的逻辑非常清晰：

预设参数部分获得最初的预设参数后；分别在代内循环和代际循环中不断更新种群每一行所对应的数值和最佳适应度；使用旧行差分再叠加的形式实现变异，经过概率的判断将变异行和原始行进行交叉，最后讲所获得结果按照目标原则进行选择。

差分进化算法相对于之前的模拟退火算法和遗传算法，在预设参数阶段是通过生成随机数的形式完成，相对来说减少了参数预设的难度。

【5】代码改写

由于差分进化算法的每一组最优解都是一维数组，所以这个算法更适用于计算多参数的函数。

此处以f (x，y) = 3cos(xy) + x + y为目标函数该写代码。

这里直接给出完整代码：

import numpy as np #引入必要模块
import  matplotlib.pyplot as plt

# 设置 matplotlib 支持中文
plt.rcParams['font.family'] = 'SimHei'
# 解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False

# 定义目标函数，这里是求向量元素平方和
# 这个函数可以求一系列的x，所有自变量x会组成数组
def objective_function(params):
    #return np.sum(x ** 2)
    x, y = params
    return 3 * np.cos(x * y) + x + y

# 差分进化算法实现
def differential_evolution(objective_function, bounds, pop_size=20, max_iter=100, F=0.5, CR=0.7):
    # bounds就是初始化的自变量取值范围，dim获取bounds的行数
    dim = len(bounds)
    # 初始化种群，随机生成在边界范围内的个体
    # 取bounds数组的第一列bounds[:, 0]为最小值界限
    # 取bounds数组的第二列bounds[:, 1]为最大值界限
    # 随机数的生成规则是从均匀分布中取值，先取列，population第i列的取值上下限，是由bounds的第i行确定
    # 按照均匀分布的原则，逐列随机生成位于左闭右开区间[bounds[:, 0],bounds[:, 1])上的随机数
    # dim本来是bounds的行数，现在转化为population的列数
    # 每一列的随机数重复取pop_size次，就组成了行列为(pop_size, dim)的新数组
    population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))
    #print('population=',population)
    # 计算初始种群每个个体的适应度
    # 虽然population是一个(pop_size, dim)形式的数组，但经过计算之后，会获得pop_size个fitness
    # fitness在这被定义成一个数组，是一个一维的pop_size列数组
    fitness = np.array([objective_function(ind) for ind in population])
    # 找出初始种群中适应度最优（值最小）的个体索引
    # 这个索引只是一个整数，也就是fitness中最小值对应population中行的排序
    best_index = np.argmin(fitness)
    # 记录初始最优解
    # 按照索引直接获取最佳的解，这个解是一个行向量，一共有dim列
    best_solution = population[best_index]
    # 记录初始最优适应度
    # 按照索引直接获取最佳的适应度，这个适应度只有一个值
    best_fitness = fitness[best_index]

    best_solution_history=[]
    best_fitness_history=[]
    # 以上部分为前置工作，相当于给定一些初始值
    # 以下部分为循环计算，会和前置工作获得的参数进行比较
    # 迭代优化
    for _ in range(max_iter):
        for i in range(pop_size):
            # 随机选择三个不同个体
            # 在0到pop_size-1的范围内，逐个取值
            # choices取值的组合形式为，不包含i的顺序组合
            choices = [j for j in range(pop_size) if j != i]
            # 在上一步获得的choices中，随机且不重复地挑选3个数，命名为r1, r2, r3
            r1, r2, r3 = np.random.choice(choices, 3, replace=False)
            # 变异操作
            # 把上一步取得的三个数叠加在一起，后两个数求差再乘以一个系数，这就是差分
            # 实际操作中，是把第r2行的数据和第r3行的数据作差，然后乘以系数，再和第r1行的数据叠加
            # 所以新的mutant也是一个数组，population有dim列，mutant也有dim列
            mutant = population[r1] + F * (population[r2] - population[r3])
            # 边界处理，确保变异后的个体在边界内,最小值不会小于bounds[:, 0]，最大值不会超过bounds[:, 1]
            mutant = np.clip(mutant, bounds[:, 0], bounds[:, 1])
            # 交叉操作
            # 使用np.random.rand取dim个介于[0,1)之间的随机数
            # 当任何一个位置的随机数<CR时，这个位置会被赋值True，否则赋值False
            cross_points = np.random.rand(dim) < CR
            # 这是一个条件比较函数，对应cross_points上取True的位置，会给trial赋值mutant, 否则赋值population[i]
            # 这个population[i]会回到循环顶部的i，population在循环之前已经定义好了
            # trial也是一个一行dim列的数据
            trial = np.where(cross_points, mutant, population[i])
            # 计算试验个体的适应度
            # 上一步计算出来的trial是一个一行dim列的值，它们会被代入目标函数进行计算
            trial_fitness = objective_function(trial)
            # 选择操作
            if trial_fitness < fitness[i]:
                # 如果计算出来的数值比之前预先获得的适应度更小
                # 就将trail行全部赋值给population[i],计算出来的数值赋值给第i个适应度
                population[i] = trial
                fitness[i] = trial_fitness
                if trial_fitness < best_fitness:
                    # 如果计算所得的数值还比之前获得的最佳适应度小
                    # 就将计算所得值直接赋值给最佳适应度
                    # 将trial赋值给最佳解
                    best_fitness = trial_fitness
                    best_solution = trial
                    #best_fitness_history.append(best_fitness)
                # 由于上述的计算过程会循环pop_size次，所以获得的trial也会有(pop_size,dim)个
        best_solution_history.append(best_solution)
        best_fitness_history.append(best_fitness)
    #输出最佳解和最佳适应度
    return best_solution, best_fitness,best_solution_history,best_fitness_history

# 定义变量边界，这里是二维问题，每个变量范围是 [-5, 5]，bounds是一个二行二列的数组
bounds = np.array([[-4, 4], [-4, 4]])
# 运行差分进化算法
best_solution,best_fitness,best_solution_history,best_fitness_history = differential_evolution(objective_function, bounds)
print("最优解:", best_solution)
print("最优适应度:", best_fitness)
print(best_fitness_history)

# 绘制目标函数的三维图
x = np.linspace(bounds[0, 0], bounds[0, 1], 100)
y = np.linspace(bounds[1, 0], bounds[1, 1], 100)
X, Y = np.meshgrid(x, y)
Z = np.zeros_like(X)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        Z[i, j] = objective_function([X[i, j], Y[i, j]])

# 创建一个包含两个子图的图形
#fig = plt.figure(figsize=(12, 5))
fig = plt.figure(figsize=(12, 5))

# 第一个子图：三维图
ax1 = fig.add_subplot(131, projection='3d')
ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('目标函数值')
ax1.set_title('目标函数三维图')

# 第二个子图：三维图
ax2= fig.add_subplot(132, projection='3d')
ax2.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
ax2.set_zlabel('目标函数值')
ax2.set_title('最佳目标函数搜索过程')
# 绘制搜索过程
search_history = np.array(best_solution_history)
#ax3.scatter(search_history[:, 0], search_history[:, 1], [objective_function(sol) for sol in search_history], c='r', marker='o')
ax2.scatter(search_history[:, 0], search_history[:, 1], [objective_function(sol) for sol in search_history], c='r', marker='o')

# 第三个子图：最优适应度历史变化
ax3 = fig.add_subplot(133, projection='3d')
# 绘制搜索过程
search_history = np.array(best_solution_history)
#ax3.scatter(search_history[:, 0], search_history[:, 1], [objective_function(sol) for sol in search_history], c='r', marker='o')
ax3.scatter(search_history[:, 0], search_history[:, 1], [objective_function(sol) for sol in search_history], c='r', marker='o')
ax3.set_xlabel('X')
ax3.set_ylabel('Y')
ax3.set_zlabel('最佳目标函数值')
ax3.set_title('最佳目标函数')

print(len(best_fitness_history))
plt.tight_layout()
plt.show()

相对于之前的代码，新代码的修改主要在于输出了所有best_solution_history和best_fitness_history所代表的每一代best_solution和best_fitness。

然后使用多子图的形式将详细的搜索过程绘制了下来：

图2 多子图展示差分进化算法搜索过程

图2中间小图和右边小图的红色点本来是一样的分布格局，因为坐标轴范围的不同，所以看起来有差异。

此时控制台的输出为：

图3 控制台输出

另外一种流行的绘图方式为，绘制迭代次数和最佳目标函数关系图：

import numpy as np #引入必要模块
import  matplotlib.pyplot as plt

# 设置 matplotlib 支持中文
plt.rcParams['font.family'] = 'SimHei'
# 解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False

# 定义目标函数，这里是求向量元素平方和
# 这个函数可以求一系列的x，所有自变量x会组成数组
def objective_function(params):
    #return np.sum(x ** 2)
    x, y = params
    return 3 * np.cos(x * y) + x + y

# 差分进化算法实现
def differential_evolution(objective_function, bounds, pop_size=20, max_iter=100, F=0.5, CR=0.7):
    # bounds就是初始化的自变量取值范围，dim获取bounds的行数
    dim = len(bounds)
    # 初始化种群，随机生成在边界范围内的个体
    # 取bounds数组的第一列bounds[:, 0]为最小值界限
    # 取bounds数组的第二列bounds[:, 1]为最大值界限
    # 随机数的生成规则是从均匀分布中取值，先取列，population第i列的取值上下限，是由bounds的第i行确定
    # 按照均匀分布的原则，逐列随机生成位于左闭右开区间[bounds[:, 0],bounds[:, 1])上的随机数
    # dim本来是bounds的行数，现在转化为population的列数
    # 每一列的随机数重复取pop_size次，就组成了行列为(pop_size, dim)的新数组
    population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))
    #print('population=',population)
    # 计算初始种群每个个体的适应度
    # 虽然population是一个(pop_size, dim)形式的数组，但经过计算之后，会获得pop_size个fitness
    # fitness在这被定义成一个数组，是一个一维的pop_size列数组
    fitness = np.array([objective_function(ind) for ind in population])
    # 找出初始种群中适应度最优（值最小）的个体索引
    # 这个索引只是一个整数，也就是fitness中最小值对应population中行的排序
    best_index = np.argmin(fitness)
    # 记录初始最优解
    # 按照索引直接获取最佳的解，这个解是一个行向量，一共有dim列
    best_solution = population[best_index]
    # 记录初始最优适应度
    # 按照索引直接获取最佳的适应度，这个适应度只有一个值
    best_fitness = fitness[best_index]

    best_solution_history=[]
    best_fitness_history=[]
    # 以上部分为前置工作，相当于给定一些初始值
    # 以下部分为循环计算，会和前置工作获得的参数进行比较
    # 迭代优化
    for _ in range(max_iter):
        for i in range(pop_size):
            # 随机选择三个不同个体
            # 在0到pop_size-1的范围内，逐个取值
            # choices取值的组合形式为，不包含i的顺序组合
            choices = [j for j in range(pop_size) if j != i]
            # 在上一步获得的choices中，随机且不重复地挑选3个数，命名为r1, r2, r3
            r1, r2, r3 = np.random.choice(choices, 3, replace=False)
            # 变异操作
            # 把上一步取得的三个数叠加在一起，后两个数求差再乘以一个系数，这就是差分
            # 实际操作中，是把第r2行的数据和第r3行的数据作差，然后乘以系数，再和第r1行的数据叠加
            # 所以新的mutant也是一个数组，population有dim列，mutant也有dim列
            mutant = population[r1] + F * (population[r2] - population[r3])
            # 边界处理，确保变异后的个体在边界内,最小值不会小于bounds[:, 0]，最大值不会超过bounds[:, 1]
            mutant = np.clip(mutant, bounds[:, 0], bounds[:, 1])
            # 交叉操作
            # 使用np.random.rand取dim个介于[0,1)之间的随机数
            # 当任何一个位置的随机数<CR时，这个位置会被赋值True，否则赋值False
            cross_points = np.random.rand(dim) < CR
            # 这是一个条件比较函数，对应cross_points上取True的位置，会给trial赋值mutant, 否则赋值population[i]
            # 这个population[i]会回到循环顶部的i，population在循环之前已经定义好了
            # trial也是一个一行dim列的数据
            trial = np.where(cross_points, mutant, population[i])
            # 计算试验个体的适应度
            # 上一步计算出来的trial是一个一行dim列的值，它们会被代入目标函数进行计算
            trial_fitness = objective_function(trial)
            # 选择操作
            if trial_fitness < fitness[i]:
                # 如果计算出来的数值比之前预先获得的适应度更小
                # 就将trail行全部赋值给population[i],计算出来的数值赋值给第i个适应度
                population[i] = trial
                fitness[i] = trial_fitness
                if trial_fitness < best_fitness:
                    # 如果计算所得的数值还比之前获得的最佳适应度小
                    # 就将计算所得值直接赋值给最佳适应度
                    # 将trial赋值给最佳解
                    best_fitness = trial_fitness
                    best_solution = trial
                    #best_fitness_history.append(best_fitness)
                # 由于上述的计算过程会循环pop_size次，所以获得的trial也会有(pop_size,dim)个
        best_solution_history.append(best_solution)
        best_fitness_history.append(best_fitness)
    #输出最佳解和最佳适应度
    return best_solution, best_fitness,best_solution_history,best_fitness_history

# 定义变量边界，这里是二维问题，每个变量范围是 [-5, 5]，bounds是一个二行二列的数组
bounds = np.array([[-4, 4], [-4, 4]])
# 运行差分进化算法
best_solution,best_fitness,best_solution_history,best_fitness_history = differential_evolution(objective_function, bounds)
print("最优解:", best_solution)
print("最优适应度:", best_fitness)
print(best_fitness_history)

# 创建一个新图形
fig,ax = plt.subplots(figsize=(12, 5))
ax.plot(best_fitness_history,label='最优适应度')
ax.set_xlabel('迭代次数')
ax.set_ylabel('最优适应度')
ax.set_title('最优适应度随迭代次数的变化')
ax.legend()

print(len(best_fitness_history))
plt.tight_layout()
plt.show()

获得的图像为：

图4 迭代次数和最佳目标函数关系图 

此外也可以只绘制目标函数和最优解：

import numpy as np #引入必要模块
import  matplotlib.pyplot as plt

# 设置 matplotlib 支持中文
plt.rcParams['font.family'] = 'SimHei'
# 解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False

# 定义目标函数，这里是求向量元素平方和
# 这个函数可以求一系列的x，所有自变量x会组成数组
def objective_function(params):
    #return np.sum(x ** 2)
    x, y = params
    return 3 * np.cos(x * y) + x + y

# 差分进化算法实现
def differential_evolution(objective_function, bounds, pop_size=20, max_iter=100, F=0.5, CR=0.7):
    # bounds就是初始化的自变量取值范围，dim获取bounds的行数
    dim = len(bounds)
    # 初始化种群，随机生成在边界范围内的个体
    # 取bounds数组的第一列bounds[:, 0]为最小值界限
    # 取bounds数组的第二列bounds[:, 1]为最大值界限
    # 随机数的生成规则是从均匀分布中取值，先取列，population第i列的取值上下限，是由bounds的第i行确定
    # 按照均匀分布的原则，逐列随机生成位于左闭右开区间[bounds[:, 0],bounds[:, 1])上的随机数
    # dim本来是bounds的行数，现在转化为population的列数
    # 每一列的随机数重复取pop_size次，就组成了行列为(pop_size, dim)的新数组
    population = np.random.uniform(bounds[:, 0], bounds[:, 1], (pop_size, dim))
    #print('population=',population)
    # 计算初始种群每个个体的适应度
    # 虽然population是一个(pop_size, dim)形式的数组，但经过计算之后，会获得pop_size个fitness
    # fitness在这被定义成一个数组，是一个一维的pop_size列数组
    fitness = np.array([objective_function(ind) for ind in population])
    # 找出初始种群中适应度最优（值最小）的个体索引
    # 这个索引只是一个整数，也就是fitness中最小值对应population中行的排序
    best_index = np.argmin(fitness)
    # 记录初始最优解
    # 按照索引直接获取最佳的解，这个解是一个行向量，一共有dim列
    best_solution = population[best_index]
    # 记录初始最优适应度
    # 按照索引直接获取最佳的适应度，这个适应度只有一个值
    best_fitness = fitness[best_index]

    best_solution_history=[]
    best_fitness_history=[]
    # 以上部分为前置工作，相当于给定一些初始值
    # 以下部分为循环计算，会和前置工作获得的参数进行比较
    # 迭代优化
    for _ in range(max_iter):
        for i in range(pop_size):
            # 随机选择三个不同个体
            # 在0到pop_size-1的范围内，逐个取值
            # choices取值的组合形式为，不包含i的顺序组合
            choices = [j for j in range(pop_size) if j != i]
            # 在上一步获得的choices中，随机且不重复地挑选3个数，命名为r1, r2, r3
            r1, r2, r3 = np.random.choice(choices, 3, replace=False)
            # 变异操作
            # 把上一步取得的三个数叠加在一起，后两个数求差再乘以一个系数，这就是差分
            # 实际操作中，是把第r2行的数据和第r3行的数据作差，然后乘以系数，再和第r1行的数据叠加
            # 所以新的mutant也是一个数组，population有dim列，mutant也有dim列
            mutant = population[r1] + F * (population[r2] - population[r3])
            # 边界处理，确保变异后的个体在边界内,最小值不会小于bounds[:, 0]，最大值不会超过bounds[:, 1]
            mutant = np.clip(mutant, bounds[:, 0], bounds[:, 1])
            # 交叉操作
            # 使用np.random.rand取dim个介于[0,1)之间的随机数
            # 当任何一个位置的随机数<CR时，这个位置会被赋值True，否则赋值False
            cross_points = np.random.rand(dim) < CR
            # 这是一个条件比较函数，对应cross_points上取True的位置，会给trial赋值mutant, 否则赋值population[i]
            # 这个population[i]会回到循环顶部的i，population在循环之前已经定义好了
            # trial也是一个一行dim列的数据
            trial = np.where(cross_points, mutant, population[i])
            # 计算试验个体的适应度
            # 上一步计算出来的trial是一个一行dim列的值，它们会被代入目标函数进行计算
            trial_fitness = objective_function(trial)
            # 选择操作
            if trial_fitness < fitness[i]:
                # 如果计算出来的数值比之前预先获得的适应度更小
                # 就将trail行全部赋值给population[i],计算出来的数值赋值给第i个适应度
                population[i] = trial
                fitness[i] = trial_fitness
                if trial_fitness < best_fitness:
                    # 如果计算所得的数值还比之前获得的最佳适应度小
                    # 就将计算所得值直接赋值给最佳适应度
                    # 将trial赋值给最佳解
                    best_fitness = trial_fitness
                    best_solution = trial
                    #best_fitness_history.append(best_fitness)
                # 由于上述的计算过程会循环pop_size次，所以获得的trial也会有(pop_size,dim)个
        best_solution_history.append(best_solution)
        best_fitness_history.append(best_fitness)
    #输出最佳解和最佳适应度
    return best_solution, best_fitness,best_solution_history,best_fitness_history

# 定义变量边界，这里是二维问题，每个变量范围是 [-5, 5]，bounds是一个二行二列的数组
bounds = np.array([[-4, 4], [-4, 4]])
# 运行差分进化算法
best_solution,best_fitness,best_solution_history,best_fitness_history = differential_evolution(objective_function, bounds)
print("最优解:", best_solution)
print("最优适应度:", best_fitness)
print(best_fitness_history)

# 绘制目标函数的三维图
x = np.linspace(bounds[0, 0], bounds[0, 1], 100)
y = np.linspace(bounds[1, 0], bounds[1, 1], 100)
X, Y = np.meshgrid(x, y)
Z = np.zeros_like(X)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        Z[i, j] = objective_function([X[i, j], Y[i, j]])

#绘制三维目标函数和最优解
fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
ax.scatter(best_solution[0], best_solution[1],best_fitness, c='r',s=100, marker='o')
ax.set_title('目标函数及最优解')
ax.set_xlabel('X')
ax.set_ylabel('Y')
plt.show()

图5 目标函数和最优解

至此，实际运算效果证明了差分进化算法的准确性。

【6】总结

学习了差分进化算的基本原理，掌握了使用python应用差分进化算法求解方程极限值的技巧。