Physics-Based Freely Scalable Continuum Deformation for UAS Traffic Coordination

TL;DR

This paper introduces a physics-inspired continuum deformation method for UAS traffic management, extending 1-D traffic models to 2-D airspace, with boundary control and vehicle clustering for resilient coordination.

Contribution

It presents a novel higher-dimensional traffic coordination framework for UAS, incorporating boundary control, resilience to failures, and vehicle clustering techniques.

Findings

Effective boundary control algorithm enhances resilience.

Model successfully extends to 2-D airspace without predefined paths.

Clustering improves microscopic coordination of UAS.

Abstract

This paper develops a novel physics-inspired traffic coordination approach and applies it to Unmanned Aircraft System (UAS) traffic management. We extend available physics-inspired approaches previously applied to 1-D traffic flow on highways and urban streets to support models of traffic coordination in higher dimension airspace for cases where no predefined paths exist. The paper considers airspace as a finite control volume while UAS coordination, treated as continuum deformation, is controlled at the airspace boundaries. By partitioning airspace into planned and unplanned spaces, the paper models nominal coordination in the planned airspace as the solution of a partial differential equation with spatiotemporal parameters. This paper also improves resilience to vehicle failures with a resilient boundary control algorithm to update the geometry of the planned space when UAS problems threaten safe coordination in existing navigable airspace channels. To support UAS coordination at the microscopic level, we propose clustering vehicles based on vehicle performance limits. UAS clusters, with each UAS treated as a particle of a virtual rigid body, use leader-follower containment to acquire the macroscopic desired trajectory.

Figures

Tables

Table 1. Table I: 2-D ideal fluid flow Potential functions

Pattern	Potential Function
Uniform	${\Phi }_U({𝐫}_{xy},{𝐔}_{\infty ,f}\left(t\right))={𝐔}_{\infty ,f}\cdot {𝐫}_{xy}$
Source	${\Phi }_{\mathrm{So}}({𝐫}_{xy},{𝐫}_{0,i,f})={\Delta }_i\ln \Vert {𝐫}_{xy}-{𝐫}_{0,i,f}\Vert$
Sink	${\Phi }_{\mathrm{Si}}({𝐫}_{xy},{𝐫}_{0,i,f})=-{\Delta }_i\ln \Vert {𝐫}_{xy}-{𝐫}_{0,i,f}\Vert$
Doublet	${\Phi }_D({𝐫}_{xy},{𝐫}_{0,i,f})=\frac{{\Delta }_i\left({𝐫}_{xy}-{𝐫}_{0,i,f}\right)\cdot {𝐧}_f}{\Vert {𝐫}_{xy}-{𝐫}_{0,i,f}\Vert }$

Table 2. Table II: 2-D ideal fluid flow Potential functions

Pattern	Stream Function
Uniform	${\Psi }_U({𝐫}_{xy},{𝐔}_{\infty ,f}\left(t\right))=u_{\infty ,f}\left(y\cos {\theta }_{0,f}-x\sin {\theta }_{0,f}\right)$
Source	${\Psi }_{\mathrm{So}}({𝐫}_{xy},{𝐫}_{0,i,f})={\Delta }_i{\tan }^{-1}\left(\frac{y-y_{0,i,f}}{x-x_{0,i,f}}\right)$
Sink	${\Psi }_{\mathrm{Si}}({𝐫}_{xy},{𝐫}_{0,i,f})=-{\Delta }_i{\tan }^{-1}\left(\frac{y-y_{0,i,f}}{x-x_{0,i,f}}\right)$
Doublet	${\Psi }_D({𝐫}_{xy},{𝐫}_{0,i,f})=\frac{-{\Delta }_i\left[\sin {\theta }_{0,f}\left(x-x_{0,i,f}\right)+\cos {\theta }_{0,f}\left(y-y_{0,i,f}\right)\right]}{{\left(x-x_{0,i,f}\right)}^2+{\left(y-y_{0,i,f}\right)}^2}$

Table 3. Table III: Communication weights w i , j , h 1 subscript 𝑤 𝑖 𝑗 subscript ℎ 1 w_{i,j,h_{1}} , w i , j , h 2 subscript 𝑤 𝑖 𝑗 subscript ℎ 2 w_{i,j,h_{2}} , and w i , j , h 3 subscript 𝑤 𝑖 𝑗 subscript ℎ 3 w_{i,j,h_{3}} ( i ∈ 𝒮 , j ∈ 𝒱 F , i , 𝒩 j , i = { h 1 , h 2 , h 3 } formulae-sequence 𝑖 𝒮 formulae-sequence 𝑗 subscript 𝒱 𝐹 𝑖 subscript 𝒩 𝑗 𝑖 subscript ℎ 1 subscript ℎ 2 subscript ℎ 3 i\in\mathcal{S},\leavevmode\nobreak\ j\in\mathcal{V}_{F,i},\leavevmode\nobreak\ \mathcal{N}_{j,i}=\{h_{1},h_{2},h_{3}\} ).

$j$	$h_1$	$h_2$	$h_3$	$w_{i,j,h_1}$	$w_{i,j,h_2}$	$w_{i,j,h_3}$
$4$	$1$	$7$	$10$	$0.50$	$0.25$	$0.25$
$5$	$2$	$8$	$9$	$0.50$	$0.25$	$0.25$
$6$	$3$	$9$	$10$	$0.50$	$0.25$	$0.25$
$7$	$4$	$8$	$10$	$0.40$	$0.30$	$0.40$
$8$	$5$	$7$	$9$	$0.29$	$0.35$	$0.36$
$9$	$5$	$6$	$8$	$0.31$	$0.4$	$0.29$
$10$	$4$	$6$	$7$	$0.45$	$0.25$	$0.30$

Equations

$$\mathcal{C}=\left(\partial \mathcal{P}\bigcup \mathcal{P}\right)\setminus \mathcal{U}$$

$$\begin{array}{rl}\mathbf{r}=& x{\hat{\mathbf{e}}}_x+y{\hat{\mathbf{e}}}_y+z{\hat{\mathbf{e}}}_z\\ \mathbf{V}=& v_x{\hat{\mathbf{e}}}_x+v_y{\hat{\mathbf{e}}}_y+v_z{\hat{\mathbf{e}}}_z\end{array},$$

$$\mathbf{r}\in \mathcal{C},\text{leavevmode}t\ge 0,\mathbf{V}\left(\mathbf{r},t\right)=K\left(\mathbf{r},t\right)▽{\Phi }_a\left(\mathbf{r},t\right),$$

$$\mathbf{r}\in \mathcal{C},\text{leavevmode}t\ge 0,\frac{d{\Phi }_a}{dt}=▽\cdot \left(▽{\Phi }_a\right)=▽\cdot \left(K\left(\mathbf{r},t\right)\mathbf{V}\left(\mathbf{r},t\right)\right),$$

$$▽\left(▽\Phi \right)=0$$

$$\begin{array}{rl}& {\int }_{\mathcal{C}\setminus \left(\partial \mathcal{U}\bigcup \partial \mathcal{P}\right)}\left(▽\cdot \left(▽{\Phi }_a\right)\right)\mathrm{d}V={\int }_{\partial \mathcal{P}}\left(▽{\Phi }_a\right)\cdot {\hat{\mathbf{n}}}_p\mathrm{d}A\\ -& \sum _{j=1}^{n_u}\left({\int }_{\partial {\mathcal{U}}_j}\text{cancelto}0\left(▽{\Phi }_a\right)\cdot {\hat{\mathbf{n}}}_{u_j}\mathrm{d}A\right)={\int }_{\partial \mathcal{P}}\left(▽{\Phi }_a\right)\cdot {\hat{\mathbf{n}}}_p\mathrm{d}A.\end{array}$$

$$\mathbf{r}\in \partial \mathcal{P},q\left(\mathbf{r},t\right)=▽\Phi \cdot {\hat{\mathbf{n}}}_p.$$

$$\forall t,{\int }_{\partial \mathcal{P}}q\left(\mathbf{r},t\right)\mathrm{d}A={\int }_{\partial \mathcal{P}}\left(▽\Phi \right)\cdot {\hat{\mathbf{n}}}_p\mathrm{d}A=0.$$

$$\forall t,\text{leavevmode}\mathbf{r}\in \mathcal{C},▽\cdot \left(▽\Phi \right)=▽\cdot \left(▽\cdot (K(\mathbf{r},t)\mathbf{V}(\mathbf{r},t)\right))=0.$$

$$-\mathbf{L}\left(t\right)\mathbf{\Phi }\left(t\right)+\mathbf{q}\left(t\right)=0,$$

$$\hat{\mathbf{L}}=\mathbf{SL}{\mathbf{S}}^T=\left[\begin{array}{ccc}{\hat{\mathbf{L}}}_{cb,cb}& {\hat{\mathbf{L}}}_{cb,ci}& {\hat{\mathbf{L}}}_{cb,u}\\ {\hat{\mathbf{L}}}_{cu,cb}& {\hat{\mathbf{L}}}_{cu,ci}& {\hat{\mathbf{L}}}_{cu,u}\\ {\hat{\mathbf{L}}}_{u,cb}& {\hat{\mathbf{L}}}_{u,ci}& {\hat{\mathbf{L}}}_{u,u}\end{array}\right],$$

$${\mathbf{S}}_{hq}=\{\begin{array}{ll}1& h=1,\cdots ,m_{cb}\wedge \left(j_h\in {\mathcal{V}}_{cb}\right)\wedge j_h=q\\ 1& h=m_{cb}+1,\cdots ,m_{cb}+m_{ci}\wedge \left(l_h\in {\mathcal{V}}_{ci}\right)\wedge l_h=q\\ 1& h=m_{cb}+m_{ci}+1,\cdots ,m\wedge \left(i_h\in {\mathcal{V}}_u\right)\wedge i_h=q\\ 0& \mathrm{else}.\end{array}$$

$$-\hat{\mathbf{L}}(t)\hat{\Phi }(t)=\hat{\mathbf{q}}(t).$$

$$\left[\begin{array}[]{c|c}-\hat{\mathbf{L}}_{cb,cb}&\mathbf{0}\\ \hline\cr\mathbf{B}_{c}&\mathbf{A_{c}}\end{array}\right]=-\left[\begin{array}[]{c|ccc}\hat{\mathbf{L}}_{cb,cb}&\mathbf{0}&\mathbf{0}\\ \hline\cr\hat{\mathbf{L}}_{ci,cb}&\hat{\mathbf{L}}_{ci,ci}&\hat{\mathbf{L}}_{ci,u}\\ \mathbf{0}&\mathbf{0}&\mathbf{I}_{m_{u}}\\ \end{array}\right]$$

$${\mathbf{A}}_c{\mathbf{\Phi }}_I+{\mathbf{B}}_c{\hat{\mathbf{\Phi }}}_{cb}=0,$$

$$\frac{\partial {\Phi }_f}{\partial x}=\frac{\partial {\Psi }_f}{\partial y}=K_f(x,y,t)u_f(x,y,t)$$

$$\frac{\partial {\Phi }_f}{\partial y}=-\frac{\partial {\Psi }_f}{\partial x}=K_f(x,y,t)v_f(x,y,t),$$

$${\mathbf{r}}_f={\mathbf{r}}_{xy}+z_f(x,y){\hat{\mathbf{e}}}_z$$

$$\begin{array}{rl}{\Phi }_f=& {\Phi }_{\mathrm{U}}\left({\mathbf{r}}_{xy},{\mathbf{U}}_{\infty ,f}(t)\right)+{\Xi }_f\sum _{i=1}^{n_u}[{\gamma }_{i,f}{\Phi }_{\mathrm{D}}\left({\mathbf{r}}_{xy},{\mathbf{r}}_{0,i,f}\left(t\right)\right)\\ +& \left({\Phi }_{\mathrm{Si}}\left({\mathbf{r}}_{xy},{\mathbf{r}}_{0,i,f}\left(t\right)\right)+{\Phi }_{\mathrm{So}}\left(\mathbf{r},{\mathbf{r}}_{0,i,f}\left(t\right)\right)\right)]\end{array}$$

$$\begin{array}{rl}{\Psi }_f=& {\Psi }_{\mathrm{U}}\left({\mathbf{r}}_{x,y},{\mathbf{U}}_{\infty ,f}(t)\right)+{\Xi }_f\sum _{i=1}^{n_u}[{\gamma }_{i,f}{\Psi }_{\mathrm{D}}\left({\mathbf{r}}_{xy},{\mathbf{r}}_{0,i,f}\right)\\ +& \left(1-{\gamma }_{i,f}\right)\left({\Psi }_{\mathrm{Si}}\left({\mathbf{r}}_{xy},{\mathbf{r}}_{0,i,f}\right)+{\Psi }_{\mathrm{So}}\left({\mathbf{r}}_{xy},{\mathbf{r}}_{0,i,f}\right)\right)],\end{array}$$

$$\partial\mathcal{U}=\left\{\mathbf{r}\in\mathcal{P}\big{|}\Psi_{f}\left(\mathbf{r}_{xy},t\right)=0,\leavevmode\nobreak\ f=1,\cdots,n_{s}\right\}.$$

$${\varphi }_f=u_{\infty ,f}\left(1+\frac{{\Delta }_f}{x^2+y^2}\right)x$$

$${\psi }_f=u_{\infty ,f}\left(1-\frac{{\Delta }_f}{x^2+y^2}\right)y$$

$$f=1,\cdots ,n_s,\text{leavevmode}\mathbf{r}\in \mathcal{C},\text{leavevmode}t\ge 0,\frac{d{\Psi }_{f,a}}{dt}=▽\cdot \left(▽{\Psi }_{f,a}\right)$$

$$f=1,\cdots ,n_s,\text{leavevmode}\mathbf{r}\in \partial \mathcal{P},\text{leavevmode}t\ge 0,\frac{\partial {\Psi }_{f,a}}{\partial t}=U_p(\mathbf{r},t)$$

$$f=1,\cdots ,n_s,\text{leavevmode}\mathbf{r}\in \partial \mathcal{U},\text{leavevmode}t\ge 0,\frac{\partial {\Psi }_{f,a}}{\partial t}=U_u(\mathbf{r},t).$$

$$f=1,\cdots ,n_s,\text{leavevmode}\mathbf{r}\in \mathcal{C},\text{leavevmode}t\ge 0,▽\cdot \left(▽{\Psi }_f\right)=0$$

$$f=1,\cdots ,n_s,\text{leavevmode}\mathbf{r}\in \partial \mathcal{P},\text{leavevmode}t\ge 0,{\Psi }_f=\mathrm{given},$$

$$f=1,\cdots ,n_s,\text{leavevmode}\mathbf{r}\in \partial \mathcal{U},\text{leavevmode}t\ge 0,{\Psi }_f=0.$$

$$f=1,\cdots ,n_s,\text{leavevmode}\mathbf{r}\in \mathcal{C},\text{leavevmode}t\ge 0,\frac{d{E}_f}{dt}=▽\cdot \left(▽E_f\right).$$

$$f=1,\cdots ,n_s,\text{leavevmode}\mathbf{r}\in \partial \mathcal{P},\text{leavevmode}t\ge 0,U_p(\mathbf{r},t)=\frac{{\Psi }_{f,a}}{dt}-k_p\left(▽E_f\cdot {\mathbf{n}}_p\right)$$

$$\begin{array}{rl}& f=1,\cdots ,n_s,\text{leavevmode}j=1,\cdots ,n_u,\text{leavevmode}\mathbf{r}\in \partial {\mathcal{U}}_j,\text{leavevmode}t\ge 0,\\ & U_{u_j}(\mathbf{r},t)=-k_{u_j}\left(▽E_f\cdot {\mathbf{n}}_{u_j}\right)\end{array}$$

$${\dot{\mathbf{\Phi }}}_I^a={\mathbf{A}}_c{\mathbf{\Phi }}_I^a+{\mathbf{B}}_c{\hat{\mathbf{\Phi }}}_{cb}^a$$

Full text

I Introduction

Multi-agent coordination has been widely investigated over the past two decades. Centralized and decentralized coordination approaches typically assume the total number of agents is fixed in a given motion space. This assumption must be relaxed in emerging multi-agent coordination applications such as Unmanned Aerial System (UAS) Traffic Management (UTM) in which the capacity of a finite airspace volume or sector must vary over time, particularly to achieve the high-density operational throughput levels anticipated for urban regions. This paper introduces the notion of free scalability for UTM in which the number of vehicles and the specific coordination solution are governed by controlling vehicle inflow and outflow at the boundaries of a motion finite space. This macroscopic modeling approach is inspired by the enroute airspace sector used for commercial aviation today adapted to an autonomous UTM paradigm of the future.

This paper models air traffic management as multi-agent coordination, a topic widely studied in the literature. Formation and cooperative control of a multi-agent system can enhance resilience to failure [1], improve efficiency, and reduce mission cost [2]. Applications include surveillance [3], air traffic management [4], formation flight [5], and connected vehicle control [6]. Virtual structures [7, 8], consensus [9, 10, 11], containment control [12, 13, 14], partial-differential equation (PDE) based approaches [15, 16, 17], continuum deformation [18, 19], and graph rigidity methods [20, 21] are available multi-agent coordination approaches. While the virtual structure is a centralized approach, the others are typically implemented in a decentralized fashion. In addition, containment control, PDE-based, and continuum deformation methods apply leader-follower coordination in which leaders move independently and followers acquire the desired coordination through strictly local communication.

The above coordination approaches assume the total number of agents is constant over time within a designated volume or sector. This is a limiting assumption for a traffic system in which vehicles can dynamically leave or enter a traffic control sector. This limitation can be addressed by modeling traffic flow as a spatiotemporal system governed by a partial differential equation [22, 23, 24, 25]. Ref. [23] uses mass conversation to model traffic flow near intersections. Ref. [26] models swarm coordination as an optimal mass transport problem governed by the continuity PDE. A coupled PDE-ODE approach is developed in [24] to analyze the underlying influence of large and slow vehicles on large-scale traffic for highways. Multi-lane traffic flow is modeled by a first-order mass conservation PDE in Ref. [27]. In Refs. [28, 29, 25] optimal control is integrated with physics-based models to optimize highway traffic flow.

Available physically-inspired approaches mostly study 1-D traffic flow in highways and urban areas with mapped roads, pathways, and highways. In contrast, this paper develops an approach to low-altitude UAS coordination (LAUC) in a free airspace with no predefined paths. Autonomous collision avoidance is one of the main LAUC challenges. Each UAS must have sense and avoid (SAA) [30, 31, 32, 33, 34] capabilities to avoid collision. SAA using laser/lidar [33] and machine vision [35] technologies is widely studied in the literature. Ref. [30] develops a vision-based SAA algorithm to avoid collision through trajectory prediction. SAA using Boolean decision logic (BDL) is proposed in [32]. SAA using geofence boundary violation detection and turn-back have been proposed in Refs. [36, 37, 38] to avoid collision by satisfying keep-in airspace permission constraints.

I-A Contribution and Outline

This paper develops a continuum deformation framework for traffic coordination (CDFTC) in a finite motion space and applies CDFTC to LAUC. CDFTC applies principles of continuum mechanics [39] to achieve freely scalable UAS coordination in a finite airspace volume. CDFTC integrates the following two main layers: (i) Macroscopic coordination using Eulerian continuum mechanics and (ii) Microscopic coordination of clustered vehicles treated as particles of a virtual rigid body. By dividing the motion space into planned and unplanned spaces, the macroscopic layer assigns UAS nominal coordination over the planned space in a centralized fashion. To this end, macroscopic coordination is defined as the solution of a parabolic spatiotemporal partial differential equation (PDE). In the microscopic layer, CDFTC applied to LAUC allows a UAS to move either individually or as part of a group. By treating agents as particles of a continuum (or deformable body), a leader-follower containment approach is developed to coordinate clusters of UAS in a decentralized fashion.

This paper applies the principles of continuum mechanics to manage UAS coordination. To our knowledge, this is the first paper to propose a physics-based macroscopic or meta-level coordination law that assures each UAS cluster can safely transit local planned airspace channels via decentralized containment control. By using an Eulerian description of continuum mechanics, time and space allocated to each UAS can be efficiently managed. This paper also supports UAS heterogeneity in both macroscopic and microscopic levels whereas most of the cooperative control literature assumes all vehicles have the same performance. By shaping the spatiotemporal parameters of the macroscopic coordination PDE, we show how UAS with different nominal operation speeds can efficiently and safely occupy a finite airspace volume or sector. At the microscopic level, the paper treats UAS as particles of [math]-D (single UAS), $1$-D, $2$-D, and $3$-D virtual rigid bodies as appropriate to acquire the desired macroscopic coordination in a decentralized fashion. For this purpose, we apply leader-follower containment theory.

This paper is organized as follows. A problem statement is presented in Section II. Functionality of the CDFTC macroscopic and microscopic layers is described in Sections III and IV, respectively. Simulation results presented in Section V are followed by concluding remarks in Section VI.

II Problem Statement

The paper applies principles of continuum mechanics to safely coordinate a large-scale MAS in a finite airspace. A finite airspace sector requiring LAUC is defined by open set $\mathcal{P}$ with boundary $\partial\mathcal{P}$. Assume that $\mathcal{P}$ contains a finite number of obstacles defined by closed set unplanned set $\mathcal{U}$. Then,

[TABLE]

defines the planned subset of the airspace. Note that $\mathcal{U}=\{\mathcal{U}_{1},\cdots,\mathcal{U}_{n_{u}}\}$ defines obstacle zones $n_{u}$ excluded from $\mathcal{P}$. A schematic of planned-unplanned airspace partitioning is shown in Fig. 1. As shown in Fig. 1, $\mathbf{n}_{u_{i}}(\mathbf{r},t)$ ($i=1,\cdots,n_{u}$) denotes the outward unit normal vector on the unplanned airspace boundary while $\mathbf{n}_{p}(\mathbf{r},t)$ is the outward unit normal vector on airspace boundary $\partial\mathcal{P}$. Furthermore, $\mathbf{t}_{u_{i}}(\mathbf{r},t)$ ($i=1,\cdots,n_{u}$) is the normal tangent vector at unplanned boundary $\partial\mathcal{U}_{i}$ ($i=1,\cdots,n_{u}$).

For macroscopic coordination, we use an Eulerian description of continuum mechanics to define nominal coordination over the planned airspace and manage time interval and space allocation to UAS across the planned airspace in a computationally-effective manner. In practice, this maximizes airspace usability by scaling services provided to UAS while they exist in the planned airspace. This paper models UAS macroscopic coordination as the solution of a PDE with spatiotemporal parameters. To enhance resilience of the macroscopic coordination against failure in the planned airspace, we offer a novel boundary control approach to update geometry of the planned airspace in a safe manner.

For microscopic coordination, UAS clusters are formed based on identifying compatible nominal operational speeds of all aerial vehicles operating in the region. UAS clusters are then treated as a particles of virtual rigid bodies with 3-D rotations and translations in the planned airspace. A leader-follower containment approach is used by each cluster to track the nominal coordination assigned at the macroscopic level in a decentralized fashion.

III Macroscopic Coordination

For macroscopic coordination, an Eulerian description of continuum mechanics is used to define a spatio-temporal UAS cluster velocity field over $\mathcal{C}$, represented as $\mathbf{V}=\mathbf{V}\left(\mathbf{r},t\right)$ where $\mathbf{r}$ and $t$ are three-dimensional position and time, respectively. Defining a Cartesian coordinate system with unit bases $\hat{\mathbf{e}}_{x}$, $\hat{\mathbf{e}}_{y}$, and $\hat{\mathbf{e}}_{z}$, position and velocity are expressed as follows:

[TABLE]

where $x$, $y$, and $z$ are position components and $v_{x}$, $v_{y}$, and $v_{z}$ are velocity components.

Governing Equation: Defining an (actual) potential function $\Phi_{a}=\Phi_{a}(\mathbf{r},t)$ where $\mathbf{r}\in\mathcal{P}$ denotes position in the planned space and $t$ denotes time. Nominal velocity of the field is given by

[TABLE]

where $K(\mathbf{r},t):(\mathcal{C},\mathbb{R}_{+})\rightarrow\mathbb{R}_{+}$ is the cost and $\Phi_{a}(\mathbf{r},t)$ is governed by the following partial differential equation (PDE):

[TABLE]

where $\bigtriangledown=\left[\frac{\partial}{\partial x}\leavevmode\nobreak\ \frac{\partial}{\partial y}\leavevmode\nobreak\ \frac{\partial}{\partial z}\right]^{T}$ is the gradient vector defined with respect to the Cartesian coordinate system.

III-A Reference Macroscopic Coordination

Let $\Phi=\Phi(\mathbf{r},t)$ represent a reference potential field. A reference macroscopic coordination is the solution of

[TABLE]

This paper applies the divergence theorem to relate the volume integral $\bigtriangledown\cdot\left(\bigtriangledown\Phi\right)$ over $\mathcal{C}\setminus\left(\partial\mathcal{U}\bigcup\partial\mathcal{P}\right)$ and the surface integral of $\bigtriangledown\Phi$ over the borders $\partial\mathcal{U}$ and $\partial\mathcal{P}$ by

[TABLE]

In (4), $\partial\mathcal{U}_{j}$ defines a boundary surface enclosing the unplanned airspace $\mathcal{U}_{j}$, $\partial\mathcal{U}=\{\partial\mathcal{U}_{1},\cdots,\partial\mathcal{U}_{n_{u}}\}$, $\mathbf{n}_{u_{j}}$ is the outward unit vector normal to the boundary surface $\partial\mathcal{U}_{j}$, and $\mathbf{n}_{p}$ is the outward unit vector normal to the boundary of the airspace. Note that $\left(\bigtriangledown\Phi_{a}\right)\cdot\hat{\mathbf{n}}_{o_{j}}$ vanishes because normal velocity component is zero at any point on the obstacle boundaries. Define boundary flux

[TABLE]

A key assumption of this paper is that

[TABLE]

To assign a nominal coordination, we assume that the left-hand side of Eq. (6) vanishes at any time $t$. This leads the nominal macroscopic coordination to be assigned as the solution of the Laplace Equation:

[TABLE]

Notice that $\Phi$ defines the steady state solution $\Phi=\Phi(\mathbf{r})$ if: (i) $q=q(\mathbf{r})$ is spatially varying and time-invariant at any point on the outer boundary $\partial\mathcal{P}$, and (ii) the geometry of unplanned space $\mathcal{U}$ does not change. Otherwise, nominal coordination can vary with time at a point $\mathbf{r}\in\mathcal{C}$, but $\Phi$ is still a solution of Laplace equation (3). In other words, both spatially-varying and spatiotemporally-varying solutions of nominal (reference) coordination $\Phi$ satisfy the Laplace Equation.

Next, Section III develops numerical and analytic solutions to the macroscopic coordination governing PDE. Microscopic coordination of each UAV/agent cluster is then discussed in Section IV. Collectively these formulations offer numerical and analytic approaches to solve the governing coordination equation given in (7).

III-A1 Numerical Solution

This paper applies the finite difference (FD) method to assign $\Phi$ over $\mathcal{C}$. Let a connected network with $m$ nodes, defined by $\mathcal{V}=\{1,2,\cdots,m\}$, be distributed over airspace $\mathcal{P}\bigcup\partial\mathcal{P}$. It is assumed that the geometry of airspace outer border $\partial\mathcal{P}$ is fixed. Therefore, the node set $\mathcal{V}$ is time-invariant, e.g. $\left|\mathcal{V}\right|=m$ is unchanged. Let $\mathcal{V}_{u}=\{i_{1},\cdots,i_{m_{u}}\}$ define those grid nodes distributed over unplanned airspace $\mathcal{U}$ and $\mathcal{V}_{c}=\mathcal{V}\setminus\mathcal{V}_{u}$ define the remaining grid nodes distributed over planned airspace $\mathcal{C}$. $\mathcal{V}_{c}$ can be expressed as $\mathcal{V}_{c}=\mathcal{V}_{cb}\bigcup\mathcal{V}_{ci}$, where

$\mathcal{V}_{cb}=\left\{j_{1},\cdots,j_{m_{cb}}\right\}$ defines the outer boundary control nodes distributed over $\partial\mathcal{P}$, 2. 2.

$\mathcal{V}_{ci}=\left\{l_{1},\cdots,l_{m_{ci}}\right\}$ defines the interior control nodes distributed over $\mathcal{C}\setminus\left(\partial\mathcal{U}\bigcup\partial\mathcal{P}\right)$.

Because unplanned airspace can be potentially dynamic, $m_{u}(t)=\left|\mathcal{V}_{u}\right|$ can change with time. Therefore, $m_{ci}(t)=\left|\mathcal{V}_{ci}\right|$ can be time-varying as well. However, $m_{cb}+m_{ci}(t)+m_{u}(t)=m$ remains unchanged at any time $t$.

Nodal $\Phi$ values, defined by ${\mathbf{\Phi}}=[\Phi_{i}]\in\mathbb{R}^{m\times 1}$ are found by solving

[TABLE]

where $\Phi_{i}\left(t\right)=\Phi\left(x_{i},y_{i},z_{i},t\right)$ ($i\in\mathcal{V}$), ${\mathbf{q}}\left(t\right)=[q_{i}]\in\mathbb{R}^{m\times 1}$ is the flux vector, and $\mathbf{L}\left(t\right)\in\mathbb{R}^{m\times m}$ is the positive semi-definite Laplacian matrix.

Define $\hat{{\mathbf{\Phi}}}=\mathbf{S}{\mathbf{q}}=\begin{bmatrix}\hat{\mathbf{\Phi}}_{cb}^{T}&\hat{\mathbf{\Phi}}_{ci}^{T}&\hat{\mathbf{\Phi}}_{u}^{T}\end{bmatrix}^{T}$,

$\hat{\mathbf{q}}=\mathbf{S}\mathbf{q}=\begin{bmatrix}\hat{\mathbf{q}}_{cb}^{T}&\hat{\mathbf{q}}_{ci}^{T}&\hat{\mathbf{q}}_{u}^{T}\end{bmatrix}^{T}$, and

[TABLE]

where superscript “$T$” denotes matrix transpose and orthogonal matrix $\mathbf{S}=[\mathbf{S}_{hq}]\in\mathbb{R}^{m\times m}$, called the interface matrix, is given by

[TABLE]

Note that $\hat{\mathbf{q}}_{ci}=\mathbf{0}\in\mathbb{R}^{m_{ci}\times 1}$ and $\hat{\mathbf{q}}_{u}=\mathbf{0}\in\mathbb{R}^{m_{u}\times 1}$. Let $\mathbf{S}(t)=\begin{bmatrix}\mathbf{S}_{cb}\\ \mathbf{S}_{cu}(t)\\ \mathbf{S}_{ci}(t)\end{bmatrix}$, where $\mathbf{S}_{cb}\in\mathbb{R}^{m_{cb}\times m}$, $\mathbf{S}_{cu}(t)\in\mathbb{R}^{m_{cu}\times m}$, and $\mathbf{S}_{ci}(t)\in\mathbb{R}^{m_{ci}\times m}$. By substituting $\mathbf{L}=\mathbf{S}^{T}\hat{\mathbf{L}}{\mathbf{S}}$, $\mathbf{q}=\mathbf{S}^{T}\hat{\mathbf{q}}$, $\mathbf{{\mathbf{\Phi}}}=\mathbf{S}^{T}\hat{\mathbf{{\mathbf{\Phi}}}}$, Eq. (8) is converted to $\hat{\mathbf{L}}\hat{{\Phi}}=\hat{{\mathbf{q}}}$, or

[TABLE]

We assume that $\hat{\mathbf{\Phi}}_{cb}$ is given, e.g. $\Phi$ values of the outer boundary nodes are not impacted by the $\Phi$ values of the in-neighbor nodes. Therefore, $\hat{\mathbf{L}}_{cb,cu}=\mathbf{0}$, $\hat{\mathbf{L}}_{cb,ci}=\mathbf{0}$, and $\hat{\mathbf{L}}_{cb,u}=\mathbf{0}$. Furthermore, $\hat{\mathbf{\Phi}}_{u}=\mathbf{0}$ at any time $t$. Therefore, the unplanned nodes defined by $\mathcal{V}_{u}$ are independent: $\hat{\mathbf{L}}_{u,cb}=\mathbf{0}$, $\hat{\mathbf{L}}_{u,ci}=\mathbf{0}$, and $\hat{\mathbf{L}}_{u,cu}=\mathbf{0}$. Consequently,

[TABLE]

where $\hat{\mathbf{L}}_{cb,cb}$ is diagonal and positive definite, $\mathbf{I}_{m_{u}}\in\mathbb{R}^{m_{u}\times m_{u}}$ is the identity matrix, and $\mathbf{A}_{c}$ has the following properties:

Diagonal entries of matrix $\mathbf{A}_{c}$ are all negative. 2. 2.

Off-diagonal entries of matrix $\mathbf{A}_{c}$ are all non-negative. 3. 3.

The sum of the entries are non positive for every row of matrix $\mathbf{A}_{c}$.

Theorem III.1

If the network distributed over the airspace is connected, matrix $\mathbf{A}_{c}$ is Hurwitz.

Steady-state coordination Eq. (8) can be rewritten as

[TABLE]

where $\mathbf{\Phi}_{I}=\begin{bmatrix}\hat{\mathbf{\Phi}}_{ci};\hat{\mathbf{\Phi}}_{u}\end{bmatrix}\in\mathbb{R}^{\left(m-m_{cb}\right)\times 1}$.

III-A2 Analytic Solution

For the analytic solution, the airspace is decomposed into $n_{s}$ motion surfaces (or floors) defined by $\Gamma(x,y,z,t)=\Gamma_{f}$ where $\Gamma_{f}$ is a constant for each $f=1,\cdots,n_{s}$. A position on floor $f$ can be expressed as $\mathbf{r}_{f}=\mathbf{r}\left(\Phi_{f},\Psi_{f},\Gamma_{f},t\right)$, where $\Phi_{f}=\Phi_{f}\left(x,y,t\right)$ is the potential function of floor $f$ satisfying the governing PDE in (3), and $\Psi_{f}=\Psi_{f}\left(x,y,t\right)$ is a stream function.

Remark 1: The level curves $\Phi_{f}\left(x,y,t\right)=\rm{constant}$ and $\Psi_{f}\left(x,y,t\right)=\rm{constant}$ perpendicularly cross at the intersection point in the image plane $x-y$. Using the Cauchy Riemann theorem [40], $\Psi_{f}$ and $\Phi_{f}$ are related as follows:

[TABLE]

where $u_{f}$ and $v_{f}$ are the nominal velocity components on floor $f$ projectes on the $x-y$ plane. Note that $\Phi_{f}$ and $\Psi_{f}$ both satisfy Laplace equations $\bigtriangledown^{2}\Phi_{f}=0$ and $\bigtriangledown^{2}\Psi_{f}=0$.

Let position vector $\mathbf{r}_{f}$ be expressed as

[TABLE]

where $\mathbf{r}_{xy}=x\hat{\mathbf{e}}_{x}+y\hat{\mathbf{e}}_{y}$ is the projection of $\mathbf{r}$ on the $x-y$ plane (image plane), and $\hat{\mathbf{e}}_{x}$, $\hat{\mathbf{e}}_{y}$, and $\hat{\mathbf{e}}_{z}$ are the bases of the Cartesian coordinate system. Macroscopic coordination is treated as ideal fluid flow or irrotational flow [41] defined in the $x-y$ image plane per the Fig. 2 schematic. “Uniform Flow”, “Source”, “Sink”, and “Doublet” are ideal flow patterns considered in this paper. Potential and stream functions $\Phi$ and $\Psi$ of these ideal flow patterns are given in Table I. Define unplanned airspace characteristic point $\mathbf{r}_{0,i,f}=[x_{0,i,f}\leavevmode\nobreak\ y_{0,i,f}]\in\mathcal{U}_{i}$ ($i=1,\cdots,n_{u}$). Unplanned airspace $\mathcal{U}_{i}$ can be excluded from $\mathcal{P}$ by combining either a uniform flow and doublet or a uniform flow, a sink and a source. Mathematically speaking, we define

[TABLE]

where $\mathbf{U}_{\infty,f}=u_{\infty,f}\mathbf{n}_{f}$ is the traffic flow direction projected on the horizontal $x-y$ plane and $\mathbf{n}_{f}(t)=\cos\theta_{0,f}(t)\hat{\mathbf{e}}_{x}+\sin\theta_{0,f}(t)\hat{\mathbf{e}}_{y}$ is a unit vector. $\gamma_{i,f}$ and $\Xi_{f}$ are binary parameters with values [math] or $1$. If floor $f$ contains no unplanned space, then $\Xi_{f}=0$. Otherwise, $\Xi_{f}=1$. If $\gamma_{i,f}=1$, unplanned airspace $\mathcal{U}_{i}$ is excluded by combining uniform flow and doublet flow; otherwise, $\mathcal{U}_{i}$ is excluded by combining uniform flow, a sink, and a source.

The boundary of the unplanned airspace is defined by

[TABLE]

In the example shown in Fig. 2, unplanned airspace $\mathcal{U}$ is separated from the planned space by combining uniform and doublet flows. Assuming $n_{s}=1$, $n_{u}=1$, $\theta_{0,f}=0$, and $\mathbf{r}_{1,f}=\mathbf{0}$,

[TABLE]

[TABLE]

define the potential and stream fields over projection plane $x-y$, respectively. $\Psi_{f}=0$ defines the boundary of the unplanned airspace excluding an obstacle by a circle with radius $\sqrt{\frac{\Delta_{f}}{u_{\infty,f}}}$.

Remark 2: $u_{\infty,f}$, $\Delta_{i,f}$, and $\mathbf{r}_{0,i,f}$ are design parameters to appropriately specify the geometry of unplanned airspace $\mathcal{U}_{i}\in\mathbf{\mathcal{missing}}{U}$. By increasing $u_{\infty,f}$ the size of the unplanned air space is decreased on floor $f$. On the other hand, the unplanned airspace $\mathcal{U}_{i}$ is enlarged when $\Delta_{i}$ is increased. Unplanned airspace $\mathcal{U}_{i}$ can be displaced in airspace $\mathcal{P}$ by updating $\mathbf{r}_{0,i,f}$.

III-B Resilient Macroscopic Coordination

To realize a robust and efficient macroscopic coordination strategy, it is desirable to minimize the nominal velocity variation imposed by unplanned airspace displacement and deformation when geometry of the unplanned airspace changes with time. To this end, an analytic boundary control solution will be designed in Section III-B1 to handle unpredicted (pop-up) unplanned airspace sectors cut out from the planned airspace. Furthermore, a numerical boundary control will be developed in Section III-B2 to enhance macroscopic coordination in cases where an unplanned airspace volume can either be stationary or move predictably over time.

III-B1 Analytic Boundary Control

We define the resilient macroscopic coordination over the floor $f$ as the solution of the following governing PDE:

[TABLE]

subject to the boundary condition

[TABLE]

where $\psi_{f,a}$ is the actual stream function defined over planned space $\mathcal{C}$. Note that $U_{p}(\mathbf{r},t)$ and $U_{u}(\mathbf{r},t)$ are the boundary controls actuated at boundaries $\mathcal{\partial P}$ and $\partial\mathcal{U}$, respectively.

Proposition 1

Let $\mathbf{n}_{u_{j}}$ and $\mathbf{t}_{u_{j}}$ denote normal and tangent unit vectors at unplanned boundary $\partial\mathcal{U}_{j}$ per Fig. 1. $\psi_{f}(\mathbf{r},t)$ is constant at any point $\mathbf{r}$ on the boundary of the planned airspace.

Considering Proposition 1, reference macroscopic coordination can be obtained as the solution

[TABLE]

subject to boundary condition

[TABLE]

Theorem III.2

Let $E_{f}(\mathbf{r},t)=\psi_{f,a}(\mathbf{r},t)-\psi_{f}(\mathbf{r},t)$ be the distributed error function defining the difference between actual and reference stream values over the planned space ($\mathbf{r}\in\mathcal{C},\leavevmode\nobreak\ t\geq 0$). $E_{f}(\mathbf{r},t)$ governed by

[TABLE]

${E_{f,a}}(\mathbf{r},t)$* asymptotically tends to [math], if the boundary controls $U_{p}$ and $U_{u}$ are chosen as follows:*

[TABLE]

where $k_{p}>0$ and $k_{u_{j}}>0$ are constant control gains.

III-B2 Numerical Boundary Control

Discretizing governing equation (2) yields:

[TABLE]

where ${\mathbf{\Phi}}_{I}^{a}\in\mathbb{R}^{\left(m-m_{cb}\right)\times 1}$ defines actual coordination given actual boundary input $\hat{\mathbf{\Phi}}_{cb}^{a}\in\mathbb{R}^{m_{B}\times 1}$ where superscript $a$ denotes actual. Define $\mathbf{E}={\mathbf{\Phi}}_{I}-{\mathbf{\Phi}}_{I}^{a}$ and $\mathbf{U}=\hat{\mathbf{\Phi}}_{cb}-\hat{\mathbf{\Phi}}_{cb}^{a}$. By subtracting Eq. (12) from Eq. (24), the error dynamics is obtained:

[TABLE]

The paper designs macroscopic control as a linear quadratic regulator (LQR) minimizing

[TABLE]

where $\mathbf{W}_{e}\in\mathbb{R}^{\left(m-m_{cb}\right)\times\left(m-m_{cb}\right)}$ is positive semi-definite and $\mathbf{W}_{u}\in\mathbb{R}^{m_{cb}\times m_{cb}}$ is positive-definite. The above cost function is minimized by choosing $\mathbf{U}=-\mathbf{K}_{e}\mathbf{E}$, where $\mathbf{K}_{E}=\mathbf{W}_{u}^{-1}\mathbf{B}_{c}^{T}\mathbf{P}\in\mathbb{R}^{m_{cb}\times(m-m_{cb})}$ is the LQR control gain and $\mathbf{P}_{e}\in\mathbb{R}^{\left(m-m_{cb}\right)\times\left(m-m_{cb}\right)}$ is the solution of the algebraic Ricatti equation [42]:

[TABLE]

IV Microscopic Coordination

First define the following position notations:

• •

$\mathbf{p}_{{}_{j,i}}$: Actual position of agent $j$ in cluster $i\in\mathcal{S}$

• •

$\mathbf{p}_{{}_{d,j,i}}$: Local desired position of agent $j$ in cluster $i\in\mathcal{S}$

• •

$\mathbf{p}_{{}_{RB,j,i}}$: Global desired position of agent $j$ in cluster $i\in\mathcal{S}$

• •

$\mathbf{p}_{{}_{0,j,i}}$: Material (Relative) position of agent $j$ in cluster $i\in\mathcal{S}$

• •

$\mathbf{r}_{i}$: Macroscopic desired position of cluster $i\in\mathcal{S}$

Assume $N(t)$ UAS clusters exist in $\mathcal{P}$ at time $t$. Existing clusters are defined by $\mathcal{S}(t)=\{1,\cdots,N(t)\}$. We assume that $\mathcal{V}_{i}=\{1,\cdots,n_{i}\}$ defines the UAS belonging to existing cluster $i\in\mathcal{S}$. Cluster $i$ contains either a single UAS or multiple UAS. With multiple UAS, graph $\mathcal{G}_{i}=\left(\mathcal{V}_{i},\mathcal{E}_{i}\right)$ defines interagent communication in cluster $i\in\mathcal{S}$ where $\mathcal{E}_{i}\subset\mathcal{V}_{i}\times\mathcal{V}_{i}$ specifies edges of of graph $\mathcal{G}_{i}$.

Agents (UAS) of cluster $i\in\mathcal{S}$ are considered particles of a $d_{i}$-D group contained by a $d_{i}$-D ($d_{i}=0,1,2,3$) simplex, called the leading simplex. Define $\mathcal{V}_{i}=\mathcal{V}_{L,i}\bigcup\mathcal{V}_{F,i}$ where

[TABLE]

specify identification numbers of leaders and followers, respectively. Leaders $1$ through $d_{i}+1$ are placed at vertices of the leading simplex. Therefore,

[TABLE]

Remark 3: If $d_{i}=0$, cluster $i\in\mathcal{S}$ consists of a single UAS, $\mathcal{V}_{i}=\mathcal{V}_{L,i}$, and $\mathcal{V}_{F,i}=\emptyset$.

Assume each $d_{i}$-D cluster is guided by $d_{i}+1$ leaders placed at vertices of the leading simplex. While leaders move independently, followers acquire desired coordination through local communication. The in-neighbor set of follower $j\in\mathcal{V}_{F,i}$ ($i\in\mathcal{S}$) is defined by

[TABLE]

In this paper, agents are modeled as double integrators:

[TABLE]

where $\beta_{1,i}$ and $\beta_{2,i}$ are positive control gains,

[TABLE]

$w_{i,j,h}>0$, and

[TABLE]

Given communication weights, weight matrix $\mathbf{W}_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$ is defined as follows:

[TABLE]

Matrix $\mathbf{W}_{i}$ can be partitioned as

[TABLE]

where $\mathbf{I}_{d_{i}+1}\in\mathbb{R}^{\left(d_{i}+1\right)\times\left(d_{i}+1\right)}$ is the identity matrix.

Theorem IV.1

If graph $\mathcal{G}_{i}$ defining inter-agent communication in cluster $i\in\mathcal{S}$ contains a spanning tree, and communication weights are all positive and defined as given in Eq. (30), then matrices $\mathbf{L}_{i}\in\mathbb{R}^{(n_{i}-d_{i}-1)\times(n_{i}-d_{i}-1)}$ and $\mathbf{W}_{i}$ are Hurwitz.

IV-A Microscopic Coordination and Control

Global desired position $\mathbf{p}_{{}_{RB,j,i}}$ is defined by the following affine transformation:

[TABLE]

where $\mathbf{Q}_{i}\in\mathbb{R}^{3\times 3}$ is an orthogonal rotation matrix. Agents’ material (relative) position $\mathbf{p}_{{}_{0,j,i}}$ ($j\in\mathcal{V}_{i},\leavevmode\nobreak\ i\in\mathcal{S}$) is assumed constant. Note that $\mathbf{p}_{{}_{\rm{RB},j,i}}=\mathbf{\mathbf{p}_{{}_{d,j,i}}}$ if agent $j\in\mathcal{V}_{L,i}$ is a leader. However, $\mathbf{p}_{{}_{\rm{RB},j,i}}$ differs from the local desired position $\mathbf{p}_{{}_{d,j,i}}$ expressed as the convex combination of in-neighbor agents per Eq. (29) if agent $j\in\mathcal{V}_{F,i}$ is a follower. Furthermore, macroscopic desired position $\mathbf{r}_{i}$ is assigned at the macroscopic level and updated by Eq. (2):

[TABLE]

subject to the arrival condition $\mathbf{r}_{i}(t_{0,i})=\mathbf{r}_{0,i}\in\partial\mathcal{P}$. Note that $t_{0,i}$ is the time at which cluster $i\in\mathcal{S}$ crosses border $\partial\mathcal{P}$ to enter the planned airspace at $\mathbf{r}_{0,i}\in\partial\mathcal{P}$, and $t_{f,i}$ is the time at which cluster $i$ leaves the planned airspace at border $\mathbf{r}_{f,i}\in\partial\mathcal{P}$, e.g. $\mathbf{r}_{f,i}=\mathbf{r}_{i}(t_{f,i})$ ($i\in\mathcal{S}$). Given $t_{0,i}$ and $\mathbf{r}_{0,i}$ reference desired trajectory $\mathbf{r}_{i}(t)$ is assigned over $t\in[t_{0,i},t_{f,i}]$.

Remark 4: Given arrival condition $\mathbf{r}_{i}(t_{0,i})=\mathbf{r}_{0,i}\in\partial\mathcal{P}$, desired reference trajectory $\mathbf{r}_{i}(t)$ ($t\in[t_{0,i},t_{f,i}]$) is unique.

Cluster Rotation Matrix $\mathbf{Q}_{i}$: Agents of cluster $i\in\mathcal{S}$ are treated as particles of a virtual rigid body (VRB) with a desired evolution given by rigid body motion kinematics. Defining a body frame with orthogonal unit bases $\mathbf{i}_{m,i}$, $\mathbf{j}_{m,i}$, and $\mathbf{k}_{m,i}$, $\mathbf{i}_{m,i}$, $\mathbf{j}_{m,i}$, and $\mathbf{k}_{m,i}$ are related to ground coordinate bases $\mathbf{e}_{x}$, $\mathbf{e}_{y}$, and $\mathbf{e}_{z}$ by

[TABLE]

where

[TABLE]

is an orthogonal matrix and $\theta_{1,i}$ and $\theta_{2,i}$ are independent rotation angles. This paper assumes that $\mathbf{i}_{m,i}$ is along the macroscopic desired velocity $\dot{\mathbf{r}}_{i}$ at any time $t$. Therefore,

[TABLE]

Note that $\dot{\mathbf{r}}_{i}=\mathbf{V}(\mathbf{r}_{i},t)$ has been assigned at the macroscopic level per Eq. 1. Given reference velocity $\dot{\mathbf{r}}_{i}$, rotation angles $\theta_{1,i}$ and $\theta_{2,i}$ are obtained as follows:

[TABLE]

IV-B Microscopic Collective Dynamics

Let $\mathbf{p}_{j,i}=[p_{x,j,i}\leavevmode\nobreak\ p_{y,j,i}\leavevmode\nobreak\ p_{z,j,i}]^{T}\in\mathbb{R}^{3\times 1}$ denote actual position of agent $j\in\mathcal{V}_{i}$. Let $\mathbf{p}_{\rm{RB},j,i}=[x_{\rm{RB},j,i}\leavevmode\nobreak\ y_{\rm{RB},j,i}\leavevmode\nobreak\ z_{\rm{RB},j,i}]^{T}\in\mathbb{R}^{3\times 1}$ denote global desired position of agent $j\in\mathcal{V}_{i}$. Define

[TABLE]

Assuming agent $j\in\mathcal{V}_{i}$ updates its position according to the dynamics in Eq. (28), collective dynamics of cluster $i\in\mathcal{S}$ is given by

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Note that “$\otimes$” is the Kronecker product symbol, $\mathbf{0}_{n_{i}\times(d_{i}+1)}\in\mathbb{R}^{n_{i}\times(d_{i}+1)}$ and $\mathbf{0}_{(n_{i}-d_{i}-1)\times(d_{i}+1)}\in\mathbb{R}^{(n_{i}-d_{i}-1)\times(d_{i}+1)}$ are the zero-entry matrices and $\mathbf{I}_{d_{i}+1}\in\mathbb{R}^{(d_{i}+1)\times(d_{i}+1)}$ is the identity matrix. Furthermore, $\rm{vec}(\cdot)$ is the matrix vectorization operator. For example,

[TABLE]

Because $\mathbf{W}_{i}$ is Hurwitz per Theorem IV.1, $\mathbf{A}_{\rm{SYS},i}$ is Hurwitz as well. Therefore, the microscopic collective dynamics of cluster $i\in\mathcal{S}$ is stable.

V Simulation Results

A series of simulation results are presented to illustrate the proposed hierarchical coordination strategy for LAUC. First, macroscopic coordination of clustered UAS with different nominal operating speeds are shown in Section V-A. Next, simulation results of resilient macroscopic coordination are presented in Section V-B. Section V-C simulates microscopic coordination of an individual cluster operating with a given planned airspace region.

V-A Heterogeneous Macroscopic Coordination

In this example, a finite airspace consisting of a single unplanned space ($n_{u}=1$) is considered, and the analytic approach developed in Section III-A2 is used to exclude the unplanned space. Assuming the airspace consists of a single floor ($n_{s}=1,\leavevmode\nobreak\ z_{f}=0$), navigable channels are formed by combining a doublet flow and a uniform flow. Therefore, Eq. (17) defines navigable channel geometry where $\Psi_{f}=\rm{const}$ specifies navigable channel boundaries. In Fig. 4, naviagble channel boundaries are shown by green. Choosing $u_{\infty}=40$ and $\Delta=4000$, $\psi_{f}=0$ is a closed circle of radius $20m$ excluding the unplanned airspace.

In this example, it is desired that nominal speed remains constant along every channel. Therefore, cost function $K=K_{f}$ ($f=1$) remains constant along every channel. To deal with UAS heterogeneity at the microscopic scale, UAS entering the airspace are clustered into $6$ groups based on their nominal speed. Assuming nominal coordination and speed are admitted by the existing UAS in the planned space, Fig. 5 shows the six UAS clusters at time $t=30s$.

V-B Resilient Macroscopic Coordination

An example of resilient macroscopic coordination is shown in Fig. 6. We assume that a UAS reports a failure (emergency) while moving along a safe navigable zone at $t=0$; note the failed UAS is shown by a red circle in Fig. 6 (a)). Afterward, the failed UAS is enclosed by a rectangular unplanned airspace $\mathcal{U}_{D}$ per Fig. 6 (a). Given the geometry of the unplanned airspace, the numerical method developed in Section III-A1 is used to update the desired macroscopic coordination as the solution of Eq. (12). Then, we design an LQR boundary control law to redefine the navigable channels so that the unplanned airspace is safely excluded. Hence, $\mathbf{U}=-\mathbf{K}_{e}\mathbf{E}$ is chosen, where $\mathbf{K}_{e}$ is the LQR gain obtained as the solution of the algebraic Ricatti equation [43]. Recovery of safe navigation zones is illustrated in Fig. 6 (b-d). This example illustrates how a pop-up unplanned airspace obstacle can be handled.

V-C Microscopic Coordination Results

As a third example, consider the microscopic coordination of cluster $i$ consisting of $n_{i}=10$ agents. Cluster $i$ treats each UAS as a particle of a $2$-D virtual rigid body (VRB). Three leaders with indices $1$, $2$, and $3$ guide the collective motion of the cluster. The remaining agents with indices $4$ through $10$ are followers acquiring the desired coordination through local communication. The graph shown in Fig. 7 defines inter-agent communication in cluster $i$. Followers’ communication weights are listed in Table III.

For this case study cluster $i$ is instructed through macroscopic control to move on the surface

[TABLE]

Assume the single-floor airspace consists of a single unplanned space. The unplanned airspace is excluded by combining a uniform flow and a doublet flow. Therefore

[TABLE]

define the potential and stream fields over the planned airspace, respectively. It is desired that cluster $i$ moves along the path $\mathbf{r}_{i}(t)=x_{i}(t)\hat{\mathbf{e}}_{x}+y_{i}(t)\hat{\mathbf{e}}_{y}+z_{i}(t)\hat{\mathbf{e}}_{z}$ where $x_{i}(t)$, $y_{i}(t)$, and $z_{i}(t)$ satisfy $psi_{f}(x_{i},y_{i})=375$ and $z_{i}=1000-5\times 10^{-3}\left[\left(x_{i}-25\right)^{2}+y_{i}^{2}\right]$. The macroscopic desired path $\mathbf{r}_{i}=x_{i}\hat{\mathbf{e}}_{x}+y_{i}\hat{\mathbf{e}}_{y}+z_{i}\hat{\mathbf{e}}_{z}$ is shown by a green curve in Fig. 9. Furthermore, components of macroscopic desired trajectory $\mathbf{r}_{i}$ ($x_{i}(t),\leavevmode\nobreak\ y_{i}(t),\leavevmode\nobreak\ z_{i}(t)$) are plotted by dashed red in Fig. 8 (a-c). Fig. 8 (a-c) also illustrates $x$, $y$, and $z$ components of the agents’ actual positions versus time. Continuous curves plot actual position components $p_{x,j,i}$, $p_{y,j,i}$, and $p_{x,j,i}$ versus time for every agent $j\in\mathcal{V}_{i}=\{1,\cdots,10\}$. Given $\mathbf{r}_{i}(t)$ macroscopic rotation angles $\theta_{1,i}(t)$ and $\theta_{2,i}(t)$ are computed at each time $t$ and plotted versus time in Fig. 8(d). Moreover, Fig. 8(e) plots deviation of follower $7$ ($\|\mathbf{p}_{7,i}-\mathbf{p}_{\rm{RB},7,i}\|$) versus time $t$. The formation of cluster $i$ at sample times $t=0s$, $t=40s$, $t=80s$, $t=120s$, $t=160s$, and $t=190s$ are shown in Fig. 9. As illustrated, cluster $i$ successfully tracks the macroscopic desired trajectory $\mathbf{r}_{i}$.

VI   Conclusion

The proposed hierarchical LAUC strategy uses an Eulerian description of continuum mechanics to manage the space and time allocated to UAS clusters to enable heterogeneous UAS clusters to safely coordinate motion through a complex airspace sector. To this end, nominal coordination is obtained as the solution of a governing PDE with spatiotemporal parameters. Assuming airspace is finite, this paper manages the airspace capacity through controlling UAS inflow and outflow at the airspace boundaries. UAS heterogeneity is managed at both macroscopic and microscopic levels. At the macroscopic level, planned airspace can be divided into safe navigable channels and navigable channels can be allocated to different UAS classes with different nominal flight speeds. At the microscopic level, similar UAS are treated as particles of a rigid-body cluster. A leader-follower containment method is applied by each cluster to acquire the prescribed macroscopic coordination in a decentralized fashion. The paper also studies LAUC resilience individual vehicle failure by wrapping a failed UAS with a no-fly zone defined as a ”pop-up” unplanned airspace region other UAS never enter. We offer a resilient linear quadratic regulator (LQR) boundary control strategy to update the geometry of the unplanned airspace and safely plan nominal coordination in the macroscopic layer.

This paper focuses on the efficient flow of UAS through a single LAUC airspace sector, prescribing constraints on UAS entry and exit stations based on UAS travel speed and macroscopic channel definitions. Follow-on research is required to assure UAS will be able to safely and efficiently coordinate transitions between adjacent planned airspace sectors.

Appendix A Proofs

Proof for Theorem III.1: Let $\mathbf{A}_{c}$ be expressed as

[TABLE]

where $\mathbf{G}\in\mathbb{R}^{\left(m-m_{cb}\right)\times\left(m-m_{cb}\right)}$ is diagonal and positive definite matrix with diagonal entry

[TABLE]

The diagonal entries of matrix $\hat{\mathbf{A}}_{c}$ are all $-1$ while off-diagonal entries are all non-negative. Furthermore, the sums of the row entries of $\hat{\mathbf{A}}_{c}$ are non-negative while the sum of the entries is negative in at least one row of matrix $\hat{\mathbf{A}}_{c}$. Let $\hat{\mathbf{A}}_{c}$ be expressed as $\hat{\mathbf{A}}_{c}=-\mathbf{I}+\mathbf{D}$.

The spectral radius of matrix $\mathbf{D}$, denoted $\rho(\mathbf{D})$, is less than $1$. Therefore, $-\hat{\mathbf{A}}_{c}$ is a non-singular M-matrix with eigenvalues located inside a disk of radius $\rho(\mathbf{D})<1$ centered at $-1+0\mathbf{j}$. Therefore, $\hat{\mathbf{A}}_{c}$ and $\hat{\mathbf{A}}_{c}$ are both Hurwitz.

Proof for Proposition 1: At $\mathbf{r}\in\mathcal{U}_{j}$

[TABLE]

Therefore, $\bigtriangledown\Phi$ is along the tangent unit vector $\mathbf{t}_{u_{j}}$ at any point $\mathbf{r}\in\mathcal{U}_{j}$ ($j=1,\cdots,n_{u}$). $\bigtriangledown\Psi_{f}(\mathbf{r},t)$ is normal to $\bigtriangledown\Phi_{f}(\mathbf{r},t)$ per Remark 1. Thus, $\bigtriangledown\Phi_{f}(\mathbf{r},t)$ remains [math] along level curve $\psi_{f}\left(\mathbf{r},t\right)$ which is constant at the boundary of the unplanned airspace.

Proof for Theorem III.2: Let

[TABLE]

be a Lypunpov candidate function defined over planned space $\mathcal{C}$. Then,

[TABLE]

Choosing ${U}_{p}$ and $U_{u_{j}}$ according to Eqs. (23a) and (23b) and considering the error governing equation (22), $\dot{V}$ can be converted to

[TABLE]

Note that $\dot{V}$ in Eq. (A) is negative semi-definite where $\dot{V}=0$ is the invariant set of the governing Eq. (22). Using LaSalle’s Invariance Principle, it is concluded that $E_{f}(\mathbf{r},t)$ assymptotically converges to [math] at any point $\mathbf{r}\in\mathcal{C}$ as $t\rightarrow\infty$.

Proof for Theorem IV.1: The sum of the row entries are [math] at row $d_{i}+2$ through $n_{i}$ of matrix $\mathbf{W}_{i}$ while the diagonal entries of matrix $\mathbf{L}_{i}$ are all $-1$. Because every leader is an in-neighbor agent to at least one follower, $d_{i}+1$ rows of matrix $\mathbf{\Omega}_{i}$ have at least one positive entry. Consequently, the sum of the row elements is negative in at least $d_{i}+1$ rows of matrix $\mathbf{L}_{i}$ while the remaining rows of matrix $\mathbf{L}_{i}$ are zero-sum-row. Therefore, the spectral radius of matrix $\mathbf{L}_{i}+\mathbf{I}_{n_{i}-d_{i}-1}$ is less that $1$ and matrix $\mathbf{L}_{i}$ is Hurwitz. If $\mathbf{L}_{i}$ is Hurwitz, $\mathbf{W}_{i}$ is Hurwitz as well.