JADE: Adaptive Differential Evolution with Optional External Archive

JADE

在本节中,我们提出一种新的差分进化算法JADE。该算法采用"DE/current-to-pbest"变异策略(带有可选存档机制),并通过自适应方式调整控制参数F和CR。JADE沿用了第II节公式(4)和(5)所描述的交叉与选择操作。
在这里插入图片描述
在这里插入图片描述

A. DE/Current-to-pbest

DE/rand/1是差分进化算法中最早提出的变异策略[1][2],也被认为是文献中最为成功且应用最广泛的方案[31]。然而,文献[6]指出DE/best/2可能具有优于DE/rand/1的特性,文献[32]则表明DE/best/1在多数技术问题的研究中表现更优。此外,文献[33]的作者认为融入最优解信息是有益的,并在其算法中采用了DE/current-to-best/1策略。相较于DE/rand/k类策略,诸如DE/current-to-best/k和DE/best/k这类贪婪策略通过融入最优解信息,在进化搜索中获得了更快的收敛速度。然而,这种最优解信息也可能导致种群多样性下降,进而引发早熟收敛等问题。
针对现有贪婪策略(式2和式3)收敛速度快但可靠性不足的问题,本文提出一种新型变异策略——“DE/current-to-pbest with optional archive”,作为本文提出的自适应DE算法JADE的核心基础。
在这里插入图片描述
如图1所示,在DE/current-to-pbest/1策略(无存档机制)中,变异向量按如下方式生成:
在这里插入图片描述
其中, x b e s t , g p x_\mathrm{best,g}^p xbest,gp 是从当前种群中前 100 p % 100p\% 100p% p ∈ ( 0 , 1 ] p\in(0,1] p(0,1])的精英个体中随机选取的, F i F_i Fi是与个体 x i x_i xi 相关联的变异因子,其值通过后文式(10)提出的自适应机制在每一代重新生成。DE/current-to-pbest本质上是DE/current-to-best策略的推广形式——通过随机选取前 100 p % 100p\% 100p%精英个体中的任一解,取代DE/current-to-best中单一全局最优解的角色。
新近探索的劣势解相较于当前种群,能够为进化方向提供额外的有效信息。设A为存档的劣势解集合,P为当前种群。在带存档机制的DE/current-to-pbest/1策略中,变异向量生成方式如下:
在这里插入图片描述
其中, x i , g x_{i,g} xi,g 𝑥 𝑟 1 , 𝑔 𝑥_{𝑟1,𝑔} xr1,g 𝑥 p b e s t , 𝑔 𝑥_{pbest,𝑔} xpbest,g的选取方式与式(6)中相同(均从当前种群P中选择),而 x ~ r 2 , g \tilde{\mathbf{x}}_{r2,g} x~r2,g 则是从当前种群P与存档A的并集𝑃∪𝐴中随机选取的。

存档操作的设计力求简洁,以避免产生显著的计算开销。存档初始化为空集,在每代进化过程中,根据式(5)的选择操作被淘汰的父代解会被添加至存档。若存档规模超过预设阈值(例如NP),则通过随机移除部分解以维持存档容量为NP。显然,当设置存档规模为零时,式(7)即退化为式(6)的特例。

该档案库不仅提供了关于进化方向的信息,还能够有效提升种群的多样性。此外,如后续参数自适应调整公式(10)-(12)所示,我们通过鼓励使用较大的F值来增强种群多样性。因此,尽管算法会偏向潜在最优解(可能为局部极小值)的进化方向,但所提出的DE/current-to-pbest/1算法仍能有效避免陷入局部极小。简而言之,DE/rand/1算法在无特定方向偏好的情况下,仅在相对较小的区域内进行搜索;而结合档案库机制的DE/current-to-pbest/1算法则在更大的搜索区域内开展探索,且其搜索方向会偏向于有潜力的进化方向。
JADE的伪代码如表I所示。接下来介绍F和C R的参数自适应。
在这里插入图片描述

B. Parameter Adaptation

在第g代时,每个个体xi的交叉概率CRi均独立地根据均值为 μ C R μ_{CR} μCR 、标准差为0.1的正态分布生成。
在这里插入图片描述
然后截断为[0,1]。将 S C R S_{CR} SCR表示为第g代所有成功交叉概率 C R i CR_i CRi的集合。均值 μ C R μ_{CR} μCR被初始化为0.5,并在每代结束时按如下方式更新:
在这里插入图片描述
其中c是0和1之间的正常数, m e a n A ( ⋅ ) meanA(·) meanA()是通常的算术平均值。

类似地,在第g代时,每个个体 x i x_i xi的变异因子 F i F_i Fi均独立地根据位置参数 μ F μ_F μF、尺度参数0.1的柯西分布生成。
在这里插入图片描述
F i ≥ 1 F_i≥1 Fi1时将其截断为1,若 F i ≤ 0 F_i≤0 Fi0则重新生成。令 S F S_F SF表示第g代中所有成功变异因子的集合。柯西分布的位置参数 μ F μ_F μF被初始化为0.5,并在每代结束时按如下方式更新:
在这里插入图片描述
其中,meanL(·) 表示莱默平均(Lehmer mean)。
在这里插入图片描述

C. Explanations of the Parameter Adaptation

μ C R μ_{CR} μCR的自适应调整遵循以下核心原则:较优的控制参数值能够产生具有更高生存概率的个体,因此这些参数值应当被传递到后续进化代际中。其核心机制在于记录近期成功的交叉概率值,并利用这些历史信息指导新一代 C R i CR_i CRi的生成。公式(8)中设定较小标准差的原因在于:若标准差过大,参数调整效率将显著降低——例如在标准差趋近无穷大的极端情况下,截断正态分布将与 μ C R μ_{CR} μCR的实际取值完全无关。

C R CR CR的调整相比, μ F μ_F μF的适配包含两项独特操作:首先, F i F_i Fi的生成采用截断柯西分布。与正态分布相比,柯西分布能更有效地实现变异因子的多样化,从而避免在贪婪变异策略(如DE/best、DE/current-to-best和DE/current-to-pbest)中经常出现的早熟收敛现象——这种现象通常发生在变异因子高度集中于某个特定值时。

其次, μ F μ_F μF 的适配过程通过采用Lehmer均值(如公式12所示)而非 μ C R μ_{CR} μCR 调整中使用的算术均值,赋予较大成功变异因子更高权重这种Lehmer均值有助于传播更大的变异因子,从而提升进化进度率。相反, S F S_F SF 的算术均值往往会小于变异因子的最优值,导致 μ F μ_F μF偏小并引发后期早熟收敛。 μ F μ_F μF 偏小的趋势主要源于进化搜索中成功概率与进度率之间的差异:采用较小 F i F_i Fi 的DE/current-to-pbest策略本质上类似于(1+1)进化策略(ES)[23],两者都在基向量的小邻域内生成后代。对于(1+1)ES而言,已有研究表明(在corridor and sphere 函数上严格证明[23][34]),变异方差越小通常成功概率越高。但趋近于0的变异方差显然会导致进化停滞。因此,通过赋予较大成功变异因子更高权重来实现更快进化进度,是一种简洁而高效的解决方案。

关于公式(9)和(11)中的常数c:当c=0时,系统不进行参数自适应调整;反之,成功的 C R i CR_i CRi F i F_i Fi 的"生命周期"约为1/c代——具体而言,经过1/c代后,当c趋近于零时,原有的 μ C R μ_{CR} μCR μ F μ_F μF值会按 ( 1 − c ) 1 / c (1−c)^{1/c} (1c)1/c→1/e≈37%的比例衰减。

D. Discussion of Parameter Settings

经典差分进化算法(DE)包含两个需要用户调整的控制参数F和CR。这些参数具有问题依赖性,因此通常需要繁琐的试错过程才能为特定问题找到合适的取值。与之相反,JADE算法新引入的两个参数c和p预计对不同问题具有不敏感性——这源于它们在JADE中的特定作用:c控制参数自适应速率,p决定变异策略的贪婪程度。如第五节所示,JADE通常在1/c∈[5,20](即 μ C R μ_{CR} μCR μ F μ_F μF值的生命周期为5-20代)且p∈[5%,20%]时表现最优,这意味着变异过程中我们仅考虑前5%-20%的高质量解。

Code

JADE.m

JADE算法。

%**************************************************************************************************
%Reference:  J. Zhang and A. C. Sanderson, "JADE: adaptive differential evolution
%                     with optional external archive," IEEE Trans. Evolut. Comput., vol. 13,
%                     no. 5, pp. 945-958, 2009.
%
% Note: We obtained the MATLAB source code from the authors, and did some
%           minor revisions in order to solve the 25 benchmark test functions,
%           however, the main body was not changed.
%**************************************************************************************************

clc;
clear all;
tic;

format long;
format compact;

'JADE'

% Choose the problems to be tested. Please note that for test functions F7
% and F25, the global optima are out of the initialization range. For these
% two test functions, we do not need to judge whether the variable violates
% the boundaries during the evolution after the initialization.
problemSet = [1 : 6 8 : 24];
for problemIndex = 1 : 23

    problem = problemSet(problemIndex)

    % Define the dimension of the problem
    n = 30;

    popsize = 100;

    switch problem

        case 1

            % lu: define the upper and lower bounds of the variables
            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            % Load the data for this test function
            load sphere_func_data
            A = []; M = []; a = []; alpha = []; b = [];

        case 2

            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load schwefel_102_data
            A = []; M = []; a = []; alpha = []; b = [];

        case 3

            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load high_cond_elliptic_rot_data
            A = []; a = []; alpha = []; b = [];

            if n == 2, load elliptic_M_D2,
            elseif n == 10, load elliptic_M_D10,
            elseif n == 30, load elliptic_M_D30,
            elseif n == 50, load elliptic_M_D50,
            end

        case 4

            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load schwefel_102_data
            A = []; M = []; a = []; alpha = []; b = [];

        case 5

            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load schwefel_206_data
            M = []; a = []; alpha = []; b = [];

        case 6

            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load rosenbrock_func_data
            A = []; M = []; a = []; alpha = []; b = [];

        case 7

            lu = [0 * ones(1, n); 600 * ones(1, n)];
            load griewank_func_data
            A = []; a = []; alpha = []; b = [];

            if n == 2, load griewank_M_D2,
            elseif n == 10, load griewank_M_D10,
            elseif n == 30, load griewank_M_D30,
            elseif n == 50, load griewank_M_D50,
            end

        case 8

            lu = [-32 * ones(1, n); 32 * ones(1, n)];
            load ackley_func_data
            A = []; a = []; alpha = []; b = [];

            if n == 2, load ackley_M_D2,
            elseif n == 10, load ackley_M_D10,
            elseif n == 30, load ackley_M_D30,
            elseif n == 50, load ackley_M_D50,
            end

        case 9

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load rastrigin_func_data
            A = []; M = []; a = []; alpha = []; b = [];

        case 10

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load rastrigin_func_data
            A = []; a = []; alpha = []; b = [];
            if n == 2, load rastrigin_M_D2,
            elseif n == 10, load rastrigin_M_D10,
            elseif n == 30, load rastrigin_M_D30,
            elseif n == 50, load rastrigin_M_D50,
            end

        case 11

            lu = [-0.5 * ones(1, n); 0.5 * ones(1, n)];
            load weierstrass_data
            A = []; a = []; alpha = []; b = [];
            if n == 2, load weierstrass_M_D2, ,
            elseif n == 10, load weierstrass_M_D10,
            elseif n == 30, load weierstrass_M_D30,
            elseif n == 50, load weierstrass_M_D50,
            end

        case 12

            lu = [-pi * ones(1, n); pi * ones(1, n)];
            load schwefel_213_data
            A = []; M = []; o = [];

        case 13

            lu = [-3 * ones(1, n); 1 * ones(1, n)];
            load EF8F2_func_data
            A = []; M = []; a = []; alpha = []; b = [];

        case 14

            lu = [-100 * ones(1, n); 100 * ones(1, n)];
            load E_ScafferF6_func_data
            if n == 2, load E_ScafferF6_M_D2, ,
            elseif n == 10, load E_ScafferF6_M_D10,
            elseif n == 30, load E_ScafferF6_M_D30,
            elseif n == 50, load E_ScafferF6_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 15

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func1_data
            A = []; M = []; a = []; alpha = []; b = [];

        case 16

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func1_data
            if n == 2, load hybrid_func1_M_D2,
            elseif n == 10, load hybrid_func1_M_D10,
            elseif n == 30, load hybrid_func1_M_D30,
            elseif n == 50, load hybrid_func1_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 17

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func1_data
            if n == 2, load hybrid_func1_M_D2,
            elseif n == 10, load hybrid_func1_M_D10,
            elseif n == 30, load hybrid_func1_M_D30,
            elseif n == 50, load hybrid_func1_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 18

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func2_data
            if n == 2, load hybrid_func2_M_D2,
            elseif n == 10, load hybrid_func2_M_D10,
            elseif n == 30, load hybrid_func2_M_D30,
            elseif n == 50, load hybrid_func2_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 19

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func2_data
            if n == 2, load hybrid_func2_M_D2,
            elseif n == 10, load hybrid_func2_M_D10,
            elseif n == 30, load hybrid_func2_M_D30,
            elseif n == 50, load hybrid_func2_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 20

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func2_data
            if n == 2, load hybrid_func2_M_D2,
            elseif n == 10, load hybrid_func2_M_D10,
            elseif n == 30, load hybrid_func2_M_D30,
            elseif n == 50, load hybrid_func2_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 21

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func3_data
            if n == 2, load hybrid_func3_M_D2,
            elseif n == 10, load hybrid_func3_M_D10,
            elseif n == 30, load hybrid_func3_M_D30,
            elseif n == 50, load hybrid_func3_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 22

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func3_data
            if n == 2, load hybrid_func3_HM_D2,
            elseif n == 10, load hybrid_func3_HM_D10,
            elseif n == 30, load hybrid_func3_HM_D30,
            elseif n == 50, load hybrid_func3_HM_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 23

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func3_data
            if n == 2, load hybrid_func3_M_D2,
            elseif n == 10, load hybrid_func3_M_D10,
            elseif n == 30, load hybrid_func3_M_D30,
            elseif n == 50, load hybrid_func3_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 24

            lu = [-5 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func4_data
            if n == 2, load hybrid_func4_M_D2,
            elseif n == 10, load hybrid_func4_M_D10,
            elseif n == 30, load hybrid_func4_M_D30,
            elseif n == 50, load hybrid_func4_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

        case 25

            lu = [2 * ones(1, n); 5 * ones(1, n)];
            load hybrid_func4_data
            if n == 2, load hybrid_func4_M_D2,
            elseif n == 10, load hybrid_func4_M_D10,
            elseif n == 30, load hybrid_func4_M_D30,
            elseif n == 50, load hybrid_func4_M_D50,
            end
            A = []; a = []; alpha = []; b = [];

    end

    % Record the best results
    outcome = [];

    %Main body which was developed by the authors

    time = 1; % 当前运行次数

    % The total number of runs
    totalTime = 1;

    while time <= totalTime

        rand('seed', sum(100 * clock)); % 设置随机种子

        % Initialize the main population
        popold = repmat(lu(1, :), popsize, 1) + rand(popsize, n) .* (repmat(lu(2, :) - lu(1, :), popsize, 1));

        valParents = benchmark_func(popold, problem, o, A, M, a, alpha, b); % 计算适应度值 NP*1

        c = 1/10;
        p = 0.05;

        CRm = 0.5;  % 初始化CR均值
        Fm = 0.5; % 初始化F均值

        Afactor = 1;  % 用于控制外部存档和种群大小的比例

        archive.NP = Afactor * popsize; % the maximum size of the archive 外部存档的大小
        archive.pop = zeros(0, n); % the solutions stored in te archive 劣解
        archive.funvalues = zeros(0, 1); % the function value of the archived solutions 劣解的适应度值

        %% the values and indices of the best solutions 升序排序,并得到对应的索引,最好的值在最前面
        [valBest, indBest] = sort(valParents, 'ascend');

        FES = 0;  % 函数评估次数
        while FES < n * 10000 %& min(fit)>error_value(problem)

            pop = popold; % the old population becomes the current population

            if FES > 1 && ~isempty(goodCR) && sum(goodF) > 0 % If goodF and goodCR are empty, pause the update
                CRm = (1 - c) * CRm + c * mean(goodCR);
                Fm = (1 - c) * Fm + c * sum(goodF .^ 2) / sum(goodF); % Lehmer mean
            end

            % Generate CR according to a normal distribution with mean CRm, and std 0.1
            % Generate F according to a cauchy distribution with location parameter Fm, and scale parameter 0.1
            [F, CR] = randFCR(popsize, CRm, 0.1, Fm, 0.1);  % 每个个体对应一个F和CR

            r0 = [1 : popsize];
            popAll = [pop; archive.pop];  % PUA
            [r1, r2] = gnR1R2(popsize, size(popAll, 1), r0);  % 生成参数r1和r2

            % Find the p-best solutions
            pNP = max(round(p * popsize), 2); % choose at least two best solutions
            randindex = ceil(rand(1, popsize) * pNP); % select from [1, 2, 3, ..., pNP]
            randindex = max(1, randindex); % to avoid the problem that rand = 0 and thus ceil(rand) = 0
            pbest = pop(indBest(randindex), :); % randomly choose one of the top 100p% solutions

            % == == == == == == == == == == == == == == == Mutation == == == == == == == == == == == == ==
            vi = pop + F(:, ones(1, n)) .* (pbest - pop + pop(r1, :) - popAll(r2, :));

            vi = boundConstraint(vi, pop, lu);  % 越界处理

            % == == == == = Crossover == == == == =
            mask = rand(popsize, n) > CR(:, ones(1, n)); % mask is used to indicate which elements of ui comes from the parent
            % choose one position where the element of ui doesn't come from the parent
            rows = (1 : popsize)';cols = floor(rand(popsize, 1) * n)+1;
            jrand = sub2ind([popsize n], rows, cols); % 使用 sub2ind 将子脚标转换为线性索引
            mask(jrand) = false;  % 确保每个个体至少有一个维度来自突变向量
            ui = vi; ui(mask) = pop(mask);

            valOffspring = benchmark_func(ui, problem, o, A, M, a, alpha, b); % 计算适应度值

            FES = FES + popsize;  % 函数评估次数+popsize

            % == == == == == == == == == == == == == == == Selection == == == == == == == == == == == == ==
            % I == 1: the parent is better; I == 2: the offspring is better
            [valParents, I] = min([valParents, valOffspring], [], 2);
            popold = pop;

            archive = updateArchive(archive, popold(I == 2, :), valParents(I == 2));  % 更新外部存档

            popold(I == 2, :) = ui(I == 2, :);

            goodCR = CR(I == 2);
            goodF = F(I == 2);

            [valBest indBest] = sort(valParents, 'ascend');

        end

        outcome = [outcome min(valParents)];  % 记录本次运行的最优值

        time = time + 1;

    end

    sort(outcome)
    mean(outcome)
    std(outcome)

end
toc;

randFCR.m

为种群中的每个个体生成F和CR。

function [F,CR] = randFCR(NP, CRm, CRsigma, Fm,  Fsigma)

% this function generate CR according to a normal distribution with mean "CRm" and sigma "CRsigma"
%           If CR > 1, set CR = 1. If CR < 0, set CR = 0.
% this function generate F  according to a cauchy distribution with location parameter "Fm" and scale parameter "Fsigma"
%           If F > 1, set F = 1. If F <= 0, regenrate F.
%
% Version: 1.1   Date: 11/20/2007
% Written by Jingqiao Zhang (jingqiao@gmail.com)

%% generate CR
CR = CRm + CRsigma * randn(NP, 1);
CR = min(1, max(0, CR));                % truncated to [0 1]

%% generate F
F = randCauchy(NP, 1, Fm, Fsigma);
F = min(1, F);                          % truncation

% we don't want F = 0. So, if F<=0, we regenerate F (instead of trucating it to 0)
pos = find(F <= 0);
while ~ isempty(pos)
    F(pos) = randCauchy(length(pos), 1, Fm, Fsigma);
    F = min(1, F);                      % truncation
    pos = find(F <= 0);
end

% Cauchy distribution: cauchypdf = @(x, mu, delta) 1/pi*delta./((x-mu).^2+delta^2)
function result = randCauchy(m, n, mu, delta)

% http://en.wikipedia.org/wiki/Cauchy_distribution
result = mu + delta * tan(pi * (rand(m, n) - 0.5));

updateArchive.m

更新存档并输入解决方案
步骤1:将新解决方案添加至存档
步骤2:移除重复项
步骤3:如有需要,随机删除部分解决方案以维持存档容量

function archive = updateArchive(archive, pop, funvalue)
% Update the archive with input solutions
%   Step 1: Add new solution to the archive
%   Step 2: Remove duplicate elements
%   Step 3: If necessary, randomly remove some solutions to maintain the archive size
%
% Version: 1.1   Date: 2008/04/02
% Written by Jingqiao Zhang (jingqiao@gmail.com)

if archive.NP == 0, return; end

if size(pop, 1) ~= size(funvalue,1), error('check it'); end

% Method 2: Remove duplicate elements
popAll = [archive.pop; pop ];
funvalues = [archive.funvalues; funvalue ];
[dummy IX]= unique(popAll, 'rows');
if length(IX) < size(popAll, 1) % There exist some duplicate solutions
  popAll = popAll(IX, :);
  funvalues = funvalues(IX, :);
end

if size(popAll, 1) <= archive.NP   % add all new individuals
  archive.pop = popAll;
  archive.funvalues = funvalues;
else                % randomly remove some solutions
  rndpos = randperm(size(popAll, 1)); % equivelent to "randperm";
  rndpos = rndpos(1 : archive.NP);
  
  archive.pop = popAll  (rndpos, :);
  archive.funvalues = funvalues(rndpos, :);
end

gnR1R2.m

选择两个随机个体r1和r2。

function [r1, r2] = gnR1R2(NP1, NP2, r0)

% gnA1A2 generate two column vectors r1 and r2 of size NP1 & NP2, respectively
%    r1's elements are choosen from {1, 2, ..., NP1} & r1(i) ~= r0(i)
%    r2's elements are choosen from {1, 2, ..., NP2} & r2(i) ~= r1(i) & r2(i) ~= r0(i)
%
% Call:
%    [r1 r2 ...] = gnA1A2(NP1)   % r0 is set to be (1:NP1)'
%    [r1 r2 ...] = gnA1A2(NP1, r0) % r0 should be of length NP1
%
% Version: 2.1  Date: 2008/07/01
% Written by Jingqiao Zhang (jingqiao@gmail.com)

NP0 = length(r0);

r1 = floor(rand(1, NP0) * NP1) + 1;
for i = 1 : inf
    pos = (r1 == r0);
    if sum(pos) == 0
        break;
    else % regenerate r1 if it is equal to r0
        r1(pos) = floor(rand(1, sum(pos)) * NP1) + 1;
    end
    if i > 1000, % this has never happened so far
        error('Can not genrate r1 in 1000 iterations');
    end
end

r2 = floor(rand(1, NP0) * NP2) + 1;
for i = 1 : inf
    pos = ((r2 == r1) | (r2 == r0));
    if sum(pos)==0
        break;
    else % regenerate r2 if it is equal to r0 or r1
        r2(pos) = floor(rand(1, sum(pos)) * NP2) + 1;
    end
    if i > 1000, % this has never happened so far
        error('Can not genrate r2 in 1000 iterations');
    end
end

boundConstraint.m

越界处理。

function vi = boundConstraint (vi, pop, lu)

% if the boundary constraint is violated, set the value to be the middle
% of the previous value and the bound
%
% Version: 1.1   Date: 11/20/2007
% Written by Jingqiao Zhang, jingqiao@gmail.com

[NP, D] = size(pop);  % the population size and the problem's dimension

%% check the lower bound
xl = repmat(lu(1, :), NP, 1);
pos = vi < xl;  % NP*D找到所有违反约束下界的值的位置
vi(pos) = (pop(pos) + xl(pos)) / 2;

%% check the upper bound
xu = repmat(lu(2, :), NP, 1);
pos = vi > xu;  % NP*D找到所有违反约束上界的值的位置
vi(pos) = (pop(pos) + xu(pos)) / 2;

参考

[1] Jingqiao Zhang, Sanderson A C. JADE: Adaptive differential evolution with optional external archive[J]. IEEE Transactions on Evolutionary Computation, 2009, 13(5): 945-958.

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

W_chuanqi

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值