Code covered by the BSD License

# Multi-Objective Optimization Differential Evolution Algorithm

### Gilberto Reynoso-Meza (view profile)

Runs a basic Differential Evolution (DE) algorithm for multi-objective optimization.

MODE.m
```%% MODE
% Multi-objective Evolutionary Algorithm (MOEA) based on Differential
% Evolution (DE).
% It implements a greedy selection based on pure dominance.
% DE algorithm has been introduced in:
%
% Storn, R., Price, K., 1997. Differential evolution: A simple and
% efficient heuristic for global optimization over continuous spaces.
% Journal of Global Optimization 11, 341  359.
%%
%% 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:
% http://cpoh.upv.es/en/research/software.html
% Matlab Central File Exchange
%%
%% Overall Description
% This code implements a basic multi-objective optimization algorithm based
% on Diferential Evolution (DE) algorithm.
%
% When one objective is optimized, the standard DE runs; if two or more
% objectives are optimized, the greedy selection step in DE algorithm is
% performed using a dominance relation.
%%
%
%%

function OUT=MODE(MODEDat)

Generaciones  = MODEDat.MAXGEN;    % Maximum number of generations.
Xpop          = MODEDat.XPOP;      % Population size.
Nvar          = MODEDat.NVAR;      % Number of decision variables.
Nobj          = MODEDat.NOBJ;      % Number of objectives.
Bounds        = MODEDat.FieldD;    % Optimization bounds.
Initial       = MODEDat.Initial;   % Initialization bounds.
ScalingFactor = MODEDat.Esc;       % Scaling fator in DE algorithm.
CrossOverP    = MODEDat.Pm;        % Crossover probability in DE algorithm.
mop           = MODEDat.mop;       % Cost function.

%% Initial random population
Parent = zeros(Xpop,Nvar);  % Parent population.
Mutant = zeros(Xpop,Nvar);  % Mutant population.
Child  = zeros(Xpop,Nvar);  % Child population.
FES    = 0;                 % Function Evaluation.

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

if size(MODEDat.InitialPop,1)>=1
Parent(1:size(MODEDat.InitialPop,1),:)=MODEDat.InitialPop;
end

JxParent = mop(Parent,MODEDat);
FES = FES+Xpop;

%% Evolution process

for n=1:Generaciones

for xpop=1:Xpop
rev=randperm(Xpop);

%% 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,1);
end
end

%% Crossover operator
for nvar=1:Nvar
if rand() > CrossOverP
Child(xpop,nvar) = Parent(xpop,nvar);
else
Child(xpop,nvar) = Mutant(xpop,nvar);
end
end

end

JxChild = mop(Child,MODEDat);
FES=FES+Xpop;

%% Selection
for xpop=1:Xpop
if JxChild(xpop,:) <= JxParent(xpop,:)
Parent(xpop,:) = Child(xpop,:);
JxParent(xpop,:) = JxChild(xpop,:);
end
end

PFront=JxParent;
PSet=Parent;

OUT.Xpop           = Parent;   % Population
OUT.Jpop           = JxParent; % Poopulation's Objective Vector
OUT.PSet           = PSet;     % Pareto Set
OUT.PFront         = PFront;   % Pareto Front
OUT.Param          = MODEDat;  % MODE Parameters
MODEDat.CounterGEN = n;
MODEDat.CounterFES = FES;

[OUT MODEDat]=PrinterDisplay(OUT,MODEDat); % To print results on screen

if FES>MODEDat.MAXFUNEVALS || n>MODEDat.MAXGEN
disp('Termination criteria reached.')
break;
end
end

OUT.Xpop=PSet;
OUT.Jpop=PFront;
[OUT.PFront, OUT.PSet]=DominanceFilter(PFront,PSet); %A Dominance Filter

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

disp('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')
disp('Red  asterisks : Set Calculated.')
disp('Black diamonds : Filtered Set.')
if strcmp(MODEDat.SaveResults,'yes')
disp(['Check out OUT_' datestr(now,30) ...
' variable on folder for results.'])
end
disp('+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++')

F=OUT.PFront;
for xpop=1:size(F,1)
if Nobj==1
figure(123); hold on;
plot(MODEDat.CounterGEN,log(min(F(:,1))),'dk', ...
'MarkerFaceColor','k'); grid on; hold on;
elseif Nobj==2
figure(123); hold on;
plot(F(xpop,1),F(xpop,2),'dk','MarkerFaceColor','k');...
grid on; hold on;
elseif Nobj==3
figure(123); hold on;
plot3(F(xpop,1),F(xpop,2),F(xpop,3),'dk','MarkerFaceColor','k');...
grid on; hold on;
end
end

%% Print and Display information
%
function [OUT Dat]=PrinterDisplay(OUT,Dat)

disp('------------------------------------------------')
disp(['Generation: ' num2str(Dat.CounterGEN)]);
disp(['FEs: ' num2str(Dat.CounterFES)]);
disp(['Pareto Front Size: ' mat2str(size(OUT.PFront,1))]);
disp('------------------------------------------------')

if mod(Dat.CounterGEN,1)==0
if Dat.NOBJ==3
figure(123);
plot3(OUT.PFront(:,1),OUT.PFront(:,2),OUT.PFront(:,3),'*r');
grid on;
elseif Dat.NOBJ==2
figure(123);
plot(OUT.PFront(:,1),OUT.PFront(:,2),'*r'); grid on;
elseif Dat.NOBJ==1
figure(123);
plot(Dat.CounterGEN,log(min(OUT.PFront(:,1))),'*r'); ...
grid on; hold on;
end
end

%% Dominance Filter
%
% A filter based on dominance criteria
%
function [Frente Conjunto]=DominanceFilter(F,C)

Xpop=size(F,1);
Nobj=size(F,2);
Nvar=size(C,2);
Frente=zeros(Xpop,Nobj);
Conjunto=zeros(Xpop,Nvar);
k=0;

for xpop=1:Xpop

for compara=1:Xpop
if F(xpop,:)==F(compara,:)
if xpop > compara
break;
end
else
if F(xpop,:)>=F(compara,:)
break;
end
end
end

k=k+1;
Frente(k,:)=F(xpop,:);
Conjunto(k,:)=C(xpop,:);
end
end
Frente=Frente(1:k,:);
Conjunto=Conjunto(1:k,:);

%% Release and bug report:
%
% November 2012: Initial release```