MOPSO_自适应网格机制_matlab

%CMOPSO1 主程序针对测试函数1 clear all; clc; %---------初始化基本参数----------% format long; %long 输出结果是长型，short是短型 a=[0 0;6 7]; %搜索空间界限 for i=1:2 vmax(i)=0.25*(a(2,i)-a(1,i)); %速度的最大值经验为权值的0.25 end N=100; %群体规模 R=100; %外部档案大小 w=0.4; %权重 D=2; %待优化问题的维数 Tmax=100; %最大迭代次数 Py=zeros(N,2); %存放粒子的个体极值 A=zeros(N,2); %存放每代粒子的适应值 Xa=[]; Xrv=zeros(31,2);%初始化各个格子的边界 Nrv=zeros(30,30);%初始化各个格子所含粒子的数目 %---------初始化粒子速度和位置----------% for i=1:N %矩阵，N行D列 rand均匀分布的伪随机数；randn正态分布的随机数 for j=1:D x(i,j)=rand*(a(2,j)-a(1,j))+a(1,j); %随机初始化位置 v(i,j)=0; %初始化速度 end end %---------计算两个目标函数的适应值，初始化个体极值----------% for i=1:N %初始化个体极值 Py(i,1)=f1(x(i,1:D)); %计算目标f1的适应值 Py(i,2)=f2(x(i,1:D)); %计算目标f2的适应值 Xy(i,1:D)=x(i,1:D); %记录个体极值的坐标位置 A(i,1)=Py(i,1); %初始矩阵A A(i,2)=Py(i,2); Xa=[Xa;x(i,1:D)]; %初始矩阵A对应的位置坐标 end %----------A集合成员的支配关系比较，得到当前的非支配解集C----------% [C,Xc]=ParetoC(A,Xa); %调用函数得到C,和Xc %----------产生超子空间----------% B=C;Xb=Xc; %将非劣解集放入外部档案 [Nrv,Xrv]=hypercubes(B); %---------用轮盘赌法从B中选取一个粒子，将其位置值赋给Pg----------% [cp]=roulette(Nrv,Xrv,B); Pg=Xb(cp,1:D); %----------主循环，迭代过程----------% for t=1:Tmax %迭代次数 Xa=[];C=[];Xc=[];clear('A'); if ~mod(t,5) fprintf('current iter:%d\n',t) end for i=1:N %粒子位置和速度更新 v(i,1:D)=w*v(i,1:D)+rand*(Xy(i,1:D)-x(i,1:D))+rand*(Pg-x(i,1:D)); for j=1:D if v(i,j)>vmax(j) v(i,j)=vmax(j); else if v(i,j)<-vmax(j) v(i,j)=-vmax(j); end end end x(i,1:D)=x(i,1:D)+v(i,1:D); for j=1:D if(x(i,j)>a(2,j)) x(i,j)=a(2,j); v(i,j)=-v(i,j); end %搜索空间上限 if(x(i,j)<a(1,j)) x(i,j)=a(1,j); v(i,j)=-v(i,j); end %搜索空间下限 end %----------计算得到新的目标向量集A----------% A(i,1)=f1(x(i,1:D)); A(i,2)=f2(x(i,1:D)); Xa(i,:)=x(i,1:D); %----------突变操作----------% [Xa]=mutation(Xa,t,a,N,Tmax,i); %----------更新个体极值Py(i)----------% con1=0; %存放Ai受约束条件的个数 con2=0; %存放Py受约束条件的个数 if (Xa(i,1)/6+Xa(i,2)-13/2)>0 con1=con1+1; %计算Ai受约束条件的个数 end if (Xa(i,1)/2+Xa(i,2)-15/2)>0 con1=con1+1; end if (Xa(i,1)*5+Xa(i,2)-30)>0 con1=con1+1; end if (Xy(i,1)/6+Xy(i,2)-13/2)>0 con2=con2+1; %计算Py受约束条件的个数 end if (Xy(i,1)/2+Xy(i,2)-15/2)>0 con2=con2+1; end if (Xy(i,1)*5+Xy(i,2)-30)>0 con2=con2+1; end if con1<con2 %Py受约束条件个数大于Ai,Ai支配Py Py(i,1)=A(i,1); Py(i,2)=A(i,2); Xy(i,1:D)=x(i,1:D); elseif con1>con2 %Ai受约束条件个数大于Py,Py支配Ai Py(i,1)=Py(i,1);Py(i,2)=Py(i,2);Xy(i,1:D)=Xy(i,1:D); elseif ((A(i,1)<Py(i,1))&&(A(i,2)<=Py(i,2))) %Ai支配Py，则Py更新 Py(i,1)=A(i,1); Py(i,2)=A(i,2); Xy(i,1:D)=x(i,1:D); elseif ((A(i,1)<=Py(i,1))&&(A(i,2)<Py(i,2))) %Ai支配Py，则Py更新 Py(i,1)=A(i,1); Py(i,2)=A(i,2); Xy(i,1:D)=x(i,1:D); elseif ((Py(i,1)<A(i,1))&&(Py(i,2)<=A(i,2))) %Py支配Ai，则Py保持不变 Py(i,1)=Py(i,1);Py(i,2)=Py(i,2);Xy(i,1:D)=Xy(i,1:D); elseif ((Py(i,1)<=A(i,1))&&(Py(i,2)<A(i,2))) %Py支配Ai，则Py保持不变 Py(i,1)=Py(i,1);Py(i,2)=Py(i,2);Xy(i,1:D)=Xy(i,1:D); else %Ai和Py之间没有支配与被支配关系，Py取两者中的任意值 m=randint(1,1,[1,2]); if m==1 Py(i,1)=A(i,1); Py(i,2)=A(i,2); Xy(i,1:D)=x(i,1:D); end end end [C,Xc]=ParetoC(A,Xa);%判断A中成员的支配关系得到新的当前非劣解集C %----------更新B----------% C=[C;B]; Xc=[Xc;Xb]; [B,Xb]=ParetoC(C,Xc); [Nrv,Xrv]=hypercubes(B); if size(B,1)>R Bt=[];Xbt=[]; for i=1:R [cp]=roulette(Nrv,Xrv,B);%调用轮盘赌法选择函数选择粒子 Bt=[Bt;B(cp,:)]; Xbt=[Xbt;Xb(cp,:)]; [B m]=setdiff(B,Bt,'rows'); Xb=Xb(m,:); [Nrv,Xrv]=hypercubes(B); end B=Bt; Xb=Xbt; end %----------更新Pg----------% [Nrv,Xrv]=hypercubes(B); [cp]=roulette(Nrv,Xrv,B); Pg=Xb(cp,1:D); %---------出图，每20代显示一次----------% if mod(t,1)==0 figure(1) xt=0:0.01:40; y1=-xt/6+13/2; y2=-xt/2+15/2; y3=-5*xt+30; plot(Xa(:,1),Xa(:,2),'*',xt,y1,xt,y2,xt,y3) xlabel('x');ylabel('y'); grid on axis([0 7 0 7]); title(strcat('迭代',num2str(t),'代后的粒子位置分布')); figure(2) plot(Xb(:,1),Xb(:,2),'*',xt,y1,xt,y2,xt,y3) xlabel('x');ylabel('y'); grid on axis([0 7 0 7]); title(strcat('迭代',num2str(t),'次后的Pareto粒子位置分布')); figure(3) plot(B(:,1),B(:,2),'*') xlabel('F1');ylabel('F2'); grid on title(strcat('迭代',num2str(t),'次后的Pareto非劣解')); end end