Empirical Reduced-Order Modeling for Boundary Feedback Flow Control

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 10¹² Riccati unknowns need to be calculated for a discretized flow model describing 10⁶ 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 Ω_gap be the region defined by [𝑎₁,𝑎₂]×[𝑏₁,𝑏₂]. Let Ω_full be the region defined by (𝑎₀,𝑎_end)×(𝑏₀,𝑏_end). Then the problem domain is given by Ω=Ω_full/Ω_gap. In this problem setup, Ω_gap is an obstacle.

Figure 1 

Problem geometry.

The system dynamics that act within the problem domain are described by the two-dimensional (2D) Burgers’ equation [10] 𝜕=1𝜕𝑡𝑤(𝑡,𝑥,𝑦)+∇⋅𝐹(𝑤)𝑟𝜕²𝜕𝑥²𝜕𝑤(𝑡,𝑥,𝑦)+²𝜕𝑦²𝑤(𝑡,𝑥,𝑦), (1) where the form of 𝐹(𝑤) is 𝐹(𝑤)=𝑐₁𝑤²(𝑡,𝑥,𝑦)2𝑐₂𝑤²(𝑡,𝑥,𝑦)2^𝑇. (2) In this case, the value for 𝑐₁ is equal to 1 and 𝑐₂ 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 Γ_topandΓ_bottom. 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: 𝑤𝑡,Γ_bottom=𝑢_bottom(𝑡)Ψ_bottom𝑤(𝑥),𝑡,Γ_top=𝑢_top(𝑡)Ψ_top(𝑥). (4) Here 𝑢_top(𝑡) and 𝑢_bottom(𝑡) are the control inputs on the top and bottom boundaries, respectively; the spatial functions Ψ_top(𝑥) and Ψ_bottom(𝑥) describe the spatial effect that the controls have on the top and bottom boundaries.

The boundary condition on the airflow intake side is 𝑤𝑡,Γ_in=𝑓(𝑦) (5) and it is parabolic in nature. The airflow outtake side has a Neumann boundary condition that has the form [24] 𝜕𝑤𝜕𝑥𝑡,Γ_out=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,𝑥,𝑦)=𝑤₀(𝑥,𝑦)∈𝐿²(Ω). (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] 𝑢_bottom(𝑡)=𝛽sin0.25𝑡²,𝑢_top𝑢(𝑡)=0,_bottom(𝑡)=0,𝑢_top(𝑡)=𝛽sin0.25𝑡²,𝑢_bottom(𝑡)=𝛽sin0.25𝑡²,𝑢_top(𝑡)=𝛽sin0.25𝑡², (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.

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 𝐿²([0,𝑇)×Ω). The solution can be approximated as 𝑤(𝑡,𝑥,𝑦)≈^𝑛_𝑘=1𝛼_𝑘(𝑡)𝜙_𝑘(𝑥,𝑦), (9) where the 𝛼_𝑘's are time varying coefficients that multiply the spatial functions 𝜙_𝑘's and 𝛼_𝑘∈𝐿²[0,𝑇),𝜙_𝑘(𝑥,𝑦)∈𝐿²(Ω), (10) where 𝐿²([0,𝑇)) and 𝐿²(Ω) 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 𝐿²([0,𝑇))⊗𝐿²(Ω)∶=^𝑛_𝑘=1𝛼_𝑘(𝑡)𝜙_𝑘(𝑥,𝑦),𝛼_𝑘(𝑡)∈𝐿²𝜙([0,𝑇)),_𝑘(𝑥,𝑦)∈𝐿²(Ω),∀𝑛integer (11) is dense in 𝐿²([0,𝑇)×Ω) [18].

Any basis for 𝐿²(Ω) 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]: 𝜇_𝑛∶=inf_{𝛼𝑘∈𝐿2𝜙([0,𝑇)),𝑘(𝑥,𝑦)∈𝐿2(Ω)}𝑤(𝑡,𝑥,𝑦)−^𝑛_𝑘=1𝛼_𝑘(𝑡)𝜙_𝑘(𝑥,𝑦)₂.𝜇_𝑛⟶_𝑛→∞0. (12) 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] max_{𝜑𝑖∈𝐿2(Ω)}𝑢_𝑖,𝜑_𝑖²𝜑_𝑖²₂. (13) To construct numerically the POD basis {𝜑_𝑗}, we build the correlation matrix L, of size 𝑀×𝑀, composed of the inner products of the snapshots. In 𝐿²(Ω) we have 𝐿_𝑘,𝑗=𝑆_𝑘,𝑆_𝑗,𝑆_𝑘,𝑆_𝑗=_Ω𝑆_𝑘𝑆^∗_𝑗𝑑𝑥𝑑𝑦, (14)(15) where ∗ denotes complex conjugate transpose.

A spectral decomposition of the matrix L is performed. The n largest eigenvalues {𝜆₁,𝜆₂,…,𝜆_𝑛} 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 {𝜑_𝑗} 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,for𝑘={1,2,…,𝑛}. (16) Then, the POD basis set is formed according to 𝜑_𝑘(𝑥,𝑦)=^𝑀_𝑗=1𝑣_𝑘,𝑗𝑆_𝑗(𝑥,𝑦) (17) 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 {𝛼_𝑘} [10].

Projecting (1) onto the POD basis yields [10] _Ω𝜕𝜕𝑡𝑤(𝑡,𝑥,𝑦)𝜑_𝑘=1(𝑥,𝑦)𝑑𝑥Re_𝜕Ω∇𝑤(𝑡,𝑥,𝑦)⋅𝑛𝜑_𝑘−(𝑥,𝑦)𝑑𝐴(𝑥)_Ω∇𝑤(𝑡,𝑥,𝑦)⋅∇𝜑_𝑘−(𝑥,𝑦)𝑑𝑥_𝜕Ω𝐹(𝑤)⋅𝑛𝜑_𝑘−(𝑥,𝑦)𝑑𝐴(𝑥)_Ω𝐹⋅∇𝜑_𝑘(𝑥,𝑦)𝑑𝑥, (18) where the first term on the right-hand side is ∇𝑤(𝑡,𝑥,𝑦)⋅𝑛𝜑_𝜕Ω𝑘=(𝑥,𝑦)𝑑𝐴(𝑥)_𝑎2𝑎1𝜕𝑤𝜕𝑦𝑡,𝑥,𝑏₁)𝜑_𝑘𝑥,𝑏₁𝑑𝑥−^𝑎2_𝑎1𝜕𝑤𝜕𝑦𝑡,𝑥,𝑏₂𝜑_𝑘𝑥,𝑏₂−𝑑𝑥_𝑏end𝑏0𝜕𝑤𝜕𝑥𝑡,𝑎₀,𝑦𝜑_𝑘𝑎₀,𝑦𝑑𝑦. (19) The Neumann boundary condition forces the portion of the boundary integral over 𝑏₀ to 𝑏_end along 𝑎_end to be 0, that is, ^𝑏end_𝑏0𝜕𝑤𝜕𝑥𝑡,𝑎_end,𝑦𝜑_𝑘𝑎_end,𝑦𝑑𝑦=0. (20) The second boundary integral is decomposed as follows 𝐹(𝑤)⋅𝑛𝜑_𝜕Ω𝑘=1(𝑥,𝑦)𝑑𝐴(𝑥)2^𝑏end_𝑏0𝑤𝑡,𝑎_end,𝑦²𝜑_𝑘𝑎_end,𝑦−𝑓(𝑦)²𝜑_𝑘𝑎₀,𝑦𝑑𝑦. (21) 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] 𝜕𝑤𝜕𝑦𝑡,𝑥,𝑏₁≈𝑢_bottom(𝑡)Ψ_bottom(𝑥)−𝑤𝑡,𝑥,𝑏₁−ℎℎ,𝜕𝑤𝜕𝑦𝑡,𝑥,𝑏₂≈𝑤𝑡,𝑥,𝑏₂+ℎ−𝑢_top(𝑡)Ψ_top(𝑥)ℎ,𝜕𝑤𝜕𝑥𝑡,𝑎₀≈𝑤,𝑦𝑡,𝑎₀+ℎ,𝑦−𝑓(ℎ)ℎ. (22) After substitutions 𝑤(𝑡,𝑥,𝑦) can be approximated as a linear combination of POD modes when the 𝛼_𝑘's are solved in the following system model. Then, the temporal model for the system is given by [10] 𝛼̇=𝐴𝛼+𝐵𝑢+𝑁(𝛼)+𝐹,𝛼(0)=𝛼_𝑜, (23) where 𝛼∈ℝ^𝑛 and the matrices A is 𝑛×𝑛,B is 𝑛×2,N and F are both vectors 𝑛×1. The output equation will be simply chosen to be 𝑦(𝑡)=𝛼(𝑡). (24) In this model, the dimension of the state vector 𝛼 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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

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 𝑢₁(𝑡)=sin3𝜋𝑡4,𝑢₂(𝑡)=sin3𝜋𝑡2. (25) 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.

(a)

(b)

Figure 4 

Boundary control accuracy.

The goal of model reduction is to construct another nonlinear system [10, 11] 𝛼̇_𝑟=𝐴_𝑟𝛼_𝑟+𝐵_𝑟𝑢+𝑁_𝑟𝛼_𝑟+𝐹_𝑟, (26) where 𝛼_𝑟∈ℝ^𝑟 and 𝑟<𝑛, 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] 𝛼=𝑇𝛼_𝑟,𝛼_𝑟=𝑇𝛼,𝑇𝑇=𝐼_𝑟, (27) where 𝐼_𝑟 is the 𝑟×𝑟 identity matrix, resulting in the following reduced model: 𝛼̇_𝑟=𝑇𝐴𝑇𝛼_𝑟+𝑇𝐵𝑢+𝑇𝑁𝑇𝛼_𝑟𝑇𝛼+𝑇𝐹,𝑦(𝑡)=𝐶_𝑟(𝑡). (28) 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𝑌(ℓ)𝑢(𝑘−ℓ),𝑌(0)=𝐷,𝑌(1)=𝐶𝐵,…,𝑌(𝑘)=𝐶𝐴^𝑘−1𝐵. (29)(30) 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] 𝑅_𝑦𝑢(𝑘)=^∞_ℓ=0𝑌(ℓ)𝑅_𝑢𝑢(𝑘−ℓ), (31) where the length of the data sequence is 𝑅_𝑢𝑢1(𝑘)=𝑚^𝑚−1_𝜏=0𝑢(𝜏)𝑢^𝑇𝑅(𝑘−𝜏),_𝑦𝑢1(𝑘)=𝑚^𝑚−1_𝜏=0𝑦(𝜏)𝑢^𝑇(𝑘−𝜏). (32) 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: 𝑃_𝑦𝑢1(𝑘)=𝑚^𝑚−1_𝜏=0𝑅_𝑦𝑢(𝜏)𝑒^{−𝑗(2𝜋𝑘/𝑚)𝜏},𝑃_𝑢𝑢1(𝑘)=𝑚^𝑚−1_𝜏=0𝑅_𝑢𝑢(𝜏)𝑒^{−𝑗(2𝜋𝑘/𝑚)𝜏}. (33) The ratio of the two power spectral densities is the frequency response function and is given by [26] 𝐺(𝑧_𝑘𝑃)=^𝑦𝑢(𝑘)𝑃_𝑢𝑢.(𝑘) (34) Then, the final step is to take the inverse Fourier transform to find the pulse response (Markov parameter) matrices [26] 𝑌_𝑘∶=𝑌(𝑘)=^∞_ℓ=0𝐺𝑧_𝑘𝑒^{𝑗(2𝜋ℓ/𝑚)𝑘}. (35) The Hankel matrix containing the Markov parameters is of the following form: 𝐻_𝑙𝑞=𝑌₁𝑌₂⋯𝑌_𝑞𝑌₂𝑌₃⋯𝑌_𝑞+1𝑌⋮⋮⋮⋮_𝑙𝑌_𝑙+1⋯𝑌_{𝑞+𝑙−1}. (36) The individual 𝑌_𝑘's correspond to the following sequence: 𝑌₀=𝐷,𝑌₁=𝐶𝐵,…,𝑌_𝑘=𝐶𝐴^𝑘−1𝐵. (37) 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: 𝑢₁(𝑡)=−sin0.55𝑡²,𝑢₂(𝑡)=−sin0.60𝑡². (38)

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 𝐺∶𝑥̇(𝑡)=𝐴𝑥(𝑡)+𝐵𝑢(𝑡),𝑦(𝑡)=𝐶𝑥(𝑡)+𝐷𝑢(𝑡), (39) 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 𝑊_𝑐 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]: 𝐴𝑊_𝑐+𝑊_𝑐𝐴^∗+𝐵𝐵^∗𝐴=0,_∗𝑊_𝑜+𝑊_𝑜𝐴+𝐶^∗𝐶=0, (40) 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] 𝑊_𝑜=Σ=𝑊_𝑐, (41) where the matrix Σ 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 𝑛(𝑛+1)/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] 𝑂_𝑙𝐶∶=𝐶𝐴𝐶𝐴²⋮𝐶𝐴^𝑙−1,𝑅_𝑞∶=𝐵𝐴𝐵𝐴²𝐵⋯𝐴^𝑞−1𝐵 (42) which give rise to the 𝑙-step observability and q-step controllability gramians 𝑊_𝑜𝑙∶=𝑂^∗_𝑙𝑂_𝑙,𝑊_𝑐𝑞∶=𝑅_𝑞𝑅^∗_𝑞. (43) As the numbers q and 𝑙 approach infinity, these empirical gramians approach the true gramians [15] lim_𝑙→∞𝑊_𝑜𝑙=𝑊_𝑜,lim_𝑞→∞𝑊_𝑐𝑞=𝑊_𝑐. (44) The goal is to find a balancing transformation matrix M that will approximately balance the empirical gramians, that is, 𝑊_𝑐𝑞∶=𝑀𝑊_𝑐𝑞𝑀^∗=𝑀^∗−1𝑊_𝑜𝑙𝑀⁻¹=∶𝑊_𝑜𝑙=Σ. (45) 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 𝐻_𝑙𝑞, containing the Markov parameters 𝐶𝐴^𝑘𝐵,𝑘=0,1,…, of the system in the following way: 𝐻_𝑙𝑞∶=𝑂_𝑙𝑅_𝑞=𝐶𝐵𝐶𝐴𝐵⋯𝐶𝐴^𝑞−1𝐵𝐶𝐴𝐵𝐶𝐴²𝐵⋯𝐶𝐴^𝑞𝐵⋮⋮⋱⋮𝐶𝐴^𝑙−1𝐵𝐶𝐴^𝑙𝐵⋯𝐶𝐴^{𝑙+𝑞−2}𝐵 (46) for integers l and q chosen such that [15] rank𝐻_𝑙𝑞=rank𝐻_{(𝑙+1)(𝑞+𝑗)}=𝑛,∀𝑗≥1. (47) In terms of the SVD decomposition of 𝐻_𝑙𝑞𝐻_𝑙𝑞=𝑈Σ𝑉^∗=𝑈₁𝑈₂Σ₁0𝑉00𝑉^∗1∗₂. (48) The balancing transformation M is constructed as [4] 𝑀=𝑅_𝑞𝑉₁Σ₁^−1/2. (49) A straightforward computation shows 𝑊_𝑐𝑞∶=𝑀⁻¹𝑊_𝑐𝑞𝑀^∗−1=𝑀^∗𝑊_𝑜𝑙𝑀=∶𝑊_𝑜𝑙=Σ₁. (50) Balanced truncation can be realized by the usual way; if 𝜎_𝑟≫𝜎_𝑟+1 for some r, then we can partition Σ₁ as Σ₁=Σ_𝑟00Σ_𝑟+1, (51) where Σ_𝑟=diag𝜎₁,𝜎₂,…,𝜎_𝑟,Σ_𝑟+1=diag𝜎_𝑟+1,𝜎_𝑟+2,…,𝜎_𝑛. (52) A columwise conformal partition of 𝑈₁ and 𝑉₁𝑈₁=𝑈_𝑟𝑈_𝑛−𝑟,𝑉₁=𝑉_𝑟𝑉_𝑛−𝑟 (53) yields the immersion/projection pair [4] 𝑇_𝑟=𝑅_𝑞𝑉_𝑟Σ_𝑟^−1/2,𝑇_𝑟=Σ_𝑟^−1/2𝑈^∗_𝑟𝑂_𝑙,𝑇_𝑟𝑇_𝑟=𝐼_𝑟 (54) and from which a reduced-order r-dimensional model with state matrices is deduced 𝐴_𝑟=𝑇_𝑟𝐴𝑇_𝑟,𝐵_𝑟=𝑇_𝑟𝐵,𝐶_𝑟=𝐶𝑇_𝑟. (55) The above construction only requires estimates of the Markov parameters 𝐶𝐴^𝑘𝐵,𝑘=0,1,…,𝑙+𝑞−1. (56) 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 𝛼̇=𝐴𝛼+𝐵𝑢+𝑁(𝛼)+𝐹,𝛼(0)=𝛼_𝑜 (57) which has an equilibrium in steady state, denoted by 𝛼_ss. 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 𝛼=𝑇𝛼_𝑟,𝛼_𝑟=𝑇𝛼. (58) This results in the following reduced-order model: 𝛼̇_𝑟=𝐴_𝑟𝛼_𝑟+𝐵_𝑟𝑢+𝑁_𝑟𝑇𝛼_𝑟+𝐹_𝑟,𝛼_𝑟(0)=𝑇𝛼_𝑜,𝐴_𝑟∶=𝑇𝐴_𝑟𝑇,𝐵_𝑟∶=𝑇𝐵,𝑁_𝑟∶=𝑇𝑁,𝐹_𝑟∶=𝑇𝐹. (59) If (57) has a linearization around the steady state equilibrium 𝛼_ss,𝛼̇=𝐴_ℓ𝛼+𝐵_ℓ𝑢,𝛼(0)=𝛼_𝑜, (60) where 𝐴_ℓ=𝜕𝐴𝛼+𝐵𝑢(𝑡)+𝑁(𝛼)+𝐹𝜕𝛼_{𝛼=𝛼𝑠𝑠},𝐵_ℓ=𝜕𝐴𝛼+𝐵𝑢(𝑡)+𝑁(𝛼)+𝐹𝜕𝑢_𝑢=𝐵, (61) and the reduced system linearization around the steady state equilibrium 𝛼_ss𝛼̇_𝑟=𝐴_ℓ𝑟𝛼_𝑟+𝐵_ℓ𝑟𝑢,𝛼_𝑟(0)=𝑇𝛼_𝑜, (62) where 𝐴_ℓ𝑟=𝜕𝐴_𝑟𝛼_𝑟+𝐵_𝑟𝑢(𝑡)+𝑁_𝑟𝑇𝛼_𝑟+𝐹_𝑟𝜕𝛼_𝛼=𝛼ss,𝐵_ℓ𝑟=𝜕𝐴_𝑟𝛼_𝑟+𝐵_𝑟𝑢(𝑡)+𝑁_𝑟𝑇𝛼_𝑟)+𝐹_𝑟𝜕𝑢_𝑢=𝐵_𝑟. (63) The linearization of both models about equilibrium 𝛼_ss is related by 𝐴_ℓ𝑟=𝑇𝐴_ℓ𝑇,𝐵_ℓ𝑟=𝑇𝐵_ℓ. (64) 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.

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.

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 𝑤_𝑟(𝑡,𝑥,𝑦). The figure shows that they behave similarly especially at the boundary where control is applied.

(a)

(b)

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 𝑤_ref(𝑥) 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.

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 𝛼_ref and the actual output 𝑦(𝑡), that is, 𝑒(𝑡)∶=𝛼(𝑡)−𝛼_ref. (65) The objective of the 𝐻^∞ 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 𝛼_ref to e is given by the sensitivity function 𝑇_𝑒𝛼ref defined by 𝑇_𝑒𝛼ref∶=(𝐼+𝑃𝐶)⁻¹. (66) The 𝐻^∞ control design reduces to the following optimization. Find C such that the closed-loop system is stable and minimizes 𝑇_{𝑒𝛼ref∞}=sup_{‖𝛼ref‖2≤1}‖𝑒‖₂. (67) Projecting 𝑤_ref(𝑥) 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.

(a)

(b)

Figure 10 

Flow initial condition and reference.

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.