Improved Thrust Performance Optimization Method for UAVs Based on the Adaptive Margin Control Approach

Abstract

This study proposes a strategy for improving the thrust performance of fixed-wing UAV turbine engines from the perspective of aircraft/engine integration. In the UAV engine control process, the inlet distortion caused by the angle of attack change is taken into account, the inlet distortion index is calculated in real time by predicting the angle of attack, and the influence of the inlet distortion on the engine model is analyzed mechanically. Then, the pressure ratio command is adjusted according to the new compressor surge margin requirement caused by the inlet distortion to finally improve the engine thrust performance. To verify the effectiveness of the algorithm, an adaptive disturbance rejection controller is designed for the flight control of a fixed-wing UAV to complete the simulation of horizontal acceleration. The simulation results show that, with this strategy, the UAV turbofan engine can improve the turbofan engine thrust performance by more than 8% under the safety conditions.

1. Introduction

The unmanned aerial vehicle (UAV), as a rapidly developing aircraft technology, plays an important role in many fields [1,2,3], such as large-scale display activities, aerial photography entertainment, mapping and geological exploration, agricultural spraying, police patrol, electric power inspection, etc. [4]. UAV propulsion systems have a decisive influence on the performance, cost, and reliability of UAVs. With the continuous enrichment of application scenarios, different types of power units have been developed for different requirements of UAVs in terms of flight speed, flight altitude, maneuvering overload, landing and take-off methods, range, and economic indicators [5]. Most of the power units currently installed on advanced UAVs are turbofan engines, and this status quo will not change for a long time to come. Li conducted work related to the prediction of future UAV power development trends [6], summarized the overloaded advanced typical combat UAVs and the main parameters of power, and compared the payloads of several types of UAVs under development with those of manned fighter aircraft in service. The study shows that the development of large-thrust turbofan engines is one of the future unmanned fighter power development trends. Hu analyzed the requirements of various unmanned aircraft systems for propulsion systems and the impact of propulsion system technology on the performance of UAVs [7,8]. Chen analyzed the development requirements of different types of engines for UAVs and proposed the integrated design of the full authority digital engine control (FADEC) system of a vehicle/engine/propeller for optimal performance matching output [9]. In summary, the research on the improvement of the thrust performance of future UAV turbine engines has considerable value, and this part of the work is not abundant at present. Cheng [10] proposed a turbofan engine Turbine Inter Burner (TIB) supplemental combustion and thrust augmentation thermodynamic cycle scheme based on ultra-compact combustion UCC technology, increasing the engine unit thrust significantly. Dong [11] analyzed the effect of jet pre-cooling on the high- and low-pressure shaft torque, obtained the characteristics of the effect of jet pre-cooling on the control scheme, and also designed the relevant control scheme to adjust the nozzle and buck ratio limiting strategies, thus improving the engine speed performance. Sun [12] established a piecewise linear model of the turbofan engine and used a weighted prediction method to effectively estimate the thrust of the engine and track the set expected thrust to achieve maximum compensation for the thrust loss, and the simulation results showed that the whole thrust control system can ensure the safe operation of the engine and achieve effective thrust boost.

The existing related work lacks the ability to improve powertrain performance from an aircraft/engine integration perspective, and the integrated control of UAVs and engines can seek the best match of their performance. In the aircraft/engine integration research, the stability control technology of the turbo-fan engine is essential. The influence of inlet distortion on the performance and stability of turbine engines is first considered, which is a common phenomenon that occurs widely in turbo-fan engines; the geometric asymmetry of the inlet airflow channel of the compressor, the obstruction of the local inlet airflow by the leading edge of the compressor inlet when the aircraft is flying at a large angle of attack, and the uneven distribution of the temperature and density of the compressor inlet flow field by the missile trailing combustion airflow are the main reasons [13,14]. In order to solve the problem of the unmeasurable engine surge margin in super maneuvering flight, Wang [15] proposed a modeling method of the engine surge margin and verified that the above scheme can accurately control the engine surge margin at 11~13% by the numerical simulation of maneuvering flight with a large angle of attack, which ensures the stability and high efficiency of engine operation. Liu [16] concluded that inlet distortion is an important factor affecting engine stability. To meet the inlet distortion requirements, the engine design needs to make large concessions; in other words, the inlet distortion required for the surge margin greatly reduces the performance of the design point of the compression system. This reduction in margin translates directly into improved system capabilities.

Indeed, the current strategy, which considering flight control and engine control separately, is capable of achieving the optimal performance of each subsystem. However, the overall performance cannot be optimized since the change in engine operation caused by the angle of attack is not taken into account. A large margin of surge is usually given in the control in order to ensure the safety of the compressor operation, which sacrifices a certain engine performance. Inspired by the development of aero-engine stability control technology, this paper proposes a method for improving the thrust performance of future UAVs. This method incorporates engine stability assessment into the engine control process, calculates the requirement of the compressor surge margin associated with the engine stability influence factor in real time, and evaluates the compressor stability, allowing the control system to minimize the surge margin to improve the engine performance and thus enhance the engine thrust performance.

2. Dynamical Modeling and Problem Description

In order to evaluate the UAV thrust performance improvement, UAV flight simulation under certain flight profiles is required, so the UAV dynamics model is first established.

2.1. UAV Dynamics Modeling

2.1.1. UAV Modeling Assumptions

Considering that the mass of the turbine engine-powered UAV is time-varying during flight and its structure is strong and elastic, there are also centrifugal and Gauche accelerations; thus, the aerodynamic forces acting on the UAV are affected by many variables such as its aerodynamic layout, structural elasticity, and flight state, which sometimes cannot be modeled mathematically accurately. Therefore, the primary factors should be considered when modeling the UAV, ignoring the secondary factors. Figure 1 shows a diagram of the coordinate system of the UAV body. The following assumptions need to be considered in the mathematical modeling of scaled-down UAVs.

• The mass of the UAV does not change during flight for a short period of time.
• The UAV is an ideal rigid body, and the effects of aircraft elasticity are ignored.
• The UAV is affected only by the acceleration of gravity, and the curvature of the Earth is neglected.
• Since the UAV is symmetric about the $x-y$ plane, $I_{xz}=0$ and $I_{yz}=0$.
• The change in the aircraft center of mass with fuel is neglected during the UAV.

2.1.2. Dynamics Equations

The dynamics equations for the motion of the center of mass of the UAV are established as follows.

$$\left\{\begin{array}{c}m\dot{V}=P\cos \alpha \cos \beta -X-mg\sin \theta \\ mV\dot{\theta }=P\left(\sin \alpha \cos {\gamma }_V+\cos \alpha \sin \beta \sin {\gamma }_V\right)+Y\cos {\gamma }_V-Z\sin {\gamma }_V-mg\cos \theta \\ -mV{\dot{\psi }}_V\cos \theta =P\left(\sin \alpha \sin {\gamma }_V-\cos \alpha \sin \beta \cos {\gamma }_V\right)+Y\sin {\gamma }_V+Z\cos {\gamma }_V\end{array}\right.$$

(1)

where m is the UAV mass, V is the UAV velocity, X,Y,Z are the drag lift and lateral force, respectively, P is the UAV thrust, α, β are the angle of attack and sideslip angle, θ is the flight path angle, γ_V is the speed tilt angle, and ψ_V is the trajectory deflection angle.

The dynamics equations of the UAV rotational motion are as follows

$$\begin{array}{c}{\dot{\omega }}_x=\left(-\frac{\left(I_{zz}-I_{yy}\right)I_{yy}-I_{xy}^2}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_y-\frac{\left(I_{xx}-I_{zz}+I_{yy}\right)I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_x\right)+\frac{I_{yy}}{I_{xx}{I}_{yy}-I_{xy}^2}{M}_x+\frac{I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}{M}_y\\ {\dot{\omega }}_y=\left(\frac{\left(I_{xx}-I_{zz}+I_{yy}\right)I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_y-\frac{\left(I_{xx}-I_{zz}\right)I_{xx}+I_{xy}^2}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_x\right){\omega }_z+\frac{I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}{M}_x+\frac{I_{xx}}{I_{xx}{I}_{yy}-I_{xy}^2}{M}_y\\ {\dot{\omega }}_z=-\frac{I_{yy}-I_{xx}}{I_{zz}}{\omega }_x{\omega }_y+\frac{I_{xy}}{I_{zz}}\left({\omega }_x^2-{\omega }_y^2\right)+\frac{M_z}{I_{zz}}\end{array}$$

(2)

where ${\omega }_x,{\omega }_y,{\omega }_z$ are the angular velocities of rotation of the UAV relative to the ground coordinate system, $M_x,M_y,M_z$ are the components of the moments of all external forces acting on the UAV on the center of mass in the UAV body coordinate system, and $I_{xz},I_{yy},I_{zz},I_{xy}$ are inertia tensor elements.

2.1.3. Kinematic Equations

In the modeling, it is assumed that the UAV is a rigid body and that the six-degrees-of-freedom motion of the UAV is the displacement and rotation around the center of mass. With the help of the transformation matrix of the ground inertial system and the aircraft dynamic system, the kinematic equations are obtained as follows [17,18].

$$\left\{\begin{array}{c}\dot{x}=V\cos \theta \cos {\psi }_V\\ \dot{y}=V\sin \theta \\ \dot{z}=-V\cos \theta \sin {\psi }_V\\ \dot{\vartheta }={\omega }_y\sin \gamma +{\omega }_z\cos \gamma \\ \dot{\psi }=\left({\omega }_y\cos \gamma -{\omega }_z\sin \gamma \right)/\cos \vartheta \\ \dot{\gamma }={\omega }_x-\tan \vartheta \left({\omega }_y\cos \gamma -{\omega }_z\sin \gamma \right)\end{array}\right.$$

(3)

where $x,y,z$ are the center-of-mass position components of the man-machine, $\vartheta$ is the pitch angle, and $\gamma$ is the roll angle.

2.2. UAV Turbofan Engine Modeling

The object of this study is a twin-shaft small-culvert-ratio mixed-exhaust turbofan engine containing no afterburner. The main components of the engine include the inlet, fan, compressor, burner, high-pressure turbine, low-pressure turbine, bypass, mixer, and nozzle; its structural schematic diagram is shown in Figure 2.

The aerodynamic and thermal processes during the operation of an aero-engine are very complex, and, when modeling them, the real physical processes need to be suitably simplified, and the following assumptions are usually made when modeling at the component level [19].

• Ignore the effect of the Reynolds number variation on the characteristics of the engine components; ignore the effect of atmospheric humidity on the parameters of the thermodynamic properties of the gas; the effect of rotor inertia is considered.
• The flow of air in the engine is treated as quasi one-dimensional flow, i.e., the airflow parameters are considered to be the same in the same cross-section of the engine.

The idea of component-level modeling is to first establish each component sub-model; then, build a non-linear system of equations based on the aerodynamic thermodynamic laws governing the gas flow process, combined with the assumptions of flow continuity and energy balance; finally, determine the engine operating state by solving the system of common operating equations to construct the engine model.

When a small-culvert-ratio turbofan engine enters into quasi-steady-state operation, each component needs to satisfy a series of common operating conditions in mechanical and aerodynamic aspects, including flow continuity and power balance. It is generally accepted that the inlet and outlet sections of each component meet the flow continuity condition and the two rotors meet the power balance condition, i.e., the power from the high-pressure turbine is balanced with the power consumed by the fan and booster stage, and the power from the low-pressure turbine is balanced with the power consumed by the high-pressure compressor. The dynamic engine simulation process also requires an iterative solution of the common set of operating equations, unlike the steady-state simulation, where the flow continuity balance condition is still satisfied, but the power balance condition is no longer satisfied. The high- and low-pressure rotor speeds are no longer used as initial guess variables but are determined by the speed and acceleration at the previous moment, i.e., the power balance conditions are transformed into the rotor dynamics equations.

2.2.1. Inlet

As the inlet section of the engine, the calculation of the relevant parameters of the air inlet tract is related to the height of the engine. The main parameters of the inlet cross-section—atmospheric pressure $P_H$ in $P_a$ and atmospheric temperature $T_H$ in $K$ at different altitudes (H km)—can be obtained by the atmoscoesa function in MATLAB. The inlet pressure $p_1^{*}$, temperature $T_1^{*}$, and flight speed $c_0$ of the engine are calculated as follows.

$$\left\{\begin{array}{c}{p}_1^{*}={\sigma }_i{p}_H{\left(1+\frac{k-1}{2}{M}_H^2\right)}^{\frac{k}{k-1}}\\ T_1^{*}=T_H\left(1+\frac{k-1}{2}{M}_H^2\right)\\ c_0=M_H\times a_H=M_H\sqrt{kR{T}_H}\end{array}\right.$$

(4)

in the above equation, $M_H$ is the Mach number; ${\sigma }_i$ is the total pressure recovery factor, where ${\sigma }_i={\sigma }_i\left(M_H\right)$; k is the adiabatic index; and $a_H$ is the local speed of sound.

2.2.2. Fan

Assuming the same main and bypass flow pressure ratios and efficiencies of the fan, we have:

$${\pi }_{cL}^{*}={\pi }_{cII}^{*},\hspace{1em}{\eta }_{cL}^{*}={\eta }_{cII}^{*}$$

(5)

where ${\pi }_{cL}^{*}$,${\pi }_{cII}^{*}$ represent the pressure ratios of the main and bypass flow and ${\eta }_{cL}^{*}$,${\eta }_{cII}^{*}$. represent the overall efficiencies of the main and bypass flow.

The exit section airflow parameters of the turbofan engine fan are calculated as follows.

$$\begin{gathered} p_{1.5}^{*}={\pi }_{cL}^{*}{p}_1^{*} \\ {p}_2^{*}=p_{1.5}^{*} \\ {m}_a=m_{aL,cor}\frac{p_1^{*}}{101325}\sqrt{\frac{288}{T_1^{*}}} \\ {T}_{1.5}^{*}=T_1^{*}\left(1+\frac{{\pi }_{cL}^{*}{}^r-1}{{\eta }_{cL}^{*}}\right) \\ {T}_{2II}^{*}=T_{1.5}^{*} \\ {\eta }_{cL}^{*}={\eta }_{cL}^{*}\left(n_{cL,cor},{\pi }_{cL}^{*},{\varphi }_L\right) \\ {m}_{aL,cor}=m_{aL,cor}(n_{cL,cor},{\pi }_{cL}^{*},{\varphi }_L) \end{gathered}$$

(6)

in the above equation, $r=\frac{k-1}{k}$, where $k=k\left(T_{1.5}^{*}\right)$; $m_a$ is the total air flow through the fan; $m_{aL,cor}$ is the converted air flow; ${\varphi }_L$ is the fan adjustable parameter; ${\eta }_{cL}^{*}$ is the converted speed at different temperatures. The last two equations in Formula (6) show that the specific values of $m_{aL,cor}$ and ${\eta }_{cL}^{*}$ are obtained by interpolating the variables associated with each of them.

2.2.3. Compressor

Similar to the fan characteristics, the generic characteristics of a compressor are calculated as follows.

$$\begin{gathered} m_{aH,cor}=m_{aH,cor}(n_{cH,cor},{\pi }_{cH}^{*},{\varphi }_H) \\ {\eta }_{cH}^{*}={\eta }_{cH}^{*}\left(n_{cH,cor},{\pi }_{cH}^{*},{\varphi }_H\right) \end{gathered}$$

(7)

where $m_{aH,cor}$ is the converted air flow; ${\varphi }_H$ is the compressor adjustable parameter; $n_{cH,cor}$ is the converted speed of the compressor.

If the fan control law ${\varphi }_H={\varphi }_H\left(n_{cH}^{*}\right)$ is known, then the mathematical relationship between the conversion speed and the rotational speed $n_{cH,cor}=n_H\sqrt{288/T_{1.5}^{*}}$, and we can obtain $n_{cH,cor}$ and ${\pi }_{cH}^{*}$. Then, from $n_{cH,cor},$ ${\pi }_{cH }^{*}$, and ${\varphi }_H$, we can calculate $m_{aH,cor}$ and ${\eta }_{cH}^{*}$, where the relevant parameters for the outlet section of the compressor are calculated as follows:

$$\begin{gathered} p_2^{*}={\pi }_{cH}^{*}{p}_{1.5}^{*} \\ {m}_{a2}=m_{aH,cor}\frac{p_{1.5}^{*}}{101325}\sqrt{\frac{288}{T_{1.5}^{*}}} \\ {T}_2^{*}=T_{1.5}^{*}\left(1+\frac{{\pi }_{cH}^{*}{}^r-1}{{\eta }_{cH}^{*}}\right) \end{gathered}$$

(8)

2.2.4. Burner

The burner characteristics are generally given by the component characteristics test; in general, the burner efficiency ${\eta }_b$ and the total pressure recovery coefficient ${\sigma }_b$ are the main characteristics parameters of the combustion chamber and are calculated as follows.

$$\begin{gathered} {\eta }_b={\eta }_b\left(\tilde{\alpha },p_2^{*},T_2^{*},T_3^{*}\right) \\ {\sigma }_b={\sigma }_b\left(c_b,\theta \right) \\ \tilde{\alpha }=\frac{m_{a2}}{m_f·14.8} \\ \theta =T_3^{*}/T_2^{*} \end{gathered}$$

(9)

in the above equation, $\tilde{\alpha }$ is the residual gas coefficient of the gas mixture; $c_b$ is the gas flow rate; $\theta$ is the heating ratio.

The ${\eta }_b$ is calculated from $\tilde{\alpha },p_2^{*},T_2^{*}$, and $T_3^{*}$ (the burner characteristics), and then the $T_3^{*}$ is calculated from the energy balance (by iteration). The energy balance equation is calculated as follows.

$$\frac{m_{a2}}{m_f}\left[H_u{\eta }_b+h_f\left(T_{f0}\right)-h_f\left(T_3^{*}\right)\right]-h_g\left(T_3^{*}\right)+h_a\left(T_2^{*}\right)=0$$

(10)

in the formula above, $T_{f0}$ is the fuel temperature, and $h_a$, $h_f$, and $h_g$ denote the enthalpy per unit (kg) of mass of air, fuel, and fuel after vaporization, respectively.

In summary, the combustion chamber output parameters can be calculated as follows.

$$\begin{gathered} p_3^{*}={\sigma }_b{p}_2^{*} \\ {m}_{g3}=m_{a2}+m_f \end{gathered}$$

(11)

2.2.5. High-Pressure Turbine (HPT) and Low-Pressure Turbine (LPT)

The efficiencies of the high-pressure turbine ${\eta }_{TH}^{*}$, the gas flow rate $m_{gH,cor}$, and the conversion speed $n_{TH,cor}$ are related to the expansion ratio ${\pi }_{TH}^{*}$, i.e., its operating state can be determined from ${\pi }_{TH}^{*}$ and $n_{TH,cor}$. The mathematical representation of $n_{TH,cor}$ and $m_{gH,cor}$, obtained from the characteristic curve diagram, is as follows.

$$\begin{gathered} m_{gH,cor}=m_{gH,cor}(n_{TH,cor},{\pi }_{TH}^{*}) \\ {\eta }_{TH}^{*}={\eta }_{TH}^{*}\left(n_{TH,cor},{\pi }_{TH}^{*}\right) \end{gathered}$$

(12)

where $n_{TH,cor}=n_H\sqrt{288/T_3^{*}}$, and $T_3^{*}$ is the engine turbine inlet section temperature in the design condition.

In summary, the high-pressure turbine outlet parameters are calculated as follows.

$$\begin{gathered} m_{g3.5}=m_{gH,cor}\frac{p_3^{*}}{101325}\sqrt{\frac{288}{T_3^{*}}} \\ {p}_{3.5}^{*}=\frac{p_3^{*}}{{\pi }_{TH}^{*}} \\ {T}_{3.5}^{*}=T_3^{*}\left[1-\left(1-{\pi }_{TH}^{*}{}^{-r\prime }\right){\eta }_{TH}^{*}\right] \end{gathered}$$

(13)

in the above equation, $r^{\prime }=\frac{k^{\prime }-1}{k^{\prime }}$, where $k^{\prime }=k^{\prime }\left(T_3^{*},\tilde{\alpha }\right)$ is the adiabatic coefficient of fuel oil after vaporization.

The structures of the high-pressure turbine and the low-pressure turbine are the same, as are the models, and the parameters of the low-pressure turbine outlet section can be obtained.

$$\begin{gathered} m_{g4}=m_{gL,cor}\frac{p_{3.5}^{*}}{101325}\sqrt{\frac{288}{T_{3.5}^{*}}} \\ {p}_4^{*}=\frac{p_{3.5}^{*}}{{\pi }_{TL}^{*}} \\ {T}_4^{*}=T_{3.5}^{*}\left[1-\left(1-{\pi }_{TL}^{*}{}^{-r\prime }\right){\eta }_{TL}^{*}\right] \end{gathered}$$

(14)

in the above equation, $r^{\prime }=\frac{k^{\prime }-1}{k^{\prime }}$, and $k^{\prime }=k^{\prime }\left(T_{3.5}^{*},\tilde{\alpha }\right)$ is the thermal insulation index. $m_{gL,cor}$, ${\pi }_{TL}^{*}$, and ${\eta }_{TL}^{*}$ correspond to $m_{gH,cor}$, ${\pi }_{TH}^{*}, \mathrm{and}$ ${\eta }_{TH}^{*}$, representing the characteristics of the low-pressure turbine.

2.2.6. Mixer

The gas flow rate contained in the mixing chamber is calculated as follows:

$$m_{g4}=K_g^{\prime }\frac{p_4^{*}{A}_{4I}q\left({\lambda }_4\right)}{\sqrt{T_4^{*}}}$$

(15)

Equivalently, the above equation is transformed as:

$$q\left({\lambda }_4\right)=\frac{m_{g4}\sqrt{T_4^{*}}}{p_4^{*}{A}_{4I}{K}_g^{\prime }}$$

(16)

where $K_g^{\prime }=\sqrt{\frac{k^{\prime }}{R^{\prime }}{\left(\frac{2}{k^{\prime }+1}\right)}^{\frac{k^{\prime }+1}{k^{\prime }-1}}}$, $R^{\prime }$ is the gas constant, $k^{\prime }=k^{\prime }\left(T_4^{*},\tilde{\alpha }\right)$, and $A_{4I}$ is the cross-sectional area of the inner culvert (main flow) into the mixing room.

Additionally, since $q\left({\lambda }_4\right)$ can be expressed as

$$q\left({\lambda }_4\right)={\left(\frac{k^{\prime }+1}{2}\right)}^{\frac{1}{k^{\prime }-1}}{\lambda }_4{\left(1-\frac{k^{\prime }-1}{k^{\prime }+1}{\lambda }_4{}^2\right)}^{\frac{1}{k^{\prime }-1}}$$

(17)

the above equation leads to ${\lambda }_4$, which, in turn, gives

$$\pi \left({\lambda }_4\right)={\left(1-\frac{k^{\prime }-1}{k^{\prime }+1}{\lambda }_4{}^2\right)}^{\frac{1}{k^{\prime }-1}}$$

(18)

$$f\left({\lambda }_4\right)=\left({\lambda }_4{}^2+1\right){\left(1-\frac{k^{\prime }-1}{k^{\prime }+1}{\lambda }_4{}^2\right)}^{\frac{1}{k^{\prime }-1}}$$

(19)

$$p_4=\pi \left({\lambda }_4\right)p_4^{*}$$

(20)

Again, the main flow and bypass flow pressures are equal, i.e., $p_{4II}=p_4$, so we have

$$\pi \left({\lambda }_4\right)=\frac{ p_{4II}}{p_{4II}^{*}}=\frac{p_4}{{\sigma }_{II}{p}_{2II}^{*}}$$

(21)

$${\lambda }_{4II}=\sqrt{\frac{k+1}{k-1}\left\{1-{\left[\pi \left({\lambda }_{4II}\right)\right]}^{\frac{1}{k-1}}\right\}}$$

(22)

So, there is

$$q\left({\lambda }_{4II}\right)={\left(\frac{k+1}{2}\right)}^{\frac{1}{k-1}}{\lambda }_{4II}{\left(1-\frac{k-1}{k+1}{\lambda }_{4II}{}^2\right)}^{\frac{1}{k-1}}$$

(23)

$$f\left({\lambda }_{4II}\right)=\left({\lambda }_{4II}{}^2+1\right){\left(1-\frac{k+1}{k-1}{\lambda }_{4II}{}^2\right)}^{\frac{1}{k-1}}$$

(24)

$$m_{4II}=K_q^{\prime }\frac{p_{4II}^{*}{A}_{4II}q\left({\lambda }_{4II}\right)}{\sqrt{T_{2II}^{*}}}$$

(25)

where $A_{4II}$ indicates the area of the bypass duct at the entrance to the mixing room.

The mixer outlet pressure is calculated as follows.

$$p_{cm}^{*}={\sigma }_{cm}\frac{m_{g4}{p}_4^{*}+m_{aII}{p}_{4II}^{*}}{m_{g4}+m_{aII}}$$

(26)

where ${\sigma }_{cm}$ is the total pressure recovery coefficient of the mixer.

${\lambda }_{cm}$ can be calculated from $f\left({\lambda }_{cm}\right)$, and then $\pi \left({\lambda }_{cm}\right)$ and $q\left({\lambda }_{cm}\right)$ can be calculated from the pneumatic function calculation formula. In summary, the gas flow and temperature at the outlet section of the mixing chamber are calculated as follows.

$$\begin{gathered} m_{gcm}=m_{g4}+m_a-m_{aH} \\ {T}_{cm}^{*}=\frac{c_p{T}_{2II}^{*}{m}_{aII}+c_p^{\prime }{T}_4^{*}{m}_{g4}}{c_p^{″}{m}_{gcm}} \end{gathered}$$

(27)

where $c_p=c_p\left(T_{2II}^{*}\right)$, $c_p^{\prime }=c_p^{\prime }\left(\tilde{\alpha },T_4^{*}\right)$, and $c_p^{″}=c_p^{″}\left(\tilde{\alpha },T_{cm}^{*}\right)$.

2.2.7. Nozzle

According to the different relationship between the ratio of the total pressure of the nozzle outlet $p_e^{*}$ and the total atmospheric pressure $p_H$ and the nozzle pressure ratio ${\pi }_{ecr}$, the engine nozzle can be divided into subcritical, critical, and supercritical operating conditions.

When $\frac{p_e^{*}}{p_H}<{\pi }_{ecr}$, the nozzle is operating in a sub-critical condition. At this point, the calculation of $p_e^{*}$ and ${\pi }_{ecr}$ are as follows.

$$p_e^{*}={\sigma }_e{p}_{cm}^{*}$$

(28)

$${\pi }_{ecr}=\frac{p_e^{*}}{p_{cr}}={\left(\frac{k^{\prime }+1}{2}\right)}^{\frac{k^{\prime }}{k^{\prime }-1}}$$

(29)

At this point, the exit section of the nozzle airflow velocity is less than the speed of sound, so ${\lambda }_e$< 1. The gas in the nozzle is completely expanded; at this time, the nozzle exit section pressure and the external atmospheric pressure are the same, that is, $p_e=p_H$.

From $p_e^{*}$ and $p_e$, we have

$$\pi \left({\lambda }_e\right)=\frac{p_e^{*}}{p_e}=\frac{p_e^{*}}{p_H}$$

(30)

from $\pi \left({\lambda }_e\right)$, we obtain ${\lambda }_e$ and find $q\left({\lambda }_e\right)$.

When $\frac{p_e^{*}}{p_H}={\pi }_{ecr}$, the nozzle works in the critical state; at this time, ${\lambda }_e$ = 1 and $q\left({\lambda }_e\right)$ = l. At this time, the gas in the nozzle is completely expanded, and the nozzle outlet section pressure and the external atmospheric pressure are the same, that is, $p_e=p_H$.

When $\frac{p_e^{*}}{p_H}>{\pi }_{ecr}$, the nozzle operation is in the supercritical state; at this time, ${\lambda }_e$ = 1 and $q\left({\lambda }_e\right)$ = l. At this time, the gas in the nozzle is not fully expanded, and the nozzle outlet cross-section pressure is greater than the atmospheric pressure, that is, $p_e>p_H$.

Therefore, its operating condition can be determined from the nozzle parameters (${\lambda }_e, q\left({\lambda }_e\right), p_e$), and the parameters related to the exit section of the tail nozzle are calculated as follows.

$$m_{ge}=K_q^{\prime }\frac{p_e^{*}{A}_{e}q\left({\lambda }_e\right)}{\sqrt{T_{cm}^{*}}}$$

(31)

$$c_e={\phi }_e{\lambda }_e{a}_{cr}={\phi }_e{\lambda }_e\sqrt{\frac{2k^{\prime }}{k^{\prime }+1}{R}^{\prime }{T}_{cm}^{*}}$$

(32)

In the above equation, $m_{ge}$ is the nozzle outlet flow rate, $c_e$ is the exhaust velocity, $a_{cr}$ is the critical speed of sound, and ${\phi }_e$ is the nozzle flow loss coefficient number.

2.3. Problem Description

As mentioned in the introduction, the aim of this paper is to improve the thrust performance of the UAV by reducing the surge margin reserve of the turbine engine. The real-time estimation of the surge margin is caused by many factors, including inlet distortion, and calculations of the proper pressure ratio of the fan and the compressor are carried out to satisfy the surge margin requirements through multivariable control.

From a mathematical point of view, the thrust lifting problem under this line of thought can be described as the optimization problem of engine thrust regulation control.

$$\left\{\begin{array}{cc}\begin{array}{c}\mathrm{Max}\\ \mathrm{s}.\mathrm{t}.\end{array}& \begin{array}{c}{P}_{opt}\left({\pi }_{cH}^{*},m_{aH,cor}, A_8,H, Ma,N_{fan}\right)\\ S{M}_{opt}<S{M}_{re}\end{array}\\ & T_{3.5}^{*}<T_{cr}\\ & p_2^{*}<p_{2,cr}^{*}\end{array}\right.$$

(33)

where $P_{opt}$ is the current state of the engine output thrust; in this optimization problem, the engine output thrust considered is related to the compressor pressurization ratio ${\pi }_{cH}^{*}$, the compressor converted air flow rate $m_{aH,cor}$, the fan speed $N_{fan}$, the ambient altitude $H$, and the flight Mach number $Ma$. $S{M}_{opt}$ and $S{M}_{re}$ represent the current surge margin and the current allowable surge margin, respectively. $T_{3.5}^{*}$ and $T_{cr}$ are the turbine front temperature and the threshold of the turbine front temperature protection, $p_2^{*}$ and $p_{2,cr}^{*}$ represent the current total compressor outlet pressure and the overpressure protection threshold. $SM$ is described as

$$SM={\left[\frac{{\left(\frac{{\pi }_{cH}^{*}}{m_{aH,cor}}\right)}_s-{\left(\frac{{\pi }_{cH}^{*}}{m_{aH,cor}}\right)}_o}{{\left(\frac{{\pi }_{cH}^{*}}{m_{aH,cor}}\right)}_o}\right]}_{N_c=const}\times 100\%$$

(34)

The problem is constrained by three inequalities: the first is the surge margin of the engine, that is, the compressor does not surge; the second is, in the optimization process, the turbine front temperature cannot exceed the overtemperature protection temperature while the engine is working; the third inequality constraint is that the total compressor outlet pressure of the engine cannot exceed the limit pressure.

In this optimization problem, the independent variables are ${\pi }_{cH}^{*}$,$m_{aH,cor}$,$H$, $Ma$, and $N_{fan}$, where $H$ and $Ma$ are environment variables that change with the motion of the vehicle and the surrounding environment. The value of $N_{fan}$ can be determined by ${\pi }_{cH}^{*}$, so the independent variables that can be adjusted in this optimization problem are ${\pi }_{cH}^{*}$, $m_{aH, \mathrm{and} cor}$, and the dependent variable is the thrust $P_{opt}$.

According to the engine model established in Section 2.2, for a known compressor operating state, we can obtain a set of engine operating parameters. In the engine control, ${\pi }_{cH}^{*}$, as one of the main parameters controlling the compressor operating catch state, can directly affect the working state of the compressor. The physical principle is to change the deflection angle of the guide vane and increase the work carried out by the compressor on the airflow such that the pressure of the airflow at the exit of the compressor is greater. The gas mixture entering the burner carries more energy, so it can increase the overall reasoning performance; however, if ${\pi }_{cH}^{*}$ is too large, it will lead to a reduction in $SM$, that is, it will lead to the occurrence of the gas turbulence of the compressor, so another control parameter $m_{aH,cor}$ should be adjusted at the same time to limit the appearance of dangerous situations.

3. Major Studies

3.1. Thrust Optimization Based on the Adaptive Margin Model

In the introduction, the idea of thrust promotion is introduced, which is to reduce the compressor surge margin reserve and increase the propulsion performance under the premise of ensuring the safe operation of the engine; a detailed description of this strategy is provided below.

As shown in Figure 3, the first step is to predict the angle of attack $\widehat{\alpha }$ in a short period of time in the future by building a dynamic model of the UAV, which is combined with the current mission; then, the distortion index $W_a$ is determined by the current compressor gas conversion flow $m_{aH,cor}$, and the available pressure ratio increment $\Delta {\pi }_{cH}$ is determined by the distortion index and the current surge margin value of the engine; then, the final pressure ratio correction plan is given according to the protection limit of the turbine engine, and the instruction correction quantity is transmitted to the controller to achieve the goal of thrust lifting.

3.1.1. Angle of Attack Prediction for UAV

Based on the symmetry and small perturbation approximations theory [20], the UAV dynamic model is expanded on the point of longitudinal motion equilibrium and rewritten as perturbation equation. The following equations of motion for the longitudinal disturbance of the UAV are given.

$$\begin{gathered} \frac{\mathrm{d}\Delta V}{\mathrm{d}t}-a_{11}\Delta V-a_{13}\Delta \theta -a_{14}\Delta \alpha =a_6{F}_{gx} \\ \frac{{\mathrm{d}}^2\Delta \theta }{\mathrm{d}t}-a_{21}\Delta V-a_{22}\frac{\mathrm{d}\Delta \theta }{\mathrm{d}t}-a_{24}\Delta \alpha -a_{24}^{\prime }\frac{\mathrm{d}\Delta \alpha }{\mathrm{d}t}=a_{25}\Delta {\delta }_z+a_{26}{M}_{gx} \\ \frac{\mathrm{d}\Delta \theta }{\mathrm{d}t}-a_{31}\Delta V-a_{33}\Delta \theta -a_{34}\Delta \alpha =a_{35}\Delta {\delta }_z+a_{36}{F}_{gy} \\ -\Delta \vartheta +\Delta \theta +\Delta \alpha =0 \end{gathered}$$

(35)

where $F_{gx},F_{gy}$ denote the disturbance force due to the disturbance introduced during the small disturbance linearization process and the disturbance torque $M_{gx}$. ${\delta }_z$ denotes the elevator deflection angle. $a_{11}$… $a_{36}$ are known as the power factors, which characterize the dynamics of the UAV. The specific expressions are shown in

Table 1. The symbol of the dynamic coefficient $a_{ij}$ and its expression.

Equations of Motion Number ‘i’	Coefficient of Motion Serial Number ‘j’
	1	2	3	4	5	6
1	$a_{11}=\frac{P^V-X^V}{m}$ $\left({\mathrm{s}}^{-1}\right)$	$a_{12}=0$	$a_{13}=-g\cos \theta$ $\left(\mathrm{m}·{\mathrm{s}}^{-2}\right)$	$a_{14}=-\frac{X^{\alpha }+P^{\alpha }}{m}$ $\left(\mathrm{m}·{\mathrm{s}}^{-2}\right)$	-	$a_{16}=\frac{1}{m}$ $\left({\mathrm{Kg}}^{-1}\right)$
2	$a_{21}=\frac{M_z^V}{I_{zz}}$ $\left({\mathrm{m}}^{-1}·{\mathrm{s}}^{-1}\right)$	$a_{22}=\frac{M_z^{{\omega }_z}}{I_{zz}}$ $\left({\mathrm{s}}^{-1}\right)$	$a_{23}=0$	$a_{24}=\frac{M_z^{\alpha }}{I_{zz}}\left({\mathrm{s}}^{-2}\right)$ $a_{24}^{\prime }=\frac{M_z^{\dot{\alpha }}}{I_{zz}}\left({\mathrm{s}}^{-1}\right)$	$a_{25}=\frac{M_z^{{\delta }_z}}{I_{zz}}\left({\mathrm{s}}^{-2}\right)$ $a_{25}^{\prime }=\frac{M_z^{\dot{{\delta }_z}}}{I_{zz}}\left({\mathrm{s}}^{-1}\right)$	$a_{26}=\frac{1}{I_{zz}}$ $\left({\mathrm{Kg}}^{-1}·{\mathrm{m}}^{-2}\right)$
3	$a_{31}=\frac{P^V\alpha +Y^V}{mV}$ $\left({\mathrm{m}}^{-1}\right)$	$a_{32}=0$	$a_{33}=\frac{g\sin \theta }{V}$ $\left({\mathrm{s}}^{-1}\right)$	$a_{34}=\frac{P+Y^{\alpha }}{mV}$ $\left({\mathrm{s}}^{-1}\right)$	$a_{35}=\frac{Y^{{\delta }_z}}{mV}$ $\left({\mathrm{s}}^{-1}\right)$	$a_{36}=\frac{1}{mV}$ $\left(\mathrm{s}·{\mathrm{Kg}}^{-1}·{\mathrm{m}}^{-1}\right)$

.

In the table above, $X^{\alpha },\cdots ,Y^V,\cdots ,M_z^{\dot{{\delta }_z}},\cdots ,P^V$ represent the partial derivative of parameters such as the aerodynamic torque and thrust with respect to the parameter $\alpha , \cdots ,V,\cdots , \dot{{\delta }_z},\cdots , {\omega }_z$.

In the longitudinal motion of UAV, the outputs of the airframe are $\Delta V,\Delta \vartheta ,\Delta \theta , \mathrm{and} \Delta \alpha$, while the input is $\Delta {\delta }_z$. If there is external interference, the input is $F_{gx}$, $F_{gy}$, and $M_{gz}$. In order to obtain the transfer function of UAV short-period longitudinal perturbation motion, the long-period output $\Delta \mathrm{V}$ is first ignored, and the Equation (40) are transformed by Laplace transform as

$$\left[\begin{array}{ccc}s\left(s-a_{22}\right)& 0& -\left(a_{24}^{\prime }s+a_{24}\right)\\ 0& s-a_{33}& -a_{34}\\ -1& 1& 1\end{array}\right]\left[\begin{array}{c}\Delta \vartheta \left(s\right)\\ \Delta \theta \left(s\right)\\ \Delta \alpha \left(s\right)\end{array}\right]=\left[\begin{array}{c}{a}_{25}\\ a_{35}\\ 0\end{array}\right]\Delta {\delta }_z\left(s\right)+\left[\begin{array}{c}{a}_{26}{M}_{gz}\left(s\right)\\ a_{36}{F}_{gy}\left(s\right)\\ 0\end{array}\right]$$

(36)

Then, the transfer function from the elevator to the angle of attack is obtained as

$$W_{{\delta }_z}^{\alpha }\left(s\right)=\frac{\Delta \alpha \left(s\right)}{\Delta {\delta }_z\left(s\right)}=\frac{-a_{35}{s}^2+\left(a_{35}{a}_{22}+a_{25}\right)s}{s\left[s^2+\left(a_{34}-a_{22}-a_{24}^{\prime }\right)s-\left(a_{34}{a}_{22}+a_{24}\right)\right]}$$

(37)

So far, we have obtained the transfer function from the elevator deflection command to the change in the angle of attack in the short-period longitudinal motion of the UAV, by which the current flight parameters of the UAV and the elevator command can be combined, and the angle of attack can be predicted in a short time.

3.1.2. Engine Model under the Influence of Inlet Distortion

The inlet distortion creates an uneven flow field at the engine inlet, which leads to stall boundary degradation and steady-state characteristic decay of the entire engine, causing the stall point of the entire compressor to move to the lower right [21] (Figure 4), that is to say, the inlet distortion affects the performance parameters of each component of the engine and results in the offset of the common working line and the reduction in the surge margin.

The inhomogeneous air flow caused by the inlet distortion changes the inlet total pressure and affects the thermodynamic process of the subsequent parts. The average total pressure recovery coefficient σ of the inlet cross-section is generally used to measure the inlet total pressure after the inlet distortion. The main formulas are as follows.

$${\sigma }_{AV}=\frac{{\int }_0^{2\pi }{\sigma }_r\left(\theta \right)d\theta }{2\pi }$$

(38)

in the formula, ${\sigma }_r\left(\theta \right)$ is the radial mean total pressure recovery coefficient function, and the expression is

$${\sigma }_r\left(\theta \right)=\frac{{\int }_{{\overline{r}}_{hub}}^{1}2\sigma \left(\overline{r},\theta \right)\overline{r}d\theta }{1-{\overline{r}}_{hub}^2}$$

(39)

where ${\overline{r}}_{hub}$ is the relative radius of the wheel hub and $\sigma \left(\overline{r},\theta \right)$ is the circumferential angle of $\theta$; the total pressure recovery coefficient is a function of the relative radius $\overline{r}$, and the schematic diagram of the parameters is shown in Figure 5. The shaded area is below ${\sigma }_{AV}$ in the low-pressure region ${\theta }^{-}$. ${\theta }_1$ and ${\theta }_2$ are the starting and ending points of the low-pressure region, respectively.

The recovery coefficient of the surface mean total pressure directly reflects the variation in total pressure after distortion [17]. Figure 6 shows the surface mean total pressure recovery coefficient ${\sigma }_{AV}$ and the mean total pressure ${\sigma }_0$ of the distorted low-pressure sector of an inlet at different angles of attack during subsonic and sonic travel. The angle of attack is directly related to the distortion. The larger the angle of attack, the stronger the distortion is. As can be seen from Figure 6, the total pressure recovery coefficient is closely related to the angle of attack. The ${\sigma }_{AV}$ attenuates rapidly from the 0° angle of attack to both sides, and the attenuation is faster at a negative angle of attack. It can be concluded that ${\sigma }_{AV}$ is closely related to the distortion strength and increases with the increase in distortion strength. Therefore, the inlet total pressure in the component-level model needs to be corrected according to the inlet distortion.

In order to quantify the impact of inlet distortion, it is necessary to select a distortion index that can reflect the magnitude of the distortion. The distortion index is mainly used to measure the degree of distortion quantitatively. At present, different countries use different calculation methods, so there is not a universal definition of the distortion index. Because the Russian comprehensive distortion index is widely used in engine stability assessment and experimental research, the research work in this paper is based on the definition of the distortion index. The distortion index only considers the total pressure distortion, not the total temperature distortion, and mainly consists of three parts, namely, the circumferential total pressure unevenness, the radial total pressure unevenness, and the mean surface turbulence. However, a large number of statistical data show that the influence of radial total pressure unevenness on the stability of different types of engine structures is very small, so it is not considered. Therefore, the total pressure distortion index can be expressed as

$$W=\Delta \overline{{\sigma }_0}+{\epsilon }_{AV}$$

(40)

The composite distortion index $W$ is the sum of the steady-state circumferential total pressure distortion index $\Delta \overline{{\sigma }_0}$ and the surface mean turbulence degree ${\epsilon }_{AV}$. For the convenience of description, the $W$s in this paper are all percentages. The $\Delta \overline{{\sigma }_0}$ indicates the difference between the surface mean total pressure recovery coefficient and the total pressure recovery coefficient in the low-pressure area.

$$\Delta \overline{{\sigma }_0}=1-\frac{{\sigma }_0}{{\sigma }_{AV}}$$

(41)

The mean surface turbulence ${\epsilon }_{AV}$ represents the quantitative characteristics of the total pressure fluctuation at the aerodynamic interface.

$${\epsilon }_{AV}=\frac{{\sum }_{i=1}^n{\epsilon }_i}{n}$$

(42)

${\epsilon }_i$ is the ratio of the root-mean-square value of the pulsating pressure to the mean value of the total pressure $P^{*}$.

$$\begin{gathered} {\epsilon }_i=\Delta P^{*}/{\overline{P}}^{*} \\ \Delta P^{*}=\sqrt{\frac{1}{T}\underset{0}{\overset{T}{\int }}{\left(P^{*}\left(\tau \right)-{\overline{P}}^{*}\right)}^2} \\ {\overline{P}}^{*}=\frac{1}{T}\underset{0}{\overset{T}{\int }}{P}^{*}\left(\tau \right) d\tau \end{gathered}$$

(43)

In the formula above, T is the sampling time of the pulse airflow.

In Section 2.2.1, we present the calculation of inlet and outlet flow parameters, and in Formula (4), we refer to the total inlet pressure recovery factor ${\sigma }_i$, which is used to characterize the total pressure loss of the inlet flow. The general total pressure recovery coefficient is calculated by the following formula.

$${\sigma }_i=\left\{\begin{array}{ll}1.0-0.075{\left(M_a-1\right)}^{1.35}& ,M_a>1.0\\ 1.0& ,M_a\le 1.0\end{array}\right.$$

(44)

For the total pressure recovery coefficient with inlet distortion, it is correlated not only with $M_a$ but also with $W$. The total pressure recovery coefficient in this paper is:

$${\sigma }_i=\left\{\begin{array}{ll}\left(1-0.01W\right)\left[1-0.075{\left(M_a-1\right)}^{1.35}\right]& ,M_a>1.0\\ 1-0.01W& ,M_a\le 1.0 \end{array}\right.$$

(45)

The conventional fan component calculation module obtains the conversion flow and efficiency based on the interpolation of the fan relative conversion speed, and the calculation expression is

$${\sigma }_i=\left\{\begin{array}{ll}\left(1-0.01W\right)\left[1-0.075{\left(M_a-1\right)}^{1.35}\right]& ,M_a>1.0\\ 1-0.01W& ,M_a\le 1.0 \end{array}\right.$$

(46)

Similarly, in Section 2.2.2, some modifications are made to the fan characteristics by taking into account the inlet distortion index $W$ in the fan converted air flow and the overall efficiency calculations, i.e., Formula (6) is changed to:

$$\begin{gathered} m_{aL,cor}=m_{aL,cor}(n_{cL,cor},{\pi }_{cL}^{*},{\varphi }_L,W) \\ {\eta }_{cL}^{*}={\eta }_{cL}^{*}\left(n_{cL,cor},{\pi }_{cL}^{*},{\varphi }_L,W\right) \end{gathered}$$

(47)

So far, the calculation flow of the components related to the inlet distortion has been modified and combined with the existing engine component-level modeling methods. The nonlinear engine component-level model with the inlet distortion can be formed, and the schematic diagram of its main calculation flow is shown in Figure 7.

3.2. Adaptive Disturbance Rejection Control

The basic configuration of adaptive disturbance rejection control is a nonlinear dynamic inverse (NDI) control law. The control logic is shown in Figure 8.

Taking the tilt-turn dynamics model of UAV as an example, first, according to the dynamics of UAV established in Section 2.1, the instruction design rule is that the time scale of the state variable is separated, and the control instructions are generated by dynamic inversion according to the response speed of the variable.

3.2.1. Design of Slow Loop Control Instruction

The desired angular rate command $\dot{{\gamma }_t}$ and $\dot{{\vartheta }_t}$ can be obtained from the desired attitude angular response, and the expression is as follows:

$$\begin{gathered} \dot{{\gamma }_t}=K_{\gamma }\left({\gamma }_{cmd}-\gamma \right) \\ \dot{{\vartheta }_t}=K_{\vartheta }\left({\vartheta }_{cmd}-\vartheta \right) \end{gathered}$$

(48)

The angular rate instructions can be obtained by associating with the dynamic equations. The results are as follows:

$$\begin{gathered} {\omega }_{x,cmd}=K_{\gamma }\left({\gamma }_{cmd}-\gamma \right)+\tan \vartheta \left({\omega }_y\cos \gamma -{\omega }_z\sin \gamma \right) \\ {\omega }_{z,cmd}=\frac{1}{\cos \gamma }\left(K_{\vartheta }\left({\vartheta }_{cmd}-\vartheta \right)-{\omega }_y\sin \gamma \right) \end{gathered}$$

(49)

Observing the dynamic model, the angular acceleration model can be considered as the result of the interaction of the moment part and the non-moment part, and the non-moment control parts are $a_1,a_2,a_3$.

$$\begin{gathered} a_1=-\frac{\left(I_{ZZ}-I_{yy}\right)I_{yy}-I_{xy}^2}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_y{\omega }_z-\frac{\left(I_{xx}-I_{zz}+I_{yy}\right)I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_x{\omega }_z \\ {a}_2=-\frac{\left(I_{xx}-I_{zz}\right)I_{xx}-I_{xy}^2}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_x{\omega }_z+\frac{\left(I_{xx}-I_{zz}+I_{yy}\right)I_{xy}}{I_{xx}{I}_{yy}-I_{xy}^2}{\omega }_y{\omega }_z \\ {a}_1=-\frac{I_{yy}-I_{xx}}{I_{zz}}{\omega }_x{\omega }_y+\frac{I_{xy}}{I_{zz}}\left({\omega }_x^2-{\omega }_y^2\right) \end{gathered}$$

(50)

The above equation can be rewritten as

$$\left[\begin{array}{c}\dot{{\omega }_x}\\ \dot{{\omega }_y}\\ \dot{{\omega }_z}\end{array}\right]={\left[\begin{array}{ccc}-I_{xy}& I_{yy}& 0\\ I_{xx}& -I_{xy}& 0\\ 0& 0& I_{zz}\end{array}\right]}^{-1}\left(\left[\begin{array}{c}l\\ m\\ n\end{array}\right]+\left[\begin{array}{c}{l}_{\delta }\\ m_{\delta }\\ n_{\delta }\end{array}\right]\right)+\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right]$$

(51)

The control quantities are $l_{\delta },m_{\delta },n_{\delta }$, and the control instruction is designed according to the above expected response method.

$$\left[\begin{array}{c}\dot{{\omega }_{x,t}}\\ \dot{{\omega }_{y,t}}\\ \dot{{\omega }_{z,t}}\end{array}\right]=\left[\begin{array}{ccc}{K}_{{\omega }_x}& 0& 0\\ 0& K_{{\omega }_y}& 0\\ 0& 0& K_{{\omega }_z}\end{array}\right]\left[\begin{array}{c}{\omega }_{x,cmd}-{\omega }_x\\ {\omega }_{y,cmd}-{\omega }_y\\ {\omega }_{z,cmd}-{\omega }_z\end{array}\right]$$

(52)

which is associated with the kinetic equations:

$$\left[\begin{array}{c}{l}_{\delta ,cmd}\\ m_{\delta ,cmd}\\ n_{\delta ,cmd}\end{array}\right]={\left[\begin{array}{ccc}-I_{xy}& I_{yy}& 0\\ I_{xx}& -I_{xy}& 0\\ 0& 0& I_{zz}\end{array}\right]}^{-1}\left(\left[\begin{array}{ccc}{K}_{{\omega }_x}& 0& 0\\ 0& K_{{\omega }_y}& 0\\ 0& 0& K_{{\omega }_z}\end{array}\right]\left[\begin{array}{c}{\omega }_{x,cmd}-{\omega }_x\\ {\omega }_{y,cmd}-{\omega }_y\\ {\omega }_{z,cmd}-{\omega }_z\end{array}\right]-\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right]\right)-\left[\begin{array}{c}l\\ m\\ n\end{array}\right]$$

(53)

3.2.2. Design of Adaptive Disturbance Rejection Control

In the above derivation process, the body torque of the UAV is strongly nonlinear during the flight process, and it cannot accurately rely on the existing data to calculate it during the flight process. Similarly, when the moment control command is calculated, the rudder surface efficiency cannot be obtained accurately, and the exact value of rudder surface deflection cannot be determined to realize the control moment; the original control law is improved to overcome its dependence on an accurate model. The adaptive disturbance rejection control operates by comparing the measured response with an internal model of the desired dynamics. It then adjusts the control signal based on the difference between the desired and measured behaviors. First, the desired angular rate is designed as follows:

$$\dot{{\omega }_{des}}=K_{\omega }{\omega }_{cmd}-{\omega }_{des}$$

(54)

Because of the simplified assumptions, model uncertainty, external disturbance, and so on, dynamic inversion cannot achieve its goal completely. Therefore, $d$ is introduced to describe all the uncertainty errors of the system. The system dynamics are written in a form similar to the expected behavior as follows.

$$\dot{\omega }=K_{\omega }({\omega }_{cmd}-{\omega }_{des})+u_{ad}+d$$

(55)

Among them, $u_{ad}$ brings the angular acceleration increment for the adaptive control moment, $u_{ad}=I^{-1}\Delta M$. The angular rate internal model instruction is designed as.

$$\dot{\widehat{\omega }}=K_{\omega }({\omega }_{cmd}-{\omega }_{des})+u_{ad}-K_d\left( \tilde{\omega }-\omega \right)$$

(56)

where $\widehat{\omega }$ is an output of the UAV state space model. $\tilde{\omega }=\widehat{\omega }-\omega$ is defined. The difference between the true angular rate and the internal model angular rate is obtained.

$$\dot{\tilde{\omega }}=\dot{\widehat{\omega }}-\dot{\omega }=-K_{ad}\tilde{\omega }-d$$

(57)

According to the solution of the first-order non-homogeneous differential equation, the functional relationship of $\tilde{\omega }$ with time can be solved.

$$\tilde{\omega }=-\frac{d}{K_{ad}}+\left(\tilde{\omega }\left(0\right)+\frac{d}{K_{ad}}\right)e^{-K_{ad}t}$$

(58)

From the above, we know that, in a certain period of time, $\tilde{\omega }$ converges to $-d/K_{ad}$, and in the same way, in this period of time, $\dot{\tilde{\omega }}$ converges to 0, that is, $-K_{ad}\tilde{\omega }-d$ converges to 0. Based on the actual $\omega$ differential equation, the $u_{ad}=K_{ad}\tilde{\omega }$ can be substituted into Formula (61) to obtain:

$$\dot{\omega }=K_{\omega }({\omega }_{cmd}-{\omega }_{des})$$

(59)

According to the solution of the first-order non-homogeneous linear differential equation,

$$\omega ={\omega }_{cmd}+\left(\omega \left(0\right)-{\omega }_{cmd}\right)e^{-K_{ad}t}$$

(60)

In this case, $\omega$ converges to ${\omega }_{cmd}$ within a certain time, so $u_{ad}=K_{ad}\tilde{\omega }$ can suppress the effect of error $d$.

According to the dynamic inversion of the above dynamic equations, the adaptive input is added to the NDI control torque command to generate the total torque command.

$$M=I{K}_{\omega }({\omega }_{cmd}-{\omega }_{des})+I{K}_{ad}\left(\widehat{\omega }-\omega \right)+I{\left[\begin{array}{ccc}{a}_1& a_2& a_3\end{array}\right]}^T$$

(61)

4. Simulation Results

In the preceding paragraph, we build the dynamics model of UAV and the engine model of UAV and add the influence of inlet distortion into the engine model. Then, the strategy of thrust optimization based on the adaptive margin model and the method of adaptive disturbance rejection flight control for UAV are proposed. In order to verify this, the flight control strategy and the thrust optimization strategy proposed earlier are simulated and verified successively. The parameters of a certain type of UAV are given in

Table 2. Parameters of a certain type of UAV.

Parameter	Symbol	Value
aircraft mass ($\mathrm{kg}$)	$m$	9295.44
wingspan ($\mathrm{m}$)	$b$	9.144
wing area (${\mathrm{m}}^2$)	$S$	27.87
mean aerodynamic chord ($\mathrm{m}$)	$\overline{c}$	3.45
roll moment of inertia ($\mathrm{kg}·{\mathrm{m}}^2$)	$I_{xx}$	12,874.8
yaw moment of inertia ($\mathrm{kg}·{\mathrm{m}}^2$)	$I_{yy}$	85,552.1
pitch moment of inertia ($\mathrm{kg}·{\mathrm{m}}^2$)	$I_{zz}$	75,673.6
product moment of inertia ($\mathrm{kg}·{\mathrm{m}}^2$)	$I_{xy}$	1331.4
product moment of inertia ($\mathrm{kg}·{\mathrm{m}}^2$)	$I_{xz}$	0
product moment of inertia ($\mathrm{kg}·{\mathrm{m}}^2$)	$I_{yz}$	0
c.g. location ($\mathrm{m}$)	$x_{cg}$	0.3$\overline{c}$
reference c.g. location ($\mathrm{m}$)	$x_{cgr}$	0.35$\overline{c}$

. The simulation is based on the model and strategy established in the previous paper, which are implemented in code in MATLAB, where the turbine engine model is built by Simulink of MATLAB.

4.1. Simulation Results of Thrust Optimization Based on the Adaptive Margin Model

In order to verify the optimization results of the engine thrust performance, we use the acceleration performance improvement of the UAV at the altitude of 7500 m as a way to test the engine thrust optimization, to increase the speed of the UAV from 150 m/s to 250 m/s with the altitude change, and to check the time of the engine acceleration plan as well as the thrust output of the engine. The simulation results are as follows:

Figure 9a shows the speed simulation results of the UAV for a given speed command, and we compare the speed results after the thrust performance optimization with the original case. It can be seen that the acceleration performance of the UAV after the thrust lift is obviously improved compared with the original situation. In the original acceleration plan, the UAV used 28.9 s to accelerate the speed to 249 m/s, the speed was increased to 250 m/s in 40 s, the time of acceleration to 249 m/s was advanced to 25.5 s, and the time of acceleration to 250 m/s was reduced to 29.9 s after the thrust optimization. This shows that, under the same engine acceleration control parameters, the engine thrust optimization strategy proposed in this paper can significantly improve the acceleration capability of UAV and shorten the acceleration time. It is also possible to see this by the slope of the velocity resultant curve in plot 9(a), i.e., the maximum acceleration of the UAV.

Figure 9b shows the result of UAV engine throttle command. Under the original condition, the engine throttle command needs to keep 19.3 s fully open, but after improvement, the change time is reduced to 17.5 s. The result shows that the full throttle time of UAV is shorter after the thrust performance is improved. In the acceleration control program of the UAV, the controller parameters used in the two simulated scenarios are exactly the same. Therefore, the optimized throttle command maximum output time is shorter than the original one, which can also indicate that the optimized case has an increased maximum thrust output compared to the original case.

To visualize the thrust improvement, Figure 10 shows the UAV thrust output during this acceleration; the percentage increase in the thrust of the engine at 100% throttle is in the corner of the image (take the first 17.5 s, which is 100% of the time for the engine throttle command in both cases). The right axis of the small graph is the label, and the unit is the percentage of thrust performance improvement. It can be seen that the maximum thrust output of the engine can be increased by at least 8% under the condition of maximum throttle command and also increases with the increase in the UAV speed, reaching 10.4% at 17.5 s.

Figure 11 is the result diagram of the working state of the compressor. It can be seen that the working point of the engine is closer to the surge boundary after the more reasonable order of the turbocharging ratio, that is, the surge margin reserve of the engine is reduced but still kept in safe working conditions, and the engine did not surge.

Brief summary:

• Under the same acceleration control parameters, the thrust optimization strategy proposed in this paper can significantly improve the acceleration capability of UAV and shorten the acceleration time.
• With improved thrust performance, the UAV requires a shorter full throttle time.
• At the same flat flying speed, the throttle command of the engine after thrust optimization is smaller than the original throttle command.
• At maximum throttle, the engine’s maximum thrust output can be increased by at least 8%.

4.2. Simulation Results of Adaptive Disturbance Rejection Control

In order to verify the designed adaptive disturbance rejection control strategy, the UAV is given different flight instructions in the longitudinal and rolling directions, and the advantages and disadvantages of the control method and the traditional PID control were observed. The flight environment is given under two conditions: altitude 7500 m, speed 150 m/s.

The first is the longitudinal simulation results, the initial pitch angle is $\gamma =$ 7.2°, the initial track angle $\alpha$ is 0°, the UAV is controlled in the first 5 s of the flat flight, the UAV is given a 15° pitch command at 5 s, and the UAV is made to climb quickly. Observe the response of the pitch angle and the command output of the rudder deviation.

For the PID control case in the simulation, we used the Ziegler–Nichols tuning method. This method provides a good starting point for the PID gains, which can then be refined through further tuning. Start by adjusting the proportional gain Kp until the system’s response is close to the desired response. Then, adjust the integral gain Ki until the steady-state error is reduced to an acceptable level. Finally, adjust the derivative gain Kd to reduce the overshoot and dampen the oscillations caused by the proportional and integral gains. A set of PID controller parameters with the fastest response time and the smallest overshoot is obtained by the above method.

As shown in Figure 12c, it took 17.39 s for the PID controller to make the pitch angle converge to the command, while the adaptive disturbance rejection control took 3.62 s to achieve the equivalent situation (0.05 degrees of static difference). Compared with the PID control, the pitch angle of the adaptive disturbance rejection control method converges to the given instruction more rapidly, the overshoot is smaller, and the overshoot arrival time is shorter.

Similarly, the simulation results for the roll direction are given in Figure 13. The initial tilt angle is 0°, and the control UAV will fly horizontally in the first 5 s. The UAV is given a 35° tilt command at 5 s to make the UAV roll quickly. Observe the response of the tilt angle and the command output of the rudder.

As shown in Figure 13a, The PID controller produces a large overshoot in the tilt angle (12.82% overshoot), which is obviously undesirable in this case, while the adaptive disturbance rejection control method has only 0.4% overshoot. Compared to the result of PID control, the tilt angle of the adaptive disturbance rejection control mode converges to the given instruction more quickly, the overshoot is smaller, and the overshoot arrival time is shorter. The aileron command, compared with the PID control in the case of swing, greatly reduced. This increases the potential of the UAV to perform other maneuvers. The results of both the longitudinal and rolling flight simulation show that the adaptive disturbance rejection control method proposed in this paper can achieve fast attitude maneuvers; compared with PID control, it can complete the maneuver command more accurately in a short time and improve the overshoot and rudder surface jitter caused by PID control.

5. Conclusions

We designed a strategy for exchanging the surge margin of the engine for the improvement of the thrust performance in order to meet the requirement of the development of the propulsion system for UAVs. First, the dynamic equation of UAV is established, which describes the straight motion of the center of mass and the rotation motion of the body around the center of mass. A mathematical model of the turbofan’s components is built. Then, an optimal thrust algorithm based on the adaptive margin model is designed. The internal working environment of the engine and the angle of attack in flight dynamics are considered. From the angle of aircraft/engine integration, the way to improve the performance of the engine propulsion system is studied; in this part, the UAV body transfer function is obtained from the UAV dynamic equation, and the UAV angle of attack is predicted by the transfer function. Then, the inlet and fan models are modified by adding a correction term considering inlet distortion. For further flight simulation, a flight control strategy based on adaptive disturbance rejection control is designed. In order to compensate for the external loop control law, a disturbance observer is introduced to estimate the flight disturbance of the UAV.

Finally, the feasibility of the proposed strategy is verified by numerical simulation. The numerical simulation results show that the thrust lifting strategy based on the adaptive margin model can ensure the safety of the flight and the stability of the compressor; the maximum thrust of the UAV is increased by at least 8% by reducing the surge margin reserve and improving the maneuverability of the UAV. The adaptive disturbance rejection control method can suppress the attitude fluctuation quickly and has strong robustness and adaptive ability.

Nowadays, there are more and more research works on aircraft/engine integration, and in this paper, we meet the demand of the thrust performance improvement of UAVs from this perspective, which has some implementability in engineering practice. However, the strategy still has some limitations, namely, a high accuracy requirement for the engine model. In fact, turbine engine power system modeling is a complex discipline, and the engine model has strong nonlinearity. The component-level modeling method used in this paper, at the cost of a lengthy and complex mechanism model and in exchange for a higher accuracy, sacrifices a certain amount of computing speed. We plan to continue to explore the thrust optimization method based on the adaptive margin model and, on the other hand, to seek a more accurate and faster engine modeling method.