Feedback Stabilization of the Two-Dimensional Navier-Stokes Equations by Value Function Approximation

The value function associated with an optimal control problem subject to the Navier-Stokes equations in dimension two is analyzed. Its smoothness is established around a steady state, moreover, its derivatives are shown to satisfy a Riccati equation at the order two and generalized Lyapunov equations at the higher orders. An approximation of the optimal feedback law is then derived from the Taylor expansion of the value function. A convergence rate for the resulting controls and closed-loop systems is demonstrated.


Introduction
In this work we continue our investigations of the value function associated with infinite-horizon optimal control problems of partial differential equations, that we initiated in [12,14]. We consider a stabilization problem of the Navier-Stokes equations in dimension two and focus on the regularity of the value function and its characterization as a solution to a Hamilton-Jacobi-Bellman (HJB) equation. This task has been the subject of tremendous research, for optimal control problems of a general structure, in general associated with finite-dimensional dynamical systems. The use of the notion of viscosity solutions has allowed to deal with the low regularity of the value function. In the present paper, to the contrary, we show that the value function is smooth and that the HJB equation is satisfied in the strict sense, in a neighborhood of the steady state. Moreover, we show that the derivatives of the value function, at the steady state, are solutions to an algebraic Riccati equation (for the order 2) and to linear equations, called generalized Lyapunov equations, for the higher orders. The main interest of these results is the fact that polynomial feedback laws can be derived from Taylor approximations of the value function. Moreover their efficiency can be analyzed.
From a methodological point of view, we mainly follow the techniques that we laid out for bilinear optimal control problems (such as control problems of the Fokker-Planck equation) in [14] and [12]. The Navier-Stokes control system considered here requires a different functional analytic treatment. In fact, the involved nonlinear terms must be tackled with different estimates, to guarantee, for example, the well-posedness of the closed-loop system. They also lead to different generalized Lyapunov equations. Moreover, from the point of view of open-loop control of the Navier-Stokes equation, this paper contains results on infinite-horizon optimal control which are not readily available elsewhere.
Feedback stabilization of the Navier-Stokes equations has been and still is an active topic of research. Among the numerous works, we refer to, e.g., [3,4,7,21,31], and the references therein.
The technique of approximation of the value function with a Taylor expansion dates back to [2,28], where optimal control problems associated to finite-dimensional control systems were investigated. We also quote follow-up work, for instance in [1,5,29]. For infinite-dimensional problems, we are only aware of [12,14]. In [13], the numerical solvability of the Lyapunov equations has been addressed. Model reduction techniques based on balanced truncation have been used in this reference to cope with the curse of dimensionality encountered when dealing with PDE controlled systems.
Let us next specify the problem which will be investigated in this paper. Throughout Ω ⊂ R 2 denotes a bounded domain with Lipschitz boundary Γ. Given two vector valued functions ϕ and ψ, we consider a solution (z,q) of the stationary Navier-Stokes equations −ν∆z + (z · ∇)z + ∇q = ϕ in Ω, divz = 0 in Ω, z = ψ on Γ. (1) Our goal is to find a control u such that the solution (z, q) to the instationary Navier-Stokes equations ∂z ∂t = ν∆z − (z · ∇)z − ∇q + ϕ +Bu in Ω × (0, T ), div z = 0 in Ω × (0, T ), z = ψ on Γ × (0, T ), z(0) =z + y 0 (2) is stabilized aroundz, i.e., lim t→∞ z(t) =z provided the initial perturbation y 0 is small in an appropriate sense. The control operatorB will be defined below. Throughout this work, we assume that div y 0 = 0. Our results are concerned with feedback stabilization of (2) and for this purpose, we consider new state variables (y, p) := (z, q) − (z,q) which satisfy the following generalized Navier-Stokes equations ∂y ∂t = ν∆y − (y · ∇)z − (z · ∇)y − (y · ∇)y − ∇p +Bu in Ω × (0, T ), div y = 0 in Ω × (0, T ), y = 0 on Γ × (0, T ), y(0) = y 0 . (3) The following sections are structured as follows. The problem statement and fundamental results on the state-equation on the time interval [0, ∞) are given in Section 2. Section 3 contains the existence theory of optimal controls, the adjoint equation, sensitivity analysis, and differentiability of the value function. The characterization of all higher order derivatives of the value as solutions to generalized Lyapunov equations are provided in Section 4. Section 5 contains the Taylor expansion of the value function, and estimates for convergence rates between the optimal solution and its approximation on the basis of feedback solutions obtained from derivatives of the value function. The paper closes with a very short outlook.
Notation. For Hilbert spaces V ⊂ Y with dense and compact embedding, we consider the Gelfand triple V ⊂ Y ⊂ V where V denotes the topological dual of V with respect to the pivot space Y . Given T ∈ R we consider the space For T = ∞, the space W (0, T ) will be denoted by W ∞ . For vector-valued functions f ∈ (L 2 (Ω)) 2 , we use the notation f ∈ L 2 (Ω). Elements f ∈ L 2 (Ω) will be denoted in boldface and are distinguished from real-valued functions g ∈ L 2 (Ω). Similarly, we use H 2 (Ω) for the space (H 2 (Ω)) 2 and H 1 0 (Ω) for (H 1 0 (Ω)) 2 . Given a closed, densely defined linear operator (A, D(A)) in Y , its adjoint (again considered as an operator in Y ) will be denoted with (A * , D(A * )). Let us introduce some notation that will be needed for the description of polynomial mappings. For δ ≥ 0 and a Hilbert space Y , we denote by B Y (δ) the closed ball in Y with radius δ and center 0. For k ≥ 1, we make use of the following norm: Given a Hilbert space Z, we say that T : Y k → Z is a bounded multilinear mapping (or bounded multilinear form when Z = R) if for all i ∈ {1, . . . , k} and for all (z 1 , . . . , z i−1 , z i+1 , . . . , z k ) ∈ Y k−1 , the mapping z ∈ Y → T (z 1 , . . . , z i−1 , z, z i+1 , . . . , z k ) ∈ Z is linear and The set of bounded multilinear mappings on Y k will be denoted by M(Y k , Z). For all T ∈ M(Y k , Z) and for all (z 1 , . . . , z k ) ∈ Y k , Given a bounded multilinear form T and z 2 , . . . , z k ∈ Y k−1 , we denote by T (·, z 2 , . . . , z k ) the bounded linear form z 1 ∈ Y → T (z 1 , . . . , z k ) ∈ R. It will be very often identified with its Riesz representative. Note that Bounded multilinear mappings T ∈ M(Y k , Z) are said to be symmetric if for all z 1 , . . . , z k ∈ Y k and for all permutations σ of {1, . . . , k}, Finally, given two multilinears mappings T 1 ∈ M(Y k , Z) and T 2 ∈ M(Y , Z), we denote by T 1 ⊗ T 2 the bounded multilinear form defined by Throughout the manuscript, we use M as a generic constant that might change its value between consecutive lines.

Problem formulation 2.1 Abstract Cauchy problem
In this section, we formulate system (3) as an abstract Cauchy problem on a suitable Hilbert space and, subsequently, define the stabilization problem of interest. This procedure is quite standard, see, for instance, [3,4,21,31,33] for details. We introduce the spaces Y := y ∈ L 2 (Ω) | div y = 0, y · n = 0 on Γ , It is well-known that Y is a closed subspace of L 2 (Ω). Moreover, we have the orthogonal decomposition where see, e.g., [33, page 15]. By P we denote the Leray projector P : L 2 (Ω) → Y which is the orthogonal projector in L 2 (Ω) onto Y . Following, e.g., [3], we define a trilinear form s by and a nonlinear operator F : V → V by For the bilinear mapping associated with the linearization of F , we introduce the operator The Oseen-Operator is then defined by The following well-known results (see, e.g., [3], [33,Lemma III.3.4]) concerning s and N will be used frequently throughout the paper.
Proposition 1. The following properties hold for N and s: With the previous result, we obtain similar properties for time-varying functions y, z, w.
where M is the constant given by Proposition 1.
Proof. Using Proposition 1 and Cauchy-Schwarz inequality (two times), we obtain that The two inequalities easily follow.
Forz ∈ V , we further introduce the Stokes-Oseen operator A via Considered as operator in L 2 (Ω) the adjoint A * , as operator in L 2 (Ω), can be characterized by (see, e.g., [31]) We note that as a consequence of Proposition 1, the operator A can be extended to a bounded linear operator from V to V in the following manner: Note that this extension is consistent, since by definition of the Leray projector P , we have P y, w Y = y, w Y for all y ∈ L 2 (Ω) and for all w ∈ V . Similarly, A * can be extended to a bounded linear operator from V to V . The control operator is chosen to satisfyB ∈ L(U, L 2 (Ω)). We further define B := PB ∈ L(U, Y ). The controlled state equation (3) can now be formulated as the abstract control systeṁ where the pressure p is eliminated. We can finally formulate the stabilization problem as an infinite-horizon optimal control problem: inf y∈W∞ u∈L 2 (0,∞;U ) J(y, u), subject to: e(y, u) = (0, y 0 ) (P ) Let us note that e :

Assumptions and first properties
Throughout the article we assume that the following assumptions hold true.
Assumption A2. There exists an operator K ∈ L(Y, U ) such that the semigroup e (A−BK)t is exponentially stable on Y .
Assumption A1 concerning the exponential feedback stabilizability of the Stokes-Oseen operator is well investigated. We refer e.g. to [3] where finite-dimensional feedback operators are constructed on the basis of spectral decomposition or alternatively by Riccati theory. In this case A2 can be satisfied with U = R m , for m appropriately large. Alternatively, we can rely on exact controllability results as obtained in [20]. They imply that the finite cost criterion holds. We can then rely on classical results, see, e.g., [30] which guarantee the existence of a stabilizing feedback operator.
Let us discuss some important consequences of the above definitions and assumptions.
Consequence C2. For all y 0 ∈ Y , for all f ∈ L 2 (0, ∞; V ), and for all T > 0, there exists a unique solution y ∈ W (0, T ) to the systeṁ This solution satisfies with a continuous function c. Assuming that y ∈ L 2 (0, ∞; Y ), we consider the equivalent equatioṅ where f λ ∈ L 2 (0, ∞; V ). By (18), the operator A λ generates an analytic, exponentially stable, semigroup on Y satisfying e A λ t Y ≤ e −δt for some δ > 0 independent of t ≥ 0, see [9, Theorem II. 1.2.12]. It follows that y ∈ W ∞ and there exists M λ such that with This estimate is obtained by adapting [9, Corollary II.3.2.1] and [9, Theorem II.3.2.2] from the temporal domain (0, T ) to (0, ∞), which can be achieved using the exponential stability of e A λ t .
Lemma 4. There exists a constant C > 0 such that for all δ ∈ [0, 1] and for all y and z ∈ W ∞ with y W∞ ≤ δ and z W∞ ≤ δ, it holds that Proof. We have The assertion now easily follows from Corollary 3.
The following lemma is formulated for an abstract generator A s of an analytic semigroup on Y . It will subsequently be used to address the asymptotic behavior of the nonlinear system (15).
Lemma 5. Let A s be the generator of an exponentially stable analytic semigroup e Ast on Y such that (18) holds. Let C denote the constant from Lemma 4. Then there exists a constant M s such that for all y 0 ∈ Y and f ∈ L 2 (0, ∞; V ) with has a unique solution y in W ∞ , which moreover satisfies Proof. Since the semigroup e Ast is exponentially stable on Y , it follows that for all (y 0 , g) ∈ Y × L 2 (0, ∞; V ) the systemż = A s z + g, z(0) = y 0 has a unique solution z ∈ W ∞ . Moreover, there exists a constant M s such that Without loss of generality we can assume that M s ≥ 1 2C . We claim that the constant M s is the one announced in the assertion. This will be shown by a fixed-point argument applied to the system (20). For this purpose, let us define M = {y ∈ W ∞ | y W∞ ≤ 2M s γ} and let us define the mapping Z : If there exists a unique fixed point of Z, then it is a unique solution of (20) in M. With C and M s given, we shall use Lemma 4 with δ = 2M s γ ≤ 1 2CMs ≤ 1. Together with (21), it follows that Again by (21) and Lemma 4 we obtain In other words, Z is a contraction in M and therefore, there exists a unique y ∈ M such that Z(y) = y. Regarding uniqueness in W ∞ , consider two solutions y, z ∈ W ∞ . For the difference e := y − z it then holdsė Multiplying with e and taking inner products yields Since A s satisfies an inequality of the form (18), we have where α ≥ 0 and β > 0. Using Proposition 1 and Young's inequality we further obtain Taking γ and δ sufficiently large, it holds that Since y, z ∈ W ∞ and e(0) = 0, with Gronwall's inequality, we conclude that e(t) = 0 for all t ≥ 0. Hence, y = z showing the uniqueness of the solution in W ∞ .
The following two corollaries are consequences of Lemma 4 and Lemma 5. The constant C which is employed is the one given by Lemma 4. Corollary 6. There exists a constant M K > 0 such that for all y 0 ∈ Y and for all f ∈ L 2 (0, ∞; V ) with has a unique solution y ∈ W ∞ satisfying Proof. By assumption A2, there exists K such that A − BK generates an exponentially stable, analytic semigroup on Y . The result then follows by applying Lemma 5 to the systeṁ and by defining u = −Ky.
In the following corollary, we assume without loss of generality that the constant M λ given by then y ∈ W ∞ and it holds that Proof. Since y ∈ L 2 (0, ∞; Y ), we can apply Lemma 5 to the equivalent systeṁ This shows the assertion.

Differentiability of the value function
In this section we perform a sensitivity analysis for the stabilization problem. The main purpose is to analyze the dependence of solutions to (P ) with respect to the initial condition y 0 and to show the differentiability of the associated value function, defined by J(y, u), subject to: e(y, u) = (0, y 0 ).

Existence of a solution and optimality conditions
In Lemma 8 we prove the existence of a solution (ȳ, u) to problem (P ), assuming that y 0 Y is sufficiently small. We derive then in Proposition 10 first-order necessary optimality conditions. Lemma 8. There exists δ 1 > 0 such that for all y 0 ∈ B Y (δ 1 ), problem (P ) possesses a solution (ȳ,ū). Moreover, there exists a constant M > 0 independent of y 0 such that Proof. Let us set, for the moment, where C is as in Lemma 4 and M K denotes the constant from Corollary 6. Applying this corollary (with f = 0), we obtain that for y 0 ∈ B Y (δ 1 ), there exists a control u ∈ L 2 (0, ∞; U ) with associated state y satisfying ). We can thus consider a minimizing sequence (y n , u n ) n∈N with J(y n , u n ) ≤ M 2 δ 2 1 (1 + α). We therefore have for all n ∈ N that Possibly after further reduction of δ 1 , we eventually obtain that where M λ is as in Corollary 7. It then follows that the sequence (y n ) n∈N is bounded in W ∞ with sup n∈N y n W∞ ≤ 2M λ y 0 Y . Extracting if necessary a subsequence, there exists (ȳ,ū) ∈ W ∞ × L 2 (0, ∞; U ) such that (y n , u n ) (ȳ,ū) ∈ W ∞ × L 2 (0, ∞; U ), and (ȳ,ū) satisfies (23).
Let us prove that (ȳ,ū) is feasible and optimal. For any T > 0 let us consider an arbitrary Sinceẏ n ẏ in L 2 (0, T ; V ), we can pass to the limit in the l.h.s. of the above equality. Moreover, since Ay n Aȳ ∈ L 2 (0, T ; V ), Analogously, we obtain that We also have By Lemma 2, it then follows that . Since V is compactly embedded in Y , we obtain that y n −ȳ L 2 (0,T ;Y ) −→ n→∞ 0 with the Aubin-Lions lemma. We can pass to the limit in (24) and obtain Density of H 1 (0, T ; V ) in L 2 (0, T ; V ) implies that e(ȳ,ū) = (0, y 0 ). Finally, by weak lower semicontinuity of norms it follows that J(ȳ,ū) ≤ lim inf n→∞ J(y n , u n ), which proves the optimality of (ȳ,ū).
For the derivation of the optimality system for (P ) we need the following technical lemma.
Then, for all f ∈ L 2 (0, ∞; V ) and y 0 ∈ Y , there exists a unique solution to the following system: Moreover, αū Moreover, there exists a constant M > 0, independent of (ȳ,ū), such that Proof. Let us set δ 2 = δ 1 . By Lemma 8, problem (P ) has a solution (ȳ,ū). In the first part of the proof, we derive abstract optimality conditions, by proving that the mapping e (used for formulating the constraints) has a surjective derivative. For proving the differentiability of e, we only need to consider the nonlinear term. We have F (y) = N (y, y) and we know that N is a bounded bilinear mapping from W ∞ × W ∞ to L 2 (0, ∞; V ), by Lemma 2. Thus N and F are Fréchet differentiable, and so is e, with Observe that by Corollary 3 By Lemma 8, it further holds that for some constant M independent of (r, s) and y 0 . It follows from the surjectivity of De(ȳ,ū) that there exists a unique pair (p, µ) In the second part of the proof, we derive the costate equation (25) and relation (26) from (29). As can be easily verified, J is differentiable with Taking z = 0 and letting v vary in L 2 (0, ∞; U ), we deduce from (30) and (31) that which proves (26). Let us set W 0 ∞ = {z ∈ W ∞ | z(0) = 0}. Taking now v = 0, we obtain that for all z ∈ W 0 ∞ , p,ż L 2 (0,∞;V ),L 2 (0,∞;V ) = p, Az − Gz L 2 (0,∞;V ),L 2 (0,∞;V ) + ȳ, z L 2 (0,∞;Y ) .
The time derivative of p, in the sense of distributions, can be extended to a linear form on W 0 ∞ with the formula: We immediately deduce from (32) and (23) that for all z ∈ W 0 ∞ , Since W 0 ∞ is dense in L 2 (0, ∞; V ) for the L 2 (0, ∞; V )-norm, we can extendṗ to a bounded linear form on L 2 (0, ∞; V ), i.e.ṗ can be extended to an element of L 2 (0, ∞; V ), moreover, the following bound holds true: It finally follows that p ∈ W ∞ and thus that the costate equation (25) is satisfied. It remains to bound p in L 2 (0, ∞; V ). Let r ∈ L 2 (0, ∞; V ) and let (z, v) satisfy De(ȳ,ū)(z, v) = (r, 0) and the bound (28) (with s = 0). Using the optimality condition (29), the expression (30) of DJ(ȳ,ū), estimate (28), and estimate (23) on (ȳ,ū), we obtain the following inequalities: Since r was arbitrary and since M does not depend on r, we obtain that p L 2 (0,∞;V ) ≤ M y 0 Y . Combining this estimate with (33), we finally obtain (27).

Sensitivity analysis
We define the space The well-posedness of Φ follows from the considerations on e(y, u) and the costate equation (25) that have been given in the proof of Proposition 10.
Proof. The cost function J is clearly infinitely differentiable. Since V(y 0 ) = J(Y(y 0 ), U(y 0 )), V is then the composition of infinitely differentiable mappings, which shows the assertion.

Derivatives of the value function
By standard arguments, we can derive a Hamilton-Jacobi-Bellman equation which provides an optimal feedback control based on the derivative of the value function.
All along the section, the first-order derivative DV(y 0 ) is either seen as a linear form on Y or is identified with its Riesz representative in Y . The identification is done for example in the term B * DV(y 0 ) 2 U appearing in the HJB equation below.
For deriving a Taylor series expansion of V, let us follow the approach from [2] and differentiate (38) in some direction z 1 ∈ D(A). To alleviate the calculations, we denote the variable y 0 in (38) by y. We then obtain A second differentiation in the directions (z 1 , z 2 ) ∈ D(A) 2 yields the equation Since V(0) = 0 and V(y) ≥ 0 for all y ∈ Y , it follows that DV(0) = 0. We can thus evaluate the last equation for y = 0 to obtain We recall that D 2 V(0) ∈ M(Y × Y, R) is a bounded and symmetric bilinear form on Y and thus can be represented (see, e.g., [26, Chapter 5, Section 2]) by an operator Π ∈ L(Y ) such that D 2 V(0)(y, z) = Πy, z Y , for all y, z ∈ Y.
As a consequence, we can formulate (40) as Equation (41) is the well-known algebraic operator Riccati equation which has been studied in detail in, e.g., [17,27]. From the stabilizability assumption A2, and the fact that the pair (A, id) is exponentially detectable as a consequence of (18), we conclude that (41) has a unique stabilizing solution Π ∈ L(Y ). In the discussion below, we denote by the closed-loop operator associated with the linearized stabilization problem. In particular, let us mention that A π generates an analytic exponentially stable semigroup e Aπt on Y . Hence, for trajectories of the formỹ = e A· y, y ∈ Y it follows thatỹ ∈ W ∞ . For higher order derivatives of V, we follow the exposition from [14]. For this purpose, let us briefly recall the symmetrization technique introduced there. Let i and j ∈ N, consider S i,j = σ ∈ S i+j | σ(1) < · · · < σ(i) and σ(i + 1) < · · · < σ(i + j) , where S i+j is the set of permutations of {1, . . . , i + j}. A permutation σ ∈ S i,j is uniquely defined by the subset {σ(1), . . . , σ(i)}, therefore, the cardinality of S i,j is equal to the number of subsets of cardinality i of {1, . . . , i + j}, that is to say |S i,j | = i+j i . For a multilinear mapping T of order i + j, we set The following proposition is a generalization of the Leibniz formula for the differentiation of the product of two functions.
Proposition 16. Let Z be a Hilbert space. Let f : Y → Z and g : Y → Z be two k-times continuously differentiable functions. Then, for all k ≥ 1, for all y ∈ Y and (z 1 , . . . , Proof. The proof is analogous to the one given in [14,Lemma 10] for Z = R. where the multilinear form R k : D(A) k → R is given by Proof. The proof relies on successive differentiations of (38). For a bilinear control problem, a similar result has been obtained in [14,Theorem 12]. In particular, it was shown that Obviously, for k ≥ 3, we have D k ( 1 2 y 2 Y ) = 0. Let us discuss the structure of the derivatives of the remaining terms appearing in (38). Applying Proposition 16 to the term B * DV(y) 2 U , we obtain Since V has a minimum at the origin, we have DV(0) = 0 and the terms for i = 0 and i = k vanish when evaluated in y = 0. By definition of the Sym-operator, for i = 1 we obtain As explained previously, we can represent D 2 V(0) in terms of the solution Π of the algebraic operator Riccati equation. This shows A similar relation can be derived for i = k − 1. Finally we consider the term D k (DV(y)F (y)). By Proposition 16, we get Since D 3+ F (y) = 0 for all ≥ 0, the previous equation simplifies as follows Evaluating the last expression in y = 0 yields since F (0) and DF (0) are both null. Combining (44), (45) and (46) proves the assertion.

Estimates for the velocity
In this section we analyze the polynomial feedback law u d derived from the Taylor series approximation of the value function ·, y, . . . , y).
The associated closed-loop system is given bẏ The open-loop control generated by u d (that is, the mapping t ≥ 0 → u d (y d (t)) ∈ U ) will also be denoted u d , without risk of confusion.
We begin with some local Lipschitz continuity estimates for the nonlinear part of the feedback law. For this purpose, we set for all k ≥ 3. The closed-loop system can be reformulated as follows: Lemma 18. For all k ≥ 3, there exists a constant C(k) > 0 such that for all y and z ∈ Y , Moreover, for all δ ∈ [0, 1], for all y and z ∈ W ∞ such that y W∞ ≤ δ and z W∞ ≤ δ, Proof. We have the identity The first inequality easily follows, with We also obtain that for all y and z ∈ W ∞ , The second inequality follows, since k ≥ 3 and δ ≤ 1.
The well-posedness of the closed-loop system can be now established with the same tools as those used in Lemma 5.
Proof. The existence of a solution y ∈ W ∞ , satisfying (50), can be obtained exactly as in Lemma 5. Thus we only discuss uniqueness. Let y and z denote two solutions to (47) in W ∞ . Let us set e = y − z. Arguing as in the proof of Lemma 5, one can prove the existence of M > 0 such that for all t ≥ 0. Since y and z ∈ W ∞ and e(0) = 0, we obtain with Gronwall's inequality that e = 0, which proves the uniqueness of the solution to the closed-loop system.
Theorem 20. Let d ≥ 2. There exist δ 6 > 0 and M > 0 such that for all y 0 ∈ B Y (δ 6 ), it holds that where (ȳ,ū) = (Y(y 0 ), U(y 0 )), y d is the solution of the closed-loop system (47) with initial condition y 0 , and u d is the generated open-loop control.
Proof. Let us fix δ 6 = min δ 5 , (4(C + d k=3 C(k))M 2 cls ) −1 , so that Proposition 14 and Theorem 19 apply for y 0 ∈ B Y (δ 6 ). By Taylor's theorem, see, e.g., [35,Theorem 4A], there exists δ > 0 such that for all y ∈ B Y (δ), where the remainder term R d satisfies for some constant M independent of y. Reducing if necessary δ 6 , we have that ȳ(t) Y ≤ δ for all t ≥ 0. Combining then (39) and the Taylor expansion (51), we obtain thaṫ Let us now consider the error dynamics e :=ȳ − y p . We have e(0) = 0, moreover by (49) and (52), Alternatively, e can be expressed as the solution of the systeṁ where the source term f is given by Considerδ ∈ (0, 1]. The precise value ofδ will be fixed later. By Lemma 11 and Theorem 19, we can reduce δ 6 so that max( ȳ W∞ , y d W∞ ) ≤δ. We first observe that Applying further Lemma 4 and Lemma 18, we obtain . For the solution of system (53) we thus obtain the estimate The constant M > 0 in the above estimate is independent ofδ. We can now defineδ = min 1, 1 2M . The first estimate on ȳ − y d W∞ follows.
Let us estimateū − u d . By (39) and by definition of the generated open-loop control u d , we have thatū Let us estimate the two terms of the right-hand side. It is easy to check that Using the techniques of Lemma 18 and the estimate on ȳ − y d W∞ , we also obtain that The second estimate onū − u d follows.
We have the following result, extending Theorem 20.
wherep and p d denote the pressure terms associated with (ȳ,ū) and (y d , u d ) respectively.
Proof. We have introduced in the proof of Lemma 21 the term g associated with a feasible pair (y, u). Let us denote byḡ and g d the corresponding terms associated with (ȳ,ū) and (y d , u d ). One can verify that as a consequence of Theorem 20, ḡ − g d L 2 (0,∞;H −1 (Ω)) ≤ M y 0 d Y . Proposition 22 follows then with similar calculations to those performed in the proof of Lemma 21.

A numerical example
In this section, we present numerical simulations for the two-dimensional Navier-Stokes equations and computed feedback laws of order 2 and 3. The discretization procedure and the example setups are classical and are taken from [6]. The main purpose is to show that the computation of higher order feedback laws is possible and, depending on the chosen parameters, visible differences to a Riccati-based feedback law can be observed.

Setup and discretization
We briefly summarize the numerical implementation provided in [6]. Therein a Taylor-Hood P 2 -P 1 finite element discretization for a two dimensional wake behind a cylinder is discussed. The computational domain Ω = (0, 2.2) × (0, 0.41) as well as a non uniform grid are shown in Figure  1. For all simulations, we use the Reynolds number Re := 1 ν = 90 and the parabolic inflow profile discussed in [6]. For the upper and lower end of the geometry, no slip boundary conditions are employed. The outflow is modeled by do nothing boundary conditions on the right end of the geometry. For the desired stabilization, we utilize a distributed, separable control acting in the control domain Ω c := [0.27, 0.32] × [0. 15, 0.25]. In particular, the control operator is of the form where the control shape functions w 1 , w 2 and w 3 are piecewise linear functions which are constant along the x 1 -direction. The finite element discretization is computed in FEniCS and the resulting matrices associated with the spatial semidiscretization are exported to MATLAB. As described in detail in [6], the (spatially) discrete system takes the form where E, K ∈ R nv×nv are the mass and stiffness matrices, G T ∈ R np×nv represents the discrete divergence operator, the tensor matricization H ∈ R nv×n 2 v represents the trilinear form (9) and B ∈ R nv×6 is the discrete control operator. Note that H can be constructed in such a way that H(z 1 ⊗z 2 ) = H(z 2 ⊗z 1 ) for any z 1 , z 2 ∈ R nv . The time invariant vectors f z ∈ R nv and f q ∈ R np are due to the elimination of the boundary nodes. The following results correspond to a discretization level with n v = 9356 and n p = 1289. The velocity profile of the unstable steady state solutionz shown in Figure 2 is obtained by a Picard iteration applied to the uncontrolled stationary system, i.e., system (58) withż(t) = 0 and u(t) = 0. The turbulent velocity profile to be stabilized is the result of a random perturbation of the form z(0) =z + z 2 2000 · randn(n v , 1).

Reformulation as an ODE system
System (58) is a system of differential-algebraic equations (DAEs) and hence the results from above are not readily applicable. While a thorough analysis in the framework of control of DAEs is certainly of interest, at this point we employ a reformulation initially proposed in [24] that allows to rewrite the dynamics as a set of ODEs for the velocity vector z. As in (3), we consider the shifted variables y = z −z and p = q −q, respectively. Consequently, we obtain Eẏ(t) = Ay(t) + H(y(t) ⊗ y(t)) + Bu(t) + Gp(t), where A = −K + H(z ⊗ I + I ⊗z). Let us note that the second equation implies G Tẏ (t) = 0. Following [24,Section 3], from the first equation, we thus obtain We can now eliminate the pressure from (59) using the relation ⊗ y(t)) + Bu(t)) .
In fact, as has been discussed in [8], the matrix P = P 2 as a discrete realization of the Leray projector. Since G T y = 0, we have P T y(t) = y(t) so that we can multiply the last equation by P to obtain (P EP T )ẏ(t) = (P AP T )y(t) + P HP T ⊗ P T (y(t) ⊗ y(t)) + (P B)u(t).
Finally, by means of a decomposition P = Θ Θ T r with Θ T Θ r = I we can project onto the n v − n p dimensional subspace range(P ) and arrive at the ODE system whereỹ = Θ T y(t). For the initialization, we useỹ(0) = Θ T y 0 . At this point, we emphasize that the explicit formulas yield dense matrices and thus are rather a theoretical tool. In particular, an explicit computation of H is infeasible for the problem dimension considered here. As a remedy, we work with an implementation that applies the above operations whenever a matrix vector multiplication is needed.

Computing the feedback gain
With the previous considerations in mind, we focus on the stabilization problem inf u∈L 2 (0,∞;R 6 ) J(ỹ 0 , u), subject to: e(ỹ u , u) = (0,ỹ 0 ) where We illustrate the effect of higher order feedback laws by computing the first two non trivial derivatives D 2 V(0) and D 3 V(0), respectively. For the computation of D 2 V(0) ≡ Π ∈ R (nv−np)×(nv−np) , we have to solve the algebraic matrix Riccati equation which in our case was done by means of the MATLAB function care. For the third order tensor D 3 V(0) ≡ X ∈ R (nv−np) 3 we have to solve a linear system of the form where π = vec(Π) denotes the vectorization of Π and the permutation matrix P is given by Let us emphasize that F is the discrete realization of the term R 3 in (43). In particular, the tensor F is symmetric. Note that computing a solution X to A T X = F is infeasible without using further tools such as model order reduction or tensor calculus as storing the vector X ∈ R (nv−np) 3 already requires more than 4 TB of data. As a remedy, we aim for a direct computation of the corresponding feedback gain without explicitly computing X . With this in mind, we proceed as in [13] and utilize a quadraturebased approximation that has been analyzed in [22]. From [22,Lemma 3], it follows that As shown in [22,Theorem 9], the previous integral can be well approximated by a tensor sum of the form where t j and w j are suitable quadrature points and weights and λ denotes a constant determined by the spectrum of the matrix pencil ( E, A). Combining the representation in (62), (63) and (64), we obtain the following approximation formula for the feedback gain with r = 30 in the numerical examples. By use of algebraic manipulations such as reshaping and transposition of matrices, the computation of the permutation matrix P as well as computation of the dense matricization H can be avoided. As a consequence, we obtain an approximation of K ∈ R 6(nv−np) 2 whose storage requires less than 4 GB of data. Let us point out that the above considerations do not fully break the curse of dimensionality but nevertheless allow us to compute a third order feedback law even for a spatially discretized PDE. For the simulation of the timevarying systems, we make use of the MATLAB function ode23 with the standard relative error tolerance 10 −3 . In each time step, the control law u 3 (ỹ) is obtained via where I 6 denotes the identity matrix for the control space R 6 .

Results
Below, we present a numerical comparison for two different values of α. In Figure 3, the control laws corresponding to (61) with α = 1 are shown. We observe that both feedback laws u 2 and u 3 , respectively, exhibit a similar behavior and create vortices which induce the desired control. Indeed, the control velocities in  to investigate the numerical convergence behavior as the order of the control laws increases. At the moment, however this is out of reach, and could be based on model reduction techniques in an independent numerical endeavor. In Figure 4, we observe that the amplitudes of the u 3 controls decay more rapidly than those of the u 2 controls. This is consistent with Figure 5, where we compare the dynamical behavior of u 2 2 2 and u 3 2 2 . Let us emphasize that for α = 10 −4 , for all t, the norm of the control law u 3 (t) is smaller than the one of u 2 (t). For the values of the cost functionals, we obtain J(ỹ u2 , u 2 ) = 0.9546, J(ỹ u3 , u 3 ) = 0.8432, for α = 1, J(ỹ u2 , u 2 ) = 0.0128, J(ỹ u3 , u 3 ) = 0.0125, for α = 10 −4 , which indicates that higher order feedback laws can be of interest for feedback stabilization.

Outlook
In the present paper we demonstrated that the approach that we carried out for obtaining Taylor approximations to the value function of optimal control problems related to the Fokker-Planck equation, is also applicable for optimal control of the Navier-Stokes equations in dimension two. The question arises to which extent analogous results can be obtained for dimension three and for boundary control problems. In dimension three the situation will be significantly different from that of the current paper. It will not be possible to work with weak variational solutions. Rather one has to resort to strong variational solutions, and thus one can expect at best that the value function is smooth on V rather than on Y . This leads to difficulties for the operator representations of the derivatives of the value function. Alternatively one can start by analyzing (43) as equations for abstract multilinear forms D k V(0), which are not necessarily obtained of derivatives of V. This is an approach which we plan to follow.