Abstract
We consider a path planning problem where a team of Unmanned Vehicles (UVs) is required to visit a given set of targets. The UVs are assumed to carry different sensors, and as a result, there are vehicle-target constraints that require each UV to visit a distinct subset of targets. The objective of the path planning problem is to find a path for each UV such that each target is visited at least once by some vehicle, the vehicle-target constraints are satisfied and the total distance travelled by the vehicles is a minimum. This path planning problem is a generalization of the Hamiltonian path problem and is NP-Hard. We develop a primal-dual heuristic and incorporate the heuristic in a Lagrangian relaxation procedure to find good, feasible solutions and lower bounds for the path planning problem. Computational results show that solutions whose costs are on an average within 14% of the optimum can be obtained relatively quickly for the path planning problem involving five UVs and 40 targets.
Keywords
1. Introduction
Surveillance applications involving Unmanned Vehicles (UVs) require UVs with different capabilities to gather information about a set of targets [3]. Information is often gathered using the sensors on-board the UVs while visiting the targets. In these missions, it is possible that UVs carry different sensors due to resource constraints. As a result, there may be vehicle-target constraints that require a specific vehicle to visit a subset of targets. We consider a basic path planning problem that commonly arises in these surveillance applications: Given a set of vehicles and targets, the objective of this problem is to find a path for each vehicle such that each target is visited at least once by some vehicle, the vehicle-target constraints are satisfied and the sum of the distances travelled by all the vehicles is a minimum.
The difficulty with solving this path planning problem for the UVs is mainly combinatorial. While the motion constraints of the UV complicate the problem, they are not a source of difficulty as long as the trajectory of least distance to go from an origin to a destination can be efficiently computed; it does not matter whether the vehicle is modelled as a double integrator or a Dubins [9] or a Reeds-Shepp [17] vehicle. To illustrate this point about the combinatorial nature of the problem, consider the path planning problem for just one UV modelled as a Dubins vehicle [9]. If the “optimal” heading angles were to be specified at each and every target, the problem of finding the optimal sequence of targets to be visited reduces to the notoriously hard Hamiltonian Path Problem (HPP), which is known to be NP-Hard in the literature [1]. However, if the “optimal” sequence were to be given, the “optimal” heading angles can be found using dynamic programming and can be determined arbitrarily accurately as a shortest path problem in a network [21], a problem for which efficient algorithms exist [7].
For the above reasons, we assume that the heading angle at which a vehicle must visit each target is known. We also assume that the cost of travelling for a UV from target A at an initial heading angle (ψA) to target B at a final heading angle (ψB) can be computed and is known. Therefore, finding a sequence of targets to visit for a vehicle is equivalent to specifying the path for the vehicle. Using these assumptions, our problem can be viewed as a multiple depot, multiple vehicle, asymmetric 1 HPP with additional vehicle-target constraints. In the context of UVs, we are not aware of algorithms that directly address this problem in the literature. However, there are approximation algorithms and heuristics for similar problems. For example, Doshi et al. [8] present heuristics and an approximation algorithm for a 2-depot, symmetric HPP. Approximation algorithms are presented for a multiple depot, symmetric HPP in [20]. A heuristic based on a transformation method is also presented for a multiple depot, asymmetric, heterogeneous TSP by Oberlin et al. in [15]. The combinatorial problem considered in this article is a special case of more difficult vehicle routing problems for which heuristics are presented in [4],[6],[5],[16]. In this article, we develop a primal-dual heuristic and incorporate the heuristic in a Lagrangian relaxation procedure to find good, feasible solutions and lower bounds for the path planning problem. Computational results show that solutions whose costs are within 14% of the optimum can be obtained relatively quickly for the path planning problem involving five UVs and 40 targets.
2. Problem Statement
Let there be n vehicles initially located at distinct depots denoted by d1, d2, …, dn and let D:={d1, …, dn (Refer to Figure 1). Let T denote the set of targets that needs to be visited by the vehicles. Let Fk ⊆ T represent the subset of targets that must be visited by the vehicle initially located at depot dk. C:=T\Uni=1 Fk denotes the subset of all the targets that can be visited by any vehicle. We assume
the paths for the vehicles start at their respective depots,
each target in Fk is visited at least once by the vehicle located at depot dk,
each common target is visited at least once by some vehicle, and,
the total cost of travelling the edges in all the paths is a minimum.
An illustration of a feasible solution to this MVP is shown in Figure 1.

A solution to the path planning problem. F1, F2 and F3 denote the subsets of targets that need to be visited by the vehicles at depots d1, d2 and d3 respectively.
3. Problem Formulation
In this section, we formulate the problem as an integer program. Let xij be the binary decision variable that represents the presence of edge (i, j) in the solution; that is, xij is equal to 1 if a vehicle is travelling from vertexi to vertex j and is equal to 0 otherwise. For any S ⊂ V, let ™−(S) ⊆ E denote the set of all the edges directed into S, i.e, for any (i, j) ∈ ™−(S),i ∈ V\S and j ∈ S. Similarly, let ™+(S) represent the set of all the edges directed out of S, that is, for any (i, j) ∈ ™+(S),i ∈ S and j ∈ V\S. The MVP is formulated as an integer program as follows:
subject to
Equations (1) and (2) require that the in-degree and out-degree of each vertex be at most equal to one. Constraints (3), (4) ensure that each common target is connected to at least one depot and each functionally heterogeneous target is connected to its respective depot. Specifically, equation (3) requires that there must be at least one incoming edge to every subset of common targets. Similarly, equation (4), requires that there must be at least one suitable edge directed into every subset consisting of common targets and at least one target from Fk for any k. For example, if the subset S consists of a target from F1 and some common targets, then a suitable incoming edge into S can be directed either from the depot d1, or from a common target, or from a target in F1 but not in S.
4. Lagrangian Relaxation of the MVP
A Lagrangian relaxation for the MVP can be obtained by removing the constraints that complicate the formulation in (1)-–(4), and penalizing them in the objective if they are violated. This idea was first successfully applied to a single vehicle routing problem by Held and Karp in [13]. Suppose the degree constraints in (1) and (2) are relaxed, and let λi and µi be the penalty variables corresponding to the degree constraint of vertex i in (1) and (2) respectively. Let λ denote the vector (λ1,…,λ|V|), and μ represent (μ1,…,μ|V|). Then, for an λ, µ≥ 0, we obtain a Lagrangian relaxation as follows:
subject to
The objective of the above relaxation can be re-written as follows:
where the modified cost ĉ ij = c ij + λj + μ i , and K(λ,μ) = Σj∈Vλj + Σi∈Vμi.
It is well known in the literature [11] that the optimal cost of a Lagrangian relaxation always provides a lower bound for the corresponding integer programming problem. In relation to the MVP, this result implies that φ*(λ,µ) ≤f*for any λ,µ ≥0 where f* denotes the optimal cost of the MVP. One can then obtain the best or the tightest lower bound for the MVP by solving the following problem: maxλ,µ≥0 φ* (λ,µ). This problem of maximizing the lower bound will be referred to as the Lagrangian dual problem in this article.
Suppose F denotes the set of all the solutions (a solution consists of a forest of edges; refer to Figure 2 for an illustration) that satisfy the constraints in (7)–(9). Then, the Lagrangian dual problem can be stated as
Even though the Lagrangian dual problem is a non-smooth optimization problem, note that the objective φ* (λ,µ) is a concave function of the penalty variables in λ and µ. Therefore, one can use a subgradient algorithm as in [13],[18] to solve the Lagrangian dual problem.
A crucial part of solving the Lagrangian dual problem involves solving the minimization problem in (6)–(9). This minimization problem aims to find a forest of edges such that
there is a directed path from a depot to each of the common targets,
there is a directed path from the depot dk to each of the targets in Fk for k =1,…,n, and,
the sum of the cost of the edges in this forest is a minimum.
If there is only one vehicle and there are no vehicle-target constraints, this minimization problem reduces to solving a directed spanning tree which can be computed in polynomial time. If there are multiple vehicles with additional vehicle-target constraints, it is not yet known if this minimization problem can be solved in polynomial time. However, we circumvent this difficulty by developing a primal-dual heuristic which provides a good lower bound and a feasible solution to the minimization problem. We use this feasible solution (a forest) to partition the common targets such that each common target is assigned to some vehicle. Once a subset of targets is assigned to a vehicle, we use the Lin-Kernighan Heuristic (LKH) to find a path for the vehicle. Using the above procedure, given the penalty variables in λ and µ, one can find a lower bound and a feasible solution for the MVP. In the next section, we discuss the primal-dual heuristic. Later, we present the subgradient algorithm that uses this primal-dual heuristic to solve the Lagrangian dual problem.

Feasible solutions for the Lagrangian relaxation. Notice that in (a) some common targets are connected to depots d1 and d2, whereas each common target is connected to exactly one depot in (b).
5. Primal Dual Heuristic for solving the Lagrangian Relaxation
Given the penalty variables in λ and µ, consider the basic minimization problem in the Lagrangian relaxation:
Henceforth, the above forest problem will be referred to as the primal problem, and will be used to find a lower bound and partition the targets. A dual corresponding to a linear programming relaxation (only equation (12) is relaxed to xij ≥0 ∀ (i, j) ∈ E) of the primal problem can be stated as follows:
subject to
The primal-dual heuristic follows the greedy procedure outlined by Goemans and Williamson in [12]. The basic structure of this heuristic involves maintaining a forest of edges and a feasible solution to the dual problem. The edges in this forest are candidates for the set of edges that finally appear in the output of the algorithm. Initially, the forest is empty and all the dual variables are set to zero. As the forest is empty, there may be several constraints that are violated in the primal problem initially. The primal-dual heuristic is an iterative algorithm where in each iteration, the subset of vertices S corresponding to one such violated constraint in the primal problem is chosen and the dual variable yS is increased until one of the constraints corresponding to say edge (i, j) in (13) becomes tight. A constraint corresponding to edge (i, j) becomes tight if
The algorithm then adds the edge (i, j) to the forest, and proceeds to the next iteration until the forest contains a feasible solution to the primal problem.
Once the forest becomes feasible to the primal problem, it is possible that there is more than one path from a depot to some of the targets. Therefore, in the last step of the heuristic, we eliminate some unnecessary edges in the forest while ensuring that the pruned solution is still feasible. We do this by sorting the edges in the reverse of the order in which they were added to the forest and removing them if the pruned solution is feasible to the primal problem. This pruning procedure is referred to as the reverse delete step in the literature [12], [19] and is usually used for developing good approximation algorithms for NP-Hard problems. The pseudo code of the primal-dual heuristic is presented in Algorithm 1. The key step of this primal-dual heuristic requires one to find a violated constraint (if any) in the primal problem given a forest of edges. This step is explained in the following subsection.
5.1 Identifying Violated Constraints
A forest is infeasible to the primal problem if there is no directed path from any of the depots to a common target or if there is no directed path from the depot dk to any of the targets in Fk for some k ∈ {1,…,n}. If S ⊆ C and there is no incoming edge into S, then the constraint associated with S in (10) will be violated. If S ⊆C ∪ Fk, |S ∩ Fk |≥ 1 for some k∈ {1,…,n} and there is no incoming edge into S from a vertex in C ∪ Fk ∪ {dk}, then the constraint associated with S in (11) will be violated. Henceforth, given a forest, we will refer to a set as a violated set if the forest does not satisfy the constraint associated with the set in the primal problem.
4 : count ← 0
5 : Given a forest which is infeasible, one can use an algorithm similar to the Edmonds' algorithm [10] for the directed minimum spanning tree problem to find a violated constraint (or a violated set). Our algorithm works as follows: first, find a strongly connected component S such that S ⊆C∪Fk for some k. Contract all the arcs in this strongly connected component so that S is replaced with a single pseudo node as shown in figures 3,4. For example, in Figure 3-b, {1,7,6} is a strongly connected component and is replaced with a pseudo node in Figure 3-c. Repeat this procedure until all such strongly connected components are replaced with their corresponding pseudo nodes. Let this contracted network be denoted by Gc = (Vc, Ec) where Vc and Ec represent the set of all the vertices and edges in the contracted network respectively. Mark each vertex v in the contracted network with a label
Next, we repeat the following procedure for each marked vertex v in this contracted network until we find a violated set. Let G'(v) represent the largest possible subgraph of Gc such that all the vertices in G'(v) are marked, for each vertex u(≠v) in G'(v), label(u)∈ {label(v),
Expand each of the pseudo nodes in G'(v) and let S denote all the vertices in the expanded graph obtained from G'(v). As we ensure that the label of each of the vertices in G'(v) is either

An example showing the procedure for finding the violated sets.
We output S as a violated set for the forest if label (v) =
6. Subgradient algorithm for addressing the Lagrangian dual problem
Given the penalty variables in λ and µ, the primal-dual heuristic presented in the previous section can be used to find a lower bound and a good, feasible solution for the Lagrangian relaxation problem in (6)–(9). In this section, we provide a subgradient algorithm for solving the Lagrangian dual problem.
The subgradient algorithm is an iterative procedure where the penalty variables are updated in each iteration based on an improving direction of the objective of the Lagrangian dual problem. Let [λ]k and [µ]k indicate the values of λ and µ during the kth iteration respectively. During the kth iteration, we compute a new set of penalty vectors, [λ]k+1 and [µ]k+1 using an update scheme that aims to change [λ]k and [µ]k along an improving direction defined by the subgradient as follows: let gλi and gµi denote the subgradient corresponding to the penalty variables λi and µi for i =1,…, |V| where

Another example showing the procedure for finding the violated sets.
Let g = (gλ1,…,gλ|v|,gμ1,…,gμ|v|) represent the vector of all the subgradients, and let [g]k represent this vector during the kth iteration. The new update [λ,µ]k+1 is computed as follows:
where αk is the step size defined by
In the above equation, γk is a parameter in the interval (0,2], f* denotes the optimal cost of the MVP and ‖[g]k‖ denotes the l2 –norm of [g]k. As f* may not be known, during the kth iteration, we use the cost of the best solution found for the MVP during the first k iterations as a proxy for f*. The expression (19) is referred to as the Polyak's Rule II. A common practice is to start with a value of γ0 ∈(0,2] and reduce its value by a constant factor if the dual objective does not increase in a specified number of iterations [18]. The overall procedure for solving the Lagrangian dual problem is given in Algorithm 2.
1: Choose the initial penalty variables in [λ0, µ0. Let k = 0.
2: Solve the Lagrangian relaxation for [λ,µ]k. This step would provide a lower bound LBk to φ* (λk, µk) and a forest Tk for the primal problem.
3: Use the forest Tk to assign each target to one of the vehicles. Break the ties arbitrarily if a common target is connected to more than one depot. Once the targets have been partitioned among the vehicles, use the LKH [14] to find a corresponding path for each of the vehicles. Let the sum of the cost of travelling all the paths be denoted by fk.
4: Terminate the algorithm if the duality gap (minki=1 fi - maxki=1 LBi) is less than a constant ε, or if the number of iterations, k, has reached a maximum value; otherwise, continue to the next step.
5: Compute [λ,µ]k+1 using equation (18), i.e., [λ,µ]k+1 = [λ,µ]k + αk [g]k. αk is the step size given by Polyak's Rule II as follows:
In the above rule, γ0 is chosen to be in the interval (0,2].
For
6: Let k ← k + 1 and go to step 2.
7. Computational Results
The objective of this section is to compute the quality of the solutions produced by the proposed algorithms for the MVP. In that pursuit, we conducted two sets of simulations for the MVP - in the first set, we fixed the number of vehicles and varied the number of targets; in the second set, we fixed the number of targets and varied the number of vehicles. All the simulations were run on a Dell Precision T5500 workstation (Intel Xeon E5630 processor @ 2.53 GHz, 12GB RAM).
For a given number of vehicles and targets, 50 test instances were generated. The coordinates of the targets and depots were generated randomly in a square of size 5000 m using a uniform distribution. For the simulations, we model the UV as a Dubins car that travels at a constant velocity and has a lower bound on turning radius. The minimum turning radius of all the vehicles was chosen to be equal to 250 m. For each generated target, a heading angle was selected randomly in the interval [0,2π]. The Dubins' result [9] was used to calculate the minimum distance required for a vehicle to travel between any two targets. Each vehicle was also assigned a subset containing three targets which the vehicle must visit. An upper bound on the deviation of the cost of the suboptimal solutions found by the subgradient algorithm from the optimum is defined as
where fI is the cost of the best solution and LBI is the best lower bound obtained for an instance I by the subgradient algorithm. The LKH program by Helsgaun [14] available at http://www.akira.ruc.dk/keld/research/LKH/was used to solve the HPP in the subgradient algorithm. The LKH program was run without changing any of its default settings.
While applying the subgradient algorithm to solve the Lagrangian dual problem, all the penalty variables were initially set to zero. γ0 was initially assigned a value of 0.001 in Equation (19), and its value was reduced by a factor of 2 whenever φ([λ]k, [µ]k) failed to increase in two iterations. As the iterations progresses, we can expect the best lower bound produced by the subgradient algorithm to continue to increase until it converges to a fixed value. Similarly, the cost (upper bound) of the best solution obtained by the subgradient algorithm may continue to decrease as the iterations progress. For each instance, the subgradient algorithm was allowed to run for at most 50 iterations. Figure 5 illustrates the way in which the best upper bound and the lower bound change with the number of iterations for a 30 target instance.

Variation of the best upper bound and the best lower bound obtained by the subgradient algorithm as a function of the number of iterations for an instance.
In the first set of simulations, the number of vehicles was fixed to be equal to five, and the number of targets was allowed to vary from 15 to 40. The first plot in Figure 6 shows the average deviation of the cost of the suboptimal solutions produced by the algorithm from the optimum. This figure also shows the average computation time required to implement the subgradient algorithm. For the test instances considered in the first set of simulations, the proposed algorithms were able to find solutions whose cost, on an average, is at most 12% away from the optimum.

Computational results for the MVP: number of vehicles fixed.
In the second set of simulations, the number of targets was fixed at 35 and the number of vehicles was varied from 2 to 7. The first plot in Figure 7 shows the average deviation of the cost of the suboptimal solutions produced by the algorithm from the optimum. Specifically, for the test instances considered in the second set of simulations, the proposed algorithms were able to find solutions whose cost, on an average, is at most 14% away from the optimum. In addition, the average computation time required for implementing the proposed algorithms was at most equal to a minute. Considering that the MVP with additional vehicle-target constraints is a difficult problem to solve, these results indicate that the approach proposed in this article is promising. Figures 8 and 9 show the Dubins' paths found by the proposed algorithms for a couple of test instances.

Computational results for the MVP: number of targets fixed.
8. Conclusion
A primal-dual heuristic was combined with the Lagrangian relaxation technique to find good solutions for a heterogeneous vehicle routing problem. Future work related to this problem can focus on two research areas. Firstly, the feasible solutions and the bounds found by the developed heuristics can be used in branch and bound solvers in order to obtain optimal solutions for the related heterogeneous vehicle routing problems. Secondly, one can aim to develop primal-dual approximation algorithms for heterogeneous travelling salesman problems where there are additional vehicle-target constraints. A result in this direction has already been reached with regard to a special case of a heterogeneous travelling salesman problem in [2].

An example showing the Dubins' paths found by the algorithms for an instance with two vehicles.

An example showing the Dubins' paths found by the algorithms for an instance with three vehicles.
Footnotes
10. Acknowledgements
This material is based upon work supported by the National Science Foundation under grant no. 1015066 and the Air Force Office of Scientific Research under grant no. FA9550-10-1-0392.
1
The cost of travelling for a UV from target A at an initial heading angle (ψA) to target B at a final heading angle (ψB) may be different than the cost of travelling from target B at ψB to target A at ψA.
