1.算法描述
PSO算法是一种随机的、并行的优化算法。它的优点是：不要求被优化函数具有可微、可导、连续等性质，收敛速度较快，算法简单，容易编程实现。然而，PSO算法的缺点在于：（1）对于有多个局部极值点的函数，容易陷入到局部极值点中，得不到正确的结果。造成这种现象的原因有两种，其一是由于待优化函数的性质；其二是由于微粒群算法中微粒的多样性迅速消失，造成早熟收敛。这两个因素通常密不可分地纠缠在一起。（2）由于缺乏精密搜索方法的配合，PSO算法往往不能得到精确的结果。造成这种问题的原因是PSO算法并没有很充分地利用计算过程中获得的信息，在每一步迭代中，仅仅利用了群体最优和个体最优的信息。（3）PSO算法虽然提供了全局搜索的可能，但是并不能保证收敛到全局最优点上。（4）PSO算法是一种启发式的仿生优化算法，当前还没有严格的理论基础，仅仅是通过对某种群体搜索现象的简化模拟而设计的，但并没有从原理上说明这种算法为什么有效，以及它适用的范围。因此，PSO算法一般适用于一类高维的、存在多个局部极值点而并不需要得到很高精度解的优化问题。
当前针对PSO算法开展的研究工作种类繁多，经归纳整理分为如下八个大类：（1）对PSO算法进行理论分析，试图理解其工作机理；（2）改变PSO算法的结构，试图获得性能更好的算法；（3）研究各种参数配置对PSO算法的影响；（4）研究各种拓扑结构对PSO算法的影响；（5）研究离散版本的PSO算法；（6）研究PSO算法的并行算法；（7）利用PSO算法对多种情况下的优化问题进行求解；（8）将PSO算法应用到各个不同的工程领域。以下从这八大类别着手，对PSO算法的研究现状作一梳理。由于文献太多，无法面面俱到，仅捡有代表性的加以综述。
PSO初始化为一群随机粒子（随机解）。然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个“极值（pbest和gbest）”来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。
根据要求文档，优化目标函数为：
对于公式（1）：
公式（1）中的第一部分称为记忆项，表示上次速度大小和方向的影响；
公式（1）中的第二部分称为自身认知项，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；
公式（1）中的第三部分称为群体认知项，是一个从当前点指向种群最好点的矢量，反映了粒子间的协调合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。
2.matlab算法仿真效果
matlab2022a仿真结果如下：
3.MATLAB核心程序 `%Pia物品a的惩罚单价
Pia = [1200,450,9000,4700];
%STOCa:物品a的存储量
STOCa = [43010.16,70850.84,14170.33,21255.18];
%Wa: 物品a的单位重量
Wa = [4,18,40,5];
%Ca: 物品a的单位体积
Ca = [0.054,0.054,0.2,0.01887];
%ra: 物品a的重要程度
ra = [0.35,0.3,0.2,0.15];
%vwm: 交通工具m的载重重量
VWm= [4100,14000];
%vcm: 交通工具m的载重体积
VCm= [45,1240];
%NVm: 初始可用交通工具数量
NVm= [6,3];
%TCm:交通工具m的燃油费
TCm= [2.61,0.00295];
%QCm:后期租车费用
QCm= [28800,900];
%Probs: 每种情景发生的概率
Probs=[0.0270,0.063,0.063,0.147,0.063,0.147,0.147,0.343];
Mm = 100000;
%LJ: 受灾点 j 到仓库的距离
Lj=[75.4,113,204.9,121,93.2,38.1];
%在情景s的第t阶段下用交通工具m到受灾点j的路线可通
% [1,1,1,1,1,1,1;
% 1,1,1,1,1,1,1;
% 1,1,1,1,1,1,1;
% 1,1,1,1,1,1,1;
% 1,1,1,1,1,1,1;
% 1,1,1,1,1,1,1;
% 1,1,1,1,1,1,1;
% 1,1,1,1,1,1,1;
% 1,0,1,1,0,0,1;
% 1,0,1,1,0,0,1;
% 1,0,1,1,0,0,1;
% 1,0,1,1,0,0,1;
% 1,0,1,1,0,0,1;
% 1,0,1,1,0,0,1;
% 1,0,1,1,0,0,1;
% 1,0,1,1,0,0,1;];
Rsjmt = zeros(M,S,T,J);
Rsjmt(1,1,1,1:6) = 1;
Rsjmt(1,1,2,1:6) = 1;
Rsjmt(1,2,1,1:6) = 1;
Rsjmt(1,2,2,1:6) = 1;
Rsjmt(1,3,1,1:6) = 1;
Rsjmt(1,3,2,1:6) = 1;
Rsjmt(1,4,1,1:6) = 1;
Rsjmt(1,4,2,1:6) = 1;
Rsjmt(1,5,1,1:6) = 1;
Rsjmt(1,5,2,1:6) = 1;
Rsjmt(1,6,1,1:6) = 1;
Rsjmt(1,6,2,1:6) = 1;
Rsjmt(1,7,1,1:6) = 1;
Rsjmt(1,7,2,1:6) = 1;
Rsjmt(1,8,1,1:6) = 1;
Rsjmt(1,8,2,1:6) = 1;
Rsjmt(2,1,1,1:6) = [1,0,1,1,0,0];
Rsjmt(2,1,2,1:6) = [1,0,1,1,0,0];
Rsjmt(2,2,1,1:6) = [1,0,1,1,0,0];
Rsjmt(2,2,2,1:6) = [1,0,1,1,0,1];
Rsjmt(2,3,1,1:6) = [1,0,1,1,0,0];
Rsjmt(2,3,2,1:6) = [1,0,1,1,1,0];
Rsjmt(2,4,1,1:6) = [1,0,1,1,0,0];
Rsjmt(2,4,2,1:6) = [1,0,1,1,1,1];
Rsjmt(2,5,1,1:6) = [1,0,1,1,0,0];
Rsjmt(2,5,2,1:6) = [1,1,1,1,0,0];
Rsjmt(2,6,1,1:6) = [1,0,1,1,0,0];
Rsjmt(2,6,2,1:6) = [1,1,1,1,0,1];
Rsjmt(2,7,1,1:6) = [1,0,1,1,0,0];
Rsjmt(2,7,2,1:6) = [1,1,1,1,1,0];
Rsjmt(2,8,1,1:6) = [1,0,1,1,0,0];
Rsjmt(2,8,2,1:6) = [1,1,1,1,1,1];
.......................................
%变量1的限制
for s1 = 1:S
for t1 = 1:T
for m1 = 1:M
for j1 = 1:J
for a1 = 1:A
if Xsajmt(s1,t1,m1,j1,a1,i) >= max1
Xsajmt(s1,t1,m1,j1,a1,i) = max1;
end
if Xsajmt(s1,t1,m1,j1,a1,i) <= min1
Xsajmt(s1,t1,m1,j1,a1,i) = min1;
end
end
end
end
end
end
%UDsajt
%UDsajt
%速度2设置
for s1 = 1:S
for t1 = 1:T
for j1 = 1:J
for a1 = 1:A
vb(s1,t1,j1,a1,i) = w*vb(s1,t1,j1,a1,i)+...
c1rand(1)(UDsajt_best(s1,t1,j1,a1,i)-UDsajt(s1,t1,j1,a1,i))+...
c2rand(1)(TUDsajt_best(s1,t1,j1,a1)-UDsajt(s1,t1,j1,a1,i));
end
end
end
end
%更新
for s1 = 1:S
for t1 = 1:T
for j1 = 1:J
for a1 = 1:A
UDsajt(s1,t1,j1,a1,i) = UDsajt(s1,t1,j1,a1,i) + vb(s1,t1,j1,a1,i);
end
end
end
end
%变量2的限制
for s1 = 1:S
for t1 = 1:T
for j1 = 1:J
for a1 = 1:A
if UDsajt(s1,t1,j1,a1,i) >= max2
UDsajt(s1,t1,j1,a1,i) = max2;
end
if UDsajt(s1,t1,j1,a1,i) <= min2
UDsajt(s1,t1,j1,a1,i) = min2;
end
end
end
end
end
%QVsm
%QVsm
%速度3设置
for s1 = 1:S
for m1 = 1:M
vc(s1,m1,i) = w*vc(s1,m1,i)+...
c1rand(1)(QVsm_best(s1,m1,i)-QVsm(s1,m1,i))+...
c2rand(1)(TQVsm_best(s1,m1)-QVsm(s1,m1,i));
end
end
%更新
for s1 = 1:S
for m1 = 1:M
QVsm(s1,m1,i) = QVsm(s1,m1,i) + vc(s1,m1,i);
end
end
%变量3的限制
for s1 = 1:S
for m1 = 1:M
if QVsm(s1,m1,i) >= max3
QVsm(s1,m1,i) = max3;
end
if QVsm(s1,m1,i) <= min3
QVsm(s1,m1,i) = min3;
end
end
end
%Nsjmt
%Nsjmt
%速度4设置
for s1 = 1:S
for t1 = 1:T
for m1 = 1:M
for j1 = 1:J
vd(s1,t1,m1,j1,i) = w*vd(s1,t1,m1,j1,i)+...
c1rand(1)(Nsjmt_best(s1,t1,m1,j1,i)-Nsjmt(s1,t1,m1,j1,i))+...
c2rand(1)(TNsjmt_best(s1,t1,m1,j1)-Nsjmt(s1,t1,m1,j1,i));
end
end
end
end
%更新
for s1 = 1:S
for t1 = 1:T
for m1 = 1:M
for j1 = 1:J
Nsjmt(s1,t1,m1,j1,i) = Nsjmt(s1,t1,m1,j1,i) + vd(s1,t1,m1,j1,i);
end
end
end
end
%变量4的限制
for s1 = 1:S
for t1 = 1:T
for m1 = 1:M
for j1 = 1:J
if Nsjmt(s1,t1,m1,j1,i) >= max4
Nsjmt(s1,t1,m1,j1,i) = max4;
end
if Nsjmt(s1,t1,m1,j1,i) <= min4
Nsjmt(s1,t1,m1,j1,i) = min4;
end
end
end
end
end
[BsJ,Xsajmt(:,:,:,:,:,i),UDsajt(:,:,:,:,i),QVsm(:,:,i),Nsjmt(:,:,:,:,i)] = func_fitness(Xsajmt(:,:,:,:,:,i),UDsajt(:,:,:,:,i),round(QVsm(:,:,i)),round(Nsjmt(:,:,:,:,i)));
if BsJ<BsJi(i)
BsJi(i) = BsJ;
Xsajmt_best(:,:,:,:,:,i) = Xsajmt(:,:,:,:,:,i);
UDsajt_best(:,:,:,:,i) = UDsajt(:,:,:,:,i);
QVsm_best(:,:,i) = QVsm(:,:,i) ;
Nsjmt_best(:,:,:,:,i) = Nsjmt(:,:,:,:,i);
end
if BsJi(i)<minJi
minJi = BsJi(i);
TXsajmt_best(:,:,:,:,:,i) = Xsajmt(:,:,:,:,:,i);
TUDsajt_best(:,:,:,:,i) = UDsajt(:,:,:,:,i);
TQVsm_best(:,:,i) = QVsm(:,:,i) ;
TNsjmt_best(:,:,:,:,i) = Nsjmt(:,:,:,:,i);
end
end
Jibest(t) = minJi;
end
Xsajmt0 = Xsajmt(:,:,:,:,:,end);
UDsajt0 = UDsajt(:,:,:,:,end);
QVsm0 = round(QVsm(:,:,end));
Nsjmt0 = round(Nsjmt(:,:,:,:,end));
02_039m`
