image thumbnail

Hybrid Differential Evolution Algorithm With Adaptive Crossover Mechanism

by

 

28 Nov 2012 (Updated )

An EA based on DE with adaptive crossover rate, population refresh and local search.

DE_TCR.m
%% DE_TCRparam
% Runs the DE_TCR optimization algorithm.
%%
%% Beta version 
% Copyright 2006 - 2012 - CPOH  
%
% Predictive Control and Heuristic Optimization Research Group
%      http://cpoh.upv.es
%
% ai2 Institute
%      http://www.ai2.upv.es
%
% Universitat Politcnica de Valncia - Spain.
%      http://www.upv.es
%
%%
%% Author
% Gilberto Reynoso Meza
% gilreyme@upv.es
% http://cpoh.upv.es/en/gilberto-reynoso-meza.html
% http://www.mathworks.es/matlabcentral/fileexchange/authors/289050
%%
%% For new releases and bug fixing of this Tool Set please visit:
% 1) http://cpoh.upv.es/en/research/software.html
% 2) Matlab Central File Exchange
%%
%% Overall Description
% This code implements a modified version of the single-objective 
% optimization algorithm described at:
%
% G. Reynoso; J. Sanchis; X. Blasco; Juan M. Herrero. Hybrid DE Algorithm 
% With Adaptive Crossover Operator For Solving Real-World Numerical 
% Optimization Problems. In IEEE Congress on Evolutionary Computation. 
% CEC 2011. (ISBN 978-1-4244-7833-0). New Orleans (USA). June 2011.
%
%%

function OUT=DE_TCR(DE_TCRDat)

%% Reading parameters from DE_TCRDat

% Problem to solve

Generations  = DE_TCRDat.MAXGEN;    % Generations.
Nvar         = DE_TCRDat.NVAR;      % Number of decision variables.
Bounds       = DE_TCRDat.FieldD;    % Optimization Bounds.
Initial      = DE_TCRDat.Initial;   % (Re)Initialization Bounds.
sop          = DE_TCRDat.sop;       % Cost function.

% Local Search

RateLS       = DE_TCRDat.RateLS;    % Probability for Local Search.

% Adaptive Mechanism

MaxCr = DE_TCRDat.TCR(1,3);         % Maximum Crossover rate.
MedCr = DE_TCRDat.TCR(1,2);         % Median Crossover rate.
MinCr = DE_TCRDat.TCR(1,1);         % Minimum Crossover rate.
MaxF  = DE_TCRDat.TF(1,3);          % Maximum Scaling Factor.
MedF  = DE_TCRDat.TF(1,2);          % Median Scaling Factor.
MinF  = DE_TCRDat.TF(1,1);          % Minimum Scaling Factor.

recombination = DE_TCRDat.Recomb;   % Linear recombination factor.
CRsuccess = DE_TCRDat.CRsuccess;    % Threshold for CR successes.

% Population Management

Xpop        = DE_TCRDat.Xpop;       % Population Size.
GammaVar    = DE_TCRDat.GammaVar;   % Interquartil difference for refresh.
minVarPop   = DE_TCRDat.minVarPop;  % Minimum diversity in population.
XpopRefresh = DE_TCRDat.XpopRefresh;% Population Refresh size.

%% Initialization of variables

% Adaptive Mechanism

ParamEVO = ...
    zeros(DE_TCRDat.CRsuccess,1);   % To record successful values on CR.
SuccessEvoCr = 0;                   % To recor successes on CR.

% Population

Child=zeros(Xpop,Nvar);             % Child population
Parent=zeros(Xpop,Nvar);            % Parent Population
Mutant=zeros(Xpop,Nvar);            % Mutant Population
JxParent=zeros(Xpop,1);             % Fitness of Parento Population.

%% Population for evolutionary process

PopParent=zeros(Xpop,Nvar);         % Population initialization

for xpop=1:Xpop                     % Uniform distribution
    for nvar=1:Nvar
        PopParent(xpop,nvar) = Initial(nvar,1) + ...
            (Initial(nvar,2) - Initial(nvar,1))*rand();
    end
end

JxPopParent = sop(PopParent,DE_TCRDat);
FES = Xpop;

%% Evolution process

for n=1:Generations
    
    DE_TCRDat.CounterGEN=n;
    OUT.JxPopParent=JxPopParent;
    OUT.JxParent=JxParent;
    OUT.Param=DE_TCRDat;
    OUT.TCR=[MinCr MedCr MaxCr]; 
    
    [OUT DE_TCRDat]=PrinterDisplay(OUT,DE_TCRDat); % To display results.

    
    %% Population Refreshment Mechanism
    VarPop=iqr(PopParent); % Measuring the interquartile range difference.
    
    if VarPop<=((Initial(:,2)-Initial(:,1))')/GammaVar %We need a refresh.
        LookingBest=sortrows([JxPopParent PopParent]);
        if VarPop<minVarPop
            Initial(:,1) = (median(PopParent))' - minVarPop';
            Initial(:,2) = (median(PopParent))' + minVarPop';
        else
            Initial(:,1) = (median(PopParent))' - ...
                (Initial(:,2)-Initial(:,1))/GammaVar;
            Initial(:,2) = (median(PopParent))' + ...
                (Initial(:,2)-Initial(:,1))/GammaVar;
        end
        for nvar=1:Nvar % Checking bounds
            if Initial(nvar,1) < Bounds(nvar,1)
                Initial(nvar,1) = Bounds(nvar,1);
            end
            if Initial(nvar,2) > Bounds(nvar,2)
                Initial(nvar,2) = Bounds(nvar,2);
            end
        end
        
        VVrefresh=zeros(Xpop-DE_TCRDat.XpopRefresh,Nvar);
        
        for xpop=1:Xpop-DE_TCRDat.XpopRefresh % New individuals.
            for nvar=1:Nvar
                VVrefresh(xpop,nvar) = Initial(nvar,1) + ...
                    (Initial(nvar,2) - Initial(nvar,1))*rand();
            end
        end
        
        FES=FES+Xpop;
        
        PopParent=[VVrefresh;
            LookingBest(1:XpopRefresh,2:end)];
        %We keep the best individuals + New individuals
        
        JxPopParent = sop(PopParent,DE_TCRDat); %Evaluate New Individuals.
    end
    
    ParamEvo=zeros(Xpop,4); % Parameter matrix initialization.
    
    %% Adaptive Mechanism
    
    for xpop=1:Xpop % Calculating DE parameters for each individual.
        % Scaling Factor.
        randU=rand;
        if randU<(MedF-MinF)*(2/(MaxF-MinF))*0.5
            ParamEvo(xpop,1) = ...
                MinF+sqrt(randU*(MaxF-MinF)*(MedF-MinF));
        else
            ParamEvo(xpop,1) = ...
                MaxF-sqrt((1-randU)*(MaxF-MinF)*(MaxF-MedF));
        end
        
        % CrossOver Rate.
        randU=rand;
        if randU<(MedCr-MinCr)*(2/(MaxCr-MinCr))*0.5
            ParamEvo(xpop,2) = ...
                MinCr+sqrt(randU*(MaxCr-MinCr)*(MedCr-MinCr));
        else
            ParamEvo(xpop,2) = ...
                MaxCr-sqrt((1-randU)*(MaxCr-MinCr)*(MaxCr-MedCr));
        end
    end
    
    Parent=PopParent;
    JxParent=JxPopParent;
    
    %% DE ALGORITHM
    
    for xpop=1:Xpop
        
        rev=randperm(Xpop);
        ScalingFactor=ParamEvo(xpop,1);
        
        % Mutant vector calculation.
        Mutant(xpop,:)= Parent(rev(1,1),:)+ScalingFactor*...
            (Parent(rev(1,2),:)-Parent(rev(1,3),:));
        
        for nvar=1:Nvar % Bounds are always verified.
            if Mutant(xpop,nvar)<Bounds(nvar,1)
                Mutant(xpop,nvar) = Bounds(nvar,1);
            elseif Mutant(xpop,nvar)>Bounds(nvar,2)
                Mutant(xpop,nvar)=Bounds(nvar,2);
            end
        end
        
        % Child calculation.
        
        FlagCr=0;  % Counter for Child=Mutant?
        FlagCr2=0; % Counter for Child=Paretn?
        CRrate=ParamEvo(xpop,2);
        if Nvar > 1
            for nvar=1:Nvar
                if CRrate>=0.95 % We try a Lineal Recombination.
                    Child(xpop,:)=Parent(xpop,:) + ...
                        recombination.*(Mutant(xpop,:)-Parent(xpop,:));
                    FlagCr=1;
                else % We try a binomial recombination.
                    if rand() > CRrate
                        Child(xpop,nvar) = Parent(xpop,nvar);
                        FlagCr=1;
                        FlagCr2=FlagCr2+1;
                    else
                        Child(xpop,nvar) = Mutant(xpop,nvar);
                    end
                end
            end
        end
        
        
        % If Child==Mutant, with a low CRrate,
        % we keep at least one decision variable from the Parent.
        if FlagCr==0 && CRrate<0.5
            j=randperm(Nvar);
            Child(xpop,j(1,1))=Parent(xpop,j(1,1));
        end
        
        % If Child==Mutant, with a high CRrate,
        % we force a lineal recombination.
        if FlagCr==0 && CRrate>=0.5 && CRrate<0.95
            Child(xpop,:)=Parent(xpop,:) + ...
                recombination.*(Mutant(xpop,:)-Parent(xpop,:));
        end
        
        % if Child==Parent,
        % we keep at least one decision variable from the Mutant.
        if FlagCr2==Nvar
            j=randperm(Nvar);
            Child(xpop,j(1,1))=Mutant(xpop,j(1,1));
        end
        
        for nvar=1:Nvar % Bounds are always verified.
            if Child(xpop,nvar) < Bounds(nvar,1)
                Child(xpop,nvar) = Bounds(nvar,1);
            elseif Child(xpop,nvar) > Bounds(nvar,2)
                Child(xpop,nvar) = Bounds(nvar,2);
            end
        end
        
        %% Local Search
        
        if rand>RateLS % We decide to perform a Local Search.
            [Child(xpop,:),~,FVAL] = ...
                DE_TCRDat.sopLS(Child(xpop,:),DE_TCRDat);
            FES=FES+FVAL;
            if FES>DE_TCRDat.MAXFUNEVALS || n>DE_TCRDat.MAXGEN
                disp('Optimization terminated.')
                break;
            end
        end
        
    end
    %% Selection process
    
    JxChild = sop(Child,DE_TCRDat);
    FES=FES+Xpop;
    
    if FES>DE_TCRDat.MAXFUNEVALS || n>DE_TCRDat.MAXGEN
        disp('Optimization terminated.')
        break;
    end
    
    for xpop=1:Xpop
        if JxChild(xpop,:) <= JxParent(xpop,:)
            ParamEvo(xpop,4)=1;
            SuccessEvoCr=SuccessEvoCr+1;
            ParamEVO(SuccessEvoCr,1) = ...
                ParamEvo(xpop,2); % We keep record of a success.
            Parent(xpop,:) = Child(xpop,:);
            JxParent(xpop,:) = JxChild(xpop,:);
        end
    end
        
    OUT.JxPopParent=JxPopParent;
    OUT.JxParent=JxParent;
    OUT.Param=DE_TCRDat;
    PopParent=Parent;
    JxPopParent=JxParent;
    
    DE_TCRDat.CounterFES=FES;
    
    [OUT DE_TCRDat]=PrinterDisplay(OUT,DE_TCRDat);
    
    if FES>DE_TCRDat.MAXFUNEVALS || n>DE_TCRDat.MAXGEN
        disp('Optimization terminated.')
        break;
    end
    
    %% Adaptive Mechanism: Adapting Triangular distribution for CR rate.
    
    if SuccessEvoCr>=CRsuccess
        SuccessEvoCr=0;
        DataCR=sortrows(ParamEVO,1);
        
        MaxCr=DataCR(end,1);
        MinCr=DataCR(1,1);
        MedCr=median(DataCR);
        
        if abs(MedCr-MaxCr)<0.1
            MaxCr=min(MedCr+0.1,1.0);
        end
        
        if abs(MedCr-MinCr)<0.1
            MinCr=max(MedCr-0.1,0.2);
        end
    end
end

PopParent=Parent;
JxPopParent=JxParent;
OUT.PopParent=PopParent;
OUT.JxPopParent=JxPopParent;
OUT.Param=DE_TCRDat;

if strcmp(DE_TCRDat.SaveResults,'yes')
    save(['OUT_DE_TCR_' datestr(now,30)],'OUT'); % Results are saved.
end

disp('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
if strcmp(DE_TCRDat.SaveResults,'yes')
    disp(['Check out OUT_DE_TCR' datestr(now,30) ...
        ' variable on folder for results.'])
end
disp('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')


function [OUT Dat]=PrinterDisplay(OUT,Dat)

if strcmp(Dat.SeeProgress,'yes')
    disp('------------------------------------------------')
    disp(['Generation: ' num2str(Dat.CounterGEN)]);
    disp(['Evaluations: ' num2str(Dat.CounterFES)]);
    disp(['Minimum: ' num2str(min(OUT.JxPopParent))]);
    disp(num2str(OUT.TCR))
    disp('------------------------------------------------')
end

if strcmp(Dat.Plotter,'yes')
    figure(123);
    plot(Dat.CounterGEN,(min(OUT.JxPopParent(:,1))),'*r'); ...
    grid on; hold on;
end

%%
%% Release and bug report:
%
% November 2012: Initial release
% 19th. April 2013: Bug in Re-Initialization ranges was fixed

Contact us