Modeling and Flight Control of Small UAV with Active Morphing Wings

In recent research works, morphing wings were studied as an interesting field for a small unmanned aerial vehicle (UAV). The previous studies either focused on selecting suitable material for the morphing wings or performing experimental tests on UAVs with morphing wings. Though, the dynamic modeling of active flexible morphing wings and their involved interactions with the aerodynamics of the UAV body are challenging subjects. Using such a model to control a small UAV to perform specific maneuvering is not investigated yet. The dynamic model of UAV with active morphing wings generates a multi-input multi-output (MIMO) system which rises the difficulty of the control system design. In this paper, the aeroelastic dynamic model of morphing wing activated by piezocomposite actuators is established using the finite element method and modal decomposition technique. Then, the dynamic model of the UAV is developed taking into consideration the coupling between the wing and piezocomposite actuators, as well as the dynamic properties of the morphing actuators with the aerodynamic wind disturbances. A model predictive control (MPC) is designed for the MIMO control system to perform specific flight maneuvering by tracking desired trajectories of UAV altitude and yaw angle. Additionally, the MPC achieves constrained behavior of pitch and roll angles to get satisfactory UAV motion. Also, the behaviors of the UAV control system using MPC are evaluated after adding Dryden wind turbulence to the UAV outputs. Finally, a UAV flight simulation is conducted which shows that the control system successfully rejects the applied disturbances and tracks the reference trajectories with acceptable behavior of pitch and roll angles.


Introduction
A smart structure is a well-known research topic in the mechanical engineering field, which has a long history of challenges and aspirations. Recently, aeronautic applications, especially small aircraft with morphing wings, have employed piezoelectric actuators for improving aerodynamic characteristics. Some of the required aerodynamic improvements are to reduce fuel consumption, reduce drag, and enhance flight conditions for navigation and maneuvering. The conventional piezoelectric materials have some limitations such as their very weak behavior at high mechanical stress, which leads to them cracking easily during the bonding process and operation. Moreover, its formation performance on curved surfaces is weak [1].
Therefore, a piezo-composite material is proposed to solve the challenges and problems encountered with conventional piezo ceramic material.
The macro fiber composite (MFC) piezoelectric material is widely used in smart structure applications. This MFC is composed of rectangular rods of piezoceramic material (PZT-5A) inserted between Kapton layers and interdigitated electrodes (IDEs) [2,3]. The MFC retains the merits of the previous active fiber composite (AFC), as it produces large longitudinal deformation by using reasonable voltage owing to the large value of the longitudinal piezoelectric constant, d 33 . MFC features a uniform and repeatable fabrication [4]. A full description of the fabrication process and properties of MFC can be found in [5][6][7][8].
Wang et al. [9,10] used an MFC piezoelectric actuator for active shape control of the morphing wing. The interaction among structural dynamics, unsteady aerodynamics, and MFC actuation was investigated to achieve active dynamic morphing for flexible wings. A bimorph configuration of MFCs on a flexible cantilever beam was modeled with 45° fiber orientation to obtain both bending and torsional deformations for active shape control. The detailed modeling of the MFC using the finite element method (FEM) was not presented in these works, which leaves open questions without answers.
Another work for designing and testing micro aerial vehicles using morphing wings controlled by MCF piezoelectric actuators was presented by Kochersberger et al. [11]. The authors stated that modeling MFC in simulating a morphing wing is challenging. Therefore, the authors followed a technique developed by Gustafson [12], which assumes that the piezoelectric expansion coefficient is equivalent to the thermal expansion coefficient. This thermal expansion coefficient was calculated from the piezoelectric constants of MFC.
Other trials are conducted in the field of morphing wings for small UAVs. Di Luca et al. [13] presented a design of morphing wing using bioinspired artificial feathers that can reshape the wing substrate to achieve specific aerodynamic requirements. The morphing wings were used on a small drone for rolling control. Their work was evaluated by simulations, wind tunnel tests, and outdoor flight. Similar work was conducted by Hui et al. [14] that investigates the influence of pigeon-inspired morphing wings on the aerodynamic performance of UAVs. Moreover, a load-carrying UAV was designed and manufactured with an active morphing wing/ tail to enhance flight performance [15]. The morphing was achieved by using servo motors placed in the wing. Another similar design of a morphing wing using servo motors was presented in [16]. This work presented an analysis of the morphing wing with computational fluid dynamics (CFD) simulation as well as wind tunnel tests.
Ohanian et al. [17] designed and manufactured a morphing wing actuated by MFC piezoelectrical material. They studied the deformations employed on the wing camber from the MFC actuation and its influence on the airfoil lift coefficients. Also, Han et al. [18] investigated the aerodynamic performance of a UAV with morphing winglets that mimics the wing-tips of real birds. These morphing winglets were created using a soft composite of shape memory alloy wires and glass fibers. The resulting aerodynamic performance was evaluated by wind tunnel tests for a UAV under various angles of attack.
From the above-discussed literature, no previous work was conducted on a model-based approach to control a UAV with morphing wings for flight maneuvering. The previous research works either focused on the design and manufacture of the morphing wings with their influence on the aerodynamic characteristics or perform experimental flight tests for UAVs with morphing wings without establishing the mathematical dynamic model of the UAV. The dynamic modeling of the active flexible morphing wings is challenging due to the involved interactions between the active morphing wings and the aerodynamics of the UAV body. The order of the morphing wing model is much higher than a regular fixed-wing model. So, the complexity of the plant is much considerable.
This paper focuses on studying and establishing the complex dynamic model of a UAV equipped with active morphing wings for maneuvering control and trajectory tracking. The developed model enables in-depth investigation and analysis of the dynamics of the aeroelastic morphing wings interacting with the UAV. Also, the developed model allows applying and evaluating different control strategies other than well-established control techniques for flight maneuvering tasks.
The paper is organized as follows: after reviewing the previous research works of UAVs with morphing wings controlled by piezoelectric composites in Section 1, Section 2 presents the mathematical modeling of a UAV with morphing wings actuated by MFCs. Section 3 demonstrates the control design of the UAV with morphing wings and Section 4 discusses the generated results. Lastly, Section 5 presents the conclusions of the presented work and the suggested future works.

Mathematical Modeling of Morphing Wing
This section discusses the detailed derivation of the mathematical model of a morphing wing operated by the MFC actuator. This dynamic model of morphing wing describes the characteristics of the base wing structure and the MFC actuator as well as the influence of the aerodynamic loads. The MFC piezoelectric material has a large number of IDEs which dictates using a very fine FE mesh that requires extremely large computational costs, especially with a large number of actuators. The process becomes unfeasible if simultaneous finite element simulations are needed, as in control tasks.
The authors of this paper developed a new efficient technique for modeling the MFC actuator which was previously published in [19] and [20]. This technique presented a simplified FE model for the homogenized MFC actuator with only two electrodes at its longitudinal extremes instead of the current modeling using a physically large number of electrodes, which results in a very fine FE mesh with a high computational cost. The realization that the restrictions imposed in practice should not also be restrictions in modeling is the key point in developing the conducted technique. The new voltage proposed in this technique theoretically produces the same electric field, strain, and deformation as the physical MFC.
The simulation accuracy was verified by comparing the simplified FE model with the detailed FE model, reference value in the datasheet and previously published FE simulation of a morphing wing. Furthermore, experimental tests were conducted using a cantilever beam equipped with a physical MFC actuator. Both simulations and experimental results showed good consistency for the proposed technique of modeling MFC with low computational cost. The proposed simplified FE model facilitates the modeling of smart structures with a large number of MFC actuators, which is not feasible using the previous detailed FE model. The morphing wing model is obtained using the FEM in the current work. The new technique of modeling MFC previously published [19] and [20] are used in this derivation. So, the dynamic mathematical model of the morphing wing is formulated as where M nodal , C nodal , and K nodal are the global mass, damping, and stiffness matrices respectively. x is the nodal displacements vector and represents the generalized coordinate of the system. F is the nodal forces vector. For n degrees of freedom (DOF), the size of the matrices M nodal , C nodal , and K nodal is [n × n], while the size of the vectors x, and F is [n × 1]. Rayleigh damping is considered in this model for its benefits in further modal analysis. Thus, the global damping matrix is formulated as where α d and β d are mass and stiffness damping coefficients, respectively. With the large number of DOF in the FE model represented in Eq. (1), It is difficult to use this model in the control design phase. Therefore, the modal decomposition method is used to reduce the FE model. The reduced model is achieved by using the following well-known formula where φ is the open-loop mode shapes normalized to mass, and q is the modal coordinate. For a specific number of modes (m), the size of the mode shapes matrix is [n × m], while the modal coordinate is a vector of size [m × 1]. The mode shapes are obtained from the FE model using ANSYS software.
By substituting both Eqs. (2) and (3) in Eq. (1) and multiplying by φ T : Since the mode shapes are normalized to mass, then the rules of normalization can be used as follows: where I is the identity matrix, ω is the natural frequency of each mode, and ζ is the damping ratio of each mode. Therefore, the final shape for the mathematical model of the morphing wing in the modal coordinate is: After converting the mathematical model from the generalized coordinate with large n nodal DOF to reduced m mode shapes, the control system can be designed efficiently. The nodal force vector (F) is due to the MFC actuation (F MFC ) and aerodynamic loads (F a ) as: and with The MFC force vector is calculated from the applied voltage (u MFC ) which represent the control signal from the controller, and (B MFC ) which is the effective coefficient matrix of the MFC actuator as shown in Eq. (9). For a total number of MFC actuators (N MFC ), the size of (B MFC ) matrix is [n × N MFC ] and the size of (u MFC ) vector is [N MFC × 1]. The effective coefficient matrix describes the structural electromechanical influence of the MFC actuator on the base wing [21]. Eq. (9) calculates analytically the effective coefficient of the MFC actuator, where B x is the strain-displacement matrix, (B u ) is the electric field-potential matrix, (E) is the electric field vector, and (V Volume, MFC ) is the occupied volume of the MFC layer on the base wing. According to the simplified proposed modeling of MFC with linearity assumption; this effective coefficient can be determined numerically from the static response of the FE model by applying a constant voltage of 1 V for one actuator, 0 V for the rest of the actuators, and extracting its equivalent nodal force vector (static response of Eq. (1)); then repeat this process for every MFC actuator. The modeling of aerodynamic loads (F a ) mentioned in eq. (7) will be fully discussed in the next subsection.

Modeling of Aerodynamic Loads
The aerodynamic loads in the current work are modeled as unsteady aerodynamic lift and moment. These unsteady aerodynamic loads are formulated according to the linearized thin airfoil theory [22,23] and represented as: with where Q L and Q M are the aerodynamic lift and moment about the elastic axis of the airfoil section per unit distance at span length z, respectively. h and α are the vertical bending distance and twisting angle or angle of attack (AOA) of the airfoil section, respectively. b and a are the half chord length of the morphing wing and non-dimensional distance from mid-chord to the elastic axis, respectively. ρ and V are the air density and free stream velocity, respectively. C(k) is Theodorsen's function, and k is the reduced frequency. The distribution of the aerodynamic loads along the morphing wing is modeled using strip theory in the spanwise direction. The aerodynamic variables are illustrated in Fig. 1. Some assumptions are applied to the current model which are: • The morphing wing has a thin airfoil. • The morphing wing has a rectangular cross-section with no cambers. • The morphing wing has no ailerons or flaps. • The left and right morphing wings are identical and symmetrically attached to the fuselage.
ccording to the above assumptions, the center of pressure and center of rotation (axis of rotation) are at the same point on the center of gravity of the airfoil where the elastic axis is also located. Therefore, the non-dimensional distance (a) is equal to zero in this model, and the center of gravity is located at the mid-chord of the airfoil.
The aerodynamic loads formulated using Theodorsen's function in Eqs. (10) and (11) are used in the frequency domain where harmonic excitation is assumed to be applied on the airfoil. The time-domain response from the unsteady aerodynamic loads is needed to actively control the morphing wing. Therefore, Wagner's function is used to model the unsteady aerodynamic loads, which can replace Theodorsen's function and then use Duhamel's principle in the time domain [24,25]. The Wagner's sectional circulatory lift (L c ) is with where λ(t) is Jones' [26] approximation for the Wagner function, A 1 = 0.165, A 2 = 0.335, b 1 = 0.0455, and b 2 = 0.3.
By applying Duhamel's principle to Eq. (13) The cumbersome f ( ) term in the integral can be eliminated using integration by parts The derivative of approximated Wagner function is then substitute by Eqs. (16) and (17) in Eq. (13)

Fig. 1 Morphing wing with applied aerodynamic loads
The integral parts in Eq. (18) are difficult to be integrated numerically [9,27]. So, two variables are introduced: where j = 1, 2. The time derivative of the introduced variables is then applying integration by parts By assuming zero initial conditions; f(0) = 0; then Therefore, the main Eqs. (10) and (11)  The aerodynamic loads are applied on the structural nodes of the FE model of the morphing wing at the mid-chord of the airfoil. The FE model of the morphing wing is divided into a number of strips, where the interpolation method is used to distribute the aerodynamic loads over the airfoil section for each strip. The total aerodynamic loads can be obtained by integrating the nodal forces applied on all strips. Consequently, the aerodynamic equations can be expressed by the nodal displacement, velocity, and acceleration vectors as where M a , C a , and K a are the aerodynamic effective matrices for acceleration, velocity, and displacement, respectively vector v is [2i × 1]. It is worth mentioning that the aerodynamic loads in Eqs. (10) and (11) should be multiplied by the width of each strip along the span. Also, the aerodynamic moments (Q M ) are replaced by two equal and opposite forces applied at the extreme nodes of the chord at each strip. Therefore, the aerodynamic nodal forces (F a ) are now entirely modeled.
To reduce the order of the aerodynamic model, the modal decomposition method is used again; Eq. (3); and the modal aerodynamic model is then, the modal aerodynamic model can be simplified into The modal aerodynamic Eq. (29) is substituted in the main morphing wing model in Eq. (6). Hence, the complete mathematical model of the morphing wing is which can be simplified to The small UAV is assumed to fly at a low speed of 20 m/s with low Reynold's number and below the flutter speed. The aeroelastic system is established within the flight envelope and the flight speed does not exceed the flutter speed, so the developed system is stable. Also, the flutter speed is estimated at approximately 25 m/s.

State Space Modeling of the Aeroelastic Morphing Wing
The general form of the state space model is where X, u, and y are the state vector, input or control vector, and output vector, respectively. A ss , B ss , C ss , and D ss are the dynamics matrix, input or control matrix, output or sensor matrix, and feedthrough matrix, respectively.
The complete modal aeroelastic model of the morphing wing in Eq. (32) can be represented in state-space form as and The length of the state vector X can be calculated as 2 m + 2i. From the FE model, 7 modes and 40 strips are used, so the total states are 94.
Four outputs are generated from this model, y tip , y lift , y pitch , and y roll which are defined as the maximum tip deflection, total lift force, total pitching moment, and total rolling moment of the morphing wing, respectively. The aerodynamic outputs are calculated as follows: where c and z, and dz are the wing semi-span and distance from the center of gravity of each strip to the fixation at the fuselage, and width of each strip along the spanwise direction, respectively. Moreover, the aerodynamic outputs can be represented from the FE model in modal coordinates as All other parameters in Eqs. (41), (42), and (43) are matrices relating the aerodynamic outputs to the modal coordinate q with its derivatives and the introduced aerodynamic variables v. Hence, the four outputs can be represented in state-space form as where C tip is a vector relating the modal coordinates to the tip deflection.

A Minimal Model of the Aeroelastic Morphing Wing
The equivalent model for the right and left morphing wings can be assembled in parallel as shown in Fig. 2. Each morphing wing model has 2 inputs for the MFC actuators and 3 outputs for total lift force, total pitching moment, and total rolling moment. Each wing model has 6 transfer functions corresponding to each output with each input. Therefore, the state-space model of the aeroelastic morphing wing has a very large order with 564 states, 4 inputs, and 3 outputs. The indicated order of the morphing wing model is much higher than a regular fixed wing model due to a large number of elastic modes. These modes are required to describe the dynamics of the flexible wings which are actively shaped by MFC actuators. This large order model consumes high computational efforts, especially during repetitive simulations in the control design stage. Hence model order reduction is needed to reduce computational costs. The Hankel singular values and approximation error are plotted in Fig. 3 their frequency responses are plotted in Fig. 4, which shows great consistency.

Static and Dynamic Response of the Morphing Wing Model Due to Step Inputs
The  Fig. 5a are fixed support which represents the fuselage, ground, and active voltages for each MFC actuator. The static deformed shape of the morphing wing resulted from applying step inputs of V 1 = − 500 V and V 2 = 500 V without considering aerodynamic loads is demonstrated in Fig. 5b. Whereas, Fig. 5c shows the static deformation of the morphing wing after applying aerodynamic loads for the same voltages at the MFC actuators. The resulting aerostatic lift force is 3.0877 × 10 −1 N, the aerostatic pitching moment is 3.242 × 10 −3 N.m, and the aerostatic rolling moment is 4.512 × 10 −2 N.m.
Referring to the equations of the aerodynamic model in subsection 2.2 the parameters of the FE model have the following properties in Table 1.
The reduced morphing wing model was implemented on MATLAB and the responses of the model due to applying step inputs of voltages V 1 = − 500 V and V 2 = 500 V are  Fig. 6. The percentage differences between the steady-state response (final value) of the reduced (modal) model and the static response of the FE model are −3.7% for the maximum tip deflection, −6.5% for the total lift force, −6.5% for the total pitching moment, and − 5% for the rolling moment.

Mathematical Model of Small UAV with Morphing Wings
After modeling the aeroelastic morphing wing influenced by aerodynamics loads in the previous subsections, the complete model of a small-light UAV will be created. This complete model consists of the UAV body model including the fuselage, wings, and tail combined with the previous model of the morphing wing. A CAD sketch in Fig. 7 demonstrates the small UAV with two morphing wings equipped with two MFCs for each wing and two horizontal tails. The UAV is assumed to be in longitudinal equilibrium, where it moves horizontally with a constant velocity of 20 m/s and its equivalent weight is applied at the center of gravity of the UAV body. The overall UAV model is developed with body mass (m) and moments of inertia (I xx , I yy , I zz ). The net mass for the UAV is estimated as 0.2 Kg, while its moments of inertia are estimated based on the statistical method in [28,29]. This method supports the assumption of similarity for the non-dimensional radii of gyration between aircraft with nearly similar configurations and orientations. The moments of inertia for the current UAV can be calculated from the Eqs. (47), (48), and (49). Where c and L are wing semi-span and overall length of the UAV, respectively. R x , R y , and R z are the non-dimensional radii of gyration of aircraft influencing rolling moment, pitch moment, and yaw moment, respectively.
The UAV is a single-engine propeller-driven aircraft with a span (2c) of 400 mm and a total length (L) of 290 mm. The radii of gyration for the current design can be approximately referred to as the airplane type "Cessna 177A" which can be extracted from Table B1 in Appendix B [29]. So, the values of the radii of gyration are: R x = 0.212 , R y = 0.362 , and R z = 0.394.
The control objectives in the current work are to achieve desired altitude and yaw angle with maintaining acceptable behavior of the pitch and roll angles of the aircraft. Three outputs of the aeroelastic morphing wing model are considered, which are the lift force, pitching moment, and the rolling moment. These three outputs of the aeroelastic where, h UAV , α UAV , θ UAV , and ψ UAV are the altitude, pitch angle, roll angle, and yaw angle of the UAV body, respectively. F L, Total , M P, Total , M R, Total and M ψ, Total are the total lift force, total pitching moment, total rolling moment, and total yawing moment applied on the UAV body, respectively. The total lift force, total pitching moment, total rolling moment, and total yawing moment are calculated as: Where F Wing , M α, Wing , and M θ, Wing are the lift force, pitching moment, and rolling moment generated from the aeroelastic deformations on the morphing wings due to MFCs actuation and the change of the UAV pitch angle as well as the initial pitch angle (4°), respectively. M ψ, Wing is the yawing moment due to the roll rate generated from the wing contribution.   and roll angle, respectively. W is the equivalent weight of the UAV body. F Tail , and M α, Tail are the lift force and pitching moment generated from the horizontal tail due to the change of the UAV pitch angle as well as the initial pitch angle, respectively. It is to be noted that, Eqs. (23) and (24) are used to calculate aerodynamic lift force and pitch moment for both wing and horizontal tail. The half chord and semi-span length for the horizontal tail are assumed to be 7 mm and 60 mm, respectively. The yawing moment due to roll rate is calculated as: Where, C L is the equilibrium lift coefficient due to the initial AOA; and based on the thin airfoil theory; it can be calculated as C L = 2πα o . C D is the induced drag coefficient which can be approximately calculated as C D,UAV = 4 ( UAV ) 2 0.7AR and AR is the aspect ratio that equals the total wing span divided by the wing chord length.
The drag force and moments can be calculated as [30]: Where ρ is the air density, A D is the exposed area, and C D is the drag coefficient which equals 2 for the rectangular cross-section. w, t UAV , and L is the fuselage width, thickness, and length, respectively. b and c are the half chord and semispan length of the wing, respectively. The nonlinear equations of the drag force and moments are linearized about an average operating velocity of the UAV altitude, pitch angle, roll angle, and yaw angle. The complete model of the UAV is achieved by assembling its body dynamic model with the aeroelastic morphing wings model, dynamic model of the wings, dynamic model of the tail, and linear drag model. Also, the static models of the UAV weight, aerodynamics from the initial pitch angle, and constant parts from the linearization of the drag equations are added to the complete model of the UAV body. A Simulink block diagram for the complete model of UAV is demonstrated in Fig. 8. Consequently, the complete model of the UAV is a multi-input, multi-output model (MIMO). This model is called the plant model (G p ) in the control system design.

Sequential Steps for the Proposed Method
The sequence of steps for the conducted method could be declared by the flowchart in Fig. 9. First, the active morphing wing is developed by FE modeling on ANSYS software. Then, modal analysis is conducted to generate the mode shapes and the modal frequencies. Next, the modal decomposition method is applied, and its parameters are extracted from ANSYS. These modal parameters are imported to MATLAB software to develop the dynamic model of the active morphing wing. Then, the dynamic model of a UAV with morphing wings is developed. After that, a control system is developed to control this UAV for maneuvering and trajectory tracking.

Design of Model Predictive Control
A closed control loop system was constructed to control the UAV to perform certain flight maneuvering by tracking desired trajectory path of altitude and yaw angle with acceptable constrained behavior of pitch and roll angles. A model predictive controller (MPC) was used in this closedloop system to achieve the requirements. The MPC has two input measured variables (UAV altitude and yaw angle), two corresponding input references, and four output manipulated The basic idea of MPC is to predict the future behavior of the plant over a finite time horizon and compute a corresponding optimal control input while satisfying specific system constraints. The cost function is formulated so that the system output tracks a given reference for a prediction horizon. So, the prediction horizon should be large enough to represent the effect of changes in the manipulated variables on the control variables. However, increasing the prediction horizon needs more computational efforts.
The control horizon is the number of moves for the manipulated variables to be optimized at the control interval. The control horizon falls between 1 and the prediction horizon. So, increasing the control horizon generally enhances the control behavior but promotes slow computations. Therefore, precise tuning for the sampling time, prediction horizon, and control horizon is required.
After tuning for the MPC parameters according to guidance procedures from the MATLAB documentary, the MPC was set to a sample time of 0.001 seconds, prediction horizon of 300, and control horizon of 10. Its manipulated variables have constraints from −500 V to +500 V and rate weights of 0.001. The references for the UAV altitude and yaw angle are defined as quintic polynomial paths. The closed-loop control system was implemented on Simulink as shown in Fig. 11. The MPC reacts only to the dynamic statespace model of the UAV plant and so, the effect of the static forces inside the UAV model disturbs the control process. To solve this problem, static nominal values are added to the manipulated variables so that they can compensate for the static forces in the plant model. These static values of the manipulated variables are calculated from where, C Lift , A, and B are the matrices of the morphing wing state-space model identified in subsection 2.3. The notation " †" denotes the Pseudoinverse which is also called the Moore Penrose inverse. F static are the total static forces in the UAV model divided by 2 for the two wings, which include UAV weight, initial lift force due to initial pitch angle, and the constant part in the linear model of the drag force. The static values are found to be about −284.5 V for u 1 and u 3 , and about 284.5 V for u 2 and u 4 .

Dryden Wind Turbulence Model
To realize the actual behavior of the UAV during flight, wind turbulence should be added to the UAV model and the control design should consider and reject this disturbance during the reference tracking process. Discrete wind turbulence with Dryden velocity spectra was used in the current study [31,32]. This disturbance model creates an average disturbance for all air turbulence conditions based on the UAV altitude, longitudinal speed, and body orientations. Also, the wind disturbance model is driven by band-limited white noise with digital filter finite difference equations. The disturbance model generates turbulence velocities (u g , v g , and w g ) and turbulence angular rates (p g , q g , and r g ) in the longitudinal (x-axis), lateral (y-axis), and vertical (z-axis) directions, respectively. These The turbulence can be defined as a stochastic process with velocity spectra functions as follows: Where, L u , L v , and L w are the turbulence scale lengths, σ u , σ v , and σ w are the turbulence intensities in the x-axis, y-axis, and z-axis, respectively. ω g is the circular frequency in rad/s, which is calculated by multiplying the UAV speed (V) by the turbulence spatial frequency Ω in rad/m.
For the Low altitude disturbance model: where, h UAV is the UAV altitude and V 6m is the wind speed at 6 m which is assumed to be 0.7 m/s.

Results and Discussion
Simulations on MATLAB software are conducted to evaluate the developed control system for the small UAV with active morphing wings. The control objectives are to track quintic reference trajectories for the UAV altitude and yaw angle. The UAV should take off from the ground to reach

Fig. 11
Simulink block diagram for the whole closed-loop control system a specific height of 20 m. Also, the UAV should achieve reference yaw angles with a specific sequence starting initially from 0° to 50° at a time of 12 sec, then decreasing to 40° at a time of 16 sec, decreasing again to 30° at a time of 18 sec, then increasing to 45° at time 24 seconds. The UAV altitude and the sequence of the yaw angles agree with specific maneuvering for the UAV in a simulated open-world environment.
The UAV successfully performs certain flight maneuvering tasks by tracking desired quintic trajectories of the UAV altitude and yaw angle as shown in Fig. 12a and b. The quintic trajectories represent more realistic motions than the step trajectories for the UAV during maneuvering. The corresponding error signals for the UAV altitude and yaw angle are shown in Fig. 12c and d. These results show that the UAV tracks well the reference trajectories with a maximum error of about -0.2587 m occurring at about 7.9 m of UAV altitude and final value error of 0.05594 m, while the maximum error for the yaw angle is -0.69° occurred at about 25.86° and final value error of -0.00292°. Also, Fig. 12e shows the UAV pitch angle starts with an initial angle of 4° and oscillates within an acceptable constrained range until it saturates at about 2.4°. The response of the UAV roll angle during maneuvering is shown in Fig. 12f which is within an acceptable constrained range. Moreover, the control signals (manipulated variables) generated from the MPC are shown in Fig. 12g, which follow the piezo voltage constraints from −500 V to +500 V.
The behaviors of the UAV control system using MPC are evaluated again after adding wind turbulence to the UAV outputs. The response of the discrete Dryden wind turbulence for UAV altitude, pitch angle, roll angle, and yaw angle are displayed in Fig. 13a, b, c and d, respectively. The amplitude of the altitude disturbance increases, as expected, with increasing the UAV altitude. Also, the variations in the yaw angle disturbance are high due to the applied control actions from the MPC on the controlled yaw angle. The initial high variations in the pitch disturbance are due to the initial variations of the uncontrolled pitch angle. The UAV control system successfully rejects the applied disturbances and tracks the reference trajectories to perform the flight maneuvering task, which is demonstrated through the UAV responses in Fig. 14a and b for the UAV altitude and yaw angle, respectively. Also, the error signals for the UAV altitude and yaw angle are shown in Fig. 14c and d, respectively. These results show that the UAV tracks well the reference trajectories with a maximum error of about -0.264 m occurring at about 7.8 m of UAV altitude and a final value error of 0.2577 m, while the maximum error for the yaw angle is -0.71° occurred at about 25° and final value error of -0.00275°. The maximum errors for the UAV altitude and yaw angle after adding the wind disturbances are larger than their corresponding errors without the wind influence. The effect of the wind disturbance on the pitch and roll angles are illustrated in Fig. 14e and f, which are still within acceptable constraint ranges. Moreover, the control signals after adding the wind disturbances are illustrated in Fig. 14g, which are still within the voltage constraints of the MFC actuators.
After investigating the MPC performance, the controller successfully copes with the challenges of simultaneously tracking the UAV altitude and yaw angle trajectories by controlling only the morphing wings. The control task of achieving acceptable constraint behaviors for the pitch and roll angles added more challenges to the MPC performance.
An open-world simulated environment is constructed to validate the UAV for performing specific maneuvering tasks. The UAV flight maneuvering is demonstrated through the screenshots acquired from the UAV animation as shown in Fig. 15. The flight maneuvering task is performed successfully by the UAV with active morphing wings under applied turbulences.
The limitations of the conducted work could be pointed to the voltage constraints of the MFC actuators which are limited to the applied voltage from +1500 V to -500 V. In this paper, the generated voltages from the controller to the MFC actuators are constrained from −500 V to +500 V. Also, the morphing wings are considered to have flexible thin airfoil, so that the MFC actuators could be able to actively change the wing's shape. Furthermore, the conducted control system is based on a linear mathematical model developed from the dynamics of the UAV with the active morphing wings. The nonlinearities already existed in the conducted system are overcome by linearization about operating conditions as well as assuming constant horizontal flying speed for the UAV. In addition, the flight stability for the UAV's uncontrolled pitch and roll angles were considered by achieving a constraint behavior for them during the UAV maneuvering. To include the roll and pitch angles in the control strategy for trajectory tracking, the horizontal and vertical tails should be also controlled which are considered for future works.

Conclusions and Future Works
The aeroelastic morphing wing was modeled using a modal decomposition approach acquired from its finite element model. This model takes into consideration the unsteady aerodynamic loads developed in the time domain. Each wing model has two MFC actuators developed using the previously published technique of FE modeling of MFC. Then a MIMO dynamic model of a small UAV with active morphing wings controlled by the MFC actuators was established. The wings model has large order which is reduced to its minimal order model after checking its Henkel plot. Hence, a closed-loop control system was designed using a model predictive controller (MPC) to perform specific flight maneuvering by tracking quintic trajectories of UAV altitude and yaw angle. The results show that the UAV successfully tracks the reference trajectories with an acceptable constrained response of UAV pitch and roll angles. Moreover, Dryden wind turbulences were added to the controlled outputs, where the MPC succeeded to reject these disturbances and track the reference trajectories.
The conducted work in the paper can be further progressed as follows: • Adding morphing tails actuated by MFC actuators to the overall UAV model and considering it in the control design. • Designing and manufacturing of a small UAV with morphing wings/tails controlled by MFC actuators. Funding Open access funding provided by The Science, Technology & Innovation Funding Authority (STDF) in cooperation with The Egyptian Knowledge Bank (EKB). The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Data Availability
The data presented in this study are available on request from the corresponding author.

Declarations
Ethical Approval The article has the approval of all the authors.

Consent to Participate
All the authors gave their consent to participate in this article.

Consent to Publish
The authors gave their authorization for the publishing of this article.

Competing Interests
There is no conflict of interest or competing interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.