3.1 Definition of RS

Let $\mathcal{U}$ be a nonempty compact set, and $\mathcal{U}$ represents the collection consisting of measurable functions from $[0, + \infty )$ to $\mathcal{U}$ . We consider a control system with state $s \in {\mathbb{R}^n}$ and a control input $u \in \mathcal{U}$ . The system dynamics are determined by:

(3) \begin{align}\dot s = f(s,u).\end{align}

We assume that $f(.,.)$ is bounded and Lipschitz continuous. Then fixing the state at s at time t and given a control input $u(.) \in \mathcal{U}$ , the evolution of system (3) is a continuous trajectory $\xi _{s,t}^u(.):[t,\infty ) \to {\mathbb{R}^n}$ and $\xi _{s,t}^u(t) = s$ . With a given target set $K \subset {\mathbb{R}^n}$ and a time horizon $\bar T$ , the RS is defined as:

Table 1.

Parameters of the reference aircraft

Definition 1. Reachable set:

(4) \begin{align}\mathcal{R}(\bar T,K) = \left\{ {s|\exists u(.) \in \mathcal{U}\;\exists \tau \in [0,\bar T]\;\xi _{s,0}^u(\tau ) \in K} \right\}.\end{align}

3.2 Definition of CRS

Set a running cost function $c(.,.):{\mathbb{R}^n} \times \mathcal{U} \to \mathbb{R}$ for the system. Given a control input $u(.)$ , the performance index of the evolution start from state s at time ${t_0}$ in time interval $[{t_0},{t_1}]$ is:

(5) \begin{align}\mathcal{J}_{{t_0}}^{{t_1}}(s,u) = \int_{{t_0}}^{{t_1}} c\left( {\xi _{s,{t_0}}^u(t),u(t)} \right)dt.\end{align}

Then, with a target set $K \in {\mathbb{R}^n}$ and an admissible cost $\bar J$ , the CRS is defined as:

Definition 2. Cost-limited reachable set:

(6) \begin{align}\begin{array}{l}{{\mathcal{R}^c}(\bar J,K) = \left\{ {{s_0}|\exists u(.) \in \mathcal{U}\;\exists \tau \in [0,\infty )\;} \right.}\\[5pt]{\left. {\xi _{s,0}^u(\tau ) \in K \wedge \mathcal{J}_0^\tau ({s_0},u) \le \bar J} \right\}.}\end{array}\end{align}

Remark 1. If $c(s,u) \equiv 1$ , then $\mathcal{J}_0^t(s,u) = t$ and the CRS is degenerated into the RS. Thus, the RS can be considered as a special case of the CRS.

3.3 Method to compute CRS

Literature [Reference Liao, Wei and Lai23] suggests that the CRS can be obtained by constructing a time-varying scalar function $V(.,.):{\mathbb{R}^n} \times \mathbb{R} \to \mathbb{R}$ via the following recursion:

(7) \begin{align}\begin{array}{l}{V(s,0) = 0}\\[5pt]{}{V(s,(n + 1)\Delta t) = }\\[5pt] \left\{ {\begin{array}{c}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{V(s,n\Delta t),}\;{s \in K}\\[5pt]{\mathop {\min }\limits_{u \in \mathcal{U}} \left[ {c(s,u)\Delta t + V(F(s,u),n\Delta t)} \right],}\;{s \notin K}\end{array}} \right.\,\,,\end{array}\end{align}

where $F(s,u)$ is the time discrete form of system (3), i.e. $s(t + \Delta t) = F(s(t),u(t))$ . Denote $\mathop {\min }\nolimits_{s,u} c(s,u) = \gamma $ , then for any $\bar J \lt \gamma \,h\Delta t$ , the CRS can be characterised as:

(8) \begin{align}{\mathcal{R}^c}(\bar J,K) = \{ s \in {\mathbb{R}^n}|V(s,h\Delta t) \le \bar J\} .\end{align}

Table 2.

Parameters for the trim set computation

In order to approximate function $V(.,.)$ , it is necessary to specify a rectangular region in the state space as the computational domain and discretise it into a Cartesian grid structure. The function $V(.,.)$ is represented as a grid interpolation.

4.1 Longitudinal trim set

As mentioned earlier, the flight safety envelope is regarded as a CRS of the trim set. Therefore, the trim set needs to be calculated first. The trim set is a set of aircraft states in which the derivatives of each state of the aircraft can be kept at 0 by adjusting the control inputs, and therefore the trim set is necessarily in the $\alpha - \theta $ -plane.

Table 3.

Solver settings for the CRS computation

When calculating the trim set, the $\alpha - \theta $ -plane is first discretised into a Cartesian grid, and then each grid point is verified one by one to ensure whether it can meet the trim requirements. The settings in the computation process are shown in Table 2.

Figure 2(a) shows the longitudinal trim envelope at Mach 0.5 at sea level. The whole $\alpha - \theta $ plane is divided into $101 \times 101$ cells, the centre of each cell is a grid point. In the current research, whether a cell is contained in the longitudinal trim envelope is determined by its centre, see Fig. 2(b).

Figure 2.

Trim envelope of the reference aircraft.

4.2 Running cost function

In order to take the influence of manoeuver overload into consideration, a term related to manoeuver overload is needed to be added to the running cost function. Manoeuver overload is the ratio between the summation of aerodynamic force and thrust and the magnitude of gravity [Reference Jiang, Zhou, Gao, Yang and Jing24], i.e.:

(9) \begin{align}G = \frac{{T + L + D}}{{mg}},\end{align}

where L and D present the lift and drag, respectively. The running cost function is set as:

(10) \begin{align}c(s,u) = 1 + {\lambda _G}\parallel \!G{\parallel _2} = 1 + \frac{{{\lambda _G}}}{{mg}}\sqrt {|{F_z}{|^2} + |{F_x}{|^2}} ,\end{align}

where, ${\lambda _G}$ is the weight of manoeuver overload. ${F_x}$ and ${F_z}$ denote the combined force of the thrust and the aerodynamic force in the horizontal and vertical directions, respectively, i.e.:

(11) \begin{align}\begin{array}{l}{{F_x} = T\cos \theta - \frac{1}{2}\rho {v^2}S{C_L}\sin (\theta - \alpha ) - \frac{1}{2}\rho {v^2}S{C_D}\cos (\theta - \alpha )}\\[5pt]{}{{F_z} = T\sin \theta + \frac{1}{2}\rho {v^2}S{C_L}\cos (\theta - \alpha ) - \frac{1}{2}\rho {v^2}S{C_D}\sin (\theta - \alpha ).}\end{array}\end{align}

All of these variables can be represented with state $(v,\alpha ,q,\theta )$ and control inputs $({\delta _f},{\delta _e})$ .

4.3 Computation of the CRS

By assuming constant velocity, the longitudinal flight safety envelope is degenerated into a three-dimensional envelope of $(\alpha ,q,\theta )$ . The flight safety envelope is defined as the CRS, which are computed by using the method introduced in Section 3. The parameters for the computation are summarised in Table 3.

Figures 3, 4, 5, 6 and 7 show the variation of the flight safety envelope with different ${\lambda _G}$ s for the given admissible cost $\bar J = 1$ . In each figure, the feasible state of the aircraft lies in the region that makes the value of the scalar function V less than or equal to $\bar J$ , i.e. the inner region of the envelope. Figure 3 appears to contain two surfaces because the flight safety envelope size is too large and exceeds the boundary of the computational domain. An interesting fact to notice here is that the flight safety envelope shrinks as ${\lambda _G}$ increases, which is also intuitive.

Figure 3.

Flight safety envelope with ${\lambda _G} = 0$ .

Figure 4.

Flight safety envelope with ${\lambda _G} = 0.25$ .

Figure 5.

Flight safety envelope with ${\lambda _G} = 0.5$ .

Figure 6.

Flight safety envelope with ${\lambda _G} = 0.75$ .

Figure 7.

Flight safety envelope with ${\lambda _G} = 1$ .

Figure 8.

Variation in the proportion of grid points inside the CRS.

To quantify the impact of ${\lambda _G}$ on the flight safety envelope, the number of grid points in the RS as a percentage of the total number of grid points in the computational domain is plotted against different ${\lambda _G}$ s in Fig. 8.

It can be seen that, with the increasing ${\lambda _G}$ , the size of the flight safety envelope decreases. In order to ensure that the aircraft’s state can reach the trim set before the performance index increases to the allowable cost, the greater the value of ${\lambda _G}$ , the more the pilot must be aware that the aircraft’s state is within the flight safety envelope.

4.4 Flight control law

When the aircraft’s state lies outside the trim set, the pilot needs to go through multiple time steps to adjust the aircraft’s state to the trim set. Denote the state at the ith time step as ${s_i}$ , according to Bellman’s principle of optimality, the optimal control input for the ith time step is

(12) \begin{align}{(\delta _f^*,\delta _e^*)_i} = \arg \mathop {\min }\limits_{{\delta _f},{\delta _e}} \left[ {c({s_i},u)\Delta t + V(F({s_i},u),(h - i)\Delta t)} \right].\end{align}

To verify the computational accuracy, we take many points in the state space as initial states and verify one by one that the aircraft states from these initial states, under the control law in Equation (12), can reach the trim set before the performance index increases to the allowable cost. If it can, then this initial state is represented by a grey pixel point, and if not, then this initial state is represented by a white pixel point. Figures 9, 10, 11, 12 and 13 are some slices of the state space.

Figure 9.

Simulation results of ${\lambda _G} = 0$ .

Figure 10.

Simulation results of ${\lambda _G} = 0.25$ .

Figure 11.

Simulation results of ${\lambda _G} = 0.5$ .

Figure 12.

Simulation results of ${\lambda _G} = 0.75$ .

Figure 13.

Simulation results of ${\lambda _G} = 1$ .

In the preceding figures, the grey areas are the sets of states that can be driven into the target sets under control law in Equation (12) before the performance index growing to the admissible cost, the blue curves are the outlines of CRSs computed by the proposed method. As can be seen, the outlines of the CRSs almost coincided with that of the grey areas. This indicates that our method has a high accuracy.

4.5 Computational complexity

When the dimension of the state space is n, the number of grids is $\prod\nolimits_{i = n}^n {N_i}$ , where ${N_i}$ denotes the number of grids in the ith dimension. In addition, the time complexity of multiple linear interpolation is proportional to ${2^n}$ . Therefore, the time complexity of the method used in this paper is $O\left( {{2^n}\prod\nolimits_{i = n}^n {N_i}} \right)$ . During the computation, the value of function $V(.,.)$ at each grid point needs to be kept in memory, thus, the space complexity is $\prod\nolimits_{i = n}^n {N_i}$ .

It is worth mentioning that in the level set method, the time consumption at each time step is proportional to $n\prod\nolimits_{i = n}^n {N_i}$ [Reference Nabi, Lombaerts, Zhang, van Kampen, Chu and de Visser11], and the length of the time step is proportional to the grid size in order to ensure numerical stability [Reference Mitchell25]. Therefore, the time complexity of the level set method is $O\left( {{n^2}\prod\nolimits_{i = n}^n {N_i}} \right)$ . The space complexity of the level set method is the same as the method used in this paper.