Visual Monitoring for Multiple Points of Interest on a 2.5D Terrain using a UAV with Limited Field-of-View Constraint

TL;DR

This paper addresses the challenge of visual monitoring of multiple points on a 2.5D terrain using a UAV with limited field-of-view, proposing an approximation algorithm for planning UAV tours.

Contribution

It introduces a two-phase strategy with a constant-factor approximation algorithm for the TSPN, reducing it to a GTSP and solving via ILP, for efficient UAV monitoring.

Findings

The proposed algorithm provides near-optimal tours in varied scenarios.

Comparative evaluation shows effectiveness of the ILP and GTSP approaches.

Preliminary field experiments validate the practical applicability.

Abstract

Varying terrain conditions and limited field-of-view restricts the visibility of aerial robots while performing visual monitoring operations. In this paper, we study the multi-point monitoring problem on a 2.5D terrain using an unmanned aerial vehicle (UAV) with limited camera field-of-view. This problem is NP-Hard and hence we develop a two phase strategy to compute an approximate tour for the UAV. In the first phase, visibility regions on the flight plane are determined for each point of interest. In the second phase, a tour for the UAV to visit each visibility region is computed by casting the problem as an instance of the Traveling Salesman Problem with Neighbourhoods (TSPN). We design a constant-factor approximation algorithm for the TSPN instance. Further, we reduce the TSPN instance to an instance of the Generalized Traveling Salesman Problem (GTSP) and devise an ILP formulation to solve it. We present a comparative evaluation of solutions computed using a branch-and-cut implementation and an off-the-shelf GTSP tool -- GLNS, while varying the points of interest density, sampling resolution and camera field-of-view. We also show results from preliminary field experiments.

Tables

Table 1. TABLE I : The table shows the number of instances solved optimally by the ILP solver. Numbers in bracket represent the number of instances for which the solver could compute at least a feasible solution.

Number of instances solved by the ILP solver.
$m$	$\delta =20°$			$\delta =30°$			$\delta =40°$
	d = 20	d = 30	d = 40	d = 20	d = 30	d = 40	d = 20	d = 30	d = 40
4	20 (20)	20 (20)	9 (20)	20 (20)	19 (20)	7 (20)	20 (20)	18 (20)	13 (20)
6	14 (20)	0 (20)	0 (20)	10 (20)	0 (20)	0 (19)	7 (20)	0 (20)	0 (19)
8	0 (20)	0 (20)	0 (18)	0 (19)	0 (17)	0 (8)	0 (18)	0 (13)	0 (8)

Table 2. TABLE II : Table shows the mean percentage relative gap of GLNS solutions w.r.t. ILP solver generated lower bounds. The numbers in bracket represent the standard deviation.

Mean percentage relative gap of GLNS solver solution w.r.t. ILP generated lower bounds.
$m$	$\delta =20°$			$\delta =30°$			$\delta =40°$
	d = 20	d = 30	d = 40	d = 20	d = 30	d = 40	d = 20	d = 30	d = 40
4	0 (0)	0 (0)	4.8 (5.6)	0 (0)	1 (4.5)	12.4 (16.7)	0 (0)	2.1 (7.4)	14 (24.7)
6	4.5 (9.5)	24.7 (11.2)	30.1 (12.2)	17.8 (23.1)	44.8 (16.8)	52.3 (18.2)	28.4 (27.9)	57.3 (18)	68 (19.8)
8	31.6 (15.2)	41.7 (17.8)	48.2 (23.2)	53.1 (20.8)	68.2 (20)	85.1 (20.2)	76.2 (19.7)	85 (16.1)	94 (10.2)

Table 3. TABLE III : The table shows the relative gap for best ILP solutions w.r.t. solver generated lower bound within a maximum time limit of 900 seconds. The numbers in bracket represent the standard deviation.

Mean percentage relative gap of ILP generated solutions.
$m$	$\delta =20°$			$\delta =30°$			$\delta =40°$
	d = 20	d = 30	d = 40	d = 20	d = 30	d = 40	d = 20	d = 30	d = 40
4	0 (0)	0 (0)	5.3 (5.9)	0 (0)	1 (4.7)	13.0 (16.9)	0 (0)	2.2 (7.7)	14.2 (24.9)
6	4.7 (9.8)	25.9 (10.7)	32 (11.4)	18.2 (23.4)	45.3 (16.8)	51.3 (14.5)	28.9 (28.1)	58.5 (17.8)	68.5 (18.7)
8	33.7 (14.3)	43.7 (17.6)	45.3 (14.8)	52.3 (18.2)	65.1 (16.3)	65.7 (13)	74.99 (18.4)	79.6 (14.3)	87.5 (9.3)

Equations

$$L\ge \frac{m}{2}\alpha r.$$

$$L_{D_I}^{\ast }\ge \frac{m}{2}\alpha h\tan \delta .$$

$$L_D^{\ast }\ge L_{D_I}^{\ast }\ge \frac{m}{2}\alpha h\tan \delta .$$

$$L_d\le (1+ϵ)(L_V^{\ast }+2m_d(h-l)\tan \delta )+2m_d\pi (h-l)\tan \delta ),$$

$$L_d$$

$$\le (1+ϵ)(L_V^{\ast }+2m_d(h-l)\tan \delta )+2m_d\pi (h-l)\tan \delta .$$

$$m_d\le {\left(\frac{2h}{h-l}\right)}^{2}m.$$

$$L_d\le \left((1+ϵ)\left(1+16\frac{h}{\alpha (h-l)}{L}_V^{\ast }\right)+16\pi \frac{h}{\alpha (h-l)}\right)L_V^{\ast },$$

$$L_d$$

$$\le (1+ϵ)\left(L_V^{\ast }+2{\left(\frac{2h}{h-l}\right)}^{2}m(h-l)\tan \delta )\right)$$

$$+2{\left(\frac{2h}{h-l}\right)}^{2}m\pi (h-l)\tan \delta ,$$

$$\le (1+ϵ)\left(L_V^{\ast }+2{\left(\frac{2h}{h-l}\right)}^2\frac{2}{\alpha h\tan \delta }{L}_V^{\ast }(h-l)\tan \delta )\right)$$

$$+2{\left(\frac{2h}{h-l}\right)}^2\frac{2}{\alpha h\tan \delta }{L}_V^{\ast }\pi (h-l)\tan \delta ,$$

$$\le \left((1+ϵ)\left(1+16\frac{h}{\alpha (h-l)}\right)+16\pi \frac{h}{\alpha (h-l)}\right)L_V^{\ast }$$

$$\le O(1)L_V^{\ast }.$$

$$\min \sum _{v_i\in \mathcal{V}}\sum _{v_j\in \mathcal{V}}{c}_{ij}{y}_{ij}$$

$$\sum _{v_j\in \mathcal{V}\setminus v_i}{y}_{ji}-\sum _{v_j\in \mathcal{V}\setminus v_i}{y}_{ij}=0,\forall v_i\in \mathcal{V}$$

$$\sum _{v_i\in S_k}\sum _{v_j\in \mathcal{V}\setminus S_k}{y}_{ji}=1,\forall k\in [1\dots m]$$

$$\sum _{(i,j)\in {\delta }^{+}(s)}{y}_{ij}=1,\forall s\in \mathbb{P}(\mathcal{S})\setminus \{\mathcal{S},\varphi \}$$

$$y_{ij}\in \{0,1\}\forall v_i,v_j\in \mathcal{V}$$

$$\sum _{(i,j)\in {\delta }^{+}(\kappa )}{y}_{ij}=1.$$

Full text

I Introduction

Visual surveillance and monitoring is an important application area for aerial robots. Crop management [1], area coverage [2, 3], terrain mapping [4], structural inspection [5], and disaster management [6, 7] are some applications where aerial robots are widely used. When planning paths in such missions, it is important to take visibility obstructions into account. Landscape features such as mountains, gorges, buildings, and bridges limit the line-of-sight of the aerial robots. In addition, operative limitations such as camera field-of-view and maximum flight altitude corresponding to the image resolution and/or regulatory requirements also restrict visibility. It is imperative that such restrictions be accounted for when planning for monitoring missions.

In this paper, we address the visual monitoring problem on 2.5D terrains using a UAV while accounting for camera field-of-view and terrain imposed visibility restrictions. We present a two-phase strategy to compute a tour for an aerial robot to visually monitor a set of points located within a terrain. A naive strategy is to visit each point of interest (or a point directly above it). This strategy does not exploit the camera field-of-view and essentially assumes the most restrictive field-of-view only along the center of the camera. When considering the flight altitude, this can lead to solutions that are numerically far from the optimal by a factor of $R$, where $R$ is the radius of the camera footprint on the ground (assuming a circular camera footprint). In fact, in the special case that the field-of-view angle of the camera sensor tends to zero, the problem reduces to an instance of the Traveling Salesman Problem (TSP) that is known to be NP-hard [8]. Since our problem is a generalization of TSP, it is at least as hard as TSP.

II Related Work

Visual monitoring and surveillance using UAVs has been an active area of research over the past decade. Various lines of work have addressed aspects related to target monitoring with differential priorities [9], multi-robot surveillance [1, 3, 10, 11], persistent monitoring [10, 12, 3], mission planning while addressing robot kinematics [13, 14] and so on. There has been limited effort towards addressing the visual monitoring problem in the presence of visibility restrictions due to terrain features. Terrain visibility however is a classic problem in computational geometry [15] and graphics literature ([16] and references within). Terrain guarding [15], and watchtower problems [17] relate to computing a set of points that lie on the terrain and at an altitude, respectively, to ensure line-of-sight area coverage on terrains. The problem of line-of-sight coverage within a polygon by a watchman (or a robot), known as the Watchman Routing Problem (WRP), is also well studied in the literature ([18] and the references within). Variations of WRP, including the homogeneous [19, 20, 21] and heterogeneous [11, 10] multiple-robot versions (both within restricted sub-domains) have also been studied extensively. Recently, Maini et. al. [11, 10] modeled the coverage problem for piece-wise linear features within terrains as a variant of the n-WRP.

There is a lot of literature within the aerial robotics community on coverage path planning. Area decomposition based on camera footprint and/or obstacle-free space is a popular choice and admits a robust discretization of the area of interest [2, 22, 14, 23, 24]. Other techniques include seed-spreader algorithms [14], potential fields [22] and graph-based search algorithms [23]. However, most of the existing works on coverage path-planning assume a flat surface and do not account for altitude variance (and hence the visibility obstructions) of the ground surface. A closely related work is that of Choi et. al. [24], who address a constant resolution coverage problem that takes into account camera viewing direction and altitude to maintain the image resolution. In this work, we address a multiple-point monitoring problem using an aerial robot while explicitly accounting for visibility restrictions due to the shape of the terrain and camera field of view. Main contributions of this paper are as follows:

• •

Extraction of visibility regions and modeling the path-planning problem as an instance of TSPN

• •

Design of a constant-factor approximation algorithm to solve the class of TSPN instances encountered within the path planning problem

• •

Reduction of the TSPN instance to GTSP thus allowing the application of existing algorithmic tools for GTSP. A new ILP formulation and a branch-and-cut implementation to solve GTSP

• •

Validation of the developed methods in simulation and field experiments

The rest of the paper is organized as follows. The application scenario and a formal problem definition are developed in Section III. In Section IV, we describe a method to compute visibility regions on the constant altitude flight plane for the points of interest. In Section V we develop a constant-factor approximation algorithm for the class of TSPN instances encountered within the path planning problem. Section VI outlines an ILP formulation and a branch-and-cut implementation to solve GTSP. Sections VII and VIII discuss evaluation results in simulation and field trials, respectively. Concluding remarks and future directions are identified in Section IX.

III Problem Formulation

Consider an environment $\mathcal{E}$ and a set of points of interest, $\mathcal{P}=\{p_{i}:i\in[1,m]\}$, as shown in Figure 1(a). We represent the topographical surface within $\mathcal{E}$ using a polyhedral (or 2.5D) terrain ([25], pg 352) and model it as a triangular irregular network (TIN). Let $h$ be the constant flight altitude for UAV operation (higher than the altitude of all terrain features). We assume that the UAV is equipped with a fixed downward-facing camera having a constant focal length and a circular field-of-view (FOV). The camera casts a conical field of view on the terrain, as shown in Figure 1(b).

While our strategies are extensible to the case of varying flight altitude, as may be required for constant resolution flights, we assume a constant flight altitude for ease of exposition. We briefly discuss the extension to constant resolution (varying altitude) flight operations, later in the text (Section V). The circular FOV is also a non-binding assumption and can be easily extended to a non-circular FOV. The FOV may then be represented as a fixed view angle, $\delta$, in each direction.

We introduce the following UAV routing problem on terrains (URPT), defined as: Given an environment with a polyhedral terrain, a set of points $\mathcal{P}$ located on the terrain, a UAV comprising of a fixed down-facing camera with a view angle $\delta$ operating at a fixed altitude $h$ above ground level, plan a minimum length tour for the UAV to visually monitor all points in the set $\mathcal{P}$.

We develop a two phase solution strategy for URPT. The first phase computes a visibility region for each point of interest, $p_{i}\in\mathcal{P}$. Visibility region of a terrain point is a closed connected space on the constant altitude flight plane, within the aerial robot’s operational region that allows the robot to monitor the corresponding point and may assume a complex geometric shape due to obstructive features within the terrain (Figure 2(b)). We outline a method to compute the visibility region for each point of interest in Section IV. Second phase involves route planning for the aerial robot to visit each of the visibility regions to complete a monitoring mission. Sections V and VI develop route planning methods for the aerial robot.

IV Visibility Computation

Consider a point, $p_{i}\in\mathcal{P}$, that needs to be monitored as shown in Figure 2(a). To compute visibility region for $p_{i}$ on the terrain, we place a viewpoint at $p_{i}$ and compute terrain visibility. The farthest visible point in a given radial direction, that obstructs all points beyond itself when viewed from a specific viewpoint is called global horizon point [16] and is expressed in terms of the elevation angle $\epsilon$. The locus of all global horizon points forms the global horizon. Visibility angle $\upsilon$ in a radial direction is computed as the complement of the elevation angle of the global horizon point. The minimum of visibility angle and camera view angle defines the boundary of the visibility region in a radial direction.

Horizon computation is a fairly well studied problem in the graphics literature ([16] and references within). We employ the approximate horizon computation method developed by Stewart [26]. However, other methods in the literature may also be used since our solution approach is independent of the horizon computation algorithm used. Our selection of Stewart’s algorithm was motivated due to its ease of implementation and computational tractability. The interested reader may refer [26] for details on the algorithm.

We use Stewart’s algorithm on the terrain shown in Figure 1(a) to compute the horizon at each point of interest. The radial space is sampled in a discrete number of values, $d$111Sampling resolution, $d$, is a user-input parameter and relates to computational complexity of the approach. Its effects are discussed in more detail in Sections VI and VII-B. and the minimum of $\upsilon$ and $\delta$ is computed in each direction. Linear interpolation of the extended projections in each direction on the constant altitude plane at height $h$ is then used to compute the boundary of the visibility region as shown in Figure 2(b).

V Polynomial-Time Approximation Algorithm for Route Planning

The route planning problem is a generalization of the NP-hard Traveling Salesperson Problem  [8]. As a result, finding the optimal solution in polynomial time is impossible (unless $P=NP$). In this section, we present a polynomial time approximation algorithm for route planning. Specifically, we present a polynomial-time algorithm that finds a tour for the UAV whose length is guaranteed to be within a constant-factor of the minimum length.

The input to our algorithm is the set of visibility regions that are computed using the method described in Section IV. The problem of finding the shortest tour that visits a set of 2D regions is known as the TSP with neighborhoods (TSPN). The neighborhoods correspond to the visibility regions, in our case. TSPN is NP-hard. However, there exists polynomial-time approximation algorithms for many special cases such as when the neighborhoods are all disks of the same radii [27] and non-overlapping convex polygons [28]. These regions may not necessarily be polygonal (may contain circular arcs) or convex and can be overlapping. Nevertheless, we show how to approximate the visibility regions by possibly-overlapping disks of the same radius. We then show that this approximation still yields a tour whose length is bounded with respect to the optimal.

In the following, let $V=\{V_{i}:i\in[1,m]\}$ be the set of input visibility regions corresponding to the points of interest $\mathcal{P}$ that the robot must monitor.

Lower Bound

We start by showing a lower bound on the length of the optimal (unknown) tour. Recall that $h$ is the height of the fixed-altitude plane on which the robot is allowed to fly and $\delta$ is the FOV angle. We construct a lower-bound approximation tour for the optimal one as follows. Replace each $V_{i}$ by a disk, say $D_{i}$, whose radius is equal to $h\tan\delta$. Let $D$ denote the collection of all the disks $D_{i}$. The disk, $D_{i}$ lies in the constant altitude flight plane at the height $h$ and centered at the same $x$ and $y$ coordinates as that of $p_{i}$.

Lemma 1

The visibility region $V_{i}$ is completely contained within the disk $D_{i}$.

Proof:

Recall that $V_{i}$ is obtained by projecting a reverse cone whose apex is at $p_{i}$ on the fixed altitude plane at height $h$. Let the coordinates of $p_{i}$ be $(x_{i},y_{i},z_{i})$. Consider a reverse cone drawn centered at $(x_{i},y_{i},0)$. Further assume that this cone is not obstructed by any point on the terrain. It is clear that this cone completely contains the cone drawn at $p_{i}$. The intersection of the larger cone with the fixed altitude plane at height $h$ yields the disk $D_{i}$. Therefore, $V_{i}$ is completely contained within $D_{i}$. (In the extreme case, $V_{i}$ is the same as $D_{i}$.) ∎

Lemma 2

Let $L_{V}^{*}$ be the length of the optimal tour that visits at least one point in each visibility region, $V_{i}\in V$. Let $L_{D}^{*}$ be the length of the optimal tour that visits at least one point in each disk, $D_{i}\in D$. We have: $L_{D}^{*}\leq L_{V}^{*}$.

Proof:

From Lemma 1, we know that $V_{i}\subseteq D_{i}$. Therefore, any tour that visits at least one point in each $V_{i}$ is also a tour that visits at least one point in each $D_{i}$. As a result, the optimal tour (of length $L_{V}^{*}$) that visits at least one point in each $V_{i}$ is also a tour that visits at least one point in each $D_{i}$. However, $L_{D}^{*}$ is the length of the optimal tour that visits at least one point in each $D_{i}$. Therefore, $L_{D}^{*}\leq L_{V}^{*}$. ∎

Finally, we lower bound the length of the optimal tour that visits at least one point in each $V_{i}$. We relate the lower bound to the maximum number of non-overlapping disks. Specifically, let $D^{I}\subseteq D$ be the largest set of disks, $D_{i}$, such that no two disks overlap with each other. This can be found out greedily by constructing the maximum independent set of the disks, as shown in [27]. Let $m$ be the number of disks in $D^{I}$.

Lemma 3

Let $L_{V}^{*}$ be the length of the optimal tour that visits at least one point in each visibility region, $V_{i}$. We have: $L_{V}^{*}\geq\dfrac{m}{2}\alpha h\tan\delta$ where $\alpha=0.4786$ and when $m\geq 3$.

Proof:

From Theorem 1 in [29], we know that any tour of length $L$ that visits at least one point in $m$ disjoint disks of radius $r$ satisfies,

[TABLE]

Therefore, the optimal tour of length $L_{D_{I}}^{*}$ that visits all the $m$ disks in $D^{I}$ must satisfy:

[TABLE]

Since $D^{I}\subseteq D$, the optimal tour that visits at least one point in each disk in $D$ will have a length:

[TABLE]

From the above equation and Lemma 2, we get the desired inequality: $L_{V}^{*}\geq\dfrac{m}{2}\alpha h\tan\delta$. ∎

Upper Bound

So far, we have only presented a lower bound on the length of the shortest tour that visits each disk $D_{i}$. Note that visiting each disk in $D_{i}$ is necessary to visit $V_{i}$ but may not be sufficient. Instead, we will replace each $V_{i}$ by an inner disk, say $d_{i}$, such that it is completely contained within $V_{i}$, i.e., $d_{i}\subseteq V_{i}$. The inner disk $d_{i}$ is also centered at the same point as $D_{i}$ and all inner disks have the same radius.

Let $l$ be the maximum height of the terrain. That is, all points on the terrain are at height of $l$ or below. It is easy to see that $l<h$. We set the radius of the inner disks to be equal to $(h-l)\tan\delta$. Using a similar argument as given in Lemma 1, we can prove that $d_{i}\subseteq V_{i}$. That is, all inner disks are completely contained within the visibility regions, $V_{i}$.

Our algorithm for solving the route planning problem is to find a tour that visits at least one point in each inner disk, $d_{i}$. This can be found using the algorithm presented by Dumitrescu and Mitchell [27]. (1) Find the maximum independent set of non-overlapping inner disks, say $d^{I}$. (2) Find a $(1+\epsilon)$–approximation to the optimal TSP tour that visits the center of all disks in $d^{I}$. (3) Follow the tour found in the second step. Every time the tour enters a new disk, take a detour to follow the circumference till you reach the same point again, and then move towards the center. Note that this step adds a detour of length at most $2\pi r$ to the TSP tour, where $r$ is the radius of the disk.

Let $L_{d}$ be the length of this tour. We now provide an upper bound to the length of this tour.

Lemma 4

Let $L_{d}$ be the length of the tour found using the proposed algorithm. We have:

[TABLE]

where $m_{d}$ is the maximum number of non-overlapping inner disks, $d^{I}$.

Proof:

The length of the tour, $L_{d}$, is equal to the distance traveled to visit the centers of the disk in $d^{I}$ (Step 2) and the detours added every the center is visited (Step 3). Let $L_{TSP}^{*}$ be the length of the optimal TSP tour that visits the center of the disks in $d^{I}$. Although finding $L_{TSP}^{*}$ is NP-hard, there exists polynomial time approximation algorithms that find a tour whose length is at most $(1+\epsilon)L_{TSP}^{*}$ for any $\epsilon>0$. Therefore,

[TABLE]

The second inequality follows from the fact that we can always construct a tour that visits the center of the disks in $d^{I}$ by first finding the optimal tour that visits at least one point in each disk in $d^{I}$ (of length $L_{d}^{*}$) and then adding a detour of at most $2r$ to visit the center. That is, $L_{TSP}^{*}\leq L_{d}^{*}+2m_{d}$ and $L_{d}^{*}\leq L_{V}^{*}$. ∎

What remains to show is the relationship between $m$ and $m_{d}$. It is easy to see that $m_{d}\geq m$, that is the number of non-overlapping outer disks (in $D$) cannot be more than the number of non-overlapping inner disks ($d_{i}$). We will show that $m_{d}$ cannot be arbitrarily larger than $m$.

Lemma 5

Let $m$ be the maximum number of non-overlapping outer disks, $D_{i}$. Let $m_{d}$ be the maximum number of non-overlapping inner disks, $d_{i}$. We have:

[TABLE]

Proof:

Consider an outer disk, $D_{i}$, whose radius is equal to $h\tan\delta$. Draw another disk, say $D_{i}^{\prime}$, whose radius is equal to $2h\tan\delta$ with the same center. Any inner disk that intersects with $D_{i}$ is completely contained within $D_{i}^{\prime}$. We now bound the maximum number of inner disks that can be packed within $D_{i}^{\prime}$ without any two overlapping. One inner disk has an area of $\pi(h-l)^{2}\tan^{2}\delta$. Therefore, at most $\left(\frac{2h}{h-l}\right)^{2}$ non-overlapping inner disks are contained within $D_{i}^{\prime}$. Since there are $m$ non-overlapping outer disks, we get the desired result. ∎

We are now ready to state the main result of this section.

Theorem 1

Let $L_{d}$ be the length of the tour found using the proposed algorithm. Let $L_{V}^{*}$ be the length of the optimal tour that visits at least one point in each visibility region, $V_{i}$. We have:

[TABLE]

where $h$ is the height of the fixed-altitude plane, $l$ is the height of the tallest point on the terrain, and $\alpha=0.4786$.

Proof:

We know from Lemma 4,

[TABLE]

∎

This shows that the proposed algorithm yields a constant-factor approximation. The constant depends on two parameters, maximum height of the terrain and the height of the fixed-altitude plane, but is otherwise independent of the input (e.g., $|\mathcal{P}|$, the width of the terrain, etc.). We would like to remark that, the same approximation algorithm may also be used in the case of constant resolution imagery (variable flight altitude) missions to compute UAV tours within a constant-factor of the optimal. It is easy to see that the enclosing outer disks ($D_{i}$) and enclosed inner disks ($d_{i}$), used to compute the lower and upper bounds on the UAV tour respectively, are still valid and may be used to compute the same constant approximation factor.

VI Route Planning

The visibility regions for all points of interest computed as discussed in Section IV are given as input to the route planning stage. Unless all visibility regions are contained within one of the visibility region, in which case the containing region may as well be ignored and the problem be solved for $m-1$ visibility regions, the tour for the aerial robot must enter each visibility region at a point on the boundary of the region. This implies that we can restrict the search for the points visited by the aerial robot to the boundaries of the respective visibility regions. Similar ideas have been used by Obermeyer et al. [4] for path planning for a non-holonomic robot through a set of polygonal spaces. Hence, we consider only the points on the boundary of each visibility region to compute a tour for the aerial robot. Each visibility region contributes $d$ (sampling parameter, refer Section IV) unique points on its boundary. To address the overlapping regions case, we duplicate all points in the overlapping regions and add them to each region in the intersection.

To formalize, let $S_{i}$ be the set of sample points on the boundary of the $i^{th}$ visibility region corresponding to the $p_{i}$ point of interest and $\mathcal{S}=\{S_{i}:i\in[1,m]\}$ be the set of all such sets. Let $\mathcal{V}=\bigcup\limits_{i=1}^{m}S_{i}$ be the set of all vertices on boundary of a visibility region. $\mathcal{V}$ also includes duplicate points that lie in the intersection of visibility regions. We define the cost function, $c_{ij}:\mathcal{V}\times\mathcal{V}\rightarrow\mathbb{R}^{+}$, as the length of path for the aerial robot to go from $v_{i}$ to $v_{j}$, where $v_{i},v_{j}\in\mathcal{V}$. In this form, the problem reduces to an instance of the well-known Generalized Traveling Salesman Problem (GTSP). We employ two strategies to solve the GTSP instance. An Integer Linear Programming (ILP) formulation solved within a branch-and-cut framework and a specialized GTSP solver called GLNS [30]. We give the ILP formulation below and describe the GLNS solver settings in the next Section.

ILP Formulation

To formulate the problem as an Integer Linear Program, we define binary decision variable, $y_{ij}$, for each pair of vertices $v_{i}$ and $v_{j}$ in the set $\mathcal{V}$. $y_{ij}=1$ if the aerial robot visits $v_{i}$ and $v_{j}$ vertices in order. Let $\delta^{+}(X)$ denote the set of pairs $(i,j)$ such that $v_{i}\in X$ and $v_{j}\in\mathcal{V}\setminus{X}$ and $\mathbb{P}(X)$ denote the power set of $X$. The objective function and constraints of the ILP formulation are defined as

Objective:

[TABLE]

Degree Constraints:

[TABLE]

Sub-tour Elimination Constraints:

[TABLE]

Variable Domain:

[TABLE]

Equations (9) and (10) represent tour constraints and ensure each visibility region is visited. Equation (11) represents the set of sub-tour elimination constraints. The number of sub-tour elimination constraints grows exponentially with increase in the number of visibility regions. Therefore, we employ a branch-and-cut strategy to solve the ILP formulation. A relaxed formulation, minus the sub-tour elimination constraints is given as input to the solver. A separation algorithm (Algorithm 1) computes valid inequalities (given by Equation (13)) at runtime and adds them to the formulation to ensure feasibility of the final solution. Branch-and-cut has been observed to be an effective strategy to improve computational time in problems of similar flavor [31].

[TABLE]

VII Simulation Results

The performance of the two stage strategy is evaluated using IBM ILOG CPLEX library (version 12.7) in C++11 and GLNS solver in Julia [30]. For visualization, we use MATLAB R2017a with TIN based modeling of the terrain.

VII-A Simulation setup

To generate the simulation instances we use an environment of size 200 units$\times$200 units. A 10$\times$10 grid was superimposed on the environment and terrain altitude at each grid point was sampled randomly between 0 to 100 units. A TIN representation of the terrain was then generated by a piecewise triangular interpolation between neighboring grid points. 20 different terrains were generated using this method. The aerial robot flight altitude was fixed at 125 units (clear of all terrain features). The size of the set $\mathcal{P}$ of points to be monitored on the terrain was varied from 4 to 8 in steps of 2. The camera FOV view angle $\delta$ was varied from 20$\degree$ to 40$\degree$ in steps of 10. The value of sampling resolution parameter $d$, that determines the number of points sampled on the boundary of the terrain, was varied from 20 to 40 in steps of 10. A total of 540 instances were generated.

The total number of points in the GTSP instances were in the range of 80 to 1500.In general, the GTSP instance size increases with increase in number of points (sets) and the value of sampling parameter $d$. It also increases with increase in camera view angle, $\delta$; as this increases the overlap between visibility regions. The GLNS solver was used in the fast mode setting and allowed a maximum time of 100 seconds. The maximum time limit was not reached for any of the instances with GLNS. The ILP formulation was implemented in a branch-and-cut framework using the lazy callback functionality of IBM ILOG CPLEX library. The solver was allowed to run for a maximum time of 900 seconds for each instance.

VII-B Results

Simulation results using GLNS and ILP solvers for GTSP instances generated using visibility regions to compute tours for the aerial robot are summarized in Tables II-III. GLNS solver found a feasible solution for each instance in under 5 seconds. Table I gives the number of instances solved by the ILP solver in 900 seconds. The ILP solver was not able to find any feasible solution for a sizable number of instances. This is attributed to the time limit of 900 seconds imposed on the ILP solver. We do not report results for larger instances ($m=10$, $d=50$) for both ILP and GLNS, as the solver could not find the optimal solution for any instance with 8 sets ($m=8$). The relative gap for both GLNS and ILP solver generated solutions rises quickly with increase in the value of $m$ and $\delta$ (Tables II and III). This points to the hardness and large size of the instances, as GLNS is a widely used tool to solve GTSP instances. Sample paths generated using the two solution methods are shown in Figure 3.

VIII Field Experiments

Field experiments were conducted at IIIT-Delhi campus to validate the proposed solution approach. The operational area for the experiments is shown in Figure 4(a). A DJI Phantom 4 quadrotor was used to perform the experiments. A DEM of the area was created using PIX4D and modeled using a TIN representation. Four target points were placed in the area as shown in Figure 4(b). Visibility regions for each target were computed using the visibility computation strategy given in Section IV (Figure 5(a)). The value of the sampling resolution parameter $d$ was set to 30. Paths for the UAV were computed as solution to the corresponding GTSP instance solved using GLNS solver, as shown in Figure 5(b). Experiment footage showing the UAV flight and camera imagery may be viewed in the video attachment to the paper.

IX Conclusion and Future Work

In this paper, we proposed a two phase strategy to computer tours for a UAV to perform multi-point visual monitoring on a terrain. The path planning problem was modeled as an instance of the TSPN problem and a constant-factor polynomial time approximation algorithm was designed for the class of TSPN instances. Further, GTSP based solution methods by discretizing the visiblity region boundaries were evaluated using two solution techniques – ILP formulation implemented in a branch-and-cut framework and GLNS (widely used GTSP solver). Field experiments were also conducted to verify the applicability and effectiveness of solution methods.

A natural extension of the proposed framework is to increase the number of vehicles and optimize the use of vehicles for monitoring larger number of points of interest. Additional venues for future exploration can be (i) design of efficient heuristics to decrease the tour cost (ii) persistent monitoring problem wherein fuel limitations of the aerial robot need to be addressed and (iii) exploring regions that have non-convex terrains.