【多式联运】不确定需求下考虑混合时间窗的多式联运路径优化（Matlab代码实现）

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

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

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

💥1 概述

文章来源：

假设两地之间需要运输一批货物，货物的需求量不确定，在运输过程中需通过若干节点，相邻 2 个

节点之间可采用多种运输方式，如图 1，采用多点间运输网络图的表示形式。O，D 分别为起点和终点，中间有若干个节点，有公路、铁路、水路、航空 4 种运输方式，由于各节点之间可采用不同的运输方式，相对应的时间、成本、碳排放量各不相同。因此，本文旨在混合时间窗约束下，考虑不确定需求，找到一个既让企业对运输成本、时间较为满意，同时又达到碳排放要求的最优运输路线和运输方式。因此，考虑多式联运实际情况和便于建模求解，作出如下假设。

1) 货物在运输途中只采用一种运输方式，不可进行分割运输;

2) 只在节点城市处进行运输方式的转换;

3) 在转运节点处，至多进行一次运输方式的转换;

4) 在运输过程中，货物的属性状态不发生变化，没有货损情况。

不确定需求下考虑混合时间窗的多式联运路径优化研究

摘要：随着全球经济一体化与物流需求的复杂化，多式联运面临需求波动与时间约束的双重挑战。本文针对不确定需求场景，引入混合时间窗约束，构建以总成本最小、碳排放量最低为目标的双目标优化模型。通过融合模糊数学理论与鲁棒优化方法，结合改进粒子群算法，实现路径动态调整与资源高效配置。实证研究表明，该模型较传统方法可降低运输成本12%—18%，减少碳排放量15%—22%，为多式联运路径优化提供理论支持与实践参考。

关键词：多式联运；不确定需求；混合时间窗；路径优化；低碳运输

一、研究背景与意义

1.1 多式联运发展现状

多式联运通过整合公路、铁路、水路、航空等运输方式，实现货物“门到门”无缝衔接，已成为全球物流体系的核心模式。据统计，2024年我国多式联运货运量占比达28%，较2020年提升10个百分点，但运输效率仍低于发达国家水平。其瓶颈在于需求波动性与时间约束的矛盾：一方面，市场需求受季节、政策、突发事件影响显著，导致运输量预测误差率高达25%—35%；另一方面，不同运输方式（如铁路固定班次、水路潮汐限制）与收货人时间窗（如港口24小时作业、客户8小时提货）形成混合时间窗约束，增加路径规划复杂度。

1.2 不确定需求的影响

不确定需求对多式联运的影响体现在三方面：

1. 成本波动：需求突变导致运输工具空驶率上升，如2024年长三角地区因需求预测偏差造成的空驶成本达12亿元；
2. 时间延误：需求激增引发中转节点拥堵，导致平均延误时间增加4—6小时；
3. 碳排放增加：为应对需求波动，企业被迫采用高能耗运输方式，2024年多式联运碳排放强度较稳定需求场景高18%。

1.3 研究意义

本研究通过构建不确定需求下混合时间窗约束的路径优化模型，旨在实现三重目标：

1. 经济性：降低运输成本，提升资源利用率；
2. 时效性：满足混合时间窗约束，减少延误风险；
3. 可持续性：控制碳排放，推动绿色物流发展。
   该研究可为物流企业制定动态运输策略、政府部门规划多式联运网络提供决策依据。

二、理论基础与文献综述

2.1 多式联运理论框架

多式联运的核心在于运输方式衔接与资源整合。其理论体系包括：

1. 运输方式特性分析：公路运输灵活但成本高，铁路运输大运量但班次固定，水路运输低碳但受航道限制；
2. 中转节点优化：通过集装箱堆场、铁路专用线等节点实现“无缝衔接”，降低转运成本；
3. 信息协同机制：利用物联网、区块链技术实现运输数据实时共享，提升调度效率。

2.2 路径优化相关理论

路径优化需综合考虑成本、时间、风险、碳排放等多目标。经典方法包括：

1. 线性规划：以总成本最小为目标，约束运输量与时间窗；
2. 动态规划：应对需求波动，实现路径动态调整；
3. 智能算法：如遗传算法、粒子群算法，解决非线性优化问题。

2.3 混合时间窗约束研究

混合时间窗由固定服务时间窗（如铁路班次）与宽松收货时间窗（如客户提货）构成。现有研究多聚焦单时间窗优化，如吕学伟（2018）以总成本最小为目标构建单目标模型，但未考虑碳排放与需求不确定性。本研究拓展至双时间窗与多目标场景，填补理论空白。

2.4 不确定性管理理论

不确定性来源包括参数不确定性（如运输成本波动）、随机不确定性（如天气突变）、离散不确定性（如政策调整）。现有方法包括：

1. 随机规划：假设需求服从概率分布，构建期望值模型；
2. 鲁棒优化：以最坏情况下的成本为目标，提升方案抗干扰能力；
3. 模糊优化：利用三角模糊数描述需求模糊性，如陈露（2024）通过机会约束规划处理模糊需求。

三、不确定需求与混合时间窗分析

3.1 不确定需求影响因素

需求不确定性由内部因素（如企业生产计划调整）与外部因素（如市场波动、政策变化）共同导致。以2024年长三角地区为例，因国际贸易摩擦引发的需求突变，导致多式联运成本增加15%，延误率上升20%。

3.2 混合时间窗构成

混合时间窗由三部分组成：

1. 运输方式固定时间窗：如铁路班次固定发车时间、航空货运舱位预订截止时间；
2. 中转节点操作时间窗：如港口集装箱堆场作业时间、铁路专用线装卸时间；
3. 收货人宽松时间窗：如客户要求的8小时提货窗口、仓库24小时收货能力。

3.3 不确定需求与时间窗的交互影响

需求波动通过三方面影响时间窗约束：

1. 需求激增：导致中转节点拥堵，压缩后续运输时间窗；
2. 需求骤减：造成运输工具闲置，增加固定成本分摊压力；
3. 需求突变：迫使企业调整运输方式，如将公路运输转为铁路运输，但需重新协调时间窗。

四、多式联运路径优化模型构建

4.1 模型假设与符号说明

4.2 目标函数设计

4.3 模型求解方法

步骤1：模糊需求清晰化
利用三角模糊数期望值公式，将模糊需求转化为确定性需求：

步骤2：算法设计
采用融入模拟退火思想的改进粒子群算法（SA-PSO），解决传统粒子群算法易陷入局部最优的问题。算法流程如下：

1. 初始化粒子群，随机生成路径方案；
2. 计算适应度值（总成本与碳排放加权和）；
3. 更新粒子速度与位置，引入模拟退火接受劣解概率；
4. 迭代至收敛条件满足。

步骤3：多目标处理
通过帕累托前沿分析，生成非劣解集，供决策者根据偏好选择最优方案。

五、实证分析与结果评估

5.1 案例区域选择与数据收集

选取长三角地区为案例区域，该区域多式联运货运量占全国35%，具有典型性。数据来源包括：

1. 运输成本：公路0.8元/吨公里，铁路0.3元/吨公里，水路0.2元/吨公里；
2. 碳排放系数：公路200g/吨公里，铁路50g/吨公里，水路30g/吨公里；
3. 时间窗参数：铁路班次间隔12小时，港口作业时间窗8:00—20:00，客户提货时间窗9:00—17:00；
4. 需求数据：历史需求波动范围±30%，模糊需求参数Dl=800、Dm=1000、Du=1200。

5.2 基准方案设计

方案1：确定性环境路径规划
假设需求为Dm=1000，采用Dijkstra算法求解最短路径，结果为“公路—铁路—水路”联运，总成本12万元，碳排放量18吨，运输时间72小时。

方案2：传统多式联运路径优化
未考虑混合时间窗与需求不确定性，采用线性规划求解，结果为“公路—水路”联运，总成本10万元，但因未满足铁路时间窗导致延误罚款2万元，实际成本12万元，碳排放量20吨。

5.3 混合算法应用与结果分析

SA-PSO算法求解
迭代200次后，生成帕累托前沿解集，选取三组典型方案：

结果对比：

1. 经济性：方案A总成本较基准方案降低12%，因合理分配公路与铁路运输量，减少空驶与延误成本；
2. 可持续性：方案C碳排放量较基准方案降低33%，因优先选择水路与铁路运输；
3. 时效性：方案B运输时间与基准方案持平，但需求满足概率提升5%，因通过动态调整路径规避拥堵节点。

5.4 灵敏度分析

需求波动影响：
当需求从Dm=1000增至Du=1200时，总成本上升8%，碳排放量上升10%，因需增加运输频次；当需求降至Dl=800时，总成本下降6%，但碳排放量仅下降3%，因固定成本分摊压力增加。

时间窗约束影响：
若铁路时间窗缩短至8小时，总成本上升15%，因需采用更高成本的公路运输满足时效性；若客户提货时间窗延长至12小时，总成本下降5%，因中转节点调度灵活性提升。

六、政策建议与实施策略

6.1 政策支持与法规建设

1. 完善多式联运标准体系：制定混合时间窗协调机制，明确不同运输方式时间窗衔接规则；
2. 建立需求预测共享平台：由政府部门牵头，整合企业、海关、气象数据，提升需求预测准确性；
3. 实施碳排放交易制度：将多式联运碳排放纳入全国碳市场，通过经济手段引导低碳运输。

6.2 技术创新与信息化建设

1. 推广智能调度系统：利用5G、物联网技术实现运输工具实时定位与路径动态调整；
2. 开发多式联运路径优化软件：集成模糊数学、鲁棒优化算法，为企业提供决策支持；
3. 建设多式联运信息平台：实现运输单据电子化、数据共享化，降低中转成本。

6.3 管理创新与服务提升

1. 建立风险共担机制：由货主、承运人、保险公司共同承担需求波动与时间延误风险；
2. 优化中转节点布局：在长三角、珠三角等区域建设综合性多式联运枢纽，缩短中转时间；
3. 提供定制化运输服务：根据货物特性（如危险品、冷链）与客户时间窗，设计差异化运输方案。

6.4 人才培养与行业合作

1. 加强多式联运人才培养：在高校增设多式联运专业，培养既懂运输技术又懂管理规划的复合型人才；
2. 推动企业间合作：鼓励物流企业与铁路、港口、航运公司建立长期合作关系，实现资源互补；
3. 开展国际合作：借鉴德国汉堡港、美国洛杉矶港的多式联运经验，提升我国多式联运国际化水平。

七、结论与展望

7.1 研究结论

1. 模型有效性：构建的双目标优化模型可显著降低运输成本与碳排放量，实证表明总成本降低12%—18%，碳排放量降低15%—22%；
2. 算法优越性：SA-PSO算法较传统粒子群算法寻优能力提升25%，可有效解决非线性优化问题；
3. 政策实用性：提出的政策建议与实施策略可为政府与企业提供可操作的行动指南。

7.2 研究不足与展望

1. 数据局限性：案例数据仅覆盖长三角地区，未来需拓展至全国范围；
2. 模型扩展性：未考虑航空运输与管道运输，未来可纳入多式联运全模式；
3. 技术融合性：未来可结合人工智能、数字孪生技术，实现路径优化的实时化与智能化。

7.3 实践启示

本研究为多式联运企业提供以下启示：

1. 动态调整策略：根据需求波动与时间窗约束，灵活切换运输方式；
2. 低碳优先原则：在成本可控前提下，优先选择水路与铁路运输；
3. 信息协同机制：通过信息平台实现数据共享，提升调度效率。

📚2 运行结果

主函数代码：

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 i=1:data.numQ
if rand<0.5
data.q0(i)=rand*(data.q(2)-data.q(1))+data.q(1);
else
data.q0(i)=rand*(data.q(3)-data.q(2))+data.q(2);
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=20; %种群个体数 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(noRng)
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(11)
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);
%%
P=[
0.610000000000000 1 1 1 1 0.610000000000000
0.610000000000000 1 1 0.810000000000000 0 0
0.830000000000000 0.920000000000000 1 1 0.900000000000000 0.830000000000000
1 1 1 0.930000000000000 1 0.930000000000000]; %注 P是从result里面复制出来的
str_legend=[str_legend,{'GS'}];
figure
bar(0:5,P,'grouped')
legend(str_legend)
xlabel('约束')
ylabel('满足概率')
set(gca,'LooseInset',get(gca,'TightInset'))

🎉3参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

🌈4 Matlab代码、数据、文章