Nonlinear Incremental Control for Flexible Aircraft Trajectory Tracking and Load Alleviation
NNonlinear Incremental Control for Flexible Aircraft TrajectoryTracking and Load Alleviation
Xuerui Wang ∗ , Tigran Mkhoyan † and Roeland de Breuker ‡ Delft University of Technology, Faculty of Aerospace Engineering,Kluyverweg 1, 2629 HS Delft, The Netherlands
This paper proposes a nonlinear control architecture for flexible aircraft simultaneoustrajectory tracking and load alleviation. By exploiting the control redundancy, the gust andmaneuver loads are alleviated without degrading the rigid-body command tracking perfor-mance. The proposed control architecture contains four cascaded control loops: positioncontrol, flight path control, attitude control, and optimal multi-objective wing control. Sincethe position kinematics are not influenced by model uncertainties, the nonlinear dynamic in-version control is applied. On the contrary, the flight path dynamics are perturbed by bothmodel uncertainties and atmospheric disturbances; thus the incremental sliding mode controlis adopted. Lyapunov-based analyses show that this method can simultaneously reduce themodel dependency and the minimum possible gains of conventional sliding mode control meth-ods. Moreover, the attitude dynamics are in the strict-feedback form; thus the incrementalbackstepping sliding mode control is applied. Furthermore, a novel load reference generator isdesigned to distinguish the necessary loads for performing maneuvers from the excessive loads.The load references are realized by the inner-loop optimal wing controller, while the excessiveloads are naturalized by flaps without influencing the outer-loop tracking performance. Themerits of the proposed control architecture are verified by trajectory tracking tasks and gustload alleviation tasks in spatial von Kármán turbulence fields.
I. Introduction T he pursuit of higher efficiency has driven aircraft design evolving towards having higher aspect ratios and morelightweight structures. The resulting increase in structural flexibility makes the conventional frequency separationbetween rigid-body and structural degrees of freedom unreliable. During aggressive flight maneuvers or in the presenceof gusts, structural integrity is also challenged. Given this background, there has been increasing interest in flexibleaircraft control algorithms that can achieve multiple objectives, which include traditional pilot command-followingflight control, gust load alleviation (GLA), maneuver load alleviation (MLA), aeroelastic mode suppression, dragminimization, etc.In the literature, it is a common practice to design two independent controllers, one for rigid-body commandtracking and another for fulfilling the remaining objectives (GLA, MLA, flutter suppression, etc.) [1–3]. As reportedin [1], conflict can occur if the roll tracking command and the load alleviation command are given to the same wingtrailing-edge control surface. In order to avoid this conflict, the aerodynamic control surfaces are usually classifiedinto two sets. For example, in [1], two inboard elevators are used for pitch control, two outboard elevators are used forroll control, while all the ailerons are used for GLA. A model prediction controller is designed for MLA in [3]; theoutboard ailerons are used to achieve MLA, while the inboard ailerons together with elevator and rudder are selected forrigid-body control. A similar selection is also adopted in [2] for executing the commands provided by linear-quadraticregulators. However, the selection and isolation cannot make full use of the control input space, and can cause alarge gap between two adjacent flaps. An integrated controller that considers both rigid-body and aeroelastic controlobjectives, while utilizing the control redundancy is desirable.Another challenge faced by multi-objective flexible aircraft flight control algorithms is the balance between loadalleviation performance and rigid-body command tracking performance. For example, the simulation results in [4] ∗ Assistant Professor, Department of Aerospace Structures and Materials, Department of Control and Operations, Faculty of AerospaceEngineering, Kluyverweg 1, 2629HS Delft, the Netherlands, [email protected], AIAA Member. † Ph.D. Candidate, Department of Aerospace Structures and Materials, Faculty of Aerospace Engineering, Kluyverweg 1, 2629HS Delft, theNetherlands, [email protected], AIAA Student Member. ‡ Associate Professor, Department of Aerospace Structures and Materials, Faculty of Aerospace Engineering, Kluyverweg 1, 2629HS Delft, theNetherlands, [email protected], AIAA Senior Member. a r X i v : . [ ee ss . S Y ] J a n how that the pitch rate command tracking is poor because the MLA function is designed to reduce the wing rootbending moment [4]. Essentially, if a flexible aircraft has more than three independent control surfaces, and if not all ofthem are in the aircraft symmetrical plane, then there are infinite deflection combinations to realize the conventionalthree-axis attitude control moments. In other words, such an aircraft is over-actuated in rigid-body command trackingtasks. Therefore, it is physically realistic to simultaneously achieve desired rigid-body control moments together withother objectives (such as load alleviation), rather than making trade-offs among them.Recently, a control method that can alleviate loads without changing rigid-body tracking performance is presentedin [5, 6]. By exploiting the null space between the input and rigid-body output, the load alleviation objective is decoupledfrom the rigid-body tracking objective. However, this method has two main limitations: first, the design of the nullspace filters strongly depends on the matrix fraction description of the linear state-space model; thus the decouplingeffectiveness is questioned in the presence of model uncertainties and unmodeled nonlinearities; second, the proposedcontroller only ensures the load bounds are not violated but does not further minimize the load variations within thebounds. These two issues will be tackled in this paper.Incremental control refers to a class of nonlinear control methods that replace a part of model information by sensormeasurements. The most well-known incremental control method is the incremental nonlinear dynamic inversion (INDI)control, which was initially proposed in [7]. In the past decade, the reduced model-dependency and strong robustness ofINDI have stimulated its application to various aerospace systems (a CS-25 certified aircraft [8], quadrotors [9], launchvehicles [10], etc.). Recently, Ref. [11] presented more generalized derivations, as well as Lyapunov-based stabilityand robustness analyses for the INDI control. Theoretical analyses also show INDI is more robust than the nonlineardynamic inversion used in [12]. The INDI control law has been used to solve a free-flying flexible aircraft GLA problemin [13]. The control objectives of rigid-body motion regulation, GLA, and structural vibration suppression are achievedsimultaneously. Moreover, theoretical analyses and various time-domain simulations demonstrate the robustness ofINDI control to atmospheric disturbance, aerodynamic model uncertainties, and even actuator faults. Apart from thesemerits, the controller designed in [13] also has limitations. Because the number of controlled states is higher than thenumber of control inputs, the weight-least square method is used in the INDI loop for making trade-offs among differentvirtual control channels. As a consequence, not all the designed virtual controls can be fully executed. Although thiscontrol structure works well for the GLA problem, it would degrade aircraft tracking performance if is directly extendedto MLA problems. In view of this issue, a different control architecture will be proposed in this paper.To further expand the applicability and enhance the robustness of INDI control, the incremental sliding mode controlframework (INDI-SMC) was proposed in [14]. This hybrid framework is derived for generic multi-input/multi-outputnonlinear systems and is able to simultaneously reduce the minimum possible control/observer gains as well as themodel dependency of conventional SMC methods. Furthermore, theoretical proofs and simulations demonstrate that awide range of model uncertainties, sudden actuator faults, and structural damage can be passively resisted by INDI-SMC,without using additional fault detection and isolation modules [14]. The hybridization idea was carried forward in [15],where the incremental backstepping sliding mode control (IBSMC) is proposed for multi-input/multi-output nonlinearstrict-feedback perturbed systems. Theoretical analyses and simulations show IBSMC has less model dependency butenhanced robustness than model-based backstepping sliding mode control in the literature. These advantages makeINDI-SMC and IBSMC promising for solving flexible aircraft control problems.The main contribution of this paper is a nonlinear control architecture designed for flexible aircraft trajectorytracking and load alleviation purposes. The proposed control architecture has the following features: 1) it has enhancedrobustness against model uncertainties, external disturbances and faults despite its reduced model dependency; 2) itcontains two load reference generation algorithms to distinguish the loads that are necessary to perform maneuvers fromthe excessive loads; 3) it can neutralize the excessive loads without degrading the rigid-body tracking performance, nomatter the excessive loads are induced by maneuvers or atmospheric disturbances.Apart from the control design, the kinematics and dynamics of a free-flying flexible aircraft are also derived inthis paper. A body-fixed reference frame is used to capture both the inertial and aerodynamic couplings between therigid-body and structural degrees of freedom. Furthermore, a modular approach to conveniently make an existingclamped-wing aeroservoelastic model free-flying is presented in this paper. This modular approach contributes tobridging the flight dynamic and aeroelasticity communities.The rest of this paper is structured as follows. The flexible aircraft kinematics and dynamics are derived in Sec. II.The nonlinear control architecture is presented in Sec. III. The simulation results are presented in Sec. IV with theconclusions drawn in Sec. V. In this paper, bold symbols represent vectors and matrices.2 I. Model Descriptions
A. Flight Kinematics
An overview of the flexible aircraft model used in this paper is shown in Fig. 1. To begin with, the following
Fig. 1 Reference frames and axis definitions of the flexible aircraft. reference frames are defined: • Inertia Frame: F 𝐼 ( 𝑂 𝐼 , 𝑥 𝐼 , 𝑦 𝐼 , 𝑧 𝐼 ). Under the assumption that the Earth is flat and non-rotating, the Earth-centerEarth-fixed reference frame can be used as an inertia frame. • Body-fixed Reference Frame: F 𝐵 ( 𝑂 𝑏 , 𝑥 𝑏 , 𝑦 𝑏 , 𝑧 𝑏 ). In this paper, 𝑂 𝑏 is selected at a fixed point on the fuselage,which is not necessarily the center of mass (c.m.). This more general approach is adopted because during flightoperations, c.m. is time-varying due to structural deformations and fuel assumptions. 𝑂 𝑏 𝑥 𝑏 and 𝑂 𝑏 𝑧 𝑏 are definedin the undeformed aircraft symmetrical plane. • Wing-fixed Reference Frame: F 𝑊 ( 𝑂 𝑤 , 𝑥 𝑤 , 𝑦 𝑤 , 𝑧 𝑤 ). The left and right wings are modeled as two Euler Bernoullibeams. 𝑂 𝑤 is defined at the root of each beam, and the distance vector from 𝑂 𝑏 to 𝑂 𝑤 is denoted as 𝒓 𝑤𝑏 . 𝑂 𝑤 𝑥 𝑤 is aligned with the undeformed beam, and 𝑂 𝑤 𝑦 𝑤 is aligned with the mean-chord. • Flight Trajectory Axes: F 𝑉 ( 𝑂 𝑏 , 𝑥 𝑉 , 𝑦 𝑉 , 𝑧 𝑉 ). 𝑂 𝑏 𝑥 𝑉 is aligned with the aircraft ground velocity vector, and 𝑂 𝑏 𝑧 𝑉 pointing downwards in the plumb plane. • Aerodynamic Axes: F 𝐴 ( 𝑂 𝑏 , 𝑥 𝐴 , 𝑦 𝐴 , 𝑧 𝐴 ). 𝑂 𝑏 𝑥 𝐴 is aligned with the aircraft airspeed vector, and 𝑂 𝑏 𝑧 𝐴 pointingdownwards in the aircraft symmetrical plane.All the reference frames are right-handed. In Fig. 1, 𝐹 𝑤 represents the wing root shear force. 𝑀 𝜃 and 𝑀 𝜙 respectively represents the wing root torsion and bending moment.Define the distance vector from 𝑂 𝐼 to 𝑂 𝑏 expressed in the inertial frame as 𝑹 𝑏 = [ 𝑋, 𝑌 , 𝑍 ] T , then the translationalkinematics are (cid:164) 𝑋 (cid:164) 𝑌 (cid:164) 𝑍 = 𝑉 cos 𝜒 cos 𝛾𝑉 sin 𝜒 cos 𝛾 − 𝑉 sin 𝛾 (1)where 𝑉 is the ground velocity of the aircraft, 𝛾 is the flight path angle, and 𝜒 is the kinematic azimuth angle.The attitude of an aircraft can be described by three Euler angles. Apart from them, for evaluating the influences ofattitudes on aerodynamic forces, the following three attitude angles are more widely used: angle of attack 𝛼 , side-slipangle 𝛽 , and the kinematic bank angle 𝜇 . The attitude kinematics of an aircraft are given as follows [16]: (cid:164) 𝜇 (cid:164) 𝛼 (cid:164) 𝛽 = cos 𝛼 cos 𝛽 𝛼 sin 𝛽 𝛼 cos 𝛽 − cos 𝛼 − (cid:169)(cid:173)(cid:173)(cid:171) − 𝑪 𝐵𝑉 − (cid:164) 𝜒 sin 𝛾 (cid:164) 𝛾 (cid:164) 𝜒 cos 𝛾 + 𝑝𝑞𝑟 (cid:170)(cid:174)(cid:174)(cid:172) (2) 𝑪 𝑖 𝑗 is used to denote the direction cosine matrix from the reference frame F 𝑗 to the reference frame F 𝑖 . In thispaper, only the left and right wings are modeled as flexible beams, while the remaining aircraft components are assumed3o be rigid. B. Aeroservoelastic Model of the Wing
In this subsection, the aeroservoelastic model of the wing will be derived in the local F 𝑊 frame. The right-wingdynamics will be derived as an example, while the dynamics of the left-wing can be derived analogously. The rigid-bodymotions influence the wing dynamics in two aspects: 1) change the local angle of attack of the wing; 2) cause additionalinertial forces if rigid-body linear and angular accelerations are non-zero. Both the aerodynamic and inertial influencesof rigid-body motions will be considered in the derivations.
1. Structural Model
The wing structure is modeled as a linear dynamic Euler Bernoulli beam, which is discretized into 𝑛 𝑠 elements with 𝑛 𝑠 + 𝑤 𝑠 (downwards positive),torsion 𝜃 𝑠 (nose-up positive) and out-of plane bending 𝜙 𝑠 (bend-up positive). A distributed trailing-edge controlsurface configuration is adopted in this paper. At each node, a flap with a deflection angle 𝛽 𝑠 is connected to the beamthrough a rotational spring. 𝛽 𝑠 > 𝑿 𝑠 = [ 𝒙 T 𝑠, , ..., 𝒙 T 𝑠,𝑛 𝑠 ] T ∈ R ( 𝑛 𝑠 + ) ,where 𝒙 𝑠,𝑖 = [ 𝑤 𝑠,𝑖 , 𝜙 𝑠,𝑖 , 𝜃 𝑠,𝑖 , 𝛽 𝑠,𝑖 ] T , 𝑖 = , ..., 𝑛 𝑠 . Assume the beam is clamped at the first node, and denote 𝒙 𝑠 = [ 𝒙 T 𝑠, , ..., 𝒙 T 𝑠,𝑛 𝑠 ] T ∈ R 𝑛 𝑠 , then the corresponding dynamic equations are 𝑴 𝑠 (cid:165) 𝑿 𝑠 + 𝑲 𝑠 𝑿 𝑠 = (cid:34) 𝑴 𝑠, 𝑴 𝑠, 𝑴 T 𝑠, 𝑴 𝑠,𝑛 𝑠 (cid:35) (cid:34) (cid:165) 𝒙 𝑠, (cid:165) 𝒙 𝑠 (cid:35) + (cid:34) 𝑲 𝑠, 𝑲 𝑠, 𝑲 T 𝑠, 𝑲 𝑠,𝑛 𝑠 (cid:35) (cid:34) 𝒙 𝑠, 𝒙 𝑠 (cid:35) = (cid:34) 𝑭 root 𝑭 ext (cid:35) (3)In Eq. (3), 𝑴 𝑠 and 𝑲 𝑠 respectively represents the structure mass and stiffness matrices. The structural damping isneglected in this model. 𝑭 root is the reaction force at the clamped wing root. 𝑭 ext is the distributed external force vector,which contains gravitational and aerodynamic forces. Moreover, when the aircraft has linear and angular accelerations,the wing frame F 𝑊 becomes a non-inertial frame. Therefore, when modeling the structural dynamics in the wing frame,the inertial forces induced by rigid-body motions should be added to 𝑭 ext . This inertial coupling effect is not modeledin the conventional mean-axes method [17].
2. Inertial Forces
Consider an infinitesimal mass element dm on an arbitrary section 𝐴 𝑤 of the right wing. When the wing isundeformed, the distance vector from 𝑂 𝑤 to dm, expressed in the right wing reference frame is 𝒓 𝑤 = [ 𝑟 𝑥 , 𝑟 𝑦 , 𝑟 𝑧 ] T .Denote the transverse displacement of this wing section as 𝑤 , and its torsion angle around the 𝑂 𝑤 𝑥 𝑤 axis as 𝜃 , then theposition vector caused by deformation is 𝒓 𝑒 = 𝑟 𝑒,𝑥 𝑟 𝑒,𝑦 𝑟 𝑒,𝑧 = 𝑤 + 𝑟 𝑦 cos 𝜃𝑟 𝑦 sin 𝜃 (4)The absolute distance from the inertial frame origin 𝑂 𝐼 to 𝑂 𝑤 equals the summation of 𝑹 𝑏 (defined in the inertialframe) and 𝒓 𝑤𝑏 (defined in the body reference frame). Denote 𝑹 𝑤 as the absolute distance from 𝑂 𝐼 to dm, projected onthe right wing reference frame, then 𝑹 𝑤 = 𝑪 𝑊 𝐵 𝑪 𝐵𝐼 𝑹 𝑏 + 𝑪 𝑊 𝐵 𝒓 𝑤𝑏 + 𝒓 𝑤 + 𝒓 𝑒 (5)By differentiating the above equation, the absolute velocity of dm expressed in the right wing reference frame equals 𝑽 𝑤 = 𝑪 𝑊 𝐵 𝑽 𝑏 + 𝑪 𝑊 𝐵 { 𝝎 𝑏 × [ 𝒓 𝑤𝑏 + 𝑪 T 𝑊 𝐵 ( 𝒓 𝑤 + 𝒓 𝑒 )]} + 𝒗 𝑒 = 𝑪 𝑊 𝐵 𝑽 𝑏 + 𝑪 𝑊 𝐵 ˜ 𝝎 𝑏 𝒓 𝑤𝑏 + 𝑪 𝑊 𝐵 ˜ 𝝎 𝑏 𝑪 T 𝑊 𝐵 ( 𝒓 𝑤 + 𝒓 𝑒 ) + 𝒗 𝑒 (6)in which 𝑽 𝑏 = [ 𝑉 𝑥 , 𝑉 𝑦 , 𝑉 𝑧 ] T and 𝝎 𝑏 = [ 𝑝, 𝑞, 𝑟 ] T are the translational and rotational velocities of the body-frame, andboth of them are expressed in F 𝐵 . 𝒗 𝑒 is the relative deformation velocity of dm, which equals [ (cid:164) 𝑟 𝑒,𝑥 , (cid:164) 𝑟 𝑒,𝑦 , (cid:164) 𝑟 𝑒,𝑧 ] T . ( ˜ ·) denotes the skew-symmetric matrix of the vector (·) . 4urthermore, differentiate Eq. (6), the absolute acceleration of dm, expressed in the right wing reference frameequals 𝒂 𝑤 = 𝑪 𝑊 𝐵 𝒂 𝑓 + 𝑪 𝑊 𝐵 ˜ 𝜶 𝑓 𝒓 𝑤𝑏 + 𝑪 𝑊 𝐵 ˜ 𝜶 𝑓 𝑪 T 𝑊 𝐵 ( 𝒓 𝑤 + 𝒓 𝑒 ) + 𝑪 𝑊 𝐵 ˜ 𝝎 𝑏 𝑪 T 𝑊 𝐵 𝒗 𝑒 + 𝒂 𝑒 + 𝑪 𝑊 𝐵 ˜ 𝝎 𝑏 ˜ 𝝎 𝑏 𝒓 𝑤𝑏 + 𝑪 𝑊 𝐵 ˜ 𝝎 𝑏 ˜ 𝝎 𝑏 𝑪 T 𝑊 𝐵 ( 𝒓 𝑤 + 𝒓 𝑒 ) (7)where 𝒂 𝑏 = [ 𝑎 𝑥 , 𝑎 𝑦 , 𝑎 𝑧 ] T and 𝜶 𝑏 = [ (cid:164) 𝑝, (cid:164) 𝑞, (cid:164) 𝑟 ] T are the translational and rotational accelerations of the body-frame, andboth of them are expressed in F 𝐵 . 𝒂 𝑒 = [ (cid:165) 𝑟 𝑒,𝑥 , (cid:165) 𝑟 𝑒,𝑦 , (cid:165) 𝑟 𝑒,𝑧 ] T is the relative deformation acceleration of dm. From Eq. (7),the acceleration component of dm that is purely induced by rigid-body motion equals 𝒂 𝑅𝑤 = 𝑪 𝑊 𝐵 𝒂 𝑓 + 𝑪 𝑊 𝐵 ˜ 𝜶 𝑓 ( 𝒓 𝑤𝑏 + 𝑪 T 𝑊 𝐵 𝒓 𝑤 ) + 𝑪 𝑊 𝐵 ˜ 𝝎 𝑏 ˜ 𝝎 𝑏 ( 𝒓 𝑤𝑏 + 𝑪 T 𝑊 𝐵 𝒓 𝑤 ) (8)Integrate Eq. (8) at the wing section 𝐴 𝑤 , then the inertial force purely induced by rigid-body motion per unit lengthis computed as 𝒇 acc = − ∬ 𝐴 𝑤 𝜌 𝒂 𝑅𝑤 d 𝑟 𝑦 d 𝑟 𝑧 (9)where 𝜌 is the volume density. The inertial moment around 𝑂 𝑤 𝑥 𝑤 that is purely induced by rigid-body motion is 𝒎 acc = − ∬ 𝐴 𝑤 𝜌 [ , 𝑟 𝑦 , 𝑟 𝑧 ] T × 𝒂 𝑅𝑤 d 𝑟 𝑦 d 𝑟 𝑧 (10)The distributed inertial forces and moments are integrated to their nearest structural node. Consider the 𝑖 -th structuralnode with a spanwise location 𝑟 𝑥 = 𝑥 𝑖 , then the integrated force vector associated with its degrees of freedom is 𝑓 𝑤 𝑓 𝜙 𝑓 𝜃 𝑓 𝛽 acc ,𝑖 = ∫ 𝑥 𝑖 𝑥 𝑖 − 𝑓 acc ,𝑧 d 𝑟 𝑥 ∫ 𝑥 𝑖 𝑥 𝑖 − ( 𝑥 𝑖 − 𝑥 ) 𝑓 acc ,𝑧 − 𝑚 acc ,𝑥 d 𝑟 𝑥 ∫ 𝑥 𝑖 𝑥 𝑖 − 𝑚 acc ,𝑦 d 𝑟 𝑥 (11)which is valid for 𝑖 = , .., 𝑛 𝑠 . The inertial loads on the structural nodes are collected into a force vector 𝑭 s,acc = [[ 𝑓 𝑤 , 𝑓 𝜙 , 𝑓 𝜃 , 𝑓 𝛽 ] acc , , ..., [ 𝑓 𝑤 , 𝑓 𝜙 , 𝑓 𝜃 , 𝑓 𝛽 ] acc ,𝑛 𝑠 ] T ∈ R 𝑛 𝑠 .
3. Aerodynamic Forces
The unsteady strip theory is used to calculate the aerodynamic force 𝑭 aero , which is caused by motions, flapdeflections and external atmospheric disturbances. Discretize the wing into 𝑛 𝑎 undeformable strips, where each one ofthem has three degrees of freedom: heave 𝑤 𝑎 (downwards positive), pitching around the elastic axis 𝜃 𝑎 (nose-up positive)and flap deflection 𝛽 𝑎 (downwards positive). Referring to Theodorsen’s theory, 𝑭 aero contains instant noncirculatoryand time-dependent circulatory parts [18]. The influences of noncirculatory force can be directly considered in theaerodynamic mass, damping and stiffness matrices. On the other hand, the circulatory part was initially expressed infrequency domain [18]. For control design purposes, the circulatory aerodynamics are converted into time-domainusing indicial function approximation.Define a nondimensional time variable as 𝜏 = 𝑉𝑡 / 𝑏 , where 𝑏 is the semi-chord. Using the Duhamel’s integral, adynamic system with indicial response function 𝑓 ( 𝜏 ) = − 𝑎 𝑒 − 𝑏 𝜏 − 𝑎 𝑒 − 𝑏 𝜏 can be realized in a state-space form as (cid:34) (cid:164) 𝑧 𝑎 (cid:164) 𝑧 𝑎 (cid:35) = (cid:34) − (cid:0) 𝑉𝑏 (cid:1) 𝑏 𝑏 − (cid:0) 𝑉𝑏 (cid:1) ( 𝑏 + 𝑏 ) (cid:35) (cid:34) 𝑧 𝑎 𝑧 𝑎 (cid:35) + (cid:34) (cid:35) 𝑢𝑦 = (cid:34) ( 𝑎 + 𝑎 ) 𝑏 𝑏 (cid:18) 𝑉𝑏 (cid:19) , ( 𝑎 𝑏 + 𝑎 𝑏 ) (cid:18) 𝑉𝑏 (cid:19)(cid:35) (cid:34) 𝑧 𝑎 𝑧 𝑎 (cid:35) + ( − 𝑎 − 𝑎 ) 𝑢 (12)The circulatory lift induced by aircraft motions and flap deflections is modeled by the Wagner’s function, which hasan exponential approximation 𝜙 ( 𝜏 ) = − . 𝑒 − . 𝜏 − . 𝑒 − . 𝜏 [19]. Since this function is in the form of 𝑓 ( 𝜏 ) ,the state-space representation in Eq. (12) can be applied. The input of the system is chosen as 𝑢 = 𝛼 𝑞𝑠,𝑟 ( 𝑽 𝑓 , 𝝎 𝑏 ) + 𝛼 𝑞𝑠,𝑎𝑒 + 𝛽 𝑞𝑠 (13)5here the subscript (·) 𝑞𝑠 indicates quasi-steady. From Eq. (6), the local in-flow velocity that is purely caused byrigid-body motion equals 𝑽 𝑅𝑤 = 𝑪 𝑊 𝐵 𝑽 𝑏 + 𝑪 𝑊 𝐵 ˜ 𝝎 𝑏 ( 𝒓 𝑤𝑏 + 𝑪 T 𝑊 𝐵 𝒓 𝑤 ) . Therefore, the local quasi-steady angle of attackinduced by rigid-body translational and rotational motion equals 𝛼 𝑞𝑠,𝑟 ( 𝑽 𝑏 , 𝝎 𝑏 ) = arctan ( 𝑉 𝑅𝑤,𝑧 / 𝑉 𝑅𝑤,𝑥 ) .The angle of attack induced by aeroelastic motions is 𝛼 𝑞𝑠,𝑎𝑒 = (cid:164) 𝑤 𝑎 𝑉 + 𝜃 𝑎 + 𝑏 ( − 𝑎 ) (cid:164) 𝜃 𝑎 𝑉 . Moreover, in Eq. (13), 𝛽 𝑞𝑠 = 𝑇 𝜋 𝛽 𝑎 + 𝑇 𝜋 𝑏𝑉 (cid:164) 𝛽 𝑎 , where 𝑇 and 𝑇 can be found in [18]. Substituting Eq. (13) and the parameters of the Wagner’sfunction into Eq. (12), then the output 𝑦 gives the effective time-dependent angle of attack. As a result, the localcirculatory lift coefficient equals 𝐶 𝑆𝐹𝐿 𝛼 𝑦 . The Theodorsen’s equations are derived for thin airfoils with lift-slope equals2 𝜋 . In this paper, more general expressions are used, where the local 𝐶 𝑆𝐹𝐿 𝛼 is obtained from steady-flow analysis basedon the vortex lattice method. This way allows the considerations of the wing-tip wake roll up effects and compressibility(Prandtl-Glauert factor) in the local 𝐶 𝑆𝐹𝐿 𝛼 . Consequently, the resulting circulatory moment coefficient around the shearcenter equals 𝐶 𝑆𝐹𝐿 𝛼 𝑦 ( . + . 𝑎 ) , and the circulatory hinge moment coefficient equals − 𝐶 𝑆𝐹𝐿 𝛼 𝑦 𝑇 𝜋 (the coefficient 𝑇 isgiven in [18]).The Küssner function is used to calculate the unsteady responses of an airfoil to an unit sharp-edged gust. One of itsexponentially response approximation is 𝜓 𝑔 ( 𝜏 ) = − . 𝑒 − . 𝜏 − . 𝑒 − 𝜏 [19]. This expression is also in the form of 𝑓 ( 𝜏 ) = − 𝑎 𝑒 − 𝑏 𝜏 − 𝑎 𝑒 − 𝑏 𝜏 , thus the dynamic system in Eq. (12) can be used to calculate the responses of an airfoil toan arbitrary-spectrum atmospheric disturbance. Consider a wing section encounters a vertical gust with velocity 𝑤 𝑔 ( 𝑡 ) ,then the gust velocity projected in the wing reference frame equals 𝑽 gust = 𝑪 𝑊 𝐵 𝑪 𝐵𝐼 [ , , 𝑤 𝑔 ( 𝑡 )] T . Therefore, the inputof Eq. (12) is taken as 𝑢 = 𝛼 𝑔 = arctan (− 𝑉 gust ,𝑧 / 𝑉 ) . Substituting the parameters of the Küssner function into Eq. (12),then the output 𝑦 gives the effective time-dependent angle of attack induced by 𝑤 𝑔 ( 𝑡 ) . The resulting lift, pitchingmoment and hinge moment coefficients are given by 𝐶 𝑆𝐹𝐿 𝛼 𝑦 multiplied with 1, 0 . + . 𝑎 and − 𝑇 / 𝜋 , respectively.Based on the above discussions, four additional states are needed by each strip to model the circulatory dynamics.These additional states are also known as lag states. Define the lag state vector as 𝒛 𝑎 ∈ R 𝑛 𝑎 , and define the aerodynamicstate vector as 𝒙 𝑎 ∈ R 𝑛 𝑎 , then the following state-space system can be derived (cid:164) 𝒛 𝑎 = 𝑨 𝑧 𝒛 𝑎 + 𝑩 𝑧 𝒙 𝑎 + 𝑩 (cid:164) 𝑧 (cid:164) 𝒙 𝑎 + 𝑩 𝑧𝑟 𝜶 𝑞𝑠,𝑟 + 𝑩 𝑧𝑔 𝜶 𝑔 (14)where 𝜶 𝑞𝑠,𝑟 ∈ R 𝑛 𝑎 and 𝜶 𝑔 ∈ R 𝑛 𝑎 are the angle of attack vectors induced by rigid-body motions and gusts. Define theaerodynamic output vector as 𝒇 aero = [ 𝒇 T aero , , ..., 𝒇 T aero ,𝑛 𝑎 ] T , where each element 𝒇 aero ,𝑖 contains lift, pitching momentaround the elastic axis, and the hinge moment, then 𝒇 aero = 𝑴 𝑎 (cid:165) 𝒙 𝑎 + 𝑪 𝑎 (cid:164) 𝒙 𝑎 + 𝑲 𝑎 𝒙 𝑎 + 𝑲 𝑧 𝒛 + 𝑲 𝑧𝑟 𝜶 𝑞𝑠,𝑟 (15)In the above equation, 𝑴 𝑎 , 𝑪 𝑎 , and 𝑲 𝑎 are the aerodynamic mass, damping and stiffness matrices, and they includeboth the noncirculatory and circulatory effects. Because in the Wagner’s function, 1 − 𝑎 − 𝑎 = .
5, so the inputof 𝜶 𝑞𝑠,𝑟 has direct influence on 𝒇 aero . On contrast, 1 − 𝑎 − 𝑎 = 𝜶 𝑔 only indirectlyinfluence 𝒇 aero though the lag state 𝒛 (Eq. (14)).To establish a coupled aeroelastic model, the degrees of freedom of the aerodynamic strips and the structure nodesneed to be connected. Normally, the number of strips 𝑛 𝑎 is larger than the number of structure nodes 𝑛 𝑠 . Two linearinterpolation matrices 𝑯 𝑎𝑠 ∈ R 𝑛 𝑎 × 𝑛 𝑠 and 𝑯 𝑠𝑎 ∈ R 𝑛 𝑠 × 𝑛 𝑎 are designed based on the nearest neighbor principle. Theaerodynamic degrees of freedom are interpolated from the structural states as 𝒙 𝑎 = 𝑯 𝑎𝑠 𝒙 𝑠 . On the other hand, thedistributed aerodynamic forces are assigned to the structural nodes as 𝑭 𝑠, aero = 𝑯 𝑠𝑎 𝒇 aero .
4. Gravitational Forces
The gravitational force per unit length expressed in the right-wing reference frame is 𝒇 grav = ∬ 𝐴 𝑤 𝑪 𝑊 𝐵 𝑪 𝐵𝐼 [ , , 𝜌𝑔 ] T d 𝑟 𝑥 d 𝑟 𝑧 (16)where 𝑔 is the gravitational acceleration. The gravitational moment around 𝑂 𝑤 𝑥 𝑤 per unit length is 𝒎 grav = ∬ 𝐴 𝑤 ([ , 𝑟 𝑦 , 𝑟 𝑧 ] T + 𝒓 𝑒 ) × 𝑪 𝑊 𝐵 𝑪 𝐵𝐼 [ , , 𝜌𝑔 ] T d 𝑟 𝑥 d 𝑟 𝑧 (17)The distributed gravitational forces and moments are integrated to the nearest structural node. Consider the 𝑖 -th6tructural node with spanwise location 𝑟 𝑥 = 𝑥 𝑖 , then the integrated force vector associated with its degrees of freedom is 𝑓 𝑤 𝑓 𝜙 𝑓 𝜃 𝑓 𝛽 grav ,𝑖 = ∫ 𝑥 𝑖 𝑥 𝑖 − 𝑓 grav ,𝑧 d 𝑟 𝑥 ∫ 𝑥 𝑖 𝑥 𝑖 − ( 𝑥 𝑖 − 𝑥 ) 𝑓 grav ,𝑧 − 𝑚 grav ,𝑥 d 𝑟 𝑥 ∫ 𝑥 𝑖 𝑥 𝑖 − 𝑚 grav ,𝑦 d 𝑟 𝑥 (18)which is valid for 𝑖 = , .., 𝑛 𝑠 . The gravitational loads at structural nodes are collected into a force vector 𝑭 s,grav = [[ 𝑓 𝑤 , 𝑓 𝜙 , 𝑓 𝜃 , 𝑓 𝛽 ] grav , , ..., [ 𝑓 𝑤 , 𝑓 𝜙 , 𝑓 𝜃 , 𝑓 𝛽 ] grav ,𝑛 𝑠 ] T ∈ R 𝑛 𝑠 .
5. Actuation Forces
In this paper, flaps are attached to structural nodes via rotational springs. In Eq. (3), the flap stiffness is consideredin 𝑲 𝑠 , while the flap mass and its inertial couplings with the beam structure are modeled in 𝑴 𝑠 . The actuationmoment around the hinge is chosen as the control input, which is more physically meaningful than the flap angleitself. Denote the actuation moment for the 𝑖 -th flap as 𝑢 𝑖 , then the total actuation force vector is written as 𝑭 𝑠, act = [[ , , , 𝑢 ] , ..., [ , , , 𝑢 𝑛 𝑠 ]] T . Denote the control input vector as 𝒖 = [ 𝑢 , ..., 𝑢 𝑛 𝑠 ] T ∈ R 𝑛 𝑠 , then 𝑭 𝑠, act alsoequals 𝑯 𝑢 𝒖 , where 𝑯 𝑢 ∈ R 𝑛 𝑠 × 𝑛 𝑠 is a boolean selection matrix.
6. Wing Aeroservoelastic Dynamic Equations
After the derivations of inertial, aerodynamic, gravitational, and actuation forces, the wing aeroservoelastic model isready to be assembled. Because the origin of the wing-fixed reference frame F 𝑊 is coincide with the first structure node,when observing the wing motions in F 𝑊 , we have (cid:165) 𝒙 𝑠, = (cid:164) 𝒙 𝑠, = 𝒙 𝑠, =
0. The 𝑭 ext in Eq. (3) equals the summation of 𝑭 𝑠, aero , 𝑭 𝑠, acc , 𝑭 𝑠, grav , and 𝑭 𝑠, act . Substituting Eqs. (14) and (15) into Eq. (3), the following state-space system can bederived: (cid:165) 𝒙 𝑠 (cid:164) 𝒙 𝑠 (cid:164) 𝒛 𝑎 = 𝑴 − 𝑎𝑒 𝑯 𝑠𝑎 𝑪 𝑎 𝑯 𝑎𝑠 𝑴 − 𝑎𝑒 ( 𝑯 𝑠𝑎 𝑲 𝑎 𝑯 𝑎𝑠 − 𝑲 𝑠,𝑛 𝑠 ) 𝑴 − 𝑎𝑒 𝑯 𝑠𝑎 𝑲 𝑧 𝑰 𝑛 𝑠 × 𝑛 𝑠 𝑛 𝑠 × 𝑛 𝑠 𝑛 𝑠 × 𝑛 𝑎 𝑩 (cid:164) 𝑧 𝑯 𝑎𝑠 𝑩 𝑧 𝑯 𝑎𝑠 𝑨 𝑧 (cid:164) 𝒙 𝑠 𝒙 𝑠 𝒛 𝑎 + 𝑴 − 𝑎𝑒 𝑯 𝑢 𝑛 𝑠 × 𝑛 𝑎 × 𝒖 (19) + 𝑛 𝑠 × 𝑛 𝑎 𝑛 𝑠 × 𝑛 𝑎 𝑩 𝑧𝑔 𝜶 𝑔 + 𝑴 − 𝑎𝑒 𝑭 𝑠, grav 𝑛 𝑠 × 𝑛 𝑎 × + 𝑴 − 𝑎𝑒 𝑯 𝑠𝑎 𝑲 𝑧𝑟 𝑛 𝑠 × 𝑛 𝑎 𝑩 𝑧𝑟 𝜶 𝑞𝑠,𝑟 + 𝑴 − 𝑎𝑒 𝑭 𝑠, acc 𝑛 𝑠 × 𝑛 𝑎 × in which 𝑴 𝑎𝑒 = 𝑴 𝑠,𝑛 𝑠 − 𝑯 𝑠𝑎 𝑴 𝑎 𝑯 𝑎𝑠 . Recall Eq. (3), the wing root reaction force is calculated as 𝑭 root = 𝑴 𝑠, (cid:165) 𝒙 𝑠 + 𝑲 𝑠, 𝒙 𝑠 = 𝑴 𝑠, 𝑴 − 𝑎𝑒 𝑯 𝑠𝑎 𝑲 𝑧𝑟 𝜶 𝑞𝑠,𝑟 + 𝑴 𝑠, 𝑴 − 𝑎𝑒 ( 𝑩 𝑢 𝒖 + 𝑭 𝑠, acc + 𝑭 𝑠, grav )+ (cid:104) 𝑴 𝑠, 𝑴 − 𝑎𝑒 𝑯 𝑠𝑎 𝑪 𝑎 𝑯 𝑎𝑠 𝑴 𝑠, 𝑴 − 𝑎𝑒 ( 𝑯 𝑠𝑎 𝑲 𝑎 𝑯 𝑎𝑠 − 𝑲 𝑠,𝑛 𝑠 ) + 𝑲 𝑠, 𝑴 𝑠, 𝑴 − 𝑎𝑒 𝑯 𝑠𝑎 𝑲 𝑧 (cid:105) (cid:164) 𝒙 𝑠 𝒙 𝑠 𝒛 𝑎 (20)The above wing aeroservoelastic model is derived in a modular approach. The control command, atmosphericdisturbance, gravitational force, the angle of attack induced by rigid-body motions, and the inertial forces induced byrigid-body accelerations are modeled as separate inputs to the system. This allows convenient contribution assessmentsof different factors on the wing dynamic loads. Furthermore, this approach provides a clear interface between theconventional clamped-wing aeroservoelastic dynamics and the free-flying wing aeroservoelastic dynamics. If 𝑭 𝑠, acc and 𝜶 𝑞𝑠,𝑟 are set to zero, while the orientation matrices in 𝑭 𝑠, grav are evaluated at the trimmed condition, then Eqs. (19) and(20) degenerate to the conventional clamped-wing aeroservoelastic dynamic equations. On the other hand, no matterhow the aerodynamic and structural dynamics are modeled and coupled (finite-element method, panel method, etc.),as long as the resulting aeroelastic model can be written in state-space form, the clamped-wing model can be madefree-flying by adding 𝑭 𝑠, acc , 𝜶 𝑞𝑠,𝑟 and the direction cosine matrices.7 . Free-flying Dynamics of the Flexible Aircraft The translational dynamics of the flexible aircraft expressed in the flight trajectory axes F 𝑉 are (cid:164) 𝑉 (cid:164) 𝜒 (cid:164) 𝛾 = 𝑉 cos 𝛾
00 0 − 𝑉 𝑪 𝑉 𝐵 𝑚 ( 𝑭 𝐵 tot − ˜ 𝑺 T 𝜶 𝑏 − ˜ 𝝎 𝑏 ˜ 𝑺 T 𝝎 𝑏 + 𝒅 𝐹 ) (21)while the rotational dynamics of the aircraft expressed in the body-fixed frame F 𝐵 are 𝑱 (cid:164) 𝝎 𝑏 = − ˜ 𝑺 (cid:164) 𝑽 𝑏 − ˜ 𝑽 𝑏 ˜ 𝑺 T 𝝎 𝑏 − ˜ 𝝎 𝑏 ˜ 𝑺𝑽 𝑏 − ˜ 𝝎 𝑏 𝑱𝝎 𝑏 + 𝑴 𝐵 tot + 𝒅 𝑀 (22)Denote the distance vector from 𝑂 𝑏 to an infinitesimal mass element dm on the aircraft as 𝒓 , which is expressed in F 𝐵 . Then 𝑺 = ∫ 𝒓 d 𝑚 is the first moment of area, which is non-zero when 𝑂 𝑏 is not coincide with the center of mass.The moment of inertia is calculated as 𝑱 = ∫ ˜ 𝒓 T ˜ 𝒓 d 𝑚 . For a flexible aircraft, both 𝑺 and 𝑱 are functions of structuraldeformations. 𝑭 𝐵 tot and 𝑴 𝐵 tot are the total external forces and moments expressed in F 𝐵 , which contain aerodynamic, gravitationaland propulsion forces. 𝒅 𝐹 and 𝒅 𝑀 are functions of elastic accelerations and cross couplings between rigid-body andelastic velocities (Coriolis effects) [13, 20]. If 𝒅 𝐹 and 𝒅 𝑀 are set to zero, while 𝑺 and 𝑱 are kept constant, then therigid-body dynamics are retrieved. If in the resulting simplified dynamics, further set 𝑺 to zero, then the conventionalrigid-body dynamics where 𝑂 𝑏 is coincide with aircraft c.m. are retrieved. 𝒅 𝐹 , 𝒅 𝑀 and the dependencies of 𝑺 , 𝑱 onelastic states essentially reflect the inertial couplings between rigid-body and elastic motions.As compared to the inertial couplings, aerodynamic couplings normally have more dominant effects. This is alsothe only coupling aspect considered in the mean-axes method [17]. As discussed in subsection II.B, rigid-body motionsexcite wing aeroelastic dynamics through the input 𝜶 𝑞𝑠,𝑟 . On the other hand, 𝑭 𝐵 tot and 𝑴 𝐵 tot are contributed by 𝒇 aero (Eq. (15)), which is a function of aeroelastic states. It is noteworthy that drag is not included in 𝒇 aero , while it should beconsidered in 𝑭 𝐵 tot and 𝑴 𝐵 tot . For each aerodynamic strip, the local drag coefficient is calculated using the quadraticexpression 𝐶 𝐷 = 𝐶 𝐷 + 𝑘 𝐷 𝐶 𝐿 . Since the fuselage and tails are assumed to be rigid, the quasi-steady strip theory isused for calculating the distributed lift and drag on them. Denote the distance vector from 𝑂 𝑏 to an arbitrary strip as 𝒓 𝑖 (expressed in F 𝐵 ), then the local airspeed is derived as 𝑽 𝑎 𝑖 = 𝑽 𝑏 + ˜ 𝝎 𝑏 𝒓 𝑖 − 𝑪 𝐵𝐼 [ , , 𝑤 𝑔 ( 𝑡 )] T . During simulations, theturbulence velocity 𝑤 𝑔 ( 𝑡 ) is calculated by interpolating the spatial vertical atmospheric disturbance field at the localstrip position 𝑹 𝑏 + 𝑪 T 𝐵𝐼 𝒓 𝑖 [13]. From 𝑽 𝑎 𝑖 , the local dynamic pressure and aerodynamic angles can be calculated. III. Control Design
After presenting the flight dynamics, kinematics, and wing aeroservoelastic dynamics, the flexible aircraft trajectorytracking control architecture will be shown in this section. Cascaded control loops are widely adopted in flexible aircraftcontrol [21, 22]. The control architecture proposed in this paper contains four cascaded control loops: position control,flight path control, attitude control, and optimal multi-objective wing control. These control loops will be designed inthe following subsections:
A. Position Control
The kinematics of aircraft positions are given in Eq. (1). The control objective is to make the aircraft followa pre-designed three-dimensional trajectory [ 𝑋 ref , 𝑌 ref , 𝑍 ref ] T . Since there is no model uncertainty in the positionkinematics, the nonlinear dynamic inversion control is adopted. Design the virtual control for the lateral and verticalpositions as 𝜈 𝑌 = (cid:164) 𝑌 ref + 𝐾 𝑌 ( 𝑌 ref − 𝑌 ) , 𝜈 𝑍 = (cid:164) 𝑍 ref + 𝐾 𝑍 ( 𝑍 ref − 𝑍 ) (23)where 𝐾 𝑌 and 𝐾 𝑍 are positive proportional gains, then when (cid:164) 𝑌 = 𝜈 𝑌 , and (cid:164) 𝑍 = 𝜈 𝑍 , the lateral and vertical trackingerrors converge to zero exponentially. Replace (cid:164) 𝑌 and (cid:164) 𝑍 in Eq. (1) by their virtual controls, and then invert the resultingequations, the references for the kinematic azimuth angle and the flight path angles are obtained as 𝜒 ref = arcsin ( 𝜈 𝑌 /( 𝑉 cos 𝛾 )) , 𝛾 ref = − arcsin ( 𝜈 𝑍 / 𝑉 ) (24)8 . Flight Path Control The control objective of this loop it to make the flight path angles track their references designed in Eq. (24). Becausethe total external force vector contains aerodynamic, gravitational and propulsion forces, the aircraft translationaldynamics given in Eq. (21) can be rewritten as (cid:164) 𝑉 (cid:164) 𝜒 (cid:164) 𝛾 = 𝑉 cos 𝛾
00 0 − 𝑉 (cid:169)(cid:173)(cid:173)(cid:171) 𝑪 𝑉 𝐵 𝑚 𝑇 + 𝑪 𝑉 𝐼 𝑚 𝑚𝑔 + 𝑪 𝑉 𝐴 𝑚 − 𝐷𝐶 − 𝐿 + 𝑑 𝑉 𝑑 𝜒 𝑑 𝛾 (cid:170)(cid:174)(cid:174)(cid:172) (25)in which [ 𝑇, , ] T is the thrust vector. 𝐷 , 𝐶 and 𝐿 are the total drag, side force, and lift of the aircraft, which aredefined in the aerodynamic axes F 𝐴 . [ 𝑑 𝑉 , 𝑑 𝜒 , 𝑑 𝛾 ] T includes 𝒅 𝐹 in Eq. (21) and the cross-coupling terms caused by thenon-zero first moment of area. Substituting the expressions of the direction cosine matrices into Eq. (25), the dynamicsof the flight path angle are (cid:164) 𝛾 = − 𝑚𝑉 (− 𝑇 sin 𝛼 cos 𝜇 + 𝑚𝑔 cos 𝛾 − 𝐿 cos 𝜇 + 𝑑 𝛾 ) (26)From the physical point of view, one of the most effective ways of changing aircraft flight path is by changing thetotal lift, which is dominated by the angle of attack 𝛼 . Therefore, 𝛼 is selected as a control input to the dynamics of 𝛾 .Denote the flight path angle dynamics in Eq. (26) as 𝛾 = 𝑓 𝛾 ( 𝒙 , 𝑢 ) + 𝑑 𝛾 , where the vector 𝒙 include the rigid-body andelastic states except 𝛼 . It can be seen from Eq. (26) that 𝑓 𝛾 ( 𝒙 , 𝑢 ) + 𝑑 𝛾 is non-affine in control. Conventional nonlinearcontrol methods including feedback linearization and backstepping cannot be directly applied. In this paper, the novelincremental sliding mode control method is adopted, which not only applies to non-affine in control dynamic systems,but also has enhanced robustness to model uncertainties [14].Design the sliding surface as 𝜎 𝛾 = 𝛾 − 𝛾 ref . The control objective is to design a reference for 𝛼 , such that when 𝛼 tracks its reference through inner-loop control, 𝜎 𝛾 also converges to zero. Denote the sampling interval as Δ 𝑡 ,the incremental dynamic equation is derived by taking the first-order Taylor series expansion of Eq. (26) around thecondition at 𝑡 − Δ 𝑡 (denoted by the subscript 0) as (cid:164) 𝛾 = 𝑓 𝛾 ( 𝒙 , 𝑢 ) + 𝑑 𝛾 = (cid:164) 𝛾 + 𝜕 𝑓 𝛾 𝜕𝑢 (cid:12)(cid:12)(cid:12)(cid:12) Δ 𝑢 + 𝜕 𝑓 𝛾 𝜕 𝒙 (cid:12)(cid:12)(cid:12)(cid:12) Δ 𝒙 + Δ 𝑑 𝛾 + 𝑅 (27)in which Δ 𝒙 and Δ 𝑢 respectively represents the variations of 𝒙 and 𝑢 in one sampling time step Δ 𝑡 . Δ 𝑑 𝛾 is the variationsof 𝑑 𝛾 in Δ 𝑡 . 𝑅 is the expansion remainder, whose Lagrange form is 𝑅 = 𝜕 𝑓 𝛾 𝜕 𝒙 (cid:12)(cid:12)(cid:12)(cid:12) 𝑚 Δ 𝒙 + 𝜕 𝑓 𝛾 𝜕 𝒙 𝜕𝑢 (cid:12)(cid:12)(cid:12)(cid:12) 𝑚 Δ 𝒙 Δ 𝑢 + 𝜕 𝑓 𝛾 𝜕 𝑢 (cid:12)(cid:12)(cid:12)(cid:12) 𝑚 Δ 𝑢 (28)in which (·)| 𝑚 means evaluating (·) at a condition where 𝒙 ∈ ( 𝒙 ( 𝑡 − Δ 𝑡 ) , 𝒙 ( 𝑡 )) , 𝑢 ∈ ( 𝑢 ( 𝑡 − Δ 𝑡 ) , 𝑢 ( 𝑡 )) . The incrementalsliding mode control law for stabilizing 𝜎 𝛾 is then designed as Δ 𝑢 = ( 𝜈 𝑛 + 𝜈 𝑠 + (cid:164) 𝛾 ref − (cid:164) 𝛾 )/ ¯ 𝐺 (29)where the virtual control 𝜈 𝑛 is designed to stabilize the nominal sliding dynamics, while the virtual control 𝜈 𝑠 is arobustify virtual control for disturbance rejection. The specific expressions of these two virtual control terms will bepresented later. The real control input equals Δ 𝑢 + 𝑢 , where 𝑢 denotes the control input at 𝑡 − Δ 𝑡 . In Eq. (29), ¯ 𝐺 represents the estimation of the control effectiveness 𝜕 𝑓 𝛾 𝜕𝑢 (cid:12)(cid:12) . It can be seen from Eq. (26) that 𝐺 = 𝜕 𝑓 𝛾 𝜕𝑢 (cid:12)(cid:12)(cid:12)(cid:12) = cos 𝜇𝑚𝑉 (cid:18) 𝑇 cos 𝛼 + 𝜕𝐿𝜕𝛼 (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) (30)its estimation is calculated as ¯ 𝐺 = cos 𝜇𝑚𝑉 (cid:0) 𝑇 cos 𝛼 + 𝑞 ∞ 𝑆 𝑤 ¯ 𝐶 𝐿𝛼 (cid:1) (cid:12)(cid:12) (31)where 𝑞 ∞ is the dynamic pressure and 𝑆 𝑤 is the wing area. Substituting Eq. (29) into Eq. (27), then the closed-loopdynamics are (cid:164) 𝜎 𝛾 = (cid:164) 𝛾 − (cid:164) 𝛾 ref = (cid:164) 𝛾 + ( 𝐺 / ¯ 𝐺 )( 𝜈 𝑛 + 𝜈 𝑠 + (cid:164) 𝛾 ref − (cid:164) 𝛾 ) − (cid:164) 𝛾 ref + 𝛿 ( 𝒙 , Δ 𝑡 ) + Δ 𝑑 𝛾 (cid:44) 𝜈 𝑛 + 𝜈 𝑠 + 𝜀 𝛾 (32)9here 𝛿 ( 𝒙 , 𝑡 ) is the closed-loop value of the state variations and the expansion reminder: 𝛿 ( 𝒙 , 𝑡 ) = (cid:20) 𝜕 𝑓 𝛾 𝜕 𝒙 (cid:12)(cid:12)(cid:12)(cid:12) Δ 𝒙 + 𝑅 (cid:21) 𝑢 = 𝑢 + Δ 𝑢 (33) 𝜀 𝛾 in Eq. (32) is the lumped perturbation term, which is expressed as 𝜀 𝛾 = 𝛿 ( 𝒙 , Δ 𝑡 ) + Δ 𝑑 𝛾 + ( 𝐺 / ¯ 𝐺 − )( 𝜈 𝑛 + 𝜈 𝑠 + (cid:164) 𝛾 ref − (cid:164) 𝛾 ) (34) Assumption 1 𝛿 ( 𝒙 , Δ 𝑡 ) in Eq. (33) and Δ 𝑑 𝛾 in Eq. (27) satisfy | 𝛿 ( 𝒙 , Δ 𝑡 )| < ¯ 𝛿 , | Δ 𝑑 𝛾 | < ¯ 𝑑 𝛾 . Since 𝒙 is continuously differentiable, then lim Δ 𝑡 → (cid:107) Δ 𝒙 (cid:107) =
0. If the partial derivatives of 𝑓 𝛾 with respect to 𝒙 ,up to any order, are bounded, then in view of Eqs. (28) and (33), the absolute value of 𝛿 ( 𝒙 , Δ 𝑡 ) approaches zero as Δ 𝑡 decreases. Therefore, this assumption holds when Δ 𝑡 is sufficiently small. Proposition 1
Under Assumption 1, if < 𝐺 / ¯ 𝐺 < , then when the sampling frequency 𝑓 𝑠 is sufficient high, theuncertainty term 𝜀 𝛾 in Eq. (34) is ultimately bounded. Proof:
Substituting Eq. (29) into Eq. (27), the resulting closed-loop dynamics are (cid:164) 𝛾 = 𝜈 𝑛 + 𝜈 𝑠 + (cid:164) 𝛾 ref + 𝜀 𝛾 . Thisclosed-loop dynamics are also valid at 𝑡 − Δ 𝑡 , thus (cid:164) 𝛾 = ( 𝜈 𝑛 + 𝜈 𝑠 + (cid:164) 𝛾 ref )| + 𝜀 𝛾, . Therefore, 𝜀 𝛾 in Eq. (34) can berewritten as 𝜀 𝛾 = 𝛿 ( 𝒙 , Δ 𝑡 ) + Δ 𝑑 𝛾 + ( 𝐺 / ¯ 𝐺 − )( 𝜈 𝑛 + 𝜈 𝑠 + (cid:164) 𝛾 ref − (( 𝜈 𝑛 + 𝜈 𝑠 + (cid:164) 𝛾 ref )| + 𝜀 𝛾, )) = ( − 𝐺 / ¯ 𝐺 ) 𝜀 𝛾, − ( − 𝐺 / ¯ 𝐺 )( 𝜈 𝛾 − 𝜈 𝛾, ) + 𝛿 ( 𝒙 , Δ 𝑡 ) + Δ 𝑑 𝛾 (35)where 𝜈 𝛾 = 𝜈 𝑛 + 𝜈 𝑠 + (cid:164) 𝛾 ref . These virtual control terms are all continuous in time, thus under sufficiently high samplingfrequency, the variation of 𝜈 in Δ 𝑡 has an upper bound, i.e., | 𝜈 𝛾 − 𝜈 𝛾, | < ¯ 𝜈 𝛾 . Analogous to the proof of Theorem 1in [14], the ultimate boundedness of 𝜀 𝛾 can be proved. (cid:3) The boundedness of 𝜀 𝛾 makes it feasible to compensate for its influences using robust control. In this paper, 𝜈 𝑠 in Eq. (29) is designed using a super-twisting sliding mode observer. Design an auxiliary sliding variable as 𝑠 𝛾 = 𝜎 𝛾 − ∫ 𝜈 𝑛 d 𝑡 , then using Eq. (32), the dynamics of 𝑠 𝛾 are (cid:164) 𝑠 𝛾 = (cid:164) 𝜎 𝛾 − 𝜈 𝑛 = 𝜈 𝑠 + 𝜀 𝛾 (36)Further design the observer virtual control as 𝜈 𝑠 = − 𝜆 𝛾 | 𝑠 | . sign ( 𝑠 ) − 𝛽 𝛾 ∫ sign ( 𝑠 ) d 𝑡 (37)Denote the upper bound of 𝜀 𝛾 as ¯ 𝜀 𝛾 , and design 𝜆 𝛾 = . 𝜀 . 𝛾 , 𝛽 = . 𝜀 𝛾 , then 𝑠 = (cid:164) 𝑠 = 𝜈 𝑠 provides a real-time observation of the uncertainty term − 𝜀 𝛾 . Moreover, thenominal dynamics (cid:164) 𝜎 𝛾 = 𝜈 𝑛 are retrieved on the sliding surface regardless of uncertainties. Simply design the nominalvirtual control as 𝜈 𝑛 = − 𝐾 𝜎 𝜎 𝛾 , with 𝐾 𝜎 >
0, then 𝜎 𝛾 converges to zero exponentially. The observer gains can also bemade adaptive, which removes the per-knowledge requirement on the uncertainty bound [24]. Remark 1
The only model information required by the incremental sliding mode control in Eq. (29) is the controleffectiveness 𝐺 . In view of Eq. (31), the model parameters in 𝐺 are only the wing area 𝑆 𝑤 and the aircraft lift-slope 𝐶 𝐿𝛼 . These two parameters are easily known from aircraft overall design. Remark 2
The presented flight path control differs from the one in [25] in three aspects: 1) the residual uncertaintyterm 𝜀 𝛾 is observed and is compensated for in this paper. 2) In Eq. (25), the aerodynamic forces are expressed inthe aerodynamic axes, while in [25], they are expressed in the body axes. Consequently, the aerodynamic coefficientused in the proposed controller is only 𝐶 𝐿𝛼 , while Ref. [25] requires both 𝐶 𝑋 𝛼 and 𝐶 𝑍 𝛼 ( 𝑋, 𝑍 are the aerodynamicforces expressed in the body-fixed frame). As compared to 𝐶 𝑋 𝛼 and 𝐶 𝑍 𝛼 , the state-dependency of 𝐶 𝐿𝛼 is lower, whichsimplifies its identification process. 3) In Ref. [25], throttle control is also included in the incremental control loop.However, aircraft throttle normally has much lower bandwidth than the aerodynamic control surfaces. In view of this, aseparate throttle controller ( 𝛿 𝑇 ) for maintaining airspeed is adopted in this paper.10fter presenting the flight path angle ( 𝛾 ) control law, the reference tracking problem for the kinematic azimuth angle 𝜒 will be solved. Recall Eq. (25), the dynamics of 𝜒 are (cid:164) 𝜒 = 𝑉 cos 𝛾 ( 𝑇 sin 𝛼 cos 𝜇 + 𝐿 sin 𝜇 + 𝑑 𝜒 ) (38)In view of these dynamics, the kinematic bank angle 𝜇 is chosen as the control input. Replace (cid:164) 𝛾 in Eq. (26) and (cid:164) 𝜒 inEq. (38) by their virtual control 𝜈 𝛾 and 𝜈 𝜒 respectively, and then invert the resulting dynamics, the reference for 𝜇 isdesigned as 𝜇 ref = arctan (cid:18) 𝜈 𝜒 𝑉 cos 𝛾𝜈 𝛾 𝑉 + 𝑔 cos 𝛾 (cid:19) (39)in which 𝜈 𝛾 has been designed in the preceding texts, while the virtual control 𝜈 𝜒 is designed as 𝜈 𝜒 = (cid:164) 𝜒 ref + 𝐾 𝜒 ( 𝜒 ref − 𝜒 ) with 𝐾 𝜒 > C. Attitude Control
The flexible aircraft attitude kinematics are given in Eq. (2), while the dynamics of angular rates are given inEq. (22). The control objective of this attitude control loop is to make 𝜇, 𝛼 and 𝛽 track their references. 𝜇 ref and 𝛼 ref have been designed in subsection III.B, while the 𝛽 ref is set to zero for mitigating side force. From a physical point ofview, the attitude of an aircraft can be changed by creating control moments around 𝑂 𝑏 𝑥 𝑏 , 𝑂 𝑏 𝑦 𝑏 and 𝑂 𝑏 𝑧 𝑏 axes. Forthe aircraft configuration shown in Fig. 1, the pitch and yaw control moments can be generated by the deflections ofelevator 𝛿 𝑒 and rudder 𝛿 𝑟 . On the other hand, the difference between the left and right wing root bending moments isessentially the main source of aircraft rolling moment. Therefore, 𝛿 𝑒 , 𝛿 𝑟 and the wing root bending moment difference 𝑀 𝜙, diff = 𝑀 𝜙,𝑙 − 𝑀 𝜙,𝑟 are chosen as control input variables in this loop. The wing root bending moment can bemeasured by strain gauges at a sampling frequency around 1 kHz.Define the state vectors as 𝒙 = [ 𝜇, 𝛼, 𝛽 ] T , 𝒙 = [ 𝑝, 𝑞, 𝑟 ] T , and define the control vector as 𝒖 = [ 𝛿 𝑒 , 𝛿 𝑟 , 𝑀 𝜙, diff ] T ,then the attitude kinematics and dynamics in Eqs. (2) and (22) can be represented as (cid:164) 𝒙 = 𝒇 ( 𝒙 ) + 𝑮 ( 𝒙 ) 𝒙 (cid:164) 𝒙 = 𝒇 ( 𝒙 , 𝒙 , 𝒙 𝑒 ) + 𝑮 ( 𝒙 , 𝒙 , 𝒙 𝑒 , 𝒖 ) + 𝒅 (40)where 𝒅 = 𝑱 − 𝒅 𝑀 , while 𝒙 𝑒 = [ (cid:164) 𝒙 𝑠 , 𝒙 𝑠 , 𝒛 𝑎 ] T represents the state vector in Eq. (19). Choose the controlled output as 𝒚 = 𝒙 , then the control objective is to make 𝒚 track its reference 𝒚 ref . In view of the kinematics in Eq. (2), there isno model uncertainty in 𝒇 and 𝑮 . On the contrary, 𝒇 and 𝑮 contains model uncertainties, and are functions ofelastic states. In addition, 𝑮 in non-affine in control, which makes direct applications of feedback linearization andbackstepping infeasible.The novel incremental backstepping sliding mode control proposed in [15] is adopted to handle this output trackingproblem. Define the error variable as 𝒛 = 𝒚 − 𝒚 ref , then (cid:164) 𝒛 = 𝒇 + 𝑮 𝒙 − (cid:164) 𝒚 ref (41)Consider a candidate Lyapunov function 𝑉 ( 𝒛 ) = 𝒛 T 𝒛 . In order to make (cid:164) 𝑉 ( 𝒛 ) ≤
0, the reference for 𝒙 isdesigned as 𝒙 ref2 = 𝑮 − (− 𝒇 − 𝑲 𝒛 + (cid:164) 𝒚 ref ) (42)in which 𝑲 is a positive definite diagonal gain matrix. Further define the tracking error for 𝒙 as 𝒛 = 𝒙 − 𝒙 ref2 , thenthe resulting dynamics are (cid:164) 𝒛 = 𝒇 ( 𝒙 , 𝒙 , 𝒙 𝑒 ) + 𝑮 ( 𝒙 , 𝒙 , 𝒙 𝑒 , 𝒖 ) + 𝒅 − (cid:164) 𝒙 ref2 . Denote the augmented state vector as 𝒙 = [ 𝒙 T , 𝒙 T , 𝒙 T 𝑒 ] T . The incremental dynamics of 𝒙 are derived by taken the first-order Taylor series expansion aroundthe condition at 𝑡 − Δ 𝑡 (denoted by the subscript 0) as: (cid:164) 𝒙 = (cid:164) 𝒙 , + 𝜕 [ 𝒇 ( 𝒙 ) + 𝑮 ( 𝒙 , 𝒖 )] 𝜕 𝒖 (cid:12)(cid:12)(cid:12)(cid:12) Δ 𝒖 + 𝜕 [ 𝒇 ( 𝒙 ) + 𝑮 ( 𝒙 , 𝒖 )] 𝜕 𝒙 (cid:12)(cid:12)(cid:12)(cid:12) Δ 𝒙 + Δ 𝒅 + 𝑹 (cid:48) (43)where Δ 𝒙 = 𝒙 − 𝒙 , Δ 𝒖 = 𝒖 − 𝒖 , and Δ 𝒅 = 𝒅 − 𝒅 respectively denotes the variations of states, control inputs, anddisturbances, in one incremental time step. 𝑹 (cid:48) in Eq. (43) is the expansion remainder, whose Lagrange form is 𝑹 (cid:48) = 𝜕 [ 𝒇 ( 𝒙 ) + 𝑮 ( 𝒙 , 𝒖 )] 𝜕 𝒙 (cid:12)(cid:12)(cid:12)(cid:12) 𝑚 Δ 𝒙 + 𝜕 [ 𝒇 ( 𝒙 ) + 𝑮 ( 𝒙 , 𝒖 )] 𝜕 𝒙 𝜕 𝒖 (cid:12)(cid:12)(cid:12)(cid:12) 𝑚 Δ 𝒙 Δ 𝒖 + 𝜕 [ 𝒇 ( 𝒙 ) + 𝑮 ( 𝒙 , 𝒖 )] 𝜕 𝒖 (cid:12)(cid:12)(cid:12)(cid:12) 𝑚 Δ 𝒖 (44)11n which (·)| 𝑚 means evaluating (·) at a condition where 𝒙 ∈ ( 𝒙 ( 𝑡 − Δ 𝑡 ) , 𝒙 ( 𝑡 )) , 𝒖 ∈ ( 𝒖 ( 𝑡 − Δ 𝑡 ) , 𝒖 ( 𝑡 )) , 𝒅 ∈ ( 𝒅 ( 𝑡 − Δ 𝑡 ) , 𝒅 ( 𝑡 )) .To stabilize 𝒛 , the incremental backstepping sliding mode control input is designed as Δ 𝒖 ibsmc = ¯ 𝑮 − ( 𝝂 𝑐 + 𝝂 𝑠 − (cid:164) 𝒙 , ) (45)where ¯ 𝑮 is the estimation of the control effectiveness matrix 𝑮 = 𝜕 [ 𝒇 ( 𝒙 )+ 𝑮 ( 𝒙 , 𝒖 ) ] 𝜕 𝒖 (cid:12)(cid:12) . The virtual control 𝝂 𝑐 = − 𝑲 𝒛 + (cid:164) 𝒙 ref2 − 𝑮 T 𝒛 is used for stabilizing the nominal dynamics. 𝑲 is a positive definite diagonal gain matrix.Consider the sliding surface 𝝈 = 𝒛 = , the robustify virtual control 𝝂 𝑠 is designed as: 𝝂 𝑠 = − 𝑲 𝑠 sig ( 𝝈 ) 𝜸 𝑠 = −[ 𝐾 𝑠, | 𝜎 | 𝛾 𝑠, sign ( 𝜎 ) , 𝐾 𝑠, | 𝜎 | 𝛾 𝑠, sign ( 𝜎 ) , 𝐾 𝑠, | 𝜎 | 𝛾 𝑠, sign ( 𝜎 )] T (46)where 𝐾 𝑠,𝑖 > , 𝛾 𝑠,𝑖 ∈ ( , ) . In the conventional sliding mode control, the discontinuous sign function has to besmoothened for chattering reduction. This process is not needed here, since without any approximation, | 𝜎 𝑖 | 𝛾 𝑖 sign ( 𝜎 𝑖 ) itself is a continuous function of 𝜎 𝑖 . Substituting Eq. (45) into Eq. (43), the closed-loop dynamics are (cid:164) 𝒙 = 𝝂 𝑐 + 𝜹 ( 𝒙 , Δ 𝑡 ) + ( 𝑮 ¯ 𝑮 − − 𝑰 )( 𝝂 𝑐 − (cid:164) 𝒙 , ) + 𝑮 ¯ 𝑮 − 𝝂 𝑠 + Δ 𝒅 (cid:44) 𝝂 𝑐 + 𝑮 ¯ 𝑮 − 𝝂 𝑠 + 𝜺 ibs (47)where 𝜹 ( 𝒙 , Δ 𝑡 ) equals the summation of the closed-loop values of 𝑹 (cid:48) (Eq. (44)) and 𝜕 [ 𝒇 ( 𝒙 )+ 𝑮 ( 𝒙 , 𝒖 ) ] 𝜕 𝒙 (cid:12)(cid:12) Δ 𝒙 (Eq. (43)). Proposition 2 If (cid:107) 𝑰 − 𝑮 ¯ 𝑮 − (cid:107) ≤ ¯ 𝑏 < for all 𝑡 , and if (cid:107) 𝜹 ( 𝒙 , Δ 𝑡 )(cid:107) ≤ ¯ 𝛿 , (cid:107) Δ 𝒅 (cid:107) ≤ Δ 𝑑 , then under sufficiently highsampling frequency, 𝜺 ibs given by Eq. (47) is bounded for all 𝑡 , and is ultimately bounded by Δ 𝜈 𝑐 ¯ 𝑏 + ¯ 𝛿 + Δ 𝑑 − ¯ 𝑏 , where Δ 𝜈 𝑐 isthe upper bounds of Δ 𝝂 𝑐 . This proposition can be proved following the proofs of Theorem 1 in [15]. It puts requirement on the controleffectiveness estimation, i.e., 𝑮 ¯ 𝑮 − has to be diagonally dominated. Proposition 3 If 𝜺 ibs (Eq. (47)) is bounded, then using the control incremental in Eq. (45) , both 𝒛 and 𝒛 are ultimatelybounded. Proof:
Consider a candidate Lyapunov function 𝑉 ( 𝒛 , 𝒛 ) = ( 𝒛 T 𝒛 + 𝒛 T 𝒛 ) , using Eqs. (42) and (45), the timederivative of 𝑉 are derived as (cid:164) 𝑉 = − 𝒛 T 𝑲 𝒛 − 𝒛 T 𝑲 𝒛 + 𝒛 T ( 𝜺 ibs − 𝑮 𝑛 ¯ 𝑮 − 𝑛 𝑲 𝑠 sig ( 𝝈 ) 𝜸 )≤ − 𝒛 T 𝑲 𝒛 − 𝒛 T 𝑲 𝒛 + ∑︁ 𝑖 = (| 𝜎 𝑖 || 𝜀 ibs ,𝑖 | + ¯ 𝑏𝐾 𝑠,𝑖 | 𝜎 𝑖 | 𝛾 𝑠,𝑖 + − 𝐾 𝑠,𝑖 | 𝜎 𝑖 | 𝛾 𝑠,𝑖 + )≤ − 𝒛 T 𝑲 𝒛 − 𝒛 T 𝑲 𝒛 − ∑︁ 𝑖 = 𝜌 𝑖 | 𝜎 𝑖 | , ∀| 𝜎 𝑖 | ≥ (cid:18) 𝜌 𝑖 + | 𝜀 ibs ,𝑖 |( − ¯ 𝑏 ) 𝐾 𝑠,𝑖 (cid:19) 𝛾𝑠,𝑖 , ∀ 𝜌 𝑖 > 𝜎 𝑖 equals (cid:16) 𝜌 𝑖 +| 𝜀 ibs ,𝑖 |( − ¯ 𝑏 ) 𝐾 𝑠,𝑖 (cid:17) 𝛾𝑠,𝑖 , whose size can be madearbitrarily small by increasing 𝐾 𝑠,𝑖 and reducing 𝛾 𝑠,𝑖 . Because 𝝈 = 𝒛 , 𝑲 and 𝑲 are positive definite, then both 𝒛 and 𝒛 are ultimately bounded. (cid:3) Remark 3
Incremental backstepping sliding mode control has lower model dependency than feedback linearization andbackstepping sliding mode control. The only model information needed for implementation is the control effectiveness¯ 𝑮 . Even so, by virtue of its sensor-based nature, incremental backstepping sliding mode control actually has enhancedrobustness against model uncertainties, sudden actuator faults and structural damage [15, 23]. Under the sameperturbation circumstance, there exists a sampling frequency such that the bound of 𝜺 ibs is smaller than that of theresidual error under backstepping control. This property enables incremental backstepping sliding mode control topassively resist a wider range of perturbations with much lower sliding mode control gains [15, 23].12 . Optimal Multi-objective Wing Control In subsection III.C, the reference values for 𝛿 𝑒 , 𝛿 𝑟 , 𝑀 𝜙, diff have been designed for aircraft attitude control purposes.The references for 𝛿 𝑒 and 𝛿 𝑟 can be directly given to the actuators of elevator and rudder, whereas 𝑀 ref 𝜙, diff still needsto be achieved by aircraft trailing-edge control surfaces. The aircraft model considered in this paper has fourteendistributed flaps, which makes the achievement of 𝑀 ref 𝜙, diff an over-actuated problem. As a consequence, after 𝑀 ref 𝜙, diff is achieved for the purpose of bank angle tracking, the wing system still have control redundancy to achieve otherobjectives such as maneuver load alleviation, gust load alleviation, flutter suppression, etc. In view of this, an optimalmulti-objective wing control law will be designed in this subsection.Wing-root shear force 𝐹 𝑤 and wing-root bending moment 𝑀 𝜙 are two important wing load indicators. Theirs valueswould deviate from the trimmed condition during aircraft maneuvers, under the excitations of atmospheric disturbances,as well as under flap inputs. In a trimmed flight, when an aircraft encounters disturbances, a reasonable control objectiveis to maintain 𝐹 𝑤 and 𝑀 𝜙 at their trimmed values. However, this objective is not applicable during aircraft maneuversfor two reasons: 1) in view of Eq. (26), the wing-root shear force 𝐹 𝑤 , which is strongly coupled with the total lift ona half-wing, is the main medium for changing flight path angle 𝛾 ; 2) as discussed in subsection III.C, the differencebetween the left and right wing-root bending moments 𝑀 𝜙, diff is the main medium for aircraft roll control.Based on these discussions, it is essential to identify the loads that are necessary to achieve commanded maneuvers,while reducing any excessive load induced by either maneuvers or external disturbances. A novel reference generator isdesigned to fulfill this goal. The designed references for the left and right wing-root shear forces are 𝐹 ref 𝑤,𝑟 = 𝐹 ref 𝑤,𝑙 = 𝐹 trim 𝑤 + 𝐺 𝐹 𝑤 ( 𝑠 )( 𝛼 ref − 𝛼 trim ) (49)In Eq. (49), 𝛼 ref is the angle of attack reference designed in subsection III.B. 𝐺 𝐹 𝑤 ( 𝑠 ) represents an estimatedmapping from 𝛼 input to the open-loop shear force response of a half wing. In this paper, 𝐺 𝐹 𝑤 ( 𝑠 ) is designed as alow-pass filer, with 𝑠 represents the Laplace variable. The zero-frequency gain of 𝐺 𝐹 𝑤 ( 𝑠 ) equals 𝑞 ∞ 𝑆 𝑤 𝐶 𝐿𝛼 𝑤 /
2, where 𝑆 𝑤 is the wing area and 𝐶 𝐿𝛼 𝑤 is the lift-slope of the wing. The dynamics of 𝐺 𝐹 𝑤 ( 𝑠 ) are caused by the time-dependentcirculatory effects. Therefore, the time-constant of 𝐺 𝐹 𝑤 ( 𝑠 ) can be either calculated using Eq. (12), or be identifiedonline (the output 𝐹 𝑤 can be measured by strain gauges). If 𝛼 ref = 𝛼 trim , the shear force references are identical to theirtrimmed values, and any variations caused by disturbances or maneuvers (such as the lift drop during a sharp-roll) willbe automatically counteracted by flaps.Furthermore, the reference generation algorithm for the left and right wing-root bending moments are given inAlgorithm 1. Algorithm 1:
Reference generation for left and right wing-root bending moments 𝑀 ref 𝜙,𝑟 ← 𝑀 trim 𝜙,𝑟 − 𝑀 ref 𝜙, diff / 𝑀 ref 𝜙,𝑙 ← 𝑀 trim 𝜙,𝑙 + 𝑀 ref 𝜙, diff / if 𝑀 ref 𝜙,𝑟 > 𝑀 ref 𝜙 then 𝑀 ref 𝜙,𝑟 ← 𝑀 ref 𝜙 ; 𝑀 ref 𝜙,𝑙 ← 𝑀 ref 𝜙, diff + 𝑀 ref 𝜙 elseif 𝑀 ref 𝜙,𝑙 > 𝑀 ref 𝜙 then 𝑀 ref 𝜙,𝑟 ← 𝑀 ref 𝜙 − 𝑀 ref 𝜙, diff ; 𝑀 ref 𝜙,𝑙 ← 𝑀 ref 𝜙 end ifend if In Algorithm 1, 𝑀 ref 𝜙, diff has been designed in subsection III.C, which is the core for achieving the designed rollangle tracking performance. 𝑀 ref 𝜙 in the upper limit on the reference for wing-root bending moment 𝑀 𝜙 , which can bedesigned much smaller than the real physical limit max ( 𝑀 𝜙 ) , i.e., 𝑀 trim 𝜙,𝑙 = 𝑀 trim 𝜙,𝑟 < 𝑀 ref 𝜙 < max ( 𝑀 𝜙 ) . This design canensure a safety margin from max ( 𝑀 𝜙 ) , and is also able to constrain the variations of 𝑀 𝜙 for extending structural fatiguelife.As formulated in Algorithm 1, to achieve a moderate roll maneuver, the left and right wing take the same responsibility.However, once one of the wing-root bending moments reaches 𝑀 ref 𝜙 under a sharp-roll command, the exceeded commandwill be automatically allocated to the other half-wing. In such a way, not only the wing loads are limited, 𝑀 ref 𝜙, diff isachieved for roll-command tracking as well. Remark 4
The wing-root shear force and bending moment reference generation algorithms use very little model13nformation, and are easy to be implemented in real-time. More importantly, by exploiting the control redundancy, theload alleviation is achieved without degrading the command tracking performance.After the references 𝐹 ref 𝑤,𝑟 , 𝐹 ref 𝑤,𝑙 , 𝑀 ref 𝜙,𝑟 , 𝑀 ref 𝜙,𝑙 are derived, a controller should be designed to track these references.Taking the right wing as an example, only two independent control variables are needed for tracking 𝐹 ref 𝑤,𝑟 and 𝑀 ref 𝜙,𝑟 .Since we have seven flaps on each wing, the remaining control space can be used for flutter suppression, aeroelasticdamping enhancement, and control energy reduction. Choose the controlled output as 𝒚 = [ 𝐹 𝑤,𝑟 , 𝑀 𝜙,𝑟 ] T , then theright-wing dynamics given by Eqs. (19) and (20) are represented as (cid:164) 𝒙 𝑒 = 𝑨𝒙 𝑒 + 𝑩 𝑢 𝒖 + 𝑩 𝑔 𝜶 𝑔 + 𝑩 grav + 𝑩 𝑟 𝜶 𝑞𝑠,𝑟 + 𝑩 acc 𝒚 = 𝑪𝒙 𝑒 + 𝑫 𝑢 𝒖 + 𝑫 grav + 𝑫 𝑟 𝜶 𝑞𝑠,𝑟 + 𝑫 acc (50)Further design an augmented state vector 𝑿 = [ 𝒙 T 𝑒 , ∫ 𝒆 T 𝑦 ] T , 𝒆 𝑦 = 𝒚 − 𝒚 ref for load commands tracking, then theaugmented system dynamics are (cid:34) (cid:164) 𝒙 𝑒 𝒆 𝑦 (cid:35) = (cid:34) 𝑨 𝑪 (cid:35) (cid:34) 𝒙 𝑒 ∫ 𝒆 𝑦 (cid:35) + (cid:34) 𝑩 𝑢 𝑫 𝑢 (cid:35) 𝒖 + (cid:34) 𝑩 𝑟 𝑫 𝑟 (cid:35) 𝜶 𝑞𝑠,𝑟 + (cid:34) 𝑩 g (cid:35) 𝜶 𝑔 + (cid:34) 𝑩 grav 𝑫 grav (cid:35) + (cid:34) 𝑩 acc 𝑫 acc (cid:35) + (cid:34) − 𝒚 ref (cid:35) (51)To evaluate the robustness of the controller, the disturbance input 𝜶 𝑔 and the inertial input vector [ 𝑩 T acc , 𝑫 T acc ] T areeliminated from the model used for control design. The control performance can be further enhanced if extra knowledgeabout 𝜶 𝑔 is available (from LIDAR or a gust estimator). The local quasi-steady angle of attack induced by rigid-bodytranslational and rotational motion 𝜶 𝑞𝑠,𝑟 (Eq. (13)) is treated as a known input in wing control design. In addition, thegravitational forces are also treated as known inputs in wing control design. Formulate a cost function as 𝐽 = lim 𝑡 𝑓 →∞ ∫ 𝑡 𝑓 [ 𝑿 T 𝑸 𝑿 + 𝒖 T 𝑹𝒖 ] d 𝑡 (52)where 𝑸 and 𝑹 are positive definite diagonal matrices. Denote the augmented system dynamic matrix as 𝑨 aug , anddenote 𝑩 aug = [ 𝑩 T 𝑢 , 𝑫 T 𝑢 ] T , then the infinite time-horizon optimal control is designed as 𝒖 = 𝑲 𝑋 𝑿 + 𝑲 𝑟 (cid:32) (cid:34) 𝑩 𝑟 𝑫 𝑟 (cid:35) 𝜶 𝑞𝑠,𝑟 + (cid:34) 𝑩 grav 𝑫 grav (cid:35) + (cid:34) − 𝒚 ref (cid:35) (cid:33) (53)where 𝑲 𝑋 = − 𝑹 − 𝑩 T aug 𝑾 , 𝑲 𝑟 = − 𝑹 − 𝑩 T aug ( 𝑾 𝑩 aug 𝑹 − 𝑩 T aug − 𝑨 T aug ) − 𝑾 (54) 𝑾 in Eq. (54) is the solution of the associated Riccati equation. The designed control input can achieve the followinggoals: 1) stabilize the wing aeroelastic system (flutter suppression); 2) enhance structural damping (by assigning theclosed-loop eigenvalues); 3) track the designed load commands; 4) minimize control energy. The states needed forfeedback control can be observed by a Kalman filter [13]. The proposed control architecture is illustrated in Fig. 2. IV. Simulation Results and Discussions
In this section, the effectiveness of the proposed control architecture is evaluated. A 3D view of the planform isshown in Fig. 3. Top, back and side views of the planform are shown in Fig. 4. The wing was designed using XFLR5.The steady-state aerodynamic lift coefficients were extracted using the vortex lattice method with a total of 1967 meshelements. The aspect ratio of the wing equals 26.67. The total mass of the aircraft equals 227 kg. The aircraft momentof inertia in trim condition is 𝐼 𝑥𝑥 = . · m , 𝐼 𝑦𝑦 = . · m , 𝐼 𝑧𝑧 = . · m , and 𝐼 𝑥𝑧 = − . · m .During maneuvers, these inertia will change with structural deformations. The actuator dynamics of elevator and rudderare modeled as first-order systems, with time constant equals 0.02 s. The deflection limits of them both equal ±
20 deg.As discussed in subsection II.B, the fourteen wing flaps are connected with the main wing structure via rotationalsprings. The inertia couplings between the flap and structure are considered in the structure mass matrix. The deflectionlimits of the flaps are set as ±
30 deg. In addition, the throttle dynamics are model as a first-order system with timeconstant equals 0.2 s, which is slower than that of control surfaces. To capture the high-frequency dynamics of thewing aeroelastic dynamics, the simulation sampling frequency is chosen as 20,000 Hz. On the contrary, the samplingfrequency used by controllers are much lower. The sampling frequency used by the attitude and wing control loopsis 100 Hz, while 50 Hz is used by the position and flight path control loops. 100 Hz is a reasonable choice for theincremental control of aircraft attitude dynamics, which has been verified by flight tests on a CS-25 aircraft [8].14 ositionControl: NonlinearDynamicInversion
AttitudeControl: IncrementalBacksteppingSliding ModeControlFlight PathControl: IncrementalSliding ModeControl LoadReferenceGenerator Optimal Multi-objectiveWing Control VelocityControl
Fig. 2 Control architecture.Fig. 3 3D view of the flexible aircraft model.
XZ Y1.12m 8m (a) Top view ( 𝑥, 𝑦 ) ZYX 1.20m (b) Side view ( 𝑥, 𝑧 ) XZ Y (c) Back view ( 𝑦, 𝑧 ) Fig. 4 Top, back and side views of the flexible aircraft.A. Trim and Model Analysis
A free-flying rigid aircraft model is set up for comparative studies. This rigid aircraft shares the same geometricand static aerodynamic properties with the free-flying flexible aircraft model. Nevertheless, its structural stiffness areassumed to be infinitely high, and its aerodynamic model relies on quasi-steady strip theory. Both the rigid and flexibleaircraft are trimmed in a steady-level flight condition at 𝐻 = 𝑉 =
35 m/s (a typical operational point for glideraircraft). The trim solutions are shown in Table 1.To analyze their characteristics, the rigid and flexible aircraft are linearized around their equilibrium points, withthe resulting eigenvalues illustrated in Fig. 5. The poles of the clamped aeroelastic wing are also shown. It can beseen from Fig. 5, most of the flexible aircraft poles appear in the high-frequency range, and are in good agreementwith those of the clamped wing. Discrepancies caused by the interactions between structural and rigid-body motions15 able 1 Trim solutions for the rigid and flexible aircraft. 𝛼 ∗ [ ◦ ] 𝛿 𝑒 ∗ [ ◦ ] 𝐹 𝐸 ∗ [N]Rigid aircraft 3.467 − . − . Fig. 5 Eigenvalues of the aereoelastic clamped wing, the rigid aircraft, and the free-flying flexible aircraft.
The rigid-body and structural couplings are clearly exposed in the above analyses, which underline the necessity of amulti-objective integrated controller. From Fig. 5, the phugoid mode should be stabilized. Moreover, the damping ratioof the short-period, Dutch roll and the new aeroelastic-lateral coupled modes need to be increased. Furthermore, thecontroller should simultaneously fulfill the trajectory tracking commands and the load alleviation requirements. Theeffectiveness of the proposed control architecture will be shown in the following subsections.
B. Maneuver Load Alleviation
Four maneuver circumstances will be considered in this subsection: sudden pull-up, sharp roll, flight path control,and 3D aircraft position control. The controller aims at minimizing the tracking errors while alleviating the loads causedby maneuvers.
1. Load Alleviation in a Pull-up Maneuver
In this simulation case, the aircraft is commanded to track an 𝛼 profile, the kinematic bank angle and the side-slipangle are commanded to remain zero. As illustrated in the first subplot of Fig. 6, two smoothly combined sigmoidfunctions 𝑓 = /( + 𝑒 − ( 𝑡 − ) ) and 𝑓 = − /( + 𝑒 − ( 𝑡 − ) ) are used to compose the command. In practice, a commonchoice of command is a discontinuous step function filtered by a low-pass filter. However, at the step points, the filteredsignal is still non-differentiable. These non-differentiable points would cause spikes in control inputs. In view of this,the sigmoid function, which is differentiable up to any order, is chosen as a smooth realization of the step function in16his paper. Since this maneuver is symmetrical, only the right wing responses are shown. [ ° ] F w , r [ N ] t [s] M , r [ N m ] Reference with MLA without MLA Limit(Ref(M )) t [s] M , r [ N m ] Fig. 6 𝛼 tracking performance and the responses of wing-root shear force and bending moment. To demonstrate the effectiveness of the proposed control law, another controller without the MLA function is alsodesigned. These two controllers give the same commands to elevator and rudder. However, for the controller withoutMLA, the flap hinge moment commands are set to zero during symmetrical maneuvers. As shown in Fig. 6, both thecontrollers with and without MLA can track the given 𝛼 command, with max (| 𝛼 − 𝛼 ref |) < . ◦ . The wing-root shearforce responses of the two controllers are also similar. Nevertheless, the wing-root bending moment is effectivelyreduced by the MLA function. Recall Algorithm 1, 𝑀 ref 𝜙,𝑟 = 𝑀 trim 𝜙,𝑟 in symmetrical maneuvers ( 𝑀 ref 𝜙, diff = (| 𝑀 𝜙,𝑟 − 𝑀 trim 𝜙,𝑟 |) from 4183 N · m to 1.344 N · m (by99.97%), and reduces rms (| 𝑀 𝜙,𝑟 − 𝑀 trim 𝜙,𝑟 |) from 2387 N · m to 0.5827 N · m (by 99.98%). It can also be seen from Fig. 7,that the MLA function make the inner-board flaps deflect downwards, and the out-board flaps deflect upwards, whichshifts the pressure center towards the wing root. More importantly, by exploiting the control redundancy, the 𝛼 trackingperformance is not influenced by the MLA function. - w s , r [ m ] - w s , r [ m ] t [s] -10010 s , r [ ° ] t [s] -0.5-0.4-0.3 s , r [ ° ] Fig. 7 Responses of the wing vertical displacements and flap angles with (left) and without MLA (right).
2. Load Alleviation in a Sharp-roll Maneuver
In this test case, the aircraft is commanded to roll from 0 deg to 40 deg within two seconds. The sigmoid function 𝑓 = /( + 𝑒 − ( 𝑡 − ) ) is adopted as a smooth realization of the step function. Another controller without the MLA functionis designed for comparisons. This controller uses the same outer-loop control algorithms with the nominal controller, but17here are two differences in the optimal wing control loop: 1) the shear force command tracking is eliminated; 2) the rightand left wing always share the same responsibility in roll tracking, i.e., 𝑀 ref 𝜙,𝑟 ≡ 𝑀 trim 𝜙,𝑟 − 𝑀 ref 𝜙, diff / , 𝑀 ref 𝜙,𝑙 ≡ 𝑀 trim 𝜙,𝑙 + 𝑀 ref 𝜙, diff / 𝑀 ref 𝜙 is set at 5000 N · m, which is much lower than the real physical limit on 𝑀 𝜙 . In thesteady-level flight, the trimmed value for root bending moment equals 4914 N · m, thus | 𝑀 ref 𝜙 − 𝑀 trim 𝜙 | is less than 2% of 𝑀 trim 𝜙 . Under this strict constrain, and tracking performance and load responses are shown in Fig. 8. [ ° ] F w , r & F w , l [ N ] t [s] M , r [ N m ] t [s] M , l [ N m ] Reference with MLA without MLA Limit(Ref(M ))
Fig. 8 Responses of the kinematic bank angle, wing-root shear force and bending moment.
In view of Fig. 8, the kinematic bank angle tracking performance with and without MLA are almost identical. Bothof them can track the given command very well, with max (| 𝜇 − 𝜇 ref |) < . ◦ . The responses of the left and rightwing-root shear forces are shown in the second subplot, from which it can be seen that without MLA, 𝐹 𝑤,𝑟 and 𝐹 𝑤,𝑙 have large variations because of the 𝜶 𝑞𝑠,𝑟 induced by roll maneuvers. By contrast, the MLA function can compensatefor the shear force variations, since 𝐹 ref 𝑤,𝑟 = 𝐹 ref 𝑤,𝑙 = 𝐹 trim 𝑤 when 𝛼 ref = ( 𝐹 𝑤,𝑟 − 𝐹 ref 𝑤,𝑟 ) is reduced from 63.83 N to 6.235 N (by 90.23%). Moreover, as illustrated in the third and fourthsubplots, the MLA algorithm successfully complies with the strict limit 𝑀 ref 𝜙 = · m. Once one of 𝑀 𝜙,𝑟 , 𝑀 𝜙,𝑙 reaches 𝑀 ref 𝜙 , the exceeded command is automatically allocated to the other wing without influencing the roll trackingperformance. On the contrary, if the MLA algorithm is not used during this maneuver, max ( 𝑀 𝜙,𝑟 ) = · m, andmax ( 𝑀 𝜙,𝑙 ) = · m.The load reduction without influencing the command tracking performance is achieved by fully exploiting the inputredundancy in control design (subsection III.D). The responses of the wing bending displacements and flap deflectionsduring the sharp-roll maneuver are shown in Figs. 9 and 10. It can be observed that by actuating the distributedtrailing-edge control surfaces, the proposed controller can not only reduce loads, but constrain the maximum wingdeflections as well. In addition, the flap deflection angles remain in the ±
30 deg limits.
3. Flight Path Control With Maneuver Load Alleviation
In this simulation case, the aircraft is commanded to follow a spiral trajectory, which is shown in Figs. 11 and 12.This trajectory can be realized by giving a step command to flight path angle 𝛾 , and a ramp command to the kinematicazimuth angle 𝜒 . Smooth realizations of the step commands are given to the aircraft, as shown in Fig. 13. The upperlimit on the reference for wing-root bending moment is still set as 𝑀 ref 𝜙 = · m.As can be viewed from Figs. 11 and 12, the proposed controller can make the flexible aircraft track the givencommand, with max (| 𝑋 − 𝑋 ref |) < .
10 m, max (| 𝑌 − 𝑌 ref |) < .
11 m and max (| 𝐻 − 𝐻 ref |) < .
42 m. Fig. 13 shows18 - w s , r [ m ] - w s , r [ m ] t [s] -20-100 s , r [ ° ] t [s] -20-100 s , r [ ° ] Fig. 9 Responses of the right wing displacements and flap deflections with (left) and without MLA (right). - w s , l [ m ] - w s , l [ m ] t [s] s , l [ ° ] t [s] s , l [ ° ] Fig. 10 Responses of the left wing displacements and flap deflections with (left) and without MLA (right).Fig. 11 3D spiral trajectory tracking performance. the rigid-body state responses during this maneuver. The flight path angle increases from 0 deg to 8 deg within 3.5 s,and the maximum tracking error equals 0.19 ◦ . A ramp command with slope 5 . ◦ /s is given to the kinematic azimuthangle 𝜒 , and the aircraft can track it with maximum tracking error equals 0.026 ◦ . The aircraft velocity is maintained atthe trimmed value using a linear throttle control ( 𝛿 𝑇 ). Since the aircraft enters a stable spiral mode after 𝑡 =
25 s, theresponses of inner-loop states are shown in 𝑡 ∈ [ , ] s. The attitude commands given by the flight path control loopare also tracked, with max (| 𝑒 𝜇 |) = . ◦ , max (| 𝑒 𝛼 |) = . ◦ , and max (| 𝑒 𝛽 |) = . ◦ .19
200 400 600
X [m] Y [ m ] X [m] H [ m ] Y [m] H [ m ] Reference Response
Fig. 12 Three views of the 3D spiral trajectory tracking performance. [ ° ] [ ° ] t [s] V [ m / s ] [ ° ] [ ° ] t [s] -0.01-0.00500.0050.01 [ ° ] p [ ° / s ] q [ ° / s ] t [s] r [ ° / s ] Reference Response
Fig. 13 Responses of flight path angles, attitude angles and angular rates in an aircraft spiral maneuver.
The load responses of the aircraft during the spiral maneuver are shown in Fig. 14. When 𝑡 ∈ [ , ] s, onlysymmetrical states are excited, thus 𝑀 𝜙,𝑟 and 𝑀 𝜙,𝑙 are maintained at their trimmed values, while 𝐹 𝑤,𝑟 , 𝐹 𝑤,𝑙 are drivento realize the 𝛾 tracking command (Eq. (26)). When 𝑡 ∈ [ , ] 𝑠 , the aircraft starts entering the spiral. Similar tothe responses in the sharp-roll maneuver illustrated in Fig. 8, in case one of 𝑀 𝜙,𝑟 , 𝑀 𝜙,𝑙 reaches 𝑀 ref 𝜙 , the exceededcommand is automatically allocated to the other wing.In this spiral maneuver, the responses of the wing bending displacements and flap deflections are shown in Fig. 15.It can be seen when 𝑡 ∈ [ , ] 𝑠 , the flap deflections are the same, while when 𝑡 ∈ [ , ] 𝑠 , the flap deflections arere-allocated for roll tracking and MLA purposes.
4. Position Control With Maneuver Load Alleviation
A more complex 3D maneuver with strongly coupled longitudinal and lateral motions will be tested in this subsection.The trajectory references and responses are illustrated in Figs. 16 and 17. The proposed controller can minimize thetrajectory errors, with max (| 𝑋 − 𝑋 ref |) < . (| 𝑌 − 𝑌 ref |) < .
019 m and max (| 𝐻 − 𝐻 ref |) < .
023 m duringthis maneuver.The responses of rigid-body states are shown in Fig. 18. Although the lateral and longitudinal control commandsare coupled in this case, the proposed nonlinear controller can decouple the resulting dynamics through the control20 M [ N m ] Limit(Ref(M )) Ref(M ,r ) M ,r Ref(M ,l ) M ,l F w [ N ] Ref(F w,r ) F w,r
Ref(F w,l ) F w,l t [s] n z [ g ] Fig. 14 Responses of the wing-root bending moments, shear forces, and load factor in a spiral maneuver. - w s , r [ m ] - w s , l [ m ] t [s] -10010 s , r [ ° ] t [s] -10010 s , l [ ° ] Fig. 15 Responses of the wing bending displacements and flap deflections in a spiral maneuver.Fig. 16 The flexible aircraft 3D trajectory tracking responses. effectiveness inversion. As a result, the pitch, roll, and yaw channels are governed by their own designed tracking21
X [m] -80-60-40-200 Y [ m ] X [m] H [ m ] -60 -40 -20 0 Y [m] H [ m ] Reference Response
Fig. 17 Three views of the 3D trajectory tracking performance. error dynamics, which are then realized by the virtual controls (Eq. (45)). In Fig. 18, max (| 𝑒 𝛾 |) = . ◦ , andmax (| 𝑒 𝜒 |) = . ◦ . [ ° ] [ ° ] t [s] V [ m / s ] [ ° ] [ ° ] t [s] -0.01-0.00500.0050.01 [ ° ] p [ ° / s ] q [ ° / s ] t [s] -505 r [ ° / s ] Reference Response
Fig. 18 Responses of flight path angles, attitude angles, and angular rates in an aircraft 3D maneuver.
The load responses are shown in Fig. 19. The differences between the right and left wing root bending moments arenecessary for achieving lateral maneuvers (Sec. III.C). By virtue of the MLA algorithm, the reference limit 𝑀 ref 𝜙 isnot exceeded. The shear force commands are needed for achieving flight path angle command, and they have similarfrequency component with the load factor 𝑛 𝑧 . Figure 20 shows the wing displacements and flap deflections. During thismaneuver, the maximum wing-tip displacement equals 0.74 m, which only deviates 0.06 m from its nominal value. Theflap deflections are smooth, and are within the ±
30 deg limit.
C. Gust Load Alleviation
In this section, the gust load alleviation performance of the designed controller will be evaluated. The aircraft iscommanded to maintain at a cruise condition: 𝑉 = , 𝐻 = , 𝛾 = 𝜒 =
0. The von Kármán turbulence modelis adopted in this paper. Using two-dimensional Fourier transforms, the von Kármán spectrum can be realized in atwo-dimensional spatial domain [13]. Realizing a turbulence spectrum in spatial domain has several benefits: first of all,it allows the consideration of gust penetration effect [27]. Moreover, the span-wise gust velocity variations are modeled,22
10 20 30 40 50 60480049005000 M [ N m ] Limit(Ref(M )) Ref(M ,r ) M ,r Ref(M ,l ) M ,l F w [ N ] Ref(F w,r ) F w,r
Ref(F w,l ) F w,l t [s] n z [ g ] Fig. 19 Responses of the wing-root bending moments, shear forces, and load factor in a 3D maneuver. - w s , r [ m ] - w s , l [ m ] t [s] -505 s , r [ ° ] t [s] -505 s , l [ ° ] Fig. 20 Responses of the wing bending displacements and flap deflections in a 3D maneuver. which are important for high aspect-ratio aircraft. Furthermore, in a spatial turbulence field, the aircraft maneuvers arenot constrained in the symmetric plane; multi-dimensional maneuvers are allowed. In addition, given the aircraft groundvelocity, the spatial domain turbulence field can be easily embedded in time-domain simulations. Fig. 21 shows a severevon Kármán turbulence field with the turbulence scale length 𝐿 𝑔 =
762 m, and the turbulence intensity 𝜎 = ( 𝑛 𝑧 − 𝑛 ∗ 𝑧 ) is reduced from 0.439 g to 0.0159 g (by 96.38%), rms ( 𝑀 𝜙,𝑟 − 𝑀 ∗ 𝜙,𝑟 ) is alleviated from 2057 N · m to 16.45 N · m (by99.20%), and rms ( 𝐹 𝑤,𝑟 − 𝐹 ∗ 𝑤,𝑟 ) is diminished from 511.3 N to 38.68 N (by 92.43%). The superscript (·) ∗ denotes thetrimmed value. Moreover, the responses of angular rates are shown in Fig. 23. In the “open-loop” case, the roll rate andpitch rate are non-zero because the local angles of attack induced by the 2D spatial turbulence field are different on eachstrip. For the high-aspect ratio aircraft model used in this paper, the turbulence field has strongest influence on the rollchannel, with rms ( 𝑝 ) = . ◦ /s in the “open-loop” case, which is reduced to 0 . ◦ /s using the proposed controller23 ig. 21 A 2D severe von Kármán turbulence field ( 𝐿 𝑔 = m, 𝜎 = m/s). (reduced by 99.65%). Besides, this control architecture diminishes rms ( 𝑞 ) = from 5.31 ◦ /s to 0.134 ◦ /s (by 97.48%), andreduces rms ( 𝑟 ) = from 5.99 ◦ /s to 0.505 ◦ /s (by 91.57%). As discussed in [28], the load factor is tightly linked withpassenger ride comfort, thus by alleviating 𝑛 𝑧 , the ride comfort can also be improved. n z [ g ] M , r [ N m ] t [s] F w , r [ N ] Open-loop Closed-loop
Fig. 22 Responses of the load factor, root bending moment, and shear force in a von Kármán turbulence field.
The wing displacements in the “open-loop” and closed-loop cases are illustrated in Fig. 24. It can be seen thatthe structural motions are damped out by the proposed controller. The maximum deviation of the wing tip is reducedfrom 1.5 m to 0.81 m (by 46.00%). The flap deflections are shown in Fig. 25. As mentioned in the preceding texts,the “open-loop” case gives zero hinge moments, but the flaps can still move due to their inertia and stiffness couplingswith the main wing structure. When the control loop is closed, the flaps deflect in oppose to the turbulence profile,and modify the lift distributions at the same time. The deflection angles are within the limit of ±
30 deg in this severeturbulence field.
V. Conclusions
In this paper, a nonlinear control architecture for flexible aircraft trajectory tracking and load alleviation is proposed.To begin with, this paper derives the kinematics and dynamics of the free-flying flexible aircraft, which capture both the24 p [ ° / s ] q [ ° / s ] t [s] -10010 r [ ° / s ] Open-loop Closed-loop
Fig. 23 Angular velocity responses in a 2D von Kármán turbulence field. - w s , r [ m ] t [s] - w s , r [ m ] Fig. 24 Displacements of the seven structural nodes on the right wing (top: open-loop, bottom: closed-loop). s , r [ ° ] t [s] -20020 s , r [ ° ] Fig. 25 Deflections of the seven flaps on the right wing (top: open-loop, bottom: closed-loop). inertial and aerodynamic couplings between the rigid-body and structural degrees of freedom. The dynamic model isderived in a modular approach, which provides a convenient way to make an existing clamped-wing aeroservoelasticmodel free-flying.Based on the flexible aircraft model, a four-loop cascaded control architecture is proposed. The position control loopis designed using nonlinear dynamic inversion, which provides references to the flight path control loop. Because model25ncertainties and disturbances present, the flight path control loop adopts the incremental sliding mode control law,whose reduced model dependency and enhanced robustness are demonstrated theoretically. Moreover, the incrementalbackstepping sliding mode control is used in the attitude control loop. Based on Lyapunov methods, this loop isproven to be stable under disturbance perturbations. Furthermore, a novel load reference generator is proposed, whichdistinguishes the loads that are necessary to perform maneuvers from the excessive loads. The load references are thentracked by the inner-loop multi-objective wing control. At the same time, the excessive loads (no matter induced bymaneuvers or gusts) are naturalized by control surface deflections. Since the proposed control architecture exploits thecontrol redundancy, the load alleviation is achieved without affecting the aircraft command tracking performance.The effectiveness of the proposed control architecture is validated by numerical simulations. It is demonstrated thatthe loads during sudden pull-up and sharp roll maneuvers can be alleviated without influencing rigid-body commandtracking performance. Moreover, a simulation case that the flexible aircraft flies through a severe 2D spatial von Kármánturbulence field verifies the robustness of the controller to model uncertainties and external disturbances. Finally, theflexible aircraft is commanded to follow a 3D trajectory in a spatial turbulence field. Simulation results show thatboth the maneuver and gust loads are effectively alleviated without affecting the trajectory tracking performance. Theextension of the method to highly flexible aircraft with nonlinear beam structures will be explored in future work.
Acknowledgement
The first author would like to thank Lan Yang, Paul Lancelot, and Yang Meng for their help on the dynamic model.
References [1] Dillsaver, M., Cesnik, C. E., and Kolmanovsky, I., “Trajectory Control of Very Flexible Aircraft with Gust Disturbance,”
AIAAAtmospheric Flight Mechanics (AFM) Conference , American Institute of Aeronautics and Astronautics, Boston, MA, 2013.https://doi.org/10.2514/6.2013-4745, URL http://arc.aiaa.org/doi/10.2514/6.2013-4745.[2] Drew, M. C., Hashemi, K. E., Cramer, N. B., and Nguyen, N. T., “Multi-objective Optimal Control of the 6-DoF AeroservoelasticCommon Research Model with Aspect Ratio 13.5 Wing,”
AIAA Scitech 2019 Forum , American Institute of Aeronautics andAstronautics, San Diego, California, 2019. https://doi.org/10.2514/6.2019-0220, URL https://arc.aiaa.org/doi/10.2514/6.2019-0220.[3] Gaggero, M., and Caviglione, L., “MODEL PREDICTIVE CONTROL FOR MANEUVER LOAD ALLEVIATION INFLEXIBLE AIRLINERS,” Vol. 16, No. 1, 2019, pp. 420–432. URL https://ieeexplore.ieee.org/document/8355726/.[4] Nguyen, N. T., Ting, E., Chaparro, D., Drew, M. C., and Swei, S. S.-M., “Multi-Objective Flight Control for DragMinimization and Load Alleviation of High-Aspect Ratio Flexible Wing Aircraft,” , American Institute of Aeronautics and Astronautics, Grapevine, Texas, 2017.https://doi.org/10.2514/6.2017-1589.[5] Hansen, J. H., Duan, M., Kolmanovsky, I., and Cesnik, C. E., “Control Allocation for Maneuver and Gust Load Alleviation ofFlexible Aircraft,”
AIAA Scitech 2020 Forum , American Institute of Aeronautics and Astronautics, Orlando, Florida, 2020, pp.1–16. https://doi.org/10.2514/6.2020-1186.c1, URL https://arc.aiaa.org/doi/10.2514/6.2020-1186.c1.[6] Duan, M., Hansen, J. H., Kolmanovsky, I. V., and Cesnik, C. E., “Maneuver load alleviation of flexible aircraft through controlallocation: A case study using X-Hale,”
International Forum on Aeroelasticity and Structural Dynamics 2019, IFASD 2019 ,Savannah, Georgia, 2019, pp. 1–16.[7] Sieberling, S., Chu, Q. P., and Mulder, J. A., “Robust Flight Control Using Incremental Nonlinear Dynamic Inversion andAngular Acceleration Prediction,”
Journal of Guidance, Control, and Dynamics , Vol. 33, No. 6, 2010, pp. 1732–1742.https://doi.org/10.2514/1.49978.[8] Grondman, F., Looye, G., Kuchar, R. O., Chu, Q. P., and van Kampen, E., “Design and Flight Testing of Incremental NonlinearDynamic Inversion-based Control Laws for a Passenger Aircraft,” ,American Institute of Aeronautics and Astronautics, Kissimmee, Florida, 2018. https://doi.org/10.2514/6.2018-0385.[9] Sun, S., Wang, X., Chu, Q., and de Visser, C., “Incremental Nonlinear Fault-Tolerant Control of a Quadrotor With CompleteLoss of Two Opposing Rotors,”
IEEE Transactions on Robotics , 2020, pp. 1–15. https://doi.org/10.1109/TRO.2020.3010626,URL https://ieeexplore.ieee.org/document/9160894/.
10] Mooij, E., “Robust Control of a Conventional Aeroelastic Launch Vehicle,”
AIAA Scitech 2020 Forum , American Institute ofAeronautics and Astronautics, Orlando, Florida, 2020. https://doi.org/10.2514/6.2020-1103, URL https://arc.aiaa.org/doi/10.2514/6.2020-1103.[11] Wang, X., van Kampen, E., Chu, Q. P., and Lu, P., “Stability Analysis for Incremental Nonlinear Dynamic Inversion Control,”
Journal of Guidance, Control, and Dynamics , 2019. https://doi.org/10.2514/1.G003791.[12] Shearer, C. M., and S. Cesnik, C. E., “Trajectory Control for Very Flexible Aircraft,”
Journal of Guidance, Control, andDynamics , Vol. 31, No. 2, 2008, pp. 340–357. https://doi.org/10.2514/1.29335, URL http://arc.aiaa.org/doi/abs/10.2514/1.29335.[13] Wang, X., van Kampen, E., Chu, Q. P., and De Breuker, R., “Flexible Aircraft Gust Load Alleviation with Incremental NonlinearDynamic Inversion,”
Journal of Guidance, Control, and Dynamics , 2019, pp. 1–18. https://doi.org/10.2514/1.G003980, URLhttps://arc.aiaa.org/doi/10.2514/1.G003980.[14] Wang, X., van Kampen, E., Chu, Q. P., and Lu, P., “Incremental Sliding-Mode Fault-Tolerant Flight Control,”
Journal ofGuidance, Control, and Dynamics , Vol. 42, No. 2, 2019, pp. 244–259. https://doi.org/10.2514/1.G003497.[15] Wang, X., and van Kampen, E., “Incremental Backstepping Sliding Mode Fault-Tolerant Flight Control,”
AIAA Scitech 2019Forum , American Institute of Aeronautics and Astronautics, San Diego, California, 2019, pp. 1–23. https://doi.org/10.2514/6.2019-0110, URL https://arc.aiaa.org/doi/10.2514/6.2019-0110.[16] Stevens, B. L., and Lewis, F.,
Aircraft Control and Simulation , Wiley, New York, 1992.[17] Waszak, M. R., and Schmidt, D. K., “Flight Dynamics of Aeroelastic Vehicles,”
Journal of Aircraft , Vol. 25, No. 6, 1988, pp.563–571. https://doi.org/10.2514/3.45623.[18] Theodorsen, T., “General Theory of Aerodynamic Instability and the Mechanism of Flutter,” Tech. rep., NACA TechnicalReport No. 496, 1935.[19] Bisplinghoff, R., Ashley, H., and Halfman, R.,
Aeroelasticity , Addison-Wesley Publishing Company, 1955.[20] Meirovitch, L., and Tuzcu, I., “Integrated Approach to the Dynamics and Control of Maneuvering Flexible Aircraft,” Tech. rep.,NASA Langley Research Center, 2003.[21] Qi, P., and Zhao, X., “Flight Control for Very Flexible Aircraft Using Model-Free Adaptive Control,”
Journal of Guidance,Control, and Dynamics , Vol. 43, No. 3, 2020, pp. 608–619. https://doi.org/10.2514/1.G004761, URL https://arc.aiaa.org/doi/10.2514/1.G004761.[22] Qi, P., Zhao, X., Wang, Y., Palacios, R., and Wynn, A., “Aeroelastic and Trajectory Control of High Altitude LongEndurance Aircraft,”
IEEE Transactions on Aerospace and Electronic Systems , Vol. 54, No. 6, 2018, pp. 2992–3003.https://doi.org/10.1109/TAES.2018.2836598, URL https://ieeexplore.ieee.org/document/8359337/.[23] Wang, X., Van Kampen, E. J., and Chu, Q., “Comparison of three control structures for inducing higher-order sliding modes,” , 2019, pp. 4073–4080. https://doi.org/10.23919/ECC.2019.8795872.[24] Edwards, C., and Shtessel, Y. B., “Adaptive continuous higher order sliding mode control,”
Automatica , Vol. 65, No. March,2016, pp. 183–190. https://doi.org/10.1016/j.automatica.2015.11.038.[25] Lu, P., van Kampen, E., de Visser, C., and Chu, Q. P., “Aircraft fault-tolerant trajectory control using Incremental NonlinearDynamic Inversion,”
Control Engineering Practice , Vol. 57, No. December, 2016, pp. 126–141. https://doi.org/10.1016/j.conengprac.2016.09.010.[26] Khalil, H. K.,
Nonlinear Systems , Prentice-Hall, New Jersey, 2002.[27] Etkin, B.,
Dynamics of Atmospheric Flight , Dover Publications, Toronto, 2005. https://doi.org/10.1007/s13398-014-0173-7.2.[28] Wang, X., van Kampen, E., and Chu, Q. P., “Gust Load Alleviation and Ride Quality Improvement with Incremental NonlinearDynamic Inversion,”
AIAA Atmospheric Flight Mechanics Conference , American Institute of Aeronautics and Astronautics,Grapevine, Texas, 2017, pp. 1–21. https://doi.org/10.2514/6.2017-1400., American Institute of Aeronautics and Astronautics,Grapevine, Texas, 2017, pp. 1–21. https://doi.org/10.2514/6.2017-1400.