function [count,ps,pf]=MMGPE(Th,func_name,VRmin,VRmax,n_obj,Particle_Number,Max_Gen)
% MO_Ring_PSO_SCD: A multi-objective particle swarm optimization using ring topology for solving multimodal multi-objective optimization problems
% Dimension: n_var --- dimensions of decision space
% n_obj --- dimensions of objective space
%% Input:
% Dimension Description
% func_name 1 x length function name the name of test function
% VRmin 1 x n_var low bound of decision variable
% VRmax 1 x n_var up bound of decision variable
% n_obj 1 x 1 dimensions of objective space
% Particle_Number 1 x 1 population size
% Max_Gen 1 x 1 maximum generations
%% Output:
% Description
% ps Pareto set
% pf Pareto front
%% Reference and Contact
% Reference: [1]Caitong Yue, Boyang Qu and Jing Liang, "A Multi-objective Particle Swarm Optimizer Using Ring Topology for Solving Multimodal Multi-objective Problems", IEEE Transactions on Evolutionary Computation, 2017, DOI 10.1109/TEVC.2017.2754271.
% [2]Jing Liang, Caitong Yue, and Boyang Qu, “ Multimodal multi-objective optimization: A preliminary study”, IEEE Congress on Evolutionary Computation 2016, pp. 2454-2461, 2016.
% Contact: For any questions, please feel free to send email to zzuyuecaitong@163.com.
%% Initialize parameters
n_var=size(VRmin,2); %Obtain the dimensions of decision space
Max_FES=Max_Gen*Particle_Number; %Maximum fitness evaluations
n_PBA=5; %Maximum size of PBA. The algorithm will perform better without the size limit of PBA. But it will time consuming.
n_NBA=3*n_PBA; %Maximum size of NBA
cc=[2.05 2.05]; %Acceleration constants
iwt=0.7298; %Inertia weight
count=zeros(500,3);
%% Initialize particles' positions and velocities
mv=0.5*(VRmax-VRmin);
VRmin=repmat(VRmin,Particle_Number,1);
VRmax=repmat(VRmax,Particle_Number,1);
Vmin=repmat(-mv,Particle_Number,1);
Vmax=-Vmin;
pos=VRmin+(VRmax-VRmin).*rand(Particle_Number,n_var); %initialize the positions of the particles
vel=Vmin+2.*Vmax.*rand(Particle_Number,n_var); %initialize the velocities of the particles
Vrmin=(VRmin+abs(VRmin));
Vrmax=(VRmax+abs(VRmin));
Pos=pos+abs(VRmin);
OriginX(1)={Pos}; % 存储第一代
%% Evaluate the population
fitness=zeros(Particle_Number,n_obj);
for i=1:Particle_Number
fitness(i,:)=feval(func_name,pos(i,:));
end
fitcount=Particle_Number; % count the number of fitness evaluations
particle=[pos,fitness]; %put positions and velocities in one matrix
%% Initialize personal best archive PBA and Neighborhood best archive NBA
row_of_cell=ones(1,Particle_Number); % the number of row in each cell
col_of_cell=size(particle,2); % the number of column in each cell
PBA=mat2cell(particle,row_of_cell,col_of_cell);
NBA=PBA;
for i=2:3
%% Update NBA
for j=1:Particle_Number
% Ring topology. Particles exchange information with its closed neighbors
if j==1
tempNBA=PBA{Particle_Number,:};
tempNBA=[tempNBA;PBA{1,:}];
tempNBA=[tempNBA;PBA{2,:}];
elseif j==Particle_Number
tempNBA=PBA{Particle_Number-1,:};
tempNBA=[tempNBA;PBA{Particle_Number,:}];
tempNBA=[tempNBA;PBA{1,:}];
else
tempNBA=PBA{j-1,:};
tempNBA=[tempNBA;PBA{j,:}];
tempNBA=[tempNBA;PBA{j+1,:}];
end
NBA_j=NBA{j,1};
tempNBA=[tempNBA;NBA_j];
tempNBA=non_domination_scd_sort(tempNBA(:,1:n_var+n_obj), n_obj, n_var);
% Optional operation. Limit the size of NBA to save time.
if size(tempNBA,1)>n_NBA
NBA{j,1}=tempNBA(1:n_NBA,:);
else
NBA{j,1}=tempNBA;
end
end
for k=1:Particle_Number
% Choose the first particle in PBA_k as pbest
PBA_k=PBA{k,1}; %PBA_k contains the history positions of particle_k
pbest=PBA_k(1,:); %Choose the first one
% Choose the first particle in NBA_k as nbest
NBA_k=NBA{k,:}; %NBA_k contains the history positions of particle_k's neighborhood
nbest=NBA_k(1,:);
% Update velocities according to Eq.(5)
vel(k,:)=iwt.*vel(k,:)+cc(1).*rand(1,n_var).*(pbest(1,1:n_var)-pos(k,:))+cc(2).*rand(1,n_var).*(nbest(1,1:n_var)-pos(k,:));
% Make sure that velocities are in the setting bounds.
vel(k,:)=(vel(k,:)>mv).*mv+(vel(k,:)<=mv).*vel(k,:);
vel(k,:)=(vel(k,:)<(-mv)).*(-mv)+(vel(k,:)>=(-mv)).*vel(k,:);
% Update positions according to Eq.(4)
pos(k,:)=pos(k,:)+vel(k,:);
% Make sure that positions are in the setting bounds.
pos(k,:)=((pos(k,:)>=VRmin(1,:))&(pos(k,:)<=VRmax(1,:))).*pos(k,:)...
+(pos(k,:)<VRmin(1,:)).*(VRmin(1,:)+0.25.*(VRmax(1,:)-VRmin(1,:)).*rand(1,n_var))+(pos(k,:)>VRmax(1,:)).*(VRmax(1,:)-0.25.*(VRmax(1,:)-VRmin(1,:)).*rand(1,n_var));
% Evaluate the population
fitness(k,:)=feval(func_name,pos(k,1:n_var));
fitcount=fitcount+1;
particle(k,1:n_var+n_obj)=[pos(k,:),fitness(k,:)];
%% Update PBA
PBA_k=[PBA_k(:,1:n_var+n_obj);particle(k,:)];
PBA_k = non_domination_scd_sort(PBA_k(:,1:n_var+n_obj), n_obj, n_var);
% Optional operation. Limit the size of PBA to save time.
if size(PBA_k,1)>n_PBA
PBA{k,1}=PBA_k(1:n_PBA,:);
else
PBA{k,1}=PBA_k;
end
PBA_k=PBA{k,1};
C(k,:)=PBA_k(1,:);
end
D= C(:,1:n_var);
Pos=D+abs(VRmin);
OriginX(i)={Pos};
end
for i=4:(Max_Gen)
for j=1:Particle_Number
% Ring topology. Particles exchange information with its closed neighbors
if j==1
tempNBA=PBA{Particle_Number,:};
tempNBA=[tempNBA;PBA{1,:}];
tempNBA=[tempNBA;PBA{2,:}];
elseif j==Particle_Number
tempNBA=PBA{Particle_Number-1,:};
tempNBA=[tempNBA;PBA{Particle_Number,:}];
tempNBA=[tempNBA;PBA{1,:}];
else
tempNBA=PBA{j-1,:};
tempNBA=[tempNBA;PBA{j,:}];
tempNBA=[tempNBA;PBA{j+1,:}];
end
NBA_j=NBA{j,1};
tempNBA=[tempNBA;NBA_j];
tempNBA=non_domination_scd_sort(tempNBA(:,1:n_var+n_obj), n_obj, n_var);
% Optional operation. Limit the size of NBA to save time.
if size(tempNBA,1)>n_NBA
NBA{j,1}=tempNBA(1:n_NBA,:);
else
NBA{j,1}=tempNBA;
end
end
t=0.01-(0.05/(Max_Gen))*(i-(Max_Gen));
aa=randperm(Particle_Number);
bb=randperm(Particle_Number);
cc=randperm(Particle_Number);
for k=1:Particle_Number
NBA_ii=NBA{k,:}; %NBA_k contains the history positions of particle_k's neighborhood
nbest=NBA_ii(1,:);
Nbest=nbest(1,1:n_var)+abs(VRmin(k,:));
for jj=1: n_var
X0=[OriginX{1,1}(aa(k),jj),OriginX{1,2}(bb(k),jj),OriginX{1,3}(cc(k),jj)];
if abs(max(X0)-min(X0))<Th(jj)
count(i,1)=count(i,1)+1;
X0(3)= Nbest(jj);
Pos(k,jj)=X0(3)+t*2*abs(max(X0)-min(X0))*(rand-0.5);
elseif (max(X0)~=min(X0))&&(abs(X0(1)-X0(2))<Th(jj)||abs(X0(1)-X0(3))<Th(jj)||abs(X0(2)-X0(3))<Th(jj))
count(i,2)=count(i,2)+1;
% X0(3)= Nbest(jj);
Pos(k,jj)=(4*X0(3)+X0(2)-2*X0(1))./3;
else
count(i,3)=count(i,3)+1;
a=2*(X0(2)-X0(3))/(X0(2)+X0(3));