A multi-vehicle covering tour problem with speed optimization
Joshua T. Margolis Yongjia Song Scott J. Mason


 Interpreting Undesirable Pixels for Image Classification on Black-Box Models
    [
    Received: date / Accepted: date
=============================================================================



* Corresponding author

E­-mail addresses: jtmargo@clemson.edu (J.T. Margolis), yongjis@clemson.edu (Y. Song),  mason@clemson.edu (S.J. Mason)





The multi-vehicle covering tour problem with time windows (MCTPTW) aims to construct a set of maximal coverage routes for a fleet of vehicles that serve (observe) a secondary set of sites given a fixed time schedule, coverage requirements, and energy restrictions. The problem is formulated as a mixed-integer second-order cone programming (MISOCP) model and is an extension of both the multi-covering tour problem and the vehicle routing problem with time windows under energy constraints. 
Further, we study a special case of the proposed model and develop a labeling algorithm to solve its Lagrangian relaxation problem, which exploits the combinatorial structure exhibited by an optimal solution to the Lagrangian relaxation.









§ INTRODUCTION






Unmanned (or autonomous) vehicles are aerial, terrestrial, and aquatic vehicles that operate without human intervention; they are becoming more prevalent in both commercial and governmental applications. As the technology surrounding them continues to advance, they will continue to be utilized in new and innovative ways. In particular, creating autonomous vehicles capable of adapting to real-time situations has recently become a priority in multiple industries. Rideshare companies (e.g., Uber and Lyft), technology companies (e.g., Apple and Google), and car manufacturers (e.g., Ford and General Motors) are all introducing self-driving cars capable of operating without human assistance. 

 
Research surrounding unmanned vehicles has highlighted the increasing desire to utilize them for a variety of purposes 
as they have many advantages: 1) they can be used to assist or even replace humans performing hazardous assignments; 2) they have increased maneuverability and deployability; and 3) they have low operation and maintenance costs <cit.>. For instance, unmanned aerial vehicles (UAVs) can be deployed for agricultural purposes <cit.>, traffic monitoring and management <cit.>, forest fire detection <cit.>, product delivery <cit.>, and 
disaster relief efforts and humanitarian causes <cit.>. For an extensive survey of additional civil UAV applications, we refer the reader to 
<cit.>. In regards to military environments, numerous research papers focus on routing unmanned vehicles for reconnaissance <cit.>, target engagement <cit.>, and surveillance missions <cit.>. 



In this paper, we study unmanned vehicle routing for target surveillance with specific vehicle characteristics (e.g., battery life, payload, speed, etc.) which can be considered as a variant of the vehicle routing problem (VRP). First studied by 
<cit.>, the VRP aims to route a collection of vehicles at minimal distance through a set of locations such that each vehicle begins and ends its journey at a fixed depot. Many survey papers exist that provide excellent summaries of the research surrounding the VRP <cit.>. Variants include the capacitated VRP <cit.>, the VRP with pickup and delivery <cit.>, and the VRP with time-windows (VRPTW) <cit.>. Exact solution methodologies for the VRPTW have been explored previously in the literature, including 
1) branch­-and­-cut algorithms <cit.> and 2) branch­-cut­-and­-price algorithms <cit.>.
For a more detailed summary, we refer the reader to recent VRP review papers <cit.>.






Specifically, we focus on the multi-vehicle covering tour problem with time windows (MCTPTW), which aims to construct a set of maximal coverage routes for a fleet of vehicles that serve (observe) a secondary set of sites given fixed time window constraints, coverage requirements, and energy restrictions. Vehicles are routed through coverage areas by visiting a set of waypoints (locations); we utilize the covering tour problem (CTP) as a basis for our model. The CTP seeks to find a minimum­-length Hamiltonian cycle along a set of vertices that passes within a predetermined distance of a secondary set of nodes so that these nodes are “covered” by the cycle <cit.>. <cit.> extend the CTP to include multiple routes (MCTP) and solve the problem heuristically. <cit.> was the first to solve the MCTP exactly via a branch­-and­-price algorithm, where the resulting subproblem is a variant of the traveling salesman problem with profits and reduces to a ring star problem that is solved using a branch­-and­-cut algorithm. 
<cit.> analyze a variant of the problem with no tour length constraints and develop a branch­-and­-cut algorithm and a two-phase metaheuristic derived from evolutionary local search to solve it. <cit.> introduce stochastic demand to the MCTP and optimize with respect to multiple objectives (i.e., cost and expected uncovered demand); they reformulate the multi-objective problem using the ϵ­-constraint method and solve it via a branch­-and­-cut approach. 




Unlike the MCTP's node­-based definition of coverage, 

the close-enough arc routing problem (CEARP) 
defines coverage based upon the proximity of traversed arcs 
in relation to the location of interest <cit.>; we implement both definitions in our study. Since the CEARP is a generalization of the directed rural postman problem (DRPP), <cit.> utilize DRPP heuristics as a solution methodology after identifying arc subsets. <cit.> solve the CEARP exactly using a mixed­-integer programming formulation and cutting-plane approaches. Both <cit.> and <cit.> examine the problem from the perspective of meter reading. Furthermore, <cit.> model the CEARP as a second­-order cone and solve it exactly using a branch-and-bound algorithm; they test their approach on two- and three-dimensional instances. Coverage constraints, moreover, can be time-dependent; thus, it is important to include the speed at which an arc is traversed as a degree of freedom. 






We model the multi-vehicle routing problem under energy, time window, and covering constraints. We formulate a deterministic mixed­-integer second-order cone programming (MISOCP) model to construct a set of maximal coverage routes that service (observe) a secondary set of sites under a fixed time schedule. Our model is an extension of both the MCTP and the VRPTW under energy constraints. Our model's solution accounts for the vehicles' energy limitations and prioritizes risk­-averse routes. We implement our model by examining the surveillance of both low and high priority targets. The remainder of the paper is organized as follows: We define our model and solution approach in the Section <ref> and study a special case which we solve using Lagrangian relaxation in Section <ref>. Next, we discuss computational experiments in Section <ref> and conclude with discussions on future research directions in Section <ref>.





§ MODEL



We seek to route unmanned vehicles for the purpose of target surveillance. In particular, we consider a fleet of vehicles that are used to cover a set of stationary targets. These targets, moreover, can be covered by traveling in proximity to and by idling near them. The process of surveilling targets of interest immediately subjects the fleet to threats of detection. Supplementary problems thus aim to evade enemy threats such as radar detection or combatant encounters (see, e.g., ), and hence seek to minimize the potential risk associated with a vehicle's route. We study the tradeoff between developing routes that are both risk-averse and coverage­-intensive, while accounting for the limitations imposed on the fleet with respect to time, energy, speed, and risk.


Specifically, the model determines (i) the number of vehicles to deploy and their respective routes, (ii) the vehicle speed along every traversed arc, and (iii) the service time, i.e., the time spent surveilling while idling at a given waypoint, such that the minimum level of surveillance required for each target is achieved. Vehicle energy limits are considered and high-valued surveillance targets are prioritized. Coverage radii are defined per vehicle such that a target must pass within the coverage area for some specified time to successfully surveil the respective target. Moreover, risk radii are defined per target such that a vehicle's threat of detection is dependent upon the amount of time spent within some range of a given target.



 §.§ Notation and Model Formulation


We formulate the proposed model over a directed network G=(V∪ W,A) with a set of K vehicles. The node set V = {0,1,…, n,n+1} consists of all waypoints (locations) that can be visited, where 0 and n+1 are the insertion and extraction depots, respectively. We assume that waypoints have been previously determined and are fixed. For notational purposes, we define V^0=V\{0,n+1}. 
Let W be the set of targets and A be the set of arcs in the network. 

Routing decisions are incorporated for each arc/vehicle pair via a binary decision variable x_ijk which equals 1 if vehicle k travels along arc (i,j)∈ A and 0 otherwise. In the context of surveillance, vehicles may only visit particular waypoints during any specific period of availability. This could be due to security rounds or temporary site clearances.
Hence, a time window [a_i,b_i] is associated with each waypoint i∈ V for a given time horizon, where [a_0,b_0] = [a_n+1,b_n+1] represents the earliest departure and the latest arrival times to the insertion and extraction depots, respectively. 
The resulting set of routes can be repeated in each time interval to account for a periodic time horizon. The decision variable s_ik represents the time at which vehicle k arrives at waypoint i∈ V, whereas decision variable y_ik represents the amount of time vehicle k idles at location i∈ V^0. We do not allow surveillance to be conducted at the depots. Idling (or servicing) and waiting are defined as loitering while active and inactive, respectively. We assume that the time a vehicle spends waiting at a waypoint prior to conducting surveillance does not influence the vehicle's risk of detection nor its level of observation, and vehicles can only conduct service at a waypoint i∈ V^0 within its associated time window [a_i,b_i]. Therefore, vehicles can arrive prior to the start of the time window and wait to begin service. Moreover, let 𝐝(i,j) define the distance between nodes i∈ V and j ∈ V∪ W. We assume that each vehicle k travels at a constant speed v_ijk on each arc (i,j)∈ A, so that 𝐝(i,j) = t_ijk× v_ijk, where t_ijk is the travel time for vehicle k∈ K to traverse arc (i,j)∈ A.




  §.§.§ Risk


We assume that each target w∈ W has an associated risk radius, η̂_w, which defines a circular area around it in which a vehicle can be detected (e.g., by a radar detection system installed at the target). Due to a vehicle's ability to travel with different speeds on different arcs, we define the risk in terms of a risk index per unit time. The risk index is directly related to the risk factor, σ_w, and indirectly related to the Euclidean distance between the target and the vehicle (we assume a horizontal plane of motion). The former depends upon the target's detection attributes and is assumed to be a constant. Let [i,j] denote the line segment defined by arc (i,j)∈ A. We define [i^w,j^w]⊆[i,j] as the intersection of line segment [i,j] with the detection area associated with a target w (see Figure <ref> for an illustration).








Thus, as in <cit.>, the associated risk index for a vehicle at location h∈[i^w,j^w] with respect to target w∈ W per unit time is defined as

    r_hw = σ_w1/𝐝(h,w)^2,

where 𝐝(·,·) gives the Euclidean distance between two points (for notational convenience, we use w to denote both the index of a target and its location). Using this definition, we can examine a particular segment [i^w,j^w] and its contribution to the total route's risk. Define P_k as the route for vehicle k∈ K consisting of a sequence of arcs (i,j) traversed by the vehicle.
A vehicle's total risk accumulated from a specific target is: 

    ∑_(i,j)∈ P_k∫_T([i,j]) r_h(t)w dt + ∑_i∈ V^0r_iwy_ik,

where y_ik denotes the amount of time vehicle k idles at waypoint i∈ V^0, T(·) is the interval of time traveling along a line segment [i,j], h(t) is the vehicle's location at time t, and dt is unit time.
Each arc (i,j)∈ P_k accrues risk from target w that is proportional to the segment [i^w,j^w] within the detection area. As defined in <cit.>, the total risk associated with arc (i,j)∈ P_k⊆ A with respect to target w∈ W is:

    ∫_T[i^w,j^w] r_h(t)w dt = ∫_T[i^w,j^w]σ_w1/𝐝(h(t),w)^2 dt = σ_wθ̂_ijw/sinθ̂_ijw𝐝(i^w,w)𝐝(j^w,w)𝐝(i^w,j^w)/𝐝(i,j)t_ijk,

where θ̂_ijw represents the angle between two vectors, from w to i^w and from w to j^w, as illustrated in Figure <ref>.



















Define for each (i,j)∈ A:

    r_ijw := σ_wθ̂_ijw/sinθ̂_ijw𝐝(i^w,w)𝐝(j^w,w),

and for each i∈ V^0:

    r_iw := σ_w1/𝐝(i,w)^2    if 𝐝(i,w) ≤η̂_w,
    
    0     otherwise.

Using these concepts, the risk incurred by a given vehicle k∈ K from target w∈ W over the course of its trajectory, i.e., (<ref>), can be rewritten as:

    ∑_w∈ W(∑_(i,j)∈ Ar_ijw𝐝(i^w,j^w)/𝐝(i,j)  t_ijk +∑_i∈ V^0r_iwy_ik).





  §.§.§ Coverage


We assume that each vehicle k∈ K has an associated coverage area defined by its coverage radius, η̅_k. Ensuring that a target passes within the coverage area during deployment is equivalent to ensuring that a vehicle is routed within some range of the target. We define [i_k^w,j_k^w]⊆[i,j] as the intersection of line segment [i,j] with the coverage area associated with a vehicle k surrounding target w. Thus, we can define the level of surveillance obtained by a vehicle equivalently as the risk of detection of the stationary targets in relation to the mobile vehicle (see Figure <ref>).






Similar to the risk of detection as defined in Section <ref>, we can express the total level of surveillance on target w∈ W obtained by vehicle k∈ K while traversing arc (i,j)∈ A as 

    ρ_kθ̅^w_ijk/sinθ̅^w_ijk𝐝(i_k^w,w)𝐝(j_k^w,w)𝐝(i_k^w,j_k^w)/𝐝(i,j)t_ijk,

where ρ_k is the coverage factor associated with vehicle k∈ K representing its surveillance capabilities, and θ̅^w_ijk represents the angle between two vectors, from w to i^w_k and from w to j^w_k, as illustrated in Figure <ref>. 

Define for each (i,j)∈ A:

    c_ijwk := ρ_kθ̅^w_ijk/sinθ̅^w_ijk𝐝(i_k^w,w)𝐝(j_k^w,w),

and for each i∈ V^0:

    c_iwk := ρ_k1/𝐝(i,w)^2    if 𝐝(i,w) ≤η̅_k,
    
    0     otherwise.


Using these concepts, the total level of coverage on target w obtained by a given vehicle k∈ K throughout its trajectory is calculated by:

    ∑_w∈ W(∑_(i,j)∈ Ac_ijwk𝐝(i_k^w,j_k^w)/𝐝(i,j)  t_ijk +∑_i∈ V^0c_iwky_ik).



In our proposed model, we make the following additional assumptions. First, we assume that a vehicle can cover multiple locations simultaneously as it traverses an arc, i.e., observation times are non-competitive. Second, we allow multiple vehicles to traverse a given arc. Additionally, we assume that a target w∈ W can be surveilled via multiple passes but a vehicle k∈ K can visit any waypoint no more than once during the planning horizon. The former is due to the fact that surveillance need not be continual, while the latter is to avoid additional risk brought by multiple visits.





  §.§.§ Energy



Each vehicle has a limited energy capacity E^max_k, such as fuel, battery, etc. It is well-known that the energy consumption for a vehicle to travel along an arc (i,j)∈ A depends on the speed at which it travels on the arc. For simplicity, we assume that a vehicle's speed must remain constant over a given arc, and we allow for instantaneous change of speed between different arcs. According to <cit.>, the engine power output of a vehicle per unit time can be expressed as a function of its speed: 

    P^e_ijk    = v_ijk(m_vgf_rcos(α)+ 1/2ρ_aC_DA_fv_ijk^2+M_vgsin(α) + m_vδ∂v_ijk/∂t_ijk),
        = v_ijk(β+ γ v_ijk^2+ξ + ψ),
         = v_ijk(β+ γ v_ijk^2),

where β is the rolling resistance, γ v_ijk^2 represents aerodynamic drag, ξ accounts for the power lost due to the slope of a surface, and ψ is the acceleration of the vehicle. 
We let ξ = ψ = 0 since we assume a two-dimensional plane of motion and vehicles are at a steady-state velocity.

Let E_ijk be the energy consumption of vehicle k∈ K over arc (i,j)∈ A. Since E_ijk = P^e_ijkt_ijk, and 𝐝(i,j) = v_ijkt_ijk, it follows that 

    E_ijk    = 𝐝(i,j)β + γ𝐝(i,j) v_ijk^2.

Furthermore, we assume that various power losses due to transmission/drivetrain efficiencies, discharge efficiency, inverter efficiency, motor efficiency, and idling are negligible. We also assume that the vehicle uses no energy (i.e., the vehicle is off) when it waits or idles at a given waypoint.




  §.§.§ Formulation



A summary of the notation is provided below.

Sets:

 

V       Set of nodes that can be visited (V = {0,1,…, n,n+1}) 

V^0     V\{0,n+1}

W       Set of targets that must be covered (W={w_1,…, w_m}) 

K       Set of vehicles (K={k_1,…, k_l}) 

A       Set of arcs in G, where A=(V\{n+1}× V)\({(i,i) : i ∈ V}∪{(0,n+1)})




1pc

Decision Variables:



x_ijk     Indicates whether or not arc (i,j)∈ A is traversed (x_ijk = 1) or not (x_ijk = 0) by vehicle k∈ K

s_ik     Time vehicle k∈ K arrives at i∈ V

t_ijk     Time it takes vehicle k∈ K to traverse arc (i,j)∈ A

v_ijk     Velocity of vehicle k∈ K along arc (i,j)∈ A

y_ik     Time vehicle k∈ K idles at node i∈ V^0




1pc
Parameters:

 -1pc
p0.5inp5.5in
0, n+1       Insertion and extraction depot, respectively

𝐝(i,j)       Distance between nodes i,j∈ V∪ W

𝐝(i_k^w,j_k^w)     Distance along arc (i,j)∈ A in which vehicle k∈ K covers target w∈ W

𝐝(i^w,j^w)     Distance along arc (i,j)∈ A in which vehicle k∈ K risks detection from target w∈ W

p̅_w     Priority factor for target w∈ W

p̂_k     Priority factor for vehicle k∈ K





r_ijw     Risk index for arc (i,j)∈ A associated with target w∈ W

r_iw     Risk index for waypoint i∈ V associated with target w∈ W

c_ijwk     Coverage index for arc (i,j)∈ A associated with target w∈ W  for vehicle k∈ K

c_iwk     Coverage index for waypoint i∈ V associated with target w∈ W for vehicle k∈ K

σ_w     Risk factor for node w∈ W

ρ_k     Coverage factor for vehicle k∈ K






t_w     Minimum allowable time to survey node w∈ W

_k     Maximum speed (≥0) for vehicle k∈ K 

_k     Minimum speed (≥0) for vehicle k∈ K

E^max_k     Maximum energy capacity for vehicle k∈ K

β     Rolling resistance coefficient

γ     Aerodynamic drag coefficient

[a_i,b_i]     Time window for node i∈ V

_ij     Big-M constant where _ij = max{b_i+𝐝(i,j)/τ-a_j, 0} for all (i,j)∈ A




table-1
We optimize the objectives


    max   = ∑_w∈ Wp̅_w∑_k∈ K(∑_(i,j)∈ Ac_ijwk𝐝(i_k^w,j_k^w)/𝐝(i,j)  t_ijk +∑_i∈ V^0c_iwky_ik), 
       min   = ∑_w∈ W∑_k∈ Kp̂_k(∑_(i,j)∈ Ar_ijw𝐝(i^w,j^w)/𝐝(i,j)  t_ijk +∑_i∈ V^0r_iwy_ik),

subject to

    ∑_k∈ K(∑_(i,j)∈ Ac_ijwk𝐝(i_k^w,j_k^w)/𝐝(i,j)  t_ijk +∑_i∈ V^0c_iwky_ik) ≥ t_w,       ∀ w∈ W,
        x_ijk𝐝(i,j) ≤ t_ijkv_ijk,       ∀ (i,j)∈ A, k∈ K,
       ∑_(i,j)∈ A(𝐝(i,j)β x_ijk + γ𝐝(i,j) v_ijk^2) ≤ E^max_k,       ∀ k∈ K,
    
    
        x_ijk≤∑_j'∈ V^0:
     (0,j')∈ A x_0j'k,       ∀ (i,j)∈ A, k∈ K,
       ∑_j∈ V:
    (0,j)∈ A x_0jk = ∑_j∈ V:
    (j,n+1)∈ Ax_j(n+1)k,       ∀ k∈ K,
       ∑_j∈ V:
    (i,j)∈ A x_ijk≤ 1,       ∀ i∈ V, k ∈ K,
    
    
    
       ∑_j∈ V:
     (i,j)∈ Ax_ijk = ∑_j∈ V:
     (j,i)∈ Ax_jik,        ∀ i∈ V^0, k ∈ K,
    
    
    
    
       𝐝(i,j)/_kx_ijk≤ t_ijk≤𝐝(i,j)/_kx_ijk,       ∀ (i,j)∈ A, k∈ K,
       _kx_ijk≤ v_ijk≤_kx_ijk,       ∀ (i,j)∈ A, k∈ K,
        s_ik + (t_ijk+y_ik) - s_jk≤ (1 - x_ijk)_ij,       ∀ (i,j)∈ A, k∈ K,
        a_i∑_j∈ V:
    (i,j)∈ Ax_ijk≤ s_ik,       ∀ i∈ V^0, k∈ K,
        s_ik +y_ik≤ b_i∑_j∈ V:
    (i,j)∈ Ax_ijk,       ∀ i∈ V^0, k∈ K,
        a_0 ≤ s_ik≤ b_n+1,       ∀ i∈{0,n+1}, k∈ K,
    
    
        y_ik = 0,       ∀ i∈{0,n+1}, k∈ K,
        y_ik≥ 0,       ∀ i ∈ V^0, k∈ K,
       x_ijk∈{0,1},       ∀ (i,j)∈ A, k∈ K.




Objective (<ref>) maximizes a weighted sum of coverage times for all targets as some targets are of higher priority than others, while objective (<ref>) minimizes the fleet's risk of detection. Constraint set (<ref>) ensures that the required minimum amount of surveillance is obtained for each target. Constraints (<ref>) relate distance to steady-state speed and travel time; the equality relationship 𝐝(i,j) = v_ijkt_ijk is relaxed to an inequality without affecting the optimal value of (<ref>) to allow for a second-order cone reformulation (see Proposition <ref>). Constraints (<ref>) ensure that each vehicle does not exceed its maximum energy capacity. Constraints (<ref>) ensure that every utilized vehicle is deployed from the insertion depot 0, while constraints (<ref>) make sure that the same number of vehicles that leave the insertion depot arrive at the extraction site n+1. For a given vehicle k∈ K and waypoint i∈ V, the constraint set (<ref>) requires that no waypoint is visited more than once per vehicle. Constraints (<ref>) enforce flow balance and effectively require that vehicles may only start or end their route at a depot. Constraints (<ref>) and (<ref>) enforce the vehicle speed restrictions, while constraints (<ref>) -(<ref>) ensure schedule feasibility with respect to time window constraints. In particular, constraints (<ref>) influence vehicle speeds by limiting the amount of time they can devote toward surveilling any subset of targets. Moreover, although a vehicle can arrive at a given waypoint prior its earliest arrival time, constraints (<ref>) ensure that service (surveillance) starts within the time window. Constraints (<ref>) define y_ik as the idling time at waypoint i∈ V for each vehicle k∈ K and do not allow a vehicle to conduct surveillance at a waypoint beyond the time window; constraints  (<ref>) and  (<ref>) restrict the total idle time to be at most the time window's length. 
Constraints (<ref>) limit total time of each vehicle's route and constraints (<ref>) do not allow service to be conducted at either the insertion or extraction sites. Finally, constraint sets (<ref>) and (<ref>) are the non-negativity and integrality restrictions on the decision variables, respectively. We have assumed an identical fleet of vehicles in our experiments. That is, all vehicles have the same surveillance range (η̅), importance factor (), coverage factor (ρ), energy capacity (E^max), and boundary speeds ( and ). 


The second-order cone constraint (<ref>) will be active for at least one optimal solution. 
 


Assume for contradiction that for any optimal solution (x^*, s^*, t^*,v^*, y^*) of (<ref>), there exists some k∈ K and some (i,j)∈ A such that x^*_ijk = 1 and t_ijk^*v_ijk^* > 𝐝(i,j)x^*_ijk. In particular, note that

    _kx^*_ijk≤ v^*_ijk≤_kx^*_ijk,   ∀ (i,j)∈ A, k∈ K,

 and
 
    ∑_(i,j)∈ A𝐝(i,j)β x^*_ijk + γ𝐝(i,j) v^*_ijk^2 ≤ E^max_k.

 Moreover, let _ijk = 𝐝(i,j)/t^*_ijk such that t_ijk^*_ijk = 𝐝(i,j). It follows then that _ijk < v_ijk^* and so
 
    ∑_(i,j)∈ A𝐝(i,j)β x^*_ijk + γ𝐝(i,j) _ijk^2 ≤ E^max_k.

Further, since 

    𝐝(i,j)/τ̅x^*_ijk≤ t_ijk^*≤𝐝(i,j)/τx^*_ijk,

 we know 
 
    _kx^*_ijk≤_ijk≤_kx^*_ijk.

Thus, there exists an alternative optimal solution in which t_ijk^*_ijk = 𝐝(i,j)x^*_ijk for all (i,j)∈ A and for all k∈ K.









 §.§ Solution Methodology





Model (<ref>) is a bi-objective optimization problem where the fleet's ability to gather information is hindered by
the risk accumulated during deployment. Thus, our intention is to derive a Pareto curve which depicts the tradeoff between the two competing objectives. The efficient frontier characterizes the set of non-dominated solutions and can be generated via several methodologies (see, e.g., ). We use the ϵ-constraint method <cit.>, which can generate all possible efficient solutions in the absence of convexity so long as ϵ is chosen appropriately <cit.>. We chose to bound the total risk (objective (<ref>)) with ϵ-constraints. That is, we solve


    max    , 
    s.t.    ≤ϵ, 
       constraints (<ref>)–(<ref>),



for a given value of ϵ that admits a feasible solution under constraints (<ref>)–(<ref>). Due to the continuous nature of our objective functions, we aim to approximate the Pareto frontier by solving model (<ref>) iteratively. Thus, given the total risk associated with the n^th iteration (^n) and a pre-defined step­-size δ, we define

    ϵ^n+1←^n + δ.

We can find a sufficient number of non-dominated points by choosing δ small enough. However, we cannot guarantee we generate the entire Pareto curve. 





§ SPECIAL CASES




In this section, we consider several special cases of the proposed model and develop specialized solution approaches for these cases. Consider the following simplifying assumptions: 1) the fleet's energy supply and the time windows associated with each waypoint can be approximated by an operational deadline (T = b_n+1 - a_0), and 2) vehicles are not subject to risk of detection, i.e., the risk constraint (<ref>) is relaxed. Specifically, the above assumptions can be fulfilled by letting [a_i, b_i] = [0,∞) for all i ∈ V^0 and E^max_k = ∞ for all  k∈ K, and introducing constraints:

    ∑_(i,j)∈ At_ijk + ∑_i∈ V^0y_ik≤ T,   ∀ k∈ K.

We refer to this new model as model (). Without the constraints on the risk of detection (<ref>), model () is loosely coupled in the sense that each vehicle k∈ K can be treated independently if constraints (<ref>) are relaxed. This motivates a Lagrangian relaxation based approach to solve model (). For each w∈ W, given a Lagrangian multiplier λ_w associated with constraint (<ref>) for target w, the Lagrangian relaxation of model () is formulated as:


    max     - ∑_w ∈ Wλ_w(∑_k ∈ K(∑_(i,j) ∈ A c_ijwk𝐝(i_k^w,j_k^w)/𝐝(i,j) t_ijk + ∑_i ∈ V^0 c_iwky_ik)-t_w), 
    s.t.    constraints (<ref>)–(<ref>), (<ref>).


The Lagrangian relaxation (<ref>) becomes separable by k∈ K for each fixed vector λ of Lagrangian multipliers.  We define  to be the set of all feasible paths from node 0 to node n+1. Due to the homogeneity assumption of the fleet, define _ijw := 𝐝(i^w_k,j^w_k), c_iw := c_iwk, c_ijw := c_ijwk, for all k∈ K. Hence, given a fixed λ∈|W|, (<ref>) can be reformulated as:


    f(λ) = ∑_w∈ Wt_wλ_w +    max_P∈∑_w∈ W(p̅_w-λ_w)[∑_(i,j)∈ Pc_ijw_ijw/𝐝(i,j)  t_ij + ∑_i∈ V^0c_i wy_i],
    s.t.   𝐝(i,j)/τ̅≤ t_ij≤d(i,j)/τ,       ∀ (i,j)∈ P,
    
        y_i = 0,       ∀ i∈{0,n+1},
        y_i≥ 0,       ∀ i ∈ V^0.




The tightest Lagrangian relaxation bound, namely the Lagrangian dual bound, can be found by optimizing over the Lagrangian vector λ∈|W|. Optimal solutions to the Lagrangian relaxation (<ref>) will produce cutting planes that will be essential to solving the Lagrangian dual problem. We consider a standard level bundle method <cit.>, a variant of the cutting plane method for solving the Lagrangian dual problem. Specifically, this method keeps track of a lower bound and an upper bound to the Lagrangian dual value, denoted as LB and UB, respectively, and updates their values in an iterative fashion until they converge to the same value, which corresponds to the Lagrangian dual value. Suppose λ̂^L is the iterate solution and suppose f^lev_L is the level parameter at the L-th iteration of the cutting plane method, the Lagrangian dual master problem is given by:


    min     λ -λ̂^L^2 
    s.t.     θ≤ f^lev_L,
       θ≥∑_w∈ W[t_w - α̂_w(λ̂^ℓ)]λ_w + ∑_w∈ W_wα̂_w(λ̂^ℓ),    ∀ℓ = 1,2,…, L,
       λ≤ 0,



where cut coefficient α̂_w(λ̂^ℓ) is obtained by solving the Lagrangian relaxation (<ref>) for Lagrangian vectors λ̂^ℓ∈|W|,  ℓ = 1,2,…, L. Solving this Lagrangian dual master problem will produce a new iterate λ̂^L+1, on top of which a cut of the form (<ref>) will be generated by solving the corresponding Lagrangian relaxation problem (<ref>) with λ̂^L+1, which will then be added to the master problem (<ref>) in the next iteration. The level parameter f^lev_L is defined as f^lev_L: = ϕ LB + (1-ϕ) UB for some ϕ∈ (0,1). UB is updated each time a Lagrangian relaxation problem (<ref>) with respect to a new iterate is solved: UB ←min{UB,f(λ^ℓ)}, and LB is updated each time where the Lagrangian master problem (<ref>) is infeasible, indicating that LB can be updated to be identical to the level parameter f^lev_L. 



  
Optimal solution structures of the Lagrangian relaxation problem (<ref>). We show that an optimal solution to the Lagrangian relaxation problem (<ref>), which will be iteratively solved within the level bundle method for computing the Lagrangian dual value, exhibits certain combinatorial structures that allow us to develop specialized solution approaches. 


There exists an optimal solution to problem (<ref>) where a vehicle will idle at no more than one waypoint, i.e., among all waypoints i∈ V^0, there exists at most one waypoint v̅ where y_v̅ > 0.


Assume for contradiction that for any optimal solution to (<ref>), there exists a set I of waypoints such that y_i > 0, i∈ I, and |I| ≥ 2. Let ∈ I be such that c_ w∑_w∈ W(p̅_w-λ_w) ≥ c_iw∑_w∈ W(p̅_w-λ_w) for all i∈ I\{}. The solution could be improved by “reallocating” all the idling time to waypoint , while the travel times t_ij for all arcs (i,j) traversed by the route corresponding to the solution are fixed. Thus, the assumed solution must not have been optimal, which yields a contradiction. In the case where there exist multiple , an alternative optimal solution can be obtained by idling only at one , which also yields a contradiction.

The optimal solution structure outlined in Proposition <ref> only holds for the Lagrangian relaxation problem (<ref>), but not for the original model (). Consider the following example in which there are two targets where each has only one waypoint in range. Further, suppose the waypoints are located at the boundary of the coverage area so that no coverage can be obtained by traversing any arc, and also assume that each target requires a positive amount of coverage (see constraint set (<ref>)). In this case, to meet the coverage level a vehicle must idle at both waypoints.

Next, we develop a specialized algorithm for solving the Lagrangian relaxation problem (<ref>), which exploits the aforementioned optimal solution structures. We will start with a simpler case where T≥max_P∈∑_(i,j)∈ P𝐝(i,j)/, i.e., the operational deadline is irrelevant as any path is feasible under this time limit even traveling with the slowest possible speed. 



 §.§ Special Case I: A Sufficiently Large Operational Deadline

In this case, since the operational deadline is sufficiently large, i.e., T≥max_P∈∑_(i,j)∈ P𝐝(i,j)/, given a path P and the travel time on each arc (i,j)∈ P, if the vehicle ever idles at any waypoint, it will spend ∑_i∈ V^0y_i = T - ∑_(i,j)∈ P t_ij on exactly one waypoint along path P according to Proposition <ref>. In the following approach, we will loop through all candidate waypoints ∈ V^0 where a vehicle may benefit from idling (to gain a positive coverage), i.e., ∑_w∈ W(p̅_w-λ_w) c_ w > 0, and find the best path that traverses this waypoint. We must also consider the case where the vehicle does not idle at all. After we finish the loop, the best path among all these paths is selected.  As an overview, to solve (<ref>), we use the following procedure: (i) For each ∈ I:={i∈ V^0 |∑_w∈ W(p̅_w-λ_w) c_i w > 0}∪{0}, find the maximal covering route P^*_ among all feasible routes _ that travel through . We denote the case when no idling occurs as = 0 and define c_0w := 0, ∀ w∈ W. Arc weights represent the coverage obtained while traveling along each respective route segment, whose values are known upon fixing ; 
ii) Select the optimal path (and the associated travel and idle time along the path) from ∪_∈ I{P^*_}.

Given a fixed waypoint  where the vehicle idles, model (<ref>) becomes:


    f_(λ) = ∑_w∈ Wt_wλ_w +    max_P∈_∑_w∈ W(p̅_w-λ_w)[∑_(i,j)∈ Pc_ijw_ijw/𝐝(i,j)  t_ij + c_ wy_],
       s.t. 𝐝(i,j)/τ̅≤ t_ij≤𝐝(i,j)/τ,        ∀ (i,j)∈ P,
        y_≥ 0,       ∀ ∈ V^0.



Since y_ = T - ∑_(i,j)∈ Pt_ij, let P^*_ and t^*_ij for all (i,j)∈ P^*_ 
denote an optimal solution to model (<ref>), we can rewrite its objective function (<ref>) as:

    f_(λ) =    ∑_w∈ Wt_wλ_w +∑_w∈ W(p̅_w-λ_w)T c_ w + max_P∈_∑_(i,j)∈ P[∑_w∈ W(p̅_w-λ_w)(c_ijw_ijw/𝐝(i,j) - c_ w)t_ij],
        = ∑_w∈ Wt_wλ_w +∑_w∈ W(p̅_w-λ_w)T c_ w + ∑_(i,j)∈ P^*_[∑_w∈ W(p̅_w-λ_w)(c_ijw_ijw/𝐝(i,j) - c_ w)t^*_ijv̅],
        =  ∑_w∈ Wt_wλ_w +∑_w∈ W(p̅_w-λ_w)[T c_ w + ∑_(i,j)∈ P^*_(c_ijw_ijw/𝐝(i,j) - c_ w)t^*_ijv̅]_α̂_w(λ),
        = ∑_w∈ W[t_w - α̂_w(λ)]λ_w + ∑_w∈ W_wα̂_w(λ),

where given v̅∈ I, the optimal travel time t^*_ijv̅ along each arc (i,j)∈ P^*_ corresponds to either its lower or upper bound. Specifically, for each (i,j)∈ P^*_:

  1) If ∑_w∈ W(p̅_w-λ_w) (_ijw/𝐝(i,j)   c_ijw - c_ w) > 0, then t^*_ijv̅ = 𝐝(i,j)/;

  2) If ∑_w∈ W(p̅_w-λ_w) (_ijw/𝐝(i,j)   c_ijw - c_ w) ≤ 0, then t^*_ijv̅ = 𝐝(i,j)/.


Thus, we see that given a fixed waypoint  for which to idle, the speed at which an arc will be traversed immediately
follows, regardless of the choice of path P_. Therefore, for each arc (i,j)∈ A, we can compute the optimal adjusted coverage as:

    c^*_(i,j),v̅ := ∑_w∈ W(p̅_w-λ_w) (_ijw/𝐝(i,j)   c_ijw - c_ w)t^*_ijv̅.
  
This simple observation allows us to perform path/route optimization using a pre-optimized travel time and the associated adjusted coverage for each arc (i,j)∈ A as defined above, so that the label correcting algorithm for standard VRP (see, e.g., <cit.>) can be applied to solve model (<ref>) with some simple extensions elaborated below. We remark that this is not the case when the operational deadline T is not sufficiently large.



  §.§.§ The labeling algorithm

Algorithm <ref> is an extension of the label correcting algorithm for the VRPTW <cit.>. Specifically, each label L represents a partial path P (i.e., L=(P)) starting from 0 and ending at some node v∈ V^0. We define M(L) as the set of waypoints visited so far by L. For any fixed ∈ I, we define c_v̅(L) as the cumulative adjusted coverage associated with L, given by c_v̅(L) = ∑_(i,j)∈ P c^*_(i,j),v̅.

An optimal complete path (starting from node 0 and ending at node n+1) is obtained by performing partial enumeration of labels (e.g., partial paths) via a dynamic programming algorithm that iteratively generates labels by extending existing labels to new waypoints. As the operational deadline is non-binding for any path, a label can be extended to visit any waypoint not already visited.



An effective dominance rule is critical to the success of the labeling algorithm, which allows many partial paths to be dominated (and thus be removed) from the dynamic programming procedure. We next describe the dominance rule that we develop for special case I.

[Dominance rule for special case I]
We say that a label L^2 is dominated by another label L^1, if: (i) L^1 and L^2 end at the same waypoint; (ii) M(L^1)⊆ M(L^2); and further:


  1. If ∈ M(L^1) and ∈ M(L^2) or if ∉ M(L^1) and ∉ M(L^2), then L^1 dominates L^2 if c_v̅(L^1) ≥ c_v̅(L^2), 

  2. If ∉ M(L^1) but ∈ M(L^2), and ∈ V^0, we compute the worst case extra cost incurred by adjusting a direct application of an extension L' of L^2 to L^1 in order to accomodate the visit to . This extra cost can be computed as follows:

    c^extra_v̅(L^2) = min_i∈ V^0:
     i∉ M(L^2){c^*_(i,v̅),v̅ + c^*_(v̅,n+1),v̅ - c^*_(i,n+1),v̅}.

If c_v̅(L^1)+c^extra_v̅(L^2) ≥ c_v̅(L^2), then L^1 dominates L^2.


Note that it is not possible for ∈ M(L^1) but ∉ M(L^2) since condition (ii) must be satisfied in order for L^1 to dominate L^2. We define c^extra_0(L^2) := 0.






  1. In this case, any feasible extension L' of L^2 to a complete path is also a feasible extension of L^1 to a complete path. Also, the cumulative adjusted coverage of the complete path corresponding to L^1, c_v̅(L^1⊕ L'), is higher than that of L^2, c_v̅(L^2 ⊕ L'), since c_v̅(L^1) ≥ c_v̅(L^2) and the extensions for L^1 and L^2 are identical.


  2. In this case, L^1 is not yet a feasible path (since it has not visited v̅ yet) and so the argument for case 1 cannot be made. According to conditions (i) and (ii), any feasible extension L' of L^2 is also a feasible extension of L^1. Suppose that node i∈ V^0, i∉ M(L^2) is the last node visited by L^2 prior to n+1 for this feasible extension L'. When we apply the same extension to L^1, since L^1 has not visited v̅ yet, we will have to use arcs (i,) and (,n+1) in the extension for L^1 instead of arc (i,n+1). As we do not know which node i∈ V^0, i∉ M(L^2) is the last node visited by L^2 prior to n+1, we must consider the worst-case scenario in order to be sufficient for L^1 to dominate L^2. This worst case corresponds to the expression shown in (<ref>): if the cumulative adjusted coverage of L^1 is greater than or equal to that of L^2 in this worst case, then we can conclude that L^2 is dominated by L^1.




















































 §.§ Special Case II: A Restrictive Operational Deadline

We now consider the case where the operational time restriction can limit the path choices. In this case, there exist paths with traversal times potentially longer than the operational time limit, i.e., T < max_P∈∑_(i,j)∈ P𝐝(i,j)/. Similar to special case I, we have:

    f(λ) = ∑_w∈ Wt_wλ_w +max.75![c]
    max_P∈   ∑_(i,j)∈ P[∑_w∈ W(p̅_w-λ_w)c_ijw_ijw/𝐝(i,j)] t_ij ,
    s.t.   ∑_(i,j)∈ Pt_ij≤ T,
    s.t.   𝐝(i,j)/τ̅≤ t_ij≤𝐝(i,j)/τ,      ∀ (i,j)∈ P,
       (a)
    [c]
    max_∈ V^0   {∑_w∈ W(p̅_w-λ_w)Tc_ w +max_P∈_∑_(i,j)∈ P[∑_w∈ W(p̅_w-λ_w)(c_ijw_ijw/𝐝(i,j) - c_ w)]t_ij} ,
    s.t.   ∑_(i,j)∈ Pt_ij≤ T,
    s.t.   𝐝(i,j)/τ̅≤ t_ij≤𝐝(i,j)/τ,      ∀ (i,j)∈ P.
       (b)


In either case (a) or (b), for any fixed path P∈ or P∈_, the underlying problem is a continuous knapsack problem. Therefore, in an optimal solution, all decision variables will be set to either their lower or upper bounds, save for at most one <cit.>. Moreover, the assignment of different values to these variables depends on the ranking of their respective coverage value per unit distance. For notational convenience, we denote f_(i,j),v̅ := ∑_w∈ W(p̅_w-λ_w)(c_ijw_ijw/𝐝(i,j) - c_ w) for any fixed ∈ I, for each (i,j)∈ A.

[Optimal travel time structures on a fixed path]
Given a fixed ∈ I, let P∈_ be a fixed path that traverses , then there is an optimal solution {t^*_ij}_(i,j)∈ P to the following optimal timing problem:


    max    ∑_(i,j)∈ P f_(i,j),v̅t_ij
    s.t.   ∑_(i,j)∈ Pt_ij≤ T 
       𝐝(i,j)/τ̅≤ t_ij≤𝐝(i,j)/τ,      ∀ (i,j)∈ P,


such that ∃ a^*∈ P where t^*_a^*∈ [𝐝(a^*)/τ̅,𝐝(a^*)/τ], and:


  * For all (i,j)∈ P, (i,j)≠ a^* such that f_(i,j),v̅/𝐝(i,j)≥f_a^*,v̅/𝐝(a^*), t^*_ij = 𝐝(i,j)/τ.

  * For all (i,j)∈ P, (i,j)≠ a^* such that f_(i,j),v̅/𝐝(i,j)≤f_a^*,v̅/𝐝(a^*), t^*_ij = 𝐝(i,j)/τ̅.



We will exploit special optimal solution structures shown in Proposition <ref> to develop a novel labeling algorithm to solve (<ref>). In particular, we introduce the concept of a token, which can be used to allow the time spent on traversing an arc (i,j)∈ P to be strictly between their respective upper and lower bounds. According to Proposition <ref>, for each route we are allowed to use at most one token, and except for the arc where the token is applied, all other arcs will be traversed with either the minimum or the maximum possible travel time. Of course, the challenge in solving the optimal timing problem together with optimal path selection as in (<ref>) is that, before a route is completed, it is unclear whether or not the token will be applied to any arc, and if so, where the token will be applied. However, as indicated by Proposition <ref>, the place where the token will be applied is restricted by arcs that have been traversed by a label (partial path). Specifically, given a label L and an arc (i,j) on which we perform a label extension, we need to consider up to two different possible label extensions from L = (P) to (P ⊕ (i,j)):

  1. If f_(i,j),v̅≤ 0, then we must extend the label by traversing arc (i,j) with the minimum time. 

  2. If f_(i,j),v̅ > 0, i.e., traversing arc (i,j) with the maximum time is optimal, then we may extend the label by traversing arc (i,j) with either the minimum time or the maximum time. We will refer to these arcs as the tradeoff arcs.
We denote the set of tradeoff arcs associated with a label L as the tradeoff set (L). It is clear that we will only use a token on a tradeoff arc. For each tradeoff arc in (L), we compute its corresponding ratio f_(i,j),v̅/𝐝(i,j). Given a label (partial path) L, we keep track of the minimum ratio f_min,(L) among all tradeoff arcs within (L) that have been traveled with the maximum time, and the maximum ratio f_max,(L) among all tradeoff arcs within (L) that have been traveled with the minimum time (see Figure <ref>).


Specifically, we label a tradeoff arc (i,j)∈(L) with label u_(i,j) = 1 if its travel time is set to its maximum possible travel time, and u_(i,j) = 0 if its travel time is set to its minimum. When a label L is to be extended with a tradeoff arc (i,j), we first compute its corresponding ratio f_(i,j),v̅/𝐝(i,j), and then we compare it with the minimum and maximum ratio that we keep track of for label L to decide which label (0 or 1) will be assigned to arc (i,j).


  * If f_(i,j),v̅/𝐝(i,j)≥ f_min,(L) := min{f_(i',j'),v̅/𝐝(i',j')| (i',j') ∈(L) such that u_(i',j') = 1}, then we will label arc (i,j) with u_(i,j) = 1 and extend label L with maximum travel time on arc (i,j).



  * If f_(i,j),v̅/𝐝(i,j)≤ f_max,(L) := max{f_(i',j'),v̅/𝐝(i',j')| (i',j') ∈(L) such that u_(i',j') = 0}, then we will label arc (i,j) with u_(i,j) = 0 and extend label L with minimum travel time on arc (i,j).


  * If neither of the above two cases occurs, we will create two new labels where one is extended with the maximum travel time and the other is extended with the minimum travel time. In this case, f_min,(L) and f_max,(L) after the extension will be updated accordingly.


Note that the feasibility of each label extension is verified before we execute the extension. For any fixed ∈ I, let t_(L) be the arrival time at the last node visited by L, given by 

    t_v̅(L) = ∑_(i,j)∈ P, (i,j)∉(L)𝐝(i,j)/τ̅ + ∑_(i,j)∈ P, (i,j)∈(L)(𝐝(i,j)/τ̅(1-u_(i,j)) + 𝐝(i,j)/τu_(i,j)).

We will check:


  * If ∈ M(L), a label L ending at waypoint i∈ V^0 can be extended to j∈ V^0, where j∉ M(L), only if t_(L) + 𝐝(i,j)/τ̅(1-u_(i,j)) + 𝐝(i,j)/τu_(i,j) + 𝐝(j,n+1)/τ̅≤ T, i.e., a feasible complete path is possible if using the maximum speed on the last leg of the path, (j, n+1).  

  * If ∉ M(L), a label L ending at waypoint i∈ V^0 can be extended to j∈ V^0, where j∉ M(L), only if t_(L) + 𝐝(i,j)/τ̅(1-u_(i,j)) + 𝐝(i,j)/τu_(i,j) + 𝐝(j,v̅)/τ̅ + 𝐝(v̅,n+1)/τ̅≤ T, i.e., a feasible complete path is possible if using the maximum speed on the last two legs of the path, (j, v̅) and (v̅, n+1).


Except for the possibility of creating up to two new labels per label extension, the labeling algorithm for special case II is identical to Algorithm <ref>. Together with the cumulative adjusted coverage c_(L), the arrival time t_(L), we also store the maximum cumulative adjusted coverage c^max_(L), the minimum cumulative adjusted coverage c^min_(L), and their corresponding arrival times t^max_(L) and t^min_(L), respectively, by applying the token on either the arc corresponding to f_min,(L) or the arc corresponding to f_max,(L). We next present the dominance rules that we develop for the labeling algorithm for special case II.

[Dominance Rule for Special Case II]
We say that a label L^2 is dominated by another label L^1, if: (i) L^1 and L^2 end at the same waypoint; (ii) M(L^1)⊆ M(L^2); (iii) c_(L^1) ≥ c_(L^2), t_(L^1) ≤ t_(L^2), and at least one of these holds strictly, and furthermore, one of the following four cases will hold:


  * Case 1: c_v̅(L^1) ≥ c^max_v̅(L^2), and t_v̅(L^1) ≤ t^min_v̅(L^2);

  * Case 2: If c_v̅(L^1) ≤ c^max_v̅(L^2), and t_v̅(L^1) ≤ t_v̅^min(L^2), then t_v̅(L^1) + (c_v̅^max(L^2)-c_v̅(L^1))/f_min,v̅(L) ≤ t^max_v̅(L^2) must hold;

  * Case 3: If c_v̅(L^1) ≥ c^max_v̅(L^2), and t_v̅(L^1) ≥ t^min_v̅(L^2), then c_v̅(L^1) - f_max,v̅(L)·(t_v̅(L^1) - t^min_v̅(L^2)) ≥ c^min_v̅(L^2) must hold;

  * Case 4: If c_v̅(L^1) ≤ c^max_v̅(L^2), and t_v̅(L^1) ≥ t^min_v̅(L^2),  then both t_v̅(L^1) + (c^max_v̅(L^2)-c_v̅(L^1))/f_min,v̅(L) ≤ t^max_v̅(L^2) and c_v̅(L^1) - f_max,v̅(L)·(t_v̅(L^1) - t^min_v̅(L^2)) ≥ c^min_v̅(L^2) must hold.






   
Figure <ref> illustrates these four cases, where each label L corresponds to a piecewise linear function with two pieces. The upper segment of each function corresponds to using the token to change the traversal time of the arc corresponding to f_max,(L). As its time can vary between its upper and lower bounds, and in doing so t_(L) and c_(L) increase, the resulting set of solutions are characterized by a line segment. A similar argument is applied to the lower segment associated with changing the traversal time of the arc corresponding to f_min,(L). We remark that the slope of the lower piece is steeper than that of the upper piece due to the optimal solution structure shown in Proposition <ref>.

 
Consider any feasible extension L' of L^2. We discuss the following two situations.

First, if the token is used in the extension L' after L^2 is extended to a complete path, it is feasible to apply the same extension and use the token at the same arc for L^1, and the resulting cumulative adjusted coverage for L^1 is higher than L^2 according to conditions (i), (ii) and (iii).

Second, if the token is used on an existing arc traversed by label L^2, in order to be feasible for L^1 to apply the same extension L' of L^2, and yet yield a higher cumulative adjusted coverage, we consider the four cases discussed in Proposition <ref>. Case 1 (Fig.  <ref>) is clearly trivial, since point (t_(L^1),c_(L^1)) dominates both extreme situations for L^2 by utilizing a token, that is, c_(L^1) ≥ c^max_(L^2), and t_(L^1) ≤ t^min_(L^2). For simplicity, we only discuss case 2 in detail below, since case 3 can be shown using a similar argument, and case 4 is a combination of case 2 and case 3.

In fact, conditions imposed in case 2 ensures that the intersection of the upper piece of L^1 with c = c^max_(L^2) corresponds to a time that is less than t^max_(L^2). Therefore, both ends of the upper piece of L^2 are dominated by some point at the upper piece of L^1 (in terms of both cumulative adjusted coverage and time), indicating that for any (c'_(L^2),t'_(L^2)) combination associated with L^2 that is achievable by applying a token on L^2, there exists a combination (c'_(L^1),t'_(L^1)) associated with L^1 by applying a token on L^1 such that c'_(L^1) ≥ c'_(L^2) and t'_(L^1) ≤ t'_(L^2).





§ NUMERICAL EXPERIMENTS


We now demonstrate our model using an illustrative example (Section <ref>) and further investigate its efficacy with larger problem instances (Section <ref>). Next, Section <ref> examines solution approaches for the special cases discussed in Section <ref>. In our experiments, we set β = 1 and γ = 1 for the rolling resistance and aerodynamic drag coefficients, respectively, and fix each vehicle's speed range such that [,] = [1,10]. The model in Section <ref> was coded in AMPL and solved via Gurobi 8.0.1.
In Section <ref>, computations were performed using an 8­-thread PC running an Intel i7 5960X 3.7 GHz processor with 128 GB RAM and instances were analyzed using Gurobi 8.0.1 with a 3600 seconds time limit. For Section <ref> computations, an Intel i7-8700 3.2 GHz processor with 32 GB of RAM was used. Gurobi 8.1.1 was run in a Python 2.7 environment with a time limit of 7200 seconds. Both the level bundle method (for solving the Lagrangian dual) and the labeling algorithms (for solving the Lagrangian relaxation outlined in Section <ref>) were implemented in C⁠[4].4ex++. We use a 7200 second time limit for the level bundle method.



 §.§ Illustrative Examples



We 
consider a symmetric network composed of four waypoints ([YellowOrange,fill=YellowOrange] (0,0) circle (.6ex);) and eight targets ([ForestGreen,fill=ForestGreen] (0,0) rectangle (1ex,1ex);), where the insertion and extraction depots (11) are identical (Figure <ref>). Coverage radii ([red,thick,dash pattern=on 7pt off 2pt on 1pt off 3pt] (0,0) – (0.85,0);) and risk radii ([black, thick,densely dotted] (0,0) – (0.85,0);) are consistent for each vehicle and target, respectively. 






For the instance parameters described in Table <ref>, subscripts are omitted as values are made constant over their respective sets. Initially, we do not permit idling at any waypoint. 
























We perform sensitivity analyses on key parameters: energy capacity, time windows, risk and coverage radii, and risk and coverage levels. We show the impact each parameter has on the solution by varying its value while holding all others constant. It is important to highlight that an optimal solution will utilize the entire fleet. In fact, in the absence of restrictive parameters (e.g., time windows, energy limits, and risk thresholds), model (<ref>) will produce identical routes for each vehicle (see Figure <ref>). Decreasing vehicle energy capacity decreases accumulated surveillance since each vehicle must now shorten its route and cover different targets independently (Figure <ref>). Moreover, by restricting the target's time window for observation, vehicles are assigned different routes to maximize coverage as they must travel more quickly to maintain schedule feasibility (Figure <ref>). When permitted to service waypoints, vehicles can increase the level of observation; however, due to the incurred risk restriction, the fleet cannot idle indefinitely (Figure <ref>). Decreasing the risk threshold, though, directly reduces the fleet's level of observation (Figure <ref>). Increasing the detection area increases the accumulated risk; hence, a vehicle must sacrifice surveillance time in order to satisfy its permissible risk level. By expanding the coverage area, observation increases as a greater proportion of the route covers targets and vehicles travel through detection areas more quickly  (Figure <ref>). The speed of travel does not diminish the level of surveillance due to the targets being observable over a greater proportion of the routes. As the model aims to maximize total coverage, decreasing the minimum surveillance level has no effect on the solution, while increasing this requirement can lead to infeasibility (Figure <ref>).


   


As certain restrictions become more or less binding, we not only see routes change, but travel speed as well. Figure <ref> shows the speed of each vehicle over the duration of the trip subject to baseline control parameters (Figs. <ref> and <ref>), decreased energy capacity (Figs. <ref> and <ref>), shortened time windows (Figs. <ref> and <ref>), and permissible service (Figs. <ref> and <ref>). The time a vehicle starts service at each waypoint is denoted as a data point, whereas the time spent conducting service at a waypoint is shown as a break in the plot.



Depending upon the limiting resource, the vehicles will either utilize a shorter route, travel at a faster speed, or idle as long as possible. When permissible, vehicles will behave identically, though their paths and speeds will differentiate when resources are limiting or targets are dispersed. In particular, limiting the available energy supply causes a vehicle's total route time to decrease, but the speed it travels remains the same (Figs. <ref> and <ref>). This decrease in route time creates a corresponding reduction in surveillance. Shortening time windows yields a similar result. However, this causes vehicles to travel faster along certain route segments (Figs. <ref> and <ref>). When compared to the effect of decreased energy supply, we see that the increase in speed does not affect the objective as much as the decrease in total trip time; this is due in part to the associated coverage index per time unit. When vehicles are able to service a given waypoint, they will do so for as long as possible such that they can be routed back to the depot within the acceptable risk level (Figs. <ref> and <ref>).



 §.§ Larger Problem Instances


Model (<ref>) was coded in AMPL and validated using instances generated from the well-known Solomon VRPTW benchmarks <cit.>. In particular, we created two distinct networks (Network 1 and Network 2) composed of 20 waypoints (including the depots) and 20 targets from the R101 problem set. The location and time windows for the waypoints are set as the first 20 customers (Customer No. 0­-19) for Network 1 and the second 20 customers (Customer No. 20­-39) for Network 2. Targets are chosen randomly from the remaining 80 customers. The associated time windows for the insertion and extraction depots are redefined for each instance such that a_0 = 0 and b_n+1 = 10000. Coverage radii for each vehicle are set to 10 units and the associated risk radius per target is fixed at 5 units. All targets that are not coverable within the given network or lie directly on the arc (i,j)∈ A are removed for feasibility and numerical purposes, respectively. 

Table <ref> shows the common instance parameters for each network. Targets were randomly assigned a priority weight such that _w∈[1,5] and all vehicles were weighted equally (i.e., _k = 1, ∀ k∈ K).










 















It is important to note that Network 1 consists of instances solvable with two vehicles, whereas Network 2 requires a fleet size of at least seven vehicles. The energy limits were chosen based upon experimental results from initial testing and are non-binding for the given networks. Table <ref> summarizes the results.



All instances from Network 1 were solved to optimality within a few minutes; however, Network 2's instances terminated upon reaching the pre­-defined time limit. While most instances have an optimality gap of less than 1%, instance 2­-4 is an exception. Unlike the other Network 2 instances, the risk threshold for instance 2­-4 is not binding upon termination. In fact, no feasible solution exists with incurred risk greater than 200 units. By directing Gurobi to focus on the solution's lower bound (i.e., set mipfocus = 3) for instance 2­-4, the optimality gap decreases to 1.4% in the same 3600 second time limit. 
Further, Table <ref> confirms the number of branch­-and­-bound nodes explored by Gurobi is far greater for Network 2's instances than for Network 1's instances, which is due in part to both network structure and fleet size.

As outlined in Section <ref>, we can approximate the Pareto frontier 
to illustrate the tradeoff between our two competing objective functions (Eqs. (<ref>) and (<ref>)). Without loss of generality, we generate an approximation to the Pareto curve for instance 1­-2 in Figure <ref> by setting δ = 10^-1 (doing so for the remaining instances would show similar insights). 
By choosing a smaller value of δ, a more accurate representation of the true efficient frontier would result; however, due to the continuous nature of the risk function, it is impossible to generate the full Pareto frontier exactly. This methodology will produce dominated and/or weakly-dominated solutions (which can be removed from the frontier) when energy limits and time windows are more restrictive.





 §.§ Lagrangian Relaxation Based Approach for Special Cases

We demonstrate our specialized algorithm (Section <ref>) for solving the Lagrangian relaxation problem (<ref>) (special case I) and  (<ref>) (special case II) by comparing the Lagrangian dual bound with the solution obtained by Gurobi for three sets of instances: 1) small (|V^0| = 9, |W|≈ 10, η̅ =10), 2) medium (|V^0| = 12, |W|≈ 11, η̅ =20), and 3) large (|V^0| = 15, |W|≈ 12, η̅ =20). All instances are generated using a similar approach to those discussed in Section <ref> and utilize a single depot and one vehicle. 



  §.§.§ Special Case I

Total operational time is defined as T := |A|×max{𝐝(i,j) | (i,j)∈ A}.  After removing (<ref>), (<ref>), (<ref>),  (<ref>) and (<ref>) from model (<ref>), we set _ij← T, ∀(i,j)∈ A and b_i← T, ∀ i∈ V, 





and solve the reduced model using Gurobi 8.1.1.
On average, the Lagrangian dual bounds are 3.9% tighter than the initial Lagrangian relaxation solutions (for the small­- and medium­-sized instances), with instance s-4 seeing the largest improvement of 12.4% (Table <ref>). The initial bounds yield a gap of 4.2% from the best solution obtained by Gurobi, whereas the Lagrangian dual bounds are very tight on our test instances – less than 0.1% of Gurobi's best solution. These bounds, moreover, are obtained by applying our proposed specialized algorithm before Gurobi finishes. We can make marginal improvements to those bounds considering the time required to obtain the Lagrangian dual bound, especially as instance size grows. For large­-scale instances, we do not improve upon the initial Lagrangian relaxation bound within the time limit (an average gap of 1.6% from the best solution obtained by Gurobi) as not many iterations have finished processing due to the complexity for solving a single iteration. However, the proposed algorithm will eventually achieve a good bound (<0.1% difference) in the absence of a time limit (Table <ref>).










 
 





 







































 


  §.§.§ Special Case II

Total operational time is defined as T := |A|×max{𝐝(i,j) | (i,j)∈ A}× 0.1 for all instances, except for s-4 where T := |A|×max{𝐝(i,j) | (i,j)∈ A}× 0.2 as the former deadline causes infeasibility.  We solve the same reduced model outlined in Section <ref> using Gurobi 8.1.1 with the addition of the following constraint:

    ∑_(i,j)∈ A t_ij≤ T.



On average, the Lagrangian dual bounds are 26.0% tighter than the initial Lagrangian relaxation solutions, with instance s-4 again seeing the largest improvement of 54.1% (Table <ref>). The initial bounds yield a gap of 40.8% from the best solution obtained by Gurobi, whereas the Lagrangian dual bounds are within an average of 0.2% of the best Gurobi solution. These bounds, moreover, are obtained by applying our proposed specialized algorithm before Gurobi finishes (except for instance m-3). We do not include results for the large­-sized instances as the specialized algorithm does not complete a single iteration within the 7200 second time limit. This is expected, though, as the algorithm crucially depends on the state space, which increases rapidly as each label extension generates two distinct labels.

 


§ CONCLUSIONS AND FUTURE RESEARCH
 


We have proposed and studied the multi­-vehicle covering tour problem with time windows (MCTPTW), which concerns unmanned vehicle routing with energy, time window, and coverage constraints. Motivated by applications for homeland security and public safety, our model aims to find an optimal route for each vehicle in a fleet through a set of secured waypoints in order to surveil a set of targets. We have formulated the problem as a mixed­-integer second­-order cone program with both coverage and risk objectives and test its performance using several benchmark instances. We have considered special cases where the fleet's energy supply and the time windows associated with each waypoint can be approximated by an operational deadline, and where vehicles are not subject to risk, in order to to develop a label correcting algorithm with an innovative set of dominance rules to solve the resulting Lagrangian relaxation problem given a set of dual multipliers. 

Computational experiments examining the efficacy and efficiency of the labeling algorithm show our methodology can produce high quality bounds in the same or less time than a commercial solver. 

For future study, we plan on improving the effectiveness of the labeling algorithm for large­-scale instances. Additionally, we aim to extend the model to manage non-stationary targets or targets whose locations are probabilistic. Other research extensions include the presence of time­-sensitive (on­-demand) targets. Similarly, introducing stochastic time windows or the potential loss of vehicles/waypoints would mimic various levels of uncertainty faced by surveillance operations. These extensions could better assist mission planners to adapt to changing environments within friendly and/or enemy territory.





