UAV Surveillance Under Visibility and Dwell-Time Constraints: A Sampling-Based Approach

TL;DR

This paper presents a sampling-based framework for UAV visual surveillance planning that accounts for visibility, dwell-time, and viewing constraints, optimizing routes with a heuristic approach.

Contribution

It introduces a novel multi-objective optimization method for UAV path planning that handles complex viewing constraints and approximates the problem as a generalized TSP.

Findings

The proposed heuristic effectively solves complex surveillance routing problems.

The framework can incorporate various viewing parameters and constraints.

Numerical studies demonstrate the approach's practicality and efficiency.

Abstract

A framework is introduced for planning unmanned aerial vehicle flight paths for visual surveillance of ground targets, each having particular viewing requirements. Specifically, each target is associated with a set of imaging parameters, including a desired (i) tilt angle, (ii) azimuth, with the option of a 360-degree view, and (iii) dwell-time. Tours are sought to image the targets, while minimizing both the total mission time and the time required to reach the initial target. An epsilon-constraint scalarization is used to pose the multi-objective problem as a constrained optimization, which, through careful discretization, can be approximated as a discrete graph-search. It is shown that, in many cases, this approximation is equivalent to a generalized traveling salesperson problem. A heuristic procedure for solving the discrete approximation and recovering solutions to the full routing problem is presented, and is shown to have resolution completeness properties. Algorithms are illustrated through numerical studies.

Tables

Table 1. Table 1: Target Input Data

$T_j$	$t_j$	${\mathrm{𝙱𝚎𝚑}}_j$	${\tau }_j$	$[{\varphi }_j^{\text{A}}-{\Delta }_j^{\text{A}},{\varphi }_j^{\text{A}}+{\Delta }_j^{\text{A}}]$	$[{\varphi }_j^{\text{T}}-{\Delta }_j^{\text{T}},{\varphi }_j^{\text{T}}+{\Delta }_j^{\text{T}}]$
$T_1$	$(5000,-5000)\text{ m}$	FULL	$2$	$-$	$[\frac{\pi }{8},\frac{3\pi }{8}]$
$T_2$	$(4300,-1750)\text{ m}$	ANGLE	$1$	$[\frac{\pi }{4},\frac{3\pi }{4}]$	$[\frac{\pi }{8},\frac{3\pi }{8}]$
$T_3$	$(0,4000)\text{ m}$	FULL	$3$	-	$[\frac{\pi }{8},\frac{3\pi }{8}]$
$T_4$	$(-8000,-2000)\text{ m}$	ANY	$1$	-	$[\frac{\pi }{8},\frac{3\pi }{8}]$
$T_5$	$(-2000,8000)\text{ m}$	ANGLE	$0$	$[\frac{3\pi }{2},2\pi ]$	$[\frac{\pi }{8},\frac{3\pi }{8}]$

Table 2. Table 2: Spacing Conditions

Spacing Condition	$\delta r$	$\delta \theta$	$\delta \alpha$
1	$1000\text{ m}$	$\pi$	$\pi$
2	$500\text{ m}$	$\pi$	$\pi$
3	$500\text{ m}$	$\pi /2$	$\pi /2$
4	$250\text{ m}$	$\pi /2$	$\pi /2$
5	$250\text{ m}$	$\pi /4$	$\pi /4$
6	$125\text{ m}$	$\pi /4$	$\pi /4$
7	$125\text{ m}$	$\pi /8$	$\pi /8$

Table 3. Table 3: Target Input Data

$T_j$	$t_j$	${\mathrm{𝙱𝚎𝚑}}_j$	${\tau }_j$	$[{\varphi }_j^{\text{A}}-{\Delta }_j^{\text{A}},{\varphi }_j^{\text{A}}+{\Delta }_j^{\text{A}}]$	$[{\varphi }_j^{\text{T}}-{\Delta }_j^{\text{T}},{\varphi }_j^{\text{T}}+{\Delta }_j^{\text{T}}]$
$T_1$	$(2131.8,1026.7)\text{ m}$	ANY	$0$	$-$	$[\frac{\pi }{6},\frac{\pi }{3}]$
$T_2$	$(-13840,-5833)\text{ m}$	ANY	$1$	$-$	$[\frac{\pi }{8},\frac{3\pi }{8}]$

Equations

$$\begin{split}\text{Minimize:}\hskip 8.53581pt&\text{Closed Trajectory Time},\\ \text{Subject To:}\hskip 8.53581pt&\text{Initial Maneuver Time}\leq\epsilon,\\ &\text{Dynamic Constraints Satisfied (Sec.~{}\ref{sec:UAV_specs}), and}\\ &\text{Correct Dwell-Time Maneuvers Performed}\\ &\hskip 14.22636pt\text{at Each Target (Sec.~{}\ref{sec:target}),}\end{split}$$

$$\begin{array}{rl}\text{Minimize:}& W(v_M,v_1)+\sum _{k=1}^{M-1}W\left(v_k,v_{k+1}\right)\\ \text{Subject To:}& W\left(v_0,v_1\right)\le ϵ,\text{ and}\\ & {\text{TAR}}_{v_{k_1}}\ne {\text{TAR}}_{v_{k_2}},\text{ for any }{k}_1\ne k_2,\end{array}$$

$${\text{INL}}_{ϵ}:=\left\{v\in \bigcup _{j=1}^M{\text{DWL}}_j\right|\text{DIST}(v_0,v)\le ϵ\},$$

$$\text{cramped}{{\text{CLS}}_{ϵ}:=\bigcup _{\sigma \in \text{Perm}(M)}({\text{DWL}}_{\sigma (1)}\cap {\text{INL}}_{ϵ})\times {\text{DWL}}_{\sigma (2)}\times \cdots \times \text{DWL}}_{\sigma (M)},$$

Full text

1 Introduction

The use of autonomous mobile sensors is becoming increasingly common in civilian and military applications, since they can provide support in tasks that are too dangerous, too expensive, or too difficult for humans to perform unaided. Tasks that can benefit from autonomous sensors include search and rescue, forest fire or oil spill monitoring, surveillance and reconnaissance, transportation and logistics management, and hazardous waste cleanup [1, 2, 3]. These applications require intelligent and practical strategies to govern autonomous behavior in the presence of numerous constraints that arise in realistic missions.

This article considers a particular mobile sensor scenario that is of interest, e.g., in military operations [4], where a single fixed-wing Unmanned Aerial Vehicle (UAV) collects visual sensor data within a large geographic environment. Here, the UAV is equipped with a gimbaled camera, and must provide surveillance imagery of multiple static ground targets, each having associated imaging constraints. These pre-specified constraints include (i) a desired tilt angle with tolerances, (ii) a desired azimuth with tolerances, including the option of a 360-degree view, and (iii) the amount of time that the UAV should dwell before moving to the next target. The goal is to construct flight paths that are optimal in some sense, while simultaneously allowing each target to be imaged with the prescribed specifications.

Ideally, we seek flight paths that simultaneously minimize two metrics: (i) the time required for the UAV to image all targets and return to the initial target, and (ii) the delay between the mission onset and the time the UAV reaches its first target. The latter goal is important because sensory data often plays a crucial role in the dynamic development of mission strategies, i.e., mission critical planning often cannot progress until some sensory information is provided. This is true, for example, in reconnaissance missions where sensory data supports the operations of a manned aircraft, e.g. [5]. Thus, it is often important that initial imagery is provided quickly, even if at the expense of total tour duration.

Since these two objectives are conflicting in general, we construct solutions by constraining the second performance metric and optimizing over the first. This is a standard scalarization approach for addressing multi-objective problems, since solutions of the constrained problem allow for the construction of Pareto-optimal fronts [6]. We develop a discrete-approximation strategy for constructing high-quality solutions to the scalarized problem, naturally leading to a complete heuristic framework for the full, multi-objective mission. This approach produces reasonable flight paths within a modular framework, while also leveraging existing solution strategies to classical routing problems to ensure straightforward, effective, practical implementation.

Specifically, the contributions of this article are as follows. First, we show how the multi-objective routing problem with both visibility and dwell-time constraints can be posed as a constrained optimization problem through reasonable assumptions about UAV trajectories and dwell-time maneuvers. Indeed, we define visibility regions at each target that reflect imaging requirements, along with a set of feasible dwell-time maneuvers, which are performed within the regions, that both accommodate the dynamic constraints of the UAV and are intuitively straightforward. These constructions are then translated into a constraint set and incorporated into a precise problem statement, which represents an $\epsilon$-constraint scalarization of the multi-objective problem.

Then, we illustrate how the constrained optimization problem, having an infinite solution space, is approximated by a discrete problem, having a finite solution space. This discrete approximation implicitly considers both the time required for appropriate dwell-time maneuvers and the visibility constraints, as it is formulated through a selective sampling procedure. Namely, the UAV configuration space is carefully sampled to obtain a discrete set whose elements are each paired with a feasible dwell-time trajectory. Elements of this set form the nodes of a graph whose edge weights are defined via an augmented Dubins distance that incorporates required dwell-times. The discrete approximation is defined as an optimal path-finding problem over the resultant graph.

Next, we present a novel heuristic method for solving the discrete approximation that leverages solutions to standard Generalized Traveling Salesperson Problem (GTSP) instances. We show that our heuristic method produces feasible solutions and, in many cases, maps optimal solutions of a GTSP directly to optimal solutions of the discrete problem.

Finally, we incorporate these constructions into a complete heuristic framework to produce high-quality solutions to the full, multi-objective routing problem. We prove that the heuristic is resolution complete in a specific mathematical sense, and finish by illustrating our methods numerically.

We note that our work is most closely related to that of [7], which uses a similar sampling-based framework in order to address a Polygon-Visiting Dubins Traveling Salesperson Problem (PVDTSP). However, the framework in [7] considers a single performance metric (total tour time), and cannot incorporate non-trivial UAV dwell-time behaviors. As such, our work is, in some sense, an expansion of [7] to incorporate a more general set of geometric and temporal imaging behaviors, and accommodate an additional performance metric that is reflective of realistic mission scenarios. In fact, the more general solution framework herein reduces to that of [7] in particular cases (Remark 3).

2 Related Literature

A large amount of recent research focuses on the development of coordination strategies for UAV applications (e.g. [8, 9, 10, 11]). However, UAV research is a sub-class of a much larger body of literature that addresses algorithmic design and high-level reasoning for general mobile sensor applications [12, 13]. These applications often necessitate solutions to complex routing problems, which may involve logical, temporal, and spatial constraints, as well as environmental uncertainty. When construction of global optima is not feasible, heuristics are often used to permit the real-time construction of solutions for use within practical systems. For example, variations on the classic Traveling Salesperson Problem (TSP) [14, 15, 16] arise frequently and, since the TSP is NP-hard, global optima usually cannot be consistently found in reasonable time. However, several new insights have been developed over the last decade for the classical TSP, along with a number of variations [17, 18, 19], and sophisticated heuristic solvers, e.g. [20] can quickly construct high-quality solutions for practical use. Other problems that are common in mobile sensor applications include the construction region selection policies [21, 22], the derivation of static and dynamic coverage schemes [23, 24], the development of persistent task execution frameworks [25], and the design of load balancing strategies [26, 27].

The problem herein is loosely interpreted as a generalization of both the Polygon-Visiting Dubins Traveling Salesperson Problem (PVDTSP) [7, 28] and the Dubins Traveling Salesperson Problem with Neighborhoods (DTSPN) [29]. The PVDTSP and the DTSPN are variations on the classic Dubins TSP (DTSP) that require vehicles to visit a series of geometric regions rather than discrete points. To incorporate explicit imaging constraints within a multi-objective framework, we adopt a strategy that is, in some sense, an extension of [7], where the authors approximate solutions to a PVDTSP by discretizing regions of interest and posing the resulting problem as a Generalized Traveling Salesperson Problem (GTSP) (also called the Set TSP, Group TSP, (Finite) One-in-a-set TSP, Multiple Choice TSP, or Covering Salesperson Problem). Sampling-based strategies that approximate a continuous motion planning problem with a discrete path-finding problem over a graph are typically called sampling-based roadmap methods [30, Ch. 5]. Such methods are traditionally used for point-to-point planning among obstacles; however, they have also been used for more general Dubins path planning, e.g. [31].

Like the TSP, the GTSP is a combinatorial optimization; however, strategies exist for computing high-quality solutions. The most popular solution approach converts the GTSP into an Asymmetric Traveling Salesperson Problem (ATSP) via a Noon-Bean transform [32]. The availability of efficient solvers for the ATSP, e.g. LKH [20], make such transformations a viable GTSP solution option in practice, though they do not produce global optima in general. Other common GTSP solution approaches make use of meta-heuristics, e.g. [33].

Aside from the explicit visibility and dwell-time constraints, the problem herein differs from typical “TSP-like” problems due to the presence of multiple (potentially conflicting) performance objectives. A vast amount of literature addresses multi-objective optimization problems, some of which focuses on purely mathematical formulations while some targets solutions for particular applications [34, 35, 36]. Solution approaches typically seek to satisfy a pre-determined notion of optimality, the most common being Pareto optimality. Pareto-optimal fronts are usually difficult to characterize directly, so they are typically constructed using scalarization methods, which reduce the multi-objective problem to a single-objective problem whose optimal solutions map to Pareto-optimal solutions with respect to the original problem. The most common scalarization techniques are linear scalarization, where the cost is a linear combination of the objectives, and $\epsilon$-constraint method, where the values of all but a single objective are explicitly treated as optimization constraints [6].

Our approach to the present problem uses an $\epsilon$-constraint method that treats the initial UAV maneuver time as an explicit optimization constraint. In general, no single approach to multi-objective optimization is always superior; the appropriate method depends on the type of information available, the user’s preferences, solution requirements, and the availability of software [35]. We choose the $\epsilon$-constraint method since it (i) is a standard approach to multi-objective optimization, (ii) gives explicit flexibility in the degree to which the initial maneuver time is considered, allowing the approach to be easily tailored to specific mission scenarios, (iii) allows for a heuristic solution procedure that is naturally modular and can leverage solutions to standard routing problems, (iv) promotes straightforward approximation of Pareto-optimal fronts by varying scalarization parameters, and (v) is straightforward and intuitive, making it an attractive approach to both practitioners and theoreticians alike.

3 Problem Formulation

3.1 UAV Specifications

Consider a single fixed-wing UAV, equipped with a GPS location device and a gimbaled, omnidirectional camera. The camera is steered by a low-level controller, which is independent of the vehicle motion controller. For simplicity, we neglect the possibility of camera occlusions. The work herein focuses on high-level UAV trajectory planning, rather than low-level vehicle or camera motion control. We consider planar motion in a global, ground-plane reference frame, assuming the UAV maintains a fixed altitude $a$ and a fixed speed $s$. We model the UAV as a Dubins vehicle [37] with minimum turning radius $r$, neglecting dynamic effects caused by wind, etc. Let $v_{0}\in{}^{2}\times[0,2\pi)$ denote the UAV’s initial configuration (location, heading).

3.2 Target Specifications

Consider $M$ static targets, each with associated imaging, i.e., visibility and dwell-time, requirements. Each target $T_{j}$ is associated with the following (fixed) parameters:

$t_{j}\in{}^{2}$, the location of the target in the ground-plane reference frame, 2. 2.

$\texttt{BEH}_{j}\in\{\texttt{ANY},\texttt{ANGLE},\texttt{FULL}\}$, the required viewing behavior, where ANY indicates no preference for the azimuth of the collected images, ANGLE indicates that the target should be imaged at a specific azimuth, and FULL indicates that a 360-degree view of the target should be provided, 3. 3.

$\left[\phi_{j}^{\textup{A}}-\Delta_{j}^{\textup{A}},\phi_{j}^{\textup{A}}+\Delta_{j}^{\textup{A}}\right]\subset\real$, a range of acceptable azimuths when $\texttt{BEH}_{j}=\texttt{ANGLE}$, as measured with respect to a reference ray in the ground plane (Fig. 1, top), 4. 4.

$\left[\phi_{j}^{\textup{T}}-\Delta_{j}^{\textup{T}},\phi_{j}^{\textup{T}}+\Delta_{j}^{\textup{T}}\right]\subset(0,\frac{\pi}{2}]$, acceptable camera tilt angles as measured with respect to a plane parallel to the ground-plane (Fig. 1, bottom), and 5. 5.

$\tau_{j}\in\mathbb{Z}_{\geq 0}$, the required number of dwell-time “loops.”

Define the visibility region, $\texttt{VIS}_{j}\subset{}^{2}$, for target $T_{j}$ as the set of locations from which the UAV is able to image the target with an acceptable tilt angle and azimuth (that is, a tilt angle within the interval $\left[\phi_{j}^{\textup{T}}-\Delta_{j}^{\textup{T}},\phi_{j}^{\textup{T}}+\Delta_{j}^{\textup{T}}\right]$ and, if $\texttt{BEH}_{j}=\texttt{ANGLE}$, an azimuth within the interval $\left[\phi_{j}^{\textup{A}}-\Delta_{j}^{\textup{A}},\phi_{j}^{\textup{A}}+\Delta_{j}^{\textup{A}}\right]$; if $\texttt{BEH}_{j}\neq\texttt{ANGLE}$, then any azimuth is acceptable). Each $\texttt{VIS}_{j}$ is uniquely defined by the UAV altitude $a$, the location $t_{j}$, the behavior $\texttt{BEH}_{j}$, and the intervals $\left[\phi_{j}^{\textup{T}}-\Delta_{j}^{\textup{T}},\phi_{j}^{\textup{T}}+\Delta_{j}^{\textup{T}}\right]$, $\left[\phi_{j}^{\textup{A}}-\Delta_{j}^{\textup{A}},\phi_{j}^{\textup{A}}+\Delta_{j}^{\textup{A}}\right]$. Algorithm 1 presents the methodology for constructing visibility regions. Note that, if $\texttt{BEH}_{j}\neq\texttt{ANGLE}$, then $\texttt{VIS}_{j}$ is a full annulus centered at $t_{j}$; otherwise, $\texttt{VIS}_{j}$ is an annular sector (Fig. 2).

As such, fixing target locations, each visibility region is parameterized by two radii (lower, upper radial limits) together with two angles (lower, upper angular limits).

The variable $\tau_{j}\in\mathbb{Z}_{\geq 0}$ indicates the number of dwell-time “loops” that are required at the target $T_{j}$. If $\tau_{j}=0$, then the UAV accomplishes its task by passing over any point within $\texttt{VIS}_{j}$. If $\tau_{j}>0$, i.e., non-trivial dwell time is specified, assume the UAV images $T_{j}$ as follows: If $\texttt{BEH}_{j}=\texttt{FULL}$, the UAV makes $\tau_{j}$ full circles around the target location at some constant radius; if $\texttt{BEH}_{j}\neq\texttt{FULL}$, then the UAV selects a pivot point within $\texttt{VIS}_{j}$ and makes $\tau_{j}$ circles about the selected point at radius $r$. Each non-trivial dwell-time maneuver must be performed entirely within the appropriate visibility region. Fig. 3 shows examples of acceptable imaging behavior for various choices of $\tau_{j}$ and $\texttt{BEH}_{j}$.

For the remaining analysis, we assume that imaging parameters are chosen to ensure problem feasibility, i.e., there exists at least $1$ dwell-time maneuver at each target satisfying the aforementioned constraints.

Remark 1** (Feasibility)**

Feasibility is always ensured if tolerance parameters are sufficiently large.

3.3 Problem Statement

The goal is to construct an optimal UAV trajectory with the following characteristics: The UAV begins its tour by moving from its initial configuration to a configuration where it can begin imaging a target and, after the initial maneuver, the UAV follows a closed trajectory, along which it images each target according to the required specifications. Note that by separating the initial maneuver from the remaining closed route, we ensure that the target imaging behavior can be effectively repeated if desired. Recall the performance metrics to be minimized: (i) the time required for the UAV to traverse the closed portion of the generated trajectory (beginning and ending at the first target), and (ii) the time required for the UAV to perform its initial maneuver, i.e., move from $v_{0}$ to the starting point of the closed portion. Since the metrics are conflicting in general, a tradeoff must be made. We consider the following formulation of the route-finding problem.

Problem 1** (Optimal UAV Tour)**

Find a UAV tour (consisting of an initial maneuver, and a closed trajectory) that solves the following optimization problem:

[TABLE]

where $\epsilon\geq 0$ is a constant parameter.

Remark 2** (Scalarization)**

The optimization problem (1) is an $\epsilon$-constraint scalarization of the multi-objective problem. A typical alternative scalarization would instead account for the initial maneuver time within a linear objective: $\alpha\left(\text{Initial Maneuver Time}\right)+\beta\left(\text{Closed Trajectory Time}\right)$ where $\alpha$, $\beta$ are constant parameters. For fixed $\alpha,\beta$, this alternative formulation is equivalent to (1) for some choice of $\epsilon$. Further, parameter selection and subsequent solution of the linear alternative typically requires construction of a Pareto optimal front by solving instances of (1).

Remark 3** (Relation to [7])**

If (i) $\epsilon$ is large (the initial maneuver is inconsequential), and (ii) $\tau_{j}=0$ for all $j$ (dwell-times are trivial), then solving Problem 1 is equivalent to solving a Polygon-Visiting Dubins Traveling Salesperson Problem, as in [7]. Our methods are loosely based on [7], though we consider a more general multi-objective framework that also incorporates non-trivial dwell-times.

4 Discrete Approximation

Problem 1, which explicitly considers all target imaging constraints, has an infinite number of potential solutions and is difficult to solve. However, by carefully sampling the UAV configuration space, we can pose a finite, discrete alternative whose optimal solutions approximate those of the original problem. The discrete approximation is still a combinatorial search; however, it is closely related to standard path-finding problems, allowing us to leverage existing solvers to produce high-quality sub-optimal solutions. This section develops the discrete approximation of interest.

4.1 Configuration Space Sampling

Recalling the Dubins model, we sample the UAV configuration space to obtain a finite collection of points of the form $v\coloneqq(x,\theta)\in{}^{2}\times[0,2\pi)$. These points serve as the basis for the discrete approximation to Problem 1. Specifically, we choose a set of points that each represent the starting and ending configuration of an appropriate dwell-time “loop” at some target (valid dwell-time maneuvers each start and end at the same configuration). That is, each sampled point $v\coloneqq(x,\theta)$ has heading $\theta$ that points in a direction tangent to a valid dwell-time loop (associated with some target $\texttt{TAR}_{v}$) passing through the location $x\in\texttt{VIS}_{j}$ (Remark 4). By explicitly pairing each $v$ with its target $\texttt{TAR}_{v}$, this procedure creates a natural one-to-one mapping between the generated points and a set of feasible dwell-time maneuvers. As such, subsequent graph formulations can “disregard” dwell-time constraints by using an augmented graph distance. Fig. 4 shows examples of valid sampled sets associated with some $T_{j}$, for various $\texttt{BEH}_{j}$ and $\tau_{j}$ values. Here, the red dot is the sampled point’s location and the arrow represents its heading (distinct points can have the same planar location).

Algorithm 2 outlines the sampling process. Here, each set $\texttt{DWL}_{j}$ is defined thusly: If $\tau_{j}\neq 0$, let $\texttt{DWL}_{j}$ be the set of points $v\coloneqq(x,\theta)\in\texttt{VIS}_{j}\times[0,2\pi)$ having location $x$ that lies on the circular image of an appropriate dwell-time maneuver and heading $\theta$ that points in a direction tangent to the same circular image at $x$. If $\tau_{j}=0$, let $\texttt{DWL}_{j}\coloneqq\texttt{VIS}_{j}\times[0,2\pi)$. Notice Algorithm 2 allows multiple “copies” of the same configuration be sampled, provided each is associated with a distinct target, i.e., there may exist $v_{k_{1}},v_{k_{2}}\in V$ with identical locations and headings so long as $\texttt{TAR}_{v_{k_{1}}}\neq\texttt{TAR}_{v_{k_{2}}}$.

Remark 4** (Dwell-Time Loops)**

It is possible that some $v\in V$ is tangent to multiple, distinct dwell-time loops associated with $\texttt{TAR}_{v}$. Notice that all such loops have identical radii, i.e., the UAV requires the same amount of time to traverse each. Thus, we can assume without loss of generality that each $v$ is the starting and ending configuration of a single loop associated with $\texttt{TAR}_{v}$.

4.2 Graph Construction

Algorithm 2 returns a discrete set $V$, where each $v\in V$ represents a point in the UAV configuration space ${}^{2}\times[0,2\pi)$ that is the starting and ending point of a valid dwell-time maneuver at $\texttt{TAR}_{v}$. Recalling that the optimal Dubins path between any two configurations $v_{k_{1}},v_{k_{2}}\in{}^{2}\times[0,2\pi)$ is well defined (and easily computed) [37], we utilize Algorithm 3 to construct a weighted, directed graph $G\coloneqq(V\cup\{v_{0}\},E,W)$ that effectively discretizes the solution space of Problem 1. Here, the edge set $E$ contains directed edges connecting each pair of nodes in $V$, along with directed edges connecting the initial UAV configuration $v_{0}$ with each node in $V$. Weights are defined via an augmented distance that includes both the time required to complete the dwell-time maneuver at the source node and the time required to travel between configurations (recall that optimal Dubins paths are asymmetric in general). We are now ready to formally define the discrete approximation to Problem 1 using the graph $G$.

Problem 2** (Discrete Approximation)**

Consider the graph $G\coloneqq(V\cup\{v_{0}\},E,W)$ resulting from Algorithm 3. Find a sequence $v_{1},v_{2},\ldots,v_{M}\in V$ that solves

[TABLE]

where $v_{0}\in{}^{2}\times[0,2\pi)$ corresponds to the initial UAV configuration and $\epsilon\geq 0$ is a constant parameter.

5 UAV Tour Construction

Notice that solutions to Problem 1 can be recovered from solutions to Problem 2. Indeed, given a solution $v_{1},\ldots,v_{M}$ to (2), we recover a feasible solution to (1) by: (i) concatenating the optimal Dubins paths between adjacent nodes in the sequence (appending the path from $v_{0}$ to $v_{1}$ at the beginning and the path from $v_{M}$ with $v_{1}$ at the end) and (ii) appending the dwell-time trajectory associated to each node $v_{1},\ldots,v_{M}$. The remainder of our analysis studies the discrete approximation (Problem 2) and its relation to the continuous formulation (Problem 1).

5.1 Solving the Discrete Problem

We leverage solutions of a classic graph path-finding problem to construct solutions to (2). In particular, we propose a heuristic framework that relates solutions of (2) to those of an (asymmetric) Generalized Traveling Salesperson Problem (GTSP), which is defined for convenience here.

Problem 3** (GTSP)**

Given a complete, weighted, directed graph $\mathcal{G}\coloneqq(\mathcal{V},\mathcal{E},\mathcal{W})$, and a family of finite, non-empty subsets $\{\mathcal{V}_{j}\subseteq\mathcal{V}\}_{j\in\{1,\dots,M\}}$, find a minimum weight, closed path that visits exactly one node from each subset $\mathcal{V}_{j}$.

Remark 5** (GTSP Formulation)**

A common alternative GTSP formulation requires the closed path to visit at least one node from each $\mathcal{V}_{j}$. If edge weights satisfy a triangle inequality, then this alternative and Problem 3 are identical. In what follows, we define GTSP instances on a complete subgraph $\mathcal{G}\subseteq G$, where $G$ is the graph constructed in Algorithm 3. In this case, edge weights in $\mathcal{G}$ represent an augmented Dubins distance and, since the Dubins distance function satisfies a triangle inequality [29], edge weights in $\mathcal{G}$ also satisfy a triangle inequality. Thus, we consider Problem 3 without loss of generality.

Remark 6** (GTSP Solutions)**

The standard GTSP is NP-hard. However, practical strategies exist for quickly constructing high-quality solutions, e.g., transformation of the problem into a standard ATSP and application of a heuristic solver (see Section 2).

Note that, in general, Problem 2 is not equivalent to a GTSP, due to the constraint on the initial maneuver. We can, however, leverage GTSP solution procedures in constructing solutions to the constrained problem. Indeed, a heuristic procedure for constructing solutions to Problem 2 using the solutions to related GTSP instances is outlined in Algorithm 4.

Here, $\texttt{INL}_{\epsilon}$ denotes the set of all nodes in $V$ that can be reached from $v_{0}$ in time less than $\epsilon$. In general, the sequences produced by Algorithm 4 will not be optimal with respect to Problem 2. They will, however, be feasible. Further, if $\texttt{INL}_{\epsilon}$ has a particular structure, then the subset $\texttt{INL}_{\epsilon}^{*}$ (see Algorithm 4) can be chosen to ensure that the GTSP instance (line $5$) is equivalent to Problem 2. The following results make this discussion precise.

Theorem 1** (Feasibility)**

Algorithm 4 produces a feasible solution to Problem 2.

Proof 5.2**.**

The GTSP solution (line $5$) will contain some $v\in\texttt{INL}_{\epsilon}^{*}\subseteq\texttt{INL}_{\epsilon}$. Thus, the permutation operation in line $6$ will produce $v_{1},v_{2},\ldots,v_{M}$ with $v_{1}\in\texttt{INL}_{\epsilon}$. It follows readily that the algorithmic output is feasible with respect to Problem 2

Remark 5.3** (Feasibility).**

Theorem 1 does not require construction of an optimal GTSP solution, i.e., the result holds as long as any feasible GTSP solution is produced in line $5$.

Theorem 5.4** (Equivalence).**

Consider Algorithm 4. Suppose $\texttt{INL}_{\epsilon}$ is nonempty and that there exists an index $\hat{\jmath}\in\{1,\dots,M\}$ satisfying either:

$\texttt{INL}_{\epsilon}\subseteq\{v\in V\;|\;\texttt{TAR}_{v}=T_{\hat{\jmath}}\}$, or 2. 2.

$\{v\in V\;|\;\texttt{TAR}_{v}=T_{\hat{\jmath}}\}\subseteq\texttt{INL}_{\epsilon}$.

If $\texttt{INL}_{\epsilon}^{*}$ (line 3) is chosen as the set of all nodes in $\texttt{INL}_{\epsilon}$ that are associated with $T_{\hat{\jmath}}$, then optimal solutions of the GTSP in line 5 map to those of Problem 2 via the operation in line $6$. That is, if a globally optimal solution to the GTSP in line 5 is produced, then the output $v_{1},\ldots,v_{M}$ of Algorithm 4 is a globally optimal solution to Problem 2.

Proof 5.5**.**

Any solution $v_{1},v_{2},\ldots,v_{M}$ to Problem 2 must contain exactly one node associated to each target, where $v_{1}\in\texttt{INL}_{\epsilon}$. If $\texttt{INL}_{\epsilon}$ only contains nodes associated to $T_{\hat{\jmath}}$ (condition (i)), then no feasible solution to Problem 2 contains any $v\notin\texttt{INL}_{\epsilon}$ with $\texttt{TAR}_{v}=T_{\hat{\jmath}}$. Thus, there is no loss of generality in considering the GTSP defined over the modified graph $\mathcal{G}$ when $\texttt{INL}_{\epsilon}^{*}=\texttt{INL}_{\epsilon}$ (as opposed to defining a GTSP over the subgraph induced by $V$). The same applies when $\texttt{INL}_{\epsilon}$ satisfies condition (ii) and $\texttt{INL}_{\epsilon}^{*}$ equals the set of all nodes associated with $T_{\hat{\jmath}}$, as this implies $\mathcal{V}=V$. Since cyclic permutation of the node sequence does not affect the cost, any optimal solution to the GTSP in line $5$ can be mapped to an optimal solution of Problem 2 via the operation in line $6$.

Figure 5 shows a graphical illustration of when the theorem conditions are met. We note that the conditions required by Theorem 5.4 are typically met whenever target spacing is large.

5.2 Complete Tour Construction

Algorithm 5 is a heuristic procedure that leverages solutions to the discrete approximation (Problem 2) to construct solutions to the full routing problem (Problem 1).

Solutions produced by Algorithm 5 are not optimal in general, though they will exhibit structural characteristics that generally improve in quality (with respect to Problem 1) as the sampling granularity is made increasingly fine. Indeed, Algorithm 5 is a resolution complete in some sense, providing justification for the sampling-based approximation approach. A precise characterization of the resolution completeness properties is included as an Appendix.

6 Numerical Examples

We illustrate our algorithms through numerical examples. For each simulated mission, solutions to Problem 1 are constructed via Algorithm 5, where GTSPs are solved through transformation into an equivalent ATSP (see [32]) that is subsequently solved using the Lin-Kernighan heuristic (as implemented by LKH [20]). In all cases, the set $\texttt{INL}_{\epsilon}^{*}$ (Algorithm 4) is chosen as the set of all points in $\texttt{INL}_{\epsilon}$ associated with some single target (satisfying Theorem 5.4 conditions whenever possible). A slightly modified version of Algorithm 2 is used for sampling in which the number of samples, $N$, associated with each target is not fixed a priori, but instead is determined by creating a grid of samples within the appropriate sampling subsets. The grid spacing is determined by $3$ parameters $\delta r$, $\delta\theta$, and $\delta\alpha$, which represent, loosely, the radial location spacing, the angular location spacing, and the angular heading spacing, resp. The parameters $\delta r$, $\delta\theta$, and $\delta\alpha$ are inversely proportional to the number of samples at each target, and the sampled set is dense in the limit as spacing parameters jointly tend to [math].

6.1 Pareto-Optimal Front

The first example is a $5$ target mission with the following UAV parameters: $r=750$ m, $a=1000$ m, $s=39$ m/s, and $v_{0}=((-2500,500)\text{ m},0)\in\mathbb{R}^{2}\times[0,2\pi)$. Target parameters are shown in Table 1 (tolerances are large enough for feasibility).

The approximate Pareto-optimal front (with respect to $\epsilon$) for Problem 1 as a function of the sampling granularity is shown in Figure 6. The figure also shows illustrations of solutions produced at spacing condition 5 when $\epsilon=65$ s (left) and $\epsilon=205$ s (right). Recalling that (i) Theorem 5.4 does not hold for all $\epsilon$, and (ii) heuristic solvers for GTSPs do not guarantee global optima, to obtain the best approximation of the Pareto-optimal front, the following steps were taken to generate each curve: First, Algorithm 5 was called for a series of $\epsilon$ values, and the resulting initial maneuver/closed trajectory times were recorded. Then, in post-processing, the approximate Pareto-optimal curve was generated by selecting, for each $\epsilon$, the lowest cost route satisfying the initial maneuver constraint out of the solutions produced in the computation stage. Note that increasing $\epsilon$ corresponds to relaxing the initial maneuver constraint, and thus the cost is non-increasing in $\epsilon$. Notice also that the Pareto-optimal fronts shift toward zero as the sampling spacing is decreased.

Remark 6.6** (Pareto-optimal Fronts).**

The optimal cost value in Problem 1 is only sensitive to changes in $\epsilon$ over some finite set $\bigcup_{j=1}^{M}[\underline{\epsilon}_{j},\overline{\epsilon}_{j}]$ that depends primarily on the target and UAV locations. This structure can be exploited for more efficient computation of Pareto-optimal fronts.

6.2 Performance

The next example illustrates the performance of Algorithm 5 in comparison to an incremental, “greedy” alternative that operates as follows: Visibility region creation and configuration space sampling are done using Algorithms 1 and  2. Starting with the initial UAV configuration, each successive UAV destination is chosen by selecting the closest node (Dubins distance) associated with a target that has not yet been imaged. A valid route is constructed by appending dwell-time maneuvers and connecting the last selected configuration ($v_{M}$) with the first ($v_{1}$). We consider a $5$ target mission with the same UAV parameters and the same target locations, imaging behaviors, and tolerances as in the previous example. However, we vary the number of dwell-time loops associated with each target (assume each target requires same number of loops).

The difference between the closed trajectory times produced by the greedy method and those produced by Algorithm 5 as a function of $\epsilon$ under the spacing condition $5$ (Fig. 6) is shown in Fig. 7. Notice that the relative performance of Algorithm 5 improves as both $\epsilon$ and the number of dwell-time loops at each target are increased. In this example, the performance of the greedy search method can be made arbitrarily poor by increasing the number of dwell-time loops. This result is primarily due to targets that require a 360-degree view, since the greedy heuristic generally chooses those points lying on the perimeter of the visibility region, which are very far from the target location. Thus, increasing the number of dwell-time loops (performed at constant radius) can dramatically increase total tour times. As such, Algorithm 5 can provide a significant advantage over similar incremental planning strategies when non-trivial dwell-times are required.

6.3 Resolution Completeness

To illustrate the resolution completeness properties, consider a mission with $2$ targets whose relevant data is in Table 3. The UAV has the same altitude, velocity, and minimum turning radius as the UAVs in the previous examples, but has initial configuration $v_{0}=((0,0)\text{ m},\pi/7)\in{}^{2}\times[0,2\pi)$.

This example has been carefully constructed so that the optimal solution is easily deduced for select $\epsilon$ conditions. In particular, when $\epsilon\geq 25$ s, the optimal tour involves the UAV making an immediate left turn at minimum radius, and then visiting the remaining target for a total closed tour length of $848.62$ s. When $\epsilon=16.26$ s, the UAV is constrained to make a shorter initial maneuver, so the UAV travels straight until it hits the visibility region, and then proceeds with the remainder of the tour for an optimal closed tour length of $881.14$ s. A schematic showing optimal routes for the two $\epsilon$ conditions just described is shown in the top left portion of Fig. 9. (the red curve corresponds to the initial turn for the $\epsilon=16.26$ s case; the remainder of the route is identical). When $\epsilon\geq 25$ s, the problem instance is non-degenerate (Definition .10) and satisfies the conditions required for resolution completeness (Theorem .12). Therefore, if a reasonable sampling method and GTSP solver are used, the solutions produced by Algorithm 5 will tend toward the global optimum with finer sampling. This behavior is illustrated by the top right plot in Figure 9, which shows the relative error between the cost produced via Algorithm 5 and the globally optimal cost (difference divided by the optimum) for the sampling conditions listed in the table when $\epsilon=130$ s. Notice that the relative error tends to zero with finer discretization. Also note that, for this example, the optimal solution to the discrete approximation involved vertices located on visibility region boundaries, resulting in an insensitivity to radial grid spacing. In contrast, when $\epsilon=16.26$ s, the problem becomes degenerate, since there is a single configuration to which the UAV can travel in order to satisfy the $\epsilon$ bound. In this case, Problem 2 is infeasible for any sampling scheme that does not choose the precise configuration in question. Therefore, Algorithm 5 is not resolution complete in the sense of Theorem .12 when $\epsilon=16.26$. Note, however, that resolution completeness holds if $\epsilon$ is increased by any arbitrarily small positive amount.

7 Conclusions

A novel algorithmic framework was presented for constructing unmanned aerial vehicle trajectories for surveillance of multiple targets. Adopting reasonable constraints on optimal dwell-time behavior and UAV maneuvers, the framework works to balance mission goals, namely, the closed trajectory and the initial maneuver times, while simultaneously accommodating both viewing and dwell-time constraints associated with each target. In particular, an $\epsilon$-constraint scalarization method is applied to rigorously pose the multi-objective problem as a constrained optimization, which is then approximated by a finite, path-planning problem that results from careful sampling of the UAV configuration space. Solutions to related GTSP problems can then be utilized to construct solutions to the discrete approximation which, in many cases, map directly to globally optimal solutions of the discrete problem. These solutions were then mapped to solutions of the continuous problem. It was shown that, under certain conditions, the complete heuristic procedure is resolution complete in the sense that solution quality generally increases with increasingly fine sampling.

Avenues of future research include the expansion to the multi-vehicle case, explicit comparisons with other routing schemes (e.g., Markov chain-based schemes), and an investigation of alternative discretization strategies. In addition, incorporation of uncertain dwell-times and the explicit pairing with other facets of complex missions, e.g. operator analysis of imagery, should be explored.

Appendix: Resolution Completeness of Algorithm 5

We start with some preliminary notions. For this appendix, we assume that all angles $\theta\in[0,2\pi)$ are equivalently represented as points on the unit circle $S^{1}\coloneqq\{x\in{}^{2}\;|\;\|x\|=1\}$ via the relation $\theta\mapsto(\cos(\theta),\sin(\theta))$. As such, we assume without loss of generality that the UAV configuration space is ${}^{2}\times S^{1}$. Recall that $\texttt{DWL}_{j}$ is defined as the set of configurations from which a UAV can begin executing a feasible dwell-time maneuver at $T_{j}$. We first parameterize the set of UAV routes that can be produced by Algorithm 5.

Definition .7** (Initial Maneuver Set).**

The parametrized initial maneuver set, $\texttt{INL}_{\epsilon}$, is defined

[TABLE]

where $\texttt{DIST}(v_{0},v)$ is the time required for the UAV to traverse the optimal Dubins path from $v_{0}$ to $v$.

Note that we have re-defined $\texttt{INL}_{\epsilon}$ as the continuous analog of the discrete set of Section 4.

Definition .8** (Closed Trajectory Set).**

The parametrized closed trajectory set, $\texttt{CLS}_{\epsilon}$, is defined

[TABLE]

where $\texttt{Perm}(M)$ is the set of all permutations of $\{1,\dots,M\}$.

Notice that we can map each $\textbf{v}\coloneqq(v_{1},v_{2},\ldots,v_{M})\in\texttt{CLS}_{\epsilon}$ to a feasible solution of Problem 1 by (i) appending the initial configuration $v_{0}$, (ii) constructing optimal Dubins routes between successive nodes (connecting $v_{M}$ to $v_{1}$), and (iii) appending dwell-time maneuvers. Strictly speaking, such a mapping is not unique, since, for some $j$, there may exist $v\in\texttt{DWL}_{j}$ that is the starting configuration for multiple distinct dwell-time maneuvers (associated to the same target). However, recalling Remark 4, this ambiguity is inconsequential to the present analysis and we can assume without loss of generality that, for each $j$, the correspondence between $\texttt{DWL}_{j}$ and the set of dwell-time maneuvers at $T_{j}$ is bijective. Under this assumption, any optimal solution to Problem 1 is uniquely associated with some element in the set $\texttt{CLS}_{\epsilon}$ via the mapping just described. With this in hand, we now assign an appropriate cost to each element of $\texttt{CLS}_{\epsilon}$.

Definition .9** (Cost).**

Define the map $\texttt{LGTH}:\texttt{CLS}_{\epsilon}\rightarrow{}_{\geq}$ where $\texttt{LGTH}((v_{1},v_{2},\ldots,v_{M}))$ is the time required for the UAV to traverse the closed trajectory (beginning at $v_{1}$) that sequentially touches and performs the dwell-time maneuver associated with all vertices $v_{1},\ldots,v_{M}$.

For each $T_{j}$, the length of an appropriate dwell-time maneuver varies continuously as a function of the maneuver’s starting configuration (element in $\texttt{DWL}_{j}$). As such, results in [7] imply that (i) the map LGTH is well-defined and continuous except on a finite set of $(3M-1)$-dimensional smooth surfaces embedded in $({}^{2}\times S^{1})^{M}$, and (ii) in the limit from one side of each discontinuity surface, LGTH is continuous up to and on the surface. The presence of discontinuity sets necessitates the following definition. Here, $\texttt{LGTH}^{*}\coloneqq\min\;\texttt{LGTH}$.

Definition .10** (Degeneracy).**

An instance of Problem 1 is called degenerate if, for every sequence $\textbf{v}_{1},\textbf{v}_{2},\ldots\in\texttt{CLS}_{\epsilon}$ such that $\texttt{LGTH}(\textbf{v}_{k})\to\texttt{LGTH}^{*}$ and $\textbf{v}^{*}\coloneqq\lim_{k\to\infty}\textbf{v}_{k}$ exists, the limit $\textbf{v}^{*}$ either (i) is not contained in $\texttt{CLS}_{\epsilon}$, (ii) belongs to a discontinuity set of LGTH, or (iii) has a first component (corresponding to the starting point of the closed trajectory) that is isolated in the set $\text{INL}_{\epsilon}$.

Finally, we introduce dense sampling procedures.

Definition .11** (Dense Sampling Procedure).**

Suppose that, for each $N\in\mathbb{N}$, a set $\mathcal{A}\subset{}^{2}\times S^{1}$ is sampled $N$ times to form a discrete subset $A_{N}$, i.e. $|A_{N}|=N$ for each $N$. Such a procedure is called a dense sampling procedure if, for any $a\in\mathcal{A}$ and any open neighborhood $U$ of $a$, there exists $\hat{N}\in\mathbb{N}$ such that $A_{N}\cap U$ is non-empty for all $N>\hat{N}$.

We are now ready to rigorously characterize resolution completeness of Algorithm 5.

Theorem .12** (Resolution Completeness).**

If Problem 1 is not degenerate, and the set $\texttt{INL}_{\epsilon}$ is such that: (i) $\texttt{INL}_{\epsilon}\neq\emptyset$, and (ii) there exists $\hat{\jmath}\in\{1,\dots,M\}$ such that either $\texttt{INL}_{\epsilon}\subseteq\texttt{DWL}_{\hat{\jmath}}$ or $\texttt{DWL}_{\hat{\jmath}}\subseteq\texttt{INL}_{\epsilon}$, then Algorithm 5 is resolution complete in the following sense:

Form a sequence $\cramped{\{(\texttt{INL\_MNVR}_{N},\texttt{CLS\_TRAJ}_{N})\}_{N\in\mathbb{N}}}$, where each element $(\texttt{INL\_MNVR}_{N},\texttt{CLS\_TRAJ}_{N})$ represents the output of a call to Algorithm 5 when: (i) the number of discrete samples at each target is $N$, (ii) $\texttt{INL}_{\epsilon}^{*}$ (Algorithm 4, line $3$) is chosen to satisfy the conditions of Theorem 5.4, and (iii) an optimal GTSP solution is found (Algorithm 4, line $5$). Then, $(\texttt{INL\_MNVR}_{N},\texttt{CLS\_TRAJ}_{N})$ is a feasible solution to Problem 1 for each $N$, and if a dense sampling procedure is used to generate the discrete node sets at each target, then the length of the tours in the sequence $\{\texttt{CLS\_TRAJ}_{N}\}$ approaches the length of an optimal solution to Problem 1 as $N\to\infty$.

Proof .13**.**

First note that the notion of resolution completeness described here is well-defined under the specified assumptions. Indeed, if the (non-discrete) set $\texttt{INL}_{\epsilon}$ satisfies the necessary assumptions, then the discrete analog always has the properties required for Theorem 5.4, and $\texttt{INL}_{\epsilon}^{*}$ can be chosen to exploit the GTSP equivalence when constructing the sequence $\{(\texttt{INL\_MNVR}_{N},\texttt{CLS\_TRAJ}_{N})\}_{N\in\mathbb{N}}$.

Feasibility of each $\texttt{INL\_MNVR}_{N}$ follows from Theorem 1, and the fact that feasible solutions to Problem 2 map to feasible solutions to Problem 1. Let $\textbf{v}_{1},\textbf{v}_{2},\ldots\in\texttt{CLS}_{\epsilon}$ be a sequence such that $\texttt{LGTH}(\textbf{v}_{i})\to\texttt{LGTH}^{*}$. Note that $\texttt{CLS}_{\epsilon}\subseteq({}^{2}\times S^{1})^{M}$ can be represented as a bounded subset in $({}^{2}\times S^{1})^{M}\subseteq{}^{4M}$. Invoking the Bolzano-Weierstrass theorem and compactness of $S^{1}$, any infinite sequence in $\texttt{CLS}_{\epsilon}$ contains a subsequence that converges to some point in $({}^{2}\times S^{1})^{M}$. Thus, recalling that Problem 1 is non-degenerate, we can assume without loss of generality that $\textbf{v}^{*}=\lim_{i\to\infty}\textbf{v}_{i}$ exists, is contained in $\texttt{CLS}_{\epsilon}$, does not belong to a discontinuity sets of LGTH, and has a first component that is not isolated in $\texttt{INL}_{\epsilon}$. Since the function LGTH is continuous at the point $\textbf{v}^{*}$, for any $\delta>0$, there exists an open neighborhood $U\subset\{(\texttt{DWL}_{\sigma(1)}\cap\texttt{INL}_{\epsilon})\times\cdots\times\texttt{DWL}_{\sigma(M)}\;|\;\sigma\text{ permutes the set }\{1,\dots,M\}\}$ of $\textbf{v}^{*}$, within which LGTH takes values within $\delta$ of $\texttt{LGTH}^{*}$. Since the sampling procedure is dense, for some $\hat{N}$, there will be a discrete node placed inside of the set $U$ for all $N\geq\hat{N}$. Since the GTSP is solved exactly, it follows that $|\texttt{LGTH}(\texttt{CLS\_TRAJ}_{N})-\texttt{LGTH}^{*}|<\delta$ for $N>\hat{N}$, proving convergence.

Theorem .12 states that, under a non-degeneracy assumption about the problem structure, reasonable implementations of Algorithm 5 will produce solutions that approximate optimal solutions to the full, multi-objective routing problem. Indeed, the resolution completeness property ensures that Problem 2 will be appropriately reflective of Problem 1.

A few final comments are in order. First, note that degenerate problems occur only in very select situations: a degenerate problem instance is made non-degenerate by perturbing target locations by an arbitrarily small amount (example degenerate problem instances for particular cases are shown in [7]). Thus, the non-degeneracy condition in Theorem .12 is not typically restrictive. Second, if $\texttt{INL}_{\epsilon}$ does not satisfy the conditions of Theorem .12, then resolution completeness does not hold in general since Algorithm 4 is not guaranteed to find an optimal solution to Problem 2 (Theorem 5.4). Nevertheless, the quality of solutions produced by Algorithm 5 generally increase as sampling granularity is made increasingly fine. Finally, it is not generally possible to find optimal solutions to GTSPs. Therefore, the utility of Theorem .12 is its ability to provide intuition about qualitative solution behavior, and to ensure that the discrete method herein is an appropriate approximation.

Acknowledgments

The authors thank Karl J. Obermeyer of Obermeyer Labs for his insight regarding vehicle routing.