# Empirical Reduced-Order Modeling for  Boundary Feedback Flow Control

Authors: Seddik M. Djouadi; R. Chris Camphouse; James H. Myatt
Journal: Journal of Control Science and Engineering
 (2008)

---

## 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. AnH∞ feedback flow controller is designed based on the reduced model to achieve tracking and then applied to the full-order model with excellent performance.


---

## Body

## 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, roughly1012 Riccati unknowns need to be calculated for a
discretized flow model describing 106 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 Section2, 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 H∞ 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 H∞ 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 Figure1. 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 [a1,a2]×[b1,b2]. Let Ωfull be the region defined by (a0,aend)×(b0,bend). 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)∂∂tw(t,x,y)+∇⋅F(w)=1r(∂2∂x2w(t,x,y)+∂2∂y2w(t,x,y)),  where the form of F(w) is (2)F(w)=[c1w2(t,x,y)2c2w2(t,x,y)2]T.  In this case, the value for c1 is equal to 1 and c2 is equal to 0. The value r 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 f(x) on a domain boundary [24]. The form of the boundary condition is (3)w(t,x,y)=f(t,x,y)∀(x,y)∈∂Ω.  The boundary conditions on the top
and bottom are described by the following: (4)w(t,Γbottom)=ubottom(t)Ψbottom(x),w(t,Γtop)=utop(t)Ψtop(x).  Here utop(t) and ubottom(t) are the control inputs on the top and bottom
boundaries, respectively;
the spatial functions Ψtop(x) and Ψbottom(x) describe the spatial effect that the controls
have on the top and bottom boundaries.The boundary condition on the airflow
intake side is(5)w(t,Γin)=f(y)  and it is parabolic in nature. The
airflow outtake side has a Neumann boundary condition that has the form [24] (6)∂∂xw(t,Γout)=0.  On all of the remaining boundaries of Ω,w(t,x,y) is set equal to 0 for all values of t. Finally, the initial conditions for
the interior are given by  (7)w(0,x,y)=w0(x,y)∈L2(Ω). 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] (8)ubottom(t)=β⁢sin⁡(0.25t2),utop(t)=0,ubottom(t)=0,utop(t)=β⁢sin⁡(0.25t2),ubottom(t)=β⁢sin⁡(0.25t2),utop(t)=β⁢sin⁡(0.25t2),       where the values for β are −3,−2 and −1 and the range for t 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{Sk(x,y)}k=1M [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 spaceL2([0,T)×Ω).
The solution can be approximated as (9)w(t,x,y)≈∑k=1nαk(t)ϕk(x,y),  where the αk's are time varying coefficients that multiply
the spatial functions ϕk's and (10)αk∈L2([0,T)),ϕk(x,y)∈L2(Ω),  where L2([0,T)) and L2(Ω) are the standard Hilbert spaces of absolutely
square integrable functions defined, respectively, on the time interval [0,T) and spatial domain Ω.
The approximation (9) can be as accurate as desired since the tensor space (11)L2([0,T))⊗L2(Ω):={∑k=1nαk(t)ϕk(x,y),αk(t)∈L2([0,T)),ϕk(x,y)∈L2(Ω),∀n⁢integer}  is dense in L2([0,T)×Ω) [18].Any basis forL2(Ω) can be used to construct the approximation of
the solution w(t,x,y).
Here we use the POD basis {ϕk} since it is optimal in the following sense [2]: (12)μn:=inf⁡αk∈L2([0,T)),ϕk(x,y)∈L2(Ω)‖w(t,x,y)−∑k=1nαk(t)ϕk(x,y)‖2.μn⟶n→∞0.  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φi is chosen to maximize the average projection
of the member ui onto φi suitably normalized [25] (13)maxφi∈L2(Ω)|〈ui,φi〉|2‖φi‖22.  To construct numerically the POD
basis {φj},
we build the correlation matrix L, of
size M×M, composed of the inner products
of the snapshots. In L2(Ω) we have (14)Lk,j=〈Sk,Sj〉,(15)〈Sk,Sj〉=∫ΩSkSj*dxdy,  where * denotes complex conjugate transpose.A spectral
decomposition of the matrixL is
performed. The n largest eigenvalues {λ1,λ2,…,λn} of the matrix L are found and placed in descending order. Then the set of
eigenvectors are identified to be {v1,v2,…,vn}.The resulting orthonormal POD basis{φj} 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] (16)λk‖vk‖2=1,fork={1,2,…,n}.  Then, the POD basis set is formed according to (17)φk(x,y)=∑j=1Mvk,jSj(x,y)  with vk,j being the jth component of the eigenvector vk.
Solving (17) gives   nφk'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{αk} [10].Projecting (1) onto the POD basis yields
[10] (18)∫Ω∂∂tw(t,x,y)φk(x,y)dx=1Re(∫∂Ω(∇w(t,x,y)⋅n)φk(x,y)dA(x)−∫Ω∇w(t,x,y)⋅∇φk(x,y)dx)−(∫∂Ω(F(w)⋅n)φk(x,y)dA(x)−∫ΩF⋅∇φk(x,y)dx),   where the first term on the right-hand
side is (19)∫∂Ω(∇w(t,x,y)⋅n)φk(x,y)dA(x)=∫a1a2∂∂yw(t,x,b1)φk(x,b1)dx−∫a1a2∂∂yw(t,x,b2)φk(x,b2)dx−∫b0bend∂∂xw(t,a0,y)φk(a0,y)dy.  The Neumann boundary condition forces the
portion of the boundary integral over b0 to bend along aend to be 0, that is, (20)∫b0bend∂∂xw(t,aend,y)φk(aend,y)dy=0. The second boundary integral is decomposed
as follows (21)∫∂Ω(F(w)⋅n)φk(x,y)dA(x)=12∫b0bend(w(t,aend,y)2φk(aend,y)−f(y)2φk(a0,y))dy.    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] (22)∂∂yw(t,x,b1)≈ubottom(t)Ψbottom(x)−w(t,x,b1−h)h,∂∂yw(t,x,b2)≈w(t,x,b2+h)−utop(t)Ψtop(x)h,∂∂xw(t,a0,y)≈w(t,a0+h,y)−f(h)h.   After substitutions w(t,x,y) can be approximated as a linear combination of
POD modes when the αk's are solved in the following system model.
Then, the temporal model for the system is given by [10] (23)α˙=Aα+Bu+N(α)+F,α(0)=αo,  where α∈ℝn and the matrices A is n×n,B is n×2,N and F are both vectors n×1. The output equation will be simply chosen to be (24)y(t)=α(t).   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.First eight POD modes.
(a)(b)(c)(d)(e)(f)(g)(h)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(25)u1(t)=sin(3πt4),u2(t)=sin(3πt2).   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.Boundary control accuracy.
(a)(b)The goal of model reduction is to
construct another nonlinear system [10, 11] (26)α˙r=Arαr+Bru+Nr(αr)+Fr,  where αr∈ℝr and r<n,
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)α=Tαr,αr=Tα,TT=Ir,  where Ir is the r×r identity matrix, resulting in the following
reduced model: (28)α˙r=TATαr+TBu+TN(Tαr)+TF,y(t)=CTαr(t).  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(29)y(k)=∑ℓ=0∞Y(ℓ)u(k−ℓ),(30)Y(0)=D,Y(1)=CB,…,Y(k)=CAk−1B.  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] (31)Ryu(k)=∑ℓ=0∞Y(ℓ)Ruu(k−ℓ),  where the length of the data sequence
is (32)Ruu(k)=1m∑τ=0m−1u(τ)uT(k−τ),Ryu(k)=1m∑τ=0m−1y(τ)uT(k−τ).  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: (33)Pyu(k)=1m∑τ=0m−1Ryu(τ)e−j(2πk/m)τ,Puu(k)=1m∑τ=0m−1Ruu(τ)e−j(2πk/m)τ.  The ratio of the two power spectral densities is the
frequency response function and is given by [26] (34)G(zk)=Pyu(k)Puu(k).  Then, the final step is
to take the inverse Fourier transform to find the pulse response (Markov
parameter) matrices [26] (35)Yk:=Y(k)=∑ℓ=0∞G(zk)ej(2πℓ/m)k.  The Hankel matrix containing the
Markov parameters is of the following form: (36)Hlq=[Y1Y2⋯YqY2Y3⋯Yq+1⋮⋮⋮⋮YlYl+1⋯Yq+l−1].  The individual Yk's correspond to the following
sequence: (37)Y0=D,Y1=CB,…,Yk=CAk−1B.  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: (38)u1(t)=−sin⁡(0.55t2),u2(t)=−sin⁡(0.60t2).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)G:x˙(t)=Ax(t)+Bu(t),y(t)=Cx(t)+Du(t),  where x(t)∈ℝn is the state vector, u(t)∈ℝm is the input, and y(t)∈ℝp is the output.The first step in applying balanced
truncation is to compute a coordinate transformationM such that the controllability and observability gramians,
denoted Wc and Wo,
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]: (40)AWc+WcA*+BB*=0,A*Wo+WoA+C*C=0,  where A* denotes complex conjugate transpose.A balanced realization needs a
similarity transformationM such that
the transformed gramians are equal and satisfy [16, 17, 19] (41)W^o=Σ=W^c,  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 withn states, the controllability and observability matrices are n×n symmetric matrices and therefore solving each
one of them involves
finding n(n+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 thel-step observability and q-step controllability matrices [15] (42)Ol:=[CCACA2⋮CAl−1],Rq:=[BABA2B⋯Aq−1B] which
give rise to the l-step observability
and q-step controllability gramians (43)Wol:=Ol*Ol,Wcq:=RqRq*.   As
the numbers q and l approach infinity, these empirical gramians approach the true
gramians [15] (44)liml→∞Wol=Wo,limq→∞Wcq=Wc.  The
goal is to find a balancing transformation matrix M that will approximately balance the empirical gramians, that is, (45)W^cq:=MWcqM*=(M*)−1WolM−1=:W^ol=Σ.  The
matrix M can then be applied back to
the original system model to produce an approximately balanced realization.The product of thel-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 Hlq,
containing the Markov parameters CAkB,k=0,1,…,
of the system in the following way: (46)Hlq:=OlRq=[CBCAB⋯CAq−1BCABCA2B⋯CAqB⋮⋮⋱⋮CAl−1BCAlB⋯CAl+q−2B]  for
integers l and q chosen such that [15] (47)rank(Hlq)=rank(H(l+1)(q+j))=n,∀j≥1.  In
terms of the SVD decomposition of Hlq(48)Hlq=UΣV*=[U1U2][Σ1000][V1*V2*].   The
balancing transformation M is
constructed as [4] (49)M=RqV1Σ1−1/2.  A
straightforward computation shows (50)W^cq:=M−1WcqM*−1=M*WolM=:W^ol=Σ1.  Balanced
truncation can be realized by
the usual way; if σr≫σr+1 for some r, then we can partition Σ1 as (51)Σ1=[Σr00Σr+1],  where (52)Σr=diag(σ1,σ2,…,σr),Σr+1=diag(σr+1,σr+2,…,σn).  A
columwise conformal partition of U1 and V1(53)U1=[UrUn−r],V1=[VrVn−r]  yields
the immersion/projection pair [4] (54)Tr=RqVrΣr−1/2,Tr=Σr−1/2Ur*Ol,TrTr=Ir  and
from which a reduced-order r-dimensional
model with state matrices is deduced (55)Ar=TrATr,Br=TrB,Cr=CTr.  The
above construction only requires estimates of the Markov parameters (56)CAkB,k=0,1,…,l+q−1.  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(57)α˙=Aα+Bu+N(α)+F,α(0)=αo  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)α=Tαr,αr=Tα.   This
results in the following reduced-order model: (59)α˙r=Arαr+Bru+Nr(Tαr)+Fr,αr(0)=Tαo,Ar:=TArT,Br:=TB,Nr:=TN,Fr:=TF.  If (57) has a linearization around the steady state equilibrium αss,(60)α˙=Aℓα+Bℓu,α(0)=αo,  where (61)Aℓ=∂(Aα+Bu(t)+N(α)+F)∂α|α=αss,Bℓ=∂(Aα+Bu(t)+N(α)+F)∂u|u=B,   and the reduced system linearization
around the steady state equilibrium αss(62)α˙r=Aℓrαr+Bℓru,αr(0)=Tαo,  where (63)Aℓr=∂(Arαr+Bru(t)+Nr(Tαr)+Fr)∂α|α=αss,Bℓr=∂(Arαr+Bru(t)+Nr(Tαr)+Fr)∂u|u=Br.  The linearization of both models about equilibrium αss  is related by (64)Aℓr=TAℓT,Bℓr=TBℓ.  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
Figure7, 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
Figure8, we compare the full-order solution w(t,x,y) of Burgers’ equation with the solution based
on the 14th-order ERA model wr(t,x,y).
The figure shows that they behave similarly especially at the boundary where
control is applied.Full- and reduced-order
models' responses.
(a)(b)


## 6. Controller Design

AnH∞ 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 wref(x) 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 H∞ is based on the 14th-order reduced model.Figure 9
Closed-loop control system for tracking.From Figure9, 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 y(t), that is, (65)e(t):=α(t)−αref.     The objective of the H∞ 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 Teαref defined by (66)Teαref:=(I+PC)−1.   The H∞ control design reduces to the following
optimization. Find C such that the closed-loop system is stable and
minimizes (67)‖Teαref‖∞=sup‖αref‖2≤1‖e‖2.   Projecting wref(x) onto the POD basis yields the tracking
coefficients for the reduced-order model. After computing the H∞ 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.Flow initial
condition and reference.
(a)(b)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 H∞ 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.



---

*Source: 154956-2009-03-24.xml*