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

粒子群算法程序

tic

D=10;%粒子群中粒子的个数

%w=0.729;%w为惯性因子

wmin=1.2;

wmax=1.4;

c1=1.49445;%正常数，成为加速因子

c2=1.49445;%正常数，成为加速因子

Loop_max=50;%最大迭代次数

%初始化粒子群

for i=1:D

X(i)=rand(1)*(-5-7)+7;

V(i)=1;

f1(i)=X(i)^2;

f2(i)=(X(i)-2)^2;

end

Loop=1;%迭代计数器

while Loop<=Loop_max%循环终止条件

%对粒子群中的每个粒子进行评价

for i=1:D

k1=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:nb1

am=find(b1(j1)==Xr(i,:));

k51(j1)=intersect(k1,am);%计算第一辆车配送城市的顺序号

G1(i)=G1(i)+g(k51(j1)+1);%计算第一辆车的配送量

end

k61=[];

k61=[0,k51,0];%定义第一辆车的配送路径

L1(i)=0;%初始化第一辆车的配送路径长度

for k11=1:nb1+1

L1(i)=L1(i)+Distance(k61(k11)+1,k61(k11+1)+1);%计算第一辆车的配送路径长度end

else%如果第一辆车配送的城市个数不大于0则

G1(i)=0;%第一辆车的配送量设为0

L1(i)=0;%第一辆车的配送路径长度设为0

end

k2=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:nb2

bm=find(b2(j2)==Xr(i,:));

k52(j2)=intersect(k2,bm);%计算第二辆车配送城市的顺序号

G2(i)=G2(i)+g(k52(j2)+1);%计算第二辆车的配送量

end

k62=[];

k62=[0,k52,0];%定义第二辆车的配送路径

L2(i)=0;%初始化第二辆车的配送路径长度

for k22=1:nb2+1

L2(i)=L2(i)+Distance(k62(k22)+1,k62(k22+1)+1);%计算第二辆车的路径长度end

else%如果第二辆车配送的城市个数不大于0则

G2(i)=0;%第二辆车的配送量设为0

L2(i)=0;%第二辆车的配送路径长度设为0

end

k3=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:nb3

cm=find(b3(j3)==Xr(i,:));

k53(j3)=intersect(k3,cm);%计算第三辆车配送城市的顺序号

G3(i)=G3(i)+g(k53(j3)+1);%计算第三辆车的配送量

end

k63=[];

k63=[0,k53,0];%定义第三辆车的配送路径

L3(i)=0;%初始化第三辆车的配送路径长度

for k33=1:nb3+1

L3(i)=L3(i)+Distance(k63(k33)+1,k63(k33+1)+1);%计算第三辆车的路径长度end

else%如果第三辆车配送的城市个数不大于0则

G3(i)=0;%第三辆车的配送量设为0

L3(i)=0;%第三辆车的配送路径长度设为0

end

L(i)=0;%初始化每个粒子对应的配送方案总路径长度

L(i)=L1(i)+L2(i)+L3(i);%计算每个粒子对应的配送方案总路径长度

if L(i)

Xvg(:)=Xv(i,:);%将粒子i设为历史最优粒子

Xrg(:)=Xr(i,:);%将粒子i设为历史最优粒子

Lg=L(i);%将粒子i的总路径长度设为最优粒子对应的配送方案的总路径长度

else

Xvg(:)=Xvg(:);%最优粒子保持不变

Xrg(:)=Xrg(:);%最优粒子保持不变

Lg=Lg;%最优粒子所对应的配送方案的总路径长度也不变

end

Limin(i)=100000;%初始化每个粒子代表的配送方案的历史最优总路径长度

if L(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,:);

end

end

%记录本次循环得到的所有粒子的位置

for i=1:D

for j=1:N

Xv1(i,j)=Xvl(i,j);%记录本次循环得到的所有粒子的位置

Xr1(i,j)=Xrl(i,j);%记录本次循环得到的所有离子的位置

end

end

LL(i)=0;%初始化每个粒子的历史最优路径总长度

for i=1:D

LL(i)=Limin(i);%对每个粒子的历史最优路径总长度进行赋值

end

%对粒子群中每个粒子进行迭代

w=wmin+(wmax-wmin)*exp((-Loop)/(Loop_max-Loop));

for i=1:D

for j=1:N

Vv(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-1