推文作者：张晨皓，东南大学在读博士生，
研究方向：时空网络建模，车辆-无人机路由优化，精确算法与分解技术；(Bilibili：交通小张爱摸鱼)
编者按
本文发表于2015年，是卡车无人机路径研究的开山之作。作者针对单车辆单无人机的配送场景，首次提出了两种带无人机的旅行商问题变体：卡车与无人机串联协作的FSTSP和并行协作的PDSTSP。文章明确给出了基于弧的混合整数线性规划（MILP）模型，并通过设计基于成本节约的大邻域搜索算法，对10至20个节点的实例进行了数值计算。其建模思路对后续一系列相关研究产生了深远影响，尤其在当前低空经济快速发展的背景下，仍具有重要的参考价值。
1.引言
优化理论和计算机科学的发展，逐渐打破了不同交通运输方式之间的壁垒。随着相关技术的军转民，无人机（UAV，unmanned aerial vehicles）与其他运输方式的联合应用逐渐进入公众视野。这种多式联运凭借其在协调统一、便捷可靠、降本增效等方面的天然优势，吸引了越来越多从业者的关注，特别是在对成本、效率高度敏感的城市物流运输和路由领域。同时，由于现实中存在自然或法律的限制，多种运输模式之间的有效协调与配合也成为无法回避的重要课题。
相较于其他模式的组合，车辆与无人机，在载重量、续航、速度、能耗污染方面具有高度互补、彼此成就的特点：在现实层面，车辆携带无人机，可以进一步扩大其飞行范围，及时对无人机进行货物的装载、更换电池或维修，无人机亦可以缩短车辆的服务时长。在理论层面，相较于卡车-拖车问题，车辆与无人机两类运输模式路网更加分离，载具具有独立动力，更具备两级车辆路由问题（2E-VRP，Two-echelon Vehicle Routing Problem）的特征。
随着研究发展，这类问题产生了丰富的变体。常被称为卡车无人机路由协同问题（TDRP，Truck-Drone collaborative Routing Problem），可根据问题中执行任务的基本单元——卡车组（truck group）的运行情况（由单个卡车组，还是多个卡车组提供服务），分为带无人机的旅行商问题（TSPD，Traveling Salesman Problem with Drone）和带无人机的车辆路由问题（VRPD，Vehicle Routing Problem with Drone），这篇文章的研究内容关注于单个卡车组，可以简单归并为TSPD。
2.问题描述
此类问题的基本描述为：给定一个有向完全图，其中节点集合由客户集合与起终点车场和构成，弧集合中各弧的权重取决于具体问题场景。由卡车与无人机构成的卡车组从起点车场出发，最终回到终点车场，期间卡车与无人机相互配合，共同服务客户集合的全部需求。常见的目标函数为最小化完工时间，或最小化运营成本。为帮助大家更好地理解该问题，下面结合现有的研究，对于前述高光部分进行简单梳理和解释。
有向完全图（建模所基于的网络）：此类问题常见的为空间客户服务网络，即网络单纯由客户节点和车场构成，一般是通过二次加工提取最短路数据之后的网络，因此一般为完全图形式；此外常见的还有交通运输网络，即最原始未加工的交通网络，由于为对最短路数据进行计算提取，一般以邻接矩阵的形式给出。本文使用的是空间客户服务网络。 客户集合（客户属性）：此类问题由于一些现实因素的考虑，例如客户存在人工签收环节/包裹重量大于无人机容量，因此该类客户对于无人机是不可访问的（ineligible），只能通过卡车进行服务。无人机可访问的（eligible）客户集合是客户集合的子集。本文考虑了不可由无人机访问的客户。 起终点车场（车场属性）：此类问题最常见的为单车场问题（也有多车场问题），终点车场只是为了在索引上区别于起点车场而引入的虚拟车场，究其本质是为了在空间网络建模中向破子圈的MTZ约束服务。本文是单车场问题。 弧集合（弧属性）：此类问题弧权重一般为距离成本、时间成本、经济成本，特性包括固定、随机和时变，在一般的建模中，对于某条弧，卡车与无人机被认为具有不同的成本，这也是两种方式差异的体现。本文使用固定的时间成本表示弧成本。 卡车组（卡车组属性）：一个卡车组由一辆卡车和一架或更多的无人机构成（ 即卡车：无人机或），在所研究的问题中，可存在一个卡车组或更多的卡车组。本文关注单卡车组，且卡车组由辆卡车和架无人机构成。 相互配合（起降、服务模式）：此类问题无人机从卡车发射，执行完任务之后由卡车回收，发射与回收可在客户节点或专用站点进行，早到者需要等待晚来者；此外，无人机可由发射它的对应卡车回收，也可由其他卡车回收；无人机在发射后，可以执行配送、取货、取送货，可访问一个客户，多个客户（multi-visit/drop），单行程，可重复在某点起降的多行程（multi-trip），且无人机飞行时间不得超过其续航。本文发射和回收发生在客户节点，由发射它的对应卡车进行，配送问题，且无人机一次只能访问一个客户，单行程。 全部需求（需求特点）：此类问题客户的需求特点主要可分为配送需求、取货需求、或取送货需求，需求量可为定值也可随机，同时也可存在时间窗约束。本文客户需求为配送需求，未考虑时间窗，需求量仅影响无人机可服务的客户集合。 最小化完工时间（目标函数）：此类问题的目标函数主要包括最小化完工时间、最小化运营成本、最小化污染排放、最大化收益、最大化客户满意度等。本文目标函数为最小化完工时间，即卡车组中卡车和无人机返回终点车场的最晚时间。
此外，越来越多的文献将包裹载重考虑在无人机续航的影响因素中，本文并未对其进行深入研究。
本文对于所关注的单车单无人机的场景，重点讨论了车辆与无人机的两种协作模式：（1）在客户点距离起点车场较远时，由于续航限制，无人机无法直接从车场起飞执行服务，而是需要被卡车载运至合适的范围再前去执行服务，这类问题在文中被定义为The flying sidekick TSP（FSTSP）；（2）在客户点距离起点车场较近时，无人机可以直接从车场起飞，不需要借助卡车，而是与其并行工作，这类问题在文中被定义为The parallel drone scheduling TSP（PDSTSP）。文章重点对FSTSP进行了建模与求解，下面分别对这两类问题进行说明。
2.1 FSTSP问题描述
当客户点距离起点车场较远时，无人机需要被卡车载运然后再去提供配送服务，两者从车场出发，最终回到车场，期间相互协作配合服务所有客户的配送需求，目标函数为最小化整体的完工时间。卡车与无人机在时间和空间层面存在同步行为（发射、回收、等待），此类问题被本文称之为The flying sidekick TSP（FSTSP）。后续的较多研究皆是以此类问题为基础，本文主要对该问题进行了详细的建模与求解。下图分别展示了（a）不带无人机的单辆卡车最优配送路线，与（b）无人机服务两个eligible客户的FSTSP场景。
上述两种情况对应的完工时间如下图所示。
2.2 PDSTSP问题描述
当客户点距离起点车场较近时，无人机和卡车可以各自从车场出发，无人机以车场作为起飞和降落的依托，无需与卡车在时间和空间层面存在同步行为（发射、回收、等待），目标函数同样为最小化整体的完工时间。此类问题被本文称之为The parallel drone scheduling TSP（PDSTSP），这类问题本质上是一个相对简单的并行机调度问题。下图分别展示了（a）不带无人机的单辆卡车最优配送路线，（b）不进行优化下，无人机和卡车各自服务职能范围内客户的配送路线，与（c）以最小化完工时间为目标，优化后的结果。
上述三种情况对应的完工时间如下图所示。
3.模型构建
3.1 模型假设
在建模前，需要解释说明本文所使用的一个概念，sortie（出击、架次）：
一个sortie通常由三个节点的索引组成：，用来描述无人机的一次“起飞-服务客户-降落”的过程，即无人机从发射点在卡车上起飞，前往客户点执行配送服务，而后返回在点被卡车回收。所有的sortie均在建模前被提前构建好，这类似于的建模思想，也类似于Carlsson（2017）文章中 和MirHassani（2013）文章中的概念，这种集成的概念也在2018年Agatz的文章中被进一步建立为基于的建模方式，而后被广泛借鉴。
模型假设如下：
每个sortie只能服务一个客户，但在无人机的服务过程中，卡车可以服务多个客户； 无人机在一个sortie中处于持续飞行的状态，无人机不得在中途临时降落以节省电量； 若无人机在某个客户节点被卡车回收，则其可能在被重新发射，但如果在某客户节点被发射，则可能无法返回点的卡车； 无人机的发射和降落发生在卡车当前服务的客户位置（此外，起点车场只能进行发射、终点车场只能进行回收），卡车所经过的弧中间位置是无法进行发射和降落的，且卡车不会重新访问一些节点来回收无人机； 除了车场之外，无人机和卡车均不得访问其他的任何非客户节点（客户服务网络）； 若无人机的架次在终点车场结束，则其退出服务不得再进行重新发射。此外，从起点车场直接进行发射的场景更加贴近PDSTSP这种并行的问题。
3.2 FSTSP问题建模
本次推文主要针对串联的FSTSP模型展开，其参数、变量分别整理如下：
其中，卡车路由决策变量为0-1变量，若卡车从节点行驶至节点，则取值为1，否则为0；
无人机路由决策变量为0-1变量，若无人机从节点发射，前往节点，而后返回节点，则取值为1，否则为0；
为本文用于破子圈约束的辅助变量，表示了被服务的节点次序，为保证优化效果，可给定一个合理的较紧值，如范围定义为；
和分别为无人机和卡车在节点的有效到达时间，区别于到达节点那一瞬间的时间，所谓“有效”是指到达节点且完成一些具体事项之后的时间，下文存在与此照应的约束；
是用于判断卡车所访问客户点和的先后而引入的0-1辅助决策变量，若先于被卡车访问，则，否则。
FSTSP的具体模型如下:
对于目标函数（1）：可以看到，该问题追求最小化卡车到达终点车场的有效时间，通过约束（12）-（15）可以实现卡车与无人机在发射与回收两个关键节点上，有效时间的对齐，因此这里的有效时间实际上是一种带有“卡车组”底色的到达时间，及最晚到达的那一个成分（卡车、无人机）的时间，其实写为 也是等价的，因为两者已经通过约束实现了对齐。
对于约束条件，小编将其概括为客户、卡车、无人机、同步四种类别的约束，如下图所示。此外，还有一些其他关于辅助决策变量取值的约束。当然，正如Drexl（2012）所言，两种模式之间的同步正是此类问题最为核心的挑战和难点，为模型带来了庞大且复杂的偶和约束，下面分类别向大家进行介绍。
客户：约束（2）为客户服务约束，要求对于任何一位客户，必须要被访问（服务），该访问要么被卡车提供，要么被无人机提供；
卡车：约束（3），（4）及（6）分别为卡车的网络流平衡约束，约束（3）要求单位的卡车流需要从起点车场流出，最终流入终点车场，期间在各个中间点满足流入和流出相等的流平衡约束；约束（5）是形式的破子圈约束，可以破除任何环路的存在，此处的取值为合理的最紧值，在辅助变量取值为的情况下，卡车路径访问的节点依次对应的会形成一组首项为，公差为的递增数列；
无人机：约束（7），（8）为无人机在发射/回收节点的流约束，分别表示无人机可以从非终点车场流出，无人机可以从非起点车场流入；约束（17），（18）为无人机访问/回收节点的时间戳约束，一种析取约束，分别表示无人机访问客户的时间戳应符合时间逻辑（不早于起飞时间戳+前往服务的飞行时长），无人机回到卡车的时间戳应符合时间逻辑（不早于访问客户的时间戳+返回的飞行时长+回收用时）；约束（19）表示若无人机选择$sortie$，则其对应的飞行时间不得超过其续航时长；,j,k>
同步：约束（9），（10）为无人机发射/回收点的卡车访问约束，分别要求对于非起点发射的sortie，车辆必须经过无人机的发射点和回收点，而对于在起点发射的sortie，车辆必须经过无人机的回收点；约束（11）为无人机sortie内卡车访问次序约束，要求在无人机sortie内，卡车一定是先访问发射点，再访问回收点；约束（12），（13）为无人机发射点的有效时间约束，要求无人机在发射节点的有效时间等于卡车在该点的有效时间，共同实现了卡车与无人机在发射节点时间上的对齐；约束（14），（15）为无人机回收点的有效时间约束，要求无人机在回收节点的有效时间等于卡车在该点的有效时间，共同实现了卡车与无人机在回收节点时间上的对齐；约束（16）为卡车的有效时间戳约束，要求卡车访问客户点的有效时间戳应当符合时间逻辑；约束（23）无人机为若干sortie之间的次序约束，对于多个无人机sortie，一定是上一个先回收，下一个再发射，任意两个sortie之间无交叉；
辅助决策变量：约束（20）-（22）定义了卡车访问客户点次序变量；约束（24）-（32）给出了一些辅助决策变量的初始值与其定义域。
讨论与思考：在这种建模框架下，对于空间层面的同步已经做了较为详尽和紧凑的约束，但是对于时间上仍较为粗略和松弛，例如上图分别在空间和时间上描述了一个卡车-无人机的路径，可以看出该模型仅仅是“锚定”了一些关键的时间点（如,,,这些时间点，要求卡车与无人机在时间上来一个“对表、校正、对齐”的效果），而对于其他一些琐碎的时间节点，如发射和降落的开始和结束时间，其实文章并未给出约束，允许在一定范围内滑动，事实上，大部分研究也未能给出准确和紧凑的描述，在时间层面对于该问题的研究相对还不够充分。
3.3 PDSTSP问题建模
相较于卡车与无人机存在同步和耦合的FSTSP，并行的PDSTSP则更简单一些，适用于较多客户距离起点车场较近（在无人机的续航范围内）的情况，下面对其模型进行相对简单的解释和说明：
模型引入了三个新的决策变量，分别为无人机访问决策变量，当客户由无人机服务时等于1，否则为0；卡车路由决策变量，当卡车从节点前往节点时等于1，否则为0；以及新的破子圈辅助决策变量，其他符号若无强调则与FSTSP的建模一致。
可以看到，PDSTSP的模型缺少了复杂的、耦合的、关于两种模式同步的约束。
目标函数（33）与FSTSP的目标函数相同，最小化最早完工时间（无人机与卡车返回终点车场的最晚时间）；而约束（34）和（35）保证了可以作为无人机和卡车返回终点车场的最晚时间；约束（36）确保每个客户被服务一次，服务可由卡车和无人机提供；约束（37）-（39）为卡车的网络流约束，要求单位的卡车流需要从起点车场流出，最终流入终点车场 ，期间在中间节点满足流平衡约束；约束（40）为卡车流的破子圈约束；约束（41）-（43）给出了决策变量和辅助决策变量的取值范围。
4.算法设计
对于FSTSP，由于此类问题的NP-hard特点，即使是看似简单的10个客户以上的问题，MILP求解器通常也需要几个小时才可以取得最优解，因此本文设计了基于的启发式算法，这实际上是一种在贪婪思想指导下的、基于单一解的、通过“移除-插入”来实现路线重新排列的大邻域搜索算法，小编将本文所使用的具体算法的分类简单整理如下。
4.1 针对FSTSP的启发式算法
本文所提出的针对FSTSP的启发式算法主函数的伪代码如下所示：
上图展示了主函数的伪代码，其中较为关键的五个子函数分别为：
：用于获得算法的初始可行解，该可行解在文中默认为针对客户集合的、只使用卡车的TSP解，包括节点的访问序列及其访问时间序列，同时也会生成初始的、用于后续插入的子路径集合，该函数我们通过蚁群算法和2-邻边交换来实现； ：针对节点的移除阶段，优先关注可以被无人机服务的（UAV-eligible）客户 ，以遍历的方式计算在卡车路径中，删除掉节点所对应的增益，该子函数的伪代码如下，主要给出了不同情况下，移除该节点的计算规则；
：针对节点的插入阶段，将被移除的节点插入到卡车子路径中，当待插入的卡车子路径中存在无人机，则计算在当前子路段中使用卡车来服务节点的最佳插入位置，服务方式，及其净增益 ，这里也会考虑卡车服务该客户对既有的无人机飞行时间的影响，该子函数的伪代码如下所示；
：与子函数3相同，针对节点的插入阶段，将被移除的节点插入到卡车子路径中，当待插入的卡车子路径中不存在无人机，则计算在当前子路段中使用无人机来服务节点的最佳插入位置，服务方式，及其净增益 ，这里主要会考虑无人机的续航，该子函数的伪代码如下所示；
：根据覆盖保存的结果, 对应地调整的插入位置和服务方式, 更新子路径集合，该子函数的伪代码如下所示。
其实这种移除-插入启发式算法，是一种破坏-修复的大邻域搜索思想，以净增益 为指标，贪婪地指导路径的重新排列，且该算法有一个启发式的默认，那就是优先选择使用无人机，认为无人机对于完工时间的优化效果要高于卡车，具体算法实现的示意图如下所示。
4.2 针对PDSTSP的启发式算法
相比于FSTSP，并行的PDSTSP的算法更加简单，算法的伪代码如下所示。
PDSTSP将旅行商问题（TSP）与并行机调度问题（PMS）相结合，旨在优化卡车与无人机对客户的交付。初始解中，将所有对无人机可访问的客户均分配给无人机服务，剩余客户由卡车处理。首先解决PMS以分配客户给无人机，然后解决TSP以确定卡车的路线。接着，通过迭代重新分配客户，用来平衡无人机与卡车的完工时间，依然使用净增益 为指标，直到无法再改进为止，关于子函数的伪代码如下所示。
5.数值实验
5.1 FSTSP的实例计算结果与分析
由前述分析可知，FSTSP启发式框架依赖于反复求解TSP的机制。因此比较了四种候选TSP求解方法。第一种是通过整数规划求解器（如Gurobi，使用“IP”表示）最优解决TSP子问题。但该种方法对于实际规模的问题可能计算效率较低。因此，还测试了三种常见且易于实现的TSP启发式算法：
第一种是Clarke-Wright提出的基于启发式的修改版。在此过程中，每个客户最初直接连接到车场。然后选择能够提供最大“”的客户对，依次连接这对客户，不断进行该步骤直到形成有效的TSP解； 第二种是“”方法，从起点车场开始，依次选择离当前位置最近的客户形成有效的TSP解； 第三种是“”方法。该过程从起点车场（初始焦点）发出光线。路线从车场开始，顺时针旋转光线。通过不断改变焦点、初始角度和光线的旋转方向，可以获得不同的解决方案。
生成了72个测试问题的集合，每个问题包含10个客户，分布在一个8英里见方的区域，车场位置随机选择为客户的x和y坐标的平均值（即重心附近）。80%到90%的客户被指定为UAV-eligible。无人机的续航选择为20、40分钟，速度为15、25、35英里/小时，采用欧氏距离计算弧成本。卡车速度设为25英里/小时，采用曼哈顿距离计算弧成本。发射和回收的时间消耗和假设为1分钟，求解结果如下图所示。
直接使用Gurobi对FSTSP进行求解的“FSTSP formulation”在72个问题上都超过了1800s的求解时长，这也导致许多测试实例中出现了负的gap，表示启发式常常找到比Gurobi更好的解。不同的TSP启发式在解的质量上存在显著差异，具体求解的细节见表1。
结果表明有效的TSP启发式有助于提升FSTSP启发式的整体表现，特别是，启发式表现良好，数值分析表明，高质量（即接近最优）的TSP解对所提FSTSP启发式框架的解质量有积极影响。
5.2 PDSTSP的实例计算结果与分析
由于“”算法在FSTSP的实例计算中表现十分有限，因此本文在PDSTSP的案例中并未选择使用该算法，而是选择使用Gurobi以及其他前述的启发式算法。每个问题限时3分钟。Gurobi能够为所有10客户问题找到最优解，但在360个20客户实例中有87个超出了时间限制，不同算法的求解结果如下图所示。
其中横轴的标签分别表示底层TSP和PMS子问题的具体求解方法，例如，“Savings/LPT”表示分别通过和来求解TSP和PMS问题，具体的求解细节图表2所示。
通过计算结果我们观察到，首先，观察到使用精确（整数规划）方法解决单个TSP和PMS子问题（标记为“IP/IP”）能够获得接近最优的结果。然而，对于20客户实例，反复求解TSP子问题所需的运行时间显著超过了Gurobi求解完整PDSTSP的时间限制。相比之下，使用启发式解决TSP子问题所需的运行时间显著更短，同时在最优性差距方面产生了具有竞争力的解。这在使用“IP”或求解PMS子问题时均成立。事实上，对于给定的TSP解法，通过获得的解质量通常与通过IP公式最优解决PMS子问题的解相同。此外，在表2中观察到，通常比IP方法在解决PMS子问题时更快。对于PDSTSP的求解结果，本文强调了TSP启发式的重要性，并表明似乎能够为PMS子问题提供接近最优的解。
5.3 无人机速度与续航的敏感性分析
该部分主要对于无人机速度和续航的组合，进行了敏感性分析。设计了90个实例，包含25、50和75个客户，考虑了不同的仓库位置。测试了九种速度与续航的组合，从每小时20英里、续航20分钟到每小时40英里、续航40分钟不等，下图展示了PDSTSP的计算结果。
实验结果显示，对于PDSTSP，速度的提升显著改善了性能，特别是对于三架无人机的测试案例。例如，在无人机飞行范围为13.33英里的场景中，速度为每小时40英里（续航20分钟）时的表现明显优于速度为每小时20英里（续航40分钟）。同样，在固定续航为20分钟的情况下，UAV以每小时40英里飞行的性能显著优于以每小时30英里飞行。下图展示了FSTSP的计算结果。
FSTSP实例的结果显示出类似的趋势。较低的速度，即便在高续航情况下，也限制了无人机可访问的客户数量。总体而言，与传统配送方式相比，PDSTSP的优势更加明显。这是因为PDSTSP对于卡车和无人机之间同步的要求较为宽松，不需要FSTSP中的高度同步。需要注意的是，在FSTSP的设定中，无人机必须在客户位置起飞或返回，而不能在中间位置，论文指出，这一限制的放松有潜力显著节省整体的配送时间。
5.4 FSTSP与PDSTSP的对比
PDSTSP和FSTSP分别针对两种不同的操作场景。PDSTSP适用于大部分客户位于无人机飞行范围内的情况，而FSTSP则适合仓库位置偏远的场景，如下图所示。
然而，在某些情况下，选择哪种建模方法并不具有明显的、绝对的判别标准。两种服务模式各有优劣，分析结果显示，虽然PDSTSP在客户可达性强时表现出色，但在客户数量较少或仓库位置偏远时，FSTSP可能更具优势。因此，未来的研究可以关注如何根据不同的操作特征来选择最优策略，包括客户分布、无人机数量、速度和续航等多个因素。这将有助于更好地理解FSTSP和PDSTSP之间的差异，为物流配送方面的优化提供更有效的解决方案。
6.FSTSP的实现代码
目前全网还没有关于这篇文章中对于FSTSP的开源代码，于是小编在这里为大家准备了复现本文FSTSP和PDSTSP的python实现源码，包括FSTSP的Gurobi建模求解及其启发式算法，PDSTSP的Gurobi建模求解。部分代码内容展示如下。全部内容已经上传至云盘，欢迎大家关注“运筹OR帷幄”公众号！
代码链接: https://pan.baidu.com/s/1VMjK0ZGSomQSZy4_P4SnMw 提取码: rc4q
也欢迎关注小编B站个人账号^_^：https://space.bilibili.com/168828224
def modelingAndSolve(data):
# 建立模型
m = Model('FSTSP')
# 模型设置：由于存在函数printSolution，因此关闭输出;以及容许误差
m.setParam('MIPGap', 0.01)
# m.setParam('OutputFlag', 0)
# 定义变量：
# Step1.建立合适的数据结构建立变量的索引，FSTSP有两个决策变量X和Y,一个MTZ辅助变量U，一个卡车时间辅助变量T，一个无人机时间辅助变量T_UAV，一个访问次序一致变量p
'______添加决策变量X______'
# 建立存储决策变量X的数据结构
X = [[[] for _ in range(0, data.nodeNum)] for _ in range(0, data.nodeNum)] # x_ij
# 根据数据结构向模型中添加对应下标的变量
for i in data.N_out_lst:
for j in data.N_in_lst:
if i != j:
X[i][j] = m.addVar(vtype=GRB.BINARY, name=f"X_{i}_{j}")
'______添加决策变量Y______'
# 建立存储决策变量Y的数据结构
Y = [[[[] for _ in range(0, data.nodeNum)] for _ in range(0,data.nodeNum)] for _ in range(0,data.nodeNum)] # Y_ijk
# 根据数据结构向模型中添加对应下标的变量
for i in data.N_out_lst:
for j in data.C_lst:
if i != j:
for k in data.N_in_lst:
if (i,j,k) in data.P_lst:
Y[i][j][k] = m.addVar(vtype=GRB.BINARY, name=f"Y_{i}_{j}_{k}")
'______添加MTZ辅助变量U______'
# 建立存储辅助变量U的数据结构
U = [[] for _ in range(0, data.nodeNum)] # U_i
# 根据数据结构向模型中添加对应下标的变量
for i in data.N_in_lst:
U[i] = m.addVar(vtype=GRB.CONTINUOUS,lb=1.0, ub=data.customerNum+2, name=f"U_{i}")
'______添加卡车时间辅助变量T______'
# 建立存储辅助变量T的数据结构
T = [[] for _ in range(0, data.nodeNum)] # T_i
# 根据数据结构向模型中添加对应下标的变量
for i in data.N_lst:
T[i] = m.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name=f"T_{i}")
'______添加无人机时间辅助变量T_UAV______'
# 建立存储辅助变量TUAV的数据结构
TUAV = [[] for _ in range(0, data.nodeNum)] # TUAV_i
# 根据数据结构向模型中添加对应下标的变量
for i in data.N_lst:
TUAV[i] = m.addVar(vtype=GRB.CONTINUOUS, lb=0.0, name=f"TUAV_{i}")
'______添加访问次序一致变量P______'
# 建立存储辅助变量P的数据结构
P = [[[] for _ in range(0, data.nodeNum)] for _ in range(0, data.nodeNum)] # P_ij
# 根据数据结构向模型中添加对应下标的变量
for i in data.N_out_lst:
for j in data.C_lst:
if j != i:
P[i][j] = m.addVar(vtype=GRB.BINARY, name=f"P_{i}_{j}")
m.update()
# var_lst = m.getVars()
# print(var_lst)
# 定义目标函数
obj = LinExpr(0) # 线性项，初始值为0，可用.addTerms(a,x) 意为将变量x增加到表达式，系数为a
obj.addTerms(1, T[data.customerNum+1])
m.setObjective(obj, sense=GRB.MINIMIZE)
# 定义约束条件:
# 1.客户点访问约束
num = 0
for j in data.C_lst:
expr = LinExpr(0)
for i in data.N_out_lst:
if i != j:
expr.addTerms(1, X[i][j])
for i in data.N_out_lst:
if i != j:
for k in data.N_in_lst:
if (i,j,k) in data.P_lst:
expr.addTerms(1, Y[i][j][k])
num += 1
m.addConstr(expr == 1, f'C1_{num}')
# 2.卡车起点车场流出约束
num = 0
expr = LinExpr(0)
for j in data.N_in_lst:
expr.addTerms(1, X[0][j])
num += 1
m.addConstr(expr == 1, f'C2_{num}')
# 3.卡车终点车场流入约束
num = 0
expr = LinExpr(0)
for i in data.N_out_lst:
expr.addTerms(1, X[i][data.customerNum+1])
num += 1
m.addConstr(expr == 1, f'C3_{num}')
# 4.卡车在客户点的流平衡约束
num = 0
for j in data.C_lst:
expr = LinExpr(0)
for i in data.N_out_lst:
if i != j:
expr.addTerms(1,X[i][j])
for k in data.N_in_lst:
if k != j:
expr.addTerms(-1,X[j][k])
num += 1
m.addConstr(expr == 0, f'C4_{num}')
# 5.卡车破子圈约束
num = 0
for i in data.C_lst:
for j in data.N_in_lst:
if i != j:
expr = LinExpr(0)
expr.addTerms(1,U[i])
expr.addTerms(-1,U[j])
expr.addTerms(data.customerNum+2, X[i][j])
num += 1
m.addConstr(expr <= data.customerNum+1, f'C5_{num}')
# 6.无人机非终点车场发射约束
num = 0
for i in data.N_out_lst:
expr = LinExpr(0)
for j in data.C_lst:
if j != i:
for k in data.N_in_lst:
if (i, j, k) in data.P_lst:
expr.addTerms(1, Y[i][j][k])
num += 1
m.addConstr(expr <= 1,f'C6_{num}')
# 7.无人机非起点车场回收约束
num = 0
for k in data.N_in_lst:
expr = LinExpr(0)
for i in data.N_out_lst:
if i != k:
for j in data.C_lst:
if (i,j,k) in data.P_lst:
expr.addTerms(1,Y[i][j][k])
num += 1
m.addConstr(expr <= 1, f"C7_{num}")
# 8.非起点架次下，无人机发射/回收节点的卡车访问约束
num = 0
for i in data.C_lst:
for j in data.C_lst:
if j != i:
for k in data.N_in_lst:
if (i,j,k) in data.P_lst:
expr = LinExpr(0)
for h in data.N_out_lst:
if h != i:
expr.addTerms(1, X[h][i])
for l in data.C_lst:
if l != k:
expr.addTerms(1, X[l][k])
expr.addTerms(-2,Y[i][j][k])
num += 1
m.addConstr(expr >= 0, f'C8_{num}')
# 9.起点架次下，无人机回收节点的卡车访问约束
num = 0
for j in data.C_lst:
for k in data.N_in_lst:
if (0,j,k) in data.P_lst:
expr = LinExpr(0)
for h in data.N_out_lst:
if h != k:
expr.addTerms(1,X[h][k])
expr.addTerms(-1,Y[0][j][k])
num += 1
m.addConstr(expr >= 0, f'C9_{num}')
# 10.无人机架次发射与回收节点卡车访问顺序约束
num = 0
for i in data.C_lst:
for k in data.N_in_lst:
if k != i:
expr = LinExpr(0)
expr.addTerms(1, U[i])
expr.addTerms(-1, U[k])
for j in data.C_lst:
if (i,j,k) in data.P_lst:
expr.addTerms(data.customerNum+2,Y[i][j][k])
num += 1
m.addConstr(expr <= data.customerNum+1, f'C10_{num}')
# 11.无人机发射节点的左侧时间窗约束
num = 0
for i in data.C_lst:
expr = LinExpr(0)
expr.addTerms(1,T[i])
expr.addTerms(-1,TUAV[i])
for j in data.C_lst:
if j != i:
for k in data.N_in_lst:
if (i,j,k) in data.P_lst:
expr.addTerms(M, Y[i][j][k])
num += 1
m.addConstr(expr <= M, f'C11_{num}')
# 12.无人机发射节点的右侧时间窗约束
num = 0
for i in data.C_lst:
expr = LinExpr(0)
expr.addTerms(1,TUAV[i])
expr.addTerms(-1,T[i])
for j in data.C_lst:
if j != i:
for k in data.N_in_lst:
if (i, j, k) in data.P_lst:
expr.addTerms(M, Y[i][j][k])
num += 1
m.addConstr(expr <= M, f'C12_{num}')
# 13.无人机回收节点的左侧时间窗约束
num = 0
for k in data.N_in_lst:
expr = LinExpr(0)
expr.addTerms(1, T[k])
expr.addTerms(-1, TUAV[k])
for i in data.N_out_lst:
if i != k:
for j in data.C_lst:
if (i, j, k) in data.P_lst:
expr.addTerms(M, Y[i][j][k])
num += 1
m.addConstr(expr <= M, f'C13_{num}')
# 14.无人机回收节点的右侧时间窗约束
num = 0
for k in data.N_in_lst:
expr = LinExpr(0)
expr.addTerms(1, TUAV[k])
expr.addTerms(-1, T[k])
for i in data.N_out_lst:
if i != k:
for j in data.C_lst:
if (i, j, k) in data.P_lst:
expr.addTerms(M, Y[i][j][k])
num += 1
m.addConstr(expr <= M, f'C14_{num}')
# 15.卡车有效到达时间约束
num = 0
for h in data.N_out_lst:
for k in data.N_in_lst:
if h != k:
expr = LinExpr(0)
expr.addTerms(1, T[h])
expr.addTerms(-1, T[k])
for l in data.C_lst:
if l != k:
for m_ in data.N_in_lst:
if (k,l,m_) in data.P_lst:
expr.addTerms(SL, Y[k][l][m_])
for i in data.N_out_lst:
if i != k :
for j in data.C_lst:
if (i,j,k) in data.P_lst:
expr.addTerms(SR,Y[i][j][k])
expr.addTerms(M,X[h][k])
num += 1
m.addConstr(expr + data.costMatrix[h][k] <= M, f'C15_{num}')# 34
# 16.无人机交付点j的时间约束
num = 0
for j in data.C_UAV_lst:
for i in data.N_out_lst:
if i != j:
expr = LinExpr(0)
expr.addTerms(1,TUAV[i])
expr.addTerms(-1,TUAV[j])
for k in data.N_in_lst:
if (i,j,k) in data.P_lst:
expr.addTerms(M,Y[i][j][k])
num += 1
m.addConstr(expr + UAV_factor * data.costMatrix[i][j] <= M, f'C16_{num}')
# 17.无人机回收点k的时间约束
num = 0
for j in data.C_UAV_lst:
for k in data.N_in_lst:
if k != j:
expr = LinExpr(0)
expr.addTerms(1,TUAV[j])
expr.addTerms(-1, TUAV[k])
for i in data.N_out_lst:
if (i,j,k) in data.P_lst:
expr.addTerms(M,Y[i][j][k])
num += 1
m.addConstr(expr + UAV_factor * data.costMatrix[j][k] + SR <= M, f'C17_{num}')
# 18.无人机架次电量约束
num = 0
for k in data.N_in_lst:
for j in data.C_lst:
if j != k:
for i in data.N_out_lst:
if (i,j,k) in data.P_lst:
expr = LinExpr(0)
expr.addTerms(1,TUAV[k])
expr.addTerms(-1, TUAV[j])
expr.addTerms(M,Y[i][j][k])
num += 1
m.addConstr(expr + UAV_factor * data.costMatrix[i][j] <= UAV_endurance + M, f'C18_{num}')
# 19.客户访问先后辅助变量的定义与建立约束1
num = 0
for i in data.C_lst:
for j in data.C_lst:
if j != i:
expr = LinExpr(0)
expr.addTerms(1, U[i])
expr.addTerms(-1, U[j])
expr.addTerms(data.customerNum+2, P[i][j])
num += 1
m.addConstr(expr >= 1, f'C19_{num}')
# 20.客户访问先后辅助变量的定义与建立约束2
num = 0
for i in data.C_lst:
for j in data.C_lst:
if j != i:
expr = LinExpr(0)
expr.addTerms(1, U[i])
expr.addTerms(-1, U[j])
expr.addTerms(data.customerNum + 2, P[i][j])
num += 1
m.addConstr(expr <= data.customerNum + 1, f'C20_{num}')
# 21.客户访问先后辅助变量的定义与建立约束3
num = 0
for i in data.C_lst:
for j in data.C_lst:
if j != i:
expr = LinExpr(0)
expr.addTerms(1, P[i][j])
expr.addTerms(1, P[j][i])
num += 1
m.addConstr(expr == 1, f'C21_{num}')
# 22.无人机架次的先后时间约束/不相交架次下先后时间约束
num = 0
for i in data.N_out_lst:
for k in data.N_in_lst:
if i != k:
for l in data.C_lst:
if (l != i) and (l != k):
expr = LinExpr(0)
expr.addTerms(1, TUAV[k])
expr.addTerms(-1, TUAV[l])
expr.addTerms(M, P[i][l])
for j in data.C_lst:
if ((i, j, k) in data.P_lst) and (j != l):
expr.addTerms(M, Y[i][j][k])
for m_ in data.C_lst:
if (m_ != i) and (m_ != k) and (m_ != l):
for n in data.N_in_lst:
if ((l, m_, n) in data.P_lst) and (n != i) and (n != k):
expr.addTerms(M, Y[l][m_][n])
num += 1
m.addConstr(expr <= 3*M, f'C22_{num}')
# 23~.辅助变量符号约束
m.addConstr(T[0] == 0, f'C23')
m.addConstr(TUAV[0] == 0, f'C24')
num = 0
for j in data.C_lst:
num += 1
m.addConstr(P[0][j] == 1, f'C25_{num}')
# 记录求解开始时间
start_time = time.time()
# 求解
m.optimize()
m.write('FSTSP.lp')
if m.status == GRB.OPTIMAL:
print("-" * 20, "求解成功", '-' * 20)
# 输出求解总用时
print(f"求解时间: {time.time() - start_time} s")
print(f"目标函数为（卡车返回车场的时间）: {m.ObjVal}")
solution = getSolution(data,m)
plotSolution(data, solution)
printSolution(data, solution)
else:
print("无解")
return m
7.结论
本论文探讨了无人机（UAV）在物流配送中最后一公里的崭新角色，开创性地提出了车辆与无人机联合路由优化的研究框架。文章针对两种配送场景，提出了FSTSP（串联的卡车无人机配送模式）和PDSTSP（并行的卡车无人机配送模式）两种问题，并给出了相应的混合整数线性规划模型，以优化无人机与传统配送卡车的路线和调度。同时，设计了简单有效的启发式算法，以适应实际应用需求。通过数值分析，论文强调了无人机在飞行速度与续航时间之间的权衡，指出这种联合配送模式能够显著提高客户订单的交付速度，降低分销成本，并减少环境影响。总体而言，本文旨在解决无人机在物流应用中的操作挑战，推动其更广泛的采用与效率的提升。
参考文献：
Carlsson J G, Song S. Coordinated logistics with a truck and a drone[J]. Management Science, 2018, 64(9): 4052-4069.
MirHassani S A, Ebrazi R. A flexible reformulation of the refueling station location problem[J]. Transportation Science, 2013, 47(4): 617-628.
Agatz N, Bouman P, Schmidt M. Optimization approaches for the traveling salesman problem with drone[J]. Transportation Science, 2018, 52(4): 965-981.
Drexl M. Synchronization in vehicle routing — a survey of VRPs with multiple synchronization constraints[J]. Transportation Science, 2012, 46(3): 297-316.
微信公众号后台回复
加群：加入全球华人OR|AI|DS社区硕博微信学术群
资料：免费获得大量运筹学相关学习资料
人才库：加入运筹精英人才库，获得独家职位推荐
电子书：免费获取平台小编独家创作的优化理论、运筹实践和数据科学电子书，持续更新中ing...
加入我们：加入「运筹OR帷幄」，参与内容创作平台运营
知识星球：加入「运筹OR帷幄」数据算法社区，免费参与每周「领读计划」、「行业inTalk」、「OR会客厅」等直播活动，与数百位签约大V进行在线交流
文章须知
推文作者：张晨皓，东南大学在读博士生
责任编辑：江镕行
微信编辑：疑疑
文章由『运筹OR帷幄』原创发布
如需转载请在公众号后台获取转载须知
推荐阅读 ：
乘积包络谱优化- Gram：一种用于滚动轴承故障诊断的增强包络分析
基于混合注意力的多小波系数融合滚动轴承剩余寿命预测方法
ReF-DDPM: 一种基于DDPM的滚动轴承故障诊断数据增强新方法
