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) if dnp(i,j)~=0 && dnp(i,j)~=firstmin && dnp(i,j) sparedis(i)=(firstmin+secondmin)/2; end end %%%%%%%%%%%%%%%%%-----------密集(稀疏)距离 ------------%%%%%%%%%%%%%%%%%%% function v = DominateRel(x,y,funcname,params) %判断x与y支配关系,返回1表示x支配y，返回-1表示y支配x,返回0表示互不支配 v=0; faix = GetFai(x,funcname,params); faiy = GetFai(y,funcname,params); if faix>faiy v=-1; elseif faiy>faix v=1; end if v~=0 %x、y构成违反约束支配关系 return; end gx = fitness(x,funcname,params);%x的目标向量 gy = fitness(y,funcname,params);%y的目标向量 len = length(gx); if sum(gx<=gy)==len%x的所有目标都比y小，x支配y v=1; elseif sum(gx>=gy)==len%y的所有目标都比x小，y支配x v=-1; end end %%%%%%%%%%%%%%%%%---------比较两粒子的相互支配关系-----------%%%%%%%%%%%%% function [gbest,nprule_out] = GetGlobalBest2(x,np,nprule) %按照一定原则随机取一个全局最优 r=size(np,2);%非劣解的行列,r:非劣解集个数，c：维数 %假定x被非劣解集np支配 IsDominated=true; Mdom=[];%支配x的集合 for i=1:r domiRel=DominateRel(np(i,:),x,funcname,params); if domiRel==1%np(i)支配x Mdom=[Mdom,i];%记下其在非劣解集中的编号 else IsDominated=false; break; end end nprule_out=nprule; if IsDominated%x被支配,从x的支配集中选出一个作为全局最优 intem=randi([1,length(Mdom)],1); gbest=np(Mdom(intem),:);nprule_out(Mdom(intem),1)=nprule_out(Mdom(intem),1)+1; else%x不被支配，从所有非劣解集中选出一个作为全局最优 nr=size(np,1); intem=randi([1,nr],1); gbest=np(intem,:);nprule_out(intem,1)=nprule_out(intem,1)+1; end end function [gbest,nprule_out] = GetGlobalBest(np,nprule,dnp_out) %随机取一个全局最优 r=size(np,1);%非劣解的行列 nprule_out=nprule; intem=1; if(randi([0,3],1)==0) if r==1 gbest = np(1,:); else sparedis = GetSpareDis(dnp_out); %[max1,max2] = find(max(max(dnp_out))==dnp_out); %intem=max1(randi([1,length(max1)],1)); n_max=find(max(sparedis)==sparedis); intem=n_max(randi([1,length(n_max)],1)); gbest = np(intem,:); end else %rr=randi([1,r],1);%随机取一个作为全局最优 %gbest = np(rr,:); tt=find(min(nprule(:,1))==nprule(:,1)); intem=tt(randi([1,length(tt)],1)); gbest = np(intem,:); end nprule_out(intem,1)=nprule_out(intem,1)+1; end %%%%%%%%%%%%%%%%%%%----------从部种群中找到全局最优-------%%%%%%%%%%%%%%% function [x,v,pbest,NP,NPRule,Dnp]=ReInit(x,v,pbest,NP,NPRule,Dnp,Nnp,D,lb,ub,unfitx,funcname,params) for i=1:10 for j=1:D x(i,j)=(ub(j)-lb(j))*rand; %随机初始化位置 v(i,j)=(ub(j)-lb(j))*rand*0.5; %随机初始化速度 pbest(i,j)=x(i,j); end end for i=1:10 faix = GetFai(x(i,:),funcname,params); if faix<=unfitx [NP,NPRule,Dnp] = compare(x(i,:),NP,NPRule,Dnp,Nnp,funcname,params); end end end %%%%%%%%%%%%%%%%%%%----------重新初始化部分粒子-------%%%%%%%%%%%%%%% function fv=fitness(x,funcname,params) %获得多目标的目标向量 fv fv=[]; switch upper(funcname) %------DTLZ1 %------ZDT1 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)); %------ZDT2 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); %------ZDT3 case 'ZDT3' n=length(x); gv=1+9*sum(x(2:n).^2)/(n-1); fv(1)=x(1); fv(2)=gv*(1-sqrt(x(1)/gv)-(x(1)/gv)*sin(10*pi*x(1))); %------ZDT4 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)); %------ZDT5 case 'ZDT5' %------ZDT4 case 'ZDT6' n=length(x); gv=1+9*(sum(x(2:n))/(n-1)).^0.25; fv(1)=1-exp(-4*x(1))*sin(6*pi*x(1)).^6; fv(2)=gv*(1-(fv(1)/gv).^2); %------CTP1 case 'CTP1' n=length(x); cx=41+sum(x(2:n).^2-10*cos(2*pi*x(2:n))); fv(1)=x(1); fv(2)=cx*exp(-fv(1)/cx); %------CTP2,CTP3,CTP4,CTP5,CTP6,CTP7 case {'CTP2','CTP3','CTP4','CTP5','CTP6','CTP7'} n=length(x); cx=41+sum(x(2:n).^2-10*cos(2*pi*x(2:n))); fv(1)=x(1); fv(2)=cx*(1-fv(1)/cx); %------CTP8,CONSTR case {'CTP8','CONSTR'} fv(1)=x(1); fv(2)=(1+x(2))/x(1); %------CTP9,SRN case {'CTP9','SRN'} fv(1)=(x(1)-2)^2+(x(2)-1)^2+2; fv(2)=9*x(1)-(x(2)-1)^2; %------CTP10,TNK,MOP-C4 case {'CTP10','TNK','MOP-C4'} fv(1)=x(1); fv(2)=x(2); %-------MOP-C1,BNH case {'MOP-C1','BNH'} fv(1)=4*x(1)^2+4*x(2)^2; fv(2)=(x(1)-5)^2+(x(2)-5)^2; %------MOP-C2,OSY case {'MOP-C2','OSY'} fv(1)=-(25*(x(1)-2)^2+(x(2)-2)^2+(x(3)-1)^2+(x(4)-4)^2+(x(5)-1)^2); fv(2)=x(1)^2+x(2)^2+x(3)^2+x(4)^2+x(5)^2+x(6)^2; %------MOP-C3 case 'MOP-C3' fv(1)=(x(1)-2)^2/2+(x(2)+1)^2/13+3; fv(2)=(x(1)+x(2)-3)^2/175+(2*x(2)-x(1))^2/17-13; fv(3)=(3*x(1)-2*x(2)+4)^2/8+(x(1)-x(2)+1)^2/27+15; end end function fai = GetFai(x,funcname,params) %构造约束函数 switch upper(funcname) %---DTLZ1,...,DTLZ4函数 case {'DTLZ1','DTLZ2','DTLZ3','DTLZ4'} fai=0; %---ZDT1,...,ZDT6函数 case {'ZDT1','ZDT2','ZDT3','ZDT4','ZDT5','ZDT6'} fai=0; %------CTP1 case 'CTP1' n=length(x); cx=41+sum(x(2:n).^2-10*cos(2*pi*x(2:n))); f1=x(1); f2=cx*exp(-f1/cx); g1=0.858*exp(-0.541*f1)-f2; g2=0.728*exp(-0.295*f1)-f2;