Flight safety warning of iced aircraft based on reachability analysis and fuzzy inference

Abstract

A flight safety warning method based on reachability analysis and fuzzy inference is proposed against aircraft icing. A nonlinear model of iced aircraft with uncertainty is proposed based on the existing research results on the effects of icing and the uncertainty in icing detection. To deal with the uncertainty caused by icing, reachability analysis is used to estimate the safe flight envelope of iced aircraft. On this basis, fuzzy inference is employed for flight safety warning which can be used to enhance the pilot’s situational awareness in icing encounters. Simulations of the GTM (Generic Transport Model) aircraft show that, the proposed method has the potential to further increase the flight crew awareness about the risk of losing control in flight under icing.

Keywords

Aircraft icing safe flight envelope reachability analysis fuzzy inference flight safety warning

Introduction

Icing alters the aerodynamic characteristics of the aircraft, which in turn alters the flight envelope. The aviation industry gives immense importance to the harmful effects of icing and developed many anti-icing and deicing techniques. However, some serious aviation accidents induced by icing are still reported occasionally. The statistics show that many severe aviation accidents can be traced back to a similar cause, that is, aircraft exceeding the safe flight envelope.¹ Icing is a major cause of aircraft exceeding the safe envelope and even loss of control in flight (LOC-I).² Under icing conditions, the safe envelope of the aircraft will shrink. If the pilot lacks the awareness of the shrunk safe envelope and still controls the aircraft according to the normal envelope, an accident would likely happen in the worst-case scenario. Therefore, it is necessary to enhance the pilot’s situational awareness of the flight envelope under icing condition. This involves icing effect estimation, safe envelope estimation, and safety warning.

Presently, most modern planes are equipped with constraint and protection systems. However, the current types of protections are static. They deter pilot inputs from exceeding certain predefined limits. However, icing significantly changes the aircraft’s aerodynamics performance and other parameters. Compared with other adverse conditions, it is more challenging to deal with the impact of icing on flight safety. Icing may occur at different locations onboard, and it has different severity. The occurrence of icing will have an adverse effect on the aircraft flight ability or performance, such as degrading aircraft performance and handling characteristics, and even leading to aircraft LOC-I.³ In addition, icing conditions are difficult to detect in flight practice.

It is a challenge to obtain the safe flight envelop for aircraft. Most methods use the parameters of the physical model to obtain the flight envelope bounds or the limitation of commands.^4,5 Various univariate or bivariate flight envelopes have been extensively studied. These envelopes typically contain flight speed, altitude, angle of attack, sideslip angle, etc.⁶ The coupling between these parameters, as well as the pilot’s operational input, is also critical for flight safety. Therefore, a multi-dimensional safe envelope should be constructed by considering flight state variables and maneuvering control variables. In Helsen et al.,¹ the safe envelope is defined as the intersection of three envelopes, namely the dynamics envelope, the structural and comfort envelope, and the environmental envelope. The dynamic envelope reflects the dynamics and kinematic limits of aircraft and is an important envelope that is closely related to flight safety. Reachability analysis is a direct safety verification method that can deal with multidimensional system models. Therefore, reachability analysis has been tried to estimate the multidimensional dynamic envelope of an aircraft.^1,7 The neural network-based methods have also been employed to achieve flight envelope bounds estimation. In Xie et al.,⁸ the safe flight envelope bounds are updated adaptively based on the safety critical information obtained by a sparse recurrent fuzzy neural network. However, these methods are more demanding on computing resources.

Estimating the safe flight envelope in icing encounter is a more difficult task. Flight envelope estimation for iced aircraft is a challenge. The Icing Contamination Envelope Protection (ICEPro) system⁹ was designed and implemented to identify degradations in airplane performance and flying qualities resulting from ice contamination and provide safe flight-envelope cues to the pilot. ICEPro uses a of a-priori information and real-time aerodynamic estimations to provide safe limits of the flight envelope during in-flight icing encounters. In Zhang et al.,¹⁰ a database-driven safe flight envelope protection method is proposed for abnormal cases like structural damage and icing occur. Based on the information given by damage classification, the flight envelopes are explicitly retrieved and processed online from the database. Therefore, this method requires estimating and storing envelopes under multiple icing conditions. This paper presents a new method to predict the safe flight envelope and perform safety monitoring in icing encounter. The proposed method introduces the effect of icing into the system as an uncertainty parameter. as long as the icing uncertain parameter can be estimated, the safe flight envelope can be given systematically, which makes it possible to estimate the safe flight envelope under icing conditions online and eliminates the need to store large amounts of envelopes beforehand.

The flight safety warning system evaluates the safety of the flight state and feeds safety information back to the pilot in some way. In Van Baelen et al.,¹¹, the risk-evolution analysis is used determine whether a flight state is safe or not, and the safety of the whole flight space is presented in the form of a cloud chart, which provides vivid and predictive information for the pilot by means of computational flight dynamics. In Van Baelen et al.,¹¹ a haptic feedback system is designed by using force feedback through the control device to provide intuitive information on the state of the aircraft relative to the flight envelope protection.

In this paper, a flight safety warning method for iced aircraft is proposed based on reachability analysis and fuzzy inference. The effect of icing on aircraft dynamics is described by an uncertain parameter. To deal with the uncertainty caused by icing, reachability analysis is used to estimate the safe flight envelope. The estimation of the envelope is obtained by solving the Hamilton-Jacobi equation, which requires an optimal solution for parameters such as control inputs and icing effects. Then the Hamilton-Jacobi equation can be efficiently solved within the framework of level set method. To improve safety of iced aircraft, a flight safety warning method is proposed based on fuzzy inference. The main input to fuzzy inference is the physically meaningful variables constructed in the flight envelope. The contributions are as follows:

(1) The use of an uncertainty parameter to describe the effects of icing on aircraft. Consequently, uncertainties due to changes in meteorological conditions, the working conditions of the anti-icing system and icing detection can be taken into account.

(2) The reachability analysis is used to obtain the safe flight envelop of iced aircraft. Consequently, the effect of the uncertainty caused by icing on the flight envelope can be theoretically handled. Compared with univariate or bivariate flight envelopes, the flight envelope obtained by the proposed method is stored using a function defined in state space, which can be used to construct variables that characterize flight safety.

(3) The designed fuzzy inference system can generate signals containing the information of flight safety. This information is not simply derived from a comparison of control input limits and actual control inputs. Instead, it is inferred from the relative relationship between the flight state and the flight envelope in the state space. The input variables of the fuzzy inference system are constructed from the safe flight envelope in the state space. These variables not only reflect the safety of the current flight state, but also have a certain degree of predictability.

Longitudinal nonlinear model of iced aircraft

Icing changes the shape of aircraft and affects aerodynamic characteristics. As a prerequisite for flight envelope estimation and flight safety warning, a dynamic model of the aircraft in icing encounter must first be established

Longitudinal nonlinear dynamics model of GTM aircraft

GTM aircraft is a general purpose transport aircraft model developed by NASA (National Aeronautics and Space Administration) to study loss of control.^12,13 Equation (1) can describe the longitudinal dynamics model.¹⁴

{\begin{matrix} \overset{\cdot}{V} & = (F_{X} \cos α + F_{Z} \sin α) / m \\ \overset{\cdot}{α} & = q + (- F_{X} \sin α + F_{Z} \cos α) / (mV) \\ \overset{\cdot}{θ} & = q \\ \overset{\cdot}{q} & = M / I_{y} \end{matrix}

(1)

Where $V$ is flight speed, $α$ is the angle of attack, $θ$ is pitch angle, $q$ is pitch rate, $m$ is aircraft mass, $I_{y}$ is the moment of pitch inertia. $F_{X}$ , $F_{Z}$ , and $M$ are given by equation (2).

{\begin{matrix} F_{X} = \bar{q} S C_{X_{T}} (α, \hat{q}, δ_{e}) + T_{X} - mg \sin θ \\ F_{Z} = \bar{q} S C_{Z_{T}} (α, \hat{q}, δ_{e}) + T_{Z} + mg \cos θ \\ M = \bar{q} S \bar{c} C_{m_{T}} (α, \hat{q}, δ_{e}) \\ + \bar{q} S C_{X_{T}} (α, \hat{q}, δ_{e}) (z_{cgref} - z_{cg}) \\ - \bar{q} S C_{Z_{T}} (α, \hat{q}, δ_{e}) (x_{cgref} - x_{cg}) \\ + T_{X} (z_{eng} - z_{cg}) - T_{Z} (x_{eng} - x_{cg}) \end{matrix}

(2)

Where $\bar{q} = 0.5 ρ V^{2}$ is the dynamic pressure, $ρ$ is the air density, $S$ is the reference area of the airfoil, $\bar{c}$ is the mean aerodynamic chord, $\hat{q} = q \bar{c} / (2 V)$ , $δ_{e}$ is the elevator angle, and $g$ is the gravitational acceleration. $C_{X_{T}} (α, \hat{q}, δ_{e})$ , $C_{Z_{T}} (α, \hat{q}, δ_{e})$ and $C_{m_{T}} (α, \hat{q}, δ_{e})$ are the total axial force coefficient on the axis, total axial force coefficient on the axis, and total pitching moment coefficient, respectively. While $T_{X}$ and $T_{Z}$ are the components of engine thrust $T$ :

{\begin{matrix} T_{X} & = k_{X} T \\ T_{Z} & = k_{Z} T \end{matrix}

(3)

Where $k_{X}, k_{Z}$ are the components of the cosine vector of the engine direction. The definition of other parameters and their values can be found in the GTM Aircraft open source simulation package, GTM DesignSim¹⁵ released by NASA.

The model used in this paper is compared with GTM DesignSim. $fx$ in Figure 1 represents the state curve produced by the present model, and gtm_design shows the longitudinal state curve of the GTM DesignSim.

Figure 1.

The comparison of state curve.

Figure 1 reveals that the model used in this paper is consistent with the longitudinal dynamic characteristics of the GTM. The difference between the two state curves is caused by the polynomial fitting of some aerodynamics coefficient tables, presented in this paper. The reasons and details for using polynomial fitting are explained in the next section. In Figure 1, the small difference in thrust is caused by changes in flight altitude and speed.

Uncertainty description of icing effects

The icing effect model proposed by Bragg et al. is widely used.¹⁶ In this model, the icing effects on aircraft aerodynamics coefficient are expressed as:

C_{(A) iced} = (1 + η_{ice} {k'}_{C_{A}}) C_{(A)} .

(4)

Where $C_{(A) iced}$ is a certain aerodynamic coefficient of iced aircraft, and $C_{(A)}$ is the corresponding aerodynamic coefficient of clean aircraft. The relationship between $C_{(A) iced}$ and $C_{(A)}$ is determined by $η_{ice} k'_{C_{A}}$ , where $η_{ice}$ reflects the severity of icing, and the parameter $k'_{C_{A}}$ is associated with the ice type and the structural parameters of aircraft. For the same type of aircraft, the $k'_{C_{A}}$ value varies under different icing conditions, reflecting the icing effects on specific aerodynamics characteristics. This parameter is usually estimated by wind tunnel tests. Therefore, it is difficult to directly study the icing effects on safe flight envelopes using the Bragg model. In equation (4), the severity of icing is regarded as an observable variable with an exact value. The idea is to obtain the severity of icing through some observation or estimation method, and then modify the aerodynamic coefficient of clean aircraft. However, for the following reasons, it is not sufficient to use exact value of the severity of icing:

(1) Aircraft icing is a dynamic development process, and meteorological conditions continue to change. In addition, the anti-icing/deicing system works intermittently. When it works, the severity of icing will be reduced, and when it does not work, the severity of icing may be increased.

(2) The severity of icing is difficult to detect and quantify accurately. The ice detection results obtained in practice are often uncertain.

(3) The effects of icing on aircraft aerodynamic characteristics are complex. The icing effect estimation model can only describe the general trend of iced aircraft aerodynamic characteristics. Therefore, there is uncertainty in the model.

Considering these reasons, we believe that it is reasonable to use the uncertainty model to describe the effects of icing. Lampton et al. proposed an estimation method of iced aircraft model based on flight test and wind tunnel test data, in which the aerodynamic coefficient is described by a range of values.³ This method takes the icing effect as a parameter of the system, which can easily estimate the aerodynamic coefficients of different aircraft under icing conditions. In this paper, the Bragg model and Lampton method are combined, and an uncertainty description method of icing effects is developed. The Bragg model can be expressed as:

C_{(A) iced} = C_{(A)} + C_{(A)} η_{ice} k'_{C_{A}} .

(5)

Where the two parameters $η_{ice}$ and $k'_{C_{A}}$ in the second part $C_{(A)} η_{ice} k'_{C_{A}}$ are difficult to obtain accurately in practice. Therefore, in this paper, the icing effects on aircraft aerodynamics are expressed as an uncertainty model as shown in equation (6).

C_{(A) iced} = C_{(A)} + C_{{(A)}_{d}} d .

(6)

Where $d$ is the uncertainty caused by icing, the exact value of which is uncertain, but it stays in the range of $d \in D$ . The $C_{{(A)}_{d}}$ is the influence coefficient of icing uncertainty. The coefficient $C_{{(A)}_{d}}$ can be estimated by the range of coefficient $C_{(A)}$ after icing.

The model in equation (6) is an attempt to describe the icing uncertainty. In the optimization problem, the uncertainty parameter $d$ will make the safety set as small as possible. This leads to a degree of conservatism that is generally acceptable for flight safety. The most significant difference between equation (6) and (5) is that, the uncertainty in the detection or estimation of icing severity is introduced in equation (6). Therefore, it is necessary to adopt a method that can deal with this icing uncertainty to estimate the safety set of iced aircraft.

The GTM aerodynamic data is given in the form of a lookup table. Each aerodynamics coefficient is composed of a series of look-up tables. By using uncertainty description, the icing effect factors are superimposed into each aerodynamics coefficient, and then the aerodynamics coefficient of the iced aircraft can be expanded as:

{\begin{matrix} C_{X_{T}} (α, \hat{q}, δ_{e}) = C_{X} (α) + C_{X} (α, δ_{e}) \\ + C_{X} (α, \hat{q}) + C_{X d} d \\ C_{Z_{T}} (α, \hat{q}, δ_{e}) = C_{Z} (α) + C_{Z} (α, δ_{e}) \\ + C_{Z} (α, \hat{q}) + C_{Z d} d \\ C_{m_{T}} (α, \hat{q}, δ_{e}) = C_{m} (α) + C_{m} (α, δ_{e}) \\ + C_{m} (α, \hat{q}) + C_{md} d \end{matrix}

(7)

Where $C_{X} (α)$ , $C_{X} (α, δ_{e})$ , and $C_{X} (α, \hat{q})$ correspond to the look-up tables of aerodynamics coefficients along the X-axis of aircraft.¹⁵ Similarly, $C_{Z} (α)$ , $C_{Z} (α, δ_{e})$ , $C_{Z} (α, \hat{q})$ , $C_{m} (α)$ , $C_{m} (α, δ_{e})$ , and $C_{m} (α, \hat{q})$ correspond to corresponding look-up tables. $C_{Xd}$ , $C_{Zd}$ , and $C_{md}$ are the uncertainty influence coefficients of icing.

The model in equation (6) separates the icing uncertainty and control inputs. Equation (7) assumes that the aerodynamic coefficients associated with the control inputs do not explicitly contain the icing uncertainty. This results in no explicit cross-term between icing uncertainty and control inputs. The models in equations (6) and (7) do not apply to cases where icing happens on the control surfaces.

Flight safety set

Safety set of nonlinear dynamics systems

Consider the nonlinear system with external inputs and uncertain parameters as shown in equation (8):

\overset{\cdot}{x} = f (x, u, d) .

(8)

Where $x \in R^{n}$ is state variables, $u \in U \subset R^{m}$ is control input, $d \in D \subset R^{q}$ is the uncertainty parameter caused by icing. Vector field $f$ is the Lipschitz continuous function. The solution trajectory of the system (8) is noted as $ξ_{f} (t; x_{0}, t_{0}, u (\cdot), d (\cdot))$ , where $x_{0}$ is the initial state, and $t_{0}$ the initial time.

For safety reasons, the state variables of the aircraft should be properly limited. Suppose the state variables of aircraft is constrained within the set $K \subset R^{n}$ . The flight envelope can be represented as a set $S$ in which there is always a feasible control $u \in U$ , guaranteeing that the flight state does not exceed the limit $K$ under icing conditions $d \in D$ . Set $S$ can be defined as equation (9):

\begin{matrix} S = V_{f} (K, [0, T]) : = {x_{0} \in R^{n} | \forall d (\cdot) \in D, \\ \exists u (\cdot) \in U, \forall τ \in [0, T], \\ ξ_{f} (τ; x_{0}, 0, u (\cdot), d (\cdot)) \in K} . \end{matrix}

(9)

Set $S$ can be transformed into an optimization control problem, which is obtained by the viscosity solution of the Hamilton-Jacobi equation.¹⁷ Let $ϕ (x, t)$ be the viscosity solution of the Hamilton-Jacobi, then:

\begin{matrix} \frac{\partial ϕ (x, t)}{\partial t} + min [0, H (x, p)] = 0 \\ ϕ (x, t_{terminal}) = ϕ_{0} (x) . \end{matrix}

(10)

Where $t_{terminal}$ is terminal time, $p = \partial ϕ (x, t) / \partial x$ is the gradient of solution function $ϕ (x, t)$ , and the geometric meaning of $H (x, p)$ is the optimal value of the inner product of the system vector field.

H (x, p) = sup_{u (\cdot) \in U} inf_{d (\cdot) \in D} p^{T} f (x, u, d) .

(11)

The function of the control variable $u$ is to keep the aircraft within the target set $K$ to ensure the safety of the aircraft. Icing uncertainty $d$ acts in the opposite way to the control variable $u$ . It will make the flight status as far as possible beyond the limit K. Therefore, in the optimization problem, the optimization of icing uncertainty $d$ is opposite to that of $u$ , which is like the pursuit-evasion problem.¹⁷ The safety set $S$ can be expressed as

S = V_{f} (K, [0, t_{terminal}]) = {x \in R^{n} | ϕ (x, 0) > 0} .

(12)

Flight safety set calculation

The key to calculating flight safety set is to determine the boundary conditions $ϕ_{0} (x)$ and seek the optimal solution of equation (11). Suppose the flight states are confined to a rectangular area $K = [V_{min}, V_{max}] \times [α_{min}, α_{max}] \times [θ_{min}, θ_{max}] \times [q_{min}, q_{max}]$ , then the boundary conditions $ϕ_{0} (x) : R^{3} \to R$ can be assigned as:

\begin{matrix} ϕ_{0} (x) = min {x_{1} - V_{min}, V_{max} - x_{1}, \\ x_{2} - α_{min}, α_{max} - x_{2}, \\ x_{3} - θ_{min}, θ_{max} - x_{3}, \\ x_{4} - q_{min}, q_{max} - x_{4}} \end{matrix}

(13)

Thus, if $ϕ_{0} (x) \geq 0$ , then $x \in K$ , otherwise $x \notin K$ . State constraints can be taken in other forms, such as hyperplane, as long as the corresponding boundary condition $ϕ_{0} (x)$ is constructed. In practice, control inputs follow their corresponding limiting ranges, like $T_{min} \leq T \leq T_{max}$ and $δ_{e_{min}} \leq δ_{e} \leq δ_{e_{max}}$ . In addition, icing uncertainty has a certain range $d_{min} \leq d \leq d_{max}$ .

In addition to the boundary conditions, the Hamilton-Jacobi equations require solving optimization problems (11) and obtaining optimal solutions $T^{*}$ , $δ_{e}^{*}$ and $d^{*}$ . The functions of $T^{*}$ and $δ_{e}^{*}$ are keeping the aircraft within $K$ as maximum as possible, while $d^{*}$ gets the aircraft away from $K$ as maximum as possible . The set $S$ obtained from the boundary conditions $ϕ_{0} (x)$ and the optimal solution $T^{*}$ , $δ_{e}^{*}$ and $d^{*}$ is the maximum control invariant set of the aircraft system corresponding to the set $K$ . Thus, for any icing uncertainty $d_{min} \leq d \leq d_{max}$ , there exists feasible control $T_{min} \leq T \leq T_{max}$ and $δ_{e_{min}} \leq δ_{e} \leq δ_{e_{max}}$ , such that the aircraft remains within the bounds $K$ . Therefore, $S$ can be regarded as the flight safety set of the aircraft,¹⁸ and the boundary of $S$ can be regarded as the safe flight envelope.

A level set algorithm can be employed to solve the viscosity solution $ϕ (x, t)$ of equation (10). In this paper, Mitchell’s level set toolbox¹⁷ was employed. The procedure of safe flight envelope calculation can be summarized in Algorithm 1.

Algorithm 1: safe flight envelope calculation
Data: The dynamics model $f (x, u, d)$ , the optimal solution $u^{}$ and $d^{}$ , the computational domain $G$ , the total evolution time $t_{terminal}$ , the target set $K$ for the initial estimate. Result: The flight safety set $S$ , the safe flight envelope $\partial S$ 1 Create the computational grid in $G$ 2 $t = t_{terminal}$ 3 $ϕ_{0} (x) \overset{(13)}{\leftarrow} K$ 4 repeat 5 $H (x, p) \leftarrow min [0, {(\nabla ϕ_{0})}^{T} f (x, u^{}, d^{})]$ 6 Determine an time increment $Δ t$ 7 $ϕ (x, t) \leftarrow \int (H (x, p), ϕ_{0} (x), Δ t)$ 8 $ϕ_{0} (x) \leftarrow ϕ (x, t)$ 9 $t \leftarrow t - Δ t$ 10 until $t = 0$ 11 $S \leftarrow {x \in R^{n} \| ϕ (x, 0) > 0}$ 12 $\partial S \leftarrow {x \in R^{n} \| ϕ (x, 0) = 0}$

Algorithm 1: safe flight envelope calculation

Data: The dynamics model

f (x, u, d)

, the optimal solution

u^{*}

and

d^{*}

, the computational domain

G

, the total evolution time

t_{terminal}

, the target set

K

for the initial estimate.
Result: The flight safety set

S

, the safe flight envelope

\partial S

1 Create the computational grid in

G

t = t_{terminal}

ϕ_{0} (x) \overset{(13)}{\leftarrow} K

4 repeat
5

H (x, p) \leftarrow min [0, {(\nabla ϕ_{0})}^{T} f (x, u^{*}, d^{*})]

6 Determine an time increment

Δ t

ϕ (x, t) \leftarrow \int (H (x, p), ϕ_{0} (x), Δ t)

ϕ_{0} (x) \leftarrow ϕ (x, t)

t \leftarrow t - Δ t

10 until

t = 0

S \leftarrow {x \in R^{n} | ϕ (x, 0) > 0}

\partial S \leftarrow {x \in R^{n} | ϕ (x, 0) = 0}

In Algorithm 1, the safety set is determined by the limits of flight state variables, the limits of control input variables and the range of the severity of icing. Once these conditions are given, The flight safety set is determined. For a specific aircraft, it is reasonable to think that limits of flight state variables and the limits of control input variables are known. Therefore, the flight safety set will vary with the given severity of icing. If the severity of icing changes dynamically during flight, then the flight safety set also changes during flight.

For the clean aircraft, $H (x, p)$ can be simplified as:

H (x, p) = sup_{u (\cdot) \in U} p^{T} f (x, u) .

(14)

From equations (1) and (2), it is clear that there is no cross-term between control variables $T$ and $δ_{e}$ in the system equation, therefore $f (x, T, δ_{e})$ can be separated as:

f (x, T, δ_{e}) = f_{x} (x) + f_{T} (x, T) + f_{δ_{e}} (x, δ_{e}) .

(15)

Thus, $H (x, p)$ is divided into three parts that can be independently optimized:

H (x, p) = H_{T} (x, p) + H_{δ_{e}} (x, p) + H_{0} (x, p) .

(16)

Where

{\begin{matrix} H_{T} (x, p) & = sup_{T \in [T_{min}, T_{max}]} p^{T} f_{T} (x, T) \\ H_{δ_{e}} (x, p) & = sup_{δ_{e} \in [δ_{e_{min}}, δ_{e_{max}}]} p^{T} f_{δ_{e}} (x, δ_{e}) \\ H_{0} (x, p) & = p^{T} f_{x} (x) \end{matrix}

(17)

$H_{0} (x, p)$ does not explicitly contain control variables and does not need to be considered when seeking the optimal solution. For $H_{T} (x, p)$ , $f_{T} (x, T)$ can be expanded as:

\begin{matrix} f_{T} (x, T) & = T [\begin{matrix} (k_{X} \cos α + k_{Z} \sin α) / m \\ (- k_{X} \sin α + k_{Z} \cos α) / (mV) \\ 0 \\ [k_{X} (z_{eng} - z_{cg}) - k_{Z} (x_{eng} - x_{cg})] / I_{y} \end{matrix}] . \\ = T K_{f_{T}} \end{matrix}

(18)

Thus, $p^{T} f_{T} (x, T)$ is linear to $T$ , and the optimal solution can be taken as:

T^{*} = {\begin{matrix} T_{min}, & if p^{T} K_{f_{T}} \leq 0 \\ T_{max}, & otherwise . \end{matrix}

(19)

For $H_{δ_{e}} (x, p)$ , $f_{δ_{e}} (x, δ_{e})$ can be expanded as:

\begin{matrix} H_{δ_{e}} (x, δ_{e}) = sup_{δ_{e} \in [δ_{e_{min}}, δ_{e_{max}}]} \\ \bar{q} S {k_{CX} C_{X} (α, δ_{e}) + k_{CZ} C_{Z} (α, δ_{e}) + k_{Cm} C_{m} (α, δ_{e})} . \end{matrix}

(20)

Where

{\begin{matrix} k_{CX} = p_{1} \cos α / m - p_{2} \sin α / (mV) \\ + p_{4} (z_{cgref} - z_{cg}) / I_{y} \\ k_{CZ} = p_{1} \sin α / m + p_{2} \cos α / (mV) \\ - p_{4} (x_{cgref} - x_{cg}) / I_{y} \\ k_{Cm} = p_{4} \bar{c} / I_{y} \end{matrix}

(21)

To facilitate the search for the optimal solution $δ_{e}^{*}$ , polynomial fitting is performed on the aerodynamics data associated with $δ_{e}$ . The aerodynamics data tables provided by the NASA GTM Simulink model cover the angle of attack in $[- 5 °, 85 °]$ and elevator in $[- 30 °, 20 °]$ . Flight safety analysis is usually performed at the angle of attack below $20 °$ . Hence, the aerodynamics coefficient is fitted within the range of angle of attack $[- 5 °, 20 °]$ in this paper. Within this range, $C_{X} (α, δ_{e})$ presents distinct nonlinearity. Therefore it is fitted as a quadratic polynomial.

\begin{matrix} C_{X} (α, δ_{e}) = P X_{00} + P X_{10} α + P X_{01} δ_{e} \\ + P X_{20} α^{2} + P X_{11} α δ_{e} + P X_{02} {δ_{e}}^{2} \end{matrix} .

(22)

Where $P X_{00} = 9.602 e - 4$ , $P X_{10} = - 5.558 e - 4$ , $P X_{01} = - 9.349 e - 3$ , $P X_{20} = - 0.03$ , $P X_{11} = 0.076$ and $P X_{02} = - 0.06987$ .

$C_{Z} (α, δ_{e})$ and $C_{m} (α, δ_{e})$ can be linearly fitted as:

C_{Z} (α, δ_{e}) = P Z_{00} + P Z_{10} α + P Z_{01} δ_{e}

(23)

Where $P Z_{00} = - 5.9 e - 3$ , $P Z_{10} = 6.7 e - 3$ and $P Z_{01} = - 0.3824$ .

C_{m} (α, δ_{e}) = P m_{00} + P m_{10} α + P m_{01} δ_{e}

(24)

Where $P m_{00} = - 0.0107$ , $P m_{10} = 0.02324$ and $P m_{01} = - 1.524$ .

Because $p^{T} f_{δ_{e}} (x, δ_{e})$ is the quadratic function of $δ_{e}$ , and it is known to have one and only one stagnation point ${\hat{δ}}_{e}$ :

\begin{matrix} {\hat{δ}}_{e} = - (P X_{11} α + P X_{01}) / 2 P X_{02} \\ - (K_{CZ} P Z_{01} + K_{Cm} P m_{01}) / 2 P X_{02} K_{CX} \end{matrix}

(25)

The optimal solution $δ_{e}^{*}$ includes three cases.

Case 1, $k_{CX} = 0$ , then:

{δ_{e}}^{*} = {\begin{matrix} {δ_{e}}_{min}, if (k_{CZ} P Z_{01} + k_{Cm} P m_{01}) \leq 0; \\ {δ_{e}}_{max}, otherwise . \end{matrix}

(26)

Case 2, $k_{CX} P X_{02} > 0$ , then:

{δ_{e}}^{*} = {\begin{matrix} {δ_{e}}_{min}, if {\hat{δ}}_{e} \geq ({δ_{e}}_{min} + {δ_{e}}_{max}) / 2; \\ {δ_{e}}_{max}, otherwise . \end{matrix}

(27)

Case 3, $k_{CX} P X_{02} < 0$ , then:

{δ_{e}}^{*} = {\begin{matrix} {δ_{e}}_{min}, if {\hat{δ}}_{e} \leq {δ_{e}}_{min}; \\ {\hat{δ}}_{e}, if {δ_{e}}_{min} < {\hat{δ}}_{e} < {δ_{e}}_{max}; \\ {δ_{e}}_{max}, if {\hat{δ}}_{e} \geq {δ_{e}}_{max} . \end{matrix}

(28)

After setting the boundary conditions and finding the optimal solutions $T^{*}$ and $δ_{e}^{*}$ , the level set toolbox can be employed to solve the Hamilton-Jacobi equation (10), and the viscosity solution $ϕ (x, t)$ can be obtained. In this paper, $V_{min} = 90 ft / s$ , $V_{max} = 120 ft / s$ , $α_{min} = 0 \deg$ , $α_{max} = 15 \deg$ , $θ_{min} = - 15 \deg$ , $θ_{max} = 30 \deg$ , $q_{min} = - 40 \deg / s$ , $q_{max} = 40 \deg / s$ , $T_{min} = 0 lbf$ , $T_{max} = 30 lbf$ , $δ_{e_{min}} = - 30 \deg$ , $δ_{e_{max}} = 20 \deg$ .

The longitudinal dynamics model of aircraft has four state variables, so the system state space is four-dimensional and the solution function $ϕ (x, 0)$ is defined in $R^{4}$ . For better visualization, the solution function is sliced along the $α$ axis, and the sliced section of $ϕ (x, 0)$ at $α = 10 °$ is shown in Figure 2.

Figure 2.

The slice of solution function $ϕ (x, 0)$ at $α = 10 °$ .

A three-dimensional subspace is obtained from the slice of $ϕ (x, 0)$ at $α = 10 °$ . Figure 2 shows an eighth of this three-dimensional subspace. From Figure 2, it can be observed that the value of the solution function $ϕ (x, 0)$ reaches its maximum in the center of the three-dimensional subspace, which is positive, while the value gradually decreases and becomes negative when it extends to the edge. The interface between positive and negative is the boundary of the flight safety set, that is, the sefe flight envelope, as shown in the black curve in Figure 2. As shown in Figure 3, the safe flight envelope is extracted and rendered into a red transparent curved surface.

Figure 3.

The slice of the safe flight envelope.

Figures 2 and 3 present the flight safety set and the envelope of the calculation up to the terminal time $t_{terminal} = 5 s$ . When calculated to terminal time $t_{terminal} = 5 s$ , state points ${x | ϕ (x, 0) \geq 0}$ within the boundary of the safety set meet the condition $H (x, p) \geq 0$ , so the safe envelope does not change anymore. This means that, for state points inside the safety set, there exist feasible control $u \in U$ which can keep the aircraft within the state limit $K$ . The terminal time $t_{terminal} = 5 s$ does not imply the stability of the system within 5 s. It means that for the state points outside the safety set, within the limit of the control variable, no matter what control is adopted, the aircraft will exceed the state limit $K$ in at most 5 s, which may lead to catastrophic consequences.¹⁹

Verification of safety set boundary

For the state points in the safety set, permitted control can always be found to keep the flight state within the target set $K$ . The optimal solution $u^{*}$ is a feasible control. A simple switching strategy can be designed, that is, when the flight state is inside the safety set, the control command is given by the pilot or the autopilot, and when the flight state reaches the envelope, it is switched to $u^{*}$ control.¹⁸ Around the envelope, the gradient $p$ of $ϕ (x, 0)$ points inside the safety set, and $u^{*}$ is optimized along the direction of the gradient. Consequently, $u^{*}$ will move the system as far as possible into the safe set. Therefore, there must be $p^{T} f (x, u^{*}) \geq 0$ near the envelope. This can be used to verify the correctness of the safety set. In the slice corresponding $α = 10 °$ , the $p^{T} f (x, u^{*})$ values of state points on the envelope are shown in Figure 4. As can be seen in Figure 4, the state points on the envelope satisfy $p^{T} f (x, u^{*}) \geq 0$ .

Figure 4.

The scatter of $p^{T} f (x, u^{*}) \geq 0$ on the envelope.

Moreover, when the system state is outside the safety set $V_{f} (K, [0, t_{terminal}])$ , the state overrun will be inevitable, and it will exceed the target set $K$ at most at time $t_{terminal}$ . This can also be used for verification to ensure that points in the safety set do not fall outside the safety incorrectly. Set the initial state to $x_{01} = [115, 10 - 10, - 10]$ . $x_{01}$ is outside the safety set, but within the target set $K$ . Take control $u^{*}$ , that is, try to keep the aircraft state in the target set $K$ . The state response curve is shown in Figure 5. It can be seen that the flight state exceeds the limit at about $0.7$ seconds. The simulation results for more points outside the safety set show that they will cause the flight state to exceed the limit $K$ in no more than $t_{terminal} = 5$ seconds. This indicates that points within the safety set are not incorrectly excluded.

Figure 5.

The state response curve of $x_{01}$ with the control $u^{*}$ .

On the other hand, for the state points within the safety set, there is always a feasible control to keep the flight state in the target set $K$ . But, if it is not controlled properly, the flight state will pass through the envelope and then exceeds the state limit. When the flight state reaches the envelope, $u^{*}$ will control the trajectory to slide on the boundary or turn back inside the safe envelope. For example, balance the aircraft to a steady state $\hat{x} = [114.9, 7.02, 10.6, 0]$ , the required engine thrust is about $9.4 lbf$ , and the elevator angle is about $- 0.7 °$ . Set the initial state to $x_{02} = [115, 7.02, - 3, 0]$ , and $x_{02}$ is still within the safety set. Under the condition that engine thrust and elevator angle are kept at the balanced values. The state response curve is shown in Figure 6.

Figure 6.

The state response curve of $x_{02}$ when the control is equal to trim value.

As shown in Figure 6, the system state was traversed from within the safety set and exceeded state limits due to lack of proper controls. If the control is switched to $u^{*}$ when the flight state reaches the boundary of the safety set, the state curve is changed, as shown in Figure 7.

Figure 7.

The state response curve of $x_{02}$ with the control $u^{*}$ .

As shown in Figure 7, $u^{*}$ effectively avoids flight state overrun. Control quantity and system real-time $ϕ (x, 0)$ value produced by enabling $u^{*}$ control are shown in Figure 8.

Figure 8.

The real-time value of control variables and $ϕ (x, 0)$ .

From Figure 8, it can be observed that $ϕ (x, 0) = 0$ at about $t = 0.5 s$ , which means the flight state reaches the envelope, at which point the control $u^{*}$ is enabled. Thereafter, the value $ϕ (x, 0) = 0$ stays close to zero, that is, the state trajectory slides on the envelope. It can be seen that $u^{*}$ maintains the flight state within the safety set to the greatest extent. However, Figure 8 reveals the deficiency of this simple switching method, that is, the control variable exhibits severe oscillation. Moreover, the flight sate slides on the envelope for a long time, and there is a risk of being disturbed and exceeding the limit. These deficiencies are due to the untimely initiation of envelope protection. The flight safety warning method proposed in this paper provides the pilot with the safety information in real-time by issuing a warning before the state reaches the envelope. With the real-time flight safety information, the pilots can conduct proper maneuvering to ensure that the aircraft does not exceed the state limits.

Damage of icing to flight safety

Flight safety set of iced aircraft

According to equations (6) and (7), there is no cross-term between icing uncertainty $d$ and control variables $T$ and $δ_{e}$ , so the effect of $d$ can be separated individually, that is, $f_{d} (x, d)$ .

\begin{matrix} f_{d} (x, d) = \bar{q} Sd K_{f_{d}} . \end{matrix}

(29)

Where

K_{f_{d}} = [\begin{matrix} (C_{Xd} \cos α + C_{Zd} \sin α) / m \\ (- C_{Xd} \sin α + C_{Xd} \cos α) / (mV) \\ 0 \\ [\bar{c} C_{md} + C_{Xd} (z_{eng} - z_{cg}) - C_{Xd} (x_{eng} - x_{cg})] / I_{y} \end{matrix}] .

(30)

Therefore, $p^{T} f_{d} (x, d)$ is linear to $d$ .

\begin{matrix} H_{d} (x, p) = inf_{d \in [d_{min}, d_{max}]} p^{T} f_{d} (x, d) . \end{matrix}

(31)

So the optimal solution of $d$ is:

d^{*} = {\begin{matrix} d_{max}, & if p^{T} K_{f_{d}} \leq 0; \\ d_{min}, & otherwise . \end{matrix}

(32)

For a specific type of aircraft, the icing effect coefficient can be estimated by wind tunnel test and flight test data. According to the variation rule of aerodynamics coefficient of GTM aircraft after icing,² the relevant parameters are estimated as: $C_{Xd} = 0.03$ , $C_{Zd} = - 0.15$ , and $C_{md} = 0.01$ . Slicing at $α = 10 °$ , the flight safe envelope of iced aircraft is shown in Figure 9.

Figure 9.

The slice of the safe flight envelope of the iced aircraft.

Analysis of icing damage

In Algorithm 1, $d$ is the uncertainty parameter caused by icing. The uncertainty parameter $d$ the opposite effect on the control variable $u$ . The function of uncertainty parameter $d$ is to make the aircraft leave the target set $K$ as far as possible. Therefore, in the optimization problem, the optimization of $d$ $(inf_{d (\cdot) \in D})$ is opposite to that of $u$ $(sup_{u (\cdot) \in U})$ . The uncertainty parameter $d$ will make $H (x, p)$ as small as possible. According to Algorithm 1, this will cause $ϕ (x, t)$ to evolve toward negative values. The safety set is described by $S = {x \in R^{n} | ϕ (x, 0) > 0}$ . This means that the uncertainty parameter $d$ will shrink the safety set.

By comparing Figures 3 and 9, it is observed that the volumes of the safety sets in Figure 9 are about 15% less than that in Figure 3. The red and the blue transparent surfaces in Figure 10 are the safe flight envelope of clean aircraft and iced aircraft, respectively. It can be known that the envelope shrinks seriously in icing encounter.

Figure 10.

Comparison of the envelope of clean aircraft and iced aircraft.

From Figure 10, it can be observed that, the high-speed part of the flight safe envelope on the $V$ axis shrinks significantly under icing conditions. The main reason is that, icing reduces the lift of the aircraft, and the angle of attack is small. To compensate for the lack of lift, the flight speed needs to be increased, so the speed exceeds the limit. The low-speed part of the flight safe envelope on the $V$ axis shows considerable shrinkage. This is mainly due to the icing that reduces the lift of the aircraft, due to the large angle of attack, in order to make up for the lack of lift, the angle of attack will be further increased, which is easy to lead to the pitch angle exceeding the limit, so the low-speed part is no longer safe. In addition, the envelope of iced aircraft shows considerable shrinkage along the $q$ axis, which is caused by the change in pitching moment after icing. Slicing by $α = 2^{o}, 5^{o}, 10^{o}, 13^{o}$ , it is found that greater $α$ leads to more severe volumetric shrinkage of safety set. This shows that aircraft safety is seriously reduced at a high angle of attack in icing encounter.

It can be found that icing seriously compressed the flight safety set. If the driver is not aware of the shrinking envelope, it is easy to make mistakes and lead to flight accidents. Therefore, it is necessary to enhance pilots’ awareness of flight safety under icing conditions and implement safety warning.

Flight safety warning of iced aircraft

Flight safety is a fuzzy concept that is described by a fuzzy set. The causal relationship between flight state variables and flight safety is complex and difficult to be described by explicit functions. Therefore, the fuzzy inference is used in this paper for the safety warning of iced aircraft.

The fuzzy comprehensive evaluation method has been used in flight safety evaluation.²⁰ In the fuzzy comprehensive evaluation method, all parameters associated with flight safety (state variables and maneuver variables) are divided into several intervals according to prior knowledge. Then, the interval of each parameter is determined according to the state of the system. Finally, the comprehensive safety evaluation results are obtained by using the composition rules. The advantage of this method is that all factors can be managed directly. However, this method is highly dependent on experience and needs to store the state classification intervals of each parameter under different icing conditions.

In contrast, the fuzzy inference method used in this paper does not directly classify flight state variables or control variables. Instead, the flight safety set is estimated based on a nonlinear model with icing uncertainties, flight state variables and control variables. Based on this safety set, the safety warning variables such as generalized distance, velocity, and acceleration are constructed. The safety warning variables are used for fuzzy inference to realize safety warning. Thus, the method proposed in this paper is based on the flight dynamics mechanism rather than completely depending on experience.

Flight safety warning variables of iced aircraft

The solution function of the Hamilton-Jacobi equation provides a measure of the distance between the state point and the safe flight envelope. $ϕ (x, 0)$ can be regarded as the generalized distance from $x$ to the boundary of the safety set, and defined as:

L = ϕ (x, 0)

(33)

$L > 0$ indicates that the aircraft state is within the safety set, $L < 0$ means that the aircraft state is outside the safety set, and $L = 0$ means that the aircraft state is on the safe flight envelope. Hence, the value of $L$ is closely related to the safety degree of the flight state.

The control variable $u$ corresponds to the input from the pilot. It is a crucial factor in determining the motion trend of the aircraft, which has a significant impact on flight safety. The control variable $u$ directly affects the velocity and acceleration of the flight state moving toward the envelope. The rate of change of $L$ with respect to time is:

\dot{L} = \frac{\partial ϕ}{\partial t} + \frac{\partial ϕ}{\partial x^{T}} \dot{x} = - \min {0, H (x, p)} + \frac{\partial ϕ}{\partial x^{T}} f (x, u, d)

(34)

For the state inside the safety set $x \in S$ , $H (x, p) \geq 0$ , so in the safety set:

\overset{\cdot}{L} = \frac{\partial ϕ}{\partial x^{T}} f (x, u, d)

(35)

Equation (35) has a clear geometric meaning, that is, the rate of change of distance $L$ is equal to the projection of the rate of change of system state $f (x, u, d)$ on the solution function gradient $\nabla ϕ$ . The acceleration of the system toward the boundary of the safety set is:

\begin{matrix} \overset{\cdot\cdot}{L} = \frac{d}{dt} (\frac{\partial ϕ}{\partial x^{T}}) f (x, u, d) + \frac{\partial ϕ}{\partial x^{T}} \frac{d}{dt} (f (x, u, d)) \\ = {\overset{\cdot}{x}}^{T} \frac{\partial \nabla ϕ}{\partial x^{T}} f (x, u, d) + \frac{\partial ϕ}{\partial x^{T}} \frac{\partial f}{\partial x^{T}} \overset{\cdot}{x} \\ = f^{T} \nabla^{2} ϕ f + \nabla ϕ \frac{\partial f}{\partial x^{T}} f \end{matrix}

(36)

Thus, for an iced aircraft, $L$ , $\overset{\cdot}{L}$ and $\overset{\cdot\cdot}{L}$ can comprehensively represent the impact of of the state variables, control variables and icing factors on flight safety. The use of $L$ , $\overset{\cdot}{L}$ and $\overset{\cdot\cdot}{L}$ to provide safety warning for iced aircraft can not only reflect the safety degree of the current system state, but also have a certain degree of predictability.

Fuzzy inference system

Fuzzy inference uses membership function and implication operator to express the relationship between factors and conclusions. In fuzzy inference, input variables should be fuzzified first. The variable value $L$ directly reflects the position relationship between the flight state and envelope. When $L$ is positive, it means that the flight state is inside the safety set. When $L$ is negative, it means that the flight state is outside the safety set. Therefore, three fuzzy sets are used in this paper to describe variable $L$ , including the inner far, near, and outer far, representing that the flight state is inside the safety set and away from the envelope, and the flight state is near the envelope, and the flight state is outside the safety set and away from the envelope.

The key to deciding the membership function is to establish relevant shape parameters. Taking the trapezoidal membership function as an example, its shape is determined by the turning point of the function curve. The parameters of the relevant membership function can be determined based on the value distribution of variables on the computing grid. The value distribution of $L$ over the entire calculation grid is shown in Figure 11. The value distribution of $L$ in the safety set is shown in Figure 12

Figure 11.

The value distribution of $L$ on the calculation grid.

Figure 12.

The value distribution of $L$ inside the safe set.

Considering the value distribution of $L$ , as shown in Figures 11 and 12, this paper takes the membership function as shown in Figure 13.

Figure 13.

The membership function for $L$ .

There is some overlap between the fuzzy sets in Figure 13. The overlap between fuzzy sets is helpful to improve the smoothness and robustness of the input and output of the inference system. Similarly, the membership functions describing variables $\overset{\cdot}{L}$ and $\overset{\cdot\cdot}{L}$ are shown in Figures 14 and 15:

Figure 14.

The membership function for $\overset{\cdot}{L}$ .

Figure 15.

The membership function for $\overset{\cdot\cdot}{L}$ .

In this paper, the inference conclusion, namely flight safety, is divided into four fuzzy sets, including safety, caution, vigilance, and danger. For safety warning, usually only a fuzzy set of inference results is required, and there is no requirement for defuzzify. Therefore, in this paper, simple triangular membership functions are used to describe the safety fuzzy sets, as shown in Figure 16.

Figure 16.

The membership function for describing flight safety.

Inference rules describe the causal relationship between input and output. The nonlinear relationship between input and output is expressed by the membership function and implication and aggregate. For the safety warning of iced aircraft, the inference rules used in this paper are as follows:

The last column of the Table 1 is the weight of the rule. The configuration of inference rules and their weights in Table 1 is rough and serves to demonstrate the feasibility of the method. In practice, the experience of pilots and aircraft designers should be integrated, and the inference rules and the weights of safety warning should be configured more delicately through wind tunnel tests and flight tests.

Table 1.

The inference rules.

No.	$L$	$\overset{\cdot}{L}$	$\overset{\cdot\cdot}{L}$	$Y$	weight
1	Outer far	—	—	Danger	1
2	Near	Negative	Negative	Danger	1
3	Near	Negative	Positive	Vigilance	0.8
4	Near	Positive	Negative	Vigilance	0.5
5	Near	Positive	Positive	Caution	0.8
6	Inner far	Negative	Negative	Caution	1
7	Inner far	Negative	Positive	Safety	0.5
8	Inner far	Positive	—	Safety	1

“—” means value is arbitrary.

Iced aircraft safety warning based on fuzzy inference

Inference rules in Table 1 are graphically expressed in MATLAB (Matrix Laboratory) using the fuzzy logic toolbox, as shown in Figure 17.

Figure 17.

The fuzzy inference rules for iced aircraft safety warning.

Each line in Figure 17 corresponds to an inference rule. Each column in the first component of the rule corresponds to an input variable. The second component of the rule, namely the inference result, is a fuzzy set. In Figure 17, the adopted implication operator is the Mamdani minimum operator, the composite operator is the maximum-minimum operator, and the center of gravity method is used for defuzzify, which forms a typical Mamdani fuzzy inference system. Based on the inference rule, the safety degree of each point in the flight state space can be inferred, and the results are shown in Figure 18.

Figure 18.

Safety assessment of the flight state space.

According to the membership function, as shown in Figure 16, the larger the value of the safety variable $Y$ shows a safer flight state. In Figure 18, black and red are dangerous areas, yellow is the area requiring vigilance, cyan is the area requiring caution, and green is a safe area. The black and the red are shown in Figure 18 correspond to the areas near the envelope in Figure 2. According to the meaning of the flight safety set, Once the flight state enters the area, it is easy to cause accidents due to the state exceeding the limit. The interior of the safety set is a safe region. It can be observed that the safety assessment results in Figure 18 are reasonable.

The values of $L$ , $\overset{\cdot}{L}$ , and $\overset{\cdot\cdot}{L}$ vary in real-time with flight state and control input. According to the real-time values of $L$ , $\overset{\cdot}{L}$ , and $\overset{\cdot\cdot}{L}$ , the inference system uses inference rules to derive the safety of aircraft for safety warning. For instance, under the control of the pilot’s rapid pulling and pushing of the driving stick, the state response curve of the aircraft is shown in Figure 19.

Figure 19.

The state response curve.

From Figure 19, it can be observed that the angle of attack, pitch angle and pitch rate of the aircraft have changed significantly, especially the angle of attack, and the magnitude of the change is close to the envelope. As a result, flight safety warning output changes drastically. The real-time safety warning output is shown in Figure 20.

Figure 20.

The output of safety warning.

Distinct colors are used to indicate the safety degree of the safety warning output in Figure 20. Green, cyan, yellow, and red represent safety, caution, vigilance, and danger, respectively. The rapid and substantial changes in flight state are reflected in the safety warning output of the system. Once the flight safety warning output is assessed as dangerous, it indicates that the pilot approaches the envelope.

The flight state curve is shown in Figure 21. At first, the plane is safely flat. From $t = 2$ , due to the pilot’s violent pull on the driving stick, the flight speed of the aircraft decreases, and the angle of attack, pitch angle and pitch rate increase. Among these, the angle of attack rapidly rises until approaches the limit. During the process, the safety warning output deteriorates continuously, from safety to dangerous. The subsequent stick pushing operation leads to the start of the recovery of the aircraft state. But due to the pushing being too fierce, the angle of attack and the pitch rate change drastically, and so, flight safety turns to a vigilance state. The flight state gradually returns to safety after the control stick returns to balance.

Figure 21.

The state curves and output of safety warning: (a) flight velocity, (b) angle of attack, (c) pitch angle, and (d) pitch rate.

The curve of angle of attack is shown in Figure 21(b). It can be found that the time of the warning output as dangerous is earlier than that when the angle of attack reaches its limit. This is because of the comprehensive use of variables $L$ , $\overset{\cdot}{L}$ , and $\overset{\cdot\cdot}{L}$ for conducting safety warning so that the information of state variables and control input are integrated effectively. The warning output not only includes the real-time information of the system state but also well predicts the motion trend of the flight state. Therefore, the safety warning method proposed in this paper has certain predictability. This makes pilots aware of approaching the envelope as early as possible.

How to apply the methods proposed in this paper to aircraft will be the focus of our subsequent research work. For the following reasons, we recognize that there are many difficulties in applying the proposed methods.

(1) The control law of the aircraft has a significant impact on flight safety. Both safety warning and the flight envelope estimating aim to keep the state within a safe range, which limits control. Introducing restrictions directly into flight control laws is a matter of great care.

(2) The flight crew is responsible for all manipulation results and should therefore have the highest control authority, and control restrictions should try not to bypass the crew.

(3) Under the existing flight control computer hardware conditions and flight control law design habits, the industry is cautious about envelope estimation methods such as the one proposed in this paper.

Therefore, we tend to use safety warning and envelope estimation to provide safety-critical information to the flight crew, rather than directly playing a role in the flight control laws. In recent work, we have experimented with incorporating safety warning and envelope estimation information into the primary flight display and head-up display. Safety warning information is transmitted to the crew through flashing text and voice. Flight crews expect a concise and intuitive way to display envelope, however the envelopes estimated in this paper are multidimensional structures described by implicit functions. How to project the estimated multidimensional envelope onto the flight display is our main challenge.

Conclusions

A flight safety warning method for iced aircraft is proposed based on reachability analysis and fuzzy inference. The effect of icing on aircraft dynamics is described by an uncertain parameter. To deal with the uncertainty caused by icing, reachability analysis is used to estimate the safe flight envelope. The estimation of the envelope is obtained by solving the Hamilton-Jacobi equation, which requires an optimal solution for parameters such as control inputs and icing effects. Under icing conditions, flight envelope shrink severely. In the example of this paper, the safety boundary shrinks by 15% in icing encounters, which will significantly increases the pilot’s workload. To improve safety of iced aircraft, a flight safety warning method is proposed based on fuzzy inference. The main input to fuzzy inference is the physically meaningful variables such as generalized distance, velocity, and acceleration. These variables not only reflect the safety of the current flight state, but also have a certain degree of predictability. However, this method has limitations. The icing uncertainty parameter is optimized toward the flight state exceeding the envelope. Therefore, the safety set obtained is conservative to some extent, and there may be false alarms as the result of fuzzy inference. The method proposed in this paper is intended to alert the crew to future unsafe or otherwise undesired conditions. Therefore, it is necessary to develop display technology for presenting safety warning information and the estimated flight envelope.

Footnotes

Acknowledgements

The authors would like to thank the editors and the anonymous reviewers for their constructive comments and suggestions that have improved the quality of this paper.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the Basic Research Program of Shaanxi Province of China (2021JQ-360).

ORCID iD

Guoqiang Yuan

References

Helsen

van Kampen

de Visser

, et al. Distance-fieldsover-grids method for aircraft envelope determination. J Guid Control Dyn 2016; 39(7): 1470–1480.

Broeren

Lee

Shah

, et al. Aerodynamic effects of simulated ice accretion on a generic transport model. In International Conference on Aircraft and Engine Icing and Ground Deicing, Chicago, IL: NASA Glenn Research Center.

Lampton

Valasek

. Prediction of icing effects on the lateral/directional stability and control of light airplanes. Aerosp Sci Technol 2012; 23(1): 305–311.

Govindarajan

de Visser

van Kampen

, et al. Optimal control framework for estimating autopilot safety margins. J Guid Control Dyn 2015; 38(7): 1197–1207.

Lombaerts

Looye

Ellerbroek

, et al. Design and piloted simulator evaluation of adaptive safe flight envelope protection algorithm. J Guid Control Dyn 2017; 40(8): 1902–1924.

Tekles

Chongvisal

Xargay

, et al. Design of a flight envelope protection system for NASA’s transport class model. J Guid Control Dyn 2017; 40: 863–877.

Stapel

de Visser

Chu

, et al. Efficient methods for flight envelope estimation through reachability analysis. In AIAA Guidance, Control, and Dynamics. Number AIAA 2016-0083 in AIAA SciTech Forum, American Institute of Aeronautics and Astronautics, pp. 1–21.

Xie

Meng

, et al. Neural-network-based dynamic flight envelop protection within integrated aircraft control system. J Guid Control Dyn 2023; 46(6): 1083–1094.

Martos

Ranaudo

Norton

, et al. Development, implementation, and pilot evaluation of a model-driven envelope protection system to mitigate the hazard of in-flight ice contamination on a twin-engine commuter aircraft. Technical Report NASA/CR-2014-218320, 2014.

10.

Zhang

Huang

Chu

, et al. Database-driven safe flight envelope protection for impaired aircraft. In AIAA Scitech 2020 Forum. Orlando, FL. DOI:10.2514/6.2020-1374.

11.

Van Baelen

Ellerbroek

(René) van Paassen

, et al. Design of a haptic feedback system for flight envelope protection. J Guid Control Dyn 2020; 43(4): 700–714.

12.

Cunningham

Cox

Murri

, et al. A piloted evaluation of damage accommodating flight control using a remotely piloted vehicle. In AIAA guidance, navigation, and control conference. Guidance, Navigation, and Control and Colocated Conferences, American Institute of Aeronautics and Astronautics, 2011.

13.

Cox

Cunningham

Jordan

. Subscale flight testing for aircraft loss of control: Accomplishments and future directions. In AIAA Guidance, Navigation, and Control Conference. American Institute of Aeronautics and Astronautics (AIAA).

14.

Kwatny

Dongmo

JET

Chang

, et al. Nonlinear analysis of aircraft loss of control. J Guid Control Dyn 2013; 36(1): 149–162.

15.

Cox

. Gtm designsim. Technical report, 2016.

16.

Bragg

Hutchison

Merret

. Effect of ice accretion on aircraft flight dynamics. In: 38th Aerospace Sciences Meeting and Exhibit. American Institute of Aeronautics and Astronautics (AIAA).

17.

Mitchell

Bayen

Tomlin

. A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games. IEEE Trans Automat Contr 2005; 50(7): 947–957.

18.

Allen

Kwatny

Bajpai

. Safe set protection and restoration for unimpaired and impaired aircraft. In AIAA Guidance, Navigation, and Control Conference. AIAA 2012-4822, American Institute of Aeronautics and Astronautics (AIAA), pp. 1–9.

19.

Lygeros

. On reachability and minimum cost optimal control. Automatica 2004; 40(6): 917–927.

20.

Pei

Xue

. Flight-safety space and cause of incident under icing conditions. J Guid Control Dyn 2017; 40(11): 2983–2990.