用粒子群算法求解多目标优化问题的Pareto解

粒子群算法程序ticD=10;%粒子群中粒子的个数%w=0.729;%w为惯性因子wmin=1.2;wmax=1.4;c1=1.49445;%正常数，成为加速因子c2=1.49445;%正常数，成为加速因子Loop_max=50;%最大迭代次数%初始化粒子群for i=1:DX(i)=rand(1)*(-5-7)+7;V(i)=1;f1(i)=X(i)^2;f2(i)=(X(i)-2)^2;endLoop=1;%迭代计数器while Loop<=Loop_max%循环终止条件%对粒子群中的每个粒子进行评价for i=1:Dk1=find(1==Xv(i,:));%找出第一辆车配送的城市编号nb1=size(k1,2);%计算第一辆车配送城市的个数if nb1>0%判断第一辆车配送城市个数是否大于0，如果大于0则a1=[Xr(i,k1(:))];%找出第一辆车配送城市顺序号b1=sort(a1);%对找出第一辆车的顺序号进行排序G1(i)=0;%初始化第一辆车的配送量k51=[];am=[];for j1=1:nb1am=find(b1(j1)==Xr(i,:));k51(j1)=intersect(k1,am);%计算第一辆车配送城市的顺序号G1(i)=G1(i)+g(k51(j1)+1);%计算第一辆车的配送量endk61=[];k61=[0,k51,0];%定义第一辆车的配送路径L1(i)=0;%初始化第一辆车的配送路径长度for k11=1:nb1+1L1(i)=L1(i)+Distance(k61(k11)+1,k61(k11+1)+1);%计算第一辆车的配送路径长度endelse%如果第一辆车配送的城市个数不大于0则G1(i)=0;%第一辆车的配送量设为0L1(i)=0;%第一辆车的配送路径长度设为0endk2=find(2==Xv(i,:));%找出第二辆车配送的城市编号nb2=size(k2,2);%计算第二辆车配送城市的个数if nb2>0%判断第二辆车配送城市个数是否大于0，如果大于0则a2=[Xr(i,k2(:))];%找出第二辆车配送城市的顺序号b2=sort(a2);%对找出的第二辆车的顺序号进行排序G2(i)=0;%初始化第二辆车的配送量k52=[];bm=[];for j2=1:nb2bm=find(b2(j2)==Xr(i,:));k52(j2)=intersect(k2,bm);%计算第二辆车配送城市的顺序号G2(i)=G2(i)+g(k52(j2)+1);%计算第二辆车的配送量endk62=[];k62=[0,k52,0];%定义第二辆车的配送路径L2(i)=0;%初始化第二辆车的配送路径长度for k22=1:nb2+1L2(i)=L2(i)+Distance(k62(k22)+1,k62(k22+1)+1);%计算第二辆车的路径长度endelse%如果第二辆车配送的城市个数不大于0则G2(i)=0;%第二辆车的配送量设为0L2(i)=0;%第二辆车的配送路径长度设为0endk3=find(3==Xv(i,:));%找出第三辆车配送的城市编号nb3=size(k3,2);%计算第三辆车配送城市的个数if nb3>0%判断第三辆车配送城市个数是否大于0，如果大于0则a3=[Xr(i,k3(:))];%找出第三辆车配送城市的顺序号b3=sort(a3);%对找出的第三辆车的顺序号进行排序G3(i)=0;%初始化第三辆车的配送量k53=[];cm=[];for j3=1:nb3cm=find(b3(j3)==Xr(i,:));k53(j3)=intersect(k3,cm);%计算第三辆车配送城市的顺序号G3(i)=G3(i)+g(k53(j3)+1);%计算第三辆车的配送量endk63=[];k63=[0,k53,0];%定义第三辆车的配送路径L3(i)=0;%初始化第三辆车的配送路径长度for k33=1:nb3+1L3(i)=L3(i)+Distance(k63(k33)+1,k63(k33+1)+1);%计算第三辆车的路径长度endelse%如果第三辆车配送的城市个数不大于0则G3(i)=0;%第三辆车的配送量设为0L3(i)=0;%第三辆车的配送路径长度设为0endL(i)=0;%初始化每个粒子对应的配送方案总路径长度L(i)=L1(i)+L2(i)+L3(i);%计算每个粒子对应的配送方案总路径长度if L(i)<Lg&&G1(i)<Q&&G2(i)<Q&&G3(i)<Q%如果第i个粒子的总路径长度优于历史最优粒子并且满足车辆容量要求Xvg(:)=Xv(i,:);%将粒子i设为历史最优粒子Xrg(:)=Xr(i,:);%将粒子i设为历史最优粒子Lg=L(i);%将粒子i的总路径长度设为最优粒子对应的配送方案的总路径长度elseXvg(:)=Xvg(:);%最优粒子保持不变Xrg(:)=Xrg(:);%最优粒子保持不变Lg=Lg;%最优粒子所对应的配送方案的总路径长度也不变endLimin(i)=100000;%初始化每个粒子代表的配送方案的历史最优总路径长度if L(i)<Limin(i)%如果本次循环得到的总路径长度优于粒子i历史最优总路径长度Limin(i)=L(i);%更新本次循环得到的总路径长度为粒子i的历史最优路径长度Xvl(i,:)=Xv(i,:);%更新本次得到的粒子i为i粒子的历史最优位置Xrl(i,:)=Xr(i,:); %更新本次得到的粒子i为i粒子的历史最优位置else%否则，保持粒子i的历史最优位置及历史最优路径长度不变Limin(i)=LL(i);Xvl(i,:)=Xv1(i,:);Xrl(i,:)=Xr1(i,:);endend%记录本次循环得到的所有粒子的位置for i=1:Dfor j=1:NXv1(i,j)=Xvl(i,j);%记录本次循环得到的所有粒子的位置Xr1(i,j)=Xrl(i,j);%记录本次循环得到的所有离子的位置endendLL(i)=0;%初始化每个粒子的历史最优路径总长度for i=1:DLL(i)=Limin(i);%对每个粒子的历史最优路径总长度进行赋值end%对粒子群中每个粒子进行迭代w=wmin+(wmax-wmin)*exp((-Loop)/(Loop_max-Loop));for i=1:Dfor j=1:NVv(i,j)=w*Vv(i,j)+c1*rand(1)*(Xvl(i,j)-Xv(i,j))+c2*rand(1)*(Xvg(1,j)-Xv(i,j));%计算位置变化率Vr(i,j)=w*Vr(i,j)+c1*rand(1)*(Xrl(i,j)-Xr(i,j))+c2*rand(1)*(Xrg(1,j)-Xr(i,j));%计算位置变化率%Vv(i,j)和Vr(i,j)进行上下限的限制if Vv(i,j)>K-1Vv(i,j)=K-1;elseif Vv(i,j)<1-KVv(i,j)=1-K;elseVv(i,j)=Vv(i,j);endendendfor i=1:Dfor j=1:NXv(i,j)=ceil(Xv(i,j)+Vv(i,j));%更新位置坐标%对Xv(i,j)进行上下限的限制if Xv(i,j)>KXv(i,j)=K;elseif Xv(i,j)<1Xv(i,j)=1;elseXv(i,j)=Xv(i,j);endXr(i,j)=Xr(i,j)+Vr(i,j);%更新位置坐标endendLoop=Loop+1;endXvg%输出粒子群中的最优粒子Xrg%输出粒子群中的最优粒子Lg%输出最优粒子所代表方案的总路径长度Loop%输出迭代的次数%计算最优粒子所代表的配送方案k1=find(1==Xvg(:));%找出第一辆车配送的城市编号k1=k1';nb1=size(k1,2);%计算第一辆车配送城市的个数if nb1>0%判断第一辆车配送城市个数是否大于0，如果大于0则a1=[Xrg(k1(:))];%找出第一辆车配送城市顺序号b1=sort(a1);%对找出第一辆车的顺序号进行排序G1=0;%初始化第一辆车的配送量k51=[];am=[];for j1=1:nb1am=find(b1(j1)==Xrg(:));k51(j1)=intersect(k1,am);%计算第一辆车配送城市的顺序号G1=G1+g(k51(j1)+1);%计算第一辆车的配送量endk61=[];k61=[0,k51,0];%定义第一辆车的配送路径L1=0;%初始化第一辆车的配送路径长度for k11=1:nb1+1L1=L1+Distance(k61(k11)+1,k61(k11+1)+1);%计算第一辆车的配送路径长度endelse%如果第一辆车配送的城市个数不大于0则G1=0;%第一辆车的配送量设为0L1=0;%第一辆车的配送路径长度设为0endk2=find(2==Xvg(:));%找出第二辆车配送的城市编号k2=k2';nb2=size(k2,2);%计算第二辆车配送城市的个数if nb2>0%判断第二辆车配送城市个数是否大于0，如果大于0则a2=[Xrg(k2(:))];%找出第二辆车配送城市的顺序号b2=sort(a2);%对找出的第二辆车的顺序号进行排序G2=0;%初始化第二辆车的配送量k52=[];bm=[];for j2=1:nb2bm=find(b2(j2)==Xrg(:));k52(j2)=intersect(k2,bm);%计算第二辆车配送城市的顺序号G2=G2+g(k52(j2)+1);%计算第二辆车的配送量endk62=[];k62=[0,k52,0];%定义第二辆车的配送路径L2=0;%初始化第二辆车的配送路径长度for k22=1:nb2+1L2=L2+Distance(k62(k22)+1,k62(k22+1)+1);%计算第二辆车的路径长度endelse%如果第二辆车配送的城市个数不大于0则G2=0;%第二辆车的配送量设为0L2=0;%第二辆车的配送路径长度设为0endk3=find(3==Xvg(:));%找出第三辆车配送的城市编号k3=k3';nb3=size(k3,2);%计算第三辆车配送城市的个数if nb3>0%判断第三辆车配送城市个数是否大于0，如果大于0则a3=[Xrg(k3(:))];%找出第三辆车配送城市的顺序号b3=sort(a3);%对找出的第三辆车的顺序号进行排序G3=0;%初始化第三辆车的配送量k53=[];cm=[];for j3=1:nb3cm=find(b3(j3)==Xrg(:));k53(j3)=intersect(k3,cm);%计算第三辆车配送城市的顺序号G3=G3+g(k53(j3)+1);%计算第三辆车的配送量endk63=[];k63=[0,k53,0];%定义第三辆车的配送路径L3=0;%初始化第三辆车的配送路径长度for k33=1:nb3+1L3=L3+Distance(k63(k33)+1,k63(k33+1)+1);%计算第三辆车的路径长度endelse%如果第三辆车配送的城市个数不大于0则G3=0;%第三辆车的配送量设为0L3=0;%第三辆车的配送路径长度设为0endk61k62k63x=City(:,1);y=City(:,2);%对各个城市进行顺序标号max_text={'0','1','2','3','4','5','6','7',};text(x+1,y+1,max_text)%画出最优粒子所代表的配送方案路径for i=1:nb1+2short1(i)=k61(i)+1;endfor i=1:nb2+2short2(i)=k62(i)+1;endfor i=1:nb3+2short3(i)=k63(i)+1;endline(x(short1),y(short1),'Marker','o')line(x(short2),y(short2),'Marker','o')line(x(short3),y(short3),'Marker','o')toc%计算程序的运行时间Time=num2str(toc)clear allticK=3;%车辆数D=200;%粒子群中粒子的个数Q=1;%每辆车的容量%w=0.729;%w为惯性因子wmin=1.2;wmax=1.4;c1=1.49445;%正常数，成为加速因子c2=1.49445;%正常数，成为加速因子Loop_max=50;%最大迭代次数%初始化城市坐标City=[18,54;22,60;58,69;71,71;83,46;91,38;24,42;18,40];n=size(City,1);%城市个数,包含中心仓库N=n-1;%发货点任务数for i=1:nfor j=1:nDistance(i,j)=sqrt((City(i,1)-City(j,1))^2+(City(i,2)-City(j,2))^2);%各城市节点之间的距离矩阵endendg=[0,0.89,0.14,0.28,0.33,0.21,0.41,0.57];%各发货点的货运量%初始化粒子群for i=1:Dfor j=1:NXv(i,j)=randi(K,1);%初始化粒子群中粒子的位置Vv(i,j)=randi(2*K-1,1)-K;%初始化粒子群中粒子的位置变化率Vr(i,j)=randi(2*N-1,1)-N;%初始化粒子群中离子的位置变化率Xvl(i,j)=Xv(i,j);%初始化粒子群中每个粒子的最优位置endendfor i=1:Da=randperm(N);for j=1:NXr(i,j)=a(j);%初始化粒子群中粒子的位置Xrl(i,j)=Xr(i,j);%初始化粒子群中每个粒子的最优位置endendLg=100000;%初始化最优粒子对应的配送方案的总路径长度Xvg=ones(1,N);%粒子群中最优的粒子Xrg=ones(1,N);%粒子群中最优的粒子Loop=1;%迭代计数器while Loop<=Loop_max%循环终止条件%对粒子群中的每个粒子进行评价for i=1:Dk1=find(1==Xv(i,:));%找出第一辆车配送的城市编号nb1=size(k1,2);%计算第一辆车配送城市的个数if nb1>0%判断第一辆车配送城市个数是否大于0，如果大于0则a1=[Xr(i,k1(:))];%找出第一辆车配送城市顺序号b1=sort(a1);%对找出第一辆车的顺序号进行排序G1(i)=0;%初始化第一辆车的配送量k51=[];am=[];for j1=1:nb1am=find(b1(j1)==Xr(i,:));k51(j1)=intersect(k1,am);%计算第一辆车配送城市的顺序号G1(i)=G1(i)+g(k51(j1)+1);%计算第一辆车的配送量endk61=[];k61=[0,k51,0];%定义第一辆车的配送路径L1(i)=0;%初始化第一辆车的配送路径长度for k11=1:nb1+1L1(i)=L1(i)+Distance(k61(k11)+1,k61(k11+1)+1);%计算第一辆车的配送路径长度endelse%如果第一辆车配送的城市个数不大于0则G1(i)=0;%第一辆车的配送量设为0L1(i)=0;%第一辆车的配送路径长度设为0endk2=find(2==Xv(i,:));%找出第二辆车配送的城市编号nb2=size(k2,2);%计算第二辆车配送城市的个数if nb2>0%判断第二辆车配送城市个数是否大于0，如果大于0则a2=[Xr(i,k2(:))];%找出第二辆车配送城市的顺序号b2=sort(a2);%对找出的第二辆车的顺序号进行排序G2(i)=0;%初始化第二辆车的配送量k52=[];bm=[];for j2=1:nb2bm=find(b2(j2)==Xr(i,:));k52(j2)=intersect(k2,bm);%计算第二辆车配送城市的顺序号G2(i)=G2(i)+g(k52(j2)+1);%计算第二辆车的配送量endk62=[];k62=[0,k52,0];%定义第二辆车的配送路径L2(i)=0;%初始化第二辆车的配送路径长度for k22=1:nb2+1L2(i)=L2(i)+Distance(k62(k22)+1,k62(k22+1)+1);%计算第二辆车的路径长度endelse%如果第二辆车配送的城市个数不大于0则G2(i)=0;%第二辆车的配送量设为0L2(i)=0;%第二辆车的配送路径长度设为0endk3=find(3==Xv(i,:));%找出第三辆车配送的城市编号nb3=size(k3,2);%计算第三辆车配送城市的个数if nb3>0%判断第三辆车配送城市个数是否大于0，如果大于0则a3=[Xr(i,k3(:))];%找出第三辆车配送城市的顺序号b3=sort(a3);%对找出的第三辆车的顺序号进行排序G3(i)=0;%初始化第三辆车的配送量k53=[];cm=[];for j3=1:nb3cm=find(b3(j3)==Xr(i,:));k53(j3)=intersect(k3,cm);%计算第三辆车配送城市的顺序号G3(i)=G3(i)+g(k53(j3)+1);%计算第三辆车的配送量endk63=[];k63=[0,k53,0];%定义第三辆车的配送路径L3(i)=0;%初始化第三辆车的配送路径长度for k33=1:nb3+1L3(i)=L3(i)+Distance(k63(k33)+1,k63(k33+1)+1);%计算第三辆车的路径长度endelse%如果第三辆车配送的城市个数不大于0则G3(i)=0;%第三辆车的配送量设为0L3(i)=0;%第三辆车的配送路径长度设为0endL(i)=0;%初始化每个粒子对应的配送方案总路径长度L(i)=L1(i)+L2(i)+L3(i);%计算每个粒子对应的配送方案总路径长度if L(i)<Lg&&G1(i)<Q&&G2(i)<Q&&G3(i)<Q%如果第i个粒子的总路径长度优于历史最优粒子并且满足车辆容量要求Xvg(:)=Xv(i,:);%将粒子i设为历史最优粒子Xrg(:)=Xr(i,:);%将粒子i设为历史最优粒子Lg=L(i);%将粒子i的总路径长度设为最优粒子对应的配送方案的总路径长度elseXvg(:)=Xvg(:);%最优粒子保持不变Xrg(:)=Xrg(:);%最优粒子保持不变Lg=Lg;%最优粒子所对应的配送方案的总路径长度也不变endLimin(i)=100000;%初始化每个粒子代表的配送方案的历史最优总路径长度if L(i)<Limin(i)%如果本次循环得到的总路径长度优于粒子i历史最优总路径长度Limin(i)=L(i);%更新本次循环得到的总路径长度为粒子i的历史最优路径长度Xvl(i,:)=Xv(i,:);%更新本次得到的粒子i为i粒子的历史最优位置Xrl(i,:)=Xr(i,:); %更新本次得到的粒子i为i粒子的历史最优位置else%否则，保持粒子i的历史最优位置及历史最优路径长度不变Limin(i)=LL(i);Xvl(i,:)=Xv1(i,:);Xrl(i,:)=Xr1(i,:);endend%记录本次循环得到的所有粒子的位置for i=1:Dfor j=1:NXv1(i,j)=Xvl(i,j);%记录本次循环得到的所有粒子的位置Xr1(i,j)=Xrl(i,j);%记录本次循环得到的所有离子的位置endendLL(i)=0;%初始化每个粒子的历史最优路径总长度for i=1:DLL(i)=Limin(i);%对每个粒子的历史最优路径总长度进行赋值end%对粒子群中每个粒子进行迭代w=wmin+(wmax-wmin)*exp((-Loop)/(Loop_max-Loop));for i=1:Dfor j=1:NVv(i,j)=w*Vv(i,j)+c1*rand(1)*(Xvl(i,j)-Xv(i,j))+c2*rand(1)*(Xvg(1,j)-Xv(i,j));%计算位置变化率Vr(i,j)=w*Vr(i,j)+c1*rand(1)*(Xrl(i,j)-Xr(i,j))+c2*rand(1)*(Xrg(1,j)-Xr(i,j));%计算位置变化率%Vv(i,j)和Vr(i,j)进行上下限的限制if Vv(i,j)>K-1Vv(i,j)=K-1;elseif Vv(i,j)<1-KVv(i,j)=1-K;elseVv(i,j)=Vv(i,j);endendendfor i=1:Dfor j=1:NXv(i,j)=ceil(Xv(i,j)+Vv(i,j));%更新位置坐标%对Xv(i,j)进行上下限的限制if Xv(i,j)>KXv(i,j)=K;elseif Xv(i,j)<1Xv(i,j)=1;elseXv(i,j)=Xv(i,j);endXr(i,j)=Xr(i,j)+Vr(i,j);%更新位置坐标endendLoop=Loop+1;endXvg%输出粒子群中的最优粒子Xrg%输出粒子群中的最优粒子Lg%输出最优粒子所代表方案的总路径长度Loop%输出迭代的次数%计算最优粒子所代表的配送方案k1=find(1==Xvg(:));%找出第一辆车配送的城市编号k1=k1';nb1=size(k1,2);%计算第一辆车配送城市的个数if nb1>0%判断第一辆车配送城市个数是否大于0，如果大于0则a1=[Xrg(k1(:))];%找出第一辆车配送城市顺序号b1=sort(a1);%对找出第一辆车的顺序号进行排序G1=0;%初始化第一辆车的配送量k51=[];am=[];for j1=1:nb1am=find(b1(j1)==Xrg(:));k51(j1)=intersect(k1,am);%计算第一辆车配送城市的顺序号G1=G1+g(k51(j1)+1);%计算第一辆车的配送量endk61=[];k61=[0,k51,0];%定义第一辆车的配送路径L1=0;%初始化第一辆车的配送路径长度for k11=1:nb1+1L1=L1+Distance(k61(k11)+1,k61(k11+1)+1);%计算第一辆车的配送路径长度endelse%如果第一辆车配送的城市个数不大于0则G1=0;%第一辆车的配送量设为0L1=0;%第一辆车的配送路径长度设为0endk2=find(2==Xvg(:));%找出第二辆车配送的城市编号k2=k2';nb2=size(k2,2);%计算第二辆车配送城市的个数if nb2>0%判断第二辆车配送城市个数是否大于0，如果大于0则a2=[Xrg(k2(:))];%找出第二辆车配送城市的顺序号b2=sort(a2);%对找出的第二辆车的顺序号进行排序G2=0;%初始化第二辆车的配送量k52=[];bm=[];for j2=1:nb2bm=find(b2(j2)==Xrg(:));k52(j2)=intersect(k2,bm);%计算第二辆车配送城市的顺序号G2=G2+g(k52(j2)+1);%计算第二辆车的配送量endk62=[];k62=[0,k52,0];%定义第二辆车的配送路径L2=0;%初始化第二辆车的配送路径长度for k22=1:nb2+1L2=L2+Distance(k62(k22)+1,k62(k22+1)+1);%计算第二辆车的路径长度endelse%如果第二辆车配送的城市个数不大于0则G2=0;%第二辆车的配送量设为0L2=0;%第二辆车的配送路径长度设为0endk3=find(3==Xvg(:));%找出第三辆车配送的城市编号k3=k3';nb3=size(k3,2);%计算第三辆车配送城市的个数if nb3>0%判断第三辆车配送城市个数是否大于0，如果大于0则a3=[Xrg(k3(:))];%找出第三辆车配送城市的顺序号b3=sort(a3);%对找出的第三辆车的顺序号进行排序G3=0;%初始化第三辆车的配送量k53=[];cm=[];for j3=1:nb3cm=find(b3(j3)==Xrg(:));k53(j3)=intersect(k3,cm);%计算第三辆车配送城市的顺序号G3=G3+g(k53(j3)+1);%计算第三辆车的配送量endk63=[];k63=[0,k53,0];%定义第三辆车的配送路径L3=0;%初始化第三辆车的配送路径长度for k33=1:nb3+1L3=L3+Distance(k63(k33)+1,k63(k33+1)+1);%计算第三辆车的路径长度endelse%如果第三辆车配送的城市个数不大于0则G3=0;%第三辆车的配送量设为0L3=0;%第三辆车的配送路径长度设为0endk61k62k63x=City(:,1);y=City(:,2);%对各个城市进行顺序标号max_text={'0','1','2','3','4','5','6','7',};text(x+1,y+1,max_text)%画出最优粒子所代表的配送方案路径for i=1:nb1+2short1(i)=k61(i)+1;endfor i=1:nb2+2short2(i)=k62(i)+1;endfor i=1:nb3+2short3(i)=k63(i)+1;endline(x(short1),y(short1),'Marker','o') line(x(short2),y(short2),'Marker','o') line(x(short3),y(short3),'Marker','o') toc%计算程序的运行时间Time=num2str(toc)clear all。