The Balanced Mode Decomposition Algorithm for Data-Driven LPV Low-Order Models of Aeroservoelastic Systems
TThe Balanced Mode Decomposition Algorithm forData-Driven LPV Low-Order Modelsof Aeroservoelastic Systems
Andrea Iannelli ∗ ,a , Urban Fasel b , Roy S. Smith a a Automatic Control Lab, Department of Electrical Engineering, ETH Z¨urich, Switzerland b Department of Mechanical Engineering, University of Washington, Seattle, USA
Abstract
A novel approach to reduced-order modeling of high-dimensional time varyingsystems, such as those encountered in aeroservoelasticity, is proposed. It lever-ages the formalism of the Dynamic Mode Decomposition technique togetherwith the concept of balanced realization. It is assumed that the only informationavailable on the system comes from input, state, and output trajectories gener-ated by numerical simulations or recorded and estimated during experiments,thus the approach is fully data-driven. The goal is to obtain an input-outputlow dimensional linear model which approximates the system across its operat-ing range. Since the dynamics of aeroservoelastic systems markedly changes inoperation (e.g. due to change in flight speed or altitude), time-varying featuresare retained in the constructed models. This is achieved by generating a LinearParameter-Varying representation made of a collection of state-consistent lineartime-invariant reduced-order models. The algorithm formulation hinges on theidea of replacing the orthogonal projection onto the Proper Orthogonal Decom-position modes, used in Dynamic Mode Decomposition-based approaches, witha balancing oblique projection constructed entirely from data. As a consequence,the input-output information captured in the lower-dimensional representationis increased compared to other projections onto subspaces of same or lower size. ∗ Corresponding author
Email address: [email protected] (Andrea Iannelli) a r X i v : . [ ee ss . S Y ] F e b oreover, a parameter-varying projection is possible while also achieving state-consistency. Therefore, it is envisaged its application for tasks such as off-lineand real-time control design, and in multi-disciplinary optimization tool chains,where typically low-order representations are employed as surrogate models.The validity of the proposed approach is demonstrated on a morphing wingfor airborne wind energy applications by comparing the performance againsttwo algorithms recently proposed in the literature. Comparisons cover bothprediction accuracy and performance in model predictive control applications. Key words:
Reduced-order modeling, aeroservoelasiticity, data-driven,balanced reduction, control systems.
1. Introduction
Data-driven approaches to extract from trajectories of high-dimensional sys-tems, parsimonious models capable of balancing accuracy of the prediction withcomplexity, are an increasingly popular research topic [1]. In fact, pioneering ante litteram contributions to the field, prompted by the goal of identifying low-order structures in complex physical problems such as turbulence, were made inthe fluid mechanics and aerodynamics communities [2]. The fundamental ideacommon to many successful approaches, developed in the wake of these earlycontributions, is to project the high-dimensional data on a lower dimensionalsubspace (also constructed from data), such that the most important features ofthe dynamics are therein preserved. A celebrated example is the Dynamic ModeDecomposition (DMD) approach [3, 4], whereby the spectrum of a low-order lin-ear dynamical model approximating the training data is obtained by leveragingthe Proper Orthogonal Decomposition (POD) [5] reduction technique. Specifi-cally, the projecting subspace provided by POD is spanned by the left singularvectors associated with the largest singular values of a data matrix gatheredfrom observations of the dynamics. The exact interpretation of the largest sin-gular values depend on the inner product used to define the data matrix. Instandard applications, where the so-called snapshot matrix, corresponding to2he correlation matrix between the dynamical states, is used, the largest singu-lar values are associated with the modes capturing most of the energy in thesystem. Thus, the projection onto the lower dimensional subspace preserves thespatial structures with the highest energy content. This criterion for choosingthe projection subspace might not always give the best results, as low-energyfeatures can have a large effect on the dynamics, e.g. in the case of non-normalsystems, 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, and thus require particularcare when the associated projection operator is computed.Despite these potential shortcomings, DMD methods have been successfullyapplied in various aerospace problems [1, 8, 9]. However, a relatively unexploredapplication 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. The standard practice to reducedimensionality is the use of well established model-based reduction technique[10], see e.g. [11] for an application. However, the increasing complexity ofthe high-fidelity solvers (often made up of distinct sub-solvers for the differentdisciplines) on one hand, and the potential advantage of recalibrating or directlysubstituting parts of the code with experimental or flight data on the other,favour the adoption of equation-free strategies. Among the possible reasons forthe lack of their application in the field, two important issues are highlightedhere.First, a common feature is the focus on internal dynamics, meant here aspartial or ordinary differential equations without external excitations and withfully observable states. The work in [12] recently extended the DMD frameworkto 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. 3econd, in aeroservoelastic applications, capturing the variation in the stabil-ity 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 [13], which are of acknowledged benefit for control relatedtasks [14, 15, 16, 17]. Unfortunately, obtaining accurate models featuring loworders is notoriously a difficult task [18], even for the well explored class ofmodel-based approaches [19, 20, 21]. One of the most common strategies is toseek low-order linear time-invariant (LTI) representations for frozen-parameterconditions (defining a grid) and then interpolate them for varying parameters.The need to work with a consistent state-space basis for the local ROMs, re-quired for a correct interpolation [22], poses a challenge for DMD-inspired data-driven approaches. State-consistency will depend indeed on the selection of theprojecting subspace. If this changes across the parameter range, as it is the casewhen one computes the POD modes at each grid point, then state-consistencywill not hold in general. Conversely, if the subspace is kept fixed for all thefrozen-parameter LTI systems, then accuracy might deteriorate since projec-tion 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 isthe proposal of a novel equation-free approach to obtain LPV low-order models,namely the Balanced Mode Decomposition (BMD) algorithm. The key ideais to use, instead of an orthogonal projection associated with one subspace(as in standard DMD), an oblique projection, which is associated with twosubspaces, namely a basis space and a test space, characterizing the range spaceand null space of the projection, respectively. Oblique projection, a strategyoften encountered in model reduction [23] and system identification [24], was alsoused for model-based reduction of LPV models in [25]. As detailed in Section3, the oblique projection proposed here can be interpreted, within the contextof DMD-type approaches, as an alternative choice to the subspace spanned bythe POD modes. It achieves two favourable properties. First, emphasis can beput on the input-output behaviour of the ROM by defining the range and null4paces of the projection as a function of the controllability and observabilityGramians. As a result, the projection will preserve the spatial structures in thedata matrices that are at the most observable and controllable, rather thosehaving higher energy content. Second, the LPV model has a consistent state-space basis across the parameter space without having to sacrifice accuracy.Indeed, one subspace (the basis space) is common to all parameters and providesthe sought after common basis, while the other subspace (the test space) isdifferent for each parameter and thus alleviates the limitation of a fixed subspaceprojection.The second contribution of the paper is to extensively compare the resultsof the BMD method with two recent extensions of DMD with control. Thefirst algorithm is the algebraic DMDc (aDMDc) [26], which extended DMDc toparameter-varying systems described by algebraic, in addition to differential,equations. Including algebraic constraints is very important when consideringstate trajectories generated by aerodynamic solvers capturing unsteady effects,such as in panel methods or unsteady vortex lattice methods [27]. The second al-gorithm is the input-output reduced-order model (IOROM) approach, proposedin [28] to construct data-driven reduced-order LPV models. Improved ways ofdefining the low-dimensional subspace such that state-consistency is achievedwhile preserving accuracy in the (orthogonal) projection were proposed therein.However, the projection operator is the same for all parameters, and is ob-tained from the POD modes as in standard DMD. An extension of IOROMto handle algebraic constraints is also developed here in order to allow for afair comparison. The algorithms are tested on a high-fidelity, fluid-structureinteraction (FSI) numerical model of an airborne wind energy (AWE) morphingwing. The FSI simulator is described in [29] and the wing was analyzed in detailin [30]. Airborne wind energy and morphing wings are paradigmatic examplesof application domains where the system’s response originates from complexinteractions across different domains, and thus could benefit from equation-free approaches. The first type of comparison investigates the accuracy of thereduced-order models to predict various outputs of the wing as the size of the5odel is decreased. In a second set of analyses, models featuring different ordersare used by a model predictive control (MPC) algorithm to track predefined tra-jectories of the airborne wind energy system with the goal of gaining insight intothe trade-off between size and performance of the different schemes. Preliminaryresults of the work were presented in [31].Figure 1 shows a conceptual representation of the proposed data-driven ROMframework. I . Aeroservoelastic system II . Fluid-structure interaction modeling III . System’s trajectories 𝑓 𝐬𝐭𝐚𝐭𝐞 𝑛+1 ,𝐢𝐧𝐩𝐮𝐭 𝑛+1 = 𝑔(𝐬𝐭𝐚𝐭𝐞 𝑛 ,𝐢𝐧𝐩𝐮𝐭 𝑛 ) Mode superposition
Snapshots
Best low-rank approximation
Morphing wing
Aerodynamic pressure
SkinSparsStringersFlexible skin
Actuator
Balanced Mode Decomposition (BMD)
Displacement & velocity IV . Data-driven reduced-order model OutputInput State
Drone u i y i x i Input, state, output trajectories(on parameters grid) Balancing basis and test spacesRegression on the (oblique projected) low-rank subspace Interpolation of grid-based LPV model
Figure 1: Overview of the algorithm and its proposed application: I . illustrative aeroeservoe-lastic testcase [30]; II . typical fluid-structure interaction problem; III . system characterizeduniquely by its states, inputs, and outputs trajectories; IV . sketch of the newly proposedBMD algorithm and comparison with two other algorithms. . Data-driven reduced-order modeling This section provides background material on the tools and concepts relevantto the reduced-order modeling algorithm proposed in this work. In Section 2.1the general data-driven low-order modeling problem is presented. Section 2.2reviews the algebraic DMD with Control (aDMDc) [26], and Section 2.3 reportson the input-output reduced-order model (IOROM) [28], that is the two ROMalgorithms from the literature used here for comparison.
The starting point is a generic discrete-time nonlinear parameter-varyingmodel which can be used to describe typical control systems, such as aeroser-voelastic systems modelled by FSI solvers, x k +1 = f ( x k , u k , ρ k ) ,y k = h ( x k , u k , ρ k ) , (1)where x ∈ R n x , u ∈ R n u , y ∈ R n y are the state, input and output, and ρ : R → R n ρ is a vector of time-varying parameters defining the operating conditions ofthe system. The problem of finding an LPV low-order approximation of (1)can be divided into two phases: first, LTI approximations for frozen values of ρ in a pre-defined grid { ρ j } n g j =1 are computed; then, an LPV model is obtainedthrough interpolation. The following discussion is concerned with the formerphase.It is assumed that for each frozen value ¯ ρ there exists an equilibrium (ortrim) point characterized by the tuple (¯ x (¯ ρ ),¯ u (¯ ρ ),¯ y (¯ ρ )) such that¯ x (¯ ρ ) = f (¯ x (¯ ρ ) , ¯ u (¯ ρ ) , ¯ ρ ) , ¯ y (¯ ρ ) = h (¯ x (¯ ρ ) , ¯ u (¯ ρ ) , ¯ ρ ) . The deviation vectors ˜ x k := x k − ¯ x (¯ ρ ), ˜ u k := u k − ¯ u (¯ ρ ), and ˜ y k := y k − ¯ y (¯ ρ )can then be used as states of an LTI approximation of the system around theequilibrium: ˜ x k +1 = A (¯ ρ )˜ x k + B (¯ ρ )˜ u k , (2a)˜ y k = C (¯ ρ )˜ x k + D (¯ ρ )˜ u k , (2b)7here ( A (¯ ρ ), B (¯ ρ ), C (¯ ρ ), D (¯ ρ )) is a state-space representation completely describ-ing the linearization about the trim point associated with ¯ ρ . The dependence onthe parameter ρ will be dropped in the remainder when clear from the context.In the data-driven setting, the only information on the system comes frominput, state, and output trajectories { x k , u k − , y k − } n s k =1 of length n s . Afterhaving subtracted from them the corresponding trim values, these trajectoriescan be used to form the following snapshot matrices X = (cid:104) x − ¯ x x − ¯ x ... x n s − − ¯ x (cid:105) ∈ R n x × n s ,X = (cid:104) x − ¯ x x − ¯ x ... x n s − ¯ x (cid:105) ∈ R n x × n s ,U = (cid:104) u − ¯ u u − ¯ u ... u n s − − ¯ u (cid:105) ∈ R n u × n s ,U = (cid:104) u − ¯ u u − ¯ u ... u n s − ¯ u (cid:105) ∈ R n u × n s ,Y = (cid:104) y − ¯ y y − ¯ y ... y n s − − ¯ y (cid:105) ∈ R n y × n s . (3)The notation [ X ; U ] will denote the operation of stacking row-wise two ma-trices X and U .The first goal is to obtain a linear time-invariant low-order approximationof (2), that is ˜ z k +1 = F ˜ z k + G ˜ u k , ˜ y k = H ˜ z k + D ˜ u k , where ˜ z ∈ R n z and n z n x . Once this is available, the family of frozen LTIsystems, or directly the signals of interest, are interpolated so that the responseof the system is available at each value ρ k for a generic time-varying trajectoryof the parameter. The algebraic Dynamic Mode Decomposition with Control (aDMDc) algo-rithm was recently proposed in [26] to extend the DMDc algorithm to systemsdescribed by algebraic-differential equations. The DMDc algorithm from [12] isfirst briefly reviewed. This algorithm seeks a data-driven approximation of the8atrices involved in the state equation (2) by means of two truncated singularvalue decompositions (SVD) of the snapshot matrices. The first one is[ X ; U ] = U Σ V (cid:62) ∼ = U r Σ r V (cid:62) r , (4)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 is unrelated to the size of the final reduced-order model,and it could be set for example by using the hard threshold criterion suggestedin [32]. The effect of choosing r on the accuracy of the model will be discussedin the result section. The second truncated SVD is computed from the snapshotmatrix X X = ˆ U ˆΣ ˆ V (cid:62) ∼ = ˆ U n z ˆΣ n z ˆ V (cid:62) n z , (5)where the columns of ˆ U n z are also called POD modes of X and are used forthe projection onto a lower dimensional space. The selection of n z fixes the sizeof the reduced-order model.An approximation of the high-order matrices appearing in (2) can be for-mulated in terms of the truncated SVD (4) and the snapshot matrix X as[ A B ] = X V (cid:62) r Σ − r U (cid:62) r . (6)Then, a low-order approximation is obtained by projecting (6) onto the set ofPOD modes by making use of (5)[ F G ] = (cid:104) ˆ U (cid:62) n z A ˆ U n z ˆ U (cid:62) n z B (cid:105) . Therefore, the low-order model obtained by DMDc is˜ z k +1 = F ˜ z k + G ˜ u k , where ˜ z ∈ R n z is the state of the low-order model and the high-order state canbe recovered by ˜ x = ˆ U n z ˜ z .The aDMDc algorithm [26] builds on the DMDc approach and addresses thepresence of algebraic constraints in the dynamic equations which might arise9hen considering unsteady aerodynamics features. Specifically, the morphingwing analyzed in [26] is described by an FSI solver that implements a 3D panelmethod with a free evolving wake inspired by the method in [27]. This leadsto 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 g ( x k +1 , u k +1 ) = f ( x k , u k , ρ k ) ,y k = h ( x k , u k , ρ k ) , (7)where g is in general a nonlinear function taking into account the dependenceof the states on the control inputs at the next time step. This dependenceresults from algebraic equations relating the doublet strengths (aerodynamicsstates) and downwash (function of the other states and the control inputs). Thiseffect is sometimes accounted for with artificial aerodynamic states by simplychanging the feedtrough matrix to the outputs. However, to correctly capturethe evolution of the states it is important to formulate the problem as stated in(7). The reader is referred to [26] for further discussion on this aspect.The proposed LTI representation of the system accounting for the algebraicconstraints due to the unsteady aerodynamics is˜ x k +1 = A ˜ x k + B ˜ u k + R ˜ u k +1 , where, as in DMDc, the objective is to find a low-order approximation for thestate equation only.The only difference with respect to DMDc is that now the first SVD decom-position is computed with respect to the snapshot matrices X , U , and U ,that is [ X ; U ; U ] = U Σ V (cid:62) ∼ = U r Σ r V (cid:62) r . And the high-order matrices are thus approximated by[
A B R ] = X V (cid:62) r Σ − r U (cid:62) r .
10 low-order approximation is then obtained by projecting (6) onto the same setof POD modes used in DMDc (5)[
F G L ] = (cid:104) ˆ U (cid:62) n z A ˆ U n z ˆ U (cid:62) n z B ˆ U (cid:62) n z R (cid:105) . This procedure results in the aDMDc low-order model˜ z k +1 = F ˜ x k + G ˜ u k + L ˜ u k +1 , (8)where the high-order state can again be obtained from ˜ x = ˆ U n z ˜ z .The approach proposed in the parametrically varying version of the aDMDcalgorithm is to use a different set of POD modes for each value of ρ in thegrid { ρ j } n g j =1 . The frozen LTI models (8) are then simulated simultaneously,the relative states are lifted to the high-order ones using the correspondingprojection matrices (e.g. ˆ U n z ( ρ j ) for the model corresponding to the j − th element in the parameter space), and the state corresponding to the desiredvalue of ρ is obtained by interpolating the high-dimensional states. A firstconsequence of this approach is that the frozen LTI models (8) do not have aconsistent basis for the state. While this has the advantage of projecting overPOD modes specifically computed for a particular value of ρ , it also requiresrunning in parallel all of the low-order models. Moreover, this algorithm doesnot provide an LPV model and thus the use of LPV robust control designstrategies is precluded [17]. While other control techniques, such as modelpredictive control, can still be successfully used (see Section 4.4), the necessityto run in parallel, multiple low-order models, is a drawback of the method whentargeting real-time applications. The input-output reduced-order model (IOROM) algorithm was proposed in[28] to compute a family of state-consistent data-driven low order LTI models(including the output equation) which can be directly parameterized by thevector ρ .Consider first the case when there is no parameter dependence (or equiva-lently, ρ is fixed). Drawing inspiration from the interpretation of DMD as linear11ynamics fitting [4], the main idea is that, given the snapshot matrices (3), thematrices ( A , B , C , D ) defining (2) can be obtained by solving the followingleast-squares problem min A,B,C,D (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) X Y − A BC D X U (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F , (9)where the subscripts F denotes the Frobenius norm of a matrix. Without ap-propriate regularization, this problem would be ill-posed for high-dimensionalsystems ( n x (cid:29) n z is performed by introducing the projection matrix Q ∈ R n x × n z , where Q (cid:62) Q = I n z , such that the orthogonal projection of ˜ x on an n z -dimensionalsubspace is given by QQ (cid:62) ˜ x . Equivalently, one can think that the original stateis approximated by ˜ x ∼ = Q ˜ z for some reduced-order state (or coefficient vector)˜ z ∈ R n z . This results in the following low-order LTI system model˜ z k +1 = ( Q (cid:62) AQ )˜ z k + ( Q (cid:62) B )˜ u k , ˜ y k = ( CQ )˜ z k + D ˜ u k . (10)The vector ˜ z = Q (cid:62) ˜ x ∈ R n z can thus be interpreted as the state of the low-rankapproximation of (2) A BC D ≈ QF Q (cid:62)
QGHQ (cid:62) D = Q I n y F GH D Q (cid:62) I n u . The projection matrix Q is constructed from the POD modes of X , that is Q = U n z , where X ∼ = U n z Σ n z V (cid:62) n z . (11)The least-squares problem giving ( F , G , H , D ) is thenmin F,G,H,D (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) X Y − Q I n y F GH D Q (cid:62) I n u X U (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F , (12)12hose solution is F GH D opt = Q (cid:62) X Y Q (cid:62) X U † , (13)where † denotes the pseudo-inverse of a matrix. It is worth noting that thereduced-order model given by the IOROM algorithm is qualitatively similarto the one associated with DMDc. The main difference (besides the outputequation, not considered in DMDc) is that the pseudo-inverse operation, whichalso amount to an SVD decomposition and thus is conceptually similar to (4),is done here directly on the projected snapshot matrices (unlike in Eq. 4 wherethis applies to X and U ). A minor difference is also that the POD modes arecomputed here with respect to X instead of X .In the parameter-varying case, the regression problem (13) is solved at eachvalue of the parameter grid { ρ j } n g j =1 by taking the corresponding snapshot ma-trices { X ( ρ j ), X ( ρ j ), U ( ρ j ), Y ( ρ j ) } n g j =1 . By always using the same projectionmatrix Q when computing the low-order models at different ρ , state-consistencyis automatically guaranteed because the orthogonal projection has the samerange space. An LPV reduced-order model is then obtained by interpolating(13) across the parameter’s range. That is˜ z k +1 = F ρ k ˜ z k + G ρ k ˜ u k + (¯ z ( ρ k ) − ¯ z ( ρ k +1 )) , ˜ y k = H ρ k ˜ z k + D ρ k ˜ u k , where ( F ρ k , G ρ k , H ρ k , D ρ k ) are obtained by interpolating the corresponding ma-trices for the value of ρ at timestep k . Note that here ˜ z k = z k − ¯ z ( ρ k ), wherethe trim point ¯ z ( ρ k ) = Q (cid:62) ¯ x ( ρ k ) can change as a function of ρ . To correctlytake into account this effect, the term (¯ z ( ρ k ) − ¯ z ( ρ k +1 )) is added [28].Since the choice of a fixed projection matrix is typically associated with lessaccuracy, two strategies are proposed in [28] to alleviate this issue. The firstconsists of using for the decomposition (11) a fat snapshot matrix X obtainedby stacking column-wise the snapshot matrices of multiple parameters. Thismatrix will feature n s n g columns, which can result in computationally expensivecalculations when this number is large. The second, less accurate but more13ractical in case several grid points are analyzed, consists of iteratively building Q by incrementally processing the snapshot data from each grid point in asimilar fashion to the Gram-Schmidt orthogonalization procedure. The formerstrategy is used here when showing results for the IOROM algorithm, togetherwith a linear interpolation of the state-space matrices.
3. Balanced Mode Decomposition with oblique projection algorithm
This section presents the technical aspects of the Balanced Mode Decompo-sition (BMD) algorithm proposed in this paper. Section 3.1 clarifies the goalsand the novelty of the contribution with respect to previous works. Section3.2 presents the algorithm and Section 3.3 details the computation of the twosubspaces defining the oblique projection. Finally, Section 3.4 presents a ver-sion of the algorithm which can handle algebraic constraints and thus allowsthe analyses in Section 4 of the morphing wing with an unsteady aerodynamicsmodel.
The main motivation for the proposal of the BMD algorithm for data-drivenLPV low-order modeling is to address two limitations of recent extensions of thecelebrated DMD method to input-output parameter-varying models. The firstone concerns the selection of Q as subspace for the projection of the higher-order dynamics, which is suboptimal as also acknowledged by the authors of[28]. 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 wherethe system’s state is in balanced coordinates [33]. This is indeed the rationalebehind balanced truncation, which consists of removing the states correspond-ing to the smallest n x − n z Hankel singular values [34]. The justification for thisis that the sum of the Hankel singular values provides a lower bound, and forsystems in balanced coordinates, an upper bound on the error of the approxi-mation achieved by removing system’s states. Even though not guaranteed to14e optimal, balanced truncation is a very effective tool in model-based orderreduction [10]. These ideas are used here to propose a new projection operatorfor the high-dimensional state.Whereas the first aspect is of general concern for POD-based projections,the second aspect addressed with the BMD algorithm is specific to the LPVsetting. Namely, how to handle state-consistency across the frozen LTI models inorder to estimate the system’s response at intermediate points in the parametergrid. In the currently available approaches, this is addressed in two possibleways. When state-consistency is not fulfilled, all reduced-order models are runin parallel by interpolating directly the high-dimensional lifted state. This is thecase of aDMDc, and while it has the advantage that the projection operatorsare parameter-dependent (i.e. at each parameter’s value one can use a differentset of POD modes), an LPV model is not available and moreover computationalefficiency might be compromised. On the other hand, a parameter-independent(i.e. fixed) projection matrix for all frozen models can be used in order toguarantee state-consistency. This is the case for the IOROM algorithm, andit has the drawback that the orthonormal basis associated with the n z mostenergetic modes will be in general different at each value of ρ .The central idea to overcome both of the aforementioned issues is to replacethe orthogonal projection employed in POD-based approaches by an obliqueprojection. Given V ∈ R n x × n z , and W ∈ R n x × n z , such that W is bi-orthogonalto V , i.e. W (cid:62) V = I n z , an oblique projection can be defined by the matrix Π = V W (cid:62) . That is, the oblique projection of ˜ x is given by V W (cid:62) ˜ x . Equivalently, onecan think that the original state is approximated as ˜ x ∼ = V ˜ z (where, as before,˜ z ∈ R n z is the reduced-order state), and the component of ˜ x that is eliminated bythe 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 (spannedby the columns of V ), such that the projection of ˜ x lies in the span of V ;and the test space (spanned by the columns of W ), such that the projection V ˜ z has zero error within it, i.e. W (cid:62) (˜ x − V ˜ z ) = 0. Technically, the high-15imensional state vector is projected along the orthogonal complement of thesubspace spanned by the columns of W onto a subspace spanned by the columnsof V . In practice, this means that what is lost by projecting ˜ x (i.e. the nullspaceof 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 (sincethis defines the state basis, and thus state-consistency will depend only on it)and a parameter-dependent W .The idea of using an oblique projection for LPV model-order reduction wasfirst proposed in [25]. Therein, the setting where a model of the system isavailable (in the form of high-order state-space matrices) is considered, andthus both the construction of V and W , and the computation of the low-ordermodel, is model-based. In the data-driven ROM literature, balancing conceptsare used in two important techniques, namely Balanced POD (BPOD) [35] andthe Eigensystem Realization Algorithm (ERA) [36]. The former is only partiallyequation-free: the controllability Gramian is computed from data, while for theobservability Gramian an adjoint simulation model is needed. Additionally, thehigh-order state-space matrices are required for the balanced projection. For thecase of ERA, a balanced model comes from impulse response simulations of themodel in the spirit of system identification algorithms from realization theory[37]. The ERA algorithm is closely related to BPOD, as it can be interpreted as adata-driven balanced truncation. An important difference is that ERA providesonly the reduced-order model and not the balancing transformation, namely theset of vectors known as balancing and adjoint modes in BPOD. These modes areessentially the counterpart of the basis and test space in BMD, respectively, andare a desirable output of a ROM algorithm as they show the most importantspatio-temporal structures in the dynamics. In an aeroservoelastic setting, thiscan provide insights into efficient design solutions.16 .2. BMD regression problem We consider first the frozen-parameter case and, by virtue of the previouslydiscussed oblique projection, propose the following low-order LTI system model˜ z k +1 = ( W (cid:62) AV )˜ z k + ( W (cid:62) B )˜ u k , ˜ y k = ( CV )˜ z k + D ˜ u k , (14)where the computation of the balancing basis V and test spaces W from system’strajectories will be detailed in Section 3.3. The vector ˜ z = W (cid:62) ˜ x ∈ R n z can thusbe interpreted as the state associated with the following low-rank approximationof (2) A BC D ≈ V F W (cid:62)
V GHW (cid:62) D = V I n y F GH D W (cid:62) I n u . (15)The matrices ( F , G , H , D ) can then be obtained with the following least-squaresproblem min F,G,H,D (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) X Y − V I n y F GH D W (cid:62) I n u X U (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F , (16)which has solution F GH D opt = W (cid:62) X Y W (cid:62) X U † . (17)To build a low-order LPV model, snapshot matrices are first collected forthe values of the parameter in the grid { ρ j } n g j =1 , and the least-squares problem(16) is solved at each grid point. Crucially, the test space W , as motivatedpreviously, is allowed to be a function of ρ . This leads to the following solutionfor the reduced-order models in the grid F ( ρ j ) G ( ρ j ) H ( ρ j ) D ( ρ j ) opt = W (cid:62) ( ρ j ) X ( ρ j ) Y ( ρ j ) W (cid:62) ( ρ j ) X ( ρ j ) U ( ρ j ) † . (18)17he BMD LPV reduced-order model is then obtained by interpolating thefrozen matrices (18) across the parameter’s range˜ z k +1 = F ρ k ˜ z k + G ρ k ˜ u k + (¯ z ( ρ k ) − ¯ z ( ρ k +1 )) , ˜ y k = H ρ k ˜ z k + D ρ k ˜ u k , (19)where, as in IOROM, ( F ρ k , G ρ k , H ρ k , D ρ k ) are obtained by interpolating the cor-responding 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 witheach ρ is in general different, and ¯ z ( ρ k ) = W (cid:62) ¯ x ( ρ k ). Note that, since V is fixed,the basis space is common to all the frozen models and thus the interpolationcan be done at the state-matrices level (as in IOROM). However, the projectionis parameter-dependent due to the use of a parameter-varying test space W ( ρ ). In order to preserve the most important features of the input-output map-ping 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 observ-ability Gramians, respectively W c and W o . This ensures that the projectionpreserves the most observable and controllable states, enabling an approximatedata-driven balanced truncation of the reduced-order LPV model.Because the approach is fully data-driven, empirical Gramians are computedfrom data matrices consisting of appropriate state trajectories using known sys-tems theoretical results. The empirical controllability Gramian can be obtained,following the definition [33, 38], by impulse response simulations (one for eachinput channel). As for the empirical observability Gramian, if the model is lin-ear and its adjoint is available, then it can be computed from impulse responsesimulations (one for each output channel) of the adjoint system, as done inbalanced POD [35]. This computation is identical to the one giving the con-trollability Gramian, but is applied to the adjoint system (for an LTI model instate-space form, this is obtained by replacing A and B by A (cid:62) and C (cid:62) ). If theabove does not hold, for example in case one has only access to the system’s18rajectories and not to the model’s matrices, the approach developed in [39],valid for nonlinear systems and not requiring an adjoint model, can be used.In this method the data matrices used for the Gramian computation consist ofstate trajectories obtained from unforced (zero input) simulations (one for eachstate) obtained by perturbing the initial condition of each state. Since these areunforced responses, when the system is sufficiently damped, it will be generallysufficient to observe only the initial time-steps and thus this calculation can beparallelized and efficiently implemented to reduce the computational time.Once W c and W o are available, a procedure based on [25, Proposition 2] isemployed to compute the test and basis spaces. This construction is reported inthe first part of the pseudocode given in Algorithm 1, which summarizes input,output, and main steps of the BMD algorithm (MATLAB notation is used formatrix operations). For a fixed value of ρ , the construction proposed in [25] isan equivalent procedure to the well known square root algorithm for balancedtruncation [38]. Indeed, it can be noted (see lines 4-10) that the subspace V is taken as a basis for span( L c ˜ U ), where L c ( L o ) is a Cholesky factor of W c ( W o ) and ˜ U consists of the first n z left singular vectors of H = L (cid:62) c L o .The singular values of H are the Hankel singular values of the system andthe SVD decomposition of H plays a fundamental role in balanced reduction[33, 38]. As for the subspace W , it can be shown that balancing is obtainedby taking W = W o V ( V (cid:62) W o V ) − . The expression in line 15 is equivalent butcomputes it with improved numerical robustness [40], by making use of a thinQR factorization (line 12). 19 lgorithm 1 Balanced Mode Decomposition with oblique projection
Input : parameter grid points { ρ j } n g j =1 ; snapshot matrices { X ( ρ j ), X ( ρ j ), U ( ρ j ), Y ( ρ j ) } n g j =1 ; empirical Gramians { W c ( ρ j ), W o ( ρ j ) } n g j =1 ; model desiredorder n z . Output : test space projection matrices { W ( ρ j ) } n g j =1 ; fixed basisspace projection matrix V ; reduced-order models at the grid points { F ( ρ j ), G ( ρ j ), H ( ρ j ), D ( ρ j ) } n g j =1 . for j = 1 , ..., n g do L c ( ρ j ) L c ( ρ j ) (cid:62) = W c ( ρ j ) Cholesky factorization of W c L o ( ρ j ) L o ( ρ j ) (cid:62) = W o ( ρ j ) Cholesky factorization of W o ( U, (cid:63), (cid:63) ) = svd( L c ( ρ j ) (cid:62) L o ( ρ j )) ˜ U = U (: , n z ) ( ¯ U , (cid:63), (cid:63) ) = svd( L c ( ρ j ) ˜ U ) ¯ Q (: , n z ( j −
1) : n z j ) = ¯ U (: , n z ) end for ( U V , (cid:63), (cid:63) ) = svd( ¯ Q ) V = U V (: , n z ) for j = 1 , ..., n g do ( ˜ Q, ˜ R ) = qr( L o ( ρ j ) (cid:62) V ) Thin QR factorization Q = ˜ Q (: , n z ) R = ˜ R (1 : n z , :) W ( ρ j ) = L o ( ρ j ) Q ( R (cid:62) ) − Time-varying test space F ( ρ j ) G ( ρ j ) H ( ρ j ) D ( ρ j ) = W (cid:62) ( ρ j ) X ( ρ j ) Y ( ρ j ) W (cid:62) ( ρ j ) X ( ρ j ) U ( ρ j ) † BMD regression problem end for
The output { F ( ρ j ), G ( ρ j ), H ( ρ j ), D ( ρ j ) } n g j =1 provided by the BMD algorithmis a grid LPV model. After an interpolation algorithm to evaluate the matrices’entries for any value of ρ inside the considered range has been chosen (linear20nterpolation will be used in this work), this model can be used for simulationand control design. Note also that recently proposed robust analysis methodsfor linear-time varying (LTV) systems [41, 42] can be applied to this model,e.g. to investigate specific aircraft manoeuvres. Indeed, by fixing a particulartrajectory 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 theaeroservoelastic modes which have been eliminated and those that have beenkept in the projection, while the parameter-independent basis space can be usedto recover at each time-step k the high-dimensional state via the transformation˜ x k = V ˜ z k .As noted in the introduction, the algorithm provides an approximate bal-anced 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[10], whereas the theoretical order reduction error bounds are only availablefor infinite-time balanced realizations. This source of error can however bemade arbitrarily small by using long enough data sequences for constructingthe Gramians. The slowest decay rate of the system’s impulse responses is thekey parameter to consider when choosing the length of the trajectory [43]. Since the BMD algorithm will be applied in Section 4 to the FSI solverdeveloped in [29], which presents the algebraic constraints described in Section2.2, an extension to handle this instance is presented here. For a fixed value of ρ , the model structure for the high-order model becomes˜ x k +1 = A ˜ x k + B ˜ u k + R ˜ u k +1 , ˜ y k = C ˜ x k + D ˜ u k + P ˜ u k +1 , where a potential effect of the algebraic constraints in the output equation isalso considered via the matrix P (in the analyses of the morphing wing this21atrix was, as expected, always found to be zero). Therefore, the low-orderapproximation becomes now A B RC D P ≈ V F W (cid:62)
V G V LHW (cid:62)
D P = V I n y F G LH D P W (cid:62) I n u
00 0 I n y . The new objective function to be minimized ismin
F,G,L,H,D,P (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) X Y − W I n y F G LH D P W (cid:62) I n u
00 0 I n y X U U (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F , and the new optimal solution is F G LH D P = W (cid:62) X Y W (cid:62) X U U † . (20)In a similar way, the IOROM algorithm has also been extended to the algebraic-differential case in order to allow its application to the test case considered inthe following.
4. Results
This section presents and discusses results obtained by applying the threealgorithms aDMDc, IOROM, and BMD to the flexible and highly camberedmorphing 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 decreasethe camber, thus replacing conventional ailerons. The reader is referred to therelated previous works for details on the wing design [44] and its investigationwith FSI tools [30].
The high-fidelity fluid-structure interaction model of the morphing wing ispresented in detail in [29, 45]. The structural model is based on a combination of22inear plate and beam elements. The external skin of the lifting surface and theVoronoi-based internal structure are modelled using plates, while the stringersare modelled with beam elements. The stiffness and mass matrices are obtainedwith the commercial software Nastran [46], and coupled via thin plate splineand inverse distance weighting with the aerodynamic model. The aerodynamicmodel is based on a 3D unsteady panel method [47]. The unsteadiness of theflow is represented by shedding at each time step a new row of vortex ringsingularities after the trailing edge. All the other wake nodes are then moved, viaa second-order RungeKutta integration scheme, using the local velocity of theflow. The aerodynamic forces on the surface are computed with the coefficientof pressure on each panel, considering the far field velocity, the induced velocityby the wing itself, and the induced velocity by all the wake panels.The state of the system x consists of the total number of structural modesof 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 u = [ α ; p ; q ; r ; F s ; F as ] , (21)where α is the angle of attack, p , q , and r are the roll, pitch, and yaw rota-tion rates and F s and F as are the (normalized) symmetric and anti-symmetricmorphing actuation inputs. Their value is associated with a camber deforma-tion and is thus related to a trailing edge deflection: specifically, the amount ofupwards (negative value) or downwards (positive value) deflection.Unless otherwise specified, we will consider as output only the first bendingmode of the wing ( n y =1), since this is usually the one associated with dynamicinstabilities and large deformations, and thus it is of particular interest for activecontrol tasks [48, 49]. The flight speed will be considered as the time-varyingparameter ( n ρ =1).The training phase, common to both the algorithms and consisting of gener-ating the snapshot matrices in Eq. (3), is carried out by exciting the system witha series of impulses deployed in random order in all input channel. Trajectoriesare of length n s =500 and are recorded with sampling time 0.006 s.23 .2. Fixed-parameter models In this first set of tests, the accuracy of the different models at fixed values ofthe flight speed V is assessed. The models are simulated using sinusoidal inputsin each input channel. The frequencies, different for each channel, are all of theorder of magnitude of the aircraft reduced frequency ω r =: V ¯ c , where ¯ c =0.29 mis the mean chord of the wing. This test is performed for 3 flight speeds in therange of operating conditions of interest, namely V =30 m/s, V =40 m/s, and V =50 m/s. To quantify the accuracy as a function of the order of the model n z ,the Euclidean norm of the error signal between the first bending mode amplitudeprovided by the high-fidelity FSI and each of the predictions obtained with thethree ROM algorithm is computed. The result is shown in Fig. 2. Figure 2: Error on the prediction of the system’s output (bending mode) for three differentvalues of the flight speed as a function of the model order.
It is clear that for all speeds BMD provides the smallest error for a verylow-order approximation of the full dynamics, as expected in view of the choiceof low-dimensional subspace where the high-dimensional data are projected. Asthe size n z of the system increases, the difference between the algorithms is lessnoticeable and, for high enough orders, the algorithms tend to give same results.24t is also noted that the results obtained with the aDMDc algorithm showedgreat sensitivity, in the range of n z displayed in Fig. 2, to the SVD truncationorder r employed in Eq. (4). Using the hard threshold criterion from [32] didnot provide good results as it resulted in a very large r (therefore the truncationincluded very low singular values deteriorating the approximation). To presenta fair comparison with the other two methods, the value of r was fine-tuned inorder for aDMCc to provide the best results possible (an effective heuristic wasto choose r = n z +30). In this way, any gap in performance cannot be ascribedto numerical accuracies depending on the selection of r . In the second set of tests, the accuracy during parameter-varying manoeuvresis 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 reduced-order models are obtained using snapshot matrices obtained gridding the flightspeed range every 2 m/s and thus using 16 different speeds ( n g =16). The same sinusoidal excitation signals used in Section 4.2 are consideredhere. In Fig. 3, the bending mode amplitude response obtained with the FSIsolver (
FSI ) is compared with the predictions of the three algorithms when theorder of the models is fixed at n z =14. All the signals are normalized by thelargest value of the bending amplitude measured in the FSI simulation.The plot confirms, also in the LPV setting, that the BMD algorithm guaran-tees 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 ituses a parameter-varying set of POD modes. However, aDMDc does not providea family of interpolated low-order models, and interpolates directly the high-order states, thus requiring parallel simulations of the low-order models. Thebetter performance of BMD, despite the fact that a part of the projection (theone related to the basis space) is constant, is ascribed to the improved selection25 igure 3: Comparison of the normalized output (bending mode) for parameter-varying simu-lations (flight speed linearly increasing from 20 m/s to 50 m/s) of models with size 14. of subspace for the projection compared to the standard POD one, common tothe other two methods. In addition to the improvement in the accuracy, theBMD 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 ingeneral real-time applications. This section investigates the accuracy of the reduced-order models for dif-ferent types of input signals. The Euclidean norm of the error signal betweenthe first bending mode amplitude provided by the high-fidelity FSI and theprediction obtained with each of the three ROM algorithms is again used asmetric to assess the quality of the approximation. Three classes of inputs areconsidered:
Sine coincides with the signal tested so far and already investigatedin [26];
Chirp excites the system by injecting in all 6 input channels definedin (21) a chirp signal with frequency linearly varying from 0.1 ω r to ω r ; PRBS excites the system by injecting in all 6 input channels a PRBS-9 sequence. This26ast input, namely a Pseudo-Random Binary Signal (PRBS) [50], is a determin-istic signal with white-noise-like properties. It is very well known in the systemidentification and experiment design fields since it has the favourable propertyof equally distributing energy across all the frequencies in the input spectrum.In this way, information on the models in different frequency ranges can be ex-tracted. Although not a common input in aeroservoelastic applications, it hasbeen used in this spirit here, since the previously adopted sets of input will onlygive information on the behaviour of the reduced models around the aircraftreduced frequency ω r . Results are shown in Fig. 4. Figure 4: Error of the prediction of the system’s output (bending mode) for speed-varyingmanoeuvres with three different types of input signals as a function of the model order.
The plots confirm the advantage in using the BMD approach when seekinga low-order model capturing parameter variations. These results are importantconsidering that they are obtained by exploring different frequency ranges ofthe system’s response. 27 .3.3. Effect of the parameter grid
In this section we analyze the effect of the flight speed grid where the reduced-order LTI models are computed, i.e. the selection of the parameter n g . Thisis an important aspect, which is known to influence both the accuracy of theLPV models and the quality of the control design based on them. Three casesare 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 =48m/s and then V =50 m/s; n g =16 which is the grid used so far. Figure 5: Error on the prediction of the system’s output (bending mode) for speed-varyingmanoeuvres with values of n g as a function of the model order. From the analyses it can be gathered that aDMDc is more robust than theother algorithms to the value of n g . In particular, both IOROM and, in a moreaccentuated way, BMD present poor performance for a few reduced order modelsin the range of n z between 30 and 40 when the flight speed grid is coarser. Thereason for this behaviour is due to the interpolation approaches employed bythe three ROM schemes. Whereas IOROM and BMD interpolate the low-orderstate-space matrices obtained at the grid points, aDMDc interpolates directly28he high-order vector states which are obtained by lifting the low-order states ˜ z from the local models (8) running in parallel. Interpolating every entry of thestate-space matrices therefore makes the choice of the grid a more delicate aspectin IOROM and BMD. The unstable behaviours resulting in very high (out of theplot) errors are indeed ascribed to numerical inaccuracies in this interpolation.It has been observed that the entries of the matrices are overall bigger as theorder n z is increased, hence justifying why these outliers take place at in theaforementioned range of model’s orders. Except for these numerical problems,BMD shows better performance even when very coarse grids are employed. The capability of the models to predict other quantities of interest, such asfor example aerodynamic coefficients depending on the system’s states, is inves-tigated. In particular, we test the accuracy when these coefficients are computeddirectly from the states. That is, the low-order states are lifted to the high-orderones, which are then used to compute the coefficients using their known rela-tionship to the states. While this is the only possible way of reconstructingthe system’s signals for aDMDc, in IOROM and BMD this can alternatively bedone by simply adding the desired quantities to the vector of outputs (as alreadydone for the bending mode amplitude before) and generating new reduced-ordermodels. This would probably be the preferred approach if the signals are usedfor control (either because they represent measurements fed to the controller orbecause they are performance measures to be optimized).Figure 6 shows the normalized lift ( C L ), pitch ( C M ), and drag ( C D ) coeffi-cients for the same constant acceleration manoeuvre considered in the previoussections and with a sinusoidal excitation. Normalization is performed, as doneearlier in Figure 3, by dividing each signal by the largest value of the corre-sponding signal in the FSI simulation.The same observations gathered earlier with respect to the trajectory ofthe bending mode are confirmed here. It is particularly interesting to observethat, even though these coefficients are not outputs of the model, and thus29 igure 6: Comparison of the normalized lift, pitch, and drag coefficients for a parameter-varying simulation with n z =14 and sinusoidal inputs. the balancing projection is not aimed directly at capturing them, the BMDalgorithm is still able to perform better than the others. Figure 7 shows thesame analyses when chirp signals are used as input to excite the model. Similarconclusions can be drawn. In this section, a control application of the morphing wing’s low-order modelsis investigated. Specifically, model predictive control (MPC) [51] is considered,given its well established use in the AWE field [52]. Two distinct referencetracking problems are examined, where pre-defined lift and first bending modeamplitude profiles are tracked while flying trajectories over a range of differ-ent flight speeds and in the presence of turbulence and gusts. The analysis ofthese manoeuvres is motivated by the interest in using active control to guar-antee a safe operation for the AWE system (with respect to some of its criticalcomponents such as the wing or tether) by keeping indicators of the structuralintegrity close to desired, and possibly pre-optimized, values. This can avoid30 igure 7: Comparison of the normalized lift, pitch, and drag coefficients for a parameter-varying simulation with n z =14 and chirp inputs. passive remedies such as reducing the load transmitted to the ground station,which in turn decreases the amount of wind energy harvested. Having effec-tive and reliable control laws to guarantee the integrity of the AWE system canrepresent an important enabler for this technology [53].In its basic form, model predictive control repeatedly solves a finite-horizonoptimal control problem of length N c subject to input and state-constraints. Ateach instant, a model of the system is employed to predict its response and thusselect the control sequence ( u i ) N c − i =0 which minimizes the cost J MPC = N c − (cid:88) k =0 (cid:16) (cid:107) ˜ y k − r k (cid:107) N + (cid:107) ˜ u k (cid:107) M + (cid:107) ∆˜ u k (cid:107) M ∆ (cid:17) , (22)where r is the reference trajectory, and for a vector x , we denote by (cid:107) x (cid:107) P the weighted l -norm ( x T P x ) . Besides the terms penalizing deviation of theoutput ˜ y from r (with the weighting matrix N ∈ R n y × n y ) and control effort(with the weighting matrix M ∈ R n u × n u ), the cost in (22) also penalizes fastchanges in the input via the term ∆˜ u k = ˜ u k − ˜ u k − (e.g. to avoid actuator ratesaturation). 31he following optimization problem will be solved to obtain the optimalinput sequence minimize ( u i ) Nc − i =0 , ( y i ) Nc − i =0 J MPC (˜ u, ˜ y ) , (23a)subject to ˜ y i = f (˜ u i , ˜ u i +1 ) , (23b)(˜ u i ) N − i =0 ∈ U , (23c)where: the cost function (23a) is defined in (22); the constraint set U (23c)consists of minimum and maximum constraints on the input; and (23b) enforcesthe dynamic constraint that relates the sequence of input to the output via aninput-output model of the system f . Precisely, f will be formulated here byusing the reduced-order aDMDc, IOROM, and BMD models. The goal of theanalyses is to compare the associated closed-loop cost J CLMPC , that is the cost(22) incurred by the true system (simulated here by the high-dimensional FSIsolver) when this is regulated by the inputs optimized solving problem (23).Since the model is used to predict the system’s output, and thus select the inputsequence, any mismatch between model and system can result in a degradationof the controller performance.The analysis considers the case where the morphing wing (Fig. 1, I) fliesa 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 flightspeed with gust length 0.5 s corresponding to a 1 degree deflection of the angleof attack α , and to turbulence generated with a Dryden filter (Fig. 8-left). Theright plot in Fig. 8 depicts the lift and the first bending mode’s amplitudeprofiles tracked by the MPC algorithm. Note that a unitary value of the firstbending mode corresponds to a wingtip displacement of 4.6 cm.The scenario considered here has only a symmetric morphing actuation in-put, i.e. ˜ u = 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 thetrailing edge in the range ± y is either the first bending mode, or the lift force generated by32 Figure 8: Speeds variations during the manoeuvres (left plot) and reference tracking profiles(right plot). the wing, depending on the case (i.e. n y =1). The control horizon is N c =10 andthe 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 theoutput deviation is increased in the latter case due to the order of magnitudedifference of the two tracked quantities (recall Fig. 8-right). Fig. 9 shows thecomparison of the costs. For the sake of clarity, the closed-loop cost J CLMPC hasbeen normalized in each case by dividing it by the corresponding value obtainedwith the BMD algorithm when n z =40.The observations gathered in the previous sections regarding the better pre-diction performance achieved with the BMD algorithm when low-order approxi-mations are considered are confirmed here in the context of control applications.Both plots show that, while for higher orders the closed-loop costs have verysimilar values, when the size of the model is decreased the BMD gives in generalthe lowest cost. Another interesting observation is that the lift tracking problemis quite robust to the use of low-order models. Indeed, the closed-loop costs are33 (a) Lift.
10 15 20 25 30 35 40024681012 (b) First bending mode.Figure 9: Normalized closed loop costs for the two MPC tracking problems as a function ofthe model order. always within two times of the lowest cost achieved at n z =40 except for thecase of aDMDc at n z =10). On the other hand, the bending tracking problemis 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 moti-vates further work on the use of low-order models for control of coupled flexiblestructures like those encountered in AWE applications. In the real-time controlsetting, it is important to stress that by using BMD and IOROM models an34rder of magnitude computational speed-up was achieved with respect to thecases where aDMDc models were used. This is because the aDMDc models havethe requirement of running several models in parallel.
5. Conclusion
The paper proposes the Balanced Mode Decomposition with oblique projec-tion algorithm, a novel data-driven algorithm for constructing low-order LPVmodels from system’s trajectories. The problem is presented and discussed andtwo recent algorithms from the literature, aDMDc and IOROM, are consideredfor comparison since they both have connections with the newly proposed ap-proach. Technical details on the BMD algorithm are given in order to clearlypoint out the innovations, and the advantages with respect to previous work.The performance of the BMD algorithm is assessed on a morphing wing forairborne wind energy applications. The results, proposed both for the fixed pa-rameter and, more extensively, for the parameter-varying case, confirm the the-oretical advantages discussed in the technical part of the paper. When seekinglow-order model representations, the BMD approach achieves, among the testedalgorithms, the lowest prediction error and best control performance when usedas model for an on-line MPC scheme. The improved accuracy is ascribed to theuse of a projecting subspace that balances the low-order states (this element is ofinterest also in a fixed-parameter setting), and to the use of a parameter-varyingprojection operator (which can thus be enriched with parameter-dependent fea-tures, instead of being fixed throughout the range). This has the advantageousfeature of being achieved while guaranteeing state-consistency.
Acknowledgements
This work is supported by the Swiss National Science Foundation undergrant no. 200021 178890. The authors wish to thank Nivethan Yogarajah andEva Ahbe for helpful discussions. 35 eclaration of competing interest
The authors declare that they have no conflict of interest.
References [1] S. L. Brunton, B. R. Noack, P. Koumoutsakos, Machine learning for fluidmechanics, Annual Review of Fluid Mechanics 1 (2) (2020) 391–421.[2] P. Holmes, J. L. Lumley, G. Berkooz, Turbulence, Coherent Structures,Dynamical Systems and Symmetry, Cambridge Monographs on Mechanics,Cambridge University Press, 1996.[3] P. Schmid, Dynamic mode decomposition of numerical and experimentaldata, Journal of Fluid Mechanics 656 (2010) 5–28.[4] J. H. Tu, C. W. Rowley, D. N. Luchtenburg, S. L. Brunton, J. N. Kutz,On dynamic mode decomposition: Theory and applications, Journal ofComputational Dynamics 52 (1) (2014) 477–508.[5] G. Berkooz, P. Holmes, J. L. Lumley, The proper orthogonal decompositionin the analysis of turbulent flows, Annual Review of Fluid Mechanics 25 (1)(1993) 539–575.[6] P. Schmid, D. Henningson, Stability and Transition in Shear Flows,Springer-Verlag New York, 2001.[7] A. Iannelli, R. S. Smith, The role of the state in model reduction withsubspace and POD-based data-driven methods, in: IEEE American ControlConference, 2021.URL http://people.ee.ethz.ch/~andreaia/pub/ACC21_state.pdf [8] S. Le Clainche, D. Rodrguez, V. Theofilis, J. Soria, Flow around ahemisphere-cylinder at high angle of attack and low Reynolds number.Part II: POD and DMD applied to reduced domains, Aerospace Scienceand Technology 44 (2015) 88 – 100.369] C. Hu, C. Yang, W. Yi, K. Hadzic, L. Xie, R. Zou, M. Zhou, Numericalinvestigation of centrifugal compressor stall with compressed dynamic modedecomposition, Aerospace Science and Technology 106 (2020).[10] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, Societyfor Industrial and Applied Mathematics, 2005.[11] C. P. Moreno, P. J. Seiler, G. J. Balas, Model reduction for aeroservoelasticsystems, Journal of Aircraft 51 (1) (2014) 280–290.[12] J. L. Proctor, S. L. Brunton, J. N. Kutz, Dynamic mode decompositionwith control, SIAM Journal on Applied Dynamical Systems 15 (1) (2016)142–161.[13] R. Toth, Modeling and identification of linear parameter-varying systems,Springer, 2010.[14] A. Packard, Gain scheduling via linear fractional transformations, Systems& Control Letters 22 (2) (1994) 79 – 92.[15] S. Bennani, D. M. C. Willemsen, C. W. Scherer, Robust control of linearparametrically varying systems with bounded rates, Journal of Guidance,Control, and Dynamics 21 (6) (1998) 916–922.[16] T. He, A. K. Al-Jiboory, G. G. Zhu, S. S.-M. Swei, W. Su, Applicationof ICC LPV control to a blended-wing-body airplane with guaranteed H ∞ performance, Aerospace Science and Technology 81 (2018) 88 – 98.[17] A. Hjartarson, P. Seiler, A. Packard, LPVTools: A Toolbox for Modeling,Analysis, and Synthesis of Parameter Varying Control Systems, in: 1stIFAC Workshop on Linear Parameter Varying Systems LPVS, 2015.[18] P. Benner, S. Gugercin, K. Willcox, A survey of projection-based modelreduction methods for parametric dynamical systems, SIAM Review 57 (4)(2015) 483–531. 3719] C. Poussot-Vassal, C. Roos, Generation of a reduced-order LPV/LFTmodel from a set of large-scale MIMO LTI flexible aircraft models, ControlEngineering Practice 20 (9) (2012) 919–930.[20] G. Ferreres, Computation of a Flexible Aircraft LPV/LFT Model UsingInterpolation, IEEE Transactions on Control Systems Technology 19 (1)(2011) 132–139.[21] A. K. Al-Jiboory, G. Zhu, S. S.-M. Swei, W. Su, N. T. Nguyen, LPVmodeling of a flexible wing aircraft using modal alignment and adaptivegridding methods, Aerospace Science and Technology 66 (2017) 92 – 102.[22] Q. Zhang, L. Ljung, R. M. Pintelon, On Local LTI Model Coherence forLPV Interpolation, IEEE Transactions on Automatic Control.[23] C. D. Villemagne, R. E. Skelton, Model reductions using a projection for-mulation, International Journal of Control 46 (6) (1987) 2141–2169.[24] P. van Overschee, L. de Moor, Subspace identification for linear systems:theory, implementation, applications, Kluwer Academic Publishers, 1996.[25] J. Theis, P. Seiler, H. Werner, LPV Model Order Reduction by Parameter-Varying Oblique Projection, IEEE Transactions on Control Systems Tech-nology 26 (3) (2018) 773–784.[26] N. Fonzi, S. L. Brunton, U. Fasel, Data-driven nonlinear aeroelastic modelsof morphing wings for control, Proceedings of the Royal Society A: Math-ematical, Physical and Engineering Sciences (2020).[27] J. Murua, R. Palacios, J. M. R. Graham, Applications of the un-steady vortex-lattice method in aircraft aeroelasticity and flight dynamics,Progress in Aerospace Sciences 55 (2012) 46 – 72.[28] J. Annoni, P. Seiler, A method to construct reduced-order parameter-varying models, International Journal of Robust and Nonlinear Control27 (4) (2017) 582–597. 3829] U. Fasel, P. Tiso, D. Keidel, G. Molinari, P. Ermanni, Reduced-OrderDynamic Model of a Morphing Airborne Wind Energy Aircraft, AIAAJournal 57 (8) (2019) 3586–3598.[30] U. Fasel, D. Keidel, G. Molinari, P. Ermanni, Aeroservoelastic Optimiza-tion of Morphing Airborne Wind Energy Wings, in: AIAA Scitech Forum,2019.[31] A. Iannelli, U. Fasel, N. Yogarajah, R. Smith, A Balanced Mode De-composition Approach for Equation-Free Reduced-Order Modeling of LPVAeroservoelastic Systems, in: AIAA SciTech Forum, 2021.[32] M. Gavish, D. L. Donoho, The optimal hard threshold for singular values is4 / √
3, IEEE Transactions on Information Theory 60 (8) (2014) 5040–5053.[33] B. Moore, Principal component analysis in linear systems: Controllabil-ity, observability, and model reduction, IEEE Transactions on AutomaticControl 26 (1) (1981) 17–32.[34] K. Glover, All optimal Hankel-norm approximations of linear multivariablesystems and their L ∞ -error bounds, International Journal of Control 39 (6)(1984) 1115–1193.[35] C. W. Rowely, Model reduction for fluids, using balanced proper orthogo-nal decomposition, International Journal of Bifurcation and Chaos 15 (03)(2005) 997–1013.[36] Z. Ma, S. Ahuja, C. W. Rowley, Reduced-order models for control of flu-ids using the eigensystem realization algorithm, Theoretical and Computa-tional Fluid Dynamics 25 (1) (2011) 233–247.[37] M. Viberg, Subspace-based methods for the identification of linear time-invariant systems, Automatica 31 (12) (1995) 1835 – 1851.[38] A. Laub, M. Heath, C. Paige, R. Ward, Computation of system balancingtransformations and other applications of simultaneous diagonalization al-gorithms, IEEE Transactions on Automatic Control 32 (2) (1987) 115–122.3939] S. Lall, J. E. Marsden, S. Glavaˇski, A subspace approach to balanced trun-cation for model reduction of nonlinear control systems, International Jour-nal of Robust and Nonlinear Control 12 (6) (2002) 519–535.[40] G. Golub, C. Van Loan, Matrix Computations, Johns Hopkins UniversityPress, 2013.[41] P. Seiler, R. Moore, C. Meissen, M. Arcak, A. Packard, Finite HorizonRobustness Analysis of LTV Systems Using Integral Quadratic Constraints,Automatica 100 (2019) 135–143.[42] A. Iannelli, P. Seiler, A. Marcos, Worst-case disturbances for time-varyingsystems with application to flexible aircraft, Journal of Guidance, Control,and Dynamics 42 (6) (2019) 1261–1271.[43] I. Markovsky, J. C. Willems, P. Rapisarda, B. L. De Moor, Algorithms fordeterministic balanced subspace identification, Automatica 41 (5) (2005)755 – 766.[44] D. Keidel, U. Fasel, L. Baumann, G. Molinari, P. Ermanni, Experimentalvalidation of a morphing wing for Airborne Wind Energy Applications,in: 8th International Conference on Adaptive Structures and Technologies,2017.[45] U. Fasel, Reduced-order aeroelastic modeling of morphing wings for opti-mization and control, Ph.D. thesis, ETH Zurich (2020).[46] W. P. Rodden, E. H. Johnson, User Guide V 68 MSC/NASTRAN Aeroe-lastic Analysis, 1994.[47] A. P. J. Katz, Low-speed aerodynamics, Cambridge University Press, 2001.[48] S. Waitman, A. Marcos, H ∞∞