交通 | 应用Benders分解方法解决多车生产路由问题

论文解读​ 曲晨辉，王飞龙

1 知识补充和文章贡献

2.1 IRP (Inventory routing problems)

IRP关注的是在一个给定的规划范围内，从一个设施到一组客户的单一产品分配。客户以给定的速度消费产品，并可将产品的库存维持在一个特定的水平上。一组相同的车辆可用于产品的配送。目标是最小化总配送成本，计算为路线成本和库存持有成本的总和，同时不造成任何客户的缺货。

库存路由问题 (简称IRP)作为物流配送系统的基础理论问题,主要研究的是为供应商制定给若干客户配送产品的计划,需要同时决策客户配送方案和车辆配送路线,在满足多项约束的情况下最终实现总开销最低的目标。

在IRP中，库存控制和路由决策必须同时做出。

2.2 PRP (production routing problem)

生产路线问题 (PRP)同时考虑了供应链中的各种决策，包括生产、库存管理和路由决策。 一般来说，在一个多时期的范围内，这些综合供应链问题要求人们确定生产和访问每个客户的时期，相应的生产和交付给客户的数量，以及车辆的详细路线计划。

2.3 IRP中的三种补货策略

在IRP问题中，主要有以下三种策略为网络中的客户配送货物：

•第一种策略称为“Order-up-to-level” (OU) 策略，车辆在每次访问某个客户时，为客户配送的货量均使得客户达到最大的库存水平
•第二种策略称为“maximum-level” (ML)策略，这种策略是对OU策略的一种放松，车辆到达客户时配送的货量可以介于客户当前货量和最大库存水平之间
•第三种策略称为“replenishment” (RP)策略，这种策略对ML策略进一步放松，允许车辆为客户配送任意非负货量

2.4 文章贡献

本文研究了多车PRP (MVPRP)，并使用OU策略为客户进行货物配送；为了精确求解MVPRP，文章基于Benders分解，提出了一种精确算法框架，将模型划分为一个主问题(master problem)和一个子问题(slave problem)，其中子问题又可以进一步划分为多个CVRP问题，并构建基于路径的模型进行求解；文章同时对MVPRP最优解的下界进行了分析，从而生成初始最优割约束，用于加快模型收敛速度。最终，借助MVPRP和VMIRP的算例验证了精确算法的有效性。

2 问题定义

为工厂附近长度为$l$的离散规划周期的所有顾客提供单一货物。在每个周期开始的时候，生产工厂可以生产最多$C$单位的产品，其中固定支出是$f$，单位产品支出是$u$。节点$i,i\in N$有LiL_iLi​单位的仓库容量以及一个库存单位支出hih_ihi​。节点$i$最开始在规划区域持有的产品数量为Ii0I_{i0}Ii0​。顾客i∈Nci\\in N_ci∈Nc​在周期$t$需要ditd_{it}dit​单位的产品，并且可以在每个周期最多被访问一次。容量为$Q$的$m$辆相同车辆在仓库可用。

其中，MVPRP包括同时决定：

•工厂何时生产和生产多少产品
•何时向每个客户交货以及交货数量
•在规划范围内的每个周期使用什么路线

目标是在规划期内使生产、库存和路由成本的总和最小化，而不造成任何客户的缺货。

在本文中，以OU政策为补充政策。此外，正如文献中通常的假设，每个周期在工厂的生产在交付前进行，在客户的交付在该时间段的开始时执行。

3 数学模型

3.1 符号定义

$T$: 时间周期$1,...,l$集合。为简化模型，引入一个虚拟的周期$l+1$，并定义T′=T∪{l+1}T'=T\\cup \\{l+1\\}T′=T∪{l+1};givtg_{ivt}givt​: 当顾客$i$在周期$t$被访问的总运输数量，并且当前访问是在周期$v$，计算如下：

eivte_{ivt}eivt​:当顾客$i$在周期$t$被访问的总库存持有成本，并且当前访问是在周期$v$，计算如下：

$\mu (i,t)$: 当顾客$i$下一次可以在不发生缺货的情况下补货的，在周期$t$后的最近的周期，计算公式为：

$\pi (i,t)$: 当顾客$i$下一次可以在不发生缺货的情况下补货的，在周期$t$前的最早的周期，计算公式为：

ptp_tpt​: 非负连续变量，表示在周期$t$的生产数量

IitI_{it}Iit​: 非负连续变量，代表在周期$t$末尾节点$i$的库存水平

yty_tyt​: 0-1变量，表示在周期$t$工厂是否进行生产

zitz_{it}zit​: 0-1变量，表示节点i,i∈NCi,i\\in N_Ci,i∈NC​是否在周期$t$被访问

sts_tst​: 非负整数变量，表示周期$t$使用的车辆数量

xijtx_{ijt}xijt​: 整数变量，可能取值在{0,1},{i,j}∈E,t∈T\\{0,1\\},\\{i,j\\}\\in E,t\\in T{0,1},{i,j}∈E,t∈T，和{0,1,2},{0,j}∈E,t∈T\\{0,1,2\\},\\{0,j\\}\\in E,t\\in T{0,1,2},{0,j}∈E,t∈T

λivt\\lambda_{ivt}λivt​: 0-1变量，表示节点$i$是否在周期$t$被访问，当前访问是周期$v$

$E(S)$: 两个点都在集合$S$中的边集，即E(S)={{i,j}∈E:i,j∈S,S⊆N}E(S)=\\{\\{i,j\\}\\in E: i,j\\in S, S\\subseteq N\\}E(S)={{i,j}∈E:i,j∈S,S⊆N}

$\delta (S)$: 与集合$S$中的一个节点相关联的边集，即δ(S)={{i,j}∈E:i∈S,j∉Sori∉S,j∈S}\\delta(S)=\\{\\{i,j\\}\\in E: i\\in S, j\\notin S~or~i\\notin S,j\\in S\\}δ(S)={{i,j}∈E:i∈S,j∈/S or i∈/S,j∈S}，$S$也可以是单独的一个节点$i$

补充解释：在比如λivt\\lambda_{ivt}λivt​中，v和t都是访问i的周期。

4 Benders分解​

4.1 基于逻辑的Benders模型重构

BMP（工厂生产、计划配送周期、配送产品的量）

BSP（进行运输，相当于$l$个周期的CVRP）

于是，给定一个解(p∗,I∗,y∗,z∗,λ∗,ω∗)(p^*,I^*,y^*,z^*,\\lambda^*,\\omega^*)(p∗,I∗,y∗,z∗,λ∗,ω∗)
有三种可能的情况：

•给定(z∗,λ∗)(z^*,\\lambda^*)(z∗,λ∗)BSP不可行，BMP中加入infeasibility cut：

•给定(z∗,λ∗)(z^*,\\lambda^*)(z∗,λ∗)BSP可行，但是ϕ(z∗,λ∗)>∑t∈Tωt∗\\phi(z^*,\\lambda^*)>\\sum_{t\\in T}\\omega_t^*ϕ(z∗,λ∗)>∑t∈T​ωt∗​ ，那么就需要加入optimality cut，在BMP中：

•其中LBRLB_RLBR​是$\varphi (z,\lambda )$的下界值。

•给定(z∗,λ∗)(z^*,\\lambda^*)(z∗,λ∗)BSP可行，但是ϕ(z∗,λ∗)=∑t∈Tωt∗=∑t∈T∑{i,j}∈Ecijxijt∗\\phi(z^*,\\lambda^*)=\\sum_{t\\in T}\\omega_t^*=\\sum_{t\\in T}\\sum_{\\{i,j\\}\\in E}c_{ij}x_{ijt}^*ϕ(z∗,λ∗)=∑t∈T​ωt∗​=∑t∈T​∑{i,j}∈E​cij​xijt∗​，(s∗,x∗)(s^*,x^*)(s∗,x∗)是BSP的最优解。那么MVPRP问题的最优值就是：

基于上述模型的精确算法如下：

•初始化。计算LBRLB_RLBR​，并令$LB=0,UB=\infty$

•当$LB$

1. 解BMP。如果BMP不可行，说明原问题不可行，就结束计算。否则，令(p∗,I∗,y∗,z∗,λ∗,ω∗)(p^*,I^*,y^*,z^*,\\lambda^*,\\omega^*)(p∗,I∗,y∗,z∗,λ∗,ω∗)为原问题最优解的成本z‾\\overline{z}z，并令LB=z‾LB=\\overline{z}LB=z

2. 解BSP。可能会出现两种情况：

   a). BSP不可行。加入infeasibility cut。

   b). BSP可行且是整数解，满足ϕ(z∗,λ∗)>∑t∈Tωt∗\\phi(z^*,\\lambda^*)>\\sum_{t\\in T}\\omega_t^*ϕ(z∗,λ∗)>∑t∈T​ωt∗​。加入optimality cut。

3. 令(p∗,I∗,y∗,z∗,λ∗,ω∗)(p^*,I^*,y^*,z^*,\\lambda^*,\\omega^*)(p∗,I∗,y∗,z∗,λ∗,ω∗)为原问题可行解的成本z‾\\overline{z}z，并令UB=min{UB,z‾}UB=min\\{UB,\\overline{z}\\}UB=min{UB,z}。

上述算法每一次迭代中，LB都代表着原问题最优值的有效下界，最后会和UB重合，即收敛。

4.2 解BSP

给定BMP问题的解(p‾,I‾,y‾,z‾,λ‾,ω‾)(\\overline{p},\\overline{I},\\overline{y},\\overline{z},\\overline{\\lambda},\\overline{\\omega})(p​,I,y​,z,λ,ω)。

BSP（进行运输，相当于l个周期的CVRP）

可以看做是CVRP(t)，可以用一个无向图Gt=(Vt,Et)G^t=(V^t,E^t)Gt=(Vt,Et)表示，其中Vt={0}∪NctV^t=\\{0\\}\\cup N_c^tVt={0}∪Nct​（就是t周期的所有需要被服务的顾客节点以及工厂节点Nct={i∈Nc:z‾it=1}N_c^t=\\{i\\in N_c:\\overline{z}_{it}=1\\}Nct​={i∈Nc​:zit​=1}，Et={{i,j}:i,j∈VtE^t=\\{\\{i,j\\}:i,j\\in V^tEt={{i,j}:i,j∈Vt，$i$就是边集。边{i,j}\\{i,j\\}{i,j}都有一个成本cijt=cijc_{ij}^t=c_{ij}cijt​=cij​，每个顾客都有一个已知非负的需求qit=giv‾itq_{it}=g_{i{\\overline{v}_i}t}qit​=givi​t​ ，其中v‾i\\overline{v}_ivi​是当λ‾iv‾it=1\\overline{\\lambda}_{i{\\overline{v}_i}t}=1λivi​t​=1时确认的。有m量相同的车，每一辆都有一个容量$Q$，在工厂可用。

该问题就包括了找到一组包含最多m条路的花费最少的路径。路由包括如下规则：
•每条路径都访问工厂
•每个顾客节点只被一条路径访问
•同一条路径被访问的顶点的需求和不能超过车辆容量Q

4.2.1 解CVRP(t)

为解决每个CVRP(t)，使用了基于集合划分模型的精确算法。

为简化计算，省略了索引$t$。令$\mathcal{R}$为所有可行路径的索引集合，aira_{ir}air​为二进制系数，表示顶点i∈Ncti\\in N_c^ti∈Nct​是否属于路径$r\in \mathcal{R}$（对于任意路径来说，a0r=1a_{0r}=1a0r​=1）。每一条路径$r\in \mathcal{R}$都有一个成本brb_rbr​，相当于路径r的TSP问题的最优值。令ξr\\xi_rξr​为二进制变量，代表当且仅当$r\in \mathcal{R}$属于最优解的时候其为1。模型如下：

CVRP(t)可能会infeasible，因为车辆可能不够。因此就要加入infeasibility cut。为检查其是否有可行解，可以通过解决装箱问题。有∣Nct∣|N_c^t|∣Nct​∣个物品，重量{qit}\\{q_{it}\\}{qit​}，最多有$m$个箱子，容量为$Q$。如果BPP是不可行的，那么就加入如下infeasibility cut：

该cut可以进一步加强。如果有这样一个解λivt≥λ‾ivt,∀i∈Nc,v=π(i,t),...,t−1\\lambda_{ivt}\\geq\\overline{\\lambda}_{ivt},\\forall i\\in N_c, v=\\pi(i,t),...,t-1λivt​≥λivt​,∀i∈Nc​,v=π(i,t),...,t−1，它也是不可行的。所以如下的infeasibility cut可以被加入BMP：

此时还有一种情况，对于一个给定的$i\\in N_c $，可以得到$g_{i{v_1}t}>g_{i{v_2}t},if~v_1$，也就是说如果顾客i当前被访问的时间是$[\\pi(i,t),\\overline{v}_i-1]$，那么在t周期内需要往i点运送的货物需求量就要大于其位于$\\overline{v}_i $的时候的量。因此如下被加强的infeasibility cut就可以得到：

如果BPP实例表明有可行解，那么就说明CVRP(t)有有限的最优解。于是可以使用基于路径枚举程序（2008, Baldacci, Christofides, and Mingozzi）和混合策略（2009, Pessoa, Uchoa, Poggi de Aragao, 2008, Pessoa, de Aragao, Uchoa）来解CVRP(t)。

4.2.2 更新BMP

CVRP(t)问题首先被检查是否可行，并加入infeasibility cut进BMP。如果至少有一个infeasibility cut被找到，程序终止，BMP重新被计算。否则，如果所有CVRP(t)都可行，那么问题就根据周期$(t=1,2,...,l)$$用精确算法解CVRP(t)，并加入optimality cut到BMP。如果有找到cut，程序终止，BMP重新计算。

如果没有cut被找到，那么CVRP(t)就都最优了。然后调用MVPRP的精确算法的step2.b.ii再检查一下是否有需要加入BMP的optimality cut。

4.3 计算路由成本的下界LBRLB_RLBR​

令fit,i∈Nc,t∈Tf_{it},i\\in N_c,t\\in Tfit​,i∈Nc​,t∈T为顾客$i$必须在周期$t$之前被访问的数目的下界，fit=⌈q‾it/min{Q,Li}⌉f_{it}=\\lceil{\\overline{q}_{it}/min\\{Q,L_i\\}}\\rceilfit​=⌈q​it​/min{Q,Li​}⌉。对于每个i∈Nc,t∈Ti\\in N_c,t\\in Ti∈Nc​,t∈T，令q‾it\\overline{q}_{it}q​it​为顾客$i$在周期$t$之前的累积需求量，计算q‾it=max{0,−Ii0+∑v=1tdiv},t=1,...,l.\\overline{q}_{it}=max\\{0,-I_{i0}+\\sum_{v=1}^td_{iv}\\},t=1,...,l.q​it​=max{0,−Ii0​+∑v=1t​div​},t=1,...,l.。另外，mLm_LmL​和mUm_UmU​是在计划区域内车辆总数的下界和上界。令q^i=∑t∈Tdit−Ii0\\hat{q}_i=\\sum_{t\\in T}d_{it}-I_{i0}q^​i​=∑t∈T​dit​−Ii0​是在计划区域内顾客i∈Nc的总产品需求量。i\\in N_c的总产品需求量。i∈Nc​的总产品需求量。mL=⌈1Q∑i∈Ncq^i⌉,mU=m⋅lm_L=\\lceil{1\\over Q}\\sum_{i\\in N_c}\\hat{q}_i\\rceil, m_U=m\\cdot lmL​=⌈Q1​∑i∈Nc​​q^​i​⌉,mU​=m⋅l$。$qi,∀i∈Ncq_i,\\forall i\\in N_cqi​,∀i∈Nc​$是在计划区域内的任意一次访问中传递给顾客$i$的产品数量的下界，$qi=mint∈T{dit}q_i=min_{t\\in T}\\{d_{it}\\}qi​=mint∈T​{dit​}

定义一条路径为图G中经过工厂0号节点的最小成本simple cycle，使用需求{qi}\\{q_{i}\\}{qi​}计算的顾客访问总需求不会超过车辆容量$Q$。令R~\\tilde{\\mathfrak{R}}R~为所有路径的索引集合，令aira_{ir}air​为一个二进制系数，表示顶点i∈Nci\\in N_ci∈Nc​是否属于路径r∈R~r\\in \\tilde{\\mathfrak{R}}r∈R~，每条路径都有一个路径成本brb_rbr​。

任意一个最优的MVPRP问题解的路径成本下界值都可以根据如下的整数规划问题的最优成本解得到：

下界LBRLB_RLBR​是根据一个bounding程序计算的，作为RF问题的线性规划松弛的一个近似最优对偶解。该程序和基于单纯形法的标准列生成问题不同，因为它用到了dual ascent启发式方法来解主问题。

5 优化BMP

这里主要是介绍如何加强BMP。

Aduyasak, Cordeau, and Jans在14年提出了很多加强F模型线性松弛界的方法，这里介绍了一个：

其中t′t't′指的是工厂为了避免缺货必须生产的最早时间，计算如下：

我们还可以在BMP中加入如下不等式：

5.1 optimality cut的初始集合

介绍两种生成初始optimality cut的方式。

令R^\\hat{\\mathfrak{R}}R^为周期$t$所有可行路径的索引集合，其中顾客$i$的需求qitq_{it}qit​属于需求的离散集合{givt:v=π(i,t),...,t−1}\\{g_{ivt}:v=\\pi(i,t),...,t-1\\}{givt​:v=π(i,t),...,t−1}，也就是说，t周期内被访问的顾客$i$的所有可能需求集合根据其潜在曾访问时期的集合{π(i,t),...,t−1}\\{\\pi(i,t),...,t-1\\}{π(i,t),...,t−1}确定。令aivra_{ivr}aivr​为二进制系数，表示当且仅当顾客$i$属于路径$r$，并且有一个等于givtg_{ivt}givt​的需求。CVRP(t)建模如下：

令wiv∈R,∀i∈Nc,v=π(i,t),...,t−1,u0≤0w_{iv}\\in R, \\forall i\\in N_c,v=\\pi(i,t),...,t-1,u_0\\leq0wiv​∈R,∀i∈Nc​,v=π(i,t),...,t−1,u0​≤0分别为约束（39）和（40）的对偶变量。CVRP(t)的线性松弛对偶模型如下所示：

5.1.1 Type 1 Cut

第一种cut是基于CVRP(t)得到的一个LP问题

令wiv≥0,∀i∈Nc,v=π(i,t),...,t−1,u0≤0w_{iv}\\ge0, \\forall i\\in N_c,v=\\pi(i,t),...,t-1,u_0\\leq0wiv​≥0,∀i∈Nc​,v=π(i,t),...,t−1,u0​≤0分别为约束（41）和（42）的对偶变量。令(w∗,u0∗)(w^*,u_0^*)(w∗,u0∗​)为（T1)的对偶问题的可行解，基于proposition1，可以得到如下有效的optimality cut：

5.1.2 Type 2 Cut

第二个cut是基于从F模型中得到的LP问题。令Rt^\\hat{\\mathfrak{R}_t}Rt​^​为周期t所有可行路径的索引集合，其中顾客$i$的需求qitq_{it}qit​属于需求的离散集合{givt:v=π(i,t),...,t−1}\\{g_{ivt}:v=\\pi(i,t),...,t-1\\}{givt​:v=π(i,t),...,t−1}，也就是说，t周期内被访问的顾客$i$的所有可能需求集合根据其潜在曾访问时期的集合{π(i,t),...,t−1}\\{\\pi(i,t),...,t-1\\}{π(i,t),...,t−1}确定。令aivrra_{ivrr}aivrr​为二进制系数，表示当且仅当顾客$i$在周期t属于路径$r$，并且有一个等于givtg_{ivt}givt​的需求。使用二进制变量ξrt\\xi_{rt}ξrt​的模型如下所示，ξrt\\xi_{rt}ξrt​表示周期t时路径r是否在解中。

令wivt∈R,∀i∈Nc,v=π(i,t),...,t−1,u0≤0w_{ivt}\\in R, \\forall i\\in N_c,v=\\pi(i,t),...,t-1,u_0\\leq0wivt​∈R,∀i∈Nc​,v=π(i,t),...,t−1,u0​≤0分别为约束（44）和（45）的对偶变量。关于ξrt\\xi_{rt}ξrt​的对偶约束为：

令$(w^*,u_0*) $为上述模型对偶可行解的变量。基于proposition1，如下不等式被加入：

6. 计算结果分析

文章将提出的算法(EXM)和Adulyasak, Cordeau, and Jans (2014)提出的三种Branch-and-cut算法在针对MVPRP算例和VMIPR算例的求解效果上进行了比较，Adulyasak, Cordeau, and Jans (2014)提出的三种Branch-and-cut算法如下所示：

• Veh-Ind：在单线程上运行无车辆索引的模型
• Non-Veh-Ind：在8个线程上运行无车辆索引的模型
• 8c-Veh-Ind：在8个线程上运行带车辆索引的模型

文章限定EXM算法最长求解时间为2h，详细实验结果数据可见原文，实验结果总结如下：

• EXM相对于Adulyasak, Cordeau, and Jans (2014)中最优的算法8c-Veh-Ind，能多求解21个算例到最优
• EXM最大求解的算例规模为包含三个时间周期的40个客户的算例，相对于8c-Veh-Ind的最大求解规模多了15个客户
• 对于没有求解到最优的算例，EXM收敛到的下界解和上界解之间的gap值相对较小

文章同样对EXM算法中不同组成成分的有效性进行了验证：

• 使用type I和type I割，使用强化版本的不可行和线性割来初始化EXM对于EXM算法的求解性能影响很大
• 如果不适用线性割了，则EXM算法调用CPLEX IP求解器的次数将会极大增加，从而减慢算法效率

对于VMIPR算例的实验结果总结如下：

• 从求解算例的数量角度，EXM相对于Non-Veh-Ind算法求解的算例多，但相对于Veh-Ind和8c-Veh-Ind算法求解的算例要少一些
• EXM相对于Non-Veh-Ind可以求解额外的四个算例到最优解，相对于8c-Veh-Ind，可以获得8个算例的更优的上界
• EXM对于规划周期数量增长的适应度要更强，对于30个客户，6个周期的算例求解到了新的上界和下界。

7 总结和未来研究方向

文章针对MVPRP (使用OU补货策略)和VMIPR基于Benders分解提出了一种精确算法。通过算例实验证明了文章提出的精确算法相对于SOTA的算法在求解MVPRP的效果上要更优，最多求解3个周期下，使用3辆车辆，40个客户的算例；同时，文章提出的算法对于VMIPR算例，最多可以求解6个周期下，使用4辆车辆，30个客户的算例。

文章提出的算法基于分解的思想，可以较好地处理路径内的约束，后续研究可以通过设计更加复杂的模型和算法，在PRP中考虑时间窗和距离限制约束。

参考文献

Zhang, Z., Luo, Z., Baldacci, R., & Lim, A. (2021). A Benders decomposition approach for the multivehicle production routing problem with order-up-to-level policy. Transportation Science, 55(1), 160-178.