%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% An improved multi-objective particle swarm optimization algorithm including multiple test functions
% Modifying some parameters in the program will better solve some functions.
%% Principal function
function []=main()
Varin = load('mydata.mat');%The parameters of the constraint function are imported.
ZDTFV=cell(1,50); %// Create cell arrays
ZDT=zeros(1,50); %//0 matrix
funcname = 'ZDT1';
times = 10;%It is equivalent to running ten programs independently.
M = 100;%MOPSO The number of iterations in
for i=1:times %//Loop 10 times, do the following iterations
tic; %//Timing begins
%[np,nprule,dnp,fv,goals,pbest]=ParticleSwarmOpt(funcname,100,200,2.0,1.0,0.5,M,30,Varin);%--ZDT3 zeros(1,9)-5->zeros(1,29)
[np,nprule,dnp,fv,goals,pbest]=ParticleSwarmOpt(funcname,100,200,2.0,1.0,0.4,M,10,[0,zeros(1,9)],[1,zeros(1,9)+5],Varin);%--ZDT4
elapsedTime=toc; %//End of Timing
ZDTFV(i)={fv};
ZDT(i)=elapsedTime;
display(strcat('The number of iterations is',num2str(i)));
end
zdt1fv=cell2mat(ZDTFV');
for i =1:times
display(strcat(i,':'));
disp(ZDT(i));%That's the time.
end
disp(zdt1fv);
disp('Next, the updated fitness value is:');
zdt1fv=GetLeastFunctionValue(zdt1fv);
disp(zdt1fv);
figure(9)
plot(zdt1fv(:,1),zdt1fv(:,2),'k*');
%Set the font form and size of the coordinate axis as follows
xlabel('$f_1$','interpreter','latex','FontSize',25);
ylabel('$f_2$','interpreter','latex','FontSize',25);
set(gca,'FontName','Times New Roman','FontSize',25)%Set the size of the coordinate axis font and scale,get current axes Returns the handle value of the current axis object
%if(strcmp(funcname,'ZDT3'))
axis([0 1 0 1]);
%else
% axis([20555 20790 190.4 190.7]);
%end
% The convergence is obtained as follows
% Take the truth evenly first Pareto The point on the optimal solution, and then the value of its two objective functions
p_true=0:0.002:1;
pf_true1=p_true;
pf_true2=1-p_true.^2;
r=size(zdt1fv,1);
for i=1:r%For every non-inferior solution,
for j=1:501
d(i,j)=sqrt((zdt1fv(i,1)-pf_true1(j))^2+(zdt1fv(i,2)-pf_true2(j))^2);
end
end
%Here's the right one d The minimum value of each row in the table is the first. i The Minimum Distance Between Points and Points
for i=1:r
dmin(i)=min(d(i,:));
end
Cmean=mean(dmin);
%Cfangcha=var(dmin)% The square of the difference between the number and the mean divided by the number-1，This is sample variance.
disp('Sample mean:');
disp(Cmean);
Cvariance=var(dmin,1);% The square of the difference between the number and the mean divided by the number,This is the definition of difference in mathematics.
disp('Sample variance:');
disp(Cvariance)
% Diversity delta
% First pair zdt1fv Sort in ascending order of the first column, i.e. from left to right in abscissa (the first target value)
zdt1fv=sortrows(zdt1fv);%End in ascending order of the first column, as follows df(Distance from the leading edge of the extreme solution on the left)and dl(Distance from the leading edge of the extreme solution on the rightmost side)
df=sqrt((zdt1fv(1,1)-0)^2+(zdt1fv(1,2)-1)^2);
r=length(zdt1fv);
dl=sqrt((zdt1fv(r,1)-1)^2+(zdt1fv(r,2)-0)^2);
for i=1:r-1
%The first i Number One and Number One i+1 The distance between the frontiers of solutions is d(i)
c(i)=sqrt((zdt1fv(i,1)-zdt1fv(i+1,1))^2+(zdt1fv(i,2)-zdt1fv(i+1,2))^2);
end
%The following requirements d The mean of
meanNum=mean(c);
%Calculation of substitution formula delta Value
%The part of the summation number
sum=0;
for i=1:r-1
sum=sum+abs(c(i)-meanNum);
end
delta=(df+dl+sum)/(df+dl+(r-1)*meanNum);
disp('The diversity is as follows:');
disp(delta);
end
%% MOPSO Function Definition
%function [np,nprule,dnp,fv,goals,pbest] = ParticleSwarmOpt(funcname,N,Nnp,cmax,cmin,w,M,D,Varin)
function [np,nprule,dnp,fv,goals,pbest] = ParticleSwarmOpt(funcname,N,Nnp,cmax,cmin,w,M,D,lb,ub,Varin)
%The objective function to be optimized:fitness
%Power plant constraint function: fitness2
%Internal population(Number of particles): N
%External population(non-inferior set):Nnp
%Fitness parameters
%Learning factor 1:cmax
%Learning factor 2:cmin
%Inertial weight:w
%Maximum number of iterations: M
%The dimension of the problem: D
%The Value of Independent Variables when the Object Function is Minimum:xm
%Minimum value of objective function:fv
%Number of iterations:cv
%Non-inferior inspection:flag
%Adaptive parameter:unifit:1->0.1
format long;
unifit = 1;
flag = 0;
NP=[];%non-inferior set
Dnp=[];%Distance of Non-inferior Solution Set
params = struct('isfmopso',true,'istargetdis',false,'stopatborder',true);%ZTD2->isfmopso(false->true)A change ZTD3 The problem should be true
%x0=lb+(ub-lb).*rand([1,D]);
%T=size(fitness(x0,funcname),2);
T = 2;
goals=zeros(M,N,T);%Note down N A particle M Subiteration T Dimensional Target Change
% %----Individuals Initializing the Population--------/////// Step 1///////////////////////////////////
% %x(1,:)=x0;
% %v(1,:)=(ub-lb).*rand([1,D])*0.5;
x = zeros(N,D);
v = zeros(N,D);
% for i=1:N
% for j=1:D
% x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %Random Initialization Location
% v(i,j)=(ub(j)-lb(j))*rand*0.5; %Random initialization speed
% end
% end
% %----Computation of target vectors----------
% %---speed control
% %vmax=(ub-lb)*0.5;
%vmin = -vmax;
%----Individuals Initializing the Population--------/////// Step 1///////////////////////////////////
for i=1:N
for j=1:D
x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %Random Initialization Location
v(i,j)=(ub(j)-lb(j))*rand*0.5; %Random initialization speed
end
end
%----Computation of target vectors----------
%---speed control
vmax=(ub-lb)*0.5;
vmin= -vmax;
%-----Find the Initial NP-----------//////// Step 2///////////////////////////////////
NP(1,:)=x(1,:);%The first default addition (non-inferior solution set), the particle is a row of uncertain columns, the number of columns represents the number of decision variables, that is, the dimension of the problem.
NPRule=[0,0,0];%Non-inferior set parameters
Dnp(1,1)=0;
for i=2:N
[NP,NPRule,Dnp,flag] = compare(flag,x(i,:),NP,NPRule,Dnp,Nnp,funcname,params,Varin);
end
%-----Best Place of Initial Self------/////// Step 3////////////////////////////////////
pbest = x;%Self-optimum solution
%-----In determining the target square for each particle-------//Step 4///////////////////////////
%------Enter the main cycle and iterate according to the formula.------------
for t=1:M
c = cmax - (cmax - cmin)*t/M;
w1=w-(w-0.3)*t/M;
for i=1:N
%-----Achieving Global Optimality-------///// Step 5/////////////////////////////////////////////
[gbest,NPRule] = GetGlobalBest(NP,NPRule,Dnp);
v(i,:)=w1*v(i,:)+c*rand*(pbest(i,:)-x(i,:))+c*rand*(gbest-x(i,:));
for j=1:D
if v(i,j)>vmax(j)
v(i,j)=vmax(j);
elseif v(i,j)<vmin(j)
v(i,j)=vmin(j);
end
end
x(i,:)=x(i,:)+v(i,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------Take measures to prevent particles from flying out of space----------//// Step 7/////////////
%Velocity position clamping
if(params.stopatborder)%Particles stochastically reside at the boundary
if x(i,1)>ub(1)
x(i,1)=ub(1);
v(i,1)=-v(i,1);
end
if x(i,1)<lb(1)
x(i,1)=lb(1)+(ub(1)-lb(1))*rand; %Random Initialization Location
v(i,1)=(ub(1)-lb(1))*rand*0.5;
end
for j=2:D
if x(i,j)>ub(j)
if(randi([0,2],1)==0)%Change 0->1
x(i,j)=ub(j);
v(i,j)=-v(i,j);
else
x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %Random Initialization Location
v(i,j)=(ub(j)-lb(j))*rand*1.5;
end
end
if x(i,j)<lb(j)
if(randi([0,0],1)==0)
x(i,j)=lb(j);
v(i,j)=-v(i,j)*unifit;
if(unifit>0.2)
unifit = unifit -0.1;
end
else
x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %Random Initialization Location
v(i,j)=(ub(j)-lb(j))*rand*0.5;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------------Target vector for each particle-----------------//Step 8///////////////
goals(t,i,:)=fitness(x(i,:),funcname,Varin);
%----------------Adjust oneself---------------------------//Step 9/////////////////
domiRel = DominateRel(pbest(i,:),x(i,:),funcname,params,Varin);%x,y Dominance
if domiRel==1%pbest Dominant New Solution
continue;
else
if domiRel==-1%New Dissolution of Dominance pbest
pbest(i,:) = x(i,:);
elseif(rand*2<1)%New Solutions and pbest Mutual non-domination
pbest(i,:) = x(i,:);
end
%-----------------Yes NP Update and maintain-----------------//Step 10////////////////
[NP,NPRule,Dnp,flag] = compare(flag,x(i,:),NP,NPRule,Dnp,Nnp,funcname,params,Varin);
if flag==1%In order to overcome the problem that the algorithm is easy to fall into local optimum, a non-inferior checking and monitoring mechanism is introduced.
[NP,flag,x,v] = fresh(NP,flag,x,v);
end
end
end
end
np = NP;%Non-inferior solution
nprule=NPRule;
dnp = Dnp;%The Distance Between Non-inferior Solutions
r=size(np,1);
fv=zeros(r,T);
for i=1:r
fv(i,:)=fitness(np(i,:),funcname,Varin);
end
end
%%%%%%%%%%%%%%%--------------End of main function--------------%%%%%%%%%%%%%%%%%
%% Maintaining Particles to External Populations
function [np_out,nprule_out,dnp_out,flag] = compare(flag,x,np,nprule,dnp,nnp,funcname,params,Varin)
%np:Existing non-inferior solutions
%x:Quantities to be compared
Nnp = nnp;%Non-inferior solution set space
r=size(np,1);%Number of Non-inferior Solutions
np_out=np;%Non-inferior solution duplicate
nprule_out = nprule;
dnp_out = dnp;%The Distance between Sets of Non-inferior Solutions
if r==0
return;
end
for i=r:-1:1
domiRel=DominateRel(x,np(i,:),funcname,params,Varin);
if domiRel==1 %NP(i)cover x Control
np_out(i,:)=[];%Elimination of non-inferior solutions
nprule_out(i,:)=[];
dnp_out(i,:)=[];
if ~isempty(dnp_out)
dnp_out(:,i)=[];
end
elseif domiRel==-1 %x cover NP(i)Control,Return no longer compares
return;
end
end
r1=size(np_out,1);%Existing ranks of non-inferior solutions
np_out(r1+1,:)=x;%If all non-dominant set particles are dominant or incomparable, then NP Joining in x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nprule_out(r1+1,:)=[0,0,0];
if r1==0
dnp_out=0;
end
for j=1:r1
dnp_out(r1+1,j)=GetDistance(np_out(j,:),x,funcname,params);
dnp_out(j,r1+1)=dnp_out(r1+1,j);
end
if r1>=Nnp %Reaching the Population Limit of Non-inferior Solutions
%---------Remove one of the smallest dense distances-------
densedis = GetDenseDis(dnp_out);
n_min = find(min(densedis)==densedis);%Find out the smallest density distance
tempIndex = randi([1,length(n_min)],1);
if min(densedis)==0
flag = 1;
end
np_out(n_min(tempIndex),:)=[];%Elimination of non-inferior solutions
nprule_out(n_min(tempIndex),:)=[];
dnp_out(n_min(tempIndex),:)=[];
if ~isempty(dnp_out)
dnp_out(:,n_min(tempIndex))=[];
end
end
end
%% Find the distance between two vectors
function dis=GetDistance(x,y,funcname,params)
if(params.istargetdis)
gx=fitness(x,funcname,Varin);
gy=fitness(y,funcname,Varin);
gxy=(gx-gy).^2;
dis=sqrt(sum(gxy(:)));
else
g=x-y;
dis=sum(sum(g.^2));
end
end
%% Dense distance (nearest distance)
function densedis = GetDenseDis(dnp)
[r,c] = size(dnp);
densedis=zeros(1,r);
for i=1:r
firstmin=Inf;
for j=1:c
if dnp(i,j)~=0 && dnp(i,j)<firstmin
firstmin = dnp(i,j);
end
end
densedis(i)=firstmin;
end
end
%% Sparse Distance (Second Close Distance)
function sparedis = GetSpareDis(dnp)
[r,c] = size(dnp);
sparedis=zeros(1,r);
for i=1:r
firstmin=Inf;
secondmin=Inf;
for j=1:c
if dnp(i,j)~=0 && dnp(i,j)<firstmin
firstmin = dnp(i,j);
end
if dnp(i,j)~=0 && dnp(i,j)~=firstmin && dnp(i,j)<secondmin
secondmin = dnp(i,j);
end
end
sparedis(i)=(firstmin+secondmin)/2;
end
end
%% Comparing the Dominance of Two Particles
function v = DominateRel(x,y,funcname,~,Varin)
%judge x and y Dominance,Return 1 representation x Control y，Return-1 Express y Control x,Returning 0 means no control over each other
v=0;
gx = fitness(x,funcname,Varin);%x Target vector
gy = fitness(y,funcname,Varin);%y Target vector
len = length(gx);
if sum(gx<=gy)==len%x All goals are better than y Small, x Control y
v=1;
elseif sum(gx>=gy)==len%y All goals are better than x Small, y Control x
v=-1;
end
end
%% Random global optimization
function [gbest,nprule_out] = GetGlobalBest(np,nprule,dnp_out)
r=size(np,1);%Ranks of non-inferior solutions
nprule_out=nprule;
intem=1;
if(round(rand)==0)
if r==1
gbest = np(1,:);
else
sparedis = GetSpareDis(dnp_out);
if(round(rand)==0)
n_max=find(max(sparedis)==sparedis);
intem=n_max(round(rand*(length(n_max)-1)+1));
gbest = np(intem,:);
else %Adding an operation to find the smallest
sparedis = GetSpareDis(dnp_out);
n_min=find(min(sparedis)==sparedis);
intem=n_min(round(rand*(length(n_min)-1)+1));
gbest = np(intem,:);
end
end
else
tt=find(min(nprule(:,1))==nprule(:,1)); %Choose one randomly as the global optimum and look at the lowest number of times that have been selected.
intem=tt(round(rand*(length(tt)-1)+1));
gbest = np(intem,:);
end
nprule_out(intem,1)=nprule_out(intem,1)+1;%For the selected particle, the nprule(The same number of rows np，The column number is 3) corresponds to this
%The value of the first column of the row is added by 1, because this column records the number of times the particle has been selected for future reference as to whether it should be selected as the global extremum again.
end
%% Power Plant Constraint Function
function fv = fitness2(x,~,Varin)
res1=0;
res2=0;
for i=1:6
res1 = res1+Varin.a(i)*x(i)*x(i)+Varin.b(i)*x(i)+Varin.c(i);
res2 = res2+Varin.Ea(i)*x(i)*x(i)+Varin.Eb(i)*x(i)+Varin.Ec(i)+Varin.Ed(i)*exp(Varin.Ee(i)*x(i));
end
fv(1) = res1;
fv(2) = res2;
end
%% ZDT1,ZDT2,ZDT3 Test function
function fv=fitness(x,funcname,~)
%Obtaining the Object Vector of Multiple Objects fv
fv=[];
switch upper(funcname)
case 'ZDT1'
n=length(x);
gv=1+9*sum(x(2:n))/(n-1);
fv(1)=x(1);
fv(2)=gv*(1-sqrt(x(1)/gv));
case 'ZDT2'
n=length(x);
gv=1+9*sum(x(2:n))/(n-1);
fv(1)=x(1);
fv(2)=gv*(1-(x(1)/gv).^2);
case 'ZDT3'
n=length(x);
gv=1+9*sum(x(2:n))/(n-1);
fv(1)=x(1);
fv(2)=gv*(1-sqrt(x(1)/gv)-(x(1)/gv)*sin(10*pi*x(1)));
case 'ZDT4'
n=length(x);
gv=1+10*(n-1)+sum(x(2:n).^2-10*cos(4*pi*x(2:n)));
fv(1)=x(1);
fv(2)=gv*(1-sqrt(x(1)./gv));
end
end
%% Excluding Non-dominant Sets from External Population
function fvout=GetLeastFunctionValue(fvin)
fvout=fvin;
n=size(fvout,1);
i=1;
while(i<=n)
j=i+1;
isdominated=false;%Judgment Conditions
while(j<=n)
a=fvout(i,:);b=fvout(j,:);
if((a(1)<b(1)&&a(2)<=b(2))||(a(1)<=b(1)&&a(2)<b(2)))%b cover a Dominated
fvout(j,:)=[];n=n-1;
else
if((b(1)<a(1)&&b(2)<=a(2))||(b(1)<=a(1)&&b(2)<a(2)))%a cover b Dominated
isdominated=true;
end
j=j+1;
end
end
if isdominated
fvout(i,:)=[];n=n-1;
else
i=i+1;
end
end
end
%% flag For 1,It shows that the result may fall into local optimum, so add 25.%Fresh particles
function [NP,flag,x,v] = fresh(NP,~,x,v)
r=size(NP,1);
flag =1;
for i=2:r
if(randi([0,3],1)==0)
for j=1:D
x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %Random Initialization Location
v(i,j)=(ub(j)-lb(j))*rand*0.5; %Random initialization speed
end
end
end
end