Journal of Control Science and Engineering
Volume 2008 (2008), Article ID 154956, 11 pages
doi:10.1155/2008/154956
Research Article

Empirical Reduced-Order Modeling for Boundary Feedback Flow Control

Seddik M. Djouadi,1 R. Chris Camphouse,2 and James H. Myatt3

1Department of Electrical Engineering and Computer Science, University of Tennessee, 1508 Middle Drive, Knoxville, TN 37996, USA
2Performance Assessment and Decision Analysis Department, Carlsbad Program Group, Sandia National Laboratories, 4100 National Parks Highway MS 1395, Carlsbad, NM 88220, USA
3Air Force Research Laboratory, Wright-Patterson Air Force Base, Dayton, OH 45433, USA

Received 27 December 2007; Revised 4 November 2008; Accepted 24 December 2008

Academic Editor: Onur Toker

Copyright © 2008 Seddik M. Djouadi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

This paper deals with the practical and theoretical implications of model reduction for aerodynamic flow-based control problems. Various aspects of model reduction are discussed that apply to partial differential equation- (PDE-) based models in general. Specifically, the proper orthogonal decomposition (POD) of a high dimension system as well as frequency domain identification methods are discussed for initial model construction. Projections on the POD basis give a nonlinear Galerkin model. Then, a model reduction method based on empirical balanced truncation is developed and applied to the Galerkin model. The rationale for doing so is that linear subspace approximations to exact submanifolds associated with nonlinear controllability and observability require only standard matrix manipulations utilizing simulation/experimental data. The proposed method uses a chirp signal as input to produce the output in the eigensystem realization algorithm (ERA). This method estimates the system's Markov parameters that accurately reproduce the output. Balanced truncation is used to show that model reduction is still effective on ERA produced approximated systems. The method is applied to a prototype convective flow on obstacle geometry. An 𝐻 ∞ feedback flow controller is designed based on the reduced model to achieve tracking and then applied to the full-order model with excellent performance.

1. Introduction

Recently there has been significant interest in model reduction for the purpose of control design [1–9]. One such application of reduced-order modeling is control design in the context of aerodynamic flow. Aerodynamic flow control is a research area of great interest to the air force and the fluid mechanics community. Presently considerable research efforts are working with feedback control law design for systems described by PDEs that need a very large number of states to accurately simulate their characteristics. However, recent advances in the design of actuators and sensors can be leveraged for better system control only if the control design methods provide a reliable low-order controller [10]. Additionally, simulation, and experimental diagnostics are making applications such as the suppression of acoustic tones in cavities, separation control for high lift, and trajectory control without the need to move hinged surfaces a possibility [11]. However, these applications require the integration of feedback control because of the need for robustness to flight condition and vehicle attitude, precision tracking, overcoming low-fidelity models, or moving a system away from a stable solution or limit cycle as efficiently as possible [11].

The traditional systematic development of feedback control laws for these systems subject to a large number of states is currently an intractable problem. Feedback control strategies offer the possibility to improve performance and reduce control power through the control of unstable structures in the flow field. Reduced models are important for the design of feedback control laws, which rely on models that capture the relevant dynamics of the input-output system and are amenable to control design.

Unfortunately, it is very difficult to create models that capture the relevant dynamics of the input-output system. For example, computational fluid dynamics simulations can provide good solutions to a discretized version of the Navier-Stokes equation [1, 12]. However, accurate simulations for simple shapes such as two-dimensional airfoils, or complex shapes, such as a full vehicle, require several thousands to millions of states. Therefore, the simulation results are not directly useful for control design [11]. Complexity in the model is a legitimate need. The large number of states is necessary to capture important flow features that occur at extremely small spatial scales. Although these small flow features might seem insignificant, if they are not captured, it is not possible to analyze if they are necessary in securing the closed-loop system’s overall stability [10].

POD has been extensively investigated in distributed parameters systems due to its order reduction capability [7–9, 13, 14], and balanced truncation, which is a simple yet efficient model reduction technique widely used in reducing model orders of high-order linear systems [15–17]. POD models of only a few dozen states can often accurately capture the input-output behavior of systems that have full-order system models of thousands of states [11]. In addition to using the POD method in conjunction with model reduction techniques, the idea of using empirical gramians is growing in popularity for use in an approximate balanced truncation [4, 18–20]. Further, some work has been done on finding nonlinear empirical Gramians for balanced truncation [20, 21].

In fluid flow configurations it is not uncommon for discretized flow models to describe thousands to millions of state variables, for example, if one uses a linear quadratic regulator (LQR) control formulation, roughly 1 0 1 2 Riccati unknowns need to be calculated for a discretized flow model describing 1 0 6 states. Existing computing power and computational algorithms are not capable of solving an LQR problem of such large dimension. For dynamical models that are very large scale, such as those describing fluid flow configurations, it is apparent that the order of the system must be reduced prior to control law design [22]. This prevents us from using closed-loop model reduction strategy wherein the system is part of a closed-loop system with a controller [23].

In the area of fluid mechanics controls must often be fixed to the boundary of the problem geometry. For example, control of flow separation over an airfoil requires that actuation and sensing be done on the airfoil surface [11]. The problem geometry used for this project is one example of a case where control is restricted to the boundaries by physical necessity.

The paper is organized as follows. In Section 2, we introduce a prototype flow problem geometry that is used to apply the proposed order reduction techniques. Section 3 introduces the eigensystem realization algorithm (ERA) to identify the Markov parameters of the system, and as a product of the empirical gramians. Section 4 introduces empirical balanced truncation. This method is based on approximate (empirical) controllability and observability gramians and uses only a single simulation/experimental test. In Section 5, empirical balanced truncation and the ERA algorithm are applied to the Galerkin model, and numerical results are provided to show the effectiveness of the proposed method. In Section 6, an 𝐻 ∞ controller based on the empirical reduced model and which achieves tracking is discussed. The responses of the controlled closed loop on the full-order model are presented and show that the 𝐻 ∞ controller achieves good tracking performances despite being designed on a much lower model than the original model. Section 7 contains concluding remarks.

2. Problem Geometry

The specific problem geometry considered is shown in Figure 1. The idea and methods presented here could be modified to apply to a different geometry or obstacle shape. The problem statement with its corresponding boundary conditions and governing equations was taken from [10]. A realistic example of this geometry in an aerodynamic application would be a payload hatch open during flight with actuator control only on the boundary. Let Ω g a p be the region defined by [ 𝑎 1 , 𝑎 2 ] × [ 𝑏 1 , 𝑏 2 ] . Let Ω f u l l be the region defined by ( 𝑎 0 , 𝑎 e n d ) × ( 𝑏 0 , 𝑏 e n d ) . Then the problem domain is given by Ω = Ω f u l l / Ω g a p . In this problem setup, Ω g a p is an obstacle.

154956.fig.001
Figure 1: Problem geometry.

The system dynamics that act within the problem domain are described by the two-dimensional (2D) Burgers’ equation [10] 𝜕 = 1 𝜕 𝑡 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) + ∇ ⋅ 𝐹 ( 𝑤 ) 𝑟 ( 𝜕 2 𝜕 𝑥 2 𝜕 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) + 2 𝜕 𝑦 2 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) ) , ( 1 ) where the form of 𝐹 ( 𝑤 ) is 𝑐 𝐹 ( 𝑤 ) = 1 𝑤 2 ( 𝑡 , 𝑥 , 𝑦 ) 2 𝑐 2 𝑤 2 ( 𝑡 , 𝑥 , 𝑦 ) 2 𝑇 . ( 2 ) In this case, the value for 𝑐 1 is equal to 1 and 𝑐 2 is equal to 0. The value 𝑟 is similar to the Reynolds number used in the Navier-Stokes equation. This parameter controls how much nonlinearity is present in the problem. The value used is 300, a small “Reynolds number,” but it still allows for the nonlinearity to show in the problem.

Dirichlet boundary conditions located on the obstacle top and bottom are denoted by Γ t o p a n d Γ b o t t o m . A Dirichlet boundary condition is a first-type boundary condition that specifies the values of the solution defined by 𝑓 ( 𝑥 ) on a domain boundary [24]. The form of the boundary condition is 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) = 𝑓 ( 𝑡 , 𝑥 , 𝑦 ) ∀ 𝑥 , 𝑦 ∈ 𝜕 Ω . ( 3 ) The boundary conditions on the top and bottom are described by the following: 𝑤 𝑡 , Γ b o t t o m = 𝑢 b o t t o m ( 𝑡 ) Ψ b o t t o m 𝑤 ( 𝑥 ) , 𝑡 , Γ t o p = 𝑢 t o p ( 𝑡 ) Ψ t o p ( 𝑥 ) . ( 4 ) Here 𝑢 t o p ( 𝑡 ) and 𝑢 b o t t o m ( 𝑡 ) are the control inputs on the top and bottom boundaries, respectively; the spatial functions Ψ t o p ( 𝑥 ) and Ψ b o t t o m ( 𝑥 ) describe the spatial effect that the controls have on the top and bottom boundaries.

The boundary condition on the airflow intake side is 𝑤 𝑡 , Γ i n = 𝑓 ( 𝑦 ) ( 5 ) and it is parabolic in nature. The airflow outtake side has a Neumann boundary condition that has the form [24] 𝜕 𝑤 𝜕 𝑥 𝑡 , Γ o u t = 0 . ( 6 ) On all of the remaining boundaries of Ω , 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) is set equal to 0 for all values of t. Finally, the initial conditions for the interior are given by 𝑤 ( 0 , 𝑥 , 𝑦 ) = 𝑤 0 ( 𝑥 , 𝑦 ) ∈ 𝐿 2 ( Ω ) . ( 7 ) A numerical solution was found by simulation using a uniformly spaced grid. The resulting system model contains a little more than 2000 states [11]. The POD model construction is based on the total energy captured. A condition that 99.9% of system energy must be captured was used for determining how many system modes were retained. This condition was met by a 40 POD basis. Although this is a major reduction from the numerical solution, it will be shown that important system dynamics can be retained with even lower state number system models.

The general approach of POD is to construct a series of solution “snapshots.” These snapshots are generated by numerical simulations of the governing system equation(s) with a variety of input equations [6, 25]. The snapshots are needed to generate the correlation matrix that is used to find the POD basis. It is important to choose relevant input signals for the numerically simulated system. Further, these inputs should be similar to the expected inputs of the real system. The inputs used for system identification are of the form [10] 𝑢 b o t t o m ( 𝑡 ) = 𝛽 s i n 0 . 2 5 𝑡 2 , 𝑢 t o p 𝑢 ( 𝑡 ) = 0 , b o t t o m ( 𝑡 ) = 0 , 𝑢 t o p ( 𝑡 ) = 𝛽 s i n 0 . 2 5 𝑡 2 , 𝑢 b o t t o m ( 𝑡 ) = 𝛽 s i n 0 . 2 5 𝑡 2 , 𝑢 t o p ( 𝑡 ) = 𝛽 s i n 0 . 2 5 𝑡 2 , ( 8 ) where the values for 𝛽 are − 3 , − 2 and − 1 and the range for 𝑡 is 0 to 10 seconds with a sample every 50 milliseconds. The squelch signal for all three values of 𝛽 is shown in Figure 2.

154956.fig.002
Figure 2: Test inputs used to generate the snapshots.

The numerical simulation was performed to create the ensemble of solution snapshots { 𝑆 𝑘 ( 𝑥 , 𝑦 ) } 𝑀 𝑘 = 1 [10]. The value for M (the number of snapshots) must be greater than the number of modes that one will choose for the approximated system model. For a good representation M should be much larger than the desired size for the POD basis [10].

The solution to the PDE is assumed be finite energy, that is, belongs to the Hilbert space 𝐿 2 ( [ 0 , 𝑇 ) × Ω ) . The solution can be approximated as 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) ≈ 𝑛 𝑘 = 1 𝛼 𝑘 ( 𝑡 ) 𝜙 𝑘 ( 𝑥 , 𝑦 ) , ( 9 ) where the 𝛼 𝑘 's are time varying coefficients that multiply the spatial functions 𝜙 𝑘 's and 𝛼 𝑘 ∈ 𝐿 2 [ 0 , 𝑇 ) , 𝜙 𝑘 ( 𝑥 , 𝑦 ) ∈ 𝐿 2 ( Ω ) , ( 1 0 ) where 𝐿 2 ( [ 0 , 𝑇 ) ) and 𝐿 2 ( Ω ) are the standard Hilbert spaces of absolutely square integrable functions defined, respectively, on the time interval [ 0 , 𝑇 ) and spatial domain Ω . The approximation (9) can be as accurate as desired since the tensor space 𝐿 2 [ 0 , 𝑇 ) ⊗ 𝐿 2 ( Ω ) ∶ = { 𝑛 𝑘 = 1 𝛼 𝑘 ( 𝑡 ) 𝜙 𝑘 ( 𝑥 , 𝑦 ) , 𝛼 𝑘 ( 𝑡 ) ∈ 𝐿 2 , 𝜙 [ 0 , 𝑇 ) 𝑘 ( 𝑥 , 𝑦 ) ∈ 𝐿 2 ( Ω ) , ∀ 𝑛 i n t e g e r } ( 1 1 ) is dense in 𝐿 2 ( [ 0 , 𝑇 ) × Ω ) [18].

Any basis for 𝐿 2 ( Ω ) can be used to construct the approximation of the solution 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) . Here we use the POD basis { 𝜙 𝑘 } since it is optimal in the following sense [2]: 𝜇 𝑛 ∶ = i n f 𝛼 𝑘 ∈ 𝐿 2 𝜙 ( [ 0 , 𝑇 ) ) , 𝑘 ( 𝑥 , 𝑦 ) ∈ 𝐿 2 ( Ω ) ‖ 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) − 𝑛 𝑘 = 1 𝛼 𝑘 ( 𝑡 ) 𝜙 𝑘 ( 𝑥 , 𝑦 ) ‖ 2 . 𝜇 𝑛 → 𝑛 → ∞ 0 . ( 1 2 ) In [2], it is shown that solving the above optimization reduces to the usual POD optimality in the average kinetic energy sense (15) discussed below [25].

The POD basis function 𝜑 𝑖 is chosen to maximize the average projection of the member 𝑢 𝑖 onto 𝜑 𝑖 suitably normalized [25] m a x 𝜑 𝑖 ∈ 𝐿 2 ( Ω ) | | | 𝑢 𝑖 , 𝜑 𝑖 | | | 2 ‖ ‖ 𝜑 𝑖 ‖ ‖ 2 2 . ( 1 3 ) To construct numerically the POD basis { 𝜑 𝑗 } , we build the correlation matrix L, of size 𝑀 × 𝑀 , composed of the inner products of the snapshots. In 𝐿 2 ( Ω ) we have 𝐿 𝑘 , 𝑗 = 𝑆 𝑘 , 𝑆 𝑗 𝑆 , ( 1 4 ) 𝑘 , 𝑆 𝑗 = Ω 𝑆 𝑘 𝑆 ∗ 𝑗 𝑑 𝑥 𝑑 𝑦 , ( 1 5 ) ∗ where { 𝜆 1 , 𝜆 2 , … , 𝜆 𝑛 } denotes complex conjugate transpose.

A spectral decomposition of the matrix L is performed. The n largest eigenvalues { 𝑣 1 , 𝑣 2 , … , 𝑣 𝑛 } . of the matrix L are found and placed in descending order. Then the set of eigenvectors are identified to be { 𝜑 𝑗 }

The resulting orthonormal POD basis 𝜆 𝑘 ‖ ‖ 𝑣 𝑘 ‖ ‖ 2 = 1 , f o r 𝑘 = 1 , 2 , … , 𝑛 . ( 1 6 ) of dimension n can be constructed using the information found from the correlation matrix L. First, the eigenvectors of L are weighted by their corresponding eigenvalues and normalized according to [10] 𝜑 𝑘 ( 𝑥 , 𝑦 ) = 𝑀 𝑗 = 1 𝑣 𝑘 , 𝑗 𝑆 𝑗 ( 𝑥 , 𝑦 ) ( 1 7 ) Then, the POD basis set is formed according to 𝑣 𝑘 , 𝑗 with 𝑣 𝑘 . being the jth component of the eigenvector 𝑛 𝜑 𝑘 Solving (17) gives { 𝛼 𝑘 } 's, which constitute the POD basis of dimension n.

The governing equation is projected onto the POD basis. The projection is accomplished via a Galerkin type projection and results in a system of ordinary differential equations (ODEs). The Galerkin projection results in only a weak solution to the PDE. However, this weak solution with finite difference approximations of the boundary conditions eventually leads to a nonlinear temporal model for the temporal or POD coefficients Ω 𝜕 𝜕 𝑡 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) 𝜑 𝑘 = 1 ( 𝑥 , 𝑦 ) 𝑑 𝑥 R e 𝜕 Ω 𝜑 ∇ 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) ⋅ 𝑛 𝑘 − ( 𝑥 , 𝑦 ) 𝑑 𝐴 ( 𝑥 ) Ω ∇ 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) ⋅ ∇ 𝜑 𝑘 − ( 𝑥 , 𝑦 ) 𝑑 𝑥 𝜕 Ω 𝜑 𝐹 ( 𝑤 ) ⋅ 𝑛 𝑘 − ( 𝑥 , 𝑦 ) 𝑑 𝐴 ( 𝑥 ) Ω 𝐹 ⋅ ∇ 𝜑 𝑘 , ( 𝑥 , 𝑦 ) 𝑑 𝑥 ( 1 8 ) [10].

Projecting (1) onto the POD basis yields [10] 𝜕 Ω 𝜑 ∇ 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) ⋅ 𝑛 𝑘 = ( 𝑥 , 𝑦 ) 𝑑 𝐴 ( 𝑥 ) 𝑎 2 𝑎 1 𝜕 𝑤 𝜕 𝑦 𝑡 , 𝑥 , 𝑏 1 ) 𝜑 𝑘 𝑥 , 𝑏 1 − 𝑑 𝑥 𝑎 2 𝑎 1 𝜕 𝑤 𝜕 𝑦 𝑡 , 𝑥 , 𝑏 2 𝜑 𝑘 𝑥 , 𝑏 2 − 𝑑 𝑥 𝑏 e n d 𝑏 0 𝜕 𝑤 𝜕 𝑥 𝑡 , 𝑎 0 𝜑 , 𝑦 𝑘 𝑎 0 , 𝑦 𝑑 𝑦 . ( 1 9 ) where the first term on the right-hand side is 𝑏 0 The Neumann boundary condition forces the portion of the boundary integral over 𝑏 e n d to 𝑎 e n d along 𝑏 e n d 𝑏 0 𝜕 𝑤 𝜕 𝑥 𝑡 , 𝑎 e n d 𝜑 , 𝑦 𝑘 𝑎 e n d , 𝑦 𝑑 𝑦 = 0 . ( 2 0 ) to be 0, that is, 𝜕 Ω 𝜑 𝐹 ( 𝑤 ) ⋅ 𝑛 𝑘 = 1 ( 𝑥 , 𝑦 ) 𝑑 𝐴 ( 𝑥 ) 2 𝑏 e n d 𝑏 0 𝑤 𝑡 , 𝑎 e n d , 𝑦 2 𝜑 𝑘 𝑎 e n d , 𝑦 − 𝑓 ( 𝑦 ) 2 𝜑 𝑘 𝑎 0 , 𝑦 𝑑 𝑦 . ( 2 1 ) The second boundary integral is decomposed as follows 𝜕 𝑤 𝜕 𝑦 𝑡 , 𝑥 , 𝑏 1 ≈ 𝑢 b o t t o m ( 𝑡 ) Ψ b o t t o m ( 𝑥 ) − 𝑤 𝑡 , 𝑥 , 𝑏 1 − ℎ ℎ , 𝜕 𝑤 𝜕 𝑦 𝑡 , 𝑥 , 𝑏 2 ≈ 𝑤 𝑡 , 𝑥 , 𝑏 2 + ℎ − 𝑢 t o p ( 𝑡 ) Ψ t o p ( 𝑥 ) ℎ , 𝜕 𝑤 𝜕 𝑥 𝑡 , 𝑎 0 ≈ 𝑤 , 𝑦 𝑡 , 𝑎 0 + ℎ , 𝑦 − 𝑓 ( ℎ ) ℎ . ( 2 2 ) This solution does not explicitly include the control inputs or boundary condition information into the governing equation. In order to do so an approximation of the partial derivatives is carried out including the control inputs and the boundary data. If h denotes the step size between the points on the uniform Cartesian grid used for the finite-difference solution, then we have [10] 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) After substitutions 𝛼 𝑘 can be approximated as a linear combination of POD modes when the ̇ 𝛼 = 𝐴 𝛼 + 𝐵 𝑢 + 𝑁 ( 𝛼 ) + 𝐹 , 𝛼 ( 0 ) = 𝛼 𝑜 , ( 2 3 ) 's are solved in the following system model. Then, the temporal model for the system is given by [10] 𝛼 ∈ ℝ 𝑛 where 𝑛 × 𝑛 , and the matrices A is 𝑛 × 2 , B is 𝑛 × 1 . N and F are both vectors 𝑦 ( 𝑡 ) = 𝛼 ( 𝑡 ) . ( 2 4 ) The output equation will be simply chosen to be 𝛼 In this model, the dimension of the state vector 𝑢 1 ( 𝑡 ) = s i n 3 𝜋 𝑡 4 , 𝑢 2 ( 𝑡 ) = s i n 3 𝜋 𝑡 2 . ( 2 5 ) is 40 which corresponds to 40 POD modes. The first 8 POD modes corresponding to the first 8 temporal coefficients are shown in Figure 3. The first model corresponds to the baseline mode and the remaining modes to actuated modes, that is, the modes due only to the input.

fig3
Figure 3: First eight POD modes.

To assess the validity of the POD model the following test inputs, which are different from the inputs used to generate the snapshots, are applied at the boundary ̇ 𝛼 𝑟 = 𝐴 𝑟 𝛼 𝑟 + 𝐵 𝑟 𝑢 + 𝑁 𝑟 𝛼 𝑟 + 𝐹 𝑟 , ( 2 6 ) In Figure 4, dashed lines denote the linear combination of POD modes restricted to the boundary. Solid lines denote the boundary test inputs. As can be seen in Figure 4, there is very good agreement between the boundary conditions specified for the full-order system and the linear combination of POD modes restricted to the boundary.

fig4
Figure 4: Boundary control accuracy.

The goal of model reduction is to construct another nonlinear system [10, 11] 𝛼 𝑟 ∈ ℝ 𝑟 where 𝑟 < 𝑛 , and 𝛼 = 𝑇 𝛼 𝑟 , 𝛼 𝑟 = 𝑇 𝛼 , 𝑇 𝑇 = 𝐼 𝑟 , ( 2 7 ) such that the behavior of the two systems is similar for states in some region of the state space. The reduced model is derived via the construction of an immersion/projection pair [21] 𝐼 𝑟 where 𝑟 × 𝑟 is the ̇ 𝛼 𝑟 = 𝑇 𝐴 𝑇 𝛼 𝑟 + 𝑇 𝐵 𝑢 + 𝑇 𝑁 𝑇 𝛼 𝑟 + 𝑇 𝐹 , 𝑦 ( 𝑡 ) = 𝐶 𝑇 𝛼 𝑟 ( 𝑡 ) . ( 2 8 ) identity matrix, resulting in the following reduced model: 𝑦 ( 𝑘 ) = ∞ ℓ = 0 𝑌 ( ℓ ) 𝑢 ( 𝑘 − ℓ ) , ( 2 9 ) 𝑌 ( 0 ) = 𝐷 , 𝑌 ( 1 ) = 𝐶 𝐵 , … , 𝑌 ( 𝑘 ) = 𝐶 𝐴 𝑘 − 1 𝐵 . ( 3 0 ) This is carried out by developing an empirical balanced truncation algorithm which is based on simulation input-output measurements of the nonlinear Galerkin model. To do so we need first to introduce balanced truncation model reduction for linear time-invariant systems.

3. Eigensystem Realization Algorithm (ERA)

Several frequency domain identification techniques are used in practice to identify the model parameters. One such method is the eigensystem realization algorithm (ERA) technique [26]. The ERA-based system realization model is created directly from empirical data and frequency domain characteristics of transfer functions. This method is applied to discrete time versions of system models.

A basic relationship between the Markov parameters and the input and output relationship in discrete time is 𝑅 𝑦 𝑢 ( 𝑘 ) = ∞ ℓ = 0 𝑌 ( ℓ ) 𝑅 𝑢 𝑢 ( 𝑘 − ℓ ) , ( 3 1 ) 𝑅 𝑢 𝑢 1 ( 𝑘 ) = 𝑚 𝑚 − 1 𝜏 = 0 𝑢 ( 𝜏 ) 𝑢 𝑇 𝑅 ( 𝑘 − 𝜏 ) , 𝑦 𝑢 1 ( 𝑘 ) = 𝑚 𝑚 − 1 𝜏 = 0 𝑦 ( 𝜏 ) 𝑢 𝑇 ( 𝑘 − 𝜏 ) . ( 3 2 ) An alternative form to (29) can be created not using the actual outputs and inputs but replacing the output term by the cross-correlation between the inputs and the corresponding outputs [26] 𝑃 𝑦 𝑢 1 ( 𝑘 ) = 𝑚 𝑚 − 1 𝜏 = 0 𝑅 𝑦 𝑢 ( 𝜏 ) 𝑒 − 𝑗 ( 2 𝜋 𝑘 / 𝑚 ) 𝜏 , 𝑃 𝑢 𝑢 1 ( 𝑘 ) = 𝑚 𝑚 − 1 𝜏 = 0 𝑅 𝑢 𝑢 ( 𝜏 ) 𝑒 − 𝑗 ( 2 𝜋 𝑘 / 𝑚 ) 𝜏 . ( 3 3 ) where the length of the data sequence is 𝐺 ( 𝑧 𝑘 𝑃 ) = 𝑦 𝑢 ( 𝑘 ) 𝑃 𝑢 𝑢 ( 𝑘 ) . ( 3 4 ) The basic process for finding the Markov parameters starts using the ratio of the power spectral density of the cross-correlation between the inputs and outputs and the power spectral density of the autocorrelation between the input signals. These power spectral densities are given by the following: 𝑌 𝑘 ∶ = 𝑌 ( 𝑘 ) = ∞ ℓ = 0 𝐺 𝑧 𝑘 𝑒 𝑗 2 𝜋 ℓ / 𝑚 𝑘 . ( 3 5 ) The ratio of the two power spectral densities is the frequency response function and is given by [26] 𝐻 𝑙 𝑞 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑌 1 𝑌 2 ⋯ 𝑌 𝑞 𝑌 2 𝑌 3 ⋯ 𝑌 𝑞 + 1 𝑌 ⋮ ⋮ ⋮ ⋮ 𝑙 𝑌 𝑙 + 1 ⋯ 𝑌 𝑞 + 𝑙 − 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 3 6 ) Then, the final step is to take the inverse Fourier transform to find the pulse response (Markov parameter) matrices [26] 𝑌 𝑘 The Hankel matrix containing the Markov parameters is of the following form: 𝑌 0 = 𝐷 , 𝑌 1 = 𝐶 𝐵 , … , 𝑌 𝑘 = 𝐶 𝐴 𝑘 − 1 𝐵 . ( 3 7 ) The individual 𝑢 1 ( 𝑡 ) = − s i n 0 . 5 5 𝑡 2 , 𝑢 2 ( 𝑡 ) = − s i n 0 . 6 0 𝑡 2 . ( 3 8 ) 's correspond to the following sequence: 𝐺 ∶ ̇ 𝑥 ( 𝑡 ) = 𝐴 𝑥 ( 𝑡 ) + 𝐵 𝑢 ( 𝑡 ) , 𝑦 ( 𝑡 ) = 𝐶 𝑥 ( 𝑡 ) + 𝐷 𝑢 ( 𝑡 ) , ( 3 9 ) In some cases, the input data for the ERA method might be provided by an experiment on a real system. However, in this paper a unique approach of using the Galerkin model in the place of the real system was used to generate the empirical data. The full-order system model was created using finite-difference methods. Recall that the control inputs were explicitly placed in the boundary conditions because the control inputs do not show up explicitly in two-dimensional Burgers' equation. However, the weak Galerkin model results in a nonlinear state space model that simplifies the relationship between the input and outputs. The chirp signals used for the excitation of the Galerkin model are of the following form and are shown in Figure 5: 𝑥 ( 𝑡 ) ∈ ℝ 𝑛

154956.fig.005
Figure 5: Excitation inputs for ERA method.

4. Empirical Balanced Truncation

The dynamics of (finite dimensional) linear time-invariant (LTI) systems are governed by a state space model of the form 𝑢 ( 𝑡 ) ∈ ℝ 𝑚 where 𝑦 ( 𝑡 ) ∈ ℝ 𝑝 is the state vector, 𝑀 is the input, and 𝑊 𝑐 is the output.

The first step in applying balanced truncation is to compute a coordinate transformation 𝑊 𝑜 , such that the controllability and observability gramians, denoted 𝐴 𝑊 𝑐 + 𝑊 𝑐 𝐴 ∗ + 𝐵 𝐵 ∗ 𝐴 = 0 , ∗ 𝑊 𝑜 + 𝑊 𝑜 𝐴 + 𝐶 ∗ 𝐶 = 0 , ( 4 0 ) and 𝐴 ∗ respectively, of the system are equal and diagonal. Assuming that (39) is stable, these gramians are solutions of the following Lyapunov equations [16, 17, 19]: 𝑊 𝑜 𝑊 = Σ = 𝑐 , ( 4 1 ) where Σ denotes complex conjugate transpose.

A balanced realization needs a similarity transformation M such that the transformed gramians are equal and satisfy [16, 17, 19] 𝑛 × 𝑛 where the matrix 𝑛 ( 𝑛 + 1 ) / 2 is a diagonal matrix containing constants in monotonically decreasing order.

Balanced model reduction requires the knowledge of the controllability and observability gramians. The latter are obtained by solving Lyapunov equations, which is prohibitive for large-scale systems. For a system with n states, the controllability and observability matrices are 𝑙 symmetric matrices and therefore solving each one of them involves finding 𝑂 𝑙 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐶 ∶ = 𝐶 𝐴 𝐶 𝐴 2 ⋮ 𝐶 𝐴 𝑙 − 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , 𝑅 𝑞 ∶ = 𝐵 𝐴 𝐵 𝐴 2 𝐵 ⋯ 𝐴 𝑞 − 1 𝐵 ( 4 2 ) unknowns. An alternative is to develop a balanced truncation algorithm based on empirical gramians, which are constructed solely from a single simulation using a sufficiently rich input. In the next paragraph we discuss how this is carried out.

The computation cost to solve large Lyapunov equations for the controllability and observability gramians prompts us to propose a balanced truncation algorithm, based on empirical gramians constructed from input-output data measurements. To this end, let us first introduce the 𝑙 -step observability and q-step controllability matrices [15] 𝑊 𝑜 𝑙 ∶ = 𝑂 ∗ 𝑙 𝑂 𝑙 , 𝑊 𝑐 𝑞 ∶ = 𝑅 𝑞 𝑅 ∗ 𝑞 . ( 4 3 ) which give rise to the 𝑙 -step observability and q-step controllability gramians l i m 𝑙 → ∞ 𝑊 𝑜 𝑙 = 𝑊 𝑜 , l i m 𝑞 → ∞ 𝑊 𝑐 𝑞 = 𝑊 𝑐 . ( 4 4 ) As the numbers q and 𝑊 𝑐 𝑞 ∶ = 𝑀 𝑊 𝑐 𝑞 𝑀 ∗ = 𝑀 ∗ − 1 𝑊 𝑜 𝑙 𝑀 − 1 𝑊 = ∶ 𝑜 𝑙 = Σ . ( 4 5 ) approach infinity, these empirical gramians approach the true gramians [15] 𝐻 𝑙 𝑞 , The goal is to find a balancing transformation matrix M that will approximately balance the empirical gramians, that is, 𝐶 𝐴 𝑘 𝐵 , 𝑘 = 0 , 1 , … , The matrix M can then be applied back to the original system model to produce an approximately balanced realization.

The product of the l-step controllability and the q-step observability matrices gives a Hankel matrix, simply a matrix that has the ith column identical to the ith row, denoted 𝐻 𝑙 𝑞 ∶ = 𝑂 𝑙 𝑅 𝑞 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐶 𝐵 𝐶 𝐴 𝐵 ⋯ 𝐶 𝐴 𝑞 − 1 𝐵 𝐶 𝐴 𝐵 𝐶 𝐴 2 𝐵 ⋯ 𝐶 𝐴 𝑞 𝐵 ⋮ ⋮ ⋱ ⋮ 𝐶 𝐴 𝑙 − 1 𝐵 𝐶 𝐴 𝑙 𝐵 ⋯ 𝐶 𝐴 𝑙 + 𝑞 − 2 𝐵 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ( 4 6 ) containing the Markov parameters 𝐻 r a n k 𝑙 𝑞 𝐻 = r a n k ( 𝑙 + 1 ) ( 𝑞 + 𝑗 ) = 𝑛 , ∀ 𝑗 ≥ 1 . ( 4 7 ) of the system in the following way: 𝐻 𝑙 𝑞 for integers l and q chosen such that [15] 𝐻 𝑙 𝑞 = 𝑈 Σ 𝑉 ∗ = 𝑈 1 𝑈 2 Σ 1 0 𝑉 0 0 ∗ 1 𝑉 ∗ 2 . ( 4 8 ) In terms of the SVD decomposition of 𝑀 = 𝑅 𝑞 𝑉 1 Σ 1 − 1 / 2 . ( 4 9 ) 𝑊 𝑐 𝑞 ∶ = 𝑀 − 1 𝑊 𝑐 𝑞 𝑀 ∗ − 1 = 𝑀 ∗ 𝑊 𝑜 𝑙 𝑀 𝑊 = ∶ 𝑜 𝑙 = Σ 1 . ( 5 0 ) The balancing transformation M is constructed as [4] 𝜎 𝑟 ≫ 𝜎 𝑟 + 1 A straightforward computation shows Σ 1 Balanced truncation can be realized by the usual way; if Σ 1 = Σ 𝑟 0 0 Σ 𝑟 + 1 , ( 5 1 ) for some r, then we can partition Σ 𝑟 𝜎 = d i a g 1 , 𝜎 2 , … , 𝜎 𝑟 , Σ 𝑟 + 1 𝜎 = d i a g 𝑟 + 1 , 𝜎 𝑟 + 2 , … , 𝜎 𝑛 . ( 5 2 ) as 𝑈 1 where 𝑉 1 A columwise conformal partition of 𝑈 1 = 𝑈 𝑟 𝑈 𝑛 − 𝑟 , 𝑉 1 = 𝑉 𝑟 𝑉 𝑛 − 𝑟 ( 5 3 ) and 𝑇 𝑟 = 𝑅 𝑞 𝑉 𝑟 Σ 𝑟 − 1 / 2 , 𝑇 𝑟 = Σ 𝑟 − 1 / 2 𝑈 ∗ 𝑟 𝑂 𝑙 , 𝑇 𝑟 𝑇 𝑟 = 𝐼 𝑟 ( 5 4 ) 𝐴 𝑟 = 𝑇 𝑟 𝐴 𝑇 𝑟 , 𝐵 𝑟 = 𝑇 𝑟 𝐵 , 𝐶 𝑟 𝑇 = 𝐶 𝑟 . ( 5 5 ) yields the immersion/projection pair [4] 𝐶 𝐴 𝑘 𝐵 , 𝑘 = 0 , 1 , … , 𝑙 + 𝑞 − 1 . ( 5 6 ) and from which a reduced-order r-dimensional model with state matrices is deduced ̇ 𝛼 = 𝐴 𝛼 + 𝐵 𝑢 + 𝑁 ( 𝛼 ) + 𝐹 , 𝛼 ( 0 ) = 𝛼 𝑜 ( 5 7 ) The above construction only requires estimates of the Markov parameters 𝛼 s s . The Markov parameters can be computed from a single simulation in which a sufficiently rich input signal is applied and the output responses are collected. In the next section, the discrete Fourier transform (DFT) is used to map time domain data into spectral densities from which frequency response estimates are calculated using the ERA [26].

5. Application to the Galerkin Model

The empirical balanced truncation based on linear systems is applied to the Galerkin model 𝛼 = 𝑇 𝛼 𝑟 , 𝛼 𝑟 = 𝑇 𝛼 . ( 5 8 ) which has an equilibrium in steady state, denoted by ̇ 𝛼 𝑟 = 𝐴 𝑟 𝛼 𝑟 + 𝐵 𝑟 𝑢 + 𝑁 𝑟 𝑇 𝛼 𝑟 + 𝐹 𝑟 , 𝛼 𝑟 ( 0 ) = 𝑇 𝛼 𝑜 , 𝐴 𝑟 ∶ = 𝑇 𝐴 𝑟 𝑇 , 𝐵 𝑟 ∶ = 𝑇 𝐵 , 𝑁 𝑟 ∶ = 𝑇 𝑁 , 𝐹 𝑟 ∶ = 𝑇 𝐹 . ( 5 9 ) The rationale for doing so is that linear subspace approximations to exact submanifolds associated with nonlinear controllability and observability require only standard matrix manipulations utilizing simulation/experimental data as explained in [4, 20, 21]. The computational advantages of the scheme presented here carry over directly to the nonlinear setting.

The reduced-order model is derived as discussed through the construction of the immersion/projection nonlinear system pair 𝛼 s s , This results in the following reduced-order model: ̇ 𝛼 = 𝐴 ℓ 𝛼 + 𝐵 ℓ 𝑢 , 𝛼 ( 0 ) = 𝛼 𝑜 , ( 6 0 ) If (57) has a linearization around the steady state equilibrium 𝐴 ℓ = 𝜕 𝐴 𝛼 + 𝐵 𝑢 ( 𝑡 ) + 𝑁 ( 𝛼 ) + 𝐹 | | | 𝜕 𝛼 𝛼 = 𝛼 𝑠 𝑠 , 𝐵 ℓ = 𝜕 𝐴 𝛼 + 𝐵 𝑢 ( 𝑡 ) + 𝑁 ( 𝛼 ) + 𝐹 | | | 𝜕 𝑢 𝑢 = 𝐵 , ( 6 1 ) 𝛼 s s where ̇ 𝛼 𝑟 = 𝐴 ℓ 𝑟 𝛼 𝑟 + 𝐵 ℓ 𝑟 𝑢 , 𝛼 𝑟 ( 0 ) = 𝑇 𝛼 𝑜 , ( 6 2 ) and the reduced system linearization around the steady state equilibrium 𝐴 ℓ 𝑟 = 𝜕 𝐴 𝑟 𝛼 𝑟 + 𝐵 𝑟 𝑢 ( 𝑡 ) + 𝑁 𝑟 𝑇 𝛼 𝑟 + 𝐹 𝑟 | 𝜕 𝛼 𝛼 = 𝛼 s s , 𝐵 ℓ 𝑟 = 𝜕 𝐴 𝑟 𝛼 𝑟 + 𝐵 𝑟 𝑢 ( 𝑡 ) + 𝑁 𝑟 𝑇 𝛼 𝑟 ) + 𝐹 𝑟 | 𝜕 𝑢 𝑢 = 𝐵 𝑟 . ( 6 3 ) 𝛼 s s where 𝐴 ℓ 𝑟 = 𝑇 𝐴 ℓ 𝑇 , 𝐵 ℓ 𝑟 = 𝑇 𝐵 ℓ . ( 6 4 ) The linearization of both models about equilibrium 𝑤 ( 𝑡 , 𝑥 , 𝑦 ) is related by 𝑤 𝑟 ( 𝑡 , 𝑥 , 𝑦 ) . Empirical balanced truncation applied to the 40th-order Galerkin model resulted in 14th-order reduced models. The first 8 temporal coefficients of the 14th-order reduced model and 2000th full-order model are plotted in Figure 6. Figure 6 shows good agreement between the temporal coefficients.

154956.fig.006
Figure 6: Projected and POD model coefficients.

In Figure 7, we compare the Hankel singular values of the 2000th full-order linearized and reduced 14th-order empirical models. As expected the Hankel singular values corresponding to the reduced-order model are smaller than the full-order model; nevertheless, the figure shows that they are close.

154956.fig.007
Figure 7: Comparison of Hankel singular values.

In Figure 8, we compare the full-order solution 𝐻 ∞ of Burgers’ equation with the solution based on the 14th-order ERA model 𝑤 r e f ( 𝑥 ) The figure shows that they behave similarly especially at the boundary where control is applied.

fig8
Figure 8: Full- and reduced-order models' responses.

6. Controller Design

An 𝐻 ∞ controller was designed based on the linearized 14th reduced model and applied to the full-order model using Matlab. The performance was to achieve tracking a fixed reference signal 𝛼 r e f specified for the full-order model. The tracking problem is depicted in Figure 9, where C is the controller and P the plant. The computation of the 𝑦 ( 𝑡 ) , is based on the 14th-order reduced model.

154956.fig.009
Figure 9: Closed-loop control system for tracking.

From Figure 9, for tracking purposes the controlled output is chosen to be the error signal e which is defined to be the difference between the reference 𝑒 ( 𝑡 ) ∶ = 𝛼 ( 𝑡 ) − 𝛼 r e f . ( 6 5 ) and the actual output 𝐻 ∞ that is, 𝛼 r e f The objective of the 𝑇 𝑒 𝛼 r e f controller C is to stabilize the closed-loop system and minimize the error e. From Figure 9, in terms of transfer function matrices of P and C, the transfer matrix from 𝑇 𝑒 𝛼 r e f ∶ = ( 𝐼 + 𝑃 𝐶 ) − 1 . ( 6 6 ) to e is given by the sensitivity function 𝐻 ∞ defined by ‖ ‖ 𝑇 𝑒 𝛼 r e f ‖ ‖ ∞ ‖ ‖ = s u p 𝛼 r e f ‖ ‖ 2 ≤ 1 ‖ ‖ 𝑒 ‖ ‖ 2 . ( 6 7 ) The 𝑤 r e f ( 𝑥 ) control design reduces to the following optimization. Find C such that the closed-loop system is stable and minimizes 𝐻 ∞ Projecting 𝐻 ∞ onto the POD basis yields the tracking coefficients for the reduced-order model. After computing the controller we close the loop on the original 2000th full-order model. The controller is only 14th order since based on the 14th-order reduced model. The projected reference onto the POD basis initial condition and reference are shown in Figure 10. The controlled flow with the action of the boundary controller is shown in Figure 11. The figure shows good tracking performance.

fig10
Figure 10: Flow initial condition and reference.
154956.fig.011
Figure 11: Controlled flow.

7. Conclusion

Empirical balanced truncation has been considered in conjunction with POD as an approach for deriving reduced-order models and applied to 2D Burgers’ equation. Like POD, empirical balanced truncation is based on simulation/experimental data and can be implemented via standard matrix computations. Improvements to the scheme originally proposed in [4, 20] have been presented that lead to reduced data requirements that may become significant for applications such as aerodynamic flow control. Essentially, the balancing transformation is constructed via Markov parameters that can be identified from measurements collected in a single experiment/simulation. The approach has been applied with favorable results to 2D Burgers’ equation, a partial differential equation in two spatial dimensions that possesses features comparable to the Navier-Stokes equations governing fluid flow. An feedback flow controller was designed based on the empirical reduced model to achieve flow tracking. The closed loop on the full-order model shows good flow tracking performance.

Acknowledgments

This work was supported in part by an Air Force Summer Faculty Fellowship Program, Air Force Office of Scientific Research under Contract AF-FA9550-08-1-0450, and National Science Foundation under Contract NSF-CMMI-0825921.

References

  1. K. Willcox and J. Peraire, “Balanced model reduction via the proper orthogonal decomposition,” AIAA Journal, vol. 40, no. 11, pp. 2323–2330, 2002.
  2. S. M. Djouadi, R. C. Camphouse, and J. H. Myatt, “Optimal order reduction for the the two-dimensional Burgers' equation,” in Proceedings of the 46th IEEE Conference on Decision and Control (CDC '07), pp. 3507–3512, New Orleans, La, USA, December 2007.
  3. M. Ilak and C. W. Rowley, “Reduced-order modeling of channel flow using traveling POD and balanced POD,” in Proceedings of the 3rd AIAA Flow Control Conference, vol. 2, pp. 1–11, San Francisco, Calif, USA, June 2006.
  4. D. A. Lawrence, “Empirical model reduction for active closed-loop flow control,” Final Report, Faculty Summer Research Opportunity, ARL, State College, Pa, USA, 2004.
  5. F. Leibfritx and S. Volkwein, “Numerical feedback controller design for PDE systems using model reduction: techniques and case studies,” in Real-Time PDE-Constrained Optimization, SIAM, Philadelphia, Pa, USA, 2006.
  6. C. W. Rowley, T. Colonius, and R. M. Murray, “Model reduction for compressible flows using POD and Galerkin projection,” Physica D, vol. 189, no. 1-2, pp. 115–129, 2004.
  7. M. N. Glauser, M. J. Young, H. Higuchi, C. E. Tinney, and H. Carlson, “POD based experimental flow control on a NACA-4412 airfoil,” in Proceedings of the 42nd AIAA Aerospace Sciences Meeting and Exhibit, pp. 4861–4870, Reno, Nev, USA, January 2004, paper no. 2004-0575.
  8. H. A. Carlson, M. Glauser, and R. Roveda, “Models for controlling airfoil lift and drag,” in Proceedings of the 42nd AIAA Aerospace Sciences Meeting and Exhibit, pp. 5627–5638, Reno, Nev, USA, January 2004, paper no. 2004-0579.
  9. K. Cohen, S. Siegel, T. McLaughlin, and J. Myatt, “Proper orthogonal decomposition modeling of a controlled Ginzburg-Landau cylinder wake model,” in Proceedings of the 21st International Communications Satellite Systems Conference and Exhibit, Reno, Nev, USA, January 2003, paper no. 2003-2405.
  10. R. C. Camphouse and J. H. Myactt, “Reduced order modelling and boundary feedback control of nonlinear convection,” in Proceedings of the AIAA Guidance, Navigation, and Control Conference, vol. 1, pp. 289–297, San Francisco, Calif, USA, August 2005, paper no. 2005-5844.
  11. R. C. Camphouse, S. M. Djouadi, and J. H. Myatt, “Feedback control for aerodynamics,” in Proceedings of the 45th IEEE Conference on Decision and Control (CDC '06), pp. 1–14, San-Diego, Calif, USA, December 2006.
  12. W. R. Graham, J. Peraire, and K. Y. Tang, “Optimal control of vortex shedding using low-order models—part I: open-loop model development,” International Journal for Numerical Methods in Engineering, vol. 44, no. 7, pp. 945–972, 1999.
  13. H. Banks, R. Del Rosario, and R. Smith, “Reduced order model feedback control design: numerical implantation in a thin shell model,” Center for Research in Scientific Computation, North Carolina State University, Raleigh, NC, USA, June 1998.
  14. E. Caraballo, M. Samimy, and J. DeBonis, “Low dimensional modeling of flow for closed-loop flow control,” in Proceedings of the AIAA Aerospace Sciences Meeting, Reno, Nev, USA, January 2003, paper no. 2003-0059.
  15. W. J. Rugh, Linear System Theory, Prentice-Hall, Upper Saddle River, NJ, USA, 2nd edition, 1996.
  16. K. Zhou, J. C. Doyle, and K. Glover, Robust and Optimal Control, Prentice-Hall, Upper Saddle River, NJ, USA, 1996.
  17. K. Zhou and J. C. Doyle, Essentials of Robust Control, Prentice-Hall, Upper Saddle River, NJ, USA, 1998.
  18. A. Defant and K. Floret, Tensor Norms and Operator Ideals, Elsevier Science, Amsterdam, The Netherlands, 1993.
  19. B. C. Moore, “Principal component analysis in linear systems: controllability, observability, and model reduction,” IEEE Transactions on Automatic Control, vol. 26, no. 1, pp. 17–32, 1981.
  20. J. Hahn and T. F. Edgar, “Reduction of nonlinear models using balancing of empirical gramians and Galerkin projections,” in Proceedings of the American Control Conference, vol. 4, pp. 2864–2868, Chicago, Ill, USA, June 2000.
  21. M. Condon and R. Ivanov, “Empirical balanced truncation of nonlinear systems,” Journal of Nonlinear Science, vol. 14, no. 5, pp. 405–414, 2004.
  22. C. A. Antoulas, Approximation of Large-Scale Dynamical Systems, Advances in Design and Control, SIAM, Philadelphia, Pa, USA, 2006.
  23. V. Sreeram, “Frequency response error bounds for time-weighted balanced truncation,” in Proceedings of the 41st IEEE Conference on Decision and Control, vol. 3, pp. 3330–3331, Las Vegas, Nev, USA, December 2002.
  24. T. J. Chung, Computational Fluid Dynamics, Cambridge University Press, Cambridge, UK, 2002.
  25. P. Holmes, J. L. Lumley, and G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, Cambridge, UK, 1996.
  26. J. Juang, Applied System Identification, Prentice-Hall, Upper Saddle River, NJ, USA, 1994.