Data-driven nonlinear aeroelastic models of morphing wings for control

Accurate and efficient aeroelastic models are critically important for enabling the optimization and control of highly flexible aerospace structures, which are expected to become pervasive in future transportation and energy systems. Advanced materials and morphing wing technologies are resulting in next-generation aeroelastic systems that are characterized by highly-coupled and nonlinear interactions between the aerodynamic and structural dynamics. In this work, we leverage emerging data-driven modeling techniques to develop highly accurate and tractable reduced-order aeroelastic models that are valid over a wide range of operating conditions and are suitable for control. In particular, we develop two extensions to the recent dynamic mode decomposition with control (DMDc) algorithm to make it suitable for flexible aeroelastic systems: 1) we introduce a formulation to handle algebraic equations, and 2) we develop an interpolation scheme to smoothly connect several linear DMDc models developed in different operating regimes. Thus, the innovation lies in accurately modeling the nonlinearities of the coupled aerostructural dynamics over multiple operating regimes, not restricting the validity of the model to a narrow region around a linearization point. We demonstrate this approach on a high-fidelity, three-dimensional numerical model of an airborne wind energy (AWE) system, although the methods are generally applicable to any highly coupled aeroelastic system or dynamical system operating over multiple operating regimes. Our proposed modeling framework results in real-time prediction of nonlinear unsteady aeroelastic responses of flexible aerospace structures, and we demonstrate the enhanced model performance for model predictive control. Thus, the proposed architecture may help enable the widespread adoption of next-generation morphing wing technologies.


Introduction
It is expected that highly flexible aeroelastic structures will become ubiquitous in future transportation and energy systems, enabled by advanced materials and emerging morphing wing technologies [1]. Indeed, more responsive and deformable aerodynamic surfaces may have transformative impact in efficiency and maneuverability, as demonstrated by the incredible performance of biological flight systems [2][3][4][5]. These flexible aerospace structures pose a considerable modeling challenge, as they involve highly coupled and nonlinear interactions between the aerodynamic and structural dynamics. Existing aeroelastic models are frequently either linearized about a single operating condition or involve expensive high-fidelity numerical simulations to resolve the relevant spatial and temporal scales [6,7]. In this work, we develop a flexible data-driven modeling architecture, based on the recent dynamic mode decomposition with control (DMDc) [8], to model highly coupled and nonlinear aeroelastic dynamics over the entire flight envelope. This work is motivated by a particularly compelling application in airborne wind energy, although the methods discussed are generally applicable to highly coupled aeroelastic systems. To achieve these flexible and tractable models, we introduce several innovations to the DMD architecture that may be used more broadly. Finally, we demonstrate the efficacy of the resulting models for model predictive control (MPC).
The development and deployment of innovative renewable energy technologies is essential to reduce greenhouse gas emissions. Airborne wind energy (AWE) is a promising technology that extracts power from high-altitude winds using tethered drones. Currently, these drones are equipped with conventional rigid-wings, relying on hinged control surfaces [9,10]. Replacing the rigid-wings with shape-adaptable or so-called morphing wings has the potential to improve power production by allowing the drone to smoothly adapt to changing flight conditions, thus enabling optimal performance over the full operational regime [11]. Such morphing wings are inherently flexible, leading to a tight coupling of their aerodynamic, structural, and rigid-body responses [12,13]. This coupling, combined with the vast flight regime such AWE drones operate in, makes modeling and controlling these systems a considerable challenge.
Ground-based power-generator AWE systems, in particular, aim to extract power by periodically reeling-out (traction phase) and reeling-in (retraction phase) the tether. Thus, the drone constantly changes between two distinct operating modes. In the traction phase, maximum power is produced by operating the drone at high flight speeds, large incidence angle, and high-lift forces; whereas in the retraction phase, the load on the tether is minimized by decreasing the incidence angle to reduce the required reel-in power. The drone is therefore required to operate both at high and low wing loading over a wide range of wind speeds, while simultaneously following a desired trajectory, thus creating the need for highly adaptable drones. To design, analyze, and control such a high-dimensional, nonlinear dynamical system, efficient numerical models are needed. Opposed to earlier work on flutter [14][15][16][17][18] or recent work on high aspect ratio wings [19][20][21], these models need to be valid over a larger set of operational conditions and flight speeds.
A compelling new family of methods capable to approach this problem are emerging data-driven modeling techniques enabling the characterization of such high-dimensional, nonlinear dynamical systems [22][23][24]. Dynamic mode decomposition (DMD) is a particularly promising technique, enabling the discovery of dynamical systems from high-dimensional data by decomposing complex dynamics into simple representations based on spatiotemporal coherent structures [22,[25][26][27]. By using measurements of the system, DMD extracts the dynamics without the need to know the underlying equations, opposed to physics-based models build on established mathematical and physical laws. DMD was first introduced by Schmid [25] and has since gained traction for modeling systems exhibiting nonlinearities. DMD is strongly connected to the Koopman operator, which is an infinitedimensional linear operator representing nonlinear dynamical systems [26]. Therefore, DMD is a promising candidate to model the nonlinear dynamics inherent to morphing AWE drones.
In this work, we present an efficient and accurate unsteady aeroelastic reduced-order model (ROM) for flexible structures, applicable to AWE morphing wing drones. The method is outlined in figure 1. Two innovations are introduced. First, an extension of the DMD method is developed for algebraic differential equations to generate a reduced-order unsteady aerodynamic and aeroelastic model. To generate the ROM, the algorithm sequentially applies a mode superposition method on a detailed 3-D structural finite element model [28], followed by using the extended algebraic DMD method on a coupled structural and unsteady 3-D panel method [29,30]. Second, an interpolation scheme smoothly connecting multiple ROMs valid at different flight velocities is introduced. This allows accurate faster-than-real-time prediction of the nonlinearities of the coupled aerostructuraldynamics of flexible AWE drones over the full nonlinear flight regime.   [31,32]), exploiting camber-morphing for roll-control [13,33,34]. II. AWE system modeling, consisting of the coupled flight dynamic and aeroelastic models. III. Reduced-order model training phase, consisting of a. extracting the most important structural and morphing modes, b. training the model from snapshots of the states x of the aeroelastic model excited through impulses in the inputs u for a set of flight speeds. IV. Data-driven reduced-order modeling method, consisting of a. generating the parametric algebraic DMD with control (aDMDc) model as in § 3 (depicted is the doublet distribution on the wing for a number of aDMDc modes) and b. applying the parametric ROM to calculate the aerodynamic forces and moments F aer o generated by the morphing wing for a specific input u.
The paper is organized as follows. In Sec. 2, the state-of-the-art in modeling flexible aerospace structures, reduced-order modeling, DMD, and MPC is reviewed. In Sec. 3, the proposed modeling approach is introduced, specifically the extension of DMD for algebraic differential equations and to interpolate between local linear models. In Sec. 4, a number of numerical examples, including a NACA0012 rigid-wing, a morphing AWE wing, and an MPC test case are discussed, highlighting the applicability of the ROM to aerodynamics-only, aeroelastic, and control problems. In Sec. 5, the results and the potential applicability for modeling and controls of general aerospace structures and other dynamic systems is discussed.

Modeling of flexible aerospace structures
The interaction between structural deformations and aerodynamic forces has long been recognised in the field of aeronautic as being of paramount importance. The first wind tunnels at NACA/NASA were specifically dedicated to aeroelastic studies. Early flight suffered from aeroelastic issues, and as flight speed increased, it was not possible to neglect these effects [35]. Now, it is common practise to consider aeroelastic effects early in the design, to avoid expensive redesigns. Due to the different equations for the structure and for the aerodynamics, the two problems are typically modeled with separate techniques, to be coupled later with an appropriate scheme [30]. Splines are usually used [36] to interpolate the structural displacements onto the aerodynamic grid, and the aerodynamic forces onto the structural nodes. Nonlinear beam models are commonly used to describe the characteristic of flexible structures with a dominant spatial dimension [12,19,20]. The sectional properties of such structures are usually pre-calculated along their dominant direction. However, more refined models are required in the case of morphing and geometrically complex wings that exhibit flexibility in the chordwise direction and, therefore, strongly interact with the aerodynamics. Detailed finite element models, based on both beam and shell elements, are therefore used to accurately represent the characteristics of such structures [29]. To increase computational efficiency, model reduction based on modal decomposition techniques are often applied [28,29].
In terms of the aerodynamics, different models have been applied to fluid-structure interaction (FSI) problems, depending on the tradeoff between computational cost and accuracy of the simulation. Early studies considered simple 2D geometries with analytical models for the aerodynamics, based on unsteady potential flow theory [15]. These methods rapidly evolved, and extensions to 3D problems, based on strip theory, are still used today [37]. However, it is now common to use the doublet lattice method (DLM) [38] for the unsteady aerodynamic generalised forces, and its steady counterpart the vortex lattice method (VLM) [39,40], for aeroelastic analysis. Compared to computational fluid dynamics based on the Navier-Stokes equations, these methods are significantly more efficient, and they provide results that are accurate enough for the early design stages. These methods do not require discretising the volume surrounding a body, but instead reduce the problem to an equivalent formulation on the boundaries of the domain, so that only the wing surface must be discretised. Therefore, these methods do not suffer from issues related to mesh deformations [41]. On the other hand, they do not represent viscous effects and are not suited for transonic applications, with the exception of Morino's method [42].
Panel methods are mainly divided into frequency-domain and time-domain methods. Representing unsteady aerodynamics in the frequency domain is useful for flutter predictions. A linear state space model is usually preferred for response analysis and control design [43][44][45]. In the field of AWE, the benchmark problem for the method in this paper, researchers have applied several timedomain models, ranging from lifting line methods [46] to quasi-steady approximations of a 3D panel method, based on source and doublet distributions [11]. However, none of these approaches considered the full unsteadiness of the flow.
Murua et al. [47] proposed a promising approach based on the unsteady vortex lattice method (UVLM). The UVLM is the direct extension of VLM in the unsteady time domain. In this specific case, linearising the nonlinear equations for lifting surfaces and aerodynamics, the state-space form is easily obtained. Additionally, it is straightforward to include the flexible-body dynamics in the simulation [20]. In practise, the full equations of the UVLM can be summarised as: where f , g and h are general nonlinear functions, x is the vector of states, composed of structural and aerodynamic nodes positions, u is the vector of inputs (e.g. the flight condition) and y is the vector of outputs (the aerodynamic forces). The superscript indicates the time step. Note that these equations depend on the control inputs at the next time step as a result of the impulsive part of the aerodynamics. Usually, when artificial aerodynamic states are added and a state-space form is obtained, this is reflected in the feed-trough matrix to the outputs. However, if the correct evolution of the states must be considered, the immediate effect of the inputs must be modeled. Under the assumption of small displacements around the trim condition, thus the deformed equilibrium, the equations can be linearized, resulting in a state space form that can be coupled with the structural equations: The limitations of UVLM are related to the approximation of the surface as an infinitely thin sheet. If the camber effect is important, as in the context of morphing structures or flexible airfoils, this requires more detailed models, such as the unsteady 3D panel method based on both sources and doublets [30]. In this case, the equations for the aerodynamics are: where A b is the aerodynamic influence matrix of the body, A w is the wake influence matrix, µ is a vector containing the doublet strengths, and w is the downwash at the control points. The wake node positions are then evolved with a simple advection equation, using the local velocity. These equations are similar to 1a, and the computation of the forces are governed by a similar expression. Both the UVLM and the 3D panel methods involve algebraic-differential equations. It is essential that our ROMs respect this structure.

Reduced-order aeroelastic models
The methods above provide accurate models of aeroelastic effects. However, faster models are often necessary for optimisation and control, even at the expense of some fidelity. This tradeoff between accuracy and efficiency has motivated reduced-order models, which model the behaviour of the system with as few states as possible. There is a wide variety of model reduction techniques in the literature [23]. Both data-driven and semi-analytical approaches are common. The latter is exemplified by the modal reduction of structural dynamics: an eigenvalue problem is solved and a reduced set of orthogonal modes are used to describe the state of the system. Other analytical methods project the governing equations onto a set of data-driven modes [48,49], for example obtained via the proper orthogonal decomposition (POD); it maybe possible to develop parametric models, although this procedure may be cumbersome.
In recent years, data-driven approaches have become increasingly powerful and widely adopted. Many data-driven modeling techniques may be categorized as system identification, where a model for the input-output behavior is constructed based on data. Time domain techniques are especially common, such as the aerodynamic impulse response (AIR) [50][51][52][53], where the system is perturbed with an impulsive input, and the output response may be used to predict the response to future input maneuvers via convolution. Although these techniques are typically linearized, nonlinear kernels may also be employed in a Volterra series [51,54]. State-space realizations are becoming increasingly common, especially for control applications [43,45]. The eigensystem realization algorithm [55] was developed specifically to model the structural response of aerospace structures, although the resulting model is formulated in terms of a nonphysical state that is difficult to interpret. However, when the entire physical state can be observed, it is possible to construct state-space models in terms of a reduced state that is related to the physical domain. The dynamic mode decomposition (DMD) [22,[25][26][27] and the sparse identification of nonlinear dynamics (SINDy) [56][57][58] both result in physically interpretable models of the dynamics, and they have been extended to input-output systems for control [8,59]. DMD results in a linear time-invariant (LTI) system, which is valid in a particular operating regime. There are several approaches to interpolate LTI systems to obtain a linear parameter varying (LPV) system, as a function of the flight condition [60][61][62][63][64]. However, it is generally challenging to interpolate models; indeed, the direct interpolation of the matrices only works for small systems [65], but for large systems the eigenvalues of the system become unstable [66]. If the transfer function is interpolated, assuming oscillatory movement, this problem is solved, but there are still limitations related to the change of an oscillatory couple of poles, to two real-valued poles [60,67]. In general, there is no universal solution for the interpolation of LTIs, and this work will develop an interpolation scheme specifically for DMD models.

Dynamic mode decomposition with control
In this work, we will develop an interpolation scheme to combine local models that are obtained via dynamic mode decomposition with control (DMDc) [8]. There are several advantages to DMDc over other linear state-space modeling techniques, primarily the formulation of the models in terms of physically interpretable modes [22]. Further, it has been shown by Kaiser et al. [59] that DMDc models may be nearly as effective as the full nonlinear dynamics for model predictive control, even in chaotic systems.
In DMDc, it is assumed that the entire state and all control inputs may be observed to train the model, and the evolution of this state may be expressed as a linear system where A and B are unknown constant matrices, x is the state vector, and u is the input vector. The unknown matrices A and B are solved through a regression procedure, based on time-series measurement data: The original dynamics may be written in terms of these data matrices as and Ω = [X T Υ T ] T . The matrix of known terms Ω can be approximated by a singular value decomposition (SVD) as Ω = UΣV T . This expression may be truncated, taking into account only the vectors associated with the largest elements of Σ. This gives the approximation where the subscript t denotes truncated matrices. It is then possible to solve for the unknown terms in Ψ via the pseudo-inverse of Ω: It is possible to build a reduced-order model by projecting the dynamics onto the leading modesÛ t from the SVD of The resultingΨ matrix can then be split intoÃ andB, so that: where x =Û t q, and q is the vector of mode amplitudes.

Model predictive control
In this section, we introduce the MPC problem, outline its potential benefits, and discuss the importance of developing efficient ROMs to allow real-time control. MPC is an effective model-based control, which has revolutionized the industrial control landscape [68], as it enables the control of strongly nonlinear systems with constraints, which are difficult to handle using traditional linear control approaches [69][70][71]. MPC computes a control sequence u(x j ) = {u j+1 , ..., u j+m c }, given a measurement x j , by solving a constrained optimization over a receding horizon T c = m c ∆t. At each time step, the optimization is repeated, updating the control sequence over the control horizon and applying the first control action u j+1 to the system. The optimal control sequence u(x j ) is obtained by minimizing a cost function J over a prediction horizon T p = m p ∆t ě T c . The cost function is subject to the discrete-time dynamics and other constraints. The cost function thus penalizes deviations of state from the reference trajectory, the control expenditure, and the rate of change of the control signal, with each term weighted by the matrices Q, R u , and R ∆u , respectively. To enable this optimization loop to run in real-time on a flight controller, MPC relies on efficient models and high-performance computing. In AWE and general flight systems, the control must rapidly respond to disturbances on short time-scales, as time delays from sensors, signal transduction, or processing can destroy the robustness of feedback control, putting limitations on the achievable performance [72]. Currently, MPC is not well suited for control of such complex, high-dimensional system with fast timescales, especially with power-constrained or computationally limited hardware available on state-of-the-art AWE drones or lightweight flexible aircraft. The short timescales associated with agile flight in a complex unsteady fluid environment make it challenging to arrive at control decisions with the small latency required for robust performance. Gust disturbances are very rapid, and computing and estimating the dynamics of the unsteady fluid requires considerable computational resources. Therefore, instead of a detailed model-based feedback control strategy that spans all relevant timescales, it is essential to develop ROMs that balance accuracy and efficiency. The benefit of MPC lies in simple and intuitive tuning and the ability to control complex phenomena, including systems with time delays, non-minimum phase dynamics, and instability. It is also simple to include known constraints and multiple operating conditions, and it provides the flexibility to formulate and tailor a control objective. The major challenge of MPC lies in the development of a suitable model via existing system identification or model reduction techniques [59,73], which may require expensive and time-consuming data collection and computations. Nonlinear models based on machine learning, such as neural networks, are increasingly used due to advances in computing power, and recently deep reinforcement learning has been combined with MPC [74,75], yielding impressive results in the large-data and high-performance computing limit. However, many modern techniques in machine learning (e.g., neural networks) rely on access to massive data sets, have limited ability to generalize, do not readily incorporate known physical constraints, and require expensive and time-consuming computations. Instead, Kaiser et al. [59] showed that a simple linear model obtained via dynamic mode decomposition with control (DMDc) [8] performs nearly as well with MPC as a full nonlinear model, and may be trained in a surprisingly short amount of time.
Currently, several studies have applied MPC to AWE, although they rely on simplified models [76][77][78][79][80][81]. For flexible aerospace structures [82][83][84][85][86], most control strategies rely on linearised models valid only around a single steady state trim position or on successive linearisation of the underlying system. These models do not allow for optimal control actions over a large range of flight speeds. Thus, there is a need to develop fast and accurate ROMs covering a large range of flight speeds, which can be used for MPC.

Methods
The models and methods used for the current research are presented in this section. The full structural model is based on a linear structure, composed of beams and shells, while the aerodynamics are modeled using an unsteady 3D panel method. The coupling is obtained via a thin plate spline (TPS) [87] and inverse distance weighting (IDW) interpolation [88] and the integration in time is performed using the Newmark method [89] with a tight coupling between structural part and aerodynamic model. The full model is described in detail in the Appendix.
Taking inspiration from the approach of [47], we develop a time-domain ROM. However, instead of a mathematically formal linearisation, we leverage the data-driven DMDc [8]. This method has been applied to a wide variety of problems, ranging from fluid dynamics to disease modeling [22,23], and it has shown promise in the model predictive control of complex systems [59]. In the following, we develop extensions to DMDc to handle algebraic-differential systems and to interpolate between multiple operating regimes, extending the applicability of these methods to aeroelastic systems over a large operating regime.

Algebraic dynamic mode decomposition method aDMDc
Here we extend the DMDc algorithm to handle algebraic-differential equations, as are found in aeroelastic problems: f (x n+1 , u n+1 ) = g(x n , u n ).
We seek a linear system of the form: q n+1 =Ãq n +Bu n +Fu n+1 .
Therefore, we will modify the standard DMDc to include the F matrix. The modification of the system is straightforward. In addition to the standard DMDc data matrices, we include the matrix of inputs shifted in time: Again, we define: Thus, it is possible to write the dynamics as: In order to solve for Ψ, we compute the SVD of Ω: and invert Ω via the pseudo-inverse, as with DMDc. As before, we truncate the SVD using the optimal hard threshold criteria of Gavish and Donoho [90]. We then obtain: The matrices U 1 , U 2 , and U 3 are obtained by splitting the truncated SVD U t according to the structure of Ω in Eq. (14), as in DMDc, andÛ t are the leading SVD modes of X 1 «Û tΣtV T t . This procedure, results in a data-driven approximation of the system in Eq. (12). It is possible to go between the reduced state q and the full physical state x via the modesÛ t .
The procedure to generate the ROM is summarized in Algorithm 1.

aDMDc for flexible aerospace structures
In this research, we will apply the aDMDc algorithm from above to aeroelastic systems. Instead of developing separate ROMs for the aerodynamic and structural models, we design a single monolithic ROM for the coupled aeroelastic system. This approach has the benefit of reducing the dimension of the external inputs to the system. If we developed two isolated models, then all of the structural modes would be inputs to the aerodynamic system. Instead, in our coupled model, only the terms influencing the flight condition of the body are inputs to the system: where α is the angle of attack, p, q, and r are the roll, pitch, and yaw rotation rates and F i is a general normalised actuation force (with values from 0 to 1) required for the morphing wing actuation inputs. The sideslip angle is not included as an input because the effect was considered negligible, but it can easily be added. Before applying aDMDc, we conduct a modal analysis to identify the most important structural modes required to describe a general motion. The ROM is then developed in terms of these reduced structural modes and the aerodynamics. The state x of the system is then: where η contains the amplitudes of the structural modes that are retained and µ contains the doublet strengths. Only the doublets required to compute the aerodynamic forces are retained, thus only those on the lifting surface and not those of the wake. The data-driven method is trained using impulsive inputs, although it is possible to use other input sequences. If the total forces on the body are of interest, for example if the aeroelastic system must be included in a multibody framework, or we want to study the stability and control derivatives, the distribution of coefficient of pressure can be computed with the full nonlinear expression, described in more detail in the appendix. Indeed, the relation between doublet distribution and pressure distribution is a nonlinear relation and, in order to retain accuracy, it can be used to obtain the forces. The procedure to generate the ROM is summarized in Algorithm 2.

Parametric aDMDc
The method above is able to reproduce the entire aeroelastic system at one flight velocity. Thus, we will construct a separate aDMDc ROM for each velocity condition. For our AWE application, it is necessary to interpolate multiple ROMs across a wide range of flight conditions. There are several valid approaches to interpolate multiple LTI systems into an LPV system. Due to the large state dimension, we opt not to interpolate the high-dimensional model, as it is likely to result in unstable models. Instead, because the low-dimensional aDMDc model state q may be used to reconstruct the high-dimensional state x, we run multiple local ROMs independently and then interpolate the high-dimensional state using a spline method. This is a distinct advantage of the DMD approach over other subspace identification methods. This approach bypasses the issues with standard interpolation schemes, related to the assumption of slow varying scheduling parameter [61]. A graphical explanation of this procedure is shown in figure 2.

Results
In this section, several applications of the proposed method are outlined. First, we present a simple example illustrating the need for the modified aDMDc algorithm. Next, we demonstrate the approach to model the aerodynamics on a rigid NACA0012 wing. We then extend this approach to challenging aeroelastic morphing wing, both at a single flight condition and interpolated across several flight conditions. Finally, the aDMDc model is used for MPC, demonstrating the ability to effectively control a nonlinear system with constraints.

Generic differential-algebraic system of equations
The addition of the F matrix is important to retain the stability of the system. This is better understood with the following example, assuming known dynamics of the form: Applying the standard DMDc algorithm, with only the A and B matrices, we obtain: The resulting dynamics are not only incorrect, but they are unstable. In contrast, the aDMDc method, results in the correct, stable dynamics of the system: Note that this system consists of only two states, thus there are no truncation or projections involved in the procedure. However, this first example shows the necessity of adding the F matrix to correctly capture the dynamics of the algebraic-differential system.

NACA0012 rigid wing
We now apply aDMDc to model the aerodynamics of a rigid, unswept, untappered, planar NACA 0012 wing with chord c = 0.4 m and aspect ratio AR = 10. The wing surface is divided into trapezoidal panels, with 98 panels along the span and 38 along the chord. The discretisation is defined so that no panels are present with high aspect ratio and the results are mesh independent.
To initialize the simulation, the wing is set into motion with a constant speed of V 8 = 50 m/s, a constant angle of attack of α 0 = 2˝and no rotation rates. When the steady state is reached, the training phase for the ROM is initiated. In this rigid airfoil case, the system state only consists of the doublet strengths on the surface of the body. A set of impulsive inputs in the α, p, q, and r are provided in a random order. The impulse amplitudes are ∆α = 8˝and ∆p = ∆q = ∆r = 30˝/s, which are large enough to excite the dynamics, but small enough that the system may be considered as a linearisation around a single flight condition. Both positive and negative impulses are used.
Algorithm 1 is used to obtain a ROM, using the snapshots matrices of the states and inputs recorded in the training phase. For a first comparison, this ROM is used to reproduce the behaviour of the system during the training phase. Afterwards, two different testing phases are performed. In these phases, the ROM and the full model are compared. The first test set consists of a sinusoidal roll rate with reduced frequency of 0.1 and amplitude ∆p = 18˝/s and a sinusoidal change in the angle of attack with reduced frequency of 0.2 and amplitude of ∆α = 4˝. The second test set consists of a random sequence of impulsive inputs, with the same amplitude as the training phase, but with different order. This case is challenging, as all frequencies are excited by the impulses, fully exploring the limitations of the ROM.
The aerodynamic forces and moment coefficients acting on the wing, for all the three phases are shown in figure 3. The agreement between the full and reduced-order model is excellent. Importantly, ROM is three orders of magnitude more efficient than the full model, as highlighted in Tab. 1, which shows the relative errors and the speed-up factor for the two testing phases.

Morphing AWE wing
In this example, we use aDMDc to model a coupled aeroelastic system, given by a flexible and highly cambered wing. The wing is the result of studies performed in the context of the ftero project at ETH Zurich [31,32]. The wing is composed of a CFRP skin and an internal structure based on a Voronoi tessellation. The trailing edge of both the right and left sides of the wing can morph, increasing or decreasing the camber, thus replacing conventional ailerons. For flexible structures, it is first necessary to reduce the number of states characterising the wing structure, before training the aDMDc ROM. The full finite element morphing wing model is generated using the commercial software Nastran. The mass and stiffness matrix of the structure is extracted and a modal decomposition is performed. We retain only the first 8 structural modes. Two additional modes are added corresponding to the deformation modes due to morphing. Afterwards, the aerodynamic mesh is generated, as in the previous example by discretizing the surface with trapezoidal panels. We then construct the ROM following Algorithm 1. In this case, two additional inputs are included to actuate the morphing surfaces, and the state of the system consists of the doublet strengths on the surface of the body and the structural modal amplitudes.
The linearisation point is the same as before: V 8 = 50 m/s, α 0 = 0˝and p 0 = q 0 = r 0 = 2˝/s. The same set of impulses for α, p, q and r is used, with the same amplitudes, resulting in a single ROM that is sufficient at the given velocity. In addition, the morphing wing actuation is impulsively fully actuated during the training phase. Two testing phases are considered, as before, consisting of sinusoidal and impulsive inputs.
The aeroelastic responses are shown for all three phases in figure 4. In all cases, the agreement between the full and reduced-order model is excellent. Table 2 summarizes the speed-up factor between the full model and the ROM and the relative errors in the predicted response.

Parametric reduced-order aeroelastic model
Large changes in operating conditions generally pose a challenge for ROMs. Indeed, when the operating condition changes enough, the ROM loses its validity and another training phase must be performed. In this example, we demonstrate the ability to interpolate several ROMs that are identified at different speeds. The result is a highly flexible model that seamlessly covers a large range of flight conditions and results in high performance for continuously varying operating conditions. We only need to interpolate models that are identified for different flow velocities, as the other parameter variations, such as angle of attack, are captured well by a single ROM for a given velocity. A set of five ROMs are generated for the ftero wing, following the scheme presented above, for velocities in the range V 8 = 35´80m/s. These ROMs are then used to predict the behaviour of the system using the reference set of inputs, with a linear velocity sweep across the entire range. The same inputs are also used for the full model, and the results are presented and compared in figure 5.
Again, there is excellent agreement between the full-order and reduced-order models. Furthermore, even after the interpolation is introduced where multiple systems are solved in parallel, the ROMs are still more than three orders of magnitude faster than the full model. These results are summarized in Tab. 3.   Gust load alleviation MPC. Maintaining the wing stress level below a critical threshold is crucial to guarantee structural integrity of the wing. The stress level in the wing is directly proportional to the modal amplitudes of the wing structure, especially to the first bending mode [82]. Therefore, the MPC objective in this first example is to keep the first bending mode amplitude close to zero, ensuring structural integrity of the wing at gust encounter. This goal is achieved by actuating the morphing wings symmetrically, thereby adjusting the lift and the wing root bending moment. The first bending mode is available for feedback in this study, which could be estimated for example with strain gauges placed at the wing root. Flying at a constant velocity of V = 60m/s and with the first bending mode being a physical state of the aDMDc model, the nonlinear calculation of the aerodynamic forces and moments can be omitted, resulting in linear model. The system time step is ∆t sys = 0.006s and the MPC time step is ∆t model = 0.018s, representing approximately a state-of-theart flight control system actuator update frequency [91]. The weighting matrices are Q = 10, R u = 0, and R ∆u = 1, the actuation is limited to u P [´1, 1], ∆u P [´0. 18, 0.18], and the control and prediction horizon are set to m p = m c = 10. The wing encounters a series of three gusts with constant amplitude ∆α = 1.5 o and increasing gust length t g,1 = 0.25s, t g,2 = 0.5s, and t g,3 = 1.0s. Figure 6 shows the gust induced angle of attack in the upper plot, the modal amplitude deviation of the first bending mode calculated with the high fidelity model in the middle plot, and the actuation of the morphing wing in the lower plot. First, the MPC is switched on (blue continuous line) and second, the MPC is switched off (red dashed line). The expected delay of the MPC is clearly visible, especially prominent for short gust lengths, which could be dealt with using additional gust-induced angle of attack sensors [82]. Nevertheless, the performance of the controller is promising, achieving modal amplitude peak reductions of ∆q 0.25s = 51.2% ∆q 0.5s = 66.5%, and ∆q 1.0s = 79.5%. Lift force tracking MPC. In the second test case, the MPC objective is to track a specified lift coefficient. We thus use the full aDMDc ROM with the nonlinear calculation of the lift. The lift is directly available for feedback in this study, which could be measured on the drone for example by a five-hole probe [92]. Tracking a specific lift is similarly achieved by symmetrically actuating the morphing wing to reduce or increase the lift coefficient. The MPC parameters are the same as in the first example, except the weighting matrices are set to Q = 1 and R ∆u = 10000. Figure 7 shows the desired normalized lift coefficient in the top plot, the error defined as the deviation of the actual from the desired normalized lift assessed by the high fidelity model in the middle plot, and the morphing actuation in the bottom plot. Again, the MPC is first switched on (blue continuous line), and then switched off (red dashed line). The MPC performs well, even when the commanded lift exceeds the attainable lift and the actuation is saturated.

AWE trajectory and gust load alleviation MPC.
The third test case is AWE specific, where the structural integrity of the tether is maintained by keeping the tether force below a critical maximum force. Currently, the tether force is controlled by changing the ground station reel out speed; however, this approach is too slow to mitigate force peaks introduced by gusts on the wing. Therefore, either the maximum allowable tether force must be reduced, leading to lower power production, or active load alleviation must be achieved on the wing. In this example, we use MPC to keep the lift coefficient on the wing close to a desired steady state value. The drone flies an AWE specific trajectory, with flight velocities ranging between V = [40, 80]m/s, shown in figure 8. Additionally, the wing encounters a gust at the maximum flight speed, with t g = 0.5s and ∆α = 1 o . In this case, an AWE specific lift coefficient trajectory is defined, with a higher lift coefficient in the traction phase at high flight speeds and lower lift coefficient at low flight speeds in the retraction phase. This is achieved by both changing the drone's angle of attack for large and slow lift coefficient changes (assumed to be controlled by an extended fixed wing drone controller using the drone elevator [29]) and by symmetrically actuating the morphing wing for fast gust load alleviation. To accurately control the system over the large range of flight speeds, the parametric aDMDc-based ROM is incorporated in the MPC. As the objective of the MPC is similar to case 2, the MPC parameters are the same as in the previous example.
The results of this test case are shown in figure 8. We compare the MPC performance with the parametric aDMDc model (blue continuous line), and the MPC performance with a non-parametric aDMDc ROM obtained at the mean flight speed and assuming a quadratic depencency of the lift force on the flight velocity (red dashed line). Thus, nonlinearities introduced by the flexibility of the morphing wing and by the reduced frequency of the unsteady aerodynamics are neglected. The results highlight the importance of using a parametric ROM covering the full operational regime to precisely track a desired lift coefficient when dealing with large flight speed ranges specific to AWE. The non-parametric ROM predicts the lift coefficient inaccurately and therefore saturates the actuator. Therefore, it not only fails to track the correct lift coefficient but also prevents any gust load alleviation. This could potentially lead to tether rupture and loss of drone, further motivating the parametric ROM.

Discussion and conclusion
In this work, we present a data-driven reduced-order modeling framework for flexible aerospace structures that is valid over a large range of operating conditions. Our models are based on the recent DMDc algorithm, and we introduce two methodological extensions to this algorithm to broaden its applicability. First, we generalize the formulation to handle algebraic equations, and second, we develop an interpolation scheme to smoothly connect several models valid in different operating regimes. The method is demonstrated on several test cases, ranging from a generic differentialalgebraic system of equations to a high fidelity, three-dimensional numerical model of a morphing wing AWE system. The reduced-order models are faster than real time and are able to accurately predict nonlinear, unsteady aeroelastic responses. Furthermore, the model is incorporated into an MPC framework, demonstrating the potential to perform load control of morphing wings in AWE typical flight conditions, characterized by large changes to the load and flight speed.
The results of this study show that our extension to the DMDc algorithm to handle algebraic equations is vital for unsteady aerodynamic and aeroelastic systems, not only to capture the correct dynamics but especially to retain stability of the system. The ROM shows excellent agreement with the full-order model for unsteady aerodynamic lift, drag, and moment coefficient prediction, with relative errors below 0.5%. For the coupled aeroelastic system, the agreement is similarly close, even for large random inputs, with maximum errors for lift and drag below 2% and for the moment coefficients below 6%. Further, we show that it is possible to smoothly interpolate between several aDMDc ROMs, resulting in a parametric ROM that is valid over a large range of operating conditions and that agrees closely with the full model. The errors are below 2% for the lift and drag and below 5% for the moment coefficients. In all cases, the ROMs are more than three orders of magnitude faster than the full model, enabling use for MPC. On three test cases, the model performance is evaluated, showing excellent results when performing load control of the morphing wing AWE system. Especially for the AWE specific flight trajectory with gust rejection, the parametric ROM is able to fully mitigate the gust induced load peak, which is not achievable by a single non-parametric ROM.
Although the data-driven modeling procedure introduced in this work is generally applicable, it was specifically motivated by AWE applications. The resulting ROMs will improve the prediction performance of current AWE models, which are mostly based on simplified quasi-steady aerodynamic models. This will lead to better performing AWE drones with enhanced power production capabilities and potentially reduce the levelized cost of energy, crucial in adopting the emerging AWE technology. Furthermore, incorporating the ROM into existing AWE-specific MPC frameworks with power production objectives [81], and enhancing them with additional structural stress level constraints, could potentially lead to drastic improvements in control performance and therefore power production. Apart from AWE specific applications of the model, this procedure may be applied to any flexible aerospace structure operating over multiple operating regimes, and more generally to any parametric dynamical system across multiple operating conditions.
In the standard unsteady formulation, during the integration in time, the state of the system continuously increases as new rows of the wake sheet are added. This increases the computational time for the integration at each time step and, for this reason, it was decided to truncate the wake. The truncation has been performed far from the body, and a detailed study on the influence of this on the generalised forces has been conducted. It was concluded that, if the truncation point is far enough, the influence on the force can be neglected.