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 [Reference Jung and Kim1, Reference Lee2]. Therefore, it is especially important to avoid flight accidents. Some studies [Reference Airplanes3] 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 [Reference Kwatny, Dongmo, Chang, Bajpai, Yasar and Belcastro4]. 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 [Reference Lyu, Gu, Zhou, Li, Shen and Zhang5, Reference Ignatyev, Khrabrov, Kortukova, Alieva, Sidoryuk and Bazhenov6], flight experiments and flight simulator-based methods [Reference Goman, Khramtsovsky and Kolesnikov7, 8] or the combination of computational fluid dynamics (CFD) and extended Kalman filter [Reference Olejnik, Dziubinski and Kiszkowiak9, Reference Menon, Sengupta, Vaddi, Yang and Kwan10]. A more effective method to estimate the flight safety envelope is reachability analysis [Reference Nabi, Lombaerts, Zhang, van Kampen, Chu and de Visser11, Reference Harno and Kim12]. 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 [Reference Lygeros13, Reference Mitchell, Bayen and Tomlin14]. 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 [Reference Nabi, Lombaerts, Zhang, van Kampen, Chu and de Visser11, Reference Zhang, de Visser and Chu15–Reference Liu, Wang, Quan, Du and Yang17].
In the past few decades, several numerical methods have been proposed to RSs, the most commonly used of which is the level set method [Reference Lygeros13, Reference Mitchell, Bayen and Tomlin14, 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 [Reference Lygeros13, Reference Mitchell, Bayen and Tomlin14]. On the basis of this principle, a number of sophisticated toolkits have been developed [Reference Mitchell19, Reference Bansal, Chen, Tanabe and Tomlin20], which have greatly facilitated the research process of flight safety envelope estimation [Reference Nabi, Lombaerts, Zhang, van Kampen, Chu and de Visser11, Reference Zhang, de Visser and Chu15–Reference Liu, Wang, Quan, Du and Yang17].
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 [Reference Biernacki, Lewkowicz, Zielinski and Wojtkowiak21]. Some research data suggest that a long-lasting manoeuver overload has a significant effect on visual response time [Reference Truszczynski, Wojtkowiak, Lewkowicz, Biernacki and Kowalczuk22]. 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 [Reference Liao, Wei and Lai23] 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. [Reference Liao, Wei and Lai23] 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:
\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:
\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:
\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:
\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:
\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:
\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:
\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:
\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 [Reference Jiang, Zhou, Gao, Yang and Jing24], i.e.:
\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:
\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.:
\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
\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.
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.
Acknowledgements
The authors gratefully acknowledge support from the Aeronautical Science Foundation of China (Grant No. 20182852021).
