快递的同城运输是指一个城市内分拨中心与大量网点间的货物运输。目前很多快递公司的同城运输都比较混乱，缺乏管理，参与运输的车辆数很多，装载率低，导致运输成本较高。交货车辆数多，还会造成分拨中心的拥堵，从而影响快递的时效。如何从全局出发，统一调度车辆，将各网点的货物运到分拨中心，是降低运输成本的关键性问题。

快递的同城运输问题本质上是一个车辆路径问题(vehicle routing problem, VRP)，最早由Dantzig等^[1]在1959年首次提出。几十年来，经过国外学者们的大量研究^[2-3]，VRP已扩展出多种相关问题。21世纪以来，国内学者也逐渐开始研究VRP。张智海等^[4]将并行遗传算法应用于带时间窗的车辆路径问题研究；张文博等^[5]将遗传算法和模拟退火算法相结合，求解了带动态时间窗的车辆路径问题；Luo等^[6]采用了改进的混合蛙跳算法来求解多车场车辆路径问题；王晓博等^[7]研究了多车型混合装卸车辆路径问题，并设计混合遗传算法求解。

与一般的VRP相比，城市VRP需要考虑更多时间、环境因素^[8]。Malandraki等^[9]、Donati等^[10]研究了有时间依赖的车辆路径问题，这类问题将两点间的行驶时间视为随时间动态变化的。Gendreau等^[11]、Fleischmann等^[12]研究了实时动态调度的车辆路径问题。Chen等^[13]综合考虑了时间依赖性以及动态调度。Mamasis等^[14]、Minis等^[15]考虑了车辆故障的突发情况；李宏^[16]研究了城市冷链物流，考虑了时间窗、交通状况、产品腐坏等因素。

本文以快递同城运输为研究对象，除了容量约束、时间窗约束等常规条件外，还考虑了车辆限行情况以及错峰交货时间。首先建立了以最小化运输成本为目标函数的混合整数规划模型，然后依次提出了启发式算法以及列生成算法，最后通过实例对3个方法进行比较，分析各自的特点并得出结论。

1 问题描述

记网点集合为 $ N = \{ 1,2,\cdots,n\} $ ，分拨中心为0，车型集合为 $\{ 1,2,\cdots,\alpha \} $ ，每种车型的货车数量充足。车辆从任意一个网点装货出发，依次前往若干网点装货，最后到分拨中心卸货。车辆的容量约束只考虑体积。每个网点都有开始装货的时间窗和装货时长，车辆到达某网点的时间不能晚于其时间窗的右区间，但早于左区间到达时，可以先在原地等待。由于网点都希望发货量大一点，所以他们的时间窗普遍比较晚，而且重合度很高，窗口很窄。为了错峰到达分拨中心，避免发生拥堵，每个网点有各自独立的最晚交货时间。此外，某些网点所在区域对车型有限行约束，即禁止某些大型车通行。该优化问题的目标函数为运输成本，包括车辆使用的固定成本以及油耗成本。

在实际中，有一些网点的货量很大，一辆车不足以装下所有货物。此时优先用该区域能通行的最大型车辆装满直接运到分拨中心，直到剩下的货量能用一辆车装下。经过这样的预处理后，可以假设每个网点只能被一辆车服务。令车辆数为 $m = n \times \alpha $ 辆，每种车型的车辆数均为 $n$ ，这样能保证各车型的数量都是足够的。记车辆集合为 $M = \{ 1,2,\cdots,m\} $ ，令第 $i(i \in M)$ 辆车的车型为 $\left\lceil {i/n} \right\rceil $ ， $\left\lceil {} \right\rceil $ 表示向上取整，从而可以按车型给出车辆集合 $M$ 中每辆车的相关参数。

2 数学模型

根据问题描述，给定输入参数，确定决策变量，建立以最小化成本为目标的混合整数规划模型。

2.1 模型参数

${v_i}$ ：网点 $i$ 的货量， $\forall i \in N$ ；

${Q_k}$ ：车辆 $k$ 的最大容积， $\forall k \in M$ ；

$[t_i^1,t_i^2]$ ：网点 $i$ 的开始装货时间窗， $\forall i \in N$ ；

${l_i}$ ：网点 $i$ 的装货时长， $\forall i \in N$ ；

${T_i}$ ：网点 $i$ 到分拨中心的最晚交货时间， $\forall i \in N$ ；

${b_{ik}}$ ：当车辆 $k$ 在网点 $i$ 处被限行时，其值为1，否则为0， $\forall i \in N,k \in M$ ；

${t_{ij}}$ ：点 $i$ 到点 $j$ 的行驶时间， $\forall i,j \in N \cup \{ 0\} $ ；

${d_{ij}}$ ：点 $i$ 到点 $j$ 的行驶距离， $\forall i,j \in N \cup \{ 0\} $ ；

$c_k^1$ ：使用车辆 $k$ 的固定成本， $\forall k \in M$ ；

$c_k^2$ ：车辆 $k$ 的单公里油耗成本， $\forall k \in M$ ；

$R$ ：一个足够大的常数。

2.2 决策变量

${Z_k}$ ：当车辆 $k$ 被使用时，其值为1，否则为0， $\forall k \in M$ ；

${y_{ik}}$ ：当网点 $i$ 由车辆 $k$ 服务时，其值为1，否则为0， $\forall i \in N,k \in M$ ；

${x_{ijk}}$ ：当车辆 $k$ 从点 $i$ 到点 $j$ 时，其值为1，否则为0， $\forall i \in N,j \in N \cup \{ 0\} ,k \in M$ ；

${T_{ik}}$ ：车辆 $k$ 在网点 $i$ 开始装货的时间， $\forall i \in N,$ $k \in M $ ；

${T_{0k}}$ ：车辆 $k$ 到达分拨中心的时间， $\forall k \in M$ 。

2.3 模型建立

$\qquad\min \;\left(\sum\limits_{k = 1}^m {c_k^1{{\textit{z}}_k} + \sum\limits_{k = 1}^m {c_k^2\sum\limits_{i = 1}^n {\sum\limits_{j = 0}^n {{d_{ij}}{x_{ijk}}} } } }\right) {\text{。}}$

(1)

$ \begin{split} \qquad & \;{\rm{s}}{\rm{.t}}{\rm{.}}\\ &\sum\limits_{j = 0}^n {{x_{ijk}}} = {y_{ik}}, \forall i \in N,k \in M; \end{split} $

(2)

$ \qquad \sum\limits_{i = 1}^n {{x_{ijk}}} {\text{≤}} {y_{jk}}, \forall j \in N,k \in M;$

(3)

$ \qquad\sum\limits_{k = 1}^m {{y_{ik}}} = 1, \forall i \in N;$

(4)

$\qquad \sum\limits_{i = 1}^n {{x_{i0k}}} = {Z_k}, \forall k \in M;$

(5)

$ \qquad {y_{ik}}{\text{≤}} 1 - {b_{ik}}, \forall i \in N,k \in M;$

(6)

$ \qquad\sum\limits_{i = 1}^n {{v_i}{y_{ik}}} {\text{≤}} {Q_k}, \forall k \in M;$

(7)

$ \qquad t_i^1 {\text{≤}} {T_{ik}} {\text{≤}} t_i^2, \forall i \in N,k \in M;$

(8)

$ {T_{jk}} {\text{≥}} {T_{ik}} + {t_{ij}} + {l_i} + R({x_{ijk}} - 1) ,\;\; \forall i \in N,j \in N \cup \{ 0\} ,k \in M;$

(9)

$ \qquad {T_{0k}}{\text{≤}}{T_i} + R(1 - {y_{ik}}), \qquad \forall i \in N,k \in M;$

(10)

$ \qquad {x_{ijk}},{y_{ik}},{{\textit{z}}_k} \in \{ 0,1\},\;\;\;\; \forall j \in N \cup \{ 0\} ,i \in N,k \in M{\text{。}}$

(11)

式(1)表示最小化成本的目标函数，其中第1项为固定成本，第2项为油耗成本；式(2)表示当网点 $i$ 由车辆 $k$ 服务时，车辆 $k$ 必定会从网点 $i$ 到其他某一点；式(3)表示当网点 $j$ 由车辆 $k$ 服务时，车辆 $k$ 可能会从某一网点到网点 $j$ ，或者网点 $j$ 正是车辆 $k$ 的起始点；式(4)表示每个网点被且仅被一辆车服务；式(5)表示被使用的车辆一定会从某一网点前往分拨中心；式(6)表示若车辆 $k$ 在网点 $i$ 处被限行，则网点 $i$ 不能由车辆 $k$ 服务；式(7)表示每辆车所服务的网点的总货量不能超过其最大容积；式(8)表示每个网点开始装货的时间窗约束；式(9)表示前后2个点开始装货(对于分拨中心则是卸货)的时间关系；式(10)表示每辆车到达分拨中心的时间约束，用于错峰；式(11)表示决策变量类型。

用CPLEX等求解器求解该模型，可以得到关于模型的最优解。但是对于NP-hard问题，只要规模稍大，就难以在短时间内得到较好的解。

3 启发式算法

结合快递业务的相关经验设计一个构造性启发式算法，能在短时间内给出一个较好的解。该算法基于2个基本思想。

1) 优先使用大型车。首先，所用车辆数越少，总成本也越低，而使用大型车能够减少总车辆数。其次，当车的容量有较多富余时，考虑能否缩小车型。因为，一方面对单辆车来说，车型越小，成本越低；另一方面，缩小车型意味着可以访问一些对大车型限行的网点。

2) 对于每一辆车要访问的网点集，都是从一个网点开始，以与已确定的网点集的行驶距离之和作为邻域搜索优先指标选择新加入的网点并检验其可行性，从而不断扩张这个网点集。由同一辆车服务的网点间的行驶距离越短，越容易满足时间窗约束，成本也会越低。

记单辆车的访问网点集 $A$ (或点 $a$ )的总货量为 ${V_A}$ (或 ${V_a}$ )，选择的车型为 ${r_A}$ ， $\left| A \right|$ 表示集合 $A$ 中的元素个数，车型 $r$ 的最大容积为 ${H_r}$ 。该算法的具体步骤如下。

步骤1 　将所有网点按与分拨中心的行驶距离降序排列，记 $S$ 为排序后的有序集合。

步骤2 　记 $S$ 中的第一个点为 $a$ ，将 $a$ 从 $S$ 中删去，记在 $a$ 处可通行的车型集合为 $\{ {r_1},{r_2},\cdots,{r_t}\} $ (按最大容积降序排列)。令 $i \leftarrow 1,A \leftarrow \{ a\} ,{r_A} \leftarrow {r_i}$ 。

步骤3 　从集合 $S$ 中取出所有未对车型 ${r_A}$ 限行的点，组成集合 $S'$ ，若 $S' = \varnothing $ ，则转到步骤8，否则继续。

步骤4 　令 $S'' \leftarrow \{ s \in S'|{V_s} {\text{≤}} {H_{{r_A}}} - {V_A}\} $ ，若 $S'' = \varnothing $ ，则转到步骤8，否则继续。

步骤5 　对集合 $S''$ 中的点排序，使得下列关系成立。其中 ${S''_i}$ 代表 $S''$ 中的第 $i$ 个点， ${\rm{dist}}({p_i},{p_j})$ 代表从点 ${p_i}$ 到点 ${p_j}$ 的行驶距离。

$\begin{split} &\qquad\sum\limits_{p \in A} {{\rm{dist}}({{S_1^{''}}},p)} {\text{≤}}\sum\limits_{p \in A} {{\rm{dist}}({{S_2^{''}}},p)} {\text{≤}} \cdots {\text{≤}} \\ &\sum\limits_{p \in A} {{\rm{dist}}({{S_i^{''}}},p)} {\text{≤}} \cdots {\text{≤}} \sum\limits_{p \in A} {{\rm{dist}}({{S_{\left| {S''} \right|}^{''}}},p)}{\text{。}} \end{split}$

步骤6 　记 $S''$ 中的第一个点为 $b$ ，将 $b$ 从 $S''$ 中删去，检验集合 $\{ b\} \cup A$ 是否存在能满足时间窗约束的访问序列，若存在，则令 $A \leftarrow \{ b\} \cup A$ ，将 $b$ 从 $S$ 和 $S'$ 中删去，更新 ${V_A}$ 并转到步骤4，否则，继续。

步骤7 　若 $S'' \ne \varnothing $ ，则回到步骤6，否则，继续。

步骤8 　如果 $i {\text{≤}} t - 1,{V_A} {\text{≤}} {H_{{r_{i + 1}}}}$ 且车型 ${r_{i + 1}}$ 在 $A$ 中的所有网点处均能通行，那么令 $i \leftarrow i + 1,{r_A} \leftarrow {r_i}$ 并转到步骤3，否则继续。

步骤9 　对集合 $A$ 中的点排序，使得访问序列 ${A_1},{A_2}, \cdots ,{A_{\left| A \right|}},0$ 满足时间窗约束，且总行驶距离最小，从而得到一辆车的调度安排 $({A^*},{r_A})$ ， ${A^*}$ 表示最优的网点访问排序。

步骤10 　若 $S \ne \varnothing$ ，则转到步骤2，否则，算法终止，得到最终的解。

在检验时间窗(步骤6)及确定最优访问顺序(步骤9)时，可以直接枚举全排列。在问题描述中已经提到，时间窗的约束非常严格，所以每辆车能访问的网点数不会太多，一般不超过4个。这是快递同城运输中固有的特点。因此即使采用穷举，计算量也不会太大。

4 列生成

定义 $\omega = (A,r)$ 为路径−车型对， $r$ 表示车型， $A$ 是一个网点序列，代表要访问的顺序。令 $ \varOmega $ 为所有可行 $ \omega $ 的集合。

4.1 集分割模型

对于 $ {\omega _j} \in \varOmega $ ，令 ${c_j}$ 表示其总成本。引入指示参数 ${a_{ij}}$ ，当 ${\omega _j}$ 会访问网点 $i$ 时，令 ${a_{ij}} = 1$ ，否则为0， $i \in N$ 。决策变量设为 ${x_j}$ ，当 $ {\omega _j} \in \varOmega $ 被采用时， ${x_j} = 1$ ，否则为0。那么可以建立如下模型。

$\qquad\min \; \sum\limits_{{\omega _j} \in \varOmega } {{c_j}{x_j}} {\text{。}}$

(12)

$ \begin{split} & \qquad\;{\rm{s}}{\rm{.t}}{\rm{.}} \\ &\qquad \sum\limits_{{\omega _j} \in \varOmega } {{a_{ij}}{x_j}} = 1, \; \forall i \in N; \end{split} $

(13)

$\qquad {x_j} \in \{ 0,1\}, \; \forall {\omega _j} \in \varOmega {\text{。}}$

(14)

式(12)为成本目标函数；式(13)表示每个网点都被且仅被访问一次；式(14)表示变量类型。

当该模型的变量数较大时，求解速度比较慢，需要用到列生成技术。在单纯形法中，对于变量较多的模型，最优解中多数变量都是值为0的非基变量，只有小部分值不为0的基变量。列生成的基本思想就是不断生成能优化目标值的列，即寻找检验数为负的变量(对于最小化问题)。在这里， ${\omega _j}$ 就相当于列生成中的“列”，它蕴含了 ${c_j}$ 、 ${a_{ij}}$ 等信息。

4.2 限制性主问题

利用前面介绍的启发式算法生成一个初始解，并转化为初始 $\omega $ 集合 $ {\varOmega _0} $ ，再将决策变量线性松弛，从而能得到如下限制性主问题(restricted master problem, RMP)。

$\qquad\min \; \sum\limits_{{\omega _j} \in {\varOmega _0}} {{c_j}{x_j}}{\text{。}} $

$ \begin{split} &\qquad\;{\rm{s}}{\rm{.t}}{\rm{.}}\\ &\qquad \sum\limits_{{\omega _j} \in {\varOmega _0}} {{a_{ij}}{x_j}} = 1, \; \forall i \in N{\text{；}} \end{split} $

$\qquad 0 {\text{≤}} {x_j} {\text{≤}} 1 ,\quad \forall {\omega _j} \in {\varOmega _0}{\text{。}}$

求解上述模型后，能得到对偶变量 $\pi = ({\pi _1}, $ ${\pi _2}, \cdots ,$ ${\pi _n}) = {c_B}{B^{ - 1}} $ ， $B$ 表示最优解中的基， ${c_B}$ 表示目标函数中对应的系数。

4.3 子问题

以单纯形法中的检验数为目标函数，确定列生成中的子问题。

$\qquad\min \; {\textit{z}} = {c_j} - \sum\limits_{i = 1}^n {{\pi _i}{a_{ij}}}{\text{。}} $

$ \begin{split} &\qquad{\rm{s}}{\rm{.t}}{\rm{.}} \\ &\qquad{\omega _j} \in \varOmega {\text{。}} \end{split} $

当最优目标值 ${{\textit{z}}^*} {\text{＜}} 0$ 时，说明子问题产生的列能减小RMP的目标值，因此要把这个列加到RMP中，然后继续求解，循环迭代；否则，说明已经没有更好的列。若此时RMP的列集合为 $ \overline \varOmega $ ，再求解 $ \overline \varOmega $ 对应的集分割模型即可。由于原问题是整数规划，故最后得到的解是近似最优解，若原问题是线性规划，那么得到的解即为最优解。

子问题是列生成算法中的关键。上述子问题是一个带资源约束的初等最短路问题，Feillet等^[17]对此提出了一个基于动态规划的精确算法。该算法的原理是通过对图中各个点标号来建立点与点之间的联系，并不断扩展这些标签以检查资源约束，从而获得符合要求的 $\omega $ 。为了减小计算量，通常不会每次都求到最优解，只要计算出足够多目标值为负的可行 $\omega $ 即可。只有在最后一次迭代时要求出目标值最小的可行 $\omega $ ，此时要遍历所有可行 $\omega $ 并确认它们的目标值均不为负。

由于本问题的约束非常严格， $ \varOmega $ 的规模不会太大，所以本文事先求出所有的可行 $\omega $ ，其优点是能避免在子问题中对相同 $\omega $ 的可行性进行重复检验。另外，还可进一步缩小 $ \varOmega $ 的规模，剔除一些最优解中必然不会采用的 $\omega $ 。由于 $\omega = (A,r)$ ，故考虑分别通过 $A$ 和 $r$ 两部分来缩减 $\omega $ 的数量。

1) 对于网点集合 $S$ ，可能有多种访问序列 ${A_i}$ 能与车型 $r$ 构成可行的 $({A_i},r)$ ，此时只保留访问顺序最优的 $({A^*},r)$ ，即总行驶距离最短的，并舍弃其他的 $({A_i},r)$ 。

2) 对于网点集 $S$ 的最优访问序列 ${A^*}$ ，可能有多个车型 ${r_i}$ 能与其构成可行的 $({A^*},{r_i})$ ，此时只保留车型最优的 $({A^*},{r^{\rm{*}}})$ ，即成本最低的，并舍弃其他的 $({A^*},{r_i})$ 。

遵循这2个原则，能大大缩小 $ \varOmega $ 的规模。记 $ {W_i} = \{ \omega = (A,r) \in \varOmega |\left| A \right| = i\} $ ，可通过如下算法求出 $ \varOmega $ 。

步骤1 　初始化 ${W_0} = \{ (\varnothing ,r)\} ,k = 0$ ， $r$ 可以是任意车型；

步骤2 　以 ${W_k}$ 为基础，对 $\forall \omega = (A,r) \in {W_k}$ ，考虑往 $A$ 中加入任意一个网点，检验可行性，若可行则确定最优访问序列及车型，从而得到 ${W_{k + 1}}$ ，令 $k \leftarrow k + 1$ ；

步骤3 　若 ${W_k} \ne \varnothing $ ，则返回步骤2，否则算法终止，得到 $ \varOmega = \bigcup\limits_{i = 1}^{k - 1} {{W_i}} $ 。

得到 $ \varOmega $ 后，在每次求解子问题时，只要遍历 $\forall {\omega _j} \in \Omega $ ，就可以快速计算出目标值为负的 $\omega $ ，而无需检验可行性。

4.4 算法实现

引入以下符号。

$e$ ：一个较小的常数，表示检验数允许的绝对误差；

$t$ ：常数，表示每次迭代中生成列的数量；

${\rm Heur}({\rm{ }}) $ ：表示前面介绍的启发式算法，返回对应的初始子集 $ {\varOmega _0} $ ；

${ \rm RV}({\rm{ }}) $ ：表示前面介绍的计算全集 $\varOmega $ 的算法，返回 $ \varOmega $ ；

${\rm RMP}({\varOmega _1}) $ ：用来求解列集为 $ {\varOmega _1} $ 的RMP，返回对偶变量 $\pi $ ；

$C(\pi ,\omega )$ ：返回子问题中对偶变量为 $\pi $ ，列为 $\omega $ 时的目标值；

${\rm OP}({\varOmega _1})$ ：求解列集为 $ {\varOmega _1} $ 的原问题(集分割模型)。

列生成算法过程如下。

1) $S \leftarrow {\rm Heur}()$

2) $ P \leftarrow {\rm RV}({\rm{ }}) $

3) ${\rm{repeat}}$

4) 　　 ${\rm{ }}\pi \leftarrow {\rm RMP}(S)$

5) 　　 ${\rm{ }}i = 0$

6) 　　 ${\rm{ }}T{\rm{ = }}\varnothing $

7) 　　 ${\rm{ for\;all }}\;\omega \in P$

8)　　 ${\rm{ if }}\;C(\pi ,\omega ) {\text{≤}} - e$

9)　　　　　 ${\rm{ }}i \leftarrow i + 1$

10) 　　　　 $T \leftarrow T \cup {\rm{\{ }}\omega {\rm{\} }}$

11)　　　　 ${\rm{ if }}\;i = = t$

12)　　　　 ${\rm{break}}$

13) 　　 ${\rm{ }}S \leftarrow S \cup T$

14) ${\rm{until }}\;i = 0$

15) ${\rm OP}(S)$

5 实验及结论

实验的运行环境为Intel(R) Core(TM) i5-5200U CPU @ 2.20 GHz 2.19 GHz。实验采用某快递公司在上海的真实数据，网点有120个，车型有4种。任意两点间的行驶时间、行驶距离通过高德API来调取，加上分拨中心，一共有12 100(121×100)条数据。设计12组实例，网点数量分别为10, 20, ···, 120。除了最后一组，每一组都是从这120个网点中随机取点，然后分别用数学模型、启发式算法和列生成3种方法进行求解。对每一组中的每一个方法，均求解5次并对成本和运行时间取平均值作为最后的结果。成本保留两位小数，运行时间保留三位小数。数学模型均调用CPLEX进行求解，设定最长运行时间为600 s，若到600 s还未求解完则直接输出当前解。

以网点数为30的组为例，用列生成求解所得的调度方案如表1所示，一共需要16辆车，平均装载率为0.90。

统计各组计算完成后的成本和运行时间，实验结果如表2所示，据此分析各算法的特点。数学模型求解时，即使只有10个点也无法在600 s内得到最优解，而且会占用大量内存；启发式算法耗时极短，且不会随着规模的提高而显著增加，其优化结果与列生成相比，平均相对误差为8.86%；列生成在每组中的成本都是最低的，例如第6组中，其成本仅为数学模型的56%，启发式算法的92%，耗时也不算长，在120个网点的规模下，只需175 s就能得到近似最优解。

实验结果显示，列生成算法能有效解决快递业的同城运输问题。如果问题规模继续扩大，可以给列生成设定最长运行时间或迭代次数，达到条件时输出当前解，由于列生成以启发式算法作为初始解，所以有上界保证。

本文采用的行驶时间数据是通过高德API提前获取的，而不是实时的。但事实上，城市内的交通状况与时段有很大的关系，部分路段在早晚高峰期会比较拥堵，使得实际行驶时间变长。在后续的相关研究中，可以考虑增加交通因素。