【多式联运】基于粒子群结合遗传算法实现陆海空多式联运问题附matlab代码

1 简介

物流运输方式由公路、水路、空运及管道等 3 种方式组成，3 种运输方式在技术上、经济上各有长短，都有适宜的使用范围，每种运输方式单独运用很难实现节约资源、降本增效。随着我国经济不断发展以及布局网络技术的不断深化，多式联运通过把传统的、单一的运输方式进行择优组合，充分利用了各个运输方式现有的设施设备，实现了运输过程中的资源整合，有利于运输过程中的可持续发展及达成规模经济中降本增效的目的，同时提高了物流行业竞争力。特别是通过公铁水多式联运路径优化，构建以运输时间最少、运输线路距离最短、运输成本最低的公铁水多式联运模式，对于物流企业节约资源、降本增效意义重大。本文提出了一个以粒子群结合遗传算法为主框架的解决方案,用来求解多式联运的路径规划问题。本文从运输需求内容、运输过程、应用场景等角度对多式联运在军事运输中的应用进行分析,定义多式联运路径规划问题,建立分别以时问最短、路线最短、成本最低为目标的多式联运路径规划模型。

2 部分代码

clear clc close all tic %% 用importdata这个函数来读取文件 shuju= xlsread('shuju.xlsx', 'Sheet1'); bl=0; cap=100; %车辆最大装载量 %% 提取数据信息 zuobiao=shuju(:,2:3); %所有点的坐标x和y customer=zuobiao(2:end,:); %顾客坐标 cusnum=size(customer,1); %顾客数 v_num=8; %车辆最多使用数目 demands=shuju(2:end,4); %需求量 a=shuju(2:end,5); %顾客时间窗开始时间[a[i],b[i]] b=shuju(2:end,6); %顾客时间窗结束时间[a[i],b[i]] s=shuju(2:end,7); %客户点的服务时间 %% 距离矩阵 h=pdist(zuobiao); lldist=squareform(h).*1.5; %路路距离矩阵 htdist=squareform(h); %飞机距离矩阵 hydist=squareform(h).*1.2; %%火车距离矩阵 %% 遗传算法参数设置 alpha=100000; %违反的容量约束的惩罚函数系数 belta=0.5;%违反时间窗约束的惩罚函数系数 belta2=0.5; chesu=[1,5,2]; NIND=300; %种群大小 MAXGEN=500; %迭代次数 Pc=0.9; %交叉概率 Pm=0.05; %变异概率 GGAP=0.9; %代沟(Generation gap) N=cusnum+v_num-1; %染色体长度=顾客数目+车辆最多使用数目-1 %% 粒子群参数 lx=3; w=1; %惯性因子 wdamp=0.99; %惯性因子衰减率 c1=1.5; %个体学习因子 c2=2.0; %全局学习因子 XvMin=1; %Xv下限 XvMax=lx; %Xv上限 VvMin=-(lx-1); %Vv下限 VvMax=lx-1; %Vv上限 %% 初始化种群 Chrom=InitPopCW(NIND,N,cusnum,a,demands,cap,XvMin,XvMax,VvMin,VvMax); %构造初始解 %% 输出随机解的路线和总距离 disp('初始种群中的一个随机值:') [VC,NV,PD]=decode(Chrom{1},cusnum); % disp(['总距离：',num2str(TD)]); disp(['车辆使用数目：',num2str(NV)]); disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~') %% 优化 gen=1; figure; hold on;box on xlim([0,MAXGEN]) title('优化过程') xlabel('代数') ylabel('最优值') ObjV=calObj(Chrom,cusnum,cap,demands,a,b,s,lldist,htdist,hydist,alpha,belta,belta2,chesu,bl); %计算种群目标函数值 preObjV=min(ObjV); [~,minInd]=min(ObjV); gbest_pos=Chrom{minInd,1}(2,:); %假设第一个粒子位置为全局最优位置 gbest_obj=preObjV; %第一个粒子位置的目标函数值 pbest_pos=cell(NIND,1); %初始化各个粒子的个体最优位置 pbest_obj=ObjV; %初始化各个粒子的个体最优的目标函数值 for i=1:NIND particle=Chrom{i,1}; %第i个粒子 position=particle(2,:); %第i个粒子的位置 pbest_pos{i,1}=position; %初始化这个粒子的个体最优 if pbest_obj(i,1)<gbest_obj %更新初始种群中的全局最优粒子 gbest_obj=pbest_obj(i,1); gbest_pos=position; end end %% while gen<=MAXGEN %% 计算适应度 ObjV=calObj(Chrom,cusnum,cap,demands,a,b,s,lldist,htdist,hydist,alpha,belta,belta2,chesu,bl); %计算种群目标函数值 line([gen-1,gen],[preObjV,min(ObjV)]);pause(0.0001)%画图最优函数 preObjV=min(ObjV); FitnV=Fitness(ObjV); %% 选择 SelCh=Select(Chrom,FitnV,GGAP); %% OX交叉操作 SelCh=Recombin(SelCh,Pc); %% 变异 SelCh=Mutate(SelCh,Pm); %% 重插入子代的新种群 Chrom=Reins(Chrom,SelCh,ObjV); ObjV=calObj(Chrom,cusnum,cap,demands,a,b,s,lldist,htdist,hydist,alpha,belta,belta2,chesu,bl); %% 更新各个粒子的位置和速度 Chrom1=Chrom; for i=1:NIND particle=Chrom{i,1}; %第i个粒子 position=particle(2,:); %第i个粒子的位置 Xv=position(1,:); velocity=particle(3,:); %第i个粒子的速度 Vv=velocity(1,:); %% 更新速度 velocity=w*velocity+ +c1*rand([1,N]).*(pbest_pos{i,1}-position)... +c2*rand([1,N]).*(gbest_pos-position); %% 速度越界处理 velocity(1,:)=max(velocity(1,:),VvMin); velocity(1,:)=min(velocity(1,:),VvMax); %% 更新位置 position=position+velocity; position(1,:)=ceil(position(1,:)); %对Xv向上取整 %% 速度镜像影响 IsOutside=(position(1,:)<XvMin | position(1,:)>XvMax ); velocity(IsOutside)=-velocity(IsOutside); %% 位置越界处理 position(1,:)=max(position(1,:),XvMin); position(1,:)=min(position(1,:),XvMax); Chrom1{i,1}(2,:)=position; end ObjV1=calObj(Chrom1,cusnum,cap,demands,a,b,s,lldist,htdist,hydist,alpha,belta,belta2,chesu,bl); %% 重插入子代的新种群 Chrom=Copy_of_Reins(Chrom,Chrom1,ObjV1,NIND,ObjV); %% 打印当前最优解 ObjV=calObj(Chrom,cusnum,cap,demands,a,b,s,lldist,htdist,hydist,alpha,belta,belta2,chesu,bl); %计算种群目标函数值 [~,minInd]=min(ObjV); for i=1:NIND if min(ObjV)<ObjV(i) particle=Chrom{minInd,1}; position=particle(2,:); pbest_pos{i,1}=position; pbest_obj(i,1)=ObjV(i); %% 更新全局最优 if preObjV<gbest_obj gbest_pos=pbest_pos{i,1}; gbest_obj=pbest_obj(i,1); end end end disp(['第',num2str(gen),'代最优解:']) fprintf('\n') %% 更新迭代次数 gen=gen+1 ; end %% 画出最优解的路线图 ObjV=calObj(Chrom,cusnum,cap,demands,a,b,s,lldist,htdist,hydist,alpha,belta,belta2,chesu,bl); %计算种群目标函数值 [minObjV,minInd]=min(ObjV); %% 输出最优解的路线和总距离 disp('最优解:') bestChrom=Chrom{minInd(1)}; [VC,NV,PD]=decode(bestChrom,cusnum); disp('-------------------------------------------------------------') [cost]=costFuction(VC,NV,PD,a,b,s,lldist,htdist,hydist,demands,cap,alpha,belta,belta2,chesu,bl); disp(['总成本：',num2str(cost)]); disp(['路路：',num2str(1)]); disp(['航天：',num2str(2)]); disp(['海运：',num2str(3)]); %% 画出最终路线图 draw_Best(VC,zuobiao,PD);

% toc

function [A,c]=cdzwz(dist,route,dis) dis=0.9*dis; n=length(route); p_l=0; DL=1; c=[]; A=[]; for i=1:n if i==1 p_l=p_l+dist(1,route(i)+1); DL=DL-(dist(1,route(i)+1)/dis); if DL<0 c=[c,i]; DL=1; end else p_l=p_l+dist(route(i-1)+1,route(i)+1); DL=DL-dist(route(i-1)+1,route(i)+1)/dis; if DL<0 [a,b]=min(dist(route(i-1)+1,end-9:end)); DL=DL+dist(route(i-1)+1,route(i)+1)/dis-a/dis; if DL<0 c=[c,i-1]; DL=1; [a,b]=min(dist(route(i-1)+1,end-9:end)); A=[A,b+40]; DL=DL-dist(route(i-1)+1,b+41)/dis; else c=[c,i]; DL=1; [a,b]=min(dist(route(i)+1,end-9:end)); A=[A,b+40]; DL=DL-dist(route(i)+1,b+41)/dis; end end    end end p_l=p_l+dist(route(end)+1,1); DL=DL-dist(route(end)+1,1)/dis; if DL<0 c=[c,n]; DL=1; [a,b]=min(dist(route(i)+1,end-9:end)); A=[A,b+40]; DL=DL-dist(route(i)+1,b+41)/dis; end end
function [A,c]=cdzwz(dist,route,dis)
dis=0.9*dis;
n=length(route);
p_l=0;
DL=1;
c=[];
A=[];
for i=1:n
if i==1
p_l=p_l+dist(1,route(i)+1);
DL=DL-(dist(1,route(i)+1)/dis);
if DL<0
c=[c,i];
DL=1;
end
else
p_l=p_l+dist(route(i-1)+1,route(i)+1);
DL=DL-dist(route(i-1)+1,route(i)+1)/dis;
if DL<0
[a,b]=min(dist(route(i-1)+1,end-9:end));
DL=DL+dist(route(i-1)+1,route(i)+1)/dis-a/dis;
if DL<0
c=[c,i-1];
DL=1;
[a,b]=min(dist(route(i-1)+1,end-9:end));
A=[A,b+40];
DL=DL-dist(route(i-1)+1,b+41)/dis;
else
c=[c,i];
DL=1;
[a,b]=min(dist(route(i)+1,end-9:end));
A=[A,b+40];
DL=DL-dist(route(i)+1,b+41)/dis;
end
end

end
end
p_l=p_l+dist(route(end)+1,1);
DL=DL-dist(route(end)+1,1)/dis;
if DL<0
c=[c,n];
DL=1;
[a,b]=min(dist(route(i)+1,end-9:end));
A=[A,b+40];
DL=DL-dist(route(i)+1,b+41)/dis;
end
end
function [A,c]=cdzwz(dist,route,dis)
dis=0.9*dis;
n=length(route);
p_l=0;
DL=1;
c=[];
A=[];
for i=1:n
if i==1
p_l=p_l+dist(1,route(i)+1);
DL=DL-(dist(1,route(i)+1)/dis);
if DL<0
c=[c,i];
DL=1;
end
else
p_l=p_l+dist(route(i-1)+1,route(i)+1);
DL=DL-dist(route(i-1)+1,route(i)+1)/dis;
if DL<0
[a,b]=min(dist(route(i-1)+1,end-9:end));
DL=DL+dist(route(i-1)+1,route(i)+1)/dis-a/dis;
if DL<0
c=[c,i-1];
DL=1;
[a,b]=min(dist(route(i-1)+1,end-9:end));
A=[A,b+40];
DL=DL-dist(route(i-1)+1,b+41)/dis;
else
c=[c,i];
DL=1;
[a,b]=min(dist(route(i)+1,end-9:end));
A=[A,b+40];
DL=DL-dist(route(i)+1,b+41)/dis;
end
end

end
end
p_l=p_l+dist(route(end)+1,1);
DL=DL-dist(route(end)+1,1)/dis;
if DL<0
c=[c,n];
DL=1;
[a,b]=min(dist(route(i)+1,end-9:end));
A=[A,b+40];
DL=DL-dist(route(i)+1,b+41)/dis;
end
end

%% 根据vehicles_customer整理出final_vehicles_customer，将vehicles_customer中空的数组移除
function [ final_VC,VU ] = deal_vehicles_customer( VC )
vecnum=size(VC,1); %车辆数
final_VC={}; %整理后的vehicles_customer
count=1; %计数器
for i=1:vecnum
par_seq=VC{i}; %每辆车所经过的顾客
%如果该辆车所经过顾客的数量不为0，则将其所经过的顾客数组添加到final_vehicles_customer中
if ~isempty(par_seq)
final_VC{count}=par_seq;
count=count+1;
end
end
%% 为了容易看，将上述生成的1行多列的final_vehicles_customer转置了，变成多行1列的了
final_VC=final_VC';
VU=size(final_VC,1); %所使用的车辆数
end

3 仿真结果

4 参考文献

[1]丁建伟. 基于遗传算法的多式联运应急管理研究[J]. 电脑开发与应用, 2009, 22(1):3.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。