The Balanced Mode Decomposition Algorithm for Data-Driven LPV Low-Order Models of Aeroservoelastic Systems

A novel approach to reduced-order modeling of high-dimensional time varying systems is proposed. It leverages the formalism of the Dynamic Mode Decomposition technique together with the concept of balanced realization. It is assumed that the only information available on the system comes from input, state, and output trajectories generated by numerical simulations or recorded and estimated during experiments, thus the approach is fully data-driven. The goal is to obtain an input-output low dimensional linear model which approximates the system across its operating range. Since the dynamics of aeroservoelastic systems markedly changes in operation (e.g. due to change in flight speed or altitude), time-varying features are retained in the constructed models. This is achieved by generating a Linear Parameter-Varying representation made of a collection of state-consistent linear time-invariant reduced-order models. The algorithm formulation hinges on the idea of replacing the orthogonal projection onto the Proper Orthogonal Decomposition modes, used in Dynamic Mode Decomposition-based approaches, with a balancing oblique projection constructed entirely from data. As a consequence, the input-output information captured in the lower-dimensional representation is increased compared to other projections onto subspaces of same or lower size. Moreover, a parameter-varying projection is possible while also achieving state-consistency. The validity of the proposed approach is demonstrated on a morphing wing for airborne wind energy applications by comparing the performance against two algorithms recently proposed in the literature. Comparisons cover both prediction accuracy and performance in model predictive control applications.

systems, which can be found in some fluid dynamics problems [6]. Moreover, as recently shown in [7], projections onto POD modes are not uniquely defined, due to the arbitrariness of the definition of the state. These findings reinforced the known fact that the quality of the approximation highly depends on the choice of inner product and thus care is required when the projection operator is computed.
Despite these potential shortcomings, POD-and DMD-based methods have been successfully applied in various aerospace and control flow problems [1,8,9,10,11]. However, a relatively unexplored application domain of data-driven (or equation-free) reduced-order modeling (ROM) is aeroservoelasticity, where the coupling among multiple disciplines (e.g. aerodynamics, structural dynamics) and components of the system (e.g. wing, actuators) often results in high-order models. Expect for a few notable recent exceptions, e.g. the works in [12,13] where nonlinear ROMs have been developed, the standard practice to reduce dimensionality is the use of well established model-based reduction technique [14]. See the work in [15] for an application and further references on this line of research. However, the increasing complexity of the high-fidelity solvers (often made up of distinct sub-solvers for the different disciplines) on one hand, and the potential advantage of recalibrating or directly substituting parts of the code with experimental or flight data on the other, favour the adoption of equationfree strategies. Among the possible reasons for the lack of their application in the field, two important issues are highlighted here.
First, a common feature of the majority of the available approaches is the focus on internal dynamics, meant here as partial or ordinary differential equations without external excitations and with fully observable states. The work in [16] recently extended the DMD framework to controlled systems (DMDc), but the key steps of the algorithm (specifically, the selection of the projecting subspace) do not substantially change. That is, emphasis is not put on preserving the input-output behaviour of the system, which is crucial for control systems.
Second, in aeroservoelastic applications, capturing the variation in the stability and response of the system as the operating conditions change is paramount. This can be done, for example, using the so-called Linear Parameter-Varying (LPV) representation [17], which are of acknowledged benefit for control related tasks [18,19,20,21]. Unfortunately, obtaining accurate models featuring low orders is notoriously a difficult task [22], even for the well explored class of model-based approaches [23,24,25]. One of the most common strategies is to seek low-order linear time-invariant (LTI) representations for frozen-parameter conditions (defining a grid) and then interpolate them for intermediate values of the parameters. Since states, and thus state-space models, are defined up to a nonsingular (similarity) transformation, a correct interpolation requires that the state of each frozen model is defined with respect to the same basis. The need to work with a consistent state-space basis for the local ROMs, required for a correct interpolation [26], poses a challenge for DMD-inspired data-driven approaches. State-consistency will depend indeed on the selection of the projecting subspace. If this changes across the parameter range, as it is the case when one computes the POD modes at each grid point, then state-consistency will not hold in general. Conversely, if the subspace is kept fixed for all the frozen-parameter LTI systems, then accuracy might deteriorate since projection will no longer take place onto the optimal (from an energy point of view) subspace for the considered parameter.
Motivated by the discussion above, the main contribution of this paper is the proposal of a novel equation-free approach to obtain LPV low-order models, namely the Balanced Mode Decomposition (BMD) algorithm. The key idea is to use, instead of an orthogonal projection associated with one subspace (as in standard DMD), an oblique projection, which is associated with two subspaces, namely a basis space and a test space, characterizing the range space and null space of the projection, respectively. Oblique projection, often encountered in model reduction [27] and system identification [28], was also used for modelbased reduction of LPV models in [29]. As detailed in Section 3, the oblique projection proposed here can be interpreted, within the context of DMD-type approaches, as an alternative choice to the subspace spanned by the POD modes, and it is instrumental to achieve two favourable properties. The first is that emphasis can be put on the input-output behaviour of the ROM by defining the range and null spaces of the projection as a function of the controllability and observability Gramians. These objects are well known in the context of model order reduction of linear systems, as they are the main ingredients to perform balanced truncation [14]. This technique consists of transforming the system in balanced coordinates and then removing the states associated with the lowest degree of controllability and observability. In the same spirit, range and null spaces of the proposed oblique projection are defined so that the identified model is (approximately) in balanced coordinates, and thus projection onto lower-dimensional subspaces will preserve the structures in the data matrices that are most observable and controllable. The second favourable property is that the LPV model has a consistent state-space basis [26] across the parameter space without having to sacrifice accuracy. This results from the fact that one subspace (the basis space) is common to all parameters and thus provides a common basis. At the same time, the other subspace (the test space) has no influence on the state's basis, and thus can be chosen different for each parameter in order to alleviate the limitations of a fixed subspace projection.
The second contribution of the paper is to extensively compare the results of the BMD method with two recent extensions of DMD with control. The first algorithm is the algebraic DMDc (aDMDc) [30], which extended DMDc to parameter-varying systems described by algebraic, in addition to differential, equations. Including algebraic constraints is very important, for example, when considering state trajectories generated by aerodynamic solvers capturing unsteady effects, such as in panel methods or unsteady vortex lattice methods [31]. The second algorithm is the input-output reduced-order model (IOROM) approach, proposed in [32] to construct data-driven reduced-order LPV models. Improved ways of defining the low-dimensional subspace such that stateconsistency is achieved while preserving accuracy in the (orthogonal) projection were proposed therein. However, the projection operator is the same for all parameters, and is obtained from the POD modes as in standard DMD. An extension of IOROM to handle algebraic constraints is also developed here in order to allow for a fair comparison. The algorithms are tested on a highfidelity, fluid-structure interaction (FSI) numerical model of an airborne wind energy (AWE) morphing wing. The FSI simulator is described in [33] and the wing was analyzed in detail in [34]. Airborne wind energy and morphing wings are paradigmatic examples of application domains where the system's response originates from complex interactions across different domains, and thus could benefit from equation-free approaches. The first type of comparison investigates the accuracy of the reduced-order models to predict various outputs of the wing as the size of the model is decreased. In a second set of analyses, models featuring different orders are used by a model predictive control (MPC) algorithm to track predefined trajectories of the airborne wind energy system with the goal of gaining insight into the trade-off between size and performance, the latter evaluated by simulating the FSI solver in closed-loop with the MPC controller.
Preliminary results of the work were presented in [35]. Figure 1 shows a conceptual representation of the proposed data-driven ROM framework.

Data-driven reduced-order modeling
This section provides background material on the tools and concepts relevant to the reduced-order modeling algorithm proposed in this work. In Section 2.1 the general data-driven low-order modeling problem is presented. Section 2.2 reviews the algebraic DMD with Control (aDMDc) [30], and Section 2.3 reports on the input-output reduced-order model (IOROM) [32]. These are the two ROM algorithms from the literature used for comparison in this paper.

Problem statement and preliminaries
The starting point is a generic discrete-time nonlinear parameter-varying model which can be used to describe typical control systems, such as aeroservoelastic systems modelled by FSI solvers  where x ∈ R nx , u ∈ R nu , y ∈ R ny are the state, input and output, and ρ : N → R nρ is a vector of time-varying parameters defining the operating conditions of the system. The problem of finding an LPV low-order approximation of (1) can be divided into two phases: first, local LTI approximations for frozen values of ρ in a pre-defined grid {ρ j } ng j=1 are computed; then, an LPV model is obtained through interpolation. The following discussion is concerned with the former phase.
It is assumed that for each frozen valueρ there exists an equilibrium (or trim) point characterized by the tuple (x(ρ),ū(ρ),ȳ(ρ)) such that can then be used as states of an LTI approximation of the system around the where (A(ρ),B(ρ),C(ρ),D(ρ)) is a state-space representation completely describing the linearization about the trim point associated withρ. The dependence of local (i.e. related to LTI approximations for frozen values of ρ) quantities on the parameter ρ will be dropped in the remainder. It is understood that they depend on the particular value of the parameter considered. Instead, the subscript ρ will be used when discussing LPV models.
In the data-driven setting, the only information on the system comes from input, state, and output trajectories {x k , u k−1 , y k−1 } ns k=1 of length n s . After having subtracted from them the corresponding trim values, these trajectories can be used to form the following snapshot matrices The notation [X 0 ; U 0 ] will denote the operation of stacking row-wise two matrices X 0 and U 0 .
The first goal is to obtain a linear time-invariant low-order approximation of (2), that isz wherez ∈ R nz and n z n x . Once this is available, the family of frozen LTI systems, or directly the signals of interest, are interpolated so that the response of the system is available at each value ρ k for a generic time-varying trajectory of the parameter.

Algebraic Dynamic Mode Decomposition with Control algorithm
The algebraic Dynamic Mode Decomposition with Control (aDMDc) algorithm was recently proposed in [30] to extend the DMDc algorithm to systems described by algebraic-differential equations. The DMDc algorithm from [16] is first briefly reviewed. This algorithm seeks a data-driven approximation of the matrices involved in the state equation (2) by means of two truncated singular value decompositions (SVD) of the snapshot matrices. The first one is where the subscript r denotes a truncation of the SVD decomposition of order r (obtained by keeping only the r largest singular values in the decomposition).
Note that the value of r does not define the size of the final reduced-order model, and it could be set for example by using the hard threshold criterion suggested in [36]. The effect of choosing r on the accuracy of the model will be discussed in the result section. The second truncated SVD is computed from the snapshot matrix X 1 where the columns ofÛ nz are also called POD modes of X 1 and are used for the projection onto a lower dimensional space. The selection of n z defines the size of the reduced-order model. The thresholds used in (4) and (5) should be chosen such that r > n z [16].
An approximation of the high-order matrices appearing in (2) can be formulated in terms of the truncated SVD (4) and the snapshot matrix X 1 as Then, a low-order approximation is obtained by projecting (6) onto the set of POD modes by making use of (5) Therefore, the low-order model obtained by DMDc (which only includes the state equation) isz wherez ∈ R nz is the state of the low-order model and the high-order state can be recovered byx =Û nzz .
The aDMDc algorithm [30] builds on the DMDc approach and addresses the presence of algebraic constraints in the dynamic equations which might arise when considering unsteady aerodynamics features. Specifically, the morphing wing analyzed in [30] is described by an FSI solver that implements a 3D panel method with a free evolving wake inspired by the method in [31]. This leads to a dependence of the states' evolution on the inputs at the next time step.
Therefore, a slightly different starting point from the general one presented in (1) has to be considered, namely where g is in general a nonlinear function taking into account the dependence of the states on the control inputs at the next time step. This dependence results from algebraic equations relating the doublet strengths (aerodynamics states) and downwash (function of the other states and the control inputs). This effect is sometimes accounted for with artificial aerodynamic states by simply changing the feedthrough matrix to the outputs. However, to correctly capture the evolution of the states it is important to formulate the problem as stated in (7). The reader is referred to [30] for further discussion on this aspect.
The proposed LTI representation of the system accounting for the algebraic constraints due to the unsteady aerodynamics is where, as in DMDc, the objective is to find a low-order approximation for the state equation only.
The only difference with respect to DMDc is that now the first SVD decomposition is computed with respect to the snapshot matrices X 0 , U 0 , and U 1 , that is And the high-order matrices are thus approximated by A low-order approximation is then obtained by projecting (6) onto the same set of POD modes used in DMDc (5) [F G L] = Û nz AÛ nzÛ nz BÛ nz R .
This procedure results in the aDMDc low-order model where the high-order state can again be obtained fromx =Û nzz .
The approach proposed in the parametrically varying version of the aDMDc algorithm is to use a different set of POD modes for each value of ρ in the grid {ρ j } ng j=1 . The frozen LTI models (8) are then simulated simultaneously, the relative states are lifted to the high-order ones using the corresponding projection matrices (e.g.Û nz (ρ j ) for the model corresponding to the jth element in the parameter space), and the state corresponding to the desired value of ρ is obtained by interpolating the high-dimensional states. A first consequence of this approach is that the frozen LTI models (8) do not have a consistent basis for the state, because the basis of the state is determined by the POD modes, which change across the parameter grid. While this has the advantage of projecting over POD modes specifically computed for a particular value of ρ, it also requires running in parallel all of the low-order models. Moreover, this algorithm does not provide an LPV model and thus the use of LPV robust control design strategies is precluded [21]. While other control techniques, such as model predictive control, can still be successfully used (see Section 4.4), the necessity to run in parallel, multiple low-order models, is a drawback of the method when targeting real-time applications.

Input-output Reduced-Order Model algorithm
The input-output reduced-order model (IOROM) algorithm was proposed in [32] to compute a family of state-consistent data-driven low order LTI statespace models (including the output equation) which can be directly parameterized by the vector ρ.
Consider first the case when there is no parameter dependence (or equivalently, ρ is fixed). Drawing inspiration from the interpretation of DMD as linear dynamics fitting [4], the main idea is that, given the snapshot matrices (3), the matrices (A, B, C, D) defining (2) can be obtained by solving the following least-squares problem where the subscript F denotes the Frobenius norm of a matrix. Without appropriate regularization, this problem would be ill-posed for high-dimensional systems (n x 1). Most importantly, even if (9) was solved accurately, it would not provide a low dimensional representation of the system. For these reasons, an orthogonal projection of the state onto a low dimensional subspace of dimension n z is performed by introducing the projection matrix Q ∈ R nx×nz , where Q Q = I nz , such that the orthogonal projection ofx on an n z -dimensional subspace is given by QQ x. Equivalently, one can think that the original state is approximated byx ∼ = Qz for some reduced-order state (or coefficient vector) z ∈ R nz . This results in the following low-order state-space model The vectorz = Q x ∈ R nz can thus be interpreted as the state of the low-rank approximation of (2) The projection matrix Q is constructed from the POD modes of X 0 , that is The least-squares problem giving (F ,G,H,D) is then whose solution is where † denotes the pseudo-inverse of a matrix. It is worth noting that the reduced-order model given by the IOROM algorithm is qualitatively similar to the one associated with DMDc. The main difference (besides the output equation, not considered in DMDc) is that the pseudo-inverse operation, which also amounts to an SVD decomposition and thus is conceptually similar to (4), is done here directly on the projected snapshot matrices. This is different than what is done in DMDc, where the SVD decomposition (4) is applied to X 0 and U 0 . A minor difference is also that the POD modes are computed here with respect to X 0 instead of X 1 .
In the parameter-varying case, the regression problem (13) is solved at each value of the parameter grid {ρ j } ng j=1 by taking the corresponding snapshot ma- . By always using the same projection matrix Q when computing the low-order models at different ρ, state-consistency is automatically guaranteed because the orthogonal projection has the same range space. An LPV reduced-order model is then obtained by interpolating (13) across the parameter's range. That is are obtained by interpolating the corresponding matrices for the value of ρ at timestep k. Note that the low-order state isz k = ) is added to correctly take into account this effect [32].
Since the choice of a fixed projection matrix is typically associated with less accuracy, two strategies are proposed in [32] to alleviate this issue. The first consists of using in the decomposition (11) a fat snapshot matrix X 0 obtained by stacking column-wise the snapshot matrices of multiple parameters. This matrix will feature n s n g columns, which can result in computationally expensive calculations when this number is large. The second, less accurate but more practical in case several grid points are analyzed, consists of iteratively building Q by incrementally processing the snapshot data from each grid point in a similar fashion to the Gram-Schmidt orthogonalization procedure. The former strategy is used here when showing results for the IOROM algorithm, together with a linear interpolation of the state-space matrices.

Balanced Mode Decomposition with oblique projection algorithm
This section presents the technical aspects of the Balanced Mode Decomposition (BMD) algorithm proposed in this paper. Section 3.1 clarifies the goals and the novelty of the contribution with respect to previous works. Section 3.2 presents the algorithm and Section 3.3 details the computation of the subspaces defining the oblique projection. Finally, Section 3.4 presents a version of the algorithm which can handle algebraic constraints and thus allows the analyses in Section 4 of the morphing wing with an unsteady aerodynamics model.

Novelty and connections with prior work
The main motivation for the proposal of the BMD algorithm for data-driven LPV low-order modeling is to address two limitations of recent extensions of the celebrated DMD method to input-output parameter-varying models. The first one concerns the use of Q (i.e. the subspace spanned by the most energetic POD modes according to the standard choice of inner product as unweighted scalar product) for the projection of the higher-order dynamics, which is suboptimal as also acknowledged by the authors of [32]. In the input-output context, a subspace typically providing lower input-output errors with respect to the others having same size n z is the one where the system's state is in balanced coordinates [37]. This is indeed the rationale behind balanced truncation, which consists of removing the states corresponding to the smallest n x − n z Hankel singular values [38]. The justification for this is that the sum of the Hankel singular values provides a lower bound, and for systems in balanced coordinates, an upper bound on the error of the approximation achieved by removing system's states. Even though not guaranteed to be optimal, balanced truncation is a very effective tool in model-based order reduction [14,39]. These ideas are used here to propose a new projection operator for the high-dimensional state.
Whereas the aspect mentioned above is relevant also for the case where the sought model is an LTI (i.e. when there is no parameter dependence), the second one is specific to the LPV setting. Precisely, the second limitation addressed by the BMD algorithm concerns how to handle state-consistency across the frozen models in order to estimate the system's response at intermediate points in the parameter grid. In the currently available approaches, this is addressed in two possible ways. When state-consistency is not fulfilled, all reduced-order models are run in parallel by interpolating directly the high-dimensional lifted state. This is the case of aDMDc, and while it has the advantage that the projection operators are parameter-dependent (i.e. at each parameter's value one can use a different set of POD modes), an LPV model is not available and moreover computational efficiency might be compromised. On the other hand, a parameter-independent projection matrix for all frozen models can be used in order to guarantee state-consistency. This is the case for the IOROM algorithm, and it has the drawback that, whereas the orthonormal basis associated with the n z most energetic modes will be in general different at each value of ρ, a fixed one that is common to all parameters is used.
The central idea to overcome both of the aforementioned issues is to replace the orthogonal projection employed in standard POD-based approaches by an oblique projection. Given V ∈ R nx×nz and W ∈ R nx×nz , such that W is bi-orthogonal to V , i.e. W V = I nz , the oblique projection ofx is given by Πx, where Π = V W . As a result, the original state is approximated asx ∼ = Vz (where, as before,z ∈ R nz is the reduced-order state), and the component ofx that is eliminated by the projection is in the nullspace of Π.
As opposed to the orthogonal projection, which is characterized by a single subspace (the one spanned by the columns of Q), the oblique projection is defined by two subspaces: the basis space (spanned by the columns of V ), such that the projection ofx lies in the span of V ; and the test space (spanned by the columns of W ), such that the projection Vz has zero error within it, i.e.
W (x − Vz) = 0. Technically, the high-dimensional state vector is projected along the orthogonal complement of the subspace spanned by the columns of W onto a subspace spanned by the columns of V . In practice, this means that what is lost by projectingx (i.e. the nullspace of the projection) is orthogonal to W , and the state basis only depends on V . The two issues discussed above are then addressed by: computing V and W from the empirical controllability and observability Gramians of the system (which leads to a model-free balanced truncation); employing a fixed V and a parameter-dependent W . Since V by definition defines the basis of the vector space where the state of each model is defined, this basis will be common to all the local state-space models, and thus state-consistency is guaranteed.
The idea of using an oblique projection for LPV model-order reduction was first proposed in [29]. Therein, the setting where a model of the system is available (in the form of high-order state-space models) is considered, and thus both the construction of V and W , and the computation of the low-order model, is model-based. In the data-driven ROM literature, balancing concepts are used in two important techniques, namely Balanced POD (BPOD) [40] and the Eigensystem Realization Algorithm (ERA) [41]. The former is only partially equation-free: the controllability Gramian is computed from data, while for the observability Gramian an adjoint simulation model is needed. Additionally, the high-order state-space matrices are required for the balanced projection. For the case of ERA, a balanced model comes from impulse response simulations of the model in the spirit of system identification algorithms from realization theory [42]. The ERA algorithm is closely related to BPOD, as it can be interpreted as a data-driven balanced truncation. An important difference is that ERA provides only the reduced-order model and not the balancing transformation, namely the set of vectors known as balancing and adjoint modes in BPOD. These modes are the counterpart of the basis and test space in BMD, respectively, and are a desirable output of a ROM algorithm as they show the most important spatiotemporal structures in the dynamics. In an aeroservoelastic setting, this can provide insights into efficient design solutions. It is recalled that BPOD can be interpreted as a special case of POD when impulse responses are used to build the snaphost matrices and the observability Gramian is used as inner product [40]. Conceptually (because in practice the algorithm formulation is articulated in a different way), it can be helpful to think of BMD as a version of DMD that makes use of this special case of POD.

BMD regression problem
We consider first the frozen-parameter case and, by virtue of the previously discussed oblique projection, propose the following low-order LTI system model where the computation of the balancing basis V and test spaces W from system's trajectories will be detailed in Section 3.3. The vectorz = W x ∈ R nz can thus be interpreted as the state associated with the following low-rank approximation of (2) The matrices (F ,G,H,D) can then be obtained with the following least-squares which has solution To build a low-order LPV model, snapshot matrices are first collected for the values of the parameter in the grid {ρ j } ng j=1 , and the least-squares problem (16) is solved at each grid point. Crucially, the test space W is allowed to be a function of ρ. This leads to the following solution for the reduced-order models in the grid where the dependence on ρ of the local quantities has been here explicitly reported in order to clearly point out that all the objects involved in the construction of the low-order state-space models are parameter-varying.
The BMD LPV reduced-order model is then obtained by interpolating the frozen matrices (18) across the parameter's rangẽ where, as in IOROM, (F ρ k ,G ρ k ,H ρ k ,D ρ k ) are obtained by interpolating the corresponding matrices for the value of ρ at timestep k, and the term (z(ρ k ) − z(ρ k+1 )) takes into account the fact that the equilibrium point associated with each ρ is in general different, andz(ρ k ) = W x(ρ k ). Note that, since V is fixed, the basis space is common to all the frozen models and thus the interpolation can be done at the state-matrices level (as in IOROM). However, the projection is parameter-dependent due to the use of a parameter-varying test space W (ρ).

Basis and test spaces construction
In order to preserve the most important features of the input-output mapping of the system when projecting into the lower order subspace of dimension n z , the matrices V and W are computed from the controllability and observability Gramians, respectively W c and W o . This ensures that the projection preserves the most observable and controllable states, enabling an approximate data-driven balanced truncation of the reduced-order LPV model.
Because the approach is fully data-driven, empirical Gramians are computed from data matrices consisting of appropriate state trajectories using known systems theoretical results. The empirical controllability Gramian can be obtained, following the definition [37,43], by impulse response simulations (one for each input channel). As for the empirical observability Gramian, if the model is linear and its adjoint is available, then it can be computed from impulse response simulations (one for each output channel) of the adjoint system, as done in balanced POD [40]. This computation is identical to the one giving the controllability Gramian, but is applied to the adjoint system (for an LTI model in state-space form, this is obtained by replacing A and B by A and C ). If the above does not hold, for example in case one has only access to the system's trajectories and not to the model's matrices, the approach developed in [44], valid also for nonlinear systems and not requiring an adjoint model, can be used. In this method the data matrices used for the Gramian computation consist of state trajectories obtained from unforced (zero input) simulations (one for each state) obtained by perturbing the initial condition of each state. Since these are unforced responses, when the system is sufficiently damped, it will be generally sufficient to observe only the initial time-steps and thus this calculation can be parallelized and efficiently implemented to reduce the computational time.
Once W c and W o are available, a procedure based on [29, Proposition 2] is employed to compute the test and basis spaces. This construction is reported in the first part of the pseudocode given in Algorithm 1, which summarizes input, output, and main steps of the BMD algorithm (MATLAB notation is used for matrix operations and rows/columns selection). For a fixed value of ρ, the construction proposed in [29] is an equivalent procedure, but more numerically stable for large-scale systems, to the well known square root algorithm for balanced truncation [43]. Indeed, it can be noted (see lines 5-11) that the subspace

Output:
test space projection matrices {W (ρ j )} ng j=1 ; fixed basis space projection matrix V ; reduced-order models at the grid points Cholesky factorization of W c 3: R =R(1 : n z , :) 16:

18: end for
The output {F (ρ j ),G(ρ j ),H(ρ j ),D(ρ j )} ng j=1 provided by the BMD algorithm is a grid LPV model. After an interpolation algorithm to evaluate the matrices' entries for any value of ρ inside the considered range has been chosen, this model can be used for simulation and control design. Note also that recently proposed robust analysis methods for linear-time varying (LTV) systems [46,47] can be applied to this model, e.g. to investigate specific aircraft manoeuvres. Indeed, by fixing a particular trajectory for ρ the LPV system is transformed into an LTV one. Moreover, the parameter-varying test space W (ρ j ) can be useful to gain insights into the aeroservoelastic modes which have been eliminated and those that have been kept in the projection, while the parameter-independent basis space can be used to recover at each time-step k the high-dimensional state via the transformationx k = Vz k .
As noted in the introduction, the algorithm provides an approximate balanced truncation. Approximation is related to the use of empirical Gramians, which are only finite-time approximations of the true ones (for this reason, also called finite-time Gramians) since their computation is trajectory-based.
As a result, they only provide in principle a finite-time balanced realization [14], whereas the theoretical order reduction error bounds are only available for infinite-time balanced realizations. This source of error can however be made arbitrarily small by using long enough data sequences for constructing the Gramians. The slowest decay rate of the system's impulse responses is the key parameter to consider when choosing the length of the trajectory [48].

Extension to handle algebraic constraints
Since the BMD algorithm will be applied in Section 4 to a wing described by an FSI solver which implements the algebraic constraints described in Section 2.2, an extension to handle this instance is presented here. For a fixed value of ρ, the model structure for the high-order model becomes where a potential effect of the algebraic constraints in the output equation is also considered via the matrix P (in the analyses of the morphing wing this matrix was, as expected, always found to be zero). Therefore, the low-order approximation (for each value of ρ in the parameter grid) becomes The new objective function to be minimized is

and the new optimal solution is
The BMD LPV reduced-order model with algebraic constraints is then, in analogy to (19), given bỹ z k+1 = F ρ kz k + G ρ kũ k + L ρ kũ k+1 + (z(ρ k ) −z(ρ k+1 )), The IOROM algorithm has also been extended to the algebraic-differential case in order to allow its application to the test case considered in Section 4.
This can be done in a similar way to what has been shown above for BMD.
Specifically, starting from (12), the new least-squares problem for the IOROM reduced-order model is which has the following solution Eq. (20), when the interpolated matrices are taken from (21), provides the IOROM LPV reduced-order model with algebraic constraints.

Results
This section presents and discusses results obtained by applying the three algorithms aDMDc, IOROM, and BMD to the flexible and highly cambered morphing wing depicted in Fig. 1. The wing is made of composite material, and the trailing edges are able to morph and, by doing so, to increase or decrease the camber, thus replacing conventional ailerons. The reader is referred to the related previous works for details on the wing design [49] and its investigation with fluid-structure interaction tools [34].

Summary of the wing's FSI model
The high-fidelity FSI model of the morphing wing is presented in detail in [33,50]. The structural model is based on a combination of linear plate and beam elements. The external skin of the lifting surface and the Voronoi-based internal structure are modelled using plates, while the stringers are modelled with beam elements. The stiffness and mass matrices are obtained with the commercial software Nastran [51]. From these, a modal decomposition is performed to extract the structural modes of the wing, which are coupled via thin plate spline and inverse distance weighting [52] with the aerodynamic model.
The aerodynamic model is based on a 3D unsteady panel method [53]. The state of the system x consists of the total number of structural modes of the wing and the doublet strengths (from the 3D panel method solver), with n x =618. The input vector u of size n u =6 is given by where α is the angle of attack, p, q, and r are the roll, pitch, and yaw rotation rates, and F s and F as are the (normalized) symmetric and anti-symmetric morphing actuation inputs. Their value is associated with a camber deformation and is thus related to a trailing edge deflection: specifically, the amount of upwards (negative value) or downwards (positive value) deflection.
As for the output channels, we will consider both single output and multiple output models. Emphasis will be given to the case where the output is the first bending mode of the wing (n y =1), since this is usually the one associated with dynamic instabilities and large deformations, and thus it is of particular interest for active control tasks [54,55]. The flight speed will be considered as the time-varying parameter (n ρ =1).
The training phase, common to both the algorithms and consisting of generating the snapshot matrices in Eq. (3), is carried out by exciting the system with a series of impulses deployed in random order in all input channels, following the same procedure adopted in [30]. The amplitude of these signals has been chosen so that nonlinear effects due to the wake's evolution are not excited. The same precaution is used for the computation of the empirical gramians. Trajectories are of length n s =500 and are recorded with sampling time 0.006 s.

Fixed-parameter models
In this first set of tests, the accuracy of the different models at fixed values of the flight speed V is assessed by means of sinusoidal inputs. The frequencies, different for each channel, are expressed in terms of the aircraft reduced frequency f r =: V c , wherec=0.29 m is the mean chord of the wing. The first, third, and fifth input channels are excited with sinusoids 1 5π f r , 1 10π f r , and 1 20π f r , respectively, while the other channels are set to zero. For reference, the first bending mode of the wing has a frequency of approximately 10 Hz [33]. This  It is clear that for all speeds BMD provides the smallest error for a very low-order approximation of the full dynamics, as expected in view of the choice of low-dimensional subspace where the high-dimensional data are projected. As the size n z of the system increases, the difference between the algorithms is less noticeable and, for high enough orders, the algorithms tend to give same results.
It is also noted that the results obtained with the aDMDc algorithm showed great sensitivity, in the range of n z displayed in Fig. 2, to the SVD truncation order r employed in Eq. (4). Using the hard threshold criterion from [36] did not provide good results as it resulted in a very large r (therefore the truncation included very low singular values deteriorating the approximation). Since finetuning the value of r to optimize the results obtained with aDMCc would have required trying several values of r for each different value of n z , this was not pursued here. Instead, the heuristic choice r=n z +10 was implemented and proved to provide reasonable results. Even though, because of this possibly suboptimal choice, the resulting gap in performance with the other two methods observed in Fig. 2 can be also ascribed to numerical inaccuracies associated with the decomposition (4), the need to optimally choose r can be considered as a disadvantage of the aDMDc problem formulation. It is finally observed that, in the extensive analyses performed, aDMCc typically showed worse performance than BMD irrespective of the choice of r.

Parameter-varying models
In the second set of tests, the accuracy during parameter-varying manoeuvres is tested. A manoeuvre of 3 s where the flight speed linearly increases from V =20 m/s to V =50 m/s is analyzed. Unless otherwise specified, the reducedorder models are obtained using snapshot matrices obtained gridding the flight speed range every 2 m/s and thus using 16 different speeds (n g =16). A linear interpolation will be used to evaluate quantities for values of ρ that are not in the grid.

Sinusoidal excitation
The same sinusoidal input signals used in Section 4.2 are considered here.
In Fig. 3, the bending mode amplitude response obtained with the FSI solver (FSI ) is compared with the predictions of the three algorithms when the order of the models is fixed at n z =14. The goal is to compare the prediction accuracy of the algorithms when a very low number of states (in comparison with the order of the original system) is employed. It is noted that the same observations can be gathered when values of n z in the same range are considered. All the signals are normalized by the largest value of the bending amplitude measured in the FSI simulation, which for the analyzed case was 0.7. Note that a unitary value of the first bending mode corresponds to a wingtip displacement of 4.6 cm. In this case, since only one value of n z was considered, the aDMDc model was here obtained by fine tuning the threshold value r in order to provide the best results. The plot confirms, also in the LPV setting, that the BMD algorithm guarantees the smallest error when a low-rank approximation of the system is desired.
In this simulation, aDMDc outperforms IOROM, possibly due to the fact that it uses a parameter-varying set of POD modes. However, aDMDc does not provide a family of interpolated low-order models, and interpolates directly the highorder states, thus requiring parallel simulations of the low-order models. The better performance of BMD, despite the fact that a part of the projection (the one related to the basis space) is constant, is ascribed to the improved selection of subspace for the projection compared to the standard POD one, common to the other two methods. In addition to the improvement in the accuracy, the BMD algorithm is also capable of providing, like the IOROM algorithm, a fam-ily of consistent LTI models with the advantages for LPV control design and in general real-time applications.

Effect of the input signals
This section investigates the accuracy of the reduced-order models for different types of input signals. The Euclidean norm of the error signal between the first bending mode amplitude provided by the high-fidelity FSI and the prediction obtained with each of the three ROM algorithms is again used as metric to assess the quality of the approximation. Three classes of inputs are considered: Sine coincides with the signal tested so far and already investigated in [30]; Chirp excites the system by injecting in all 6 input channels defined in (22) a chirp signal with frequency linearly varying from 1 50π f r to 2 15π f r ; PRBS excites the system by injecting in all 6 input channels a PRBS-9 sequence. This last input, namely a Pseudo-Random Binary Signal (PRBS) [56], is a deterministic signal with white-noise-like properties. It is very well known in the system identification and experiment design fields since it has the favourable property of equally distributing energy across all the frequencies in the input spectrum.
In this way, information on the models in different frequency ranges can be extracted. Although not a common input in aeroservoelastic applications, it has been used in this spirit here, since the previously adopted sets of input will only give information on the behaviour of the reduced models around the aircraft reduced frequency f r . Results are shown in Fig. 4.
The plots confirm the advantage in using the BMD approach when seeking a low-order model capturing parameter variations. These results are important considering that they are obtained by exploring different frequency ranges of the system's response.
An interesting aspect observed in Fig. 4 is that none of the algorithms exhibit a monotonic improvement of the model's accuracy (measured here by the relative prediction's error) as a function of the systems order. Indeed, in a very few cases, a small deterioration can be observed between two consecutive values of n z , before the curve keeps decreasing as n z increases. Even though this might seem surprising at first glance, there is no theoretical guarantee that such a monotonic improvement is achieved in this type of algorithms. The reason is that adding one mode to the low dimensional subspace where the dynamics is described might (in principle) deteriorate the model's approximation, if that mode alone does not add meaningful information. This is for example the case when pairs of modes describe relevant features of the system (e.g. modes of vibration), and so only when both of them are used in the projection there is an improvement. In practice, numerical errors due to the SVD truncation and the interpolation of the models across the parameter's grid can also have an effect on these results. However, it is noted that the number of occurrences and the entity of this deterioration are limited.

Effect of the parameter grid
In this section we analyze the effect of the flight speed grid where the reducedorder LTI models are computed, i.e. the selection of the parameter n g . This is an important aspect, which is known to influence both the accuracy of the LPV models and the quality of the control design based on them. Three cases are compared in Fig. 5: n g =4 where the grid includes one plant every 10 m/s; n g =8 where the grid includes one plant every 4 m/s from V =20 m/s to V =48 m/s and then V =50 m/s; n g =16 which is the grid used so far. PRBS inputs are used to excite the system in these tests. From the analyses it can be gathered that aDMDc is more robust than the other algorithms to the value of n g . In particular, both IOROM and BMD present poor performance for a few reduced order models in the range of n z between 30 and 40 when the flight speed grid is coarser. The reason for this behaviour is due to the interpolation approaches employed by the three ROM schemes. Whereas IOROM and BMD interpolate the low-order state-space matrices obtained at the grid points, aDMDc interpolates directly the high-order vector states which are obtained by lifting the low-order statesz from the local models (8) running in parallel. Interpolating every entry of the state-space ma-trices therefore makes the choice of the grid a more delicate aspect in IOROM and BMD. The unstable behaviours resulting in very high (out of the plot) errors are indeed ascribed to numerical inaccuracies in this interpolation. It has been observed that the entries of the matrices are overall bigger as the order n z is increased, hence justifying why these outliers take place within the aforementioned range of model's orders. While there does not seem to be a fundamental reason to explain it, it is apparent that the sensitivity of BMD to the coarseness of the parameter grid is more accentuated. A possible explanation is that, because at each grid point the computation of the empirical Gramians is required, dealing with coarse grid exacerbates the numerical inaccuracies associated with the projections on low-order models. Improved interpolation schemes, not considered in this work, could be employed to ameliorate this issue. Except for these isolated numerical problems, BMD shows better performance even when very coarse grids are employed.

Prediction capability for other signals
The capability of the models to predict other quantities of interest, such as for example aerodynamic coefficients depending on the system's states, is investigated. In particular, we test the accuracy when these coefficients are added to the vector of output (this is done for BMD) or computed directly from the states (this is done for IOROM and aDMDc). In the latter case, the loworder states are lifted to the high-order ones, which are then used to compute the coefficients using their known relationship to the states. While this is the only possible way of reconstructing the system's signals for aDMDc, in IOROM and BMD this can alternatively be done by simply adding the desired quantities to the vector of outputs. This would probably be the preferred approach if the signals are used for control (either because they represent measurements fed to the controller or because they are performance measures to be optimized). The different choice done here for BMD and IOROM is for the sake of exploring different models, and results showed that whether the signals were computed from output channels or retrieved from the states had a very minor impact on the predictive accuracy.  The same observations gathered earlier with respect to the trajectory of the bending mode are confirmed here. It is particularly interesting to observe that, even though these coefficients are not outputs of the model, and thus the balancing projection is not aimed directly at capturing them, the BMD algorithm is still able to perform better than the others. Figure 7 shows the same analyses when chirp signals are used as input to excite the model. Similar conclusions can be drawn.

Reduced-order models for Model Predictive Control
In this section, a control application of the morphing wing's low-order models is investigated. Specifically, model predictive control (MPC) [57] is considered, given its well established use in the AWE field [58]. Two distinct reference tracking problems are examined, where pre-defined lift and first bending mode amplitude profiles are tracked while flying trajectories over a range of different flight speeds and in the presence of turbulence and gusts. The analysis of these manoeuvres is motivated by the interest in using active control to guarantee a safe operation for the AWE system (with respect to some of its critical components such as the wing or tether) by keeping indicators of the structural integrity close to desired, and possibly pre-optimized, values. This can avoid passive remedies such as reducing the load transmitted to the ground station, which in turn decreases the amount of wind energy harvested. Having effective and reliable control laws to guarantee the integrity of the AWE system can represent an important enabler for this technology [59].
In its basic form, model predictive control repeatedly solves a finite-horizon optimal control problem of length N c subject to input and state-constraints. At each instant, a model of the system is employed to predict its response and thus select the control sequence (u i ) Nc−1 i=0 which minimizes the cost where r is the reference trajectory, and for a vector x, we denote by x P the weighted l 2 -norm (x T P x) 1 2 . Besides the terms penalizing deviation of the outputỹ from r (with the weighting matrix N ∈ R ny×ny ) and control effort (with the weighting matrix M ∈ R nu×nu ), the cost in (23) also penalizes fast changes in the input via the term ∆ũ k =ũ k −ũ k−1 (e.g. to avoid actuator rate saturation).
The following optimization problem will be solved to obtain the optimal input sequence where: the cost function (24a) is defined in (23); the constraint set U (24c) enforces minimum and maximum values for the input; and (24b) enforces the dynamic constraint that relates the sequence of input to the output via an inputoutput model of the system f . Precisely, f will be formulated here by using the reduced-order aDMDc, IOROM, and BMD models. The goal of the analyses is to compare the associated closed-loop cost J CL MPC , that is the cost (23) incurred by the true system (simulated here by the high-dimensional FSI solver) when this is regulated by the inputs optimized solving problem (24). Since the model is used to predict the system's output, and thus select the input sequence, any mismatch between model and system can result in a degradation of the controller performance.
The analysis considers the case where the morphing wing (Fig. 1, I) flies a trajectory with flight velocities ranging between V =27 m/s and V =50 m/s ( Fig. 8-left). Additionally, the wing is subject to a gust at the maximum flight speed with gust length 0.5 s corresponding to a 1 degree deflection of the angle of attack α, and to turbulence generated with a Dryden filter (Fig. 8-left).
The right plot in Fig. 8   The scenario considered here has only a symmetric morphing actuation input, i.e.ũ = F s . This normalized input is constrained to be in the interval [-3, 3] at each time-step. This is associated with allowable deflections of the trailing edge in the range ± 9 mm. Problem (24) is solved using the MATLAB implementation provided in [60], where the application of MPC with models obtained via DMDc was investigated.
The outputỹ is either the first bending mode, or the lift force generated by the wing, depending on the case. The control horizon is N c =10 and the weights used in the cost are: N =1300, M =10, and M ∆ =0.1 (lift tracking) and N =13000, M =10, and M ∆ =0.1 (bending tracking). The penalty on the output deviation is increased in the latter case due to the difference in magnitude of the two tracked quantities (recall Fig. 8-right). Figure 9 shows the comparison of the closed-loop cost J CL MPC resulting from closing the true plant (simulated by means of the high-dimensional FSI solver) with the MPC controller generated using the low-order models. For the sake of clarity, the closed-loop cost J CL MPC has been normalized in each case by dividing it by the corresponding value obtained with the BMD algorithm when n z =40.
The observations gathered in the previous sections regarding the better prediction performance achieved with the BMD algorithm when low-order approximations are considered are confirmed here in the context of control applications.
Both plots show that, while for higher orders the closed-loop costs have very similar values, when the size of the model is decreased the BMD gives in general the lowest cost. Another interesting observation is that the lift tracking problem is quite robust to the use of low-order models. Indeed, the closed-loop costs are always within two times of the lowest cost achieved at n z =40 except for the case of aDMDc at n z =10). On the other hand, the bending tracking problem is shown to be more challenging when low-order representations are employed.
Whereas no attempt to further optimize the MPC problem tuning was made (all the design parameters were kept the same independent of n z ), this motivates further work on the use of low-order models for control of coupled flexible structures like those encountered in AWE applications. In the real-time control setting, it is important to stress that by using BMD and IOROM models an order of magnitude computational speed-up was achieved with respect to the cases where aDMDc models were used. This is because the aDMDc models have the requirement of running several models in parallel.

Conclusion
The paper proposes the Balanced Mode Decomposition with oblique projection algorithm, a novel data-driven algorithm for constructing low-order LPV models from system's trajectories. Two recent algorithms from the literature, aDMDc and IOROM, are considered for comparison since they both have connections with the newly proposed approach. Technical details on the BMD algorithm are given in order to clearly point out the innovations, and the advantages with respect to previous work. The performance of the BMD algo-rithm is assessed on a morphing wing for airborne wind energy applications.
The results, proposed both for the fixed parameter and, more extensively, for the parameter-varying case, confirm the theoretical advantages discussed in the technical part of the paper. When seeking low-order model representations, the BMD approach achieves generally, among the tested algorithms, the lowest prediction error and best control performance when used as model for an on-line MPC scheme. The improved accuracy is ascribed to the use of a projecting subspace that balances the low-order states (this element is of interest also in a fixed-parameter setting), and to the use of a parameter-varying projection operator (which can thus be enriched with parameter-dependent features, instead of being fixed throughout the range). This has the advantageous feature of being achieved while guaranteeing state-consistency. Owing to these appealing features, it is envisaged the application of BMD for tasks such as off-line and real-time control design, and in multi-disciplinary optimization tool chains, where typically low-order representations are employed as surrogate models.