Predictive Control of Trajectory Tracking for Flapping-Wing Aircraft Based on Linear Active Disturbance Rejection

Abstract

This article discusses the problem of controlling the trajectory of a flapping-wing aircraft in the face of external disturbances. As the applications for flapping-wing aircraft have diversified, the external disturbances to which the system is exposed have become more complex. Existing control methods have difficulty with effectively counteracting these disturbances. Therefore, this paper suggests a control method that combines linear active disturbance rejection with model predictive control to solve the tracking problem under disturbances, improve the system’s disturbance rejection capability, and ensure the accuracy of trajectory tracking. First, a linear active disturbance controller (LADRC) is developed for the position system to monitor and compensate for internal uncertainties and environmental disturbances in a timely manner. Secondly, the attitude control system is equipped with a model predictive controller (MPC) to effectively determine the optimal control variables and achieve stable attitude tracking. The method is evaluated through simulation studies to assess its performance in tracking a reference trajectory in the presence of disturbances. The findings demonstrate that the approach can accurately track the reference trajectory even when the system is subject to sinusoidal disturbances. This indicates that the method exhibits robustness and practicality.

1. Introduction

In recent years, the development of drone technology has led to the emergence of flapping-wing aircraft [1]. These aircraft mimic the wing movement principles of insects and offer a highly efficient and adaptable alternative to traditional fixed-wing drones. Flapping-wing aircraft demonstrate superior flight efficiency, agility, and performance in complex environments by emulating the flight capabilities of natural organisms. The low noise emission, compact design, and suitability for novel applications have made flapping-wing aircraft a focal point of research interest. Biomimetics is trying to draw inspiration from the wing structure, feather design, flight posture adjustment, and energy consumption of birds to develop more efficient and stable aircraft. In bird flight, their lift is mainly caused by aerodynamic forces generated by forward flight speed, and their flight direction is changed by adjusting flaps and tails. The main challenge of the flapping-wing aircraft is the design of the trajectory tracking controller [2]. However, the unique flight mechanics and nonlinear dynamic characteristics of flapping-wing aircraft mean that conventional linear control methods often fail to meet the requirements for trajectory tracking. To address this, extensive research has been conducted on control algorithms for flapping-wing aircraft, including PIDs, linear quadratic regulators, neural network control, and model-based predictive control.

Due to the challenges in modelling flapping-wing aircraft (FWA), the predominant control method remains the PID algorithm [3,4,5]. G. M. Li [6] proposed an FWA autonomous control architecture based on PID control, which enables accurate adjustment of FWA pitch, yaw, and roll angles. Wei Liao [7] designed attitude position control, altitude control, turn control, heading control, and semi-autonomous tracking flight control using PID control based on negative feedback. Jiang et al. [8] employed a cascaded PID algorithm to control the attitude of flapping-wing aircraft. Furthermore, it incorporates the L1 guidance law into the horizontal navigation task of the flapping-wing aircraft, enabling the aircraft to accomplish trajectory tracking missions. However, PID control lacks adaptability and cannot be applied to systems with variable dynamics.

Compared to the PID control method, the LQR control method is more flexible and accurate, particularly in controlling complex nonlinear and time-varying systems. The LQR method designs an objective function, which is a quadratic function of state and control inputs, to obtain the optimal control law. In the face of external perturbations and modelling uncertainties, LQR is highly adaptive. Kim [9] introduced an aerodynamic model that considers leading edge vortices interacting with the wings and tail. The best control gain matrices were calculated for a range of flight speeds using the LQR method. Gao Hejia [10] obtained a state-space model for the hovering of a flapping-wing aircraft under the Huangtai scenario through system identification. An LQR optimal controller was designed, and experimental results showed that adjusting the wing angles could achieve position and attitude trajectory tracking control for the aircraft. However, LQR control is not robust to parameter changes.

In order to improve the robustness of the system to parameter changes. Gao et al. [11] addressed the issues of dynamic uncertainty and actuator failure in flapping-wing aircraft systems by proposing a neural network algorithm to approximate the system’s dynamic uncertainty and actuator failure, effectively improving the robustness of the system. Shuai et al. [12] addressed the tracking control problem of robot systems with actuator failures. A neural network integral sliding mode fault-tolerant controller was designed to estimate the robot’s actuator failures. This approach effectively tracked the desired trajectory and enhanced the system’s robustness and fault tolerance. In addition, Shi et al. [13] proposed nonlinear $H_{\infty }$ attitude control to overcome the effects of time-varying parameters and unknown disturbances on the system. Fatehi et al. [14] used the u-analysis method to study the robust control of the flapping wing and predicted the wing flutter phenomenon in the presence of multiple uncertainties.

Previous studies have focused on controlling FWA. However, due to the uncertainty of the natural environment, such as wind disturbance, these control methods may not be immune enough. Compared to the previously mentioned advanced control algorithms, model predictive control and linear active disturbance rejection [15,16] exhibit greater control immunity. In their work, Z. Q. Gao [17] proposed LADRC that has significantly fewer parameters than the nonlinear ADRC controller. This makes it more convenient for parameter adjustment and algorithm implementation while maintaining immunity. On the other hand, MPC [18] emphasises the functionality of the model. It does not focus on its inner structure. Unlike traditional control methods, MPC optimization is limited to a moving time window. This optimization strategy is suitable for trajectory tracking systems with dynamic changes and uncertainties, enhancing the real-time control, robustness, and disturbance rejection capabilities of FWA during flight.

This paper proposes a model predictive control method based on linear active disturbance rejection (LADRC-MPC) to address the strong nonlinear characteristics of flapping-wing aircraft and their susceptibility to internal and external uncertainty disturbances. The proposed method combines the advantages of LADRC and MPC. LADRC is used to estimate internal and external disturbances in real-time, while MPC enables stable tracking of attitude control variables. The integration of these two methods increases the robustness and adaptability of the control system. Simulation results verify the validity of the suggested control method. The findings indicate that LADRC-MPC outperforms conventional control techniques by achieving more accurate and stable trajectory tracking amidst external disturbances and uncertainties.

The paper presents the structure of the trajectory tracking control algorithm as follows: In Section 2, a dynamic model of a flapping-wing aircraft is presented, and uncertain external disturbances are considered. Section 3 presents the control design and mathematical formulations. Section 4 describes two different trajectories for simulation and analysis of experimental results. Section 5 presents conclusions on the work discussed above.

2. Modelling of Flapping-Wing Aircraft

In this section, considering the flapping-wing aircraft as a rigid body, the mathematical model of the system is derived through dynamic analysis using Euler’s equations of motion. A flapping-wing aircraft is a driving system with six degrees of freedom. Figure 1 illustrates a schematic diagram of a flapping-wing aircraft with an Earth-referenced coordinate system.

Compared to fixed-wing and quadrotor aircraft, flapping wings generate the necessary lift and thrust for flight through wing vibration. This complexity in the flow field presents a significant challenge for modelling flapping-wing aircraft. For the convenience of model construction, we make the following assumptions [19].

Assumption 1.

The aerodynamic force and resistance of the motor can be ignored.

Assumption 2.

The angle of attack and pitch angle of the wings are assumed to be small, allowing for linearization of the equations of motion.

In the volume coordinate system of the flapping-wing aircraft, the force experienced by its linear motion and the moment experienced by its rotational motion can be expressed as:

$$\begin{array}{c}\hfill {\left[\begin{array}{ccc}{F}_x\hfill & F_y\hfill & F_z\hfill \end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{ccc}m\dot{u}\hfill & m\dot{v}\hfill & m\dot{w}\hfill \end{array}\right]}^{\mathrm{T}}\\ \hfill {\left[\begin{array}{ccc}{M}_x\hfill & M_y\hfill & M_z\hfill \end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{ccc}I\dot{p}\hfill & I\dot{q}\hfill & I\dot{r}\hfill \end{array}\right]}^{\mathrm{T}}\end{array}$$

(1)

where m represents the mass of the aircraft, and I is the inertia matrix.

$$\begin{array}{c}\hfill I=\left[\begin{array}{ccc}{I}_x& I_{xy}& I_{zx}\\ I_{xy}& I_y& I_{yz}\\ I_{zx}& I_{yz}& I_z\end{array}\right]\end{array}$$

(2)

Taking into account the symmetry of the inertial moments of the flapping-wing aircraft, $I_{\mathrm{x}\mathrm{y}}=I_{\mathrm{y}\mathrm{z}}=I_{\mathrm{zx}}=0$.

The flapping-wing aircraft system can be represented in the body frame as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}m\dot{V}=F-m(\omega +V)\hfill \\ I\dot{\omega }=M-\omega \times I\omega \hfill \end{array}\right.\end{array}$$

(3)

where F denotes the total force acting on the aircraft, M is the total torque applied to the aircraft, $V={\left[\begin{array}{ccc}u\hfill & v\hfill & w\hfill \end{array}\right]}^{\mathrm{T}}$ represents the velocity matrix of the aircraft, and $\omega ={\left[\begin{array}{ccc}p\hfill & q\hfill & r\hfill \end{array}\right]}^{\mathrm{T}}$ is the matrix of angular velocities of the flapping-wing aircraft around each axis.

The lift generated by the wings of a flapping aircraft can be expressed as:

$$F_T=\left[\begin{array}{c}{F}_{Tx}\\ F_{Ty}\\ F_{Tz}\end{array}\right]=\left[\begin{array}{c}2{\int }_0^b(\rho /2)C_L{V}^{2}C\left(y\right)dy·cos(\alpha +\phi )cos\varphi -(\rho /2)C_D{V}^{2}C\left(y\right)dy·sin(\alpha +\phi )cos\varphi \\ 0\\ 2{\int }_0^b(\rho /2)C_Z{V}^{2}C\left(y\right)dy·sin(\alpha +\phi )-(\rho /2)C_D{V}^{2}C\left(y\right)dy·cos(\alpha +\phi )\end{array}\right]$$

(4)

where $C_L$ is the lift coefficient of the flapping-wing aircraft wings, $C_D$ is the drag coefficient of the flapping-wing aircraft wings, b is the length of one side’s wing extension, and V is the body velocity.

The drag experienced by the fuselage of a flapping-wing aircraft can be expressed as:

$$f=\left[\begin{array}{c}{f}_x\hfill \\ f_y\hfill \\ f_z\hfill \end{array}\right]=\left[\begin{array}{c}(\rho /2)v^2{S}_b\left(-c_{x}cos\alpha cos\beta +c_{y}cos\alpha sin\beta -c_{z}sin\alpha \right)\\ {(\rho /2)}^2{S}_b\left(c_{x}sin\beta +c_{z}cos\beta \right)\\ (\rho /2)v_b^2{S}_b\left(-c_{x}sin\alpha cos\beta +c_{y}sin\alpha sin\beta +c_{z}cos\alpha \right)\end{array}\right]$$

(5)

The lift generated by the tail wings of a flapping-wing aircraft can be expressed as:

$$F_t=\left[\begin{array}{c}{F}_{tx}\hfill \\ F_{ty}\hfill \\ F_{tz}\hfill \end{array}\right]=\left[\begin{array}{c}cos{\alpha }_{0}cos\beta ·F_{Vx}+cos{\alpha }_{0}sin\beta ·F_{Vy}-sin{\alpha }_0·F_{Vz}\hfill \\ -sin\beta ·F_{Vx}+cos\beta ·F_{Vy}\hfill \\ sin{\alpha }_{0}cos\beta ·F_{Vx}+sin{\alpha }_{0}sin\beta ·F_{Vy}+cos{\alpha }_0·F_{Vz}\hfill \end{array}\right]$$

(6)

where $F_{Vx},F_{Vy},F_{Vz}$ represents the aerodynamic force on the tail wing in the velocity coordinate system.

The moments of each axis of the flapping-wing aircraft in the framework system of coordinates can be given as follows:

$$\left\{\begin{array}{c}{M}_x=M_{bx}+M_{wx}+M_{tx}\\ M_y=M_{by}+M_{wy}+M_{tx}\\ M_z=M_{bz}+M_{wz}+M_{tx}\end{array}\right.$$

(7)

where $M_b,M_w,M_t$ represents the fuselage damping torque, wing aerodynamic torque, and tail aerodynamic torque, respectively.

During actual flight, the position and direction of the flapping-wing aircraft are often described using Euler angles, which include roll, pitch, and yaw angles to represent the aircraft’s attitude in space. From the conversion of the inertial reference system to the object frame of reference, the equation for the conversion between the Euler angle ${\left[\begin{array}{ccc}\dot{\varphi }\hfill & \dot{\theta }\hfill & \dot{\psi }\hfill \end{array}\right]}^{\mathrm{T}}$ and the angular velocity ${\left[\begin{array}{ccc}p\hfill & q\hfill & r\hfill \end{array}\right]}^{\mathrm{T}}$ components can be derived as follows:

$$\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]=\left[\begin{array}{ccc}1& sin\varphi tan\theta & cos\varphi tan\theta \\ 0& cos\varphi & -sin\varphi \\ 0& sin\varphi /cos\theta & cos\varphi /cos\theta \end{array}\right]\left[\begin{array}{c}p\hfill \\ q\hfill \\ r\hfill \end{array}\right]$$

(8)

From Equations (1)–(3), it can be inferred that the nonlinear equations of the flapping-wing aircraft system can be written as:

$$\left\{\begin{array}{c}\dot{u}=\left(F_x-mgsin\theta \right)/m+vr-wq\hfill \\ \dot{v}=\left(F_y-mgcos\theta cos\varphi \right)/m+pw-ru\hfill \\ \dot{w}=\left(F_z-mgcos\theta sin\varphi \right)/m+qu-pv\hfill \\ \dot{p}=\frac{1}{I_x}\left(\left(I_y-I_z\right)qr+M_x\right)\hfill \\ \dot{q}=\frac{1}{I_y}\left(\left(I_z-I_x\right)pr+M_y\right)\hfill \\ \dot{r}=\frac{1}{I_z}\left(\left(I_x-I_y\right)pq+M_z\right)\hfill \end{array}\right.$$

(9)

The physical parameters employed in the aforementioned modelling and simulation are presented in Table 1.

External atmospheric disturbances, such as gusts, are disturbances during outdoor flights of flapping-wing aircraft and can impact the trajectory tracking stability of the aircraft. The model of the flapping-wing aircraft considering disturbances $o={\left[\begin{array}{cc}{o}_1\hfill & o_2\hfill \end{array}\right]}^{\mathrm{T}}$ is represented as:

$$\left\{\begin{array}{c}m\dot{V}=F-m(\omega +V)+o_1\hfill \\ I\dot{\omega }=M-\omega \times I\omega +o_2\hfill \end{array}\right.$$

(10)

where $o_1$ represents the force disturbance modelling, and $o_2$ represents the torque disturbance modelling. These disturbances include both internal disturbances within the system and external random disturbances from the environment.

3. Control System Design

This article adopts a dual closed-loop control scheme for the trajectory control of flapping-wing aircraft. The attitude loop is mainly responsible for controlling the aircraft attitude, including pitch, roll, and yaw angles, to maintain a stable flight attitude. The position loop is responsible for controlling the position of the aircraft, including altitude and horizontal position, to achieve the targetted position control. The design is mainly based on LADRC position control law and MPC attitude control law. Firstly, we calculate $F_x,F_y,F_z$ based on the reference trajectory and position controller, and then $M_x,M_y,M_z$ is obtained by bringing the given expected angle ${\psi }_{ref}$, ${\theta }_{ref},{\varphi }_{ref}$ into the attitude model prediction controller; finally, $F_x,F_y,F_z,M_x,M_y,M_z$ is substituted into the flapping-wing aircraft model, and the actual position and attitude state variables of the flapping-wing aircraft are fed back to their respective controllers to form a closed-loop control system. Figure 2 depicts the flapping-wing aircraft control system diagram.

3.1. Position Controller

This article uses second-order LADRC to design a position controller, which consists of three parts: tracking differentiator (TD), linear extended state observer (LESO), and linear error feedback control law (LESF). An outer loop position control system is responsible for controlling the position of an object in an inertial reference frame. The translation motion equations can describe how the position and velocity of an object in space change over time, providing a basis for achieving position control in the inertial reference frame. The translation motion equations of the flapping-wing aircraft can be derived from Equation (11) and expressed as Equation (12).

$$\begin{array}{c}\hfill a=R\times F/m=R\times \left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]/m\end{array}$$

(11)

where $R=\left[\begin{array}{ccc}cos\psi cos\theta & sin\theta cos\psi sin\varphi -cos\varphi sin\psi & sin\theta cos\psi cos\varphi +sin\varphi sin\psi \\ sin\psi cos\theta & sin\theta sin\psi sin\varphi +cos\varphi cos\psi & sin\theta sin\psi cos\varphi -sin\varphi cos\psi \\ -sin\theta & sin\varphi cos\theta & cos\theta cos\varphi \end{array}\right]$

$$\left\{\begin{array}{c}\ddot{x}=\left(F_{x}cos\psi cos\theta +F_y(sin\theta cos\psi sin\varphi -cos\varphi sin\psi )+F_z(cos\varphi sin\theta cos\psi +sin\varphi sin\psi )\right)/m\hfill \\ \ddot{y}=\left(F_{x}sin\psi cos\theta +F_y(sin\theta sin\psi sin\varphi +cos\varphi cos\psi )+F_z(cos\varphi sin\theta cos\psi -sin\varphi cos\psi )\right)/m\hfill \\ \ddot{z}=\left(F_{y}sin\varphi cos\theta -F_{x}sin\theta +F_{z}cos\varphi cos\theta \right)/m-g\hfill \end{array}\right.$$

(12)

From the translational dynamics model of flapping-wing aircraft, it can be seen that there is strong coupling between the roll, pitch, and yaw channels of the flapping-wing aircraft, making it difficult to design the controller. However, the LESO in LADRC can effectively solve this problem. The coupled and unmodelled parts in each channel are treated as total disturbances and observed by LESO for real-time feedback compensation, thus achieving decoupled control of each channel. The original nonlinear and strongly coupled control system is transformed into a purely linear system. Based on the small-angle assumption, the translational equation of motion of the flapping-wing aircraft is given as follows:

$$\left\{\begin{array}{c}\ddot{x}=\left(F_{x}cos\psi cos\theta +F_{z}cos\varphi sin\theta cos\psi \right)/m\hfill \\ \ddot{y}=\left(F_{y}cos\varphi cos\psi -F_{z}sin\varphi cos\psi \right)/m\hfill \\ \ddot{z}=F_{z}cos\varphi cos\theta /m-g\hfill \end{array}\right.$$

(13)

By defining the virtual control quantity $U_x,U_y,U_z$, Equation (13) can be derived as follows:

$$\left\{\begin{array}{c}{F}_x=\frac{U_{z}m-F_{z}cos\varphi sin\theta cos\psi }{cos\psi cos\theta }\hfill \\ F_y=\frac{U_{z}m+F_{z}sin\varphi cos\psi }{cos\varphi cos\psi }\hfill \\ F_z=\frac{\left(U_z+g\right)m}{cos\varphi cos\theta }\hfill \end{array}\right.$$

(14)

The coupling between each position channel can be viewed by LADRC technology as an indeterminate internal disturbance that is included in the overall disturbance. Real-time estimation of the current unmodelled dynamic disturbances is possible thanks to the LESO in the LADRC controller. This eliminates the need to decouple all the channels in the system, as each can be controlled independently.

3.1.1. Z-Channel Position Controller Design

In this autonomous channel, our goal is to obtain the control variable $U_z$ from the reference input z. Once we have successfully controlled this variable, we can use Equation (14) to calculate the corresponding control variable $F_z$.

The algorithm for the tracking differentiator (TD) can be designed as:

$$\left\{\begin{array}{c}{\dot{v}}_1=v_2\hfill \\ {\dot{v}}_2=flan\left(v_1-v_d,v_2,R,h\right)\hfill \end{array}\right.$$

(15)

where $v_1$ is the estimated value of $v_d$, $v_2$ is the derivative of $v_1$, R is the velocity factor that determines the convergent velocity, h is the filter factor, and flan() is the optimally synthesised fast control law.

Secondly, definition error $e=z_1-y$.

The linear extended state observer is the core part of the linear self-immunity control, and its function is to estimate the system state and the total perturbation affecting the output in real time according to the system inputs and outputs, and the estimated perturbation is called the expanded state. The linear extended state observer is designed as:

$$\left\{\begin{array}{c}{\dot{z}}_1=z_2+{\beta }_{1}e\hfill \\ {\dot{z}}_2=z_3+{\beta }_{2}e+b_{0}u\hfill \\ {\dot{z}}_3={\beta }_{3}e\hfill \end{array}\right.$$

(16)

where $z_1,z_2$ represent the observed values of $v_1,v_2$, $z_3$ represents the total disturbance of the altitude channel, ${\beta }_1,{\beta }_2,{\beta }_3$ represent the observer gain, and y represents the output of the altitude channel. The value of ${\beta }_{\mathit{i}}$ can be selected as: ${\beta }_1=3w_0,{\beta }_2=3w_0^2,{\beta }_3=w_0^3$; $w_0$ is the observer bandwidth.

Finally, the control law can be designed as:

$$\left\{\begin{array}{c}{U}_z=k_p\left(v_1-z_1\right)+k_d\left(v_2-z_2\right)\hfill \\ F_z=\frac{\left(U_z+g\right)m}{cos\varphi cos\theta }\hfill \end{array}\right.$$

(17)

where $k_p,k_d$ is the gain parameter of the controller in the altitude channel: $k_p=w_c^2,k_d=2w_c$. The controller bandwidth, $w_c$, is equal to 0.25 times the observer bandwidth, $w_0$.

3.1.2. X-Channel Position Controller Design

In this channel, the virtual control quantity $U_x$ is defined, and the control law for $x-\theta$ can be designed based on Equations (14)–(17).

$$\left\{\begin{array}{c}{\dot{v}}_{x1}=v_{x2}\hfill \\ {\dot{v}}_{x2}=flan\left(v_{x1}-v_d,v_{x2},R,h\right)\hfill \\ e=z_{x1}-y_x\hfill \\ {\dot{z}}_{x1}=z_{x2}+{\beta }_{1}e\hfill \\ {\dot{z}}_{x2}=z_{x3}+{\beta }_{2}e+b_{0}u\hfill \\ {\dot{z}}_{x3}={\beta }_{3}e\hfill \\ U_x=k_p\left(v_{x1}-z_{x1}\right)+k_d\left(v_{x2}-z_{x2}\right)\hfill \\ F_x=\frac{U_{z}m-F_{z}cos\varphi sin\theta cos\psi }{cos\psi cos\theta }\hfill \end{array}\right.$$

(18)

where $z_{x1},z_{x2}$ represents the observed value of $v_{x1},v_{x2}$, $z_{x3}$ represents the observed value of the total disturbance in the channel, ${\beta }_1,{\beta }_2,{\beta }_3$ represents the observer gain, $y_x$ represents the output of the channel, and $k_p,k_d$ represent the gain parameters of the controller in the channel.

3.1.3. Y-Channel Position Controller Design

In this channel, the virtual control quantity $U_y$ is defined, and the control law for $y-\varphi$ can be designed based on Equations (14)–(17).

$$\left\{\begin{array}{c}{\dot{v}}_{y1}=v_{y2}\hfill \\ {\dot{v}}_{y2}=flan\left(v_{y1}-v_d,v_{y2},R,h\right)\hfill \\ e=z_{y1}-y_y\hfill \\ {\dot{z}}_{y1}=z_{y2}+{\beta }_{1}e\hfill \\ {\dot{z}}_{y2}=z_{y3}+{\beta }_{2}e+b_{0}u\hfill \\ {\dot{z}}_{y3}={\beta }_{3}e\hfill \\ U_y=k_p\left(v_{y1}-z_{y1}\right)+k_d\left(v_{y2}-z_{y2}\right)\hfill \\ F_y=\frac{U_{z}m+F_{z}sin\varphi cos\psi }{cos\varphi cos\psi }\hfill \end{array}\right.$$

(19)

where $z_{y1},z_{y2}$ represents the observed value of $v_{y1},v_{y2}$, $z_{y3}$ represents the observed value of the total disturbance in the channel, ${\beta }_1,{\beta }_2,{\beta }_3$ represents the observer gain, $y_y$ represents the output of the channel, and $k_p,k_d$ represent the gain parameters of the controller in the channel.

3.2. Attitude Controller

To maintain the stability of the flapping-wing aircraft’s attitude, the attitude controller computes the best control inputs at each time step using model predictive control. The three primary components of the MPC attitude controller are roll optimization, feedback correction, and the predictive model. The MPC attitude controller computes the system output in the prediction time domain by first combining the prediction model with the currently measured data. The first item in the control chain is used as the system’s real control input, and obtaining the system’s input in the control time domain involves solving the optimisation problem, taking into account the system’s constraints. Based on Equation (9), the attitude dynamics equation for the flapping-wing aircraft rotating about its centre of mass can be represented as follows:

$$\left\{\begin{array}{c}\ddot{\varphi }=\frac{I_y-I_z}{I_x}\dot{\theta }\dot{\psi }+\frac{M_x}{I_x}\hfill \\ \ddot{\theta }=\frac{I_z-I_x}{I_y}\dot{\varphi }\dot{\psi }+\frac{M_y}{I_y}\hfill \\ \ddot{\phi }=\frac{I_x-I_y}{I_z}\dot{\varphi }\dot{\theta }+\frac{M_z}{I_z}\hfill \end{array}\right.$$

(20)

By substituting the above nonlinear state equations into the state-space equations, we have the following:

$$\left\{\begin{array}{c}\dot{x}=Ax+Bu\hfill \\ y=Cx+Du\hfill \end{array}\right.$$

(21)

where $x\in {\mathbb{R}}^6$ is the state vector, $U\in {\mathbb{R}}^3$ is the input vector, and $y\in {\mathbb{R}}^3$ is the output vector. Using the Taylor expansion to expand the nonlinear expression of the system near the expected value and then using the Lagrange equation, the state matrix and control matrix of the system can be obtained as follows:

$$A=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 0& \frac{\dot{\psi }}{2}\frac{I_y-I_z}{I_x}& 0& \frac{\dot{\theta }}{2}\frac{I_y-I_z}{I_x}\\ 0& 0& 0& 1& 0& 0\\ 0& \frac{\dot{\psi }}{2}\frac{I_z-I_x}{I_y}& 0& 0& 0& \frac{\dot{\varphi }}{2}\frac{I_z-I_x}{I_y}\\ 0& 0& 0& 0& 0& 1\\ 0& \frac{\dot{\theta }}{2}\frac{I_x-I_y}{I_z}& 0& \frac{\dot{\varphi }}{2}\frac{I_x-I_y}{I_z}& 0& 0\end{array}\right],B=\left[\begin{array}{ccc}0& 0& 0\\ \frac{1}{I_x}& 0& 0\\ 0& 0& 0\\ 0& \frac{1}{I_x}& 0\\ 0& 0& 0\\ 0& 0& \frac{1}{I_x}\end{array}\right]$$

$$C=\left[\begin{array}{cccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \end{array}\right],D=\left[\begin{array}{ccc}0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill \end{array}\right]$$

Discretization of the above state-space equation gives the following discrete state equation:

$$\left\{\begin{array}{c}{x}_{k+1}=A_d{x}_k+B_d{u}_k\hfill \\ y_k=C_d{x}_k\hfill \end{array}\right.$$

(22)

where matrices $A_d=I+A{t}_s,B_d=B{t}_s$, and $C_d=C$; $t_s$ is the sampling time.

To ensure that the changes in the control input are minimized and the control action is smooth, we construct an augmented state-space equation:

Set up ${\xi }_k=\left[\begin{array}{c}{x}_k\hfill \\ u_{k-1}\hfill \end{array}\right]$.

Rewrite Equation (22) as:

$$\left\{\begin{array}{c}{\xi }_{k+1}=\tilde{A}{\xi }_k+\tilde{B}\Delta u_k\hfill \\ {\gamma }_k=\tilde{C}{\xi }_k\hfill \end{array}\right.$$

(23)

where ${\gamma }_k$ is the output of the state equation; $\Delta u_k$ is the system control increment at the current time t; $\tilde{A}=\left[\begin{array}{cc}{A}_d& B_d\\ O& I\end{array}\right]$; $\tilde{B}={\left[\begin{array}{cc}{B}_d\hfill & I\hfill \end{array}\right]}^{\mathrm{T}}$; $\tilde{C}=\left[\begin{array}{cc}{C}_d\hfill & 0\hfill \end{array}\right]$; I is the identity matrix.

The predictive model of the controller can be represented as (N is the time domain of the control input; $N_p$ is the predicted output time domain): $Y=\Psi {\xi }_k+\Theta \Delta u_k$

$$Y=\left[\begin{array}{c}\gamma (t+1\mid t)\\ \gamma (t+2\mid t)\\ ⋮\\ \gamma (t+N\mid t)\\ ⋮\\ \gamma \left(t+N_p\mid t\right)\end{array}\right],\Psi =\left[\begin{array}{c}\tilde{C}\tilde{A}\\ \tilde{C}{\tilde{A}}^2\\ ⋮\\ \tilde{C}{\tilde{A}}^N\\ ⋮\\ \tilde{C}{\tilde{A}}^{N_p}\end{array}\right],\Delta u_k=\left[\begin{array}{c}\Delta u(t\mid t)\\ \Delta u(t+1\mid t)\\ ⋮\\ \Delta u(t+N\mid t)\end{array}\right]$$

$$\Theta =\left[\begin{array}{cccc}\tilde{C}\tilde{B}& 0& \cdots & 0\\ \tilde{C}\tilde{A}\tilde{B}& \tilde{C}\tilde{B}& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮\\ \tilde{C}{\tilde{A}}^N\tilde{B}& \tilde{C}{\tilde{A}}^{N-1}\tilde{B}& \cdots & \tilde{C}\tilde{A}\tilde{B}\\ ⋮& ⋮& \cdots & ⋮\\ \tilde{C}{\tilde{A}}^{N_p-1}\tilde{B}& \tilde{C}{\tilde{A}}^{N_p-2}\tilde{B}& \cdots & \tilde{C}{\tilde{A}}^{N_p-N-1}\tilde{B}\end{array}\right]$$

The designed objective function is:

$$J=\sum _{j=1}^{N_p}{∥\gamma (t+j\mid t)-{\gamma }_d(t+j\mid t)∥}^{2}Q+\sum _{j=1}^N{\parallel \Delta u(t+j\mid t)\parallel }^{2}R$$

(24)

where $\gamma$ represents the predicted output, ${\gamma }_d$ the expected output, $\Delta u$ is the increment of the system input, and Q and R are weighted matrices. The first term in the equation represents the tracking ability for a given trajectory; the second term ensures that the control signal does not change too much.

For actual control objects, there are various constraints, and this article mainly considers the following three aspects of constraints:

$$\begin{array}{c}\Delta u_{min}\le \Delta u(t+j)\le \Delta u_{max}\hfill \\ u_{min}\le u(t+j)\le u_{max}\hfill \\ \gamma (t+N+1\mid t)={\gamma }_d(t+N+1)\hfill \end{array}$$

where $\Delta u(t+j)$ represents the control input increment; $u(t+j)$ is the control input; $\gamma (t+N+1\mid t)$ is the predicted output; ${\gamma }_d(t+N+1\mid t)$ is the expected output; $\Delta u_{min},\Delta u_{max}$ are the minimum and maximum values of the control increment, respectively; $u_{min}=u_{max}$ and represent the minimum and maximum values of the control input, respectively.

This section presents the design of the position and attitude controllers. The attitude controller employs the model predictive control approach to guarantee optimal control performance within the internal loop while continuously optimising system performance. Furthermore, the position controller is based on the linear active disturbance response controller to enhance noise immunity.

4. Experimental Simulation and Analysis

In this section, a series of typical simulation experiments were performed using the MATLAB platform to validate the trajectory tracking and disturbance rejection abilities of the flapping-wing aircraft. The aim was to verify whether the proposed control strategy could achieve the desired tracking effect. Specifically, the altitude position tracking test is carried out in a two-dimensional space, and spiral trajectory tracking was conducted in a three-dimensional space. Through experimental simulations, tracking of trajectories was obtained under both disturbed and undisturbed conditions. The simulation results are discussed in Section 3.1 and Section 3.2, respectively.

4.1. Two-Dimensional Height Position Simulation Verification

In this section, the 2D altitude tracking performance of the flapping-wing aircraft is first verified. The physical parameters of the controller parameters are shown in Table 2.

Simulation 1: The purpose of this test is to verify the control capability of the proposed control strategy. In the absence of disturbances, the input signal is a fixed height $Z_d=5$ m. In the simulation experiment, the initial height of the flapping aircraft is 0 m, and the trajectory tracking curve of the flapping aircraft is shown in Figure 3. The test results confirm that the control method has highly efficient and reliable performance for position tracking of the flapping aircraft.

4.2. Three-Dimensional Curve Trajectory Tracking Simulation Verification

To further validate the control performance of the suggested control strategy, this section conducts a simulation verification of three-dimensional curve trajectory tracking. The parameters of the flapping-wing aircraft used in the simulation experiment are the same as in the previous section. The purpose of this test is to evaluate the control performance of the proposed control strategy in three-dimensional space. In the simulation experiment, the reference trajectory of the flapping-wing aircraft is a three-dimensional space curve. The starting position of the flapping-wing aircraft is ${\left[\begin{array}{ccc}x\left(0\right)\hfill & y\left(0\right)\hfill & z\left(0\right)\hfill \end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{ccc}0\hfill & -1\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}}$, and the starting Euler angles are $\left[\begin{array}{ccc}\varphi \left(0\right)\hfill & \theta \left(0\right)\hfill & \psi \left(0\right)\hfill \end{array}\right]={\left[\begin{array}{ccc}{0}^{\circ }\hfill & 0^{\circ }\hfill & 0^{\circ }\hfill \end{array}\right]}^{\mathrm{T}}$. Additionally, external position translation disturbance $o_1$ and attitude rotation disturbance $o_2$ are introduced to simulate uncertainties in the natural environment during the flight of the flapping-wing aircraft, such as gust disturbances. The replacement of environmental gust disturbances with time-varying sinusoidal disturbances is expressed as follows:

$$\begin{array}{c}{o}_1={\left[\begin{array}{ccc}0.3sin\left(\pi \mathrm{t}\right)\hfill & 0.2sin\left(\pi \mathrm{t}\right)\hfill & 0.1sin\left(\pi \mathrm{t}\right)\hfill \end{array}\right]}^{\mathrm{T}}\mathrm{m}/{\mathrm{s}}^2\hfill \\ o_2={\left[\begin{array}{ccc}0.03sin\left(\pi \mathrm{t}\right)\hfill & 0.02sin\left(\pi \mathrm{t}\right)\hfill & 0.01sin\left(\pi \mathrm{t}\right)\hfill \end{array}\right]}^{\mathrm{T}}\mathrm{rad}/{\mathrm{s}}^2\hfill \end{array}$$

The results of the flapping-wing aircraft trajectory tracking are presented in Figure 4. This figure clearly demonstrates the superior control performance of the proposed controller. When attempting to track a desired trajectory, it takes a longer period of time for PID control to reach a state of equilibrium, and the tracking error is significantly larger. In comparison, LADRC-MPC control can track the reference trajectory more quickly and exhibits a significantly stronger tracking effect than PID control, and it is accompanied by a smaller tracking error. Furthermore, the flapping-wing aircraft is able to track the trajectory effectively in the presence of disturbances, thereby demonstrating the robust nature of LADRC-MPC control.

Figure 5 shows the position trajectories for the x-, y- and z-axes. It is clear that every time the flapping-wing aircraft rotates, there is a brief overshoot in the y-direction. However, the controller corrects this deviation in time to accurately track the reference position. It can be seen from Figure 6 that after the addition of disturbance, the tracking error of the system increases and the tracking performance decreases, but the tracking error can still be maintained within a small area, indicating that the proposed method has good anti-interference ability.

The tracking curves of the attitude angles are shown in Figure 7 and Figure 8, which show that the attitude angles remain stable during three-dimensional trajectory tracking. Comparing the results, it can be seen that when the sinusoidal perturbation is introduced, the error increases slightly compared to the case without perturbation, but the flapping-wing aircraft still maintains stable flight. The control torque is shown in Figure 9, where green represents the control input without disturbance and blue represents the control input under disturbance.

In order to enhance the authenticity of the simulation by adding sensor noise interference, the real environmental conditions encountered by flapping-wing aircraft were simulated. During the simulation process, Gaussian white noise with a standard deviation of 0.5 was used instead of sensor noise. Figure 10 clearly demonstrates the trajectory tracking performance of the flapping-wing aerial vehicle under sensor noise. After adding sensor noise disturbance, the trajectory tracking of flapping-wing aircraft shows a reduced oscillation frequency while maintaining robust performance in trajectory tracking. Figure 11 shows the position error with and without sensor noise. The robustness of trajectory tracking algorithm under sensor noise conditions was verified through simulation experiments.

This section conducted height tracking and spiral trajectory tracking tests on the flapping-wing aircraft system using the proposed controller. The test results demonstrate the good control performance of the controller.

5. Conclusions

A model predictive control (MPC) scheme based on linear active disturbance rejection control (LADRC) is proposed as a solution to the challenge of accurate trajectory tracking for flapping-wing aircraft in the presence of external disturbances. This approach integrates the predictive capabilities of MPC with the disturbance rejection features of LADRC, with the objective of enhancing the aircraft’s ability to follow a predetermined path accurately despite the influence of unpredictable environmental factors. The scheme employs LADRC to design an outer loop position controller, thereby circumventing the necessity for an accurate mathematical model in order to achieve control of the position system. The position system controller considers the coupling and external disturbances between channels as the total disturbance. The LESO is employed to estimate and compensate for disturbances. The introduction of TD enables the system to cope more effectively with external disturbances and parameter changes, thereby enhancing its overall robustness. An attitude controller based on MPC is designed to improve the controller’s robustness to disturbances and to remove tracking errors caused by parameter deviations. The results of the simulation experiments demonstrate that the control method employed in this paper is an effective means of achieving trajectory tracking. However, the limitation of this article is that it lacks physical experimental verification. The objective is to apply this scheme to practical projects and to make improvements based on experimental results. This represents a future research direction.