% Plotting and verbose
if(size(POS_fit,2)==2)
h_fig = figure(1);
h_par = plot(POS_fit(:,1),POS_fit(:,2),'or'); hold on;
h_rep = plot(REP.pos_fit(:,1),REP.pos_fit(:,2),'ok'); hold on;
try
set(gca,'xtick',REP.hypercube_limits(:,1)','ytick',REP.hypercube_limits(:,2)');
axis([min(REP.hypercube_limits(:,1)) max(REP.hypercube_limits(:,1)) ...
min(REP.hypercube_limits(:,2)) max(REP.hypercube_limits(:,2))]);
grid on; xlabel('f1'); ylabel('f2');
end
drawnow;
end
if(size(POS_fit,2)==3)
h_fig = figure(1);
h_par = plot3(POS_fit(:,1),POS_fit(:,2),POS_fit(:,3),'or'); hold on;
h_rep = plot3(REP.pos_fit(:,1),REP.pos_fit(:,2),REP.pos_fit(:,3),'ok'); hold on;
try
set(gca,'xtick',REP.hypercube_limits(:,1)','ytick',REP.hypercube_limits(:,2)','ztick',REP.hypercube_limits(:,3)');
axis([min(REP.hypercube_limits(:,1)) max(REP.hypercube_limits(:,1)) ...
min(REP.hypercube_limits(:,2)) max(REP.hypercube_limits(:,2))]);
end
grid on; xlabel('f1'); ylabel('f2'); zlabel('f3');
drawnow;
axis square;
end
display(['Generation #0 - Repository size: ' num2str(size(REP.pos,1))]);
% Main MPSO loop
stopCondition = false;
while ~stopCondition
% Select leader
h = selectLeader(REP);
% Update speeds and positions
VEL = W.*VEL + C1*rand(Np,nVar).*(PBEST-POS) ...
+ C2*rand(Np,nVar).*(repmat(REP.pos(h,:),Np,1)-POS);
POS = POS + VEL;
% Perform mutation
POS = mutation(POS,gen,maxgen,Np,var_max,var_min,nVar,u_mut);
% Check boundaries
[POS,VEL] = checkBoundaries(POS,VEL,maxvel,var_max,var_min);
% Evaluate the population
POS_fit = fun(POS);
% Update the repository
REP = updateRepository(REP,POS,POS_fit,ngrid);
if(size(REP.pos,1)>Nr)
REP = deleteFromRepository(REP,size(REP.pos,1)-Nr,ngrid);
end
% Update the best positions found so far for each particle
pos_best = dominates(POS_fit, PBEST_fit);
best_pos = ~dominates(PBEST_fit, POS_fit);
best_pos(rand(Np,1)>=0.5) = 0;
if(sum(pos_best)>1)
PBEST_fit(pos_best,:) = POS_fit(pos_best,:);
PBEST(pos_best,:) = POS(pos_best,:);
end
if(sum(best_pos)>1)
PBEST_fit(best_pos,:) = POS_fit(best_pos,:);
PBEST(best_pos,:) = POS(best_pos,:);
end
% Plotting and verbose
if(size(POS_fit,2)==2)
figure(h_fig); delete(h_par); delete(h_rep);
h_par = plot(POS_fit(:,1),POS_fit(:,2),'or'); hold on;
h_rep = plot(REP.pos_fit(:,1),REP.pos_fit(:,2),'ok'); hold on;
try
set(gca,'xtick',REP.hypercube_limits(:,1)','ytick',REP.hypercube_limits(:,2)');
axis([min(REP.hypercube_limits(:,1)) max(REP.hypercube_limits(:,1)) ...
min(REP.hypercube_limits(:,2)) max(REP.hypercube_limits(:,2))]);
end
if(isfield(MultiObj,'truePF'))
try delete(h_pf); end
h_pf = plot(MultiObj.truePF(:,1),MultiObj.truePF(:,2),'.','color',0.8.*ones(1,3)); hold on;
end
grid on; xlabel('f1'); ylabel('f2');
drawnow;
axis square;
end
if(size(POS_fit,2)==3)
figure(h_fig); delete(h_par); delete(h_rep);
h_par = plot3(POS_fit(:,1),POS_fit(:,2),POS_fit(:,3),'or'); hold on;
h_rep = plot3(REP.pos_fit(:,1),REP.pos_fit(:,2),REP.pos_fit(:,3),'ok'); hold on;
try
set(gca,'xtick',REP.hypercube_limits(:,1)','ytick',REP.hypercube_limits(:,2)','ztick',REP.hypercube_limits(:,3)');
axis([min(REP.hypercube_limits(:,1)) max(REP.hypercube_limits(:,1)) ...
min(REP.hypercube_limits(:,2)) max(REP.hypercube_limits(:,2)) ...
min(REP.hypercube_limits(:,3)) max(REP.hypercube_limits(:,3))]);
end
if(isfield(MultiObj,'truePF'))
try delete(h_pf); end
h_pf = plot3(MultiObj.truePF(:,1),MultiObj.truePF(:,2),MultiObj.truePF(:,3),'.','color',0.8.*ones(1,3)); hold on;
end
grid on; xlabel('f1'); ylabel('f2'); zlabel('f3');
drawnow;
axis square;
end
display(['Generation #' num2str(gen) ' - Repository size: ' num2str(size(REP.pos,1))]);
% Update generation and check for termination
gen = gen + 1;
if(gen>maxgen), stopCondition = true; end
end
hold off;
end% Function that updates the repository given a new population and its
% fitness
function REP = updateRepository(REP,POS,POS_fit,ngrid)
% Domination between particles
DOMINATED = checkDomination(POS_fit);
REP.pos = [REP.pos; POS(~DOMINATED,:)];
REP.pos_fit= [REP.pos_fit; POS_fit(~DOMINATED,:)];
% Domination between nondominated particles and the last repository
DOMINATED = checkDomination(REP.pos_fit);
REP.pos_fit= REP.pos_fit(~DOMINATED,:);
REP.pos = REP.pos(~DOMINATED,:);
% Updating the grid
REP = updateGrid(REP,ngrid);
end% Function that corrects the positions and velocities of the particles that
% exceed the boundaries
function [POS,VEL] = checkBoundaries(POS,VEL,maxvel,var_max,var_min)
% Useful matrices
Np = size(POS,1);
MAXLIM = repmat(var_max(:)',Np,1);
MINLIM = repmat(var_min(:)',Np,1);
MAXVEL = repmat(maxvel(:)',Np,1);
MINVEL = repmat(-maxvel(:)',Np,1);
% Correct positions and velocities
VEL(VEL>MAXVEL) = MAXVEL(VEL>MAXVEL);
VEL(VEL<MINVEL) = MINVEL(VEL<MINVEL);
VEL(POS>MAXLIM) = (-1).*VEL(POS>MAXLIM);
POS(POS>MAXLIM) = MAXLIM(POS>MAXLIM);
VEL(POS<MINLIM) = (-1).*VEL(POS<MINLIM);
POS(POS<MINLIM) = MINLIM(POS<MINLIM);
end% Function for checking the domination between the population. It
% returns a vector that indicates if each particle is dominated (1) or not
function dom_vector = checkDomination(fitness)
Np = size(fitness,1);
dom_vector = zeros(Np,1);
all_perm = nchoosek(1:Np,2); % Possible permutations
all_perm = [all_perm; [all_perm(:,2) all_perm(:,1)]];
d = dominates(fitness(all_perm(:,1),:),fitness(all_perm(:,2),:));
dominated_particles = unique(all_perm(d==1,2));
dom_vector(dominated_particles) = 1;
end% Function that returns 1 if x dominates y and 0 otherwise
function d = dominates(x,y)
d = all(x<=y,2) & any(x<y,2);
end
