Design of Ice Tolerance Flight Envelope Protection Control System for UAV Based on LSTM Neural Network for Detecting Icing Severity

Abstract

Icing on an unmanned aerial vehicle (UAV) can degrade aerodynamic performance, reduce flight capabilities, impair maneuverability and stability, and significantly impact flight safety. At present, most flight control methods for icing-affected aircraft adopt a conservative control strategy, in which small control inputs are used to keep the aircraft’s angle of attack and other state variables within a limited range. However, this approach restricts the flight performance of icing aircraft. To address this issue, this paper innovatively proposes a design method of an ice tolerance flight envelope protection control system for a UAV on the base of icing severity detection using a long short-term memory (LSTM) neural network. First, the icing severity is detected using an LSTM neural network without requiring control surface excitation. It relies solely on the aircraft’s historical flight data to detect the icing severity. Second, by modifying the fuzzy risk level boundaries of the icing aircraft flight parameters, a nonlinear mapping relationship is established between the tracking command risk level, the UAV flight control command magnitude, and the icing severity. This provides a safe range of tracking commands for guiding the aircraft out of the icing region. Finally, the ice tolerance flight envelope protection control law is developed, using a nonlinear dynamic inverse controller (NDIC) as the inner loop and a nonlinear model predictive controller (NMPC) as the outer loop. This approach ensures boundary protection for state variables such as the angle of attack and roll angle while simultaneously enhancing the robustness of the flight control system. The effectiveness and superiority of the method proposed in this paper are verified for the example aircraft through mathematical simulation.

1. Introduction

Aircraft icing results in a decrease in aerodynamic performance, characterized by a reduction in lift, an increase in drag, a decrease in the stall angle, and an increase in aerodynamic nonlinearity under the same flight conditions [1,2,3]. Changes in the aircraft’s aerodynamic characteristics cause deviations from the original flight control design point, leading to a decline in maneuverability and stability, which significantly affects flight safety.

In response to these issues, manned aircraft employ the following solutions. 1. Anti-icing systems are used, such as heating elements to melt the ice or mechanical means to remove it [4]. 2. The pilot monitors the aircraft’s icing severity through ice detection sensors [5] and adjusts the flight strategy to ensure that the aircraft remains within a safe operational range, reducing the probability of accidents.

Manned aircraft operating along fixed flight paths can predict the probability of icing occurrences in advance, thereby mitigating the associated risks by adjusting their routes accordingly. However, UAVs operate in expansive mission areas under complex meteorological conditions, and owing to constraints in cost and weight, they are unable to predict icing risk along flight paths or detect the icing severity. Furthermore, anti-icing systems are typically heavy, energy-intensive, and costly, which is why most UAVs are not equipped with such systems. In addition, UAVs typically operate in autonomous flight mode. Owing to the slow nature of the icing process, when flight data indicate abnormalities, the icing severity is often already significant, which leads to significant degradation in aerodynamic performance, thereby reducing the UAV’s ability to fulfil its mission requirements and potentially compromising flight safety. Therefore, it is essential to design an icing severity detection method for UAVs and to reconfigure the flight control law in real time on the basis of the detection results to ensure flight safety.

Icing detection is a critical step in ensuring flight safety. Existing methods for icing detection can be categorized into direct and indirect detection techniques [5]. Direct ice detection relies on sensors, such as capacitive, ultrasonic, and image-based detection methods [6,7]. Indirect ice detection is accomplished by detecting changes in the aircraft’s aerodynamic performance, primarily on the base of the aircraft’s motion and performance models, thereby enabling the detection of the icing severity. Reference [8] enables the detection of the icing severity by calculating the total energy derivative of the aircraft, thereby calculating the change in the drag coefficient. However, this method relies on an accurate performance model. Currently, the primary method for icing detection is observer-based recognition of the icing severity [9,10,11]. However, this method relies on active excitation and models that can accurately describe aircraft dynamics under nominal conditions. Icing causes a reduction in the aircraft’s stall angle of attack, thereby enhancing aerodynamic nonlinearity, and inappropriate excitation signals may induce instability in the aircraft’s motion, posing a risk to flight safety. Therefore, detecting the icing severity without the use of excitation signals is particularly important.

The calculation of safe flight control commands for icing aircraft is crucial for ensuring flight safety. References [12,13] design fuzzy risk boundaries for each flight parameter on the basis of the aircraft’s aerodynamic characteristics and, subsequently, obtain the time–domain response curves of the aircraft to input command signals through mathematical simulations. The risk assessment levels for different tracking commands are determined on the base of preset fuzzy risk boundaries. This method can assess the risk level associated with the aircraft’s command tracking and provides a certain degree of prediction for flight risk. However, the fuzzy risk boundaries it establishes are not adjusted for icing aircraft, leading to potential applicability issues. Reference [14] determined a range of safe state variables for icing aircraft through reachable set analysis. If the aircraft remains within this range, a control strategy exists that enables the aircraft to achieve the desired state. However, the method is computationally intensive and time-consuming, which makes it challenging to meet real-time computational requirements.

The design of control law for icing aircraft is a crucial step in ensuring flight safety. For icing aircraft, applying a large control input may cause the aircraft’s angle of attack and other state variables to fall within dangerously large or small ranges, thereby jeopardizing flight safety. In Reference [15], on the basis of the stall angle of attack, maximum roll angle, and minimum velocity under varying degrees of icing, the NDIC is employed to calculate the maximum available deflection angles of each control surface, which are then used as the limiting values for each control surface, thereby establishing the flight envelope protection control system. However, the effectiveness of the NDIC depends on the accuracy of aerodynamic modeling. If the aerodynamic model is not sufficiently accurate, the limit values of the control surfaces calculated at different times may vary significantly. Using the limit position during oscillations as the input for the protection control system can lead to oscillatory behavior or even divergence in the aircraft’s response. The model predictive control (MPC) approach has been widely used in recent years, and the MPC controller generates control commands for the current time by solving a constrained optimization problem for a future time horizon [16,17]. Reference [18] proposed an ice tolerance flight envelope protection control system based on variable-weight multiple-model predictive control. The effectiveness of MPC in flight envelope protection control applications is demonstrated. However, the method requires repetitive iterative design at predefined operating points, which entails significant computational effort.

This paper proposes a design method for UAV ice tolerance flight envelope protection control based on LSTM neural network icing severity detection to address the issues outlined above; an LSTM neural network is employed to detect the icing severity without the need for excitation signals on the basis of historical aircraft flight data. Second, the existing fuzzy risk boundaries for the aircraft’s state and maneuver quantities are corrected on the basis of the aerodynamic model of the icing UAV. The risk level topography of the UAV’s three-dimensional trajectory commands is then obtained through mathematical simulations, enabling the determination of the range of safe tracking commands for the UAV. Next, an ice tolerance flight envelope protection control law is designed on the basis of NMPC. The NDIC is applied to design the inner loop control law, whereas the NMPC is used to design the outer loop control law, thereby ensuring the protection of key state variables such as the angle of attack and roll angle. Finally, the effectiveness and superiority of the method are demonstrated through a mathematical simulation-based calculation.

The novelty of the research methodology presented in this paper can be summarized as follows:

• The icing severity detection method based on an LSTM neural network does not require active excitation and achieves rapid and accurate icing severity detection based on the UAV’s historical flight data while ensuring flight safety. Additionally, the offline-trained neural network has low computational requirements, making it suitable for UAV implementation.
• By correcting the fuzzy risk boundaries of UAV flight parameters in icing conditions, a three-dimensional risk level topography of trajectory commands is generated. This enables the determination of a range of safe tracking commands for the UAV.
• The integration of NMPC and NDIC simplifies the design process for icing-tolerant flight envelope protection. The control law effectively safeguards critical state variables, such as the angle of attack and roll angle, thereby ensuring flight safety for UAVs operating in icing conditions.

2. Flight Dynamics Modeling of the Icing UAV

2.1. General Parameters and Aerodynamic Modeling of the Icing UAV

In this paper, a specific type of UAV is used as a case study, with the three-view diagram and the overall parameters of the UAV presented in Figure 1 and

Table 1. General parameters of the UAV.

Parameters	Description	Value
m kg	mass	280
Sref (m2)	wing reference area	7
bref m	reference aircraft wingspan	9.8
cref m	reference aerodynamic chord length	0.7
Tmax KN	maximum thrust	2.5
Ixz (kg∙m2)	xz-axis product of inertia	2.65
Ix (kg∙m2)	x-axis moment of inertia	380
Iy (kg∙m2)	y-axis moment of inertia	150
Iz (kg∙m2)	z-axis moment of inertia	460

, respectively.

The aerodynamic force and moment coefficients of the UAV are calculated as shown in Equation (1) [19,20]:

$$\left\{\begin{array}{l}{C}_D=C_{D_0}+C_{D_{\alpha }}\cdot \Delta \alpha +C_{D_{\delta e}}\cdot \Delta \delta e\hfill \\ C_L=C_{L_0}+C_{L_{\alpha }}\cdot \Delta \alpha +C_{L_{\dot{\alpha }}}\cdot {\displaystyle \frac{\dot{\alpha }\overline{c}}{2V}}+C_{L_q}\cdot {\displaystyle \frac{q\overline{c}}{2V}}+C_{L_{\delta e}}\cdot \Delta \delta e\hfill \\ C_Y=C_{Y_{\beta }}\cdot \Delta \beta +C_{Y_p}\cdot {\displaystyle \frac{pb}{2V}}+C_{Y_r}\cdot {\displaystyle \frac{rb}{2V}}+C_{Y_{\delta r}}\cdot \Delta \delta r\hfill \\ C_m=C_{m_0}+C_{m_{\alpha }}\cdot \Delta \alpha +C_{m_{\dot{\alpha }}}\cdot {\displaystyle \frac{\dot{\alpha }\overline{c}}{2V}}+C_{m_q}\cdot {\displaystyle \frac{q\overline{c}}{2V}}+C_{m_{\delta e}}\cdot \Delta \delta e\hfill \\ C_l=C_{l_{\beta }}\cdot \Delta \beta +C_{l_p}\cdot {\displaystyle \frac{pb}{2V}}+C_{l_r}\cdot {\displaystyle \frac{rb}{2V}}+C_{l_{\delta a}}\cdot \Delta \delta a+C_{l_{\delta r}}\cdot \Delta \delta r\hfill \\ C_n=C_{n_{\beta }}\cdot \Delta \beta +C_{n_p}\cdot {\displaystyle \frac{pb}{2V}}+C_{n_r}\cdot {\displaystyle \frac{rb}{2V}}+C_{n_{\delta a}}\cdot \Delta \delta a+C_{n_{\delta r}}\cdot \Delta \delta r\hfill \end{array}\right.$$

(1)

where α and β are the heading and sideslip angles of the aircraft; p, q, and r are the roll, pitch, and yaw angular velocities, respectively; $\overline{c}$ is the average aerodynamic chord length of the aircraft; b is the wingspan of the aircraft; and V is the flight speed of the aircraft. δ_e, δ_a, and δ_r are the deflections of the elevator, aileron, and rudder, respectively. Equation (1) shows that the total aerodynamic force of an aircraft consists of three main components: (1) the basic aerodynamic force/moment coefficients C_i,basic such as C_L0, C_D0, and C_m0; (2) the control surface aerodynamic force/moment coefficients C_i,surface such as C_Lδe, C_mδe, and C_Yδr; and (3) the dynamic aerodynamic force/moment coefficients C_i,dynamic, such as C_mα, C_Yβ, C_lr, and C_lp.

Aircraft icing causes changes in aerodynamic derivatives. Bragg proposed a method for the aerodynamic modeling of icing with corrections to the aerodynamic derivatives on the basis of the clean configuration of the aircraft [21].

$$C_{A\mathrm{i}}=\left[1+\eta \left(t\right)\times K_{CA}\right]C_A$$

(2)

In Equation (2), C_A is the aerodynamic derivative value before icing and C_Ai is the aerodynamic derivative value after icing. The severity of aircraft icing, η(t), is a parameter that characterizes the severity of icing. It is related to the icing meteorological conditions and aircraft-related parameters and is a function of the icing time. K_CA is the coefficient icing factor (CIF), which reflects the extent of the impact of icing on each aerodynamic derivative of the aircraft. The value of the K_CA is closely related to the aircraft aerodynamic layout. By calculating K_CA, the aerodynamic correction method’s applicability to different aircrafts can be ensured. Additionally, as shown in Figure 1, the research object in this paper has a layout similar to that of a conventional aircraft, making this method effective in characterizing the aerodynamic properties of the icing UAV.

The above correction method effectively captures the main trends of different aerodynamic derivatives after icing, and its simplicity makes it commonly used in icing aircraft modeling. However, the constant value of K_CA in this correction method cannot accurately capture the changes in aerodynamic characteristics caused by icing at larger angles of attack, nor can it provide a reliable basis for designing protective control law for the flight envelope. Reference [22] combined wind tunnel experimental data and proposed a nonlinear coefficient icing factor (Variable CIF, VCIF) that varies with the angle of attack. As an example, the variation curve of the longitudinal aerodynamic coefficient icing factor for the UAV with respect to the angle of attack is shown in Figure 2.

Figure 2 shows that the variation in each coefficient icing factor is small when the angle of attack is in the range of −5° to 5°. However, when the angle of attack exceeds 5°, the coefficient icing factor becomes a nonlinear function of the angle of attack rather than remaining constant. Therefore, when analyzing the ice tolerance flight envelope protection problem at larger angles of attack, the constant coefficient icing factor K_CA in Equation (2) should be corrected to a nonlinear coefficient icing factor K_CA(α), which varies with the angle of attack. This correction more accurately captures the effect of icing on the aerodynamic characteristics of the aircraft at higher angles of attack. The improved aerodynamic modeling method for nonlinear icing is presented in Equation (3).

$$C_{A\mathrm{i}}=\left[1+\eta \left(t\right)\times K_{CA}(\alpha )\right]C_A$$

(3)

According to Equation (3), the aerodynamic coefficients for different icing severity and angles of attack can be calculated. Taking the longitudinal aerodynamic characteristics of an icing aircraft as an example, the curves of the longitudinal aerodynamic coefficients of the icing UAV, namely, C_L, C_D, and C_m, are shown in Figure 3.

As shown in Figure 3, icing causes only small aerodynamic changes in the low angle of attack region (0° to 5°), but leads to significant aerodynamic changes in the high angle of attack region (5° to 15°). This behavior more accurately reflects the phenomena of reduced lift, a rapid increase in drag, and an upwards pitching moment of the iced aircraft at larger angles of attack.

2.2. UAV Flight Dynamics Modeling

The nonlinear flight dynamics equations for the UAV can be expressed as Equations (4)–(7).

$$m\left[\begin{array}{c}\dot{u}+qw-rv\\ \dot{v}+ru-pw\\ \dot{w}+pv-qu\end{array}\right]=\left[\begin{array}{l}X\\ Y\\ Z\end{array}\right]+mg\left[\begin{array}{l}2(q_1{q}_3-q_0{q}_2)\\ 2(q_2{q}_3+q_0{q}_1)\\ q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right]$$

(4)

$$\left[\begin{array}{l}L\\ M\\ N\end{array}\right]=\left[\begin{array}{l}{I}_x\dot{p}+(I_z-I_y)qr-I_{zx}(pq+\dot{r})\\ I_y\dot{q}+(I_x-I_z)rp-I_{zx}(p^2-r^2)\\ I_z\dot{r}+(I_y-I_x)pq-I_{zx}(qr+\dot{p})\end{array}\right]$$

(5)

$$\left[\begin{array}{l}\dot{x}\\ \dot{y}\\ \dot{z}\end{array}\right]=\left[\begin{array}{ccc}1-2(q_2^2+q_3^2)& 2(q_1{q}_2-q_0{q}_3)& 2(q_1{q}_3+q_0{q}_2)\\ 2(q_1{q}_2+q_0{q}_3)& 1-2(q_1^2+q_3^2)& 2(q_2{q}_3-q_0{q}_1)\\ 2(q_1{q}_3-q_0{q}_3)& 2(q_2{q}_3+q_0{q}_1)& 1-2(q_1^2+q_2^2)\end{array}\right]\left[\begin{array}{l}u\\ v\\ w\end{array}\right]$$

(6)

$$\left[\begin{array}{l}{\dot{q}}_0\\ {\dot{q}}_1\\ {\dot{q}}_2\\ {\dot{q}}_3\end{array}\right]={\displaystyle \frac{1}{2}}\left[\begin{array}{cccc}0& -p& -q& -r\\ p& 0& r& -q\\ q& -r& 0& p\\ r& q& -p& 0\end{array}\right]\left[\begin{array}{l}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]$$

(7)

In Equations (4)–(7), m represents the mass of the UAV; g represents the gravitational acceleration; u, v, and w represent the axial, lateral, and normal velocity components in the body-fixed coordinate system, respectively; q₀, q₁, q₂, and q₃ represent the quaternions used to compute the attitude; X, Y, and Z represent the axial, lateral, and normal aerodynamic forces in the body axis, respectively; L, M, and N represent the roll, pitch, and yaw aerodynamic moments, respectively; I_x, I_y, I_z, and I_zx represent the moments of inertia and products of inertia; and x, y, and z represent the axial, lateral, and normal positions in the ground coordinate system, respectively.

The equation for the aircraft speed, angle of attack, and sideslip angle is shown in Equation (8).

$$\left\{\begin{array}{l}V=\sqrt{u^2+v^2+w^2}\\ \alpha =\arctan (w/u)\\ \beta =\arcsin (v/V)\end{array}\right.$$

(8)

The relationship between the quaternion and the attitude angles of the vehicle is shown in Equation (9).

$$\left[\begin{array}{l}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]=\left[\begin{array}{l}\cos {\displaystyle \frac{\varphi }{2}}\cos {\displaystyle \frac{\theta }{2}}\cos {\displaystyle \frac{\psi }{2}}+\sin {\displaystyle \frac{\varphi }{2}}\sin {\displaystyle \frac{\theta }{2}}\sin {\displaystyle \frac{\psi }{2}}\\ \sin {\displaystyle \frac{\varphi }{2}}\cos {\displaystyle \frac{\theta }{2}}\cos {\displaystyle \frac{\psi }{2}}-\cos {\displaystyle \frac{\varphi }{2}}\sin {\displaystyle \frac{\theta }{2}}\sin {\displaystyle \frac{\psi }{2}}\\ \cos {\displaystyle \frac{\varphi }{2}}\sin {\displaystyle \frac{\theta }{2}}\cos {\displaystyle \frac{\psi }{2}}+\sin {\displaystyle \frac{\varphi }{2}}\cos {\displaystyle \frac{\theta }{2}}\sin {\displaystyle \frac{\psi }{2}}\\ \cos {\displaystyle \frac{\varphi }{2}}\cos {\displaystyle \frac{\theta }{2}}\sin {\displaystyle \frac{\psi }{2}}-\sin {\displaystyle \frac{\varphi }{2}}\sin {\displaystyle \frac{\theta }{2}}\cos {\displaystyle \frac{\psi }{2}}\end{array}\right]$$

(9)

In Equation (9), φ, θ, and ψ are the roll, pitch, and yaw angles in the body–axis system, respectively.

The UAV control surface model adopts a first-order inertia system, as shown in Equation (10).

$${\delta }_i={\delta }_{i,cmd}{\displaystyle \frac{1}{T_{i}s+1}}$$

(10)

In Equation (10), I = a, e, and r represent the aileron, elevator, and rudder, respectively. T is the control surface response time, and the subscript “cmd” denotes the control surface command. The specific parameters of the control surface model are presented in

Table 2. UAV control surface parameters.

Control Surface	δe/°	δa/°	δr/°
Range	[−30, 30]	[−30, 30]	[−45, 45]
Response Time	30	30	30
Response Rate Limit	[−60, 60]

.

The mathematical simulation process for the icing UAV is as follows. First, the aerodynamic coefficients C_Lα, C_mα, etc., on the right side of Equation (1) under the current icing severity are calculated using the aerodynamic correction model in Equation (3). Next, the aerodynamic forces and moments acting on the UAV are determined based on its flight state, including the angle of attack, velocity, and other state variables. From this, the X, Y, and Z, and moments L, M, and N in the body–axis system, are calculated. Finally, the UAV’s state response at the current moment is calculated using the flight dynamics in Equations (4)–(7).

3. A Method for Detecting Icing Severity Based on an LSTM Neural Network

3.1. Principles of Neural-Network-Based Detection Methods

Traditional methods for detecting aircraft icing severity rely on identifying aerodynamic parameters through the active excitation of control surfaces and the formulation of identification equations. However, these approaches face significant challenges: the design of appropriate excitation signals and identification equations is complex, and improper signal selection can lead to substantial identification errors while posing risks to flight safety. Additionally, solving the identification equations often requires iterative optimization, imposing high computational demands that may exceed the UAV’s onboard capabilities.

In response to these limitations, this paper presents an innovative icing severity detection method using an LSTM neural network. This method eliminates the need for active excitation and achieves rapid and accurate icing severity detection by extracting information from historical UAV flight data. Furthermore, the use of an offline-trained neural network reduces computational requirements, making the method highly suitable for real-time application in UAV operations.

Aircraft icing is a nonlinear time-varying process. As shown in Figure 3, the aerodynamic data of the icing UAV exhibit a more significant deviation than those of the clean configuration. For the same flight conditions and identical control surface deflections, the state variable responses of the UAV differ because of variations in the aerodynamic characteristics of the aircraft. On the basis of the above analysis, the differences in the UAV’s state variable responses can be captured by neural networks, enabling the detection of the icing severity.

When detecting the icing severity of the UAV, the state and control input parameters from the previous period are used as inputs to the neural network, with the output being the current icing severity of the UAV. The input and output parameters selected for the neural network are presented in Equation (11):

$$\left\{\begin{array}{l}S=[\begin{array}{ccccccccc}V& \varphi & \theta & \alpha & \beta & {\delta }_e& {\delta }_r& {\delta }_a& {\delta }_p\end{array}]\\ IN=[S^{t-n}{S}^{t-n+1}\begin{array}{c}\dots \end{array}{S}^{t-1}{S}^t]\\ OUT=[{\eta }^t]\end{array}\right.$$

(11)

where S denotes the input state parameters selected by the neural network, including the flight velocity V; the attitude angles ϕ and θ; the airflow angles α and β; and the maneuvering variables δ_e, δ_r, δ_a, and δ_p. IN denotes the neural network input, which consists of the state parameters S from t−n to the moment t. OUT denotes the neural network output, which is the icing severity η at moment t.

3.2. Construction of the Neural Network Structure and Parameter Setting

The LSTM neural network [23] is a specialized type of recurrent neural network (RNN) that addresses the issue of long-term dependencies by remembering information over extended periods. The LSTM unit, illustrated in Figure 4, controls the flow of information adaptively through the use of forget gates, input gates, and output gates, allowing it to capture long-term dependencies in the data. As shown in Equation (11), the neural network input data are inherently temporal. To effectively capture the time-dependent information within the flight data, an LSTM layer is employed in the network architecture.

The prediction network constructed in reference [24] achieves a balance between prediction accuracy and computational complexity, and a similar neural network structure is adopted in this paper, as shown in Figure 5. The selection process for the network’s layers and neurons is as follows. 1. The initial number of layers and neurons per layer is determined based on the complexity of the problem. 2. The network structure is refined by evaluating factors such as training time, number of iterations, detection accuracy, and generalization performance. Through iterative training and optimization, the final configuration of the network’s layers and neurons are established for the specific requirements of this study.

This approach ensures an optimal trade-off between computational efficiency and detection performance.

As shown in Figure 5, the neural network architecture consists of an input layer, full connected layers, activation layers, LSTM layers, dropout layers, and an output layer. Among them, the fully connected layer performs feature extraction and fusion by adjusting the neuron weights; the activation layer uses the rectified linear unit (ReLU) function to introduce nonlinearities, enabling the network to handle complex relationships; the LSTM layer captures temporal dependencies in UAV flight data; and the loss layer helps reduce overfitting while simultaneously improving the neural network’s generalization ability. The parameters of the neural network structure are listed in

Table 3. Neural network structure parameters.

Name Numerical Value	Name Numerical Value
Number of Layers	14
Input Layer’s Nodes	9
LSTM Layer’s Number	3
LSTM Layer’s Nodes	300
Full Connect’s Number	4
Full Connect’s Nodes	200
Output Layer’s Nodes	1

.

As shown in Table 3, the neural network consists of nine input nodes and one output node, with a total of 14 layers. Among these, there are three LSTM layers, each containing 300 neurons. The fully connected layers consist of four layers, with the number of neurons in the final fully connected layer set to one to connect to the output node. The other fully connected layers each contain 200 neurons.

The parameters for the neural network training process are shown in

Table 4. Neural network training parameters.

Name Numerical Value	Name Numerical Value
Max Epochs	500
Mini BatchSize	128
Gradiet Threshold	0.5
Initial LearnRate	1 × 10−3
LearnRate DropPeriod	50
LearnRate DropFactor	0.8

.

As shown in Table 4, the neural network is trained for 500 iterations, with 128 datasets used in each iteration. The initial learning rate is set to 1 × 10⁻³, and every 50 iterations, the learning rate is reduced by 20% to prevent convergence to a local optimum.

3.3. Data Preparation for Neural Network Training and Testing in Icing Severity Detection

On the basis of the icing UAV flight dynamics model established in Section 2, the corresponding flight control law is formulated and used to generate data for neural network training and testing. The PID flight control law is designed for the example UAV, and the structure of the longitudinal control law is shown in Figure 6. The longitudinal control law uses a pitch angle configuration for stabilizing flight, which is divided into inner and outer loops. The inner loop acts as a stabilization link, providing feedback on the UAV’s pitch angle velocity and angle of attack to improve the longitudinal modal characteristics of the UAV. The outer loop uses the pitch angle as a control command.

The structure of the lateral-directional control law is shown in Figure 7. The lateral heading control law consists of both roll axis and yaw axis channels. The inner loop of the roll axis channel acts as a damper, feeding back the UAV’s roll angular velocity to improve the roll axis damping characteristics, whereas the outer loop uses the roll angle as a control command. In the yaw axis channel, the inner loop functions as a stabilization link, feeding back the UAV’s yaw angular velocity and sideslip angle to enhance the heading modal characteristics. The outer loop, in this case, uses the sideslip angle as a control command.

The closed-loop characteristics of the UAV can be adjusted by modifying the parameters of the inner and outer loops in the control law structure in Figure 6 and Figure 7, thereby generating flight simulation data corresponding to different closed-loop characteristics. These data can then be used to train the neural network, improving its generalization performance across various flight control laws.

The data used for neural network training were generated through mathematical simulations. Before each simulation, the UAV’s flight altitude, speed, icing severity, and tracking commands were randomly assigned, with the variation ranges of the corresponding parameters shown in Equation (12). The flight control law parameters designed in Figure 6 and Figure 7 are also subjected to a random 10% variation to ensure that the neural network can capture the changes in the UAV’s aerodynamic characteristics caused by icing.

$$\left\{\begin{array}{l}V\in [110\mathrm{m}/\mathrm{s},150\mathrm{m}/\mathrm{s}]\\ H\in [6000\mathrm{m},7000\mathrm{m}]\\ V_{cmd}\in [-15\mathrm{m}/\mathrm{s},20\mathrm{m}/\mathrm{s}]\\ {\theta }_{cmd}\in [-{15}^{\circ },{15}^{\circ }]\\ {\varphi }_{cmd}\in [-{20}^{\circ },{20}^{\circ }]\\ {\beta }_{cmd}=0^{\circ }\\ \eta \in [0,1]\end{array}\right.$$

(12)

To ensure the safety of icing UAV flights and to quickly reconstruct the flight control law, it is essential to rapidly detect the icing severity on the basis of the flight status. As indicated in Reference [9], the aircraft icing severity curve changes relatively slowly. Therefore, in the mathematical simulation calculation, the UAV’s icing severity is assumed to increase by 0.002 per second, starting from the initial icing severity as given in Equation (12).

A total of 5000 sets of mathematical simulations were conducted, and the simulation data were segmented into fixed lengths for training and validation according to Equation (11). In total, 70% of the data were used as the training set for model learning and parameter tuning, whereas the remaining 30% were used as the testing set to evaluate the model’s performance and generalization ability.

3.4. Neural Network Training Results and Effectiveness Evaluation

On the basis of the flight data obtained from the simulation, the neural network was trained, and the root mean square error (RMSE) and mean absolute error (MAE) were used as evaluation metrics during the training process. The Equations for calculating the RMSE and MAE are shown in Equations (13) and (14).

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{({\eta }_i-{\eta }_{true})}^2}}$$

(13)

$$MAE={\displaystyle \frac{{\displaystyle \sum _{i=1}^n|{\eta }_i-{\eta }_{true}|}}{n}}$$

(14)

In Equations (13) and (14), η_i represents the icing severity predicted by the neural network for each training dataset, η_true represents the actual icing severity for that dataset, and n represents the total number of training sets. The variation in the evaluation parameters during the neural network training process is shown in Figure 8.

As shown in Figure 8, the neural network exhibits clear convergence, and after 500 iterations, it reaches convergence with the MAE below 0.002 and the RMSE below 0.02, indicating satisfactory training results.

The trained neural network is used to detect the icing severity in the test set, which contains a total of 1500 datasets. The mean error and standard error of the icing severity detection are calculated for detection times of 3 s, 5 s, and 7 s, with the results presented in

Table 5. Neural network detection errors at various detection times.

Detection Time	Mean Error	Standard Error
3 s	0.0119	0.0022
5 s	0.0105	0.0012
7 s	0.0087	0.00087

.

As shown in Table 5, the mean error and standard deviation of the neural network detection gradually decrease as the detection time increases. At 7 s, the mean error is less than 0.01, and the standard deviation is less than 0.0001, indicating high accuracy of the detection results.

Moreover, the detection error rate for each dataset in the test set at different detection times is calculated using Equation (15).

$${\eta }_{error}(\%)={\displaystyle \frac{|{\eta }_{lstm}-{\eta }_{true}|}{{\eta }_{true}}}$$

(15)

In Equation (15), η_error represents the icing severity detection error rate, η_lstm is the icing severity detected by the neural network, and η_true is the true icing severity of the dataset. The detection error rate was categorized into four intervals: within 2%, between 2% and 5%, between 5% and 10%, and greater than 10%. The distribution of the detection error rates is shown in Figure 9.

As shown in Figure 9, the icing severity detection error rate gradually decreases with increasing time. The detection value gradually converges to the vicinity of the set value, and the data within the range of detection error rate <2% increase, whereas the data in other intervals gradually decrease. When the detection time is 7 s, the detection error rate for more than 90% of the data is less than 2%. This shows that the accuracy of the icing severity detection method based on the LSTM neural network is high.

To evaluate the generalization performance of the neural network, it is tested using simulation data that differ from the training dataset. The mathematical simulation parameters are shown in Equation (16), with expanded ranges for the initial altitude, speed, tracking commands, and icing severity compared with those in Equation (12). The flight control law parameters shown in Figure 6 and Figure 7 were also subjected to a 20% parameter perturbation to evaluate the sensitivity of the neural network to variations in aerodynamic characteristics. For each icing severity, 100 mathematical simulations with random tracking commands and initial states were conducted to obtain validation data.

$$\left\{\begin{array}{l}V\in [90\mathrm{m}/\mathrm{s},160\mathrm{m}/\mathrm{s}]\\ H\in [5500\mathrm{m},7500\mathrm{m}]\\ V_{cmd}\in [-20\mathrm{m}/\mathrm{s},30\mathrm{m}/\mathrm{s}]\\ {\theta }_{cmd}\in [-{20}^{\circ },{20}^{\circ }]\\ {\varphi }_{cmd}\in [-{30}^{\circ },{30}^{\circ }]\\ {\beta }_{cmd}=0^{\circ }\\ \eta \in \{0.1,0.4,0.7,1,1.3\}\end{array}\right.$$

(16)

The trained neural network is applied to detect the icing severity on the newly acquired data, and the results of the neural network’s detection error for the icing severity are shown in Figure 10 and

Table 6. Statistical parameters of detection errors at different icing severities.

η	Mean Error	Standard Error
0.1	−0.0033	0.0095
0.4	−0.0012	0.0094
0.7	−0.0001	0.0074
1	−0.0015	0.0077
1.3	0.0008	0.0061

.

Figure 10 and Table 6 show that, for the newly acquired data, the mean error in icing severity detection remains below 0.002, whereas the standard deviation does not exceed 0.01. The neural network effectively detects the icing severity on the basis of historical flight data, even when the initial speed, altitude, and tracking commands are expanded beyond their original ranges. Moreover, the neural network accurately detected the icing severity for previously untouched flight data with an icing severity of 1.3, demonstrating its ability to extract aerodynamic degradation caused by icing from aircraft data. This demonstrates that the icing severity detection method based on the LSTM neural network achieves higher accuracy and exhibits superior generalization performance.

4. Fuzzy-Risk-Boundary-Based Method for Calculating the Risk Level for Commands

4.1. Methodology for Calculating the Risk Level of Commands

After UAV icing, inappropriate input tracking commands can generate large maneuver volumes, which, in turn, push the flight state into higher-risk intervals, jeopardizing flight safety. Therefore, it is necessary to calculate the risk level of different tracking commands on the basis of the changes in the aerodynamic characteristics of the UAV after icing, from which the command value with a lower risk value is selected as the tracking command for the mission.

The steps of the method are as follows. 1. The risk boundaries of the flight state variables and maneuvering variables are determined on the basis of the aerodynamic characteristics of the UAV. 2. The tracking command signal is selected according to the flight mission, and its calculation range is determined. 3. Through mathematical simulation, the UAV time–domain response curve is obtained, the risk value of each flight parameter is weighted to calculate the risk level of different tracking commands, and a three-dimensional risk level topology map for the commands is generated. This method effectively evaluates the risk level of UAV tracking commands and offers predictive insights into flight safety.

In the flight process, the risk boundary values for each state and control variable are not fixed. Therefore, the risk state of each flight parameter is classified fuzzily in the form of risk intervals. Moreover, the flight risk level of a UAV is jointly determined by its various flight parameters, and the UAV’s flight is considered safe when all the flight parameters are within the lower range of risk values. When any of the flight parameters deviate from the safe range, the flight safety of the UAV is at risk. Therefore, important flight parameters were selected as assessment indicators, and the risk levels of these indicators were evaluated to determine the overall flight risk level of the UAV.

The evaluation metrics selected in this paper include the airspeed V, angle of attack α, sideslip angle β, roll angle φ, attitude angle θ, vertical velocity H_dot, normal overload n_z, elevator deflection δ_e, aileron deflection δ_a, and rudder deflection δ_r. The risk boundaries of each state variable for the example UAV are classified on the basis of the fuzzy risk boundary setting method from Reference [13], as shown in

Table 7. Fuzzy risk boundaries for flight parameters of clean-configuration UAV (η = 0).

Parameter	Unit	Break Point
		e		c		a		b		d		f
V	m/s	45		50		55		125		140		154
α	deg	−10		−2.5		0		8		14		18
β	deg	−30		−15		−10		10		15		30
ϕ	deg	−35		−25		−15		15		25		35
θ	deg	−40		−10		−5		16		25		55
Hdot	m/s	−25		−7		−4.5		15		25		40
nz	-	−2.5		0.25		0.5		1.5		1.75		3.75
δe	deg	−29.9		−23		−19		9		13		19.9
δa	deg	−29.9		−23		−18		18		23		29.9
δr	deg	−44.9		−38		−30		30		38		44.9

.

The color stripes in Table 7 represent the risk level of the flight parameters, with green, yellow, red, and black from the center to the edges indicating a progressively increasing risk level. The weight design within the range of different risk class intervals is shown in Equation (17):

$$R_i\left(x_i,a,b,c,d,e,f\right)=\left\{\begin{array}{l}1\begin{array}{c} a<x_i\le b\end{array}\\ 2\begin{array}{c} c<x_i\le a,b<x_i\le d\end{array}\\ 3\begin{array}{c} e<x_i\le c,d<x_i\le f\end{array}\\ 4\begin{array}{c} x_i\le e,x_i>f\end{array}\end{array}\right.$$

(17)

where x_i represents different flight parameters, R_i (i = 1, 2,…, 10) represents the risk level of each flight parameter, and a, b,…, f represent the corresponding boundary values of the different risk level intervals, corresponding to the six boundary values for each flight parameter in Table 7. If the combined risk level change during UAV tracking commands is characterized by R, then R is calculated as shown in Equation (18):

$$R={\displaystyle \sum _{i=1}^{10}{P}_i{R}_i}$$

(18)

where P_i is the weight representing the degree of contribution of the risk level of different flight parameters to the total risk level, which is calculated via the entropy weight and grey correlation algorithm proposed in Reference [13], which is more objective and adaptable.

4.2. Fuzzy Risk Boundary Correction Method

The fuzzy risk boundaries designed for clean configurations may not be applicable because of the significant changes in the aerodynamic characteristics of UAV after icing. Therefore, this paper proposes an innovative method to correct fuzzy risk boundaries on the basis of the UAV’s icing severity and flight dynamics model, enabling the more accurate calculation of the risk level of tracking commands.

As mentioned in Section 2, the value of the aerodynamic derivative in Equation (1) changes after icing. Therefore, to determine the fuzzy risk level boundaries for different flight parameters under UAV icing conditions, the overall force and moment characteristics of the UAV need to be analyzed as a benchmark. Taking the angle of attack risk boundary correction calculation for an icing UAV as an example, the aerodynamic coefficients involved in the angle of attack in Equation (1) include the lift coefficient C_L, the drag coefficient C_D, and the pitching moment coefficient C_m. The method for determining the boundary limit value of the angle of attack risk from these three perspectives is outlined below.

The lift, drag, and pitching moment experienced by the example UAV in a trimmed condition are taken as a baseline. The angle of attack at which the lift, drag, and pitching moment after icing do not deviate beyond a specified margin from the baseline value is selected as the limiting value for the angle of attack boundary protection. The margins are selected on the basis of the boundary limits for different risk classes. Selecting the angle of attack limit values in this manner ensures that the iced UAV experiences an acceptable amount of lift loss, drag increase, and pitching moment variation within the angle of attack limit region. With this method, the engine thrust required at the same level of flight speed remains essentially unchanged, and the static stability can be adjusted to be nearly the same under both icing and clean conditions by modifying the control law parameters. After the three limiting values for the angle of attack are determined, the smallest limiting value is selected through a comparison, which represents the risk boundary of the angle of attack that ultimately satisfies the three design requirements for the lift, drag, and pitching moment. Similarly, risk boundaries for other flight state parameters can be determined via this approach.

The risk boundaries for control surface deflections, such as elevator, rudder, and aileron, are determined in a similar manner. On the basis of the force and moment characteristics related to these control inputs in Equation (1), appropriate safety margins are defined relative to the UAV’s trimmed baseline values. These margins are then used to establish the risk boundaries for the maneuvering variables. Unlike flight state quantities, maneuvering variables are inherently constrained by the aircraft’s design limitations. Therefore, additional considerations must be accounted for when determining the risk boundary values to ensure that they meet all the design requirements.

Using the fuzzy risk boundary correction method based on the UAV flight dynamics model described above, the fuzzy risk boundaries for the flight parameters of the UAV under an icing severity of 1 were determined, as shown in

Table 8. Fuzzy risk boundaries for flight parameters of the icing UAV (η = 1).

Parameter	Unit	Break Point
		e		c		a		b		d		f
V	m/s	47.5		51		58		92		99.45		106.35
α	deg	−2.1		1.4		2.6		6.4		9.3		11.2
β	deg	−21.3		−7.1		−4.7		4.7		7.1		21.3
ϕ	deg	−24.3		−14.9		−9.75		9.75		14.9		24.35
θ	deg	−32.0		−7.3		−3.2		14.0		21.4		46.1
Hdot	m/s	−25		−7		−4.5		15		25		40
nz	-	−2.46		0.23		0.51		1.49		1.74		3.72
δe	deg	−14.2		−10.9		−9.0		4.3		6.2		9.4
δa	deg	−11.7		−9.0		−7.0		7.0		9.0		11.7
δr	deg	−25.4		−21.4		−17.0		17.0		21.4		25.4

.

A comparison of Table 7 and Table 8 reveals that, after the aircraft experiences icing, the risk intervals for various flight parameters become narrower, indicating an increased level of flight risk.

4.3. Calculation of the Risk Level Topology Map for 3D Commands

On the basis of the data in Table 7 and Table 8, the risk levels corresponding to different tracking commands and icing severities can be determined. This allows for identifying the range of safe commands necessary to accomplish the flight mission while ensuring flight safety.

Reference [25] revealed that aircraft icing has specific environmental requirements, meaning that parameters such as the liquid water content (LWC) and mean volume drop diameter (MVD) need to be within a certain range, so aircraft icing occurs only within a specific range. Therefore, in this paper, the UAV tracking command is selected as the trajectory command to ensure that the UAV follows a stable route, allowing it to exit the icing area, reduce the icing severity, and thereby lower the flight risk. On the basis of the above analyses, V, γ, and χ are selected as the flight control commands of the UAV for calculating the risk level topology map. The tracking commands and icing severity are defined as shown in Equation (19).

$$\left\{\begin{array}{l}{V}_{cmd}\in [70\mathrm{m}/\mathrm{s},120\mathrm{m}/\mathrm{s}]\\ {\gamma }_{cmd}\in [-{15}^{\circ },{20}^{\circ }]\\ {\chi }_{cmd}\in [-{65}^{\circ },{65}^{\circ }]\\ \eta =\begin{array}{ccc}0& or& 1\end{array}\end{array}\right.$$

(19)

The time–domain response curves for each flight parameter, corresponding to different commands and icing severities, are obtained through mathematical simulation. Firstly, the fuzzy risk boundaries applicable to the UAV are determined on the basis of the icing severity. The risk level of the flight parameters is then determined on the basis of their values at different moments, which, in turn, allows for the calculation of the risk weight for each flight parameter via Equation (17). Afterwards, the risk level weights of each flight parameter are considered comprehensively, and the overall risk level of the UAV at different moments is determined via Equation (18). Finally, the combined risk level at different moments is weighted to determine the risk level of the tracking command. The results of the risk level calculations for different tracking commands and icing severities are presented in Figure 11. Figure 11a shows the topology map of the 3D command risk level for the clean configuration, whereas Figure 11b illustrates the topology map of the 3D command risk level for an icing severity of 1. The risk level is categorized into four levels: green, yellow, red, and black. As the color intensifies, the risk level increases accordingly.

Figure 11a shows that, for the same trajectory angle command, the command risk level decreases as the speed tracking command increases. The reason for this is that, as the speed increases, the required change in the angle of attack to track the same trajectory angle command decreases. Consequently, each flight parameter remains within the low-risk interval, leading to a lower risk level for the tracking command. As the velocity tracking command decreases, to track the same trajectory angle command, the UAV requires a larger control surface deflection to maintain flight state quantities, such as the angle of attack, within the high-risk interval to generate the aerodynamic forces required to complete the mission, resulting in a higher risk level for the tracking command. When tracking the same speed command, as the tracking trajectory angle command decreases, the risk level of the tracking command also decreases. The reason for this is that, as the tracking trajectory angle command decreases, the increment of lift required by the UAV decreases. With the same dynamic pressure conditions, the increment of the angle of attack required also decreases, leading to a reduction in the variation of control surface deflection, which, in turn, lowers the risk value of the flight parameters. The tracking command risk level increases as the tracking trajectory angle command is further reduced to a larger negative value. The reason for this is that, to track a large negative trajectory angle command, the UAV will need to control the elevator to produce a significant positive deviation, causing the UAV’s angle of attack to continuously decrease, potentially reaching a negative value and entering the left end of the high-risk interval in Table 7, which leads to an increase in the overall risk level of the UAV.

A comparison of Figure 11a,b reveals that the black range in Figure 11b is significantly larger than that in Figure 11a, indicating that icing on the UAV leads to a substantial increase in the risk value of the tracking command. The effect of UAV icing on the available trajectory angle is more pronounced, with greater contraction on the upper boundary and less contraction on the lower boundary. This is primarily due to the deterioration of the UAV’s aerodynamic characteristics caused by icing, which leads to increased drag and reduced lift under the same flight conditions. The upper boundary trajectory angle is positive, and the maximum thrust can no longer overcome the additional drag induced by icing during climb. In contrast, the lower boundary trajectory angle is negative, allowing the potential energy from gliding to offset some of the increased drag.

On the basis of the above analysis, the risk level of the corresponding command can be calculated considering both the icing severity and the value of the tracking command. This approach provides the UAV with a safe range of tracking commands, ensuring that the UAV completes its flight tasks safely.

5. NMPC-Based Ice Tolerance Flight Envelope Protection Control Law Design

Under icing conditions, the aerodynamic characteristics of the UAV deteriorate, leading to a subsequent decline in flight performance. The control law design should ensure the fast, accurate, and stable tracking of command signals while prioritizing flight safety. In this paper, an NMPC-based ice tolerance flight envelope protection control law is designed following the aforementioned principles.

5.1. Design of the Inner Loop Structure of the Control Law Based on NDIC

NDIC is a superior control method for achieving the decoupled control of each state variable without the need for complex gain adjustments while providing clear physical insights into the system’s behavior [26]. Therefore, NDIC is employed to design the inner loop control law for the ice tolerance flight envelope protection of the UAV. The parameters of the inner loop control law are adjusted on the basis of the icing severity detection result discussed in Section 3, thereby optimizing the performance of the controller. The basic structure of the inner loop flight control law is shown in Figure 12.

As shown in Figure 12, the inner loop structure of the ice tolerance flight envelope protection control law consists of two loops. The inner loop control variables are p, q, and r, whereas the middle loop control variables are α, β, and μ, with μ representing the roll angle around the velocity vector axis. Ideally, the relationship between the control command and the actual response of the variable follows a first-order inertial link. The relationship between the command value and the actual response is expressed in Equation (20):

$$X={\displaystyle \frac{{\omega }_x}{s+{\omega }_x}}{X}_c$$

(20)

where X_c is the command for the airflow angle or attitude angular velocity and X is the actual response of the airflow angle or attitude angular velocity. ω_x represents the tracking response dynamic characteristic band, which includes ω_α, ω_β, ω_μ, ω_p, ω_q, and ω_r.

5.2. Design of the Outer Loop Control Law Based on NMPC

Aircraft icing is a nonlinear and time-varying process. As the icing severity increases, the nonlinearity of UAV aerodynamic performance becomes more pronounced, and the safety margins for each flight parameter progressively diminish. The state quantities that are related mainly to flight safety during flight are the angle of attack, roll angle, and sideslip angle. Therefore, the protection of the above variables is particularly important.

The analysis of Figure 3 reveals that the aerodynamic characteristics of the icing UAV exhibit minimal differences within the range of small angles of attack. In this range, icing has a limited impact on flight performance. However, beyond a certain angle of attack threshold, the aerodynamic characteristics change significantly. Therefore, the angle of attack for the icing UAV should not be too large during command tracking.

The NMPC generates a control sequence for the current moment by solving a constrained optimization problem over a set of future time steps. While ensuring the effectiveness of command tracking, it also ensures that the system parameters remain within the specified limit boundaries. Therefore, the NMPC offers significant advantages in terms of flight envelope protection.

On the basis of the above analysis, this paper proposes a flight envelope protection control method based on NMPC. The NMPC controller is applied to design the outermost loop of the UAV control system. The nonlinear equation of state for the outer loop of the UAV is shown in Equation (21) [27], where γ and χ are the flight path angle and flight path heading angle of the UAV and T, L, and D are the thrust, lift, and drag of the UAV, respectively. The state and output quantities are V, γ, and χ, while the control variables are T, α, and μ. L and D are calculated in real time by solving the flight dynamics Equation (4) using the current flight state parameters. This approach is independent of the aerodynamic model of the icing UAV, thereby effectively ensuring the accuracy of the NMPC.

$$\left\{\begin{array}{l}\dot{V}={\displaystyle \frac{1}{m}}(T\cos \alpha -D-Mg\sin \gamma )\\ \dot{\gamma }={\displaystyle \frac{1}{mV}}(T\sin \alpha \cos \mu +L\cos \mu -Mg\cos \gamma )\\ \dot{\chi }={\displaystyle \frac{1}{mV\cos \gamma }}(T\sin \alpha \sin \mu +L\sin \mu )\end{array}\right.$$

(21)

The structure of the UAV outer loop NMPC controller is constructed according to Equation (21), as shown in Figure 13.

As shown in Figure 13, the UAV determines the time-varying limiting boundaries for the angle of attack and velocity roll angle on the basis of the icing severity detection values. The optimization objective function is constructed by considering the input tracking command values V_cmd, γ_cmd, and χ_cmd, the designed output quantities, and the control variable weights. The UAV control commands δ_p, α, and μ in the predicted time domain are obtained via optimization solving on the basis of the UAV outer loop mathematical model. These control commands are then used as the inner loop inputs for the UAV within the control time domain range. The NMPC controller continuously repeats the above process, using rolling optimization to obtain the UAV inner loop control commands. This approach ensures effective command tracking while simultaneously protecting the key state variables.

Figure 12 and Figure 13 are combined to establish the overall structure of the UAV ice tolerance flight envelope protection control law, as shown in Figure 14.

The parameters to be set during the design process of the NMPC controller include the following: error weights, control weights, control rate of change weights, prediction time horizon, control time horizon, and control boundaries. The error weights should be designed to ensure that the UAV can track the V, γ, and χ commands both quickly and accurately [28,29,30]. The control weights are applied to the UAV’s δ_p, α, and μ commands to prevent the inner loop control commands from becoming too large or too small. The control command changing rate weights regulate the rate of change in the inner loop commands to prevent excessively rapid changes, which could cause instability and lead to the dispersion of the UAV’s motion. The design of the prediction and control time domains involves trade-offs between system stability, dynamic performance, predictability, and computational complexity [30,31,32].

According to Section 3, Section 4 and Section 5, an ice tolerance flight envelope protection control system for the UAV can be established, as illustrated in Figure 15.

As shown in Figure 15, the working principle of the ice tolerance flight envelope protection control system is as follows: first, the icing severity is detected by the neural network on the basis of the historical flight data of the UAV. Then, the risk levels of different trajectory tracking commands for that icing level are calculated, and the safe tracking commands from them are selected as the UAV tracking command inputs. Finally, the UAV control law is switched to the NMPC flight envelope protection control law, and the parameters of the NMPC flight envelope protection control law are adjusted according to the results of the icing severity detection, completing the control law reconstruction and ensuring flight safety.

6. Integrated Mathematical Flight Simulation Verification

6.1. Design of Control-Law-Related Parameters

NDIC is selected for comparison with the ice tolerance flight envelope protection control system based on the LSTM neural network designed in this paper. The outer loop of the NDIC adopts the trajectory command configuration, whereas the inner loop structure remains the same as that shown in Figure 12. The outer loop control variables are V, γ, and χ. The relationship between the control command and the actual response of the variables remains a first-order inertial link, where ω_V, ω_γ, and ω_χ represent the dynamic characteristic bands of the tracking response for the outer loop commands. The structure of the NDIC is shown in Figure 16, with the corresponding reference model parameters listed in

Table 9. Trajectory command configuration nonlinear dynamic inverse control law reference model parameters.

Parameter	ωV	ωγ	ωχ	ωα	ωβ	ωμ	ωp	ωq	ωr
Value (rad/s)	0.1	0.1	0.1	1	1	1	5	5	5

.

The NMPC controller design parameters are shown in

Table 10. NMPC controller design parameter results.

Design Parameters	Symbol	Design Value
Error weights	Q	diag{8,15,4}
Control weights	R	diag{1,5,2}
Control rate of change weights	S	diag{1,5,1}
Prediction time domain	p	4 s
Control time domain	m	1 s
Control boundary	U	function

, where the error weight variables are ordered as V, γ, and χ, and the control-related weights are ordered as δ_p, α, and μ. The limiting boundaries of the control variables α and μ are adjusted in real time on the basis of the icing severity detection values.

During the flight of an icing UAV, the aerodynamic characteristics are influenced primarily by the angle of attack and sideslip angle. In the design of the icing UAV flight envelope protection control law, the variable that requires particular focus for protection is the angle of attack. Therefore, the weight values related to the angle of attack in Table 10 are larger than the other control variable weight values, which effectively prevents large angle of attack commands and excessive rates of change in the angle of attack. Moreover, to ensure that the UAV can exit the icing area along a stable route, the primary output variable of interest is the flight path angle γ. The weight for the flight path angle in Table 10 is set higher than those for the other two output variables. The limiting boundaries for the angle of attack and roll angle in Table 10 are calculated in real time on the icing severity detected by the neural network. Additionally, to mitigate the risk of loss of control due to large roll angles during UAV flight, limit boundaries for the roll angle are imposed according to fuzzy risk level boundaries [33].

6.2. Mathematical Simulation Task Design and Analysis

The initial state of the UAV is set with a velocity of V = 94 m/s, an altitude of H = 6000 m, and an icing severity of η = 0.9. The UAV remains in level flight and detects the icing severity via a neural network. When the icing severity detection value increases slowly over time, the UAV is considered to have completed the icing severity detection, and the control law parameters are adjusted on the basis of the detection value. Then, the safe tracking command for exiting the icing area at this icing severity is selected on the basis of the interpolation from Figure 11, verifying the effectiveness of the ice tolerance flight envelope protection control method, which uses the LSTM neural network to detect the icing severity. The research presented in this paper was conducted using MATLAB R2023b. The mathematical simulation results are shown in Figure 17.

The analysis in Figure 17b,c reveals that the UAV maintains a level flight state during the initial 5 s, whereas the neural network detects the icing severity on the basis of the UAV’s flight data. The analysis in Figure 17i reveals that there is some error between the detection value and the reference value during the first 5 s of icing severity detection. This is due to the limited amount of flight data input to the neural network at the beginning of the detection, which leads to some error in the network’s detection results. As time progresses, the amount of flight data input into the neural network increases, enhancing the amount of information captured by the network. The neural network’s detection value converges to the reference value within approximately 5 s, after which it increases gradually with time. At this point, the neural network is considered to have completed the icing severity detection.

The control law is reconstructed on the basis of the results of the neural network detection, and its parameters are adjusted accordingly. The inner loop α and μ commands in Figure 14 are generated by the NMPC controller on the basis of the outer loop trajectory commands, whereas the sideslip angle command is maintained at zero to prevent the generation of undesired side forces, yaw, and roll moments. The angle of attack limit boundary is determined by considering two factors: (1) the stall angle of attack, which is calculated on the basis of the detected icing severity, and (2) the angle of attack limit value, which is derived from the fuzzy risk boundary. Considering the two factors mentioned above, the boundary value for the angle of attack in the simulation task of detaching from the icing area is set to 9°. From the aerodynamic characteristics of the icing UAV shown in Figure 3, it can be observed that icing causes the stall angle of attack to decrease gradually, enhancing aerodynamic nonlinearity. The UAV’s aerodynamic characteristics remain relatively stable when the angle of attack is within the 9° range. Therefore, setting the limiting value of the angle of attack to 9° during the process of detaching from the icing area is reasonable. The roll angle boundary value is set by the fuzzy risk boundary with a maximum value of 15°.

Then, on the basis of the fuzzy-risk-boundary-based command risk level calculation method established in Section 4, the risk level for different trajectory tracking commands at this icing level can be calculated. If the selected command is too small, it will not fully utilize the UAV’s flight performance, despite having a low risk level, and will increase the time the UAV spends in the icing area, thereby increasing the flight risk. Therefore, commands with a risk level within the level two risk boundary are selected for the tracking command library, from which commands are chosen as shown in Equation (22).

$$\left\{\begin{array}{l}V=100\mathrm{m}/\mathrm{s}\\ \gamma ={10}^{\circ }\\ \chi ={10}^{\circ }\end{array}\right.$$

(22)

From the analyses in Figure 17a–c, it can be observed that, in the case of the icing UAV using the NDIC, although the UAV ultimately completes the commanded tracking, the NDIC generates a large angle of attack response to ensure the rapidity of tracking, as shown in Figure 17d, with the maximum angle of attack reaching 12°. Figure 3 shows that the deterioration of aerodynamic characteristics in icing UAV becomes more pronounced at larger angles of attack, which has a greater impact on flight safety. Moreover, since the NDIC relies on an accurate aerodynamic model, the gap between the aerodynamic characteristics of the icing UAV and the clean configuration becomes large at higher angles of attack, resulting in poor control performance of the NDIC. Figure 17b,d,h shows that the trajectory angle, angle of attack, and elevator responses all oscillate. If the design parameters of the control law reference model are too large, it may lead to the loss of control of the UAV, thereby seriously endangering flight safety.

In the case of icing UAV using NMPC controllers as trajectory loop controllers, the NMPC controllers solve the optimization problem for a future period on the basis of designed parameters such as error weights, control weights, control command rate of change weights, and limiting boundaries. Under the premise of meeting the design requirements, the current angle of attack, throttle, and roll angle commands are generated, and the prediction model is corrected on the basis of the actual and predicted response values of the UAV. With the angle of attack restricted, the lift generated by increasing the angle of attack is limited, making it difficult to quickly complete the flight path angle tracking task. To ensure rapid tracking, the NMPC controller increases the flight speed by adjusting the throttle command, as shown in Figure 17k, thereby generating more lift and thrust. The tracking of the flight path angle is achieved quickly and accurately while ensuring that the angle of attack does not exceed the limiting boundary, as shown in Figure 17b,d. The tracking of the flight path heading angle also ensures both rapidity and accuracy while meeting the design requirements, as shown in Figure 17c. After the flight path angles are tracked, the NMPC controller reduces the throttle command to track the flight speed.

Figure 17 shows that the ice tolerance flight envelope protection control system based on the LSTM neural network for detecting icing severity designed in this paper is highly effective. It enables quick, accurate, and stable command tracking while ensuring flight envelope protection, thereby making full use of the aircraft’s performance under the premise of ensuring flight safety.

It is important to note that the ice-tolerant flight envelope protection control system design method based on an LSTM neural network for detecting icing severity proposed in this paper relies on the established aerodynamic model of the icing UAV. If the aerodynamic model of the icing UAV is inaccurate, the effectiveness of both the detection and control methods may be compromised. Therefore, developing an accurate aerodynamic model of icing UAVs through wind tunnel experiments and other methods is a crucial area for further investigation.

7. Conclusions

(1)

A detection method based on the LSTM neural network is established to address the icing severity detection problem for an icing UAV. The method does not depend on active excitation and identifies changes in the UAV’s aerodynamic characteristics after icing on the basis of its historical flight data over a certain period, thereby enabling the detection of the icing severity. The method demonstrated high accuracy and can be effectively applied to command risk level calculations and the design of control law reconstruction.

(2)

The existing fuzzy risk boundaries of flight parameters for a clean-configuration UAV were revised and recalculated on the basis of the UAV’s aerodynamic model and the changes in its aerodynamic characteristics after icing. Thus, a fuzzy-risk-boundary-based command risk level calculation method for an icing UAV was established. According to the flight task, the tracking command was selected as the trajectory command, and 3D command risk level topology map with clean configuration and icing severity 1 was calculated and its change rule was analyzed, thereby providing a library of safe tracking commands for UAVs to perform their tasks.

(3)

To address the icing UAV ice tolerance flight control problem, an ice tolerance flight envelope protection control law based on the NMPC controller was developed. Compared with the NDIC, the control system designed in this paper can achieve fast, accurate, and stable command tracking while also providing protection for critical state variables.

In this paper, the proposed ice tolerance flight envelope protection control system based on an LSTM neural network for detecting icing severity is simple, in principle, with high detection accuracy, and the designed control law can give full play to the UAV flight performance while guaranteeing flight safety.