MOPSO.rar_MOPSO多目标_mopso_mopso matlab_多目标_多目标 优化

function [np,nprule,dnp,fv,goals,pbest] = ParticleSwarmOpt(fitness,unfitx,N,Nnp,cmax,cmin,w,M,D,lb,ub,vars) %待优化的目标函数:fitness %约束容忍度unfitx=0.01,逐步降到0 %内部种群(粒子数目)：N %外部种群(非劣解集):Nnp %学习因子1:cmax %学习因子2:cmin %惯性权重:w %最大迭代次数：M %问题的维数：D %目标函数取最小值时的自变量值:xm %目标函数的最小值:fv %迭代次数:cv format long; NP=[];%非劣解集 Dnp=[];%非劣解集距离 pbest=[];%自身最优解 gbest=[];%全局最优解 goals=[];%记下粒子目标变化 NPRule=[];%非劣解集参数 vmax=[]; vmin=[]; Loc=vars{1}; Mvol=vars{2}; [lr,lc]=size(Loc); %----初始化种群的个体--------///////第1步 x(1,:)=lb; % x(1,:)=ones(1,D)*2; % x(2,:)=ub; v(1,:)=(ub-lb).*rand([1,D])*0.5; % v(2,:)=(ub-lb).*rand([1,D])*0.5; for i=2:N % for j=1:D % x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %随机初始化位置 % v(i,j)=(ub(j)-lb(j))*rand*0.5; %随机初始化速度 % end for li=1:lr for lj=1:lc j=Loc(li,lj); if(j>0) sx=0; %求这一列x的和 for k=1:li-1 if(Loc(k,lj)>0) sx=sx+x(i,Loc(k,lj)); end end x(i,j)=round(lb(j)+(ub(j)-lb(j)-sx)*rand);%随机初始化位置 v(i,j)=(ub(j)-lb(j))*rand*0.5; %随机初始化速度 end end end end %----计算目标向量---------- %---速度控制 vmax=(ub-lb)*0.5; vmin= -vmax; %-----求出初始NP-----------////////第2步 NP(1,:)=x(1,:);%第一个默认加入NP NPRule=[0,0,0]; Dnp(1,1)=0; for i=2:N faix = GetFai(x(i,:),vars); if faix<=unfitx [NP,NPRule,Dnp] = compare(x(i,:),NP,NPRule,Dnp,Nnp,vars); end end %-----初始自身最好位置------///////第3步 pbest = x; %-----在确定每个粒子所对就的目标方格-------//第4步 %------进入主要循环，按照公式依次迭代------------ for t=1:M if mod(t,100)==0 unfitx = 0.01-0.01*(t+200)/M; if unfitx <0 unfitx =0 ; end % [x,v,pbest,NP,NPRule,Dnp]=ReInit(x,v,pbest,NP,NPRule,Dnp,Nnp,D,lb,ub,unfitx); end c = cmax - (cmax - cmin)*t/M; w1=w-(w-0.3)*t/M; %c = cmax; %c = cmax - (cmax - cmin)*mod(t,51)/50; %w1=w-(w-0.3)*mod(t,51)/50; %-----获得全局最优-------/////第5步 %if mod(t,3)==1 %[gbest,NPRule] = GetGlobalBest(NP,NPRule,Dnp); %end for i=1:N %-------------------更新粒子的位置和速度----------////第6步 %v(i,:)=w*v(i,:)+cmin*rand*(pbest(i,:)-x(i,:))+cmax*rand*(gbest-x(i,:)); %[gbest,NPRule] = GetGlobalBest2(x(i,:),NP,NPRule);%%%%%12-17 [gbest,NPRule] = GetGlobalBest(NP,NPRule,Dnp); %w1=Inf; %while(w1>1.2 || w1<0) % w1=0.8+randn*0.8; %end 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 %-------------------采取措施，避免粒子飞出空间----------////第7步 %速度位置钳制 %sita=0.5; % for j=1:D % if x(i,j)+v(i,j)>ub(j) % v(i,j)=rand*(ub(j)-x(i,j)); % x(i,j)=x(i,j)+v(i,j); % elseif x(i,j)+v(i,j)<lb(j) % v(i,j)=rand*(x(i,j)-lb(j)); % x(i,j)=x(i,j)-v(i,j); % else % x(i,j)=x(i,j)+v(i,j); % end % if x(i,j)>ub(j) || x(i,j)<lb(j) % x(i,j) % end % end x(i,:)=x(i,:)+v(i,:); for j=1:D if x(i,j)>ub(j) if(randi([0,0],1)==0) x(i,j)=ub(j); v(i,j)=-v(i,j); else x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %随机初始化位置 v(i,j)=(ub(j)-lb(j))*rand*0.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); else x(i,j)=lb(j)+(ub(j)-lb(j))*rand; %随机初始化位置 v(i,j)=(ub(j)-lb(j))*rand*0.5; end end end x(i,:)=round(x(i,:)); %------------------每个粒子的目标向量-----------------////第8步 goals(t,i,:)=fitness(x(i,:),vars); %----------------调整自身---------------------------//第10步 domiRel = DominateRel(pbest(i,:),x(i,:),vars);%x,y的支配关系 if domiRel==1%pbest支配新解 continue; else if domiRel==-1%新解支配pbest pbest(i,:) = x(i,:); elseif(randi([0,1],1)==0)%新解与pbest互相不支配 pbest(i,:) = x(i,:); end %-----------------对NP进行更新和维护-----------------//第9步 faix = GetFai(x(i,:),vars); if faix<=unfitx [NP,NPRule,Dnp] = compare(x(i,:),NP,NPRule,Dnp,Nnp,vars); end end end end np = NP;%非劣解 nprule=NPRule; dnp = Dnp;%非劣解之间的距离 [r,c]=size(np); for i=1:r fv(i,:)=fitness(np(i,:),vars); end function [np_out,nprule_out,dnp_out] = compare(x,np,nprule,dnp,nnp,vars) %np:现有非劣解 %x:需要比较的量 Nnp = nnp;%非劣解集空间 [r,c]=size(np);%非劣解的行列 np_out=np;%非劣解复本 nprule_out = nprule; dnp_out = dnp;%非劣解集点之间距离 if r==0 return; end for i=r:-1:1 domiRel=DominateRel(x,np(i,:),vars); if domiRel==1 %NP(i)被x支配 np_out(i,:)=[];%非劣解剔除该解 nprule_out(i,:)=[]; dnp_out(i,:)=[]; if ~isempty(dnp_out) dnp_out(:,i)=[]; end elseif domiRel==-1 %x被NP(i)支配,返回不再比较 return; end end [r1,c1]=size(np_out);%现有非劣解的行列 np_out(r1+1,:)=x;%与所有非支配集粒子比较均占优或不可比较，则NP中加入x faix = GetFai(x,vars); if faix > 0 nprule_out(r1+1,:)=[0,faix,0] else nprule_out(r1+1,:)=[0,0,0]; end if r1==0 dnp_out=0; end for j=1:r1 dnp_out(r1+1,j)=GetDistance(np_out(j,:),x); dnp_out(j,r1+1)=dnp_out(r1+1,j); end if r1>=Nnp %未达到非劣解种群极限 %---------移除密集距离最小的一个------- densedis = GetDenseDis(dnp_out); n_min = find(min(densedis)==densedis);%找出密度距离最小的一个 tempIndex = randi([1,length(n_min)],1); np_out(n_min(tempIndex),:)=[];%非劣解剔除该解 nprule_out(n_min(tempIndex),:)=[]; dnp_out(n_min(tempIndex),:)=[]; if ~isempty(dnp_out) dnp_out(:,n_min(tempIndex))=[]; end end function dis=GetDistance(x,y) %求两向量之间的距离 g=x-y; dis=sum(sum(g.^2)); %gx=fitness(x); %gy=fitness(y); %gxy=(gx-gy).^2; %dis=sqrt(sum(gxy(:))); function densedis = GetDenseDis(dnp) %密集距离 [r,c] = size(dnp); 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 % densedis(i)=(firstmin+secondmin)/2; densedis(i)=firstmin; end function sparedis = GetSpareDis(dnp) %稀疏距离 [r,c] = size(dnp); 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 function v = DominateRel(x,y,vars) %判断x与y支配关系,返回1表示x支配y，返回-1表示y支配x,返回0表示互不支配 v=0; faix = GetFai(x,vars); faiy = GetFai(y,vars); if faix>faiy v=-1; elseif faiy>faix v=1; end if v~=0 %x、y构成违反约束支配关系 return; end gx = fitness(x,vars);%x的目标向量 gy = fitness(y,vars);%