Approximation Algorithm for a Heterogeneous Vehicle Routing Problem

1. Introduction

Recent advances in small sensing devices and unmanned vehicles (UVs) have provided an efficient way of collecting useful information in civil and military applications such as crop monitoring [1, 2], forest temperature monitoring [3] and border surveillance [4]. These applications frequently require the UVs to travel large swathes of land, and collect data both from onboard sensors and small sensing devices (also referred to motes) deployed in the ground. This article develops novel algorithms required for efficiently deploying heterogeneous UV networks in monitoring applications.

This article addresses the following resource allocation problem called the multiple heterogeneous UV problem (MHUVP) that arises in monitoring applications: Given a team of small heterogeneous UVs located at distinct depots (or initial locations), a set of target locations and the cost of traveling between any two locations for each UV, the problem aims to find a sequence of targets for each UV such that each target is visited by a UV, the UVs return to their respective depots after visiting the targets, and the sum of the travel cost of all the UVs is minimized. Determining the sequence of targets to visit for an UV will in turn specify the path for the UV. The travel cost of an UV is defined as the length of the path traveled by the UV. A collection of UVs is heterogeneous if the cost of traveling between any two target locations depend on the type of the UV. The travel costs for each UV are assumed to be asymmetric, i.e., the cost of traveling from location A to location B may be different from the cost of traveling from location B to location A. Heterogeneity and the asymmetric travel costs are common in surveillance applications [5] where the UVs have a bound on its minimum turning radius [6].

In the simplest of the settings, where there is only one UV, this problem is a generalization of the traveling salesman problem (TSP) and is NP-Hard [7]. In the presence of multiple UVs [8], there are two main issues that complicate this problem. First, the assignment of targets to the UVs is unknown. Second, the UVs are heterogeneous. Even though special cases of this problem have recently received attention in the literature, most of the known work has focused on multiple-homogenous UV problems [9, 6, 10, 11, 12] or problems with symmetric travel costs [13, 14, 15, 16]. Approximation algorithms for multiple vehicles with symmetric, heterogeneous costs were presented in [14, 16, 17]. Currently, there are no algorithms with provable performance guarantees for multiple vehicles with asymmetric, heterogeneous costs. The proposed algorithm in this short note aims to fill this gap.

The main focus of this article is on developing approximate solutions with a priori guarantees for solving the MHUVP. A priori guarantees are provided formally through the development of an approximation algorithm. An α - approximation algorithm is an algorithm which runs in polynomial time and finds a feasible solution whose cost is at most a specified factor (α) away from the optimal cost for every instance of the given problem. This factor is referred to as the approximation factor of the algorithm. This approximation factor provides a theoretical upper bound on the quality of the solution produced by an algorithm for any instance of the problem. These upper bounds are known a priori, i.e., they are known even before one implements the approximation algorithm for some specific instances of the problem. For these reasons, the bound provided by the approximation factor is conservative. The solutions produced by the approximation algorithms can be either used as initial solutions for other improvement heuristics or modified to improve the quality of the resulting solutions.

There is currently no approximation algorithm for the MHUVP. For a single UV case with asymmetric travel costs, Frieze et al. [18] presented a well known [log₂(p)]-approximation algorithm where p is the number of vertices. Even though there have been improvements on the [log₂(p)] -approximation factor [19], there is currently no constant factor approximation algorithm for the asymmetric, single UV case. For the multiple UV case when the travel costs are heterogeneous but symmetric, authors in [13] developed a 1.5 m approximation algorithm where m denotes the number of UVs. This result uses the 3/2 -approximation algorithm for the symmetric, single TSP [20].

The main contribution of this article is in developing a m[log₂(n + 1)] -approximation algorithm for the MHUVP where n denotes the number of targets. This is the first approximation result for a path planning problem involving heterogeneous UVs with asymmetric travel costs. This approximation algorithm first partitions the targets among UVs by solving a Linear Programming relaxation and uses the Frieze et al. algorithm [18] to find a path for each of the UVs. We then perform a simple modification of this approximation algorithm to develop a heuristic. This modification replaces the Frieze et al. algorithm with the Lin-Kernighan heuristic [21] to find solutions of superior quality. The proposed algorithms are implemented for hundreds of typical instances [22] uniformly generated in a square area with up to four UVs and 40 targets to corroborate the performance of the proposed algorithms.

This article is organized as follows: we state and mathematically formulate the problem in the next two sections. Section 4 presents the main algorithm. The proof of its approximation ratio is shown in section 5. The modified heuristic is presented in section 6. The simulation results along with conclusions are presented in sections 7 and 8 respectively.

2. Problem Statement

Let T ={1,…,n} be the set of vertices that denotes all the target locations (or targets) and D = {d₁,d₂,…,d_m} be the set of vertices that correspond to the initial depots of the m UVs. We refer to k^th UV as the vehicle which starts from depot d_k. Let V_k = {d_k}∪T denote the set of all the vertices corresponding to the k^th UV. Let (u,v) denote the directed edge from vertex u to vertex v. Let E_k be the set of all the directed edges joining any two vertices in V_k. For any S ⊂ V_k, let δ k + ( S ) represent the set of all the edges (u,v) ∈ E_k such that u ∈ S and v ∈ V_k \ S. Similarly, let δ k − ( S ) denote the set of all the edges (u,v) ∈ E_k such that u ∈ V_k \ S and v ∈ S. The cost of traveling an edge e = (u,v) ∈ E_k is denoted by c o s t e k . All the travel costs are assumed to satisfy the triangle inequality, i.e., c o s t ( u , v ) k ≤ c o s t ( u , w ) k + c o s t ( w , v ) k for all u,v,w ∈ V_k,k = 1,…,m. A path (Path_k) for the k^th UV is defined by the sequence of vertices visited by the UV, i.e., Path_k:=(d_k,u₁,u₂,…,u_ik,d_k) where u₁,…,u_ik ∈T. The cost of traveling Path_k is defined as c o s t ( P a t h k ) : = c o s t ( d k , u 1 ) k + ∑ j = 1 i k − 1 c o s t ( u j , u j + 1 ) k + c o s t ( u i k , d k ) k . The objective of the MHUVP is to find a path for each UV such that each target is visited by some UV and the sum of the travel costs of all the UVs is minimized.

3. Problem Formulation

In this section, we formulate the MHUVP as an integer linear program. This formulation will be used to develop the approximation algorithm. Let x_e^k denote the binary variable that decides whether edge e is present in the path of the k^th UV. x_e^k = 1 if and only if edge e is traveled by the k^th UV and x_e^k =0 otherwise. Let ψ_i^k be the binary variable used to assign target i to the k^th UV. ψ_i^k = 1 if target i is assigned to the k^th UV and ψ_i^k = 0 otherwise. Let x : = { x e k , e ∈ E k , k = 1 , ⋯ , m } and ψ = { ψ i k : i ∈ V k , k = 1 , ⋯ , m } denote all the decision variables. The MHUVP can be formulated as follows:

C o p t = min ∑ k = 1 m ∑ e ∈ E k c o s t e k ​ x e k

subject to

∑ k = 1 m ψ i k = 1 ∀ i ∈ T ,

(1)

∑ e ∈ δ k + ( { i } ) x e k = ψ i k ∀ i ∈ T , k = 1 , ⋯ , m ,

(2)

∑ e ∈ δ k + ( { d k } ) x e k ≥ ψ i k ∀ i ∈ T , k = 1 , ⋯ , m ,

(3)

∑ e ∈ δ k − ( { i } ) x e k = ∑ e ∈ δ k + ( { i } ) x e k ∀ i ∈ V k , k = 1 , ⋯ , m ,

(4)

∑ e ∈ δ k + ( S ) x e ≥ ψ i k ∀ S ⊆ T , i ∈ S , k = 1 , ⋯ , m ,

(5)

∑ e ∈ δ k − ( S ) x e ≥ ψ i k ∀ S ⊆ T , i ∈ S , k = 1 , ⋯ , m ,

(6)

ψ i k ∈ { 0,1 } ∀ i ∈ T , k = 1 , ⋯ , m ,

(7)

x e k ∈ { 0,1 } ∀ e ∈ E k , k = 1 , ⋯ , m .

(8)

The constraints in (1) state that each target must be assigned to exactly one UV. The out-degree constraints of each target are stated in (2). The constraints in (3) state that the out-degree of a depot must be at least equal to one if any target is assigned to the UV corresponding to the depot. The constraints in (4) state that the in-degree and out-degree of any vertex must be equal. If target i is assigned to the k^th UV, then the number of edges in E_k leaving any subset of targets containing i must be at least equal to 1. This connectivity constraint is stated in (5). Similarly, if target i is assigned to the k^th UV, constraints in (6) state that the number of edges in E_k entering any subset of targets containing i must be at least equal to 1. The constraints in (5),(6) are generally referred to as the cut constraints in the literature.

4. Approximation Algorithm for Solving the MHUVP

The approximation algorithm (Approx) first partitions the targets by solving a linear programming (LP) relaxation of the integer program and then uses the Frieze et al. algorithm [18] to find a path for each UV. The steps involved in this algorithm are as follows:

Solve the following LP relaxation of the integer program.

C L P * = min ∑ k = 1 m ∑ e ∈ E k c o s t e k ​ x e k

(9)

subject to

∑ k = 1 m ψ i k = 1 ∀ i ∈ T

(10)

∑ e ∈ δ k − ( { i } ) x e k = ∑ e ∈ δ k + ( { i } ) x e k ∀ i ∈ V k , k = 1 , ⋯ , m

(11)

∑ e ∈ δ k ( S ) x e ≥ ψ i k ∀ S ⊆ T , i ∈ S , k = 1 , ⋯ , m

(12)

∑ e ∈ δ k − ( S ) x e ≥ ψ i k ∀ S ⊆ T , i ∈ S , k = 1 , ⋯ , m

(13)

ψ i k ≥ 0 ∀ i ∈ T , k = 1 , ⋯ , m

(14)

x e k ≥ 0 ∀ e ∈ E k , k = 1 , ⋯ , m

(15)

Let the linear program above be denoted by LP* and its optimal solution be (x*, ψ*).

Partitioning the targets: ψ i k * denotes the optimal fraction of the target assigned to the k^th UV in LP*. Assign target i to the k ^ t h UV that corresponds to the largest fraction, i.e., ψ i k ^ * = max l = 1 m ψ i l * . Break ties arbitrarily. At the end of this step, each target will be assigned to exactly one UV. For k = 1, …, m, let T _k be the set of targets assigned to the k^th UV. Note that for any target i ∈ T k , ψ i k * ≥ 1 m (This fact will be later used in the proof of the approximation ratio).

Path synthesis: For k = 1,…,m, apply the Frieze et al. algorithm to the vertices in T_k ∪ {d_k} to obtain a path for the k^th UV.

5. Proof of the Approximation Ratio

First, we show that Approx runs in polynomial time. Second, we show the bound on the deviation of the solution produced by Approx from the optimum.

Lemma 1Approx runs in polynomial time.

Proof. The main steps in Approx include solving LP* and implementing the Frieze et al. algorithm for each of the UVs. It is already known that Frieze et al. algorithm runs in polynomial time [18]. Therefore, the proof of the lemma would be complete if LP* is solvable in polynomial time. The difficulty involved in solving LP* mainly rests on addressing the cut constraints (12),(13) as all the other constraints are polynomial in size of the problem. For all i ∈ T and k = 1,…,m, using the max-flow, min-cut theorem, the cut constraints in (12) can be equivalently replaced with flow constraints which require at least a shipment of ψ i k units of commodity from target i to depot d_k. Therefore, for target i and depot d_k, the cut constraints in (12) can be replaced with the following set of flow constraints that are polynomial in size:

f u v i k ≤ x ( u , v ) k ​ ∀ e = ( u , v ) ∈ V k ,

(16)

∑ v ∈ T ( f u v i k − f v u i k ) = { ψ i k u = i , 0 u ∉ { i , d k } , − ψ i k u = d k ,

(17)

f u v i k ≥ 0 ​ ​ ​ ∀ u , v ∈ V k .

(18)

In the above constraints, f u v i k denotes the amount of commodity shipped from target i to depot d_k flowing through edge (u,v). These flow constraints can be expressed for each target and depot, and therefore, all the cut constraints in (12) can be expressed equivalently using a polynomial number of flow constraints. The same idea can also be applied to the cut constraints in (13). Therefore, LP* can be solved in polynomial time [23]. Hence proved.

The following is the main result of this article:

Theorem 5.1 The approximation ratio of Approx is m[log₂(n + 1)].

Proof. For a single UV, asymmetric TSP (ATSP) with n targets and one depot, Williamson [24] showed that the cost of the solution produced by the Frieze et al. algorithm is at most log₂(n + 1) times the optimal cost of the Held-Karp relaxation of the ATSP. In the context of our problem, this result implies that c o s t ( P a t h k ) ≤ log 2 ( n + 1 ) C r e l a x k where cost(Path_k) denotes the cost of the path found for the k^th UV and C r e l a x k represents the optimal Held-Karp relaxation cost. If at least one target is assigned to the k^th UV (|T _k | ≥1), this Held-Karp relaxation cost is defined as follows:

C r e l a x k = min ∑ e ∈ E k c o s t e k ​ x e k ∑ e ∈ δ k + ( { i } ) x e k = ∑ e ∈ δ k − ( { i } ) x e k = 0 ∀ i ∈ T \ T k ,

(19)

∑ e ∈ δ k + ( { i } ) x e k = ∑ e ∈ δ k − ( { i } ) x e k = 1 ∀ i ∈ T k ∪ { d k } ,

(20)

∑ e ∈ δ k + ( S ) x e ≥ 1 ∀ S ⊂ T k ∪ { d k } ,

(21)

x e k ≥ 0 ∀ e ∈ E k .

(22)

As the travel costs satisfy the triangle inequality, Nguyen [25] has shown that the degree constraints in the above formulation can be relaxed without changing the optimal cost of the Held-Karp relaxation as follows:

∑ e ∈ δ k − ( { i } ) x e k = ∑ e ∈ δ k + ( { i } ) x e k ∀ i ∈ V k , ∑ e ∈ δ k − ( { i } ) x e k = ∑ e ∈ δ k + ( { i } ) x e k ∀ i ∈ V k ,

(23)

∑ e ∈ δ k + ( S ) x e ≥ 1 ∀ S ⊂ T k ∪ { d k } ,

(24)

x e k ≥ 0 ∀ e ∈ E k .

(25)

Now, we prove that the sum of the optimal cost of the Held-Karp relaxations, ∑ k = 1 m C r e l a x k , can be bounded by at most m times the optimal cost of LP* defined in (9)–(15), i.e., ∑ k = 1 m C r e l a x k ≤ m C L P * . Let (x*,ψ*) be an optimal solution of LP*. Then, construct a solution x ˜ e k : = m x e k * ​ ∀ e ∈ V k , k = 1 , ⋯ , m . x ˜ e k is a feasible solution to the linear program defined in (23)-(25) due to the following reasons:

For ​ all ​ ​ i ∈ V k ​ ​ ​ ​ ​ ​ ∑ e ∈ δ k − ( { i } ) x ˜ e k = m ∑ e ∈ δ k + ( { i } ) x e k * ​ ​ ​ = m ∑ e ∈ δ k − ( { i } ) x e k * ( from ​ ​ equation 11 ) ​ ​ ​ = ∑ e ∈ δ k − ( { i } ) x ˜ e k .

For ​ all ​ ​ S ⊆ T k , | S | ≥ 1 ​ , ​ ​ ​ ​ ​ ​ ​ ​ ​ ∑ e ∈ δ k + ( S ) x ˜ e k = m ∑ e ∈ δ k + ( S ) x e k * ≥ m ψ i k * ​ ​ ∀ i ∈ S ( from ​ ​ equation ​ ​ 12 ) ≥ 1. ( A s ​ ​ ψ i k * ≥ 1 m for ​ ​ any ​ ​ i ∈ T k )

For ​ all ​ ​ S ⊂ T k ∪ { d k } , | S ∩ { d k } | ≥ 1 , ​ ​ ​ ​ ​ ​ ∑ e ∈ δ k + ( S ) x ˜ e k = ∑ e ∈ δ k − ( V k \ S ) x ˜ e k = m ∑ e ∈ δ k − ( V k \ S ) x e k * ≥ m ψ i k * ​ ​ for ​ ​ all i ∈ V k \ S ( from ​ ​ equation ​ ​ 13 ) ≥ 1. ( As ​ ​ | T k \ S | ≥ 1 , ​ ​ there ​ ​ is ​ ​ at ​ ​ least ​ ​ one target ​ ​ i ∈ V k \ S ​ ​ such ​ ​ that ​ ​ ψ i k * ≥ 1 m ) ,

Therefore, for k = 1 , ⋯ , m , x ˜ e k is a feasible solution for the linear program in (23)–(25). Consequently, ∑ k = 1 m C r e l a x k ≤ ∑ e ∈ E k c o s t e k ​ x ˜ e k = m ∑ e ∈ E k c o s t e k ​ x e k * = m C L P * . By putting together all the above results, we have:

∑ k = 1 m c o s t ( P a t h k ) ≤ ⌈ log 2 ( n + 1 ) ⌉ ∑ k = 1 m C r e l a x k ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ≤ m ⌈ log 2 ( n + 1 ) ⌉ C L P * ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ≤ m ⌈ log 2 ( n + 1 ) ⌉ C o p t .

6. Heuristic

A modification of the approximation algorithm (referred to as Approx_lkh) can also be easily obtained for improving the quality of the solutions produced by Approx. In this modification, LP* is solved and the partitioning of the targets is performed as in Approx. After partitioning, the Lin-Kernigan Heuristic (LKH) [21] is used instead of the Frieze et al. algorithm to find a path for each UV. The LKH is considered to be one of the best algorithms for solving a single vehicle TSP. Computational results in the next section show that combining the critical partitioning step of Approx and LKH leads to finding high quality solutions for MHUVP.

7. Simulation results

Approx and Approx_lkh was implemented for hundreds of problem instances involving unmanned aerial vehicles (UAVs). The LKH program by Helsgaun [26] was used to solve the ATSP. The LKH program was run without changing any of its default settings.

For the simulations, the number of UAVs was varied from two to four and the number of targets were varied from 20 to 40. All the targets were randomly generated in a square of area five × 5 km ² using a uniform distribution. For the purposes of simulating heterogeneity, the minimum turning radius (r) of each UAV was chosen uniformly from the interval [100,150] metres. The approach or the heading angle at each target was fixed and randomly chosen from a uniform distribution on the interval [0,2π]. Dubins' result [27] was used to calculate the minimum distance required to travel between any two targets. These travel distances are asymmetric and satisfy the triangle inequality.

For a given number of UVs (m) and targets (n), 50 instances were randomly generated. For a given problem instance I, the bound on the a posteriori guarantee provided by an algorithm is defined as C s o l f I / C L P * I where C s o l f I is the cost of the feasible solution found by the algorithm and C L P * I is the optimal cost of the LP* defined in (9). All the simulations were run on a Dell Precision T5500 workstation (Intel Xeon E5630 processor @ 2.53GHz, 12GB RAM). The average CPU time required by the algorithms for an instance was in the order of a minute.

The average a posteriori guarantee of the solutions found by Approx and Approx_lkh are shown in tables 1, 2, 3 along with the theoretical approximation guarantees. The average computation times of the algorithms are shown in table 4. There are two conclusions from these results. First, quality solutions that are, on an average, within 1% of the optimum can be found by using Approx_lkh. Second, even though the theoretical a priori guarantee of the approximation algorithm is quite large, the simulation results show that the a posteriori guarantee of the solutions found by Approx was approximately 1.50. These results imply that Approx finds solutions with bounds that are significantly better than the guarantees indicated by the approximation factor. Sample solutions using Approx and Approx_lkh are shown in figure 1. We commonly observed the following over most of the simulations: The paths using Approx criss-crossed more corroborating its weak performance with respect to the quality of solutions produced by Approx. As there were no limits on the number of targets visited by any vehicle, the vehicle with a least turning radius usually visits most of the targets unless the depot location of a vehicle with a larger turning radius is closer to some of the targets. Future work can address this issue by including capacity constraints for each vehicle. As expected, the Dubins distance between any two adjacent targets in a tour (for any algorithm) is closer to its corresponding Euclidean distance if the targets are reasonably far apart. The complications typically appear when targets are closely located.

OPEN IN VIEWER

Table 1.

Comparison of the theoretical and simulation results for m=2

No. of Targets	Theoretical upper bound (m[log2(n + 1)])	A Posteriori Bound using Approx	A Posteriori Bound using Approxlkh
20	8.7846	1.4942	1.0058
25	9.4009	1.4516	1.0075
30	9.908	1.4896	1.0061
35	10.3399	1.5251	1.0066
40	10.7151	1.4830	1.0089

OPEN IN VIEWER

Table 2.

Comparison of the theoretical and simulation results for m=3

No. of Targets	Theoretical Upper Bound (m[log2(n + 1)])	A Posteriori Bound using Approx	A Posteriori Bound using Approxlkh
20	13.1770	1.471	1.0045
25	14.1013	1.4825	1.0068
30	14.8626	1.4577	1.009
35	15.5098	1.4897	1.008
40	16.0727	1.4843	1.0102

OPEN IN VIEWER

Table 3.

Comparison of the theoretical and simulation results for m=4

No. of Targets	Theoretical upper bound (m[log2(n + 1)])	A Posteriori Bound using Approx	A Posteriori Bound using Approxlkh
20	17.5693	1.4469	1.0046
25	18.8018	1.4769	1.0077
30	19.8168	1.4477	1.0072
35	20.6797	1.5049	1.0085
40	21.4302	1.4973	1.0092

OPEN IN VIEWER

Table 4.

Average computation times (in secs) for the proposed algorithms

	2 vehicles		3 vehicles		vehicles
No. of Targets	Approx	Approxlkh	Approx	Approxlkh	Approx	Approxlkh
20	1.648	1.571	2.998	3.121	3.930	3.871
25	5.189	5.157	8.373	8.465	11.013	11.12
30	12.468	12.320	19.690	19.875	23.968	24.351
35	27.478	27.469	40.025	40.567	45.508	45.795
40	52.339	53.776	71.348	72.240	86.713	86.908

OPEN IN VIEWER

Figure 1.

An example of the UAV paths obtained by Approx and Approx_lkh respectively. Two heterogeneous vehicles were involved in these simulations.

8. Conclusions

An approximation algorithm and a heuristic was presented for a multiple, heterogeneous, unmanned vehicle path planning problem. Simulation results were also presented to show that the proposed approach can be used to find high quality solutions to the path planning problem. The approximation factor of the proposed algorithm is m[log₂(n + 1)]. Future work can aim to improve this approximation factor by exploiting the monotonicity property present in the traveling costs of the UAVs, i.e., without loss of generality, one can always order the vehicles such that the travel costs between any two target locations increases monotonically with the index of the vehicles: c o s t ( i , j ) 1 ≤ c o s t ( i , j ) 2 ⋯ ≤ c o s t ( i , j ) m . Recently, this property has been used to improve the approximation factor for the symmetric case in [17].