Feedback Control of a Thermal Fluid Based on a Reduced Order Observer ? Weiwei Hu ∗ John R. Singler ∗∗ Yangwen Zhang ∗∗ ∗

Institute for Mathematics and its Applications, University of Minnesota, Minneapolis, MN 55455 USA (email: [email protected]). ∗∗ Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409 (email: [email protected], [email protected]) Abstract: We discuss the problem of designing a feedback control law based on a reduced order observer, which locally stabilizes a two dimensional thermal fluid modeled by the Boussinesq approximation. We consider mixed boundary control for the Boussinesq equations in an open bounded and connected domain. In particular, the controllers are finite dimensional and act on a portion of the boundary through Neumann/Robin boundary conditions. A linear Luenberger observer is constructed based on point observations of the linearized Boussinesq equations. The current setting of the system leads to a problem with unbounded control inputs and outputs. Linear Quadratic Gaussian (LQG) balanced truncation is employed to obtain the reduced order model for the linearized system. The feedback law can be obtained by solving an extended Kalman filter problem. The numerical results show that the nonlinear system coupled with the reduced order observer through the feedback law is locally exponentially stable. 1. INTRODUCTION We consider the problem of feedback stabilization of a two dimensional thermal fluid. The transport of thermal energy in a viscous incompressible fluid can be modeled by the Boussinesq approximation, which couples the NavierStokes equations with the convection-diffusion equation for the temprature of the fluid. Feedback control for fluid flows is a very active area and has been widely studied; see, e.g., Choi et al. (1993); Burns et al. (1998); Wang (2003); Barbu et al. (2006); Lee and Choi (2006); Raymond (2006, 2007); Raymond and Thevenet (2010); Badra (2012); B¨ ansch and Benner (2012); Nguyen and Raymond (2015); Brunton and Noack (2015). Our current work focuses on the low order stabilizing feedback control design that is implementable in real time. Recent work in Burns et al. (2016) considered the LQR control design for a two dimensional Boussinesq equations. It is considered that the control inputs are finite dimensional and act on a portion of the boundary through Neumann/Robin boundary conditions. Dirichlet boundary conditions are imposed on the rest of the boundary. The standard Riccati-based feedback law can be derived. Numerical experiments show that the Riccati-based feedback law locally exponentially stabilizes an unstable steady state solution to the nonlinear Boussinesq system. However, full state feedback control is not practical for most flow control applications. In this work, we continue to use the setup in Burns et al. (2016), but consider the problem of stabilizing a possible ? W. Hu was supported in part by a USC Zumberge Individual Research Grant. J. Singler and Y. Zhang were supported in part by National Science Foundation grant DMS-1217122.

unstable steady state solution to the nonlinear Boussinesq equations by a feedback control law based on a reduced order observer. In particular, point observations are used for the output measurement and a linear Luenberger observer design is employed for the the state estimation of the linearized Boussinesq system. The current setting naturally leads to a problem with unbounded control inputs and outputs. To obtain a reduced order observer, we need an effective reduced order model for the unbounded input-output system. For model reduction of linear systems, a well-known class of methods with excellent properties are the balanced truncation algorithms (Antoulas, 2005; Zhou et al., 1996). The fundamental algorithm in this class, (standard) balanced truncation, is only applicable to exponentially stable systems. For unstable systems with no eigenvalues on the imaginary axis, a generalization of balanced truncation was introduced in Zhou et al. (1999). Ideas from this work have recently been used for generalized balanced truncation model reduction of unstable systems derived from a spatial discretization of a linear PDE system; see, e.g., Ahuja and Rowley (2010); Benner et al. (2016); Flinois et al. (2015) and the references therein. Another balanced truncation approach that is directly applicable to unstable systems is LQG balanced truncation. This method is derived from the algebraic Riccati equations arising in LQG feedback control, and therefore LQG balanced truncation has been frequently used for model reduction in feedback control applications. This method has been applied to compute reduced order controllers for many PDE systems; see, e.g., Batten and Evans (2010); Benner and Heiland (2015); Breiten and Kunisch (2015, 2016); Evans (2003); Singler and Batten (2009). We use

this approach here to develop a linear reduced order controller for the nonlinear Boussinesq system. We note that a reduced order stabilizing feedback controller for the Navier-Stokes equations was designed using LQG balanced truncation in Benner and Heiland (2015); however, boundary control was not considered in that work. The main contribution of this work is the investigation of the performance of a low order LQG balanced feedback controller for the Boussinesq equations with boundary control. 2. THE MODEL Let Ω be an open bounded and connected domain with a Lipchitz boundary Γ. The Boussinesq approximation is given by 1 Gr ¯ 2 θ + fv , ∂t v + v · ∇v = div (∇v + (∇v)T ) − ∇p + e Re Re (1) div v = 0, (2) 1 ∂t θ + v · ∇θ = ∆θ + fθ , (3) ReP r where v(x, t) is the velocity, p(x, t) is the pressure, θ(x, t) is the fluid temperature, Re is the Reynolds number, Gr is the Grashof number, P r is the Prandtl number, and ¯ = [0, 1]T is the gravitational force direction. We assume e fv is a time independent external body force and fθ is a time independent heat source density. Consider a 2D domain shown in Figure 1. Assume that the controlled

Inlet

ΓI

ΓO

ΓH Radiant heating strip

1 ∂θ |Γ = 0, ReP r ∂n O m X 1 ∂θ |Γ = (bθHi |ΓH )(x)uθHi (t), ReP r ∂n H i=1

(6) (7) θ|ΓD = 0,

(8) where T (v, p) is the fluid Cauchy stress tensor defined by 1 (∇v + (∇v)T ) − qI. T (v, q) = Re For uvi = uθj = 0, let (ve , pe , θe ) be a steady-state (equilibrium) solution to equations (1)–(3). Notice that for large Reynolds numbers or strong external body forces, the steady-state solution can be unstable. We introduce the new variables w = v − ve , T = θ − θe and q = p − pe . Then the translated system is given by ∂w 1 = ∆w − w · ∇ve − ve · ∇w − w · ∇w ∂t Re Gr ¯ 2 T, − ∇q + e (9) Re ∇ · w =0, (10) ∂T 1 = ∆T − w · ∇θe − ve · ∇T − w · ∇T. (11) ∂t ReP r Let x(t) = [w(t), T (t)]T . Then the controlled translated equations (1)–(7) can be rewritten as ˙ x(t) = Ae x(t) + Bu(t) + F(x), (12) where Ae is the translated linearized system operator, B is the boundary control operator and F is the nonlinear mapping. The details for the formulation of (12) can be found in Burns et al. (2016). By the linearized theory of hydrodynamic stability in Sattinger (1973), the stability of (ve , pe , θe ) is determined by the spectrum of Ae associated with the boundary conditions. (13) T

Fig. 1. 2D domain airflow is coming in through the inlet ΓI , which is a subset of the boundary Γ, with Robin boundary control for both velocity and temperature. The airflow exits at the outlet ΓO with stress-free fluid and natural (or unforced) convective flux boundary conditions. In addition, there is a radiant heating strip, denoted by ΓH , on the floor with Neumann boundary control for temperature. We impose no slip boundary conditions for the velocity on Γf = Γ \ (ΓI ∪ ΓO ) and zero Dirichlet boundary condition for temperature on ΓD = Γ \ (ΓI ∪ ΓO ∪ ΓH ). Here the boundaries ΓI , ΓO and ΓH are disjoint. The boundary conditions can be formulated as follows. m X (T (v, q) · n + αv)|ΓI = (bvi |ΓI )(x)uvi (t), (4) i=1

v|Γf = 0,

X 1 ∂θ + βθ)|ΓI = (bθIi |ΓI )(x)uθIi (t), ReP r ∂n i=1

Linearizing system (12) yields ˙ z(t) = Ae z(t) + Bu(t),

Outlet

T (v, q) · n|ΓO = 0,

m

(

(5)

in a Hilbert space H, where z(t) = [w(t), T (t)] . The control space is U = R3m . Moreover, consider the output measurement of the linearized system (13) by point observations y(t) = Cz(t) = [z(ξ1 , t), z(ξ2 , t), . . . , z(ξn , t)]T ∈ R3n , (14) where C = [δ(x−ξ1 ), δ(x−ξ2 ), . . . , δ(x−ξn )]T is the point observation operator and ξi ∈ Ω, i = 1, 2 . . . , n, are the points of observation. Next we apply the linear Luenberger observer design to system (13)–(14). Consider z˙ c = Ae zc + Bu + L(Czc − y), (15) where L is called the filtering operator. Note that Ae generates an analytic C0 -semigroup {T (t)}t≥0 on H. The resolvent set of Ae contains a sector. Thus, there are at most a finite number of eigenvalues of Ae in the right complex half-plane {λ ∈ C : Reλ ≥ 0}. Therefore, there exists a real number λ0 ∈ ρ(Ae ), sufficiently large, such that λ0 I − Ae is a strictly positive operator and the fractional powers (λ0 I − Ae )σ are well defined for 0 ≤ σ ≤ 1. In addition, since B is a Neumann/Robin type boundary

1

control operator, (λ0 −Ae )−( 4 +) B ∈ L (R3m ; H) for some ¯  > 0 by the trace theorem. Also, since H 2γ (Ω) ⊂ C(Ω) −( 12 +) n for 2γ > d/2, we have C(λ0 − Ae ) ∈ L (H; R ) for some  > 0 and d = 2. By using Banach fixed point theorem, we can prove that if there exist operators K and L such that Ae + BK and Ae + LC generate analytic C0 -semigroups, which are exponentially stable, then the full nonlinear Boussinesq system (12) coupled with the linear Luenberger observer (15) through the feedback law u = Kzc (16) is locally exponentially stable on an appropriate state space. The detailed proof will be provided in a future paper. The objective of this paper is to construct a reduced order Luenberger observer such that the coupled nonlinear system      ˙ x(t) Ae BKr x(t) = z˙ rc (t) −Lr C Are + Br Kr + Lr Cr zrc (t)   F(x(t)) + . (17) 0

(7) Set Ar + Br Kr + Lr Cr = Sr (A + BK + LC)Tr , Kr = KTr and Lr = Sr L. Although the matrices Π and P are positive definite in theory, in practice they can have low numerical rank and and the numerical Cholesky factorization in step 2 above may fail. Here is a simple way to obtain the factors. First, compute Π and P and the gains K = B T Π and L = P C T . Then rewrite the Riccati equations as Lyapunov equations: (A − BK)T Π + Π(A − BK) + K T K + C T C = 0, (A − LC)P + P (A − LC)T + LLT + BB T = 0. Next, we use the Matlab Cholesky factor Lyapunov equation solver lyapchol to compute Rc and Rf . (This is equivalent to taking one step of a Newton iteration for each Riccati equation (Kleinman, 1968).) A short Matlab code to compute the LQG balancing transformation is below: function [R,S,mu,K,L] = LQGbal(A,B,C) % Computes a vector of LQG characteristic values mu, % LQG balancing transformation matrices R and S, % and the feedback gain matrices K and L % (R, S, and mu are not truncated)

is locally asymptotically stable. 3. LQG BALANCED TRUNCATION In this section, we discuss LQG balanced truncation for the input-output system Σ(Ae , B, C) defined by (13)– (14). The finite element method is used to approximate Σ(Ae , B, C) in space. In particular, we use Taylor-Hood finite element pair for the simulation of Navier-Stokes equations, where linear elements are used for pressure and quadratic elements are used for velocity. Quadratic elements are also used for temperature with respect to the same triangulation. Moreover, the penalty method is used to relax the incompressibility condition so that the pressure can be solved in terms of velocity. For more details on the finite element method used here, see Hu (2012). Let A, B and C be the matrix approximations of Ae , B, C. We apply LQG balanced truncation to derive the reduced order observer for the coupled system (17) linearized about zero. We sketch the process as follows. (1) Compute the solutions Π and P of the normalized control and filter algebraic Riccati equations AT Π + ΠA − ΠBB T Π + C T C = 0, and AP + P AT − P C T CP + BB T = 0. (2) Compute the Cholesky factorizations Π = RcT Rc and P = RfT Rf . (3) Compute the singular value decomposition of Rc RfT such that Rc RfT = U M V T , where U and V are unitary and M is diagonal with the LQG characteristic values µ1 ≥ µ2 ≥ · · · ≥ 0 along the diagonal in decreasing order. −1/2 −1/2 (4) Form the matrix M −1/2 = diag(µ1 , . . . , µn ). (5) Form the balancing transformation T = RfT V M −1/2 and its left inverse S = M −1/2 U T Rc (6) Choose r and truncate: Let Tr be the matrix consisting of the first r columns of T , and let Sr be the matrix consisting of the first r rows of S.

% Compute Riccati solutions and gains Pi = are(A,B*B',C'*C); K = B'*Pi; P = are(A',C'*C,B*B'); L = P*C'; % Take one Newton step and factor Riccati solutions Rc = lyapchol((A - B*K)',[K' C']); Rf = lyapchol(A - L*C,[L B]); % Compute LQG characteristic values, transformation [U,Mu,V] = svd(Rc*Rf'); mu = diag(Mu); Mumh = diag(mu.ˆ(-1/2)); R = Rf'*V*Mumh; S = Mumh*U'*Rc;

4. NUMERICAL SIMULATIONS Numerical simulations are conducted in this section to demonstrate our method introduced previously. Let the domain of the 2D model be Ω = [0, 1] × [0, 1]. The inlet, outlet and the radiant heating strip are defined by ΓI = {x = 1, 0.7 ≤ y ≤ 0.9}, ΓO = {x = 0, 0.1 ≤ y ≤ 0.4}, and ΓH = {y = 0, 0.4 ≤ x ≤ 0.6}. In the numerical experiments, we set Re = 100 and P r = 0.7 for air. We assume that the contribution of buoyancy is Gr 10 relatively strong by setting Re 2 = 9 > 1. Furthermore, we consider the external body force for velocity given by fv1 (x, y) = fv2 (x, y) = 0 and the heat source density given by fθ (x, y) = 7 sin(2πx) cos(2πy), which is shown in Figure (2a). In this setting the airflow is only driven by the heat. We use a mesh size of h = 0.05 and the penalty parameter ε = 10−5 in the simulations. We set the parameters α = β = 1 in the boundary conditions, which is sufficient to demonstrate the main idea. For these parameters, we compute the steady-state solution (ve1 , ve2 , θe ) using Newton’s iterative method. In particular, we take the steady-state solution associated with fθ0 (x, y) = 5 sin(2πx) cos(2πy) as the initial point for the iteration. We denote it by (v01 , v02 , θ0 ). As indicated in Figure (2a), the cold air in the upper left corner descends and the warm air in the lower right corner rises. This results in the velocity field shown in Figure (2b), circulating in a counter clockwise fashion, which in turn causes the temperature distribution to advect counter clockwise as we

can see in Figure (2c). Moreover, the steady state solution (ve1 , ve2 , θe ) is unstable since Ahε has a pair of simple complex conjugate eigenvalues λo = 0.15214 ± 0.84910i in the positive right-half plane. The dimension of Ahε in this case is 4624 × 4624.

14 13 12 11 10

1

9

6

0.9

8

0.8

4 7

0.7 2

6

0

5

−2

4 0

0.6 0.5 0.4 0.3

20

30

40

50

−4

0.2 0.1 0 0

L2−norm of the nonlinear open−loop solution T 10

−6 0.2

0.4

0.6

0.8

Fig. 3a. L2 -norm of the uncontrolled temperature T

1

2.8

Fig. 2a. Heat Source fθ

2.6 2.4

1 2.2

0.9

2

y

0.8 0.7

1.8

0.6

1.6

0.5

1.4

0.4 1.2

L2−norm of the nonlinear open−loop solution w1

0.3 1 0

0.2

10

20

30

40

50

0.1 0 −0.2

0

0.2

0.4

0.6

0.8

1

1.2

x

Fig. 3b. L2 -norm of the uncontrolled horizontal velocity w1

Fig. 2b. Vector field of steady state velocity ve 2.8 2.6

1 2.5 0.9

2.4

2 0.8 1.5 0.7

2.2

1

0.6

0.5

0.5

0

2 1.8

0.4

−0.5

0.3

−1

0.2

−1.5

1.6 1.4

−2

0.1

−2.5 0 0

0.2

0.4

0.6

0.8

1.2

L2−norm of the nonlinear open−loop solution w

1 1 0

10

20

30

2

40

50

Fig. 2c. Steady state temperature θe For the nonlinear open-loop and closed-loop simulations, we continue to use (v01 , v02 , θ0 ) as the initial condition. Moreover, backward Euler discretization is used for approximating the system over the time interval [0, 50] with a time step size ∆t = 0.05. We present the L2 -norms of the open-loop temperature and velocity in Figures (3a)– (3c). We can see that the nonlinear open-loop oscillates and does not converge. Since the unstable eigenvalues are simple, one control input vector b = [bv , bθI , bθH ]T can stabilize the linearized system if hb, ReΦiΓ 6= 0 and hb, ImΦiΓ 6= 0, (18) T

for Φ ∈ Ker(λO1 − Ahe ) (see Theorem 3.2 in Burns et al. (2016)). In this work, we use the same boundary input functions as in Burns et al. (2016). We take

Fig. 3c. L2 -norm of the uncontrolled vertical velocity w2 (bv |ΓI )(y) = (e

−

0.0001 [(0.7−y)(0.9−y)]2

(bθI |ΓI )(y) = 0.2e and

, 0)T ,

−0.00001 [(0.7−y)(0.9−y)]2

,

−0.00001

bθH |ΓH )(x) = 0.4e [(0.4−x)(0.6−x)]2 . It can be verified numerically that (18) holds and therefore (Ahe , Bh ) is stabilizable. The control space in this case is U = R3 . In the current setting, there are nine functional gains (Burns et al. (2016)). To find the optimal sensor locations, we first compute the best rank-1 approximation of the nine functional gains using the singular value decomposition ¯ y). Let ρ(x, y) = |k(x, ¯ y)| (SVD) and denote it by k(x, and the plots about ρ(x, y) are shown in Figures 4–5. As

we can see in Figure 4, the peak regions of the functional gains represent the major system response compared with the rest of the flat regions. Therefore, it is plausible to focus the observations on those regions.

1

0.8

y

0.6

0.4

15 0.2

10 0 0

0.2

0.4

0.2

0.4

5

x

0.6

0.8

1

0.6

0.8

1

Fig. 6. Three sensors

0 1 1 0.5

1

0.5 0 0

y

x

0.8

¯ y)| Fig. 4. Density function ρ(x, y) = |k(x, y

0.6

0.4 1

0.2

0.8

0 0

x

y

0.6

0.4

Fig. 7. Six sensors

0.2

0 0

0.2

0.4

x

0.6

0.8

closed-loop system with a four dimensional reduced order observer decay to zero exponentially as shown in Figure 9. Note that the full order observer is 4624 dimensional.

1

Fig. 5. Flat view of the density function

18

Here we use the Centroidal Voronoi Tessellations (Du et al., 1999) to find the centroids associated with ρ(x, y). The sensors will be placed in the locations of the centroids (Faulds and King, 2000). Notice that it is possible for the system undetectable if the sensors are located on the zeros of the unstable eigenfunctions of Ae . For a given number n and density function ρ(x), there exists a Voronoi tessellation Ω = ∪ni=1 ωi . Let ξi∗ denote the centroids of the corresponding Voronoi regions ωi , for i = 1, 2, . . . , n. Then ξi∗ is given by R xρ(x) dx ∗ , i = 1, 2, . . . , n. ξi = Rωi ρ(x) dx ωi We have the configurations for three and six sensor locations shown in Figures 6–7 and the sensors are marked by black dots associated with different colored Voronoi regions. We used the six sensor locations shown in Figure 7 to construct the Luenberger observer in the rest of our simulations. Their locations are given in Table 1 below.

16

14

12

10

8

6 0

L2−norm of the overall response of the nonlinear open−loop 10

20

30

40

50

Fig. 8. L2 -norm of the overall response of the nonlinear open-loop 18

L2−norm of the nonlinear closed−loop response L2−norm of the reduced order observer response

16 14 12

xc yc

0.58522 0.18955

0.45388 0.81047

0.22420 0.19358

0.78962 0.81111

0.29502 0.051416

0.22611 0.42120

Table 1. Sensor locations

10 8 6 4

Finally, we compare the overall responses of the nonlinear open-loop system with the coupled nonlinear closed-loop system and the reduced order observer. Figure 8 provides the L2 -norm of the overall response of the nonlinear openloop over the time interval [0, 50], which oscillates. However, the L2 -norms of the overall responses of the nonlinear

2 0 0

10

20

30

40

50

Fig. 9. L2 -norms of the overall responses of the nonlinear closed-loop and the reduced order observer

5. CONCLUSION We discussed the feedback control problem based on a reduced order observer for the Boussinesq system. We also suggested a systematic way to locate the sensors for the point observations. Since the system is unstable and has unbounded inputs and outputs, LQG balanced truncation is employed to obtain a reduced order observer. Numerical simulation based on a coarse mesh shows that the feedback control based on a low order observer is rather effective to stabilize the unstable steady state solution to the nonlinear Boussineq system. Computations on finer meshes will be considered in our future paper and the rigorous proofs of the numerical simulations will be addressed. REFERENCES Ahuja, S. and Rowley, C. (2010). Feedback control of unstable steady states of flow past a flat plate using reduced-order estimators. Journal of Fluid Mechanics, 645, 447–478. Antoulas, A.C. (2005). Approximation of Large-Scale Dynamical Systems. SIAM, Philadelphia, PA. Badra, M. (2012). Abstract settings for stabilization of nonlinear parabolic system with a Riccati-based strategy. Application to Navier-Stokes and Boussinesq equations with Neumann or Dirichlet control. Discrete and Continuous Dynamical Systems-Series A, 252(09), 5042–5075. B¨ ansch, E. and Benner, P. (2012). Stabilization of incompressible flow problems by Riccati-based feedback. In Constrained Optimization and Optimal Control for Partial Differential Equations, 5–20. Springer. Barbu, V., Lasiecka, I., and Triggiani, R. (2006). Tangential Boundary Stabilization of Navier-Stokes Equations. Mem. Amer. Math. Soc., 181(852). Batten, B.A. and Evans, K.A. (2010). Reduced-order compensators via balancing and central control design for a structural control problem. Internat. J. Control, 83(3), 563–574. Benner, P. and Heiland, J. (2015). LQG-balanced truncation low-order controller for stabilization of laminar flows. In Active Flow and Combustion Control 2014, 365–379. Springer. Benner, P., Saak, J., and Uddin, M.M. (2016). Balancing based model reduction for structured index-2 unstable descriptor systems with application to flow control. Numerical Algebra, Control and Optimization, 6(1), 1– 20. Breiten, T. and Kunisch, K. (2015). Feedback stabilization of the Schl¨ ogl model by LQG-balanced truncation. In Proceedings of the European Control Conference, 1171– 1176. Breiten, T. and Kunisch, K. (2016). Compensator design for the monodomain equations. ESAIM: Control, Optimisation and Calculus of Variations. To appear. Available at http://math.unigraz.at/mobis/publications/SFB-Report-2014-017.pdf. Brunton, S.L. and Noack, B.R. (2015). Closed-loop turbulence control: Progress and challenges. Applied Mechanics Reviews, 050901. Burns, J.A., King, B.B., and Rubio, D. (1998). Feedback control of a thermal fluid using state estimation. Int. J. Comput. Fluid Dyn., 11(1-2), 93–112.

Burns, J., He, X., and Hu, W. (2016). Feedback stabilization of a thermal fluid system with mixed boundary control, to appear. Computers & Mathematics with Applications. Choi, H., Temam, R., Moin, P., and Kim, J. (1993). Feedback control for unsteady flow and its application to the stochastic Burgers equation. Journal of Fluid Mechanics, 253, 509–543. Du, Q., Faber, V., and Gunzburger, M. (1999). Centroidal Voronoi tessellations: applications and algorithms. SIAM Review, 41(4), 637–676. Evans, K.A. (2003). Reduced Order Controllers for Distributed Parameter Systems. Ph.D. thesis, Virginia Polytechnic Institute and State University. Faulds, A. and King, B. (2000). Sensor location in feedback control of partial differential equation systems. In Proceedings of the 2000 IEEE International Conference on Control Applications, 536–541. Flinois, T.L., Morgans, A.S., and Schmid, P.J. (2015). Projection-free approximate balanced truncation of large unstable systems. Physical Review E, 92(023012). Hu, W. (2012). Approximation and control of the Boussinesq equations with application to control of energy effcient building systems. Ph.D. thesis, Virginia Polytechnic Institute and State University, Blacksburg, Virginia. Kleinman, D.L. (1968). On an iterative technique for Riccati equation computations. IEEE Trans. Autom. Control, AC-13, 114–115. Lee, H.C. and Choi, Y. (2006). Analysis and approximation of linear feedback control problems for the Boussinesq equations. Computers & Mathematics with Applications, 51(5), 829–848. Nguyen, P.A. and Raymond, J.P. (2015). Boundary stabilization of the Navier–Stokes equations in the case of mixed boundary conditions. SIAM Journal on Control and Optim., 53(5), 3006–3039. Raymond, J.P. (2006). Feedback boundary stabilization of the two-dimensional Navier-Stokes equations. SIAM J. Control Optim., 45(3), 790–828. Raymond, J.P. (2007). Feedback boundary stabilization of the three-dimensional incompressible Navier–Stokes equations. Journal de math´ematiques pures et appliqu´ees, 87(6), 627–669. Raymond, J.P. and Thevenet, L. (2010). Boundary feedback stabilization of the two dimensional Navier-Stokes equations with finite dimensional controllers. Discrete and Continuous Dynamical Systems-Series A, 27(3), 1159–1187. Sattinger, D.H. (1973). The mathematical problems of hydrodynamic stability. Topics in Stability and Bifurcation Theory. Lecture Notes in Mathematics, 309, 125–140. Singler, J.R. and Batten, B.A. (2009). A comparison of balanced truncation methods for closed loop systems. In Proceedings of the American Control Conference, 820– 825. Wang, G. (2003). Stabilization of the Boussinesq equation via internal feedback controls. Nonlinear Analysis: Theory, Methods & Applications, 52(2), 485–506. Zhou, K., Doyle, J.C., and Glover, K. (1996). Robust and Optimal Control. Prentice-Hall. Zhou, K., Salomon, G., and Wu, E. (1999). Balanced realization and model reduction for unstable systems. Internat. J. Robust Nonlinear Control, 9(3), 183–198.