Self-Scheduled LPV Control of Asymmetric Variable-Span Morphing UAV

Abstract

In this study, a novel framework for the flight control of a morphing unmanned aerial vehicle (UAV) based on linear parameter-varying (LPV) methods is proposed. A high-fidelity nonlinear model and LPV model of an asymmetric variable-span morphing UAV were obtained using the NASA generic transport model. The left and right wing span variation ratios were decomposed into symmetric and asymmetric morphing parameters, which were then used as the scheduling parameter and the control input, respectively. LPV-based control augmentation systems were designed to track the normal acceleration, angle of sideslip, and roll rate commands. The span morphing strategy was investigated considering the effects of morphing on various factors to aid the intended maneuver. Autopilots were designed using LPV methods to track commands for airspeed, altitude, angle of sideslip, and roll angle. A nonlinear guidance law was coupled with the autopilots for three-dimensional trajectory tracking. A numerical simulation was performed to demonstrate the effectiveness of the proposed scheme.

1. Introduction

An advanced aerial platform known as a morphing unmanned aerial vehicle (UAV) is capable of controlled, large-scale shape transformations while in flight, which improves efficiency, adaptability, and performance across a range of flight scenarios and missions. Modern UAV missions are becoming increasingly complex [1,2], making it difficult for fixed-shape aircraft to perform their assigned tasks effectively and efficiently. Consequently, the demand for adaptable, morphing aircraft has gradually increased. Since the first manned flight, aircraft designers have been intrigued by the concept of variable-geometry aircraft. Previously, the ability to change shape was used only to modify flight characteristics. However, recent advancements in materials, sensors, and actuators have sparked interest in the development of morphing aircraft [3]. In the development process, complex trade-off analyses and multidisciplinary optimization are frequently employed. Effective actuator drive systems also permit morphing for direct flight control. As new concepts for morphing aircraft emerge, the demand for efficient control methods is increasing [4]. However, control system design frameworks continue to fail to maximize all possible benefits of morphing.

Conventional fixed-wing aircraft are designed for optimal flight performance within the intended flight conditions and specified missions at the expense of performance in other flight regimes. In contrast, the ability to change configurations allows the morphing aircraft to fly nearly optimally in a broader range of flight conditions. In particular, morphing wings can be used to generate aerodynamic control forces by controlling the airflow directly around a vehicle. This additional degree of freedom can be used (i) to accomplish a movement that is normally impossible with conventional control surfaces alone, (ii) to increase agility and maneuverability when paired with conventional control surfaces, or (iii) to partially replace conventional control surfaces in the event of actuator failure. During flight, morphing configurations can be optimized in real time to maximize system-level benefits such as maneuverability, energy efficiency, survivability, etc. Numerous studies have explored the benefits of morphing, but few have focused on the control system that enables its full use.

It is difficult to use morphing parameters as control inputs in designing the control system of a morphing aircraft for two primary reasons. First, conventional control surfaces generate control torque for the roll, pitch, and yaw axes using long moment arms and minimal aerodynamic forces. However, wing morphing frequently generates forces and moments along multiple axes simultaneously. In addition, there is typically a significant nonlinear relationship between morphing parameters and the aerodynamic forces and torques that are generated. In a large-scale shape change, the location of the center of mass (CM), the moment of inertia (MOI), and the product of inertia (POI) can also vary considerably. Second, the majority of morphing actuators have a slower dynamic response to the command than conventional hydraulically actuated control surfaces. The low bandwidth of morphing actuators may significantly diminish the overall performance of control systems. Using the morphing parameters as control inputs necessitates careful consideration of the bandwidth gap.

If morphing parameters are considered exogenous parameters, the designed controller should be able to change its characteristics according to the changing morphing parameters. Since the dynamic characteristics change greatly depending on the morphing shape, it is generally difficult to ensure a sufficient level of stability and performance with a fixed controller. For example, control systems with traditional gain scheduling techniques cannot generally guarantee closed-loop stability when the scheduling parameter changes rapidly. In addition, interpolation or blending between point controllers involves a trial-and-error procedure with little theoretical guidance. Furthermore, in the region of transition between design points, robustness and performance guarantees of specific operating points are lost. Nevertheless, assuming that the parameters change slowly nullifies the fundamental advantages of morphing.

Historically, gain scheduling has been a common design technique for nonlinear flight control system design. A gain self-scheduled linear parameter-varying (LPV) control method has recently attracted interest as a suitable control method for morphing aircraft. Consequently, the majority of early studies employing LPV architecture involved flight control systems. The LPV framework enables the flight controller of a morphing aircraft to be designed with theoretical guarantees of robustness and performance across a wide variety of operating conditions and configurations. Incorporating the nonlinearity of the plant into scheduling parameters is also an effective method for addressing this nonlinearity. Because the LPV framework shares the basic control interconnection structure with a standard multivariable controller, such as a ${\mathcal{H}}_{\infty }$ controller, different actuator bandwidths can be considered during the design process. In addition, the stability guarantee for arbitrarily fast parameter variations allows the morphing parameters to be used freely to the maximum without compromising stability.

The Wright Flyer, the first powered airplane, can be considered a morphing aircraft because it deforms its wings without a separate control surface [5]. Through the 1980s, various types of variable-geometry aircraft were developed primarily for military purposes. Then, beginning in the 1990s, modern morphing aircraft emerged, and morphing technology has recently been applied to small unmanned aerial vehicles. This trend is primarily attributable to the development of new materials, actuators, and sensors, as well as the fact that modern aircraft must perform significantly more complicated tasks than in the past [6].

Classical gain scheduling techniques have been widely adopted for systems with a wide range of operations. Control design methods based on linear time-invariant (LTI) models such as classical control, loop-shaping, ${\mathcal{H}}_2$ (LQG) optimal control, ${\mathcal{H}}_{\infty }$ suboptimal control, mixed-sensitivity control, and mu-synthesis can be applied. Various linear and classical control methods have been applied to control morphing aircraft, including proportional–integral–derivative (PID) control [7]. Many studies have adopted conventional gain scheduling techniques or switching methods [8,9]. Learning-based or data-driven approaches (reinforcement learning, neural-network-based adaptive control) [10,11,12] and adaptive approaches [12,13,14] have also been adopted to address the time-varying nature of morphing aircraft. Control allocation methods have been adopted to address an overactuated model of morphing aircraft from incorporating morphing parameters into a control input [15]. Various nonlinear control methods have also been applied, including sliding mode control [16,17], nonlinear dynamic inversion [18], backstepping control [19], and disturbance-observer-based control [20,21,22]. LPV and robust control methods are also extensively applied [23]. Many LPV techniques, such as switching LPV control, model predictive control based on the LPV model, data-driven strategies, and LPV control with scheduling uncertain parameters, have been applied to aerospace domain problems [24,25,26,27,28]. Fuzzy gain scheduling techniques have also been studied [29,30].

Remark 1.

When comparing to other methods such as gain-scheduled ${\mathcal{H}}_{\infty }$ control, nonlinear dynamic inversion control, and learning-based control, the LPV control method adopted in this study provides a rigorous stability guarantee under arbitrarily fast morphing, which requires moderate online computation for practical implementation. Note that, for an LPV controller to work properly, the LPV model must be able to successfully mimic the nonlinear model.

In this study, LTI models of longitudinal and lateral-directional dynamics were obtained by linearizing the nonlinear model across all flight conditions and shape transformations. Interpolation was used to parameterize the family of pointwise LTI models to obtain the LPV models. The symmetric component of span variation was considered as a scheduling parameter alongside airspeed and altitude, while the asymmetric component was considered as a lateral-directional axis control input alongside ailerons and rudders. CAS and autopilots were designed for manual and automatic flight, respectively. The symmetric morphing parameters were used to assist the morphing aircraft in executing the commanded maneuver and improving the flight performance across a variety of flight conditions, while the asymmetric morphing parameter served as an additional control input for the lateral-directional axis. Numerical simulation of various maneuvers demonstrates that the maneuverability of a morphing aircraft can be improved by adjusting the symmetric morphing parameter appropriately. In addition, the proposed control system successfully follows the flight trajectory under rapid variation in flight conditions and extremely fast morphing, whereas the baseline controller with the same parameter settings rapidly loses closed-loop stability.

This study’s contribution can be divided into three parts. First, a high-fidelity model of an asymmetric variable-span morphing UAV was developed, from which, nonlinear and LPV morphing aircraft models were derived. The nominal model is presented first, followed by the damage model, from which, the morphing model was derived. Through Jacobian linearization, a pointwise LTI model was derived from the trim condition. By introducing scheduling parameters, an LPV model was constructed, which was then analyzed in the frequency domain and time domain to gain physical insight into the controller design. Second, the design of the control augmentation system is based on LPV methods that utilize variable span morphing parameters as the control input and scheduling parameters. Based on LPV techniques, a control augmentation system was designed for the morphing-assisted maneuvers of morphing UAVs. A normal acceleration CAS was used for controlling longitudinal motion, whereas a roll rate CAS was used for controlling lateral motion. The CAS employed symmetric and asymmetric morphing to enhance agility and maneuverability, respectively. The effectiveness of the proposed scheme is demonstrated through numerical simulations of push-over and pull-up maneuvers and a high-g turn. Finally, an autopilot for the morphing-assisted flight of morphing UAVs based on LPV methods was designed. Autopilots for airspeed and altitude were designed to control longitudinal motion. A roll angle autopilot was designed for controlling lateral motion. The autopilot employed symmetric and asymmetric morphing to enhance agility and maneuverability, respectively. The effectiveness of the proposed scheme is demonstrated using numerical simulations.

• A high-fidelity model of an asymmetric variable-span morphing UAV was developed, from which, nonlinear and LPV morphing aircraft models were derived.
• The design of the control augmentation system is based on LPV methods that utilize variable span morphing parameters as the control input and scheduling parameters.
• An autopilot for morphing-assisted flight of morphing UAVs based on LPV methods was designed.

Remark 2.

This paper is focused on application rather than the theory itself. Many practical considerations are provided in the manuscript throughout the LPV modeling and control design procedure, which is vital in the control system design for morphing aircraft. In particular, the capability of symmetric and asymmetric morphing is fully exploited to maximize the performance of the morphing aircraft. Moreover, it is demonstrated in the simulation results that morphing can be used to either assist the maneuver or flight control.

This paper is organized as follows. In Section 2, the LPV model of morphing aircraft is derived. In Section 3, LPV-based CAS is designed, and in Section 4, LPV-based autopilot is designed. In Section 5, a summary of the main results of this study and suggestions for future work are provided.

2. Asymmetric Variable-Span Morphing UAV Model

The NASA GTM shown in Figure 1 was selected as the baseline model for the development of a morphing UAV model in this study. The GTM simulation software is accessible via an open-source agreement with NASA [31]. GTM’s high-fidelity models of dynamics, aerodynamics, sensors, and actuators are ideal for not only analyzing the effects of morphing in flight but also demonstrating that the proposed technique can achieve the desired performance in a real system.

In the GTM, there are six damage cases, so a nonlinear variable-span morphing model can be derived for the scenario where 25% of the left wingtip is lost. It is possible to model the effects of variable-span morphing by assuming that the wingtips can telescope continuously rather than being lost. Therefore, the mass, CM, MOI, POI, and aerodynamic coefficients depend on the morphing parameters. It is possible to generate LTI models that decouple longitudinal and lateral-direction motions by linearizing the acquired nonlinear parameter variation model of the morphing UAV at each operational point. By examining the LTI model in the time and frequency domains, it is possible to acquire the necessary knowledge for controller design. By interpolating the LTI model family, the LPV model can be constructed. It can be demonstrated that the LPV model adequately captures the nonlinearity of the original system.

2.1. Nonlinear Model of a Morphing UAV

A nonlinear model of a morphing UAV was derived. First, the motion equations of the nominal model, aerodynamic model, sensor model, and actuator model are described in detail. In the case of asymmetric variable-span morphing, the CM, inertia matrix, and aerodynamic coefficients were then obtained as functions of the morphing parameters.

2.1.1. Nominal Model of a Baseline Model

Before presenting the morphing aircraft model, the construction process for the baseline model is described. The motion equations of a morphing UAV are identical to those of a conventional fixed-wing aircraft. The following states were utilized directly in numerical simulations.

Flat-Earth and constant gravitational acceleration were assumed in navigation equations. The geographic coordinate system was used to generate the vehicle’s position in the north-east-down (NED) coordinate system based on the world geodetic system 1984 (WGS84), which was used only for visualizing the vehicle’s spatial trajectory. The geodetic coordinate navigation equations can be represented as

$$\begin{array}{cc}\hfill \dot{\phi }& ={\displaystyle \frac{V_N}{R_0+h}}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \dot{\lambda }& ={\displaystyle \frac{V_{E}secL}{R_0+h}}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \dot{h}& =-V_D\hfill \end{array}$$

(3)

$$\left[\begin{array}{c}{V}_N\\ V_E\\ V_D\end{array}\right]=\left[\begin{array}{ccccc}c\theta c\psi & & -c\varphi s\psi +s\varphi s\theta c\psi & & s\varphi s\psi +c\varphi s\theta c\psi \\ c\theta s\psi & & c\varphi c\psi +s\varphi s\theta s\psi & & -s\varphi c\psi +c\varphi s\theta s\psi \\ -s\theta & & s\varphi c\theta & & c\varphi c\theta \end{array}\right]\left[\begin{array}{c}U\\ V\\ W\end{array}\right]$$

(4)

where $\phi$ is the geodetic latitude, $\lambda$ is the longitude, h is the height above ground, $R_0$ is the Earth’s radius, $V_N$, $V_E$, and $V_D$ are the north, east, and down components of the velocity vector in the local geographic coordinate system, respectively, U, V, and W are the x-axis, y-axis, and z-axis components of the velocity vector in the body-fixed coordinate system, respectively, $\varphi$, $\theta$, and $\psi$ are the Euler angles (roll, pitch, and yaw angles) of the vehicle body axes relative to the geographic system, respectively, and s and c stand for the sine and cosine, respectively. The forward-right-down system was utilized for the body-fixed coordinate system, whose origin is fixed at the vehicle’s CM.

Assuming constant gravitational acceleration, the force equations can be represented as [32]

$$\begin{array}{cc}\hfill \dot{U}& =RV-QW+{\displaystyle \frac{X_G+X_A+X_T}{\mathrm{m}}}\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \dot{V}& =-RU+PW+{\displaystyle \frac{Y_G+Y_A+Y_T}{\mathrm{m}}}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \dot{W}& =QU-PV+{\displaystyle \frac{Z_G+Z_A+Z_T}{\mathrm{m}}}\hfill \end{array}$$

(7)

where

$$\begin{array}{cc}\hfill X_G& =-gsin\theta \hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill Y_G& =gsin\varphi cos\theta \hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill Z_G& =gcos\varphi cos\theta \hfill \end{array}$$

(10)

and P, Q, and R are the x-axis, y-axis, and z-axis components of the angular velocity vector in the body-fixed coordinate system, respectively, $X_G$, $X_A$, $X_T$, $Y_G$, $Y_A$, $Y_T$, $Z_G$, $Z_A$, and $Z_T$ are the x-axis, y-axis, and z-axis components of the gravity, aerodynamic force, and engine thrust vectors in the body-fixed coordinate system, respectively, g is the gravitational acceleration, and m is the vehicle mass.

The kinematic equations of the Euler angle dynamics can be represented as

$$\begin{array}{cc}\hfill \dot{\varphi }& =P+tan\theta (Qsin\varphi +Rcos\varphi )\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \dot{\theta }& =Qcos\varphi -Rsin\varphi \hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \dot{\psi }& ={\displaystyle \frac{(Qsin\varphi +Rcos\varphi )}{cos\theta }}\hfill \end{array}$$

(13)

The moment equations can be represented as

$$\left[\begin{array}{c}P\\ Q\\ R\end{array}\right]=J^{-1}\left(\left[\begin{array}{c}{l}_A+l_T\\ m_A+m_T\\ n_A+n_T\end{array}\right]-\left[\begin{array}{c}P\\ Q\\ R\end{array}\right]\times \left(J\left[\begin{array}{c}P\\ Q\\ R\end{array}\right]\right)\right)$$

(14)

where $l_A$, $l_T$, $m_A$, $m_T$, $n_A$, and $n_T$ are the x-axis, y-axis, and z-axis components of the aerodynamic torque vector and the torque vector generated by the engine about the CM in the body-fixed coordinate system, respectively, and J is the inertia matrix of the vehicle with respect to the body-fixed coordinate system. Because the motion equations are formulated with respect to the CM, gravitational moments are always zero.

NASA GTM aerodynamic data were utilized directly, with the aerodynamic coefficients provided as data lookup tables of aerodynamic angles. Dependency on altitude is accounted for by the atmospheric density term, which is derived from the 1976 US standard atmosphere. The dynamic pressure term includes airspeed dependence. The GTM is a medium-sized (57.71 lb) unmanned aerial vehicle that flies in the subsonic region (Mach 0.09–0.18), with compressibility effects perhaps being negligible below Mach 0.3. Therefore, it may be claimed that the GTM is capable of simulating with appropriate accuracy.

The aerodynamic force components can be represented as

$$\begin{array}{cc}\hfill X_A& =\overline{q}S{C}_X\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill Y_A& =\overline{q}S{C}_Y\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill Z_A& =\overline{q}S{C}_Z\hfill \end{array}$$

(17)

where $\overline{q}=0.5\rho V_T^2$ is the dynamic pressure, S is the planform area of the main wing, $\rho$ is the atmospheric density, $V_T$ is the true airspeed (TAS), and $C_X$, $C_Y$, and $C_Z$ are body-axis aerodynamic force coefficients.

The aerodynamic torque components can be represented as

$$\left[\begin{array}{c}{l}_A\\ m_A\\ n_A\end{array}\right]=\left[\begin{array}{c}\overline{q}Sb{C}_l\\ \overline{q}S\overline{c}{C}_m\\ \overline{q}Sb{C}_n\end{array}\right]+({\mathbf{r}}_{cp}-{\mathbf{r}}_{cm})\times \left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\end{array}\right]$$

(18)

where b is the wing span, $\overline{c}$ is the mean aerodynamic chord (MAC), $C_l$, $C_m$, and $C_n$ are body-axis aerodynamic moment coefficients, respectively, and ${\mathbf{r}}_{cp}$ and ${\mathbf{r}}_{cm}$ are the position vectors of the center of pressure (25% MAC) and the CM (21.99% MAC) in the ARS.

The coefficient vector $\mathbf{C}={\left[C_X{C}_Y{C}_Z{C}_l{C}_m{C}_n\right]}^T$ consists of the basic airframe components ${\mathbf{C}}_{\mathrm{ba}}$, control surface components ${\mathbf{C}}_{\mathrm{cs}}$, dynamic derivative components ${\mathbf{C}}_{\mathrm{dd}}$, and stall-rolling moment asymmetry components ${\mathbf{C}}_{\mathrm{rma}}$. The basic airframe component is given in the form of a data lookup table that depends on the angle of attack (AOA) and the angle of sideslip (AOS). The control surface components also depend on the corresponding control surface deflections. The aerodynamic coefficients are built up as [31]

$$\mathbf{C}={\mathbf{C}}_{\mathrm{ba}}+\Delta {\mathbf{C}}_{\mathrm{cs}}+\Delta {\mathbf{C}}_{\mathrm{dd}}+\Delta {\mathbf{C}}_{\mathrm{rma}}$$

(19)

Five sensor systems are included in the GTM.

• Air data system (ADS)
• MIDG II GPS/INS (MIDG)
• MAG3 6-DOF analog IMU with a triaxial magnetometer, accelerometer, and gyro (MAG3)
• Surface potentiometers (Pots)
• Engine control unit (ECU)

The ADS provides the barometric altitude, AOA, AOS, dynamic pressure, TAS, and temperature ratio. At a sampling rate of 50 Hz, the MIDG provides Euler angles, geodetic latitude, longitude, altitude, and geographical velocity. MAG3 offers 3-axis acceleration and angular velocity. The pots provide control surface deflections, and the ECU provides data on the left and right engines’ internal operations. Sensor measurements are corrupted by appropriate dynamics, scale factor error, bias, white noise, saturation, truncation error, and transport delay.

The control surface models for elevators, ailerons, rudders, spoilers, and flaps include the rate limit, saturation, transport delay, and an optional dead zone with the following first-order dynamics and 5 Hz bandwidth (BW).

$$G_{cs}\left(s\right)={\displaystyle \frac{10\pi }{s+10\pi }}$$

(20)

Two Jetcat P70 engines were modeled, including throttle-to-RPM nonlinearity (as a cubic polynomial), RPM-to-thrust nonlinearity (as a cubic polynomial), a ram drag model, the effect of atmospheric pressure, and saturation with the first-order lag model shown below.

$$G_{eng}\left(s\right)={\displaystyle \frac{-0.1474s+0.7314}{s^2+1.336s+0.7314}}$$

(21)

The fuel flow rate is a parabolic function of the engine speed.The gyroscopic effect of the rotating rotor, off-axis thrust, torque generated by a thrust that is not aligned with the CM, and differential thrust is also accounted for.

2.1.2. Morphing UAV Model

Based on the nominal UAV model, a morphing UAV model was derived. The GTM includes six damage cases:

• Rudder Off;
• Vertical Tail Off;
• Left Outboard Flap Off;
• Left Wingtip Off;
• Left Elevator Off;
• Left Stabilizer Off.

Except for the fourth case, it corresponds to the loss of a control surface. Case 4 involves the loss of the outboard (approximately 25% of the semispan) left wingtip. The primary aerodynamic effect is a rolling moment bias that causes the left wing to drop (-$DeltadCl$), and this effect grows as alpha increases. Due to the location of the ailerons on the wingtip, the original model cannot use ailerons. In this study, however, it was assumed that the wingtip can telescope and the aileron can still function as it was installed slightly further inboard.

The effect of the roll axis moment can be maintained if the size of the aileron is increased proportionally to its proximity to the centerline. In this manner, the magnitude of the force generated by the ailerons increases; however, when the ailerons are used in opposite directions, the majority of the force is canceled out. It was assumed that the right wing could also telescope, similar to the left wing. Using the changes in mass, inertia matrix, and aerodynamic coefficient when the tip of the left wing is removed, it is possible to model the mass, inertia matrix, and aerodynamic coefficient when the left and right wings are adjusted independently. In this study, the left and right morphing parameters $etal$ and $etar$ were defined, representing the length variation ratio of each wing between −25% and 0%. The variables for symmetric and asymmetric morphing were then defined as follows:

$$\begin{array}{cc}\hfill {\eta }_s& ={\displaystyle \frac{{\eta }_l+{\eta }_r}{2}}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\eta }_a& ={\displaystyle \frac{{\eta }_l-{\eta }_r}{2}}\hfill \end{array}$$

(23)

It was also assumed that the differences in aerodynamic coefficients between the shortest and longest variable spans can be approximated by linear interpolation with acceptable error. Note that the CM and inertia matrices can be obtained precisely using the model data. It can be demonstrated that the shift in the center of gravity is a linear function of the morphing variables, whereas the MOIs and POIs are parabolic functions.

Variations in the mass ${\Delta }_d\mathrm{m}$, CM ${\Delta }_d{\mathbf{r}}_{cm}$, inertia matrix ${\Delta }_{d}J$, and aerodynamic coefficients ${\Delta }_d\mathbf{C}$ are known in the damaged configuration. Variations in the mass $\Delta \mathrm{m}$, CM $\Delta {\mathbf{r}}_{cm}$, inertia matrix $\Delta J$, and aerodynamic coefficients $\Delta \mathbf{C}$ can be represented as functions of the morphing parameters in the morphed configuration. The quasistatic assumption was adopted in this study, and the unsteady effect during the morphing process was not modeled because the original model does not include this effect.

The aircraft can be divided into three components: the moving left wingtip (subscript l), stationary main body (subscript b), and moving right wingtip (subscript r). The CM variations in the left and right wingtips in the morphed configuration can be represented as

$$\Delta {\mathbf{r}}_l=-{\displaystyle \frac{{\eta }_l}{25}}\left[\begin{array}{c}{\displaystyle \frac{btan\Lambda }{8}}\\ {\displaystyle \frac{b}{8}}\\ {\displaystyle \frac{btan\Gamma }{8cos\Lambda }}\end{array}\right],\Delta {\mathbf{r}}_r=-{\displaystyle \frac{{\eta }_r}{25}}\left[\begin{array}{c}{\displaystyle \frac{btan\Lambda }{8}}\\ -{\displaystyle \frac{b}{8}}\\ {\displaystyle \frac{btan\Gamma }{8cos\Lambda }}\end{array}\right]$$

(24)

where $\Lambda$ is the leading edge sweep angle and $\Gamma$ is the wing dihedral angle. The gross CM variation in the morphed configuration can be represented as

$$\Delta {\mathbf{r}}_{cm}=\left({\displaystyle \frac{{\mathrm{m}}_l}{\mathrm{m}}}\right)\Delta {\mathbf{r}}_l+\left({\displaystyle \frac{{\mathrm{m}}_r}{\mathrm{m}}}\right)\Delta {\mathbf{r}}_r$$

(25)

where ${\mathrm{m}}_l$ and ${\mathrm{m}}_r$ are the masses of the left and right wingtips, respectively. Finally, the morphed gross CM can be represented as

$${\mathbf{r}}_{cm}={\overline{\mathbf{r}}}_{cm}+\Delta {\mathbf{r}}_{cm}$$

(26)

Note that

$${\mathrm{m}}_l={\mathrm{m}}_r={\Delta }_d\mathrm{m}$$

(27)

and

$$\mathrm{m}={\mathrm{m}}_l+{\mathrm{m}}_b+{\mathrm{m}}_r$$

(28)

where ${\mathrm{m}}_b$ is the main body mass.

The nominal gross CM can be represented as

$$\mathrm{m}{\overline{\mathbf{r}}}_{cm}={\mathrm{m}}_l{\overline{\mathbf{r}}}_l+{\mathrm{m}}_b{\overline{\mathbf{r}}}_b+{\mathrm{m}}_r{\overline{\mathbf{r}}}_r$$

(29)

where ${\overline{\mathbf{r}}}_l$, ${\overline{\mathbf{r}}}_b$, and ${\overline{\mathbf{r}}}_l$ are the nominal CM of the left wingtip, main body, and right wingtip, respectively. Note that

$${\overline{\mathbf{r}}}_r=\left[\begin{array}{ccccc}1& & 0& & 0\\ 0& & -1& & 0\\ 0& & 0& & 1\end{array}\right]{\overline{\mathbf{r}}}_l$$

(30)

The damaged gross CM can be represented as

$$(\mathrm{m}-{\mathrm{m}}_l)({\overline{\mathbf{r}}}_{cm}+{\Delta }_d{\mathbf{r}}_{cm})={\mathrm{m}}_b{\overline{\mathbf{r}}}_b+{\mathrm{m}}_r{\overline{\mathbf{r}}}_r$$

(31)

Subtracting Equation (31) from Equation (29) yields ${\overline{\mathbf{r}}}_l$ as

$${\mathrm{m}}_l{\overline{\mathbf{r}}}_{cm}-(\mathrm{m}-{\mathrm{m}}_l){\Delta }_d{\mathbf{r}}_{cm}={\mathrm{m}}_l{\overline{\mathbf{r}}}_l$$

(32)

which can be rewritten as

$${\overline{\mathbf{r}}}_l={\overline{\mathbf{r}}}_{cm}-\left({\displaystyle \frac{\mathrm{m}}{{\mathrm{m}}_l}}-1\right){\Delta }_d{\mathbf{r}}_{cm}$$

(33)

Now, ${\overline{\mathbf{r}}}_r$ can be obtained by Equation (30). Then, ${\overline{\mathbf{r}}}_b$ can be calculated by either Equation (29) or Equation (31). The components’ morphed CM can be represented as

$${\mathbf{r}}_i={\overline{\mathbf{r}}}_i+\Delta {\mathbf{r}}_i,i=l,b,r$$

(34)

By applying the parallel axis theorem, also known as the Huygens–Steiner theorem, the nominal and morphed gross inertia matrix can be represented as

$$\overline{J}={\overline{J}}_l-{\mathrm{m}}_l{[{\overline{\mathbf{r}}}_l-{\overline{\mathbf{r}}}_{cm}]}^2+{\overline{J}}_b-{\mathrm{m}}_b{[{\overline{\mathbf{r}}}_b-{\overline{\mathbf{r}}}_{cm}]}^2+{\overline{J}}_r-{\mathrm{m}}_r{[{\overline{\mathbf{r}}}_r-{\overline{\mathbf{r}}}_{cm}]}^2$$

(35)

$$J={\overline{J}}_l-{\mathrm{m}}_l{[{\mathbf{r}}_l-{\mathbf{r}}_{cm}]}^2+{\overline{J}}_b-{\mathrm{m}}_b{[{\mathbf{r}}_b-{\mathbf{r}}_{cm}]}^2+{\overline{J}}_r-{\mathrm{m}}_r{[{\mathbf{r}}_r-{\mathbf{r}}_{cm}]}^2$$

(36)

where ${\overline{J}}_l$, ${\overline{J}}_b$, and ${\overline{J}}_r$ are the inertia matrices about their own centers of mass. Note that $\left[\mathbf{r}\right]$ is a skew-symmetric matrix associated with $\mathbf{r}={\left[xyz\right]}^T$ as

$$\left[\mathbf{r}\right]=\left[\begin{array}{ccccc}0& & -z& & y\\ z& & 0& & -x\\ -y& & x& & 0\end{array}\right]$$

(37)

The following equation can be useful:

$$-{\left[\mathbf{r}\right]}^2=\left[\begin{array}{ccccc}{y}^2+z^2& & -xy& & -xz\\ -yx& & x^2+z^2& & -yz\\ -zx& & -xy& & x^2+y^2\end{array}\right]=\mathrm{tr}\left(\mathbf{r}{\mathbf{r}}^T\right)I_3-\mathbf{r}{\mathbf{r}}^T$$

(38)

Finally, the inertia matrix variation due to morphing can be obtained by subtracting Equation (35) from Equation (36) as

$$\begin{array}{cc}\hfill \Delta J=& {\mathrm{m}}_l{[{\overline{\mathbf{r}}}_l-{\overline{\mathbf{r}}}_{cm}]}^2+{\mathrm{m}}_b{[{\overline{\mathbf{r}}}_b-{\overline{\mathbf{r}}}_{cm}]}^2+{\mathrm{m}}_r{[{\overline{\mathbf{r}}}_r-{\overline{\mathbf{r}}}_{cm}]}^2\hfill \\ \hfill -& {\mathrm{m}}_l{[{\mathbf{r}}_l-{\mathbf{r}}_{cm}]}^2-{\mathrm{m}}_b{[{\mathbf{r}}_b-{\mathbf{r}}_{cm}]}^2-{\mathrm{m}}_r{[{\mathbf{r}}_r-{\mathbf{r}}_{cm}]}^2\hfill \end{array}$$

(39)

The parameter variations due to ${\eta }_l$ and ${\eta }_s$ are shown in Figure 2 and Figure 3. The x-axis moment of inertia can be varied up to 23%, which means that the roll rate can be increased up to 30% by retracting the span due to angular momentum conservation. The CM can be shifted up to 1.3% MAC off the plane of symmetry by asymmetric morphing. Note that the CM moving back and forth up to 0.23% MAC due to a nonzero sweep may affect longitudinal stability, where the nominal static margin is 3.01% MAC.

In the variable-span morphing wing, the lift coefficient is nearly linearly proportional to the span [33] because the area of the lifting surface is linearly proportional to the span. The total drag in the subsonic flow consists of the profile drag and the induced drag. The drag coefficient for the induced drag of a high-aspect lightly swept wing in subsonic flow can be modeled as follows [34]:

$$C_{D_i}={\displaystyle \frac{C_L^2}{\pi eAR}}$$

(40)

where e is the efficiency factor, which is close to unity, and $AR$ is the aspect ratio defined as follows:

$$AR={\displaystyle \frac{b^2}{S}}$$

Therefore, the induced drag also varies linearly with respect to ${\eta }_s$. When an aircraft’s wing is not stalled, the aircraft’s parasite drag is nearly entirely composed of skin friction [32]. The quantity of skin friction drag depends on the aircraft’s wetted area. Therefore, the parasite drag and, in turn, the total drag can also be assumed to linearly vary with respect to ${\eta }_s$.

Similarly, the coefficients for the pitching moment and crosswind force can be assumed to vary linearly. In contrast to the longitudinal and force coefficients, the rolling and yawing moment coefficients vary in a parabola with a span due to the linear variation in the moment arm. However, if the moving portion of the span is relatively short compared to the entire wing, the effect of variation at the moment arm may be negligible. Therefore, the morphing-induced increase in the rolling and yawing moment coefficients can also be assumed to be proportional to the morphing parameter.

The aerodynamic coefficient variation due to left span morphing can be represented as

$${\Delta }_l\mathbf{C}(\alpha ,\beta )=-\left({\displaystyle \frac{{\eta }_l}{25}}\right){\Delta }_d\mathbf{C}(\alpha ,\beta )$$

(41)

where ${\Delta }_d\mathbf{C}$ is the aerodynamic coefficient variation due to left wingtip loss. Note that

$${\Delta }_r\mathbf{C}(\alpha ,\beta )=\mathrm{diag}(1,-1,1,-1,1,-1){\Delta }_l\mathbf{C}(\alpha ,-\beta )$$

(42)

Finally, the aerodynamic coefficient variation in the morphed configuration can be represented as

$$\Delta \mathbf{C}(\alpha ,\beta )={\Delta }_l\mathbf{C}(\alpha ,\beta )+{\Delta }_r\mathbf{C}(\alpha ,\beta )$$

(43)

When $\beta \approx 0$ in Equation (42),

$$\begin{array}{cc}\hfill {\Delta }_r{\mathbf{C}}_{lon}& ={\Delta }_l{\mathbf{C}}_{lon}\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\Delta }_r{\mathbf{C}}_{lat}& =-{\Delta }_l{\mathbf{C}}_{lat}\hfill \end{array}$$

(45)

where ${\mathbf{C}}_{lon}={\left[C_X{C}_Z{C}_m\right]}^T$ and ${\mathbf{C}}_{lat}={\left[C_Y{C}_l{C}_n\right]}^T$. Then,

$$\begin{array}{cc}\hfill \Delta {\mathbf{C}}_{lon}& =-2\left({\displaystyle \frac{{\eta }_s}{25}}\right){\Delta }_d\mathbf{C}\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill \Delta {\mathbf{C}}_{lat}& =-2\left({\displaystyle \frac{{\eta }_a}{25}}\right){\Delta }_d\mathbf{C}\hfill \end{array}$$

(47)

Note that the longitudinal motion is affected by symmetric morphing, whereas the lateral-directional motion is affected by asymmetric morphing. If the configuration where ${\eta }_l={\eta }_r=-12.5$% is selected as a nominal configuration, the maximum variation in the aerodynamic coefficients becomes $\pm {\Delta }_d\mathbf{C}$. The morphing parameters are suitable for the control input because they appear affinely in the forces and moments. The increments are nearly linearly proportional in the prestall region where $\alpha \in [-510]$ deg, while the conventional control surface increments are only lightly affected by the AOA, as shown in Figure 4. Consequently, the control effectiveness of asymmetric morphing in the roll axis is low when the AOA is small. However, the effectiveness of asymmetric morphing becomes even greater than the ailerons in the poststall region ($\alpha >10$ deg), as shown in Figure 4. Furthermore, symmetric morphing can be used to delay stall in low-speed flights or high-g maneuvers by increasing the lift coefficient, which, in turn, decreases the AOA toward the prestall region while generating the same lift.

The nonlinear model used in the following simulation scenarios is shown in Figure 5.

2.2. Derivation of an LPV Model of a Morphing UAV

2.2.1. Trim Analysis and Scheduling Parameter Selection

In this study, the scheduling parameters were airspeed and altitude. In addition, the parameter for symmetric morphing was included in the scheduling parameters because it has significant effects on longitudinal aerodynamics, which affects the overall flight performance. The asymmetric morphing parameter has a minimal effect on the longitudinal motion, but it can be utilized as a control input for the lateral motion.

Note that a rectangular parameter grid is necessary for the LPV controller design. However, the trimmable region of the GTM appears triangular in the TAS versus altitude diagram. The service ceiling is approximately 30,000 feet, where the trimmable TAS converges to 100 knots. Therefore, a new speed-related parameter must be introduced so that the parameter domain becomes rectangular. In this study, the following synthetic airspeed was introduced [35]:

$$V_S=V_0-{\displaystyle \frac{h_0}{h_0-h}}(V_0-V_T),h<h_0$$

(48)

where $V_0$ = 100 kt and $h_0$ = 30,000 ft. Note that Equation (48) becomes singular when $h=h_0$. However, the normal operating altitude of medium-size UAVs such as GTM is much lower than 30,000 ft. Then, the flight envelope becomes a rectangle, as shown in Figure 6.

Now, the scheduling parameter vector $\mathsf{\rho }\triangleq {\left[Vh{\eta }_s\right]}^T$ with three parameters is considered. The considered parameter grid is

• $V_S\in [60,70,80,90,100,110,120]$ in kt;
• $h\in [0,5000,\mathrm{10,000}]$ in ft;
• ${\eta }_s\in [-25,-12.5,0]$ in %.

As a result of the trim analysis, the true airspeed that can satisfy the trim condition at sea level is between approximately 60 kt and 120 kt. This trim airspeed range tends to narrow on both sides as the altitude increases, which appears constant even when the altitude changes if the synthetic airspeed is introduced. A 10 kt airspeed interval was set to properly reflect the nonlinear dependence of the aerodynamic force and torque. Atmospheric density changes as altitude changes, leading to nonlinear changes in dynamic pressure. The altitude was set at 5000 ft intervals because the nonlinearity is not strong in the atmospheric density curve. In the case of the morphing shape, the models at both ends were given from the nominal and damaged models, and the intermediate shape point was added to the parameter domain.

Remark 3.

The range of speed was set within the range where the GTM could maintain its trim state. The range of altitude is the normal operational altitudes. The range of the symmetric morphing parameter is the range provided by the damage model. If the number of the grid points is increased, the computation time may become excessive or a feasible solution may not be found. The interval of each variable was set to an appropriate value that captures enough nonlinearity but does not overly increase the total number of the grid points. Whether it is dense enough can be verified by comparing the simulation results for the LPV model and the nonlinear model or by investigating the results of the nonlinear simulation.

2.2.2. Lpv Modeling

The motion equations can be decoupled into equations describing the longitudinal motion and lateral-directional motion. The LTI models for the decoupled motions can be obtained at each parameter grid through Jacobian linearization. The wind axis was chosen in this study because of its advantages in decoupling.

The wind-axis longitudinal motion equations can be expressed as

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{V}}_T\\ \dot{\alpha }\\ \dot{q}\\ \dot{\theta }\end{array}\right]& =\left[\begin{array}{ccccccc}{X}_V+X_{T_V}cos({\alpha }_0+{\alpha }_T)& & X_{\alpha }& & 0& & -gcos{\gamma }_0\\ {\displaystyle \frac{Z_V-Z_{T_V}sin({\alpha }_0+{\alpha }_T)}{V_0}}& & {\displaystyle \frac{Z_{\alpha }}{V_0}}& & 1+{\displaystyle \frac{Z_q}{V_0}}& & -{\displaystyle \frac{gsin{\gamma }_0}{V_0}}\\ M_V+M_{T_V}& & M_{\alpha }+M_{T_{\alpha }}& & M_q& & 0\\ 0& & 0& & 1& & 0\end{array}\right]\left[\begin{array}{c}{V}_T\\ \alpha \\ q\\ \theta \end{array}\right]\hfill \\ & +\left[\begin{array}{ccc}{X}_{{\delta }_t}cos({\alpha }_0+{\alpha }_T)& & X_{{\delta }_e}\\ -{\displaystyle \frac{X_{{\delta }_t}sin({\alpha }_0+{\alpha }_T)}{V_0}}& & {\displaystyle \frac{Z_{{\delta }_e}}{V_0}}\\ M_{{\delta }_t}& & M_{{\delta }_e}\\ 0& & 0\end{array}\right]\left[\begin{array}{c}{\delta }_t\\ {\delta }_e\end{array}\right]\hfill \end{array}$$

(49)

The short-period mode dynamics can be extracted by eliminating the phugoid mode dynamics as

$$\left[\begin{array}{c}\dot{\alpha }\\ \dot{q}\end{array}\right]=\left[\begin{array}{ccc}{\displaystyle \frac{Z_{\alpha }}{V_0}}& & 1+{\displaystyle \frac{Z_q}{V_0}}\\ M_{\alpha }+M_{T_{\alpha }}& & M_q\end{array}\right]\left[\begin{array}{c}\alpha \\ q\end{array}\right]+\left[\begin{array}{c}{X}_{{\delta }_e}\\ {\displaystyle \frac{Z_{{\delta }_e}}{V_0}}\end{array}\right]{\delta }_e$$

(50)

The short-period mode dynamic was used for a longitudinal CAS design. The z-axis acceleration can be augmented to the measurement equation as

$$\left[\begin{array}{c}\alpha \\ q\\ a_z\end{array}\right]=\left[\begin{array}{ccc}1& & 0\\ 0& & 1\\ Z_{\alpha }& & 0\end{array}\right]\left[\begin{array}{c}\alpha \\ q\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ Z_{{\delta }_e}\end{array}\right]{\delta }_e$$

(51)

Note that the transfer function from ${\delta }_e$ to $a_z$ has a nonminimum phase (NMP) zero. In this case, the initial response to a step input may have the opposite sign to the final response. The altitude dynamics can be augmented to the longitudinal dynamics as

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{V}}_T\\ \dot{\alpha }\\ \dot{q}\\ \dot{\theta }\\ \dot{h}\end{array}\right]& =\left[\begin{array}{ccccccccc}{X}_V+X_{T_V}cos({\alpha }_0+{\alpha }_T)& & X_{\alpha }& & 0& & -gcos{\gamma }_0& & 0\\ {\displaystyle \frac{Z_V-Z_{T_V}sin({\alpha }_0+{\alpha }_T)}{V_0}}& & {\displaystyle \frac{Z_{\alpha }}{V_0}}& & 1+{\displaystyle \frac{Z_q}{V_0}}& & -{\displaystyle \frac{gsin{\gamma }_0}{V_0}}& & 0\\ M_V+M_{T_V}& & M_{\alpha }+M_{T_{\alpha }}& & M_q& & 0& & 0\\ 0& & 0& & 1& & 0& & 0\\ 0& & -V_0& & 0& & V_0& & 0\end{array}\right]\left[\begin{array}{c}{V}_T\\ \alpha \\ q\\ \theta \\ h\end{array}\right]\hfill \\ & +\left[\begin{array}{ccc}{X}_{{\delta }_t}cos({\alpha }_0+{\alpha }_T)& & X_{{\delta }_e}\\ -{\displaystyle \frac{X_{{\delta }_t}sin({\alpha }_0+{\alpha }_T)}{V_0}}& & {\displaystyle \frac{Z_{{\delta }_e}}{V_0}}\\ M_{{\delta }_t}& & M_{{\delta }_e}\\ 0& & 0\\ 0& & 0\end{array}\right]\left[\begin{array}{c}{\delta }_t\\ {\delta }_e\end{array}\right]\hfill \end{array}$$

(52)

The augmented dynamics were used for the longitudinal autopilot design.

For the control of the lateral-directional motion, the asymmetric morphing parameter was included in the control input. The lateral-directional dynamics, where the differential thrust is not considered, can be represented as

$$\left[\begin{array}{c}\dot{\beta }\\ \dot{p}\\ \dot{r}\\ \dot{\varphi }\end{array}\right]=\left[\begin{array}{ccccccc}{\displaystyle \frac{Y_{\beta }}{V_0}}& & {\displaystyle \frac{Y_p}{V_0}}& & {\displaystyle \frac{Y_r}{V_0}}-1& & {\displaystyle \frac{gcos{\theta }_0}{V_0}}\\ L_{\beta }& & L_p& & L_r& & 0\\ N_{\beta }& & N_p& & N_r& & 0\\ 0& & {\displaystyle \frac{cos{\gamma }_0}{cos{\theta }_0}}& & {\displaystyle \frac{sin{\gamma }_0}{cos{\theta }_0}}& & 0\end{array}\right]\left[\begin{array}{c}\beta \\ p\\ r\\ \varphi \end{array}\right]+\left[\begin{array}{ccccc}{\displaystyle \frac{Y_{{\delta }_a}}{V_0}}& & {\displaystyle \frac{Y_{{\delta }_r}}{V_0}}& & {\displaystyle \frac{Y_{{\eta }_a}}{V_0}}\\ L_{{\delta }_a}& & L_{{\delta }_r}& & L_{{\eta }_a}\\ N_{{\delta }_a}& & N_{{\delta }_r}& & N_{{\eta }_a}\\ 0& & 0& & 0\end{array}\right]\left[\begin{array}{c}{\delta }_a\\ {\delta }_r\\ {\eta }_a\end{array}\right]$$

(53)

Equation (53) was used for the lateral-directional autopilot design. When the airspeed is sufficiently high, the roll angle dynamics become negligible. In this case, the roll angle dynamics can be eliminated as

$$\left[\begin{array}{c}\dot{\beta }\\ \dot{p}\\ \dot{r}\end{array}\right]=\left[\begin{array}{ccccc}{\displaystyle \frac{Y_{\beta }}{V_0}}& & {\displaystyle \frac{Y_p}{V_0}}& & {\displaystyle \frac{Y_r}{V_0}}-1\\ L_{\beta }& & L_p& & L_r\\ N_{\beta }& & N_p& & N_r\end{array}\right]\left[\begin{array}{c}\beta \\ p\\ r\end{array}\right]+\left[\begin{array}{ccccc}{\displaystyle \frac{Y_{{\delta }_a}}{V_0}}& & {\displaystyle \frac{Y_{{\delta }_r}}{V_0}}& & {\displaystyle \frac{Y_{{\eta }_a}}{V_0}}\\ L_{{\delta }_a}& & L_{{\delta }_r}& & L_{{\eta }_a}\\ N_{{\delta }_a}& & N_{{\delta }_r}& & N_{{\eta }_a}\end{array}\right]\left[\begin{array}{c}{\delta }_a\\ {\delta }_r\\ {\eta }_a\end{array}\right]$$

(54)

Equation (54) was used for the lateral-directional CAS design. By applying interpolation to the family of LTI models, an LPV model can be obtained.

3. Cas Design Based on the LPV Method for Morphing-Assisted Maneuvers

Depending on the level of autonomy of the UAV, a human pilot is involved to some extent in controlling the UAV. In general, the speed of the rotation modes determines a UAV’s responsiveness to maneuvering commands. If the rotation modes are unstable or lightly damped, it is typically difficult for a human pilot to manually control the UAV. Therefore, a stability augmentation system (SAS) is required to ensure that these modes have the desired dynamic characteristics. In addition to stabilizing the mode, a control augmentation system is intended in order to provide a specific response to the command. Even though slow flight modes such as phugoid and spiral can be manually controlled, an automatic control system is necessary to relieve the pilot of hand flying because it is undesirable for a pilot to pay constant attention. An autopilot is an automatic control system that provides both pilot relief and specialized functions such as path following and automatic landing [32].

In this study, the control system design framework shown in Figure 7 was considered. The flight control system (FCS) depends on the control mode, guidance commands, and scheduling parameters. The CAS modes and autopilot modes for the longitudinal and lateral-directional channels were considered for the control modes. The guidance command was determined according to the intended form of flight. The controller gains of the FCS were scheduled on the morphing parameters and the flight conditions in a manner that ensures stability, and the actuator commands were computed where the morphing actuator is also included in the control inputs. The morphing system governs the morphing configuration in a manner that satisfies the control command, aids the commanded maneuver, and provides the desired dynamic characteristics based on flight conditions. Conventional control surfaces are usually designed to provide aerodynamic forces and moments primarily for the intended channel with minimal cross-coupling effects and impacts on the mass properties. However, the morphing configuration change affects the characteristics of the airframe dynamics differently than the conventional control surfaces.

In this section, CASs are designed for morphing-assisted maneuvers of morphing UAVs based on LPV methods to provide a rigorous stability guarantee under arbitrary morphing. A longitudinal CAS is designed to track the normal acceleration command in Section 3.1. A lateral-directional CAS is designed to track the AOS and roll rate commands in Section 3.2. The CASs utilize the symmetric morphing parameter for an improved performance and the asymmetric morphing parameter for flight control. The symmetric morphing strategy is discussed in Section 3.3. In Section 3.4, numerical simulation is performed for the push-over and pull-up and the high-g turn to demonstrate the effectiveness of the proposed scheme. The proposed gain self-scheduled flight control system is compared with the gain-scheduled ${\mathcal{H}}_{\infty }$ controller.

3.1. Longitudinal CAS Design for Normal Acceleration Control

The design of a normal acceleration CAS is based on LPV methods for controlling longitudinal motion. When the pilot must maneuver the UAV to its performance limits, high-performance UAVs require a specialized CAS. In this situation, normal acceleration is an appropriate variable for controlling the pitch axis. The accelerometer output has a component proportional to the AOA, allowing the short-period unstable mode to be stabilized. In addition, the accelerometer is typically less noisy and more reliable than the AOA sensor.

Note that the transfer function between the elevator deflection and normal acceleration has a nonminimum phase (NMP) zero. The initial response to a negative step elevator deflection will be a negative normal acceleration, followed by the expected positive normal acceleration. When the elevator is deflected upward to produce a positive normal acceleration, the tail experiences an increase in downward force. Consequently, the center of mass may momentarily fall, causing normal acceleration to decrease before increasing again. Therefore, the NMP behavior must be taken into account when designing the normal acceleration control system. The remaining degree of freedom in the throttle setting can be manually adjusted, or an autothrottle can be employed.

The closed-loop interconnection structure shown in Figure 8 was considered for the design of the longitudinal CAS. The 1-degree-of-freedom controller was considered. The open-loop plant P was composed of the state equation, Equation (50), and the output equation, Equation (51). The control design problem was cast in the model-matching framework. The control objective was to minimize the weighted normal acceleration model matching error and the weighted control effort against the command, noise, and disturbance.

The external input vector w was assumed to be a broadband ${\mathcal{L}}_2$ signal. For design convenience, the ${\mathcal{L}}_2$-norm of the external inputs can be regarded as unity. In this case, the control objective becomes making the performance output less than unity, and the weighting filters were designed accordingly. Note that the frequency-dependent weights can be designed as real, rational, and proper transfer matrices whose elements may depend on the scheduling parameter; that is, the performance objectives can also be scheduled on the scheduling parameter. The open-loop plant P was the LPV model obtained in Section 2.2. The actuator model $P_a$ was designed to reflect the realistic responses of the control surfaces to the control input commands. The reference weight $W_r$ corresponds to the maximum expected guidance command. The model weight $W_m$ corresponds to the ideal model response to the unit step guidance command. The performance weight $W_p$ shapes the relative importance of the model matching error throughout the frequency range. The noise weight $W_n$ corresponds to the expected sensor noise level across the frequency range. The control weight $W_c$ shapes the control effort penalty to the actuator command in the undesirable frequency range. The disturbance weight $D_a$ corresponds to the expected disturbance to the actuator output. The time delay $T_d$ corresponds to the expected transport delay.

The elevator was modeled as a first-order system with a bandwidth of 5 Hz. Except for some highly agile maneuvers, in most cases, the vertical acceleration command does not exceed 2 g, and, therefore, the reference weight was set to 2 g. The ideal model for the normal acceleration response to the command was set to have a natural frequency of 2.5 rad/s and a damping ratio of 0.8 to allow for some overshoot. Performance weights were set to ensure the model matching error does not exceed five times the normal acceleration measurement noise level at low frequencies. The sensor measurement error was modeled as white Gaussian noise, and the standard deviation was set to match the sensor specifications of the GTM. In the elevator model, a disturbance of 0.01 deg was added below 0.5 rad/s to account for errors due to the dead zone, saturation, and rate limit. A first-order Padé approximant corresponding to a transport delay of 30 ms was applied to all sensor measurements. The Laplace transform of a time delay of T is $e^{-sT}$, and the exponential transfer function can be approximated by a rational transfer function using Padé approximation formulas. Given the order of a rational function, the Padé approximant is known as the best approximation. High-order Padé approximations produce transfer functions with clustered poles. In general, Padé approximations with high order ($N>10$) are not preferred because their poles are sensitive to perturbations.

The generalized open-loop plant has 11 states and depends on three scheduling parameters. The LPV controller was synthesized to satisfy performance objectives according to the performance specifications defined in the generalized open-loop plant. The synthesized LPV controller guarantees that the closed-loop system is quadratically stable and that the ${\mathcal{L}}_2$ gain from the external input $\mathbf{w}\in {\mathcal{L}}_2$ to the performance output $\mathbf{z}$ is less than $\gamma \ge 0$ under arbitrary time-varying scheduling parameter $\mathsf{\rho }\in \mathcal{P}$. The LPV synthesis problem was solved twice. In the first iteration, an optimal solution that minimizes the induced ${\mathcal{L}}_2$-norm of $F_l(G,K)$ was obtained. In the second iteration, there is a suboptimal solution whose $\gamma$ is at most 20% greater than the $\gamma$ obtained in the first iteration for better numerical conditioning. As a result of LPV synthesis, suboptimal $\gamma =2.1505$ was obtained, which is 10.9% greater than the optimal $\gamma =1.9387$. The ${\mathcal{H}}_{\infty }$ controller was synthesized with respect to the pointwise LTI models, whose $\gamma$ ranges from 0.2006 to 0.2258. Note that the ${\mathcal{H}}_2$ controller can also be synthesized in this way. The synthesis time is largest in the LPV controller and smallest in the ${\mathcal{H}}_2$ controller.

3.2. Lateral-Directional CAS Design for Turn Coordination and Roll Rate Control

Regarding lateral-directional control, the most prevalent control augmentation system is the roll-rate command system. This mechanism can be constructed to roll the aircraft about its own velocity vector rather than the body axis. The roll-yaw stability augmentation system is adequate for the majority of aircraft, but a more refined lateral-directional CAS is required for aircraft that must maneuver rapidly at a high angle of attack (AOA). At high alpha, the lateral aerodynamic control surfaces have a tendency to cause the aircraft to roll about its longitudinal axis, which can lead to extremely undesirable phenomena such as the kinematic coupling of alpha and beta. The primary purpose of a roll is to initiate a turn, which is accomplished by utilizing the AOA to generate the lift that will ultimately produce the required centripetal acceleration.

The interconnection shown in Figure 9 was considered for the control design. The aileron, rudder, and morphing actuator were modeled as first-order systems, where the bandwidth was set to 5 Hz for the conventional control surfaces and 0.5 Hz for the morphing actuator to address a relatively slow response of the morphing actuators. Except for some highly agile maneuvers, in most cases, the AOS and the roll rate command do not exceed 3 deg and 10 deg/s, respectively. The ideal model for the AOS response to the command was set to have a natural frequency of 2 rad/s and a damping ratio of 0.8 to allow for some overshoot. The ideal model for the roll rate response to the command was set to have a natural frequency of 10 rad/s and a damping ratio of 0.8 to allow for faster convergence. Performance weights were set to ensure that the model matching error does not exceed two times the corresponding measurement noise level at low frequencies. The sensor measurement error was modeled as white Gaussian noise, and the standard deviation was set to match the sensor specifications of the GTM. In the actuator models, a disturbance of 0.01 deg was added below 0.5 rad/s to account for errors due to the dead zone, saturation, rate limit, and additional uncertainties arising from complex morphing mechanisms. A first-order Padé approximant corresponding to a transport delay of 30 ms was applied to all sensor measurements. The LPV controller and the ${\mathcal{H}}_{\infty }$ controller were synthesized.

3.3. Span Morphing Strategy

3.3.1. Effects of Span Morphing

In aircraft design, the wing span is one of the most crucial geometrical elements. It is directly related to the lift-induced drag of the wing: as the span increases, an increasing proportion of the inboard wing behaves as if in two-dimensional flow conditions [5]. Due to the increased distance between the vortices and the wing sections at the root, the downwash caused by the wingtip vortices is reduced in the inboard region of the wing. If the wing is swept, significant variations in the CM and CP can be anticipated, altering the trim and stability conditions. With span morphing, both the wing span and lifting area are altered. Therefore, the wing’s AR, which affects the induced drag, and the wetted area, which affects the profile drag, are both varied. For increasing $C_L$ values, configurations with a larger span may offer not only more lift but also less drag due to lower AOA requirements and less induced drag. Consequently, wing span primarily influences the total wing drag due to its effect on induced drag. To maximize the lift-to-drag ratio of the wing throughout the flight envelope, one can anticipate greater aerodynamic gains when operating at high $C_L$ values, i.e., at slower flight speeds.

Additionally, the effect on maneuverability must be taken into account. An increase in span has further effects on stability and control. The first and most obvious result is an increase in rolling MOI. Because the rolling MOI is proportional to the square of the distance between the mass elements and the longitudinal axis of the aircraft, there is a quadratic rise with the span. Therefore, if it is not possible to reduce the wing’s structural mass toward the tips, even if the device generating the rolling moment stays close to the wingtip, large reductions in the achievable roll rate can be anticipated. Increased the aerodynamic damping of the rolling motion due to a higher wingtip speed also contributes to a loss in maneuverability, which increases rolling stability. If the wing is swept, an increase in the span is likely to induce a change in the aircraft’s center of gravity, which will become more noticeable as the wing structural mass grows in relation to the aircraft’s mass. The position of the wing’s aerodynamic center is also modified. In turn, the trim drag and static stability margin are affected by changes in the center of gravity (CG) and aerodynamic center (AC). Similarly, if the span is raised asymmetrically, the CG and AC will move off the symmetry plane of the airplane, and there will be a stronger coupling between the longitudinal and lateral movements as the aircraft becomes more asymmetric [5].

3.3.2. Criteria for Span Variation

For improved agility and maneuverability, the symmetric morphing parameter should be increased when the normal acceleration command is large, and the symmetric morphing parameter should be decreased when the roll rate command is large.

3.4. Nonlinear Simulation of Morphing-Assisted Maneuvers

The GTM program was used to simulate a morphing aircraft’s flight dynamics. It implements general motion equations for rigid bodies for the vehicle dynamics and derives aerodynamic forces using a standard coefficient expansion performed as table lookups. The dynamics of the servo actuator and sensor bandwidth and errors were also included.

3.4.1. Push-Over and Pull-Up

First, a push-over and pull-up maneuver was considered. A normal acceleration command of 0g is engaged for the first 5 s; then, the 4g command is engaged. The throttle is fixed to 90% throughout the flight. Figure 10 and Figure 11 show the result using the LPV controller, and Figure 12 and Figure 13 show the result using the ${\mathcal{H}}_{\infty }$ controller. In case 1, the aircraft performs a push-over with the span shortened and then pulls up with the span extended, and, in case 2, the span is controlled in the opposite way. Both aircraft perform well before entering the stall, but the configuration profile in case 1 is superior in maintaining a longer pull-up. Additionally, the normal acceleration tracking performance of the ${\mathcal{H}}_{\infty }$ controller degrades when the morphing configuration changes rapidly, while the LPV controller shows a faster convergence to the command. However, the LPV-based CAS results in a larger roll angle because the CAS tries to control only the roll rate instead of the roll angle, which can be improved by implementing a roll angle autopilot. The simulation results for the ${\mathcal{H}}_2$ controller are omitted because the results are very similar to those for the ${\mathcal{H}}_{\infty }$ controller.

3.4.2. High-G Turn

Next, a high-g turn was considered. A roll rate command of 10 deg/s is engaged for the first 8 s to reach an 80 deg roll angle; then, a 0 deg/s command is engaged. A normal acceleration command of 1g is engaged for the first 8 s; then, the 5 g command is engaged. The throttle is fixed to the maximum throughout the flight. In case 1, the aircraft rotates on the roll axis with the span shortened, then increases the span and increases the normal acceleration to make a high-g turn. In case 2, the span is controlled in the opposite way. Figure 14, Figure 15, Figure 16 and Figure 17 shows that both aircraft perform well before entering the stall, but one configuration profile in case 1 is superior in maintaining a longer turn. Note that the LPV controller is quicker in convergence, but the ${\mathcal{H}}_{\infty }$ controller is better at maintaining the command longer where the configuration is not changed.

4. Autopilot Design Based on LPV Methods for Morphing-Assisted Flights

When the flight management system, including autopilots, is activated, a UAV can complete an entire mission without manual remote control. The majority of flying quality requirements have no direct bearing on the autopilot design. Pilot relief autopilot modes require the autopilot to meet standards for steady-state error and disturbance rejection while placing less importance on the dynamic response. To prevent distracting or potentially hazardous transient motions, careful consideration must also be given to how the autopilot is turned on and off.

In this section, an autopilot is designed for the morphing-assisted flight of morphing UAVs based on LPV methods. A longitudinal autopilot was designed to track the true airspeed and altitude commands. A lateral-directional autopilot was designed to track the AOS and roll angle commands. The autopilot utilizes symmetric and asymmetric morphing for improved agility and maneuverability, respectively. Numerical simulation is performed to demonstrate the effectiveness of the proposed scheme.

4.1. Longitudinal Autopilot Design for Airspeed and Altitude Control

An airspeed and altitude autopilot was designed to control longitudinal motion. Altitude hold is a vital pilot relief option that allows an aircraft to be maintained at a predetermined altitude in accordance with air traffic control requirements in an air route corridor. Typically, the speed-hold autopilot is utilized during climb and descent. During a climb, the throttles may be set to a relatively high level of power, and the elevator will receive a speed input to maintain a constant speed. While the speed will vary based on altitude, keeping the speed constant will maximize fuel efficiency. In addition, the descent will be conducted at a constant speed, with the throttles set to near idle. For effective cruising, both the throttle and the elevator will be used to maintain speed and altitude at the cruising altitude.

In longitudinal aircraft dynamics, there is significant coupling between the two control inputs (engine throttle and elevator deflection angle) and the two main outputs (speed and altitude). The interconnection shown in Figure 18 was considered for the control design.

The engine lag was modeled as a second-order system with NMP zero. The reference weight was set to 4 kt for the true airspeed and 10 ft for the altitude. The ideal models for the true airspeed and altitude response to the command were set to have a natural frequency of 2 rad/s and a damping ratio of 0.8 for a gentle response. Performance weights were set to ensure that the model matching error does not exceed two times the measurement noise levels at low frequencies. The sensor measurement error was modeled as white Gaussian noise, and the standard deviation was set to match the sensor specifications of the GTM. In the actuator model, a disturbance of 0.01 deg was added below 0.5 rad/s to account for various errors.

4.2. Lateral-Directional Autopilot Design for Turn Coordination and Roll Angle Control

A roll angle autopilot was designed based on LPV methods to control lateral-directional motion. The purpose of the autopilot is to keep the wings level. To achieve a coordinated turning motion, additional control systems must be used to regulate the sideslip and pitch rate if the aircraft is maintained at an angle other than at the wings level. In turn, the pitch rate command will determine whether the aircraft gains or loses altitude. If the roll reference can be altered, the aircraft can be directed in any direction with single control. These control systems can provide inner loops for various autopilots, allowing an aircraft to maintain a constant compass heading or follow a radio navigation beam despite crosswinds.

There is significant coupling between the control inputs and the main outputs in the lateral-direction aircraft dynamics. The interconnection shown in Figure 19 was considered for the control design.

The aileron, rudder, and morphing actuator were modeled in the same way as in the CAS design. The reference weight was set to 30 deg for the roll angle command, and this condition can be forced in the guidance command. The ideal model for the AOS response to the command was set to have a natural frequency of 2 rad/s and a damping ratio of 0.8 to allow for some overshoot. The ideal model for the roll rate response to the command was set to have a natural frequency of 5 rad/s and a damping ratio of 0.8 to allow for faster convergence. Performance weights were set to ensure that the model matching error does not exceed two times the corresponding measurement noise level at low frequencies. The sensor measurement error was modeled as white Gaussian noise, and the standard deviation was set to match the sensor specifications of the GTM. In the actuator model, a disturbance of 0.01 deg was added below 0.5 rad/s to account for errors due to the dead zone, saturation, rate limit, and additional uncertainties arising from complex morphing mechanisms.

4.3. Nonlinear Guidance for Trajectory Tracking

A lookahead distance was used to compute the desired course angle. The lateral acceleration command was generated from the guidance law [36].

$$a_{s_{cmd}}=2{\displaystyle \frac{V^2}{d}}sin\sigma$$

(55)

To convert the lateral acceleration command into the bank angle command, the relation $a_s\approx g\varphi$ derived from the aircraft coordinate turn was utilized. A heading controller was used to compute the required roll angle.

$${\varphi }_c={tan}^{-1}\left(K_{\psi }{\displaystyle \frac{V}{g}}({\psi }_c-\psi )\right)$$

(56)

$$\dot{\psi }+K_{\psi }\psi =K_{\psi }{\psi }_c$$

(57)

$${\displaystyle \frac{\psi \left(s\right)}{{\psi }_c\left(s\right)}}={\displaystyle \frac{K_{\psi }}{s+K_{\psi }}}$$

(58)

The heading controller gain was set to $K_{\psi }=3.9$.

The optimal span morphing parameter can be determined as a function of flight conditions (airspeed and altitude). Furthermore, the span morphing parameter can be increased to attenuate the altitude drop when a large roll angle command is engaged. In this study, the span morphing parameter deviation from the setpoint was obtained in proportion to the roll angle command.

4.4. Nonlinear Simulation of Morphing-Assisted Flights

4.4.1. Waypoint Following at Low Altitude

In the waypoint following at low altitude, the morphing UAV was commanded to follow four waypoints placed 3000 ft apart at the same altitude. The resulting trajectory and state response are shown in Figure 20, Figure 21 and Figure 22. Notably, the ${\mathcal{H}}_{\infty }$ controller exhibits larger oscillation in both the longitudinal and lateral-directional motions, possibly due to rapid parameter variations.

4.4.2. Circular Trajectory Tracking at High Altitude

In the circular trajectory tracking at high altitude, the morphing UAV was commanded to follow a horizontal orbit with a radius of 2000 ft. The resulting trajectory and state response are shown in Figure 23, Figure 24 and Figure 25. The ${\mathcal{H}}_{\infty }$ controller exhibits undesirable chattering phenomena in both the longitudinal and lateral-directional motions as a result of aggressive high gains.

4.4.3. Helical Ascent under Fast Morphing

In the helical ascent under fast morphing, the morphing UAV was commanded to follow a helical path with a radius of 2000 ft and a climb rate of 30 ft/s, which corresponds to a flight path angle of approximately 9.18 deg. The resulting trajectory and state response are shown in Figure 26, Figure 27 and Figure 28. Extremely fast morphing was used to test the marginal performance of the designed controllers. Note that the oscillation of the ${\mathcal{H}}_{\infty }$ controller gradually grows and eventually diverges, whereas the LPV controller maintains the tracking performance.

4.4.4. Spiral Descent with Morphing Scheduling

In spiral descent with morphing scheduling, the morphing UAV was commanded to follow a helical path with a radius of 2000 ft and a rate of descent of −20 ft/s, which corresponds to an approximate flight path angle of −8.98 deg. The symmetric morphing configuration was scheduled on the flight conditions so that the aerodynamic performance is improved. The resulting trajectory and state response are shown in Figure 29, Figure 30 and Figure 31. Note that the oscillation of the ${\mathcal{H}}_{\infty }$ controller quickly loses stability even when the scheduling parameter undergoes gentle changes, whereas the LPV controller maintains the tracking performance.

5. Conclusions

A novel framework is proposed for the flight control of morphing unmanned aerial vehicles (UAVs). The proposed scheme reflects the benefits of symmetric span morphing in the longitudinal performance while utilizing asymmetric span morphing for roll control. The control system was designed based on linear parameter-varying (LPV) methods that naturally suit the parameter-varying nature of morphing UAVs. First, a high-fidelity model of an asymmetric variable-span morphing UAV was derived from the NASA generic transport model. Second, the control augmentation system was designed based on LPV methods to track normal acceleration and roll rate commands while maintaining a small angle of sideslip. Finally, the autopilot was designed based on LPV methods to track airspeed, altitude, and roll angle commands while maintaining a small angle of sideslip.

Remark 4.

The presented LPV controller was obtained by solving complex optimization problems in the design process. However, once the design was carried out, the amount of online calculation was very small because the resulting controller has the form of a gain-scheduled linear controller. Therefore, the amount of online calculation is not a problem in practice.

Future research could focus on a technique for determining the optimal shape in real time. The optimal shape can also be determined offline by considering flight conditions only. Therefore, the online optimum shape determiner must be able to appropriately consider both current flight conditions and external guidance commands. When a large amount of lift is needed, for instance, a longer span can be extended even under the same flight conditions. When a high-speed dash or rapid rotation of a roll axis is needed, however, a short span can be advantageous. When only one morphing parameter is considered as a scheduling parameter, the problem can be easily solved using a straightforward line search. However, with two or more morphing parameters, it may be more challenging to solve the problem.