0 引言

随着交通运输和经济社会发展的深度融合，多式联运以其独特的优势，逐渐成为促进交通运输高质量发展的核心策略，备受业界关注。在研究多式联运的相关问题时，充分考虑碳排放因素显得尤为重要。然而，多式联运过程面临着诸多挑战。一方面，货物需求的快速增长与当前多式联运路径优化能力的不足存在明显的矛盾；另一方面，综合交通网络布局结构不够合理，交通拥堵时有发生，路网阻抗对多式联运过程产生一定的影响。因此，研究路网阻抗背景下的考虑经济成本和碳排放因素的多式联运路径鲁棒优化问题，已变得尤为迫切。

在多式联运路径优化领域，国内外学者已就多式联运网络展开了广泛而深入的研究。吴晓黎^[1]构建了以运输成本和时间为双目标的优化模型，并运用改进的Dijkstra算法进行求解。刘丹^[2]提出综合考虑可持续运输时间、费用和碳排放的多目标多式联运优化模型。张艳萍^[3]建立系统总运营成本、总运输时间和运输车辆驾驶员心理成本的多目标混合整数动态规划模型。邓红星^[4]构建了一个以碳排放、运输费用和运输时间最小化为目标的多式联运多目标优化模型。吴志勇^[5]构建的模型不仅考虑运输费用，而且考虑碳排放和时间窗成本，设计鲸鱼优化算法求解模型。Resat^[6]设计了适合土耳其的多式联运模型。Yin^[7]构建了以总成本、运输时间和碳排放最小化为目标的0-1规划模型。Shoukat^[8]运用混合整数规划模型来解决成本最小化和二氧化碳排放问题。然而，在这些多目标多式联运路径优化的研究中，均未充分考虑到路网阻抗对运输时间和碳排放可能产生的实际影响。

此外，也有一些学者针对不确定条件下的多式联运路径优化问题进行了深入的研究。户佐安^[9]针对运输时间及需求不确定的多式联运路径优化问题，采用基于情景的鲁棒优化方法优化模型。杨喆^[10]在模糊需求和模糊时间的背景下，构建了旨在实现低碳与低成本的多式联运路径优化模型。刘松^[11]考虑需求和运输环境的双重不确定性，提出了应急物资多式联运可靠路径的优化模型。张旭^[12]针对需求与碳交易价格的不确定性，建立了混合鲁棒随机优化模型，以优化多式联运路径。陈泊梨^[13]运用随机规划理论对货物运输速度值、准时送达概率阈值及转运时间值进行了估算。Demir^[14]提出了一个综合考虑旅行时间不确定性的绿色多式联运服务网络设计模型。Wang^[15]针对需求不确定性，提出了一个公铁联运双目标模型，该模型同时考虑了期望值准则和临界值准则。Zhang^[16]借助轮辐网络的概念，构建了一个不确定多式联运优化模型。Kalinina^[17]提出了一个不确定交货期下的多式联运多目标机会约束规划模型，以应对交货期的不确定性。Rouky^[18]深入研究了不确定条件下鲁棒最优化理论在多式联运问题中的应用，探讨了其在应对不确定性方面的潜力。Tian^[19]针对不确定条件下的多式联运问题，提出了一种广义区间模糊混合整数规划模型，以应对运输过程中的不确定性因素。然而，当前关于考虑需求不确定的冷藏集装箱多式联运路径优化问题的相关研究仍相对较少，这一领域仍有待进一步探索。

综上所述，许多学者在低碳多式联运及不确定多式联运方面作了诸多研究，但同时考虑路网阻抗和需求不确定的冷藏集装箱多式联运路径问题研究仍然鲜有出现。本研究在路网阻抗背景下，构建考虑运输成本和碳排放成本的多式联运优化模型，采用改进的樽海鞘算法优化模型, 讨论鲁棒优化对不确定需求的决策方案的影响，以及不同遗憾系数下的成本与最优解的关联性，为企业制订合理的韧性多式联运方案提供理论支撑。同时，验证优化后的樽海鞘算法的有效性，从而进一步拓宽樽海鞘算法的应用范畴。

1 模型描述 1.1 问题描述

在客户可接受的时间窗内将一批货物从起点运输到终点，需要考虑如何充分利用不同运输方式的优势，包括公路运输、铁路运输和水路运输。同时，需要考虑交通拥堵、运输时间、运输需求的不确定性，以及集装箱的性质等综合因素。在此基础上，设计具有经济性、环境可持续性和鲁棒性的多式联运运输路线和运输方式，以实现最大程度降低总成本的目标，并提高货物运输效率。

1.2 研究假设

研究假设如下所示。

(1) 货物不可拆分，仅可整箱运输。

(2) 相邻节点间至多只有一种相同的运输方式。

(3) 货物中转只能发生在节点处。

(4) 节点处的货物中转最多只发生1次。

(5) 路网阻抗只发生在公路运输，排除铁路运输和水路运输交通拥堵的影响。

(6) 货物中转的等待时间忽略不计。

(7) 不考虑铁路运输和水路运输的时间班点。

(8) 运输需求波动范围已知。

1.3 符号说明

设定模型所需的集合、参数和决策变量，如表 1所示。

2 模型构建 2.1 成本模型

本研究综合考虑了运输成本、转运成本、时间成本、货损成本及冷藏成本等多个维度。

(1) 运输成本

节点i和节点j之间的运输成本为单位运输基价、运货量、运输距离的乘积。

$ C_{\mathrm{y}}=\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} C_{i j}^k \cdot q_{i j} \cdot d_{i j}^k \cdot x_{i j }^k。$

(1)

(2) 转运成本

节点i处的转运成本为单位转运成本和转运量的乘积。

$ C_z=\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} C_i^{k l} \cdot q_{i j} \cdot y_{i }^{k l}。$

(2)

(3) 时间成本

总时间T包括节点间的运输时间和节点处的转运时间。铁路和水路有专线运输，不存在拥堵情况，运输时间为距离与速度的比值。

$ T=\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} t_{i j}^k \cdot x_{i j}^k+\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} t_i^{k l} \cdot q_{i j} \cdot y_i^{k l}, $

(3)

$ t_{i j}^k=\frac{d_{i j}^k}{v_k} \cdot x_{i j }^k。$

(4)

因路网阻抗，公路运输会有交通拥堵，采用路阻函数来计算公路运输时间^[20]。

$ t_{i j}^k=t_{i j}^{k 0} \cdot\left[1+\alpha \cdot\left(\frac{f}{c}\right)^\beta\right] \cdot x_{i j}^k, $

(5)

式中t_ij^k0为没有交通拥堵的时间。

时间成本的计算公式为：

$ \begin{gathered} C_{\mathrm{t}}=P_{\mathrm{e}} \cdot q_{i j} \cdot \max \left\{T_{\min }-T, 0\right\}+ \\ P_1 \cdot q_{i j} \cdot \max \left\{T-T_{\max }, 0\right\}。\end{gathered} $

(6)

(4) 货损成本和冷藏成本

货损成本为单位货损成本与运输时间的乘积，冷藏成本为单位冷藏成本与运输时间的乘积。

$ C_{\mathrm{h}}=h_1 T+h_2 T $

(7)

2.2 碳排放模型

在冷藏集装箱多式联运路径优化问题中，碳排放主要源自运输过程、转运过程及冷藏过程3个环节。

(1) 运输碳排放量

运输过程的碳排放量为碳排放系数与燃油消耗的乘积。

$ E_y=\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} E_{i j}^k \cdot x_{i j}^k=\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} \theta^k \cdot e_{i j}^k \cdot x_{i j}^k。$

(8)

公路运输的燃油消耗采用综合模式排放模型计算^[21]。

$ F_{\mathrm{r}}=\varepsilon\left(\delta N V_{\mathrm{r}}+p / \eta\right) / \gamma。$

(9)

货车牵引功率和辅助功率计算公式为：

$ p=p_{\text {tract }} / \lambda+p_{\text {acc }}, $

(10)

$ \begin{gathered} p_{\text {tract }}=\left(M v_{\mathrm{a}}+M g \sin \theta+0.5 C_{\mathrm{d}} \rho S v^2+\right. \\ \left.M g C_{\mathrm{r}} \cos \theta\right) v / 1\;000。\end{gathered} $

(11)

$ \begin{gathered} e_{i j}^k=F_{\mathrm{r}} \cdot \frac{d_{i j}}{v}=\frac{\varepsilon d_{i j} \delta N V_{\mathrm{r}}}{\gamma v}+\frac{\varepsilon d_{i j} p_{\mathrm{acc}}}{\gamma v \eta}+ \\ \frac{\varepsilon d_{i j}\left(M v_{\mathrm{a}}+M g \sin \theta+0.5 C_{\mathrm{d}} \rho S v^2+M g C_{\mathrm{r}} \cos \theta\right)}{\gamma \eta \lambda 1\;000} 。\end{gathered} $

(12)

速度变化对尾气排放具有显著影响，特别是加速和减速过程对尾气排放的影响更为显著，这一影响超过了单纯考虑恒定速度时的情况。然而，在实际路况下，车辆配送的动态性使得精确模拟变得相当困难。因此，为了提升汽车油耗理论模型的实用性和简便性，本研究忽略了道路坡度、路面状况、司机驾驶行为因素、克服重力、加减速度等能耗，加速坡度及辅助功率均取值为0，简化后公路运输燃油消耗模型如式(13)所示。

$ e_{i j}^k=\frac{\varepsilon d_{i j} \delta N V_{\mathrm{r}}}{\gamma v}+\frac{\varepsilon d_{i j} 0.5 C_{\mathrm{d}} \rho S v^2}{1\;000 \gamma \eta \lambda}+\frac{\varepsilon d_{i j} M g C_{\mathrm{r}}}{1\;000 \gamma \eta \lambda} 。$

(13)

采用国际政府间气候变化专门委员会(IPCC)发布的《国家温室气体排放指南》中的基于行驶距离的碳排放计算方法，准确表示铁路运输的燃油消耗情况，确保碳排放量的估算更为科学、合理。

$ e_{i j}^k=e^k \cdot d_{i j} \cdot q。$

(14)

多式联运中的水路运输燃油消耗不仅受到船舶主机功率的影响，还与辅助设备的功率、船舶的实际行驶速度及设计速度密切相关。

$ e_{i j}^k=\left[p_\alpha\left(\frac{v_{i j}}{v_0}\right)^3+p_{\mathrm{b}}\right] \cdot \frac{d_{i j}}{24 v_0} 。$

(15)

(2) 转运碳排放量

转运过程中的碳排放量为转运量与单位转运碳排放量的乘积。

$ E_{\mathrm{z}}=\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} E_i^{k l} \cdot y_i^{k l}=\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} e_{k l} \cdot q_{i j} \cdot y_i^{k l} 。$

(16)

(3) 冷藏碳排放量

冷藏过程中的碳排放量为单位冷藏碳排放量、运货量、运输距离的乘积。

$ E_{\mathrm{h}}=\sum\limits_{i \in M} \sum\limits_{j \in M} \sum\limits_{k \in K} e_{\mathrm{h}} \cdot q_{i j} \cdot d_{i j}^k \cdot x_{i j}^k。$

(17)

企业的碳排放成本为碳排放量及碳税税率相乘得出。

$ C_{\mathrm{e}}=\left(E_{\mathrm{y}}+E_{\mathrm{z}}+E_{\mathrm{h}}\right) \cdot \gamma_{0}。$

(18)

2.3 模型建立

建立冷藏集装箱多式联运路径成本模型：式(19)为总成本；式(20)为相邻节点间相同的运输方式最多只有一种；式(21)规定了节点处的货物转换次数上限为一次；式(22) 确保了节点处的中转方式与前后运输方式相匹配；式(23)为决策变量，表示如果节点i处中转方式从k转换成l，y_i^kl为1，否则为0；式(24)为决策变量，表示如果货物从节点i以运输方式k运输到节点j，x_{i, j}^k为1，否则为0。

$ \min C=C_{\mathrm{y}}+C_{\mathrm{z}}+C_{\mathrm{t}}+C_{\mathrm{h}}+C_{\mathrm{e}}, $

(19)

$ \sum\limits_k x_{i, j}^k \leqslant 1, $

(20)

$ \sum\limits_k \sum\limits_l y_i^{k l} \leqslant 1, $

(21)

$ \sum\limits_{h, i \in M} x_{h, i}^k+\sum\limits_{i, j \in M} x_{i, j}^l \geqslant 2 y_i^{k, l}, $

(22)

$ y_i^{k l} \in\{0, 1\}, $

(23)

$ x_{i, j}^k \in\{0, 1\} 。$

(24)

2.4 鲁棒优化模型

为解决需求不确定问题，采用基于情景分析的鲁棒优化方法，构建了针对需求不确定性的冷藏集装箱多式联运优化模型。该模型旨在通过情景模拟和鲁棒性设计，提高多式联运路径规划的灵活性和适应性，以应对市场需求的波动变化。

$ \begin{gathered} \min C=\sum\limits_{s=1}^S p_{\mathrm{s}} \cdot C_{\mathrm{s}}(x)= \\ \sum\limits_{s=1}^S p_{\mathrm{s}} \cdot\left(C_{\mathrm{y}}+C_{\mathrm{z}}+C_{\mathrm{t}}+C_{\mathrm{h}}+C_{\mathrm{e}}\right), \end{gathered} $

(25)

式中C_s为运输情景s中对应的运输费用。

模型需要满足如下约束条件，同时式(20)~(24)成立。

$ C_{\mathrm{s}}(x) \leqslant(1+\beta) C_{\mathrm{s}}^*, $

(26)

$ \sum\limits_{s=1}^S p_{\mathrm{s}}=1, $

(27)

式中C_s^*为运输情景s中最优解对应的最小运输费用。

3 改进樽海鞘算法

樽海鞘算法是一种模拟樽海鞘群体链式觅食行为的元启发式优化算法，通过领导者和追随者的协作机制来求解复杂优化问题。为在规定时间内求解冷藏集装箱多式联运最优方案，设计了改进的樽海鞘算法对模型进行求解。首先，对算法的编码方式进行了改进，以更好地适应离散路径优化问题的特点；其次，引入自适应种群动态分配数量结构；最后，引入了混沌映射机制，对领导者和追随者的位置更新公式进行改进，求解步骤如图 1所示。

3.1 编码方式

任意樽海鞘个体编码包含运输路径节点的编码区域及不同运输方式的编码区域。n为运输节点个数，n－1为运输方式总数。

$ X=\left[x_1, x_2, \cdots, x_n, x_{n+1}, \cdots, x_{2 n-1}\right] 。$

(28)

3.2 自适应动态种群分配结构

传统樽海鞘算法由1个领导者和n个追随者组成，固定的种群结构会导致全局搜索能力不足，同时局部搜索容易陷入局部最优。为解决这一问题，本研究提出自适应动态种群分配结构。使领导者和追随者的种群数量比例随着迭代次数增多而变化，前期大量领导者进行全局探索，后期大量追随者进行细致的局部搜索，从而优化了算法的整体精确度。

$ \begin{gathered} r=r_1-\left(r_1-r_2\right) \cdot \\ \left\{b \cdot \tan \left[\frac{\mathtt{π} }{4}\left(1-\frac{t}{t_{\max }}\right)\right]-e \cdot R\right\}, \end{gathered} $

(29)

式中，r为领导者种群数量比例，因此1－r为追随者种群数量比例；r₁，r₂，b，e为控制参数，分别取值0.70，0.30，0.75，0.20；t为当前迭代次数；t_max为最大迭代次数；R为随机函数。

3.3 领导者位置更新公式

传统的樽海鞘算法依赖领导者引导整个种群移动，但这种单一的领导者搜索策略限制了樽海鞘链的搜索能力。为减少搜索的盲目性并增强算法的全局搜索能力，在领导者位置更新策略中引入共享机制。同时，为了增加领导者种群的多样性，还引入了混沌映射方法。本研究采用Tent映射。

$ r_k=\left\{\begin{array}{l} r_{k-1} / \lambda, r_{k-1} \in(0, \lambda] \\ \left(1-r_{k-1}\right) /(1-\lambda), r_{k-1} \in(\lambda, 1] \end{array}, \right. $

(30)

式中，r为Tent映射混沌数；λ为参数，若η取值0.5，系统容易演化为周期性质，从而失去混沌状态。为了解决这一问题，在Tent映射基础上引入随机变量，以维持系统的混沌特性。

$ r_k^*=\left\{\begin{array}{l} r_{k-1} / \lambda+R \quad r_{k-1} \in(0, \lambda] \\ \left(1-r_{k-1}\right) /(1-\lambda)+R \quad r_{k-1} \in(\lambda, 1] \end{array}。\right. $

(31)

领导者位置更新如式(32)所示。当混沌数小于0.5时，当前个体与食物源进行部分映射交叉操作(PMX); 当混沌数不小于0.5时，当前个体与上次迭代的当前个体进行部分映射交叉操作。交叉操作如图 2所示。

$ F_i(t+1)=\left\{\begin{array}{l} \operatorname{PMX}\left[F_i(t),\; F_\alpha(t)\right] \quad r_k^*<0.5 \\ \operatorname{PMX}\left[F_i(t), \;F_i(t-1)\right] \quad r_k^* \geqslant 0.5 \end{array}, \right. $

(32)

式中F_i和F_α分别为当前个体和食物源。

3.4 追随者位置更新公式

传统樽海鞘算法中，追随者仅与自身及相邻个体进行信息交换，这可能导致算法陷入局部最优解。为了克服这一局限，加强了追随者与当前最优个体之间的信息交流，从而有效避免算法过早收敛。同时，采用了Bernoulli映射方法，增加种群的多样性。

$ r_{k+1}=\left\{\begin{array}{ll} r_k /(1-\delta), & r_k \in(0, 1-\delta] \\ \left(r_k-1+\delta\right) / \delta, & r_k \in(1-\delta, 1] \end{array}, \right. $

(33)

式中δ为混沌参数，一般取值为0.4。追随者位置更新如式(34)所示。当δ＜0.5时，采用当前个体与相邻个体进行部分映射交叉操作，当δ≥0.5时，采用当前个体与最优个体进行部分映射交叉操作。

$ F_i(t+1)=\left\{\begin{array}{ll} \mathrm{PMX}\left[F_i(t), \;F_{i-1}(t)\right], \;\;r_{k+1}<0.5 \\ \mathrm{PMX}\left[F_i(t), \;F_i^*(t)\right], \;\;r_{k+1} \geqslant 0.5 \\ \end{array}, \right. $

(34)

式中*代表最优。

4 算例分析 4.1 算例数据

为验证改进樽海鞘算法求解所构建模型的有效性，以中国某地区多式联运网络为例，构建网络结构如图 3所示。

一批不确定数量的货物从始发地运输到目的地，需要在60~90 h内完成，否则提前到达产生的单位仓储费用为50元/(h·TEU)，单位惩罚成本为95元/(h·TEU)。冷藏集装箱的单位冷藏成本为12元/(h·TEU)，单位货损成本为6元/(h·TEU)，单位冷藏碳排放量为0.2 kg/h，碳税单价为20元/kg。实际交通量和设计交通量比值范围为[0, 3]。不同运输方式的速度和运输基价数据如表 2所示。不同运输方式转换数据如表 3所示。

4.2 模型求解及结果分析

所有算法均采用MatlabR2020a编程实现，算法中的计算参数设置为种群规模N=30，最大迭代次数T_max=200。构建5种不同运输情景，各货运量(集装箱数量)分别为15，17，19，21，23个。

得到的鲁棒优化解为6.874 5万元。改进樽海鞘算法求解模型的收敛情况如图 4所示。情景最优解和鲁棒最优解对应的总成本对比结果如表 4所示。

由表 4可知：(1)鲁棒优化路线对应的各情景的到达时间均在客户允许的范围内，这表明该鲁棒优化路线具有可行性。(2)鲁棒优化路线对应的总成本不小于确定情景对应的总成本，这说明鲁棒路线虽然稳定，抗风险能力强，但是企业要额外付出更多的成本。(3)成本偏差均小于设定的遗憾系数，这表明该鲁棒最优解满足企业客户的需求，所构建的模型具有鲁棒性。

4.3 改进冷藏集装箱多式联运模型的检验

对优化模型的先进性进行验证，结果如表 5所示。在未考虑路网阻抗的情境下，模型倾向于选择成本较低的路径，这可能导致在高峰时段遇到延误。而优化后的模型则能够识别并优先选择成本稍高但能规避拥堵的路线，以确保货物能够准时送达。通过对比2种模型的测试结果可知，考虑路网阻抗的改进模型在准时性、服务可靠性方面均优于传统模型，这对于提升客户满意度至关重要。

对优化模型的鲁棒性进行验证。鲁棒优化解与遗憾系数具有关联性，如图 5所示。当遗憾系数较小时，客户的要求较高，成本与碳排放期望值较大，模型确保了可行性，而在低遗憾系数时，智能避免无解情况。随着系数的增加，成本与碳排放期望值均呈现减小的趋势，并在达到最优值后稳定提供鲁棒解。这体现了优化模型在适应性和优化效率上的先进性。

4.4 改进樽海鞘算法的检验

为了验证改进算法中各个优化步骤的有效性，进行了不同优化策略算法的对比试验。在其他条件不变的情况下，将本研究改进算法、未引入Tent映射共享机制的领导者位置更新策略、未引入自适应动态种群分配结构和未引入Bernoulli映射共享机制的追随者位置更新策略分别定义为PMSSA，PMSSA-1，PMSSA-2，PMSSA-3。算法收敛曲线如图 6所示。

由图 6可知，不同改进樽海鞘算法策略求解目标函数所得结果有显著差异。本研究所提出的改进樽海鞘算法在前期展现出卓越的全局搜索能力，而到了后期，其收敛速度得到了显著的提升。这说明引入混沌映射与交叉算子、自适应动态种群分配结构策略对提升樽海鞘算法寻优能力具有显著优势。改进领导者位置更新策略增强算法的全局搜索能力，改进追随者位置更新策略使算法跳出局部最优，防止过早收敛。

5 结论

冷藏集装箱多式联运与普通多式联运存在显著差异，主要体现在其运输全过程都处于严格的低温条件下，其产生的运输成本和碳排放量远高于普通多式联运。为达到经济效益和低碳效益最大化目标，研究冷藏集装箱多式联运路径优化问题显得尤为重要。本研究在考虑需求不确定和路网阻抗等因素的情况下，研究冷藏集装箱多式联运路径优化问题，建立运输总成本优化模型，改进传统樽海鞘算法，并应用于某多式联运网络。

研究结果表明：(1)本研究的鲁棒优化解能够达到客户的要求，具有鲁棒性，使用所构建的模型能够得到较为满意的解决方案。(2)成本期望值随着遗憾系数的增大而减小，客户可根据对鲁棒的需求来选择合适的运输路线和运输方式。(3)在传统樽海鞘算法中引入领导者与追随者位置更新策略和自适应动态种群分配结构能显著提高算法的寻优能力。(4)所构建的改进模型在服务准时性、鲁棒性方面均优于传统模型，不仅提升了物流效率，还确保了货物交付的高可靠性。