# Differential Evolution Optimizer

MODE.m

% multi-objective differential evolution optimizer

% last updated on 10/1/2017

clear;

global lbv; global ubv; global lbv_out; global ubv_out;

nvp = 100; % number of vectors in population

nit = 100; % number of iterations

idp = 4;

% # of variables, # of objective functions, and # of constraints

prbp = {30, 2, 0; … % ZDT1 problem

30, 2, 0; … % ZDT2 problem

30, 2, 0; … % ZDT3 problem

2, 2, 2}; … % TNK problem

switch (idp)

case 1

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 2

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 3

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 4

lbv = zeros(prbp{idp,1},1); ubv = [pi;pi];

case 5 %Project 2a

lbv = [2.6; 0.7; 17; 7.3; 7.3; 2.9; 5.0]; ubv = [3.6; 0.8; 28; 8.3; 8.3; 3.9; 5.5];

case 6 %Project 2b

lbv = [2.6; 0.7; 17; 7.3; 7.3; 2.9; 5.0]; ubv = [3.6; 0.8; 28; 8.3; 8.3; 3.9; 5.5];

lbv_out = 7; ubv_out = 8;

end

% create structures for design variables, objectives, and constraints

v_pop = struct(‘pop’,zeros(nvp,prbp{idp,1})); % structure for population vectors (variables)

obj_cons = struct(‘objf’,zeros(nvp,prbp{idp,2}),’consf’,zeros(nvp,1)); % structure for objective functions and constraints

vp = v_pop; tmpvp = v_pop; % parent and offspring population vectors

obc = obj_cons; tmpobc = obj_cons; % parent and offspring population vectors objective and constraint function values

vp.pop = rand(nvp,prbp{idp,1}); % generate initial population

for i=1:nvp, vp.pop(i,:) = vp.pop(i,:)’.*(ubv-lbv)+lbv; end

% evaluate objective functions and constraints

for i=1:nvp, [obc.objf(i,:) obc.consf(i,1)] = feval(‘optim_pr’,vp.pop(i,:),idp); end

% initialize algorithmic parameters

fm = 0.5; % mutation scale factor

cr = 0.2; % crossover rate

% arrays and vectors used during the Pareto ranking procedure of the nondominated solutions

fcount=zeros(2*nvp,1); dmflag = zeros(2*nvp,1); dupflag = zeros(2*nvp,1);

% start algorithmic iterations

for it=1:nit, it

% update the population using the “rand/1/bin” algorithm

for i=1:nvp

jrand = ceil(rand(1)*prbp{idp,1});

vr1 = ceil(rand(1)*nvp);

vr2 = ceil(rand(1)*nvp);

while(vr2==vr1), vr2 = ceil(rand(1)*nvp); end

vr3 = ceil(rand(1)*nvp);

while(vr3==vr1 || vr3==vr2) vr3 = ceil(rand(1)*nvp); end

for j=1:prbp{idp,1}

if(rand(1)<=cr || j==jrand)

tmpvp.pop(i,j) = vp.pop(vr1,j)+fm*(vp.pop(vr2,j)-vp.pop(vr3,j)); % differential mutation

else

tmpvp.pop(i,j) = vp.pop(i,j);

end

end

end

% check if offspring population vector values are within the boundary constraints

for i=1:nvp

for j=1:prbp{idp,1}

if(tmpvp.pop(i,j)<lbv(j))

tmpvp.pop(i,j) = lbv(j);

elseif(tmpvp.pop(i,j)>ubv(j))

tmpvp.pop(i,j) = ubv(j);

end

end

end

% compute objective and constraint functions of offspring

for i=1:nvp, [tmpobc.objf(i,:), tmpobc.consf(i,1)] = feval(‘optim_pr’,tmpvp.pop(i,:),idp); end

% combine parent and offspring populations

apv = vertcat(vp.pop,tmpvp.pop); apf = vertcat(obc.objf,tmpobc.objf); apc = vertcat(obc.consf,tmpobc.consf);

% initialize counters

ndfcount = 0; % number of nondominated fronts

ndcount = 0; % number of nondominated solutions

dcount = 2*nvp; % number of remaining solutions

tndcount = 0;

ndvp =[]; ndobf = []; ndobc = [];

% Pareto ranking procedure

while(ndcount<nvp) % add solutions and fronts until enough solutions have been found to form new parent population

dmflag(:,1) = 0; dupflag(:,1) = 0;

for j=1:dcount

if(dmflag(j,1)==0)

for jj=1:dcount

if(dmflag(jj,1)==0 && jj~=j)

if(apc(jj,1)<apc(j,1))

dmflag(j,1) = 1; % jth solution is dominated

break;

elseif(apc(jj,1)>apc(j,1))

dmflag(jj,1) = 1; % jjth solution is dominated

elseif(apc(j,1)==apc(jj,1)) % the two solutions have the same amount of constraint violation

ind1 = 0; ind2 = 0;

for ij=1:prbp{idp,2}

if(apf(j,ij)<apf(jj,ij))

ind1 = ind1 + 1;

elseif(apf(j,ij)>apf(jj,ij))

ind2 = ind2 + 1;

end

end

if(ind1>0 && ind2==0)

dmflag(jj,1) = 1;  % jjth solution is dominated

elseif(ind2>0 && ind1==0)

dmflag(j,1) = 1;  % jth solution is dominated

elseif(ind1==0 && ind2==0)

if(apv(jj,:)==apv(j,:))

dmflag(jj,1) = 1; % jjth solution is identical to jth solution

dupflag(jj,1) = 1; % do not include jjth solution in the new population

end

end

end

end

end

end

end

ndfcount = ndfcount + 1; % increase number of fronts by 1

fcount(ndfcount,1) = 0; % number of solutions in current front

tdcount = 0; % temporary variable to store number of dominated solutions

tempv = []; tempf = []; tempc = [];

tndcount = ndcount; % temporary variable to store the number of nondominated solutions in the previous front

for j=1:dcount

% store nondominated solutions of current front

if(dmflag(j,1)==0)

ndcount =  ndcount + 1;

ndvp(ndcount,:) = apv(j,:); ndobf(ndcount,:) = apf(j,:); ndobc(ndcount,1) = apc(j,1);

end

% store dominated solutions

if(dmflag(j,1)==1 && dupflag(j,1)==0)

tdcount = tdcount + 1;

tempv(tdcount,:) = apv(j,:); tempf(tdcount,:) =  apf(j,:); tempc(tdcount,1) =  apc(j,1);

end

end

fcount(ndfcount,1) = ndcount – tndcount;

dcount = tdcount;

apv = []; apf = []; apc = [];

apv = tempv; apf = tempf; apc = tempc;

end

% generate new population

if(ndcount>nvp)

saccp = fcount(ndfcount,1)-(ndcount-nvp); % number of solutions from the last front that will be kept

% compute the modified crowding distance metric to keep the solutions that reside in the least crowded regions

mcd = zeros(ndcount,1); indx = zeros(ndcount,1);

norml = max(ndobf)-min(ndobf); % find range of objective function values in current set of nondominated solutions

for i=1:prbp{idp,2}

if(norml(i)==0), norml(i) = 1; end

[~, indx] = sort(ndobf(1:ndcount,i));

for j=1:ndcount

if(j==1)

mcd(indx(j)) = mcd(indx(j)) + ((ndobf(indx(j+1),i)-ndobf(indx(j),i))/norml(i))^2;

elseif(j==ndcount)

mcd(indx(j)) = mcd(indx(j)) + ((ndobf(indx(j),i)-ndobf(indx(j-1),i))/norml(i))^2;

else

mcd(indx(j)) = mcd(indx(j)) + (ndobf(indx(j),i)-ndobf(indx(j-1),i))*(ndobf(indx(j+1),i)-ndobf(indx(j),i))/norml(i);

end

end

end

[dum, indx] = sort(mcd); indx = flipud(indx);

if(ndfcount==1) % only one nondominated front

for j=1:nvp

vp.pop(j,:) = ndvp(indx(j),:); obc.objf(j,:) = ndobf(indx(j),:); obc.consf(j,:) = ndobc(indx(j),1);

end

else % more than one nondominated fronts

vp.pop(1:nvp-saccp,:) = ndvp(1:nvp-saccp,:); obc.objf(1:nvp-saccp,:) = ndobf(1:nvp-saccp,:); obc.consf(1:nvp-saccp,1) = ndobc(1:nvp-saccp,1);

sindx = find(indx>(ndcount-fcount(ndfcount,1)));

for j=1:saccp

vp.pop(nvp-saccp+j,:) = ndvp(indx(sindx(j),1),:); obc.objf(nvp-saccp+j,:) = ndobf(indx(sindx(j),1),:); obc.consf(nvp-saccp+j,1) = ndobc(indx(sindx(j),1),1);

end

end

else

vp.pop = ndvp(1:nvp,:); obc.objf = ndobf(1:nvp,:); obc.consf = ndobc(1:nvp,1);

end

if(pfcount>0) % plot if the true Pareto front is available

plot(pf(:,1),pf(:,2),’rx’,obc.objf(:,1),obc.objf(:,2),’bo’,obc.objf(1:min(fcount(1,1),nvp),1),obc.objf(1:min(fcount(1,1),nvp),2),’gs’);

xlabel(‘f_1’);

ylabel(‘f_2′);

drawnow();

else % plot if the true Pareto front is not available

plot(obc.objf(:,1),obc.objf(:,2),’o’,obc.objf(1:min(fcount(1,1),nvp),1),obc.objf(1:min(fcount(1,1),nvp),2),’+’);

xlabel(‘f_1’);

ylabel(‘f_2’);

drawnow();

end

end

if(pfcount>0) feval(‘metrics’,obc.objf(1:min(fcount(1,1),nvp),:),pf); end % evaluate metrics

csvwrite(‘functions’,obc.objf(1:min(fcount(1,1),nvp),:)); % print objective function values

csvwrite(‘variables’,vp.pop(1:min(fcount(1,1),nvp),:)); % print decision variable values

MOPSO.m

% multi-objective differential evolution optimizer

% last updated on 10/1/2017

clear;

global lbv; global ubv; global lbv_out; global ubv_out;

nvp = 100; % number of particles in population

nit = 300; % number of iterations

idp = 4;

% # of variables, # of objective functions, and # of constraints

prbp = {30, 2, 0; … % ZDT1 problem

30, 2, 0; … % ZDT2 problem

30, 2, 0; … % ZDT3 problem

2, 2, 2}; … % TNK problem

switch (idp)

case 1

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 2

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 3

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 4

lbv = zeros(prbp{idp,1},1); ubv = [pi;pi];

case 5 %Project 2a

lbv = [2.6; 0.7; 17; 7.3; 7.3; 2.9; 5.0]; ubv = [3.6; 0.8; 28; 8.3; 8.3; 3.9; 5.5];

case 6 %Project 2b

lbv = [2.6; 0.7; 17; 7.3; 7.3; 2.9; 5.0]; ubv = [3.6; 0.8; 28; 8.3; 8.3; 3.9; 5.5];

lbv_out = 7; ubv_out = 8;

end

% create structures for design variables, objectives, and constraints

v_pop = struct(‘pop’,zeros(nvp,prbp{idp,1})); % structure for population vectors (variables)

obj_cons = struct(‘objf’,zeros(nvp,prbp{idp,2}),’consf’,zeros(nvp,1)); % structure for objective functions and constraints

vp = v_pop; archvp = v_pop; % parent and offspring population vectors

obc = obj_cons; archobc = obj_cons; % parent and offspring population vectors objective and constraint function values

vp.pop = rand(nvp,prbp{idp,1}); % generate initial population

for i=1:nvp, vp.pop(i,:) = vp.pop(i,:)’.*(ubv-lbv)+lbv; end

% evaluate objective functions and constraints

for i=1:nvp, [obc.objf(i,:), obc.consf(i,1)] = feval(‘optim_pr’,vp.pop(i,:),idp); end

archvp = vp; archobc = obc;

% initialize algorithmic parameters

c1 = 1.5; c2 = 1.5; w = 0.8; v_clm = 5; % initialize algorithmic parameters

vc = v_pop; vc.pop = zeros(nvp,prbp{idp,1}); % create and initialize array for the velocity vector

% Velocity clamping

range = ubv – lbv; % find range of decision variables

Vmax = v_clm*range;

Vmax_v = (Vmax*ones(1,nvp))’;

% arrays and vectors used during the Pareto ranking procedure of the nondominated solutions

fcount=zeros(2*nvp,1); dmflag = zeros(2*nvp,1); dupflag = zeros(2*nvp,1);

% start algorithmic iterations

for it=1:nit, it

% update the population using the original PSO

for i=1:nvp % generate trial vectors

vr1 = ceil(rand(1)*nvp);

vr2 = ceil(rand(1)*nvp); while vr2==vr1, vr2 = ceil(rand(1)*nvp); end

vc.pop(i,:) = w*vc.pop(i,:) + c1*rand(1,prbp{idp,1}).*(archvp.pop(vr1,:)-vp.pop(i,:)) + c2*rand(1,prbp{idp,1}).*(archvp.pop(vr2,:)-vp.pop(i,:));

end

vp.pop = vp.pop + vc.pop;

for i=1:nvp % check boundary values

for j=1:prbp{idp,1}

if vp.pop(i,j)<lbv(j), vp.pop(i,j) = lbv(j); end

if vp.pop(i,j)>ubv(j), vp.pop(i,j) = ubv(j); end

end

end

% compute objective and constraint

for i=1:nvp, [obc.objf(i,:), obc.consf(i,:)] = feval(‘optim_pr’,vp.pop(i,:),idp); end % evaluate objective functions and constraints

% combine archive and current particle positions

apv = vertcat(vp.pop,archvp.pop); apf = vertcat(obc.objf,archobc.objf); apc = vertcat(obc.consf,archobc.consf);

% initialize counters

ndfcount = 0; % number of nondominated fronts

ndcount = 0; % number of nondominated solutions

dcount = 2*nvp; % number of remaining solutions

tndcount = 0;

ndvp =[]; ndobf = []; ndobc = [];

% Pareto ranking procedure

while(ndcount<nvp) % add solutions and fronts until enough solutions have been found to form new parent population

dmflag(:,1) = 0; dupflag(:,1) = 0;

for j=1:dcount

if(dmflag(j,1)==0)

for jj=1:dcount

if(dmflag(jj,1)==0 && jj~=j)

if(apc(jj,1)<apc(j,1))

dmflag(j,1) = 1; % jth solution is dominated

break;

elseif(apc(jj,1)>apc(j,1))

dmflag(jj,1) = 1; % jjth solution is dominated

elseif(apc(j,1)==apc(jj,1)) % the two solutions have the same amount of constraint violation

ind1 = 0; ind2 = 0;

for ij=1:prbp{idp,2}

if(apf(j,ij)<apf(jj,ij))

ind1 = ind1 + 1;

elseif(apf(j,ij)>apf(jj,ij))

ind2 = ind2 + 1;

end

end

if(ind1>0 && ind2==0)

dmflag(jj,1) = 1;  % jjth solution is dominated

elseif(ind2>0 && ind1==0)

dmflag(j,1) = 1;  % jth solution is dominated

elseif(ind1==0 && ind2==0)

if(apv(jj,:)==apv(j,:))

dmflag(jj,1) = 1; % jjth solution is identical to jth solution

dupflag(jj,1) = 1; % do not include jjth solution in the new population

end

end

end

end

end

end

end

ndfcount = ndfcount + 1; % increase number of fronts by 1

fcount(ndfcount,1) = 0; % number of solutions in current front

tdcount = 0; % temporary variable to store number of dominated solutions

tempv = []; tempf = []; tempc = [];

tndcount = ndcount; % temporary variable to store the number of nondominated solutions in the previous front

for j=1:dcount

% store nondominated solutions of current front

if(dmflag(j,1)==0)

ndcount =  ndcount + 1;

ndvp(ndcount,:) = apv(j,:); ndobf(ndcount,:) = apf(j,:); ndobc(ndcount,1) = apc(j,1);

end

% store dominated solutions

if(dmflag(j,1)==1 && dupflag(j,1)==0)

tdcount = tdcount + 1;

tempv(tdcount,:) = apv(j,:); tempf(tdcount,:) =  apf(j,:); tempc(tdcount,1) =  apc(j,1);

end

end

fcount(ndfcount,1) = ndcount – tndcount;

dcount = tdcount;

apv = []; apf = []; apc = [];

apv = tempv; apf = tempf; apc = tempc;

end

% update archive

if(ndcount>nvp)

saccp = fcount(ndfcount,1)-(ndcount-nvp); % number of solutions from the last front that will be kept

% compute the modified crowding distance metric to keep the solutions that reside in the least crowded regions

mcd = zeros(ndcount,1); indx = zeros(ndcount,1);

norml = max(ndobf)-min(ndobf); % find range of objective function values in current set of nondominated solutions

for i=1:prbp{idp,2}

if(norml(i)==0), norml(i) = 1; end

[~, indx] = sort(ndobf(1:ndcount,i));

for j=1:ndcount

if(j==1)

mcd(indx(j)) = mcd(indx(j)) + ((ndobf(indx(j+1),i)-ndobf(indx(j),i))/norml(i))^2;

elseif(j==ndcount)

mcd(indx(j)) = mcd(indx(j)) + ((ndobf(indx(j),i)-ndobf(indx(j-1),i))/norml(i))^2;

else

mcd(indx(j)) = mcd(indx(j)) + (ndobf(indx(j),i)-ndobf(indx(j-1),i))*(ndobf(indx(j+1),i)-ndobf(indx(j),i))/norml(i);

end

end

end

[dum, indx] = sort(mcd); indx = flipud(indx);

if(ndfcount==1) % only one nondominated front

for j=1:nvp

archvp.pop(j,:) = ndvp(indx(j),:); archobc.objf(j,:) = ndobf(indx(j),:); archobc.consf(j,:) = ndobc(indx(j),1);

end

else % more than one nondominated fronts

archvp.pop(1:nvp-saccp,:) = ndvp(1:nvp-saccp,:); archobc.objf(1:nvp-saccp,:) = ndobf(1:nvp-saccp,:); archobc.consf(1:nvp-saccp,1) = ndobc(1:nvp-saccp,1);

sindx = find(indx>(ndcount-fcount(ndfcount,1)));

for j=1:saccp

archvp.pop(nvp-saccp+j,:) = ndvp(indx(sindx(j),1),:); archobc.objf(nvp-saccp+j,:) = ndobf(indx(sindx(j),1),:); archobc.consf(nvp-saccp+j,1) = ndobc(indx(sindx(j),1),1);

end

end

else

archvp.pop = ndvp(1:nvp,:); archobc.objf = ndobf(1:nvp,:); archobc.consf = ndobc(1:nvp,1);

end

if(pfcount>0) % plot if the true Pareto front is available

plot(pf(:,1),pf(:,2),’rx’,archobc.objf(:,1),archobc.objf(:,2),’bo’,archobc.objf(1:min(fcount(1,1),nvp),1),archobc.objf(1:min(fcount(1,1),nvp),2),’gs’);

xlabel(‘f_1’);

ylabel(‘f_2′);

drawnow();

else % plot if the true Pareto front is not available

plot(archobc.objf(:,1),archobc.objf(:,2),’o’,archobc.objf(1:min(fcount(1,1),nvp),1),archobc.objf(1:min(fcount(1,1),nvp),2),’+’);

xlabel(‘f_1’);

ylabel(‘f_2’);

drawnow();

end

end

if(pfcount>0) feval(‘metrics’,archobc.objf(1:min(fcount(1,1),nvp),:),pf); end % evaluate metrics

csvwrite(‘functions’,archobc.objf(1:min(fcount(1,1),nvp),:)); % print objective function values

csvwrite(‘variables’,archvp.pop(1:min(fcount(1,1),nvp),:)); % print decision variable values

Solution

metrics.m

function metrics(ofv,pf)

% compute the generational distance (gd) and the inverted generational distance (igd)

gd = 0; mindpf = zeros(size(ofv,1),1);

nrml = max(pf)-min(pf); % normalize using the current range of values in each objective function

for ji=1:size(ofv,1)

for j=1:size(pf,1)

pfd = sqrt(sum(((ofv(ji,:)-pf(j,:))./nrml).^2));

if(pfd<mindpf(ji)||j==1), mindpf(ji) = pfd; end

end

gd = gd + mindpf(ji,1)^2;

end

gd = sqrt(gd)/size(ofv,1)

igd=0; mindpf = zeros(size(pf,1),1);

for ji=1:size(pf,1)

for j=1:size(ofv,1)

pfd = sqrt(sum(((ofv(j,:)-pf(ji,:))./nrml).^2));

if(pfd<mindpf(ji)||j==1), mindpf(ji,1) = pfd; end

end

igd = igd + mindpf(ji,1)^2;

end

igd = sqrt(igd)/size(pf,1)

end

MODE.m

% multi-objective differential evolution optimizer

% last updated on 10/1/2017

clear;

global lbv; global ubv; global lbv_out; global ubv_out;

nvp = 100; % number of vectors in population

nit = 1000; % number of iterations

idp = 5;

% # of variables, # of objective functions, and # of constraints

prbp = {   30, 2,  0; … % ZDT1 problem

30, 2,  0; … % ZDT2 problem

30, 2,  0; … % ZDT3 problem

2, 2,  2; … % TNK problem

7, 2, 11; … % Speed Reducer

4, 2, 4};     % Ravindran-Reklaitis

switch (idp)

case 1

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 2

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 3

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 4

lbv = zeros(prbp{idp,1},1); ubv = [pi;pi];

case 5 %Speed Reducer

lbv = [2.6; 0.7; 17; 7.3; 7.3; 2.9; 5.0];

ubv = [3.6; 0.8; 28; 8.3; 8.3; 3.9; 5.5];

pf = []; pfcount = size(pf,1);

case 6 %Ravindran-Reklaitis

lbv = [0.1 0.1 0.1 0.1]’;

ubv = [2 10 10 2]’;

pf = []; pfcount = size(pf,1);

end

% create structures for design variables, objectives, and constraints

v_pop = struct(‘pop’,zeros(nvp,prbp{idp,1})); % structure for population vectors (variables)

obj_cons = struct(‘objf’,zeros(nvp,prbp{idp,2}),’consf’,zeros(nvp,1)); % structure for objective functions and constraints

vp = v_pop; tmpvp = v_pop; % parent and offspring population vectors

obc = obj_cons; tmpobc = obj_cons; % parent and offspring population vectors objective and constraint function values

vp.pop = rand(nvp,prbp{idp,1}); % generate initial population

for i=1:nvp, vp.pop(i,:) = vp.pop(i,:)’.*(ubv-lbv)+lbv; end

% evaluate objective functions and constraints

for i=1:nvp, [obc.objf(i,:), obc.consf(i,1)] = feval(‘optim_pr’,vp.pop(i,:),idp); end

% initialize algorithmic parameters

fm = 0.5; % mutation scale factor

cr = 0.2; % crossover rate

% arrays and vectors used during the Pareto ranking procedure of the nondominated solutions

fcount=zeros(2*nvp,1); dmflag = zeros(2*nvp,1); dupflag = zeros(2*nvp,1);

% start algorithmic iterations

for it=1:nit, it

% update the population using the “rand/1/bin” algorithm

for i=1:nvp

jrand = ceil(rand(1)*prbp{idp,1});

vr1 = ceil(rand(1)*nvp);

vr2 = ceil(rand(1)*nvp);

while(vr2==vr1), vr2 = ceil(rand(1)*nvp); end

vr3 = ceil(rand(1)*nvp);

while(vr3==vr1 || vr3==vr2) vr3 = ceil(rand(1)*nvp); end

for j=1:prbp{idp,1}

if(rand(1)<=cr || j==jrand)

tmpvp.pop(i,j) = vp.pop(vr1,j)+fm*(vp.pop(vr2,j)-vp.pop(vr3,j)); % differential mutation

else

tmpvp.pop(i,j) = vp.pop(i,j);

end

end

end

% check if offspring population vector values are within the boundary constraints

for i=1:nvp

for j=1:prbp{idp,1}

if(tmpvp.pop(i,j)<lbv(j))

tmpvp.pop(i,j) = lbv(j);

elseif(tmpvp.pop(i,j)>ubv(j))

tmpvp.pop(i,j) = ubv(j);

end

end

end

% compute objective and constraint functions of offspring

for i=1:nvp, [tmpobc.objf(i,:), tmpobc.consf(i,1)] = feval(‘optim_pr’,tmpvp.pop(i,:),idp); end

% combine parent and offspring populations

apv = vertcat(vp.pop,tmpvp.pop); apf = vertcat(obc.objf,tmpobc.objf); apc = vertcat(obc.consf,tmpobc.consf);

% initialize counters

ndfcount = 0; % number of nondominated fronts

ndcount = 0; % number of nondominated solutions

dcount = 2*nvp; % number of remaining solutions

tndcount = 0;

ndvp =[]; ndobf = []; ndobc = [];

% Pareto ranking procedure

while(ndcount<nvp) % add solutions and fronts until enough solutions have been found to form new parent population

dmflag(:,1) = 0; dupflag(:,1) = 0;

for j=1:dcount

if(dmflag(j,1)==0)

for jj=1:dcount

if(dmflag(jj,1)==0 && jj~=j)

if(apc(jj,1)<apc(j,1))

dmflag(j,1) = 1; % jth solution is dominated

break;

elseif(apc(jj,1)>apc(j,1))

dmflag(jj,1) = 1; % jjth solution is dominated

elseif(apc(j,1)==apc(jj,1)) % the two solutions have the same amount of constraint violation

ind1 = 0; ind2 = 0;

for ij=1:prbp{idp,2}

if(apf(j,ij)<apf(jj,ij))

ind1 = ind1 + 1;

elseif(apf(j,ij)>apf(jj,ij))

ind2 = ind2 + 1;

end

end

if(ind1>0 && ind2==0)

dmflag(jj,1) = 1;  % jjth solution is dominated

elseif(ind2>0 && ind1==0)

dmflag(j,1) = 1;  % jth solution is dominated

elseif(ind1==0 && ind2==0)

if(apv(jj,:)==apv(j,:))

dmflag(jj,1) = 1; % jjth solution is identical to jth solution

dupflag(jj,1) = 1; % do not include jjth solution in the new population

end

end

end

end

end

end

end

ndfcount = ndfcount + 1; % increase number of fronts by 1

fcount(ndfcount,1) = 0; % number of solutions in current front

tdcount = 0; % temporary variable to store number of dominated solutions

tempv = []; tempf = []; tempc = [];

tndcount = ndcount; % temporary variable to store the number of nondominated solutions in the previous front

for j=1:dcount

% store nondominated solutions of current front

if(dmflag(j,1)==0)

ndcount =  ndcount + 1;

ndvp(ndcount,:) = apv(j,:); ndobf(ndcount,:) = apf(j,:); ndobc(ndcount,1) = apc(j,1);

end

% store dominated solutions

if(dmflag(j,1)==1 && dupflag(j,1)==0)

tdcount = tdcount + 1;

tempv(tdcount,:) = apv(j,:); tempf(tdcount,:) =  apf(j,:); tempc(tdcount,1) =  apc(j,1);

end

end

fcount(ndfcount,1) = ndcount – tndcount;

dcount = tdcount;

apv = []; apf = []; apc = [];

apv = tempv; apf = tempf; apc = tempc;

end

% generate new population

if(ndcount>nvp)

saccp = fcount(ndfcount,1)-(ndcount-nvp); % number of solutions from the last front that will be kept

% compute the modified crowding distance metric to keep the solutions that reside in the least crowded regions

mcd = zeros(ndcount,1); indx = zeros(ndcount,1);

norml = max(ndobf)-min(ndobf); % find range of objective function values in current set of nondominated solutions

for i=1:prbp{idp,2}

if(norml(i)==0), norml(i) = 1; end

[~, indx] = sort(ndobf(1:ndcount,i));

for j=1:ndcount

if(j==1)

mcd(indx(j)) = mcd(indx(j)) + ((ndobf(indx(j+1),i)-ndobf(indx(j),i))/norml(i))^2;

elseif(j==ndcount)

mcd(indx(j)) = mcd(indx(j)) + ((ndobf(indx(j),i)-ndobf(indx(j-1),i))/norml(i))^2;

else

mcd(indx(j)) = mcd(indx(j)) + (ndobf(indx(j),i)-ndobf(indx(j-1),i))*(ndobf(indx(j+1),i)-ndobf(indx(j),i))/norml(i);

end

end

end

[dum, indx] = sort(mcd); indx = flipud(indx);

if(ndfcount==1) % only one nondominated front

for j=1:nvp

vp.pop(j,:) = ndvp(indx(j),:); obc.objf(j,:) = ndobf(indx(j),:); obc.consf(j,:) = ndobc(indx(j),1);

end

else % more than one nondominated fronts

vp.pop(1:nvp-saccp,:) = ndvp(1:nvp-saccp,:); obc.objf(1:nvp-saccp,:) = ndobf(1:nvp-saccp,:); obc.consf(1:nvp-saccp,1) = ndobc(1:nvp-saccp,1);

sindx = find(indx>(ndcount-fcount(ndfcount,1)));

for j=1:saccp

vp.pop(nvp-saccp+j,:) = ndvp(indx(sindx(j),1),:); obc.objf(nvp-saccp+j,:) = ndobf(indx(sindx(j),1),:); obc.consf(nvp-saccp+j,1) = ndobc(indx(sindx(j),1),1);

end

end

else

vp.pop = ndvp(1:nvp,:); obc.objf = ndobf(1:nvp,:); obc.consf = ndobc(1:nvp,1);

end

if(pfcount>0) % plot if the true Pareto front is available

plot(pf(:,1),pf(:,2),’rx’,obc.objf(:,1),obc.objf(:,2),’bo’,obc.objf(1:min(fcount(1,1),nvp),1),obc.objf(1:min(fcount(1,1),nvp),2),’gs’);

xlabel(‘f_1’);

ylabel(‘f_2′);

drawnow();

else % plot if the true Pareto front is not available

plot(obc.objf(:,1),obc.objf(:,2),’o’,obc.objf(1:min(fcount(1,1),nvp),1),obc.objf(1:min(fcount(1,1),nvp),2),’+’);

xlabel(‘f_1’);

ylabel(‘f_2’);

drawnow();

end

end

if(pfcount>0) feval(‘metrics’,obc.objf(1:min(fcount(1,1),nvp),:),pf); end % evaluate metrics

csvwrite(‘functions.csv’,obc.objf(1:min(fcount(1,1),nvp),:)); % print objective function values

csvwrite(‘variables.csv’,vp.pop(1:min(fcount(1,1),nvp),:)); % print decision variable values

MOPSO.m

% multi-objective differential evolution optimizer

% last updated on 10/1/2017

clear;

global lbv; global ubv; global lbv_out; global ubv_out;

nvp = 100; % number of particles in population

nit = 1000; % number of iterations

idp = 5;

% # of variables, # of objective functions, and # of constraints

prbp = {   30, 2,  0; … % ZDT1 problem

30, 2,  0; … % ZDT2 problem

30, 2,  0; … % ZDT3 problem

2, 2,  2; … % TNK problem

7, 2, 11; … % Speed Reducer

4, 2, 4};     % Ravindran-Reklaitis

switch (idp)

case 1

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 2

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 3

lbv = zeros(prbp{idp,1},1); ubv = ones(prbp{idp,1},1);

case 4

lbv = zeros(prbp{idp,1},1); ubv = [pi;pi];

case 5 %Speed Reducer

lbv = [2.6; 0.7; 17; 7.3; 7.3; 2.9; 5.0];

ubv = [3.6; 0.8; 28; 8.3; 8.3; 3.9; 5.5];

pf = []; pfcount = size(pf,1);

case 6 %Ravindran-Reklaitis

lbv = [0.1 0.1 0.1 0.1]’;

ubv = [1 10 10 2]’;

pf = []; pfcount = size(pf,1);

end

% create structures for design variables, objectives, and constraints

v_pop = struct(‘pop’,zeros(nvp,prbp{idp,1})); % structure for population vectors (variables)

obj_cons = struct(‘objf’,zeros(nvp,prbp{idp,2}),’consf’,zeros(nvp,1)); % structure for objective functions and constraints

vp = v_pop; archvp = v_pop; % parent and offspring population vectors

obc = obj_cons; archobc = obj_cons; % parent and offspring population vectors objective and constraint function values

vp.pop = rand(nvp,prbp{idp,1}); % generate initial population

for i=1:nvp, vp.pop(i,:) = vp.pop(i,:)’.*(ubv-lbv)+lbv; end

% evaluate objective functions and constraints

for i=1:nvp, [obc.objf(i,:), obc.consf(i,1)] = feval(‘optim_pr’,vp.pop(i,:),idp); end

archvp = vp; archobc = obc;

% initialize algorithmic parameters

c1 = 1.5; c2 = 1.5; w = 0.8; v_clm = 5; % initialize algorithmic parameters

vc = v_pop; vc.pop = zeros(nvp,prbp{idp,1}); % create and initialize array for the velocity vector

% Velocity clamping

range = ubv – lbv; % find range of decision variables

Vmax = v_clm*range;

Vmax_v = (Vmax*ones(1,nvp))’;

% arrays and vectors used during the Pareto ranking procedure of the nondominated solutions

fcount=zeros(2*nvp,1); dmflag = zeros(2*nvp,1); dupflag = zeros(2*nvp,1);

% start algorithmic iterations

for it=1:nit, it

% update the population using the original PSO

for i=1:nvp % generate trial vectors

vr1 = ceil(rand(1)*nvp);

vr2 = ceil(rand(1)*nvp); while vr2==vr1, vr2 = ceil(rand(1)*nvp); end

vc.pop(i,:) = w*vc.pop(i,:) + c1*rand(1,prbp{idp,1}).*(archvp.pop(vr1,:)-vp.pop(i,:)) + c2*rand(1,prbp{idp,1}).*(archvp.pop(vr2,:)-vp.pop(i,:));

end

vp.pop = vp.pop + vc.pop;

for i=1:nvp % check boundary values

for j=1:prbp{idp,1}

if vp.pop(i,j)<lbv(j), vp.pop(i,j) = lbv(j); end

if vp.pop(i,j)>ubv(j), vp.pop(i,j) = ubv(j); end

end

end

% compute objective and constraint

for i=1:nvp, [obc.objf(i,:), obc.consf(i,:)] = feval(‘optim_pr’,vp.pop(i,:),idp); end % evaluate objective functions and constraints

% combine archive and current particle positions

apv = vertcat(vp.pop,archvp.pop); apf = vertcat(obc.objf,archobc.objf); apc = vertcat(obc.consf,archobc.consf);

% initialize counters

ndfcount = 0; % number of nondominated fronts

ndcount = 0; % number of nondominated solutions

dcount = 2*nvp; % number of remaining solutions

tndcount = 0;

ndvp =[]; ndobf = []; ndobc = [];

% Pareto ranking procedure

while(ndcount<nvp) % add solutions and fronts until enough solutions have been found to form new parent population

dmflag(:,1) = 0; dupflag(:,1) = 0;

for j=1:dcount

if(dmflag(j,1)==0)

for jj=1:dcount

if(dmflag(jj,1)==0 && jj~=j)

if(apc(jj,1)<apc(j,1))

dmflag(j,1) = 1; % jth solution is dominated

break;

elseif(apc(jj,1)>apc(j,1))

dmflag(jj,1) = 1; % jjth solution is dominated

elseif(apc(j,1)==apc(jj,1)) % the two solutions have the same amount of constraint violation

ind1 = 0; ind2 = 0;

for ij=1:prbp{idp,2}

if(apf(j,ij)<apf(jj,ij))

ind1 = ind1 + 1;

elseif(apf(j,ij)>apf(jj,ij))

ind2 = ind2 + 1;

end

end

if(ind1>0 && ind2==0)

dmflag(jj,1) = 1;  % jjth solution is dominated

elseif(ind2>0 && ind1==0)

dmflag(j,1) = 1;  % jth solution is dominated

elseif(ind1==0 && ind2==0)

if(apv(jj,:)==apv(j,:))

dmflag(jj,1) = 1; % jjth solution is identical to jth solution

dupflag(jj,1) = 1; % do not include jjth solution in the new population

end

end

end

end

end

end

end

ndfcount = ndfcount + 1; % increase number of fronts by 1

fcount(ndfcount,1) = 0; % number of solutions in current front

tdcount = 0; % temporary variable to store number of dominated solutions

tempv = []; tempf = []; tempc = [];

tndcount = ndcount; % temporary variable to store the number of nondominated solutions in the previous front

for j=1:dcount

% store nondominated solutions of current front

if(dmflag(j,1)==0)

ndcount =  ndcount + 1;

ndvp(ndcount,:) = apv(j,:); ndobf(ndcount,:) = apf(j,:); ndobc(ndcount,1) = apc(j,1);

end

% store dominated solutions

if(dmflag(j,1)==1 && dupflag(j,1)==0)

tdcount = tdcount + 1;

tempv(tdcount,:) = apv(j,:); tempf(tdcount,:) =  apf(j,:); tempc(tdcount,1) =  apc(j,1);

end

end

fcount(ndfcount,1) = ndcount – tndcount;

dcount = tdcount;

apv = []; apf = []; apc = [];

apv = tempv; apf = tempf; apc = tempc;

end

% update archive

if(ndcount>nvp)

saccp = fcount(ndfcount,1)-(ndcount-nvp); % number of solutions from the last front that will be kept

% compute the modified crowding distance metric to keep the solutions that reside in the least crowded regions

mcd = zeros(ndcount,1); indx = zeros(ndcount,1);

norml = max(ndobf)-min(ndobf); % find range of objective function values in current set of nondominated solutions

for i=1:prbp{idp,2}

if(norml(i)==0), norml(i) = 1; end

[~, indx] = sort(ndobf(1:ndcount,i));

for j=1:ndcount

if(j==1)

mcd(indx(j)) = mcd(indx(j)) + ((ndobf(indx(j+1),i)-ndobf(indx(j),i))/norml(i))^2;

elseif(j==ndcount)

mcd(indx(j)) = mcd(indx(j)) + ((ndobf(indx(j),i)-ndobf(indx(j-1),i))/norml(i))^2;

else

mcd(indx(j)) = mcd(indx(j)) + (ndobf(indx(j),i)-ndobf(indx(j-1),i))*(ndobf(indx(j+1),i)-ndobf(indx(j),i))/norml(i);

end

end

end

[dum, indx] = sort(mcd); indx = flipud(indx);

if(ndfcount==1) % only one nondominated front

for j=1:nvp

archvp.pop(j,:) = ndvp(indx(j),:); archobc.objf(j,:) = ndobf(indx(j),:); archobc.consf(j,:) = ndobc(indx(j),1);

end

else % more than one nondominated fronts

archvp.pop(1:nvp-saccp,:) = ndvp(1:nvp-saccp,:); archobc.objf(1:nvp-saccp,:) = ndobf(1:nvp-saccp,:); archobc.consf(1:nvp-saccp,1) = ndobc(1:nvp-saccp,1);

sindx = find(indx>(ndcount-fcount(ndfcount,1)));

for j=1:saccp

archvp.pop(nvp-saccp+j,:) = ndvp(indx(sindx(j),1),:); archobc.objf(nvp-saccp+j,:) = ndobf(indx(sindx(j),1),:); archobc.consf(nvp-saccp+j,1) = ndobc(indx(sindx(j),1),1);

end

end

else

archvp.pop = ndvp(1:nvp,:); archobc.objf = ndobf(1:nvp,:); archobc.consf = ndobc(1:nvp,1);

end

if(pfcount>0) % plot if the true Pareto front is available

plot(pf(:,1),pf(:,2),’rx’,archobc.objf(:,1),archobc.objf(:,2),’bo’,archobc.objf(1:min(fcount(1,1),nvp),1),archobc.objf(1:min(fcount(1,1),nvp),2),’gs’);

xlabel(‘f_1’);

ylabel(‘f_2′);

drawnow();

else % plot if the true Pareto front is not available

plot(archobc.objf(:,1),archobc.objf(:,2),’o’,archobc.objf(1:min(fcount(1,1),nvp),1),archobc.objf(1:min(fcount(1,1),nvp),2),’+’);

xlabel(‘f_1’);

ylabel(‘f_2’);

drawnow();

end

end

if(pfcount>0) feval(‘metrics’,archobc.objf(1:min(fcount(1,1),nvp),:),pf); end % evaluate metrics

csvwrite(‘functions.csv’,archobc.objf(1:min(fcount(1,1),nvp),:)); % print objective function values

csvwrite(‘variables.csv’,archvp.pop(1:min(fcount(1,1),nvp),:)); % print decision variable values

optim_pr.m

function [obf,obc] = optim_pr(vp, idp) % minimization of all objective functions

obc = 0;

switch (idp)

case 1

obf(1) = vp(1);

esum = 1+(9/29)*sum(vp(2:end),2);

obf(2) = esum*(1-sqrt(obf(1)/esum));

case 2

obf(1) = vp(1);

esum = 1+(9/29)*sum(vp(2:end),2);

obf(2) = esum*(1-(obf(1)/esum)^2);

case 3

obf(1) = vp(1);

esum = 1+(9/29)*sum(vp(2:end),2);

obf(2) = esum*(1-sqrt(obf(1)/esum)-(obf(1)/esum)*sin(10*pi*obf(1)));

case 4

obf(1) = vp(1);

obf(2) = vp(2);

obc = abs(min(vp(1)^2+vp(2)^2-1-0.1*cos(16*atan2(vp(1),vp(2))),0));

obc = obc + abs(max(((vp(1)-0.5)^2+(vp(2)-0.5)^2)/0.5-1,0));

case 5 %Speed Reducer

obf(1) = 7.854*vp(1)*vp(2)^2*(vp(3)^2/3+1.493*vp(3)-4.309)-1.508*…

vp(1)*(vp(6)^2+vp(7)^2)+7.477*(vp(6)^3+vp(7)^3)+0.7854*…

(vp(4)*vp(6)^2+vp(5)*vp(7)^2);

obf(2) = sqrt(((745*vp(4)/vp(2)/vp(3))^2 + 1.69e7)/…

(0.1*vp(6)^3));

obc = max(27/(prod(vp(1:3))*vp(2))-1, 0);

obc = obc + max(397.5/(vp(1)*vp(2)^2*vp(3)^2)-1, 0);

obc = obc + max(1.93*vp(4)^3/(vp(2)*vp(3)*vp(6)^4)-1, 0);

obc = obc + max(1.93*vp(5)^3/(vp(2)*vp(3)*vp(7)^4)-1, 0);

obc = obc + max(1/110/vp(6)^3*sqrt(((745*vp(4))/(vp(2)*vp(3)))^2+16.9e6)-1, 0);

obc = obc + max(1/85/vp(7)^3*sqrt(((745*vp(5))/(vp(2)*vp(3)))^2+157.5e6)-1, 0);

obc = obc + max(vp(2)*vp(3)/40-1, 0);

obc = obc + max(5*vp(2)/vp(1)-1, 0);

obc = obc + max(vp(1)/12/vp(2)-1, 0);

obc = obc + max((1.5*vp(6)+1.9)/vp(4)-1, 0);

obc = obc + max((1.1*vp(7)+1.9)/vp(5)-1, 0);

case 6 %Ravindran-Reklaitis

% constants

F = 6000;

L = 14;

tau_d = 13.6e3;

sigma_d = 3e4;

E = 30e6;

G = 12e6;

% anonymous functions

M = @(x) F*(L + x(2)/2);

R = @(x) sqrt(x(2)^2/4 + ((x(1) + x(3))/2)^2);

J = @(x) 2*(0.707*x(1)*x(2)*(x(2)^2/12 + …

((x(1) + x(3))/2)^2));

tau_p  = @(x) F/sqrt(2)/x(1)/x(2);

tau_pp = @(x) M(x)*R(x)/J(x);

t = @(x) x(2)/2/R(x);

tau = @(x) sqrt(tau_p(x)^2 + 2*tau_p(x)*tau_pp(x)*t(x) + …

tau_pp(x)^2);

sigma = @(x) 6*F*L/x(4)/x(3)^2;

I = @(x) 1/12*x(3)*x(4)^3;

alpha = @(x) 1/3*G*x(3)*x(4)^3;

Pc = @(x) 4.013*sqrt(E*I(x)*alpha(x))/L^2*(1 – x(3)/2/L*…

sqrt(E*I(x)/alpha(x)));

DEL = @(x) 4*F*L^3/E/x(4)/x(3)^3;

% objective functions

c3 = 0.37*0.283;

c4 = 0.12*0.283;

obf(1) = (1+c3)*vp(1)^2*vp(2) + c4*vp(3)*vp(4)*(L + vp(2)); % DEL(vp)

obf(2) = DEL(vp);

% constraints

obc = obc + min(tau_d – tau(vp), 0);

obc = obc + min(sigma_d – sigma(vp), 0);

obc = obc + min(vp(4) – vp(1), 0);

obc = obc + min(vp(2), 0);

obc = obc + min(vp(3), 0);

obc = obc + min(Pc(vp) – F, 0);

obc = obc + min(vp(1) – 0.125);

obc = obc + min(0.25 – DEL(vp));

end

end