接近最优的考虑转弯代价的覆盖路径规划

Krupke D M. Near-Optimal Coverage Path Planning with Turn Costs[J]. arXiv preprint arXiv:2310.20340, 2023.

摘要

覆盖路径规划是机器人技术中的一个基本挑战，它在航空监视、制造、清洁、检查、农业等多个领域都有广泛的应用。其主要目标是为代理设计一条轨迹，以有效地覆盖给定区域，同时最小化时间或能源消耗。现有的实用方法往往缺乏坚实的理论基础，依赖于纯粹的启发式方法，或者将问题过度抽象为简单的网格图上的旅行商问题（TSP）。此外，考虑的成本函数很少考虑转弯成本，不平均覆盖需求的奖励收集变体，或任意几何区域。

在本文中，我们描述了一系列系统方法，用于处理源自复杂多边形环境的任意网格。这种适应为计算高效覆盖路径铺平了道路，并为现实世界的机器人应用提供了坚实的理论基础。通过全面评估，我们证明该算法还表现出低最优性差距，同时有效地处理复杂环境。此外，我们展示了它在处理部分覆盖和适应异质通行成本方面的多功能性，提供了权衡覆盖质量和时间效率的灵活性。

1 引言

覆盖路径规划对于多种应用领域都是一个重要问题，比如航空监视[14]、清洁[13]、铣削[37]、割草[30]、病虫害控制[9]等等。这个问题已经获得了相当多的关注，主要是从实践的角度出发，但也涉及到一些理论成果。该问题在多个层面上都难以解决，因为它包含了NP难题和PSPACE难题，例如旅行商问题（TSP）、覆盖问题和钢琴搬运者问题。

问题的最简单理论抽象是网格图上的TSP。在这里，我们简单地在区域上放置一个网格，其单元格大小与代理的覆盖能力相匹配，然后计算其上的最短游程。因为TSP在许多应用中都会出现，它是最为广泛研究的优化问题之一，尽管它的难度已经被证明，但仍有能力很强的求解器。Concorde求解器[3]能够求解数万个顶点的实例并证明其最优性[5]，也有其他算法可以为更大的实例计算出好的解。例如，Bormann等人[13]使用Concorde来优化覆盖路径。

尽管在网格图上解决TSP旨在最小化游程长度，这是能源消耗中的一个重要因素，但这种狭隘的优化标准可能导致意想不到的后果。在多旋翼飞行器等应用中，更直的飞行路径通常更加节能[15, 39]。单纯只关注最小化覆盖游程长度的目标通常会鼓励曲折的路线，因为这种方法能够，例如，一次性覆盖两个车道。因此，这些表面上较短的游程实际上可能执行起来更加昂贵。

这个问题在带转弯成本的铣削问题（Milling with Turn Costs）中得到了解决，该问题不仅最小化路径长度，还最小化了路径在网格中执行的转弯角度总和[7]。尽管还没有捕捉到所有动态，但它作为各种情景下更现实的近似，缓解了仅关注长度最小化的不足。不幸的是，转弯成本提高了问题的复杂性，以至于不仅是本身的问题，就连循环覆盖放松也成为NP难题[24]。尽管最优可解决的问题规模从不到100个顶点[20]增加到超过1000个顶点[25]，但与经典TSP的仍然很大差距显示了计算实际动态模型的最优解的限制，即使是对于高度简化的环境。

除了复杂的动态之外，我们有时并不需要覆盖整个区域。在许多情况下，甚至无法实现真正的100%覆盖，因为工具根本无法适应每一个角落。相反，我们有一个可行的区域可以移动，以及其内一个实际上是“有价值”的更小子集。真空机器人可以在整个房间内移动，但通常有易于积垢的区域和较干净的区域，这些区域不需要每次都清洁。收割机可以沿着整个田间移动，但是作物产量可能是异质的；收割机不需要收割所有东西，而只是大部分产量。对于航空监督而言，有些区域的兴趣度比其他区域更高。此外，可能有些区域比其他区域更难通过，例如对无人机（UAV）的风场[53]，以及对地面车辆的复杂地形或倾斜[30]。

Fekete和Krupke[24, 25]为网格图上的带转弯成本的铣削问题提出了一个常数因子近似算法，该算法也能够通过跳过惩罚来处理部分覆盖。在本文中，我们将这一算法推广到从多边形环境和异质成本获得的任意网格上，这使我们能够基于理论基础计算出高效的轨迹，以适用于现实世界的应用，见图1。我们在评估中展示了算法能够计算出与最优解平均接近（10%至15%）的解，就网格表示而言。尽管对于任意网格可能失去了常数因子近似保证，本文展示了如何将一个针对正方形网格的理论算法推广到现实世界的应用中。

1.1 相关工作

规划一条工具的轨迹以覆盖一个区域，例如割草或吸尘，被称为覆盖路径规划问题（CPP）。CPP已经在不同的应用、模型（例如多机器人）、约束和目标中得到了广泛的关注，这在多个综述文章中有所体现[17, 27, 13, 14]。有多种方法，两种最突出的是：（1）将较大的区域分解成可以使用螺旋或之字形模式覆盖的简单区域（[42, 19, 18]），以及（2）在区域上应用（规则的正方形）网格，每个网格单元大致代表覆盖区域，将几何覆盖问题转换为网格图上的离散巡回问题（[12, 40, 54, 50, 39]）。在这篇文章中，我们使用了第二种方法，但是推广到可以更好地适应区域的任意网格，因为合适的网格可以极大地改善可实现的路径。当仅考虑轨迹长度时，问题变成了著名的旅行推销员问题（TSP），即使在正方形网格中也是NP难的[32]，但由于算法工程学的广泛应用，在实践中可以很好地解决。为了解决显著的动态问题，我们需要考虑转弯成本，这使得问题变得更加困难。即使以前简单的放松也变成了NP难题[24]，但是有常数因子近似解可用[6, 7, 24]。在网格图上，大约1000个顶点的实例可以求得最优解，而近似算法已经应用到多达300000个顶点的实例[25]。对于平面中的一般点，问题被称为角度度量TSP，目前只有对数级的近似解[1]。将其进一步推广到抽象图形是二次TSP问题，例如，在生物信息学中扮演重要角色[26]。这些问题中，只有不到100个顶点的实例可以预期在合理时间内解决到最优[33, 45, 2]。在实际方面，CPP在距离和转弯成本方面的不同程度模型上已被考虑，例如仅最小化转弯次数[34]，转弯角度总和[12, 43, 39]（像这篇文章一样），甚至基于模型和实验的成本函数[42, 15]。CPP中也考虑了包含异质成本函数，例如[54, 30]，以及简单路径规划[38, 52, 46]。

本文的另一个方面是能够基于某些价值分布选择性地覆盖区域。有一些论文也考虑了部分覆盖路径规划。Papachristos等人[43]以及Ellefsen、Lepikson和Albiez[23]考虑了具有距离和转弯成本的三维结构的部分检查。Jensen等人[34]和Soltero等人[51]执行没有固定半径的覆盖，但最小化（加权）兴趣点到轨迹的距离。Murtaza等人[40]计算了对区域的全覆盖，但根据概率分布优先处理子区域，以快速找到目标。Sharma等人[50]也计算了区域的全覆盖，但由于预算有限，结果是多次巡回，尝试尽可能有效地覆盖。然而，所有这些问题都与我们的问题有显著的不同。在理论方面，有允许以罚金跳过顶点的Penalty TSP和尝试在预算范围内尽可能多覆盖的Budget TSP。Ausiello等人[8]提供了此类问题的概述。

在本文中，我们做出以下贡献：

• 我们通过使用网格算法，将常规网格图中的覆盖巡回近似算法推广到更真实的多边形实例，为现实世界的应用计算出更有效的覆盖巡回路径，并以稳健的理论基础为基石。

• 我们通过使用基于行进距离和转角总和的线性组合模型，以及异质路径成本的局部乘法因子，来近似代理的动态特性。这也允许创建软障碍，应该（但不必）避免。

• 我们通过对未覆盖区域使用罚金来研究部分覆盖，这允许在覆盖质量和时间效率之间进行权衡。可以加权区域，以目标重要区域进行巡回。

• 我们通过使用大邻域搜索（LNS），局部改进巡回路径，能够将巡回路径改进几个百分点。

• 我们在超过500个实例上评估了实施的最优性差距，这些实例是半自动生成的，以模仿现实世界的场景。提供数据和代码。

我们没有保持原算法的近似因子，但我们展示了该实现仍能够利用可靠的下界，在任意网格上计算出良好的解决方案。由于缺乏真实世界的实例和模型，评估仅在合成实例上进行，这些实例是半自动生成的，试图模仿农业区域、多建筑位置和复杂建筑结构。没有执行与不受网格限制的几何模型的比较，因为难以获得强有力的下界。然而，对不同网格和网格的可行解决方案质量进行了比较，并且在评估中使用了最佳的网格策略。我们注意到，在选择网格分辨率时，关注边缘的覆盖而不是点的覆盖能够改善结果。与方形网格相比，使用六边形网格也显示出好处，尤其是在转弯成本较高的情况下。此外，重要的是要注意，并非所有的网格算法都适合解决我们特定的问题。相关研究已附在附录C中。

1.3 预备知识

给定一个图 $G = (P, E)$，其中$P \subset \mathbb{R}^2$是一组路径点，它们跨越了一个潜在的轨迹，而$E$是连接两个路径点的线段。另外，我们给出了一个价值函数$\text{val} : P \rightarrow \mathbb{R}_+$，它为每个路径点分配一个价值，和一个成本函数$\text{cost} : P^3 \rightarrow \mathbb{R}_+$，它为每三个连续的路径点$u, v, w$分配一个成本，其中$uv, vw \in E$。我们称这样的三元组为通过中间点$v$的通道。目标是找到一个路径$T=p_0,p_1,\cdots,p_{|T|-1},p_0$其中$p_ip_{i+1}\in E$对于所有$i\in\{0,\cdots,|T|-1\}$成立，使得最小化如下目标函数

$$\begin{equation} \min_{T} \underbrace {\sum_{i=0}^{|T|-1} \text{cost}(p_{i-1}, p_i, p_{i+1})}_{\text{Touring cost}} + \underbrace{\sum_{p \in P, p \notin T} \text{val}(p)}_{\text{Coverage loss}} \end{equation}$$

我们定义通过路径$uvw$的成本为两线段长度及其转角的线性组合，由$\tau \in \mathbb{R}^+$加权。它还可能乘以一个局部因子$\alpha_v$进行缩放。

$$\text{cost}(u, v, w) = \alpha_v \left( \frac{d(u, v) + d(v, w)}{2} \right) + \tau \cdot \text{turn}(u, v, w)$$

距离减半是为了避免对边重复计算成本。更深入的讨论请参见附录A。

2 泛化算法

在本节中，我们展示了如何调整Fekete和Krupke的算法 [24, 25] 来解决包括昂贵区域和有价值区域的多边形实例。更准确地说，我们展示了如何使用嵌入图来近似区域，如何使原有算法适应任意嵌入图，并且加入优化。泛化算法有七个步骤：(1) 将多边形实例转换为离散图的路径点。 (2) 在这个图中使用线性规划计算分数解。 (3) 使用分数解来选择原子条带。对于一般网格而言，这一步比对于正方形网格更加复杂。 (4) 在原子条带上进行匹配并获得循环覆盖。 (5) 改进循环覆盖。 (6) 连接循环以形成路径。 (7) 改进路径。步骤 (1)、(3)、(5) 和 (7) 与原始算法有显著不同，我们将详细描述它们。然而，如果我们在步骤 (1) 中给定一个规则的正方形网格，并且禁用局部优化步骤 (5) 和 (7)，那么算法的行为几乎与原始算法相同。结果轨迹如图2所示。

2.1 步骤1: 离散化

我们应用网格生成算法 dmsh (v0.2.17)[49] 并通过optimesh [48] 进行额外平滑处理，对多边形进行处理以获得合适的网格，如图3a所示。两个路径点之间的最优距离设定为$0.95·4\sqrt{3}·r$，其中$r$代表覆盖半径，假设工具具有圆形覆盖区域。在示例中我们设定$r=1$，但算法适用于任何$r$。在三角形网格中，距离$4/\sqrt{3}·r\approx3.31·r$导致平行线之间的距离刚好是$2·r$。由于dmsh倾向于将顶点设置得过于分散而非过于靠近，我们通过将距离减少5%来对此进行抵消。这种稀疏网格上的路径会在转弯时错过一些区域，但我们通过最小化转弯次数，并在后处理中略微放大转弯来补偿错过的覆盖区域。覆盖值是通过路径点的Voronoi单元覆盖的区域来估计的，见图3b。我们也可以使用代理在航点处的覆盖范围，但这种方法不太准确，因为覆盖主要是在沿边缘移动时发生的。获得能产生良好巡视路径的网格并非易事，需要进行大量的实验来找到一个好的网格生成算法和参数。离散化的许多方面也都在与几何操作的臭名昭著的数值问题作斗争，必须谨慎处理。也可以使用 gmsh [28] 的 Parallelograms Packing 算法来获得类似质量好的网格。gmsh 更快也更稳定，但与 dmsh 相比，质量的离群值更多。还有许多其他的网格生成算法，但它们大多数不适合我们的目的，因为它们不允许平滑的轨迹和大小相等的单元，而是专注于不同的质量。更多细节可以在附录C中找到。

2.2 第二步：线性松弛

给定图$G = (P, E)$，我们可以通过使用线性规划获得循环覆盖的分数解。我们处理覆盖航点$v\in P$的通道$uvw = wvu$，这些航点来自或指向邻近航点$u, w\in N(v)$。对于每个通道$uvw$，变量$xuvw≥0$ 表示该通道被使用的频率。此外，我们使用变量$s_v≥0$表示跳过该航点并支付其覆盖损失的费用。

$$\begin{equation} \min \sum_{v \in P} \text{val}(v) \cdot s_{v} + \sum_{u,w \in N(v)} \text{cost}(u, v, w) \cdot x_{uvw} \end{equation}$$ $$\begin{equation} \text{s.t.} \sum_{u,w \in N(v)} x_{uvw} + s_{v} \geq 1 \quad \forall v \in P \end{equation}$$ $$\begin{equation} 2 \cdot x_{uvw} + \sum_{u \in N(v),u \neq w} x_{wvu} = 2 \cdot x_{uvw} + \sum_{u \in N(w),u \neq v} x_{uvw} \quad \forall vw \in E \end{equation}$$

方程（2.2）的目标很简单，就是最小化漏覆盖的价值和巡回成本。方程（2.3）强制每个航点要么被覆盖要么被跳过，方程（2.4）确保流的一致性，即每条边应该从两边被等量使用。关于覆盖全部区域的分数解和部分覆盖的示例，请参见图4a和图4b。

2.3 第三步：原子条带

在下一步中，我们希望使用前一解的分数解作为提示，计算一个循环覆盖。如果成本仅依赖于距离，那么可以通过最小权重完美匹配有效地计算循环覆盖。为此，我们会将每个航点替换为两个顶点，并将它们与所有其他顶点通过对应距离连接起来，这个距离可以通过Dijkstra算法有效计算。为了实现部分覆盖，我们会在一个航点的两个顶点之间添加一条对应覆盖损失价值的边。最小权重完美匹配将会强制每个航点都有一条进入和一条离开的轨迹，即在一个循环中，或者只使用内部边并跳过航点。

考虑到转弯成本，循环覆盖问题变得 NP 难，但 Fekete 和 Krupke [24] 表明，我们可以使用前一步的分数解来估计我们通过航点的方向，并将相应的转弯成本移到边权重上。在方格网格中，可以证明这种技术会产生4近似解，在三角形网格中则为6近似解。这可以想象为将每个航点替换为一个$\epsilon$长度的线段，如图5所示，其方向在分数解中使用最多。我们称这些$\epsilon$长度的线段为原子条带（atomic strips）。在端点上计算最小权重完美匹配，将得到包括所有这些线段的最优循环覆盖。如果这些线段被正确选择（这是 NP 难的），那么最小权重完美匹配实际上对应于航点的最优循环覆盖。

网格使得选择这些原子条带更加复杂，因为可能有不止两个或三个合理的方向。原子条带的一个有用特性是，转弯越大，有多种方向是最优的。对于一个 U 形转弯，每个方向都是最优的。通道越笔直，好的方向就变得越重要；但通常这些情况可以从分数解中很容易判断出来。因此，将潜在方向限制为与之相邻的边的方向是明智的。我们根据分数解的通道与其拟合程度对每个方向进行加权，选择具有最高总和的方向。

将所有航点彼此连接会导致二次数量的边，其权重难以计算。Fekete 和 Krupke [25] 指出，更有效的做法是仅将航点与其邻居连接（这样权重也容易计算），并允许可选的原子条带处理潜在的必要重叠轨迹。可选的原子条带可以通过在其端点之间添加权重为零的边来实现，允许它在不增加额外成本的情况下被中和。Arkin 等人 [7] 表明，在方格网格中，每个顶点最多被访问四次，从而限制了必要可选原子条带的数量。对于三角形网格，必要的访问次数可以是线性的，如图6所示，这会破坏使用此优化时的近似因子。然而，这是一个人工构造的情况，在我们的情况中，每个航点通常只被覆盖一次或两次。另一个挑战是，可选的原子条带还必须与更长边的原始轨迹匹配，以重建实际成本。否则，通过可选的原子条带连接两个航点可能比直接连接它们更昂贵。为了解决这个问题，可以为任何方向添加一定数量的可选条带，但这也会增加计算复杂性。因此，我们将每个航点的原子条带数量限制为常数$k$，并且允许每个航点$p\in P$ 最多有一个邻居 $n\in N(p)$的原子条带。这将使辅助图的复杂度保持在$O(|P| ·k2)$。不同$k$的示例可以在图7中看到，详细的实现描述在附录B中。

2.4 第四步：匹配

我们留下了一个带有原子条带端点的加权图，我们想要计算一个最小匹配。在网格中，任何相邻航点的原子条带端点之间都有边。权重对应于两个航点之间的巡回成本，并且对应于端点的方向。此外，每个原子条带都有一条连接其两个端点的边。对于强制原子条带，权重对应于机会损失，即当不覆盖它时的分配覆盖值。对于所有其他情况，成本为零，以允许跳过它们而不产生额外成本。设$k$为航点上的最大原子条带数量，那么匹配实例中的顶点和边的数量在$O(|P| ·k2)$范围内。我们使用 Kolmogorov 的 Blossom V 算法 [35] 解决相应的最小权重完美匹配实例。作者声称最坏情况的复杂度是$O(n^3m)$，这可能是不切实际的，但在实践中，即使对于大型实例，它也显示出足够快的速度。通过匹配的端点连接原子条带，得到一组循环，见图9a，我们可以在第六步中连接它们。

2.5 第五步：局部优化

在我们继续将循环连接成单一巡回路径之前，我们可以优化循环覆盖。为此，我们选择解的一个小但昂贵部分，并通过混合整数规划计算（几乎）最优解。这可以重复多次，直到获得令人满意的解决方案，见图8。请注意，可以如[41, 25]所述，在常规方格网格中解决包含1000个顶点的许多实例达到最优解。此外，对于不规则网格，通常可以在几秒钟内解决少于100个顶点的小实例。我们将局部优化所需的顶点数表示为$t$。通过选择一个昂贵的根和选择广度优先搜索的前$t$个顶点，我们选择要优化的昂贵区域。解中一个航点的成本被表示为覆盖它的通道的成本，或者如果它未被使用，则表示为相应的机会损失。为了使选择更加健壮，我们还通过对所有直接邻居的开销进行求和来包括它们的开销。

通过简单地将分数变量替换为整数变量，2.2节中的线性规划问题就变成了相应的混合整数规划（MIP）问题。在这个MIP中，我们固定了给定解的所有变量，除了与$t+1$个选定的航点对应的变量。当然，我们根本不需要在这个MIP中包括已经固定的航点，只需将相应的常数放入方程(2.4)中即可。这确保了局部解与固定的外部解保持一致。在优化局部MIP后，我们替换解中的这部分内容，并排除了根及其邻居作为后续迭代的根选择。这是因为昂贵的部分可能已经在其局部区域内是最优的，不应再次进行优化。MIP的一个有用特性是，如果我们的（局部）解已经是（几乎）最优的，那么优化过程通常会更快。如果我们将相应的起始解提供给MIP求解器，它只需要找到一个匹配的下界。通过使用运行时间和实际的改进，可以改进下一个区域的选择，或者动态增加它。通过选择不相交的区域，这种优化方法还允许有效的并行化。然而，我们将这些优化留给未来的工作，只是针对固定的区域大小$t$执行$i$次迭代。

2.6 第六步：连接循环

现在，我们只需要连接这些循环以形成一条巡回路径。对于相邻的循环，这相当简单，只涉及最小的额外成本：只需穿过连接两个循环的每条边，并通过最便宜的边进行合并，见图9b。一个简单的优化是使用两条平行边一次，而不是一条边两次，但这在2.7节中也会自动完成。如果这些循环之间相隔较远，情况就会变得更加复杂。连接成本可能实际上超过了相应循环的巡回成本。如果循环所覆盖的区域价值不够高，我们最好只是删除该循环，见图10。

要选择任何循环，首先需要知道每个循环的价值。我们通过覆盖的航点价值之和来估算一个循环的价值。如果一个航点出现在不同的循环中，只有第一个循环获得它的价值。这可能发生在两个循环交叉并且由于转弯成本无法连接的情况。由于这种情况很少发生，所以如果航点的价值准确的话，估算的循环价值就是准确的。否则，循环的价值可能被低估，导致解的质量稍低。接下来，我们需要知道连接任意两个循环的成本。这可以通过在边图上的一种迪杰斯特拉变种来实现。在网格的边图上工作使我们能够包括不仅路径的距离，还有任意两个边之间的转弯成本。为了简化事情，我们使用了一个有向版本，其中还包括了我们通过边的方向，如图11所示。

距离和转弯成本分配给了蓝色的弧线。现在，使用边图中的出边可以简单地分配边的距离成本。如果我们将$k$设为网格中的最大度数，则辅助图中最多有$O(|P| ·k)$个顶点和$O(|P| ·k2)$条边。使用迪杰斯特拉算法，我们可以在$O(|P| ·k2log|P|)$的时间内计算出任意两个边之间的最便宜路径（忽略可能收集到的覆盖价值）。成本是对称的，所以在两个方向上都是最优的。现在仍然缺少的是合并（加倍）路径与循环的成本。对于与两个循环相邻的边的所有组合来说，检查它们的成本将会很昂贵。相反，我们可以选择其中一个循环，并将所有与之相邻的边初始化为与它的最终连接成本。现在，我们只需要使用已经计算的迪杰斯特拉算法的距离来找到到目标循环的最便宜的相邻边。有了这两个信息，我们可以在循环和它们的连接上计算一个奖励收集斯坦纳树（PCST）。得到的树对应于有价值的循环以及如何连接它们。计算最优的PCST是NP难的，但这里获得的循环覆盖通常足够小，可以使用整数规划来最优地解决。否则，可以使用Goemans和Williamson的2近似方法的实现[31]。如果存在一些零或负的连接成本，我们可以在计算PCST之前直接连接相应的循环。使用PCST而不仅仅是贪婪地连接循环，还可以集成那些单独来看不够有价值但与其他循环组合起来有价值的循环。

使用PCST，我们现在从PCST中的任意一个循环开始，通过深度优先搜索迭代地合并循环（使用Dijkstra方法计算的加倍路径）。每当我们合并两个循环时，路径会创建额外的停靠点，这些停靠点可能比最初计算的连接路径更便宜。然而，我们不需要重新计算整个Dijkstra树，而只需简单地减少相应边的成本，并让减少的成本传播。注意：在将循环与加倍路径连接时，实际上是从循环中替换通道。从这些被移除的通道起始的最短路径将变得无效。由于这种情况很少发生且可以检测到，所以只有在即将使用这样的无效最短路径时才应进行重新计算。

2.7 第七步：局部优化

在连接循环以形成巡回路径之后，连接部分通常具有很高的冗余，如图12所示。幸运的是，我们可以将第2.5节的局部优化方法扩展到连接的巡回路径。挑战在于确保在局部优化后巡回路径仍然保持连接。使用的MIP不会强制连接性，并可能再次断开巡回路径。一种朴素的方法是只接受保持连接性的局部改进，并且丢弃所有其他改进。这当然相当严格，我们可以找到一个更好的解决方案。

在MIP中进行子路径消除比例如旅行商问题更加困难：不仅所有的访问是可选的，而且两个巡回路径可以交叉而不连接。因此，仅强制两条边离开一个连接组件并不能得到期望的结果。在[41]中，我们实际上有一个对应的MIP。因为我们已经从一个巡回路径开始，并且知道我们必须连接一个内部解（在要优化的小区域内）到固定的外部解决方案，所以我们可以设计一个更简单的分离约束。有两种类型的子路径：一种完全位于区域内，另一种只部分位于区域内。只有当局部解错误地连接了外部解时，我们才会得到第二种类型的子路径的不可行解。然而，这两种类型都可以处理。我们要么希望子路径C溶解或成为连接路径的一部分。为此，子路径的顶点通道需要未使用，或者离开子路径的顶点通道已被使用。我们选择第一种类型的任意顶点通道，并要求第二种类型的总和大于它。请注意，这假定存在一个外部的固定解决方案，否则不是精确的。

让$X_A$是包含在区域$A$中并可以通过局部优化修改的顶点通道变量。这包括所有的变量$x_{uvw}$，其中$u、v、w$属于$A$。如果$u$或$w$不在$A$中，那么$v_u$ resp. $v_w$ 必须在解决方案中使用，即该边连接了可变的内部解决方案和固定的外部解决方案。所有其他变量都是固定的。让$X_C$是由子路径$C$使用的顶点通道变量。让$X^\prime_C$是与子路径$C$共享一条边但不在$X_C$中的顶点通道变量。这些是离开路径$C$的轨迹的顶点通道。现在，我们可以陈述一个约束，如果它是由$X_C$上的优化创建的，则会消除$C$。

$$\begin{equation} \sum_{x\in X^\prime_C \cap X_A}x\geq x_c\quad x_c\in X_C\cap X_A,\text{C is subtour} \end{equation}$$

有更有效的选项来连接，例如，更远的子路径，但这几乎不适用于仅优化小区域的情况。如果 MIP 在固定次数的迭代中没有为区域$A$提供连接的解决方案，我们将丢弃不可行的解决方案，不会在这次迭代中更改$A$。多次应用这种方法可以显著改善解决方案，如图13所示的示例。

3. 评估

在这一部分，我们将评估算法在一组基准实例上的性能。首先，我们评估新优化对性能的影响，然后评估在基准实例上的整体性能。我们使用联合和差异来生成我们的基准的 500 个随机实例。生成过程是受监督的，参数是手动调整的，以创建模仿复杂的农田、建筑、建筑群和其他真实世界情境的实例。有价值的区域和成本增加的区域也是通过随机放置厚多边形来选择的，可能存在重叠并且相加。这些实例的示例可以在图2中看到。这些示例的选择是随机的，因此应该反映在 500 个实例中的分布。所有实验在配有 AMD Ryzen 7 5800X（8 ×3.8 GHz）CPU 和 128 GB RAM 的 Ubuntu 工作站上运行。代码在 Python 3.8.8 中运行，并使用 Gurobi 9.1.2。

3.1 局部优化

在第一个实验中，我们评估了原始算法中未考虑的局部优化步骤。重要的问题是：（1）使用这种优化我们能够提高多少解决方案？我们必须确保提高是值得额外复杂性的。 （2）我们应该重点优化循环覆盖还是巡回路径？尽管巡回路径是最终结果，但循环覆盖的优化成本较低。 （3）迭代次数和区域大小有多大影响？运行时间随着迭代次数的增加呈线性增加，但随着区域大小的增加呈指数增加。但是，问题的NP困难性也意味着迭代不能完全替代区域。为了回答这些问题，我们计算了以下解决方案：在循环覆盖或巡回路径上进行局部优化，迭代次数为0、10、25、50、100和200次，区域大小为50个顶点。此外，我们还计算了在循环覆盖或巡回路径上进行50次局部优化的解决方案，但区域大小变化为0、10、25、50、75和100个顶点。

图14a中的结果表明，面积大小为50个顶点的优化在部分覆盖的两个步骤中都带来了明显的改进。循环覆盖上的10次迭代已经将最优性差距（与下界相比）降低了约10%。进一步的迭代失去了效果，这是可以预料到的，因为我们优先考虑了昂贵的区域，但仍然可以看到改进。虽然优化在循环覆盖上取得了成功，但在巡回路径上表现得更加出色。在这里，前10次迭代将最优性差距降低了超过20%。进一步的迭代也仍然比循环覆盖的迭代更强，但它们的改进仍然很快下降。这意味着循环覆盖已经接近最优，但将循环连接成巡回路径并不是非常高效的。巡回路径上的局部优化可以轻松找到连接解决方案中的（局部）次优部分并明显改进它们。

图14b中的结果令人惊讶地非常相似：将区域大小翻倍与将迭代次数增加四倍有类似的效果。一个区别是，对于优化循环覆盖，较大的区域比对于巡回路径更重要。仅优化包含10个顶点的小区域几乎不会改善解决方案。另一方面，对于巡回路径，这种小区域已经可以产生显著的差异。这是非常有用的信息，因为优化10个顶点非常快，仍然可以通过蛮力方法完成。因此，我们可以在短时间内对这样小的区域进行多次迭代。较大的区域仍然具有优势，大小为100的50次迭代大致与大小为50的200次迭代一样有效。

迭代次数和区域大小的运行时差异可以在图15中看到。令人惊讶的是，较大区域的运行时几乎呈线性增长（注意：迭代的x轴是指数的，但对于区域来说几乎是线性的）。然而，应该谨慎使用这些数据，因为它可能会被扭曲。实施仅针对质量进行了优化，而不是运行时。用于确保我们没有意外断开路径并需要插入约束的连通性检测尤其低效。它总是检查整个解决方案，而不仅仅是分析更改的部分，这个过程是用纯Python编写的。这使得路径变体有显著的开销，可以消除。路径变体仍然会变慢，因为解决方案经常会断开连接，并且需要使用额外的约束重新连接。对于较大的优化区域，应该开发和使用额外的约束。对于下一个实验，我们在两个步骤中都使用50个顶点的25次迭代。对于巡回路径，我们最多使用10个切割平面迭代。

为了评估算法的整体性能，我们再次在500个实例上运行算法。图16a中的图表显示了解决方案的质量如何随实例大小而发展。质量再次通过目标值与分数解提供的下界之间的差异来衡量（参见第2.2节）。我们可以看到目标值大约比分数解高出10%到15%，正如我们在先前的实验中已经看到的那样。然而，在这里我们观察到，对于较大的实例，质量会稍微下降。基于工具半径为1.0，较大的（图形）实例具有数千个顶点。这种恶化可能会趋于稳定，但数据相对嘈杂，范围太小，无法做出任何确定的假设。对于较低的转弯成本，间隙通常较小，但这并不令人惊讶，因为转弯成本使问题的组合更加复杂。这至少会影响分数解的质量，该解为我们提供了下界。从解决方案中无法确定实际解是否具有更大的最优性差距。

3.3 运行时间

本文的主要重点是解决方案的质量，但运行时也是一个重要因素。原始算法能够解决具有30万个顶点以上的实例，尽管这可能需要数小时并需要强大的工作站。本文中考虑的实例只有几千个顶点，因为实现仅针对质量进行了优化，而不是运行时。尽管相对较小，但这些实例仍然不是微不足道的，如图2所示。这些实例需要几分钟的运行时间，如图16b所示。原型的效率可以在多个地方进行改进。然而，与原始算法相比，存在内在的挑战。首先，原始算法受益于方形网格的简单性，该网格仅具有三种类型的通道。其次，它利用了基本的整数算术，而本文中的算法需要浮点算术，可能会影响收敛行为。

4. 总结

在本文中，我们展示了如何将网格图上的一种常数近似算法适应复杂多边形环境导出的任意网格上的覆盖路径问题。虽然在这个过程中可能会丧失近似系数（如果网格不恰好是完美的方形网格），但我们证明了该算法在实践中仍然具有较低的最优性差距。此外，我们展示了它在处理部分覆盖和适应异构通道成本方面的多功能性，提供了在覆盖质量和时间效率之间进行权衡的灵活性。这种改编为实际机器人应用提供了坚实的理论基础，可以计算高效的覆盖路径。潜在的未来工作包括问题的多机器人变体，其中可以使用固定数量的机器人。如果只关心总成本的总和，则可以通过仅调整连接步骤（步骤6）来扩展当前方法。如果关心各个成本，那么可以通过不仅基于线性松弛来决定方向（步骤3），而且通过将线性松弛扩展到多个机器人（实际上为每个机器人复制它），并另外决定使用哪个机器人来推广所提出的方法。一个在实践中相关但在算法上具有挑战性的变体是在给定预算的情况下最大化覆盖质量。在我们的方法中的一个问题是依赖于分数松弛，而已知它对于预算约束是弱的。然而，通过附加约束或执行一些分支步骤，潜在地可以改进线性松弛。

附录A. 高级问题定义

实际上，由于篇幅限制，本文实现的算法有一个更复杂的问题定义。在接下来的内容中，我们将讨论完整的潜在问题定义以及背后的动机，以及有关离散化的更多细节。本节内容不是理解本文所必需的，但它可以更深入地了解实现的意图。

A.1 几何模型

我们在这一节中对我们的优化方法进行评估，使用了一个简化的但仍然通用的二维几何模型。这个模型可以适用于许多现实场景，而许多规格并不是由于算法限制，而只是为了简化评估而使用的。虽然基于模拟的评估可能会产生更真实的结果，但它会变得不太通用，并且需要大量真实情况的实例，这在实际操作中很难获得。

首先让我们讨论一下我们如何建模机器人。在接下来的内容中，我们主要谈论机器人，但通常包括所有类型的工具，如铣床或无人机。我们将机器人建模为半径为$r > 0$的圆，其位置$p\in \mathbb{R^2}$由其中心点定义。机器人立即覆盖其下方的所有内容，也就是说，如果它位于位置$p$，覆盖区域$\text{Cov(p) }= {p^\prime \in \mathbb{R^2}|\quad||p - p^\prime||\leq r}$。这使得机器人具有旋转不变性，并简化了许多计算。这种圆形覆盖看起来乍一看对于割草机等应用似乎不太现实；但在一次运动中，垂直于轨迹的线的覆盖几乎与圆形相同。

环境，例如墙壁或障碍物，可能会限制机器人的移动。我们将可行区域表示为$F \subset \mathbb{R^2}$，即机器人所有可行位置的集合，并通过一个（非简单）多边形来近似表示它。在示例和评估中，我们从一个更大的多边形开始，代表着一个房间，然后通过删除距离边界太近的部分来缩小它。$F$不需要与可覆盖区域重叠，这允许我们将机器人的形状与其覆盖分开。

我们将机器人的轨迹定义为一系列闭合的航点链,$\omega_0,\omega_1,\cdots, \omega_{|T|-1}\in F$,机器人在航点之间直线移动。我们用$\text{SEGMENTS}(T)=\omega_0\omega_1,\omega_1\omega_2,\cdots,\omega_{|T|-1}\omega_0$表示相应的线段,我们要求所有轨迹中的线段$s \in \text{Segments} (T)$完全包含在可行区域$F$中。在接下来的部分中，轨迹也被称为巡回路线。多个（闭合的）轨迹的中间解决方案，这些轨迹仍然需要连接成一个巡回路线，称为循环覆盖（cycle cover），它的元素是循环或子巡回路线。

除了可行区域$F$外，我们还有有价值区域和昂贵区域。有价值区域$\mathbb{Q}= Q_0, Q_1,\cdots \subset \mathbb{R^2}$，带有权重$t(Q_i) \in \mathbb{R}^+$，表示我们想要覆盖的部分。昂贵区域$\varepsilon= E_0, E_1,\cdots \subset F$，带有权重$m(E_i) \in \mathbb{R}^+$，代表具有增加巡回成本的区域。这两种类型的区域同样被多边形近似，以简化计算。

目标是计算一个可行的巡回路线，最大化覆盖价值并最小化巡回成本。为了将覆盖最大化和成本最小化结合起来，我们将覆盖最大化转化为最小化问题。这通过考虑机会损失来实现，即错失区域的价值,

$$\min_T \text{COVERAGELOSS}_{Q,r}(T)+\text{TOURCOST}_{\varepsilon}(T)\\ s.t.\quad s\subseteq F \quad \forall s\in \text{SEGMENTS}(T)$$

这个目标相对于其他目标的优势在于，它的下界是零，这使得更好的比较成为可能。我们定义了以下覆盖损失和巡游成本。

令$C_r(T)={p\in \mathbb{R}^2|\exist s\in \text{SEGMENTS}(T),p^\prime \in s: p\in \text{COV}(p^\prime)}$表示一次巡回的覆盖区域。需要注意的是，将一个点包括在一个线段中是指该点可以位于线段上的任何位置，而不仅仅是在线段的端点上。这使得我们可以通过最大可实现的覆盖值减去实际实现的值来正式定义覆盖损失。

$$\text{COVERAGELOSS}_{Q,r}(T)=\sum_{Q\in \mathcal{Q}}|Q|\cdot t(Q)-\sum_{Q\in \mathcal{Q}}|Q\cap C_r(T)|\cdot t(Q)$$

巡游成本包括加权距离和转向角度。

$$\text{TOURCOST}_{\varepsilon}(T)=\lambda_0\cdot \text{DISTCOST}_{\varepsilon}(T)+\lambda_1\cdot \text{TURNCOST}_{\varepsilon}(T)$$

两个权重$\lambda_0$,$\lambda_1\geq 0$用于加权距离和巡游成本，并且在我们的实验中会对它们进行变化。
让$\mu_{\varepsilon}:\mathbb{R}^2\rightarrow \mathbb{R}^+$定义在工具位置上的成本乘数，这允许我们模拟由环境引起的局部成本变化。它的计算方式是：$\mu_{\varepsilon}(p)=\prod_{E\in \varepsilon,p\in E}m(E)$,距离代价定义为：

$$\text{DISTCOST}(T)=\sum_{s\in \text{SEGMENTS}(T)}\int_{p\in s}\mu_\varepsilon(p)dp$$

对于$\varepsilon = \emptyset$，这变为$\sum_{s\in Segments (T)}||s||$。而转弯成本只发生在航路点，并且也受到成本乘数的影响。

$$\text{TURNCOST}_\varepsilon(T)=\sum ^{|T|-1}_{i=0}\mu_\varepsilon(\omega_i)\cdot TURN(\omega_{i-1},\omega_i,\omega_{i+1})$$

$Turn (p_0, p_1, p_2)$ 表示在穿越$p_0\rightarrow p_1\rightarrow p_2$时在$p_1$处的转向角度。航路点的索引是取模$|T|$形成一个循环的。

A.2 离散化

在应用我们的近似技术之前，我们需要将多边形区域转化为潜在航路点的图形。然后，我们只需要在图中找到一条巡游路线，其中每个顶点都提供一定的覆盖，而巡游成本基于使用的边和边的转换来计算，而不是复杂的几何问题。最简单和最常见的策略是在可行区域上放置一个规则的正方形网格。完全包含在内的点和边将成为我们的图形。这也直接允许我们使用Fekete和Krupke的算法[25]。然而，并不总是最优的。其他选择包括使用规则的三角形网格或不规则生成的网格。为了降低计算成本，通常最好减少顶点数并使顶点的边度低。不仅计算成本会增加，顶点附近的松弛质量也会下降。图17中显示了各种不同网格的示例。

计算边缘成本和顶点处的转向成本是直截了当的，可以直接使用等式(A.2)中的定义。图形只是实际解空间的一个子集，因此我们可以预先计算图形使用的各个部分的成本。覆盖顶点的价值更加复杂。我们可以简单地分配机器人在这一点上覆盖的价值，但这很容易高估或低估真实价值。如果真正的覆盖发生在移动到和从顶点时，它可能会低估价值。如果其他顶点靠近并且覆盖区域重叠，它可能会高估价值。通常情况下，如果图中的值之和等于原始实例中的最大值，那就很好。我们可以简单地缩放所有的值来实现这一点，但由于图可能具有异构分布，使用Voronoi图是一个更好的选择。Voronoi图是计算几何学中的一个经典方法，它将区域划分为每个顶点分配最接近的区域。使用这些区域的值给我们一个等于原始值的值分配，同时也对顶点的邻域敏感。这可以在图18中看到。

在接下来的内容中，我们用$G= (P, E)$表示生成的图形，称$P$为潜在航路点。每个在$G$上的巡游都由$E$中的段组成，这些段完全包含在可行区域内，因此是原始多边形内的可行巡游。我们用$val(p)$表示分配给航路点$p\in P$的覆盖值，它对应于$p$的Voronoi单元的覆盖值。边缘$pp^\prime \in E$的距离成本由$dist(p, p^\prime ) =\int pp^\prime \mu_\varepsilon (x)dx$定义，根据等式(A.3)。通过邻居$n$和$n^\prime \in N(p)$通过$p$的转向成本由$turn(n, p, n^\prime) = \mu_\varepsilon (p)·Turn(n, p, n^\prime)$根据等式(A.4)定义。
获得一个良好的图形是一个基本问题，整个附录C都集中讨论了这个问题。

B 第3步的实施细节

在这一步中，我们将实例转化为可以使用最小权重完美匹配来计算积分循环覆盖的形式，也就是说，允许解决方案由多个巡回组成。在没有转向成本的情况下，实际上可以在多项式时间内计算最优循环覆盖，因为边的成本是独立的。有了转向成本，边的成本取决于前一边的方向，使得即使在网格图中，这个问题也是NP难的[24]。我们使用前一步的分数解来预测相应的方向，并再次使成本独立。

我们可以将这个过程想象成将每个航路点替换为一个 epsilon 长度的线段，就像图5中所示。在端点上计算最小权重的完美匹配，就像图19中所示，将得到包括所有这些线段的最优循环覆盖。连接边上的必要转向在每个连接边上都是固定的，因此可以与距离一起计算在边权重中。我们将这些 epsilon 长度的线段称为原子带（atomic strips）。跳过一个航路点的可能性可以通过在其原子带的两个端点之间添加一条边来实现，其权重为未覆盖的部分。

原子带的方向具有基本重要性：如果我们判断它们正确，那么最小权重的完美匹配实际上对应于航路点上的最优循环覆盖。如果我们错误地判断了原子带的方向，最小权重的完美匹配可能会执行昂贵的转向操作来集成它。幸运的是，在航路点处转向时，精确的方向不那么重要，因为最优方向的范围随着转向角度的增加而增加，如图20所示。对于一个 U 形转弯，每个方向都是最优的。通道越直，良好的方向就越重要；但通常这些情况很容易从分数解中判断出来。这一观察使我们能够将方向限制为相邻边的方向，即邻居的方向。在接下来的内容中，我们用相邻航路点$N(p)$来表示航路点$p\in P$的原子带的可用方向。如果航路点$p$的原子带具有方向$n\in N(p)$，则其一个端点朝向$n$。

一个航路点可能需要被多次穿越，因为我们受限于网格$G= (P, E)$中的通道。这可以通过传递边来轻松实现，即两个连续的边$uv$和$vw\in E$会自动创建一条边$uw$，其成本为两者的组合成本。然而，我们在[25]中学到，引入可选的原子带并且只允许直接连接在性能上要好得多。可选的原子带可以通过在其端点之间添加一个权重为零的边来实现，参见图21。我们称航路点的非可选原子带为主导原子带。

在一个像 Fekete 和 Krupke [25] 中的正方形网格中，如果我们为每个邻居添加一个原子带，并且将分数解中使用最多的那个声明为主导，那么我们可以获得一个4近似解。在一个三角形网格中，可能会有一些被穿越线性次数的航路点，就像图6中所示，但这是一个人工生成的实例。在我们的实例中，每个航路点通常只被覆盖一次或两次。由于每个原子带都会增加计算复杂性，我们将原子带的数量限制为一个常数$k$，并允许每个航路点$p\in P$最多有一个原子带与每个邻居$n\in N(p)$。任务是选择一个子集$A\subseteq N(p)$，其中$|A| \leq k$作为原子带，并确定其中的主导原子带。

如果$k\geq |N(p)|$，我们可以简单地选择$A=N(p)$。这允许我们在没有额外开销的情况下使用任何通道两次，因为任何航路点通道要么有两个邻居，每个都有一个最优的原子带，要么通道是一个U形转弯。如果$k < |N(v)|$，事情会变得更加复杂，因为我们想要优化三个经常相互矛盾的目标：

• 我们想要改善期望情况，即具有最高可能性的通道应尽量代价低。

• 我们想要最小化平均情况下的成本开销，即任何通道的平均开销。

• 我们想要最小化最坏情况下的成本，即最坏情况下的成本。

对于这种情况，我们的策略包括两个阶段。首先，我们根据分数解中的边使用情况选择原子带。这优化了期望情况。其次，我们通过最小化未在分数解中使用的通道的平方开销之和来填充剩余的原子带。这优化了平均情况和最坏情况（使用更高的指数会将焦点转移到最坏情况）。具体策略如算法1中所示。$FS(v,\omega)=\sum _{u\in N(v)}x_{uvw}$表示分数解中使用边$vw\in E$的情况，$OH(uvw, A)$表示如果通道$uvw$必须使用$A\subseteq N(v)$中的原子带，那么最小的开销。开销对应于额外的转向成本，用于容纳（可能不对齐的）原子带。如果一个航路点$v\in P$具有覆盖值，即$val(v)>0$，我们仍然必须选择主导带，否则会损失机会。对于非常直的通道，可能没有可以在没有额外开销的情况下使用的原子带（这也可能是由于数值问题引起的）。因此，我们提出了一种更动态的方法。

从集合A中选择主导带的过程是通过原子带的使用情况进行的，如下所示：

$$\text{SELECTDOMINANT}(v,A)=argmax_{\alpha \in A}\sum_{u,\omega \in N(v)}x_{uvw}\cdot \text{USAGE}(uvw,a)$$

让$Turn a(u, v, w )$表示如果通道$uvw$被强制使用原子带$a$时的转向角度。使用取决于强制通道使用它所引起的转向开销。如果没有额外的开销，使用率为 1.0，如果有额外的开销，使用率将以指数方式下降。$\text{USAGE}(uvw, a ) = \lambda^{(Turn_a(u,v,w )−Turn (u,v,w ))/ϕ}$，其中λ是在额外转向角度为ϕ时的使用率，示例见表1。较高的值允许更大的间隙，如果网格不规则，则这是必要的。在实验中，我们使用λ=0.25和ϕ=45°。不同k的示例见图7。