1.算法概述
粒子群优化 (PSO)算法是通过模拟鸟群觅食过程中的迁徙和群聚行为而提出的一种基于群体智能的全局随机搜索算法。PSO是将群体(swarm)中的个体看作是在D维搜索空间中没有质量和体积的粒子(particle)，每个粒子以一定的速度在解空间运动，并向自身历史最佳位置pbest和邻域历史最佳位置pbest聚集，实现对候选解的进化。
PSO从这种模型中得到启示并用于解决优化问题。PSO中，每个优化问题的潜在解都是搜索空间中的一只鸟，称之为粒子。所有的粒子都有一个由被优化的函数决定的适值( fitness value)，每个粒子还有一个速度决定它们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。
PSO初始化为一群随机粒子，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个极值来更新自己；第一个就是粒子本身所找到的最优解，这个解称为个体极值；另一个极值是整个种群目前找到的最优解，这个极值是全局极值。另外也可以不用整个种群而只是用其中一部分作为粒子的邻居，那么在所有邻居中的极值就是局部极值。
在找到这两个最优值时，粒子根据如下的公式(5)和(6)来更新自己的速度和位置：
式(5)右边由三部分组成，第一部分为“惯性”或“动量”部分，反映了粒子的运动”习惯”，代表粒子有维持自己先前速度的趋势；第二部分为”认知”部分，反映了粒子对自身历史经验的记忆或回忆，代表粒子有向自身历史最佳位置逼近的趋势；第三部分为”社会”部分，反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或邻域历史最佳位置逼近的趋势。
其中，定义的适应度函数表达式如下：
即成本函数分为四个部分：运输费用、代理人转换费用、运输方式转化费用和时间惩罚费用。
2.仿真效果预览
matlab2022a仿真结果如下：
对于代理人，结果原来应该给出了，运行完，查看MATLAB的指令窗口，如下所示：
对比两种方式，仿真对比结果如下所示（VIEW2）：
3.核心MATLAB程序
%d(i,j)表示节点i到节点j之间的运输距离,0表示两点不可达到
F = 0;
d = func_dis(F);%调用距离函数
%q(m)表示作业m的运量；
q = 1e3*[4 7 6 3 5 7 4 7];
%不同代理人不同的运输方式的单位费用
w1 = 2;
w2 = 4;
w3 = 3;
w4 = 2;
cost = func_kcost1(w1,w1,w3,w4,g,G);%初始化价格，实际在公式中，通过输入运输量来确定具体的价格
%C(s,k,i,j)表示节点i到节点j由代理人s选择k种运输方式的单位运输费用，因为存在折扣问题，所以此变量为单调递减；
for i = 1:n
for j = 1:n
for s = 1:g
for k = 1:G
C(s,k,i,j) = cost(s,k);
end
end
end
end
%x(s,k,m,i,j) = 1 表示作业m在节点i和节点j之间由代理人s采用k种运输方式代理；否则x(s,k,m,i,j)=0；
x = zeros(g,G,M,n,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%p(s,l,i)表示在节点i由代理人s转换到代理人l的转换费用；
for i = 1:n
tmp = func_kcost2(i);
p(:,:,i) = tmp;
end
%r(s,l,m,i) = 1 表示作业m在节点i由代理人s转换成代理人l；否则r(s,l,m,i)=0;
r = zeros(g,g,M,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%R(k,v,m,i) = 1 表示作业m在节点i由k种运输方式转换为v种运输方式，否则R(k,v,m,i)=0;
R = zeros(G,G,M,n);
%Z(k,v,i)表示在节点i由k种运输方式转换为v种运输方式的单位中转费用；
for i = 1:n
tmp = func_kcost3(i);
Z(:,:,i) = tmp;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%T(m)表示作业m的时间期限；
%H(m)表示作业m的惩罚值；
%f(Tm,qm)为惩罚函数；
ET = [30.39 34 27 38.24 40.1 47.06 33.83 33.82];
LT = [31.50 40 28.5 39 43.5 48 34.5 35.00];
for m = 1:M
T(m) = LT(m) - ET(m);
end
%% 222222222222222222222222222222222222222222222222222222222222222
%% 222222222222222222222222222222222222222222222222222222222222222
%% 222222222222222222222222222222222222222222222222222222222222222
%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(2*num_particle,D,D,M);
V(:,:,:,:,1) = round((xmin-xmax)+2*(xmax-xmin)*rand(2*num_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(2*i-1:2*i,:,:,:,1)=particle(2*I-1:2*I,:,:,:,1); %群最优粒子位置
end
tmps = 0;
route = zeros(n,n,M,2*num_particle);
for ii=2:itmax
ii
V(:,:,:,:,ii) = 0.729*V(:,:,:,:,ii-1)+c1*rand*(pbest(:,:,:,:,ii-1)-particle(:,:,:,:,ii-1))+...
c2*rand*(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(2*i-1:2*i,:,:,:,ii) = particle(2*i-1:2*i,:,:,:,ii);
else
pbest(2*i-1:2*i,:,:,:,ii) = pbest(2*i-1:2*i,:,:,:,ii-1);
end
end
%*************************
%下面计算惯性权重
pmin = min(fit(:,ii));
a(ii) = mean(sum(abs(fit(:,ii)-pmin)));
%下面更新gbest and gbest_value
[Cs,I] = min(pbest_value(:,ii));
gbest_value(ii)=min(Cs,gbest_value(ii-1));
for i=1:num_particle
if gbest_value(ii) == Cs
gbest(2*i-1:2*i,:,:,:,ii)=pbest(2*I-1:2*I,:,:,:,ii);
else
gbest(2*i-1:2*i,:,:,:,ii)=gbest(2*I-1:2*I,:,:,:,ii-1);
end
end
end
%x2 = zeros(g,G,M,n,n,2*num_particle);
finals_cost = gbest_value(end);
% save Simulation_Results\result.mat gbest_value itmax
02_009m