差分进化算法的变体之一——JADE算法

对于DE算法而言，随着迭代次数的增加，个体间的差异会逐渐降低，收敛速度也会随之下降，这会使得DE算法容易陷入局部最优和早熟收敛。所以很多研究者在原始经典的DE算法上寻求各种改进来提高DE算法的寻优能力、收敛速度、克服早熟收敛等。

DE算法主要涉及种群规模NP、缩放因子F、交叉概率CR这三个控制参数。原始的经典DE算法通常都是根据经验来选取一组固定的参数大小：NP∈[5D，10D]；F通常取0.5；CR∈[0，1]，通常取0.3，一般可以保证较高的寻优成功率和较快的收敛速度。而现在有很多基于控制参数改进的DE变体算法，大体可以分为适应性DE(如：jDE，JADE，SHADE等)和自适应DE(如：SPDE，DESAP，SELSDE等)。两者的共同点是，控制参数在迭代过程都会变化；不同点是，自适应DE的控制参数与种群中的个体有关且在迭代过程中会经过类似突变、交叉等操作来达到自适应的目的。

之前看了一篇关于JADE算法的文献——《JADE: Adaptive Differential Evolution with Optional External Archive》，其中讲了JADE算法的实现过程，这里相较于DE的改进主要有三点：a.用到的变异策略是Current-to-pbest；b.归档种群与不归档；c.适应性参数调整CR和F。

下面给出我自己实现的JADE的matlab代码，其中柯西分布的相关函数文件夹(CauchyFunction)是网上下载，代码分为三部分，JADE_demo.m是JADE算法的函数文件，TestFunction.m是benchmark中的13类标准测试函数组成的函数文件，JADE_RunTest.m是为了方便最后统计JADE算法run不同测试函数时的均值和方差。

JADE_demo.m：

function optValue=JADE_demo(ma,mi,gen,type)

    addpath('../CauchyFunction');

    %%初始化
    Pmax=ma;%搜索上界
    Pmin=mi;%搜索下界
    NP=100;%种群规模
    Dim=30;%个体维数
    G=0;%设置迭代器（当前迭代代数）
    uCR=0.5;%初始化交叉概率
    uF=0.5;%初始化缩放因子
    p0=0.05;
    top=p0*NP;%每代中最优的top个
    A=[];%初始归档种群为空集
    t=1;%记录归档种群A中的个体个数
    Gmax=gen;%最大迭代次数
    index=type;%测试函数类型
    
%     index=input('Input the index of TestFunction:');

    P=(Pmax-Pmin)*rand(NP,Dim)+Pmin;%随机产生初始种群个体

    %%总体大循环
    while G<Gmax
        Scr=[];%初始成功参加变异的个体的交叉概率为空集
        Sf=[];%初始成功参加变异的个体的缩放因子为空集
        n1=1;%记录Scr中的元素个数
        n2=1;%记录Sf中的元素个数

        for i=1:NP
            fitnessP(i)=TestFunction(index,P(i,:));%计算种群个体适应值
            %对第G代的每个个体计算交叉概率和缩放因子
            CR(i)=normrnd(uCR,0.1);%服从正态分布
            F(i)=cauchyrnd(uF,0.1);%服从柯西分布
            while (CR(i)>1||CR(i)<0)
                CR(i)=normrnd(uCR,0.1);
            end
            if (F(i)>1)
                F(i)=1;
            end
            while (F(i)<=0)
                F(i)=cauchyrnd(uF,0.1);
            end
        end
        [fitnessBestP,indexBestP]=min(fitnessP);
        bestP=P(indexBestP,:);

        %根据个体适应值进行排序，得到最好的前top个个体
        [~,indexSortP]=sort(fitnessP);
        for i=1:top
           bestTopP(i,:)=P(indexSortP(i),:); 
        end

        %%变异操作
        for i=1:NP   
            %从top个个体中随机选出一个作为Xpbest
            k0=randperm(top,1);
            Xpbest=bestTopP(k0,:);

            %从当前种群P中随机选出P1
            k1=randi(NP);
            P1=P(k1,:);
            while (k1==i||k1==k0)
                k1=randi(NP);
                P1=P(k1,:); 
            end

            %从P∪A中随机选出P2
            PandA=[P;A];
            [num,~]=size(PandA);
            k2=randi(num);
            P2=PandA(k2,:);
            while (k2==i||k2==k0||k2==k1)
                k2=randi(num);
                P2=PandA(k2,:); 
            end

            V(i,:)=P(i,:)+F(i).*(Xpbest-P(i,:))+F(i).*(P1-P2);   
        end

        %%交叉操作
        for i=1:NP
            jrand=randi([1,Dim]); 
            for j=1:Dim
                k3=rand;
                if(k3<=CR(i)||j==jrand)
                    U(i,j)=V(i,j);
                else
                    U(i,j)=P(i,j);
                end
            end
        end

        %%边界处理
        for i=1:NP
           for j=1:Dim
              while (U(i,j)>Pmax||U(i,j)<Pmin)
                 U(i,j)=(Pmax-Pmin)*rand+Pmin; 
              end
           end
        end

        %%选择操作
        for i=1:NP
            fitnessU(i)=TestFunction(index,U(i,:));
            if(fitnessU(i)<fitnessP(i))
                A(t,:)=P(i,:);
                P(i,:)=U(i,:);
                fitnessP(i)=fitnessU(i);
                Scr(n1)=CR(i);
                Sf(n2)=F(i);
                t=t+1;
                n1=n1+1;
                n2=n2+1;
                if(fitnessU(i)<fitnessBestP)
                   fitnessBestP=fitnessU(i);
                   bestP=U(i,:);
                end
            end
        end

        %判断归档种群A的规模是否在NP之内，若大于，则随机移除个体使其规模保持NP
        [tA,~]=size(A);
        if tA>NP
            nRem=tA-NP;
            k4=randperm(tA,nRem);
            A(k4,:)=[]; 
            [tA,~]=size(A);
            t=tA+1;
        end

        %自适应参数，更新uCR和uF
        c=0.1;
        [~,ab]=size(Scr);
        if ab~=0
            newSf=(sum(Sf.^2))/(sum(Sf));
            uCR=(1-c)*uCR+c.*mean(Scr);
            uF=(1-c)*uF+c.*newSf;
        end

        G=G+1;
        bestfitnessG(G)=fitnessBestP;
    end

    optValue=fitnessBestP;
    
    %%画图
    %plot(bestfitnessG);
%     optValue=num2str(fitnessBestP);
%     Location=num2str(bestP);
%     disp(strcat('the optimal value','=',optValue));
%     disp(strcat('the best location','=',Location));

    %%测试函数
    % function y=testFun(x)
    %     y=0;
    %     for i=2:length(x)
    %         y=y+100*(x(i)-x(i-1)^2)^2+(x(i-1)-1)^2;
    %     end
    % end

end

TestFunction.m：

function [ y ] = TestFunction( type ,x)
%TESTFUNCTION 此处显示有关此函数的摘要
%   此处显示详细说明
    if type==1
        y=sum(x.^2);
    elseif type==2
        y=sum(abs(x))+prod(abs(x));
    elseif type==3
        len = length(x);
        y = 0;
        for i=1:len
            y = y+sum(x(1:i)).^2;
        end
    elseif type==4
        y = max(abs(x));
    elseif type==5
        dimension = length(x);
        y = sum(100*(x(2:dimension)-x(1:dimension-1).^2).^2 + (1-x(1:dimension-1)).^2);
    elseif type==6
        y = sum((floor(x+0.5)).^2);
    elseif type==7
        len = length(x);
        y = 0;
        for i=1:len
           y = y+i*x(i).^4;
        end

        y = y+rand();
    elseif type==8
        D = length(x);
        y = sum(-x.*sin(sqrt(abs(x))))+418.9828872724399*D;
    elseif type==9
        y = sum(x.^2-10*cos(2*pi.*x)+10);
    elseif type==10
         D = length(x);
         y = -20*exp(-0.2*sqrt(1/D*sum(x.^2)))-exp(1/D*sum(cos(2*pi.*x)))+20+exp(1);
    elseif type==11
        len = length(x);
        y = 1/4000*sum(x.^2)-prod(cos(x(1:len)/sqrt(1:len)))+1;
    elseif type==12
        D = length(x);
        Y=funcY(x);
        u=funcU(x,10,100,4);
    
        y=1/D*pi.*( 10*(sin(pi*Y(1))).^2+sum(((Y(1:D-1)-1).^2).*(1+10.*(sin(pi.*Y(2:D))).^2))+(Y(D)-1).^2 )+sum(u);
    elseif type==13
        D = length(x);
        u = funcU(x,5,100,4);
        y = 0.1.*((sin(3*pi*x(1))).^2 +sum((x(1:D-1)-1).^2.*(1+(sin(3*pi.*(x(2:D)))).^2 )))+0.1*((x(D)-1).^2)*(1+(sin(2*pi*x(D))).^2)+sum(u);
    end
        

end

function [ y ] = funcU( x,a,k,m )
%FUNCU 此处显示有关此函数的摘要
%   此处显示详细说明 a=10 or 5, k=100,m =4
    len = length(x);
    for i=1:len
        if x(i)>a
            y(i)=k*(x(i)-a).^m;
        else if -a<=x(i)&&x(i)<=a
                y(i)=0;
            else
                y(i)=k*(-x(i)-a).^m;
            end
        end
    end
end

function [ y ] = funcY( x )
%FUNCY 此处显示有关此函数的摘要
%   此处显示详细说明
    y = 1+(x+1)/4;

end

JADE_RunTest.m：

clear all

format long;
format compact;

% max=100;%搜索上界
% min=-100;%搜索下界
runtime=50;
% type=13;

tic()
% for i=1:type
    max=input('Input the upper bound:');
    min=input('Input the lower bound:');
    gMax=input('Input the Gen:');
    indexFun=input('Input the type of TestFuntion:');
    for j=1:runtime
        bestValue(j)=JADE_demo(max,min,gMax,indexFun);
    end
    mean_bestValue=mean(bestValue);
    std_bestValue=std(bestValue);
    disp(strcat('第',num2str(indexFun),'个函数结果为：mean=',num2str(mean_bestValue),' , std=',num2str(std_bestValue)));
% end
toc()