A Hamilton–Jacobi Reachability-Based Minimum Separation Estimation of Integrated Manned and Unmanned Operation in Uncertain Environments

Abstract

This work presents a minimum separation calculation for the integrated operation of manned and unmanned aerial vehicles in an uncertain airspace environment. Different from traditional path-planning-based research, this study investigated the minimum safe separation distance from a novel perspective of reachability analysis. The proposed computational method made use of the Hamilton–Jacobi partial differential equation (HJPDE) to obtain the backward reachable tube. Firstly, this work modeled the integrated operation in the UAS traffic management scenario, particularly focusing on the uncertainties. Then, a probabilistic reachability tube computation method was derived. Next, this work calculated the safe separation distances based on reachability analysis for three scenarios: a deterministic environment, an environment with relative position uncertainty, and an environment with relative heading angle uncertainty. By calculating the reachable tubes for a given response time, the worst-case minimum safe distances from the UAV’s perspective were determined, and the quantitative patterns were summarized. The results in this work indicate that, with an increase in the risk level and under the premise of a 1 s response time, the minimum safe separation increases from 26.7 m to 30.0 m. Finally, the paper discusses the results, explaining their rationality from both mathematical and physical perspectives.

1. Introduction

1.1. Background

With the growing application of unmanned aircraft systems (UASs) [1], UASs play an important role in the civil area of logistics delivery [2], search and rescue [3], and aerial surveillance [4]. Especially with the development of urban air mobility (UAM) concepts [5], the application of UASs in urban spaces is receiving increased attention. To integrate drones and UASs into airspace, the Federal Aviation Administration commissioned the National Aeronautics and Space Administration to conduct the UAS Traffic Management (UTM) project, which involves a series of studies and flight tests [6]. Europe has also initiated a similar project, called U-Space [7]. An essential objective of UTM and U-Space is to maintain flight safety among UASs by keeping proper separation distances [8]. However, achieving a safe transition from sparsely populated UASs in a single-mode operating environment to a complex, large-scale, and mixed environment with various types of aircraft integrated into the airspace presents a significant challenge due to differing UAS performances, a range of sensors, and the absence of comprehensive operational guidelines and separation standards [9].

Manned–unmanned teaming (MUM-T) has been long applied in the military field [10]. Military environments have more urgency and pose more communication issues to overcome; hence, the precision, speed, and tireless productivity inherent in MUM-T systems make it adaptable for commercial use [11]. Recent progress in the development of UAM vehicles and infrastructures is expanding the use of MUM-T, with the emphasis it places on the use of autonomous systems, UASs, remotely piloted aircraft, and vertical takeoff and landing (VTOL) aircraft. UTM systems are designed to integrate manned aircraft vehicle (MAV) and unmanned aircraft vehicle (UAV) operations in the same airspace [12], enhancing the safety and efficiency while ensuring minimal human intervention. However, due to the different dynamic characteristics of MAVs and UAVs, the safety guarantee of an MAV and a UAV is more complex than that between two UASs.

1.2. Related Work

A variety of methods have been applied to assess the flight safety in integrated airspace [9,13]. A category of existing methods involves improving traditional path planning algorithms [14] to research the low-altitude collision avoidance problem, such as the well-known $A^{★}$ algorithm [15], the rapidly exploring random tree algorithm [16], and the model predictive control algorithm [17]. Another category of existing methods is based on traditional air traffic management collision models, such as the Reich model and Monte Carlo modeling [18]. However, these methods are often based on experience and lack mathematical and physical significance. Moreover, only a few studies have been conducted under uncertainty. In [19], the authors stated that the location uncertainty or error of an aircraft is a critical issue to ensure safety. A risk domain formulation was incorporated, and the collision avoidance problem was solved by a sampling-based method [19] and the reinforcement learning method [20].

Reachability analysis [21] presents a powerful tool for providing safety and goal satisfaction guarantees to safety-critical autonomous systems under bounded disturbances. A backward reachable set was computed numerically by solving Hamilton–Jacobi (HJ) partial differential equations (PDEs) [22]. As a tool to solve the optimal control problem, HJ reachability has been applied to several vehicle-related scenarios such as automated vehicles [23], maneuver sequence design [24], aircraft autolander [25], and aerial refueling [26]. In the realm of UTM research, Chen proposed an operation model for platoons of unmanned aerial vehicles flying on air highways [27]. In this work, each vehicle was modeled as a hybrid system, and the reachability analysis method was utilized to guarantee the safety and goal satisfaction for all mode transitions. Moreover, the uncertainty problem could be addressed with probabilistic reachability analysis by formulating a stochastic continuous system and solving the related stochastic differential equation (SDE) [28]. When considering an uncertain environment, a Brownian motion process was added to the deterministic process [29], so that probability-based safety could be calculated with numerical approaches. In [30], considering the uncertainties and randomness introduced by wind perturbations in the aerial refueling process, probabilistic reachability was introduced to solve the docking risk evaluation problem through computing the docking success probability.

The benefit of reachability analysis is that the method provides a mathematical proof and calculation of a set (called a “reachable set”) where all states within it could reach or avoid another set of states (called a “target set”) in a specified time horizon with appropriate inputs despite any disturbances. Different from Monte Carlo simulations, reachability analysis considers all the possible states [22].

1.3. Motivation

Considering the fact that few methods have been applied to address the separation problem in uncertain integrated airspace, in this work, the reachability analysis method was studied. The motivation of this work was to calculate the safe distance between the aircraft by computing the reachable set, given the dynamic characteristics of an MAV and a UAV. By analyzing the computed results of reachable sets in different scenarios, the safe separation distances and the minimum LOS resolution time were concluded. A mid-air collision can be avoided if the separation distance is guaranteed. First, the related dynamic equations of an MAV and a UAV are given, and the problem formulation is derived. Then, the improved procedure of reachability analysis is introduced, and the target set and unsafe set are illustrated according to the integrated operation scenario. Then, a case study in integrated airspace of three different scenarios including deterministic and stochastic situations is conducted. Based on the calculation results, quantity analysis is conducted, and the minimum safe separation is obtained. Finally, the discussion section explains the results from both a mathematical and physical perspective.

Our contributions are summarized as follows:

• Differing from collision avoidance path planning methods, this work calculated the minimum safe separation for the integrated operation of manned and unmanned aerial vehicles from the perspective of reachability analysis.
• In the reachability analysis, this work considered uncertainties in the relative position and relative heading angle, computing probabilistic reachable tubes by calculating the expectation of the value function for different probabilities.
• By calculating the envelope boundaries under heading angles, this work determined the minimum safe separation and analyzed the results from both mathematical and physical perspectives.

The reminder of this paper is organized as follows: Section 2 presents the problem formulation of this work, including the research objective, reachability analysis formulation, and continuous dynamic systems modeling. Section 3 presents the solution method of the problem formulation, especially the stochastic reachability computing method. In Section 4, the numerical simulation and deep analysis are presented, including the quantitative results of three different scenarios; then, a discussion of the quantified results is presented in Section 5. Finally, Section 6 concludes the paper and sketches out potential directions for future work.

2. Problem Formulation

2.1. Problem Statement

This paper focuses on the problem setting where an MAV and UAV are cruising in mid-air with the intent to avoid conflict between them. As a problem of flight safety guarantee, this work formulated the relative position of the two aircraft as a differential equation. The aircraft’s positions were determined by its state and action, and the relative distance was determined by their time-variant coordinates. Depending on whether uncertainties are taken into account, the relative differential equations take the form of an ordinary differential equation (ODE) or an SDE, namely,

$$\dot{x}=f(t,x,u),x\in {\mathbb{R}}^n$$

(1)

or

$$\dot{x}=f(t,x,u)+\sigma (x\left(t\right),t)dB\left(t\right),x\in {\mathbb{R}}^n$$

(2)

where x and u are the state and input of the system, respectively. $B\left(t\right)$ is a Brownian motion process of appropriate dimension. The drift term f and the diffusion term $\sigma \left(x\right(t),t)dB\left(t\right)$ represent the deterministic component and the probabilistic component of the system evolution, respectively. The objective of the problem is to ensure a safe separation and relative minimum LOS resolution time between the MAV and UAV, in order to avoid a loss-of-separation (LOS) event, namely,

$$\exists t\in [0,T],\parallel {\mathrm{x}}_{Mt}-{\mathrm{x}}_{Ut}\parallel ⩽D_{MU}$$

(3)

where T denotes the interest time period, and ${\mathrm{x}}_{Mt}$ and ${\mathrm{x}}_{Ut}$ are the locations of the MAV and UAV at time t, respectively. $D_{MU}$ is the minimum separation of the two aircraft to ensure safety [19].

To achieve this goal, a set of states needs to be calculated in which at least one LOS event would occur no matter what action the UAV took to avoid the unsafe situation.

2.2. Reachability Analysis Formulation

Reachability analysis has been studied since the 1990s. As introduced above, based on the characteristics of dynamic system and optimal control theory, it has been successfully applied in the safety assessment problem of vehicles. In reachability analysis, the most commonly used approach is backward reachable set calculation [22]. In this work, a backward reachable set is defined as a set that satisfies

$$\mathcal{B}\left(t\right)=\{x\in {\mathbb{R}}^n:\exists u_M\in {\mathcal{U}}_M\left(t\right),\forall u_U\in {\mathcal{U}}_U\left(t\right),\zeta (0;x,t,u_U,u_M)\in \mathcal{D}\left(t\right)\}$$

(4)

where ${\mathcal{U}}_M\left(t\right)$ and ${\mathcal{U}}_U\left(t\right)$ denote the feasible set of strategies for an MAV and a UAV, respectively. $\zeta (s;x,t,u_U,u_M)$ are the trajectories of the relative dynamic equation starting from state x at time t to state at time s under control $u_U$ and $u_M$, in which $t\in [T,0]$. $\mathcal{D}\left(t\right)$ represents the unsafe set of states, called the “target set”, which should be achieved or avoided.

Intuitively, a backward reachable set is the set of states from which the system can reach a target at exactly time 0, while a backward reachable tube represents the set of states from which the system can reach or avoid a target within a duration of $\left|t\right|$ [31]. Therefore, in this case, rather than a backward reachable set, the backward reachable tube merits further consideration, because the unsafe set should be avoided at all times during the flight period. A backward reachable tube can be defined as

$$\mathcal{T}\left(t\right)=\{x\in {\mathbb{R}}^n:\exists u_M\in {\mathcal{U}}_M\left(t\right),\forall u_U\in {\mathcal{U}}_U\left(t\right),\exists s\in [T,0],\zeta (s;x,t,u_U,u_M)\in \mathcal{D}\left(t\right)\}$$

(5)

which means regardless of the UAV’s maneuvers, there is always at least one scenario in which an LOS event would occur.

As shown in Figure 1, trajectories originating from states inside the backward reachable set may cross the target set within the time horizon, while the state at the exact time 0 is located outside the target set. So, according to the definitions, the desired set, in this work called the “safe tube”, can be denoted as $\mathcal{D}\left(t\right)={\complement }_{{\mathbb{R}}^n}\mathcal{T}\left(t\right)$.

2.3. Continuous Dynamic Systems Formulation

In this work, to simplify the problem, it is assumed that the MAV and UAV are flying within a layered airspace structure [32], as shown in Figure 2. The airspace is divided into several parts with each part having a constant altitude; hence, the flight height of the aircraft does not change during the cruising period.

The following two-dimensional kinematic model was used to characterize the dynamic properties of the two aircraft [19].

$$\left\{\begin{array}{c}\dot{x_i}=v_i·cos\left({\psi }_i\right)\hfill \\ \dot{y_i}=v_i·sin\left({\psi }_i\right)\hfill \\ \dot{{\psi }_i}=a_{{\psi }_i}\hfill \end{array}\right.i=M,U$$

(6)

in which $v_U$ and $v_M$ are the velocities of the UAV and MAV, respectively, and are considered as constants in this work. ${\psi }_U$ and ${\psi }_M$ denote the two aircrafts’ headings and are driven by angular velocities $a_{{\psi }_U}$ and $a_{{\psi }_M}$.

As shown in Figure 3, the velocities and positions of the UAV and MAV are defined within their respective body frames. Similarly, their relative positions include distances in two directions as well as the heading difference. Without loss of generality, by setting the origin of the relative position coordinate system at the center of the UAV and aligning the x-axis direction with the velocity direction of the UAV, the relative position relationship between the two aircraft can be solved in a rotating coordinate system. The relative dynamic equations are

$$\left[\begin{array}{c}\dot{x_r}\\ \dot{y_r}\\ \dot{{\psi }_r}\end{array}\right]=\left[\begin{array}{c}{a}_{{\psi }_U}{y}_r-v_U+v_{M}cos{\psi }_r\\ -a_{{\psi }_U}{x}_r+v_{M}sin{\psi }_r\\ a_{{\psi }_M}-a_{{\psi }_U}\end{array}\right]$$

(7)

in which $x_r$, $y_r$, and ${\psi }_r$ denote the relative position and heading, respectively, as illustrated in the figure. The proof of Equation (7) can be found in [33], which utilized the knowledge of Lie algebra in $SE\left(2\right)$.

3. Solution Methods

In this section, the level-set method is introduced to solve the reachability problem, starting with the analysis of deterministic scenarios, and then, incorporating uncertainty factors to analyze stochastic scenarios.

3.1. Level-Set Method

Level-set methods are a collection of numerical algorithms for solving a particular class of partial differential equations. Questions of reachability for continuous and hybrid systems can be formulated as optimal control or game theory problems. The safety separation problem addressed in this paper could be formulated as a differential game. Equation (1) can be rewritten in the following Isaacs partial differential equation form:

$$\dot{x}=f(t,x,u_U,u_M),x\in {\mathbb{R}}^n$$

(8)

Fromthe perspective of the UAV, it is necessary to consider the worst-case scenario, where the actions of the UAV drive it away from the set of unsafe states, while the actions of the MAV drive the UAV towards the set of unsafe states. So in Equation (8), the input $u_U$ is considered as the control keeping the system safe, and the input $u_M$ consists of disturbances, which are trying to make the system unsafe.

As has been proved in Figure 1, the viscosity solutions of the Hamilton–Jacobi partial differential equation (HJPDE) can be linked to the results of the reachability problem. Computation of the backwards reachable set or tube is normally encoded as a terminal-value HJPDE [34]. The level-set method assumes that the avoid set $\mathcal{T}\left(t\right)$ for the reachability problem is closed and can be represented as the zero sub-level set of a bounded and Lipschitz continuous function $l:{\mathbb{R}}^n\to \mathbb{R}$:

$$\mathcal{T}\left(t\right)=\{x:l(x)\le 0\}$$

(9)

Therefore, to solve the backward reachable tube problem in Equation (5), it is necessary to solve the following:

$$\begin{array}{cc}\hfill & \underset{u_U\in {\mathcal{U}}_U\left(t\right)}{max}\underset{u_M\in {\mathcal{U}}_M\left(t\right)}{min}\underset{t\in [T,0]}{min}l\left(x\left(t\right)\right)s.t\hfill \\ \hfill & \dot{x}=f(t,x,u_U,u_M)\hfill \\ \hfill & \mathcal{T}\left(t\right)=\{x:l(x)\le 0\}\hfill \end{array}$$

(10)

The max and min in Equation (10) refer to the optimization operation on the Lipschitz continuous function l, because from the perspective of the UAV, it is desirable for $l\left(x\right)$ to be as large as possible so that fewer x satisfy $l\left(x\right)<0$. Conversely, from the perspective of the MAV in the worst situation, it is desirable for $l\left(x\right)$ to be as small as possible so that more x will satisfy $l\left(x\right)<0$. The evolution of the backwards reachable set is accomplished by solving

$$\frac{\partial }{\partial t}\varphi (x,t)+min\left[0,H\left(x,\frac{\mathrm{d}}{\mathrm{d}x}\varphi (x,t)\right)\right]=0$$

(11)

where $\varphi (x,t)$ is an implicit surface representation of the finite time backwards reachable set with a terminal condition

$$\varphi \left(x\left(0\right),0\right)=l\left(x\left(0\right)\right)$$

(12)

H denotes the Hamiltonian, where

$$H(x,p)=\underset{u_U\in {\mathcal{U}}_U\left(t\right)}{max}\underset{u_M\in {\mathcal{U}}_M\left(t\right)}{min}{\mathbf{p}}^{T}f(t,x,u_U,u_M)$$

(13)

and

$${\mathbf{p}}^T=\nabla \varphi (x,t)$$

(14)

The solution approach was proven in [35].

3.2. Deterministic Backward Reachable Analysis

First, the simplified scenario without uncertainty was considered, and the relative backward reachable tube was computed. The relative location between the two aircraft was formulated as a deterministic equation in Equation (7).

3.2.1. Target Set

The target set in reachability analysis represents the terminal situation $l\left(x\left(0\right),0\right)$ in solving the HJPDE. In this case, the target set denotes the unsafe set when an LOS would occur. As described in Equation (3), the target set can be defined as

$$\mathcal{D}\left(t\right)=\{X:\sqrt{(x_r^2+y_r^2)}\le D_{MU}\mathrm{for}\forall {\psi }_r\}$$

(15)

3.2.2. Hamiltonian Function Calculation

Based on Equations (3) and (13), the Hamiltonian function can be derived. Let $x1=x_r,$ $x2=y_r,x3={\psi }_r$, and $u_U=a_{{\psi }_U},u_M=a_{{\psi }_M}$; then, it can be obtained as below:

$$\begin{array}{cc}\hfill H(x,p)& =\underset{u_U\in {\mathcal{U}}_U\left(t\right)}{max}\underset{u_M\in {\mathcal{U}}_M\left(t\right)}{min}{\mathbf{p}}^{T}f(t,x,u_U,u_M)\hfill \\ \hfill & =\underset{a_{{\psi }_U}\in {\mathcal{U}}_U\left(t\right)}{max}\underset{a_{{\psi }_M}\in {\mathcal{U}}_M\left(t\right)}{min}(-p_1{v}_u+p_1{v}_{M}cos{\psi }_r+p_2{v}_{M}sin{\psi }_r+\hfill \\ & p_1{a}_{{\psi }_U}{y}_r-p_2{a}_{{\psi }_U}{x}_r+p_3{a}_{{\psi }_M}-p_3{a}_{{\psi }_U})\hfill \\ \hfill & =-p_1{v}_u+p_1{v}_{M}cos{\psi }_r+p_2{v}_{M}sin{\psi }_r+\hfill \\ & a_{{\psi }_{U}max}|p_1{y}_r-p_2{x}_r-p_3|-a_{{\psi }_{M}max}\left|p_3\right|\hfill \end{array}$$

(16)

in which the subscript max represents the upper bound of the angular velocities of the two aircraft.

3.2.3. Algorithm Steps

The HJPDE in the form of Equation (11) can be solved using the level-set toolbox developed by Mitchell [34]. The Initial-State-Constrained Backward Reachable Tube (ISCBRT) algorithm proposed in this work aims to terminate the computing procedure when the bound of the backward reachable tube arrives at the given initial relative state. Based on the toolbox, the pseudocode of the algorithm is shown below.

The output of Algorithm 1 is the backward reachable tube $\mathcal{T}\left(t\right)$ and the maximum 2-norm value of distance ${max\left|\right|D\left|\right|}_2$. $\mathcal{T}\left(t\right)$ represents the unsafe set of states in which there is a possibility for an LOS event to happen in specific circumstances. $\mathcal{T}\left(t\right)$ represents the minimum safe time for the UAV to escape the LOS event, and the safe terminal time $T_{safe}$ represents the minimum response time of the UAV. So, if the minimum responding time of the UAV is given, the safe separation distance can be determined by repeating Algorithm 1.

Algorithm 1 Calculate the safe LOS resolution time based on the ISCBRT
Input: Initial relative state $X_0$, Unsafe set $\mathcal{D}\left(t\right)$, Safe terminal time $T_{safe}$
1:	Set the computing state space $\mathcal{S}$, grid size N, and time step size $\mathsf{\Delta }t$
2:	Set the relative Hamiltonian function $H(x,p)$ and partial function ${\alpha }_i\left(x\right)=\underset{p}{max}\left|\frac{\partial H\left(p\right)}{\partial p_i}\right|$ for HJPDE computing.
3:	$\mathcal{T}\left(t\right)=\mathcal{D}\left(t\right)$, $t=0$, $i=0$
4:	while $X_0\notin \mathcal{T}\left(t\right)$ or $H(x,p)\ne 0$ and $t<T_{safe}$ do
5:	Solve the HJPDE (11) at the backward computation time $t+=-i\times \mathsf{\Delta }t$
6:	Let $x=0$ to obtain $\mathcal{T}\left(t\right)=\varphi (0,t)$
7:	Let $i+=1$
8:	end while
9:	Calculate the slice isosurfaces of the specific dimensions.
10:	Calculate the maximum 2-norm value of distance ${max\left|\right|D\left|\right|}_2$ among all the isosurfaces.
Output: Backward reachable tube $\mathcal{T}\left(t\right)$, maximum 2-norm value of distance ${max\left|\right|D\left|\right|}_2$.

3.3. Stochastic Backward Reachable Analysis with Location Uncertainty

3.3.1. Solution of the Probabilistic Backwards Reachable Tube

Consider a controlled stochastic process in the form of (2), and let $\rho \in [0,1]$ and $t\le T$. Then, the backwards reachable tube under the probability of success $\rho$ can be denoted as [36]

$$\mathcal{T}{\left(t\right)}^{\rho }=\{x\in {\mathbb{R}}^n:\exists u_M\in {\mathcal{U}}_M\left(t\right),\forall u_U\in {\mathcal{U}}_U\left(t\right),\exists s\in [T,0],\mathbb{P}[\zeta (s;x,t,u_U,u_M)\in \mathcal{D}\left(t\right)]>\rho \}$$

(17)

For the initial condition states in $\mathcal{T}{\left(t\right)}^{\rho }$, the probability that there exists at least one trajectory that reaches set $\mathcal{D}\left(t\right)$ during the time horizon $\left|t\right|$ is $\rho$. To solve the probability problem, the expectation should be computed rather than being an implicit surface in a particular situation. The expected pay-off is optimized over all possible further realizations of the Wiener process [37]:

$$\varphi (x,t)=E\left[\underset{u_U}{sup}\underset{u_M}{inf}\varphi \left(x+\Delta f(x,u_U,u_M)+\xi ,t+\Delta t\right)\right]$$

(18)

In Equation (18), $\xi$ represents a process with Gaussian distribution. Figure 4 illustrates a 10-s standard Brownian process with a standard deviation $\sigma =\sqrt{\Delta t}$.

The stochastic differential game in the form of (2) can be interpreted in Itô’s sense [28].

$$\begin{array}{c}\hfill \varphi \left(x+\mathsf{\Delta }f(x,u_U,u_M)+\xi ,t+\mathsf{\Delta }t\right)\approx \\ \hfill & \varphi (x,t)+{\varphi }_t(x,t)\mathsf{\Delta }t+{\varphi }_x(x,t)\mathsf{\Delta }x+\hfill \\ \hfill & \frac{1}{2}\left({\varphi }_{xx}(x,t)\mathsf{\Delta }{x}^2+{\varphi }_{tt}(x,t)\mathsf{\Delta }{t}^2+2{\varphi }_{xt}(x,t)\mathsf{\Delta }x\mathsf{\Delta }t\right)+\hfill \\ \hfill & O\left({\mathsf{\Delta }}^3\right)\hfill \end{array}$$

(19)

Let $\varphi (x,t)$ be the viscosity solution of the stochastic HJPDE, it can be obtained that

$$\frac{\partial }{\partial t}\varphi (x,t)+min\left[0,H\left(x,\frac{\mathrm{d}}{\mathrm{d}x}\varphi (x,t)\right)\right]+\frac{1}{2}\mathrm{Tr}\left[\sigma (x,t){\sigma }^T(x,t)\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}\varphi (x,t)\right]=0$$

(20)

3.3.2. Indicator Function Design

In the solution of a probabilistic backwards reachable tube problem, the signed distance function rather than a zero-level contour of the value function ${\mathbb{1}}_{\mathbb{D}}$ [36] is utilized, because the maximum objective of the probability is 1. In the context of the avoidance problem, the aim of the solution is to drive the system away from the unsafe set. So, the terminal condition of the LS function is then equal to zero for all states that are inside the initial set $\mathcal{D}$, and all other states are equal to one, as shown in Equation (21).

$${\mathbb{1}}_{\mathbb{D}}\left(x\right)=\left\{\begin{array}{cc}\hfill 0& x\in \mathcal{D}:=\mathbf{Distance}\left(x\right)\ge 0\\ \hfill 1& x\notin \mathcal{D}:=\mathbf{Distance}\left(x\right)<0\end{array}\right.$$

(21)

The term “distance” denotes the signed distance from the unsafe target. As in Ref. [38], to overcome the difficulty cased by the discontinuity of the indicator functions when solving the viscosity solutions of HJPDE, a regularized version of ${\mathbb{1}}_{\mathbb{D}}\left(x\right)$ is considered.

$${\mathbb{1}}_{\mathbb{DR}}\left(x\right)=min(1,max(0,\frac{1}{ϵ}\mathbf{Distance}(x)\ge 0))$$

(22)

Usually, the smoothing factor $ϵ$ is set related to the grid size of the state space. The origin and regularization of the indicator function are shown in Figure 5a,b, respectively.

Combining the values defined in Equations (18) and (22), the probability backwards reachable tube can be approximated by solving Equation (20) in the form of the HJPDE. The errors in approximation were studied in [36].

4. Numerical Simulation Results

In this section, the numerical simulations of the deterministic and stochastic backward reachable tube are described, and the results are discussed.

4.1. Simulation Settings

As discussed in the above section, the key to solving the backward reachable tube problem and implementing the ISCBRT algorithm lies in solving the HJPDE using the level-set method. The simulations analyzed herein were conducted with the level-set method toolbox developed by Mitchell [28]. The correctness and flexibility of the toolbox has been proven in several research works. The main parameters in the level-set method toolbox include the grid structure, estimation time horizon, dynamic system constraints, and the related Hamilton function.

The research objective of this paper is to use Equation (8) for the deterministic environment and Equation (2) for the stochastic environment. The simulation parameters of the dynamic system are shown in

Table 1. Settings of the dynamic system parameters in the simulation.

Parameter	Description	Value
$D_{MU}$	Radius of LOS event	5 m
$V_U$	Constant velocity of UAV	5 m/s
$V_M$	Constant velocity of MAV	20 m/s
${\mathcal{U}}_U\left(t\right)$	Bound of UAV angular velocity	2 rad/s2
${\mathcal{U}}_M\left(t\right)$	Bound of MAV angular velocity	1 rad/s2

.

In Table 1, the flight dynamic differences between the UAV and MAV are considered. The UAV exhibits greater agility but slower speed, whereas the MAV possesses higher velocity but inferior agility. Note that for an MAV flying at a low altitude, especially emerging vehicles such as electric vertical takeoff and landing aircraft (eVTOL), 20 m/s could be a reasonable speed.

4.2. Deterministic Backward Reachable Analysis

The deterministic backward reachable tube was computed based on the parameters in Table 1. In this case, the state space was set as ${\mathbb{R}}^3=[-40,40]\times [-40,40]\times [0,2\pi ]$, in which ${\psi }_r$ is a periodic dimension. The grid size was set to be $\delta x=[0.3,0.3,0.06]$ in three dimensions. According to the ISCBRT algorithm described in Section 3.2.3, the simulation time horizon was set to be 3 s, as the term $\frac{\partial }{\partial t}\varphi (x,t)$ tends to be 0 as the simulation time increases. The computation result is shown in Figure 6. The simulations were conducted using Matlab R2020b with an Intel i5-13500H CPU @ 2.60 GHz and 32 GB RAM on a Windows 11 operating system.

In Figure 6, the red part denotes the backward reachable tube, and the green part represents the unsafe set in which an LOS event may occur. This implies that if the initial relative states of the two aircraft are within this red isosurface, there is a potential for entering an unsafe set within the estimated 3-s time horizon. The two-dimensional projections of the isosurface in the $x-{\psi }_r$ and $y-{\psi }_r$ directions are shown in Figure 7. It should be noted that the green part, which denotes the unsafe part, is not drawn in the projection.

4.2.1. Influences of the Heading Angle

It is noted that the tube differs depending on the initial heading angle. Figure 8a–d illustrate the slices of the heading angles of 0, ${90}^{\circ }$, ${180}^{\circ }$, and ${270}^{\circ }$. It can be seen that the tubes are extended from the unsafe set, and the extent is in the opposite direction to the initial heading angle. For example, in Figure 8a, the backwards reachable tube is extended in the negative x direction, which means that when two aircraft are heading in the same direction, if the UAV is in front of the manned aircraft ($x_r<0$), it is very likely to enter an unsafe set.

In the numerical computation, limited by the precision of the grid division, not every heading angle can obtain a closed slice surface. Figure 9 presents closed slice surfaces for all valid heading angles within the range of 0 to ${360}^{\circ }$.

It can be observed that the slices of the backwards reachable tube rotate rotate as the heading angle turns. To determine the safety separation of the MAV and UAV, the maximum safety distances of each heading angle are calculated, that is, the maximum 2-norm value of the slice surface boundaries. The minimum safety distances related to the heading angles are shown in Figure 10.

The maximum unsafe distances varied with the heading angles. The largest value was obtained at ${\psi }_r={180}^{\circ }$, with a value of $D_{min}=max{D}_s=26.7$ m. This implies that the worst-case scenario would be two aircraft flying towards each other. The safety separation should be beyond the maximum unsafe distances among all the heading angles; so, it can be concluded that the safety separation is $26.7$ m given the minimum response time $T_{safe}=3.0$ s.

4.2.2. Influences of the Safe Time

Under the premise of the research conclusions from the previous section, the effect of the minimum response time on the minimum safe distance was determined. By repeating the computing process in Section 4.2.1, the quantitative relationship was concluded, as shown in Figure 11.

The conclusions in Figure 11 clearly prove that, as the response speed increases, the minimum safe distance or separation between the two aircraft continues to decrease.

4.3. Stochastic Backward Reachable Analysis

In this section, the uncertainty of the positions and heading angles of the MAV is discussed. Whether in centralized or decentralized air traffic management systems [39], there is always uncertainty in the situational awareness of a UAV relative to an MAV. For example, in a centralized scenario, errors might stem from system latency, while in the decentralized scenario, errors could originate from onboard sensors.

4.3.1. Analysis with Position Uncertainty

It is assumed in this work that the related position in the x and y directions follows a Gaussian distribution. It can be described as

$$\begin{array}{c}\hfill \widetilde{x_r}\left(t\right)=x_r\left(t\right)+\mathsf{\Delta }{x}_r\left(t\right)\\ \hfill \widetilde{y_r}\left(t\right)=y_r\left(t\right)+\mathsf{\Delta }{y}_r\left(t\right)\end{array}$$

(23)

in which $\mathsf{\Delta }{x}_r\left(t\right)\sim \mathcal{N}(0,{\sigma }_x^2(x_r,t)),\mathsf{\Delta }{y}_r\left(t\right)\sim \mathcal{N}(0,{\sigma }_y^2(y_r,t))$. So, combined with Equation (7), Equation (2) can be rewritten as Equation (24).

$$\left[\begin{array}{c}\dot{x_r}\\ \dot{y_r}\\ \dot{{\psi }_r}\end{array}\right]=\left[\begin{array}{c}{a}_{{\psi }_U}{y}_r-v_U+v_{M}cos{\psi }_r\\ -a_{{\psi }_U}{x}_r+v_{M}sin{\psi }_r\\ a_{{\psi }_M}-a_{{\psi }_U}\end{array}\right]+\left[\begin{array}{c}{\sigma }_x(x_r,t)\\ {\sigma }_y(y_r,t)\\ 0\end{array}\right]dB\left(t\right)$$

(24)

The simulation parameters were set, as shown below in

Table 2. Settings of the uncertain parameters in the simulation.

Parameter	Value
${\sigma }_x(x_r,t)$	1 m/$\sqrt{\mathrm{s}}$
${\sigma }_y(y_r,t)$	1 m/$\sqrt{\mathrm{s}}$

, while the dynamic system parameters were the same as those in Table 1. The indicator function was designed as in Equation (22), and the smoothing parameters were $ϵ=\frac{1}{\delta x}$, in which $\delta x$ denotes the grid size of the x dimension.

The probability of the backwards reachable tube of $P=0.6$ is shown in Figure 12. The isosurfaces indicate that, starting from the initial states within this set, there is a $0.8$ probability of inevitably entering the unsafe state set within 1 s. However, it should be noted that, unlike the results under deterministic scenarios, initial states outside of this set may also have the potential to enter the unsafe set.

Similarly, the 2D projections of the isosurface in the $x-{\psi }_r$ and $y-{\psi }_r$ directions are shown in Figure 13.

It can be seen that compared to the deterministic scenarios, the shape of the backwards reachable tube under uncertain conditions is more curved and narrow. Moreover, the extent of its influence is broader, due to the positional uncertainties that increase the likelihood of the tube extending in various directions. Therefore, to conduct the computation in the limited state space, the simulation time was limited to 1 s, which was less than that in the deterministic backward reachable analysis.

Figure 14a–d show the slices with different probabilities of different heading angles. In [20], the risk domain at the risk level $\alpha$ of a random variable X was defined as a set $\mathcal{D}\in {\mathbb{R}}^d$ that satisfies

$$\mathrm{Pr}(X\in \mathcal{D})\ge 1-\alpha$$

(25)

Therefore, in Figure 14a–d, the risk domain is inside the enclosed area. As the level of risk increases, the area enclosed by the tube’s envelope decreases. In the numerical computations, the value of the value function was obtained by summing the deterministic part represented by the first two terms with the probabilistic part represented by the third term in Equation (20). So, the impact of the probability on the envelope was primarily derived from the third term.

The minimum safe distances of different risk levels related to the position uncertainty are shown in Figure 15. It can be concluded that the minimum safe separation increases with the increase in the risk level. Specifically, the largest safe separation is $30.0$ m given the minimum response time $T_{safe}=1.0$ s under a risk level of $\alpha =0.1$.

4.3.2. Analysis with Heading Angle Uncertainty

In this section, the relative heading angle ${\psi }_r$ follows a Gaussian distribution. It can be described as

$$\left[\begin{array}{c}\dot{x_r}\\ \dot{y_r}\\ \dot{{\psi }_r}\end{array}\right]=\left[\begin{array}{c}{a}_{{\psi }_U}{y}_r-v_U+v_{M}cos{\psi }_r\\ -a_{{\psi }_U}{x}_r+v_{M}sin{\psi }_r\\ a_{{\psi }_M}-a_{{\psi }_U}\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ {\sigma }_{{\psi }_r}(x_r,t)\end{array}\right]dB\left(t\right)$$

(26)

In this simulation, the standard variance in the relative angle was set as ${\sigma }_{{\psi }_r}=0.3$ rad/$\sqrt{\mathrm{s}}$. The result of the stochastic backwards reachable tube with a probability of $P=0.8$ is given below in Figure 16, and the projection figures are given in Figure 17a,b.

Similarly, the envelope of the tube varied with the initial relative heading angles. The slices of the tube with different risk levels are shown in Figure 18a–d.

It can be observed that the overall shape of the result tube is consistent, and the trend of expansion is similar to that discussed in the previous section. As the level of risk increases, the envelope of the set continues to expand. This is because the influence of the uncertain terms is relatively minor compared to the deterministic terms.

The minimum safe distances of different risk levels related to the heading angle uncertainty are shown in Figure 19. In this scenario, the minimum safe separation between two aircraft also occurs when the relative heading angle is ${180}^{\circ }$. When the risk level is at its minimum of $\alpha =0.1$, the minimum separation is 29.7 m.

5. Discussion and Future Work

5.1. Discussion

Based on the results obtained from the numerical computations described in the previous section, a comparison can be made between the minimum separations under deterministic and uncertain scenarios, as shown in

Table 3. Comparison of the minimum safety separation in deterministic and uncertain scenarios, given the response time of 1 s.

Scenario	Minimum Safety Separation
Deterministic	16.1 m
${\sigma }_x={\sigma }_x=1\mathrm{m}/\sqrt{\mathrm{s}},\alpha =0.1$	30.0 m
${\sigma }_x={\sigma }_x=1\mathrm{m}/\sqrt{\mathrm{s}},\alpha =0.2$	28.8 m
${\sigma }_x={\sigma }_x=1\mathrm{m}/\sqrt{\mathrm{s}},\alpha =0.3$	28.1 m
${\sigma }_x={\sigma }_x=1\mathrm{m}/\sqrt{\mathrm{s}},\alpha =0.4$	27.4 m
${\sigma }_{{\psi }_r}=0.3\mathrm{rad}/\sqrt{\mathrm{s}},\alpha =0.1$	29.7 m
${\sigma }_{{\psi }_r}=0.3\mathrm{rad}/\sqrt{\mathrm{s}},\alpha =0.2$	28.4 m
${\sigma }_{{\psi }_r}=0.3\mathrm{rad}/\sqrt{\mathrm{s}},\alpha =0.3$	27.5 m
${\sigma }_{{\psi }_r}=0.3\mathrm{rad}/\sqrt{\mathrm{s}},\alpha =0.4$	26.7 m

. The data for the deterministic scenario were extracted from Figure 11.

The consideration of uncertainty increases the minimum separation because measurement errors in the position or heading angle can potentially lead to initial states that are not within the unsafe set, guiding the two aircraft into a unsafe set. From the results, the overall trends across the three scenarios are similar, with the differences in details being less than $10\%$.

From the perspective of mathematics, the computation of the uncertainty reachable tube consists of two parts: the calculation of the drift term and the diffusion term, which represents the deterministic component; and the probabilistic component of the system’s evolution, respectively. Since the deterministic equations are the same, the drift terms across the three scenarios are identical, resulting in similar shapes of the reachable sets. The differences are primarily determined by the diffusion terms, which are solutions based on Gaussian distributions; thus, different probabilities lead to disparities in the outcomes.

From a physical perspective, as airborne sensors advance, the errors in the distance and heading angle become less significant compared to the deterministic kinetic characteristics of the two aircraft; hence, they do not exert a substantial impact on the overall trend.

Furthermore, it is important to note that this study is based on the worst-case scenario from the perspective of a UAV, which is assumed to avoid unsafe events as much as possible, whereas the opposite MAV is considered to have a tendency towards unsafe events. However, in actual operations, the interaction between two aircraft is not a zero-sum game. Based on situational awareness, all participants in air traffic will strive to prevent mid-air conflicts. Therefore, the results of this study are conservative.

5.2. Limitations and Future Work

Considering that this research was entirely based on numerical simulations, as future work, a hardware-based [40] actual flight test could be conducted to validate the conclusions of this work. Another limitation of this work is that both obstacles and wind are valid disruptions in manned and unmanned aircraft scenarios. In the future, obstacles and wind will be incorporated as certain or stochastic disturbance terms in the HJPDE into the reachability analysis and separation estimation research.

6. Conclusions

In this study, the problem of calculating the reachable tube of unsafe events during the integrated operation of an MAV and a UAV was established as an optimization problem of the viscosity solutions of the HJPDE. Then, a backward calculation based on the Hamilton–Jacobi partial differential equation was derived to obtain the initial unsafe set that the UAV needs to avoid, given the response events and unsafe sets. In addition, the level-set method was utilized to solve the value function of the optimal problem.

Secondly, the problem of calculating the initial unsafe set under both deterministic and uncertain scenarios was considered. Specifically, in uncertain scenarios, the relative positions and relative heading angles of the two aircraft were modeled using a Brownian model or a Gaussian distribution model, respectively. By computing the expectation of the value function in the form of an HJPDE, the unsafe envelopes for a specific risk level were obtained. Therefore, a conservatively biased minimum separation distance was obtained, leading to the quantification of the safety risk assessment results. The result shows that with the increase in the risk level, the minimum safe separation increases.

Finally, the calculation results under different uncertainties and risk levels were discussed, explaining the rationality and similarity of the backward reachable set from both mathematical and physical perspectives. The next steps of this work primarily include considering the calculation of the minimum separation distance for the integrated operation of manned and unmanned aerial vehicles in cooperative scenarios, as well as the design of UAV maneuver control laws under different risk levels.