Estimation of safe flight envelope considering manoeuver overload based on cost-limited reachable set

Abstract

Depending on magnitude and duration, any manoeuvering overload can damage the structure of an aircraft and adversely affect the pilot’s concentration and reaction time. These are all threats to flight safety. The flight safety envelope estimation method based on the classical reachable set cannot take into account the effect of manoeuvering overload. To overcome this limitation, a generalized reachable set known as a cost-limited reachable set is introduced into the computation of flight safety envelopes in this paper. It differs from the classical reachable set in that the performance index of the system can be set as the time integral of a running cost, and it can discuss the ability to reach the trim set before the performance index grows to the admissible cost. When computing the flight safety envelope, the running cost is set as a weighted sum of time consumption and manoeuver overload factor, and the flight safety envelope is defined as a cost-limited reachable set of the trim set. The flight safety envelopes and optimal control laws under the different weight of manoeuver overload factors are analyzed.

Nomenclature

$C_L$

lift coefficient

$C_D$

drag coefficient

$C_m$

pitch moment coefficient

$m$

mass of the aircraft

$I_{yy}$

inertia moment about y axis

$T$

thrust

$D$

drag

$L$

lift

$v$

velocity

$\alpha$

angle of attack

$\theta$

pitch angle

$q$

pitch angle velocity

$\delta_f,\delta_e$

flap and elevator deflections

$u$

control input

$\mathcal{R}$

reachable set

$\mathcal{R}^c$

cost-limited reachable set

$K$

target set

$G$

maneuver overload

1.0 Introduction

In the coming years, aircraft will affect our lives even more profoundly, due to their applications in a variety of commercial areas, such as aerial photography, agriculture and transportation [1, 2]. Therefore, it is especially important to avoid flight accidents. Some studies [3] statistically analysed the aircraft accidents that occurred between 2006 and 2015, and the results indicated that loss of control (LOC) causes about 23% of all catastrophic accidents.

A flight safety envelope is a collection of states that ensure the safe manoeuver of an aircraft, and LOC occurs when the aircraft’s state leaves the flight safety envelope [4]. The flight safety envelope is of great significance and its estimation has attracted a lot of attention in the academic community. Many approaches have been proposed to estimate the flight safety envelope, such as wind tunnel experiments [5, 6], flight experiments and flight simulator-based methods [7, 8] or the combination of computational fluid dynamics (CFD) and extended Kalman filter [9, 10]. A more effective method to estimate the flight safety envelope is reachability analysis [11, 12]. In reachability analysis, a subset of the state space is specified as the target set, and then the aim is to find a set of initial states, known as a reachable set (RS), such that the trajectory of the system from this set can reach the target set in a given time [13, 14]. In some studies, the trim set (the set of the trimmed states) of the aircraft is considered as the target set, and the flight safety envelope is defined as the RS of the trim set [11, 15–17].

In the past few decades, several numerical methods have been proposed to RSs, the most commonly used of which is the level set method [13, 14, 18]. The level set method computes RS by numerically solving a Hamilton-Jacobi (HJ) partial differential equation (PDE) for a time-varying scalar function, and the RS is formulated as the zero sublevel set of the solution [13, 14]. On the basis of this principle, a number of sophisticated toolkits have been developed [19, 20], which have greatly facilitated the research process of flight safety envelope estimation [11, 15–17].

Prolonged manoeuver overload can cause structural damage to the aircraft and adversely affect the pilot physically and psychologically, thus affecting the pilot’s piloting skills [21]. Some research data suggest that a long-lasting manoeuver overload has a significant effect on visual response time [22]. This effect undoubtedly degrades flight safety.

However, the classical reachability analysis lacks a mechanism to discuss the effect of manoeuver overload on the flight safety envelope. The classical definition of RS is concerned with whether a system’s trajectory can reach the target set in finite time. Fortunately, literature [23] proposes a more generalised definition of RS, which is referred to as cost-limited reachable set (CRS), and proposes a method based on dynamic programming to compute it. In the context of CRS, one can specify a running cost function for the system and then discuss whether the state of the system can reach the target set before the performance index (time integral of the running cost function) grows to the admissible cost.

In summary, the main contributions of this paper are as follows:

• For the first time, manoeuver overload is introduced into the computation of flight safety envelopes based on reachability analysis. To achieve this, manoeuver overload is considered as a part of the running cost function and the CRS is adopted as the flight safety envelope.

• In some studies about flight safety envelope, due to the limitations of the level set method, the irregular trim set is replaced by a rectangle when computing the RS. The dynamic programming-based method proposed in Ref. [23] can easily handle irregular target sets, and benefited from this method, this paper does not make any simplifications for irregular trim sets.

• In the calculation of the CRS, the running cost function is set as a weighted sum of the manoeuver overload factor and the time consumption. We quantitatively analyse the influence of the weight of the manoeuver overload factor on the flight safety envelope.

The remainder of this paper is organised as follows. In Section 2, a dynamic model of the longitudinal motion of an aircraft is presented. Section 3 presents the introduction of RS and CRS. In Section 4, the trim set is computed, and the running cost function is constructed. Then, the influences of the manoeuver overload factor on the flight safety envelope and the optimal control inputs are analysed. Lastly, a brief conclusion is presented in Section 5.

2.0 Longitudinal motion of an aircraft

In this paper, flight safety envelopes are defined as CRSs of trim sets. However, due to the curse of dimensionality, the computational effort in computing these sets is enormous. Some studies have decoupled the aircraft dynamics model into two subsystems, the lateral and the longitudinal ones, based on the symmetry assumption of the aircraft when computing the flight safety envelopes. This paper focuses on the longitudinal subsystem.

Figure 1. Ground coordinate system.

The state variables of the longitudinal dynamics model of the aircraft include velocity v, angle of approach $\alpha $ , pitch angle $\theta $ and pitch angle velocity q, as shown in Fig. 1.

The dynamic equation of the longitudinal subsystem is:

(1) \begin{align}\begin{array}{l}{\dot s = {{\left[ {\begin{array}{c}{\dot v,\dot \alpha ,\dot q,\dot \theta }\end{array}} \right]}^{\rm{T}}} = }\\ {}{\left[ {\begin{array}{c}{\dfrac{1}{m}\left( {T\cos \alpha - \dfrac{1}{2}\rho {v^2}S{C_D}(\alpha ,{\delta _f},{\delta _e}) - mg\sin (\theta - \alpha )} \right)}\\ {q + \dfrac{1}{{mv}}\left( { - T\sin \alpha - \dfrac{1}{2}\rho {v^2}S{C_L}(\alpha ,{\delta _f},{\delta _e}) + mg\cos (\theta - \alpha )} \right)}\\ {\dfrac{1}{{2{I_{yy}}}}\rho {v^2}S\bar c{C_m}(\alpha ,q,{\delta _e})}\\ q\end{array}} \right]}\\ {}{ = f(s,u),}\end{array}\end{align}

where $s = {\left[ {v,\alpha ,q,\theta } \right]^{\rm{T}}}$ and $u = {\left[ {{\delta _f},{\delta _e}} \right]^{\rm{T}}}$ are the state and control input of the aircraft, respectively. T is the thrust, ${I_{yy}}$ and m are the inertia moment about y axis and mass of the aircraft. $\rho $ is air density ${\delta _f}$ and ${\delta _e}$ are flap and elevator deflections respectively, ${C_L},{C_D},{C_m}$ are the aerodynamic coefficients. The key parts of such a dynamic model are the aerodynamic and structure parameters of the aircraft. In the current research, the velocity is maintained as a constant by adjusting the thrust, and the following equations describe the aerodynamic coefficients:

(2) \begin{align}\left\{ {\begin{array}{l}{{C_L} = {C_{L0}} + {C_{L\alpha }}\alpha + {C_{L{\delta _f}}}{\delta _f} + {C_{L{\delta _e}}}{\delta _e}}\\[5pt]{{C_D} = {C_{D0}} + {C_{D\alpha }}\alpha + {C_{D\alpha 2}}{\alpha ^2} + {C_{D{\delta _f}}}{\delta _f} + {C_{D{\delta _e}}}{\delta _e}}\\[5pt]{{C_m} = {C_{m0}} + {C_{m\alpha }}\alpha + {C_{mq}}\frac{{q\bar c}}{{2v}} + {C_{m{\delta _e}}}{\delta _e}}\end{array}.} \right.\end{align}

The parameters in Equations (1) and (2) are summarised in Table 1.

3.0 RS and CRS

In this section, we introduce the RS and the CRS, and make a comparison between the definitions of them.

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 [23] 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.0 Computation of longitudinal flight safety envelope

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 [24], 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}$ [11], and the length of the time step is proportional to the grid size in order to ensure numerical stability [25]. 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.

5.0 Conclusions

This paper discusses the influence of manoeuver overload on the flight safety envelope on the basis of reachability analysis. The longitudinal flight safety envelope is computed as a CRS of the trim set. The CRS is characterised by a sublevel set of a value function, which is computed by a method based on dynamic programming. When computing the value function, the running cost is set as a weighted sum of time consumption and overload ratio. Several cases where the value of the weight of the overload ratio is within the interval $[0,1]$ are verified. The results show that the size of the flight safety envelope decreases with increasing weight of the manoeuver overload ratio. Increasing the weight of the overload ratio from 0 to 1 decreases the size of the flight safety envelope by about 84%.

At this moment, only the longitudinal flight safety envelope is considered. Future work will extend to a full model. It should be noted that, in this case, the flight safety envelope in a higher dimensional state space should have to be computed. From the complexity analysis, it is clear that the complexity of both the method used in this paper and the level set method contain terms proportional to the number of grids, which means that the time and memory consumption of these methods grow exponentially with the dimensionality of the state space. Therefore, a more computationally efficient method and a more efficient storage method have to be developed to analyse a full aircraft model.