差分进化算法原理及其MATLAB/Python代码

1.算法简介

引用自：Storn R, Price K. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces[J]. Journal of global optimization, 1997, 11: 341-359.

今天给大家带来的是一个非常经典的智能优化算法–差分进化算法(DE)，说它是现有各类优化算法的鼻祖之一也不为过(除此之外还有PSO、GA、蚁群优化算法ACO等)。小编在本科阶段打数学建模竞赛最钟情的算法就是本文介绍的差分进化算法DE，即使它似乎没那么新颖，但小编认为其核心思想与独特的算法设计值得每一位研究优化算法的同学学习，或许会启发你的idea哦~

首先简单介绍一下本文的主角，差分进化算法（Differential Evolution，简称DE）是一种用于求解优化问题的现代启发式算法，尤其擅长于解决高维、非线性、多模态的连续优化问题。DE是由Rainer Storn和Kenneth Price在1995年提出的，它是基于群体智能的策略，类似于遗传算法，但有着自己独特的优势和简化的设计。

2.整体架构流程

在了解差分进化算法之前我建议各位可以看一下有关遗传算法(GA)的原理，因为差分进化算法是基于遗传算法(GA)的改进版，这里我放个我看着还不错的讲解链接：遗传算法原理及其matlab程序实现

相对于遗传而言，其主要改进如下：(遗传算法简称GA，差分进化简称DE)

• 编码方式：GA使用二进制编码，DE使用实数编码，DE能够直接处理连续优化问题，且无需进行编码和解码的过程。
• 变异操作：GA的变异通过随机改变染色体上的基因来实现。DE中的变异是通过加权差分向量来生成新个体，具有更强的方向性和探索性。
• 交叉操作：GA中的交叉操作类似于生物遗传学中的有性繁殖。DE中的交叉则更多地表现为参数混合，即在一定概率下选择突变向量或目标向量的参数，这一过程通常称为“差分交叉”。
• 参数控制：DE控制参数较少，而GA多，DE的调参更简单。

2.1 差分进化（Differential Evolution, DE）算法概述

差分进化（DE）是一种并行直接搜索方法，它利用NP个D维参数向量作为每一代的种群，形式如下：

$$x_i^G;i=1;2;...;NP(1)$$

其中，(G)代表代数，(NP)在整个最小化过程中保持不变。初始向量种群随机选取，应当覆盖整个参数空间。通常假定所有随机决策遵循均匀概率分布，除非另有说明。如果有初步解存在，初始种群可以通过在标称解(x_{nom};0)上加上正态分布的随机偏差来生成。

DE通过将两个种群向量之差的加权值加上第三个向量来生成新的参数向量，这一操作称为突变。突变向量的参数随后与另一个预设向量（目标向量）的参数混合，产生所谓的试用向量。这一参数混合过程常被称为交叉（crossover）。如果试用向量产生的代价函数值小于目标向量，则试用向量将在下一代中取代目标向量，这一过程称为选择。每个种群向量在一代中都必须至少一次作为目标向量，因此一代中有(NP)次竞争。

2.2 DE的基本策略具体描述如下：

2.2.1 突变

对于每个目标向量(x_i^G ; i = 1 ; 2 ; 3 ; . . . ; NP)，根据下式生成突变向量：

$$v_i^{G+1}=x_{r1}^G+F\cdot (x_{r2}^G-x_{r3}^G)(2)$$

2.2.2 交叉

为了增加扰动参数向量的多样性，引入了交叉。为此，形成试用向量：

$$u_i^{G+1}=(u_{1i}^{G+1};u_{2i}^{G+1};...;u_{Di}^{G+1})$$

其中，

$$u_j^{i;G+1}=\{\begin{array}{ll}{v}_j^{i;G+1}& \text{if }(r\text{ andb}(j)<C_R)\text{ or }j=r_{nbr}(i)\\ x_j^{i;G}& \text{if }(r\text{ andb}(j)>C_R)\text{ and }j\ne r_{nbr}(i)\end{array}$$

2.2.3 选择

为了决定是否成为下一代(G + 1)的成员，试用向量(u_i^{{G+1})与目标向量(x_i}G)采用贪心准则进行比较。如果(u_i^{{G+1})产生的代价函数值小于(x_i}G)，则(x_i^{{G+1})设置为(u_i}{G+1})；否则，保留旧值(x_i^G)。

3.伪代码

原论文中使用了基于C语言格式的伪代码：

while (count < gen_max):  # Halt after gen_max generations.
    for (i=0; i<NP; i++):  # Start loop through population.
        a=rnd_uni() * NP; while (a==i);  # Randomly pick 3 vectors, all different.
        b=rnd_uni() * NP; while ((b==i || b==a));
        c=rnd_uni() * NP; while ((c==i || c==a || c==b));
        j=rnd_uni() * D;  # Randomly pick the first parameter.

        for (k=1; k<D; k++):  # Load D parameters into trial[].
            if (rnd_uni() < CR || k==D):
                trial[j]=x1[a][j]+F*(x1[b][j]-x1[c][j]);  # Source for trial[j] is a random vector plus weighted differential.
            else:
                trial[j]=x1[i][j];  # Trial parameter comes from target vector.
            j=(j+1)%D;  # Get next parameter, modulo D.

        score=evaluate(trial);  # Evaluate trial with your function.
        if (score<=cost[i]):  # If trial[] improves on x1[i][].
            for (j=0; j<D; j++) x2[i][j]=trial[j];  # Move trial[] to secondary array and store improved cost.
        else:
            for (j=0; j<D; j++) x2[i][j]=x1[i][j];  # Otherwise, move x1[i][] to secondary array.

    for (i=0; i<NP; i++):  # After each generation.
        for (j=0; j<D; j++) x1[i][j]=x2[i][j];  # Move secondary array into primary array.
    count++;  # End of generation...increment counter.

这个伪代码展示了差分进化算法的核心循环，其中包括了三个主要步骤：突变/重组、评估/选择以及种群更新。在这个循环中，首先从当前种群中随机选择三个不同的个体（a、b、c），然后对每个维度的参数进行突变和重组，生成一个新的候选解（trial）。接着，评估候选解的适应度，并与当前个体的适应度进行比较，若候选解更好，则将其替换掉当前个体。最后，在每一代结束时，将辅助数组中的个体复制回主数组，以便开始新的一代。计数器count记录了已经完成的迭代次数，当达到最大迭代次数gen_max时，算法停止。

4. MATLAB与Python代码

4.1 MATLAB 函数

function [best_solution, best_cost] = differential_evolution(evaluate_func, D, NP, F, CR, gen_max)
    % 差分进化算法
    % 输入参数:
    %   evaluate_func: 评估函数，用于计算解的适应度
    %   D: 问题的维度
    %   NP: 种群大小
    %   F: 缩放因子
    %   CR: 交叉率
    %   gen_max: 最大迭代次数
    % 输出参数:
    %   best_solution: 最优解
    %   best_cost: 最优解的适应度值

    % 初始化种群
    x1 = rand(NP, D);  % 随机生成初始种群
    x2 = zeros(NP, D);  % 用于存储新一代种群
    cost = zeros(NP, 1);  % 存储每个个体的适应度值
    
    % 评估初始种群
    for i = 1:NP
        cost(i) = evaluate_func(x1(i, :));
    end
    
    count = 0;  % 迭代计数器
    while count < gen_max
        for i = 1:NP  % 对每个个体进行操作
            % 随机选择三个不同的向量
            a = randi(NP);
            while a == i
                a = randi(NP);
            end
            
            b = randi(NP);
            while b == i || b == a
                b = randi(NP);
            end
            
            c = randi(NP);
            while c == i || c == a || c == b
                c = randi(NP);
            end
            
            % 创建试验向量
            j = randi(D);  % 随机选择一个维度开始
            trial = x1(i, :);  % 初始化试验向量
            for k = 1:D
                if rand() < CR || k == D  % 判断是否进行交叉
                    % 差分变异
                    trial(j) = x1(a, j) + F * (x1(b, j) - x1(c, j));
                end
                j = mod(j, D) + 1;  % 移动到下一个维度
            end
            
            % 评估试验向量
            trial_cost = evaluate_func(trial);
            
            % 选择：如果试验向量更优，则替换原个体
            if trial_cost <= cost(i)
                x2(i, :) = trial;
                cost(i) = trial_cost;
            else
                x2(i, :) = x1(i, :);
            end
        end
        
        % 更新种群
        x1 = x2;
        count = count + 1;  % 迭代次数加1
    end
    
    % 找到最优解
    [best_cost, idx] = min(cost);
    best_solution = x1(idx, :);
end

4.2 Python 函数

import numpy as np

def differential_evolution(evaluate_func, D, NP, F, CR, gen_max):
    """
    差分进化算法
    
    参数:
    evaluate_func : 函数
        评估函数，用于计算解的适应度
    D : int
        问题的维度
    NP : int
        种群大小
    F : float
        缩放因子
    CR : float
        交叉率
    gen_max : int
        最大迭代次数
    
    返回:
    best_solution : numpy.ndarray
        最优解
    best_cost : float
        最优解的适应度值
    """
    
    # 初始化种群
    x1 = np.random.rand(NP, D)  # 随机生成初始种群
    x2 = np.zeros((NP, D))  # 用于存储新一代种群
    cost = np.zeros(NP)  # 存储每个个体的适应度值
    
    # 评估初始种群
    for i in range(NP):
        cost[i] = evaluate_func(x1[i])
    
    count = 0  # 迭代计数器
    while count < gen_max:
        for i in range(NP):  # 对每个个体进行操作
            # 随机选择三个不同的向量
            a, b, c = np.random.choice(np.delete(np.arange(NP), i), 3, replace=False)
            
            # 创建试验向量
            j = np.random.randint(D)  # 随机选择一个维度开始
            trial = x1[i].copy()  # 初始化试验向量
            for k in range(D):
                if np.random.rand() < CR or k == D-1:  # 判断是否进行交叉
                    # 差分变异
                    trial[j] = x1[a, j] + F * (x1[b, j] - x1[c, j])
                j = (j + 1) % D  # 移动到下一个维度
            
            # 评估试验向量
            trial_cost = evaluate_func(trial)
            
            # 选择：如果试验向量更优，则替换原个体
            if trial_cost <= cost[i]:
                x2[i] = trial
                cost[i] = trial_cost
            else:
                x2[i] = x1[i]
        
        # 更新种群
        x1 = x2.copy()
        count += 1  # 迭代次数加1
    
    # 找到最优解
    best_idx = np.argmin(cost)
    best_solution = x1[best_idx]
    best_cost = cost[best_idx]
    
    return best_solution, best_cost