考虑两阶段变质的低碳冷链物流库存–路径联合优化研究

Joint Optimization Study on Inventory-Routing for Low-Carbon Cold Chain Logistics with Two-Stage Deterioration

摘要: 针对低碳冷链物流库存–路径问题，本研究考虑冷链产品的库存和配送两阶段变质情况及碳排放，计算货损成本、碳排放成本及其它各项成本；同时在供应商库存管理模式下，分析零售商的库存状态，并计算相关成本。以系统总成本最小化为目标，建立考虑两阶段冷链产品变质的低碳库存–路径问题模型(Low-Carbon Inventory-Routing Problem, LCIRP)，并设计了一种采用大规模邻域搜索中的“破坏”和“修复”思想的两阶段遗传算法进行求解。最后，通过对比有无碳排放约束的两种模型和是否考虑库存路径联合优化的两种模式的求解结果，以及对碳交易价格的灵敏度分析，算例结果表明LCIPR模型在碳排放成本和系统总成本方面有较大的优势，能同时兼顾企业经济效益和社会环境效益，更好地实现双赢局面。

Abstract: For the low-carbon cold-chain logistics inventory-routing problem, this study considers the inventory and distribution of cold-chain products in two stages of deterioration and carbon emissions, calculates the cost of damage, carbon emissions and other costs; while in the supplier inventory management mode, it analyzes the inventory status of retailers and calculates relevant costs. Taking the minimum total system cost as the goal, a low-carbon inventory-routing problem model (LCIRP) considering two-stage cold-chain product deterioration is established, and a two-stage genetic algorithm using large-scale neighborhood search with “destruction” and “repair” ideas is designed for solution. Finally, through comparing the results of two models with and without carbon emission constraints, as well as two modes with and without considering the joint optimization of inventory and routing, and a sensitivity analysis of carbon trading prices, the example results demonstrate that the LCIPR model has significant advantages in terms of carbon emission cost and total system cost, which can balance both economic benefits and social environmental benefits of enterprises, thus better achieving a win-win situation.

1. 引言

随着人民生活水平的不断提高，冷链物流行业在近年来得到了快速发展，我国对冷链产品的需求也呈现出日益增长的趋势。然而，由于冷链产品具有易变质的特性，因此在物流运输过程中容易产生损耗。这些损耗主要集中在库存和配送两个环节之间，而这两个环节之间又存在着效益悖反的关系。因此，对冷链产品的库存路径问题进行联合优化具有重要意义[1]。与普通产品不同，冷链产品在库存和配送过程中需要进行低温处理，这需要消耗更多的燃料或电力来维持产品的新鲜度和质量，从而产生更多的碳排放，对环境造成恶劣影响。目前，低碳环保问题已经成为全球共同关注的焦点，我国作为一个能源消耗大国，也面临着严峻的碳排放挑战[2]。故本文在冷链与低碳的双重背景下，针对物流系统中的库存–路径问题(Inventory Routing Problem, IRP)进行研究，旨在通过探讨库存路径联合优化的方法，降低冷链物流系统的总成本，同时降低碳排放，为我国冷链物流行业的发展提供新的思路和方法。

库存路径问题是供应链管理中一个重要的问题，旨在解决库存与配送路径之间的协调问题。在考虑库存决策和配送路径决策的前提下，库存路径问题旨在确定各个零售商的订货方案以及配送路径方案，以实现整个供应链系统的优化。自1984年Federgruen [3]首次提出将库存和路径联合优化的思想以来，库存路径问题引起了广泛关注，越来越多的学者开始研究库存路径问题，并引入了许多现实因素，以解决更加复杂的问题。其中，Rusdiansyah [4]等将时间窗约束加入到库存路径问题当中，该论文考虑了每个零售商在特定时间段内对商品的需求，如果订单在需求时间段内未到达，将导致一定比例的需求无法满足。这种时间窗约束对于现实供应链系统来说非常重要，因为无法满足客户需求可能会导致销售损失。另外，Mirzaei [5]等研究了考虑损失成本的库存路径问题，该论文将需求的一定比例视为损失销售，并考虑了由于缺货而导致的销售损失。这种考虑损失成本的研究对于优化供应链系统中的库存和配送决策具有重要意义。

近年来，冷链物流库存–路径问题的相关研究引起许多学者的关注。黎莎等[6]认为在冷链物流中货物新鲜度和送达时间是提高服务效率，降低成本的关键，为提高冷链物流库存路径的组织协调性，该论文考虑新鲜度在流通过程中的动态变化会导致配送方案和库存分配的差异。梁承姬等[7]以配送中心为研究对象，在配送之前根据产品初始新鲜度将产品分成不同等级，采取等级向上的替换模式来提升产品整体新鲜度，并确保满足客户需求，同时针对配送和库存分配问题进行了联合优化，建立了双层规划模型进行决策，并提出改进的自适应遗传算法对其进行求解。Rohmer等[8]提出了一个冷链产品的两级库存路径问题，产品从供应商运送到中间仓库储存，再由车辆送至客户处，库存持有成本为仓库存储产生，目标是降低总运输和持有成本，并采用自适应大邻域搜索元启发式方法结合简化公式进行求解。Onggo等[9]以库存、运输、货损和缺货成本为成本函数，最小化系统的总成本为目标函数，建立了随机需求冷链库存路径模型，并提出了一种将蒙特卡罗模拟与迭代局部搜索相结合的相似启发式算法来求解。Ji等[10]针对具有时间窗约束的冷链产品库存路径问题，建立了以总成本最小化为目标函数的混合整数线性规划模型，并通过引入不确定性集，将模型进一步转化为混合整数鲁棒规划模型。

随着绿色物流的发展，国内外学者在研究库存–路径问题时进一步考虑运输过程中的碳排放。Gutierrez等[11]在低碳随机库存路径问题中引入了电动道路问题，其中混合动力汽车在道路网络中导航并有充电机会，以满足零售商对单一产品的随机需求，使用等时图表示道路网络，允许跟踪车辆电池电量，同时制定了一个数学规划启发式算法，并在现实问题实例上进行了测试，展示了不同策略的可行性。Dezhtaherian等[12]研究了二级供应链中冷链产品的低碳库存路径问题，包括分销商和零售商，产品有固定到期日期，价格随时间变化，需求固定，考虑了外包仓库和车辆来管理库存，问题以混合整数线性规划模型表述，并使用OPL软件和组合蚁群算法求解。王梦梦等[13]在将配送中心选址、配送路径和库存策略相结合的基础上引入了碳排放约束，建立了随机需求的冷链物流选址–路径–库存联合优化模型，并提出一种两阶段启发式算法来求解。Shamayleh等[14]优化了零售商的补货策略，针对冷产品的离散时变需求和碳排放，建立了三个优化模型，其目标函数分别为成本最小化、碳排放量最小化以及成本和碳排放量联合最小化，并通过确定三种模型的最优解的结构性质，提出基于动态规划的精确解算法。Bai等[15]在供应商库存管理模式下，以一个供应商和两个竞争零售商为研究对象，提出了一个考虑碳限额约束的中心化系统的优化模型。

综上所述，尽管已有学者对冷链物流库存–路径问题和低碳库存–路径问题都有所研究，但仍存在以下研究缺口：1) 在低碳冷链物流库存–路径问题方面，虽然有些文献研究了低碳冷链物流配送的车辆路径问题，但很少有文献深入探讨低碳冷链物流库存–路径问题。这是一个有待进一步研究的重要领域，因为库存和路径决策之间存在着密切的联系，低碳排放需要考虑库存和路径两方面的优化。2) 在计算冷链产品的损耗成本和整个系统中的碳排放时，许多文献只考虑了运输过程中的损耗成本和碳排放，而没有考虑到库存和配送两阶段的损耗成本和碳排放。这是一个重要的研究缺口，因为在冷链物流中，库存和配送阶段同样重要，它们对整个供应链系统的效率和碳排放都会产生影响。3) 现有的文献还缺乏对车辆使用成本、车辆油耗和碳排放、配送和库存环节的制冷系统的碳排放、时间窗惩罚成本以及库存成本的全面考虑。这些因素都对冷链物流系统的优化和碳排放的减少有着重要影响。故综合考虑这些因素，可以更全面地评估冷链物流系统的效率和碳排放。因此，本文在库存–路径问题的基础上，引入冷链产品的库存和配送两阶段变质，同时考虑库存和配送两个环节的碳排放；并考虑时间窗、车容量限制等约束条件，以系统总成本最小化作为目标函数，构建LCIRP模型。此外，本文还提出了一种改进的两阶段遗传算法，在第二阶段遗传算法中引入大规模邻域搜索操作。最后，通过算例仿真实验对比分析库存–路径联合优化和单阶段优化两种模式，为冷链物流企业的未来发展方向提供理论基础；研究碳价格变化对碳排放和企业系统总成本的影响，为政府制定碳交易相关政策提供参考。

2. 问题描述与模型建立

2.1. 问题描述

低碳冷链物流库存–路径问题可以描述为：一个冷链物流配送中心为N个零售商配送一种冷链产品，配送中心与各个零售商的位置和距离信息已知。配送中心配备K辆相同型号的冷藏车，每辆冷藏车的最大容量为M。每个零售商的需求服从正态分布，并都配备相同型号的冷冻机，每个冷冻机的最大容量为b。由于零售商每天的需求服从正态分布，配送中心每天都会为每个零售商进行一次配送，当零售商出现缺货时，会产生缺货成本；产品储存在零售商仓库中，会产生库存持有成本。考虑冷链产品在运输和储存过程中会发生变质，并且会产生碳排放，本文将其转化为成本加入到成本函数中。其目标是在满足车容量限制以及零售商需求和时间窗限制的前提下，制定配送方案，使整个系统成本最小化。

本研究同时考虑供应商库存管理模式，该模式通过更好地协调供应商和零售商之间的库存管理来减少库存成本和浪费，并有效提高物流运输效率[16]。在该模式下，冷链物流配送中心每天根据各个零售商的需求和库存状况来制定补货策略，从而确定每个零售商的配送量以及车辆的配送计划。零售商每天收到产品后才开始销售，当零售商的库存量小于当天的需求时，便会产生缺货成本；产品存储在零售商的仓库中，会产生库存持有成本；在配送和存储过程中，产品会发生变质，以及冷藏车和冷冻机会产生碳排放，从而会造成货损成本和碳排放。通过分析，该问题具有两阶段性，第一阶段在零售商的需求发生前，冷链物流配送中心根据零售商的需求和当前库存量来制定配送方案，由此产生第一阶段与配送相关的成本，第二阶段在零售商需求发生后，确定零售商的库存量和缺货量，从而确定与库存方面相关的成本，因此，本文将该问题建模为两阶段模型。

2.2. 符号定义

以下是符号、集合、变量的定义：

$E$ ：计划期的时间集合$E=\left\{1,2,\cdots ,L\right\}$

$A$ ：配送节点集合$A=\left\{0,1,2,\cdots ,N\right\}$ ，0表示配送中心

$A^{\prime }$ ：零售商集合$A^{\prime }=A\backslash \left\{0\right\}$

$B$ ：配送车辆集合$B=\left\{1,2,\cdots ,K\right\}$

$h_i$ ：零售商i的单位重量单位时间的产品库存持有成本(元/kg*天)

$p$ ：产品单位重量的价格(元/kg)

$\lambda$ ：储存过程中产品的变质率

$r_{it}$ ：零售商i第t天的需求量

$I_{it}$ ：零售商i第t天的起始库存量

$d_{ijk}^t$ ：第t天时车辆k在道路$\left(i,j\right)$ 上的行驶距离

$M$ ：最大车容量

$w_{ijk}^t$ ：第t天时车辆k在道路$\left(i,j\right)$ 上的载重

$w_0$ ：车辆的自重

${\tau }_{si}^t$ ：零售商i在第t天的服务时间

${\tau }_{ijk}^t$ ：第t天时车辆k在道路$\left(i,j\right)$ 上行驶的时间

${\tau }_{ik}^t$ ：第t天时车辆k到达节点i的时间点

${\tau }_0$ ：车辆从配送中心出发的时间

$\left[e_0,l_0\right]$ ：配送中心的时间窗

$\left[e{\tau }_i^t,l{\tau }_i^t\right]$ ：零售商i第t天的服务时间窗

$F_{ijk}$ ：车辆k从节点i行驶到节点j产生的油耗量(L)

$f_{ijk}$ ：车辆k从节点i到节点j因制冷产生的油耗量(L)

$\psi$ ：燃油排放系数

$\phi$ ：车辆的固定发车费用(元/辆)

${\sigma }_1$ ：车辆在行驶时的制冷系统油耗率(L/h)

${\sigma }_2$ ：车辆在服务时的制冷系统油耗率(L/h)

$C_f$ ：单位油耗费用(元/L)

$C_e$ ：单位碳排放费用(元/kg)

$v$ ：车辆行驶速度(km/h)

$b$ ：每台冷冻机的容量

$\eta$ ：一台冷冻机的单位时间电能消耗(kw/台*天)

$\theta$ ：一度电所造成的碳排放总量(kg/kw)

${\omega }_e$ ：早到零售商的时间窗惩罚系数(元/h)

${\omega }_l$ ：晚到零售商的时间窗惩罚系数(元/h)

$Q_{it}$ ：零售商i在第t天的配送数量

变量：

$S_{it}$ ：零售商i在第t天的实际配送数量

$T_i$ ：零售商i的配送周期

$x_{ijk}^t$ ：0~1变量，第t天车辆$k$ 是否由节点i驶向节点j，如果是则为1，否则为0

$y_{ik}^t$ ：0~1变量，第t天零售商i是否由车辆k服务，如果是则为1，否则为0

2.3. 模型建立

a) 目标函数

供应商库存管理模式下，零售商的需求和库存状况可与供应商进行共享，零售商无需进行订货，从而也无订货提前期，供应商根据零售商的共享信息制定补货策略，确定配送计划。本文建立的LCIRP模型考虑的是供应商库存管理模式下，优化目标是系统总成本最小化，其包括库存环节的库存成本$C_1$ 和缺货成本$C_2$ ，配送环节的车辆使用成本$C_4$ ，油耗成本$C_5$ 和时间窗惩罚成本$C_7$ ，以及库存和配送两阶段的货损成本$C_3$ 和碳排放成本$C_6$ 。目标函数如式(1)所示。

$Z=C_1+C_2+C_3+C_4+C_5+C_6+C_7$ (1)

b) 成本函数

1) 零售商的库存持有成本$C_1$

零售商在对产品进行储存的过程中，会产生库存持有成本，其大小与库存量和储存时间有关，库存量越大和储存时间越长，库存持有成本则会越大。由于本文考虑的需求是服从正态分布的多周期随机需求，零售商的库存是实时变化的，每周期的期初库存和期末库存都有可能不相同，因此本文采用每周期的平均库存来计算零售商的库存持有成本。库存持有成本如式(2)所示。

$C_1={\displaystyle \sum _{t=1}^L{\displaystyle \sum _{i=1}^N\frac{I_{it}+I_{i,t+1}+Q_{it}}{2}h}}$ (2)

其中期末库存等于下一周期的期初库存$I_{i,t+1}$ ，其大小由期初库存、需求量和配送量决定，如式(3)所示。

$I_{i,t+1}=I_{it}+Q_{it}-r_{it}$ (3)

2) 缺货成本$C_2$

缺货成本是指供应商无法满足零售商的需求，从而产生的成本。缺货成本由零售商缺少的产品数量决定，与期初库存、需求量、配送量相关，如果当天需求量小于当天期初库存量与配送量之和，则会发生缺货成本。缺货成本如式(4)所示。

$C_2={\displaystyle \sum _{t=1}^L{\displaystyle \sum _{i=1}^N\max \left\{r_{it}-I_{it}-Q_{it},0\right\}p}}$ (4)

3) 货损成本$C_3$

冷链产品在配送和储存的过程中，产品都会发生损坏，从而产生货损成本。本文采用文献[17]提出的新鲜度衰减函数$q\left(t\right)=q_0{\text{e}}^{-\lambda t}$ 来计算配送和库存两阶段的货损成本，其中$q\left(t\right)$ 是t时刻的新鲜度，$q_0$ 是初始新鲜度。由于供应商是根据零售商的需求和库存状况制定的配送计划，考虑到配送过程中可能会导致部分冷链产品损坏，其实际发出的冷链产品数量应该要大于零售商的计划配送量。当变质率为$\lambda$ 时，则供应商第t天对零售商i的实际配送量$S_{it}$ 与第t天零售商i的配送量$Q_{it}$ 的关系如式(5)所示。

$S_{it}=\frac{Q_{it}}{1-\lambda }$ (5)

储存过程中的货损成本$C_3^s$ 如式(6)所示。

$C_3^s={\displaystyle \sum _{t=1}^L{\displaystyle \sum _{i=1}^{N}p\frac{I_{it}+I_{i,t+1}+Q_{it}}{2}\left(1-{\text{e}}^{-24\lambda }\right)}}$ (6)

配送过程中的货损成本$C_3^d$ 如式(7)所示。

$C_3^d={\displaystyle \sum _{k=1}^K{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{t=1}^{L}p{S}_{it}\left(1-{\text{e}}^{-\lambda ({\tau }_{ik}^t-{\tau }_0)}\right)y_{ik}^t}}}={\displaystyle \sum _{k=1}^K{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{t=1}^L\frac{p{Q}_{it}}{1-\lambda }\left(1-{\text{e}}^{-\lambda ({\tau }_{ik}^t-{\tau }_0)}\right)y_{ik}^t}}}$ (7)

总货损成本为如式(8)所示。

$C_3=C_3^s+C_3^d$ (8)

4) 车辆使用成本$C_4$

车辆在配送过程中会产生车辆使用成本$C_4$ ，本文主要考虑的是车辆的固定发车费用。车辆配送过程中的使用成本如式(9)所示。

$C_4={\displaystyle \sum _{k=1}^K{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{t=1}^L{x}_{0jk}^t\phi }}}$ (9)

5) 油耗成本$C_5$ 和碳排放成本$C_6$

车辆在配送过程中的油耗主要由行驶过程中的油耗和配送过程中制冷系统所需的油耗组成，本文采用Barth等[18]提出的模型，来计算车辆在行驶过程中的油耗。而在配送过程中制冷系统产生的油耗是由行驶过程中制冷系统产生的油耗和服务期间制冷系统产生的油耗组成，一般情况下由于服务时要打开车门进行卸货，会导致车厢内温度升高，从而车辆在行驶时的油耗率要小于服务时的油耗率，即${\sigma }_1<{\sigma }_2$ [19]。其行驶过程中的油耗公式和制冷系统的油耗公式分别如(10)和(11)所示，其中${\beta }_1$ 为发动机模块系数，${\beta }_2$ 速度模块系数，${\beta }_3$ 载重模块系数。

$F_{ijk}^t=\frac{{\beta }_1{d}_{ijk}^t}{v}+{\beta }_2{d}_{ijk}^t{v}^2+{\beta }_3\left(w_0+w_{ijk}^t\right)d_{ijk}^t$ (10)

$f_{ijk}^t=\frac{{\sigma }_1{d}_{ijk}^t}{v}+{\sigma }_2{\tau }_{si}^t$ (11)

车辆在配送过程中的油耗成本如式(12)所示。

$C_5=C_f{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{j=0}^N{\displaystyle \sum _{i=0}^N{\displaystyle \sum _{t=1}^L{F}_{ijk}^t{x}_{ijk}^t}}}}+C_f{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{j=0}^N{\displaystyle \sum _{i=0}^N{\displaystyle \sum _{t=1}^L{f}_{ijk}^t{x}_{ijk}^t}}}}$ (12)

冷链产品在配送过程中和储存过程中都会产生碳排放，储存过程中的碳排放主要由冷冻机产生，本文通过计算计划期内所使用的冷冻机数量，从而根据冷冻机的数量、冷冻机单位时间的电能消耗以及每度电所造成的碳排放量来计算储存过程中的碳排放成本。根据文献[20]表明，车辆行驶过程中的碳排放与油耗量成正比，油耗量乘以燃油排放系数$\psi$ 即可得到碳排放量。

储存过程中的碳排放成本$C_6^s$ 如式(13)所示。

$C_6^s={\displaystyle \sum _{t=1}^L{\displaystyle \sum _{i=1}^N\frac{I_{it}+I_{i,t+1}+Q_{it}}{2b}\eta \theta C_e}}$ (13)

配送过程中的碳排放成本如式(14)所示。

$C_6^d=C_e{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{j=0}^N{\displaystyle \sum _{i=0}^N{\displaystyle \sum _{t=1}^L\psi F_{ijk}^t{x}_{ijk}^t}}}}+C_e{\displaystyle \sum _{k=1}^K{\displaystyle \sum _{j=0}^N{\displaystyle \sum _{i=0}^N{\displaystyle \sum _{t=1}^L\psi f_{ijk}^t{x}_{ijk}^t}}}}$ (14)

总碳排放成本如式(15)所示。

$C_6=C_6^s+C_6^d$ (15)

6) 违背时间窗惩罚成本$C_7$

由于冷链产品的易变质特性，客户对其配送时间的要求会比一般产品更加严格。因此本文通过引入时间窗惩罚成本来更好地满足客户对配送时间的要求，如果早到则会产生一个等待成本，晚到则会产生一个客户对其的惩罚成本。时间窗惩罚成本为如式(16)所示。

$C_7={\omega }_e{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{t=1}^L\max \left\{e{\tau }_i^t-{\tau }_{ik}^t,0\right\}}}+{\omega }_l{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{t=1}^L\max \left\{{\tau }_{ik}^t-l{\tau }_i^t,0\right\}}}$ (16)

c) 模型建立

第一阶段在零售商发生需求之前，确定配送量，使用的冷藏车数量，每辆车的配送路径。第一阶段模型为：

$\min Z_1=C_3^d+C_4+C_5+C_6^d+C_7$ (17)

$s.t.{\displaystyle \sum _{i=1}^N{Q}_{it}{y}_{ik}^t}\le M,\forall t\in E,\forall k\in B$ (18)

${\displaystyle \sum _{k=1}^K{y}_{ik}^t}=1,\forall i\in A^{\prime },\forall t\in E$ (19)

${\displaystyle \sum _{j=1}^N{x}_{0jk}^t}={\displaystyle \sum _{j=1}^N{x}_{j0k}^t}\le 1,\forall k\in B,\forall t\in E$ (20)

${\displaystyle \sum _{k=1}^K{\displaystyle \sum _{j=1}^N{x}_{0jk}^t}}={\displaystyle \sum _{k=1}^K{\displaystyle \sum _{i=1}^N{x}_{i0k}^t}}\le K,i\ne j,\forall t\in E$ (21)

${\tau }_{0k}^t\le l_0,\forall k\in B,\forall t\in E$ (22)

${\tau }_{jk}^t={\tau }_{ik}^t+{\tau }_{si}^t+{\tau }_{ijk}^t,\forall i,j\in A^{\prime },\forall t\in E,\forall k\in B,i\ne j$ (23)

式(17)表示配送过程中的成本最小化；式(18)表示冷藏车每次的配送量不大于车辆最大容量；式(19)表示每个零售商每次最多只能由一辆冷藏车进行配送；式(20)表示车辆每次配送只有一条路径，且都是从供应商出发，最后返回供应商；式(21)表示离开供应商的车辆数等于返回供应商的车辆数；式(22)表示车辆最后返回供应商的时间点要早于供应商最后关闭的时间点；式(23)表示车辆到达下一节点的时间点等于车辆到达上一节点的时间点、上一节点的服务时间以及车辆在两节点之间的行驶时间之和。

第二阶段在零售商需求发生之和后，确定配送量之后，发生库存持有成本和缺货成本。第二阶段模型为：

$\min Z_2=C_1+C_2+C_3^s+C_6^s$ (24)

$I_{i,t+1}=I_{it}+Q_{it}-r_{it},\forall i\in A^{\prime },\forall t\in E$ (25)

$I_{it},I_{i,t+1},Q_{it},r_{it}\ge 0$ (26)

式(24)表示储存过程中的成本最小化；式(25)表示每个零售商下一周期的起初库存等于该周期的期末库存；式(26)表示规定各库存变量的取值范围。

以系统总成本最小化作为优化目标，建立如下模型：

$\min Z=C_1+C_2+C_3+C_4+C_5+C_6+C_7$ (27)

$s.t.{\displaystyle \sum _{i=1}^N{Q}_{it}{y}_{ik}^t}\le M,\forall t\in E,\forall k\in B$ (28)

${\displaystyle \sum _{k=1}^K{y}_{ik}^t}=1,\forall i\in A^{\prime },\forall t\in E$ (29)

${\displaystyle \sum _{j=1}^N{x}_{0jk}^t}={\displaystyle \sum _{j=1}^N{x}_{j0k}^t}\le 1,\forall k\in B,\forall t\in E$ (30)

${\displaystyle \sum _{k=1}^K{\displaystyle \sum _{j=1}^N{x}_{0jk}^t}}={\displaystyle \sum _{k=1}^K{\displaystyle \sum _{i=1}^N{x}_{i0k}^t}}\le K,i\ne j,\forall t\in E$ (31)

${\tau }_{0k}^t\le l_0,\forall k\in B,\forall t\in E$ (32)

${\tau }_{jk}^t={\tau }_{ik}^t+{\tau }_{si}^t+{\tau }_{ijk}^t,\forall i,j\in A^{\prime },\forall t\in E,\forall k\in B,i\ne j$ (33)

$I_{i,t+1}=I_{it}+Q_{it}-r_{it},\forall i\in A^{\prime },\forall t\in E$ (34)

$I_{it},I_{i,t+1},Q_{it},r_{it}\ge 0$ (35)

3. 算法设计

本文构建的LCIRP模型是一个两阶段模型，包括库存阶段和路径优化阶段。对于这类模型的求解，通常需要进行分解，将其分解成两个问题进行求解。然而，这种方法实现起来比较困难，并且每次求解只能得到一个可行解，而且不能保证得到最优解。为了解决这个问题，本文提出了一种具有大规模邻域搜索操作的两阶段遗传算法(Two-Stage Genetic Algorithm with Large Neighborhood Search, TSGA-LNS)，适用于求解两阶段模型。该算法分为两个阶段。在第一阶段，本文使用遗传算法对库存阶段进行求解，并将每次迭代得到的解传递到第二阶段的遗传算法中。在第二阶段，本文继续使用遗传算法对路径优化阶段进行求解，并通过加入局部搜索操作来提高解的质量。在局部搜索操作中，本文采用了大规模邻域搜索中的“破坏”和“修复”的思想，具体来说，首先随机选择一些节点，并将它们从现有的路径中删除，以打破原有的路径结构，然后使用修复策略，以最远插入启发式方法重新构建这些节点之间的连接，以形成新的路径。通过这种方式，算法可以拓展解的空间，并提高解的质量。算法步骤如下：

Step 1. 对数据进行处理，使其能够适用于算法。

Step 2. 对库存策略使用第一阶段遗传算法，初始化算法的各变量和参数，确定第一阶段算法的种群大小，迭代次数iterations设为1，最大迭代次数设为n-iterations，周期数设为n-periods，生成初始种群population。

Step 3. 对种群population中的染色体进行选择、交叉、变异操作，得到新种群population1。并根据适应度函数得到适应度值。

Step 4. 将种群population1的解作为第二阶段遗传算法的参数，进而对路径优化策略使用第二阶段遗传算法，初始化算法的各变量和参数，确定第二阶段算法的种群大小，迭代次数gen设为1，最大迭代次数设为MAXGEN，生成初始种群Chrom，将种群Chrom进行解码，得到初始解为配送方案VC，并计算该初始解的适应度值fit。

Step 5. 对种群Chrom中的染色体进行选择、交叉、变异、局部搜索和重组操作，得到新的种群Chrom1。

Step 5.1. 选择操作采用二元竞标赛选择操作，对种群中的两个个体进行比较，选择适应度值更大的那一个个体，放到自带种群当中，最后选择出的子代种群为Selch。

Step 5.2. 将选择出的子代种群Selch进行交叉、变异操作。

Step 5.3. 为了进一步提高种群质量，对交叉变异操作后的子代种群进行局部搜索操作。

Step 5.3.1. 对子代种群Selch中的每一条染色体进行破坏操作，从顾客集合中(染色体)随机选出一个顾客，根据相关性R再次依次移除需要数量to Remove的顾客，最后输出破坏后的配送方案rfvc。

Step 5.3.2. 将破坏后的配送方案rfvc进行修复操作，采用最远插入启发式方法，将最小插入成本最大的顾客插回到当前配送方案，最后得到修复后的配送方案ReIfvc。

Step 6. 在对子代种群Selch进行局部搜索操作后，需要使用重组操作将子代种群与当前迭代起始的父代种群进行合并，最后得到新的种群Chrom。并根据适应度函数得到适应值。

Step 7. 若gen ≤ MAXGEN，则跳至Step 5。

Step 8. 将Step 3和Step 6得到的适应度值进行求和，得到最终适应度值。

Step 9. 若iterations = n-iterations，则输出解集，算法结束。

4. 算例仿真与结果讨论

4.1. 算例数据

为了验证模型和算法的有效性，本文参照文献[21]，并对冷链配送展开行业调研，最终设计了中等规模的算例仿真实验。配送中心和各个零售商的坐标、时间窗、服务时间和期初库存相关信息如表1所示。其余模型参数设置如下：h = 0.1元/kg*天，p = 5元/kg，λ = 0.005，M = 250 kg，ψ = 2.621 kg/L，φ = 200元/辆，σ₁ = 3 L/h，σ₂ = 6 L/h，C_f = 6.5元/L，C_e = 2元/kg，v = 50 km/h，b = 50 kg，η = 40 kw/天，ω_e = 15元/h，ω_l = 15元/h，θ = 0.6 kg/kw，ω₀ = 1900 kg，${\beta }_1=1.004\times {10}^{-3}$ ，${\beta }_2=3.896\times {10}^{-7}$ ，${\beta }^3=1.328\times {10}^{-8}$ 。此外，本文将计划期设置为7天，零售商每天的服务时间窗和服务时间都一致，供应商每天都会对零售商进行库存管理，确定每天的配送量，零售商每天收到产品才开始进行销售；各个零售商的需求是服从均值为50，标准差为10的正态分布。

算法参数设置如下，第一阶段算法的交叉概率$p_c=0.9$ ，第一阶段算法的变异概率$p_m=0.1$ ，第一阶段算法的种群大小${\xi }_1=100$ ，第一阶段算法的最大迭代次数maxgen = 50，第二阶段算法的交叉概率$P_c=0.9$ ，第二阶段算法的变异概率$p_m=0.05$ ，第二阶段算法的种群大小${\xi }_2=50$ ，第二阶段算法的最大迭代次数MAXGEN = 100。

Table 1. Retailers and distribution center-related information

表1. 零售商与配送中心相关信息

编号	坐标	时间窗	服务时间(h)	期初库存(kg)
0	(50, 50)	[0:00~14:00]	-	-
1	(85, 70)	[10:00~10:30]	0.3	25
2	(30, 16)	[9:30~10:00]	0.3	25
3	(18, 65)	[1:00~2:30]	0.3	25
4	(64, 38)	[9:00~9:30]	0.3	25
5	(14, 26)	[0:20~1:00]	0.3	25
6	(37, 82)	[8:00~8:40]	0.3	25
7	(70, 90)	[2:40~3:20]	0.3	25
8	(32, 44)	[3:40~4:30]	0.3	25
9	(77, 24)	[7:40~8:20]	0.3	25
10	(53, 22)	[5:00~5:30]	0.3	25
11	(22, 92)	[6:30~7:20]	0.3	25
12	(8, 46)	[8:10~9:10]	0.3	25
13	(58, 75)	[1:30~2:30]	0.3	25
14	(72, 58)	[7:10~8:10]	0.3	25
15	(90, 34)	[5:40~6:40]	0.3	25

注：编号0为配送中心，编号1~15为零售商。

4.2. LCIRP模型结果分析

本文算法用Matlab编程，在CPU为1.60 GHz，内存4 GB的计算机上运行，程序运行软件为MatlabR2019a。由于本文所考虑的需求服从正态分布的多周期随机需求，因此每次运行得到的需求量都会有所差异。为了获得更优的解，本文采取了以下实验步骤：首先运行10次，从中选择最优的结果；然后根据这次最优结果的需求量再运行10次，并对这10次结果进行比较，以获得更优的求解结果。其求解结果如下所示，各周期的配送路线方案如表2所示，各零售商每周期的配送量如表3所示。

Table 2. The distribution routes for each cycle

表2. 各周期的配送路线

周期	配送路线
1	0-3-13-7-11-6-1-0, 0-5-8-10-15-9-4-14-12-2-0
2	0-13-7-15-14-9-10-6-4-1-0, 0-5-3-12-2-0
3	0-13-7-1-0, 0-5-3-10-15-14-9-11-6-12-2-8-4-0
4	0-5-3-11-12-2-0, 0-13-10-14-15-9-4-1-7-6-0
5	0-5-3-11-6-4-0, 0-13-7-8-12-2-10-15-9-14-1-0
6	0-5-13-7-15-12-2-14-9-4-3-11-1-0
7	0-7-15-9-1-10-4-14-0, 0-5-3-11-6-8-12-0

Table 3. The distribution volume for each retailer per cycle (Unit: kg)

表3. 各零售商每周期的配送量(单位：kg)

零售商	周期1	周期2	周期3	周期4	周期5	周期6	周期7
1	39	60	42	34	38	59	77
2	39	48	64	48	74	35	0
3	37	45	69	30	32	37	42
4	44	68	61	32	35	79	76
5	65	44	76	56	46	79	48
6	52	34	42	61	36	0	31
7	34	36	60	35	61	55	44
8	53	0	68	0	38	0	30
9	47	66	74	53	64	79	30
10	36	62	35	35	69	0	31
11	36	0	56	32	52	40	56
12	41	80	65	76	69	31	51
13	31	66	35	75	46	35	0
14	59	45	49	65	34	68	33
15	40	55	67	55	44	56	73

LCIRP模型的系统总成本为13329.2元，其中库存成本为3712.3元，配送成本为6037.6元，储存过程中的碳排放成本为887.3元，配送过程中的碳排放成本为2692元。

4.3. 与不考虑碳排放约束的库存–路径模型的对比分析

本文建立的LCIRP模型考虑了碳排放约束，并将碳排放转化为成本，加入到成本最小化的目标函数中。现删除目标函数中的碳排放成本函数，将其得到的结果与本文模型得到的结果进行对比。控制其他参数不变，只移除碳排放成本函数，且算例数据和实验步骤不变，对其进行求解，两模型的求解结果如表4所示。

Table 4. Whether to consider the cost of two inventory-routing models with carbon emissions (Unit: yuan)

表4. 是否考虑碳排放的两种库存–路径模型成本(单位：元)

	总成本	库存成本	配送成本	总碳排放成本
考虑碳排放	13329.2	3712.3	6037.6	3579.3
不考虑碳排放	15375.8	4585.8	6616.3	4173.7

由上表可以看出，LCIRP模型的总成本以及各项成本都比不考虑碳排放成本的库存–路径模型要低。通过分析计算，相对于不考虑碳排放的模型的总成本，考虑碳排放的模型的总成本降低了13.3%，碳排放量降低了14.2%。由此可知，对于企业而言，不考虑碳排放时，系统总成本有所增加，不利于提升经济效益；对于社会而言，碳排放量有所增加，不利于提升环境效益。因此，政府应该鼓励企业去考虑碳排放，不仅有利于环境保护，同时也可以帮助企业控制成本。

4.4. 库存–路径两阶段优化模型与单阶段优化模型的对比分析

为比较库存–路径两阶段优化模型与单阶段优化模型的两种结果，在算例数据和实验步骤不变的前提下，本文用GA-LNS求解单阶段优化模型，单阶段优化模型是将库存路径分开考虑，两个环节相互独立，供应商将零售商的需求量作为配送量，不进行库存管理。将库存–路径两阶段优化方案称为方案1，单阶段优化方案称为方案2，方案2各周期的配送路线方案如表5所示，各零售商每周期的配送量如表6所示。

Table 5. The distribution routes for each cycle in scheme 2

表5. 方案2各周期的配送路线

周期	配送路线
1	0-10-15-9-4-14-1-0, 0-13-7-11-12-2-0, 0-5-3-8-6-0
2	0-5-8-10-14-6-0, 0-3-11-12-2-0, 0-13-7-15-9-4-1-0
3	0-8-10-15-9-4-0, 0-13-7-11-6-14-1-0, 0-5-3-12-2-0
4	0-13-7-11-12-2-0, 0-8-10-14-15-9-4-1-0, 0-5-3-6-0
5	0-5-8-12-2-0, 0-3-7-11-6-1-0, 0-13-10-15-9-14-4-0
6	0-5-3-8-12-1-0, 0-10-15-14-9-4-0, 0-13-7-11-6-2-0
7	0-8-10-15-9-1-0, 0-13-7-14-12-2-0, 0-5-3-11-6-4-0

Table 6. The distribution volume for each retailer per cycle in scheme 2 (Unit: kg)

表6. 方案2各零售商每周期的配送量(单位：kg)

零售商	周期1	周期2	周期3	周期4	周期5	周期6	周期7
1	55	68	27	59	53	37	46
2	86	78	37	80	57	49	57
3	49	65	64	64	57	38	57
4	55	60	57	47	53	42	59
5	39	42	21	64	53	42	64
6	49	48	53	53	41	50	48
7	61	61	41	51	38	39	50
8	42	54	48	61	39	50	56
9	65	51	35	43	39	74	44
10	48	59	42	36	36	55	48
11	64	53	52	66	42	57	58
12	52	38	39	51	57	76	43
13	49	31	46	32	58	41	51
14	53	44	55	57	67	48	29
15	64	39	60	51	64	30	48

由于方案2是单独考虑库存和配送路径，库存和配送路径是两个相互独立的环节，供应商不对零售商进行库存管理，而零售商的需求量则为配送量。考虑到冷链产品的易变质性质，根据公式(5)来分析计算方案2的库存环节相关成本及碳排放成本。最终得到的总成本为14542.3元，其中库存环节的总成本为3606.6元，配送环节的总成本为6994.7元，碳排放总成本为3941.0元。将方案1和方案2的结果进行对比，其结果对比如表7所示。

Table 7. Comparison of results between different schemes

表7. 不同方案的结果对比

	总成本	库存成本	配送成本	库存碳排放成本	配送碳排放成本	总碳排放成本
方案1	13329.2	3712.3	6037.6	887.3	2692	3579.3
方案2	14542.3	3606.6	6994.7	761.7	3179.3	3941.0

分析上表，可以得出以下结论：

1) 相比于单阶段优化模型，库存–路径两阶段优化模型可以更好地节省系统成本，通过对比发现，系统总成本降低了8.3%。对于冷链物流企业而言，对库存和配送进行联合优化可以帮助企业降低系统成本，是一种更为经济的管理模式。

2) 通过对比两种方案的库存环节成本和配送环节成本，可以看出方案1和方案2也有所差异，由于方案2仅根据零售商需求量进行配送，没有对其库存进行管理，而方案1通过对零售商库存进行管理，同时考虑了零售商的期初期末库存和缺货成本，所以方案2零售商的平均库存量会较小，并且不会有缺货成本，从而库存环节的成本会更低。对于配送环节的成本，由于方案1零售商会进行库存管理，零售商并不是每周期都会进行配送，从而配送环节的成本会更低。通过对比发现，两个环节成本之和还是降低了8.0%，所以即使选择库存路径联合优化方案的库存环节成本会更高，但从两个环节的总体来看，选择库存路径联合优化方案还是更加经济。

3) 对于碳排放而言，方案1由于会对零售商进行库存管理，零售商不是每周期都会进行配送，通过油耗量计算得到的碳排放量则会更低，并且冷藏车服务次数会降低，服务时的碳排放量也会更低。通过计算可得，相比于方案2，方案1的碳排放量降低了9.2%，所以选择库存路径联合优化方案可以更好地降低碳排放量，从而更好地达到环保效果。

通过以上分析，无论从经济方面还是环境方面，库存–路径两阶段优化模型都是更佳选择。对于企业而言，可以将物流活动中的库存和配送两个环节联合起来考虑，帮助企业更好地提升经济效益；对于政府而言，应该提倡库存路径联合优化模式，降低企业碳排放量，从而促进社会环境的发展。

4.5. 基于碳交易价格的灵敏度分析

Table 8. Carbon price change table

表8. 碳价格变化表

碳价格(元/kg)	总成本(元)	碳排放成本(元)	碳排放量(kg)
0.5	10868.2	932.6	1865.2
1	11990.3	1855.7	1855.7
1.5	12526.2	2717.9	1811.9
2	13329.2	3579.3	1789.7
2.5	14204.0	4450.5	1780.2
3	14657.7	5319.6	1773.2
3.5	14843.5	6190.8	1768.8

由于我国当前碳交易价格不稳定，容易产生浮动。为了了解碳交易价格对企业和社会的影响，本文对碳交易价格进行灵敏度分析，结果如表8和图1所示。可以看出，当碳价格在[0.5, 3.5]之间时，碳排放量随着碳交易价格的增长而降低，从1865.2 kg降低到1768.8 kg，降低了5.2%，此时企业的系统总成本不断增长。此外，随着碳交易价格的增长，企业的碳排放量的降低趋势也逐渐变缓。

由图1可看出，当碳价格偏低时，企业的碳排放量将失去有效约束，从而无法实现环保目标；反之，若碳价格过高，企业的运营成本也会相应增加，不利于经济发展。因此，碳价格过高或过低都无法有效约束企业的碳排放行为。值得注意的是，当碳交易价格处于[1, 2]区间时，企业将更加敏感地关注碳价格，并主动采取措施降低碳排放。然而，当碳交易价格超过2元/kg时，企业碳排放量将逐渐失去对碳交易价格的敏感度。因此，政府在制定碳交易政策时，需要综合考虑碳交易价格对企业的影响，以合理控制碳交易价格，平衡经济成本和环保效益。如何使企业自愿参与碳减排活动，是当前推动绿色环保运动亟待解决的问题。

Figure 1. Graph of the relationship between total cost and carbon emissions under different carbon trading prices

图1. 不同碳交易价格下的总成本和碳排放量关系图

5. 结束语

本文以供应商和零售商所组成的二级冷链物流网络为研究对象，建立了LCIRP模型，并设计了一种两阶段遗传算法，在算法中采用了大规模邻域搜索(LNS)中的“破坏”和“修复”的思想，通过对算例仿真结果进行分析，并得出以下研究结论：1) 考虑碳排放的库存–路径模型在降低系统总成本、约束碳排放量方面有较好的效果。2) 相比于单阶段优化模型，LCIRP模型能同时兼顾企业经济效益和社会环境效益，能更好地实现双赢局面。3) 碳交易价格对冷链物流的系统总成本和碳排放量有着很大的影响，过高过低都不能对碳排放行为有较好的约束。4) 为了确保社会环境的可持续发展，政府和企业必须携手合作，政府应积极倡导实施库存路径联合优化策略，并完善碳价格的制定，以促使企业自觉地采取措施降低碳排放。

基金项目

国家自然科学基金(71871143)；上海“科技创新行动计划”社会发展科技攻关项目(22dz1203405)；教育部人文社科基金项目(22YJC790189)；教育部人文社科基金项目(23YJC790046)。