Incremental MPC for morphing aircraft: A robust and computationally efficient optimal control approach

Abstract

Compared with conventional aircrafts, the dynamics of morphing aircraft exhibit significant nonlinearity and high uncertainty posing substantial challenges to the controller design. This paper proposed an incremental model predictive control (IMPC) approach for a morphing aircraft with the varying wingspan and sweep angle. The proposed method, formulated as the constrained optimal control problems (OCPs), is robust to dynamic uncertainties and capable of handling input and state constraints such as mechanical limits and angle of attack. Firstly, we establish a nonlinear affine model of the morphing aircraft according to the aircraft’s dynamic model. Subsequently, a time-scale separation (TSS) method is introduced according to the aircraft’s characteristics, dividing the system into fast and slow subsystems. Each subsystem is approximated as an incremental system, reducing dependency on the nominal model. In addition, the constrained OCP is cast into a quadratic programing (QP) problem using the incremental model, resulting in a linear MPC. Finally, several simulations are performed to validate the effectiveness of the IMPC method, in terms of stability, robustness, and computational efficiency.

Introduction

Background

Most conventional aircrafts are designed for specific flight conditions, leading to suboptimal aerodynamic performance across different flight scenarios and potentially compromising certain mission phases. The morphing aircraft are new concept aircraft that actively switch their shapes in response to dynamic flight environments or specific mission requirements.^1,2 They have attracted significant attention due to their potential to achieve optimal flight performance in diverse mission environments, including take-off, hover, reconnaissance, attack, and landing. Innovations such as retractable and folded wings, morphing swept wings, more complex variable geometry wings, and actively hinged wings are increasingly becoming the focal points of cutting-edge research in aeronautical engineering.

The flight performance is optimized through dynamically adjusting the aircraft’s aerodynamic configuration. The deformation, however, introduces significant variations in dynamics, such as changes in the center of gravity, variations in rotational inertia, and modifications to the wingspan and surface area are potential to cause instability.^3,4 Consequently, the morphing process introduces strong disturbances to the system dynamics and poses challenges to the controller design.

Related works

Significant advancements in controllers have been achieved for morphing aircraft. Linear control theory has been established in control engineering practice for a long time and has already shown certain advantages. Song et al.^5,6 introduced the PID control to the morphing aircraft system for the first time and integrated the Lagrange equation, pseudo-inverse control and classical PID methods to solve the deformation problems of sliding skin wings. Yan et al. ⁷ developed a parameter-tunable PID controller for different configurations of the morphing aircraft. Soneda et al. ⁸ proposed a classical linear quadratic regulator (LQR) controller to regulate the pitch motion of the morphing aircraft. Neal et al. ⁹ employed LQR control to determine the optimal aerodynamic forces required for following any flight path, theoretically proving that LQR and fuzzy logic achieve satisfactory control. Linear parameter-varying (LPV) control theory, which effectively captures the impact of parameter variations in system models, has been widely adopted in the design of attitude control systems for morphing aircraft. ¹⁰ In practical applications, the LPV approach necessitates the trim and small-disturbance linearization of nonlinear systems, which causes modeling inaccuracies.

Due to the nonlinear characteristics and strong coupling of the morphing aircraft, traditional linear control methods require small perturbation linearization and decoupling analysis of the nonlinear motion model. However, these methods exhibit significant degradation in control performance with increased deformation complexity or disturbance levels. In recent years, nonlinear control methods with modeling advantages have received significant attention from scholars, and become the primary development direction of control system design for the morphing aircraft. The nonlinear dynamic inversion (NDI) method designs control signals for the system utilizing dynamic inversion transformation and the Lyapunov theory. ¹¹ And the Backstepping methods^12,13 decompose complex high-order nonlinear systems into first-order subsystems and design virtual control signals for each subsystem. The traditional NDI and Backstepping design typically depend on a precise mathematical model of the system. Considering the high uncertainty in the morphing aircraft model, robust control techniques such as adaptive control often need to be combined with the NDI or Backstepping methods.^14,15 The sliding mode control (SMC), because of its rapid response and robustness, ¹⁶ has been widely applied in the control of morphing aircraft. For instance, Yue et al. ¹⁷ developed a sliding mode controller to address modeling challenges caused by the deformations in the wingspan. Nonetheless, the discontinuous mode switching in SMC induces chattering, which needs to be combined with filters or adaptive methods. ¹⁸

Although the adaptive technique and SMC provide robust solutions to a certain degree, they fail to take into account input and state constraints. Model predictive control (MPC), an optimization-based control system, stands out for its ability to systematically address inputs and state constraints.^19,20 The primary advantage of MPC lies in its capacity to handle these constraints in a systematic and non-conservative manner while achieving optimal control performance. This unique capability has led to its rising popularity in the field of morphing aircraft control. Dai et al. ²¹ introduced a morphing waverider MPC strategy using the barrier Lyapunov function. To enhance the robustness, the interference observer of the deformation torque is employed and the deformation torque is estimated online. Pereira et al. ²² explored two distinct MPC frameworks, conducting numerical simulations on the nonlinear model of flexible-wing aircraft and investigating the predictive model’s accuracy in the MPC design process. Nonetheless, uncertainties in models and parameters degrade the prediction accuracy, thereby impacting control precision badly. Besides, the computational demands of traditional nonlinear model predictive control (NMPC) are considerable, posing challenges for real-time control.

Method and contributions

In this paper, we developed an incremental model predictive control (IMPC) approach for morphing aircraft. First, according to the principle of timescale-scale separation (TSS) theory, a two-layer controller is designed. Then, to improve robustness against model uncertainties, we develop a robust incremental model. Next, for each layer, an incremental model predictive controller is developed through constructing an optimal control problem (OCP) with inputs and state constraints, where state predictions are produced using the discretized incremental model. Each constrained OCP of IMPC is cast into a quadratic programing (QP). The contributions of this article are summarized as follows:

• (1) Enhanced Robustness: Distinct from traditional MPC methods, IMPC does not require a concrete dynamic model of the morphing aircraft. This is because the method approximates the continuous nonlinear system as an incremental form and the state predictions are produced using the linear discretized incremental system.^23,24

• (2) Optimal Control Performance: The discrete incremental controller is designed using MPC framework without terminal components, where inputs and state constraints are expressed as inequality constraints in the OCP. Consequently, optimal control performance is attained while considering both inputs and state constraints.

• (3) High Computational Efficiency: Relying on the incremental method not only enhances the robustness of the controller but also converts complex nonlinear equations into linear ones. The development IMPC is essentially a global linear MPC and each constrained OCP is cast into a QP. Compared with traditional nonlinear MPC, ²⁵ the proposed IMPC reduces computational complexity significantly, which permits real-time control in milliseconds.

Organization

The remainder of this paper is organized as follows. First, the morphing aircraft model and control objective are introduced in Section 2. Next, the control structure is determined according to TSS, and the constrained IMPC is presented in Section 3. In addition, Section 4 provides a theoretical analysis of the recursive feasibility of IMPC. Finally, a series of simulations on a morphing aircraft verify the effectiveness of the proposed IMPC in Section 5, with conclusions presented in Section 6.

Problem formulation

Modeling

In this study, we consider the morphing aircraft Navion L-17 ²⁶ with the wingspan and wing sweep angle which both vary symmetrically. The change range of the wingspan is defined from the original wingspan to twice the wingspan, and the change range of the wing sweep angle is 0°–40°. Figures 1 and 2 illustrate the changing process of the wingspan and sweep angle of the aircraft respectively.

Figure 1.

Navion L-17 morphing wingspan.

Figure 2.

Navion L-17 morphing sweep angle.

The wingspan deformation rate ${\xi }_b$ and the sweep angle deformation rate ${\xi }_s$ are defined respectively as follows:

$${\xi }_b=\frac{b-b_{min}}{b_{max}-b_{min}},{\xi }_s=\frac{s-s_{min}}{s_{max}-s_{min}}$$

where $b$ is the wingspan, $b_{min}$ and $b_{max}$ are the wingspan under deformation free and the wingspan under maximum format respectively. $s$ is the sweep angle, $s_{min}$ and $s_{max}$ are the aircraft’s sweep angle without deformation and the sweep angle at maximum deformation respectively.

Here we only consider the longitudinal motion of the aircraft during the wing deformation process. The longitudinal dynamic equations in the deformation process are expressed as follows^10,27,28:

$$\{\begin{array}{c}\stackrel{·}{V}=\frac{1}{m}T\cos \alpha -\frac{1}{m}D-g\sin (\theta -\alpha )\hfill \\ \stackrel{·}{\alpha }=-\frac{1}{\mathrm{mV}}T\sin \alpha -\frac{1}{\mathrm{mV}}L+q+\frac{1}{V}g\cos (\theta -\alpha )\hfill \\ \stackrel{·}{\theta }=q\hfill \\ \stackrel{·}{q}=\frac{1}{I_y}M\hfill \\ \stackrel{·}{h}=V\sin (\theta -\alpha )\hfill \end{array}$$

(1)

where $V$ is the velocity, $\alpha$ is the angle of attack, $\theta$ is the angle of pitch, $q$ is the pitch rate. $h$ is the altitude, $g$ is the gravity acceleration. $m$ and $I_y$ illustrate the mass of the aircraft and the moment of inertia about the pitch axis. The thrust force $T$ , drag force $D$ , lift force $L$ and pitch moment $M$ in (1) are expressed as follows ²⁹ :

$$\{\begin{array}{c}T=T_{\delta }{\delta }_t\hfill \\ L=0.5\rho V^2{S}_{\omega }{C}_L(\alpha ,{\delta }_e,{\xi }_b,{\xi }_s)\hfill \\ D=0.5\rho V^2{S}_{\omega }{C}_D(\alpha ,{\xi }_b,{\xi }_s)\hfill \\ M=0.5\rho V^2{S}_{\omega }{c}_A{C}_m(\alpha ,{\delta }_e,{\xi }_b,{\xi }_s)\hfill \end{array}$$

(2)

where $T_{\delta }$ is the thrust coefficient of the engine, ${\delta }_t$ denotes the throttle opening, ${\delta }_e$ denotes the elevator deflection, $S_{\omega }$ is the wing reference area, $c_A$ is the average aerodynamic chord length, $\rho$ is the atmospheric density. $C_L$ is the lift coefficient, $C_D$ is the resistance coefficient, $C_m$ is the pitch moment coefficient.

Then the variable relations between aerodynamic parameters and deformation rates ${\xi }_b$ and ${\xi }_s$ are determined. Lift coefficient $C_L$ , drag coefficient $C_D$ , and pitching moment coefficient $C_M$ are defined as follows:

$$\{\begin{array}{c}{C}_L=C_{L_{\alpha =0}}+C_{L_{\alpha }}\alpha +C_{L_{{\delta }_e}}{\delta }_e\hfill \\ C_D=C_{D_{\alpha =0}}+C_{D_{\alpha }}\alpha +C_{D_{{\alpha }^2}}{\alpha }^2\hfill \\ C_M=C_{M_{\alpha =0}}+C_{M_{\alpha }}\alpha +C_{M_{{\delta }_e}}{\delta }_e\hfill \end{array}$$

(3)

The coefficient values are shown as follows, with $\mathbf{C}=(c_{\mathrm{ij}}{)}_{7\times 3}$ in Section 5.

$$\begin{array}{c}{\left[\begin{array}{ccccccc}{C}_{L_{\alpha =0}}\hfill & C_{L_{\alpha }}\hfill & C_{D_{\alpha =0}}\hfill & C_{D_{\alpha }}\hfill & C_{D_{{\alpha }^2}}\hfill & C_{M_{\alpha =0}}\hfill & C_{M_{\alpha }}\hfill \end{array}\right]}^{\mathbf{T}}\hfill \\ =\mathbf{C}{\left[\begin{array}{ccc}1\hfill & {\xi }_b\hfill & {\xi }_s\hfill \end{array}\right]}^{\mathbf{T}}\hfill \end{array}$$

(4)

Finally, the longitudinal aerodynamic model of the morphing wingspan and morphing sweep angle aircraft is reformulated into an affine nonlinear model as follows:

$$\stackrel{·}{\mathbf{x}}=\mathbf{f}(\mathbf{x})+\mathbf{g}(\mathbf{x})\mathbf{u}+\mathbf{h}(\mathbf{x})\mathbf{d}$$

(5)

where $\mathbf{x}={\left[\begin{array}{ccccc}V\hfill & \alpha \hfill & \theta \hfill & q\hfill & h\hfill \end{array}\right]}^{\mathbf{T}}$ , $\mathbf{u}={\left[\begin{array}{cc}{\delta }_e\hfill & {\delta }_t\hfill \end{array}\right]}^{\mathbf{T}}$ , $\mathbf{d}={\left[\begin{array}{cc}{\xi }_b\hfill & {\xi }_s\hfill \end{array}\right]}^{\mathbf{T}}$ .

Besides, the functions $\mathbf{f}(\mathbf{x})$ , $\mathbf{g}(\mathbf{x})$ , and $\mathbf{h}(\mathbf{x})$ are expressed as follows:

$$\mathbf{f}(\mathbf{x})=\left[\begin{array}{c}-\frac{1}{2m}\rho V^2{S}_{\omega }(c_{31}+c_{41}\alpha +c_{51}{\alpha }^2)-g\sin (\theta -\alpha )\hfill \\ -\frac{1}{2m}\rho V{S}_{\omega }(c_{11}+c_{12}\alpha )+q+\frac{1}{V}g\cos (\theta -\alpha )\hfill \\ q\hfill \\ \frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A(c_{61}+c_{71}\alpha )\hfill \\ V\sin (\theta -\alpha )\hfill \end{array}\right]$$

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{g}(\mathbf{x})=\left[\begin{array}{cc}0\hfill & \frac{1}{m}{T}_{{\delta }_t}\cos \alpha \hfill \\ \multicolumn{1}{c}{-\frac{1}{2m}\rho V{S}_{\omega }{C}_{L_{{\delta }_e}}}& -\frac{1}{\mathrm{mV}}{T}_{{\delta }_t}\sin \alpha \hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{\frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A{C}_{M_{{\delta }_e}}}& 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]}\end{array}$$

$$\begin{array}{c}\mathbf{h}(\mathbf{x})=\hfill \\ \left[\begin{array}{cc}-\frac{1}{2m}\rho V^2{S}_{\omega }(c_{32}+c_{42}\alpha +c_{52}{\alpha }^2)\hfill & -\frac{1}{2m}\rho V^2{S}_{\omega }(c_{33}+c_{43}\alpha +c_{53}{\alpha }^2)\hfill \\ \multicolumn{1}{c}{-\frac{1}{2m}\rho V{S}_{\omega }(c_{12}+c_{22}\alpha )}& -\frac{1}{2m}\rho V{S}_{\omega }(c_{13}+c_{23}\alpha )\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{\frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A(c_{62}+c_{72}\alpha )}& \frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A(c_{63}+c_{73}\alpha )\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\hfill \end{array}$$

Note that in this study, the morphing process is considered as a disturbance to the dynamics. In addition, the aerodynamic modeling error is inevitable. Thus, a robust controller is necessary.

Control objective

The objective of this study is to develop a robust controller for the aircraft with morphing wingspan and morphing sweep capabilities. Considering both state and input constraints, the controller enables the morphing aircraft to track predefined velocity and altitude while resisting model uncertainties. Simultaneously, to ensure in-flight safety, it is imperative to maintain the angle of attack within an acceptable range. Additionally, the aircraft’s actuators are restricted by the design of the power system and energy supply, thereby constraining the input to a defined range.

We aim to design a robust optimal controller which is robust while meeting the safety and physical constraints of the morphing aircraft, such as constraints on flight angle of attack, elevator angle, and throttle opening. The following box constraints are considered:

$${\alpha }_{min}\le \alpha \le {\alpha }_{max},{\delta }_{e_{min}}\le {\delta }_e\le {\delta }_{e_{max}},{\delta }_{t_{min}}\le {\delta }_t\le {\delta }_{t_{max}}$$

• Remark 1. (Aircraft Safety Constraints): The angle of attack is defined as the angle between the wing chord line and the relative airflow. Generally, an increase in the angle of attack results in an increase in lift. However, each wing features a critical angle of attack beyond which airflow separation occurs, leading to a stall. Maintaining the angle of attack below this critical threshold is crucial for stall prevention.

Controller design

Control structure

In the morphing aircraft system, the time scales of each state are quite different. According to Kawaguchi et al.,^30,31 state variables are divided into four levels according to the time scale, as shown in Table 1.

Table 1.

Time-scale separation of state variables.

Classes	System state
Layer 1 (ultrafast variable)	$q$
Layer 2 (fast variable)	$\alpha$ , $\theta$
Layer 3 (slow variable)	$V$
Layer 4 (ultraslow variable)	$h$

Considering that the position response is slow compared with the attitude response, the whole control system of aircraft is divided into a fast sub-system (System I) and a slow sub-system (System II).^32,33 The ${\theta }_d$ generated by the upper controller is tracked by controlling $\theta$ , thereby achieving altitude $h$ tracking. Consequently, a two-layer control structure is established, as illustrated in Figure 3.

Figure 3.

Control structure with separation method.

In system I, with $V_d$ and ${\theta }_d$ serving as the tracking reference, the state vector for system I is ${\mathbf{x}}_1={\left[\begin{array}{cccc}V\hfill & \alpha \hfill & \theta \hfill & q\hfill \end{array}\right]}^{\mathbf{T}}$ , where the input vector is ${\mathbf{u}}_1={\left[\begin{array}{cc}{\delta }_e\hfill & {\delta }_t\hfill \end{array}\right]}^{\mathbf{T}}$ , and the deformation input vector is ${\mathbf{d}}_1={\left[\begin{array}{cc}{\xi }_b\hfill & {\xi }_s\hfill \end{array}\right]}^{\mathbf{T}}$ , the system dynamics is expressed as follows:

$${\stackrel{·}{\mathbf{x}}}_1={\mathbf{f}}_1({\mathbf{x}}_1)+{\mathbf{g}}_1({\mathbf{x}}_1){\mathbf{u}}_1+{\mathbf{h}}_1({\mathbf{x}}_1){\mathbf{d}}_1$$

(6)

With

$${\mathbf{f}}_1({\mathbf{x}}_1)=\left[\begin{array}{c}-\frac{1}{2m}\rho V^2{S}_{\omega }(c_{31}+c_{41}\alpha +c_{51}{\alpha }^2)-g\sin (\theta -\alpha )\hfill \\ -\frac{1}{2m}\rho V{S}_{\omega }(c_{11}+c_{12}\alpha )+q+\frac{1}{V}g\cos (\theta -\alpha )\hfill \\ q\hfill \\ \frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A(c_{61}+c_{71}\alpha )\hfill \end{array}\right]$$

(7)

And

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{g}}_1({\mathbf{x}}_1)=\left[\begin{array}{cc}0\hfill & \frac{1}{m}{T}_{{\delta }_t}\cos \alpha \hfill \\ \multicolumn{1}{c}{-\frac{1}{2m}\rho V{S}_{\omega }{C}_{L{\delta }_e}}& -\frac{1}{\mathrm{mV}}{T}_{{\delta }_t}\sin \alpha \hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{\frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A{C}_{m{\delta }_e}}& 0\hfill \end{array}\right]}\end{array}$$

(8)

$$\begin{array}{c}{\mathbf{h}}_1({\mathbf{x}}_1)=\hfill \\ \left[\begin{array}{cc}-\frac{1}{2m}\rho V^2{S}_{\omega }(c_{32}+c_{42}\alpha +c_{52}{\alpha }^2)\hfill & -\frac{1}{2m}\rho V^2{S}_{\omega }(c_{33}+c_{43}\alpha +c_{53}{\alpha }^2)\hfill \\ \multicolumn{1}{c}{-\frac{1}{2m}\rho V{S}_{\omega }(c_{12}+c_{22}\alpha )}& -\frac{1}{2m}\rho V{S}_{\omega }(c_{13}+c_{23}\alpha )\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{\frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A(c_{62}+c_{72}\alpha )}& \frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A(c_{63}+c_{73}\alpha )\hfill \end{array}\right]\hfill \end{array}$$

(9)

For system II, the dynamic equation is expressed as follows:

$${\stackrel{·}{x}}_2=V\sin (u_2-\alpha )$$

(10)

where, $x_2=h$ represents the system II state variable. $u_2$ is the virtual input of system II, and the reference value $h_d$ is tracked by controlling ${\theta }_d$ , that is, $u_2={\theta }_d$ .

• Remark 2. (Time-scale Separation): Time scale separation is according to the cause and effect of each state variable in the system.^34,35 Thus, the state variables influenced by the control variables can be considered as “fast” states, and the state variables influenced by fast state variables can be treated as “slow” variables. In the aircraft model, the actuator inputs are often treated as the fastest, the body attitude states are considered slower than the actuator inputs, but can be considered faster than the position dynamics. Mathematically, a small parameter $\epsilon$ represents the ratio of the fast time scale to the slow time scale. By considering the limit as $\epsilon \to 0$ , we derive reduced models that capture the essential behavior of the system without the computational complexity of simultaneously simulating all time scales.

Incremental system design

Before the controller design, the uncertain model dynamics are approximated to enhance the robustness by the incremental system. We take system I for example. Consider the following nonlinear affine model:

$${\stackrel{·}{\mathbf{x}}}_1={\mathbf{f}}_1({\mathbf{x}}_1)+{\mathbf{g}}_1({\mathbf{x}}_1){\mathbf{u}}_1+{\mathbf{h}}_1({\mathbf{x}}_1){\mathbf{d}}_1$$

(11)

Under the current operating point $({\mathbf{x}}_{1,0},{\mathbf{u}}_{1,0})$ with a near neighbor, the formula (11) is rewritten as a first-order Taylor series expansion:

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{\mathbf{x}}}_1={\stackrel{·}{\mathbf{x}}}_{1,0}+\frac{\partial }{\partial {\mathbf{x}}_1}({\mathbf{f}}_1({\mathbf{x}}_1)+{\mathbf{g}}_1({\mathbf{x}}_1){\mathbf{u}}_1+{\mathbf{h}}_1({\mathbf{x}}_1){\mathbf{d}}_1){|}_{\genfrac{}{}{0ex}{}{{\mathbf{x}}_1={\mathbf{x}}_{1,0}}{{\mathbf{u}}_1={\mathbf{u}}_{1,0}}}\mathrm{\Delta }{\mathbf{x}}_1}\\ +\frac{\partial }{\partial {\mathbf{u}}_1}({\mathbf{f}}_1({\mathbf{x}}_1)+{\mathbf{g}}_1({\mathbf{x}}_1){\mathbf{u}}_1+{\mathbf{h}}_1({\mathbf{x}}_1){\mathbf{d}}_1){|}_{\genfrac{}{}{0ex}{}{{\mathbf{x}}_1={\mathbf{x}}_{1,0}}{{\mathbf{u}}_1={\mathbf{u}}_{1,0}}}\mathrm{\Delta }{\mathbf{u}}_1+{\epsilon }_1\hfill \\ ={\stackrel{·}{\mathbf{x}}}_{1,0}+\frac{\partial }{\partial {\mathbf{x}}_1}({\mathbf{f}}_1({\mathbf{x}}_1)+{\mathbf{g}}_1({\mathbf{x}}_1){\mathbf{u}}_1+{\mathbf{h}}_1({\mathbf{x}}_1){\mathbf{d}}_1){|}_{\genfrac{}{}{0ex}{}{{\mathbf{x}}_1={\mathbf{x}}_{1,0}}{{\mathbf{u}}_1={\mathbf{u}}_{1,0}}}\mathrm{\Delta }{\mathbf{x}}_1\hfill \\ +{\mathbf{g}}_1({\mathbf{x}}_{1,0})\mathrm{\Delta }{\mathbf{u}}_1+{\epsilon }_1\hfill \end{array}$$

(12)

where ∂ denotes a partial differential operator, ${\epsilon }_1$ represents higher-order terms, which can be neglected when the sampling period is sufficiently small. And $\mathrm{\Delta }{\mathbf{x}}_1={\mathbf{x}}_1-{\mathbf{x}}_{1,0}$ , $\mathrm{\Delta }{\mathbf{u}}_1={\mathbf{u}}_1-{\mathbf{u}}_{1,0}$ .

According to Table 1 and the TSS principle, the fast incremental terms have a much larger effect on the system than the slow incremental terms, ³⁶ so the slow incremental terms in the equation are negligible, that is, the velocity dynamic equation neglects incremental terms w.r.t. altitude and velocity, while the angle of attack and pitch dynamic equations disregard incremental terms w.r.t. altitude, velocity, angle of attack and pitch.

After performing a first-order Taylor expansion and TSS, the system I system equations are expressed in the following form:

$$\{\begin{array}{c}\stackrel{·}{V}\cong {\stackrel{·}{V}}_0+\frac{\partial }{\partial \alpha }(a_V)\mathrm{\Delta }\alpha +\frac{\partial }{\partial \theta }(a_V)\mathrm{\Delta }\theta +\frac{\partial }{\partial q}(a_V)\mathrm{\Delta }q+g_V\mathrm{\Delta }{\mathbf{u}}_1\hfill \\ \stackrel{·}{\alpha }\cong {\stackrel{·}{\alpha }}_0+\frac{\partial }{\partial q}(a_{\alpha })\mathrm{\Delta }q+g_{\alpha }\mathrm{\Delta }{\mathbf{u}}_1\hfill \\ \stackrel{·}{\theta }\cong {\stackrel{·}{\theta }}_0+\frac{\partial }{\partial q}(a_{\theta })\mathrm{\Delta }q+g_{\theta }\mathrm{\Delta }{\mathbf{u}}_1\hfill \\ \stackrel{·}{q}\cong {\stackrel{·}{q}}_0+g_q\mathrm{\Delta }{\mathbf{u}}_1\hfill \end{array}$$

(13)

Where

$$\begin{array}{c}{a}_V=[-\frac{1}{2m}\rho V^2{S}_{\omega }(c_{31}+c_{41}\alpha +c_{51}{\alpha }^2)-g\sin (\theta -\alpha )]\hfill \\ +\left[\begin{array}{cc}0\hfill & \frac{1}{m}{T}_{{\delta }_t}\cos \alpha \hfill \end{array}\right]{\mathbf{u}}_1\hfill \\ +\left[\begin{array}{c}-\frac{1}{2m}\rho V^2{S}_{\omega }(c_{32}+c_{42}\alpha +c_{52}{\alpha }^2)-\frac{1}{2m}\rho V^2{S}_{\omega }(c_{33}+c_{43}\alpha +c_{53}{\alpha }^2)\hfill \end{array}\right]{\mathbf{d}}_1\hfill \\ a_{\alpha }=[-\frac{1}{2m}\rho V{S}_{\omega }(c_{11}+c_{12}\alpha )+q+\frac{1}{V}g\cos (\theta -\alpha )]\hfill \\ +\left[\begin{array}{cc}-\frac{1}{2m}\rho V{S}_{\omega }{C}_{L{\delta }_e}\hfill & -\frac{1}{\mathrm{mV}}{T}_{{\delta }_t}\sin \alpha \hfill \end{array}\right]{\mathbf{u}}_1\hfill \\ +\left[\begin{array}{cc}-\frac{1}{2m}\rho V{S}_{\omega }(c_{12}+c_{22}\alpha )\hfill & -\frac{1}{2m}\rho V{S}_{\omega }(c_{13}+c_{23}\alpha )\hfill \end{array}\right]{\mathbf{d}}_1\hfill \\ a_{\theta }=\left[q\right]+\left[\begin{array}{cc}0\hfill & 0\hfill \end{array}\right]{\mathbf{u}}_1+\left[\begin{array}{cc}0\hfill & 0\hfill \end{array}\right]{\mathbf{d}}_1\hfill \end{array}$$

$$g_V=\left[\begin{array}{cc}0\hfill & \frac{1}{m}{T}_{{\delta }_t}\cos \alpha \hfill \end{array}\right]$$

$$g_{\alpha }=\left[\begin{array}{cc}-\frac{1}{2m}\rho V{S}_{\omega }{C}_{L{\delta }_e}\hfill & -\frac{1}{\mathrm{mV}}{T}_{{\delta }_t}\sin \alpha \hfill \end{array}\right]$$

$$g_{\theta }=\left[\begin{array}{cc}0\hfill & 0\hfill \end{array}\right]$$

$$g_q=\left[\begin{array}{cc}\frac{1}{2I_y}\rho V^2{S}_{\omega }{c}_A{C}_{m{\delta }_e}\hfill & 0\hfill \end{array}\right]$$

Subsequently, to implement IMPC in the digital environment, the Euler method is used to obtain a discrete-time version of ${\stackrel{·}{x}}_1$ and ${\stackrel{·}{x}}_{1,0}$ .

$${\stackrel{·}{x}}_1=\frac{x_1(k+1)-x_1(k)}{T_s}+\omega (k)$$

(14a)

$${\stackrel{·}{x}}_{1,0}=\frac{x_1(k)-x_1(k-1)}{T_s}+\omega (k-1)$$

(14b)

where $x_1(k)$ is the state variable $x_1$ at time instant $k$ . $T_s$ denotes the time interval for state updates. The incremental element $\mathrm{\Delta }x$ is discretely represented as the difference between successive discrete points, given by $x(k)-x(k-1)$ , $\omega (k)$ is discretization error. Since the sampling frequency is sufficiently fast, the discretization error is ignored. The discrete-time formulations (13) are expressed as follows:

$$\begin{array}{c}\multicolumn{1}{c}{V(k+1)=2V(k)-V(k-1)+T_s\frac{\partial }{\partial \alpha }(a_V)\alpha (k)}\\ -T_s\frac{\partial }{\partial \alpha }(a_V)\alpha (k-1)+\frac{\partial }{\partial \theta }(a_V)\theta (k)\hfill \\ -\frac{\partial }{\partial \theta }(a_V)\theta (k-1)+g_V{\mathbf{u}}_1\hfill \end{array}$$

(15a)

$$\begin{array}{c}\multicolumn{1}{c}{\alpha (k+1)=2\alpha (k)-\alpha (k-1)+\frac{\partial }{\partial q}(a_{\alpha })q(k)}\\ -\frac{\partial }{\partial q}(a_{\alpha })q(k-1)+g_{\alpha }{\mathbf{u}}_1\hfill \end{array}$$

(15b)

$$\begin{array}{c}\multicolumn{1}{c}{\theta (k+1)=2\theta (k)-\theta (k-1)+\frac{\partial }{\partial q}(a_{\theta })q(k)}\\ -\frac{\partial }{\partial q}(a_{\theta })q(k-1)+g_{\theta }{\mathbf{u}}_1\hfill \end{array}$$

(15c)

$$q(k+1)=2q(k)-q(k-1)+g_q{\mathbf{u}}_1$$

(15d)

Define ${\mathbf{X}}_1(k)={\left[\begin{array}{cc}{x}_1(k)\hfill & x_1(k-1)\hfill \end{array}\right]}^{\mathbf{T}}$ , the following linear system is obtained:

$${\mathbf{X}}_1(k+1)={\mathbf{A}}_1{\mathbf{X}}_1(k)+{\mathbf{B}}_1\mathrm{\Delta }{\mathbf{U}}_1(k)$$

(16)

with

$$\begin{array}{c}{\mathbf{A}}_1=\left[\begin{array}{cccccccc}2\hfill & -1\hfill & a_{11}\hfill & -a_{11}\hfill & a_{12}\hfill & a_{12}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 2\hfill & -1\hfill & 0\hfill & 0\hfill & a_{21}\hfill & -a_{21}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 2\hfill & -1\hfill & a_{31}\hfill & -a_{31}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 2\hfill & -1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right]\hfill \\ {\mathbf{B}}_1=\left[\begin{array}{cc}0\hfill & g_{12}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{g_{21}}& g_{22}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{g_{41}}& 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\hfill \end{array}$$

$$\begin{array}{c}{\mathbf{X}}_1(k)=\hfill \\ {\left[\begin{array}{cccccccc}V(k)\hfill & V(k-1)\hfill & \alpha (k)\hfill & \alpha (k-1)\hfill & \theta (k)\hfill & \theta (k-1)\hfill & q(k)\hfill & q(k-1)\hfill \end{array}\right]}^{\mathbf{T}}\hfill \end{array}$$

$$\mathrm{\Delta }{\mathbf{U}}_1(k)={\left[\begin{array}{cc}\mathrm{\Delta }{\delta }_e(k)\hfill & \mathrm{\Delta }{\delta }_t(k)\hfill \end{array}\right]}^{\mathbf{T}}$$

where,

$$\begin{array}{c}\multicolumn{1}{c}{a_{11}=-T_s\frac{1}{2m}\rho V(k{)}^2{S}_{\omega }(c_{41}+2c_{51}\alpha (k))+T_{s}g\cos (\theta (k)-\alpha (k))}\\ -\frac{1}{m}{T}_{{\delta }_t}\sin \alpha {\delta }_t(k)-\frac{1}{2m}\rho V(k{)}^2{S}_{\omega }(c_{42}+2c_{52}\alpha (k)){\xi }_b(k)\hfill \\ -\frac{1}{2m}\rho V(k{)}^2{S}_{\omega }(c_{43}+2c_{53}\alpha (k)){\xi }_s(k)\hfill \end{array}$$

$$\begin{array}{c}{a}_{12}=-T_{s}g\cos (\theta (k)-\alpha (k))\hfill \\ a_{21}=T_s\hfill \\ a_{31}=T_s\hfill \end{array}$$

$$\begin{array}{c}{g}_{12}=T_s\frac{1}{m}{T}_{{\delta }_t}\cos \alpha (k)\hfill \\ g_{21}=-T_s\frac{1}{2m}\rho V(k)S_{\omega }{C}_{L{\delta }_e}\hfill \\ g_{22}=-T_s\frac{1}{\mathrm{mV}(k)}{T}_{{\delta }_t}\sin \alpha (k)\hfill \\ g_{41}=T_s\frac{1}{2I_y}\rho V^2(k)S_{\omega }{c}_A{C}_{m{\delta }_e}\hfill \end{array}$$

Similarly, the discrete increment equation for system II is formulated as follows:

$${\mathbf{X}}_2(k+1)={\mathbf{A}}_2{\mathbf{X}}_2(k)+{\mathbf{B}}_2\mathrm{\Delta }{\mathbf{U}}_2(k)$$

(17)

where ${\mathbf{X}}_2(k)={\left[\begin{array}{cc}h(k)\hfill & h(k-1)\hfill \end{array}\right]}^{\mathbf{T}}$ , $\mathrm{\Delta }{\mathbf{U}}_2(k)=\mathrm{\Delta }{\theta }_d(k)$ , ${\mathbf{A}}_2=\left[\begin{array}{cc}2\hfill & -1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right]$ , ${\mathbf{B}}_2=\left[\begin{array}{c}{g}_2\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$ , $g_2=T_{s}V(k)\cos (\theta (k)-\alpha (k))$ .

The incremental error can be small when the sampling rate is sufficiently high. To guarantee the required approximation accuracy of incremental systems, the sampling period has to be chosen sufficiently small. The control precision is usually consistent with the sampling rate, which is generally up to 100 Hz or even higher for high-precision aircraft systems.^31,37,38 And the above discrete incremental method is feasible.

Optimal control problem

During the controller design for the morphing aircraft, it is typically desired for the aircraft to track a given trajectory rapidly and smoothly while ensuring the smoothness of system inputs and the closed-loop system stability. To guarantee that the morphing aircraft follows the desired trajectory quickly and steadily, the design of the stage cost function is particularly important. At the same time, the stage cost function not only considers the tracking error, but also take into account the energy consumption and control oscillation to achieve optimal control performance.

For this purpose, considering controller I as an example, we define the stage cost using a reference state vector ${\mathbf{X}}_{1_d}={\left[\begin{array}{cccccccc}{V}_d(k)\hfill & V_d(k-1)\hfill & {\alpha }_d(k)\hfill & {\alpha }_d(k-1)\hfill & {\theta }_d(k)\hfill & {\theta }_d(k-1)\hfill & q_d(k)\hfill & q_d(k-1)\hfill \end{array}\right]}^{\mathbf{T}}$ . We assume that the reference trajectory of the aircraft is predetermined, the tracking stage cost is as follows:

$$\ell ({\mathbf{X}}_{1_{k+i|k}},\mathrm{\Delta }{\mathbf{U}}_{1_{k+i|k}},k+i)=\Vert {\mathbf{X}}_{1_{k+i|k}}-{\mathbf{X}}_{1_d}(k+i){\Vert }_{{\mathbf{Q}}_1}^2+\Vert \mathrm{\Delta }{\mathbf{U}}_{1_{k+i|k}}{\Vert }_{{\mathbf{R}}_1}^2$$

(18)

where the notation $\Vert \begin{array}{c}·\hfill \end{array}{\Vert }_{{\mathbf{Q}}_1}$ and $\Vert \begin{array}{c}·\hfill \end{array}{\Vert }_{{\mathbf{R}}_1}$ indicates the squared weighted Euclidean norms. ${\mathbf{Q}}_1$ is the tracking error weight matrix and ${\mathbf{R}}_1$ is the input weight matrix of system I. ${·}_{k+i|k}$ denotes predictions of states and system inputs.

According to the stage cost (18), in order to achieve optimal control performance, we develop the reference tracking IMPC by formulating the following constrained OCP within a finite horizon $N\in {\mathbb{I}}_{>0}$ :

$$V_{{\mathrm{I}}_N}(\mathbf{X}(k),k)=\underset{\mathrm{\Delta }{\hat{\mathbf{U}}}_1(k)}{min}\sum _{i=0}^{N-1}\ell ({\mathbf{X}}_{1_{k+i|k}},\mathrm{\Delta }{\mathbf{U}}_{1_{k+i|k}},k+i)$$

(19a)

$$\mathrm{s}.\mathrm{t}.{\hat{\mathbf{\alpha }}}_{min}\le \hat{\mathbf{\alpha }}(k)\le {\hat{\mathbf{\alpha }}}_{max}$$

(19b)

$${\hat{\mathbf{U}}}_{1min}\le {\hat{\mathbf{U}}}_1(k)\le {\hat{\mathbf{U}}}_{1max}$$

(19c)

• Remark 4. (Optimal Performance): In equation (18), the first term on the right side of the equation reflects the tracking capability of the control system with respect to the desired trajectory, aiming to minimize the tracking error and ensure the aircraft accurately follows the intended path. The second term is employed to reduce the degree of variation in the inputs, reflecting the requirements for smoothness and energy efficiency of system inputs.

According to (16), ${\overline{\mathbf{X}}}_1(k)={\left[\begin{array}{cccc}{\mathbf{X}}_1(k+1)\hfill & {\mathbf{X}}_1(k+2)\hfill & \cdots \hfill & {\mathbf{X}}_1(k+n)\hfill \end{array}\right]}^{\mathbf{T}}$ is defined as the state prediction sequence for system I, $\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)={\left[\begin{array}{cccc}\mathrm{\Delta }{\mathbf{U}}_1(k)\hfill & \mathrm{\Delta }{\mathbf{U}}_1(k+1)\hfill & \cdots \hfill & \mathrm{\Delta }{\mathbf{U}}_1(k+n-1)\hfill \end{array}\right]}^{\mathbf{T}}$ is defined as the incremental control sequence for system I, the incremental model prediction process is expressed as:

$${\overline{\mathbf{X}}}_1(k)={\mathbf{M}}_1{\mathbf{X}}_1(k)+{\mathbf{C}}_1\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)$$

(20)

with

$${\mathbf{M}}_1=\left[\begin{array}{c}{\mathbf{A}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_1^2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathbf{A}}_1^n}\end{array}\right],{\mathbf{C}}_1=\left[\begin{array}{cccc}{\mathbf{B}}_1\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_1{\mathbf{B}}_1}& {\mathbf{B}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_1^{n-1}{\mathbf{B}}_1}& {\mathbf{A}}_1^{n-2}{\mathbf{B}}_1\hfill & \cdots \hfill & {\mathbf{B}}_1\hfill \end{array}\right]$$

In accordance with (20), the OCP is cast into the following QP problem:

$$V_{{\mathrm{I}}_N}(\mathbf{X}(k),k)=min\frac{1}{2}\mathrm{\Delta }{\overline{\mathbf{U}}}_1^{\mathbf{T}}(k){\mathbf{H}}_1\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)+{\mathbf{G}}_1\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)$$

(21a)

$$\mathrm{s}.\mathrm{t}.{\overline{\mathbf{\alpha }}}_{min}\le \overline{\mathbf{\alpha }}(k)\le {\overline{\mathbf{\alpha }}}_{max}$$

(21b)

$${\overline{\mathbf{U}}}_{1min}\le {\overline{\mathbf{U}}}_1(k)\le {\overline{\mathbf{U}}}_{1max}$$

(21c)

where ${\mathbf{H}}_1={{\mathbf{C}}_1}^{\mathbf{T}}{\mathbf{Q}}_1{\mathbf{C}}_1+{\mathbf{R}}_1$ is the Hessian matrix and is a symmetric positive-definite matrix, ${\mathbf{G}}_1={\mathbf{X}}_1^{\mathbf{T}}(k){{\mathbf{M}}_1}^{\mathbf{T}}{\mathbf{Q}}_1{\mathbf{C}}_1-{\overline{\mathbf{X}}}_{1_d}^{\mathbf{T}}(k){\mathbf{Q}}_1{\mathbf{C}}_1$ is the gradient vector. ${\overline{\mathbf{\alpha }}}_{min}$ and ${\overline{\mathbf{\alpha }}}_{max}$ are the minimum and maximum values of the aircraft’s angle of attack, with $\overline{\mathbf{\alpha }}(k)={\left[\begin{array}{cccc}\alpha (k+1)\hfill & \alpha (k+2)\hfill & \cdots \hfill & \alpha (k+n)\hfill \end{array}\right]}^{\mathbf{T}}$ , ${\overline{\mathbf{U}}}_{1min}$ and ${\overline{\mathbf{U}}}_{1max}$ are the minimum and maximum values of the system inputs.

To fit the form with the QP problem, the following transformation is required for (21b) and (21c). For the constraint (21b), with $\alpha (k+1)=\left[\begin{array}{ccc}{\mathbf{O}}_{1\times 2}\hfill & {\mathbf{I}}_{1\times 1}\hfill & {\mathbf{O}}_{1\times 5}\hfill \end{array}\right]{\mathbf{X}}_1(k)={\mathbf{L}}_1{\mathbf{X}}_1(k)$ , $\overline{\mathbf{\alpha }}(k)$ is expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{\overline{\mathbf{\alpha }}(k)=\left[\begin{array}{cccc}{\mathbf{L}}_1\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{L}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{L}}_1\hfill \end{array}\right]{\overline{\mathbf{X}}}_1(k)}\\ =\left[\begin{array}{cccc}{\mathbf{L}}_1\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{L}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{L}}_1\hfill \end{array}\right]{\mathbf{M}}_1{\mathbf{X}}_1(k)+\left[\begin{array}{cccc}{\mathbf{L}}_1\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{L}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{L}}_1\hfill \end{array}\right]{\mathbf{C}}_1\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)\hfill \end{array}$$

(22)

Combining (22) and the constraint equation (21b), we derive the following expression:

$$\left[\begin{array}{cccc}{\mathbf{L}}_1\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{L}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{L}}_1\hfill \end{array}\right]{\mathbf{C}}_1\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)\ge {\overline{\mathbf{\alpha }}}_{min}-\left[\begin{array}{cccc}{\mathbf{L}}_1\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{L}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{L}}_1\hfill \end{array}\right]{\mathbf{M}}_1{\mathbf{X}}_1(k)$$

(23a)

$$\left[\begin{array}{cccc}{\mathbf{L}}_1\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{L}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{L}}_1\hfill \end{array}\right]{\mathbf{C}}_1\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)\le {\overline{\mathbf{\alpha }}}_{max}-\left[\begin{array}{cccc}{\mathbf{L}}_1\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{L}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{L}}_1\hfill \end{array}\right]{\mathbf{M}}_1{\mathbf{X}}_1(k)$$

(23b)

Besides, for the constraint (21c), ${\overline{\mathbf{U}}}_1(k)$ is expressed as:

$${\overline{\mathbf{U}}}_1(k)=\left[\begin{array}{c}{\mathbf{I}}_{2\times 2}\hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}\end{array}\right]{\mathbf{U}}_1(k-1)+\left[\begin{array}{cccc}{\mathbf{I}}_{2\times 2}\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}& {\mathbf{I}}_{2\times 2}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}& {\mathbf{I}}_{2\times 2}\hfill & \cdots \hfill & {\mathbf{I}}_{2\times 2}\hfill \end{array}\right]\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)$$

(24)

Integrating (24) and the constraint (21c), we obtain the following expression:

$$\left[\begin{array}{cccc}{\mathbf{I}}_{2\times 2}\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}& {\mathbf{I}}_{2\times 2}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}& {\mathbf{I}}_{2\times 2}\hfill & \cdots \hfill & {\mathbf{I}}_{2\times 2}\hfill \end{array}\right]\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)\ge {\overline{\mathbf{U}}}_{1min}-\left[\begin{array}{c}{\mathbf{I}}_{2\times 2}\hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}\end{array}\right]{\mathbf{U}}_1(k-1)$$

(25a)

$$\left[\begin{array}{cccc}{\mathbf{I}}_{2\times 2}\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}& {\mathbf{I}}_{2\times 2}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}& {\mathbf{I}}_{2\times 2}\hfill & \cdots \hfill & {\mathbf{I}}_{2\times 2}\hfill \end{array}\right]\mathrm{\Delta }{\overline{\mathbf{U}}}_1(k)\le {\overline{\mathbf{U}}}_{1max}-\left[\begin{array}{c}{\mathbf{I}}_{2\times 2}\hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathbf{I}}_{2\times 2}}\end{array}\right]{\mathbf{U}}_1(k-1)$$

(25b)

Up to this stage, we have formulated the complete IMPC optimal quadratic programing problem for controller I. And we solve the above QP to obtain the current control input ${\mathbf{U}}_1(k)={\mathbf{U}}_1(k-1)+\mathrm{\Delta }{\mathbf{U}}_1(k)$ .

Similar to controller I development, we derive the QP from the controller II OCP:

$$V_{{\mathrm{II}}_N}(\mathbf{X}(k),k)=min\frac{1}{2}\mathrm{\Delta }{\overline{\mathbf{U}}}_2^{\mathbf{T}}(k){\mathbf{H}}_2\mathrm{\Delta }{\overline{\mathbf{U}}}_2(k)+{\mathbf{G}}_2\mathrm{\Delta }{\overline{\mathbf{U}}}_2(k)$$

(26a)

$$\mathrm{s}.\mathrm{t}.{\overline{\mathbf{U}}}_{2max}-{\mathbf{I}}_{N\times 1}{\mathbf{U}}_2(k-1)\le {\mathbf{L}}_2\mathrm{\Delta }{\overline{\mathbf{U}}}_2(k)\le {\overline{\mathbf{U}}}_{2max}-{\mathbf{I}}_{N\times 1}{\mathbf{U}}_2(k-1)$$

(26b)

where: $\mathrm{\Delta }{\overline{\mathbf{U}}}_2(k)={\left[\begin{array}{cccc}\mathrm{\Delta }{\mathbf{U}}_2(k)\hfill & \mathrm{\Delta }{\mathbf{U}}_2(k+1)\hfill & \cdots \hfill & \mathrm{\Delta }{\mathbf{U}}_2(k+n-1)\hfill \end{array}\right]}^{\mathbf{T}}=\mathrm{\Delta }{\overline{\mathbf{\theta }}}_d(k)$ is the incremental control sequence for system II, ${\mathbf{H}}_2={\mathbf{C}}_2^{\mathbf{T}}{\mathbf{Q}}_2{\mathbf{C}}_2+{\mathbf{R}}_2$ is the Hessian matrix of controller II, ${\mathbf{G}}_2={\mathbf{X}}_2^{\mathbf{T}}(k){\mathbf{M}}^{\mathbf{T}}{\mathbf{Q}}_2{\mathbf{C}}_2-{\overline{\mathbf{X}}}_{2_d}^{\mathbf{T}}(k){\mathbf{Q}}_2{\mathbf{C}}_2$ is the gradient vector of controller II, ${\mathbf{Q}}_2$ and ${\mathbf{R}}_2$ respectively are the tracking error weight matrix and input weight matrix, ${\mathbf{L}}_2$ is the lower triangular matrix with all elements 1. The current control input is ${\mathbf{U}}_2(k)={\mathbf{U}}_2(k-1)+\mathrm{\Delta }{\mathbf{U}}_2(k)$ . And also:

$${\mathbf{M}}_2=\left[\begin{array}{c}{\mathbf{A}}_2\hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_2^2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathbf{A}}_2^n}\end{array}\right],{\mathbf{C}}_2=\left[\begin{array}{cccc}{\mathbf{B}}_2\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_2{\mathbf{B}}_2}& {\mathbf{B}}_2\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_2^{n-1}{\mathbf{B}}_2}& {\mathbf{A}}_2^{n-2}{\mathbf{B}}_2\hfill & \cdots \hfill & {\mathbf{B}}_2\hfill \end{array}\right]$$

As a summary, Figure 4 visualizes the structure and the control signal flow of the proposed IMPC.

Figure 4.

Control structure of the IMPC. “TD” denotes the sampling time delay.

Recursive feasibility

In this section, the recursive feasibility of the optimal problem in the local area around the reachable reference trajectory is proven. We will take the constrained OCP (19) for example. Initially, we introduce the definition of the reachable reference trajectory.^20,23,39

• Definition 1. (Reachable Reference Trajectory^20,23,39): The reference trajectory ${\mathbf{X}}_{1_d}$ is deemed reachable if the variables ${\mathbf{X}}_1$ and the associated controller ${\mathbf{U}}_1$ lie within the compact constraint sets ${\overline{\mathbb{X}}}_1$ and ${\overline{\mathbb{u}}}_1$ , where ${\overline{\mathbb{X}}}_1\oplus {\mathbb{C}}_s^{4n}\subseteq {\mathbb{X}}_1$ , ${\overline{\mathbb{u}}}_1\oplus {\mathbb{C}}_r^{4n}\subseteq {\mathbb{U}}_1$ , ${\mathbb{X}}_1$ := $\{{\mathbf{X}}_1:{{\mathbf{X}}_1}_{min}\le {\mathbf{X}}_1\le {{\mathbf{X}}_1}_{max}\}$ , ${\mathbb{U}}_1:=\{{\mathbf{U}}_1:{{\mathbf{U}}_1}_{min}\le {\mathbf{U}}_1\le {{\mathbf{U}}_1}_{max}\}$ , and ⊕ denotes the Minkowski sum.

According to the definition of the Definition 1 and reference, it is deduced that a constant $V_{max}\in {\mathbb{R}}_{>0}$ is exist, ensuring the optimal solutions originating from ${\mathbf{X}}_1(k)$ ( ${\mathbf{X}}_1(k)\in \mathbb{D}$ , $\mathbb{D}:=\{{\mathbf{X}}_1,V_{{\mathrm{I}}_N}({\mathbf{X}}_1(k),k)\le V_{\max }\}$ ) satisfy ${\mathbf{X}}_{1_{k+j|k}}\in {\overline{\mathbb{X}}}_1\prime$ and ${\mathbf{U}}_{1_{k+j-1|k}}\in {\overline{\mathbb{u}}}_1\prime$ for all $i\in {\mathbb{I}}_{[1,N]}$ , where ${\overline{\mathbb{X}}}_1\prime$ and ${\overline{\mathbb{u}}}_1\prime$ represent tightened constraint sets, that is, ${\overline{\mathbb{X}}}_1\prime \oplus {\mathbb{C}}_{s\prime }^{4n}\subseteq {\mathbb{X}}_1$ , ${\overline{\mathbb{u}}}_1\prime \oplus {\mathbb{C}}_{r\prime }^{4n}\subseteq {\mathbb{U}}_1$ , with positive scalars $s\prime$ and $r\prime$ . In summary, the reachable reference trajectory is tractable and strictly resides within the compact constraint sets.

Subsequently, in Theorem 1, the recursive feasibility of predictive control is proven within the local region $\mathbb{D}$ which surrounding the reachable reference trajectory, with the primary task being the demonstration of the regional decrement of $V_{I_N}({\mathbf{X}}_1(k),k)$ .

• Theorem 1. (Recursive Feasibility): Problem (21a)–(21c) permits recursive feasibility in a local area surrounding the reachable reference trajectory.

• Proof: By selecting $V_{I_N}({\mathbf{X}}_{1^{k+1|k}}^{*},k+1)$ as the auxiliary value function, followed by a feasibility analysis of the predicted predictions without terminal components according to Kohler et al., ²⁰ Wang et al., ²³ and Grüne ³⁹ it follows that for an adequately wide prediction horizon $N$ .

$$V_{I_N}({\mathbf{X}}_{1^{k+1|k}}^{*},k+1)-V_{I_N}({\mathbf{X}}_1(k),k)\le -{\alpha }_V(V_{I_N}({\mathbf{X}}_1(k),k))$$

(27)

where ${\alpha }_V(·)\in {\kappa }_{\infty }$ , ${\alpha }_V(·)\le \mathrm{id}(·)$ . Additionally, because of the approximation error ${\epsilon }_1$ and the discretization error $\omega$ , we have $\Vert {\mathbf{X}}_{1k+1|k}^{\ast }-{\mathbf{X}}_1(k+1)\Vert \le \overline{\epsilon }$ , where $\overline{\epsilon }=({\epsilon }_1+\omega )·T_s$ is a composite error.

Then, due to the Lipschitz continuity of $V_{I_N}({\mathbf{X}}_1(k),k)$ , ²³ a constant $K_{\ell }\in {\mathbb{R}}_{>0}$ exists such that:

$$V_{I_N}({\mathbf{X}}_1(k+1),k+1)-V_{I_N}({\mathbf{X}}_{1^{k+1|k}}^{*},k+1)\le K_{\ell }\overline{\epsilon }$$

(28)

Combining (27) with (28), we obtain the property of regional decrement for $V_{I_N}({\mathbf{X}}_1(k),k)$ , that is:

$$\begin{array}{c}{V}_{I_N}({\mathbf{X}}_1(k+1),k+1)-V_{I_N}({\mathbf{X}}_1(k),k)\hfill \\ \le -{\alpha }_V(V_{I_N}({\mathbf{X}}_1(k),k))+{\epsilon }_N\hfill \end{array}$$

(29)

with ${\epsilon }_N:=K_{\ell }\overline{\epsilon }$ represents the cumulative error bound.

Finally, according to the regional decreasing property of $V_{I_N}({\mathbf{X}}_1(k),k)$ , the problem’s recursive feasibility is confirmed. Assuming that the composite error term $\overline{\epsilon }$ is small enough ( $\overline{\epsilon }\le {\alpha }_V(V_{max})/K_{\ell }$ and ${\epsilon }_N\le {\alpha }_V(V_{max})$ ) and ${\mathbf{X}}_1(k)\in \mathbb{D}$ , it holds that:

$$\begin{array}{c}{V}_{I_N}({\mathbf{X}}_1(k+1),k+1)\le (\mathrm{id}-{\alpha }_V)(V_N({\mathbf{X}}_1(k),k))+{\epsilon }_N\hfill \\ \le (\mathrm{id}-{\alpha }_V)(V_{max})+{\alpha }_V(V_{max})\hfill \\ \le V_{max}\hfill \end{array}$$

(30)

Thereby, if ${\mathbf{X}}_1(k)\in \mathbb{D}$ , then ${\mathbf{X}}_1(k+1)\in \mathbb{D}$ . By applying induction, it is proved that $V_{I_N}({\mathbf{X}}_1(k+j),k+j)$ holds for all $j\in {\mathbb{I}}_{>0}$ . Put differently, $\mathbb{D}$ is a positively invariant set. Therefore, $V_{I_N}({\mathbf{X}}_1(k),k)\le V_{max}$ is recursively maintained. Based on the definition of reachable reference trajectory, when ${\mathbf{X}}_1(k)\in \mathbb{D}$ , it satisfies the input constraints and state constraints. That is to say, the feasibility of the problem in the local region near the reachable reference trajectory is proven.

• Remark 5. (Prediction Horizon Determination): For MPC without terminal components, a long prediction horizon is essential for system feasibility. However, in the morphing aircraft systems, the sampling period is notably short. Solving an OCP with a large prediction horizon within such a short period presents significant challenges. Theoretical analyses in Kohler et al., ²⁰ Wang et al., ²³ and Grüne ³⁹ have determined the minimum prediction horizon. However, in practical applications, the prediction horizon is often overestimated. Wang et al. ²³ and Grüne ³⁹ argue that, in practice, a shorter prediction period is sufficient for ensuring local stability. Nevertheless, by analyzing the required prediction horizon in Kohler et al., ²⁰ Wang et al., ²³ we conclude that a larger prediction horizon generally leads to better control performance. Therefore, it provides guidance for adjusting the prediction horizon: within the available computational resources, it is advisable to choose a prediction horizon as large as possible to improve control performance.

Simulation

Simulation setup

In this section, the proposed IMPC method is validated through numerical simulations of the morphing aircraft.

The nonlinear longitudinal model of the morphing aircraft is derived from the detailed aerodynamic data presented in Wang et al. ⁴⁰ The effectiveness of the control method designed in this article is verified through the simulation of the Navion L-17 aircraft model. The parameters of the aircraft are as follows: the mass $m$  = 1247  $\mathrm{kg}$ , $S_{\omega }=17.1m^2$ , $c_A=1.737m$ , the moment of inertia about the vertical axis is $I_y=4067.45\mathrm{kg}·m^2$ , the atmospheric density is $\rho =1.0555\mathrm{kg}/m^3$ .

The aerodynamic parameters of the morphing aircraft are as follows: the thrust coefficient $T_{{\delta }_t}$ is 41.3 $N/\%$ ; $C_{L_{{\delta }_e}}$ is 0.3209, $C_{m_{{\delta }_e}}$ is −1.0199. The coefficients matrix $\mathbf{C}$ in the equation (4) is:

$$\mathbf{C}=[c_{\mathrm{ij}}{]}_{7\times 3}=\left[\begin{array}{ccc}0.3984\hfill & 0.3390\hfill & -0.2549\hfill \\ \multicolumn{1}{c}{6.5783}& 4.1114\hfill & -3.6219\hfill \\ \multicolumn{1}{c}{0.0280}& 0.0113\hfill & -0.0077\hfill \\ \multicolumn{1}{c}{0.2137}& 0.0583\hfill & -0.1173\hfill \\ \multicolumn{1}{c}{1.9421}& 0.0619\hfill & -0.8077\hfill \\ \multicolumn{1}{c}{0.0876}& -0.4863\hfill & -0.3455\hfill \\ \multicolumn{1}{c}{-0.4765}& -4.9762\hfill & -4.9473\hfill \end{array}\right]$$

The controller parameters are selected as follows: The prediction horizon $N$ =60, and the sampling period $T_s=5$ ms. The weighting matrices $\mathbf{Q}$ and $\mathbf{R}$ are expressed in the following form:

$${\mathbf{Q}}_i=\left[\begin{array}{cccc}{\mathbf{Q}}_{\mathrm{Li}}\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{Q}}_{\mathrm{Li}}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{Q}}_{\mathrm{Li}}\hfill \end{array}\right],{\mathbf{R}}_i=\left[\begin{array}{cccc}{\mathbf{R}}_{\mathrm{Li}}\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{R}}_{\mathrm{Li}}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{R}}_{\mathrm{Li}}\hfill \end{array}\right],(i=1,2)$$

where ${\mathbf{Q}}_{L1}$ and ${\mathbf{R}}_{L2}$ are $\mathrm{diag}\{100,0\}$ and $\mathrm{diag}\{0.01\}$ respectively, and ${\mathbf{Q}}_{L2}$ = $\mathrm{diag}\{0.05,0,0,0,0.01,0,0,0\}$ , ${\mathbf{R}}_{L2}$ = $\mathrm{diag}\{0.001,0.001\}$ .

To validate the proposed IMPC in this article, three scenarios are conducted. Firstly, to evaluate the tracking control performance and deformation stability of the IMPC method, a comparative experiment with the state-of-the-art controllers was conducted under identical targets and conditions. Besides, an experimental scenario with model uncertainty is established, and the robustness of the IMPC method under imprecise modeling is verified by comparative experiments. Additionally, the IMPC method’s computational efficiency is compared with other nonlinear model predictive control (NMPC) methods to demonstrate its superior efficiency.

Scenarios and results

Scenario 1: (Tracking feasibility)

To verify the optimal tracking performance of the proposed IMPC, the following target tracking scenario is designed for the morphing aircraft.

The flight instructions for the aircraft are depicted in Figure 5. Initially, the aircraft maintains a constant velocity of 50 m/s at an altitude of 4900 m. At $t=20$  s, the aircraft begins to accelerate, achieving a velocity of 60 m/s by $t=50$  s and is sustained for the remainder of the flight. Aircraft descent initiates at $t=30$  s, and by $t=90$  s, the aircraft reaches an altitude of 5000 m. Subsequently, the aircraft maintains this velocity and altitude for the duration of the flight.

Figure 5.

Reference of velocity and altitude.

Controllers of aircraft with the same parameters are selected for comparison, including linear parameter-varying control (LPV) ¹⁰ and adaptive sliding mode control (ASMC). ²⁷ It is assumed that the wingspan deformation rate ${\xi }_b$ and sweep angle deformation rate ${\xi }_s$ of the aircraft are fixed as 0. The tracking errors and control inputs under the above flight path are shown in Figures 6–9:

Figure 6.

Tracking errors of velocity under different controllers.

Figure 7.

Tracking errors of altitude under different controllers.

Figure 8.

Elevator deflection under different controllers.

Figure 9.

Throttle opening under different controllers.

To quantitatively analyze the tracking error of the controller, the root mean square (RMS) values of velocity and altitude tracking errors are calculated, as presented in Table 2:

Table 2.

Root mean square values of tracking errors.

State	IMPC	ASMC	LPV
$V$	0.01077	0.03853	0.14440
$h$	0.07406	0.28366	0.89757

From the simulation results, as shown in Figures 6 and 7 and Table 2, it is observed that IMPC exhibits smaller tracking errors, faster response speed and superior tracking performance compared to other controllers. In addition, Figures 8 and 9 and Table 2 show that the control input of IMPC is relatively small and exhibits less fluctuation. It achieves optimal control effect.

Subsequently, to verify the attitude stability of the morphing aircraft during the deformation process, the flight command for the aircraft was set as follows: The aircraft flies at a constant velocity of 80 m/s at an altitude of 6000 m, with the wingspan length ${\xi }_b$ undergoing continuous and even change between 10 and 40 s, and the sweep angle ${\xi }_s$ varying continuously and evenly between 50 and 80 s, as shown in Figure 10.

Figure 10.

Wingspan and sweep angle deformation rates over time.

The tracking errors of aircraft’s velocity and altitude under sweep angle and wingspan length deformation are depicted in Figures 11 and 12, with the RMS values in Table 3.

Figure 11.

Tracking errors of velocity under deformation.

Figure 12.

Tracking errors of altitude under deformation.

Table 3.

RMS values of tracking errors under deformation.

State	IMPC	ASMC	LPV
$V$	0.00006261	0.00017713	0.00077345
$h$	0.06244	0.19176	0.31814

Obviously, the controller eliminates the instability caused by the wing deformation, and ensures the constant velocity and constant altitude flight in the process of changing wingspan and sweep angle. After completing the deformation process, the aircraft is able to maintain the velocity and altitude before the deformation and enter a new equilibrium state.

Scenario 2: (Robustness)

In this section, in order to verify the robustness of the proposed method to the model uncertainty, Parameter uncertainty is introduced into the controlled object in the simulation of morphing aircraft.

There are some approximate steps in the building process of the morphing aircraft model, and the modeling parameters are also uncertain. Therefore, the aerodynamic modeling in the nonlinear model of morphing aircraft is not accurate. The controller must have the ability to ensure stability in case of modeling errors. The unpredictable aspects of the model are characterized by its parametric perturbations. In this scenario, we describe the uncertainty of the model parameters by adding perturbations to the aerodynamic coefficients as follows:

$$\{\begin{array}{c}{C}_L\prime =(1+\varphi (t))C_L\hfill \\ \multicolumn{1}{c}{C_D\prime =(1+\varphi (t))C_D}\\ \multicolumn{1}{c}{C_m\prime =(1+\varphi (t))C_m}\\ \multicolumn{1}{c}{T_{\delta }\prime =(1+\varphi (t))T_{\delta }}\end{array}$$

where, $C\prime$ denotes the perturbed value, and $C$ represents the nominal value. $\varphi (t)$ follows a normal distribution function $N(0,{\sigma }^2)$ . The method proposed in this article and the comparative method utilizes identical sequences of perturbation samples. Comparative results are shown in Figures 13–16 and the RMS values in Table 4:

Figure 13.

Tracking errors of velocity with $\sigma =5\%$ .

Figure 14.

Tracking errors of altitude with $\sigma =5\%$ .

Figure 15.

Tracking errors of velocity with $\sigma =20\%$ .

Figure 16.

Tracking errors of altitude with $\sigma =20\%$ .

Table 4.

RMS values with model perturbations.

State	$\sigma =5\%$		$\sigma =20\%$
	IMPC	ASMC	IMPC	ASMC
$V$	0.0010083	0.0021932	0.0021852	0.0059785
$h$	0.0380934	0.2250289	0.0738087	0.4203995

As shown in Figures 13–16 and Table 4, introducing interference to the model significantly impacts the tracking performance of the ASMC controller. The tracking fluctuation at $\sigma =20\%$ is greater than that at $\sigma =5\%$ . In this case, IMPC effectively resists the interference caused by model uncertainty due to its independence from concrete mathematical models.

In addition, to validate the robustness of the proposed method under constant uncertainty, the following simulation was designed: within the aforementioned simulation framework, a constant uncertainty of 20% was applied to the aerodynamic coefficients at 50 s.

$$\{\begin{array}{c}{C}_{L^{\prime }}=(1+20\%)C_L\\ C_{D^{\prime }}=(1+20\%)C_D\\ C_{m^{\prime }}=(1+20\%)C_m\\ T_{\delta }^{\prime }=(1+20\%)T_{\delta }\end{array}$$

The simulation results are shown in Figures 17 and 18.

Figure 17.

Tracking errors of altitude with constant uncertainty.

Figure 18.

Tracking errors of velocity with constant uncertainty.

From the results, it can be observed that the IMPC approach exhibits relatively small errors under the constant uncertainty and does not produce steady-state errors due to model uncertainty.

Scenario 3: (Computational efficiency)

In this section, we statistically analyze the computational efficiency of the controller in Scenario 1. Traditional NMPC uses nonlinear models for prediction, thus the cost function to be solved is usually a non-convex function, with the existence of local optimum. Solving the global optimal value of predictive control requires more complex computational methods. The process of using models for prediction is a recursive process, and the predictive function is a nested function, which brings more complex computation to the non-convex optimization process of NMPC.

Solving the non-convex optimization process of NMPC using sequential quadratic programing (SQP) and particle swarm optimization (PSO) ⁴¹ is a widely used method, and this section conducts simulation experiments on the morphing aircraft using these methods. The calculation time of the model predictive control solution for each cycle is statistically analyzed, and the sampling period of each controller is set to 0.005 s. The statistical results are shown in Figure 19, which were obtained under the same simulation computer. The simulation is created by an AMD™ (Ryzen 7 5800H @3.2 GHz) CPU computer.

Figure 19.

Solution time statistics.

In order to quantitatively analyze the calculation effect, the mean and standard deviation of the single-step solution time are calculated, and the results are shown in Table 5:

Table 5.

Mean and standard deviation for solution time (s).

statistics	IMPC	SQP-NMPC	PSO-NMPC
Mean	0.00222	0.04729	0.14051
Standard deviation	0.00084	0.00119	0.00333

As illustrated in Figure 19 and Table 5, the average IMPC solution time is much smaller than that of SQP-NMPC and PSO-NMPC, indicating that the average IMPC solution time is relatively short. In the context of the sampling period (5 ms), real-time control is guaranteed. Furthermore, the standard deviation of the IMPC solution time is also small, indicating that the IMPC solution time remains consistent at each step, thus demonstrating its computational stability.

Conclusion

In this article, an Incremental Model Predictive Control approach is proposed for the trajectory tracking of the morphing aircraft. To mitigate the dependence on concrete mathematical models, incremental system approximations are employed with TSS for the equations of system dynamics. It improves the robustness of the controller in terms of high tracking accuracy. Utilizing the approximated discrete linear system, the IMPC is formulated to comply with state and input constraints, thereby obviating the need for a concrete mathematical model. In addition, IMPC results in a linear MPC. Compared with nonlinear MPC, the computational complexity of IMPC dramatically decreases. Finally, the efficacy of the proposed IMPC is validated through a series of simulations. Simulation results demonstrate that the IMPC achieves optimal control performance under state and input constraints, eliminates the dependence on a concrete mathematical model, and enhances the computational efficiency.