基于模糊需求和模糊运输时间的多式联运路径优化(Matlab代码实现)

程序猿鑫

于 2024-04-23 09:15:34 发布

阅读量705

点赞数 10

文章标签： matlab 算法人工智能

本文链接：https://blog.csdn.net/Ke_Yan_She/article/details/138107880

版权

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码、数据

💥1 概述

文献来源：

摘要：考虑不确定性的模糊多式联运路径优化研究，可以在满足运输方案经济环保双重要求的同时，增强运输方案的鲁棒性，提高企业的抗风险能力。本文建立了模糊需求和模糊运输时间下低碳低成本多式联运路径优化模型，针对连续型元启发式算法无法直接求解离散型组合优化模型的问题，设计了基于优先级的通用编码方式；在此基础上，为进一步提高算法的求解质量，提出了带启发式因子的特殊解码方式；并且提出了一种带邻域搜索策略的自适应差分进化算法。结果表明，改进算法获得的最终方案在蒙特卡罗采样的大多数场景下满足约束，方案稳定性强，目标值最低。

关键词：

不确定优化;模糊;局部搜索;差分进化算法;蒙特卡罗采样;

随着确定性多式联运问题研究的成熟，研究者聚焦于更贴合实际的不确定多式联运问题。柳伍生等

人建立了基于模糊机会约束的绿色多式联运路径鲁棒优化模型，使用GS对实际算例进行求解，验证所建模型的鲁棒性[11]。李立等人考虑枢纽拥堵状况下路径选择的不确定性，建立绿色多式联运鲁棒优化模型。基于对偶理论将模型等价转化为混合整数线性规划模型，并使用Matlab/Yalmip 进行编程，调用Gurobi进行求解[12]。张旭等人针对需求与碳交易价格双重不确定的情景，建立混多式联运路径合鲁棒随机优化模型，在蒙特卡罗采样的基础上，使用灾变自适应遗传算法进行有效性检验[13]。Li 等人考虑碳排放多任务的多式联运路径优化模型与算法，设计分步方法结合遗传算法组合求解[14]。Peng等人研究了带时间表的不确定多式联运路线问题，设计了一种基于数据驱动的高效转换方法结合NSGAII算法求解算例[15]。

综上，诸多国内外学者对多式联运进行了研究，但仍存在部分局限性。首先，建立多式联运模型时

所考虑的时间窗一般被划分为两种，分别是硬性约束时间窗和软性惩罚时间窗；部分研究仅仅考虑了其中一种形式的时间窗，但实际运输中，不同节点对时间窗的要求是不同的。其次，大多数多式联运模型中的所有参数均为确定值，但在实际运输中，部分因素是不确定的，例如在货物的实际运输中，由于季节性需求、突然补货等因素，运输量会发生不确定的变动。对于考虑到多式联运中的不确定因素的研究，其中大部分论文的求解方法主要是GS或其他类似的求解器，但求解器无法求解大规模问题；另一些论文使用了元启发式算法，但其主要是离散类算法如蚁群、遗传算法等，连续型元启发式算法如粒子群和差分进化等几乎没有被使用，此外这些文章并没有详细展示元启发式算法和不确定优化问题耦合点。鉴于上述研究的局限性，本文主要从以下几个方面进行考虑：首先，将货物的运输量和各子路径不同运输方式的运输时间视为不确定量，用三角模糊数表示，研究模糊需求和模糊运输时间下的低碳低成本多式联运路径优化问题。之后在模型中引入混合时间窗，货物必须在硬时间窗内抵达，超出软时间窗范围，则须额外支付提前或者延迟费用。此外考虑到中国碳排放政策的特点和要求，本文将碳排放作为约束加入到模型。最后，本文基于蒙特卡罗采样结合改进差分进化算法对不确定多式联运优化模型进行求解，详细阐述了元启发式算法求解不确定优化模型的全部流程，即不确定优化模型转化，编码设计、解码设计，蒙特卡罗采样，以及适应度函数设计，为元启发式算法求解不确定优化模型提供一种标准的范式。

考虑不确定性的模糊多式联运路径优化研究，可以在满足运输方案经济环保双重要求的同时，增强运输方案的鲁棒性，提高企业的抗风险能力。本文建立了模糊需求和模糊运输时间下低碳低成本多式联运路径优化模型，针对连续型元启发式算法无法直接求解离散型组合优化模型的问题，设计了基于优先级的通用编码方式；在此基础上，为进一步提高算法的求解质量，提出了带启发式因子的特殊解码方式。

2024.4.4 之前的代码缺少不确定时间的处理和启发因子，已经添加上了

📚2 运行结果

主函数代码：

%使用AFO算法以及其他GA和PSO算法求解不确定多式联运路径优化问题。
%同时和MATLAB自带的全局优化搜索器进行对比
%Matlab版本需要2021以后的版本
% 考虑不确定性的模糊多式联运路径优化研究，可以在满足运输方案经济环保双重要求的同时，
%增强运输方案的鲁棒性，提高企业的抗风险能力。
%本文建立了模糊需求和模糊运输时间下低碳低成本多式联运路径优化模型，针对连续型元启发式算法无法直接求解离散型组合优化模型的问题，设计了基于优先级的通用编码方式；在此基础上，为进一步提高算法的求解质量，提出了带启发式因子的特殊解码方式。
% 2024.4.4 之前的代码缺少不确定时间的处理和启发因子，已经添加上了
clc;
clear all;
close all;
warning off
%%
noRng=1;
rng('default')
rng(noRng)
%%
global data
data.numN=15; %节点数量
data.Cap_Ts=xlsread("节点处的最大中转运输能力.xlsx");
temp=round(rand(data.numN,3)*10+15)*10;
data.Cap_Ts(:,2:end)=temp;
data.Windows=xlsread("节点的时间窗.xlsx");
data.D=xlsread("节点间距离.xlsx");
data.Cap_Tp=xlsread("节点间最大运输能力.xlsx");
data.T=data.D;
data.v=[76,60,30];
for i=1:length(data.Cap_Tp(:,1))
no1=data.Cap_Tp(i,1);
no2=data.Cap_Tp(i,2);
for j=1:3
if isnan(data.Cap_Tp(i,2+j))
data.D(i,2+j)=nan;
end
end
data.T(i,[3,6,9])=round(data.D(i,3:5)./data.v/1.2);
data.T(i,[4,7,10])=round(data.D(i,3:5)./data.v);
data.T(i,[5,8,11])=round(data.D(i,3:5)./data.v/0.8);
end
data.Windows(:,3)=data.Windows(:,3).*(1-rand(size(data.Windows(:,3))));
data.Windows(:,4)=data.Windows(:,4).*(1+rand(size(data.Windows(:,3))))+10;
data.Windows(:,3:4)=max(0,round(data.Windows(:,3:4)+randn(size(data.Windows(:,3:4)))));
data.Windows(:,5)=max(data.Windows(:,4));
data.CT=[0,3.09,5.23; % 转换成本
3.09,0,26.62;
5.23,26.62,0];
data.TT=[0,1,1; %转换时间
1,0,2;
1,2,0];
data.ET=[0,1.56,6; % 转换碳排放
1.56,0,3.12;
6,3.12,0];
data.q=[120,150,180];

data.E0=[0.796,0.028,0.04];
data.CW=[30,50];
data.S=1;
data.E=15;
data.alpha=0.8;
data.beta1=0.8;
data.beta2=0.8;
data.beta3=0.8;
data.C0=[0.3,0.2,0.1]; %三种运输方式的运输成本
data.weight=[1,1];
data.maxB=100000;
data.maxE=21000;
%%
%%
data.numQ=100;
for ii=1:data.numQ
if rand<0.5
data.q0(ii)=rand*(data.q(2)-data.q(1))+data.q(1);
else
data.q0(ii)=rand*(data.q(3)-data.q(2))+data.q(2);
end
for i=1:length(data.Cap_Tp(:,1))
for j=1:3
if rand<0.5
data.T0{ii}(i,j)=data.T(i,j*3)+rand*(data.T(i,j*3+1)-data.T(i,j*3));
else
data.T0{ii}(i,j)=data.T(i,j*3+1)+rand*(data.T(i,j*3+2)-data.T(i,j*3+1));
end
end
end
end
%%
G=graph(data.D(:,1),data.D(:,2),data.D(:,1)*0+1);
figure
plot(G)
set(gca,'LooseInset',get(gca,'TightInset'))
%%
lb=0;
ub=1;
dim=length(data.D(:,1))*3;
option.lb=lb;
option.ub=ub;
option.dim=dim;
if length(option.lb)==1
option.lb=ones(1,option.dim)*option.lb;
option.ub=ones(1,option.dim)*option.ub;
end
option.fobj=@aimFcn_1;
option.showIter=0;

%% 算法参数设置 Parameters
% 基本参数
option.numAgent=50; %种群个体数 size of population
option.maxIteration=100; %最大迭代次数 maximum number of interation
%% 遗传算法
option.p1_GA=0.7;
option.p2_GA=0.1;
%% 粒子群
option.w_pso=0.1;
option.c1_pso=1.2;
option.c2_pso=1.2;
%% AFO

option.v_lb=-(option.ub-option.lb)/4;
option.v_ub=(option.ub-option.lb)/4;
option.w2=0.5; %weight of Moving strategy III
option.w4=1;%weight of Moving strategy III
option.w5=1;%weight of Moving strategy III
option.pe=0.01; % rate to judge Premature convergence
option.gap0=ceil(sqrt(option.maxIteration*2))+1;
option.gapMin=5; % min gap
option.dec=2; % dec of gap
option.L=10; % Catastrophe
%% DE
option.F=0.5;
option.CR=0.5;
%%
str_legend=[{'GA'},{'PSO'},{'AFO'}];
aimFcn=[{@GA},{@PSO},{@AFO3}];
%% 初始化
rng(1)
x=ones(option.numAgent,option.dim);
y=ones(option.numAgent,1);
for i=1:option.numAgent
x(i,:)=rand(size(option.lb)).*(option.ub-option.lb)+option.lb;
y(i)=option.fobj(x(i,:),option,data);
end
%% 使用算法求解
bestX=x;
for i=1:length(aimFcn)
rng(noRng)
tic
[bestY(i,:),bestX(i,:),recording(i)]=aimFcn{i}(x,y,option,data);
tt(i)=toc;
end
%% 绘制迭代曲线
figure
hold on
for i=1:length(aimFcn)

if i>1
plot((recording(i).bestFit),'LineWidth',2)
else
plot((recording(i).bestFit),'--','LineWidth',2)
end
end
legend(str_legend)
xlabel('评价次数（*100）')
ylabel('适应度函数值')
set(gca,'LooseInset',get(gca,'TightInset'))
%% 计算结果
for i=1:length(str_legend)
str=[str_legend{i},'优化后方案'];
[~,result(i)]=option.fobj(bestX(i,:),option,data);

% drawPC(result(i),data,str)
end
%%
rng(1)
x0=[0.227391667368465 0.879816756512174 0.758978072405287 0.510838515791228 0.192188536212261 0.651441269146713 0.418963422194475 0.552435663022093 0.692178644766480 0.525121273337296 0.901443303774014 0.870433999021102 0.727907836575027 0.635555150728710 0.377948489141470 0.0665747380900594 0.380375322185706 0.271295139005737 0.556390072117235 0.836324594294146 0.412650509862244 0.435408044462144 0.339160221282954 0.517673083186858 0.305139543163461 0.759153928221489 0.661651047625365 0.695663261206050 0.579862243307374 0.0203768757943083 0.567587976562174 0.417560315118765 0.194865572678659 0.0355023446924871 0.0656090882042756 0.283822647397696 0.409304642722500 0.797371466667136 0.139734109717042 0.943567670211017 0.736340144724420 0.899302180486718 0.0958418793238567 0.494947220485218 0.377043083335839 0.392090901825893 0.872224588221799 0.225018973369553 0.658496649609525 0.985772242893824 0.312242936944674 0.0176752228038184 0.151305501546866 0.117920480497746 0.137344184189911 0.929080746036836 0.737056036783124 0.264583371593473 0.935709972806256 0.930542122124921 0.458371373498295 0.0216966409045397 0.0615124906604629 0.569072523079761 0.912715979600965 0.684016483416011 0.0872976004423956 0.854123989826840 0.699891878132857 0.603233918882651 0.626266466189881 0.764357628595205 0.633979476955922 0.162893341780099 0.277411062836412 0.839224723775726 0.908550538167801 0.890586782754062 0.466391769502905 0.488448750263230 0.365236675941588 0.730292329677324 0.748971571306197 0.706903738180738];
gs = GlobalSearch;
problem = createOptimProblem('fmincon','x0',x0,...
'objective',@(x)aimFcn_1(x,option,data),'lb',x0*0,'ub',x0*0+1);
x = run(gs,problem);
[fit1,result(i+1)]=aimFcn_1(x,option,data);