A New Framework for Model Reduction of Complex Nonlinear Dynamical Systems

Considerable progress in computing technology in the past decades did not alleviate difficulty inherent in simulating complex dynamical systems. Reduced order models (ROMs) can be used to unburden these systems of redundant computations. While a variety of methods have been developed for reduced order modeling, they cannot be used for parametric study of nonlinear and complex systems, wherein we constantly change the parameters, input values, and energy levels. Parametric study is essential to determine the dynamics of complex systems. Only robust and persistent reduced order models, which remain stable with these changes, can be used for parametric study. In this dissertation, we develop a framework which measures the robustness and persistency of reduced order models. The framework quantifies the changes in the reduced models and singles out the most robust and persistent ones. The main advantage of this methodology is that it is applicable to the majority of databased model reduction methods. The approach begins with specifying a range of system’s initial states, parameters, and inputs for the parametric study. The data is collected from simulations of the system with the parameters chosen randomly within that range. The dominant structures of data are then identified using the multivariate analysis methods such as proper orthogonal decomposition (POD) and smooth orthogonal decomposition (SOD). The framework identifies the persistent and robust structures and combines them to obtain the models suitable for parametric study within the specified range. Our aim is to investigate the fidelity of the framework for persistent model order reduction of large and complex dynamical systems. The framework is validated using several numerical examples including a large linear system and two complex nonlinear systems with material and geometrical nonlinearities. While the method is used for identifying the robust subspaces obtained from both POD and SOD, the results show that SOD outperforms POD in terms of stability, accuracy, speed, and robustness. Also, showing that SOD-based ROMs are robust, we no longer need to simulate full-scale models for many parameters. We only need to do few simulations using the full-scale model to build ROMs. In addition, we extend the application of the proposed approach to model order reduction of nonlinear control systems. We use SOD to identify the dynamically relevant modal structures of the control system. The identified SOD subspaces are used to develop persistent ROMs. Performance of the resultant SOD-based ROM is compared with POD-based ROM by evaluating their robustness to the changes in system’s energy level. Results show that SOD-based ROMs are valid for a relatively wider range of the nonlinear control system’s energy when compared with POD-based models. In additions, the SOD-based ROMs show considerably faster computation time compared to the POD-based ROMs of the same order. For the considered dynamic system, SOD provides more effective reduction in dimension and complexity compared to POD.

study is essential to determine the dynamics of complex systems. Only robust and persistent reduced order models, which remain stable with these changes, can be used for parametric study.
In this dissertation, we develop a framework which measures the robustness and persistency of reduced order models. The framework quantifies the changes in the reduced models and singles out the most robust and persistent ones. The main advantage of this methodology is that it is applicable to the majority of databased model reduction methods. The approach begins with specifying a range of system's initial states, parameters, and inputs for the parametric study. The data is collected from simulations of the system with the parameters chosen randomly within that range. The dominant structures of data are then identified using the multivariate analysis methods such as proper orthogonal decomposition (POD) and smooth orthogonal decomposition (SOD). The framework identifies the persistent and robust structures and combines them to obtain the models suitable for parametric study within the specified range.
Our aim is to investigate the fidelity of the framework for persistent model order reduction of large and complex dynamical systems. The framework is validated using several numerical examples including a large linear system and two complex nonlinear systems with material and geometrical nonlinearities. While the method is used for identifying the robust subspaces obtained from both POD and SOD, the results show that SOD outperforms POD in terms of stability, accuracy, speed, and robustness. Also, showing that SOD-based ROMs are robust, we no longer need to simulate full-scale models for many parameters. We only need to do few simulations using the full-scale model to build ROMs.
In addition, we extend the application of the proposed approach to model order reduction of nonlinear control systems. We use SOD to identify the dynamically relevant modal structures of the control system. The identified SOD subspaces are used to develop persistent ROMs. Performance of the resultant SOD-based ROM is compared with POD-based ROM by evaluating their robustness to the changes in system's energy level. Results show that SOD-based ROMs are valid for a relatively wider range of the nonlinear control system's energy when compared with POD-based models. In additions, the SOD-based ROMs show considerably faster computation time compared to the POD-based ROMs of the same order.
For the considered dynamic system, SOD provides more effective reduction in dimension and complexity compared to POD.

ACKNOWLEDGMENTS
For all researchers who paved the way before me upon whose shoulders I stand.
I would like to thank the many people who helped and supported me throughout my journey in the graduate school, and those who were alongside me in discovering beautiful Rhode Island.
I would like to express my deepest gratitude to my supervisor, Professor David Chelidze, for his unwavering support and mentorship. I very much enjoyed learning from him new insights into the world of complex dynamics. For this, I will be always grateful to him.
My sincere thanks must also go to the Dr. Richard Vaccaro, Dr. David Taggart, Dr. Hongyan Yuan and Dr. Jason Dahl for serving on my committee. In particular, I took two classes in the field of control with Dr. Vaccaro, which was an important motivation for extending this research to control systems. Considerable progress in computing technology in the past few decades did not alleviate difficulty inherent in simulating complex dynamical systems. Examples of such systems are large-scale finite difference/element, multi-body dynamics, or geometrically nonlinear models, and molecular dynamics simulations [1,2,3,4,5,6]. A reduced order model (ROM) for these systems can be used to significantly reduce redundant computations and data storage requirements [7]. In particular, persistent ROMs, which are robust to the changes in initial conditions, system parameters, and loading conditions, can be used in parametric studies that are prohibitive when using a full-scale model.
While a variety of methods for model order reduction (MOR) have been developed, very few of them provide persistent ROMs. Often emphasis is only on the accuracy of the ROMs and their ability to capture the dynamics of the fullscale models for a fixed set of parameters, and operating and loading conditions. However, the importance of the robustness of a ROM to the changes in those parameters is often not accentuated. We consider a ROM to be persistent if it is robust to changes in a full-scale model's energy, forcing, and parameters. Databased reduced order modeling with no persistency is of limited scope; ROMs built on the data generated from the simulations of a full-scale model can only be used for simulating the same exact configuration of the model. This ROM might still be of great utility if we can study the long-time dynamics of a system (e.g., protein folding), but cannot be utilized in parametric studies, wherein we repeatedly change the parameters, input values, and energy levels.
In this dissertation, we present a new framework for obtaining persistent ROMs, which are valid within a defined range of the system's energy, which is imposed by changing the input parameters. Our framework can be applied to all data-based MOR methods. We make use of data from simulations or experiments to develop the ROMs. Our goal is to ensure that the obtained persistent ROMs are robust and can be used for simulating the system with any chosen parameter from the defined range.

Background and Prior Work
Persistent MOR for linear systems has not attracted extensive research focus since the linear modal structure is not dependent on the energy of a system. As a result, if a ROM is properly developed for one energy level of a linear system, it should also be valid for the other energy levels. The methodologies for MOR for these systems are mostly projection based, where the linear subspaces used in the projection can be related to the modal space which is span by linear normal modes (LNMs). For example, the modes identified using proper orthogonal decomposition (POD) (also known as singular value decomposition, principal component analysis, or Karhunen-Loève expansion) [8,9,10,11,12] approximate the LNMs for systems with uniform mass distribution [13]. In [14], a method for combining POD and a Hessian-based model reduction approach is proposed. Other popular methodologies for the MOR of linear systems include the Galerkin reduction using linear normal modes (LNMs) [15,16], Krylov subspace projections [17], Hankel norm approximations [18,19], and truncated balance realizations [20,21].
A recent survey published by Benner et al [22] is a good read for reviewing the projection-based model reduction methods for parametric studies. According to this work, which is mostly focused on linear systems, the model order reduction theory in the case of nonlinear systems is much less developed. In another work [23], parametric ROMs are developed for thermal modeling of electric motors.
The simulation time for these models are reduced by a factor of 500, however, the considered models are linear. In [24], the authors suggest an interpolation methodology which is applied to the POD basis. Using this methodology, they estimate aeroelastic damping ratio coefficients to within 10% accuracy. A parametric model reduction approach proposed in [14,25], is applied to a convection-diffusion model which is parametrized by the initial conditions. In these cases, the authors state the limitation of the method regarding the lack of an explicit connection to the transient observability gramian for nonlinear cases.
Nonlinearity, the integral part of complex dynamical systems, makes the development of persistent ROMs a much harder problem. Many approaches for nonlinear MOR are based on the extending methodologies used for linear MOR.
For example, linearization about an equilibrium point was used for the reduction of weakly nonlinear systems [26,27]. Many other approaches are derived from POD [8,28,11,29,30,9], and some from balanced truncation [31,32,33]. Other approaches include neural networks [34], Volterra theory [35], and inertial manifold approximation [36]. More recently, a method called Proper Generalized Decomposition (PGD) has been developed as a generalization of POD in order to construct a priori ROM [37,38,39,40]. This method has a potential for solving multidimensional problems since it does not require any knowledge of the solution [39,41].
The interested reader can find a review on PGD-based MOR techniques in [42].
In summary, a majority of the methodologies commonly used for MOR of nonlinear systems can be categorized into two groups. In the first group, nonlinear normal modes (NNMs) or their approximations [43,44,45,46,47,48,49] are used. In the second group, combined with the Galerkin projection, linear subspaces obtained from spatiotemporal decompositions such as POD and smooth orthogonal decomposition (SOD) are utilized [2,50,10,30,9,51]. Linear subspaces are of considerable current interest because they are computationally tractable and do not neglect the nonlinearity of the original vector-field [8], while, in general, the calculation of NNMs is difficult [52,53,54,55]. Also, MOR based on NNMs suffers from another major drawback related to changes in the NNMs with the variation in system's level of energy [56,53]. The dependence of the NNMs on the energy level causes an insufficient robustness of the corresponding NNM-based ROMs to the changes in the system's energy level. Thus, NNM-based ROMs cannot be considered truly persistent.

Our Approach to Persistent Reduced Order Modeling
Our approach is based on identifying robust subspaces which embed NNMs and do not change drastically as the system changes its energy level. Note that linear subspaces used for MOR are to be identified in such a way that the active NNMs are embedded in them [13]. These subspaces may still change as the NNMs change with the system's level of energy [56]. However, depending on the decomposition method, some particular subspaces may be robust to variations such as changes in initial conditions, external excitations, energy levels, or systems parameters. Our hypothesis is that while an individual NNM may change with energy, a linear subspace embedding this mode may not undergo any considerable change.
Identifying such linear subspaces would enable us to obtain the persistent ROMs that are robust to a relatively wide range of system parameters and operating conditions.
The new framework for persistent MOR of large, complex systems based on the concepts of subspace robustness and dynamical consistency is investigated [57,58,2]. Subspace robustness characterizes how a linear subspace changes under different conditions of the system, which can be used for complex systems to identify the subspace characteristics that lead to a persistent MOR. Dynamical consistency evaluates the deterministic properties of the full-scale system's trajectory projection onto the corresponding linear subspace. It indicates the ability of the identified subspace to potentially-but not necessarily-result in a stable and accurate ROM.
The utility of our framework will be initially evaluated by applying it to the POD subspaces since they are widely used for MOR. POD's drawback for deterministic systems is that it only considers the statistical (i.e., spatial) characteristics of the data [59]. It only prioritizes the maximal variances in the multivariate data and may disregard important dynamical features that have small variances.
Changing the energy level of a system may drastically alter dynamic features that previously had small variances, which will not be reflected in the identified POD modal structure. Therefore, POD, while providing an optimal reduction-in the least squares sense-for a system with fixed set of parameters and forcing, might not be a suitable choice for the persistent MOR of complex systems. For example, a nonlinear Euler-Bernoulli beam exhibits small-amplitude longitudinal oscillations as the input energy level increases. These oscillations may not affect the dynamics structure identified by POD. The subspace obtained from SOD, which was first used in 2005 for vibration mode identification [59], will also be considered within our framework. SOD can be viewed as an extension of POD, which acquires the ability to separate multivariate data based on inherent characteristic frequencies.
In other words, it not only considers the spatial statistics, but also looks at the temporal characteristics of data. Thus, SOD subspaces are likely to be less sensitive to the changes in the energy and properties of the system, and may provide for the persistent MOR.
The focus of this study is on complex, nonlinear dynamical systems. However, a lightly damped linear system will be considered first. The rationale behind this consideration is twofold: (1) the assertion that POD recovers LNMs for systems with uniform mass distribution [13] has been only tested on fairly low-dimensional systems, with fairly long time series; and (2)  In addition, we can also evaluate the ability of these methods to provide robust subspace identification for a system that actually possesses this robustness in all LNMs.
Following the example of the linear systems, MOR of three large-scale, complex nonlinear systems will be studied as the main subject of this work. POD and SOD will be used for multivariate analyses of the associated ill-conditioned data matrices from these systems. The POD-and SOD-spanned subspaces will be tested using the framework to identify the robust subspaces for persistent ROM development. The resultant ROMs subjected to different energy levels will be simulated using several numerical examples. The validity of the results will be investigated in terms of the stability and accuracy of the ROMs.
This dissertation is organized as follows. In the current chapter, the procedure for projection-based nonlinear model reduction is reviewed. Multivariate analysis methods using POD and SOD are reviewed and demonstrated using geometrical interpretations. Chapter 2 describes the developed framework for persistent MOR. In Chapter 2, we obtain the metrics for subspace robustness and dynamical consistency, which are used as the basics of persistent MOR framework. In Chapter 3, we apply the framework to a large-scale linear system. Chapter 4 focuses on obtaining persistent ROMs for nonlinear systems. Chapter 5 intoduces separated multivariate analysis for model order reduction. In Chapter 6, we extend the idea of persistent MOR using SOD to control systems. Chapter 7 concludes this dissertation, highlights the main finding, and suggests future work.

Nonlinear Model Reduction
The full state-space model of a deterministic, nonlinear dynamical system has the following general form:ẏ where y ∈ R 2n is a dynamic state variable, n ∈ N is the number of the system's degrees of freedom, f : R 2n × R → R is some nonlinear flow, and t is time.
As mentioned earlier, there are two approaches used for model order reduction of a dynamical system. The first and the most common one, which is the underlying idea of this dissertation, is based on projecting the dynamical system onto a lower dimensional linear subspace of the state space to yield a reduced order model. In the other approach, the state variable of the system is mapped onto some lowerdimensional nonlinear manifold using a nonlinear coordinate transformation [60,61,62,63,64,65,66,67,68,69,70,71,72,73]. In the next section, both approaches are reviewed.

Model Reduction using Galerkin Projection
There are several methods to identify a linear subspace for dynamical systems.
In case of linear systems, LNMs [74,69] are suitable as the basis. A basis can be also identified using multivariate analysis. The state variable trajectory data can be arranged in the matrix Y = [y 1 , y 2 , . . . , y 2n ]. Multivariate analysis methods, as it will be outlined in section 1.3, will be applied to this data matrix to identify a basis for model reduction. As an output for these analysis, the dominance of each basis vector or mode is given by its corresponding principal value. The most dominant kdimensional basis vectors are arranged in the matrix P k = [e 1 , e 2 , . . . , e k ] ∈ R 2n×k .
We neglected the non-dominant (or submissive) basis vectors, thus, k ≤ n. The reduced state variable is obtained using the following coordinate transformation: The coordinate transformation is plugged into Eq. (1) to yield: Multiplying both sides by P † k one obtains: where (·) † indicates the pseudoinverse of (·). Eq. (4) is the reduced order model for the full-scale model described by Eq. (1), and can be stated in terms of new vector-valued function g:q = g(q, t).
Eq. (5) provides the k-dimensional ROM described in the state-space form. For the ROM, k differential equations are needed to be simulated versus n equations in the full-scale model. Since k ≤ n, the reduced-order model is expected to be faster. The simulation results of the reduced-order model can be collected in a snapshot matrix Q = [q 1 , q 2 , . . . , q k ]. Eq. (2) transforms a k-dimensional point q to an n-dimensional point y. The snapshot matrixŶ can be obtained using a similar transformation:Ŷ The computational cost of this transformation to obtain the full-scale model's snapshot data is negligible compared to the total simulation time.
• Given a POD is looking for the best possible approximation in the least squares sense. scalar field Nodal measurements are the columns of matrix Y Figure 1: A scalar field

Multivariate Analysis
Multivariate analysis is based on the statistical principle of multivariate statistics, which involves the process of simultaneously analyzing multiple independent variables using matrix algebra [75,76]. It is being used as a method to identify the modal structure of dynamical systems. The extracted modes from the multivariate analysis can be used for MOR. This section begins with the description of POD, and SOD as an extension to POD. This includes the mathematical formulations as well as geometrical interpretations for both methods, which are provided for finite dimensional cases.
POD and SOD are applied to the recorded measurements of a scalar field. Scalar field, by definition, associates a scalar value to every point in a space, and is coordinate-independent, meaning that any two observer will agree on the value of the measurement. For example, for a beam shown in Fig. 1, the scalar field consists of position and velocity of the nodal points. The measurements can be taken from the scalar field using sensors, or obtained by computations. The recorded measurements form the data matrices for the multivariate analysis.

Proper and Smooth Orthogonal Decomposition for Finite-Dimensional Cases
The state variable measurements of the full-scale system are recorded to form position and velocity data matrices X ∈ R r×n and V ∈ R r×n , respectively. X is composed of r snapshots of n position state variables. Similarly, V is composed of r snapshots of n velocity state variables. Thus, data matrix Y, which we will refer to as full data matrix throughout this dissertation, is formed by combining X and V together, i.e., Y = [X V].
The time derivative of X is V. To obtain a time derivative of V or an acceleration data matrix K, we can use a full model of our dynamical system, Eq. (1).
Alternatively, it can be approximated by K ≈ DV, where D is the matrix form of some differential operator such as forward difference given as Therefore, an ensemble of time derivative of Y will beẎ = [V K]. Provided that Y andẎ are zero mean, the corresponding auto-covariance matrices can be formed by

POD
In POD, we are looking for a basis vector φ ∈ R 2n such that a projection of the data matrix onto this vector has maximal variance. The description of POD translates into the following constrained maximization problem: Plugging Eq. (8) into the above problem, and defining λ(φ) =λ (φ) r−1 , one obtains: Using the problem side constraint, we can rewrite Eq. (9) as follows: We take the first derivative of λ(φ) with respect to φ and equate it to zero, in order to maximize λ(φ): As a result, Using Eq. (9), Eq. (12) is simplified to obtain the solution to the POD problem in terms of the solution to the eigenvalue problem of the auto-covariance matrix Σ yy : where λ k are proper orthogonal values (POVs), φ k ∈ R 2n are proper orthogonal modes (POMs), and proper orthogonal coordinates (POCs) are columns of Q = YΦ, in which Φ = [φ 1 , φ 2 , . . . , φ 2n ] ∈ R 2n×2n . POVs are ordered such that λ 1 ≥ λ 2 ≥ . . . ≥ λ 2n , and reflect the variances in Y data along the corresponding POMs.

Geometric Interpretation of POD
In order to illustrate POD using a simple example, let us consider a scalar field Y that consists of the mean-shifted measurements y 1 (t) and y 2 (t) of a two-degreeof-freedom system. Plotting these data points results in a zero-mean cloud of data shown in Fig. 2. We aim to obtain two POMs as the solution of an optimization (maximization) problem for the two-dimensional case. The norm of the projection of the data onto POMs must be the maxima. The i-th data point y i = (y 1 (t), y 2 (t)) is specified by a red dot. An arbitrary vector φ and its direction are also specified. the (global) maximum will be obtained as vector φ approaches φ 1 , the first POM.
The maximization problem has also another (local) maximum at φ = φ 2 , the second POM. However, the second POM is the trivial solution since it is imposed by the eigenvalue problem in Eq. (13) that φ 1 and φ 2 are orthogonal to each other.
Also, due to the side constraint of the problem φ 1 and φ 2 are normal vectors, and as a result, orthonormal.
It is apparent that the data points have the maximum variance along the first Associated with each POM is a POV, denoted by λ k , which is related to the norm of the data projection onto φ k . Thus, it reflects the variances in Y data along the vector φ k . The greatest POV comes with the first POM along which the data variances is maximum. The second greatest POV comes with the second POM along which the variances are (locally) maximum, and so on. Therefore, each POV represents the amplitude dominance of its corresponding mode.
It is common to refer to POVs as the energy of the modes in engineering context. In fluid mechanics with velocity measurements, POV can be related to the kinetic energy. Chatterjee [77] discusses that in structural dynamics problems with both position and velocity measurements, thinking of POVs as modal energies is not correct. However, throughout this dissertation they are referred as energies since they are an implicit combination of potential and kinetic energy of the system.

SOD
In SOD, we are looking for a basis vector ψ ∈ R 2n such that a projection of the data matrix onto this vector has both maximal variance and minimal roughness (i.e., maximal smoothness). Roughness can be defined as the L 2 norm of the rate of change of data. Thus, the roughness of a one-dimensional scalar filedẎψ is equal to Ẏ ψ . This description of SOD is translated to the following mathematical form: which can be stated as maximizing the following function: We rewrite the above equation in the following form: Substituting Eq. (8) in Eq. (15), one obtain: In order to maximize λ(ψ), we set the first derivative equal to zero: As a result, Eq. (17) can be simplified using Eq. (16) Eq. (19) is the generalized eigenvalue problem of the matrix pairs Σ yy and Σẏẏ which yields the solution to the SOD problem. In this equation, scalars λ k are smooth orthogonal values (SOVs), and vectors ψ k ∈ R 2n are smooth projection modes (SPMs). A matrix that contains all the SPMs has the form Ψ = [ψ 1 , ψ 2 , . . . , ψ 2n ] ∈ R 2n×2n , and a matrix that contains all the SOVs has the form Λ = diag([λ 1 , λ 2 , . . . , λ 2n ]) ∈ R 2n×2n . Using these definitions, Eq. (19) can be summarized into the following matrix form: A matrix of smooth orthogonal modes (SOMs), Φ can be defined to satisfy the following biorthogonality condition: where I is identity matrix, and as a result Φ = Ψ −T . Smooth orthogonal coordinates (SOCs) are defined as the projection of the data onto SPMs given by: While the modes obtained from SOD are not orthogonal, its coordinates are orthogonal. A proof, using Eq. (21), follows that: One of the properties of the generalized eigenvalue decomposition given by Eq. (19), is that SPMs are orthogonal with respect to Σẏẏ. Thus, SOCs are orthogonal and the corollary is proved.
The degree of coordinates' smoothness is described by the magnitude of the corresponding SOV. Thus, the greater in magnitude the SOV, the smoother in time is the corresponding coordinate. It should be noted that if we were to replace Σẏẏ with the identity matrix, the formulation will yield POD.

Geometric Interpretation of SOD
We consider a scalar field Y consisting of the mean-shifted measurements y 1 (t) and y 2 (t) of a two-degree-of-freedom system with the sampling rate of ∆t = 1.
Plotting these data points results in a zero-mean cloud of data shown in Fig. 3.
We aim to obtain two SOMs, φ 1 and φ 2 , as the solution of an optimization (maximization) problem for the two-dimensional case. We indicate two consecutive data points y i = (y 1 (t), y 2 (t)) i and y i+1 = (y 1 (t), y 2 (t)) i+1 in the figure. The first derivative of the i-th data point or the vector of data evolution is approximated by We refer to this vector as the velocity vector and depict it between data (i) and (i + 1) in the figure. The projection of this vector onto φ 1 is also shown.
We assume that φ 1 , the first SOM, is wandering in the 2D space of the data.
By definition of SOD, we aim to maximize the norm of the projections of data onto this vector. However, we also aim to minimize the norm of the projection Associated with each SOM is a SOV, denoted by λ k , which reflects the ratio of variances in Y data to variances in their first derivativesẎ along ψ k . The greatest SOV comes with the first SOM along which the ratio is maximum. Compare this to the first POM along which only the variance of data is maximum. The second greatest SOV comes with the second SOM along which the ratio is (locally) maximum, and so on. Therefore, each SOV represents the dominance of its corresponding mode in terms of overall spatial variation and temporal smoothness of the coordinate. 1 SOCs are orthogonal to each other In a sense data Y comes from the consecutive mapping of a system's state onto another state using a vector-valued function (flow) f . POD only considers the spatial geometric consequences of the mapping and neglects the temporal flow under which the states have undergone. SOD considers both the states and the flows in terms of overall spatial variation and temporal smoothness of the coordinate.

POD and SOD based model order reduction
For POD based model reduction, we form a matrix of most dominant POMs, , using the solution to Eq. (13). Here we have a particular case of the general procedure described in Section 1.2.1, wherein the reduced state variable is obtained using the coordinate transformation, y = Φ k q.
The coordinate transformation is plugged into Eq. (1) to yield: The k-dimensional subspace Φ k is orthogonal. Thus, we multiply both sides by For SOD based model reduction, we obtain the matrix of SPMs, Ψ and the matrix of SOMs, Φ as the solution to Eq. (19). Since Φ and Ψ are biorthogonal, their k-dimensional representatives, Φ k and Ψ k are also biorthogonal. We use the coordinate transform y = Φ k q as SOD reduced state variable in Eq. (1). This yields: Multiplying both sides by Ψ T k , one obtains: [17] P. Feldmann and R. W. Freund, "Efficient linear circuit analysis by padé approximation via the lanczos process," Computer-Aided Design of Integrated Circuits and Systems, IEEE Transactions on, vol. 14, no. 5, pp. 639-649, 1995.
[38] M. Chevreuil and A. Nouy, "Model order reduction based on proper generalized decomposition for the propagation of uncertainties in structural dynamics," International Journal for Numerical Methods in Engineering, vol. 89, no. 2, pp. 241-268, 2012.
[39] F. Chinesta, A. Ammar, and E. Cueto, "Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models," Archives of Computational methods in Engineering, vol. 17, no. 4, pp. 327-350, 2010.
[41] P. Ladevèze and L. Chamoin, "On the verification of model reduction methods based on the proper generalized decomposition," Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 23, pp. 2032-2047.

CHAPTER 2
Persistent Reduced Order Models The core of the Galerkin projection-based model reduction is the selection of a suitable linear basis. The linear basis spans a linear subspace for the state space and can be LNMs, POMs, SOMs, etc. In common approaches, the number of modes (i.e., the dimension of the linear subspace) is gradually increased until the ROM adequately captures the full-scale model's dynamics. In this chapter, we propose a systematic approach for selecting the linear subspaces for model reduction.
This dissertation focuses on ROMs based on linear subspace of the system's full phase space or its nonlinear extensions. In both cases, this subspace should satisfy two basic requirements to provide useful ROMs: (1) it needs to embed or capture the active NNM manifold, and (2) this embedding needs to be robust to the changes in initial conditions, system parameters, and forcing functions. Here, we evaluate a subspace selection criterion based on the two new concepts of dynamical consistency and subspace robustness.
The appropriate subspace for model reduction can be selected based on newly developed criteria [1,2,3]. These criteria quantify two concepts: dynamical consistency-which demonstrates how well the linear subspace embeds the nonlinear manifold, and subspace robustness-which explains the sensitivity of the subspace to changes in system's level of energy.

Subspace Robustness
We would require that the selected model reduction subspace be robust with respect to the variations in the data used for its estimation. For example, all LNM subspaces are robust since they are unique and not data based, but POMs can vary Figure 4: Schematic of a nonlinear dynamical system with respect to initial conditions, system parameters, and forcing function [4]. For example, POMs capturing maximal energy for a deterministic steady state motion will be generally different from POMs for a stochastically excited system.
For models based on multivariate data analysis, the corresponding subspaces are identified using a finite set of simulated (or experimental) trajectory points of Eq. (1). Since the flow f is nonlinear, the simulated trajectories can be sensitive to the initial conditions and parameters. The response will also depend on the type of forcing used during simulation/testing. The model reduction subspace should not be sensitive to these variations that are expected in practice; otherwise the corresponding ROM might be invalid for some initial conditions, or for perturbed system and forcing parameters.

Robustness of Modal Subspaces
In order to quantify the subspace robustness, we provide a model of nonlinear dynamical system adopted from [5]. The system, shown in Fig. 4, consists of 5 weightless links with the length of 2l which are connected to each other by torsional springs and dampers. Springs and dampers are not drawn for the sake of clarity.
The coordinate θ i measures the angular position of the i-th link as shown in the figure.
A nonlinear system can exhibit different behaviors based on its level of energy, which include both approximately linear behavior near the stable equilibrium points and nonlinear behavior far from those equilibrium points. Our system shows similar behavior. This will be discussed in detail in Chapter 6. Closer to the equilibria the system is described by LNMs, while as we get farther the system evolves on the NNM manifold, which may also change shape as system energy changes.
Therefore, as energy increases not only the angle of the linear subspace that we get from multivariate analysis of the data changes, but we may also need a higher dimensional subspace to capture the NNM of the system. Different data sets from the system simulations with different inputs or initial conditions have different energy level. Therefore, their extracted modal matrices and the corresponding lower-dimensional subspaces may be different.
One of the ways of preparing data sets for multivariate analysis is subjecting the simulated system to random forcing. In order to illustrate the changes in the modal structure, we excite our nonlinear system by the white noise with a chosen cut-off frequency. We expect that as we increase the forcing amplitude, the higher frequencies in the system's response come into account. As a result the modal structure of the system, indicated by the corresponding subspaces, need to be altered to account for higher frequencies.
We need a metric that measures the difference in the modal structure of two different data sets which have different energy levels. One possibility is to measure the minimal angle between their corresponding subspaces using the following definition.
Definition: The minimal angle for two nonzero subspaces P 1 , P 2 ∈ R k is defined to be the number 0 ≤ θ ≤ π 2 that satisfies: For example, we generate data sets with different energy levels by changing the initial condition of the unforced links system. This way, we can control the energy level of the system. The initial angular position and velocity of all links except the first one are set to zero. The initial conditions for the first link is selected from the range −5 ≤ θ 1 (0) ≤ 5 and −2 ≤θ 1 (0) ≤ 2. The data set for each individual selection of θ 1 (0) andθ 1 (0) is simulated and recorded. POD and SOD are applied to each data set to extract the corresponding modal matrices P. Using the minimal angle between two subspaces 1 , we can estimate the changes in the k-dimensional subspaces of the estimated modal matrices for different data sets. We observe that we obtain different modal subspaces for different energy levels of the systems which are imposed by changing initial conditions or external forcing.
One of the goals of MOR in our work is to obtain a global subspace which is suitable for a range of variations in the energy level of a system under investigation.
The conventional method for proper subspace identification for MOR is based on selecting those subspaces which capture most of the system's energy. However, this method would not assure that the subspace is suitable for ROM for an energy-varied system. Therefore, a new metric is required to measure if the obtained subspace is robust or not to the variations in systems' energy. In the following, we discuss a metric to measure the robustness of different subspaces with respect to each other.

Metric for Subspace Robustness
We can change the system's subspaces obtained from multivariate analysis by changing its energy level in two ways: (1) changing the initial conditions of an unforced; and (2) changing the external forcing. For example, we can vary the external forcing by changing its frequency content and/or forcing amplitude.
Regardless of how we change the system's energy, we do s simulations or experiments and assemble the corresponding data matrices. We apply the intended multivariate analysis to the data and obtain s different modal spaces, P 1 , P 2 , . . . , P s corresponding to each simulation. The k-dimensional subspaces P i k and P j k of the modal space are considered linearly dependent if the minimal angle between them, denoted by θ ij , is equal to zero. On the other hand they are said to be linearly independent, if θ ij = π 2 . Each subspace P k consists of k dominant modes. While these k individual modes can be totally different between two data sets, the subspace spanned by them can still be linearly dependent. For example, we need two LNMs to span a plane containing a damped linear oscillator degree-of-freedom in the 2ndimensional vector space of a system. However, these modes are not unique-their linear combination would also span the same plane, which means that as the modes of system change with its energy level, they can still span the same subspace. Here, we propose a subspace robustness metric which determines if the MOR subspace is robust for a range of energy levels. The metric is a quantification of changes in the subspaces for the range of energies. For the subspace robustness close to one we can argue that the subspace is robust to the changes in energy level.
In case of s simulations it is difficult to simply use the angles between all the subspaces to develop a metric for subspace robustness. Here we propose to use singular values of all combined subspaces. Let us assume that k columns of matrix P i k span the k-dimensional subspace P i k . We look at the vectors spanning the subspaces as data which live in the 2n-dimensional space and apply the singular value decomposition to find the principal directions within the data. We form the subspace robustness data matrix S by arranging the subspaces in the following order: . (27) From singular value decomposition of matrix S, we obtain 2n direction vectors φ i in the 2n-dimensional space of data. The standard deviation of subspace data along vector φ i is given by σ i = Sφ i . We define r k = k i=1 σ i φ i to be the extension vector of the subspace data in the k-dimensional space. Then Ker(r k ) = 2n i=k+1 σ i φ i is the extensiuon vector in the null space of the k-dimensional subspace. Thus, the total extension vector in the 2n-dimensional space is r 2n = r k +Ker(r k ). The magnitude of the kernel extension vector, Ker(r k ) , measures the leak of the data into the null space of the k-dimensional space. We compare this magnitude to that of the k-dimensional extension vector, r k . Therefore, the leak into higher dimensional space is evaluated by the angle of extension vectors in the k-dimensional space and its kernel as follows: We define a lower bound for α k by taking the assumption that all the vectors spanning the subspaces are equally distributed in the space. In this case all singular values of matrix S are equal, i.e., σ i = σ. Thus, a lower bound for the k-dimensional subspace,ᾱ k , is:ᾱ Usingᾱ k we map the angleᾱ k ≤ α k ≤ π 2 to 0 to 1 to define γ k as follows: which we call the subspace robustness of the k-dimensional subspace. Alternatively, if we don't consider a lower bound, we can use the following equation [1,2,3]: Singular value decomposition is applied to the whole data to obtain three components of the extension vectors shown in the figure. As an example, r 2 = σ 1 φ 1 + σ 2 φ 2 is the two-dimensional covariance vector of data. Ker(r 2 ) = σ 3 φ 3 is the kernel covariance vector. We calculate the angle between the two-dimensional subspace and its kernel using Eq. (28):

Geometric Interpretation of Subspace Robustness
A lower bound for two dimensional subspace of a three-dimensional space is defined via Eq. (29):ᾱ . Now we can determine the robustness of our two-dimensional subspace via Eq. (30).

Dynamical Consistency
Unfolding of an attractor used in delay coordinate embedding [6] is the underlying idea of dynamical consistency. It can be determined by the premise behind a method of false nearest neighbors [7]. A linear subspace used for reduced order modeling is said to be dynamically consistent if the resultant trajectories are deterministic and smooth. We quantify the dynamical consistency of a reduced-order, k-dimensional flow g which is defined in Eq. (5) as the projection of the original flow f onto a k-dimensional subspace. We estimate the number of false nearest neighbor (FNN) points since they are caused by folding of or intersections in a trajectory. The data points consisting the phase space trajectory can be neighbor because they are temporally close to each other ( i.e, q i and q i+1 as two consecutive data points), or due to the geometric structure of the flow. Another reason for being neighbor is the folding of trajectories in a phase space. In the latter case they are considered as FNNs and will possibly separate if the dimension of g is increased to k + 1.
The metric for dynamical consistency is defined as the ratio of the number of false nearest neighbors (FNN) over the total number of nearest neighbor pairs in a particular k-dimensional subspace: where N k fnn is the estimated number of FNNs in k-dimensional subspace due to projection, and N nn is the total number of nearest neighbor pairs used in the estimation. If ζ k is close to unity, then that k-dimensional subspace is dynamically consistent.
The nearest neighbor search for each test point is accomplished by utilizing a kd-search algorithm [8]. N k f nn is estimated by comparing the distance between the temporally uncorrelated nearest neighbors in a k-dimensional space to the distance between the same points in the (k + 1)-dimensional space. If the change in the distance is one order of magnitude larger than the original k-dimensional distance, then these points are denoted as FNNs in k-dimensional space.

CHAPTER 3
Large-Scale Linear Systems In this chapter, we consider a large-scale linear system with uniform mass distribution for applying the persistent MOR framework. It is shown that in this case, given enough data from long time series, POD is able to recover LNMs [1], and SOD does not require uniform mass distribution for convergence to the LNMs [2]. Yet we need to test the performance of both methods for a large scale systems with limited data in different loading scenarios.

Derivation of the full-scale model
The linear system under investigation is an n-degree-of-freedom mass-springdamper system shown in Fig. 8, where n blocks of masses are connected in series to each other as well as both sides of the support by linear dampers and springs.
The masses can vibrate in x-direction with no friction. The system is described by the following governing differential equations: ] T as the vector of 2n state variables, the full state-space model of the system can be obtained where f e (t) = [0, . . . , 0] 1×n , [1, . . . , 1] 1×n T ; 0 ∈ R n×n is a zero matrix; I ∈ R n×n is an identity matrix; and M, K, and C are n × n mass, stiffness, and damping matrices, respectively.

Application of Subspace Robustness and Dynamical Consistency
The MOR objective is to develop persistent ROMs for harmonically excited systems considered here. The linear system will be excited by a force with a frequency close to its first natural frequency. Modal subspaces for model reduction can be obtained from different types of excitations, including both harmonic and random. With random forcing, we are more likely to explore nearly all the statespace of a dynamical system and excite its all dominant frequencies. However, the particular forcing function has to be carefully selected, especially for the systems that have combined slow and fast dynamics. This is to limit the contamination of the identified modes by forcing that can obscure the true modal structure of the system.
The white noise is used to excite the system because there is no relatively fast dynamics in the presence of a slow dynamics. We perform 12 independent was done for a total time equal to 100 cycles of a harmonic forcing, with the frequency equal to 110% of the first natural frequency of the linear system. With the chosen parameters, the total simulation time was equal to 709.8 s. We recorded 100 data samples in each cycle of applied external forcing. Therefore, a total of 10,000 data points were recorded from each simulation.
In each case, POD and SOD were used to extract the modes out of each data set using the procedure explained in Chapter 1, Section 1.3. The first k dominant modes identified from each simulation independently, spanning 12 k-dimensional subspaces, were concatenated into the matrix S as explained in Section 2.1.2.
Singular value decomposition was applied to matrix S in order to extract the singular modes and the corresponding singular values. Using singular values and Eq. (31), the robustness of the k-dimensional subspaces were evaluated for each model and decomposition scheme. deterministic trajectories, each corresponding to different forcing amplitudes, and then averaged out. In case both subspace robustness and dynamical consistencies of the extracted modes were close to unity, we considered them as suitable for persistent MOR.

Numerical Example
The parameters of the linear system were fixed as follows: The obtained POMs and SOMs for this system are used to get the subspace robustness and dynamical consistency of the ROM subspaces. For the randomly driven linear system, as depicted in Fig. 9, the SOD subspace robustness metric reaches and stays close to unity for k ≥ 3. POD subspace robustness is close to unity at k = 2 and k = 3. It drops at k = 4 and again reaches unity at k = 25 in a non-monotonic manner.
The POD subspace robustness does not behave monotonically and does not stay close to unity once it reaches it. Therefore, in the cases considered, POMs cannot approximate fixed LNMs in a robust manner, except maybe few lower modes. This shows that POMs are not robust under constraints of limited timehistory of high-dimensional data, which makes it unreliable for persistent MOR.
In contrast, the SOD subspace robustness monotonically increases, reaches unity for a low dimension, and does not fluctuate thereafter. Figure 9 also shows the dynamical consistency of POD-and SOD-based sub-spaces for the randomly driven linear system. For both POD and SOD, the dynamical consistencies are similar reaching unity at k = 2 for SOD and k = 5 for POD. This means that the projection of the linear system's deterministic trajectories onto the five-dimensional POD-based or the two-dimensional SOD-based dominant subspaces have no singular point or intersection with itself-or they do not violate the uniqueness of the deterministic evolution.
The linear system subjected to harmonic excitation was simulated using the POD-and SOD-based ROMs via Eq. (58). The phase portraits for the vibrations of the thirtieth mass are depicted in Fig. 10 and Fig. 11. The ROM simulations results shows a very good visual correspondence to the full-scale system using both POD and SOD. Both methods are able to capture the dynamics in two-and threedimensional ROMs and none of them outperformed the other irrespective of the robustness of the corresponding subspaces.
A question arises as to why some relatively non-robust POD subspace-based ROMs like four-dimensional model still correlate with the full-scale model. It should be noted that the subspace robustness metric is of more importance for lower dimensions since they possess most of the energy of the system. The twoand three-dimensional POD subspaces for the linear system are robust and capture most of the system's energy, thus, they provide for good ROMs of the system.
Increasing the dimension of the ROM reduces the robustness of the associated POD subspace to 0.85 but it does not affect its accuracy or stability. This is mainly due to the fact that the fourth POM does not capture enough associated energy to have sizable effect on the corresponding ROM. This also explains why robustness of MOR based on POD has not been much of research concerns for linear systems. As shown in Fig. 12, POD captures most of the energy in the very first few modes, which are robust to the changes in the energy of the system. Therefore, any suitably developed POD-based ROM could probably account for other similar conditions.

CHAPTER 4 Complex Nonlinear Dynamical Systems
Complexity in dynamical systems can arise for different reasons. In nonlinear dynamical systems, it is related to the size, nonlinearity, weak coupling between the DOFs resulting in simultaneous presence of slow and fast dynamics, or a combination of them. In this chapter, we study three different complex dynamical systems: a system with nonlinear spring coupling, a nonlinear Euler-Bernoulli beam, and a mass-spring-grid system.

System with Nonlinear Spring Coupling
The first nonlinear system used here is obtained by adding a nonlinear spring to the linear system as shown in Fig. 13. In this case, the complexity is caused by the large size as well as the material nonlinearity.
The system dynamics is described by the following full state-space equations of motion:ż where everything is the same as the linear system described in Chapter 3 except nonlinear term f nonlinear (z) ∈ R 2n , which has only one nonzero element: [f nonlinear (z)] 2n = −αz 3 n .

Reduced Order Models
As we discussed, with a random forcing we are more likely to explore nearly all the state-space of the system and excite all dominant frequencies. There is no relatively fast dynamics with the presence of slow dynamics. Thus, we use white noise to excite the system.
We applied 12 different external stochastic excitations with different energy Figure 13: Schematic of the system with nonlinear spring coupling level, imposed by changing the forcing amplitude, in order to obtain data for multivariate analysis. Each simulation was done for a total time equal to 100 cycles of a harmonic forcing with the frequency equal to 110% of the first natural frequency of the corresponding linearized system. The total simulation time was equal to 495.1 sec. We recorded 100 data samples in each cycle of applied external forcing.
Therefore, a total of 10,000 data points were recorded from each simulation.
POD and SOD were used to extract the modes out of each data set and matrix S was formed as explained in Section 2.1.2, the similar procedure that we used for the linear system. Subspace Robustness was calculated using Eq. (31) and the dynamical consistency using Eq. (34). Similar to the linear system, the dynamical consistencies were obtained for five deterministic trajectories, each corresponding to different forcing amplitudes, and then averaged out.
ROMs for nonlinear systems are expected to be more sensitive to the robustness of the corresponding MOR linear subspaces. Therefore, in case the lowdimensional subspaces have good robustness, non-robust higher dimensional subspace may destabilize the numerical scheme for the model or at least adversely affect its accuracy.
For investigation of the system with nonlinear spring coupling, the number of DOFs was set to 60 and the other parameters were fixed as follows: m = 1 kg, k = 3600 N/m, c = 720 Nm/s, This results in a rich dynamic response with two stable and one unstable static equilibrium points. The bifurcation diagram, shown in Fig. 14, is plotted for the amplitude of the 30th mass, obtained from the full scale model of the system.
The system is subjected to a harmonic forcing with the frequency Ω = 9.97 Hz, which is 110% of the first natural frequency of the corresponding linearized system.
The forcing amplitude is used as a bifurcation parameter, which increases with an increment size of 0.002. For each forcing amplitude, the first 25 cycles of the system's response are neglected in order to remove the effect of transient.
For the steady-state response, x 30 is recorded at a fixed phase, the end of each cycle of the applied harmonic force, for a total of 25 cycles. Before increasing the forcing amplitude by one increment, this process is repeated three times with different initial conditions to ensure that all branches of the bifurcation diagram are obtained. Our particular aim for the persistent ROM is to reproduce these bifurcation results, which will demonstrate robustness of the ROM to a range of forcing amplitudes or different input energy levels.
The subspace robustness and dynamical consistency for this system are depicted in Fig. 15. The robustness for the SOD subspaces reaches unity at k = 4, while for POD it does not happen until the very end. POD subspace robustness is fluctuating and sometimes getting worse as the subspace dimension increases.
These fluctuations are of greater importance for lower dimensional subspaces since most of the system's response energy is captured in these subspaces. The dynamical consistency for both methods is similar and reaches unity at k = 2. At k = 5, however, the dynamical consistency of the SOD method slightly drops, which may affect the accuracy of the corresponding ROM. While two-and three-dimensional POD subspace robustness are relatively close to unity, they do not account for a significant portion of the system's total energy to provide stable ROMs. The robustness of the four-dimensional POD subspace is low, which causes the diverging results of the corresponding ROM simulation. The five-dimensional POD subspace has better robustness, and the simulations showed that it provides stable ROM, yet is not robust enough to accurately reproduce the bifurcation diagram. The six-dimensional POD-based ROM has better robustness and captures more energy, thus, it results in stable and accurate simulations.
One-through three-dimensional SOD subspaces do not result in stable ROMs because their subspace robustness is relatively low and also they do not capture enough energy of the system. Four-and higher-dimensional SOD subspaces are robust and provide persistent ROMs capable of reproducing the bifurcation diagram The lowest dimensional ROM which provides accurate and robust results is four-dimensional for SOD and six-dimensional for POD. The corresponding bifurcation diagrams are shown in Fig. 16 and Fig. 17. These ROMs are more than fifty times faster in simulation than the full scale model. Thus, we used a four-time finer increment size for the forcing amplitude to provide more details in the bifurcation diagrams of the system. Visual Comparison of these diagrams with the reference diagram shown in Fig. 14 indicates that there is a close match between those of the six-dimensional POD and the full scale model. For the SOD, the bifurcation diagram is little shrunk around q 0 = 0.33. While this bifurcation diagram is not strictly accurate, it still provides a reliable qualitative description the full-scale system's dynamics. In addition, the results for the six-dimensional SOD are as good as the six-dimensional POD, while the test shows that the four-dimensional POD is not even stable, which is consistent with the significant drop of its subspace robustness at k = 4.
In Fig. 18, six-dimensional POD and four-dimensional SOD ROMs are compared to the full-scale model driven by the harmonic forcing with Ω = 9.97 Hz and as compared to k = 4 and 6. This may be explained by the drop in the dynamical consistency of SOD for k = 5, which was shown in Fig. 15.

Geometrically Nonlinear Systems 4.2.1 Nonlinear Euler-Bernoulli Beam
The mathematical model of a nonlinear Euler-Bernoulli beam shown in Fig. 20 is described in classical books devoted to spatial objects dynamics [1,2]; here only the governing differential equations are given: where ρ is density, E is Young's modulus of elasticity, A is the area of beam section, I is moment of inertia, d 1 and d 2 are damping coefficients, and q(x, t) is the transverse forcing. Also, u(x, t) and w(x, t) denote for axial and transverse vibrations, respectively. The beam under investigation is fixed-fixed at its ends with the following boundary conditions that are attached to Eq. (38): In this case, both transverse and longitudinal displacements at beam's ends, as well as tangents to the slope at its ends are equal to zero. Besides, the following initial conditions are considered for the beam: are taken into account, which introduces geometrical nonlinearity to the system.
There is no exact analytical solution for the nonlinear system; however, a finitedimensional state-space model of the system can be solved numerically using finite difference method.
The full-dimensional state-space model of the system in matrix form can be obtained using finite difference method (FDM). The space coordinate of the beam, as shown in Fig. 21, is meshed using M nodes. The parameters u i and w i indicate axial and transverse displacements of the i-th node. The space derivatives in Eq.
(38) are substituted by their second order finite difference approximation O (∆x 2 ) which are given as follows: where ∆x is space mesh size. After simplifying, this leads to the following systems of ordinary differential equations: where α = d 1 ρA , β = d 2 ρA , and G i and H i are nonlinear function of axial and transverse nodal displacements and time.
The boundary and initial conditions can be treated using two fictitious space layers, which are illustrated in Fig. 21. Therefore, using finite difference approximation one obtain: By introducing the state-variables as z i = w i , z M +i = u i , y 2M +i =ẇ i for i = in which g(z, t) = [0 1×2M , H 1 (z, t), . . . , H M (z, t) , G 1 (z, t) , . . . , G M (z, t)] T ∈ R n is the vector of nonlinear functions. The matrix T 2n×2n has the following form: where 0 ∈ R M ×M is a zero matrix and I ∈ R M ×M is an identity matrix. Eq.
(44) represents the full-dimensional mathematical model of the system and is to be solved by numerical methods.

Numerical Example for Nonlinear Euler-Bernoulli Beam
The performances of POD-and SOD-based ROMs will be evaluated in terms of stability and accuracy of the ROMs. For evaluating the accuracy phase portraits vibrations, but small variances for those corresponding to axial vibrations. PODbased MOR only considers the variances while SOD-based MOR not only considers the spatial variances but also looks at the temporal characteristics of the trajectories in the full-dimensional state-space model. Thus, we expect SOD to capture the dynamic of the system in a lower-dimensional subspace. In next section, the results of both POD-and SOD-based MOR are discussed.

ROMs for Nonlinear Euler-Bernoulli Beam
To identify the modal structure using both POD and SOD, the beam is excited by white noise for 3.14 s. This guarantees that nearly all the state space is explored and all dominant frequencies of the beam are covered. The state-variable data from the beam vibration are calculated using the numerical integration. A total number of 2600000 snapshots are recorded with the rate of ∆t of the numerical scheme and stored in the data matrix Y. POD and SOD are applied to the data matrix to extract the modes as described in Section 1.3.
A set of thirty randomly driven trajectory data is used to estimate the subspaces robustness using the aforementioned procedure. Fig. 24 illustrates the sub- The dimension of the lowest dimensional SOD-based ROM that is stable is 47. explicit nature of the FDM, which was used for solving both full-scale and reduced order models. However, another possibility, which could explain this problem, is the weak coupling between axial and transverse displacement of the nodes. This results in small magnitude of the data corresponding to state variables of the axial nodes in the full dimensional state space. Therefore, it might be difficult for both POD and SOD to identify the corresponding subspaces of these data properly.
Nevertheless, SOD was shown to have a better performance in identifying these subspaces.

Mass-Spring-Grid System
Another geometrically nonlinear system that we consider will be called massspring-grid system throughout this dissertation. The schematic of this system is shown in Fig. 28. The system is allowed to vibrate in both x and y directions as a grid of equidistant masses, dampers, and springs. Each spring is assumed to have a corresponding damper acting in parallel, which are not shown for clarity.
The springs are pinned to the masses or the walls allowing for their motions as the masses displace in both directions. As we will see, unlike the nonlinear Euler- a y x Figure 28: Schematic of the mass-spring-grid system with geometric nonlinearity Bernoulli beam, the space coordinates of the mass-spring-grid system can be model analytically. This removes any numerical instability related to meshing the space coordinates. Therefore, studying this system allows us to better investigate the performance of the "persistent MOR" methodology. In case the system is forced only in x-direction, Eq. (35) is sufficient to describe it. However, any small deviation from x-directional oscillation will cause geometric nonlinearity in the motion. For the purpose of this section, we only excite this system in y-direction. The governing differential equations for the i-th mass are as follows: where a is the initial distance between the masses, l is the free length of the springs, l i = (a + ∆x i ) 2 + (∆y i ) 2 , and ∆x i and ∆y i and are given as follows: We should note that here only linear damping is considered for the system.
The state-space vector to model this system is a vector of 4n variables defined Thus,Eq. (46) and Eq. (47) can be rewritten as follows: where 0, I, M, K, and C are n × n zero, identity, mass, stiffness, and damping matrices, respectively.

Reduced Order Models
For the mass-spring-grid system consisting of twenty masses, specifying the following parameters will result in a rich dynamical behavior: The subspace robustness for POD-based MOR of this system is close to unity for k = 1, . . . , 4 and drops to 0.6 at k = 5, as shown in Fig. 29. For SOD, subspace robustness starts at near 0.85 for k = 1, monotonically increases with the increase in the dimension and saturates at 1 near k = 20. Also, three-and higherdimensional POD-based ROMs are dynamically consistent, while for SOD five and higher dimensional subspaces are dynamically consistent. The importance of subspace robustness and dynamical consistency metrics for identifying the optimal MOR subspace is reflected in Fig. 30, where POD-based ROMs lose their stability as subspace robustness drops after k = 5. While five-dimensional ROM is still stable, for the k = 6, 7, or 8 it loses its stability. The importance of monotonically increasing subspace robustness for SOD-based ROMs is illustrated in Fig. 31.
These ROMs become and remain stable as the robustness metric approaches unity and stays there.

Introduction
In the previous chapter we showed that SOD outperformed POD in obtaining low-dimensional ROMs. However, the data matrices associated with many dynamical systems are ill-conditioned. This may destabilize the low-dimensional ROMs or reduce their performance. For instance, in Section 4.2.1 we observed that some low-dimensional POD-and SOD-based ROMs were not stable for the nonlinear Euler-Bernoulli beam [1].
The data matrices from the full state-space models are composed of measured data corresponding to position and velocity state variables. Position and velocity data have a lower condition number, i.e. they can form better conditioned matrices separately. Therefore, we propose to perform the multivariate analysis separately to identify the so-called position and velocity modes. Following the identification of these modes, we can put them together to span the modal subspace for MOR.
In this chapter the mathematical formulation of the separated POD-as well as the separated SOD-based MOR methods are provided. We apply these methods to the nonlinear spring coupling system described in Section 4.1.

Separated Multivariate Analysis for Reduced Order Models
The full data matrices in many dynamical systems are ill-conditioned. Therefore, extracting the modes out of these data using the aforementioned multivariate analysis methods may result in noisy modes. These modes have many unnecessary sign changes and may not account for the true modal structure of the system. Here we propose to do the multivariate analysis on position and velocity data matrices separately, since they have a lower condition number. We consider the position and the velocity data matrices X and V, as described in section 1.3.1. The separated multivariate analysis as its name suggests will be done on X and V separately. The solution of separated POD analysis on position data matrix is given by: Likewise, the solution of separated POD analysis on velocity data matrix is given by: In the above equations, Σ xx and Σ vv are autocovariance matrices of position and velocity data. Also, the superscript x and v denote that the corresponding scalar or vector is related to the position and velocity data matrices, respectively. The After choosing l dominant position modes given in the matrix form Φ x l as well as r dominant velocity modes given in the matrix form Φ v r , they can be put together to, at last, yield the modal matrix: where k = l + r is the dimension of the modal matrix.
For separated SOD analysis, we need to solve the following generalized eigenvalue problems on position and velocity data and their time derivatives. For position data, one has: similarly for velocity data, one obtain: in which Σẋẋ and Σvv are autocovariance matrices of velocity and acceleration  Figure 34: Subspace robustness metric for separated SOD data, respectively. Similar to POD, the modal matrix Ψ will be composed of l-dimensional Ψ x l and r-dimensional Ψ v r . Thus, one get: (56)

Numerical Example
As a numerical example, we consider the system with nonlinear spring coupling This will result in a rich dynamic response with two equilibrium points. Subjected to harmonic forcing with Ω = 9.97 Hz, which is close to the frequency of the first linear normal mode, the bifurcation diagram is plotted using the full scale model of the system as shown in Fig. 32. We aim to reproduce these results using ROMs. While any type of data obtained from the system oscillations could be used for multivariate analysis, the data from random excitation would be the best choice for MOR subspace identification. With random forcing, we are more likely to explore nearly all phase space of the system and awake most of its dominant modes. Therefore, to get the data matrices, the white noise is applied to the system for a total time of 272 seconds. The response of the system is obtained using Runge-Kutta method with the sampling time of 2.178 × 10 −3 sec. Thus, a total of 125000 data will be recorded to form position and velocity data matrices. These data will be processed using POD, SOD, separated POD and separated SOD as the multivariate methods.

Results
We excited the system by the white noise with 12 different amplitudes to get 12 independent cases. For each case, we record the data matrices and apply separated multivariate analysis. The resultant modes were mixed using Eq. (53) and Eq. (56). The obtained modes are then used for ROMs.
For optimal selection of the modes, we apply the subspace robustness metric to position and velocity modal matrices obtained from separated POD and SOD.
The results are depicted in Fig. 33 and Fig. 34, for POD and SOD, respectively. a reasonable accuracy, the four-dimensional full SOD-based model has a slightly better performance than the seven-dimensional full POD-based model. This is consistent with the previous work done on SOD-based MOR [1,2], where SODbased ROMs were shown to capture the dynamics of the nonlinear systems using lower dimensional models.
An interesting result is the effect of separating the data matrices on stabilizing the low-dimensional POD-based ROMs. Now not only the separated POD-based ROMs are stable for all dimensions, but also the performance of the models is improved. This is illustrated in Fig. 35 the results. In some cases the results do not change much and in some cases, e.g. q 0 = 0.12, the performance of the ROM significantly increases. In this particular case, the ROM is able to distinguish similar, neighbor trajectories.
With a more than 50 times faster ROM, we are able to reproduce the bifurcation diagrams in more details. Fig. 37 shows the bifurcation diagrams obtained for the four methods. It can be observed that there is a slight improvement in the diagram of the separated POD-based MOR compared to that of the full POD-based MOR. For instance, there is better match for the forcing amplitude in the range of 0.2 to 0.25 for the separated POD diagram with that of the full-scale model in Fig. 32. Although this improvement is not significant, it shows the improving nature of the separated POD-based reduced order modeling.
The bifurcation diagram of four-dimensional full SOD-based ROM is consis-

CHAPTER 6
Persistent ROMs for control Systems

Introduction
A high-fidelity mathematical model is essential to control a complex nonlinear dynamical system. These models are often high-dimensional, which means that complex differential equations are needed to describe them. Therefore, in many cases, they may not be computationally tractable. This makes the real-time control difficult to implement. A ROM of a complex system can result in a computationally tractable, accurate model for the control system [1].
Computationally complex dynamical systems usually evolve on a lowerdimensional curved nonlinear manifold embedded in a higher dimensional state space of the system. Geometric structures of nonlinear manifolds have not been extensively incorporated in nonlinear control theory since identification of highdimensional manifold is difficult [2,3,4]. Also, even if we overcome this problem, the stability and accuracy of the reduced model is still guaranteed only for a small range of operating conditions or modal parameters [4].
In this chapter, we use SOD [5,6,7] as a new tool for MOR for nonlinear control systems [8]. Our method is categorized under Galerkin projection based reduced order modeling which projects the high-dimensional nonlinear system onto an appropriate linear subspace to yield a lower-dimensional system. We also evaluate the persistency of the identified linear subspaces. Since persistent linear subspaces are robust to the changes in system's operating conditions, they expand regions within the system's state space in which the ROMs are valid. We aim to obtain a persistent ROM which allows the control system to globally operate within a region of interest.
Projection onto the linear subspace does not negate the nonlinearity of the Figure 38: Schematic of the nonlinear control system original system [9]. While the resultant ROM for the control system is still nonlinear, its corresponding state is lower dimensional which makes the control system computationally manageable.
For nonlinear control systems, however, we examine the output of the persistent ROM for a given input in comparison to the output of the full-scale control model. For the input we use a set of impulse functions as random input. This approach has two advantages: (1) under random input it would be difficult to stay in a limited region of the space; and (2) random input imitates the non-deterministic impulses generated by the control scheme as inputs to the system.
For the purpose of this work we consider the model presented in [1]. We describe and apply SOD as a new reduced order modeling method for nonlinear control systems. We also formalize the subspace robustness as a metric to identify the persistent subspaces for reduced order control models in such a way that they are globally valid for a range of the system's energy. Finally, the developed methodology in this chapter will be tested using numerical simulations of a nonlinear control system.

Background and Prior Work
Within the realm of complex dynamical systems, reduced order modeling is being extensively used to reduce the redundant computations and data storage requiremenst [10,7,11,12,13,14,15]. The methodologies for MOR of complex dynamical systems were discussed in Chapter 1.
The research on MOR of control systems is extensive. It includes well understood, and established theories and methodologies for reduction of linear control systems. Examples of these methods are POD, used for instance to design control systems for PDEs [16,17] and optimal control of fluids [18], Hankel norm approximation [19,20,21], and balanced truncation [22] which was proposed by Moore [23]. The reader may review other methods for MOR for linear control system in Refs. [22,24,25].
Model reduction of nonlinear control systems is not as well understood as for linear systems. For example, POD is being frequently used [26], however, it suffers from some limitations that are discussed in [27]: POD-based models are very sensitive to the data used [9] and may become unstable even near stable equilibrium points [28]. Another method is balanced truncation which is developed for nonlinear control system in two distinct approaches: one is based on energy function used in the works by Scherpen [29,30,31,32] and the other is proposed by Lall based on empirical balanced truncation [1].

Model Reduction Using Galerkin Projection
We consider a nonlinear control system in the form: where y ∈ R 2n is state vector of the system, n is number of degrees-of-freedom, t is time, f : R 2n ×R p → R 2n is a nonlinear flow function describing the dynamics of the system, u(t) ∈ R p is the input to the system, and z(t) ∈ R w is the system output or the state vector which is based on the desired observation, h : R 2n → R w . The goal of the control system is to control the output z(t), however, if the system is largescale or highly nonlinear, we will aim to obtain a reduced order nonlinear control model. A reduced order control model is easier to implement and is essential for a real-time and accurate control.
Galerkin projection based MOR methods are based on transforming the 2ndimensional state vector y to a k-dimensional state vector q, given that k < 2n.
The transformation is performed by a full-rank projection matrix P k ∈ R 2n×k in the form q = P † k y, with (.) † defined as the pseudoinverse of (.), to yield the reduced order model:q (t) = P † k f (P k q(t), u(t)) , Matrix P represents a description of the modal space of a dynamical system.
Matrix P k is the k-dimensional modal subspace formed by k dominant modes of the modal space. While it can be analytically obtained for linear dynamical systems using linear normal modes theory, another method to obtain P, regardless of system's linearity or nonlinearity, is using multivariate analysis of its response.
Multivariate analysis is applied to the data matrices from the full model simulations or experiments. In this work, all the data matrices are obtained from simulations.
We present an example of a nonlinear control system derived from the work by Lall et al. [1] in which they developed the balanced truncation method for nonlinear control systems. In next section, we obtain persistent reduced-order control models for this system. In this section, we model the system adopted from [1]. The system, shown in We obtain the governing differential equation of the system using the Lagrange's equation: where V and T are potential and kinetic energy, and F i is the generalized forcing term. Now we consider y = [θ 1 , . . . , θ 5 ,θ 1 , . . . ,θ 5 ] T to be the state vector. By substituting the state vector in the equations of motion, we obtain its state space form: in which M y(t) is the time-varying mass matrix and L is the matrix of the linear terms. Both are given in Appendix A. Also, f n is the vector of the nonlinear terms and u(t) is the single input to the system. The output of the system is defined as the horizontal position of the tip of the 5th link and is to be controlled.
We simulate Eq. (60) as a full-scale model of the control system using harmonic excitation, u(t) = f 0 sin ωt.

Reduced Order Nonlinear Control System
In order to construct ROM, we first randomly or stochastically drive the fullscale model to collect the required data from s different simulations. We use multivariate analysis to obtain the modal structure from each simulation. Then we apply the subspace robustness to the modal structures to select the dimension of the persistent subspace that can be used for the global reduced model, as described in Section 2.1, Chapter 2. Using the obtained subspace we construct the model While any record of the system states can be used as data for multivariate analysis, we use random excitation as the system input and collect the response of the system in the data matrices. This way we ensure that all neighbors of data points within the space of the system has been covered and that the modal structure we obtain from the analysis of data will be a better representation of the important dynamical characteristics of the system. Since we aim to build a relatively global reduced order control system which is valid for a range of energy levels, we do 12 simulations with different energy levels. To impose the changes in the energy, we only change the amplitude of the excitation while keeping the frequency content similar for all cases.
The link system has a linear modal frequency range up to 3 Hz. We limit the frequency of the random excitation to 5 Hz to assure that all linear modes are covered while data are not contaminated by noise. We select 12 equally distributed choices of the random forcing amplitude from the range of 0.1 ≤ q 0 ≤ 3. We excite the link system by the random forcing to obtain 12 data matrices Y 1 , Y 2 , . . . , Y 12 .
We identify the modal structure of each data set using POD and SOD. We calculate the subspace robustness of POD and SOD modes using Eq. (30). Fig. 40 shows the subspace robustness of POD and SOD for each dimension. The POD subspace robustness for k = 1 is very close to unity which means that the first dominant POMs from all the simulations are linearly dependent. The POD subspace robustness is also close to one for k = 7, 8 and 10. On the other hand, the SOD subspace robustness is always close to one. A subspace robustness closer to one suggests few changes occur in subspaces from different simulation. This means that there is less leakage to the higher dimensional subspaces and the subspace is persistent to changes in system's energy level. Therefore, SOD subspaces, are more persistent compared to those of POD.
Following the identification of dimension for which the subspaces are robust and persistent, in order to obtain the global reduced order control model, we combine all the data matrices together to obtain a large response matrix, Y, as follows: We extract the corresponding POMs and SOMs, as the modal space given by P, and its k-dimensional representation of the k dominant modes given by P k . In case k is the dimension of persistent subspace, we expect P k via Eq. (58) to result in a persistent ROM within the range of energies of the nonlinear control system.
Please note that for POD, POMs (denoted by Φ) are orthogonal and thus, P k = Φ k and P † k = Φ T k . For SOD, however, SOMs and SPMs are bi-orthogonal (Φ T Ψ = I), thus, P k = Φ k and P † k = Ψ T k .
Also, from matrix Y we can extract POVs and SOVs to measure the dominance of the modes. Fig. 41 depicts the POVs and SOVs. We look for the drops in their values in order to identify the low-dimensional control models. There is no significant drop in the POVs for lower k values as we observe that they gradually decrease. The POV after k = 8 drops more drastically. However, SOVs come in pairs and the drops are distinguishable. A clear drops occur at k = 2, k = 4, and k = 6. Yet, we do not expect a good control model for k = 2 from SOD since the higher dimensional modes still have a significant SOV.
The full scale nonlinear control system will be controlled by a sequence of unit inputs. The proper choice of input merely depends on the design on the controller and the control method. Therefore, a good ROM for nonlinear control system is expected to mimic the output of the full scale model excited by a random input since we have no further knowledge about the specific controller.
We generate a filtered random input with the frequency content up to 5 Hz.
We excite both full-scale and ROM control systems by this input and compare their outputs, which are in this case the horizontal positions of the 5th link. For SOD, all the ROMs except for the three-and five-dimensional are stable, although the lowest dimensional ROM which provides good results is four-dimensional. In Fig. 42 This confirms the results of the subspace robustness metric for the POD subspace.
In Fig. 44, we show the computation speed of the reduced control models and compare it to the full scale model of the control system. For both POD and SOD, the computation speeds of the unstable models are estimated by interpolation.
We observe that the eight-dimensional POD model computation time is close to the full scale control model, while its performance is not as good. Nine-and tendimensional models are even slower than the full-scale model. We note that the ten-dimensional POD model is just a POD realization of the full-scale model with the same dimension. On the other hand, the four-dimensional SOD control model is more than 6 times faster than the full-scale model of the control system.
We also notice that the computation time of the SOD models, unlike POD, increases almost linearly. More interestingly, even a 10-dimensional SOD model, which has the same dimension as the full-scale model, is about twice faster, while it provides a perfect tracking of the output. We did not expect to get these results, however, at this point we speculate that since we used MATLAB ode45 solver for numerical integration which adjusts the sampling rate based on the smoothness of the flow, and since SOD provides a smoother realization of the full-scale model of the control system, it results in speeding up the integrations. We will further investigate this effect in our future work.
[3] R. J. Kuether, M. R. Brake, and M. S. Allen, "Evaluating convergence of reduced order models using nonlinear normal modes," in Model Validation  In particular, SOD-based ROMs outperformed the POD-based ones in terms of the stability and robustness of the model. Also, the obtained persistent ROMs could be successfully used to study the dynamics of computationally expensive complex models for a relatively wide range of parameters and conditions. In Chapter 5, we proposed the separated multivariate analysis using POD and SOD to deal with the rank-deficiency of nonlinear dynamical systems. Following the mathematical formulation of the separated POD-and SOD-based MOR, a nonlinear mass-spring-damper was studied. The full-scale state-space model of the system was obtained and the data matrices were extracted. The performance of four methods for MOR including POD, SOD, separated POD and separated SOD was investigated in terms of stability and accuracy. The results confirmed that SOD was able to capture the dynamics of the system in a lower dimensional space.
Separated multivariate analysis was shown to improve the accuracy and stability of the POD-and SOD-based MOR.
In Chapter 6, a new approaches for MOR of nonlinear control systems was presented. An example of a system with five inverted links was used to examine our approach. The modal subspaces which were identified using projection based reduced order modeling methods were shown to be dependent on the system's energy. The subspace robustness metric was proposed to obtain robust and persistent reduced order control models. These models were aimed to be valid for a range of the system's energy. The developed metric was used to evaluate for POD-and SOD-based subspaces. POD subspaces were shown persistent only for the high dimensional models. SOD subspaces were persistent for all the dimensions. The resultant reduced order control models were tested using different random inputs.
Low-dimensional POD-based ROMs were not stable and the high dimensional ones were not as accurate as the low-dimensional SOD ROMs. A four-dimensional SOD ROM closely tracked the output of the nonlinear control system to different random inputs. These results were consistent with the subspace robustness metric. The accurate SOD ROMs were shown to be six times faster than the full-scale model. These ROMs outperformed the best POD ROM, which was not significantly faster than the full-scale control system. Also, we showed that the smoothing effect of SOD may speed up the full-scale model simulations, as we observed that the 10-dimensional full-scale SOD model was as accurate as, but two times faster than the original full-scale system.