1.算法概述
粒子群优化算法(PSO)是一种进化计算技术(evolutionary computation)，1995 年由Eberhart 博士和kennedy 博士提出，源于对鸟群捕食的行为研究 。该算法最初是受到飞鸟集群活动的规律性启发，进而利用群体智能建立的一个简化模型。粒子群算法在对动物集群活动行为观察基础上，利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得最优解。
PSO 过程:
第1 步: 种群随机初始化。
第2 步: 对种群内的每一个个体计算适应值(fitness value) , 适应值与最优解的距离直接有关。
第3 步: 种群根据适应值进行复制。
第4 步: 如果终止条件满足, 则停止; 否则转到第2 步。
PSO 算法中每个优化问题的解都有是搜索空间中的一只鸟, 称为粒子。与其他进化计算技术不同的是群体中的每个粒子可以记忆自己到过的最优位置, 并能感知邻近群体已达到的最优位置, 每个粒子能够根据自身到过的最优位置和邻近群体已到过的最优位置来更新自己, 然后粒子们不断地追随当前的最优粒子在解空间搜索。
有M个作业（运输任务），S个代理人，k种运输方式，n个节点（城市）。每个作业都有时间限制，由第四方物流对各作业、各代理人、各种运输方式进行整合。各作业在任意节点上可由任意代理人进行代理，即在各节点可进行代理人之间的代理转换；各作业在任意节点之间只可选择一种运输方式，代理人在节点之间的运输能力不同，根据不同运量提供的代理价格折扣不同。模型为求解出各作业的运输路线及在各节点上选择由哪个代理人选择哪种运输方式。
成本分为四个部分：运输费用、代理人转换费用、运输方式转化费用和时间惩罚费用。
2.仿真效果预览
matlab2022a仿真结果如下：
3.核心MATLAB程序 `..............................................................
%PSO
%x(s,k,m,i,j) = 1 表示作业m在节点i和节点j之间由代理人s采用k种运输方式代理；否则x(s,k,m,i,j)=0；
%r(s,l,m,i) = 1 表示作业m在节点i由代理人s转换成代理人l；否则r(s,l,m,i)=0;
%R(k,v,m,i) = 1 表示作业m在节点i由k种运输方式转换为v种运输方式，否则R(k,v,m,i)=0;
%由于算法较为复杂，这里无法直接将所有因素考虑，这里采用分级优化，即对性能影响最大的因素进行优化，再给予优化结果进行次级因素优化
%确定路线
%确定路线
%确定路线
%初始化x，r，R,初始化的值是随便设置的
for i = 1:n
for j = 1:n
if d(i,j) ~= 0 & d(i,j) ~= F
x(:,:,:,i,j) = 1;
r(:,:,:,i) = 1;
R(:,:,:,i) = 1;
else
x(:,:,:,i,j) = 0;
r(:,:,:,i) = 0;
R(:,:,:,i) = 0;
end
end
end
All_cost = fitness(M,n,g,G,C,q,d,p,Z,T,LT,ET,R,r,x);
%下面开始PSO优化
itmax = 300;%进化代数，就是预设的迭代次数。
W(1) = 0.729;% 粒子先前速度保持。惯性权重
a(1) = 0.316;% 用于计算W。
c1 = 2; %认知部分 加速系数
c2 = 2; %社会部分 加速系数
xmax = 1;
xmin = 0;
ii = 1;
num_particle = 100;
D = size(d,1);
particle = zeros(2*num_particle,D,D,M,itmax);
particle(:,:,:,:,1) = xmin+(xmax-xmin)rand(2num_particle,D,D,M);
V(:,:,:,:,1) = round((xmin-xmax)+2*(xmax-xmin)rand(2num_particle,D,D,M));
fit = zeros(num_particle,itmax);% 用于存储粒子的适应值
pbest = zeros(2*num_particle,D,D,M,itmax); % 用于存储粒子的位置
x2 = zeros(g,G,M,n,n,2*num_particle);
for m = 1:M
for i = 1:n
for j = 1:n
for nn = 1 : 2*num_particle
x2(:,:,m,i,j,nn) = particle(nn,i,j,m,1);
end
end
end
end
x_tmp = zeros(g,G,M,n,n);
for nn = 1 : num_particle
x_tmp = x2(:,:,:,:,:,nn);
fit(nn,1) = fitness(M,n,g,G,C,q,d,p,Z,T,LT,ET,R,r,x_tmp);
end
%*********************************************************
pbest(:,:,:,:,1) = particle(:,:,:,:,1);
pbest_value(:,1) = fit(:,1); %个体最优值
[Cs,I] = min(pbest_value(:,1));
gbest_value(1) = Cs; % 群最优值
for i=1:num_particle
gbest(2i-1:2i,:,:,:,1)=particle(2I-1:2I,:,:,:,1); %群最优粒子位置
end
tmps = 0;
route = zeros(n,n,M,2*num_particle);
for ii=2:itmax
ii
V(:,:,:,:,ii) = 0.729V(:,:,:,:,ii-1)+c1rand*(pbest(:,:,:,:,ii-1)-particle(:,:,:,:,ii-1))+...
c2rand(gbest(:,:,:,:,ii-1)-particle(:,:,:,:,ii-1));
V(:,:,:,:,ii) = min(V(:,:,:,:,ii),xmax-xmin);
V(:,:,:,:,ii) = max(V(:,:,:,:,ii),xmin-xmax);
particle(:,:,:,:,ii) = particle(:,:,:,:,ii-1)+V(:,:,:,:,ii);
particle(:,:,:,:,ii) = min(particle(:,:,:,:,ii),xmax);
particle(:,:,:,:,ii) = max(particle(:,:,:,:,ii),xmin);
for m = 1:M
for i = 1:n
for j = 1:n
for nn = 1 : 2*num_particle
if d(i,j) > 0
x2(:,:,m,i,j,nn) = double(particle(nn,i,j,m,ii)>0.5);%对于优化结果，只取0或者1
else
x2(:,:,m,i,j,nn) = 0;%对于优化结果，只取0或者1
end
end
end
end
end
for m = 1:M
for i = 1:n
for j = 1:n
for nn = 1 : 2*num_particle
if d(i,j) > 0
route(i,j,m,nn) = particle(nn,i,j,m,ii);
else
route(i,j,m,nn) = 0;
end
end
end
end
end
for nn = 1 : num_particle
x_tmp = x2(:,:,:,:,:,nn);
fit(:,ii) = fitness(M,n,g,G,C,q,d,p,Z,T,LT,ET,R,r,x_tmp);
end
%下面更新 pbest and pbest_value
pbest_value(:,ii)=min(pbest_value(:,ii-1),fit(:,ii));
for i=1:num_particle
if pbest_value(i,ii) == fit(i,ii)
pbest(2i-1:2i,:,:,:,ii) = particle(2i-1:2i,:,:,:,ii);
else
pbest(2i-1:2i,:,:,:,ii) = pbest(2i-1:2i,:,:,:,ii-1);
end
end
........................................................
02-006m`
