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

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.


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. 1Icing is a major cause of aircraft exceeding the safe envelope and even loss of control in flight (LOC-I). 2 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 worstcase 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. 3In 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,5Various univariate or bivariate flight envelopes have been extensively studied.These envelopes typically contain flight speed, altitude, angle of attack, sideslip angle, etc. 6 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., 1 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,7The neural networkbased methods have also been employed to achieve flight envelope bounds estimation.In Xie et al., 8 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 9 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., 10 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., 11 , 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., 11 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.
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,13quation ( 1) can describe the longitudinal dynamics model. 14 Where V is flight speed, a is the angle of attack, u 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).
Where q = 0:5rV 2 is the dynamic pressure, r is the air density, S is the reference area of the airfoil, c is the mean aerodynamic chord, q = q c= 2V ð Þ, d e is the elevator angle, and g is the gravitational acceleration.
Þ 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: 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 15 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 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. 16In this model, the icing effects on aircraft aerodynamics coefficient are expressed as: 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 h ice k 0 C A , where h ice reflects the severity of icing, and the parameter k 0 C A is associated with the ice type and the structural parameters of aircraft.For the same type of aircraft, the k 0 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. 3This 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: Where the two parameters h ice and k 0 C A in the second part C (A) h ice k 0 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).
Where d is the uncertainty caused by icing, the exact value of which is uncertain, but it stays in the range of 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: Where C X (a), C X (a, d e ), and C X (a, q) correspond to the look-up tables of aerodynamics coefficients along the X-axis of aircraft. 15Similarly, C Z (a), C Z (a, d e ), C Z (a, q), C m (a), C m (a, d e ), and C m (a, 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.

Safety set of nonlinear dynamics systems
Consider the nonlinear system with external inputs and uncertain parameters as shown in equation ( 8): Where 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 j f (t; x 0 , t 0 , u( Á ), d( Á )), 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 & R n .The flight envelope can be represented as a set S in which there is always a feasible control u 2 U, guaranteeing that the flight state does not exceed the limit K under icing conditions d 2 D. Set S can be defined as equation (9): Set S can be transformed into an optimization control problem, which is obtained by the viscosity solution of the Hamilton-Jacobi equation. 17Let f(x, t) be the viscosity solution of the Hamilton-Jacobi, then: Where t terminal is terminal time, p = ∂f(x, t)=∂x is the gradient of solution function f(x, t), and the geometric meaning of H(x, p) is the optimal value of the inner product of the system vector field.
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. 17The safety set S can be expressed as Flight safety set calculation The key to calculating flight safety set is to determine the boundary conditions f 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 3 ½a min , a max 3½u min , u max 3½q min , q max , then the boundary conditions f 0 (x) : R 3 !R can be assigned as: x 2 À a min , a max À x 2 , Thus, if f 0 (x)50, then x 2 K, otherwise x 6 2 K. State constraints can be taken in other forms, such as hyperplane, as long as the corresponding boundary condition f 0 (x) is constructed.In practice, control inputs follow their corresponding limiting ranges, like T min 4T4T max and d e min 4d e 4d e max .In addition, icing uncertainty has a certain range d min 4d4d max .
In addition to the boundary conditions, the Hamilton-Jacobi equations require solving optimization problems (11) and obtaining optimal solutions T Ã , d Ã e and d Ã .The functions of T Ã and d Ã 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 f 0 (x) and the optimal solution T Ã , d Ã 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 4d4d max , there exists feasible control T min 4T4T max and d e min 4d e 4d e max , such that the aircraft remains within the bounds K. Therefore, S can be regarded as the flight safety set of the aircraft, 18 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 f(x, t) of equation (10).In this paper, Mitchell's level set toolbox 17 was employed.The procedure of safe flight envelope calculation can be summarized in Algorithm 1.
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: From equations ( 1) and ( 2), it is clear that there is no cross-term between control variables T and d e in the system equation, therefore f(x, T, d e ) can be separated as: Thus, H(x, p) is divided into three parts that can be independently optimized: Where 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: Thus, p T f T (x, T) is linear to T, and the optimal solution can be taken as: 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 ∂S 1 Create the computational grid in G 2 t = t terminal 3 f 0 (x) (13) For H d e (x, p), f d e (x, d e ) can be expanded as: Where To facilitate the search for the optimal solution d Ã e , polynomial fitting is performed on the aerodynamics data associated with d e .The aerodynamics data tables provided by the NASA GTM Simulink model cover the angle of attack in ½À58, 858 and elevator in ½À308, 208.Flight safety analysis is usually performed at the angle of attack below 208.Hence, the aerodynamics coefficient is fitted within the range of angle of attack ½À58, 208 in this paper.Within this range, C X (a, d e ) presents distinct nonlinearity.Therefore it is fitted as a quadratic polynomial.
Because p T f d e (x, d e ) is the quadratic function of d e , and it is known to have one and only one stagnation point de : The optimal solution d Ã e includes three cases.
Case 1, k CX = 0, then: Case 2, k CX PX 02 .0, then: Case 3, k CX PX 02 \ 0, then: After setting the boundary conditions and finding the optimal solutions T Ã and d Ã e , the level set toolbox can be employed to solve the Hamilton-Jacobi equation (10), and the viscosity solution f(x, t) can be obtained.In this paper, V min = 90ft=s, V max = 120ft=s, a min = 0 deg, a max = 15 deg, u min = À 15 deg, u max = 30 deg, q min = À 40deg=s, q max = 40deg=s, T min = 0 lbf, The longitudinal dynamics model of aircraft has four state variables, so the system state space is fourdimensional and the solution function f(x, 0) is defined in R 4 .For better visualization, the solution function is sliced along the a axis, and the sliced section of f(x, 0) at a = 108 is shown in Figure 2.
A three-dimensional subspace is obtained from the slice of f(x, 0) at a = 108.Figure 2 shows an eighth of this three-dimensional subspace.From Figure 2, it can be observed that the value of the solution function f(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.
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 fxjf(x, 0)50g within the boundary of the safety set meet the condition H(x, p)50, so the safe envelope does not change anymore.This means that, for state points inside the safety set, there exist feasible control u 2 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. 19

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. 18Around the envelope, the gradient p of f(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 Ã )50 near the envelope.This can be used to verify the correctness of the safety set.In the slice corresponding a = 108, 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 Ã )50.
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.
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 x = ½114:9, 7:02, 10:6, 0, the required engine thrust is about 9:4lbf, and the elevator angle is about À0:78.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.
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.
As shown in Figure 7, u Ã effectively avoids flight state overrun.Control quantity and system real-time f(x, 0) value produced by enabling u Ã control are shown in Figure 8.
From Figure 8, it can be observed that f(x, 0) = 0 at about t = 0:5s, which means the flight state reaches the envelope, at which point the control u Ã is enabled.Thereafter, the value f(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.

Flight safety set of iced aircraft
According to equations ( 6) and ( 7), there is no crossterm between icing uncertainty d and control variables T and d e , so the effect of d can be separated individually, that is, f d (x, d). Where Therefore, p T f d (x, d) is linear to d.
So the optimal solution of d is: 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, 2 the relevant parameters are estimated as: C Xd = 0:03, C Zd = À 0:15, and C md = 0:01.Slicing at a = 108, the flight safe envelope of iced aircraft is shown in Figure 9.

Analysis of icing damage
In .The uncertainty parameter d will make H(x, p) as small as possible.
According to Algorithm 1, this will cause f(x, t) to evolve toward negative values.The safety set is described by S = x 2 R n jf(x, 0) .0 f g .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.
From Figure 10, it can be observed that, the highspeed 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 a = 2 o , 5 o , 10 o , 13 o , it is found that greater a 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. 20In 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.f(x, 0) can be regarded as the generalized distance from x to the boundary of the safety set, and defined as: 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: For the state inside the safety set x 2 S, H(x, p)50, so in the safety set: 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 rf.The acceleration of the system toward the boundary of the safety set is: Thus, for an iced aircraft, L, _ L and € L can comprehensively represent the impact of of the state variables, control variables and icing factors on flight safety.The use of L, _ L and € 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 Considering the value distribution of L, as shown in Figures 11 and 12, this paper takes the membership function as shown in Figure 13.
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 _ L and € L are shown in Figures 14 and 15: 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.
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.

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.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.
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       rapid pulling and pushing of the driving stick, the state response curve of the aircraft is shown in Figure 19.
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.
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.
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, _ L, and € 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.

Figure 1 .
Figure 1.The comparison of state curve.

Figure 3 .
Figure 3.The slice of the safe flight envelope.

Figure 4 .
Figure 4.The scatter of p T f (x, u Ã )50 on the envelope.

Figure 5 .
Figure 5.The state response curve of x 01 with the control u Ã .

Figure 6 .
Figure 6.The state response curve of x 02 when the control is equal to trim value.

Figure 7 .Figure 8 .
Figure 7.The state response curve of x 02 with the control u Ã .Figure8.The real-time value of control variables and f(x, 0).
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(Á)2D is opposite to that of u sup u(Á)2U

Figure 10 .
Figure 10.Comparison of the envelope of clean aircraft and iced aircraft.

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

Figure 11 .
Figure 11.The value distribution of L on the calculation grid.

Figure 12 .
Figure 12.The value distribution of L inside the safe set.Figure 15.The membership function for € L.

Figure 15 .
Figure 12.The value distribution of L inside the safe set.Figure 15.The membership function for € L.

Figure 14 .
Figure 14.The membership function for _ L.

Figure 13 .
Figure 13.The membership function for L.
are reasonable.The values of L, _ L, and € L vary in real-time with flight state and control input.According to the realtime values of L, _L, and € 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

Figure 16 .
Figure 16.The membership function for describing flight safety.

Figure 18 .
Figure 18.Safety assessment of the flight state space.

Figure 19 .
Figure 19.The state response curve.

Figure 17 .
Figure 17.The fuzzy inference rules for iced aircraft safety warning.

Figure 20 .
Figure 20.The output of safety warning.

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

Table 1 .
The inference rules.