A Hybrid Approach to Persistent Coverage in Stochastic Environments
AA Hybrid Approach to Persistent Coverage in StochasticEnvironments (cid:63)
William Bentz a , Dimitra Panagou a a Department of Aerospace Engineering, University of Michigan, 1320 Beal Ave, Ann Arbor, MI, 48109, USA
Abstract
This paper considers the persistent coverage of a 2-D manifold that has been embedded in 3-D space. The manifold is subjectto continual impact by intruders which travel at constant velocities along arbitrarily oriented straight-line trajectories. Thetrajectories of intruders are estimated online with an extended Kalman filter and their predicted impact points contributenormally distributed decay terms to the coverage level. A formal hybrid control strategy is presented that allows for power-constrained 3-D free-flyer agents to persistently monitor the domain, track and intercept intruders, and periodically deployfrom and return to a single charging station on the manifold. Guarantees on intruder interception with respect to agent powerlifespans are formally proven. The efficacy of the algorithm is demonstrated through simulation.
Key words:
Coverage Control; Collision avoidance; Autonomous mobile robots; Cooperative control; Multi-agent systems .
The advent of inexpensive autonomous research plat-forms has spurred recent interest in teams of mobilesensors collaborating on complex surveillance and mon-itoring tasks. Coverage control problems have been par-ticularly popular due to their numerous applications:e.g., environmental monitoring (Smith et al. 2011), bat-tlefield surveillance (Bokareva et al. 2006), lawn mow-ing and vacuuming, search and rescue (Murphy et al.2008), and hull inspections (Choset & Kortenkamp 1999,Hollinger et al. 2013). The latter application is activelysupported by NASA whose work on the Mini AER-Cam paves the way for a future of extravehicular robotic(EVR) free flyers performing independent visual inspec-tions of spacecraft exterior areas of interest (Fredrickson (cid:63)
The authors would like to acknowledge the support ofthe Automotive Research Center (ARC) in accordancewith Cooperative Agreement W56HZV-14-2-0001 U.S. ArmyTARDEC in Warren, MI and the support by an Early Ca-reer Faculty grant from NASAs Space Technology ResearchGrants Program. The material in this paper was partiallypresented at the 56th IEEE Conference on Decision and Con-trol, December 12-15, 2017, Melbourne, Australia and the2018 American Control Conference, June 27-29, 2018, Mil-waukee, Wisconsin. See Section 1.2 for a comparison betweenthe present work and the conference versions.
Email addresses: [email protected] (William Bentz), [email protected] (Dimitra Panagou). et al. 2004). Free flyer visual inspection is the primarymotivating example for our work.Coverage is often partitioned into three classes of prob-lems: static, dynamic, and persistent. Static coverageproblems (e.g., area coverage, k-coverage and point cov-erage) often explore the optimal arrangement of sensornodes in a network and the agents tend to immobilizeafter this arrangement has been achieved (Cortes et al.2004). Dynamic coverage problems involve the activeexploration of a domain. Agents typically must sweeptheir sensors over all points of a domain until some de-sired level of coverage has been achieved (Hussein & Sti-panovi´c 2007, Liu et al. 2013, Stipanovi´c et al. 2013).Persistent coverage is often similar to dynamic coveragewith the addition of information decay within the envi-ronment: i.e., agents are required to continually returnto areas of interest in order to restore a deterioratingcoverage level.The term ”persistent coverage” appears as early asHokayem et al. (2007) where agents must cover allpoints in a 2-D convex polygonal domain every T (cid:63) timeunits. This was accomplished with the design of concen-tric polygonal trajectories with agents following closedpaths in steady state. The work in Song et al. (2013)is similar but also introduces a linear coverage decayrate for specific points of interest. In this paper, aswell as Smith et al. (2012), controller design is akin toregulating the velocity along paths generated offline toincrease observation time at select points of interest. As Preprint submitted to Automatica 7 August 2019 a r X i v : . [ c s . M A ] A ug he decay rates are known and time invariant, optimalspeed control is computed via linear programming.Palacios-Gas´os et. al have published multiple worksrecently on persistent coverage (Palacios-Gas´os et al.2016 a , b , 2017) which build specifically upon Smithet al. (2012). While the earlier work assumed both theexistence and knowledge of an optimal path to coverall points of interest, Palacios-Gas´os et al. (2016 b ) usestechniques from discrete optimization and linear pro-gramming to iteratively compute this path. The effectis that if the coverage decay rate of a specific point ofinterest is found to be insufficient to justify the transittime required to service it, then the point may be re-moved online from the path of the robot. Prior works,i.e., Smith et al. (2012), Song et al. (2013), would haveinstead driven the robot to quickly pass through thepoint. Similar techniques are used in Mitchell et al.(2015) which also considers that agents must periodi-cally return to refueling depots.In H¨ubel et al. (2008) and Song et al. (2011), the de-sired coverage level of the domain is maintained withdensity maps that yield additional observation time atselect areas of interest. In Song et al. (2011) the mapsare time-invariant while H¨ubel et al. (2008) considerstime-varying density maps that may be designed aroundmoving points of interest (e.g., aerial surveillance tar-gets). However, the latter work only uses density mapsin the derivation of control laws and not in the differen-tial equations governing coverage level evolution.Common themes through all of these persistent cover-age works are convex 2-D domains, predictable environ-ments, and simplified sensing and dynamic models foragents. Coverage surfaces embedded in R are consid-ered in Cheng et al. (2008); however, this work is closerto that of Hokayem et al. (2007) in that agents also followpreplanned trajectories without considering spatially-dependent coverage decay maps.Monitoring of stochastic environments is presented inPasqualetti et al. (2014), Yu et al. (2015), outside ofthe strict persistent coverage formulation. In Yu et al.(2015), the authors consider that agents must observeevents at multiple points of interest and the precise ar-rival times of events are unknown a priori . Arrival timestatistics are used to inform a multi-objective schedul-ing protocol that results in fixed cyclic servicing policies.In Pasqualetti et al. (2014), the environment containssmart intruders which actively attempt to evade a cam-era surveillance network. Camera motion is restricted toa single pan axis and thus the the system model is es-sentially that of a 1-D pursuit evasion problem. In this paper, we present (i) a formal hybrid controlstrategy for multi-agent persistent coverage of non-planar convex surfaces embedded in R that does notmake overly simplifying assumptions with respect to agent dynamic and sensing models, (ii) guarantees onagent interception of stochastic intruders, and (iii) anenergy-aware agent deployment and scheduling proto-col. To the best of our knowledge, we are the first topresent a formal hybrid approach to persistent coverage.Although many of the cited works use density functionsto encode points of interest, the difference in our ap-proach is both subtle and powerful. The density functionin H¨ubel et al. (2008) evolves subject to the motion ofan intruder and informs the control laws; however, it hasno effect upon the dynamics of the coverage level. Thus,only the intruder’s current location has any influence onthe motion of the agents, and the time-history of theintruder’s trajectory is forgotten. This necessitates thatagents must travel faster than targets in order to coverpoints associated with peaks in the density function be-fore they vanish. Works that do include a density func-tion in the coverage level evolution, e.g Palacios-Gas´oset al. (2016 a , b , 2017), tend to have fixed decay rates thatcannot respond or adapt to a changing environment. Incontrast, our algorithm utilizes a time-varying densityfunction, which is estimated online via extended KalmanFilter, to directly encode coverage decay over the surfacearound the predicted impact points of intruders. Thisencodes a memory effect which drives some agents tofollow coverage gradients towards areas that have previ-ously been or will soon be impacted by intruders.Our agents operate with finite resources and are requiredto periodically return to a refueling station while ob-serving stochastic events at locations and times that arenot known a priori . This approach is different from re-lated works, such as Yu et al. (2015) and Mitchell et al.(2015), where the locations of events are fixed and theauthors are concerned with optimal servicing routes be-tween these known stations.This hybrid system is a successor to our previous worksin Bentz & Panagou (2017, 2018). In Bentz & Panagou(2017), we derived the first of our hybrid modes (i.e.,local coverage mode) and our intruder state estimator.However, the agents had no power constraints and theapproach was unable to provide any formal guaranteeson intruder interception without additional operatingmodes. In Bentz & Panagou (2018), we derived theseadditional modes to present a hybrid approach to per-sistent coverage. Agents were now scheduled to inter-cept intruders and followed path-length optimal trajec-tories. However, formal guarantees on intruder intercep-tion were still limited to cases in which no collision-avoidance deadlocks had occurred. Furthermore, agentsdid not make effective use of local coverage mode as theywould often travel to the predicted impact points of in-truders and then remain stationary until the moment ofimpact thus contributing to a rising coverage error.This work extends the interception guarantee of Bentz& Panagou (2018) to an arbitrary number of collisionavoidance maneuvers and presents an entirely newmethod of collision avoidance over the prior works. It2lso reformulates numerous guard conditions within theautomaton to allow agents to explore actively aroundthe predicted impact points of intruders. Furthermore,this contribution revises our sensing function defini-tion utilized in the prior works which suffered from asingularity at the sensing cone vertex.This paper is organized as follows: Section 2 describesthe agents sensing and kinematic models and providesan overview of our hybrid control strategy, Section 2.2presents the trajectory estimator for particle intrudersand defines our coverage decay rate map, Sections 3-6describe each hybrid mode in detail, Section 7 verifiesthe algorithm in simulation, Section 8 summarizes ourcontributions and Section 9 is an appendix containingthe formal definition of our hybrid automaton. Consider a network of spherical autonomous agents in-dexed i ∈ { , ..., N } , of radius r i , whose motion is sub-ject to 3-D rigid body kinematics (Beard 2008): ˙ x i ˙ y i ˙ z i = cos Θ i cos Ψ i sin Φ i sin Θ i cos Ψ i − cos Φ i sin Ψ i cos Θ i sin Ψ i sin Φ i sin Θ i sin Ψ i + cos Φ i cos Ψ i − sin Θ i sin Φ i cos Θ i cos Φ i sin Θ i cos Ψ i + sin Φ i sin Ψ i cos Φ i sin Θ i sin Ψ i − sin Φ i cos Ψ i cos Φ i cos Θ i u i v i w i , (1) ˙Φ i ˙Θ i ˙Ψ i = i tan Θ i cos Φ i tan Θ i i − sin Φ i i sec Θ i cos Φ i sec Θ i q i r i s i , (2)where p i = [ x i y i z i ] T is the position vector and Ω i =[Φ i Θ i Ψ i ] T is the vector of 3-2-1 Euler angles takenwith respect to a global Cartesian coordinate frame G with origin O . The linear velocities [ u i v i w i ] T and an-gular velocities [ q i r i s i ] T are both presented in thebody fixed frame B i with origin p i . The state vector ofagent i is defined as ˜ q i = [ p Ti Ω Ti ] T . In the sequel, the ro-tation matrices of (1) and (2) shall be denoted R and R respectively. The agents travel within a stationary do-main, D ⊂ R . Their task is to survey a two-dimensionalmanifold, C ⊂ D , known as our surface of interest. Forthe purpose of this work we assume that the surface isan ellipsoid of revolution; however, it should be notedthat the coverage laws, as well as the collision avoidancestrategy, can be easily adapted for any convex surface.The ellipsoid has semi-major axis x C ,r and semi-minoraxis z C ,r aligned with the global coordinate axes ˆ x G andˆ z G respectively with center at O . The circumflex (i.e.,hat) symbols denote unit vectors.Each agent, i , is equipped with a forward facing sensorwhose spherical sector footprint shall be referred to as S i . This model, though intended to be generic, is sim-ilar to conical camera models presented in other workson dynamic coverage (see Xie & Zhang (2013)). Ourmodel differs in terms of its heterogeneity, i.e. S i providesanisotropic sensing data that degrade in quality towardsthe periphery of the footprint and changes with respectto distance from the sensor. Degradation over distanceis not monotonically decreasing but instead contains apeak located near the vertex of S i as in Hexsel et al.(2013). This is motivated by the fact that the probabil-ity of event detection by a camera decreases when ei-ther very far from or very close to the lens. Anisotropicsensing is encoded through the definition of the sensingconstraint functions for each agent i : c i = β i R − (˜ x − x i ) − (˜ y − y i ) − (˜ z − z i ) , (3a) c i = α i − φ i , (3b)for β i = min { , η i (cid:16) (˜ x − x i ) + (˜ y − y i ) − (˜ z − z i ) (cid:17) } with real constant η i >> R is the sensing range,˜ p i = [˜ x ˜ y ˜ z ] T is the position of a point within S i withrespect to G , α i is the angle between the periphery andcenterline of the spherical sector (the ˆ x B i axis), and φ i is the angle between r ˜ p i /p i = ˜ p i − p i (resolved in G by construction) and the ˆ x B i axis given as the inversecosine of the dot product of ˆ r ˜ p i /p i and ˆ x B i resolvedin G : φ i = arccos (cid:0) ˆ r ˜ p i /p i · ˆ x B i | G (cid:1) . Note that: ˆ r ˜ p i /p i = √ (˜ x − x i ) +(˜ y − y i ) +(˜ z − z i ) [(˜ x − x i ) (˜ y − y i ) (˜ z − z i )] T ,and ˆ x B i | G is determined by multiplying R by [1 0 0] T :ˆ x B i | G = (cid:104) cos Ψ i cos Θ i sin Ψ i cos Θ i − sin Θ i (cid:105) T . Agent i is thus capable of detecting objects that lie withinan angle of 2 α i > x B i axis and a range of R >
0. The model for agent i is depicted in Fig. 1. Fig. 1. Agent i is modeled as a sphere of radius r i and has aforward facing sensor footprint, S i . Sensing constraint func-tions c ki , ∀ k ∈ { , } , encode a decay in sensing quality alongthe depth and towards the periphery of S i . Let us denote max { , c ki } = C ki . One can define thesensing function that represents the quality of informa-tion available at each point over the sensing domain as:3 i (˜ q i , ˜ p ) = (cid:40) C i C i C i + C i , if card (cid:0) ¯ C i (cid:1) < ∧ r ˜ p i /p i > , otherwise, (4)where ¯ C i is the set of zero elements in C ki . S i (˜ q i , ˜ p ) takesa value of zero outside of S i . Note that S i (˜ q i , ˜ p ) is de-fined over all of D and thus has static bounds. S i (˜ q i , ˜ p )is continuous in ˜ p while taking a value of zero along ∂ S i .In verifying this continuity, it is important to note that S i (˜ q i , ˜ p ) approaches zero from within S i in the limit thateither card (cid:0) ¯ C i (cid:1) = 2 or r ˜ p i /p i = 0 are satisfied. The for-mer condition may be verified by taking a limit of thefirst piecewise definition of (4) as C i and C i tend tozero. The latter condition results from our definitions of β i and η i which encode that S i (˜ q i , ˜ p ) drops off rapidlywhen in very close proximity to the vertex of S i . Increas-ing η i >> p i . Define the coverage level providedby agent i at time t as: Q i ( t, ˜ p ) = (cid:90) t S i (˜ q i ( τ ) , ˜ p ) C (˜ p ) dτ, (5)where C is defined as: C (˜ p ) = (cid:26) , ∀ ˜ p ∈ C ;0 , ∀ ˜ p / ∈ C , and en-codes that the accumulation of sensing information onlyoccurs along our surface of interest, C .As the agents cover C , a set of N p high-speed particleintruders denoted k ∈ { , ..., N p } , each of which trav-els in an arbitrary direction at constant velocity, passthrough the domain. The particles are assumed to beuncontrolled and cannot deviate from their initial tra-jectories. No assumptions are made with respect to thesource of the particles or whether they are intelligentlygenerated. Each particle shall have an associated mapdecay term, Λ k ( τ, ˜ p ), which is defined later in Section2.3. We may now define the global coverage level: Q ( t, ˜ p ) = N (cid:88) i =1 Q i ( t, ˜ p ) − N p (cid:88) k =1 (cid:90) t Λ k ( τ, ˜ p ) C (˜ p ) dτ. (6)In this work, coverage refers to the accumulation of sens-ing data over time. Points, ˜ p , are said to be sufficientlycovered when Q ( t, ˜ p ) ≥ C (cid:63) . The goal is to derive a hy-brid control strategy which persistently sweeps S i across C while emphasizing surveillance within some bound ofthe predicted impact points of particles k ∈ { , ..., N p } on C . More specifically, we establish theoretical guaran-tees on the worst case path length from any agent toany arbitrary impact point thus guaranteeing intercep-tion for prescribed bounds on intruder speed, detectionrange, and agent velocity. This must be done while avoid-ing collisions. Let us define collision and interception. Definition 1
Agent i is said to have intercepted particle k if i is within a ε bound of the estimated impact pointof k for a finite interval of time leading up to the impact. Agent i shall spend this duration of time sweeping thearea in local coverage thus gathering information.Note that intruders are unaffected by agents and shallalways impact the surface and then disappear. This doesnot damage the agent which is free to resume other tasksupon conclusion of interception. Definition 2
Agent i avoids collision so long as (cid:107) p i ( t ) − p j ( t ) (cid:107) > r i + r j , ∀ t ≥ , ∀ j (cid:54) = i ∈ { , ..., N } and (cid:107) n i (cid:107) > r i where the vector n i has direction normal to C andlength equal to the Euclidean distance of its intersectionpoint on C to p i .Agents operate with finite power resources and arerequired to return every T (cid:63) time units to a fueling sta-tion denoted F . Thus, a scheduling protocol is derivedwhereby agents periodically deploy from F to coverwithin assigned partitions of C . These partitions arebounded by latitude lines and are sorted by geodesicdistance from F with agents deploying to the partitionfurthest from F and then transferring between adjacentpartitions every T (cid:63) N time units as their power resourcedwindles requiring a return to F within T (cid:63) time unitsafter deployment. This partitioning scheme also has thebenefit of ensuring that the network of agents is welldistributed across C with agents nominally assigned tointercept intruders with predicted impact points withintheir own partition.Agent i is capable of localizing itself in G and detectingwhether there exists j such that (cid:107) p i ( t ) − p j ( t ) (cid:107) ≤ R .Furthermore, agents i and j can communicate their de-ployment times to one another. A centralized network isrequired to publish the current coverage level Q ( t, ˜ p ) to all agents and to estimate the trajectories of intrudersusing an omnidirectional range sensor whose measure-ments are fed through an extended Kalman filter. Com-putation of Q ( t, ˜ p ) is contingent upon continuous trans-mission of agent state ˜ q i to the centralized network. Thecentralized network assigns each intruder to an unas-signed agent at closest latitude to the predicted impactpoint. It also transmits detection time as well as esti-mated location and time of impact to the agent. We assume that the omnidirectional range sensor (e.g.,LiDAR) is co-located with O and provides measure-ments of each particle’s position in spherical coordinates.We also assume that particle detection and state esti-mate initialization occur while the distance of the parti-cle from O is greater than or equal to R det + x C ,r where R det is a lower bound on distance from detection to im-pact. We define the model for the motion of particle k : In practice, it is not necessary to publish Q ( t, ˜ p ) , ∀ ˜ p ∈ D toevery agent. For agent i to compute its local coverage controlsignal, it need only values for Q ( t, ˜ p ) within a closed ball ofradius R due to the fact that S i (˜ q i , ˜ p ) = 0 , ∀ ˜ p / ∈ ¯ B R p i ( t ).This substantially reduces the communication overhead. q k ( t ) = (cid:34) × I × × × (cid:35) ˜ q k ( t ) , (7)˜ z k ( t ) = (cid:112) x k + y k + z k atan2 ( y k , x k )arccos (cid:18) z k √ x k + y k + z k (cid:19) + (cid:15), (8)where ˜ q k = [ x k , y k , z k , ˙ x k , ˙ y k , ˙ z k ] T and ˜ z k = [ ρ k , θ k , ψ k ] T are the Cartesian state and spherical coordinate mea-surement vectors of particle k resolved in G . We as-sume that particle speed is upper bounded such that (cid:112) ˙ x k + ˙ y k + ˙ z k ≤ U intmax . ρ k , θ k , and ψ k are the range,azimuthal angle, and polar angle of k respectively. In thesequel, the matrix in (8) shall be denoted ˜ h ( x k , y k , z k ).Assume that the measurement noise, (cid:15) , is zero-meanGaussian and has covariance R = diag (cid:16) σ ρ , σ θ , σ ψ (cid:17) .This system models high-speed particles incident upona surface with negligible drag (e.g., micrometeoroidsimpacting a spacecraft hull); thus, it is reasonable toomit the process noise. The state and covariance esti-mates, ˆ˜ q and P k , are computed with a continuous-timeextended Kalman filter which is initialized upon particle k ’s detection at time t dk . At any time t , we define our decay rate map for parti-cle k in terms of its predicted position and covarianceevolution over a horizon T H,k ( t ). As the particles are as-sumed to travel at fixed velocities , the predicted valuesfor Cartesian position ˜ p (cid:48) k ( t + τ ) and associated covari-ance P k ( t + τ ) are defined as ˜ p (cid:48) k ( t + τ ) = G ( τ ) ˆ˜ q k ( t ),and P (cid:48) k ( t + τ ) = G ( τ ) P k ( t ) G ( τ ) T respectively where G ( τ ) = [ I × τ I × ] and ˆ˜ q k ( t ) is our current estimatefor ˜ q k ( t ). We define the decay rate map associated withparticle k as the integral of our predicted normal distri-bution N (˜ p (cid:48) k ( t + τ ) , P (cid:48) k ( t + τ )) through horizon T H,k : Λ k ( t, ˜ p ) = (cid:90) T H,k ( t )0 λ k N (cid:0) ˜ p (cid:48) k ( t + τ ) , P (cid:48) k ( t + τ ) (cid:1) dτ, (9) where λ k > λ k < S i . For t < t dk , define Λ k ( t, ˜ p ) =0 , ∀ ˜ p ∈ D . Our formulation for (9) essentially takes anormal distribution for the position of particle k at time t and cumulatively propagates it forward in time up toour horizon T H,k ( t ). The horizon is lower-bounded by anestimate of the remaining time until impact of particle k on C . This may be computed using ˜ q k ( t ) along with The guarantee of intruder interception presented in thiswork can be extended to intruders with time-varying veloc-ities that are bounded by U intmax . However, it is still requiredthat intruders follow straight line trajectories such that thenetwork may estimate fixed impact points. the surface geometry. With this design, Q ( t, ˜ p ) decaysalong the predicted trajectory of k with tapering omni-directional decay rates spreading out from the predictedpath. This design lends itself naturally to our local cov-erage formulation, which is gradient following in nature,in that the agents may follow these tapered decay pathstowards the predicted impact points on our surface ofinterest. Fig. 2. Agent i operates in accordance with this automaton.For clarity, elements of the reset map and brief descriptionsof transitions are included. The coverage strategy for agent i is represented by thehybrid automaton in Fig. 2. Rigorous definitions of allentities of the automaton, including the guard conditionsand reset maps, are included in the appendix. Note thateach agent operates in accordance with its own automa-ton and thus an arbitrary number of agents may be inany operating mode at any given time. Before proceed-ing, we provide a brief overview of each mode. • Local Coverage: This mode governs the active ex-ploration of our surface of interest C . When active,the agent continuously seeks to orient and trans-late S i such that S i intersects portions of C with5 lower coverage level. This is conceptually similarto following the gradient of the coverage error. Anagent currently assigned to an intruder may operatein Local Coverage while within an ε bound of thethe intruder’s predicted impact point. Any agentnot currently assigned to an intruder shall operatein Local Coverage assuming that it is within its as-signed latitude partition. Operation in Local Cov-erage is always concurrent with agent assignment tothe lowest concentric surface (see Surface Transferbelow) and transition to Local Coverage can occurfrom any mode aside from Return to Base Mode. • Particle Intercept: In this mode, an intruder is as-signed to agent i and i is guided along its assignedsurface to the predicted impact point of the in-truder. After intruder assignment occurs, the agentwill nominally remain in Particle Intercept Modeuntil after the intruder impacts C ; however, theagent may temporarily leave the mode before im-pact to avoid collision through Surface TransferMode or to explore in Local Coverage within a ε bound of the the intruder’s predicted impact point. • Partition Transfer: This mode is defined for agentsthat are not currently assigned to an intruder andits purpose is to ensure that the agents are spa-tially distributed across the entire surface area of C .Activation of this mode will guide agent i along alongitudinally-oriented geodesic trajectory until itsposition satisfies a set of latitude constraints, i.e.,agents travel to the southernmost latitude partitionupon deployment and transition through progres-sively northern partitions as their fuel is depleted.Transition to this mode can occur from any othermode. The partitioning scheme is shown in Fig. 6. • Surface Transfer: This mode’s primary purpose is toensure that agents avoid collision with two distinctcases resulting in its activation. In the first case,two or more agents have violated a safe-proximitycondition. The mode removes select agents from thedeadlock by guiding them along vectors normal to C to a higher-altitude ellipsoidal surface concentricwith C . An agent trajectory is then temporarily con-fined to this newly assigned surface until it reachesthe surface projection of its destination. The sec-ond case occurs under the condition that the agenthas arrived at the projection of its destination ona higher-altitude surface. The mode is activated toreturn the agent to the innermost surface. Transi-tion to this mode can occur from any other modeaside from Return to Base as agents in the lattermode always take priority in a deadlock. The sur-face transfer geometry is illustrated in Fig. 4. • Return to Base: The final mode is activated whenthe time since an agent’s deployment has surpassedsome threshold. It guides the power-critical agentalong the optimal trajectory to the refueling sta-tion. After charging, the agent is redeployed. Agentdeployments occur one at a time with a fixed period.A power critical agent in Particle Intercept Mode or Surface Transfer Mode shall first complete its taskof intercepting the assigned intruder or transfer-ring surfaces before transitioning to Return to BaseMode. If an agent is designated as power-criticalwhile in Partition Transfer Mode it shall immedi-ately abandon its task and transition to Local Cov-erage which shall result in instantaneous transitionto Return to Base Mode. Theoretical guaranteeson successful return to base with respect to agentpower lifespan in accordance with our automatonis presented in Theorem 2 of Section 5.
Local coverage constitutes the first of five hybrid modesin our automaton. This mode is gradient following innature and commands agent i to always seek to orientand translate S i such that the volume of uncovered spaceintersecting S i is increased. In this way, it emphasizesactive exploration of the domain by agents that are notcurrently assigned to either monitor intruders or relocatewithin the domain. The control laws are designed suchthat agent motion in local coverage shall tend to reducethe rate of growth of the global coverage error. Definethe global coverage error with respect to C (cid:63) as: E ( t ) = (cid:90) D h ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) d ˜ p, (10)where h ( w ) = (max { , w } ) with first derivative h (cid:48) = dhdw = 3(max { , w } ) and second derivative h (cid:48)(cid:48) = d hdw =6(max { , w } ). Our local coverage control laws are de-rived via differentiation of (10). This is included in theAppendix in the interest of space. The result is the se-lection of the following control strategy: u loci = k u a i ( t, Q ( t, ˜ p )) (cid:112) a i + a i + a i + ˆ x B i · ρ l,i , (11a) v loci = k v a i ( t, Q ( t, ˜ p )) (cid:112) a i + a i + a i + ˆ y B i · ρ l,i , (11b) w loci = k w a i ( t, Q ( t, ˜ p )) (cid:112) a i + a i + a i + ˆ z B i · ρ l,i , (11c) r loci = ¯ r i sat (cid:0) k r a i ( t, Q ( t, ˜ p ))¯ r i (cid:1) + ˆ y B i · ρ a,i , (11d) s loci = ¯ s i sat (cid:0) k s a i ( t, Q ( t, ˜ p ))¯ s i (cid:1) + ˆ z B i · ρ a,i , (11e) where: ρ l,i = − ln (cid:18) γR − r i ( (cid:107) n i (cid:107) − r i ) (cid:19) R − ˆ n i , (12) ρ a,i = ξ R − n i · ˆ z G ) − Θ i atan2 ( − ˆ n i · ˆ y G , − ˆ n i · ˆ x G ) − Ψ i . (13) l,i is a collision avoidance term with respect to the sur-face of interest. It takes a value of zero when agent i ’snormalized distance from C is γR for γ ∈ (0 ,
1] and islogarithmically repulsive and attractive from the surfacewhen the distance is decreased or increased respectively.Small values for γ tend to direct the agent to travel closerto the surface. This coincides with a smaller cross sectionof S i intersecting the surface but is also typically associ-ated with a higher quality of sensing. A larger choice for γ will direct the agent to fly at a higher altitude with re-spect to the surface and thus the area covered by S i willtend to be broader with a decreased quality of sensing. ρ a,i , for ξ <<
1, encodes that the agents should tendto align ˆ x B i with − ˆ n i if the coverage terms associatedwith r i and s i have become sufficiently small. The phys-ical meaning of ρ a,i is to direct S i back onto C if it hasreached a configuration in which it no longer intersects C . See Fig. 3 for further details. Fig. 3. As agent i explores C , ρ l,i is parallel to n i for (cid:107) n i (cid:107) < γR , antiparallel to n i for (cid:107) n i (cid:107) > γR , and the zerovector otherwise. This term prevents collision of i with C andprevents i from flying away from C . ρ a,i tends to direct S i onto C . ¯ r i and ¯ s i are saturation limits for the coverage angularvelocity inputs to the system. k u , k v and k w are tun-ing gains which are chosen to satisfy (cid:112) k u + k v + k w ≤ U agtmax . As ρ (cid:96),i is normal to the surface, it can be shownthat U agtmax is an upper bound to agent velocity tangen-tial to C . Assuming that particle k is embedded within the surfaceupon impact, its position shall intersect C at most onetime. We define particle k ’s estimated impact time as t ck = min (cid:18) t ∈ R + | ( ˆ x k + ˙ˆ x k t ) x C ,r + ( ˆ y k + ˙ˆ y k t ) x C ,r + ( ˆ z k + ˙ˆ z k t ) z C ,r = 1 (cid:19) ,with estimated impact point ˜ p (cid:48) k ( t ck ) = G ( t ck − t ) ˆ˜ q k ( t ).Upon detection, particle k is assigned to a free agent i with the minimum distance from the estimatedpoint of impact along the ˆ z G direction. We definea new index, i k , as the index of the agent assignedto intercept particle k at destination p id = ˜ p (cid:48) k ( t ck ): i k = argmin i ∈{ ,...,N }| i p (cid:54) =1 ,f i (cid:54) =1 (cid:107) ¯ z (cid:48) k ( t ck ) − z i ( t dk ) (cid:107) . Note that ¯ z (cid:48) k ( t ck ) is the z component of ¯ p (cid:48) k ( t ck ) and f i ∈ { , } is a particle assignment flag for agent i defined as 0when the agent is free (i.e., not currently assigned aparticle). i p ∈ { , ..., N } , the power index of agent i ,shall be fully described in Section 5; however, it shouldbe noted that the definition of i k implies that there areat most N − γR normal to C in thenominal case that they are not maneuvering to avoidcollision. We define an ellipsoid of revolution, C , whichis concentric with C and has the property that each semi-principal axis is γR longer than its associated counter-part in C , i.e., x C ,r = x C ,r + γR , and z C ,r = z C ,r + γR .The nominal trajectories of i are attractive to C .Agents maneuvering to avoid collision shall transferto additional concentric ellipsoidal surfaces each sepa-rated by a distance of R . These surfaces are denoted C , C , ..., C N − with associated semi-principal axes x C ,r = x C ,r + ( γ + 1) R and z C ,r = z C ,r + ( γ + 1) R , x C ,r = x C ,r + ( γ + 2) R and z C ,r = z C ,r + ( γ + 2) R ,etc. Surface assignment and transfer scheduling in colli-sion avoidance mode is described in full detail in Section6 and the geometry is illustrated in Fig. 4. Fig. 4. Three agents enter a deadlock in (a). The green agent,which has the greatest time since deployment, is prioritizedto continue on C and the red and blue agents are each trans-ferred to C before entering a second deadlock in (b). The blueagent, which has the second greatest time since deployment,is prioritized to continue on C and the red agent is trans-ferred to C before continuing along geodesic to proj C ˜p (cid:48) k ( t ck )in (c). Red agent transfers back to C directly above predictedimpact point of particle k in (d). Note that surface transfertrajectories are always normal to C µ i , ∀ µ i ∈ { , ..., N − } . When agent i has been assigned to intercept particle k , f i is set to 1 and it is said to have transitioned into par-ticle intercept mode. In this mode, agent i shall nomi-nally follow the optimal trajectory along C to within a ε bound of the projection of point ˜ p (cid:48) k ( t ck ) onto C (de-noted proj C ˜p (cid:48) k ( t ck )). The agent shall then transition to7ocal coverage to actively explore within this ε bounduntil t > t ck at which time f i is set to 0. If local coverageguides the agent out of the ε bound, particle interceptmode will again guide the agent back inside the bound.The optimal trajectory is referred to as a geodesic andits computation may be executed in an iterative man-ner. Specifically, we use Vincenty’s formulae as presentedin Vincenty (1975 a ). For cases involving nearly antipo-dal points in which the standard inverse method doesnot converge, we use Vincenty’s supplemental algorithmpresented in Vincenty (1975 b ).As an input, Vincenty’s algorithm requires an ellipsoid ofrevolution along with two points, current and desired po-sition, on that surface. The algorithm returns a headingangle measured clockwise from North. This heading an-gle shall be referred to as χ i . We now define the headingunit vector ˆ ν i which lies in a plane tangent to the surfaceat p i . It may be computed by rotating the North-pointingvector at p i clockwise by an angle of χ i within the tan-gent plane. For our implementation of Vincenty’s algo-rithm, we input the following: C µ i for surface assignmentindex µ i ∈ { , ..., N − } , p i , and proj C µ i ˜p (cid:48) k ( t ck ). The po-sition controller used to guide agent i to proj C µ i ˜p (cid:48) k ( t ck )is composed of two terms: one which commands veloc-ity tangential to C µ i along ˆ ν i and one logarithmic termwhich commands velocity normal to C µ i in order to con-strain the geodesic trajectory of i to C µ i . The particleintercept mode position control law is: u pimi v pimi w pimi = U agtmax R − ˆ ν i − ln (cid:16) γ + µ i ) R − r i ( (cid:107) n i (cid:107) − r i ) (cid:17) ˆ n i (cid:107) ˆ ν i − ln (cid:16) γ + µ i ) R − r i ( (cid:107) n i (cid:107) − r i ) (cid:17) ˆ n i (cid:107) . (14) As agent i travels along the geodesic, it is desirable thatit should point S i towards C . Therefore, the orientationcontroller for particle intercept mode is similar to thatof Section 3: q pimi r pimi s pimi = R − n i · ˆ z G ) − Θ i atan2 ( − ˆ n i · ˆ y G , − ˆ n i · ˆ x G ) − Ψ i , (15)which is essentially a proportional controller that tendsto align ˆ x B i with − ˆ n i . As (14) commands the vehi-cle to follow the optimal length path along C µ i to proj C µ i ˜p (cid:48) k ( t ck ), we can establish a few guarantees onsystem performance. To simplify notation, define: g C N − = ∞ (cid:88) n =1 (cid:18) (2 n − n n ! (cid:19) (cid:16) x C N − ,r − z C N − ,r x C N − ,r + z C N − ,r (cid:17) n (2 n − , (16)and g C is defined similarly in terms of the semi-principalaxes of C . Lemma 1
Let us assume that agent i has been assignedto particle k with f i := 1. Given an arbitrary agentposition p i ( t dk ) and an arbitrary predicted impact pointfor the intruder ˜ p (cid:48) k ( t ck ), there exists an upper bound tothe maximum path length until interception: P max ≤ πx C N − ,r + π (cid:0) x C N − ,r + z C N − ,r (cid:1) g C N − + 2 ( N − R . PROOF.
At the moment that f i := 1 we have thatagent i transitions to Particle Intercept Mode. Un-der the condition that the agent has not yet comewithin proximity of the predicted impact point, i.e., (cid:107) p i − proj C µ i ˜p (cid:48) k ( t ck ) (cid:107) > ε , we have that only G ( ζ i , ζ i )and G ( ζ i , ζ i ) are defined (see Appendix). These twotransitions occur sequentially for each deadlock eventthat agent i encounters as it approaches proj C µ i ˜p (cid:48) k ( t ck ).In any given deadlock arrangement, one agent remainson its current surface without ascending to a higher one.This implies that µ i = 1 for at most N − C . Furthermore, this implies that µ i = 2for at most N − µ i = N for at most zero agents. The worst cast surface assign-ment that can be incurred during sequential cycles of(( ζ i , ζ i ) , ( ζ i , i µ i = N − C µ i shall always be less than thegeodesic path length between the same two points pro-jected onto surface C µ i +1 , we may bound the geodesicportion of the trajectory by one that is constrained en-tirely to C N − . We denote this term P geo . As any twopoints on C µ i can be connected by a path of constantlatitude P lat followed by a path of constant longitude P long , we have that: P geo ≤ P lat + P long . (17)For two generic points on C N − , we have that: P lat ≤ πx C N − ,r , (18) P long ≤ π (cid:0) x C N − ,r + z C N − ,r (cid:1) g C N − . (19)where the bound on P lat is half of the circumferenceof the ellipsoid of revolution C N − about its equa-tor and the bound on P long is half of the perime-ter of the revolved ellipse. The infinite series expres-sion term, denoted g C N − in (19), is first presentedin Ivory (1798). The remaining portion of the pathlength is simply the straight line segments connecting C to C N − and back again. This length is precisely2 ( N − R . Thus, P max = P geo + 2 ( N − R as illus-trated in Fig. 5. Invoking (18) and (19) gives us P max ≤ πx C N − ,r + π (cid:0) x C N − ,r + z C N − ,r (cid:1) g C N − + 2 ( N − R .This concludes the proof. Theorem 1
Assuming that the bounds on intruder ve-locity and range from detection to impact satisfy R det U intmax > ig. 5. The longest possible path from p i ( t dk ) to ˜ p (cid:48) k ( t ck ),taken by agent i assigned to intercept particle k , is illustratedabove. We denote this path as P max and it may be upperbounded as established in Lemma 1. P max U agtmax , where P max may be bounded via Lemma 1, agent i shall reach proj C ˜p (cid:48) k ( t ck ) before t ck . PROOF.
Given the fact that agents in the particle in-tercept and surface transfer modes travel at speed U agtmax ,we have that the time t ik required to travel from p i ( t dk )to proj C ˜p (cid:48) k ( t ck ) must satisfy: t ik ≤ P max U agtmax . Given that t ck − t dk ≥ R det U intmax , agent i reaching proj C ˜p (cid:48) k ( t ck ) be-fore t ck implies that R det U intmax > t ik . This is guaranteed if R det U intmax > P max U agtmax This concludes the proof.
Remark 1
At any given time, there are at most N − N − πx C N − ,r + π ( x C N − ,r + z C N − ,r ) g C N − +2( N − RU agtmax . As this is a persistent coverage protocol, which operatesindefinitely, it is necessary to establish an agent deploy-ment and scheduling protocol that realistically consid-ers the agents’ finite power and/or propulsive resources.Our strategy is to periodically deploy agents from a fu-eling station F which we assume to be located at theNorth pole of C , i.e., at the point [0 0 z C ,r ] T . Define T (cid:63) as the power lifespan of each agent in the network.Given T (cid:63) and N , we define our deployment and schedul-ing protocol such that one agent is deployed from F ev-ery T (cid:63) N seconds. The initial deployment is that of agent i = 1 at t = 0 seconds with agent i = 2 following at t = T (cid:63) N . This continues indefinitely with the second de-ployment of agent i = 1 occurring at t = T (cid:63) seconds.In order to adequately distribute agents across C , it isdesirable to partition the domain and assign agents tomonitor separate regions. Specifically, partitioning thedomain by latitude, rather than longitude, ensures thatagents are poised to intercept particles without the need for frequent crossings of the equator which tend to be as-sociated with larger values of P geo on an oblate spheroid.Define the power index of agent i as i p ( t ) = 1 +mod (cid:0) i − − (cid:4) tNT (cid:63) (cid:5) , N (cid:1) where the first argument ofour modulo operation is the dividend and the secondargument is the divisor. The lower-bracketed delimitersrepresent the floored division operation. Upon deploy-ment from F , agent i has power index i p = N and thisindex is reduced by one every T (cid:63) N seconds until i p = 1,i.e., agent i is the power critical agent. Note that no twoagents may share the same power index as a result ofour periodic deployment and scheduling protocol.Latitude partitions are characterized by a static upperbound in ˆ z G denoted ¯ z i p − and a static lower bound¯ z i p − . Rather than dynamically sizing partitions relativeto agent power resources, we divide partitions such that N − C toexplore. This choice maximizes the coverage of any indi-vidual partition as the allocation of a larger partition toa recently deployed agent would result in less effectivecoverage of that partition. Agents are sorted by their re-maining power and transfer between partitions that areprogressively closer to F as their power resource expires.Define the surface area of our ellipsoid of revolution C as: A C = 2 πx C ,r (cid:18) − z C ,r x C ,r (cid:19)(cid:18)(cid:114) − z C ,r x C ,r (cid:19) artanh (cid:32)(cid:115) − z C ,r x C ,r (cid:33) . (20) The agent with i p = 2 is assigned to monitor the par-tition characterized by upper bound at north pole of C ,i.e., ¯ z = z C ,r . The lower bound ¯ z may be computed bydividing (20) by ( N − z ,and then numerically solving for ¯ z : A C N − (cid:90) ¯ z z C ,r π (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) x C ,r − x C ,r ˜ z z C ,r z x C ,r x C ,r (cid:16) z C ,r − z C ,r ˜ z (cid:17) d ˜ z. (21) One may then iteratively solve for the remaining boundsfor increasing values of i p up to i p = N − A C N − (cid:90) ¯ zip − zip − π (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) x C ,r − x C ,r ˜ z z C ,r z x C ,r x C ,r (cid:16) z C ,r − z C ,r ˜ z (cid:17) d ˜ z. (22) The final computation of (22) for i p = N is not neces-sary as ¯ z N − is the south pole of C , i.e., ¯ z N − = − z C ,r ,although this may be shown through numerical com-putation as well. Our partitioning strategy for the casewhere N = 4 is presented in Fig. 6.Note that no partition has been assigned to the agentfor which i p = 1. This is the power critical agent andit shall have flag f i := 1 at the instant i p := 1. Thepower critical agent cannot be assigned a new particle tointercept after i p := 1 as this opens the possibility that9 ig. 6. Our domain partitioning scheme for C is illustratedabove. Agents with i p ∈ { , , } are indicated with blue,green and black S i respectively. Their partitions are sepa-rated by latitude lines upper bounded at ¯ z i p − and lowerbounded at ¯ z i p − . The power critical has S i indicated in red. particle assignment could occur near the end of the T (cid:63) N time window during which time the agent with i p = 1should be transitioning back to F to exchange its powersource. The power critical agent will instead spend themajority of this time window in local coverage modeassisting the other agents in gathering information. Itcan only be tasked with intercepting a particle if thisassignment had occurred previously when i p = 2. Inthis scenario, the agent should be capable of interceptingparticle k and then transitioning back to F so long asa bound is established on the length of our deploymentscheduling window T (cid:63) N . Theorem 2
If agent power lifespan T (cid:63) satisfies T (cid:63) N ≥ t ck − t dk + π U agtmax ( x C ,r + z C ,r ) g C , ∀ k then the agentwith i p = 1 shall always be capable of reaching F within T (cid:63) N of the time at which i p := 1. PROOF.
Consider the worst-case scenario in whichthe agent with i p = 2 is assigned to intercept par-ticle k at the instant before i p := 1. It’s remainingflight time is currently T (cid:63) N . The time required to in-tercept the particle is t ck − t dk , after which our con-trol strategy dictates that the agent will follow ageodesic trajectory to F . As F lies at the north poleof C , this will be a trajectory of constant longitudewhich may be upper bounded by a length half theperimeter of our revolved ellipsoid: π ( x C ,r + z C ,r ) g C by definition. As the agent is controlled by (14)with a North-pointing ˆ ν i , it will proceed along thisgeodesic at speed U agtmax . Thus the time required tocomplete this trajectory is π U agtmax ( x C ,r + z C ,r ) g C and we may bound our deployment window: T (cid:63) N ≥ t ck − t dk + π U agtmax ( x C ,r + z C ,r ) g C , ∀ k . This concludesthe proof. Remark 2
The appropriate design method for this surveillance system is to first ensure that the time fromdetection to impact of any arbitrary particle, t ck − t dk ,as governed by the omnidirectional range sensor sat-isfies Theorem 1. One must subsequently ensure thatpower lifespan T (cid:63) , for all agents, satisfies Theorem 2. If an agent with i p ∈ { , ..., N } lies outside of its pre-scribed partition, and we have i f = 0, then the agentshall enter partition transfer mode. This mode uses thesame geodesic position and orientation controllers (14)and (15) with the destination position set to the point: p id = [ x id y id z id ] T = x C ,r cos (cid:16) arcsin (cid:16) z id z C ,r (cid:17)(cid:17) cos (cid:18) atan2 (cid:0) y i ( t ) , x i ( t ) (cid:1)(cid:19) y C ,r cos (cid:16) arcsin (cid:16) z id z C ,r (cid:17)(cid:17) sin (cid:18) atan2 (cid:0) y i ( t ) , x i ( t ) (cid:1)(cid:19) ¯ z i p − , if z i < ¯ z i p − ; or ¯ z i p − , if z i > ¯ z i p − , i.e., the closest point along the agent’s current longitudewhich lies on the boundary of its assigned partition.The return to base mode is similar to partition transfermode but is defined for the agent with i p = 1. This modeis activated when the time since agent i ’s last deploymentfrom F , denoted t i F ≥ T (cid:63) − π U agtmax ( x C ,r + z C ,r ) g C asestablished in Theorem 2. The control strategy is thesame as partition transfer mode with the desired positionset to F . Control laws for partition transfer mode andreturn to base shall be denoted with superscripts ptm and rtb respectively. The primary purpose of surface transfer mode is to en-code collision avoidance and it can be transitioned intofrom any other mode aside from the return to base mode.This mode is triggered for agent i , assigned to surface C µ i , when we have the condition that (cid:107) p i − p j (cid:107) ≤ R for i (cid:54) = j . Denote ˜ j = i ∪ j as the set of agents satisfy-ing this condition. Agents in ˜ j are ranked by t ˜ j F . Oneagent, denoted i pr , whose value for t ˜ j F is highest, i.e., i pr = argmax ˜ j (cid:16) t ˜ j F (cid:17) is permitted to proceed. The re-maining agents increment their surface assignment in-dices, µ i , by one and transition to surface transfer mode.This mode controls the agents to follow ˆ n i until theyhave transferred to their newly assigned concentric sur-face at a height R above the previous. Note that in gen-eral, convexity of surface C is required to ensure thatintersections of n i and n j , ∀ i (cid:54) = j, lie within the interiorspace that is bounded by the surface. The surface trans-fer position control strategy is: u stmi v stmi w stmi = U agtmax R − ln (cid:16) γ + µ i ) R − r i ( (cid:107) n i (cid:107) − r i ) (cid:17) ˆ n i (cid:107) ln (cid:16) γ + µ i ) R − r i ( (cid:107) n i (cid:107) − r i ) (cid:17) ˆ n i (cid:107) . (23)10s the agents ascend to a point at which R does not in-tersect C , sensing information is not gathered in avoid-ance mode and thus the avoidance orientation control issimply [ q stmi r stmi s stmi ] T = [0 0 0] T .Agents are said to have converged upon their newly as-signed surface when | ln (cid:107) n i (cid:107)− r i ( γ + µ i ) R − r i | < ε . At this point,each agent shall transition back to its prior mode as de-scribed in the following two scenarios.(1) If agent i had been in either particle intercept orpartition transfer mode before the deadlock, it shallresume that mode and continue along a geodesictrajectory on the newly assigned surface until itreaches the projection of its destination. At thispoint, the condition that (cid:107) p i − proj C µ i p id (cid:107) ≤ ε triggers a reset µ i := 0 concurrent with a transi-tion back to surface transfer mode thus allowingthe agent to transfer back to C . The agent thenresumes coverage of C in its prior mode. For addi-tional details on flag conditions in these transitions,see guards G ( ζ i , ζ i ) , G ( ζ i , ζ i ) , G ( ζ i , ζ i ), and G ( ζ i , ζ i ) of our hybrid automaton as presented inthe appendix.(2) If agent i had been in local coverage mode beforethe deadlock, it shall then transition back to localcoverage mode concurrent with reset µ i := 0. Thistransition is dependent upon the conditions that f i = 0 and that the agent is operating within its as-signed partition. The agent shall oscillate betweenlocal coverage and surface transfer at an altitude of R above C until i pr has moved along C to resolvethe deadlock. At this point, the local coverage con-troller shall attract agent i back to the surface.While similar work on multi-agent systems often invokeavoidance barrier functions to encode collision avoid-ance, such as in Panagou et al. (2016), it may be impos-sible to bound the time that agents spend avoiding oneanother in these maneuvers—especially when the algo-rithm is scaled to many agents. In contrast, our tech-nique results in an explicit bound on path length to anintruder as was proven in Lemma 1. With an additionalassumption on the size of agents, we can establish a guar-antee on collision avoidance for agents in surface transfermode. Theorem 3
For agents { i, j } ∈ ˜ j , the condition thatmin( R ˜ j ) > r i + 2 r j implies that i does not collide with j . PROOF.
Consider the case in which i (cid:54) = i pr and j (cid:54) = i pr . Both agents operate in accordance with (23) andfollow trajectories along ˆ n i and ˆ n j respectively. Bothunit vectors are normal to surface C µ i , an ellipsoid ofrevolution, and thus diverge from one another away from C µ i . Agents i and j shall enter surface transfer mode atan instant when (cid:107) p i − p j (cid:107) ≥ min (cid:16) R ˜ j (cid:17) and their distance shall tend to increase under (23). Thus min( R ˜ j ) > r i + r j and subsequently min( R ˜ j ) > r i + 2 r j imply that theyavoid collision.Consider the case in which i = i pr and thus j (cid:54) = i pr . Inthe instant that j transitions to surface transfer modewe have that (cid:107) p i − p j (cid:107) ≥ min (cid:16) R ˜ j (cid:17) . Thus the distancefor i to travel until collision is greater than or equal tomin (cid:16) R ˜ j (cid:17) − r i − r j . This straight line path for i is aconservative estimate as the true path is curved. Colli-sion will be avoided if agent j , whose path is normal tothe surface, may cover a distance r i + r j before i coversmin (cid:16) R ˜ j (cid:17) − r i − r j . As j moves at speed U agtmax and i ’stangential speed is upper bounded by U agtmax , this con-dition is satisfied if min (cid:16) R ˜ j (cid:17) − r i − r j > r i + r j . Thismay equivalently be written as min (cid:16) R ˜ j (cid:17) > r i + 2 r j .These arguments apply to the case in which j = i pr and i (cid:54) = i pr as well. This concludes the proof. A simulation was performed in MATLAB to verify theefficacy of the algorithm. Four agents are deployed tomonitor the surface of an ellipsoid of revolution, C , whoseradius in the xy -plane is 80 and whose radius in the z -plane is 20. For each agent, R = 10, r i = 1, α i = 30 ◦ , k u = 1, k v = 5, k w = 1, k r = 0 . k s = 0 .
1, ¯ r i = 0 . s i = 0 . U agtmax = 6, and T (cid:63) = 792. Upon initializationof the simulation, C was set to a fully covered level of C (cid:63) = 20 which would begin decaying upon detection ofthe first intruder k ∈ { , ..., } at t = 600 seconds. Thefour agents were deployed from F sequentially at times t = 0, t = T (cid:63) , t = T (cid:63) , and t = T (cid:63) seconds respectively.Upon deployment, each agent was initialized in local cov-erage mode with Φ i = 0, Θ i = π , and Ψ i = 0. Intruderstraveled in random directions with U intmax = 0 .
7, thoughwere still constrained to always impact the surface, andwere generated every 35 seconds beginning at t = T (cid:63) seconds. The detection system had a lower bound onrange R det = 80, decay rate parameter λ k = 0 .
05, andmeasurement variances σ ρ = 0 . σ θ = 0 .
25 deg ,and σ ψ = 0 .
25 deg respectively Agents were able to suc-cessfully intercept all particles along their geodesic tra-jectories while actively avoiding collision over the entireduration of the attack (see Fig. 7 and Fig. 8); however,it should be noted that one avoidance anomaly occurredbefore the initial intruder was generated during the in-terval of t = 418 − R = 10 and our simulation time step size was 1, itis clear that this anomaly occurred due to a selectionof U agtmax that was too large relative to the time step. Ina continuous time implementation, a transition to sur-face transfer mode would have occurred between t = 418and t = 419 thus preventing collision. Aside from thisanomaly, the simulation parameters adequately approx-imated the continuous time agent kinematics.11 ig. 7. Agent i = 3, indicated with green S i , is on a collision course with agent i = 1, indicated with blue S i , during theinterval from t = 2125 to t = 2130. At t = 2135, agent i = 1 has transitioned to surface transfer mode and is following atrajectory normal to the surface while agent i = 4, indicated with black S i , follows a collision course through t = 2160. Agent i = 4 transitions to surface transfer mode as well leading to the conditions that µ = 2 and µ = 1, i.e., agent i = 1 is assignedto the second tier of avoidance surfaces at a higher altitude than i = 4 as illustrated at t = 2180. Both agents proceed alongtheir respective C µ i towards their destination with i = 1 having arrived and transferred back to C before t = 2195. Note thatagent trajectories for t ≥ i = 2 follows its geodesic trajectory to thepredicted impact point of particle k over time lapse (a)-(d).The true trajectory of the particle is indicated in red andthe estimated trajectory in green. The coverage error on C , normalized with respect to themaximum error in which all of C takes a value of zerofor Q ( t, ˜ p ), as well as the minimum inter-agent distanceover time are presented in Fig. 9. The error tends to spikeupon particle detections with agents effectively curtail-ing these spikes as they cover around the vicinity of pre-dicted impact points in local coverage mode. Two par-ticularly large spikes occur at t = 3225 and t = 5180respectively. These anomalies are each associated withparticle impacts occurring close to the equator of the el-lipsoid where even small values of σ θ and σ ψ result inan estimated particle trajectory that does not initially intersect C thus delaying an agent assignment. In bothcases, the estimated trajectory did eventually intersect C with enough time to allow for agent interception. How-ever, this delay in assignment significantly reduced thetime the agent spent exploring in the vicinity of the pre-dicted impact point thus contributing to a noticeablerise in the coverage error. One potential solution to thisproblem would be to prescribe some boundary toleranceto our surface C thus loosening our definition of an im-pacting particle for the sake of measurement uncertainty. Fig. 9. The coverage error, normalized to the maximum pos-sible value, is presented. Anomalies are observed at t = 3225and t = 5180 respectively. The minimum distance betweenany two agents at any given time is also presented with ananomaly observed at t = 419. ig. 10. A typical agent’s hybrid modes are presented overtime. Abbreviations from top to bottom refer to surfacetransfer mode, partition transfer mode, particle interceptmode, return to base, and local coverage mode respectively. Agent i = 1’s operating modes with respect to time arepresented in Fig. 10. It should be noted that the mostfrequent transition out of particle intercept mode is tolocal coverage mode. This corresponds to an agent ar-riving at the estimated impact point of a particle andthen surveying the local area up until the moment of im-pact. As the agent surveys it tends to hit the ε proxim-ity boundary to the impact point thus requiring a shortoperation in particle intercept mode to direct the agentback within the ε boundary.To demonstrate scalability, we have made an additionalsimulation with 100 agents available online at: https://1drv.ms/f/s!AsiVOlIEkwNEgX2o1eV2hJ_bbaQU . In this paper, we presented a hybrid formulation forthe persistent coverage problem in an environment sub-ject to stochastic intruders. This formulation was moti-vated in part by extravehicular applications of the NASAMini AERCam. Agents operated with finite power re-sources and were required to periodically return to a re-fueling station while patrolling assigned latitude parti-tions along the surface of an ellipsoid. Formal guaran-tees were established on the ability of agents to inter-cept all intruders and the efficacy of the algorithm wasdemonstrated in simulation. This approach succeeds ourprevious work in Bentz & Panagou (2017) and Bentz& Panagou (2018) by extending the guarantees on in-truder interception to an arbitrary number of collisionavoidance maneuvers. It also removes singularities in thesensing function definition and redefines the guard con-ditions in a manner that supports a more effective useof local coverage around the vicinity of intruder impactpoints.
Our local coverage control laws are derived via differ-entiation of (10) of which we seek to reduce the rateof growth. It is a volume integral, so a few mathemat-ical preliminaries are required. Recall the generalizedtransport theorem (GTT) (Slattery 1999): ddt (cid:82) R ( s ) f dV = (cid:82) R ( s ) ∂f∂t dV + (cid:82) S ( s ) f v ( s ) · n dA , where f is any scalar-,vector-, or tensor-valued function of position and time, S ( s ) is the boundary of the volume R ( s ) over which f isintegrated, n is the unit vector normal to the boundary,and v ( s ) is the velocity of the boundary. V and A referto volume and area respectively. Invoking GTT allowsfor differentiation of (10) with respect to time:˙ E ( t ) = (cid:90) D h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) (cid:18) − ∂Q ( t, ˜ p ) ∂t (cid:19) d ˜ p + (cid:90) ∂ D (cid:0) h ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) (cid:1) v ( s ) · n dA, (24)where ∂ D is the boundary of D . D is time invariant andthus v ( s ) = 0. (24) reduces to ˙ E ( t ) = (cid:82) D h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) (cid:16) − ∂Q ( t, ˜ p ) ∂t (cid:17) d ˜ p , which expands to: ˙ E ( t ) = − (cid:90) D h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) (cid:18) N (cid:88) i =1 S i (˜ q i ( t ) , ˜ p ) C (˜ p ) − N p (cid:88) k =1 Λ k ( t, ˜ p ) C (˜ p ) (cid:19) d ˜ p = N (cid:88) i =1 (cid:90) D − h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) S i (˜ q i ( t ) , ˜ p ) C (˜ p ) d ˜ p (cid:124) (cid:123)(cid:122) (cid:125) =ˆ e i ( t ) − N p (cid:88) k =1 (cid:90) D − h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p ))Λ k ( t, ˜ p ) C (˜ p ) d ˜ p (cid:124) (cid:123)(cid:122) (cid:125) =˜ e k ( t ) = N (cid:88) i =1 ˆ e i ( t ) − N p (cid:88) k =1 ˜ e k ( t ) . (25) ˆ e i ( t ) is the rate of change of the coverage error due tothe motion of the agents while ˜ e k ( t ) is the rate of changeof the coverage error due to a contrived information de-cay surrounding the predicted impact point of particle k on C . Our strategy is to control the agents’ kinemat-ics, recovered in the derivative of ˆ e i ( t ), to decrease (25).Note that we do not presume that our local coveragestrategy provides any additional bounds on (10). Nor dowe provide guarantees on the rate of growth of this con-trived quantity. Curtailing the growth of the coverageerror simply imparts the desired effect of active explo-ration in the vicinity of impact points into our system.Using this strategy, the agents actively seek to increasetheir rate of coverage by rotating and/or translating S i to be encompass the most uncovered space in the localvicinity.Taking the derivative of ˆ e i ( t ) with respect to time yields:13 xpand ddt ( S i (˜ q i ( t ) , ˜ p )): ddt ( S i (˜ q i ( t ) , ˜ p )) = ∂S i ∂x i ˙ x i ( t ) + ∂S i ∂y i ˙ y i ( t ) + ∂S i ∂z i ˙ z i ( t ) + ∂S i ∂ Ψ i ˙Ψ i ( t ) + ∂S i ∂ Θ i ˙Θ i ( t ) = (cid:18) ∂S i ∂x i cos Θ cos Ψ + ∂S i ∂y i cos Θ sin Ψ − ∂S i ∂z i sin Θ (cid:19) u i ( t )+ (cid:18) ∂S i ∂x i (sin Φ sin Θ cos Ψ − cos Φ sin Ψ) + ∂S i ∂y i (sin Φ sin Θ sin Ψ + cos Φ cos Ψ) + ∂S i ∂z i sin Φ cos Θ (cid:19) v i ( t )+ (cid:18) ∂S i ∂x i (cos Φ sin Θ cos Ψ + sin Φ sin Ψ) + ∂S i ∂y i (cos Φ sin Θ sin Ψ − sin Φ cos Ψ) + ∂S i ∂z i cos Φ cos Θ (cid:19) w i ( t )+ (cid:18) ∂S i ∂ Ψ i sin Φ sec Θ + ∂S i ∂ Θ i cos Φ (cid:19) r i ( t ) + (cid:18) ∂S i ∂ Ψ i cos Φ sec Θ − ∂S i ∂ Θ i sin Φ (cid:19) s i ( t ) . (27)Now introduce the following definitions: a i ( t, Q ( t, ˜ p )) = (cid:90) Di h (cid:48)(cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) S i (˜ q i ( t ) , ˜ p ) C (˜ p ) ∂Q ( t, ˜ p ) ∂t d ˜ p, (28) a i ( t, Q ( t, ˜ p )) = (cid:90) Di h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) C (˜ p ) (cid:18) ∂S i ∂x i cos Θ cos Ψ + ∂S i ∂y i cos Θ sin Ψ − ∂S i ∂z i sin Θ (cid:19) d ˜ p, (29) a i ( t, Q ( t, ˜ p )) = (cid:90) Di h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) C (˜ p ) (cid:18) ∂S i ∂x i (sin Φ sin Θ cos Ψ − cos Φ sin Ψ) + ∂S i ∂y i (sin Φ sin Θ sin Ψ + cos Φ cos Ψ) + ∂S i ∂z i sin Φ cos Θ (cid:19) d ˜ p, (30) a i ( t, Q ( t, ˜ p )) = (cid:90) Di h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) C (˜ p ) (cid:18) ∂S i ∂x i (cos Φ sin Θ cos Ψ + sin Φ sin Ψ) + ∂S i ∂y i (cos Φ sin Θ sin Ψ − sin Φ cos Ψ) + ∂S i ∂z i cos Φ cos Θ (cid:19) d ˜ p, (31) a i ( t, Q ( t, ˜ p )) = (cid:90) Di h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) C (˜ p ) (cid:18) ∂S i ∂ Ψ i sin Φ sec Θ + ∂S i ∂ Θ i cos Φ (cid:19) d ˜ p, (32) a i ( t, Q ( t, ˜ p )) = (cid:90) Di h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) C (˜ p ) (cid:18) ∂S i ∂ Ψ i cos Φ sec Θ − ∂S i ∂ Θ i sin Φ (cid:19) d ˜ p. (33)One can then rewrite (26) as:˙ˆ e i ( t ) = a i ( t, Q ( t, ˜ p )) − u i ( t ) a i ( t, Q ( t, ˜ p )) − v i ( t ) a i ( t, Q ( t, ˜ p )) − w i ( t ) a i ( t, Q ( t, ˜ p )) − r i ( t ) a i ( t, Q ( t, ˜ p )) − s i ( t ) a i ( t, Q ( t, ˜ p )) . (34) ˙ˆ e i ( t ) = (cid:90) D i (cid:18) h (cid:48)(cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) S i (˜ q i ( t ) , ˜ p ) C (˜ p ) ∂Q ( t, ˜ p ) ∂t − h (cid:48) ( C (cid:63) C (˜ p ) − Q ( t, ˜ p )) ddt ( S i (˜ q i ( t ) , ˜ p )) C (˜ p ) (cid:19) d ˜ p. (26) The sensing footprint is independent of Φ i assuming thatthe centerline of the spherical sector is aligned with theˆ x B i axis. ddt ( S i (˜ q i ( t ) , ˜ p )) is expanded in (27) and throughthe definitions in (28-33) one may restate (26) as (34). Ifone were to command zero inputs to this system, it be-comes clear that a i ( t, Q ( t, ˜ p )) may be physically inter-preted as the rate at which the coverage rate is reducingdue to information saturation at any particular positionand orientation of the sensing footprint, S i . As the foot-print remains stationary, there are diminishing returnson the value of newly acquired information. Thus, theadditional terms in (34) allow for the coverage rate to beincreased by mobilizing the sensor. One strategy is thatof (11). To provide a compact notation in this section, define ´ f i = x i ( x C ,r + r i ) + y i ( y C ,r + r i ) + z i ( z C ,r + r i ) . The coverage strategyfor agent i is represented by the hybrid automaton in Fig.2, described by the following entities (Lygeros 2004): • A set of discrete states: Z i = { ζ i , ζ i , ζ i , ζ i , ζ i } , • A set of continuous states: ˜ q i = { x i , y i , z i , Φ i , Θ i , Ψ i } , • A vector field: f ( ζ i , ˜ q i ) = R (cid:2) u loci v loci w loci r loci s loci (cid:3) T ,f ( ζ i , ˜ q i ) = R (cid:2) u rtbi v rtbi w rtbi q rtbi r rtbi s rtbi (cid:3) T ,f ( ζ i , ˜ q i ) = R (cid:104) u pimi v pimi w pimi q pimi r pimi s pimi (cid:105) T ,f ( ζ i , ˜ q i ) = R (cid:2) u ptmi v ptmi w ptmi q ptmi r ptmi s ptmi (cid:3) T ,f ( ζ i , ˜ q i ) = R [ u stmi v stmi w stmi T where R = (cid:34) R R (cid:35) , • A set of initial states: { ζ i } × { ˜ q i ∈ R | p i = F∧ Φ i ∈ [ − π, + π ] ∧ Θ i ∈ (cid:2) − π , + π (cid:3) ∧ Ψ i ∈ [ − π, + π ] } , • A domain:
Dom ( ζ i ) = { ˜ q i ∈ R | ´ f i ≥ ∧ (cid:0) i p ∈ { , ..., N } = ⇒ ¯ z i p − ≤ z i ≤ ¯ z i p − (cid:1) } ,Dom ( ζ i ) = { ˜ q i ∈ R | ´ f i ≥ } ,Dom ( ζ i ) = { ˜ q i ∈ R | ´ f i ≥ } ,Dom ( ζ i ) = { ˜ q i ∈ R | ´ f i ≥ ∧ (cid:0) i p ∈ { , ..., N } = ⇒ z i < ¯ z i p − ∨ z i > ¯ z i p − (cid:1) } ,Dom ( ζ i ) = { ˜ q i ∈ R | ´ f i ≥ } , • A set of edges: E = { ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , ( ζ i , ζ i ) , } , • A set of guard conditions:14 ( ζ i , ζ i ) = { i p = 1 ∧ t i F ≥ T (cid:63) − πg C U agtmax (cid:0) x C ,r + z C ,r (cid:1) } ,G ( ζ i , ζ i ) = { ( ∃ k | i = i k ) ∧ (cid:16) (cid:107) p i − proj C µ i p id (cid:107) > ε (cid:17) } ,G ( ζ i , ζ i ) = { i p (cid:54) = 1 ∧ (cid:0) z i < ¯ z i p − ∨ z i > ¯ z i p − (cid:1) } ,G ( ζ i , ζ i ) = {(cid:107) p i − p j (cid:107) ≤ R ∧ i pr (cid:54) = argmax ˜ j (cid:16) t ˜ j F (cid:17) } ,G ( ζ i , ζ i ) = {(cid:107) p i − F(cid:107) ≤ ε ∧ t i F = T (cid:63) } ,G ( ζ i , ζ i ) = { (cid:0) t ≥ t ck ∧ (cid:0)(cid:0) i p ∈ { , ..., N } ∧ ¯ z i p − ≤ z i ≤ ¯ z i p − (cid:1) ∨ (cid:16) i p = 1 ∧ t i F < T (cid:63) − πg C U agtmax (cid:0) x C ,r + z C ,r (cid:1)(cid:17)(cid:17)(cid:17) ∨ (cid:16) t < t ck ∧ (cid:107) p i − proj C µ i p id (cid:107) ≤ ε ∧ µ i = 0 (cid:17) } ,G ( ζ i , ζ i ) = {(cid:107) p i − proj C µ i p id (cid:107) ≤ ε ∧ µ i = 0 ∧ t ≥ t ck ∧ i p = 1 ∧ t i F ≥ T (cid:63) − πg C U agtmax (cid:0) x C ,r + z C ,r (cid:1) } G ( ζ i , ζ i ) = {(cid:107) p i − proj C µ i p id (cid:107) ≤ ε ∧ µ i = 0 ∧ t ≥ t ck ∧ i p ∈ { , ..., N } ∧ (cid:0) z i < ¯ z i p − ∨ z i > ¯ z i p − (cid:1) } ,G ( ζ i , ζ i ) = G ( ζ i , ζ i ) ∨ {(cid:107) p i − proj C µ i p id (cid:107) ≤ ε ∧ (cid:16) (cid:107) p i − p j (cid:107) > R, ∀ j ∨ i pr = argmax ˜ j (cid:16) t ˜ j F (cid:17)(cid:17) ∧ µ i > } ,G ( ζ i , ζ i ) = { i p = 1 ∨ (cid:0) i p (cid:54) = 1 ∧ ¯ z i p − ≤ z i ≤ ¯ z i p − (cid:1) } ,G ( ζ i , ζ i ) = G ( ζ i , ζ i ) ,G ( ζ i , ζ i ) = G ( ζ i , ζ i ) ∨ {(cid:107) p i − proj C µ i p id (cid:107) ≤ ε ∧ (cid:16) (cid:107) p i − p j (cid:107) > R, ∀ j ∨ i pr = argmax ˜ j (cid:16) t ˜ j F (cid:17)(cid:17) ∧ µ i > } ,G ( ζ i , ζ i ) = { f i = 0 ∧ | ln (cid:16) (cid:107) n i (cid:107)− r i ( γ + µ i ) R − r i (cid:17) | < ε ∧ (cid:0) i p = 1 ∨ (cid:0) i p ∈ { , ..., N } ∧ ¯ z i p − ≤ z i ≤ ¯ z i p − (cid:1)(cid:1) } ,G ( ζ i , ζ i ) = { f i = 1 ∧ | ln (cid:16) (cid:107) n i (cid:107)− r i ( γ + µ i ) R − r i (cid:17) | < ε } ,G ( ζ i , ζ i ) = { f i = 0 ∧ | ln (cid:16) (cid:107) n i (cid:107)− r i ( γ + µ i ) R − r i (cid:17) | < ε ∧ (cid:0) i p (cid:54) = 1 ∧ (cid:0) z i < ¯ z i p − ∨ z i > ¯ z i p − (cid:1)(cid:1) } . • Additional parameters include a clock set: C = { t i F } , aflag: f i ∈ { , } , an assignment index µ i = { , ..., N − } and, • A reset map: R ( ζ i , ζ i , f i ) = { } , R ( ζ i , ζ i , µ i ) = { µ i + 1 } , R ( ζ i , ζ i , t i F ) = { } , R ( ζ i , ζ i , f i ) = { t ≥ t ck ; 1 otherwise } , R ( ζ i , ζ i , f i ) = { } ,R ( ζ i , ζ i , f i ) = { } , R ( ζ i , ζ i , µ i ) = { (cid:107) p i − proj C µ i p id (cid:107)≤ ε ∧ (cid:16) (cid:107) p i − p j (cid:107) > R, ∀ j ∨ i pr = argmax ˜ j (cid:16) t ˜ j F (cid:17)(cid:17) ∧ µ i > µ i + 1 otherwise } , R ( ζ i , ζ i , f i ) = { } ,R ( ζ i , ζ i , µ i ) = { (cid:107) p i − proj C µ i p id (cid:107) ≤ ε ∧ (cid:16) (cid:107) p i − p j (cid:107) > R, ∀ j ∨ i pr = argmax ˜ j (cid:16) t ˜ j F (cid:17)(cid:17) ∧ µ i > µ i + 1 otherwise } , R ( ζ i , ζ i , µ i ) = { } , and continuousstates do not reset between transitions. References
Beard, R. W. (2008), Quadrotor dynamics and control.lecture notes. http://scholarsarchive.byu.edu/cgi/viewcontent.cgi?article=2324&context=facpub .Bentz, W. & Panagou, D. (2017), Persistent coverage ofa two-dimensional manifold subject to time-varyingdisturbances, in ‘Proc. of the 56th IEEE Conferenceon Decision and Control’, Melbourne, Australia.Bentz, W. & Panagou, D. (2018), Energy-aware persis-tent coverage and intruder interception in 3D dynamicenvironments, in ‘Proc. of the 2018 American ControlConference’, Milwaukee, WI.Bokareva, T., Hu, W., Kanhere, S., Ristic, B., Gordon,N., Bessell, T., Rutten, M. & Jha, S. (2006), Wireless sensor networks for battlefield surveillance, in ‘Pro-ceedings of the land warfare conference’, pp. 1–8.Cheng, P., Keller, J. & Kumar, V. (2008), Time-optimalUAV trajectory planning for 3D urban structure cov-erage, in ‘Proc. of the 2008 IEEE/RSJ InternationalConference on Intelligent Robots and Systems’, Nice,France, pp. 2750–2757.Choset, H. & Kortenkamp, D. (1999), ‘Path planningand control for free-flying inspection robot in space’, Journal of Aerospace Engineering (2), 74–81.Cortes, J., Mart´ınez, S., Karatas, T. & Bullo, F. (2004),‘Coverage control for mobile sensing networks’, IEEETrans. on Robotics and Automation (2), 243–255.Fredrickson, S. E., Duran, S., Howard, N. & Wa-genknecht, J. D. (2004), Application of the mini AER-cam free flyer for orbital inspection, in ‘Defense andSecurity’, International Society for Optics and Pho-tonics, pp. 26–35.Hexsel, B., Chakraborty, N. & Sycara, K. (2013), Dis-tributed coverage control for mobile anisotropic sen-sor networks, Technical Report CMU-RI-TR-13-01,Robotics Institute, Pittsburgh, PA.Hokayem, P., Stipanov´ıc, D. & Spong, M. (2007), Onpersistent coverage control, in ‘Proc. of the 46th IEEEConference on Decision and Control’, New Orleans,LA, USA, pp. 6130–6135.Hollinger, G. A., Englot, B., Hover, F. S., Mitra, U. &Sukhatme, G. S. (2013), ‘Active planning for under-water inspection and the benefit of adaptivity’, TheInternational Journal of Robotics Research (1), 3–18.H¨ubel, N., Hirche, S., Gusrialdi, A., Hatanaka, T., Fu-jita, M. & Sawodny, O. (2008), Coverage control withinformation decay in dynamic environments, in ‘Proc.of the 17th IFAC World Congress’, Seoul, South Ko-rea, pp. 4180–4185.Hussein, I. I. & Stipanovi´c, D. M. (2007), ‘Effective cov-erage control for mobile sensor networks with guar-anteed collision avoidance’, IEEE Trans. on ControlSystems Technology (4), 642–657.Ivory, J. (1798), ‘VIII. A new series for the rectificationof the ellipsis; together with some observations on theevolution of the formula ( a + b − ab cos θ ) n ’, Trans-actions of the Royal Society of Edinburgh (2), 177190.Liu, B., Dousse, O., Nain, P. & Towsley, D. (2013),‘Dynamic coverage of mobile sensor networks’, IEEETransactions on Parallel and Distributed systems (2), 301–311.Lygeros, J. (2004), Lecture notes on hybrid sys-tems. https://robotics.eecs.berkeley.edu/~sastry/ee291e/lygeros.pdf .Mitchell, D., Corah, M., Chakraborty, N., Sycara, K. &Michael, N. (2015), Multi-robot long-term persistentcoverage with fuel constrained robots, in ‘Proc. of the2015 IEEE International Conference on Robotics andAutomation’, IEEE, pp. 1093–1099.Murphy, R. R., Tadokoro, S., Nardi, D., Jacoff, A., Fior-ini, P., Choset, H. & Erkmen, A. M. (2008), Search andrescue robotics, in ‘Springer Handbook of Robotics’,15pringer Berlin Heidelberg, pp. 1151–1173.Palacios-Gas´os, J. M., Montijano, E., Sagues, C. &Llorente, S. (2016 a ), Multi-robot persistent coverageusing branch and bound, in ‘Proc. of the 2016 Amer-ican Control Conference’, pp. 5697–5702.Palacios-Gas´os, J. M., Montijano, E., Sag¨u´es, C. &Llorente, S. (2016 b ), Multi-robot persistent coveragewith optimal times, in ‘Proc. of the 55th IEEE Confer-ence on Decision and Control’, IEEE, pp. 3511–3517.Palacios-Gas´os, J. M., Talebpour, Z., Montijano, E.,Sag¨u´es, C. & Martinoli, A. (2017), Optimal path plan-ning and coverage control for multi-robot persistentcoverage in environments with obstacles, in ‘Proc. ofthe 2017 IEEE International Conference on Roboticsand Automation’, pp. 1321–1327.Panagou, D., Stipanovi´c, D. M. & Voulgaris, P. G.(2016), ‘Distributed coordination control for multi-robot networks using lyapunov-like barrier functions’, IEEE Transactions on Automatic Control, to appear (3), 617–632.Pasqualetti, F., Zanella, F., Peters, J. R., Spindler, M.,Carli, R. & Bullo, F. (2014), ‘Camera network coordi-nation for intruder detection’, IEEE Transactions onControl Systems Technology (5), 1669–1683.Slattery, J. C. (1999), Advanced Transport Phenomena ,Cambridge UP.Smith, R. N., Schwager, M., Smith, S. L., Jones, B. H.,Rus, D. & Sukhatme, G. S. (2011), ‘Persistent oceanmonitoring with underwater gliders: Adapting sam-pling resolution’,
Journal of Field Robotics (5), 714–741.Smith, S. L., Schwager, M. & Rus, D. (2012), ‘Persis-tent robotic tasks: Monitoring and sweeping in chang-ing environments’, IEEE Transactions on Robotics (2), 410–426.Song, C., Feng, G., Fan, Y. & Wang, Y. (2011), ‘Decen-tralized adaptive awareness coverage control for multi-agent networks’, Automatica (12), 2749 – 2756.Song, C., Liu, L., Feng, G., Wang, Y. & Gao, Q. (2013),‘Persistent awareness coverage control for mobile sen-sor networks’, Automatica (6), 1867–1873.Stipanovi´c, D. M., Valicka, C., Tomlin, C. J. & Bewley,T. R. (2013), ‘Safe and reliable coverage control’, Nu-merical Algebra, Control and Optimization , 31–48.Vincenty, T. (1975 a ), ‘Direct and inverse solutions ofgeodesics on the ellipsoid with application of nestedequations’, Survey review (176), 88–93.Vincenty, T. (1975 b ), Geodetic inverse solution betweenantipodal points. Scanned by Charles Karney from thecopy in R.H. Rapp’s library at Ohio State University.The report is a work of the U.S. Government and sois in the public domain.Xie, L. & Zhang, X. (2013), 3D clustering-based cam-era wireless sensor networks for maximizing lifespanwith minimum coverage rate constraint, in ‘Proc. ofthe 2013 IEEE Global Communications Conference(GLOBECOM)’, pp. 298–303.Yu, J., Karaman, S. & Rus, D. (2015), ‘Persistent mon-itoring of events with stochastic arrivals at multiple stations’, IEEE Transactions on Robotics31