Model Reduction Framework with a New Take on Active Subspaces for Optimization Problems with Linearized Fluid-Structure Interaction Constraints
RReceived: Added at production Revised: Added at production Accepted: Added at productionDOI: xxx/xxxx
ARTICLE TYPE
Model Reduction Framework with a New Take on ActiveSubspaces for Optimization Problems with LinearizedFluid-Structure Interaction Constraints
Gabriele Boncoraglio* | Charbel Farhat | Charbel Bou-Mosleh Department of Aeronautics andAstronautics, Stanford University, DurandBuilding, 496 Lomita Mall, Stanford, CA94305-4035, USA Department of Mechanical Engineering,Stanford University, Building 530, 440Escondido Mall, Stanford, CA 94305-3030,USA Institute for Computational andMathematical Engineering, StanfordUniversity, Huang Engineering Center, 475Via Ortega, Suite 060, Stanford, CA94305-4042, USA Department of Mechanical Engineering,Notre Dame University-Louaize, ZoukMosbeh, Lebanon
Correspondence *Gabriele Boncoraglio , Department ofAeronautics and Astronautics, StanfordUniversity, Durand Building, 496 LomitaMall, Stanford, CA, 94305-4035.Email: [email protected]
Summary
In this paper, a new take on the concept of an active subspace for reducing the dimen-sion of the design parameter space in a multidisciplinary analysis and optimization(MDAO) problem is proposed. The new approach is intertwined with the concepts ofadaptive parameter sampling, projection-based model order reduction, and a databaseof linear, projection-based reduced-order models equipped with interpolation onmatrix manifolds, in order to construct an efficient computational framework forMDAO. The framework is fully developed for MDAO problems with linearized fluid-structure interaction constraints. It is applied to the aeroelastic tailoring, under flutterconstraints, of two different flight systems: a flexible configuration of NASA’s Com-mon Research Model; and NASA’s Aeroelastic Research Wing
KEYWORDS: active subspace, constrained optimization, fluid-structure interaction, gradient-based optimization, inter-polation on matrix manifolds, model reduction, parameter sampling
Multidisciplinary analysis and optimization (MDAO) problems with fluid-structure interaction (FSI) constraints arise in manyengineering applications including design and optimal control. Solving such problems requires repeated evaluations of theFSI constraints when using, for example, the nested analysis and design (NAND) approach. Typically, each evaluation of FSIconstraints incurs the solution of some coupled computational fluid dynamics (CFD) - computational structural dynamics (CSD)problems. In an industrial setting, these CFD-CSD problems are usually large-scale when high-fidelity simulations are required.Furthermore, the number of solutions of coupled CFD-CSD problems increases when the design parameter space is high-dimensional. This is because the number of evaluations of the objective function and constraints incurred by an MDAO problemgenerally increases with the dimension 𝑁 of . For all these reasons, solving highly parameterized, FSI-constrained, MDAOproblems can be very computationally intensive. Consequently, the main objective of this paper is to present an innovative,comprehensive, computational framework for efficiently solving such problems when the FSI constraints are linearized in thespatial, temporal, or frequency domain. Such linearizations are common when the FSI constraints pertain to stability or controlproblems. a r X i v : . [ m a t h . NA ] F e b Boncoraglio
ET AL
The computational cost of a linearized, coupled FSI problem can be reduced by substituting the underlying coupled high-dimensional model (HDM) with a less computationally intensive surrogate model such as a linear, coupled, projection-based,reduced-order model (PROM). Such a linear PROM can be local (or pointwise), or global . In the context of this paper, the labellocal describes a PROM that was trained at a single parameter point 𝝁 ∈ ⊂ ℝ 𝑁 of the design parameter space . On theother hand, the label global refers to a PROM that was trained at multiple points in so that it is robust and accurate within alarge region of this parameter space.The overhead cost – also known as the offline cost – associated with the construction of a local PROM is lower than thatassociated with the construction of a global PROM, given that it involves training at a single parameter point. Unfortunately,a local PROM does not usually perform well away from the training point. Therefore, it is not suitable for the solution of anMDAO problem unless it is reconstructed everytime a new parameter point is queried, which is inefficient. Alternatively, aglobal PROM may be suitable for MDAO problems, from the robustness and accuracy viewpoints. However, such a PROM isbound to be of a larger dimension than needed at any parameter point visited online by the optimization trajectory, given that itwas trained offline at multiple points of the parameter space; this also hinders computational efficiency.In order to overcome the inefficiency issues highlighted above, different approaches have been proposed in the literature. Inthe context of linear PROMs, two of them are noteworthy. First, the adaptive approach consisting in progressively updating aninitial PROM by continuously training it at a subset of the parameter points previously visited by the optimization trajectory .In this approach, which breaks the traditional offline-online framework of model order reduction, the reduced optimizationproblems are equipped with a nonlinear trust-region based on a residual error indicator to keep the optimization trajectory in aregion of where the PROM is accurate. Hence, whereas this approach is equally applicable to linear PROMs, it was designedprimarily for nonlinear PROMs. The second approach consists in carefully sampling the parameter space at a small set ofpoints, constructing at each sampled parameter point a local PROM, and storing in a database all constructed local PROMs aswell as some related quantities such as the underlying reduced-order bases (ROBs) . This approach is natural for linear PROMsas once the database has been constructed offline, a linear PROM – as well as its sensitivities – can be constructed in real-time at any unsampled parameter point queried by the optimization trajectory using interpolation on matrix manifolds . Thisapproach has been recently demonstrated for MDAO problems characterized by a small-dimensional design parameter space –say 𝑁 ≤ – and shown to deliver excellent speedup factors. For high-dimensional design parameter spaces however – say 𝑁 ≥ , building a database of local PROMs suffers from the curse of dimensionality and rapidly becomes impractical if notinfeasible. For example, sampling only 2 parameter points along each dimension of a design parameter space of dimension 10leads to the construction of = 1 , PROMs! Adaptive parameter sampling procedures – also known as greedy sampling(or simply greedy) procedures – and particularly those equipped with a saturation technique attenuate somehow this issue butinsufficiently. Alternatively, this paper addresses this issue by focusing on a PROM computational framework for MDAO wherethe concept adaptive parameter sampling, that of a database of linear PROMs, and a new take on the concept of active subspaces(AS) are intertwined.The remainder of this paper is organized as follows. Section 2 formulates the problem of interest in sufficient details to keepthis paper as self-contained as possible. Section 3 reviews and intertwines the concepts of AS, projection-based model orderreduction (PMOR), and a database of linear FSI PROMs equipped with interpolation on matrix manifolds, in the context ofMDAO problems with linearized FSI constraints. It also proposes an effective parameter sampling method for constructing anAS. Section 4 illustrates the proposed PROM-based computational framework for MDAO with the aeroelastic tailoring, underflutter constraints, of two different systems: a flexible configuration of NASA’s Common Research Model (CRM); and NASA’sAeroelastic Research Wing ARW-2. The main focus of this paper is on the acceleration of the solution of MDAO problems under a high-dimensional, linearized,FSI constraint, as well as other linear or nonlinear constraints, including partial differential equation (PDE)-based constraints. oncoraglio
ET AL Such problems can be written as min q , 𝐳 , 𝝁 𝑓 ( q , 𝐳 , 𝝁 ) s.t. 𝐜 ( q , 𝐳 , 𝝁 )) ≤ 𝐿 ( q , 𝝁 ) = 𝑁𝐿 ( 𝐳 , 𝝁 ) = (1)where: • 𝑓 ( ⋅ ) represents a scalar objective function to minimize. • q ∈ ℝ 𝑁 𝑞 denotes a first high-dimensional vector of 𝑁 𝑞 semi-discrete or discrete fluid and structural state variables involvedin a high-dimensional, parametric, linearized, FSI constraint. • 𝐳 ∈ ℝ 𝑁 𝑧 denotes a second high-dimensional vector of 𝑁 𝑧 semi-discrete or discrete state variables involved in a high-dimensional, parametric, PDE-based constraint. • 𝝁 ∈ ⊂ ℝ 𝑁 denotes a vector of 𝑁 design optimization parameters and is assumed to be relatively large – say, 𝑁 ≥ .Each component of this vector may represent some structural material property or an aerodynamic shape parameter. • 𝐜 ( ⋅ ) represents a general set of linear and nonlinear algebraic constraints. • 𝐑 𝐿 ( q , 𝝁 ) is a high-dimensional, semi-discretized or discretized, parametric PDE that is linear in q , and the subscript 𝐿 designates the linear aspect. • 𝐑 𝑁𝐿 ( 𝐳 , 𝝁 ) is a high-dimensional, semi-discretized or discretized, parametric PDE that is nonlinear in 𝐳 , and the subscript 𝑁𝐿 designates the nonlinear aspect.Specifically, the focus is set on the gradient-based solution of the MDAO problem (1) using the NAND approach. Two methodologies for solving MDAO problems such as (1) can be considered: Simultaneous Analysis and Design (SAND)and Nested Analysis and Design (NAND). In the SAND approach, both vectors of discrete state variables 𝐪 and 𝐳 as well asthe vector of design parameters 𝝁 are considered optimization parameters for problem (1). In this case, the MDAO problemhas potentially millions of optimization parameters given the high-dimensionality of 𝐑 𝐿 ( q , 𝝁 ) and 𝐑 𝑁𝐿 ( 𝐳 , 𝝁 ) . In the NANDapproach however, the only optimization parameters are the components of 𝝁 : in this case, 𝐪 and 𝐳 are recognized as implicitfunctions of 𝝁 due to the constraints 𝐑 𝐿 ( q , 𝝁 ) = and 𝐑 𝑁𝐿 ( 𝐳 , 𝝁 ) = . It follows that when storage is an issue, the NANDapproach is preferrable. Furthermore, the NAND approach has the additional benefit or enabling the reuse of available solversfor the PDE-based constraints. For these reasons, it is here the method of choice.Hence, problem (1) is rewritten asmin q , 𝐳 , 𝝁 𝑓 ( q , 𝐳 , 𝝁 ) s.t. 𝐜 ( q , 𝐳 , 𝝁 ) ≤ 𝐿 ( q , 𝝁 ) = 𝑁𝐿 ( 𝐳 , 𝝁 ) = → min 𝝁 ∈ 𝑓 ( q ( 𝝁 ) , 𝐳 ( 𝝁 ) , 𝝁 ) s.t. 𝐜 ( q ( 𝝁 ) , 𝐳 ( 𝝁 ) , 𝝁 ) ≤ (2)Given at iteration 𝑖 of the NAND approach a queried parameter point 𝝁 𝑖 ∈ , 𝐑 𝐿 ( q , 𝝁 ) = and 𝐑 𝑁𝐿 ( 𝐳 , 𝝁 ) = are solved toobtain 𝐪 ( 𝝁 𝑖 ) and 𝐳 ( 𝝁 𝑖 ) , and then the objection function 𝑓 and the set of algebraic constraints represented by 𝐜 are evaluated. Thisprocess can be written as 𝝁 𝑖 → { 𝐑 𝐿 ( 𝐪 ( 𝝁 𝑖 ) , 𝝁 𝑖 ) = 𝑁𝐿 ( 𝐳 ( 𝝁 𝑖 ) , 𝝁 𝑖 ) = → { 𝐪 = 𝐪 ( 𝝁 𝑖 ) 𝐳 = 𝐳 ( 𝝁 𝑖 ) Unfortunately, the number of solutions of the PDE-based constraints generally increases when the dimension of the designparameter space is increased, which is a disadvantage of the NAND approach and a further motivation for model reducion. Boncoraglio
ET AL
The three-field framework pioneered in is adopted for modeling a coupled FSI problem such as that represented here by 𝐑 𝐿 ( q , 𝝁 ) = . In this framework, no limiting assumption is made about the behavior of the fluid subsytem whose governingequations of equilibrium are usually formulated in the Arbitrary Lagrangian Eulerian (ALE) setting and discretized by a CFDapproach. Similarly, no limiting assumption is made about the behavior of the structural subsystem whose governing equationsof equilibrium are formulated in the Lagrangian setting and discretized by a finite element (FE) method. The CFD mesh ofthe fluid subsystem is viewed as a third subsystem that is assimilated with a fictitious discrete structural system , or the FEdiscretization of a fictitious deformable continuous body – so that it may move and/or deform as needed. In either case, it isreferred to here as the pseudo-structural subsystem.Most importantly, the three-field (fluid, fluid mesh, and structure) framework couples the fluid and structural subsystemsby simply enforcing the kinematic transmission conditions expressing the slip or no-slip wall boundary conditions of the fluidsubsystem, and the force transmission conditions expressing equilibrium at the fluid/structure interface. It is overviewed belowwith emphasis on the linearization of its semi-discrete form. The three-field framework for computational FSI outlined above can be summarized by the following three equations ⎧⎪⎪⎪⎨⎪⎪⎪⎩ 𝜕 ( ) 𝜕𝑡 + ∇ ⋅ ( ( ) − 𝜕𝑥𝜕𝑡 ) − ∇ ⋅ ( ) = 0 𝜌 𝜕 𝜕 𝑡 − 𝑑𝑖𝑣 ( 𝜎 ( 𝜖 ( ))) − = 0 ̃𝜌 𝜕 𝑥𝜕 𝑡 − 𝑑𝑖𝑣 ( ̃𝐸 ∶ ̃𝜖 ( 𝑥 ) ) = 0 (3)where: the first equation is the Navier-Stokes equation written in ALE conservative form – and can be augmented with turbulencemodeling; the second equation expresses the dynamic equilibrium of the structural subsystem; the third equation models the fluidmesh as a fictitious, dynamic, deformable body; and all dependences on the parameter point 𝝁 have not been stated explicitly inorder to keep the notation as simple as possible.In the first of equations (3), denotes the vector of conservative fluid variables, 𝑥 denotes the position vector of the fluidmesh, is the determinant of the jacobian of 𝑥 with respect to the reference configuration of the fluid mesh, 𝑡 denotes time, ∇ denotes the gradient with respect to 𝑥 , ⋅ designates the standard dot product, denotes the vector of ALE convective fluxes, and denotes the vector of diffusive fluxes.In the second of equations (3), denotes the structural displacement vector, 𝑑𝑖𝑣 designates the divergence operator, 𝜌 denotesthe material density, 𝜎 denotes the stress tensor, 𝜖 denotes the strain tensor, and is the vector of body forces acting on thestructure.In the third of equations (3), the symbols ̃ and ̃ designate the fictitious aspect of the deformable body with which the fluidmesh is assimilated.As already stated, the first two equations described above are coupled by the kinematic and equilibrium transmission con-ditions which have simple expressions that can be found, for example, in : therefore, they are not repeated here. Similarly, thesecond and third sets of equations are coupled by the continuity of the structural displacement and velocity fields across thefluid/structure interface.Finally, each of equations (3) is equipped with its own boundary and initial conditions: these are not explicitly stated here forthe sake of simplicity. oncoraglio ET AL In this work, the fluid subsystem is semi-discretized by a finite volume (FV) method: however, it can be equally approximatedby a FE method. On the other hand, the structural and pseudo-structural subsystems are semi-discretized by the FE method. Inthis case, the coupled FSI system (3) is transformed into ⎧⎪⎨⎪⎩ ̇̂ ( A ( x ) w ) + F ( w , x , ̇ x ) − F ′ ( w , x ) = ̈ u + f 𝑖𝑛𝑡 ( u , ̇ u ) − f 𝑒𝑥𝑡 ( u , w ) = ̃ Kx − ̃ K 𝑐 u = (4)In the first of the above semi-discrete equations: 𝐀 ∈ ℝ 𝑁 𝑓 × 𝑁 𝑓 denotes the diagonal matrix storing the cell volumes of the FVsemi-discretization and 𝑁 𝑓 denotes the dimension of the semi-discrete fluid HDM; 𝐅 ∈ ℝ 𝑁 𝑓 and 𝐅 ′ ∈ ℝ 𝑁 𝑓 denote the vectorsof semi-discrete convective and diffusive fluxes, respectively; 𝐱 ∈ ℝ 𝑁 𝑚 denotes the semi-discrete position vector of the fluidmesh and 𝑁 𝑚 denotes the dimension of the associated HDM; 𝐰 ∈ ℝ 𝑁 𝑓 denotes the vector of semi-discrete conservative statevariables of the fluid subsystem; and ̇ () denotes the first derivative with respect to time.In the second of equations (4): 𝐌 ∈ ℝ 𝑁 𝑠 × 𝑁 𝑠 denotes the usual mass matrix and 𝑁 𝑠 denotes the dimension of the semi-discretestructural HDM; 𝐟 𝑖𝑛𝑡 ∈ ℝ 𝑁 𝑠 and 𝐟 𝑒𝑥𝑡 ∈ ℝ 𝑁 𝑠 denote the vectors of semi-discrete internal and external forces, respectively; 𝐮 ∈ ℝ 𝑁 𝑠 denotes the vector of semi-discrete structural state variables (displacements and rotations); and ̈ () denotes the secondderivative with respect to time.The third of the above equation is a quasi-static semi-discretization of the third of equations (3). In this equation, ̃ 𝐊 denotesthe pseudo-structural stiffness matrix of the fluid mesh, and ̃ K 𝑐 is a transfer matrix that describes the effect of a structural motionon the fluid mesh motion. Here, attention is focused on those FSI constraints where for all practical purposes, the fluid can be assumed to be inviscid( 𝐅 ′ = in (4)), and the structure undergoes small rotations and deformations. These assumptions are realistic for many FSIapplications such as those pertaining to flutter, stability, and control. In this case, equations (4) can be linearized about a static,fluid-structure equilibrium point ( 𝐰 , 𝐱 , 𝐮 ) characterized by ̇ 𝐰 = , ̇ 𝐱 = , and ̇ 𝐮 = , following the approach described in .Furthermore, 𝐱 can be eliminated from the three-way coupled system using the relationship 𝐱 = ̃ K −1 ̃ K 𝑐 u = 𝐓𝐮 , which followsfrom the third of equations (4). The resulting linearized system of equations can be written as { 𝐀 ̇̄ 𝐰 + 𝐇 ̄ 𝐰 + ( E + C ) 𝐓 ̇̄ u + B 𝐓 ̄ u = 𝐌 ̈̄ 𝐮 + 𝐃 ̇̄ 𝐮 + 𝐊 ̄ 𝐮 = 𝐏 ̄ 𝐰 (5)In the first of the above equations: ̄ 𝐰 , ̄ 𝐮 , and ̇̄ 𝐮 denote perturbations of 𝐰 , 𝐮 , and ̇ 𝐮 around the equilibrium point ( 𝐰 , 𝐱 , 𝐮 ) ,respectively; 𝐀 = 𝐀 ( 𝐱 ) , 𝐇 = 𝜕 F 𝜕 𝐰 ( 𝐰 , 𝐱 , ̇ 𝐱 ) , 𝐄 = 𝜕 A 𝜕 x ( 𝐱 ) 𝐰 , 𝐂 = 𝜕 F 𝜕 ̇ 𝐱 ( 𝐰 , 𝐱 , ̇ 𝐱 ) , and 𝐁 = 𝜕 F 𝜕 x ( 𝐰 , 𝐱 , ̇ 𝐱 ) . The matrices 𝐇 , 𝐄 , 𝐂 , and 𝐁 arise from the first-order expansion in Taylor series of the semi-discrete fluid equation around the aforementionedstatic equilibrium point. In this expansion, the terms 𝐀 ̇ 𝐰 , 𝐄 ̇ 𝐱 , ( 𝜕 𝐀 𝜕 𝐱 ̄ 𝐱 ) |||| ̇ 𝐰 , ( 𝜕 𝐄 𝜕 𝐰 ̄ 𝐰 ) |||| ̇ 𝐱 , and 𝐅 ( 𝐰 , 𝐱 , ̇ 𝐱 ) vanish because ̇ 𝐰 = and ̇ 𝐱 = .In the second of equations (5), 𝐃 = 𝜕 f 𝑖𝑛𝑡 𝜕 ̇ 𝐮 ( 𝐮 , ̇ 𝐮 ) , 𝐊 = 𝜕 f 𝑖𝑛𝑡 𝜕 𝐮 ( 𝐮 , ̇ 𝐮 ) − 𝜕 f 𝑒𝑥𝑡 𝜕 𝐮 ( 𝐮 , 𝐰 ) , and 𝐏 = 𝜕 f 𝑒𝑥𝑡 𝜕 𝐰 ( 𝐮 , 𝐰 ) . The specificform of this equation results from the linearization about a fluid-structure equilibrium point, and therefore a point characterizedby − 𝐌 ̈ 𝐮 − 𝐟 𝑖𝑛𝑡 ( 𝐮 , ̇ 𝐮 ) + 𝐟 𝑒𝑥𝑡 ( 𝐮 , 𝐰 ) = .In the remainder of this paper, the symbols ̄ and are dropped in order to keep the notation as simple as possible. Hence, thelinearized FSI system (5) is rewritten as follows { 𝐀 ̇ 𝐰 + 𝐇𝐰 + 𝐑 ̇ 𝐮 + 𝐆𝐮 = 𝐌 ̈ 𝐮 + 𝐃 ̇ 𝐮 + 𝐊𝐮 = 𝐏𝐰 (6)where: 𝐀 ∈ ℝ 𝑁 𝑓 × 𝑁 𝑓 , 𝐇 ∈ ℝ 𝑁 𝑓 × 𝑁 𝑓 , 𝐑 = ( 𝐄 + 𝐂 ) 𝐓 ∈ ℝ 𝑁 𝑓 × 𝑁 𝑠 , 𝐆 = 𝐁𝐓 ∈ ℝ 𝑁 𝑓 × 𝑁 𝑠 , 𝐰 ∈ ℝ 𝑁 𝑓 , 𝐮 ∈ ℝ 𝑁 𝑠 ; and 𝐌 ∈ ℝ 𝑁 𝑠 × 𝑁 𝑠 , 𝐃 ∈ ℝ 𝑁 𝑠 × 𝑁 𝑠 , 𝐊 ∈ ℝ 𝑁 𝑠 × 𝑁 𝑠 , and 𝐏 ∈ ℝ 𝑁 𝑠 × 𝑁 𝑓 . Boncoraglio
ET AL
Let A = ⎡⎢⎢⎢⎣ A 0 𝑁 𝑓 ,𝑁 𝑠 𝑁 𝑓 ,𝑁 𝑠 𝑁 𝑠 ,𝑁 𝑓 M 0 𝑁 𝑠 ,𝑁 𝑠 𝑁 𝑠 ,𝑁 𝑓 𝑁 𝑠 ,𝑁 𝑠 M ⎤⎥⎥⎥⎦ ∈ ℝ 𝑁 𝑞 × 𝑁 𝑞 , B = ⎡⎢⎢⎣ H R G − P D K0 𝑁 𝑠 ,𝑁 𝑓 − M 0 𝑁 𝑠 ,𝑁 𝑠 ⎤⎥⎥⎦ ∈ ℝ 𝑁 𝑞 × 𝑁 𝑞 , and q = ⎡⎢⎢⎣ w ̇ uu ⎤⎥⎥⎦ ∈ ℝ 𝑁 𝑞 (7)where 𝑁 𝑞 = 𝑁 𝑓 + 2 𝑁 𝑠 The linearized FSI system (6) can be rewritten in more compact form as ̃ 𝐑 𝐿 ( q ( 𝝁 ) , 𝝁 ) = A ( 𝝁 ) ̇ q ( 𝝁 ) + B ( 𝝁 ) q ( 𝝁 ) = (8)where all dependences on the parameter point 𝝁 have been recalled. Solving FSI-constrained MDAO problems in a high-dimensional design parameter space can be computationally intensive,even when the FSI constraints are linearized. Here, the computational expense associated with the enforcement of the parametric,linearized, FSI constraint ̃ 𝐑 𝐿 ( q ( 𝝁 ) , 𝝁 ) = (8) is reduced by substituting this constraint with a linear, 𝝁 -parametric, FSI PROM,and treating the dependence of this PROM on the parameter point 𝝁 using the concept of a database 𝐵 of local PROMs .When is relatively high-dimensional, the feasibility of this concept is achieved by significantly reducing the computationalcost associated with the offline construction of the database 𝐵 , as follows. First, a low-dimensional space of design parametergeneralized coordinates ⊂ ℝ 𝑛 is generated – with 𝑛 ≪ 𝑁 – using a new take on the concept of an AS . Then, issampled using a greedy procedure, and an FSI PROM is constructed at each design parameter point 𝝁 𝑖 ∈ associated witheach generalized coordinates parameter point 𝝁 𝑟 𝑖 ∈ sampled in , and stored in the database 𝐵 . Finally, at each unsampledparameter point 𝝁 ⋆ ∈ queried by the chosen optimization procedure, a local PROM counterpart of (8) is constructed inreal-time by interpolation on matrix manifolds , and applied to the real-time solution of the reduced-order counterpart of thelinearized, FSI constraint (8).The key elements of the approach outlined above for reducing the computational cost associated with the enforcement of thehigh-dimensional, linearized, FSI constraint (8) in order to accelerate the solution of the MDAO problem (2) are discussed below. In a gradient-based iterative procedure (or algorithm) for the solution of unconstrained optimization problems, the vector ofoptimization parameters 𝝁 is typically updated at each iteration 𝑘 as follows 𝝁 𝑘 +1 = 𝝁 𝑘 − 𝜂 𝑘 𝐂 𝑘 ∇ 𝑓 ( 𝝁 𝑘 ) (9)where 𝜂 𝑘 is a scalar parameter often referred to as the step size, 𝐂 𝑘 ∈ ℝ 𝑁 × 𝑁 is an algorithm-dependent matrix, and ∇ 𝑓 ( 𝝁 𝑘 ) ∈ ℝ 𝑁 is the gradient of the scalar objective function 𝑓 with respect to the vector of optimization parameters 𝝁 evaluated at 𝝁 𝑘 .From (9), it follows that 𝝁 𝑘 +1 is computed in the subspace generated by the set of gradient vectors 𝑘 = { ∇ 𝑓 ( 𝝁 ) , … , ∇ 𝑓 ( 𝝁 𝑘 ) } .Following the standard ideas of model reduction, the AS technique seeks to approximate each iterate 𝝁 𝑘 in a low-rank matrixrepresentation of 𝑘 −1 . This can be written as 𝝁 𝑘 ≈ 𝐕 𝜇 𝝁 𝑘𝑟 (10)where 𝐕 𝜇 ∈ ℝ 𝑁 × 𝑛 is a basis of dimension 𝑛 ≪ 𝑁 , and 𝝁 𝑟 ∈ ℝ 𝑛 is the vector of generalized coordinates of 𝝁 – orin AS parlance, the vector of active parameters. Specifically, 𝐕 𝜇 is considered to be a global ROB for 𝝁 , and therefore theapproximation (10) is assumed to be accurate in the entire design parameter space . Typically, this ROB is constructed usingthe proper orthogonal decomposition (POD) method of snapshots or its singular value decomposition (SVD)-based variant.This calls for sampling the gradients of 𝑓 ( ⋅ ) at various parameter points 𝝁 𝑖 ∈ , collecting the sampled values in a snapshotmatrix 𝐒 = [ ∇ 𝑓 ( 𝝁 ) , … , ∇ 𝑓 ( 𝝁 𝑁 𝐒 ) ∈ 𝑅 𝑁 × 𝑁 𝐒 ] (11) oncoraglio ET AL where 𝑁 𝐒 denotes the number of aforementioned sampled gradients, and compressing 𝐒 using SVD in order to obtain 𝐕 𝜇 . Thisprocess is summarized in Algorithm 1. Algorithm 1
Compression of the matrix of gradient snapshots.
Input:
Snapshot matrix 𝐒 ∈ ℝ 𝑁 × 𝑁 𝐒 , and tolerance 𝜀 . Output: 𝑛 , and ROB 𝐕 𝜇 of dimension 𝑛 . Compute the thin SVD of 𝐒 : 𝐒 = 𝐔𝚺𝐃 , where 𝐔 = [ 𝐮 𝐮 … 𝐮 𝑟𝑎𝑛𝑘 ] , 𝑟𝑎𝑛𝑘 denotes the rank of 𝐒 ,and 𝚺 = ⎡⎢⎢⎢⎢⎢⎣ 𝜎 𝜎 ⋱ 𝜎 𝑟𝑎𝑛𝑘 −1
00 0 … 0 𝜎 𝑟𝑎𝑛𝑘 ⎤⎥⎥⎥⎥⎥⎦ , where 𝜎 ≥ 𝜎 ≥ … ≥ 𝜎 𝑟𝑎𝑛𝑘 > 𝑛 = minimum integer for which 𝑟𝑎𝑛𝑘 ∑ 𝑗 = 𝑛 +1 𝜎 𝑗𝑟𝑎𝑛𝑘 ∑ 𝑗 =1 𝜎 𝑗 ≤ 𝜀 𝐕 𝜇 = [ 𝐮 𝐮 … 𝐮 𝑛 ] The sampling underlying the computation of the snapshot matrix (11) is usually performed using a nonadaptive algorithm.For example, the Latin Hypercube Sampling (LHS) method can be used for this purpose. Typically, the number of samples 𝑁 𝑔 is logarithmically increased with the dimension 𝑁 of the design parameter space . For example, reference recommends 𝑁 𝐒 = 𝛼𝛽 log 𝑁 (12)where 𝛼 is a free parameter referred to as the oversampling factor, and 𝛽 is another free parameter. While the values of 𝛼 and 𝛽 are problem-dependent, reference suggests a value between 2 and 10 for 𝛼 , and a value larger than 𝑛 + 1 for 𝛽 .Substituting the AS approximation (10) into (2) leads to the following MDAO problem which features a reduced number ofoptimization parameters min 𝝁 𝑟 ∈ 𝑓 ( q ( 𝐕 𝜇 𝝁 𝑟 ) , 𝐳 ( 𝐕 𝜇 𝝁 𝑟 ) , 𝐕 𝜇 𝝁 𝑟 ) s.t. 𝐜 ( q ( 𝐕 𝜇 𝝁 𝑟 ) , 𝐳 ( 𝐕 𝜇 𝝁 𝑟 ) , 𝐕 𝜇 𝝁 𝑟 ) ≤ → min 𝝁 𝑟 ∈ ̃𝑓 ( q ( 𝝁 𝑟 ) , 𝐳 ( 𝝁 𝑟 ) , 𝝁 𝑟 ) s.t. ̃ 𝐜 ( q ( 𝝁 𝑟 ) , 𝐳 ( 𝝁 𝑟 ) , 𝝁 𝑟 ) ≤ (13)The main advantage of solving the MDAO problem (13), which is based on the AS approximation (10), instead of solvingthe original MDAO problem (2) is two-fold: • First, and most importantly, mitigating the curse of dimensionality associated with the construction of a database of local,linear FSI PROMs, when is high-dimensional, since in this case parameter sampling is performed in a space of designparameter generalized coordinates of much lower dimension than the design parameter space on which it is based. • Second, lower the overall computational cost of the iterative optimization procedure by working with a smaller numberof optimization parameters, without necessarily sacrificing the optimality of the computed solution.As outlined above however, the construction of the AS suffers from two main drawbacks: • The guideline (12) for estimating the number 𝑁 𝐒 of gradients to sample at various parameter points 𝝁 𝑖 ∈ and collectin the snapshot matrix 𝐒 is an adhoc estimate. Furthermore, because 𝛼 and 𝛽 are free, problem-dependent algorithmicparameters, this estimate may still call for sampling a very large number of parameter points in , which defeats the mainpurpose of the AS concept. • Using a nonadaptive algorithm such as LHS for sampling the gradients of 𝑓 ( ⋅ ) does not necessarily lead to an appropriateAS for the solution of the MDAO problem (13), particularly when 𝑁 𝐒 itself is limited by the curse of dimensionality. For Boncoraglio
ET AL example, it is shown in Section 4 that for one of the two applications considered in this paper, the solution of the MDAOproblem (13) based on an AS constructed as outlined above is not an optimal solution of the orginal MDAO problem (2).For this reason, a new take on the concept of an AS for the solution on an MDAO problem such as (13) is proposed next.
In order to address the drawbacks of the original concept of an AS for the solution of a highly parameterized MDAO problemhighlighted in Section 3.1.1, an alternative approach for defining and constructing an AS is proposed here. This approach isbased on the observation that when solving a constrained
MDAO problem such as (2), the update of each iterate 𝝁 𝑘 performedby an optimization procedure can be written as 𝝁 𝑘 +1 = 𝝁 𝑘 + Δ 𝝁 𝑘 (14)where Δ 𝝁 𝑘 ∈ ℝ 𝑁 is determined by the optimization procedure at its 𝑘 -th iteration – and therefore is optimization-procedure-dependent. Note that expression (14) is more general than its counterpart (9). Specifically: • When solving an unconstrained optimization problem using a gradient-based procedure, Δ 𝝁 𝑘 = − 𝜂 𝑘 𝐂 𝑘 ∇ 𝑓 ( 𝝁 𝑘 ) , wherethe scalar 𝜂 𝑘 and the matrix 𝐂 𝑘 are algorithm-dependent. For example, 𝐂 𝑘 = 𝐈 for the steepest descent algorithm, 𝐂 𝑘 =∇ 𝑓 ( 𝝁 𝑘 ) for Newton’s method, and 𝐂 𝑘 = 𝐁 𝑘 for the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, where 𝐁 𝑘 isan approximation of the Hessian ∇ 𝑓 ( 𝝁 𝑘 ) . When using instead a non gradient-based optimization procedure such as agenetic algorithm, Δ 𝝁 𝑘 is determined by some stochastic process and/or a probability model. • When solving a constrained optimization problem, Δ 𝝁 𝑘 depends on whether the optimization procedure is of the interiorpoint or active set type. For example, if the optimization problem is solved using a sequential linear programming (SLP) orsequential quadratic programming (SQP) algorithm, Δ 𝝁 𝑘 is determined by the iterative solution of a sequence of problemsformulated for the purpose of satisfying the Karush-Kuhn-Tucker (KKT) conditions.Hence, depending on the chosen optimization procedure, Δ 𝝁 𝑘 can have different forms pertaining to the gradient or Hessianof 𝑓 ( ⋅ ) . For this reason, an alternative concept of an AS for the solution of a highly parameterized MDAO problem that accountsfor how Δ 𝝁 𝑘 is computed by the chosen optimization procedure is proposed here. Specifically, this concept is based on thespecific form of the increment Δ 𝝁 rather than on ∇ 𝑓 ( 𝝁 ) . Consequently, this alternative concept of an AS calls for constructingthe ROB 𝐕 𝜇 by collecting and compressing the snapshot matrix 𝕊 = [ 𝝁 , Δ 𝝁 , … , Δ 𝝁 ( 𝑁 𝕊 −2 ) ] ∈ 𝑅 𝑁 × 𝑁 𝕊 (15)where 𝝁 is the same initial guess as for the solution of the original MDAO problem (2), and 𝑁 𝕊 denotes the total number ofsnapshots collected in 𝕊 .Furthermore, since in the context of this work the main purpose of the AS is to mitigate the curse of dimensionality facedby an adaptive procedure such as a greedy procedure (for example, see ? ) for sampling the design parameter space , it is notcompelling to construct 𝕊 (15) by sampling the same parameter space using a nonadaptive algorithm! For this reason, it isproposed here to construct 𝕊 (15) by following instead the trajectory of the solution of the dramatically less computationallyintensive version of the original MDAO problem (2), where the FSI constraint (8) is dropped – that is, by solving the auxiliaryoptimization problem min 𝝁 ∈ 𝑓 ( 𝐳 ( 𝝁 ) , 𝝁 ) s.t. 𝐜 ( 𝐳 ( 𝝁 ) , 𝝁 ) ≤ (16)using the same optimization procedure and same initial guess 𝝁 as for the solution of the original MDAO problem (2), andcollecting in 𝕊 (15) the snapshots { Δ 𝝁 𝑘 } 𝑁 𝕊 −2 𝑘 =0 .If the aforementioned optimization procedure does not involve the computation of the true Hessian of the objective function,the cost of sampling Δ 𝝁 𝑘 = 𝝁 𝑘 +1 − 𝝁 𝑘 is roughly equal to that of sampling ∇ 𝑓 ( 𝝁 𝑘 ) . If on the other hand the optimizationprocedure involves the computation of ∇ 𝑓 ( 𝝁 𝑘 ) , the additional cost associated with sampling Δ 𝝁 𝑘 instead of ∇ 𝑓 ( 𝝁 𝑘 ) is thatincurred by, for example, enforcing the KKT conditions if the formulation of problem (16) incorporates constraints, or perform-ing the matrix-vector multiplication 𝐂 𝑘 ∇ 𝑓 ( 𝝁 𝑘 ) if it does not. In either case, the additional cost is affordable and offsetted by thefollowing advantages of the proposed concept of an AS for the solution of highly parameterized MDAO problems: • The sampling procedure associated with its construction is adaptive, whereas that associated with the original concept ofan AS outlined in Section 3.1.1 is typically nonadaptive (as in the case of the LHS method highlighted above): therefore oncoraglio
ET AL for the same number of samples, the sampling procedure associated with the proposed concept of an AS can be expectedto deliver a more effective ROB 𝐕 𝜇 than its counterpart associated with the original concept of an AS. • The number of snapshots needed for the construction of an effective AS of the type advocated here is automaticallydetermined by the problem-independent convergence criterion governing the solution of the auxiliary optimizationproblem (16). On the other hand, that needed for the construction of an effective AS of the original type reviewed inSection 3.1.1 is governed by the adhoc and problem-dependent criterion (12), and therefore is determined in practice bya far less convenient trial and error approach.
PMOR reduces the order (dimension, or size) of an HDM such as the linearized, FSI constraint (8) by: • Performing subspace approximations, which in this case can be written as 𝐰 = 𝐕 𝑤 𝐰 𝑟 , 𝐮 = 𝐕 𝑢 𝐮 𝑟 ⇒ 𝐪 = 𝐕 𝑞 𝐪 𝑟 = ⎡⎢⎢⎢⎣ 𝐕 𝑤 𝑁 𝑓 ,𝑛 𝑠 𝑁 𝑓 ,𝑛 𝑠 𝑁 𝑠 ,𝑛 𝑓 𝐕 𝑢 𝑁 𝑠 ,𝑛 𝑠 𝑁 𝑠 ,𝑛 𝑓 𝑁 𝑠 ,𝑛 𝑠 𝐕 𝑢 ⎤⎥⎥⎥⎦ ⎡⎢⎢⎣ w 𝑟 ̇ u 𝑟 u 𝑟 ⎤⎥⎥⎦ (17)where 𝐕 𝑤 ∈ ℝ 𝑁 𝑓 × 𝑛 𝑓 is a right fluid ROB with 𝑛 𝑓 ≪ 𝑁 𝑓 , 𝐕 𝑢 ∈ ℝ 𝑁 𝑠 × 𝑛 𝑠 is a right structural ROB with 𝑛 𝑠 ≪ 𝑁 𝑠 – andtherefore 𝐕 𝑞 ∈ ℝ ( 𝑁 𝑞 × 𝑛 𝑞 is a fluid-structure ROB with 𝑛 𝑞 ≪ 𝑁 𝑞 – 𝐰 𝑟 ∈ ℝ 𝑛 𝑓 , 𝐮 𝑟 ∈ ℝ 𝑛 𝑠 , and 𝐪 𝑟 ∈ ℝ 𝑛 𝑞 where 𝑛 𝑞 = 𝑛 𝑓 + 2 𝑛 𝑠 ≪ 𝑁 𝑞 Substituting the subspace approximation 𝐪 = 𝐕 𝑞 𝐪 𝑟 (17) in the parametric, linearized, FSI constraint (8) and using the ASapproximation (10) leads to A ( 𝝁 𝑟 ) 𝐕 𝑞 ̇ q 𝑟 ( 𝝁 𝑟 ) + B ( 𝝁 𝑟 ) 𝐕 𝑞 q 𝑟 ( 𝝁 𝑟 ) = (18)which is an overdetermined system of equations. • Performing a Galerkin or Petrov-Galerkin projection on each of the three block equations embedded in system (18), inorder to transform it into a square system. Typically, a Galerkin projection – which does not require the computation of aleft ROB – is justified for the structural subsystem. However, a Petrov-Galerkin projection may be preferred for the fluidsubsystem for numerical stability reasons . It follows that in general, the left fluid-structure ROB can be written as 𝐖 𝑞 = ⎡⎢⎢⎢⎣ 𝐖 𝑤 𝑁 𝑓 ,𝑛 𝑠 𝑁 𝑓 ,𝑛 𝑠 𝑁 𝑠 ,𝑛 𝑓 𝐕 𝑢 𝑁 𝑠 ,𝑛 𝑠 𝑁 𝑠 ,𝑛 𝑓 𝑁 𝑠 ,𝑛 𝑠 𝐕 𝑢 ⎤⎥⎥⎥⎦ (19)where 𝐖 𝑤 ∈ ℝ 𝑁 𝑓 × 𝑛 𝑓 is a left fluid ROB. Projecting next the overdetermined system (18) onto the subspace representedby the left fluid-structure ROB 𝐖 𝑞 (19) transforms this system into the square counterpart ( 𝐖 𝑇𝑞 A ( 𝝁 𝑟 ) 𝐕 𝑞 ) ̇ q 𝑟 ( 𝝁 𝑟 ) + ( 𝐖 𝑇𝑞 B ( 𝝁 𝑟 ) 𝐕 𝑞 ) q 𝑟 ( 𝝁 𝑟 ) = where the superscript 𝑇 designates the transpose of a quantity. The above reduced system is a parametric, linear, FSIPROM. Using (7), (10), (17), and (19), it can be rewritten in a more compact form as follows ̃ 𝐑 𝐿 𝑟 ( q ( 𝝁 𝑟 ) , 𝝁 𝑟 ) = A 𝑟 ( 𝝁 𝑟 ) ̇ q 𝑟 ( 𝝁 𝑟 ) + B 𝑟 ( 𝝁 𝑟 ) q 𝑟 ( 𝝁 𝑟 ) = where A 𝑟 ( 𝝁 𝑟 ) = 𝐖 𝑇𝑞 A ( 𝝁 𝑟 ) 𝐕 𝑞 = ⎡⎢⎢⎢⎣ 𝐖 𝑇𝑤 𝐀 ( 𝝁 𝑟 ) 𝐕 𝑤 𝑛 𝑓 ,𝑛 𝑠 𝑛 𝑓 ,𝑛 𝑠 𝑛 𝑠 ,𝑛 𝑓 𝐕 𝑇𝑢 𝐌 ( 𝝁 𝑟 ) 𝐕 𝑢 𝑛 𝑠 ,𝑛 𝑠 𝑛 𝑠 ,𝑛 𝑓 𝑛 𝑠 ,𝑛 𝑠 𝐕 𝑇𝑢 𝐌 ( 𝝁 𝑟 ) 𝐕 𝑢 ⎤⎥⎥⎥⎦ ∈ ℝ 𝑛 𝑞 × 𝑛 𝑞 (20)and B 𝑟 ( 𝝁 𝑟 ) = 𝐖 𝑇𝑞 B ( 𝝁 𝑟 ) 𝐕 𝑞 = ⎡⎢⎢⎣ 𝐖 𝑇𝑤 𝐇 ( 𝝁 𝑟 ) 𝐕 𝑤 𝐖 𝑇𝑤 𝐑 ( 𝝁 𝑟 ) 𝐕 𝑢 𝐖 𝑇𝑤 𝐆 ( 𝝁 𝑟 ) 𝐕 𝑢 − 𝐕 𝑇𝑢 𝐏 ( 𝝁 𝑟 ) 𝐕 𝑤 𝐕 𝑇𝑢 𝐃 ( 𝝁 𝑟 ) 𝐕 𝑢 𝐕 𝑇𝑢 𝐊 ( 𝝁 𝑟 ) 𝐕 𝑢 𝑛 𝑠 ,𝑛 𝑓 − 𝐕 𝑇𝑢 𝐌 ( 𝝁 ) 𝐕 𝑢 𝑛 𝑠 ,𝑛 𝑠 ⎤⎥⎥⎦ ∈ ℝ 𝑛 𝑞 × 𝑛 𝑞 The parametric, linearized, FSI constraint considered this work is a flutter constraint, as flutter remains one of the mostimportant considerations in aircraft design. Consequently, for a given 𝝁 𝑟 ∈ – and therefore, a given parameter point 𝝁 = 𝐕 𝝁 𝝁 𝑟 ∈ – the fluid and structural ROBs are computed as follows: Boncoraglio
ET AL • The right fluid ROB 𝐕 𝑤 is computed as described in . For this purpose, and only for this purpose, the first of the linearizedequations (5) is rewritten in the frequency domain by assuming a periodic solution of the form ̄ 𝐰 = ̄ 𝐰 𝑎 𝑒 𝐼𝜅𝑡 and a periodicexcitation of the form ̄ 𝐮 = ̄ 𝐮 𝑎 𝑒 𝐼𝜅𝑡 , where the subscript 𝑎 designates the amplitude, 𝐼 denotes the pure imaginary numbersatisfying 𝐼 = −1 , 𝜅 denotes the (aeroelastic) reduced frequency of interest, and 𝑡 denotes as usual time. Next, a few –say 𝑁 𝑚 – ground-based natural mode shapes of the structural subsystem are considered, together with a reduced frequencyband of interest. For each considered ground-based mode shape ̄ 𝐮 𝑎 𝑚 , 𝑚 = 1 , … 𝑁 𝑚 , a sweep is performed in the reducedfrequency band of interest, which for this purpose is sampled into 𝑁 𝑙 points to obtain 𝑁 𝑚 𝑁 𝑙 excitation inputs of the form ̄ 𝐮 𝑎 𝑚 𝑒 𝐼𝜅 𝑙 𝑡 . For each of these excitation inputs, a fluid amplitude solution ̄ 𝐰 𝑎 𝑚 ( 𝑘 𝑙 ) is computed by prescribing ̄ 𝐮 = ̄ 𝐮 𝑎 𝑚 𝑒 𝐼𝜅 𝑙 𝑡 inthe aforementioned frequency domain version of the first of the linearized equations (5) and solving this equation. Next, thereal and imaginary parts of each computed solution ̄ 𝐰 𝑎 𝑚 ( 𝑘 𝑙 ) are collected in a matrix of fluid solution snapshots. Finally,the 𝑁 𝑚 𝑁 𝑙 collected fluid snapshots are compressed using SVD to construct the desired fluid ROB 𝐕 𝑤 (for example, seeAlgorithm 1). In summary, this procedure for constructing 𝐕 𝑤 amounts to training the fluid ROB for different deformedconfigurations of the structural subsystem flapping at multiple frequencies, in a relevant frequency band. • The left fluid ROB 𝐖 𝑤 is computed in two steps as follows. First, this ROB is temporarily set to 𝐖 𝑤 = 𝐕 𝑤 and theeigenvalues of the reduced fluid matrix 𝐖 𝑇𝑤 𝐇𝐕 𝑤 (see (20)) are computed in real-time. If all eigenvalues of this reducedmatrix turn out to be positive, 𝐖 𝑇𝑤 𝐇𝐕 𝑤 is stable and the fluid subsystem is reduced using a Galerkin projection ( 𝐖 𝑤 = 𝐕 𝑤 ). On the other hand, if any eigenvalue of 𝐖 𝑇𝑤 𝐇𝐕 𝑤 turns out to be negative, 𝐖 𝑇𝑤 𝐇𝐕 𝑤 = 𝐕 𝑇𝑤 𝐇𝐕 𝑤 is unstable. In thiscase, 𝐖 𝑤 is updated by adding to it 1 or 2 columns that are constructed using the procedure described in , which provablyguarantees the stability of 𝐖 𝑇𝑤 𝐇𝐕 𝑤 , and the fluid is reduced by a Petrov-Galerkin projection ( 𝐖 𝑤 ≠ 𝐕 𝑤 ) . • The right structural ROB 𝐕 𝑢 is chosen as a collection of low frequency, ground-based, natural mode shapes, which isjustified by the fact that a flutter response is typically dominated by these modes. From (20), it follows that if these modeshapes are mass-orthonormalized, A 𝑟 ( 𝝁 𝑟 ) = 𝐖 𝑇𝑞 A ( 𝝁 𝑟 ) 𝐕 𝑞 = ⎡⎢⎢⎢⎣ 𝐖 𝑇𝑤 𝐀 ( 𝝁 𝑟 ) 𝐕 𝑤 𝑛 𝑓 ,𝑛 𝑠 𝑛 𝑓 ,𝑛 𝑠 𝑛 𝑠 ,𝑛 𝑓 𝐈 𝑛 𝑠 ,𝑛 𝑠 𝑛 𝑠 ,𝑛 𝑠 𝑛 𝑠 ,𝑛 𝑓 𝑛 𝑠 ,𝑛 𝑠 𝐈 𝑛 𝑠 ,𝑛 𝑠 ⎤⎥⎥⎥⎦ ∈ ℝ 𝑛 𝑞 × 𝑛 𝑞 (21)and B 𝑟 ( 𝝁 𝑟 ) = 𝐖 𝑇𝑞 B ( 𝝁 𝑟 ) 𝐕 𝑞 = ⎡⎢⎢⎣ 𝐖 𝑇𝑤 𝐇 ( 𝝁 𝑟 ) 𝐕 𝑤 𝐖 𝑇𝑤 𝐑 ( 𝝁 𝑟 ) 𝐕 𝑢 𝐖 𝑇𝑤 𝐆 ( 𝝁 𝑟 ) 𝐕 𝑢 − 𝐕 𝑇𝑢 𝐏 ( 𝝁 𝑟 ) 𝐕 𝑤 𝐕 𝑇𝑢 𝐃 ( 𝝁 𝑟 ) 𝐕 𝑢 𝛀 𝑟 ( 𝝁 𝑟 ) 𝑛 𝑠 ,𝑛 𝑓 − 𝐈 𝑛 𝑠 ,𝑛 𝑠 𝑛 𝑠 ,𝑛 𝑠 ⎤⎥⎥⎦ ∈ ℝ 𝑛 𝑞 × 𝑛 𝑞 where 𝛀 𝑟 is the diagonal matrix storing the squares of the natural angular frequencies associated with the chosen collectionof low frequency, ground-based, natural mode shapes.From (20), it also follows that the parametric, linear, FSI PROM associated with the parametric, linearized, FSI constraint (8)can be described by the following tuple of FSI reduced-order operators 𝐹 𝑆𝐼𝑟 ( 𝝁 𝑟 ) = { A 𝑟 ( 𝝁 𝑟 ) , B 𝑟 ( 𝝁 𝑟 ) } = { 𝐖 𝑇𝑞 A ( 𝝁 𝑟 ) 𝐕 𝑞 , 𝐖 𝑇𝑞 B ( 𝝁 𝑟 ) 𝐕 𝑞 } (22)However, the solution of the MDAO problem (13) using a gradient-based optimization algorithm and the parametric, PROMcounterpart (20) of the parametric, FSI constraint (8) also requires the computation of the sensitivities 𝜕 A 𝑟 𝜕𝜇 ( 𝝁 𝑟 ) and 𝜕 B 𝑟 𝜕𝜇 ( 𝝁 𝑟 ) at each point 𝝁 𝑟 𝑖 ∈ (and therefore parameter point 𝝁 𝑖 = 𝐕 𝝁 𝝁 𝑟 𝑖 ∈ ) visited by the optimization algorithm. For this reason,the extended tuple 𝐹 𝑆𝐼𝑟 ( 𝝁 𝑟 ) = { A 𝑟 ( 𝝁 𝑟 ) , B 𝑟 ( 𝝁 𝑟 ) , { 𝜕 A 𝑟 ( 𝝁 𝑟 ) 𝜕𝜇 [ 𝑗 ] , 𝜕 B 𝑟 ( 𝝁 𝑟 ) 𝜕𝜇 [ 𝑗 ] } 𝑁 𝑗 =1 } (23)where 𝝁 [ 𝑗 ] denotes the 𝑗 -th component of the parameter point 𝝁 ∈ ⊂ ℝ 𝑁 is also introduced. The linear, FSI PROM (20) is 𝝁 -parametric. As already stated, the approach chosen here for treating the parameter depen-dence of such a PROM consists in: constructing offline a database 𝐵 of local tuples of FSI reduced-order operators 𝑟 ( 𝝁 𝑟 𝑖 ) = { A 𝑟 ( 𝝁 𝑟 𝑖 ) , B 𝑟 ( 𝝁 𝑟 𝑖 ) } , 𝑖 = 1 , … , 𝑁 𝐵 ; and interpolating these online on appropriate matrix manifolds to generate in real-time the oncoraglio ET AL tuple 𝑟 ( 𝝁 𝑟 ⋆ ) at each parameter point 𝝁 ⋆ = 𝐕 𝝁 𝝁 𝑟 ⋆ queried by the optimization algorithm, but where 𝑟 ( 𝝁 𝑟 ⋆ ) is not available in 𝐵 . This approach requires first addressing two important issues: • Selecting the number 𝑁 𝐵 and locations in of the points 𝝁 𝑟 𝑖 – and therefore locations in of the parameter points 𝝁 𝑖 = 𝐕 𝝁 𝝁 𝑟 𝑖 – where to construct the local tuples of reduced-order operators defining the local, linear PROMs of interest. • Ensuring that all constructed tuples stored in 𝐵 are consistent in the sense defined in and examplified as well asexplained below.The first issue highlighted above defines the parameter sampling problem. It is common to many local and global approachesfor treating the parameter dependence of a PROM – that is, for training a PROM in (here) or in (in general). Among allelements of the approach described in this paper for reducing the computational cost associated with the enforcement of thehigh-dimensional, linearized, FSI constraint (8) in order to accelerate the solution of the MDAO problem (13) (here) or (2) (ingeneral), it is the most vulnerable to the curse of dimensionality – specifically, the dimension 𝑁 of . This issue is treated inSection 3.4.The second issue highlighted above can be illustrated as follows. Consider, for example, the linear, structural PROM associatedwith the linear, FSI PROM (20) and describable by the sub-tuple of 𝐹 𝑆𝐼𝑟 ( 𝝁 𝑟 ) (22) 𝑆𝑟 ( 𝝁 𝑟 ) = { 𝐌 𝑟 ( 𝝁 𝑟 ) , 𝐃 𝑟 ( 𝝁 𝑟 ) , 𝐊 𝑟 ( 𝝁 𝑟 ) , 𝐏 𝑟 ( 𝝁 𝑟 ) } where 𝐌 𝑟 ( 𝝁 𝑟 ) = 𝐕 𝑇𝑢 ( 𝝁 𝑟 ) 𝐌 ( 𝝁 𝑟 ) 𝐕 𝑢 ( 𝝁 𝑟 ) , 𝐃 𝑟 ( 𝝁 𝑟 ) = 𝐕 𝑇𝑢 ( 𝝁 𝑟 ) 𝐃 ( 𝝁 𝑟 ) 𝐕 𝑢 ( 𝝁 𝑟 ) 𝐊 𝑟 ( 𝝁 𝑟 ) = 𝐕 𝑇𝑢 ( 𝝁 𝑟 ) 𝐊 ( 𝝁 𝑟 ) 𝐕 𝑢 ( 𝝁 𝑟 ) , 𝐏 𝑟 ( 𝝁 𝑟 ) = 𝐕 𝑇𝑢 ( 𝝁 𝑟 ) 𝐏 ( 𝝁 𝑟 ) 𝐕 𝑤 ( 𝝁 𝑟 ) and the local aspect of the right structural ROB is emphasized by the notation 𝐕 𝑢 ( 𝝁 𝑟 ) . Furthermore, consider the case wherethis ROB is made of low frequency, ground-based, natural mode shapes of the structural subsystem, and the design parameterspace includes material properties such as the structural Young modulus and/or structural density. It is known that in this case(for example, see ), the ordering of the natural mode shapes is material property dependent. Consequently, it is possible that for 𝝁 𝑟 = 𝝁 𝑟 ( 𝝁 = 𝐕 𝝁 𝝁 𝑟 = 𝝁 ), the first 3 columns of 𝐕 𝑢 ( 𝝁 𝑟 ) correspond to – for example – the first bending mode, second-bendingmode, and first torsion mode of the structural subsystem, respectively, while for 𝝁 𝑟 = 𝝁 𝑟 ( 𝝁 = 𝐕 𝝁 𝝁 𝑟 = 𝝁 ), they correspond tothe first bending mode, first torsion mode, and second bending mode, respectively. In this case, 𝑆𝑟 ( 𝝁 𝑟 ) and 𝑆𝑟 ( 𝝁 𝑟 ) are said tobe inconsistent, because their underlying ROBs 𝐕 𝑢 ( 𝝁 𝑟 ) and 𝐕 𝑢 ( 𝝁 𝑟 ) are inconsistent in the sense that algebraic manipulationsperformed on these ROBs – and therefore, on their associated structural PROMs – such as linear combinations and interpolationsare not physically meaningful, unless the columns of 𝐕 𝑢 ( 𝝁 𝑟 ) and 𝐕 ( 𝝁 𝑟 ) are first re-ordered to be physically consistent (forexample, first bending, second bending, and then first torsion, etc., for both 𝐕 𝑢 ( 𝝁 𝑟 ) and 𝐕 ( 𝝁 𝑟 ) ). Note however that the orderingof the ROBs is only one source of inconsistency among many others.To enforce consistency between the pre-computed tuples 𝐹 𝑆𝐼𝑟 ( 𝝁 𝑟 𝑖 ) = { A 𝑟 ( 𝝁 𝑟 𝑖 ) , B 𝑟 ( 𝝁 𝑟 𝑖 ) } = { 𝐖 𝑇𝑞 ( 𝝁 𝑟 𝑖 ) A ( 𝝁 𝑟 𝑖 ) 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) , 𝐖 𝑇𝑞 ( 𝝁 𝑟 𝑖 ) B ( 𝝁 𝑟 𝑖 ) 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) } , 𝑖 = 1 , … , 𝑁 , it is first noted that any ROB such as 𝐕 𝑞 remains a ROBrepresenting the same subspace approximation after post-multiplication by an orthogonal matrix 𝐐 ∈ ℝ 𝑛 𝑞 × 𝑛 𝑞 – that is, after“rotation” of this ROB. In other words, 𝐕 𝑞 and 𝐕 𝑞 𝐐 represent the same subspace when 𝐐 is an orthogonal matrix (for example, 𝐐 𝑇 𝐐 = 𝐈 𝑛 𝑞 ,𝑛 𝑞 ). This illustrates the following facts: • There is no unique ROB associated with a given subspace approximation. • Two different ROBs such as 𝐕 𝑞 and 𝐕 𝑞 𝐐 associated with the same subspace approximation define two different generalizedcoordinates systems for performing the same approximation.It follows that one methodology for enforcing consistency between several pre-computed PROMs is to ensure that theirunderlying ROBs are constructed for the same generalized coordinates system . This methodology consists in selecting oneof the pre-computed ROBs as a reference basis, for example, 𝐕 𝑞 ( 𝝁 𝑟 𝑟𝑒𝑓 ) = 𝐕 𝑞 ( 𝝁 𝑟 ) , and finding for each other pre-computed ROB 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) and associated tuple the orthogonal transformation matrix 𝐐 𝑖 that solves the following minimization problemmin Q ∈ ( 𝑛 𝑞 ) ‖ 𝐕 𝑞 ( 𝝁 𝑟 𝑟𝑒𝑓 ) − 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) 𝐐 ‖ 𝐹 (24)where ( 𝑛 𝑞 ) denotes the set of orthogonal matrices of size 𝑛 𝑞 , and the subscript 𝐹 designates the Frobenius norm.There are two reasons for 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) to be different from 𝐕 𝑞 ( 𝝁 𝑟 𝑟𝑒𝑓 ) : 𝝁 𝑟 𝑖 ≠ 𝝁 𝑟 𝑟𝑒𝑓 ; and 𝐕 𝑞 ( 𝝁 𝑟 𝑟𝑒𝑓 ) and 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) may define two differentgeneralized coordinates systems. By solving the minimization problem (24), the orthogonal matrix 𝐐 𝑖 is determined so that 𝝁 𝑟 𝑖 ≠ 𝝁 𝑟 𝑟𝑒𝑓 is the only reason why 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) 𝐐 𝑖 ≠ 𝐕 𝑞 ( 𝝁 𝑟 𝑟𝑒𝑓 ) , and therefore the 2 ROBs are consistent in the sense defined above. Boncoraglio
ET AL
Problem (24) is known as the orthogonal Procustes problem . Most importantly, it has an analytical solution that can becomputed in real-time using Algorithm 2. This solution, 𝐐 𝑖 , is also a congruence transformation matrix that can be used toenforce the consistency of each tuple 𝐹 𝑆𝐼𝑟 ( 𝝁 𝑟 𝑖 ) , 𝑖 = 1 , … , 𝑁 𝐵 , with the reference tuple 𝐹 𝑆𝐼𝑟 ( 𝝁 𝑟 𝑟𝑒𝑓 ) by transforming this tupleinto 𝐹 𝑆𝐼𝑟 ( 𝝁 𝑟 𝑖 ) = { A 𝐶𝑟 ( 𝝁 𝑟 𝑖 ) , B 𝐶𝑟 ( 𝝁 𝑟 𝑖 ) } = { 𝐐 𝑇𝑖 𝐕 𝑇𝑞 ( 𝝁 𝑟 𝑖 ) A ( 𝝁 𝑟 𝑖 ) 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) 𝐐 𝑖 , 𝐐 𝑇𝑖 𝐕 𝑇𝑞 ( 𝝁 𝑟 𝑖 ) B 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) 𝐐 𝑖 } = { 𝐐 𝑇𝑖 A 𝑟 ( 𝝁 𝑟 𝑖 ) 𝐐 𝑖 , 𝐐 𝑇𝑖 B 𝑟 ( 𝝁 𝑟 𝑖 ) 𝐐 𝑖 } (25)Hence, after 𝐵 is constructed, and each time it is updated, its consistency is enforced by selecting arbitrarily a referencebasis and applying Algorithm 2 to perform congruence transformations on its pre-computed tuples of FSI reduced-operators. Algorithm 2
Congruence transformation of a tuple of reduced-order operators for enforcing consistency.
Input:
Reference point 𝝁 𝑟 𝑟𝑒𝑓 , associated fluid-structure ROB 𝐕 𝑞 ( 𝝁 𝑟 𝑟𝑒𝑓 ) , and another ROB 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) where 𝝁 𝑟 𝑖 ≠ 𝝁 𝑟 𝑟𝑒𝑓 . Output:
Optimal congruence transformation matrix 𝐐 𝑖 . Compute SVD of 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) 𝑇 𝐕 𝑞 ( 𝝁 𝑟 𝑟𝑒𝑓 ) = 𝐔𝚺𝐙 𝑇 Compute 𝐐 𝑖 = 𝐔𝐙 𝑇 When the design parameter space is high-dimensional, the offline pre-computation of a database of linear, FSI PROMs maynot be feasible. This, because it may require sampling a large number of parameter points 𝝁 𝑖 ∈ where to pre-computelinear FSI PROMs, in order to enable the accurate interpolation at a queried but unsampled parameter point of a tuple of FSIreduced-operators such as (22). In this case, the new take on the AS approach described in Section 3.1.2 can be used to find alower-dimensional subspace in which to solve the FSI-constrained MDAO problem (13). After this subspace is determined andrepresented by the ROB 𝐕 𝜇 as explained in Section 3.1, parameter sampling can be performed in the lower-dimensional spaceof generalized coordinates . For this purpose, two different approaches can be pursued: • A priori sampling.
This appropach samples points in the target space either randomly, or according to a nonadaptive,pre-designed scheme. Examples include the full factorial sampling, random sampling, and LHS methods. Such samplingmethods are not optimal, because they have no explicit awareness of where the interpolant will be inaccurate. For thisreason, they may require a larger than necessary number of sampled points in order to deliver the expected accuracy atinterpolation time: as such, they may lead to unaffordable databases of linear PROMs when 𝑁 is relatively high. • Adaptive sampling . This alternative approach is typically iterative. It requires the availability of an error estimator orindicator for the PROM. It samples points in regions of the target space where the current instance of the database isassessed by the error estimator/indicator to be inaccurate. Consequently, it avoids over-sampling: in this sense, it reducesthe size 𝑁 𝐵 of the needed database and makes its construction affordable. For this reason, this parameter samplingapproach is chosen here for sampling .Specifically, the indirect sampling of the design parameter space is performed here by directly sampling the parameterspace using an iterative, adaptive sampling algorithm equipped with: • A pre-selected set of candidate sample points Ξ 𝑟 that is large enough to faithfully represent . • A residual-based error indicator 𝑒 ( 𝝁 𝑟 , 𝐵 ) .This algorithm can be described as a greedy procedure that samples at each 𝑖 -th iteration the parameter point 𝝁 𝑟 𝑖 ∈ Ξ 𝑟 wheresome norm of 𝑒 ( 𝝁 𝑟 , 𝐵 ) is maximized, and builds at this point the tuple of FSI reduced-operators 𝐹 𝑆𝐼𝑟 ( 𝝁 𝑟 𝑖 ) (22). Hence, thisalgorithm leads to an incremental construction of 𝐵 that can be written as 𝑖 = { 𝑖 −1 , 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) , 𝑟 ( 𝝁 𝑟 𝑖 ) } , where 𝝁 𝑟 𝑖 = arg max 𝝁 𝑟 ∈ Ξ 𝑟 𝑒 ( 𝝁 𝑟 , 𝑖 −1 ) (26) 𝑖 is the instance of the database of linear FSI PROMs 𝐵 at iterations 𝑖 , and the storage in the database 𝐵 of eachpre-computed local ROB 𝐕 𝑞 ( 𝝁 𝑟 𝑖 ) is justified below. oncoraglio ET AL In general, the MDAO problem of interest may contain bounding constraints on the design optimization parameters that form 𝝁 – for example, 𝐜 ( q ( 𝐕 𝜇 𝝁 𝑟 ) , 𝐳 ( 𝐕 𝜇 𝝁 𝑟 ) , 𝐕 𝜇 𝝁 𝑟 ) ≤ in (13) may contain bounding constraints on 𝝁 . In this case, each sampledparameter point of the form 𝝁 𝑖 = 𝐕 𝜇 𝝁 𝑟 𝑖 must satisfy the feasibility constraint 𝝁 lb ≤ 𝐕 𝜇 𝝁 𝑟 𝑖 ≤ 𝝁 ub (27)where 𝝁 lb ∈ ℝ 𝑁 and 𝝁 ub ∈ ℝ 𝑁 are the lower and upper bounds on 𝝁 , respectively. The case where no upper or lower boundingconstraint on 𝝁 is specified can be simply accommodated by setting the entries of 𝝁 ub and/or 𝝁 lb , as needed, to ±∞ (very largenumbers in practice).It follows that all candidate parameter points 𝝁 𝑟 𝑖 ∈ Ξ 𝑟 must also satisfy the constraint (27). This can be achieved by con-structing Ξ 𝑟 using Algorithm 3 which samples points in that are feasible in by solving the feasibility problem (28). It can beinitialized using a single design parameter point in the center of the AS, or a small number of design parameter points randomlychosen in the AS. Algorithm 3
Construction of a set of candidate parameter points Ξ 𝑟 using a feasible active subspace algorithm. Input:
Bounding vectors 𝝁 lb and 𝝁 ub , scalar constant values 𝑐 and 𝑐 , cardinal 𝑁 of Ξ 𝑟 , and AS ROB 𝐕 𝜇 . Output:
Set of feasible candidate points in the AS, Ξ 𝑟 , and set of feasible points in the full space , Ξ . Sample 𝑁 points { 𝝁 𝑟 𝑖 ∈ } 𝑁𝑖 =1 , using a uniform tensor product design of experiment 𝑐 𝑑𝑖𝑎𝑔 ( 𝑠𝑖𝑔𝑛 ( 𝐕 𝜇 ) 𝑇 𝐕 𝜇 ) ≤ 𝝁 𝑟 𝑖 ≤ 𝑐 𝑑𝑖𝑎𝑔 ( 𝑠𝑖𝑔𝑛 ( 𝐕 𝜇 ) 𝑇 𝐕 𝜇 ) , ∀ 𝑖 = 1 , … , 𝑁 or 𝑐 ≤ 𝝁 𝑟 𝑖 ≤ 𝑐 , ∀ 𝑖 = 1 , … , 𝑁 for 𝑖 = 1 ∶ 𝑁 do Solve for 𝝁 𝑟 𝑖 the feasibility problem 𝝁 𝑖⋆ = min 𝝁 ∈ ℝ 𝑁𝜇 𝑇 𝝁 s.t. 𝐕 𝑇𝜇 𝝁 = 𝝁 𝑟 𝑖 𝝁 lb ≤ 𝝁 ≤ 𝝁 ub (28) If the feasibility problem has a nontrivial solution, update Ξ 𝑟 and ΞΞ 𝑟 = { Ξ 𝑟 ∪ 𝝁 𝑟 𝑖 } , Ξ = { Ξ ∪ 𝝁 ⋆𝑖 } end for As for the residual-based error indicator, which is needed for sampling at each iteration 𝑖 the parameter point 𝝁 𝑟 𝑖 =arg max 𝝁 𝑟 ∈ Ξ 𝑟 𝑒 ( 𝝁 𝑟 , 𝑖 −1 ) , a practical and cost-effective choice in the context of the solution of the MDAO problem (13) is 𝑒 ( 𝝁 𝑟 , 𝐵 ) = ‖‖‖‖ ̃ 𝐑 𝐿 ( ̃ 𝐕 𝑞 ( 𝐕 𝜇 𝝁 𝑟 ) 𝐪 𝑟 ( 𝝁 𝑟 ) )‖‖‖‖ (29)Here, ̃ 𝐑 𝐿 is the residual (8) associated with the linearized FSI system (6), 𝐕 𝑞 ( 𝐕 𝜇 𝝁 𝑟 ) 𝐪 𝑟 ( 𝝁 𝑟 ) is the reconstructed, high-dimensional, linearized FSI solution at the reconstructed design parameter point 𝐕 𝜇 𝝁 𝑟 , and the tilde notation designates that thisROB is not available in the content of the database 𝐵 𝑖 −1 that is available at the beginning of the 𝑖 -th iteration of the greedysampling procedure. Since ̃ 𝐕 𝑞 ( 𝐕 𝜇 𝝁 𝑟 𝑖 ) is expected to satisfy the same orthogonality procedure satisfied by every previouslysampled ROB 𝐕 𝑞 ( 𝐕 𝜇 𝝁 𝑟 𝑗 ) , 𝑗 = 1 , … , 𝑖 −1 , it can be rigorously approximated by interpolation on a Grassmann manifold 𝐺 𝑁 𝑞 ,𝑛 𝑞 .Alternatively, it can be approximated faster using constant extrapolation from the nearest, previously sampled parameter point 𝝁 𝑟 𝑗 . In either case, this justifies the storage in 𝐵 of the pre-computed ROBs noted in (26).The evaluation of ̃ 𝐑 𝐿 at the reconstructed FSI solution ̃ 𝐕 𝑞 ( 𝐕 𝜇 𝝁 𝑟 ) 𝐪 𝑟 ( 𝝁 𝑟 ) is typically inexpensive compared to the solutionof the problem ̃ 𝐑 𝐿 ( q ( 𝝁 ) , 𝝁 ) = A ( 𝝁 ) ̇ q ( 𝝁 ) + B ( 𝝁 ) q ( 𝝁 ) = . Nevertheless, since Ξ 𝑟 must be large enough to faithfully represent , the evaluation of the error indicator (29) at a large number of candidate parameter points may still be overwhelming for someapplications. In this case, the error indicator (29) can be applied only to a subset of Ξ 𝑟 that is re-selected at the beginning ofeach iteration of the greedy sampling procedure. Boncoraglio
ET AL
After a database 𝐵 of consistent tuples of FSI reduced-operators (25) – which is simply referred to in the remainder of this paperas a consistent database 𝐵 – is constructed, the linear FSI problem (8) can be solved in real-time at each unsampled parameterpoint 𝝁 𝑟 ⋆ ∈ (and therefore 𝝁 ⋆ ∈ ) queried by the optimization procedure chosen for solving the MDAO problem (13) intwo steps as follows: • Interpolate in real-time the content of the consistent database 𝐵 on appropriate matrix manifolds to compute a linear,FSI PROM of the form given in (22) at the queried parameter point 𝝁 ⋆ = 𝐕 𝜇 𝝁 𝑟 ⋆ . For this purpose: note that the contentof 𝐵 , which consists of tuples of FSI reduced-operators of the form given in (25) and ROBs (see (26)), can be organizedblock-by-block according to (21) and (17), respectively; and therefore, perform matrix interpolation block-by-block on anappropriate matrix manifold . • Apply the interpolated tuple 𝐹 𝑆𝐼𝑟 ( 𝝁 ⋆ ) to solve in real-time the linearized, FSI constraint (8).Here, the appropriate matrix manifold can be chosen as follows: • The manifold of invertible real matrices of size 𝑛 𝑓 , 𝐺𝐿 ( 𝑛 𝑓 , ℝ ) , for the submatrices of the form 𝐀 𝑟 ( 𝝁 𝑟 ) = 𝐖 𝑇𝑤 𝐀 ( 𝝁 𝑟 ) 𝐕 𝑤 ∈ ℝ 𝑛 𝑓 × 𝑛 𝑓 of a matrix of the form A 𝑟 ( 𝝁 𝑟 ) (21), and the submatrices of the form 𝐇 𝑟 ( 𝝁 𝑟 ) = 𝐖 𝑇𝑤 𝐇 ( 𝝁 𝑟 ) 𝐕 𝑤 ∈ ℝ 𝑛 𝑓 × 𝑛 𝑓 of a matrixof the form B 𝑟 ( 𝝁 𝑟 ) (21). • The manifold of symmetric positive definite (SPD) matrices of size 𝑛 𝑠 , 𝑆𝑃 𝐷 ( 𝑛 𝑠 ) , for the submatrices of the form 𝐃 𝑟 ( 𝝁 𝑟 ) = 𝐕 𝑇𝑢 𝐃 ( 𝝁 𝑟 ) 𝐕 𝑢 ∈ ℝ 𝑛 𝑠 × 𝑛 𝑠 and 𝐊 𝑟 ( 𝝁 𝑟 ) = 𝐕 𝑇𝑢 𝐊 ( 𝝁 𝑟 ) 𝐕 𝑢 = 𝛀 𝑟 ( 𝝁 𝑟 ) ∈ ℝ 𝑛 𝑠 × 𝑛 𝑠 of a matrix of the form B 𝑟 ( 𝝁 𝑟 ) (21). • The manifold of 𝑛 𝑓 × 𝑛 𝑠 real matrices, ℝ 𝑛 𝑓 × 𝑛 𝑠 , for the submatrices of the form 𝐑 𝑟 ( 𝝁 𝑟 ) = 𝐖 𝑇𝑤 𝐑 ( 𝝁 𝑟 ) 𝐕 𝑢 ∈ ℝ 𝑛 𝑓 × 𝑛 𝑠 and thesubmatrices of the form 𝐆 𝑟 ( 𝝁 𝑟 ) = 𝐖 𝑇𝑤 𝐆 ( 𝝁 𝑟 ) 𝐕 𝑢 ∈ ℝ 𝑛 𝑓 × 𝑛 𝑠 of a matrix of the form B 𝑟 ( 𝝁 𝑟 ) (21). • The manifold of 𝑛 𝑠 × 𝑛 𝑓 real matrices, ℝ 𝑛 𝑠 × 𝑛 𝑓 , for the submatrices of the form 𝐏 𝑟 ( 𝝁 𝑟 ) = 𝐕 𝑇𝑢 𝐏 ( 𝝁 𝑟 ) 𝐕 𝑤 ∈ ℝ 𝑛 𝑠 × 𝑛 𝑓 of a matrixof the form B 𝑟 ( 𝝁 𝑟 ) (21).As for the interpolation procedure, it can be summarized in 3 steps as follows (see ):1. Choose a reference block 𝐗 𝑟𝑒𝑓 = 𝐗 ( 𝝁 𝑟 𝑟𝑒𝑓 ) among the matrix blocks 𝐗 𝑗 = 𝐗 ( 𝝁 𝑟 𝑗 ) , 𝑗 = 1 , … , 𝑁 𝐵 , to be interpolated onthe matrix manifold .2. Apply the logarithm map to transfer each matrix block 𝐗 𝑗 to the tangent space to at 𝐗 𝑟𝑒𝑓 , 𝐗 – that is, compute thelogarithm of each matrix 𝐗 𝑗 , 𝚪 𝑗 = 𝚪 ( 𝝁 𝑟 𝑗 ) = Log 𝐗 𝑟𝑒𝑓 𝐗 𝑗 , 𝑗 = 1 , … , 𝑁 𝐵 .3. In 𝐗 , interpolate in the AS the computed matrix logarithms 𝚪 𝑗 = 𝚪 ( 𝝁 𝑟 𝑗 ) , 𝑗 = 1 , … , 𝑁 𝐵 , entry-by-entry, usingweighted sums of radial basis functions applied directly to the parameter points 𝝁 𝑟 𝑗 , 𝑗 = 1 , … , 𝑁 𝐵 . Let 𝚪 ⋆ = 𝚪 ( 𝝁 𝑟 ⋆ ) denote the result of this interpolation.4. Apply the exponential map to 𝚪 ⋆ in order to bring this matrix to the manifold – that is, compute the exponential of thematrix 𝚪 ⋆ , 𝐗 ⋆ = 𝐗 ( 𝝁 𝑟 ⋆ ) = Exp 𝐗 𝑟𝑒𝑓 𝚪 ⋆ .A related interpolation method proposed in can be used to interpolate the sensitivity matrices defining the extendedtuples (23) without having to pre-compute and store in 𝐵 the additional quantities { 𝜕 A 𝑟 ( 𝝁 𝑟 ) 𝜕𝜇 [ 𝑗 ] } 𝑁 𝑗 =1 and { 𝜕 B 𝑟 ( 𝝁 𝑟 ) 𝜕𝜇 [ 𝑗 ] } 𝑁 𝑗 =1 char-acterizing such extended tuples – which are needed for the fast solution of the MDAO problem (13) using a gradient-basedoptimization method.Table 1 gives the expressions of the logarithm and exponential maps for the various matrix manifolds identified above.Algorithm 4 summarizes the method of interpolation in the AS on a matrix manifold. In this section, the computational framework described in this paper – which intertwines the proposed alternative AS conceptwith adaptive parameter sampling, PMOR, and the concept of a database of linear PROMs equipped with interpolation on matrix oncoraglio
ET AL TABLE 1
Logarithm and exponential maps for some matrix manifolds of interest.Manifold 𝑚 × 𝑛 𝐺𝐿 ( 𝑛, ℝ ) 𝑆𝑃 𝐷 ( 𝑛 ) Log 𝐗 ( 𝐘 ) 𝐘 − 𝐗 log ( 𝐘𝐗 −1 ) log ( 𝐗 −1∕2 𝐘𝐗 −1∕2 ) Exp 𝐗 ( 𝚪 ) 𝐗 + 𝚪 exp ( 𝚪 ) 𝐗 𝐗 exp ( 𝚪 ) 𝐗 Algorithm 4
Interpolation in the AS on a matrix manifold . Input: 𝑁 PROM matrices 𝐗 𝑗 = 𝐗 ( 𝝁 𝑟 𝑗 ) belonging to , 𝝁 𝑟 𝑗 ∈ , 𝑗 = 1 , … , 𝑁 𝐵 , and queried but unsampled parameterpoint 𝝁 𝑟 ⋆ ∈ . Output:
Interpolated PROM matrix 𝐗 ⋆ = 𝐗 ( 𝝁 𝑟 ⋆ ) . Choose a reference PROM matrix 𝐗 𝑟𝑒𝑓 = 𝐗 ( 𝝁 𝑟 𝑟𝑒𝑓 ) for 𝑗 = 1 , … , 𝑁 do Compute 𝚪 𝑗 = Log 𝐗 𝑟𝑒𝑓 𝐗 𝑗 end for Interpolate in the AS the matrices 𝚪 𝑗 , 𝑗 = 1 , … , 𝑁 𝐵 , using weighted sums of radial basis functions to obtain 𝚪 ⋆ = 𝚪 ( 𝝁 𝑟 ⋆ ) . Compute 𝐗 ⋆ = 𝐗 ( 𝝁 𝑟 ⋆ ) = Exp 𝐗 𝑟𝑒𝑓 𝚪 ⋆ .manifolds – is applied to the solution of MDAO problems with flutter constraints. For this purpose, two different aeroelasticsystems are considered: a flexible configuration of NASA’s CRM; and NASA’s ARW-2. In both cases, air is modeled as a perfectgas, and the flow is assumed to be inviscid.All HDM, PMOR, and PROM computations discussed below are performed in double precision floating point arithmetic ona Linux cluster, using AERO Suite . Specifically, all fluid and structural ROBs and all aeroelastic PROMs are computed asdescribed in Section 3.2.In , the authors considered a computational framework for linear PMOR that is related to that presented in this paper, butwithout any AS concept. They applied it to MDAO problems similar to those considered here, but with fewer optimizationparameters. They discussed in detail the computational benefits of the framework and reported three orders of magnitude speedupfactors due to PMOR. Because the main difference between the computational framework described in and its counterpartpresented in this paper is the incorporation in the process of the AS concept presented in Section 3.1.2 in order to enablecomputational feasibility for a larger number of optimization parameters, only the speedup factors due to the considered ASconcept are reported herein. These speedup factors can be simply multiplied by those due to PMOR in the absence of an ASto assess the potential of the overall computational framework described in this paper for accelerating the solution of MDAOproblems with linearized FSI constraints. In aeronautics, flutter constraints are typically formulated in terms of a lower bound on the modal damping ratios of the lineardynamical system (18), or more practically here, its reduced-order counterpart (20). The latter linear, aeroelastic PROM can bere-written as ̇ q 𝑟 ( 𝝁 𝑟 ) + N 𝑟 ( 𝝁 𝑟 ) q 𝑟 ( 𝝁 𝑟 ) = , where N 𝑟 ( 𝝁 𝑟 ) = A −1 𝑟 ( 𝝁 𝑟 ) B 𝑟 ( ̃ 𝝁 𝑟 ) (30)and A 𝑟 ( 𝝁 𝑟 ) and B 𝑟 ( 𝝁 𝑟 ) are given in (21).The eigenvalue problem associated with (30) is N 𝑟 ( 𝝁 𝑟 ) ̂ 𝐪 𝑗 ( 𝝁 𝑟 ) = 𝜆 𝑗 ( 𝝁 𝑟 ) ̂ 𝐪 𝑗 ( 𝝁 𝑟 ) , 𝑗 = 1 , … , 𝑛 𝑞 (31)where 𝜆 𝑗 ( 𝝁 𝑟 ) ∈ ℂ denotes the 𝑗 -th eigenvector and ̂ 𝐪 𝑗 ( 𝝁 𝑟 ) ∈ ℂ 𝑛 𝑞 denotes its associated eigenvector (note that if an eigenpair ( 𝜆 𝑗 ( 𝝁 𝑟 ) , ̂ 𝐪 𝑗 ( 𝝁 𝑟 ) ) is complex-valued, its complex conjugate is also an eigenpair of (31)). Boncoraglio
ET AL
Then, the modal damping ratios of the linear, aeroelastic PROM (20) are given by 𝜁 𝑗 = − 𝜆 𝑅𝑗 √( 𝜆 𝑅𝑗 ) + ( 𝜆 𝐼𝑗 ) , 𝑗 = 1 , … , 𝑛 𝑞 where 𝜆 𝑅𝑗 and 𝜆 𝐼𝑗 denote the real and imaginary parts of the complex eigenvalue 𝜆 𝑗 , respectively, and the flutter constraints aretypically formulated as 𝜁 𝑗 ( 𝝁 𝑟 ) ≥ 𝜁 lb , 𝑗 = 1 , … , 𝑛 𝑞 where 𝜁 lb is a regulation-specified lower bound. First, the following sixth-dimensional MDAO problem for a flexible configuration of NASA’s CRM, whose geometrical andmaterial descriptions are summarized in Table 2, is consideredmaximize 𝝁 ∈ ⊂ ℝ 𝐿 ( 𝝁 ) 𝐷 ( 𝝁 ) subject to 𝜎 VM ( 𝝁 ) ≤ 𝜎 ub 𝝁 lb ≤ 𝝁 ≤ 𝝁 ub 𝜻 ( 𝝁 ) ≥ 𝜁 lb (32)where: • 𝐿 ( 𝝁 ) and 𝐷 ( 𝝁 ) denote the parametric lift and drag, respectively, associated with the following flight conditions: – Level flight at the altitude of ℎ = 10 , m, where the free-stream pressure and density are 𝑃 ∞ = 23835 . Pa and 𝜌 ∞ = 0 . kg ∕ m , respectively. – Free-stream Mach number 𝑀 ∞ = 0 . , angle of attack AoA = 1 . ◦ , and angle of sideslip AoS = 0 ◦ . • 𝜎 VM ( 𝝁 ) denotes the parametric von Mises stress at any point of the structural model of the CRM, and its upper bound isset to 𝜎 ub = 1 .
331 × 10 Pa. • 𝝁 lb and 𝝁 ub are lower and upper bounds for the vector of optimization parameters, respectively, and define box constraintsfor 𝝁 . • 𝜻 ( 𝝁 ) is the parametric vector of damping ratios 𝜁 𝑗 , 𝜁 lb > is its regulation-specified lower bound, ∈ ℝ and ∈ ℝ are vectors of ones and zeros, respectively, the lower bound coefficient 𝜁 lb is set to 𝜁 lb = 4 .
75 × 10 −3 , and 𝜻 ( 𝝁 ) ≥ 𝜁 lb defines a flutter avoidance constraint.Because the focus is set here on level flight at a zero angle of sideslip, the considered CRM configuration is a symmetric one.Hence, only half of its geometry is modeled.The computational fluid domain is chosen to be a hemisphere with a symmetry plane along the middle of the fuselage. Itis discretized by a three-dimensional (3D), unstructured, body-fitted CFD mesh with , grid points (see Figure 1). Theparametric lift and drag are obtained by postprocessing the computed flow solution.The parametric von Mises stress field is computed using the FE structural representation of half of the flexible CRM configu-ration shown in Figure 2, which has 58,888 degrees of freedom (dofs). In this representation, the FE structural modeling of thewing is that described in ? . The other components of the FE structural representation shown in Figure 2 were developed by theauthors and incorporated in the final FE structural model to facilitate two-way coupled aeroelastic computations.The flutter avoidance constraint is the most CPU intensive constraint in the formulation of the MDAO problem (32). It isapproximated using a database of linear FSI PROMs.The six design optimization parameters stored in 𝝁 are organized in two groups: three optimization parameters { 𝝁 [1] , 𝝁 [2] , 𝝁 [3]} that are shape design parameters, and three optimization parameters { 𝝁 [4] , 𝝁 [5] , 𝝁 [6]} that are structuraldesign parameters.The three shape design parameters are chosen as follows: 𝝁 [1] and 𝝁 [2] to control the dihedral angle of the wing; and 𝝁 [3] to control its sweep angle. All shape changes dictated by the optimization procedures are effected using the open-source 3D oncoraglio ET AL TABLE 2
Geometrical and material descriptions of a flexible CRM configuration.Item ValueGeometryWing span 58.76 mroot chord 11.92 mtip chord 2.736 mMaterialsSkin, except flaps(various composite materials)Stiffeners (aluminum) E . Pa 𝜌 .
78 × 10 kg ∕ m 𝜈 FIGURE 1
CRM: hemispherical computational fluid domain with a symmetry plane (left); and partial view of its discretizationby an inviscid, unstructured, body-fitted mesh (right).computer graphics software Blender . Figure 3 illustrates the concept of lattice-based deformations of Blender to effect shapechanges , and Figure 4 graphically depicts the effects on the shape of the wing of the specified lower and upper bounds for thedihedral and sweep angles.The three structural design parameters are chosen as three thickness increments for three different groups of stiffeners (seeFigure 5). Each increment is defined as a percentage of the initial thickness to which it is applied, and is constrained to have amagnitude less than 10% (box constraint). Boncoraglio
ET AL
FIGURE 2
CRM – FE structural respresentations of: wing (top); fuselage and control surfaces (middle); and skin (bottom). oncoraglio
ET AL FIGURE 3
CRM – parameters controling: dihedral angle of the wing (left); and sweep angle of the wing (right).
FIGURE 4
CRM – wing shapes associated with: upper (green) and lower (red) bounds for the dihedral angle (left); and upper(green) and lower (red) bounds for the sweep angle (right).
FIGURE 5
CRM – organization of the stiffeners of the wing in three different groups represented by three different colors. Boncoraglio
ET AL
Three databases of FSI PROMs are built for the solution of the MDAO problem (32): • 𝐵 , which is constructed by directly sampling the design parameter space using a greedy procedure similar to thatdescribed in Section 3.4, but applied directly to . The greedy procedure samples in this case 𝑁 𝐵 = 119 feasible designparameter points. • 𝐵 , which is constructed by building first the lower-dimensional space of design parameter generalized coordinates cl using the classical AS method described in Section 3.1.1 and its associated empirical formula (12), then indirectlysampling by directly sampling . In this case, the free parameters of formula (12) are chosen as 𝛼 = 10 and 𝛽 = 6 ,which yields 𝑁 𝐒 = 47 . Hence, 47 gradients of the objective function are computed at 47 feasible parameter points that arerandomly sampled in using the LHS method. The computation of each gradient consumes 0.3 hour wall-clock time ona Linux cluster with 240 processors. Next, the computed snapshots are compressed into a ROB 𝐕 𝜇 of dimension 𝑛 cl = 3 (see below). • 𝐵 , which is constructed by building first al using the alternative AS method described in Section 3.1.2, then indirectlysampling by directly sampling , using the greedy procedure described in Section 3.4. In this case, the snapshots of theform Δ 𝝁 are adaptively computed by solving the auxiliary optimization problemmaximize 𝝁 ∈ ⊂ ℝ 𝐿 ( 𝝁 ) 𝐷 ( 𝝁 ) subject to 𝜎 VM ( 𝝁 ) ≤ 𝜎 ub 𝝁 lb ≤ 𝝁 ≤ 𝝁 ub (33)which is obtained by removing from the MDAO problem (32) the CPU intensive FSI constraint. This auxiliary problemis solved in 16 iterations using the sequential least-squares programming method SLSQP, which uses the Han-Powellquasi-Newton method with a BFGS update and an 𝐿 test function in the step-length algorithm. Hence, the solution of theoptimization problem (33) generates 17 snapshots. The computation of each of the snapshots consumes about 0.34 hourwall-clock time – roughly the same amount of time as the computation of a gradient snapshot in the classical AS method– on the same Linux cluster. Next, the snapshots are compressed into a ROB 𝐕 𝜇 of dimension 𝑛 al = 3 (see below).Figure 6 plots the distributions of the singular values of both matrices 𝐒 (11) and 𝕊 (15). For the snapshot matrix 𝐒 associatedwith the classical AS method, there is no clear cutoff singular value for constructing a ROB 𝐕 𝜇 , as the singular values are shownto continuously decay. For the snapshot matrix 𝕊 associated with the alternative AS method however, Figure 6-right shows thatthe last three singular values are much smaller than the first three ones. This suggests constructing a ROB 𝐕 𝜇 as the set of thefirst three columns of the matrix 𝐔 arising from the SVD of 𝕊 – that is, constructing a ROB 𝐕 𝜇 of dimension 𝑛 al = 3 . Hence,the dimension of 𝐒 is also set to 𝑛 cl = 3 , in order to enable various meaningful performance comparisons.The concept of indirectly sampling by directly sampling the lower-dimensional space using the feasible AS algorithm(Algorithm 3) is illustrated in Figure 7, for the case of the problem considered herein and 𝐵 . First, a set of 1,000 preliminarypoints is generated in al using a uniform tensor product in this parameter space. Next, Algorithm 3 is applied to select amongthese points a set Ξ 𝑟 of 323 feasible candidate points (see Figure 7-left), and determine the associated set of feasible points Ξ ∈ (see Figure 7–right). Then, the greedy algorithm described in Section 3.4 is directly applied to Ξ 𝑟 – and indirectly to Ξ – to perform the final parameter sampling and adaptively construct the database of FSI PROMs 𝐵 . Note that the feasibleparameter points in Ξ are not uniformly spread in , but only in some regions of this design parameter space that depend onthe ROB 𝐕 𝜇 (see Figure 7). The plot in the 𝑖 -th diagonal of Figure 7-left shows the distribution of the 𝑖 -th component 𝝁 𝑟 [ 𝑖 ] ofthe candidate parameter points in Ξ 𝑟 . Its counterpart in Figure 7-right shows the distribution of the 𝑖 -th component 𝝁 [ 𝑖 ] of thecorresponding candidate parameter points in Ξ . The plot in the ( 𝑖, 𝑗 ) off-diagonal of Figure 7-right shows the distribution of thecontribution of the 𝑗 -th generalized coordinate 𝝁 𝑟 [ 𝑗 ] to the 𝑖 -th component 𝝁 [ 𝑖 ] of the candidate parameter point 𝝁 in Ξ .In each constructed database, each pre-computed linear FSI PROM is built as described in Section 3.2, with 𝑛 𝑓 = 100 and 𝑛 𝑠 = 6 . Hence, each pre-computed linear FSI PROM has the dimension 𝑛 𝑞 = 𝑛 𝑓 + 2 𝑛 𝑠 = 112 .Table (3) reports on the performance of the offline phase for problem (32) of the computational framework for MDAOdescribed in this paper, depending on whether an AS is used or not, and which method is chosen to construct it. The followingobservations are noteworthy: oncoraglio ET AL FIGURE 6
Singular values of: the snapshot matrix 𝐒 computed using the classical AS method (left); and the counterpart matrix 𝕊 computed using the alternative AS method (right). FIGURE 7
Set Ξ 𝑟 ∈ al containing 323 feasible candidate points determined by the feasible AS algorithm as well as thealternative approach for constructing an AS (left), and associated set Ξ ∈ containing the same number of feasible candidatepoints (right). • Even when using the same greedy procedure equipped with the same convergence tolerance, the database of linear FSIPROMs 𝐵 ends up being less populated than 𝐵 . This is because in this case, 𝐵 is constructed to be as accurate as 𝐵 , but in a smaller parameter space. Similarly, 𝐵 is constructed to be as accurate as 𝐵 – albeit using a differentconstruction procedure – but in a smaller parameter space. • The computational cost of the greedy procedure is roughly proportional to the number of sampled parameter points, whichitself appears to be proportional to the dimension of the sampled parameter space. • Furthermore, even after accounting for the computational overhead associated with constructing a representative basis 𝐕 𝜇 , an AS speeds up the construction of a database of PROMs. It remains to asses next however, the performance of eachdatabase constructed using an AS relative to the performance of 𝐵 , which is built without an AS. Boncoraglio
ET AL
TABLE 3
CRM: performance results for the offline phaseMethod 𝑁 𝑁 𝐒 or 𝑁 𝕊 Wall-clock time 𝑁 𝐵 Wall-clock time Wall-clock time Speed-upcomputing 𝐕 𝜇 greedy procedure total offline relative to no ASw/o AS 6 N.A. N.A. 119 120.87 hrs 120.87 hrs N.A.Classical AS 3 47 14.06 hrs 61 61.95 hrs 76.01 hrs 1.59Alternative AS 3 17 5.83 hrs 62 62.97 hrs 68.80 hrs 1.76 Figure 8 and Figure 9 display the convergence histories of the objective function and constraints of the MDAO problem (32),respectively. They highlight that, at least for this CRM MDAO problem, the classical AS method leads to an ineffective optimiza-tion strategy. Indeed, both figures reveal that in the case of the classical AS method, the PROM-based computational frameworkdescribed in this paper fails to improve the objective function beyond its initial value. On the other hand, the same computationalframework equipped with the alternative AS method leads to almost the same optimal objective function and constraints as itscounterpart framework equipped with direct sampling of the entire design parameter space, and in almost the same number ofiterations. In particular, note that during the first few iterations, the alternative AS method leads to the same values of the opti-mized objective function and constraints as directly sampling the full space . This is because during the first few iterations,the flutter constraint is inactive, the auxiliary MDAO problem (33) is in this case identical to the original MDAO problem (32),and thereore the basis 𝐕 𝜇 based on the snapshot matrix 𝕊 (15) is optimal.Table 4) shows the optimal solution of the MDAO problem (32) obtained using the three different PROM databases 𝐵 , 𝐵 , and 𝐵 . The reader can observe that the optimal solution obtained using the PROM database 𝐵 – that is, the classicalAS method for sampling – is far from its counterpart obtained using the PROM database 𝐵 – that is, by directly sampling thedesign parameter space . On the other hand, the local optimal solution obtained using the PROM database 𝐵 – that is, theaternative AS method for sampling – is relatively close to its counterpart obtained using the database of linear FSI PROMs 𝐵 . FIGURE 8
CRM: convergence histories of the objective function. oncoraglio
ET AL FIGURE 9
CRM: convergence histories of the constraints.
TABLE 4
CRM: locally optimal computed solutions 𝝁 [1] 𝝁 [2] 𝝁 [3] 𝝁 [4] 𝝁 [5] 𝝁 [6] Initial configuration 0.0 0.0 0.0 0.0 0.0 0.0Specified lower bound -0.5 -0.5 -0.5 -0.05 -0.1 -0.1Specified upper bound 2.0 2.0 2.0 0.1 0.1 0.1Optimal configuration using 𝐵 -0.50 -0.16 1.47 0.07 0.05 0.06Optimal configuration using 𝐵 . −4 −1 . −5 −2 . −3 −4 . −2 −3 . −3 . −3 Optimal configuration using 𝐵 -0.32 -0.15 1.35 0.10 0.05 0.08In summary, for this MDAO problem with six design optimization parmeters only, the alternative subspace method speedsup the offline computations by almost a factor two, while leading to a reasonably close local optimum solution.Next, another MDAO problem with a larger number of design optimization parameters is considered. Next, the following fifteenth-dimensional MDAO problem for NASA’s ARW-2, whose geometrical and material descriptionsare summarized in Table 5, is considered maximize 𝝁 ∈ ⊂ ℝ 𝐿 ( 𝝁 ) 𝐷 ( 𝝁 ) subject to 𝑊 ( 𝝁 ) ≤ 𝑊 ub 𝜎 VM ( 𝝁 ) ≤ 𝜎 ub 𝝁 lb ≤ 𝝁 ≤ 𝝁 ub 𝜻 ( 𝝁 ) ≥ 𝜁 lb (34)where: • 𝐿 ( 𝝁 ) and 𝐷 ( 𝝁 ) denote as before the parametric lift and drag, respectively, associated here with the following flightconditions: Boncoraglio
ET AL – Level flight at the altitude of ℎ = 4 , ft, where the free-stream pressure and density are 𝑃 ∞ = 12 . lb/in and 𝜌 ∞ = 1 . −7 lb ⋅ s ∕ in , respectively. – Free-stream Mach number 𝑀 ∞ = 0 . , angle of attack AoA = 0 ◦ , and angle of sideslip AoS = 0 ◦ . • 𝑊 ( 𝝁 ) denotes the parametric weight of the ARW-2, and 𝑊 ub is its upper bound set to 𝑊 ub = 400 lbs. • The upper bound for the parametric von Mises stress field is set to 𝜎 ub = 2 . psi. • ∈ ℝ and ∈ ℝ are vectors of ones and zeros, respectively, and the lower bound coefficient 𝜁 lb is set to 𝜁 lb = 3 . −4 . TABLE 5
Geometrical and material descriptions of the ARW-2.Item ValueGeometryWing span 104.9 inroot chord 40.2 intip chord 12.5 inMaterialsSkin, except flaps(various composite materials)Stiffeners (aluminum) E .
03 × 10 psi 𝜌 . −4 𝑙𝑏𝑓 ⋅ 𝑠 ∕ in 𝜈 , grid points (see Figure 10). As in the previous example, the parametric lift and drag are obtainedby postprocessing the computed flow solution. The detailed FE structural model of the wing includes representations of thespars, ribs, hinges, and control surfaces (see Figure 11). It contains a total of , dofs. The parametric weight and von Misesstress field are computed using this FE representation. The flutter avoidance constraint, which is also in this case the most CPUintensive constraint, is approximated using a database of linear FSI PROMs.The 15 design optimization parameters stored in 𝝁 are organized in two groups: eight shapedesign parameters { 𝝁 [1] , 𝝁 [2] , 𝝁 [3] , 𝝁 [4] , 𝝁 [5] , 𝝁 [6] , 𝝁 [7] , 𝝁 [8]} , and seven structural design parameters { 𝝁 [9] , 𝝁 [10] , 𝝁 [11] , 𝝁 [12] , 𝝁 [13] , 𝝁 [14] , 𝝁 [15]} . The eight shape design parameters are chosen as follows: 𝜇 acts on thestretching of the wingspan; 𝜇 and 𝜇 modify the twist angle of the wing; 𝜇 and 𝜇 modify its dihedral angle; 𝜇 acts on itsbacksweep angle; 𝜇 modifies the wingspan taper; and 𝜇 controls the chord stretching of the wing.Again, all shape changes dictated by the optimization procedures are effected using Blender .Figures 12–14 graphically depict the effects on the shape of the ARW-2 of the specified lower and upper bounds for thewingspan stretching, twist angle, dihedral angle, backsweep angle, wingspan taper, and chord stretching of this wing.The seven structural design parameters are chosen as seven thickness increments for seven different groups of stiffeners (seeFigure 15). As before, each increment is defined as a percentage of the initial thickness to which it is applied, and is constrainedto have a magnitude less than 10% (box constraint). oncoraglio ET AL FIGURE 10
ARW-2: cylindrical computational fluid domain with a symmetry plane (left); and discretization by an inviscid,unstructured, body-fitted mesh (right).
FIGURE 11
ARW-2 – FE respresentations of: spars and ribs (top); and skin (bottom). Boncoraglio
ET AL
FIGURE 12
ARW-2 – wing shapes associated with: upper (green) and lower (red) bounds on the wingspan stretching (top);and upper (green) and lower (red) bounds on the twist angle (bottom). oncoraglio
ET AL FIGURE 13
ARW-2 – wing shapes associated with: upper (green) and lower (red) bounds on the dihedral angle (top); and upper(green) and lower (red) bounds on the sweep angle (bottom). Boncoraglio
ET AL
FIGURE 14
ARW-2 – wing shapes associated with: upper (green) and lower (red) bounds on the wingspan taper (top); andupper (green) and lower (red) bounds on the chord stretching (bottom).
FIGURE 15
ARW-2 – organization of the stiffeners of the wing in seven different groups represented by seven different colors. oncoraglio
ET AL FIGURE 16
Singular values of the snapshot matrix 𝕊 computed using the alternative AS method. Given that for the MDAO problem (32) the PROM-based computational framework equipped with the classical AS methodfor indirectly sampling the design parameter space failed to improve the objective function beyond its initial value, only twodatabases of FSI PROMs are considered here for the solution of the MDAO problem (34): • 𝐵 , which is constructed as for the MDAO problem (32) (see Section 4.2.1). However, the greedy procedure samplesin this case 𝑁 𝐵 = 415 feasible design parameter points. • 𝐵 , which is constructed as for the MDAO problem (32) (see Section 4.2.1). In this case however, the snapshots of theform Δ 𝝁 are adaptively computed by solving the auxiliary optimization problemmaximize 𝝁 ∈ ⊂ ℝ 𝐿 ( 𝝁 ) 𝐷 ( 𝝁 ) subject to 𝑊 ( 𝝁 ) ≤ 𝑊 ub 𝜎 VM ( 𝝁 ) ≤ 𝜎 ub 𝝁 lb ≤ 𝝁 ≤ 𝝁 ub which is obtained by removing from the MDAO problem (34) the CPU intensive FSI constraint. This auxiliary problemis solved in 19 iterations using the SLSQP method. The 19 iterations generate 20 snapshots. The computation of each ofthe snapshots consumes about 5.3 minutes wall-clock time on a Linux cluster with 36 processors. Next, the snapshots arecompressed into a ROB 𝐕 𝜇 of dimension 𝑛 al = 7 (see below).Figure 16 plots the distributions of the singular values of the snapshot matrix 𝕊 (15). It shows that the last eight singularvalues of this matrix are much smaller than the first seven ones. This suggests constructing a ROB 𝐕 𝜇 as the set of the first sevencolumns of the matrix 𝐔 arising from the SVD of 𝕊 – that is, constructing a ROB 𝐕 𝜇 of dimension 𝑛 al = 7 .Figure 17 graphically depicts the characteristics of the sampling performed for constructing the database 𝐵 . First, a set of2,197 preliminary points is generated in al using a uniform tensor product in this parameter space. Next, Algorithm 3 is appliedto select among these points a set Ξ 𝑟 of 429 feasible candidate points (see Figure 17-left), and determine the associated set offeasible points Ξ ∈ (see Figure 17–right). Then, the greedy algorithm described in Section 3.4 is directly applied to Ξ 𝑟 – andindirectly to Ξ – to perform the final parameter sampling and adaptively construct the database of FSI PROMs 𝐵 . Each plotin the 𝑖 -th diagonal of Figure 17-left, 𝑖 -th diagonal of Figure 17-right, ( 𝑖, 𝑗 ) off-diagonal of Figure 17-left, and ( 𝑖, 𝑗 ) off-diagonalof Figure 17-right has the same meaning as in Figure 7 of Section 4.2.1. Boncoraglio
ET AL
FIGURE 17
Set Ξ 𝑟 ∈ al containing 429 feasible candidate points determined by the feasible AS algorithm as well as thealternative approach for constructing an AS (left), and associated set Ξ ∈ containing the same number of feasible candidatepoints (right).In each constructed database, each pre-computed linear FSI PROM is built as described in Section 3.2, with 𝑛 𝑓 = 100 and 𝑛 𝑠 = 6 . Hence, each pre-computed linear FSI PROM has the dimension 𝑛 𝑞 = 𝑛 𝑓 + 2 𝑛 𝑠 = 112 .Table 6 reports on the performance of the offline phase for problem (34) of the computational framework for MDAO describedin this paper, depending on whether an AS is used or not. The same observations made in the previous example can also bemade in this case. Additionally, the following comments are noteworthy: • Whereas the advocates of the classical AS method recommend a number of samples 𝑁 𝐒 that grows logarithmically withthe dimension of the design parameter space 𝑁 (see (12)), the number of samples 𝑁 𝕊 obtained for both applicationsexplored in this paper using the greedy procedure does not obey the empirical formula (12). • As 𝑁 is increased, the speed-up relative to the case where no AS is used increases. oncoraglio ET AL TABLE 6
ARW-2: performance results for the offline phaseMethod 𝑁 𝑁 𝕊 Wall-clock time 𝑁 𝐵 Wall-clock time Wall-clock time Speed-upcomputing 𝐕 𝜇 greedy procedure total offline relative to no ASw/o AS 15 N.A. N.A. 415 93.2 hrs 93.2 hrs N.A.Alternative AS 7 20 1.77 hrs 203 38.5 hrs 40.27 hrs 2.31 Figure 18 and Figure 19 display the convergence histories of the objective function and constraints of the MDAO problem (34),respectively. They highlight that, as in the case of the CRM MDAO problem discussed in Section 4.2, the alternative AS methodleads in this case to almost the same optimal objective function and constraints as the direct sampling of , and in roughlythe same number of iterations. Again, the alternative AS method leads during the first few iterations to values of the optimizedobjective function and constraints that are similar to those obtained by directly sampling the full space . This, despite the factthat the flutter constraint is active early on during the iterations, and therefore the ROB 𝐕 𝜇 based on the snapshot matrix 𝕊 (15)is by design sub-optimal; nevertheless, this ROB leads to good results. FIGURE 18
ARW-2: convergence histories of the objective function.Table 7 compares the local optimal solutions of the MDAO problem (34) obtained using the two different PROM databases 𝐵 and 𝐵 . The reader can observe that the optimal solution computed using the PROM database 𝐵 – that is, the aternativeAS method for sampling – is sufficiently close to its counterpart computed using the database 𝐵 . Since both databases oflinear FSI PROMs lead to the same optimal value of the objective function (see Figure 18), the two solutions they deliver aredifferent because they correspond to different local optimal points.In summary, for this MDAO problem with fifteen design optimization parmeters, the alternative subspace method speeds upthe offline computations by a factor bigger than two, while leading to a reasonably close local optimum solution. Boncoraglio
ET AL
FIGURE 19
ARW-2: convergence histories of the constraints.
TABLE 7
ARW-2: locally optimal computed solutions 𝝁 [1] 𝝁 [2] 𝝁 [3] 𝝁 [4] 𝝁 [5] 𝝁 [6] 𝝁 [7] 𝝁 [8] 𝝁 [9] 𝝁 [10] 𝝁 [11] 𝝁 [12] 𝝁 [13] 𝝁 [14] 𝝁 [15] Initial 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0Lower bound -2 -1 -0.1 -4 -4 -4 -4 -2 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1Upper bound 2 1 0.1 4 4 4 4 2 0.1 0.1 0.1 0.1 0.1 0.1 0.1Optimal using 𝐵 𝐵 A new take on the concept of an active subspace is presented for mitigating the curse of dimensionality associated with sam-pling a design parameter space for the purpose of model reduction in multidisciplinary analysis and optimization (MDAO).Instead of performing an empirical, a priori sampling of the gradient of the objective function of interest in order to construct areduced-order basis for the optimization parameters, the alternative approach presented in this paper samples instead the solutiontrajectory of an economical version of the MDAO problem of interest in incremental form using an adaptive, greedy procedureguided by an efficient, residual-based error indicator. The new approach is intertwined with the concepts of projection-basedmodel order reduction and a database of linear, projection-based reduced-order models (PROMs) equipped with interpolationon matrix manifolds, in order to construct a complete and efficient computational framework for MDAO problems incorporat-ing a computationally intensive linearized fluid-structure interaction constraint. The framework is illustrated using two different oncoraglio
ET AL MDAO problems associated with two different aeroelastic systems governed by flutter constraints: a flexible configuration ofNASA’s Common Research Model; and NASA’s Aeroelastic Research Wing
ACKNOWLEDGMENTS
The authors acknowledge partial support by the Office of Naval Research under Grant N00014-17-1-2749, partial support bythe Boeing Company under Contract Sponsor Ref. 134824, and partial support by a research grant from the King AbdulazizCity for Science and Technology (KACST). This document does not necessarily reflect the position of these institutions, and noofficial endorsement should be inferred.
References
1. Zahr MJ, Farhat C. Progressive construction of a parametric reduced-order model for PDE-constrained optimization.
International Journal for Numerical Methods in Engineering
AIAA Journal
SIAM Journal onScientific Computing
ESAIM: Mathematical Modelling and Numerical Analysis
International Journal for Numerical Methods in Fluids
Computers& structures
Transactions of FAMENA
Computer Methods in Applied Mechanics and Engineering
Computer Methods in Applied Mechanics and Engineering
Quarterly of Applied Mathematics
Technometrics arXiv preprint arXiv:1408.0545 Boncoraglio
ET AL
14. Amsallem D, Farhat C. Stabilization of projection-based reduced-order models.
International Journal for NumericalMethods in Engineering
ComputerMethods in Applied Mechanics and Engineering
Journal of Computational Physics
Psychometrika
Journal of Computational Physics, (submitted) .19. Farhat C, Geuzaine P, Brown G. Application of a three-field nonlinear fluid–structure formulation to the prediction of theaeroelastic parameters of an F-16 fighter.
Computers & Fluids
AIAA Journal