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

最新推荐文章于 2025-12-19 16:12:26 发布

原创最新推荐文章于 2025-12-19 16:12:26 发布 · 233 阅读

10 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #算法 #人工智能 #支持向量机

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

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

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

💥1 概述

文献来源：

基于模糊需求和模糊运输时间的多式联运路径优化研究

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

关键词：

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

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

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

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

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

聚焦于存在不确定性因素的模糊多式联运路径优化研究，具有显著意义。此类研究不仅能够在满足运输方案经济性与环保性双重标准的前提下，有效增强运输方案的鲁棒性，还能显著提升企业应对风险的能力。

本文致力于构建一个在模糊需求与模糊运输时间情境下的低碳低成本多式联运路径优化模型。鉴于连续型元启发式算法难以直接应用于离散型组合优化模型的求解难题，本研究创新性地设计了一种基于优先级的通用编码方式。在此基础上，为进一步优化算法的求解效果，又提出了一种融合启发式因子的特殊解码方式。

对不确定时间的处理机制以及启发因子的应用，目前这些问题均已得到妥善解决并完成添加。

一、多式联运路径优化的基本概念与挑战

多式联运路径优化旨在通过整合铁路、公路、水路、航空等多种运输方式，实现运输成本、时间、碳排放等目标的最优组合。其核心问题包括运输节点选择、运输方式组合、路径序列规划以及不确定性处理。研究显示，现有模型常面临以下挑战：

不确定性建模：需求波动（如季节性变化、突发事件）和运输时间模糊性（如天气、交通拥堵）需通过模糊数学工具表达。
多目标冲突：成本最小化、时间最短化、碳排放最低化等目标需权衡。
复杂约束：货物换装成本、时间窗限制、节点容量限制等需综合考虑。

二、模糊需求建模方法

三角模糊数
采用 (dl,dm,du)表示模糊需求，其中 dmdm 为最可能值，dl和 du为上下界。例如，冷链物流中客户需求可表示为 D~=(d~1,d~2,d~3)，结合概率分布描述不同市场情景。
梯形模糊数
适用于更广泛的模糊区间，如运输订单需求表示为 [ef,ed]，通过隶属度函数量化预测误差。
模糊随机变量
结合概率与模糊性，如需求上升、稳定、下降的情景分别赋予概率和模糊数。
模糊机会约束规划
将模糊约束转化为置信水平下的确定性约束，例如货物运输量需满足 Cr(∑xij≥d~)≥α。

三、模糊运输时间建模方法

三角模糊数
运输时间 tij=(tl,tm,tu)tij=(tl,tm,tu)，其中 tmtm 为最可能时间，隶属度函数呈三角形分布。例如鲜果配送中，运输时间受天气影响，可能值区间为 (1,3,5)(1,3,5) 小时。
模糊预约时间窗
替代传统硬时间窗，通过凸模糊数表示客户时间偏好，如 μt(ti)μt(ti) 表示满意度函数。
动态模糊模拟
结合蒙特卡洛方法生成模糊时间样本，评估不同路径方案的鲁棒性。

四、考虑不确定性的扩展模型

鲁棒优化模型
引入盒状不确定性集处理运输时间和货运率的不确定性，通过双化方法转化为可求解的混合整数规划。
随机期望值模型
假设需求服从正态分布 Q∼N(μ,σ2)，将模糊目标转化为期望值优化。
模糊机会约束规划模型
结合置信水平 αα，将模糊约束转化为确定性约束，适用于应急物资运输。
绿色路径优化模型
在碳排放交易政策下，以低碳为目标，结合模糊需求与时间窗约束。

五、多目标优化算法设计

模糊期望值法
将模糊目标函数转换为确定性目标，例如总成本通过加权法整合多目标。
改进进化算法
- NSGA-III：基于参考点的非支配排序，适用于多目标优化。
- 自适应差分进化算法：引入邻域搜索策略，提升离散组合问题的求解效率。
混合智能算法
结合模糊模拟与遗传算法，生成初始解后通过局部搜索优化。
模糊机会约束规划与博弈论结合
在多目标模型中动态调整权重，平衡经济性与时间可靠性。

六、实际案例应用

集装箱多式联运案例
- 模糊需求建模：使用三角模糊数描述运量不确定性，结合Lingo软件求解路径方案。
- 运输时间模糊性：考虑节点作业时间窗，通过灵敏度分析验证客户满意度与成本的权衡。
冷链物流案例
- 模糊随机需求：定义 D~=(d~1,d~2,d~3)，结合概率分布优化库存与路径。
公路-铁路联运案例
- 模糊软时间窗：最大化服务水平，通过帕累托前沿实现经济性与及时性权衡。

七、未来研究方向

精细模糊建模：探索区间二型模糊数，增强对复杂不确定性的刻画。
实时动态优化：结合物联网数据，实现运输路径的实时调整。
多政策协同：整合碳税、碳交易等政策，提升绿色运输效果。
跨学科技术融合：结合机器学习预测模糊参数，优化算法效率。

八、结论

模糊需求与模糊运输时间的多式联运路径优化研究，通过三角模糊数、机会约束规划、进化算法等工具，有效应对现实中的不确定性。实际应用中需结合具体场景选择建模方法，并通过混合算法平衡多目标冲突。未来需进一步探索动态优化与智能化技术的融合，以实现更高效的物流决策支持。

📚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);