m基于PSO粒子群优化的第四方物流的作业整合算法matlab仿真,对比有代理人和无代理人两种模式下最低运输费用、代理人转换费用、运输方式转化费用和时间惩罚费用

目录

1.算法概述

2.仿真效果预览

3.核心MATLAB程序

4.完整MATLAB程序

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

4.完整MATLAB程序

matlab源码说明_我爱C编程的博客-CSDN博客

V