\section{Introduction}

To describe the behavior of linear control systems, a number of strict definitions have been proposed,
such as controllability, observability, etc., and many mature analytical methods have been proposed for analyzing
these properties \cite{b1}. These methods have been widely applied in engineering for solving stability and reliability problems.
However, in nonlinear control systems, properties similar to those mentioned above are difficult to be defined and analyzed \cite{b2}.

Reachability, viability and invariance analyses are important methods to study the behavior of nonlinear control systems.
In reachability, viability and invariance problems, one specifies a target set in the state
space and a time horizon and then seeks to find the reachable, viable and invariant sets of the target set \cite{a001}.
However, computing these sets is a challenging task. Due to the diversity of system states and control inputs,
verifying each state one by one always consumes a lot of computation time \cite{a002,i1}.
This makes the approaches based on simulation insufficient, and thus
formal verification methods are needed. Unfortunately, except for some solvable problems \cite{a003},
to exactly compute these sets is generally intractable. Therefore,
the common practice is to compute their approximate expressions.
For the computation of reachable set,
numerous techniques have been proposed in the past decades, which can be divided into two categories:
the Lagrangian methods \cite{a004} and the methods based on state space discretization.
Lagrangian methods mainly include ellipsoidal methods \cite{i3,a006,i4}, polyhedral methods \cite{i2,a005,a007} and support
vector machines \cite{a004},
which are able to solve
high-dimensional problems but have strict requirements on the form of dynamic systems and are primarily used to solve linear problems.
The main advantage of the state space discretization-based approach is
its less stringent requirements on the form of the dynamical system,
and it can solve nonlinear problems.
This category of methods mainly includes level set method \cite{a001,a002,a008} and distance fields over grids (DFOG)
method \cite{a0010,a0011,a0012},
among which the level set method is more widely used in engineering.

More importantly, the level set method creatively characterizes the reachable, viable and invariant sets in a unified form.
In level set method, the Hamilton-Jacobi (HJ) partial differential equation (PDE) about a value function is numerically solved,
and the reachable, viable and invariant sets can then be expressed as zero level sets of solutions of the aforementioned equation \cite{a001}.
Different types of sets can be obtained by setting different terminal conditions and switching the MIN or MAX
operator during the computation.
Based on this principle, some mature toolboxes have been developed and have greatly promoted the development of
computational techniques for reachable \cite{tr1,web1,web2,web3},
viable and invariant sets and solved many practical engineering problems,
such as flight control systems \cite{i5,a009,a0013}, ground traffic management systems \cite{i6,a0014}, air
traffic management systems \cite{a0015}, etc.

Although the level set method has developed into a mature, reliable and commonly used tool over the years, it still has some shortcomings.
\begin{itemize}
	\item Level set method involves solving the viscous solution of a HJ PDE, during this process, an optimization problem
	referred to as dissipation function needs to be solved \cite{tr1}. To our best knowledge, although this optimization problem can be
	reduced in affine nonlinear systems, there is no efficient method to solve it in more general systems.
	Therefore, the application of level set methods is mainly limited to affine nonlinear systems.
	\item In level set method, the terminal condition of the HJ PDE is always set as a signed
	distance function of the target set. However, it is difficult to construct such a function for an irregular target set.
	In some studies, the irregular target set is replaced by a rectangle when construct the required signed distance function \cite{a009,a0013}.
\end{itemize}
To overcome these shortcomings, this paper proposes a novel unified framework for solving reachability,
viability, and invariance problems. In this framework, the reachable, viable and invariant sets
can be described as some level sets (non-zero level sets) of different value functions.
This approach differs from the level set approach in that the value functions in this framework are obtained
not by solving the HJ PDE but by recursion and interpolation, and such a mechanism has the following advantages:
\begin{itemize}
	\item Since there is no PDE solution involved, the computation of dissipation function is effectively avoided.
	Therefore, the proposed method can be applied to non-affine nonlinear systems.
	\item As it is unnecessary to construct a signed distance function of the target set as the terminal condition of the HJ PDE,
	this method can easily handle irregular target sets.
\end{itemize}

The structure of this paper is as follows. Section II
introduces the definitions of reachable, viable and invariant sets and briefly describes the level set method.
Section III describes the proposed method in detail.
Some numerical examples are given in Section IV and Section V to illustrate the
validity of the proposed method. The results are summarized
in Section VI.

\section{Background}
\subsection{Definitions of reachable, viable and invariant sets}
Consider a continuous time control system with fully observable state:
\begin{align}
	\label{sys1}
	\dot{s}=f(s,u)
\end{align}
where $s\in \mathbb{R}^n$ is the system state, $u\in\mathcal{U}$ is the control input ($\mathcal{U}$ is a closed set).
The function $f(.,.):\mathbb{R}^n\times \mathcal{U}\to \mathbb{R}^n$ is Lipschitz continuous and bounded.
Let $\mathscr{U}$ denote the set of Lebesgue measurable functions from the time interval $[0,\infty)$ to $\mathcal{U}$.
Then, given the initial state $s_0$ at time $t_0$ and $u(.)\in \mathscr{U}$, the evolution of system (\ref{sys1}) in time
interval $[t_0,t_1]$ can be expressed as a continuous trajectory $\phi_{t_0}^{t_1}(.,s_0,u):[t_0,t_1]\to \mathbb{R}^n$
and $\phi_{t_0}^{t_1}(t_0,s_0,u)=s_0$. Given a target set $K$ and a time horizon $T$,
literature \cite{a001} raises three questions of interest:
\begin{itemize}
	\item \textbf{Reachability} Does there exist a state $s_0$ and a $u(.)\in \mathscr{U}$ and a $t\in [0,T]$ such that the trajectory
	satisfies $\phi_{0}^{T}(t,s_0,u) \in K$?
	\item \textbf{Viability} Does there exist a $s_0$ and a $u(.)\in \mathscr{U}$ for which the trajectory $\phi_{0}^{T}(.,s_0,u)$
	satisfies $\phi_{0}^{T}(t,s_0,u) \in K$ for any $t\in [0,T]$?
	\item \textbf{Invariance} Does there exist a $s_0$ such that the trajectories $\phi_{0}^{T}(.,s_0,u)$ for all $u(.)\in \mathscr{U}$
	satisfies $\phi_{0}^{T}(t,s_0,u)\in K$ for any $t\in [0,T]$?
\end{itemize}
One would like to characterize some sets of initial states at time $0$ to answer yes to these questions. These sets
are referred to as reachable, viable and invariant sets, respectively. Formally, the reachable,
viable and invariant sets can be expressed as the following definitions:
\begin{definition}
	\textbf{Maximal reachable set}
	\begin{align}
		\begin{split}
		&\mathcal{R}_{\mathrm{max}}(K,T)\\
		&=\left\{s_0| \exists u(.)\in\mathscr{U}\ \exists t\in [0,T]\ \phi_{0}^{T}(t,s_0,u)\in K  \right\}
		\end{split}
	\end{align}
\end{definition}

\begin{definition}
	\textbf{Viable set}
	\begin{align}
		\begin{split}
		&\mathcal{V}(K,T)\\
		&=\left\{s_0| \exists u(.)\in\mathscr{U}\ \forall t\in [0,T]\ \phi_{0}^{T}(t,s_0,u)\in K  \right\}
		\end{split}
	\end{align}
\end{definition}

\begin{definition}
	\textbf{Invariant set}
	\begin{align}
		\begin{split}
		&\mathcal{I}(K,T)\\
		&=\left\{s_0| \forall u(.)\in\mathscr{U}\ \forall t\in [0,T]\ \phi_{0}^{T}(t,s_0,u)\in K  \right\}
		\end{split}
	\end{align}
\end{definition}

\subsection{Level set method}
In level set method, the following HJ PDEs about $V_1$ and $V_2$ are numerically solved:
\begin{align}
	\begin{split}
		&\left\{\begin{array}{l}
		\displaystyle{ \frac{\partial V_1}{\partial t}(s,t)+\min \left[0,\max_{u\in \mathcal{U}} \frac{\partial V_1}{\partial s}(s,t) f(s,u) \right]=0}\\
		\text{s.t.  } V_1(s,T)=l(s)
		\end{array}\right.\\
		&\left\{\begin{array}{l}
			\displaystyle{ \frac{\partial V_2}{\partial t}(s,t)+\min \left[0,\min_{u\in \mathcal{U}} \frac{\partial V_2}{\partial s}(s,t) f(s,u) \right]=0}\\
			\text{s.t.  } V_2(s,T)=l(s)
		\end{array}\right.
	\end{split}
\end{align}
where $l(.):\mathbb{R}^n \to \mathbb{R}$ is a Lipschitz continuous function and satisfies:
\begin{align}
	K=\left\{s\in\mathbb{R}^n|l(s)\geq 0 \right\}
\end{align}
A common choice for the function $l(.)$ is the signed distance to the target set $K$.
Then, the viable and invariant sets are denoted as zero level sets of two value functions, i.e.:
\begin{align}
	\begin{split}
		\mathcal{V}_K(T)=\left\{s_0| V_1(s_0,0) \geq 0 \right\}\\
		\mathcal{I}_K(T)=\left\{s_0| V_2(s_0,0) \geq 0 \right\}
	\end{split}
\end{align}
The analytical form of the solution of the aforementioned equations is difficult to obtain,
therefore, the two value functions are characterized as viscosity solutions to these PDEs.
Literature \cite{a001} points out that the invariance and reachability problems
are duals of one another and
\begin{align}
	\label{rlt1}
	\mathcal{R}_{\text{max}}(K,T)=\complement_{\mathbb{R}^n} \left[\mathcal{I} ({\complement_{\mathbb{R}^n } K },T)\right]
\end{align}
Thus, the level set method represents the reachable, viable and invariant sets in a unified form.

To compute these sets by using the level set toolboxes, the user needs to specify the following parameters:
\begin{itemize}
	\item \textbf{Computational domain}: a cubic region in the state space.
	\item \textbf{Grids number}: the number of grid nodes in each dimension of the computational domain.
	\item \textbf{Terminal condition of the HJ PDE}: the signed distance to the target set.
	\item \textbf{Hamiltonian}:
	\begin{align}
		\begin{split}
		&\mathcal{H}(s,p_1)=\max_{u\in\mathcal{U}} \left[p_1^\mathrm{T} f(s,u)  \right] \text{ or }\\
		&\mathcal{H}(s,p_2)=\min_{u\in\mathcal{U}} \left[p_2^\mathrm{T} f(s,u)  \right]
		\end{split}
	\end{align}
	where $\displaystyle{p_1= \frac{\partial V_1}{\partial s}(s,t)},\displaystyle{p_2= \frac{\partial V_2}{\partial s}(s,t)}$.
	\item \textbf{Dissipation function}:
	\begin{align}
		\begin{split}
		&\mathcal{P}_i(s)=\max_{p_1} \left| \frac{\partial}{\partial p_1^i}\mathcal{H}(s,p_1) \right|
		\text{ or }\\
		&\mathcal{P}_i(s)=\max_{p_2} \left| \frac{\partial}{\partial p_2^i} \mathcal{H}(s,p_2) \right|
		\end{split}
	\end{align}
	where $p_1^i$ and $p_2^i$ are the $i$th component of $p_1$ and $p_2$, respectively.
\end{itemize}

If $\mathcal{U}=[u^1_{\text{min}},u^1_{\text{max}}]\times [u^2_{\text{min}},u^2_{\text{max}}]\times ...$ and
the nonlinear system (\ref{sys1}) is an affine nonlinear system, i.e. the nonlinear system can be separated
into the following form:
\begin{align}
	f(s,u)=f_s(s)+g(s)u
\end{align}
Then the Hamiltonians can be simply written as:
\begin{align}
	\begin{split}
		&\mathcal{H}(s,p_1)=p_1^\mathrm{T}f_s(s)+\sum_i \max\left[p_1^i g^i(s) u^i_{\text{min}},p_1^i g^i(s) u^i_{\text{max}} \right]\\
		&\text{or}\\
		&\mathcal{H}(s,p_2)=p_2^\mathrm{T}f_s(s)+\sum_i \min\left[p_2^i g^i(s) u^i_{\text{min}},p_2^i g^i(s) u^i_{\text{max}} \right]
	\end{split}
\end{align}
and the dissipation function can be overapproximated by the following equations:
\begin{align}
	\mathcal{P}_i=|f^i_s(s)|+\max\left[ |g^i(s)u^i_{\text{min}}|,|g^i(s)u^i_{\text{max}}|   \right]
\end{align}
where $f^i_s(.)$ and $g^i(.)$ are the $i$th component of $f_s(.)$ and $g(.)$, respectively.

However, for a non-affine nonlinear system, it is difficult to construct Hamiltonians and dissipation function
in the above forms. In particular, for dissipation function, there is no effective approach to compute it efficiently.

\section{Recursive Value Function Method}
In this section, we propose a method that can avoid the above-mentioned difficulties,
namely recursive value function method.
We first introduce a definition, minimal reachable set, before formally presenting this method.
\begin{definition}
	\textbf{Minimal reachable set}
	\begin{align}
		\begin{split}
		&\mathcal{R}_{\mathrm{min}}(K,T)\\
		&=\left\{s_0| \forall u(.)\in\mathscr{U}\ \exists t\in [0,T]\ \phi_{0}^{T}(t,s_0,u)\in K  \right\}
		\end{split}
	\end{align}
\end{definition}
Informally, $\mathcal{R}_{\mathrm{min}}(K,T)$ is a set of states from which trajectories can reach the target set
for any $u(.)\in \mathscr{U}$ within time $T$.
Similar to the relationship between the maximal reachable set and the invariant set,
the viable set and the minimal reachable set are also dual to each other, i.e.:
\begin{align}
	\label{rlt2}
	\mathcal{V}(K,T)=\complement_{\mathbb{R}^n} \left[\mathcal{R}_{\mathrm{min}} \left(\complement_{\mathbb{R}^n} K,T\right)\right]
\end{align}
Therefore, the invariant set can be characterized by computing the maximum reachable set,
while the viable set can be represented by computing the minimum reachable set.
To avoid repetitive descriptions, this paper focuses on the computation of the maximal and minimal reachable sets.

\subsection{Characterizations of minimal and maximal reachable sets}
Consider a case where the control input $u(.)$ is chosen to avoid or delay as much as possible
the entry of the system state into the target set.
If a trajectory $\phi_{0}^{T}(.,s_0,u)$ can still enter the target set within the time horizon $T$ in this case,
the initial state of this trajectory $s_0$ must belong to the minimum reachable set.
Thus, a value function $W_1(.):\mathbb{R}^n \to \mathbb{R}$ can be constructed as Equ. (\ref{maxvalue}).

\begin{align}
	\label{maxvalue}
	W_1(s_0)=\begin{cases}
		\left\{\begin{array}{rl}
			\displaystyle{\max_{u(.)}} &t_f\\
			\text{s.t.}& \dot{s}(t)=f(s(t),u(t))\ \forall t\in [0,t_f]\\
				 \ &s(0)=s_0\\
				 \ &u(t)\in\mathcal{U} \ \forall t\in [0,t_f]\\
				 \ &s(t_f)\in K
			\end{array}\right\},\\ \quad \quad \  \text{if } \forall u(.)\in \mathscr{U}\ \exists t\in [0,\infty)\ \phi_{0}^{\infty}(t,s_0,u)\in K\\
		\infty, \quad \text{otherwise}
	\end{cases}
\end{align}

Then the minimal reachable set can be characterized by:
\begin{align}
	\mathcal{R}_{\mathrm{min}}(K,T)=\left\{s_0|W_1(s_0)\leq T  \right\}
\end{align}

Consider another case where the control input $u(.)$ is chosen to enable the system state
to reach the target set in as short a time as possible.
If a trajectory can enter the target set within time horizon $T$ in this case,
then the initial state of this trajectory must belong to the maximal reachable set.
Thus, a value function $W_2(.):\mathbb{R}^n \to \mathbb{R}$ can be constructed as Equ. (\ref{minvalue}).

\begin{align}
		\label{minvalue}
		W_2(s_0)=\begin{cases}
			\left\{\begin{array}{rl}
				\displaystyle{\min_{u(.)}} &t_f\\
				\text{s.t.}& \dot{s}(t)=f(s(t),u(t))\ \forall t\in [0,t_f]\\
				 \ &s(0)=s_0\\
				 \ &u(t)\in\mathcal{U} \ \forall t\in [0,t_f]\\
				 \ &s(t_f)\in K
			\end{array}\right\},\\ \quad \quad \ \ \text{if } \exists u(.)\in \mathscr{U}\ \exists t\in [0,\infty)\ \phi_{0}^{\infty}(t,s_0,u)\in K\\
		\infty, \quad \text{otherwise}
		\end{cases}
\end{align}

Then the maximal reachable set can be characterized by:
\begin{align}
	\mathcal{R}_{\mathrm{max}}(K,T)=\left\{s_0|W_2(s_0)\leq T  \right\}
\end{align}

The optimal control problems in Equ. (\ref{maxvalue}) and Equ. (\ref{minvalue})
are optimal control problems with terminal state constraint.
In order to convert them into optimal control problems without terminal state constraints
to facilitate the construction of the recursive formulas of the value functions,
we shall need to define a modified dynamic system:
\begin{align}
	\label{sys2}
	\dot{s}= \bar{f}(s,u)=\begin{cases}
		f(s,u), &s\notin K\\
		\mathbf{0}, &s\in K
	\end{cases}
\end{align}
and a modified running cost function:
\begin{align}
	\label{runningcost}
	\bar{c}(s) =\begin{cases}
		1, &s\notin K\\
		0, &s\in K
	\end{cases}
\end{align}
Given the state $s_0$ at time $t_0$ and $u(.)\in \mathscr{U}$, the evolution of the modified system (\ref{sys2}) in
time interval $[t_0,t_1]$
can also be expressed as a continuous trajectory $\bar{\phi}_{t_0}^{t_1}(.,s_0,u):[t_0,t_1] \to \mathbb{R}^n$.

\begin{remark}
	The trajectory $\bar{\phi}_{t_0}^{t_1}(.,s_0,u)$ is the same as trajectory $\phi_{t_0}^{t_1}(.,s_0,u)$
	as long as it evolves outside the target set $K$. When the trajectory $\bar{\phi}_{t_0}^{t_1}(.,s_0,u)$
	touches the border of $K$, then it stays on the border under the dynamics (\ref{sys2}).
\end{remark}

Then, two modified value functions can be defined via optimal control problems without terminal state constraints:
\begin{align}
	\overline{W}_1^{\bar{T}}(s_0)=\left\{\begin{array}{rl}
	\displaystyle{\max_{u(.)}} & \displaystyle{ \int_0^{\bar{T}} \bar{c}(s(t)) \ dt}\\
	\text{s.t.} & \dot{s}(t)=\bar{f}(s(t),u(t)) \ \forall t \in [0,{\bar{T}}] \\
	\   & s(0)=s_0 \\
	\   & u(t)\in\mathcal{U} \ \forall t\in [0,{\bar{T}}]
	\end{array}\right.
\end{align}

\begin{align}
	\overline{W}_2^{\bar{T}}(s_0)=\left\{\begin{array}{rl}
	\displaystyle{\min_{u(.)}} & \displaystyle{ \int_0^{\bar{T}} \bar{c}(s(t)) \ dt}\\
	\text{s.t.} & \dot{s}(t)=\bar{f}(s(t),u(t)) \ \forall t \in [0,{\bar{T}}] \\
	\   & s(0)=s_0 \\
	\   & u(t)\in\mathcal{U} \ \forall t\in [0,{\bar{T}}]
	\end{array}\right.
\end{align}

Then, we may deduce the following results:
\begin{thm}
	For any ${\tau}<\bar{T}$ and any $s\in\left\{s_0| \overline{W}^{\bar{T}}_1(s_0)\leq {\tau} \right\}$,
	we have $\overline{W}^{\bar{T}}_1(s)=W_1(s)$.
\end{thm}

\begin{prf}

	It is obvious that the following equation holds for any $s\in \mathbb{R}^n$:
	\begin{align}
		\overline{W}^{\bar{T}}_1(s)\leq \int_0^{\bar{T}} \max_{s'\in\mathbb{R}^n} \bar{c}(s')\ dt=\int_0^{\bar{T}} 1 dt={\bar{T}}
	\end{align}
	In addition, $\overline{W}^{\bar{T}}_1(s_0)=\bar{T}$ if and only if there exists a $u(.)\in\mathscr{U}$
	such that the trajectory $\bar{\phi}_0^{\bar{T}}(.,s_{0},u)$ cannot reach the target set $K$
	within time $\bar{T}$. Therefore $\overline{W}^{\bar{T}}_1(s_0)\leq {\tau}<\bar{T}$ means that
	for any $u(.)\in \mathscr{U}$, the trajectory $\bar{\phi}_{0}^{\bar{T}}(.,s_0,u)$
	can reach the target set $K$ at a certain time $t_f\in [0,{\bar{T}}]$, i.e. $\bar{\phi}_{0}^{\bar{T}}(t_f,s_{0},u)\in K$.
	In this case, we have:
	\begin{align}
		\begin{split}
		\overline{W}^{\bar{T}}_1(s_0)&=\max_{u(.)\in \mathscr{U}} \int_0^{\bar{T}} \bar{c}(\bar{\phi}_{0}^{\bar{T}}(t,s_{0},u)) dt\\
		&=\left\{\begin{array}{rl}
			\displaystyle{\max_{u(.)\in \mathscr{U}}} &\displaystyle{ \left[ \int_0^{t_f} 1 dt + \int^{\bar{T}}_{t_f} 0 dt \right]}\\
			\text{s.t.} & \bar{\phi}_{0}^{\bar{T}}(t_f,s_{0},u)\in K
		\end{array}\right. \\
		&=\left\{\begin{array}{rl}
			\displaystyle{\max_{u(.)\in \mathscr{U}}} &t_f\\
			\text{s.t.} & \phi_{0}^{\bar{T}}(t_f,s_{0},u)\in K
		\end{array}\right.\\
		&=W_1(s_0) \ \ \ \ \ \square
		\end{split}
	\end{align}
\end{prf}

\begin{thm}
	For any ${\tau}<\bar{T}$ and any $s\in\left\{s_0| \overline{W}^{\bar{T}}_2(s_0)\leq {\tau} \right\}$,
	we have $\overline{W}^{\bar{T}}_2(s)=W_2(s)$.
\end{thm}

\begin{prf}

	It is obvious that the following equation holds for any $s\in \mathbb{R}^n$:
	\begin{align}
		\overline{W}^{\bar{T}}_2(s)\leq \int_0^{\bar{T}} \max_{s'\in\mathbb{R}^n} \bar{c}(s')\ dt=\int_0^{\bar{T}} 1 dt={\bar{T}}
	\end{align}
	In addition, $\overline{W}^{\bar{T}}_2(s_0)=\bar{T}$ if and only if for any $u(.)\in\mathscr{U}$,
	the trajectory $\bar{\phi}_0^{\bar{T}}(.,s_{0},u)$ cannot reach the target set $K$
	within time $\bar{T}$. Therefore $\overline{W}^{\bar{T}}_2(s_0)\leq {\tau}<\bar{T}$ means that
	there exists a $u(.)\in \mathscr{U}$, such that the trajectory $\bar{\phi}_{0}^{\bar{T}}(.,s_0,u)$
	can reach the target set $K$ at a certain time $t_f\in [0,{\bar{T}}]$, i.e. $\bar{\phi}_{0}^{\bar{T}}(t_f,s_{0},u)\in K$.
	In this case, we have:
	\begin{align}
		\begin{split}
		\overline{W}^{\bar{T}}_2(s_0)&=\max_{u(.)\in \mathscr{U}} \int_0^{\bar{T}} \bar{c}(\bar{\phi}_{0}^{\bar{T}}(t,s_{0},u)) dt\\
		&=\left\{\begin{array}{rl}
			\displaystyle{\max_{u(.)\in \mathscr{U}}} &\displaystyle{ \left[ \int_0^{t_f} 1 dt + \int^{\bar{T}}_{t_f} 0 dt \right]}\\
			\text{s.t.} & \bar{\phi}_{0}^{\bar{T}}(t_f,s_{0},u)\in K
		\end{array}\right. \\
		&=\left\{\begin{array}{rl}
			\displaystyle{\max_{u(.)\in \mathscr{U}}} &t_f\\
			\text{s.t.} & \phi_{0}^{\bar{T}}(t_f,s_{0},u)\in K
		\end{array}\right.\\
		&=W_1(s_0) \ \ \ \ \ \square
		\end{split}
	\end{align}
\end{prf}
Theorem 1 indicates that, given the modified value functions $\overline{W}^{\bar{T}}_1(.)$,
due to the equivalence between $\overline{W}^{\bar{T}}_1(.)$ and $W_1(.)$ in the region of
$\left\{s_0|\overline{W}^{\bar{T}}_1(s_0) \leq \tau  \right\} $ for any $\tau \in [0,\bar{T})$,
for any $T<\bar{T}$, the minimal reachable set $\mathcal{R}_{\mathrm{min}}(K,T)$ can be characterized by:
\begin{align}
	\label{minreach}
	\begin{split}
	\mathcal{R}_{\mathrm{min}}(K,T)=&\left\{s_0|W_1(s_0)\leq T  \right\}\\
	=&\left\{s_0|\overline{W}^{\bar{T}}_1(s_0)\leq T  \right\}
	\end{split}
\end{align}
Similarly, for any $T<\bar{T}$, the maximal reachable set $\mathcal{R}_{\mathrm{max}}(K,T)$ can be characterized by:
\begin{align}
	\label{maxreach}
	\begin{split}
	\mathcal{R}_{\mathrm{max}}(K,T)=&\left\{s_0|W_2(s_0)\leq T  \right\}\\
	=&\left\{s_0|\overline{W}^{\bar{T}}_2(s_0)\leq T  \right\}
	\end{split}
\end{align}

\subsection{Recursive Formulas of the Modified Value Functions}
According to Equ. (\ref{minreach}) and Equ. (\ref{maxreach}),
The crucial point in characterizing the minimal and maximal reachable sets is to compute the modified value functions.

For any $\tau\in [0,\bar{T}]$, the following equations hold:
\begin{align}
	\label{recur1}
	\begin{split}
		&\overline{W}^{\bar{T}}_1(s_0)=\max_{u(.)\in \mathscr{U}} \int_0^{\bar{T}} \bar{c}\left( \bar{\phi}_0^{\bar{T}}(t,s_0,u)  \right) dt \\
		=&\max_{u(.)\in \mathscr{U}} \left[\int_0^{\tau} \bar{c}\left( \bar{\phi}_0^{\bar{T}}(t,s_0,u)  \right) dt + \right. \\
		&\left. \max_{u'(.)\in \mathscr{U}} \int_\tau^{\bar{T}} \bar{c}\left( \bar{\phi}_0^{\bar{T}}(t,s_0',u')  \right)  \right] \\
		=&\max_{u(.)\in \mathscr{U}} \left[\int_0^{\tau} \bar{c}\left( \bar{\phi}_0^{\bar{T}}(t,s_0,u)  \right) dt +
		\overline{W}^{\bar{T}-\tau}_1(s_0') \right]
	\end{split}
\end{align}
and
\begin{align}
	\label{recur2}
	\begin{split}
		&\overline{W}^{\bar{T}}_2(s_0)=\min_{u(.)\in \mathscr{U}} \int_0^{\bar{T}} \bar{c}\left( \bar{\phi}_0^{\bar{T}}(t,s_0,u)  \right) dt \\
		=&\min_{u(.)\in \mathscr{U}} \left[\int_0^{\tau} \bar{c}\left( \bar{\phi}_0^{\bar{T}}(t,s_0,u)  \right) dt + \right. \\
		&\left. \min_{u'(.)\in \mathscr{U}} \int_\tau^{\bar{T}} \bar{c}\left( \bar{\phi}_0^{\bar{T}}(t,s_0',u')  \right)  \right] \\
		=&\min_{u(.)\in \mathscr{U}} \left[\int_0^{\tau} \bar{c}\left( \bar{\phi}_0^{\bar{T}}(t,s_0,u)  \right) dt +
		\overline{W}^{\bar{T}-\tau}_2(s_0') \right]
	\end{split}
\end{align}
where $s_0'=\bar{\phi}_0^{\bar{T}}(\tau,s_0,u)$.

If $\Delta t$ is small enough, $u(.)$ can be regarded as a constant in interval $[k\Delta t,(k+1)\Delta t]$ for $k\in \mathbb{N}$.
And system (\ref{sys2}) can be converted into the discretized form:
\begin{align}
	\begin{split}
		s(t+\Delta t)&=F(s(t),u(t))\\
		&=\begin{cases}
			s(t)+f(s(t),u(t))\Delta t, &s(t)\notin K\\
			s(t), &s(t)\in K
		\end{cases}
	\end{split}
\end{align}
Replace $\bar{T}$ and $\tau$ with $(k+1)\Delta t$ and $\Delta t$ in Equ. (\ref{recur1}) and Equ. (\ref{recur2}), then we have:
\begin{align}
	\label{recur3}
	\begin{split}
		&\overline{W}^{0}_1(s_0)=0\\
		&\overline{W}^{(k+1)\Delta t}_1(s_0)\approx \max_{u\in\mathcal{U}} \left[\bar{c}(s_0)\Delta t+ \overline{W}^{k\Delta t}_1 \left( F(s_0,u) \right)  \right]
	\end{split}
\end{align}
and
\begin{align}
	\label{recur4}
	\begin{split}
		&\overline{W}^{0}_2(s_0)=0\\
		&\overline{W}^{(k+1)\Delta t}_2(s_0)\approx \min_{u\in\mathcal{U}} \left[\bar{c}(s_0)\Delta t+ \overline{W}^{k\Delta t}_2 \left( F(s_0,u) \right)  \right]
	\end{split}
\end{align}
Equ. (\ref{recur3}) and Equ. (\ref{recur4}) are the recursive formulas of the modified value functions.

\subsection{Approximation of the Modified Value Functions}
Normally, for any $k\in \mathbb{Z}^+$, it is difficult to compute the analytic forms of
$\overline{W}^{k\Delta t}_1(.)$ and $\overline{W}^{k\Delta t}_2(.)$.
This subsection introduce a method based on interpolation to approximate these functions.
Take a three-dimensional system as an example, where $s=\left[x,y,z\right]^\mathrm{T}$.
The proposed method can be divided into the following steps:
\begin{itemize}
	\item[(1)] Select a proper rectangular computational domain
	$\Omega=\left[x_d,x_u \right]\times \left[y_d,y_u \right] \times \left[z_d,z_u \right]$,
	and $K\subset \Omega$. And discretize $\Omega$ into a $N_x\times N_y \times N_z$ Cartesian grid structure,
	then the grid points can be expressed as a set:

	\begin{align}
		\mathcal{X}=\left\{ \left[ x_d+ i_x \delta x,y_d+ i_y \delta y,z_d+  i_z \delta z
		 \right]^\mathrm{T}  \bigg|\right. i_x=0, \nonumber \\
		...,N_x-1, i_y=0,...,N_y-1,i_z=0,...,N_z-1 \bigg\}
	\end{align}
	where $\displaystyle{\delta x=\frac{ x_u-x_d }{ N_x-1 },\delta y=\frac{ y_u-y_d }{ N_y-1 },\delta z=\frac{ z_u-z_d }{ N_z-1 }}$.
	Let $\mathcal{X}[i_x][i_y][i_z]$ be an element in $\mathcal{X}$, and
	\begin{align}
		\mathcal{X}[i_x][i_y][i_z]=\left[ x_d+i_x\delta x,y_d+i_y \delta y,z_d+i_z \delta z \right]^\mathrm{T}
	\end{align}
	\item[(2)] Let $ \mathcal{W}$ be a $ N_x\times N_y \times N_z $ array and $ \mathcal{V}[i_x][i_y][i_z] $ denotes the
	value function of state $\mathcal{X}[i_x][i_y][i_z]$. Initialize $ \mathcal{W}$ as $\mathbf{0}$.
	\item[(3)] Let $ \mathcal{W}' $ be a $ N_x\times N_y \times N_z $ array and for $i_x=0,...,N_x-1$, $i_y=0,...,N_y-1$,
	$i_z=0,...,N_z-1$, compute $\mathcal{W}'[i_x][i_y][i_z]$ according to Equ. (\ref{algequ1}).
	\begin{align}
		\label{algequ1}
		\begin{split}
		&\mathcal{W}'[i_x][i_y][i_z]= \\
		&\begin{cases}
		\max\limits_{u\in \mathcal{U}}  \left[
		\hat{c}\left( \mathcal{X}[i_x][i_y][i_z] \right) \Delta t +
		\widehat{W}\left( F(\mathcal{X}[i_x][i_y][i_z],u)  \right) \right], \\
		\quad \quad   \text{for computation of } \overline{W}^{\bar{T}}_1(.)\\

		\min\limits_{u\in \mathcal{U}}  \left[
		\hat{c}\left( \mathcal{X}[i_x][i_y][i_z] \right) \Delta t +
		\widehat{W}\left( F(\mathcal{X}[i_x][i_y][i_z],u)  \right) \right], \\
		\quad \quad  \text{for computation of } \overline{W}^{\bar{T}}_2(.)
		\end{cases}
	\end{split}
	\end{align}
	where $\widehat{W}(.)$ is a trilinear interpolation \cite{web4} function using $\mathcal{W}$. Then copy $\mathcal{W}'$ to $\mathcal{W}$.
	\item[(4)] After $k$ times of repetitions of Step (3), the interpolation function $\widehat{W}(.)$ is an
	approximation of the value function
	$\overline{W}^{k\Delta t}_1(.)$ or $\overline{W}^{k\Delta t}_2(.)$.
	\item[(5)] Finally, if $k\Delta$ is greater than the given time horizon $T$,
	the minimal and maximal reachable sets can be characterized by:
	\begin{align}
		\begin{split}
		\mathcal{R}_{\mathrm{min}}(K,T)=\left\{s_0 \in \Omega| \overline{W}^{k\Delta t}_1(s_0)\leq T \right\}\\
		\mathcal{R}_{\mathrm{max}}(K,T)=\left\{s_0 \in \Omega| \overline{W}^{k\Delta t}_2(s_0)\leq T \right\}
		\end{split}
	\end{align}
\end{itemize}

The pseudocode of these steps is shown in Algorithm 1.

\begin{algorithm*}[htbp]
    \caption{Computation of modified value functions}
	\begin{algorithmic}[1]
	  \State \textbf{Input:} Recursion number $m$, number of grid points $N_x,N_y,N_z$, computational domain
	  $\Omega=\left[x_d,x_u \right]\times \left[y_d,y_u \right] \times \left[z_d,z_u \right]$;

	  \State Construct the modified system in discretized form:
	  \begin{align*}
		\begin{split}
			s(t+\Delta t)=F(s(t),u(t))
			=\begin{cases}
				s(t)+f(s(t),u(t))\Delta t, &s(t)\notin K\\
				s(t), &s(t)\in K
			\end{cases}
		\end{split}
	\end{align*}
	\State Construct the modified running cost:
	\begin{align*}
		\bar{c}(s) =\begin{cases}
			1, &s\notin K\\
			0, &s\in K
		\end{cases}
	\end{align*}

	  \State Discretize $\Omega$ into a $N_x\times N_y \times N_z$ Cartesian grid structure, and
	  store the coordinates of the grid points to a three dimensional array $\mathcal{X}$;
	  \State Let $\mathcal{W}$ and $\mathcal{W}'$ be two $N_x\times N_y \times N_z$ arrays;
	  \For{$i_x\leftarrow 0,...,N_x-1 $}
		\For{$i_y\leftarrow 0,...,N_y-1 $}
			\For{$i_z\leftarrow 0,...,N_z-1 $}
				  \State $\mathcal{W}[i_x][i_y][i_z]\leftarrow 0$;
				  \State $\mathcal{X}[i_x][i_y][i_z]\leftarrow \left[ x_d+i_x\delta x,y_d+i_y \delta y,z_d+i_z \delta z \right]^\mathrm{T}$;
			\EndFor
        \EndFor
	  \EndFor
	  \\

	  \For{$k \leftarrow 1,...,m$}
	  	\State Construct the trilinear interpolation function $\widehat{W}(.)$ using $\mathcal{W}$;
        \For{$i_x\leftarrow 0,...,N_x-1 $}
		  \For{$i_y\leftarrow 0,...,N_y-1 $}
		    \For{$i_z\leftarrow 0,...,N_z-1 $}

		\State
		\begin{align*}
			\mathcal{W}'[i_x][i_y][i_z]\leftarrow \begin{cases}
			\max\limits_{u\in \mathcal{U}}  \left[
			\hat{c}\left( \mathcal{X}[i_x][i_y][i_z] \right) \Delta t +
			\widehat{W}\left( F(\mathcal{X}[i_x][i_y][i_z],u)  \right) \right], & \text{for computation of } \overline{W}^{\bar{T}}_1(.)\\

			\min\limits_{u\in \mathcal{U}}  \left[
			\hat{c}\left( \mathcal{X}[i_x][i_y][i_z] \right) \Delta t +
			\widehat{W}\left( F(\mathcal{X}[i_x][i_y][i_z],u)  \right) \right], & \text{for computation of } \overline{W}^{\bar{T}}_2(.)
			\end{cases}
		\end{align*}
			\EndFor
          \EndFor
        \EndFor
        \State Copy $\mathcal{W}'$ to $\mathcal{W}$;
      \EndFor
      \State Return $\widehat{W}(.)$;
    \end{algorithmic}
\end{algorithm*}

In the following two sections, we illustrate the accuracy and generality of the proposed method with some numerical examples.
In order to verify the accuracy, we provide some examples on affine nonlinear systems with regular target sets,
and compare the computational results with those of the level set method.
To demonstrate generality, we provide some examples that are difficult to be solved by level set method,
where the dynamical systems are non-affine nonlinear and the target sets are irregular.

\section{Longitudinal flight control examples}
Consider the dynamics of the longitudinal motion of an aircraft.
Longitudinal motion is the subsystem in the aircraft's plane of
symmetry ($x-z$ plane) and is characterized by the following state:
velocity $v$, attack angle $\alpha$, pitch rate $q$, and pitch angle $\theta$.
These variables are described in the ground coordinate system, see Fig. \ref{figlong}.
\begin{figure}[H]
	\centering
	\includegraphics[width=6cm]{fig/corr.pdf}
	\caption{Longitudinal motion of an aircraft.}
	\label{figlong}
\end{figure}

The dynamic equation of the longitudinal motion can be written as Equ. (\ref{longair}) \cite{a009}:
\begin{align}
	\label{longair}
	\begin{split}
		&\dot{s}=\left[  \dot{v}, \dot{\alpha} , \dot{q} , \dot{\theta}   \right]^\mathrm{T}=\\
		&\left[
		\begin{array}{c}
			\displaystyle{\frac{1}{m} \left(P\cos \alpha-\frac{1}{2}\rho v^2 S C_D(\alpha,\delta_e)-mg\sin(\theta-\alpha) \right)}\\
			\displaystyle{q+\frac{1}{mv} \left(-P\sin \alpha-\frac{1}{2}\rho v^2 S C_L(\alpha,\delta_e)+mg\cos(\theta-\alpha) \right)}\\
			\displaystyle{\frac{1}{2I_{yy}}\rho v^2 S \bar{b} C_m(\alpha,q,\delta_e) }\\
			\displaystyle{q }
		\end{array}\right]
	\end{split}
\end{align}

where $s=\left[v,\alpha,q,\theta  \right]^\mathrm{T}$ is system state, $\delta_e \in [-0.3,0.3]$ is the elevator deflection
and is regarded as the control input,
and $P$ is the thrust, $m$ is the mass of the aircraft, $\rho$ is the air density, $I_{y}$ is the mass moment of inertia about $y$ axis,
$C_L,C_D,C_m$ are the aerodynamic derivatives.
In the following examples of this section, it is assumed that
the thrust $P$ is varied to keep the velocity $v$ be a constant, that is:
\begin{align}
	\begin{split}
	&\dot{v}=0 \Longrightarrow \\
	&P=\frac{1}{\cos \alpha} \left[\frac{1}{2} \rho v^2 S  C_D(\alpha,\delta_e) +mg \sin (\theta-\alpha)\right]
	\end{split}
\end{align}

Under this assumption, the state space degenerates to a three-dimensional space,
i.e., $s=\left[\alpha,q,\theta  \right]^\mathrm{T}$.
The aerodynamic model is constructed as the following equations:
\begin{align}
    \label{longaero}
    \left\{\begin{array}{l}
        \displaystyle{C_L=C_{L0}+C_{L\alpha}\alpha+C_{L\delta_e}\delta_e}\\
        \displaystyle{C_D=C_{D0}+C_{D\alpha}\alpha+C_{D\alpha 2}\alpha^2+C_{D\delta_e}\delta_e}\\
        \displaystyle{C_m=C_{m0}+C_{m\alpha}\alpha+C_{mq}\frac{q\bar{b}}{2v}+C_{m\delta_e}\delta_e}
    \end{array}\right.
\end{align}
The parameters in Equ. (\ref{longair}) and Equ. (\ref{longaero}) are summarized
in Table \ref{table1}.
\begin{table}[H]
    \centering
    \caption{Physical parameters of the reference aircraft}
    \label{table1}
    \begin{tabular}{ll|ll}
    \hline
    \textbf{Parameter}    & \textbf{Value}    &  \textbf{Parameter}                & \textbf{Value} \\ \hline
    $m$       &  235717           &  $I_{y}$   &  22428285  \\
    $S$   &  524             &  $\bar{b}$                    &  6.32 \\
    $\rho$     &  1.293        &  $C_{L0}$                                     &  0.1 \\
    $C_{L\alpha}$          &  2.4         &       $v $       &  200   \\
    $C_{L \delta e}$      &  0.2        &  $C_{D0}$                                     &  0.00108 \\
    $C_{D\alpha}$         &  0.01         &  $C_{D\alpha 2}$      &  0.6  \\
	$C_{m\delta e}$           &  -1.2          &  $C_{D \delta e}$    &  0.05 \\
    $C_{m0}$                                     &  0.04         &  $C_{m \alpha}$        &  -0.2 \\
    $C_{mq}$        &  -1         &            &   \\ \hline
    \end{tabular}
\end{table}

\subsection{Maximal reachable set}
In this example, the task is to drive the state trajectory into
the target set $K=[-0.05,0.05]\times [-0.1,0.1] \times [-0.1,0.1]$ in the given time horizon $T=1$ as much as possible.
Thus the maximal reachable needs to be computed.
The solver settings of our method for the maximal reachable set computation are:
\begin{itemize}
	\item \textbf{Computational domain}: $\Omega=[-0.4,0.4] \times [-0.75,0.75] \times [-0.75,0.75]$.
	\item \textbf{Number of grid points}:  $101\times 101\times 101$.
	\item \textbf{Time step size}: $\Delta t=0.01$.
	\item \textbf{Number of recursions}: $101$ ($101 \times \Delta t > T$).
\end{itemize}
The maximal reachable set can be characterized by a $1$ level set of the value function $\overline{W}_2^{101\Delta t}(.)$:
\begin{align}
	\mathcal{R}_{\text{min}}(K,T)=\left\{s_0\left|\overline{W}_2^{101\Delta t}(s_0)\leq 1 \right. \right\}
\end{align}

This example can also be solved by the level set method.
In level set method, the maximal reachable set can be characterized by the complement of the invariant set of the target set's complement,
therefore, the signed distance function as the terminal condition of the HJ PDE can be set as:
\begin{align}
	\begin{split}
	V_2(s,T)=l(s)=\max ( -0.05-\alpha,\alpha-0.05,\\ -0.1-q, q-0.1,-0.1-\theta,\theta-0.1   )
	\end{split}
\end{align}
The Hamiltonian can be written as:
\begin{align}
	\begin{split}
	\mathcal{H}(s,p_2)=&p_2^1 \bigg(q+\frac{1}{mv} \big(-\tan \alpha (\bar{D}+mg\sin (\theta-\alpha))\\
	 &+\bar{L} +mg \cos(\theta-\alpha)   \big)   \bigg)+p_2^2 \frac{\bar{M}}{I_{y}}+p_2^3 q\\
	 &-\delta_e |p_2^1 k_{e1}+p_2^2 k_{e2} |
	\end{split}
\end{align}
and the dissipation function can be simplified as:
\begin{align}
	\begin{split}
		\mathcal{P}_1=&\bigg| q+\frac{1}{mv} \big(-\tan \alpha (\bar{D}+mg\sin (\theta-\alpha))+\bar{L} \\
		&+ mg \cos(\theta-\alpha)   \big) \bigg|+ |k_{e1}|\delta e\\
		\mathcal{P}_2=&\left|\frac{\bar{M}}{I_{y}} \right| + |k_{e2}|\delta e\\
		\mathcal{P}_3=&|q|
	\end{split}
\end{align}
where
\begin{align*}
	\begin{split}
		&\bar{D}=\frac{1}{2}\rho v^2 S \left(C_{D0}+C_{D\alpha}\alpha+C_{D\alpha 2}\alpha^2 \right)\\
		&\bar{L}=\frac{1}{2}\rho v^2 S \left(C_{L0}+C_{L\alpha}\alpha \right)\\
		&\bar{M}=\frac{1}{2}\rho v^2 S \bar{b} \left( C_{m0}+C_{m\alpha}\alpha+C_{mq}\frac{q\bar{b}}{2v}   \right)\\
		&k_{e1}= \frac{\rho v^2 S \left(C_{D\delta e} \tan \alpha - C_{m\alpha }  \right)}{2mv} \\
		&k_{e2}= -\frac{1}{2I_{yy}} \rho v^2 S \bar{b} C_{m\delta e}
	\end{split}
\end{align*}
The setups of the computational domain and the grid points' number in level set method are the same as those in our method.
The maximal reachable set can be represented by the zero level set of the value function $V_2(.,T)$,i.e.:
\begin{align}
	\mathcal{R}_{\text{min}}(K,T)=\left\{s_0\left|V_2(s_0,T)\leq 0 \right. \right\}
\end{align}
The comparison between the result of the proposed method and level set method is shown in Fig. \ref{figexp1}.
It can be seen that, the envelope of the maximal reachable set computed by the proposed method almost
coincides with that computed by the level set method.
This indicates that the proposed method has a high accuracy in computing the maximal reachable set.
\begin{figure}[H]
	\centering
	\includegraphics[width=8.85cm]{fig/figexp1.png}
	\caption{Computation results of the maximum reachable set (the result of our method is expressed in blue
	and the result of the level set method is expressed in red).}
	\label{figexp1}
\end{figure}

\subsection{Viable set}
In this example, the task is to maintain the state inside the target set $K=[-0.35,0.35]\times[-0.35,0.35] \times[-0.35,0.35]$
for more than the given time horizon $T=1.5$. This is a problem about the computation of viable set.
In our method, the viable set is characterized by the complement of the minimal reachable set of the target set's complement.
The solver settings of our method for the minimal reachable set computation are:
\begin{itemize}
	\item \textbf{Computational domain}: $\Omega=[-0.4,0.4] \times [-0.4,0.4] \times [-0.4,0.4]$.
	\item \textbf{Number of grid points}:  $101\times 101\times 101$.
	\item \textbf{Time step size}: $\Delta t=0.01$.
	\item \textbf{Number of recursions}: $151$ ($151 \times \Delta t > T$).
\end{itemize}
The viable set can be characterized by a $1.5$ level set of the value function $\overline{W}_1^{151\Delta t}(.)$:
\begin{align}
	\mathcal{V}(K,T)=\left\{s_0\left|\overline{W}_1^{151\Delta t}(s_0)\geq 1.5 \right. \right\}
\end{align}

In level set method,
the signed distance function as the terminal condition of the HJ PDE can be set as:
\begin{align}
	\begin{split}
	V_1(s,T)=l(s)=\min ( \alpha+0.35,0.35-\alpha,\\ q+0.35,0.35-q,\theta+0.35,0.35-\theta   )
	\end{split}
\end{align}
The Hamiltonian can be written as:
\begin{align}
	\begin{split}
	\mathcal{H}(s,p_1)=&p_1^1 \bigg(q+\frac{1}{mv} \big(-\tan \alpha (\bar{D}+mg\sin (\theta-\alpha))\\
	 &+\bar{L} +mg \cos(\theta-\alpha)   \big)   \bigg)+p_1^2 \frac{\bar{M}}{I_{y}}+p_1^3 q\\
	 &+\delta_e |p_1^1 k_{e1}+p_1^2 k_{e2} |
	\end{split}
\end{align}
The dissipation function is set the same as in the previous example.
The viable set can be represented by the zero level set of the value function $V_1(.,T)$,i.e.:
\begin{align}
	\mathcal{V}(K,T)=\left\{s_0\left|V_1(s_0,T)\geq 0 \right. \right\}
\end{align}

The comparison between the result of the proposed method and level set method is shown in Fig. \ref{figexp2}.
Fig. indicates that the viable set characterized by the proposed method is the same as that characterized by level set method.
This illustrates that the proposed method also has a high accuracy in computing the viable set.
\begin{figure}[H]
	\centering
	\includegraphics[width=8.85cm]{fig/figexp2.png}
	\caption{Computation results of the viable set (the result of our method is expressed in blue
	and the result of the level set method is expressed in red).}
	\label{figexp2}
\end{figure}

\section{Aircraft ground motion examples}
Consider a simplified aircraft ground motion dynamics system:
\begin{align}
	\dot{s}=\left[
		\begin{array}{c}
			\dot{v}_x\\\dot{v}_y\\\dot{r}
		\end{array}
	\right]=
	\left[\begin{array}{c}
		\displaystyle{rv_y+\frac{F_x}{m}}\\
		\displaystyle{-rv_x+\frac{F_y}{m}}\\
		\displaystyle{\frac{M_z}{I_{z}}}
	\end{array}\right]
\end{align}
where $s=[v_x,v_y,r]^\mathrm{T}$ is system state.
$v_x$ and $v_y$ are the longitudinal and lateral velocities, respectively, $r$ is the yaw angular velocity,
$F_x$ and $F_y$ are the longitudinal and lateral resultant forces, $M_z$ is the resultant moment.
The expressions for $F_x,F_y,M_z$ are as follows:
\begin{align}
	\begin{cases}
		F_x=-D \cos \beta - Y \sin \beta -Q_n-Q_{ ml }-Q_{mr}+P\\
		F_y=Y \cos \beta -D\sin\beta-F_n-F_{ml}-F_{mr}\\
		\displaystyle{M_z=n+\left(F_{ml}+F_{mr} \right)a_m-F_n a_n +\left(  Q_{ ml }+Q_{mr} \right)\frac{b_w}{2}}
	\end{cases}
\end{align}
where $D=\frac{1}{2} \rho v^2 S C_{Y\beta} \beta$ is the drag, $Y=\frac{1}{2} \rho v^2 S C_{D0}$ is the aerodynamic lateral force,
$F_n,F_{ml},F_{mr}$ are the ground lateral forces, $Q_n,Q_{ ml },Q_{mr}$ are the ground longitudinal forces.
See Fig. \ref{figground} for the specific meanings of these variables.
\begin{figure}[H]
	\centering
	\includegraphics[width=8cm]{fig/TU.pdf}
	\caption{Aircraft ground motion}
	\label{figground}
\end{figure}
Assume that the ground friction coefficient $\mu=0.04$, then
these variables can be derived by solving the force and moment balance equations \cite{a0016}.
It is worth pointing out that the expression of the lateral force of the nose-wheel is:
\begin{align}
	\label{eqfn}
	\begin{split}
	F_n=-\mu_d \sin \left\{ \mu_c \arctan \left(\mu_b \tan \beta_n\right) \right. - \\
	 \left\{ \mu_e \left[\mu_b \tan \beta_n-\arctan \left(\mu_b \tan \beta_n\right)  \right] \right\}P_n
	\end{split}
\end{align}
and
\begin{align}
	\label{eqbn}
	\beta_n=\arctan \frac{ (v+ra_n)\cos \theta_W-u\sin \theta_W  }{ u\cos \theta_W+(v+ra_n)\sin \theta_W  }
\end{align}
where $\theta_W\in [-0.15,0.15]$ is the nose-wheel deflection and is considered as the control input.
See literatures \cite{a0017,a0018} for the meaning of each symbol in Equ. (\ref{eqfn}) and Equ. (\ref{eqbn}).
In summary, the aircraft ground motion is a highly nonlinear system, that it cannot be translated into an affine nonlinear form.
The parameters of the reference aircraft are summarized in Table \ref{table2}.
\begin{table}[H]
    \centering
    \caption{Physical parameters of the reference aircraft}
    \label{table2}
    \begin{tabular}{ll|ll}
    \hline
    \textbf{Parameter}    & \textbf{Value}    &  \textbf{Parameter}                & \textbf{Value} \\ \hline
    $m$       &  104915.9           &  $I_{z}$   &  10504308.1  \\
    $S$   &  249.9             &  $\bar{b}$                    &  12.1 \\
    $\rho$     &  1.293        &  ${b}_w$                  &  23.8 \\
    $a_n$          &  17.9         &       $a_m$       &  2.3   \\
    $b_W$      &  6.9        &  $\mu_d$       &  0.1014 \\
    $\mu_b$         &  -10.11         &  $\mu_c$      &  1.438  \\
	$\mu_e$           &  -0.8507         &  $C_{D0}$    &  0.061 \\
    $C_{Y\beta}$      &  -1.4            &  $C_{n\beta}$        &  0.2 \\
	$C_{nr}$        &  -1.5         &   $H$    &   \\
	$C_{L0}$        &  -0.053         &        &    \\\hline
    \end{tabular}
\end{table}

\subsection{Maximal reachable set}
In this example, we let the target set $K$ be a set with irregular shape,
which is regarded as a combination of some cubic cells,
and whether a cell is
contained in the target set is determined by its center, see Fig. \ref{figirr}.
\begin{figure}[H]
	\centering
	\includegraphics[width=8.85cm]{fig/figirrtarget.png}
	\caption{Irregular target set}
	\label{figirr}
\end{figure}
The time horizon $T=2$.
The solver settings of our method in this example are:
\begin{itemize}
	\item \textbf{Computational domain}: $\Omega=[20,50] \times [-20,20] \times [-0.5,0.5]$.
	\item \textbf{Number of grid points}:  $101\times 101\times 101$.
	\item \textbf{Time step size}: $\Delta t=0.01$.
	\item \textbf{Number of recursions}: $201$.
\end{itemize}
The maximal reachable set is shown in See Fig. \ref{figexp3}.:
\begin{figure}[H]
	\centering
	\includegraphics[width=8.85cm]{fig/figexp3.png}
	\caption{Computation result of the maximal reachable set.}
	\label{figexp3}
\end{figure}

\subsection{Viable set}
In this example, we aim to keep the aircraft state in the target set $K=[22,48]\times[-18,18] \times[-0.4,0.4]$
for more than the given time horizon $T=2$. The solver settings are the same as those in the previous example.
The computation result is shown in Fig. \ref{figexp4}.
\begin{figure}[H]
	\centering
	\includegraphics[width=8.85cm]{fig/figexp4.png}
	\caption{Computation result of the viable set.}
	\label{figexp4}
\end{figure}

\section{Conclusions}
The works presented in this paper propose a unified framework for addressing
reachability, viability and invariance problems. In this framework,
the reachable, viable and invariant sets are characterized by some non-zero level sets
of different value functions which are approximated by recursion and interpolation.
Such a mechanism can effectively avoid the construction of
dissipation function, which is an optimization problem in level set method
and is difficult to be solved in non-affine nonlinear systems.
Therefore, compared to the level set method, the generality of the proposed method is relatively better.
It also avoids the construction of the signed distance function of the target set
which act as the terminal condition of the HJ PDE, such that it can
handle the irregular target sets.

The primary weakness of the proposed method,
and many others, is the exponential growth of memory and computational cost as the system dimension increases \cite{a008,i7}.
Some approaches have been proposed to mitigate these costs, such as
projecting the reachable/viable/invariant
set of a high-dimensional system into a collection of lower
dimensional subspaces \cite{a002,i8,i9} and
exploiting the structure of the
system using the principle of timescale separation and
solving these problems in a sequential manner \cite{i10,i11,i12}.

To improve the interpolation algorithm to reduce the required quantity of grids might be another effective approach
to overcome the curse of dimensionality. These will be considered in our future work.

\section*{Acknowledgements}
\par The authors gratefully acknowledge support from National Defense Outstanding Youth Science Foundation (Grant No. 2018-JCJQ-ZQ-053),
and Central University Basic Scientific Research Operating Expenses Special Fund Project Support (Grant No. NF2018001).
Also, the authors would like to thank
the anonymous reviewers, associate editor, and editor for
their valuable and constructive comments and suggestions.

\bibliographystyle{unsrt}