A Hybrid Approach to Persistent Coverage in Stochastic Environments

TL;DR

This paper presents a hybrid control strategy enabling power-constrained 3D agents to persistently monitor a 2D manifold, intercept intruders, and manage energy, with formal guarantees and simulation validation.

Contribution

It introduces a novel hybrid control approach for persistent coverage and intruder interception in stochastic environments with formal performance guarantees.

Findings

Successful simulation demonstrations of the algorithm.

Formal proofs of intruder interception guarantees.

Effective management of agent power constraints.

Abstract

This paper considers the persistent coverage of a 2-D manifold that has been embedded in 3-D space. The manifold is subject to continual impact by intruders which travel at constant velocities along arbitrarily oriented straight-line trajectories. The trajectories of intruders are estimated online with an extended Kalman filter and their predicted impact points contribute normally 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 deploy from and return to a single charging station on the manifold. Guarantees on intruder interception with respect to agent power lifespans are formally proven. The efficacy of the algorithm is demonstrated through simulation.

Equations

$$\begin{gathered} \left[\begin{array}{c}\dot{x_i}\\ \dot{y_i}\\ \dot{z_i}\end{array}\right]=[\begin{array}{cc}\cos {\Theta }_i\cos {\Psi }_i& \sin {\Phi }_i\sin {\Theta }_i\cos {\Psi }_i-\cos {\Phi }_i\sin {\Psi }_i\\ \cos {\Theta }_i\sin {\Psi }_i& \sin {\Phi }_i\sin {\Theta }_i\sin {\Psi }_i+\cos {\Phi }_i\cos {\Psi }_i\\ -\sin {\Theta }_i& \sin {\Phi }_i\cos {\Theta }_i\end{array} \\ \begin{array}{c}\cos {\Phi }_i\sin {\Theta }_i\cos {\Psi }_i+\sin {\Phi }_i\sin {\Psi }_i\\ \cos {\Phi }_i\sin {\Theta }_i\sin {\Psi }_i-\sin {\Phi }_i\cos {\Psi }_i\\ \cos {\Phi }_i\cos {\Theta }_i\end{array}]\left[\begin{array}{c}{u}_i\\ v_i\\ w_i\end{array}\right], \end{gathered}$$

$$\left[\begin{array}{c}\dot{{\Phi }_i}\\ \dot{{\Theta }_i}\\ \dot{{\Psi }_i}\end{array}\right]=\left[\begin{array}{ccc}1& \sin {\Phi }_i\tan {\Theta }_i& \cos {\Phi }_i\tan {\Theta }_i\\ 0& \cos {\Phi }_i& -\sin {\Phi }_i\\ 0& \sin {\Phi }_i\sec {\Theta }_i& \cos {\Phi }_i\sec {\Theta }_i\text{par}\end{array}\right]\left[\begin{array}{c}{q}_i\\ r_i\\ s_i\end{array}\right],$$

$$c_{1i}$$

$$c_{2i}$$

$$S_i({\tilde{q}}_i,\tilde{p})=\{\begin{array}{ll}\text{hfill}\frac{C_{1i}{C}_{2i}}{C_{1i}+C_{2i}},\text{hfill}& \text{if card(Cˉi)<2∧rp~i/pi>0;}\\ \text{hfill}0,\text{hfill}& \text{otherwise,}\end{array}$$

$$Q_i(t,\tilde{p})={\int }_0^t{S}_i\left({\tilde{q}}_i(\tau ),\tilde{p}\right)C\left(\tilde{p}\right)d\tau ,$$

$$Q(t,\tilde{p})=\sum _{i=1}^N{Q}_i(t,\tilde{p})-\sum _{k=1}^{N_p}{\int }_0^t{\Lambda }_k\left(\tau ,\tilde{p}\right)C\left(\tilde{p}\right)d\tau .$$

$${\dot{\tilde{q}}}_k\left(t\right)$$

$${\tilde{z}}_k\left(t\right)$$

$$\displaystyle\Lambda_{k}\left(t,\tilde{p}\right)=\int_{0}^{T_{H,k}(t)}\lambda_{k}\mathcal{N}\bigl{(}\tilde{p}^{\prime}_{k}\left(t+\tau\right),\mathbf{P^{\prime}}_{k}\left(t+\tau\right)\bigr{)}\,d\tau,$$

$$E(t)=\underset{\mathcal{D}}{\int }h\left(C^{\star }C\left(\tilde{p}\right)-Q(t,\tilde{p})\right)d\tilde{p},$$

$$u_i^{loc}$$

$$v_i^{loc}$$

$$w_i^{loc}$$

$$r_i^{loc}$$

$$s_i^{loc}$$

$${\rho }_{l,i}$$

$${\rho }_{a,i}$$

$$\left[\begin{array}{c}{u}_i^{pim}\\ v_i^{pim}\\ w_i^{pim}\end{array}\right]=U_{max}^{agt}{\mathcal{R}}_1^{-1}\frac{{\hat{\nu }}_i-\ln \left(\frac{1}{\left(\gamma +{\mu }_i\right)R-\text{mathcalligra}{r}_i}\left(\Vert n_i\Vert -\text{mathcalligra}{r}_i\right)\right){\hat{n}}_i}{\Vert {\hat{\nu }}_i-\ln \left(\frac{1}{\left(\gamma +{\mu }_i\right)R-\text{mathcalligra}{r}_i}\left(\Vert n_i\Vert -\text{mathcalligra}{r}_i\right)\right){\hat{n}}_i\Vert }.$$

$$\left[\begin{array}{c}{q}_i^{pim}\\ r_i^{pim}\\ s_i^{pim}\end{array}\right]={\mathcal{R}}_2^{-1}\left[\begin{array}{c}0\\ \arcsin \left({\hat{n}}_i\cdot {\hat{z}}_{\mathcal{G}}\right)-{\Theta }_i\\ atan22\left(-{\hat{n}}_i\cdot {\hat{y}}_{\mathcal{G}},-{\hat{n}}_i\cdot {\hat{x}}_{\mathcal{G}}\right)-{\Psi }_i\end{array}\right],$$

$$g_{{\mathcal{C}}_{N-1}}=\left[1+\sum _{n=1}^{\infty }{\left(\frac{\left(2n-1\right)!!}{2^{n}n!}\right)}^2\frac{{\left(\frac{x_{{\mathcal{C}}_{N-1},r}-z_{{\mathcal{C}}_{N-1},r}}{x_{{\mathcal{C}}_{N-1},r}+z_{{\mathcal{C}}_{N-1},r}}\right)}^{2n}}{{\left(2n-1\right)}^2}\right],$$

$${\mathcal{P}}_{geo}\le {\mathcal{P}}_{lat}+{\mathcal{P}}_{long}.$$

$${\mathcal{P}}_{lat}$$

$${\mathcal{P}}_{long}$$

$$A_{\mathcal{C}}=2\pi x_{\mathcal{C},r}^2\left(1+\frac{1+\left(1-\frac{z_{\mathcal{C},r}^2}{x_{\mathcal{C},r}^2}\right)}{\left(\sqrt{1-\frac{z_{\mathcal{C},r}^2}{x_{\mathcal{C},r}^2}}\right)}\mathrm{artanh}\left(\sqrt{1-\frac{z_{\mathcal{C},r}^2}{x_{\mathcal{C},r}^2}}\right)\right).$$

$$\frac{A_{\mathcal{C}}}{N-1}={\int }_{z_{\mathcal{C},r}}^{{\bar{z}}_1}2\pi \sqrt{\left(x_{\mathcal{C},r}^2-\frac{x_{\mathcal{C},r}^2{\tilde{z}}^2}{z_{\mathcal{C},r}^2}\right)\left(1+\frac{{\tilde{z}}^2{x}_{\mathcal{C},r}^4}{x_{\mathcal{C},r}^2\left(z_{\mathcal{C},r}^4-z_{\mathcal{C},r}^2{\tilde{z}}^2\right)}\right)}d\tilde{z}.$$

$$\frac{A_{\mathcal{C}}}{N-1}={\int }_{{\bar{z}}_{i_p-2}}^{{\bar{z}}_{i_p-1}}2\pi \sqrt{\left(x_{\mathcal{C},r}^2-\frac{x_{\mathcal{C},r}^2{\tilde{z}}^2}{z_{\mathcal{C},r}^2}\right)\left(1+\frac{{\tilde{z}}^2{x}_{\mathcal{C},r}^4}{x_{\mathcal{C},r}^2\left(z_{\mathcal{C},r}^4-z_{\mathcal{C},r}^2{\tilde{z}}^2\right)}\right)}d\tilde{z}.$$

$$p_{id}={\left[x_{id}{y}_{id}{z}_{id}\right]}^T=$$

$$\displaystyle\begin{bmatrix}x_{\mathcal{C},r}\cos\left(\arcsin\left(\frac{z_{id}}{z_{\mathcal{C},r}}\right)\right)\cos\biggl{(}\operatorname{atan2}2\bigl{(}y_{i}(t),x_{i}(t)\bigr{)}\biggr{)}\\ y_{\mathcal{C},r}\cos\left(\arcsin\left(\frac{z_{id}}{z_{\mathcal{C},r}}\right)\right)\sin\biggl{(}\operatorname{atan2}2\bigl{(}y_{i}(t),x_{i}(t)\bigr{)}\biggr{)}\\ \bar{z}_{i_{p}-1},\text{if $z_{i}<\bar{z}_{i_{p}-1}$};\text{or }\bar{z}_{i_{p}-2},\text{if $z_{i}>\bar{z}_{i_{p}-2}$}\end{bmatrix},$$

$$\left[\begin{array}{c}{u}_i^{stm}\\ v_i^{stm}\\ w_i^{stm}\end{array}\right]=U_{max}^{agt}{\mathcal{R}}_1^{-1}\frac{\ln \left(\frac{1}{\left(\gamma +{\mu }_i\right)R-\text{mathcalligra}{r}_i}\left(\Vert n_i\Vert -\text{mathcalligra}{r}_i\right)\right){\hat{n}}_i}{\Vert \ln \left(\frac{1}{\left(\gamma +{\mu }_i\right)R-\text{mathcalligra}{r}_i}\left(\Vert n_i\Vert -\text{mathcalligra}{r}_i\right)\right){\hat{n}}_i\Vert }.$$

$$\begin{array}{rl}& \frac{d}{dt}(S_i({\tilde{q}}_i(t),\tilde{p}))=\frac{\partial S_i}{\partial x_i}{\dot{x}}_i(t)+\frac{\partial S_i}{\partial y_i}{\dot{y}}_i(t)+\frac{\partial S_i}{\partial z_i}{\dot{z}}_i(t)+\frac{\partial S_i}{\partial {\Psi }_i}{\dot{\Psi }}_i(t)+\frac{\partial S_i}{\partial {\Theta }_i}{\dot{\Theta }}_i(t)=\left(\frac{\partial S_i}{\partial x_i}\cos \Theta \cos \Psi +\frac{\partial S_i}{\partial y_i}\cos \Theta \sin \Psi -\frac{\partial S_i}{\partial z_i}\sin \Theta \right)u_i(t)\\ & +(\frac{\partial S_i}{\partial x_i}(\sin \Phi \sin \Theta \cos \Psi -\cos \Phi \sin \Psi )+\frac{\partial S_i}{\partial y_i}(\sin \Phi \sin \Theta \sin \Psi +\cos \Phi \cos \Psi )+\frac{\partial S_i}{\partial z_i}\sin \Phi \cos \Theta )v_i(t)\\ & +(\frac{\partial S_i}{\partial x_i}(\cos \Phi \sin \Theta \cos \Psi +\sin \Phi \sin \Psi )+\frac{\partial S_i}{\partial y_i}(\cos \Phi \sin \Theta \sin \Psi -\sin \Phi \cos \Psi )+\frac{\partial S_i}{\partial z_i}\cos \Phi \cos \Theta )w_i(t)\\ & +\left(\frac{\partial S_i}{\partial {\Psi }_i}\sin \Phi \sec \Theta +\frac{\partial S_i}{\partial {\Theta }_i}\cos \Phi \right)r_i(t)+\left(\frac{\partial S_i}{\partial {\Psi }_i}\cos \Phi \sec \Theta -\frac{\partial S_i}{\partial {\Theta }_i}\sin \Phi \right)s_i(t).\end{array}$$

$$a_{i0}(t,Q(t,\tilde{p}))$$

$$a_{i1}(t,Q(t,\tilde{p}))$$

$$a_{i2}(t,Q(t,\tilde{p}))$$

$$a_{i3}(t,Q(t,\tilde{p}))$$

$$a_{i4}(t,Q(t,\tilde{p}))$$

$$a_{i5}(t,Q(t,\tilde{p}))$$

Full text

††thanks: The authors would like to acknowledge the support of the Automotive Research Center (ARC) in accordance with Cooperative Agreement W56HZV-14-2-0001 U.S. Army TARDEC in Warren, MI and the support by an Early Career Faculty grant from NASA’s Space Technology Research Grants Program. The material in this paper was partially presented at the 56th IEEE Conference on Decision and Control, December 12-15, 2017, Melbourne, Australia and the 2018 American Control Conference, June 27-29, 2018, Milwaukee, Wisconsin. See Section 1.2 for a comparison between the present work and the conference versions.

,

1 Introduction

1.1 Background

The advent of inexpensive autonomous research platforms has spurred recent interest in teams of mobile sensors collaborating on complex surveillance and monitoring tasks. Coverage control problems have been particularly popular due to their numerous applications: e.g., environmental monitoring (Smith et al. 2011), battlefield surveillance (Bokareva et al. 2006), lawn mowing and vacuuming, search and rescue (Murphy et al. 2008), and hull inspections (Choset & Kortenkamp 1999, Hollinger et al. 2013). The latter application is actively supported by NASA whose work on the Mini AERCam paves the way for a future of extravehicular robotic (EVR) free flyers performing independent visual inspections of spacecraft exterior areas of interest (Fredrickson et al. 2004). Free flyer visual inspection is the primary motivating example for our work.

Coverage is often partitioned into three classes of problems: static, dynamic, and persistent. Static coverage problems (e.g., area coverage, k-coverage and point coverage) often explore the optimal arrangement of sensor nodes in a network and the agents tend to immobilize after this arrangement has been achieved (Cortes et al. 2004). Dynamic coverage problems involve the active exploration of a domain. Agents typically must sweep their sensors over all points of a domain until some desired level of coverage has been achieved (Hussein & Stipanović 2007, Liu et al. 2013, Stipanović et al. 2013). Persistent coverage is often similar to dynamic coverage with the addition of information decay within the environment: i.e., agents are required to continually return to areas of interest in order to restore a deteriorating coverage level.

The term ”persistent coverage” appears as early as Hokayem et al. (2007) where agents must cover all points in a 2-D convex polygonal domain every $T^{\star}$ time units. This was accomplished with the design of concentric polygonal trajectories with agents following closed paths in steady state. The work in Song et al. (2013) is similar but also introduces a linear coverage decay rate for specific points of interest. In this paper, as well as Smith et al. (2012), controller design is akin to regulating the velocity along paths generated offline to increase observation time at select points of interest. As the decay rates are known and time invariant, optimal speed control is computed via linear programming.

Palacios-Gasós et. al have published multiple works recently on persistent coverage (Palacios-Gasós et al. 2016b, a, 2017) which build specifically upon Smith et al. (2012). While the earlier work assumed both the existence and knowledge of an optimal path to cover all points of interest, Palacios-Gasós et al. (2016b) uses techniques from discrete optimization and linear programming to iteratively compute this path. The effect is that if the coverage decay rate of a specific point of interest is found to be insufficient to justify the transit time required to service it, then the point may be removed online from the path of the robot. Prior works, i.e., Smith et al. (2012), Song et al. (2013), would have instead driven the robot to quickly pass through the point. Similar techniques are used in Mitchell et al. (2015) which also considers that agents must periodically return to refueling depots.

In Hübel et al. (2008) and Song et al. (2011), the desired coverage level of the domain is maintained with density maps that yield additional observation time at select areas of interest. In Song et al. (2011) the maps are time-invariant while Hübel et al. (2008) considers time-varying density maps that may be designed around moving points of interest (e.g., aerial surveillance targets). However, the latter work only uses density maps in the derivation of control laws and not in the differential equations governing coverage level evolution.

Common themes through all of these persistent coverage works are convex 2-D domains, predictable environments, and simplified sensing and dynamic models for agents. Coverage surfaces embedded in $\mathbb{R}^{3}$ are considered in Cheng et al. (2008); however, this work is closer to that of Hokayem et al. (2007) in that agents also follow preplanned trajectories without considering spatially-dependent coverage decay maps.

Monitoring of stochastic environments is presented in Yu et al. (2015), Pasqualetti et al. (2014), outside of the strict persistent coverage formulation. In Yu et al. (2015), the authors consider that agents must observe events at multiple points of interest and the precise arrival times of events are unknown a priori. Arrival time statistics are used to inform a multi-objective scheduling protocol that results in fixed cyclic servicing policies. In Pasqualetti et al. (2014), the environment contains smart intruders which actively attempt to evade a camera surveillance network. Camera motion is restricted to a single pan axis and thus the the system model is essentially that of a 1-D pursuit evasion problem.

1.2 Contribution

In this paper, we present (i) a formal hybrid control strategy for multi-agent persistent coverage of non-planar convex surfaces embedded in $\mathbb{R}^{3}$ that does not make overly simplifying assumptions with respect to agent dynamic and sensing models, (ii) guarantees on agent interception of stochastic intruders, and (iii) an energy-aware agent deployment and scheduling protocol. To the best of our knowledge, we are the first to present a formal hybrid approach to persistent coverage.

Although many of the cited works use density functions to encode points of interest, the difference in our approach is both subtle and powerful. The density function in Hübel et al. (2008) evolves subject to the motion of an intruder and informs the control laws; however, it has no effect upon the dynamics of the coverage level. Thus, only the intruder’s current location has any influence on the motion of the agents, and the time-history of the intruder’s trajectory is forgotten. This necessitates that agents must travel faster than targets in order to cover points associated with peaks in the density function before they vanish. Works that do include a density function in the coverage level evolution, e.g Palacios-Gasós et al. (2016b, a, 2017), tend to have fixed decay rates that cannot respond or adapt to a changing environment. In contrast, our algorithm utilizes a time-varying density function, which is estimated online via extended Kalman Filter, to directly encode coverage decay over the surface around the predicted impact points of intruders. This encodes a memory effect which drives some agents to follow coverage gradients towards areas that have previously been or will soon be impacted by intruders.

Our agents operate with finite resources and are required to periodically return to a refueling station while observing stochastic events at locations and times that are not known a priori. This approach is different from related works, such as Yu et al. (2015) and Mitchell et al. (2015), where the locations of events are fixed and the authors are concerned with optimal servicing routes between these known stations.

This hybrid system is a successor to our previous works in 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 the approach was unable to provide any formal guarantees on intruder interception without additional operating modes. In Bentz & Panagou (2018), we derived these additional modes to present a hybrid approach to persistent coverage. Agents were now scheduled to intercept intruders and followed path-length optimal trajectories. However, formal guarantees on intruder interception were still limited to cases in which no collision-avoidance deadlocks had occurred. Furthermore, agents did not make effective use of local coverage mode as they would often travel to the predicted impact points of intruders and then remain stationary until the moment of impact thus contributing to a rising coverage error.

This work extends the interception guarantee of Bentz & Panagou (2018) to an arbitrary number of collision avoidance maneuvers and presents an entirely new method of collision avoidance over the prior works. It also reformulates numerous guard conditions within the automaton to allow agents to explore actively around the predicted impact points of intruders. Furthermore, this contribution revises our sensing function definition utilized in the prior works which suffered from a singularity at the sensing cone vertex.

This paper is organized as follows: Section 2 describes the agents sensing and kinematic models and provides an overview of our hybrid control strategy, Section 2.2 presents the trajectory estimator for particle intruders and defines our coverage decay rate map, Sections 3-6 describe each hybrid mode in detail, Section 7 verifies the algorithm in simulation, Section 8 summarizes our contributions and Section 9 is an appendix containing the formal definition of our hybrid automaton.

2 Problem Formulation

2.1 Agent Modeling

Consider a network of spherical autonomous agents indexed $i\in\{1,...,N\}$, of radius $\mathcalligra{r}_{i}$, whose motion is subject to 3-D rigid body kinematics (Beard 2008):

[TABLE]

[TABLE]

where $p_{i}=\left[x_{i}\;\;y_{i}\;\;z_{i}\right]^{T}$ is the position vector and $\Omega_{i}=\left[\Phi_{i}\;\;\Theta_{i}\;\;\Psi_{i}\right]^{T}$ is the vector of 3-2-1 Euler angles taken with respect to a global Cartesian coordinate frame $\mathcal{G}$ with origin $\mathcal{O}$. The linear velocities $\left[u_{i}\;\;v_{i}\;\;w_{i}\right]^{T}$ and angular velocities $\left[q_{i}\;\;r_{i}\;\;s_{i}\right]^{T}$ are both presented in the body fixed frame $\mathcal{B}_{i}$ with origin $p_{i}$. The state vector of agent $i$ is defined as $\tilde{q}_{i}=[p_{i}^{T}\;\;\Omega_{i}^{T}]^{T}$. In the sequel, the rotation matrices of (1) and (2) shall be denoted $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ respectively. The agents travel within a stationary domain, $\mathcal{D}\subset\mathbb{R}^{3}$. Their task is to survey a two-dimensional manifold, $\mathcal{C}\subset\mathcal{D}$, known as our surface of interest. For the purpose of this work we assume that the surface is an ellipsoid of revolution; however, it should be noted that the coverage laws, as well as the collision avoidance strategy, can be easily adapted for any convex surface. The ellipsoid has semi-major axis $x_{\mathcal{C},r}$ and semi-minor axis $z_{\mathcal{C},r}$ aligned with the global coordinate axes $\hat{x}_{\mathcal{G}}$ and $\hat{z}_{\mathcal{G}}$ respectively with center at $\mathcal{O}$. The circumflex (i.e., hat) symbols denote unit vectors.

Each agent, $i$, is equipped with a forward facing sensor whose spherical sector footprint shall be referred to as $\mathcal{S}_{i}$. This model, though intended to be generic, is similar to conical camera models presented in other works on dynamic coverage (see Xie & Zhang (2013)). Our model differs in terms of its heterogeneity, i.e. $\mathcal{S}_{i}$ provides anisotropic sensing data that degrade in quality towards the periphery of the footprint and changes with respect to distance from the sensor. Degradation over distance is not monotonically decreasing but instead contains a peak located near the vertex of $\mathcal{S}_{i}$ as in Hexsel et al. (2013). This is motivated by the fact that the probability of event detection by a camera decreases when either very far from or very close to the lens. Anisotropic sensing is encoded through the definition of the sensing constraint functions for each agent $i$:

[TABLE]

for $\beta_{i}=\min\{1,\eta_{i}\left(\left(\tilde{x}-x_{i}\right)^{2}+\left(\tilde{y}-y_{i}\right)^{2}-\left(\tilde{z}-z_{i}\right)^{2}\right)\}$ with real constant $\eta_{i}>>1$. $R$ is the sensing range, $\tilde{p}_{i}=\left[\tilde{x}\;\;\tilde{y}\;\;\tilde{z}\right]^{T}$ is the position of a point within $\mathcal{S}_{i}$ with respect to $\mathcal{G}$, $\alpha_{i}$ is the angle between the periphery and centerline of the spherical sector (the $\hat{x}_{\mathcal{B}_{i}}$ axis), and $\phi_{i}$ is the angle between $r_{\tilde{p}_{i}/p_{i}}=\tilde{p}_{i}-p_{i}$ (resolved in $\mathcal{G}$ by construction) and the $\hat{x}_{\mathcal{B}_{i}}$ axis given as the inverse cosine of the dot product of $\hat{r}_{\tilde{p}_{i}/p_{i}}$ and $\hat{x}_{\mathcal{B}_{i}}$ resolved in $\mathcal{G}$: $\phi_{i}=\arccos\left(\hat{r}_{\tilde{p}_{i}/p_{i}}\cdot\hat{x}_{\mathcal{B}_{i}}|_{\mathcal{G}}\right)$. Note that: $\hat{r}_{\tilde{p}_{i}/p_{i}}=\frac{1}{\sqrt{(\tilde{x}-x_{i})^{2}+(\tilde{y}-y_{i})^{2}+(\tilde{z}-z_{i})^{2}}}\left[(\tilde{x}-x_{i})\;\;(\tilde{y}-y_{i})\;\;(\tilde{z}-z_{i})\right]^{T}$, and $\hat{x}_{\mathcal{B}_{i}}|_{\mathcal{G}}$ is determined by multiplying $\mathcal{R}_{1}$ by $\left[1\;\;0\;\;0\right]^{T}$: $\hat{x}_{\mathcal{B}_{i}}|_{\mathcal{G}}=\begin{bmatrix}\cos\Psi_{i}\cos\Theta_{i}\;\;\sin\Psi_{i}\cos\Theta_{i}\;\;-\sin\Theta_{i}\end{bmatrix}^{T}$. Agent $i$ is thus capable of detecting objects that lie within an angle of $2\alpha_{i}>0$ about the $\hat{x}_{\mathcal{B}_{i}}$ axis and a range of $R>0$. The model for agent $i$ is depicted in Fig. 1.

Let us denote $\max\{0,c_{ki}\}=C_{ki}$. One can define the sensing function that represents the quality of information available at each point over the sensing domain as:

[TABLE]

where $\bar{C}_{i}$ is the set of zero elements in $C_{ki}$. $S_{i}(\tilde{q}_{i},\tilde{p})$ takes a value of zero outside of $\mathcal{S}_{i}$. Note that $S_{i}(\tilde{q}_{i},\tilde{p})$ is defined over all of $\mathcal{D}$ and thus has static bounds. $S_{i}(\tilde{q}_{i},\tilde{p})$ is continuous in $\tilde{p}$ while taking a value of zero along $\partial\mathcal{S}_{i}$. In verifying this continuity, it is important to note that $S_{i}(\tilde{q}_{i},\tilde{p})$ approaches zero from within $\mathcal{S}_{i}$ in the limit that either $card\left(\bar{C}_{i}\right)=2$ or $r_{\tilde{p}_{i}/p_{i}}=0$ are satisfied. The former condition may be verified by taking a limit of the first piecewise definition of (4) as $C_{1i}$ and $C_{2i}$ tend to zero. The latter condition results from our definitions of $\beta_{i}$ and $\eta_{i}$ which encode that $S_{i}(\tilde{q}_{i},\tilde{p})$ drops off rapidly when in very close proximity to the vertex of $\mathcal{S}_{i}$. Increasing $\eta_{i}>>1$ has the effect of shifting the sensing drop off point closer to $p_{i}$. Define the coverage level provided by agent $i$ at time $t$ as:

[TABLE]

where $C$ is defined as: C\left(\tilde{p}\right)=\begin{cases}1,\hfill&\text{ \forall\tilde{p}\in\mathcal{C};}\\ \hfill 0,\hfill&\text{\forall\tilde{p}\notin\mathcal{C},}\end{cases} and encodes that the accumulation of sensing information only occurs along our surface of interest, $\mathcal{C}$.

As the agents cover $\mathcal{C}$, a set of $N_{p}$ high-speed particle intruders denoted $k\in\{1,...,N_{p}\}$, each of which travels in an arbitrary direction at constant velocity, pass through the domain. The particles are assumed to be uncontrolled and cannot deviate from their initial trajectories. No assumptions are made with respect to the source of the particles or whether they are intelligently generated. Each particle shall have an associated map decay term, $\Lambda_{k}\left(\tau,\tilde{p}\right)$, which is defined later in Section 2.3. We may now define the global coverage level:

[TABLE]

In this work, coverage refers to the accumulation of sensing data over time. Points, $\tilde{p}$, are said to be sufficiently covered when $Q\left(t,\tilde{p}\right)\geq C^{\star}$. The goal is to derive a hybrid control strategy which persistently sweeps $\mathcal{S}_{i}$ across $\mathcal{C}$ while emphasizing surveillance within some bound of the predicted impact points of particles $k\in\{1,...,N_{p}\}$ on $\mathcal{C}$. More specifically, we establish theoretical guarantees on the worst case path length from any agent to any arbitrary impact point thus guaranteeing interception for prescribed bounds on intruder speed, detection range, and agent velocity. This must be done while avoiding collisions. Let us define collision and interception.

Definition 1**.**

Agent $i$ is said to have intercepted particle $k$ if $i$ is within a $\varepsilon_{1}$ bound of the estimated impact point of $k$ for a finite interval of time leading up to the impact. Agent $i$ shall spend this duration of time sweeping the area in local coverage thus gathering information. **

Note that intruders are unaffected by agents and shall always impact the surface and then disappear. This does not damage the agent which is free to resume other tasks upon conclusion of interception.

Definition 2**.**

Agent $i$ avoids collision so long as $\lVert p_{i}(t)-p_{j}(t)\rVert>\mathcalligra{r}_{i}+\mathcalligra{r}_{j},\;\forall t\geq 0,\;\forall j\neq i\in\{1,...,N\}$ and $\lVert n_{i}\rVert>\mathcalligra{r}_{i}$ where the vector $n_{i}$ has direction normal to $\mathcal{C}$ and length equal to the Euclidean distance of its intersection point on $\mathcal{C}$ to $p_{i}$. **

Agents operate with finite power resources and are required to return every $T^{\star}$ time units to a fueling station denoted $\mathcal{F}$. Thus, a scheduling protocol is derived whereby agents periodically deploy from $\mathcal{F}$ to cover within assigned partitions of $\mathcal{C}$. These partitions are bounded by latitude lines and are sorted by geodesic distance from $\mathcal{F}$ with agents deploying to the partition furthest from $\mathcal{F}$ and then transferring between adjacent partitions every $\frac{T^{\star}}{N}$ time units as their power resource dwindles requiring a return to $\mathcal{F}$ within $T^{\star}$ time units after deployment. This partitioning scheme also has the benefit of ensuring that the network of agents is well distributed across $\mathcal{C}$ with agents nominally assigned to intercept intruders with predicted impact points within their own partition.

Agent $i$ is capable of localizing itself in $\mathcal{G}$ and detecting whether there exists $j$ such that $\lVert p_{i}(t)-p_{j}(t)\rVert\leq R$. Furthermore, agents $i$ and $j$ can communicate their deployment times to one another. A centralized network is required to publish the current coverage level $Q(t,\tilde{p})$111In practice, it is not necessary to publish $Q(t,\tilde{p}),\forall\tilde{p}\in\mathcal{D}$ to every agent. For agent $i$ to compute its local coverage control signal, it need only values for $Q(t,\tilde{p})$ within a closed ball of radius $R$ due to the fact that $S_{i}(\tilde{q}_{i},\tilde{p})=0,\forall\tilde{p}\notin\bar{B}_{R}p_{i}(t)$. This substantially reduces the communication overhead. to all agents and to estimate the trajectories of intruders using an omnidirectional range sensor whose measurements are fed through an extended Kalman filter. Computation of $Q(t,\tilde{p})$ is contingent upon continuous transmission of agent state $\tilde{q}_{i}$ to the centralized network. The centralized network assigns each intruder to an unassigned agent at closest latitude to the predicted impact point. It also transmits detection time as well as estimated location and time of impact to the agent.

2.2 Intruder Modeling

We assume that the omnidirectional range sensor (e.g., LiDAR) is co-located with $\mathcal{O}$ and provides measurements of each particle’s position in spherical coordinates. We also assume that particle detection and state estimate initialization occur while the distance of the particle from $\mathcal{O}$ is greater than or equal to $R_{det}+x_{\mathcal{C},r}$ where $R_{det}$ is a lower bound on distance from detection to impact. We define the model for the motion of particle $k$:

[TABLE]

where $\tilde{q}_{k}=\left[x_{k},\,y_{k},\,z_{k},\,\dot{x}_{k},\,\dot{y}_{k},\,\dot{z}_{k}\right]^{T}$ and $\tilde{z}_{k}=\left[\rho_{k},\,\theta_{k},\,\psi_{k}\right]^{T}$ are the Cartesian state and spherical coordinate measurement vectors of particle $k$ resolved in $\mathcal{G}$. We assume that particle speed is upper bounded such that $\sqrt{\dot{x}_{k}^{2}+\dot{y}_{k}^{2}+\dot{z}_{k}^{2}}\leq U_{max}^{int}$. $\rho_{k}$, $\theta_{k}$, and $\psi_{k}$ are the range, azimuthal angle, and polar angle of $k$ respectively. In the sequel, the matrix in (8) shall be denoted $\tilde{h}\left(x_{k},\,y_{k},\,z_{k}\right)$. Assume that the measurement noise, $\epsilon$, is zero-mean Gaussian and has covariance $\mathbf{R}=diag\left(\sigma_{\rho}^{2},\sigma_{\theta}^{2},\sigma_{\psi}^{2}\right)$. This system models high-speed particles incident upon a surface with negligible drag (e.g., micrometeoroids impacting a spacecraft hull); thus, it is reasonable to omit the process noise. The state and covariance estimates, $\hat{\tilde{q}}$ and $\mathbf{P}_{k}$, are computed with a continuous-time extended Kalman filter which is initialized upon particle $k$’s detection at time $t_{dk}$.

2.3 Information Decay

At any time $t$, we define our decay rate map for particle $k$ in terms of its predicted position and covariance evolution over a horizon $T_{H,k}(t)$. As the particles are assumed to travel at fixed velocities222The guarantee of intruder interception presented in this work can be extended to intruders with time-varying velocities that are bounded by $U^{int}_{max}$. However, it is still required that intruders follow straight line trajectories such that the network may estimate fixed impact points., the predicted values for Cartesian position $\tilde{p}^{\prime}_{k}\left(t+\tau\right)$ and associated covariance $\mathbf{P}_{k}\left(t+\tau\right)$ are defined as $\tilde{p}^{\prime}_{k}\left(t+\tau\right)=G\left(\tau\right)\hat{\tilde{q}}_{k}\left(t\right)$, and $\mathbf{P^{\prime}}_{k}\left(t+\tau\right)=G\left(\tau\right)\mathbf{P}_{k}\left(t\right)G\left(\tau\right)^{T}$ respectively where $G\left(\tau\right)=\left[\mathbb{I}_{3\times 3}\;\;\tau\mathbb{I}_{3\times 3}\right]$ and $\hat{\tilde{q}}_{k}\left(t\right)$ is our current estimate for $\tilde{q}_{k}\left(t\right)$. We define the decay rate map associated with particle $k$ as the integral of our predicted normal distribution $\mathcal{N}\left(\tilde{p}^{\prime}_{k}\left(t+\tau\right),\mathbf{P^{\prime}}_{k}\left(t+\tau\right)\right)$ through horizon $T_{H,k}$:

[TABLE]

where $\lambda_{k}>0$ is a tuning parameter for the decay rate. We recommend choosing $\lambda_{k}<1$ as this typically allows for the rate of coverage to exceed the decay rate over points intersecting $\mathcal{S}_{i}$. For $t<t_{dk}$, define $\Lambda_{k}\left(t,\tilde{p}\right)=0,\,\forall\tilde{p}\in\mathcal{D}$. Our formulation for (9) essentially takes a normal distribution for the position of particle $k$ at time $t$ and cumulatively propagates it forward in time up to our horizon $T_{H,k}(t)$. The horizon is lower-bounded by an estimate of the remaining time until impact of particle $k$ on $\mathcal{C}$. This may be computed using $\tilde{q}_{k}\left(t\right)$ along with the surface geometry. With this design, $Q(t,\tilde{p})$ decays along the predicted trajectory of $k$ with tapering omnidirectional decay rates spreading out from the predicted path. This design lends itself naturally to our local coverage formulation, which is gradient following in nature, in that the agents may follow these tapered decay paths towards the predicted impact points on our surface of interest.

2.4 Algorithmic Overview

The coverage strategy for agent $i$ is represented by the hybrid automaton in Fig. 2. Rigorous definitions of all entities of the automaton, including the guard conditions and reset maps, are included in the appendix. Note that each agent operates in accordance with its own automaton and thus an arbitrary number of agents may be in any operating mode at any given time. Before proceeding, we provide a brief overview of each mode.

• •

Local Coverage: This mode governs the active exploration of our surface of interest $\mathcal{C}$. When active, the agent continuously seeks to orient and translate $\mathcal{S}_{i}$ such that $\mathcal{S}_{i}$ intersects portions of $\mathcal{C}$ with a lower coverage level. This is conceptually similar to following the gradient of the coverage error. An agent currently assigned to an intruder may operate in Local Coverage while within an $\varepsilon_{1}$ bound of the the intruder’s predicted impact point. Any agent not currently assigned to an intruder shall operate in Local Coverage assuming that it is within its assigned latitude partition. Operation in Local Coverage is always concurrent with agent assignment to the lowest concentric surface (see Surface Transfer below) and transition to Local Coverage can occur from any mode aside from Return to Base Mode.

• •

Particle Intercept: In this mode, an intruder is assigned to agent $i$ and $i$ is guided along its assigned surface to the predicted impact point of the intruder. After intruder assignment occurs, the agent will nominally remain in Particle Intercept Mode until after the intruder impacts $\mathcal{C}$; however, the agent may temporarily leave the mode before impact to avoid collision through Surface Transfer Mode or to explore in Local Coverage within a $\varepsilon_{1}$ bound of the the intruder’s predicted impact point.

• •

Partition Transfer: This mode is defined for agents that are not currently assigned to an intruder and its purpose is to ensure that the agents are spatially distributed across the entire surface area of $\mathcal{C}$. Activation of this mode will guide agent $i$ along a longitudinally-oriented geodesic trajectory until its position satisfies a set of latitude constraints, i.e., agents travel to the southernmost latitude partition upon deployment and transition through progressively northern partitions as their fuel is depleted. Transition to this mode can occur from any other mode. The partitioning scheme is shown in Fig. 6.

• •

Surface Transfer: This mode’s primary purpose is to ensure that agents avoid collision with two distinct cases resulting in its activation. In the first case, two or more agents have violated a safe-proximity condition. The mode removes select agents from the deadlock by guiding them along vectors normal to $\mathcal{C}$ to a higher-altitude ellipsoidal surface concentric with $\mathcal{C}$. An agent trajectory is then temporarily confined to this newly assigned surface until it reaches the surface projection of its destination. The second case occurs under the condition that the agent has arrived at the projection of its destination on a higher-altitude surface. The mode is activated to return the agent to the innermost surface. Transition to this mode can occur from any other mode aside from Return to Base as agents in the latter mode always take priority in a deadlock. The surface transfer geometry is illustrated in Fig. 4.

• •

Return to Base: The final mode is activated when the time since an agent’s deployment has surpassed some threshold. It guides the power-critical agent along the optimal trajectory to the refueling station. After charging, the agent is redeployed. Agent deployments 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 task of intercepting the assigned intruder or transferring surfaces before transitioning to Return to Base Mode. If an agent is designated as power-critical while in Partition Transfer Mode it shall immediately abandon its task and transition to Local Coverage which shall result in instantaneous transition to Return to Base Mode. Theoretical guarantees on successful return to base with respect to agent power lifespan in accordance with our automaton is presented in Theorem 2 of Section 5.

3 Local Coverage Mode

Local coverage constitutes the first of five hybrid modes in our automaton. This mode is gradient following in nature and commands agent $i$ to always seek to orient and translate $\mathcal{S}_{i}$ such that the volume of uncovered space intersecting $\mathcal{S}_{i}$ is increased. In this way, it emphasizes active exploration of the domain by agents that are not currently assigned to either monitor intruders or relocate within the domain. The control laws are designed such that agent motion in local coverage shall tend to reduce the rate of growth of the global coverage error. Define the global coverage error with respect to $C^{\star}$ as:

[TABLE]

where $h(w)=(\max\{0,w\})^{3}$ with first derivative $h^{\prime}=\frac{dh}{dw}=3(\max\{0,w\})^{2}$ and second derivative $h^{\prime\prime}=\frac{d^{2}h}{dw^{2}}=6(\max\{0,w\})$. Our local coverage control laws are derived via differentiation of (10). This is included in the Appendix in the interest of space. The result is the selection of the following control strategy:

[TABLE]

where:

[TABLE]

$\rho_{l,i}$ is a collision avoidance term with respect to the surface of interest. It takes a value of zero when agent $i$’s normalized distance from $\mathcal{C}$ is $\gamma R$ for $\gamma\in(0,1]$ and is logarithmically repulsive and attractive from the surface when the distance is decreased or increased respectively. Small values for $\gamma$ tend to direct the agent to travel closer to the surface. This coincides with a smaller cross section of $\mathcal{S}_{i}$ intersecting the surface but is also typically associated with a higher quality of sensing. A larger choice for $\gamma$ will direct the agent to fly at a higher altitude with respect to the surface and thus the area covered by $\mathcal{S}_{i}$ will tend to be broader with a decreased quality of sensing. $\rho_{a,i}$, for $\xi<<1$, encodes that the agents should tend to align $\hat{x}_{\mathcal{B}_{i}}$ with $-\hat{n}_{i}$ if the coverage terms associated with $r_{i}$ and $s_{i}$ have become sufficiently small. The physical meaning of $\rho_{a,i}$ is to direct $\mathcal{S}_{i}$ back onto $\mathcal{C}$ if it has reached a configuration in which it no longer intersects $\mathcal{C}$. See Fig. 3 for further details.

$\bar{r_{i}}$ and $\bar{s_{i}}$ are saturation limits for the coverage angular velocity inputs to the system. $k_{u}$, $k_{v}$ and $k_{w}$ are tuning gains which are chosen to satisfy $\sqrt{k_{u}^{2}+k_{v}^{2}+k_{w}^{2}}\leq U_{max}^{agt}$. As $\rho_{\ell,i}$ is normal to the surface, it can be shown that $U_{max}^{agt}$ is an upper bound to agent velocity tangential to $\mathcal{C}$.

4 Particle Intercept Mode

Assuming that particle $k$ is embedded within the surface upon impact, its position shall intersect $\mathcal{C}$ at most one time. We define particle $k$’s estimated impact time as $t_{ck}=\min{\left(t\in\mathbb{R}^{+}\mid\frac{\left(\hat{x}_{k}+\dot{\hat{x}}_{k}t\right)^{2}}{x_{\mathcal{C},r}^{2}}+\frac{\left(\hat{y}_{k}+\dot{\hat{y}}_{k}t\right)^{2}}{x_{\mathcal{C},r}^{2}}+\frac{\left(\hat{z}_{k}+\dot{\hat{z}}_{k}t\right)^{2}}{z_{\mathcal{C},r}^{2}}=1\right)}$, with estimated impact point $\tilde{p}^{\prime}_{k}\left(t_{ck}\right)=G\left(t_{ck}-t\right)\hat{\tilde{q}}_{k}\left(t\right)$. Upon detection, particle $k$ is assigned to a free agent $i$ with the minimum distance from the estimated point of impact along the $\hat{z}_{\mathcal{G}}$ direction. We define a new index, $i_{k}$, as the index of the agent assigned to intercept particle $k$ at destination $p_{id}=\tilde{p}^{\prime}_{k}\left(t_{ck}\right)$: $i_{k}=\operatorname*{\arg\!\min}_{i\in\{1,...,N\}\mid i_{p}\neq 1,f_{i}\neq 1}\lVert\bar{z}_{k}^{\prime}(t_{ck})-z_{i}(t_{dk})\rVert$. Note that $\bar{z}_{k}^{\prime}(t_{ck})$ is the $z$ component of $\bar{p}_{k}^{\prime}(t_{ck})$ and $f_{i}\in\{0,1\}$ is a particle assignment flag for agent $i$ defined as [math] when the agent is free (i.e., not currently assigned a particle). $i_{p}\in\{1,...,N\}$, the power index of agent $i$, shall be fully described in Section 5; however, it should be noted that the definition of $i_{k}$ implies that there are at most $N-1$ agents available for particle interception at any given time.

As with our local coverage strategy, it is assumed that agents shall maintain a distance $\gamma R$ normal to $\mathcal{C}$ in the nominal case that they are not maneuvering to avoid collision. We define an ellipsoid of revolution, $\mathcal{C}_{0}$, which is concentric with $\mathcal{C}$ and has the property that each semi-principal axis is $\gamma R$ longer than its associated counterpart in $\mathcal{C}$, i.e., $x_{\mathcal{C}_{0},r}=x_{\mathcal{C},r}+\gamma R$, and $z_{\mathcal{C}_{0},r}=z_{\mathcal{C},r}+\gamma R$. The nominal trajectories of $i$ are attractive to $\mathcal{C}_{0}$. Agents maneuvering to avoid collision shall transfer to additional concentric ellipsoidal surfaces each separated by a distance of $R$. These surfaces are denoted $\mathcal{C}_{1},\mathcal{C}_{2},...,\mathcal{C}_{N-1}$ with associated semi-principal axes $x_{\mathcal{C}_{1},r}=x_{\mathcal{C},r}+\left(\gamma+1\right)R$ and $z_{\mathcal{C}_{1},r}=z_{\mathcal{C},r}+\left(\gamma+1\right)R$, $x_{\mathcal{C}_{2},r}=x_{\mathcal{C},r}+\left(\gamma+2\right)R$ and $z_{\mathcal{C}_{2},r}=z_{\mathcal{C},r}+\left(\gamma+2\right)R$, etc. Surface assignment and transfer scheduling in collision avoidance mode is described in full detail in Section 6 and the geometry is illustrated in Fig. 4.

When agent $i$ has been assigned to intercept particle $k$, $f_{i}$ is set to $1$ and it is said to have transitioned into particle intercept mode. In this mode, agent $i$ shall nominally follow the optimal trajectory along $\mathcal{C}_{0}$ to within a $\varepsilon_{1}$ bound of the projection of point $\tilde{p}_{k}^{\prime}(t_{ck})$ onto $\mathcal{C}_{0}$ (denoted $\textit{proj}_{\mathit{\mathcal{C}_{0}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}$). The agent shall then transition to local coverage to actively explore within this $\varepsilon_{1}$ bound until $t>t_{ck}$ at which time $f_{i}$ is set to [math]. If local coverage guides the agent out of the $\varepsilon_{1}$ bound, particle intercept mode will again guide the agent back inside the bound. The optimal trajectory is referred to as a geodesic and its computation may be executed in an iterative manner. Specifically, we use Vincenty’s formulae as presented in Vincenty (1975a). For cases involving nearly antipodal points in which the standard inverse method does not converge, we use Vincenty’s supplemental algorithm presented in Vincenty (1975b).

As an input, Vincenty’s algorithm requires an ellipsoid of revolution along with two points, current and desired position, on that surface. The algorithm returns a heading angle measured clockwise from North. This heading angle shall be referred to as $\chi_{i}$. We now define the heading unit vector $\hat{\nu}_{i}$ which lies in a plane tangent to the surface at $p_{i}$. It may be computed by rotating the North-pointing vector at $p_{i}$ clockwise by an angle of $\chi_{i}$ within the tangent plane. For our implementation of Vincenty’s algorithm, we input the following: $\mathcal{C}_{\mu_{i}}$ for surface assignment index $\mu_{i}\in\{0,...,N-1\}$, $p_{i}$, and $\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}$. The position controller used to guide agent $i$ to $\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}$ is composed of two terms: one which commands velocity tangential to $\mathcal{C}_{\mu_{i}}$ along $\hat{\nu}_{i}$ and one logarithmic term which commands velocity normal to $\mathcal{C}_{\mu_{i}}$ in order to constrain the geodesic trajectory of $i$ to $\mathcal{C}_{\mu_{i}}$. The particle intercept mode position control law is:

[TABLE]

As agent $i$ travels along the geodesic, it is desirable that it should point $\mathcal{S}_{i}$ towards $\mathcal{C}$. Therefore, the orientation controller for particle intercept mode is similar to that of Section 3:

[TABLE]

which is essentially a proportional controller that tends to align $\hat{x}_{\mathcal{B}_{i}}$ with $-\hat{n}_{i}$. As (14) commands the vehicle to follow the optimal length path along $\mathcal{C}_{\mu_{i}}$ to $\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}$, we can establish a few guarantees on system performance. To simplify notation, define:

[TABLE]

and $g_{\mathcal{C}_{0}}$ is defined similarly in terms of the semi-principal axes of $\mathcal{C}_{0}$.

Lemma 1**.**

Let us assume that agent $i$ has been assigned to particle $k$ with $f_{i}:=1$. Given an arbitrary agent position $p_{i}(t_{dk})$ and an arbitrary predicted impact point for the intruder $\tilde{p}_{k}^{\prime}(t_{ck})$, there exists an upper bound to the maximum path length until interception: $\mathcal{P}_{max}\leq\pi x_{\mathcal{C}_{N-1},r}+\frac{\pi}{2}\left(x_{\mathcal{C}_{N-1},r}+z_{\mathcal{C}_{N-1},r}\right)g_{\mathcal{C}_{N-1}}+2\left(N-1\right)R$.* {pf} At the moment that $f_{i}:=1$ we have that agent $i$ transitions to Particle Intercept Mode. Under the condition that the agent has not yet come within proximity of the predicted impact point, i.e., $\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}\rVert>\varepsilon_{1}$, we have that only $G\left(\zeta_{i2},\zeta_{i4}\right)$ and $G\left(\zeta_{i4},\zeta_{i2}\right)$ are defined (see Appendix). These two transitions occur sequentially for each deadlock event that agent $i$ encounters as it approaches $\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}$.*

In any given deadlock arrangement, one agent remains on its current surface without ascending to a higher one. This implies that $\mu_{i}=1$ for at most $N-1$ agents as the remaining agent would be unable to encounter a deadlock on $\mathcal{C}_{0}$. Furthermore, this implies that $\mu_{i}=2$ for at most $N-2$ agents etc. until we have $\mu_{i}=N$ for at most zero agents. The worst cast surface assignment that can be incurred during sequential cycles of $\left(\left(\zeta_{i2},\zeta_{i4}\right),\left(\zeta_{i4},{i2}\right)\right)$ would be $\mu_{i}=N-1$.**

As the geodesic path length between any two points projected onto surface $\mathcal{C}_{\mu_{i}}$ shall always be less than the geodesic path length between the same two points projected onto surface $\mathcal{C}_{\mu_{i}+1}$, we may bound the geodesic portion of the trajectory by one that is constrained entirely to $\mathcal{C}_{N-1}$. We denote this term $\mathcal{P}_{geo}$. As any two points on $\mathcal{C}_{\mu_{i}}$ can be connected by a path of constant latitude $\mathcal{P}_{lat}$ followed by a path of constant longitude $\mathcal{P}_{long}$, we have that:**

[TABLE]

For two generic points on $\mathcal{C}_{N-1}$, we have that:**

[TABLE]

where the bound on $\mathcal{P}_{lat}$ is half of the circumference of the ellipsoid of revolution $\mathcal{C}_{N-1}$ about its equator and the bound on $\mathcal{P}_{long}$ is half of the perimeter of the revolved ellipse. The infinite series expression term, denoted $g_{\mathcal{C}_{N-1}}$ in (19), is first presented in Ivory (1798). The remaining portion of the path length is simply the straight line segments connecting $\mathcal{C}_{0}$ to $\mathcal{C}_{N-1}$ and back again. This length is precisely $2\left(N-1\right)R$. Thus, $\mathcal{P}_{max}=\mathcal{P}_{geo}+2\left(N-1\right)R$ as illustrated in Fig. 5. Invoking (18) and (19) gives us $\mathcal{P}_{max}\leq\pi x_{\mathcal{C}_{N-1},r}+\frac{\pi}{2}\left(x_{\mathcal{C}_{N-1},r}+z_{\mathcal{C}_{N-1},r}\right)g_{\mathcal{C}_{N-1}}+2\left(N-1\right)R$. This concludes the proof.**

Theorem 1**.**

Assuming that the bounds on intruder velocity and range from detection to impact satisfy $\frac{R_{det}}{U_{max}^{int}}>\frac{\mathcal{P}_{max}}{U_{max}^{agt}}$, where $\mathcal{P}_{max}$ may be bounded via Lemma 1, agent $i$ shall reach $\textit{proj}_{\mathit{\mathcal{C}_{0}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}$ before $t_{ck}$.* {pf} Given the fact that agents in the particle intercept and surface transfer modes travel at speed $U_{max}^{agt}$, we have that the time $t_{ik}$ required to travel from $p_{i}\left(t_{dk}\right)$ to $\textit{proj}_{\mathit{\mathcal{C}_{0}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}$ must satisfy: $t_{ik}\leq\frac{\mathcal{P}_{max}}{U_{max}^{agt}}$. Given that $t_{ck}-t_{dk}\geq\frac{R_{det}}{U_{max}^{int}}$, agent $i$ reaching $\textit{proj}_{\mathit{\mathcal{C}_{0}}}\mathit{\tilde{p}_{k}^{\prime}(t_{ck})}$ before $t_{ck}$ implies that $\frac{R_{det}}{U_{max}^{int}}>t_{ik}$. This is guaranteed if $\frac{R_{det}}{U_{max}^{int}}>\frac{\mathcal{P}_{max}}{U_{max}^{agt}}$ This concludes the proof.*

Remark 1**.**

At any given time, there are at most $N-1$ agents available to intercept particles. Thus, satisfaction of Theorem 1 implies that the network is capable of intercepting all particles so long as a maximum of $N-1$ particle impacts occur in any moving time window of $\frac{\pi x_{\mathcal{C}_{N-1},r}+\frac{\pi}{2}\left(x_{\mathcal{C}_{N-1},r}+z_{\mathcal{C}_{N-1},r}\right)g_{\mathcal{C}_{N-1}}+2\left(N-1\right)R}{U_{max}^{agt}}$.**

5 Energy-aware Scheduling Protocol

5.1 Domain Partitioning

As this is a persistent coverage protocol, which operates indefinitely, it is necessary to establish an agent deployment and scheduling protocol that realistically considers the agents’ finite power and/or propulsive resources. Our strategy is to periodically deploy agents from a fueling station $\mathcal{F}$ which we assume to be located at the North pole of $\mathcal{C}_{0}$, i.e., at the point $\left[0\;\;0\;\;z_{\mathcal{C}_{0},r}\right]^{T}$. Define $T^{\star}$ as the power lifespan of each agent in the network. Given $T^{\star}$ and $N$, we define our deployment and scheduling protocol such that one agent is deployed from $\mathcal{F}$ every $\frac{T^{\star}}{N}$ seconds. The initial deployment is that of agent $i=1$ at $t=0$ seconds with agent $i=2$ following at $t=\frac{T^{\star}}{N}$. This continues indefinitely with the second deployment of agent $i=1$ occurring at $t=T^{\star}$ seconds.

In order to adequately distribute agents across $\mathcal{C}$, it is desirable to partition the domain and assign agents to monitor separate regions. Specifically, partitioning the domain by latitude, rather than longitude, ensures that agents are poised to intercept particles without the need for frequent crossings of the equator which tend to be associated with larger values of $\mathcal{P}_{geo}$ on an oblate spheroid.

Define the power index of agent $i$ as $i_{p}(t)=1+\mod{\left(i-2-\left\lfloor\frac{tN}{T^{\star}}\right\rfloor,N\right)}$ where the first argument of our modulo operation is the dividend and the second argument is the divisor. The lower-bracketed delimiters represent the floored division operation. Upon deployment from $\mathcal{F}$, agent $i$ has power index $i_{p}=N$ and this index is reduced by one every $\frac{T^{\star}}{N}$ seconds until $i_{p}=1$, i.e., agent $i$ is the power critical agent. Note that no two agents may share the same power index as a result of our periodic deployment and scheduling protocol.

Latitude partitions are characterized by a static upper bound in $\hat{z}_{\mathcal{G}}$ denoted $\bar{z}_{i_{p}-2}$ and a static lower bound $\bar{z}_{i_{p}-1}$. Rather than dynamically sizing partitions relative to agent power resources, we divide partitions such that $N-1$ agents are assigned equal surface areas of $\mathcal{C}$ to explore. This choice maximizes the coverage of any individual partition as the allocation of a larger partition to a recently deployed agent would result in less effective coverage of that partition. Agents are sorted by their remaining power and transfer between partitions that are progressively closer to $\mathcal{F}$ as their power resource expires. Define the surface area of our ellipsoid of revolution $\mathcal{C}$ as:

[TABLE]

The agent with $i_{p}=2$ is assigned to monitor the partition characterized by upper bound at north pole of $\mathcal{C}$, i.e., $\bar{z}_{0}=z_{\mathcal{C},r}$. The lower bound $\bar{z}_{1}$ may be computed by dividing (20) by $\left(N-1\right)$, equating with the integral of ellipse cross sectional circumferences parametrized by $\tilde{z}$, and then numerically solving for $\bar{z}_{1}$:

[TABLE]

One may then iteratively solve for the remaining bounds for increasing values of $i_{p}$ up to $i_{p}=N-1$:

[TABLE]

The final computation of (22) for $i_{p}=N$ is not necessary as $\bar{z}_{N-1}$ is the south pole of $\mathcal{C}$, i.e., $\bar{z}_{N-1}=-z_{\mathcal{C},r}$, although this may be shown through numerical computation as well. Our partitioning strategy for the case where $N=4$ is presented in Fig. 6.

Note that no partition has been assigned to the agent for which $i_{p}=1$. This is the power critical agent and it shall have flag $f_{i}:=1$ at the instant $i_{p}:=1$. The power critical agent cannot be assigned a new particle to intercept after $i_{p}:=1$ as this opens the possibility that particle assignment could occur near the end of the $\frac{T^{\star}}{N}$ time window during which time the agent with $i_{p}=1$ should be transitioning back to $\mathcal{F}$ to exchange its power source. The power critical agent will instead spend the majority of this time window in local coverage mode assisting the other agents in gathering information. It can only be tasked with intercepting a particle if this assignment had occurred previously when $i_{p}=2$. In this scenario, the agent should be capable of intercepting particle $k$ and then transitioning back to $\mathcal{F}$ so long as a bound is established on the length of our deployment scheduling window $\frac{T^{\star}}{N}$.

Theorem 2**.**

If agent power lifespan $T^{\star}$ satisfies $\frac{T^{\star}}{N}\geq t_{ck}-t_{dk}+\frac{\pi}{2U_{max}^{agt}}\left(x_{\mathcal{C}_{0},r}+z_{\mathcal{C}_{0},r}\right)g_{\mathcal{C}_{0}},\;\forall k$ then the agent with $i_{p}=1$ shall always be capable of reaching $\mathcal{F}$ within $\frac{T^{\star}}{N}$ of the time at which $i_{p}:=1$.**

{pf}

Consider the worst-case scenario in which the agent with $i_{p}=2$ is assigned to intercept particle $k$ at the instant before $i_{p}:=1$. It’s remaining flight time is currently $\frac{T^{\star}}{N}$. The time required to intercept the particle is $t_{ck}-t_{dk}$, after which our control strategy dictates that the agent will follow a geodesic trajectory to $\mathcal{F}$. As $\mathcal{F}$ lies at the north pole of $\mathcal{C}_{0}$, this will be a trajectory of constant longitude which may be upper bounded by a length half the perimeter of our revolved ellipsoid: $\frac{\pi}{2}\left(x_{\mathcal{C}_{0},r}+z_{\mathcal{C}_{0},r}\right)g_{\mathcal{C}_{0}}$ by definition. As the agent is controlled by (14) with a North-pointing $\hat{\nu}_{i}$, it will proceed along this geodesic at speed $U_{max}^{agt}$. Thus the time required to complete this trajectory is $\frac{\pi}{2U_{max}^{agt}}\left(x_{\mathcal{C}_{0},r}+z_{\mathcal{C}_{0},r}\right)g_{\mathcal{C}_{0}}$ and we may bound our deployment window: $\frac{T^{\star}}{N}\geq t_{ck}-t_{dk}+\frac{\pi}{2U_{max}^{agt}}\left(x_{\mathcal{C}_{0},r}+z_{\mathcal{C}_{0},r}\right)g_{\mathcal{C}_{0}},\;\forall k$. This concludes the proof.**

Remark 2**.**

The appropriate design method for this surveillance system is to first ensure that the time from detection to impact of any arbitrary particle, $t_{ck}-t_{dk}$, as governed by the omnidirectional range sensor satisfies Theorem 1. One must subsequently ensure that power lifespan $T^{\star}$, for all agents, satisfies Theorem 2.**

5.2 Partition Transfer and Return to Base

If an agent with $i_{p}\in\{2,...,N\}$ lies outside of its prescribed partition, and we have $i_{f}=0$, then the agent shall enter partition transfer mode. This mode uses the same geodesic position and orientation controllers (14) and (15) with the destination position set to the point:

[TABLE]

i.e., the closest point along the agent’s current longitude which lies on the boundary of its assigned partition.

The return to base mode is similar to partition transfer mode but is defined for the agent with $i_{p}=1$. This mode is activated when the time since agent $i$’s last deployment from $\mathcal{F}$, denoted $t_{i\mathcal{F}}\geq T^{\star}-\frac{\pi}{2U_{max}^{agt}}\left(x_{\mathcal{C}_{0},r}+z_{\mathcal{C}_{0},r}\right)g_{\mathcal{C}_{0}}$ as established in Theorem 2. The control strategy is the same as partition transfer mode with the desired position set to $\mathcal{F}$. Control laws for partition transfer mode and return to base shall be denoted with superscripts $ptm$ and $rtb$ respectively.

6 Surface Transfer Mode

The primary purpose of surface transfer mode is to encode collision avoidance and it can be transitioned into from any other mode aside from the return to base mode. This mode is triggered for agent $i$, assigned to surface $\mathcal{C}_{\mu_{i}}$, when we have the condition that $\lVert p_{i}-p_{j}\rVert\leq R$ for $i\neq j$. Denote $\tilde{j}=i\cup j$ as the set of agents satisfying this condition. Agents in $\tilde{j}$ are ranked by $t_{\tilde{j}\mathcal{F}}$. One agent, denoted $i_{pr}$, whose value for $t_{\tilde{j}\mathcal{F}}$ is highest, i.e., $i_{pr}=\operatorname*{\arg\!\max}_{\tilde{j}}\left(t_{\tilde{j}\mathcal{F}}\right)$ is permitted to proceed. The remaining agents increment their surface assignment indices, $\mu_{i}$, by one and transition to surface transfer mode. This mode controls the agents to follow $\hat{n}_{i}$ until they have transferred to their newly assigned concentric surface at a height $R$ above the previous. Note that in general, convexity of surface $\mathcal{C}$ is required to ensure that intersections of $n_{i}$ and $n_{j},\,\forall i\neq j,$ lie within the interior space that is bounded by the surface. The surface transfer position control strategy is:

[TABLE]

As the agents ascend to a point at which $R$ does not intersect $\mathcal{C}$, sensing information is not gathered in avoidance mode and thus the avoidance orientation control is simply $\left[q_{i}^{stm}\;\;r_{i}^{stm}\;\;s_{i}^{stm}\right]^{T}=\left[0\;\;0\;\;0\right]^{T}$.

Agents are said to have converged upon their newly assigned surface when $\lvert\ln{\frac{\lVert n_{i}\rVert-\mathcalligra{r}_{i}}{\left(\gamma+\mu_{i}\right)R-\mathcalligra{r}_{i}}}\rvert<\varepsilon_{2}$. At this point, each agent shall transition back to its prior mode as described in the following two scenarios.

If agent $i$ had been in either particle intercept or partition transfer mode before the deadlock, it shall resume that mode and continue along a geodesic trajectory on the newly assigned surface until it reaches the projection of its destination. At this point, the condition that $\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert\leq\varepsilon_{1}$ triggers a reset $\mu_{i}:=0$ concurrent with a transition back to surface transfer mode thus allowing the agent to transfer back to $\mathcal{C}_{0}$. The agent then resumes coverage of $\mathcal{C}_{0}$ in its prior mode. For additional details on flag conditions in these transitions, see guards $G\left(\zeta_{i2},\zeta_{i4}\right),G\left(\zeta_{i4},\zeta_{i2}\right),G\left(\zeta_{i3},\zeta_{i4}\right)$, and $G\left(\zeta_{i4},\zeta_{i3}\right)$ of our hybrid automaton as presented in the appendix. 2. 2.

If agent $i$ had been in local coverage mode before the deadlock, it shall then transition back to local coverage mode concurrent with reset $\mu_{i}:=0$. This transition is dependent upon the conditions that $f_{i}=0$ and that the agent is operating within its assigned partition. The agent shall oscillate between local coverage and surface transfer at an altitude of $R$ above $\mathcal{C}_{0}$ until $i_{pr}$ has moved along $\mathcal{C}_{0}$ to resolve the deadlock. At this point, the local coverage controller shall attract agent $i$ back to the surface.

While similar work on multi-agent systems often invoke avoidance barrier functions to encode collision avoidance, such as in Panagou et al. (2016), it may be impossible to bound the time that agents spend avoiding one another in these maneuvers—especially when the algorithm is scaled to many agents. In contrast, our technique results in an explicit bound on path length to an intruder as was proven in Lemma 1. With an additional assumption on the size of agents, we can establish a guarantee on collision avoidance for agents in surface transfer mode.

Theorem 3**.**

For agents $\{i,j\}\in\tilde{j}$, the condition that $\min(R_{\tilde{j}})>2\mathcalligra{r}_{i}+2\mathcalligra{r}_{j}$ implies that $i$ does not collide with $j$.* {pf} Consider the case in which $i\neq i_{pr}$ and $j\neq i_{pr}$. Both agents operate in accordance with (23) and follow trajectories along $\hat{n}_{i}$ and $\hat{n}_{j}$ respectively. Both unit vectors are normal to surface $\mathcal{C}_{\mu_{i}}$, an ellipsoid of revolution, and thus diverge from one another away from $\mathcal{C}_{\mu_{i}}$. Agents $i$ and $j$ shall enter surface transfer mode at an instant when $\lVert p_{i}-p_{j}\rVert\geq\min\left(R_{\tilde{j}}\right)$ and their distance shall tend to increase under (23). Thus $\min(R_{\tilde{j}})>\mathcalligra{r}_{i}+\mathcalligra{r}_{j}$ and subsequently $\min(R_{\tilde{j}})>2\mathcalligra{r}_{i}+2\mathcalligra{r}_{j}$ imply that they avoid collision.*

Consider the case in which $i=i_{pr}$ and thus $j\neq i_{pr}$. In the instant that $j$ transitions to surface transfer mode we have that $\lVert p_{i}-p_{j}\rVert\geq\min\left(R_{\tilde{j}}\right)$. Thus the distance for $i$ to travel until collision is greater than or equal to $\min\left(R_{\tilde{j}}\right)-\mathcalligra{r}_{i}-\mathcalligra{r}_{j}$. This straight line path for $i$ is a conservative estimate as the true path is curved. Collision will be avoided if agent $j$, whose path is normal to the surface, may cover a distance $\mathcalligra{r}_{i}+\mathcalligra{r}_{j}$ before $i$ covers $\min\left(R_{\tilde{j}}\right)-\mathcalligra{r}_{i}-\mathcalligra{r}_{j}$. As $j$ moves at speed $U_{max}^{agt}$ and $i$’s tangential speed is upper bounded by $U_{max}^{agt}$, this condition is satisfied if $\min\left(R_{\tilde{j}}\right)-\mathcalligra{r}_{i}-\mathcalligra{r}_{j}>\mathcalligra{r}_{i}+\mathcalligra{r}_{j}$. This may equivalently be written as $\min\left(R_{\tilde{j}}\right)>2\mathcalligra{r}_{i}+2\mathcalligra{r}_{j}$. These arguments apply to the case in which $j=i_{pr}$ and $i\neq i_{pr}$ as well. This concludes the proof.**

7 Simulations

A simulation was performed in MATLAB to verify the efficacy of the algorithm. Four agents are deployed to monitor the surface of an ellipsoid of revolution, $\mathcal{C}$, whose radius in the $xy$-plane is 80 and whose radius in the $z$-plane is 20. For each agent, $R=10$, $\mathcalligra{r}_{i}=1$, $\alpha_{i}=30^{\circ}$, $k_{u}=1$, $k_{v}=5$, $k_{w}=1$, $k_{r}=0.1$, $k_{s}=0.1$, $\bar{r}_{i}=0.4$, $\bar{s}_{i}=0.4$, $U_{max}^{agt}=6$, and $T^{\star}=792$. Upon initialization of the simulation, $\mathcal{C}$ was set to a fully covered level of $C^{\star}=20$ which would begin decaying upon detection of the first intruder $k\in\{1,...,4\}$ at $t=600$ seconds. The four agents were deployed from $\mathcal{F}$ sequentially at times $t=0$, $t=\frac{T^{\star}}{4}$, $t=\frac{2T^{\star}}{4}$, and $t=\frac{T^{\star}}{4}$ seconds respectively. Upon deployment, each agent was initialized in local coverage mode with $\Phi_{i}=0$, $\Theta_{i}=\frac{\pi}{2}$, and $\Psi_{i}=0$. Intruders traveled in random directions with $U_{max}^{int}=0.7$, though were still constrained to always impact the surface, and were generated every 35 seconds beginning at $t=\frac{3T^{\star}}{4}$ seconds. The detection system had a lower bound on range $R_{det}=80$, decay rate parameter $\lambda_{k}=0.05$, and measurement variances $\sigma_{\rho}^{2}=0.0625$, $\sigma_{\theta}^{2}=0.25\textnormal{ deg}^{2}$, and $\sigma_{\psi}^{2}=0.25\textnormal{ deg}^{2}$ respectively

Agents were able to successfully intercept all particles along their geodesic trajectories while actively avoiding collision over the entire duration of the attack (see Fig. 7 and Fig. 8); however, it should be noted that one avoidance anomaly occurred before the initial intruder was generated during the interval of $t=418-420$. As our sensing range for agents was $R=10$ and our simulation time step size was 1, it is clear that this anomaly occurred due to a selection of $U_{max}^{agt}$ that was too large relative to the time step. In a continuous time implementation, a transition to surface transfer mode would have occurred between $t=418$ and $t=419$ thus preventing collision. Aside from this anomaly, the simulation parameters adequately approximated the continuous time agent kinematics.

The coverage error on $\mathcal{C}$, normalized with respect to the maximum error in which all of $\mathcal{C}$ takes a value of zero for $Q\left(t,\tilde{p}\right)$, as well as the minimum inter-agent distance over time are presented in Fig. 9. The error tends to spike upon particle detections with agents effectively curtailing these spikes as they cover around the vicinity of predicted impact points in local coverage mode. Two particularly large spikes occur at $t=3225$ and $t=5180$ respectively. These anomalies are each associated with particle impacts occurring close to the equator of the ellipsoid where even small values of $\sigma_{\theta}^{2}$ and $\sigma_{\psi}^{2}$ result in an estimated particle trajectory that does not initially intersect $\mathcal{C}$ thus delaying an agent assignment. In both cases, the estimated trajectory did eventually intersect $\mathcal{C}$ with enough time to allow for agent interception. However, this delay in assignment significantly reduced the time the agent spent exploring in the vicinity of the predicted impact point thus contributing to a noticeable rise in the coverage error. One potential solution to this problem would be to prescribe some boundary tolerance to our surface $\mathcal{C}$ thus loosening our definition of an impacting particle for the sake of measurement uncertainty.

Agent $i=1$’s operating modes with respect to time are presented in Fig. 10. It should be noted that the most frequent transition out of particle intercept mode is to local coverage mode. This corresponds to an agent arriving at the estimated impact point of a particle and then surveying the local area up until the moment of impact. As the agent surveys it tends to hit the $\varepsilon_{1}$ proximity boundary to the impact point thus requiring a short operation in particle intercept mode to direct the agent back within the $\varepsilon_{1}$ boundary.

To demonstrate scalability, we have made an additional simulation with 100 agents available online at: https://1drv.ms/f/s!AsiVOlIEkwNEgX2o1eV2hJ_bbaQU.

8 Conclusions

In this paper, we presented a hybrid formulation for the persistent coverage problem in an environment subject to stochastic intruders. This formulation was motivated in part by extravehicular applications of the NASA Mini AERCam. Agents operated with finite power resources and were required to periodically return to a refueling station while patrolling assigned latitude partitions along the surface of an ellipsoid. Formal guarantees were established on the ability of agents to intercept all intruders and the efficacy of the algorithm was demonstrated in simulation. This approach succeeds our previous work in Bentz & Panagou (2017) and Bentz & Panagou (2018) by extending the guarantees on intruder interception to an arbitrary number of collision avoidance maneuvers. It also removes singularities in the sensing function definition and redefines the guard conditions in a manner that supports a more effective use of local coverage around the vicinity of intruder impact points.

9 Appendix

9.1 Additional Derivations for Local Coverage Strategy

Our local coverage control laws are derived via differentiation of (10) of which we seek to reduce the rate of growth. It is a volume integral, so a few mathematical preliminaries are required. Recall the generalized transport theorem (GTT) (Slattery 1999): $\frac{d}{dt}\int\limits_{R_{(s)}}fdV=\int\limits_{R_{(s)}}\frac{\partial f}{\partial t}\,dV+\int\limits_{S_{(s)}}f\mathbf{v_{(s)}}\cdot\mathbf{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$ is integrated, $\mathbf{n}$ is the unit vector normal to the boundary, and $\mathbf{v_{(s)}}$ is the velocity of the boundary. $V$ and $A$ refer to volume and area respectively. Invoking GTT allows for differentiation of (10) with respect to time:

[TABLE]

where $\partial\mathcal{D}$ is the boundary of $\mathcal{D}$. $\mathcal{D}$ is time invariant and thus $\mathbf{v_{(s)}}=0$. (24) reduces to $\dot{E}(t)=\int\limits_{D}\ h^{\prime}(C^{\star}C\left(\tilde{p}\right)-Q(t,\tilde{p}))\left(\frac{-\partial Q(t,\tilde{p})}{\partial t}\right)\,d\tilde{p}$, which expands to:

[TABLE]

$\hat{e}_{i}(t)$ is the rate of change of the coverage error due to the motion of the agents while $\tilde{e}_{k}(t)$ is the rate of change of the coverage error due to a contrived information decay surrounding the predicted impact point of particle $k$ on $\mathcal{C}$. Our strategy is to control the agents’ kinematics, recovered in the derivative of $\hat{e}_{i}(t)$, to decrease (25). Note that we do not presume that our local coverage strategy provides any additional bounds on (10). Nor do we provide guarantees on the rate of growth of this contrived quantity. Curtailing the growth of the coverage error simply imparts the desired effect of active exploration in the vicinity of impact points into our system. Using this strategy, the agents actively seek to increase their rate of coverage by rotating and/or translating $\mathcal{S}_{i}$ to be encompass the most uncovered space in the local vicinity.

Taking the derivative of $\hat{e}_{i}(t)$ with respect to time yields:

[TABLE]

The sensing footprint is independent of $\Phi_{i}$ assuming that the centerline of the spherical sector is aligned with the $\hat{x}_{\mathcal{B}_{i}}$ axis. $\frac{d}{dt}(S_{i}(\tilde{q}_{i}(t),\tilde{p}))$ is expanded in (LABEL:eq:10) and through the definitions in (28-33) one may restate (26) as (34).

If one were to command zero inputs to this system, it becomes clear that $a_{i0}(t,Q(t,\tilde{p}))$ may be physically interpreted as the rate at which the coverage rate is reducing due to information saturation at any particular position and orientation of the sensing footprint, $\mathcal{S}_{i}$. As the footprint remains stationary, there are diminishing returns on the value of newly acquired information. Thus, the additional terms in (34) allow for the coverage rate to be increased by mobilizing the sensor. One strategy is that of (11).

9.2 Formal Hybrid Formulation

To provide a compact notation in this section, define $\acute{f}_{i}=\frac{x_{i}^{2}}{\left(x_{\mathcal{C},r}+\mathcalligra{r}_{i}\right)^{2}}+\frac{y_{i}^{2}}{\left(y_{\mathcal{C},r}+\mathcalligra{r}_{i}\right)^{2}}+\frac{z_{i}^{2}}{\left(z_{\mathcal{C},r}+\mathcalligra{r}_{i}\right)^{2}}$. The coverage strategy for 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}=\{\zeta_{i0},\zeta_{i1},\zeta_{i2},\zeta_{i3},\zeta_{i4}\},$

• •

A set of continuous states: $\tilde{q}_{i}=\{x_{i},y_{i},z_{i},\Phi_{i},\Theta_{i},\Psi_{i}\},$

• •

A vector field: $\newline f(\zeta_{i0},\tilde{q}_{i})=\mathcal{R}\left[u_{i}^{loc}\;v_{i}^{loc}\;w_{i}^{loc}\;0\;r_{i}^{loc}\;s_{i}^{loc}\;\right]^{T},\newline f(\zeta_{i1},\tilde{q}_{i})=\mathcal{R}\left[u_{i}^{rtb}\;v_{i}^{rtb}\;w_{i}^{rtb}\;q_{i}^{rtb}\;r_{i}^{rtb}\;s_{i}^{rtb}\;\right]^{T},\newline f(\zeta_{i2},\tilde{q}_{i})=\mathcal{R}\left[u_{i}^{pim}\;v_{i}^{pim}\;w_{i}^{pim}\;q_{i}^{pim}\;r_{i}^{pim}\;s_{i}^{pim}\;\right]^{T},\newline f(\zeta_{i3},\tilde{q}_{i})=\mathcal{R}\left[u_{i}^{ptm}\;v_{i}^{ptm}\;w_{i}^{ptm}\;q_{i}^{ptm}\;r_{i}^{ptm}\;s_{i}^{ptm}\;\right]^{T},\newline f(\zeta_{i4},\tilde{q}_{i})=\mathcal{R}\left[u_{i}^{stm}\;v_{i}^{stm}\;w_{i}^{stm}\;0\;0\;0\;\right]^{T}\textnormal{where}\;\mathcal{R}=\begin{bmatrix}\mathcal{R}_{1}&0\\ 0&\mathcal{R}_{2}\end{bmatrix},$

• •

A set of initial states: $\{\zeta_{i3}\}\times\{\tilde{q}_{i}\in\mathbb{R}^{6}\mid p_{i}=\mathcal{F}\\ \land\,\Phi_{i}\in\left[-\pi,+\pi\right]\,\land\,\Theta_{i}\in\left[\frac{-\pi}{2},\frac{+\pi}{2}\right]\,\land\,\Psi_{i}\in\left[-\pi,+\pi\right]\},$

• •

A domain: $Dom\left(\zeta_{i0}\right)=\{\tilde{q}_{i}\in\mathbb{R}^{6}\mid\acute{f}_{i}\geq 1\land\\ \left(i_{p}\in\{2,...,N\}\implies\bar{z}_{{i_{p}}-1}\leq z_{i}\leq\bar{z}_{{i_{p}}-2}\right)\},\\ Dom\left(\zeta_{i1}\right)=\{\tilde{q}_{i}\in\mathbb{R}^{6}\mid\acute{f}_{i}\geq 1\},\\ Dom\left(\zeta_{i2}\right)=\{\tilde{q}_{i}\in\mathbb{R}^{6}\mid\acute{f}_{i}\geq 1\},\\ Dom\left(\zeta_{i3}\right)=\{\tilde{q}_{i}\in\mathbb{R}^{6}\mid\acute{f}_{i}\geq 1\land\\ \left(i_{p}\in\{2,...,N\}\implies z_{i}<\bar{z}_{{i_{p}}-1}\lor z_{i}>\bar{z}_{{i_{p}}-2}\right)\},\\ Dom\left(\zeta_{i4}\right)=\{\tilde{q}_{i}\in\mathbb{R}^{6}\mid\acute{f}_{i}\geq 1\},$

• •

A set of edges: $E=\{\left(\zeta_{i0},\zeta_{i1}\right),\left(\zeta_{i0},\zeta_{i2}\right),\left(\zeta_{i0},\zeta_{i3}\right),\\ \left(\zeta_{i0},\zeta_{i4}\right),\left(\zeta_{i1},\zeta_{i3}\right),\left(\zeta_{i2},\zeta_{i0}\right),\left(\zeta_{i2},\zeta_{i1}\right),\left(\zeta_{i2},\zeta_{i3}\right),\\ \left(\zeta_{i2},\zeta_{i4}\right),\left(\zeta_{i3},\zeta_{i0}\right),\left(\zeta_{i3},\zeta_{i2}\right),\left(\zeta_{i3},\zeta_{i4}\right),\left(\zeta_{i4},\zeta_{i0}\right),\\ \left(\zeta_{i4},\zeta_{i2}\right),\left(\zeta_{i4},\zeta_{i3}\right),\},$

• •

A set of guard conditions:\\ G\left(\zeta_{i0},\zeta_{i1}\right)=\{i_{p}=1\land t_{i\mathcal{F}}\geq T^{\star}-\frac{\pi g_{\mathcal{C}_{0}}}{2U_{max}^{agt}}\left(x_{\mathcal{C}_{0},r}+z_{\mathcal{C}_{0},r}\right)\},\\ G\left(\zeta_{i0},\zeta_{i2}\right)=\{\left(\exists k\mid i=i_{k}\right)\land\left(\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert>\varepsilon_{1}\right)\},\\ G\left(\zeta_{i0},\zeta_{i3}\right)=\{i_{p}\neq 1\land\left(z_{i}<\bar{z}_{{i_{p}}-1}\lor z_{i}>\bar{z}_{{i_{p}}-2}\right)\},\\ G\left(\zeta_{i0},\zeta_{i4}\right)=\{\lVert p_{i}-p_{j}\rVert\leq R\land i_{pr}\neq\operatorname*{\arg\!\max}_{\tilde{j}}\left(t_{\tilde{j}\mathcal{F}}\right)\},\\ G\left(\zeta_{i1},\zeta_{i3}\right)=\{\lVert p_{i}-\mathcal{F}\rVert\leq\varepsilon_{1}\land t_{i\mathcal{F}}=T^{\star}\},\\ G\left(\zeta_{i2},\zeta_{i0}\right)=\{\bigl{(}t\geq t_{ck}\land\\ \left.\bigl{(}\left(i_{p}\in\{2,...,N\}\land\bar{z}_{{i_{p}}-1}\leq z_{i}\leq\bar{z}_{{i_{p}}-2}\right)\lor\right.\\ \left.\left.\left(i_{p}=1\land t_{i\mathcal{F}}<T^{\star}-\frac{\pi g_{\mathcal{C}_{0}}}{2U_{max}^{agt}}\left(x_{\mathcal{C}_{0},r}+z_{\mathcal{C}_{0},r}\right)\right)\right)\right)\lor\\ \left(t<t_{ck}\land\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert\leq\varepsilon_{1}\land\mu_{i}=0\right)\},\\ G\left(\zeta_{i2},\zeta_{i1}\right)=\{\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert\leq\varepsilon_{1}\land\mu_{i}=0\land t\geq t_{ck}\land i_{p}=1\land t_{i\mathcal{F}}\geq T^{\star}-\frac{\pi g_{\mathcal{C}_{0}}}{2U_{max}^{agt}}\left(x_{\mathcal{C}_{0},r}+z_{\mathcal{C}_{0},r}\right)\}\\ G\left(\zeta_{i2},\zeta_{i3}\right)=\{\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert\leq\varepsilon_{1}\land\mu_{i}=0\land t\geq t_{ck}\land i_{p}\in\{2,...,N\}\land\left(z_{i}<\bar{z}_{{i_{p}}-1}\lor z_{i}>\bar{z}_{{i_{p}}-2}\right)\},\\ G\left(\zeta_{i2},\zeta_{i4}\right)=G\left(\zeta_{i0},\zeta_{i4}\right)\lor\{\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert\leq\varepsilon_{1}\land\left(\lVert p_{i}-p_{j}\rVert>R,\forall j\lor i_{pr}=\operatorname*{\arg\!\max}_{\tilde{j}}\left(t_{\tilde{j}\mathcal{F}}\right)\right)\land\mu_{i}>0\},\\ G\left(\zeta_{i3},\zeta_{i0}\right)=\{i_{p}=1\lor\left(i_{p}\neq 1\land\bar{z}_{{i_{p}}-1}\leq z_{i}\leq\bar{z}_{{i_{p}}-2}\right)\},\\ G\left(\zeta_{i3},\zeta_{i2}\right)=G\left(\zeta_{i0},\zeta_{i2}\right),\\ G\left(\zeta_{i3},\zeta_{i4}\right)=G\left(\zeta_{i0},\zeta_{i4}\right)\lor\{\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert\leq\varepsilon_{1}\land\left(\lVert p_{i}-p_{j}\rVert>R,\forall j\lor i_{pr}=\operatorname*{\arg\!\max}_{\tilde{j}}\left(t_{\tilde{j}\mathcal{F}}\right)\right)\land\mu_{i}>0\},\\ G\left(\zeta_{i4},\zeta_{i0}\right)=\{f_{i}=0\land\lvert\ln{\left(\frac{\lVert n_{i}\rVert-\mathcalligra{r}_{i}}{\left(\gamma+\mu_{i}\right)R-\mathcalligra{r}_{i}}\right)}\rvert<\varepsilon_{2}\land\left(i_{p}=1\lor\\ \left(i_{p}\in\{2,...,N\}\land\bar{z}_{{i_{p}}-1}\leq z_{i}\leq\bar{z}_{{i_{p}}-2}\right)\right)\},\\ G\left(\zeta_{i4},\zeta_{i2}\right)=\{f_{i}=1\land\lvert\ln{\left(\frac{\lVert n_{i}\rVert-\mathcalligra{r}_{i}}{\left(\gamma+\mu_{i}\right)R-\mathcalligra{r}_{i}}\right)}\rvert<\varepsilon_{2}\},\\ G\left(\zeta_{i4},\zeta_{i3}\right)=\{f_{i}=0\land\lvert\ln{\left(\frac{\lVert n_{i}\rVert-\mathcalligra{r}_{i}}{\left(\gamma+\mu_{i}\right)R-\mathcalligra{r}_{i}}\right)}\rvert<\varepsilon_{2}\land\left(i_{p}\neq 1\land\left(z_{i}<\bar{z}_{{i_{p}}-1}\lor z_{i}>\bar{z}_{{i_{p}}-2}\right)\right)\}.

• •

Additional parameters include a clock set: $C=\{t_{i\mathcal{F}}\},$ a flag: $f_{i}\in\{0,1\}$, an assignment index $\mu_{i}=\{0,...,N-1\}$ and,

• •

A reset map: $R\left(\zeta_{i0},\zeta_{i2},f_{i}\right)=\{1\},R\left(\zeta_{i0},\zeta_{i4},\mu_{i}\right)=\\ \{\mu_{i}+1\},R\left(\zeta_{i1},\zeta_{i3},t_{i\mathcal{F}}\right)=\{0\},R\left(\zeta_{i2},\zeta_{i0},f_{i}\right)=\\ \{0\textnormal{ if }t\geq t_{ck};1\textnormal{ otherwise}\},R\left(\zeta_{i2},\zeta_{i1},f_{i}\right)=\{0\},\\ R\left(\zeta_{i2},\zeta_{i3},f_{i}\right)=\{0\},R\left(\zeta_{i2},\zeta_{i4},\mu_{i}\right)=\{0\textnormal{ if }\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert\\ \leq\varepsilon_{1}\land\left(\lVert p_{i}-p_{j}\rVert>R,\forall j\lor i_{pr}=\operatorname*{\arg\!\max}_{\tilde{j}}\left(t_{\tilde{j}\mathcal{F}}\right)\right)\land\mu_{i}>0;\mu_{i}+1\textnormal{ otherwise}\},R\left(\zeta_{i3},\zeta_{i2},f_{i}\right)=\{1\},\\ R\left(\zeta_{i3},\zeta_{i4},\mu_{i}\right)=\{0\textnormal{ if }\lVert p_{i}-\textit{proj}_{\mathit{\mathcal{C}_{\mu_{i}}}}\mathit{p_{id}}\rVert\leq\varepsilon_{1}\land\\ \left(\lVert p_{i}-p_{j}\rVert>R,\forall j\lor i_{pr}=\operatorname*{\arg\!\max}_{\tilde{j}}\left(t_{\tilde{j}\mathcal{F}}\right)\right)\land\mu_{i}>0;\\ \mu_{i}+1\textnormal{ otherwise}\},R\left(\zeta_{i4},\zeta_{i0},\mu_{i}\right)=\{0\},$ and continuous states do not reset between transitions.