Keywords: modèles réduits, échantillonage optimal, contrôle

Modèle d'ordre reduit : échantillonage optimal et strategies de contrôle lineaire Résumé : Nous proposons une méthode d'échantillonnage optimale pour construire un modèle d'ordre réduit construit avec la Décomposition Orthogonales aux valeurs Propre (POD) qui soit robuste en fonction des variation des paramètres d'entrée. Cette méthode a été appliquée sur un cas test qui est l'écoulement confiné autour d'un cylindre de section carré lorsque le nombre de Reynolds varie. Par ailleurs nous examinons également le lien entre les modes instables et la modélisation POD. Un exemple de contrôle basé sur une approche linearisée est présenté.

Introduction

In fluid mechanics one of the most popular method to get a reduced-order model is the Proper Orthogonal Decomposition (POD) originally introduced by Lumley in 1967 [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF] in turbulence context. The main drawback for flow control is that the POD basis is not optimal to represent a flow generated with different system parameters with respect to those used to build the basis. To get rid of this problem, different strategies can be employed. The first one is to update the POD basis during the optimization. One method is to use trust region method (TRPOD see [START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by trust-region methods and pod reduced-order models[END_REF]), another is to calibrate over several dynamical cases [START_REF] Weller | Feedback flow control design by low-order modeling for laminar wake stabilization[END_REF]. Yet another method is to build a robust POD basis that can be used all along the optimization process. This kind of POD basis can be generated using chirp excitation [START_REF] Bergmann | Optimal rotary control of the cylinder wake using proper orthogonal decomposition reduced-order model[END_REF] or using an appropriate sampling of the input parameter space.

In this spirit, the first part of this study is devoted to the construction of a robust model that can be used for control without updating of the POD basis. The idea is to sample in an efficient way the input parameter subspace. Two classes of sampling methods are commonly used: the one shot method and the iterative one. In the one shot method the sampling is obtained by partitioning the range of variation of the input parameter space. The partitions can be found using different strategies as, for instance, the uniform distribution, the orthogonal sampling, the Sobol algorithm etc... An alternative strategy to the classical partition strategies is the Centroidal Voronoi Tessellations (CVT, see [START_REF] Du | Centroidal voronoi tessellations: Applications and algorithms[END_REF][START_REF] Burkardt | Centroidal Voronoi Tessellation-Based Reduced-Order Modeling of Complex Systems[END_REF][START_REF] Burkardt | Centroidal Voronoi Tessellation-Based Reduced-Order Modeling of Complex Systems[END_REF]) which leads to an efficient partition. This kind of tessellations can be efficiently computed using the Lloyd algorithm [START_REF] Du | Convergence of the lloyd algorithm for computing centroidal voronoi tessellations[END_REF]. The main drawback of the one shot strategy is that the number of sampling points has to be fixed a-priori and, in the case of CVT, the final configuration is strongly dependent on the initial condition. Thus, an a-priori analysis of the density function used to compute the centroidal tessellation is necessary to determine the proper refinement when sampling the range of variation of the input parameter. The other class of methods consist to add sampling points in an iterative way. Thus, we can choose the degree of accuracy by fixing a stopping criterion. One efficient iterative method is based on Greedy sampling (see [START_REF] Bui-Thanh | Model reduction for large-scale systems with high-dimensional parametric input space[END_REF]). In Greedy sampling, the new value of the input parameter to sample is chosen on the maximum of the density function, i.e. where the error or the residual given by the POD basis is larger. In this study we propose a new approach that couples Constrained CVT and Greedy methods.

In the second part of this study the capability of a linearized low order model of the flow is assessed in a control optimisation. A linear model of the flow can be used when the target solution is a steady state. In this case the objective of the controller consists in stabilizing a steady state of the system. Thus, small oscillations of the system around this target state are well represented by a linear model. At the same time, designing the controller using a linear model involves standard techniques and is simpler than using a non-linear model. Moreover, it is also interesting to explore the capabilities of reduced-order models in estimating unstable modes in the linear stability analysis of a flow since this aspect is typically very demanding in terms of computational costs. Indeed, this analysis requires codes simulating the linearized flow equations and, possibly, generating the matrix of the linearized system, which is not always possible when working with complex tools as those typically used in engineering applications. For this reason, the starting point of the present analysis is just the availability of a nonlinear code for simulating the Navier-Stokes equations which cannot be linearized. The reduced order model of the linearized flow equations is built using only this tool.

On the other hand, the use of a non-linear reduced-order model for flow control, although more expensive and complex, allows more general control strategies (i.e. minimization of general cost functions, different control targets etc...). In [START_REF] Weller | Feedback flow control design by low-order modeling for laminar wake stabilization[END_REF] a control strategy based on a non-linear model is reported. In that reference it is also shown that the proposed strategy, when used for the particular objective of stabilizing a steady state for the system, has a clear behavior in terms of the spectrum of the linearized Navier-Stokes operator around the target flow.

Flow configuration

In this study the two dimensional confined square cylinder wake flow (figure 1(a)) is chosen as a prototype of separated flow. The Navier-Stokes equations write:

∂ u u u ∂t + (u u u • ∇ ∇ ∇)u u u = -∇ ∇ ∇p + 1 Re ∆u u u, ( 1a
) ∇ ∇ ∇ • u u u = 0, (1b)
where u u u = (u, v) T and p denote respectively the velocity and pressure fields. The parameter Re = U ∞ L/ν denotes the Reynolds number, with U ∞ = u(0, H/2) the maximal inflow velocity of the incoming Poiseuille flow, L the length of the side of the square cylinder and ν the kinematic viscosity. Otherwise, we used the same numerical methods and parameters as those introduced in [START_REF] Galletti | Low-order modelling of laminar flow regimes past a confined square cylinder[END_REF], i.e. the blockage ratio β = L/H is equal to 1/8 and the domain Ω is (-10L, 22L) × (-4L, 4L). For control purposes we placed two jets in opposite phase on the upper and lower faces of the cylinder, as shown in figure 1(b). Following the modeling of the actuators in [START_REF] Weller | Robust model identification of actuated vortex wakes[END_REF][START_REF] Weller | Feedback flow control design by low-order modeling for laminar wake stabilization[END_REF] the boundary conditions on the surface of the cylinder on the jets areas Γ c are modified:

v(x x x,t) = c(t), x x x ∈ Γ c
As it was already pointed out, we will present two different reduced order modeling (ROM) based on Proper Orthogonal Decomposition (POD).

INRIA 2 Reynolds dependent pressure extended reduced order model based on proper orthogonal decomposition

Proper Orthogonal Decomposition

The Proper Orthogonal Decomposition (POD) was first introduced in turbulence by Lumley [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF] in 1967 as an unbiased definition of the coherent structures widely known to exist in a turbulent flow. A comprehensive review of the POD can be found in Refs [START_REF] Aubry | The dynamics of coherent structures in the wall region of a turbulent boundary layer[END_REF][START_REF] Sirovich | Turbulence and the dynamics of coherent structures. Parts I,II and III[END_REF][START_REF] Cordier | Lecture series 2002-04 on post-processing of experimental and numerical data[END_REF]. The POD, also known as Karhunen-Loève decomposition, principal component analysis or empirical eigenfunctions method, consists of looking for the deterministic function Φ Φ Φ(x x x) that is most similar in an average sense to the realizations U U U(x x x,t). For instance, the realizations U U U(x x x,t) can be velocity fields, pressure fields, temperature fields, etc. Since in this study the data are issued from numerical simulations, the method to compute POD modes introduced by Sirovich [19] is adopted (see [START_REF] Cordier | Lecture series 2002-04 on post-processing of experimental and numerical data[END_REF] for justifications). In this case, the constrained optimization problem reduces to the following Fredholm integral eigenvalue problem:

T 0 C(t,t ′ )a n (t ′ ) dt ′ = λ n a n (t) (2)
where the temporal correlation tensor C(t,t ′ ) is defined by:

C(t,t ′ ) = 1 T U U U(x x x,t),U U U(x x x,t ′ ) Ω . (3)
The inner product (., .) Ω between two fields U U U and V V V is computed as:

(U U U,V V V ) Ω = Ω U U U •V V V dx x x = Ω n c ∑ i=1 U i V i dx x x,
where U i represents the i th component of the vector U U U with dimension n c . The eigenvalues λ n (n = 1, 2, . . . ) determined in (2) are all real and positive and form a decreasing and convergent series. Each eigenvalue represents the contribution of the corresponding mode Φ Φ Φ n to the information content of the original data. Note that if U U U are the velocity fields, the information content reduces to the kinetic energy.

In Eq. ( 2), a n are the time-dependent POD eigenfunctions of order n. These modes form an orthogonal set, satisfying the condition:

1 T T 0 a n (t)a m (t) dt = λ n δ nm . (4)
The associated eigenvectors Φ Φ Φ n (also called empirical eigenfunctions) form a complete orthogonal set and are normalized, so that they verify

(Φ Φ Φ n , Φ Φ Φ m ) Ω = δ nm .
The spatial basis functions Φ i n can then be calculated from the realizations U i and the coefficients a n with:

Φ i n (x x x) = 1 T λ n T 0 U i (x x x,t)a n (t) dt. ( 5
)
Since the POD eigenfunctions can be represented as linear combinations of the realizations, they inherit all the properties of the original data. For instance, the eigenfunctions are divergence free for incompressible flows. Moreover, the eigenfunctions verify the boundary conditions of the numerical simulation used to determine the flow realizations.

The set of POD modes {Φ Φ Φ n } N r n=1 is complete in the sense that any realization U U U(x x x,t) contained in the original data set, can be expanded with arbitrary accuracy (as a function of N r ≥ 1) in the eigenfunctions as

U U U(x x x,t) ≃ U U U [1,•••,N r ] (x x x,t) = N r ∑ n=1 a n (t)Φ Φ Φ n (x x x). ( 6
)
Hereafter, we consider that the ensemble used to determine the POD modes consists of N t flow realizations (called time snapshots)

U U U(x x x,t i ), x x x ∈ Ω, taken at t i ∈ [0, T ] , i = 1, • • • , N t .
The energetic optimality of the POD basis functions suggests that only a very small number of POD modes may be necessary to describe efficiently any flow realization of the input data i.e. N r ≪ N s . In practice, N r is usually determined as the smallest integer M such that the Relative Information Content,

RIC(M) = ∑ M i=1 λ i / ∑ N s
i=1 λ i , is greater than a predefined percentage of energy, δ .

Reynolds adaptive pressure extended reduced order model

In many practical applications for incompressible flows, reduced order basis are usually computed form the velocity fields. We thus have U U U(x x x, t) ≡ u u u(x x x, t), and decomposition (6) becomes:

u u u(x x x,t) ≃ N r ∑ n=1 a n (t)φ φ φ n (x x x), ( 7
)
where φ φ φ n denote the velocity POD basis functions. A low dimensional dynamical system is obtained via a Galerkin projection of the Navier-Stokes equations [START_REF] Aubry | The dynamics of coherent structures in the wall region of a turbulent boundary layer[END_REF]. The Galerkin projection formally writes:

φ φ φ i , ∂ u u u ∂t + (u u u • • • ∇ ∇ ∇)u u u Ω = -(φ φ φ i , ∇ ∇ ∇p) Ω + φ φ φ i , 1 Re ∆u u u Ω . (8)
All the term in system (8) can be evalueted from velocity POD basis functions, except the boundary pressure term (φ φ φ i , ∇ ∇ ∇p) δ Ω . It is demonstrated that the contribution of the pressure term vanishes in many wall bounded flows (since homogeneous boundary condition for the POD basis functions can be obtained after some algebra manipulations). However, Noack [START_REF] Noack | The need for a pressure-term representation in empirical galerkin models of incompressible shear flows[END_REF] proved that neglecting the pressure term for convectively unstable shear layers (as the mixing layer or the wake flow) can lead to large amplitude errors in the Galerkin model. Following the idea in [START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by trust-region methods and pod reduced-order models[END_REF] the pressure term can be easily computed using p = p (see decomposition (9b)). An important key issue is that, knowing the pressure field, it is possible to evaluate the Navier-Stokes residuals that can be as an error estimator. This estimator can then be used to perform a robustness improvement procedure as described in §3.2. Thus, we use a global basis for both the velocity INRIA and pressure fields (see [START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by trust-region methods and pod reduced-order models[END_REF]). The exact flow fields u u u and p are then approximated by:

u u u(x x x, t) = N r ∑ i=1 a i (t)φ φ φ i (x x x) (9a) p(x x x, t) = N r ∑ i=1 a i (t)ψ i (x x x). ( 9b
)
The velocity and the pressure basis functions, φ φ φ i and ψ ψ ψ i respectively, are determined from the temporal correlation tensor (3) computed using U U U(x x x, t) = (u u u(x x x, t), p(x x x, t)) T . The basis functions φ φ φ i and ψ ψ ψ i are determined as Φ Φ Φ(x x x, t) = (φ φ φ (x x x, t), ψ(x x x, t)) T , Φ Φ Φ(x x x, t) being obtained from [START_REF] Bui-Thanh | Model reduction for large-scale systems with high-dimensional parametric input space[END_REF].

The substitution of equations ( 9) in the Navier-Stokes momentum equations (1a) and a Galerkin projection lead to the Reduced Order Model:

N r ∑ j=1 L i j da j dt = N r ∑ j=1 B i j a j + N r ∑ j=1 N r ∑ k=1 C i jk a j a k , (10)
with initial conditions

a i (0) = (U U U(x x x, 0), Φ Φ Φ i (x x x)) Ω i = 1, • • • , N r , (11)
where the coefficients1 L i j , B i j and C i jk are given by:

L i j = + φ φ φ i , φ φ φ j Ω , (12a)
B i j = -φ φ φ i , 1 Re ∆φ φ φ j -∇ ∇ ∇ψ j Ω , (12b)
C i jk = -φ φ φ i , φ φ φ j • ∇ ∇ ∇ φ φ φ k Ω . (12c)
In order to build a Reynolds adaptive Reduced Order Model we split the viscous terms in the model from B i j ; this leads to the final model, which has the classical form with extra terms:

N r ∑ j=1 L i j da j dt = N r ∑ j=1 B Re i j a j + N r ∑ j=1 B p i j a j + N r ∑ j=1 N r ∑ k=1 C i jk a j a k , (13)
where :

B Re i j = -φ φ φ i , 1 Re ∆φ φ φ j Ω , (14a)
B p i j = + (φ φ φ i , ∇ ∇ ∇ψ j ) Ω . (14b)
Every dynamic associated with a Reynolds number belonging in a predefined interval can be approximated with more or less efficiency using model [START_REF] Galletti | Accurate model reduction of transient and forced flows[END_REF]. For simplicity reasons, we considered only the Reynolds number as the system input parameter, but all the concepts introduced here for the Reynolds number can be easily extended to other system parameters, as for instance for high dimensional control space as done in Ref. [START_REF] Bui-Thanh | Model reduction for large-scale systems with high-dimensional parametric input space[END_REF]. As discussed in several papers [START_REF] Galletti | Accurate model reduction of transient and forced flows[END_REF][START_REF] Bergmann | Optimal rotary control of the cylinder wake using proper orthogonal decomposition reduced-order model[END_REF][START_REF] Couplet | Calibrated reduced-order pod-galerkin system for fluid flow modelling[END_REF], the initial value problem (13) can be inaccurate or even unstable, since it may not take into account enough information about the flow dynamics. Indeed, although a number N r of modes can be sufficient to capture most of the flow energy, the neglected modes continue to play an important role in the flow dynamics through their interaction with the resolved ones. In order to avoid this issue we perform a calibration procedure of the low order model.

Calibration procedure

In order to build a robust order model we applied the calibration technique described in [START_REF] Weller | Robust model identification of actuated vortex wakes[END_REF] and resumed in this section for a Reynolds dependent reduced order model. Following that approach, the system coefficients are adjusted to fit the solution of system (13) to several dynamics. This avoids the effects of the neglected modes and it allows to build a stable model that reproduces accurately the actual dynamics used to calculate the POD basis.

Setting:

X i = {L -1 i j } j=1•••N r • {B Re i j } j=1•••N r , {L -1 i j } j=1•••N r • {B p i j } j=1•••N r , {L -1 i j } j=1•••N r • {C i jk } j,k=1•••N r T and f f f (a a a(t), Re) = {a j (t)} j=1•••N r , {a j (t)a k (t)} j,k=1•••N r , 1 Re
the equation in ( 13) can be written in the compact form:

ȧi (t) = f f f (a a a(t), Re) • X i
We consider a data base that includes a simulation obtained with Reynolds number Re = Re to calculate the POD basis. The system coefficients obtained by Galerkin projection of the NSE on the low-dimensional subspace are denoted X. The calibration procedure then consists in choosing X as the solution of:

min X N r ∑ i=1 T 0 ȧi (t) -f f f ( â a a(t), Re) • X i 2 dt + α N r ∑ r=1 X i -X i 2 ( 15
)
where

âi (t) = U U U(•,t), Φ Φ Φ r
and where α is the Tikhonov regularization parameter and is chosen ≪ 1. In this calibration procedure only the terms of B Re i j are calibrated. This is due to the assumption INRIA that the errors in the Galerkin model are due mainly to the fact that it neglects the small scales and therefore a large part of the viscous effects. B Re i j indeed results from the projection of the viscous term of the Navier-Stokes equations. The same procedure can be performed if the POD basis is calculated over N dynamics with different Reynolds numbers Re 1 , . . . , Re N . The new stable model is then calculated to fit the dynamics at various Reynolds numbers. Thus, the system (15) to solve becomes:

min X N r ∑ i=1 N ∑ ℓ=1 T 0 ȧℓ i (t) -f f f ( â a a ℓ (t), Re ℓ ) • X i 2 dt + α N r ∑ r=1 X i -X i 2 ( 16
)
where

âℓ i (t) = U U U ℓ (•,t), Φ Φ Φ r with U U U ℓ (•,t
) the snapshots at the instant t calculated with Re = Re ℓ .

In section ( §3.4) an example of the improvement on the behaviour of the reduced order model given by the calibration procedure is summarized.

Improvement of the POD ROM robustness

The aim of this section is to improve the robustness of reduced order models based on proper orthogonal decomposition (POD ROM). As we have just seen ( §2), two steps are involved to build a POD ROM: the first one is to compute a POD basis ( §2.1), and the second one is to derive a reduced order model by projecting the governing equations onto the subspace generated by the POD basis. The POD basis gives an optimal representation, in terms of kinetic energy, of the database of snapshots used to build the basis itself and generated by the system. However, when the input parameters vary, the basis becomes inaccurate, as it is the case in control problems (see [START_REF] Prabhu | The influence of control on Proper Orthogonal Decomposition of wall-bounded turbulent flows[END_REF][START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by trust-region methods and pod reduced-order models[END_REF]). The focus of this section is to improve the representation capabilities of a POD basis of a given flow when the Reynolds number (input parameter of the system) varies in a given range, so as to provide a single ROM that is efficient for the considered range.

As already stated, all the concepts introduced in this study can be easily extended to other system parameters, as, for instance, to a set (even large) of control parameters. The Reynolds number space (here, only a discrete interval) under consideration is denoted I = [Re L , Re R ], where we chose Re R = 180 and Re L = 40 or Re L = 70, depending on the considered case. Reynolds numbers Re L = 70 and Re R = 180 correspond approximately to the lower and higher bound for the 2D periodic regime for the flow around a confined square cylinder. The case with Re L = 40 is considered to investigate if the reduced order model is robust enough to predict a (Hopf) bifurcation of the system (that occurs at Re c ≈ 65 in this case). Numerically, I is discretized with ∆Re = 5, and it is denoted as I h .

In order to improve the functional subspace, we want to enrich the database in an one-shot way by adding some sets of snapshots at different Reynolds numbers Re i ∈ I . Let U [Re 1 ,...,Re N ] be the database composed by N sets of snapshots taken independently at Re 1 , . . . , Re N , where N is a free parameter depending on the desired accuracy of the POD basis. The projection of the numerical solution of the Navier-Stokes equations U U U(x x x, t) onto the N r retained POD modes is:

U U U [N r ] (x x x, t) = N r ∑ n=1 a n (t)Φ Φ Φ n (x x x) (17)
In the following, we will always use N r = 31. The number of basis functions is arbitrarily chosen enough large because it will be kept all along this study, even when they are computed using a database collected from N > 1 different Reynolds numbers Re i . Since we will always use N r = 31 we simply note that

U U U(x x x, t) ≡ U U U [N r ] (x x x, t) = N r ∑ n=1 a n (t)Φ Φ Φ n (x x x). ( 18
)
The temporal coefficients a n (t) can be evaluated in two ways:

• by projecting the numerical solution of the NSE onto the functional subspace:

a n (t) = Ω U U U(x x x, t)Φ Φ Φ n (x x x)dx x x, (19)
• by resolving a reduced order model.

In what follows we denoted U U U DNS and U U U ROM the fields computed with temporal coefficient obtained by projection and by model prediction respectively. Without loss of generality, U U U can be either U U U DNS or U U U ROM In order to test our criterion to improve the POD basis, described in §3.3, we will consider three initial bases: The three cases are summarized in Fig. 2. In this study we will always arbitrarily consider N t = 200.

Effect of the Reynolds number variations onto the projection error

In the following, the reconstruction capabilities of a given POD basis is estimated when the Reynolds number varies in the interval this is to compare, at each Re ∈ I h , the numerical solution U U U(x x x, t) of the Navier-Stokes equations to the POD reconstruction U U U(x x x, t) computed using a POD basis that corresponds to a given database U [Re 1 ,...,Re N ] . The numerical solution of the Navier-Stokes equations can be formally written as:

I = [Re L , Re R ]. A natural
U U U(x x x, t) = U U U(x x x, t) + U U U ′ (x x x, t), (20)
where U U U ′ (x x x, t) denotes the missing scales, i.e. the error made restricting the solution to the first N r basis functions

U U U ′ (x x x, t) = U U U(x x x, t) -U U U(x x x, t). ( 21
)
We defined the average of the L 2 norm over a temporal horizon T for missing scales [START_REF] Weller | Robust model identification of actuated vortex wakes[END_REF] by :

U U U ′ 2 = T U U U ′ (x x x, t) 2 dt. ( 22
)
Since we have to compute the numerical solution U U U(x x x, t) of the Navier-Stokes equation onto I h to achieve such comparison, the POD output flow fields U U U(x x x, t) can easily be computed with the projected coefficients a n (t) evaluated from [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. Parts I,II and III[END_REF]. This error indicates how the description capability of the POD basis changes due to variations of the Reynolds number (system parameter). In what follows, the temporal horizon T is taken to be equal to three vortex shedding periods (that depends on Re). Figure 3 shows the evolution of the error U U U ′ 2 versus the Reynolds number for the three initial databases described above. For all cases, we can see that the error is very small at Re i inside the POD database (in the order of 10 -7 ), and then it growths when the value of the Reynolds number moves away from Re i . This traduces the fact that the POD basis computed from a database collected from given dynamics is not able to give a good representation of flows that present other dynamics. We can see in figure 3 that model C seems to be more robust than models A and B. Indeed, the maximal error is smaller for model C. The reason is that the POD basis for case C is computed from 2 different dynamics and the other cases from only one. The aim is then to determine a sampling of K new sets of snapshots {Re i } N+K i=N+1 ∈ I h , to compute the most robust POD basis.

A residuals based error estimator

Since the evaluation of the error U U U ′ (x x x, t) involved the computation of the numerical solutions U U U(x x x, t) of the Navier-Sokes equations for each Re ∈ I h , the evaluation of criteria based on the error ( 22) are demanding from a computational viewpoint. It is then interesting to find an accurate estimation of the error (22). To this purpose, we introduced the average of the L 2 norm, over the same temporal horizon T , of the residuals of the Navier-Stokes operator R, evaluated using flow fields projected onto the POD basis U U U(x x x, t) :

R( U U U) 2 = T R( U U U(x x x, t)) 2 dt. ( 23
)
A comparison between the error U U U ′ 2 and the residuals R( U U U DNS ) 2 over I h is performed in figure 4. It is interesting to note that these two quantities show a similar behavior for all the considered test cases, especially after the Hopf bifurcation at Re ≈ 65. Indeed, the ratio

R( U U U DNS ) 2 / U U U ′
2 is approximately a constant over I h for all test cases. The residuals R( U U U DNS ) 2 is thus a good estimator of the error U U U ′ 2 . However using the residual R( U U U DNS ) 2 is not tractable in practice since it requires the computation of the numerical Navier-Stokes solution U U U. The idea is then to approximate the projection residuals R( U U U DNS ) with the prediction ones R( U U U ROM ).

A comparison between the projection and prediction residuals over I h is performed in figure 5 for the cases A, B and C. We recall that the models are calibrated on the initial set of dynamics and are then integrated using each Reynolds number Re ∈ I h . Projection and prediction residuals show a closed correlation for all the considered cases except for Reynolds number below the bifurcation. However, the predicted residuals are close to the error (compare figures 5 and 4). It is noticeable that the discontinuity in the residuals evolution marks the capability of the model to predict the exact position of the bifurcation at Re ≈ 65 without any knowledge about the critical Reynolds number of the dynamical bifurcation. Finally, the predicted residual R( U U U ROM ) 2 is a good estimator of the error U U U ′ 2 and can thus be used as a criterion to sample the input parameter space (here, I h ).

A residual based sampling method

As described in sec. 1, two classes of sampling methods are commonly used: the one shot and the iterative methods. In this study we propose an approach that couples the ideas of the two classes. In particular we will present a one shot method, a Constrained Centroidal Voronoi Tessellation, based on the Greedy ideas. The method proposed here is based on the residuals of the Navier-Stokes operator predicted using the flow fields predicted by the POD reduced order model. This allows to reduce the computational costs in respect of using the reconstruction error U U U ′ . We consider an initial database U [Re 1 ,...,Re N ] composed by N × N t snapshots collected at [Re 1 , . . . , Re N ]. Since we want a robust POD basis, we look for a sampling {Re i } M i=1 ∈ I M h such that the database U [Re 1 ,...,Re M ] produces models leading to reduction (or minimization in the optimal case) of the error evaluated over the whole subspace I h , where M has to be fixed as a function of the desired robustness. As already pointed out, the residuals of the Navier-Stokes equations R( U U U ROM ) can be easily calculated by integration of the calibrated ROM (13) for all Reynolds numbers in the discretized space Re ∈ I h . Thus, the density function used in the Constrained Centroidal Voronoi Tessellation [START_REF] Du | Constrained centroidal voronoi tessellations for surfaces[END_REF] is the Navier-Stokes residuals predicted by the POD model. We perform a Constrained Centroidal Voronoi Tessellation procedure starting from a random subspace Re z M 0 z=N+1 ∈ I M 0 -N , with M 0 > M. The initial Reynolds numbers [Re 1 , . . . , Re N ] are frozen while the new points are computed as being the centroids of the tessellation elements with respect to density function R( U U U) 2 . We exclude point k > N with INRIA the smaller average density function over the k th tessellation. This is done following Greedy method in order to refined where the density function reaches higher values. The size of the sampling is then M 1 = M 0 -1. This is an iterative process, and while M i > M we recompute a new Degenerated CCVT and exclude a new point k > N.

The final configuration M i = M is weakly dependent on the initial configuration for M 0 ≫ M. The Greedy Degenerated CCVT is summed up below, where the goal is to find a K-dimensional sampling to add at the N-dimensional initial sampling. 0. Random sampling with dimension K 0 > K. (the N first points are frozen).

At iteration i, start sampling process with dimension

M i = K i + N • Perform a Constrained Centroidal Voronoi Tessellation • if M i = M stop • if M i > M go to 2
2. Identify and exclude point k > N of the element with minimum integral

• M i+1 = M i -1. Increment i = i + 1, then go to 1
The sampling method presented above can be easily transposed for input parameter subspaces with dimension greater than one. The use of the residuals as error estimation leads to negligible computational costs, even for high dimensional input parameter spaces, as for instance active control space.

Sampling at different Reynolds numbers

The sampling technique described in the previous section was applied in order to improve the robustness of low order models for the three initial test cases described in sec. 3. The three models are calibrated over the dynamics of the databases as explained in sec. 2.3.

In order to display the effect of the calibration on the behaviour of the reduced order model, we calculate the residuals predicted both by (i) the calibrated and (ii) uncalibrated model for the case A. The two model are integrated with a variation of the input parameter Re over the whole interval of [Re L , Re R ]. The improvement of the dynamical behaviour is clearly visible if we calculate the average difference between the predicted residuals and the projected ones. In Tab. 1 the differences between the residuals predicted by the two models (calibrated and non calibrated) and the projected residuals are shown. Note that calibration leads to an accuracy gain of about 54% on the difference between the predicted and the projected residuals.

We notice that the effect of the calibration is progressively more important as the dimension of the space of the input parameters is increased.

In particular, as shown in [START_REF] Weller | Robust model identification of actuated vortex wakes[END_REF], when a space of parameters of an active control is considered and the pressure is not included in the basis, the calibrated model outperforms the uncalibrated one concerning robustness and accuracy. Table 1: Difference between the predicted residuals and the projected ones using Non Calibrated and Calibrated ROM for the test case A. ∆Res = (Res ROM -Res DNS ).

In order to increase the robustness we chose to add K = 2 new sampling points in Re. Starting with K 0 = 6 initial random Reynolds numbers, after four iterations the method gives the final sets of Reynolds numbers : The average error and the standard deviation evaluated over the whole subspace I are respectively defined by:

• Case A: I N f =
E = 1 Re R -Re L I U U U ′ (Re) 2 dRe, (24)
R = I ( U U U ′ (Re) 2 -E) 2 dRe. ( 25
)
While the error E measures the accuracy of the POD ROM, the standard deviation R measures its robustness.

Both projection U U U ′

2 prediction U U U ′ ROM 2 errors are plotted. The POD ROM prediction errors are close to the DNS projection ones. This proves the good behaviour of each calibrated reduced order model, showing that they are able to predict the system dynamics with a negligible error. The standard deviation has been evaluated for the POD models build using the sampling points found with both the Greedy Degenerated CCVT and the constraint uniform sampling CUS strategies.

For a given quantity F we define the percentage of the relative difference by ∆F = 100(F CUS -F CCV T )/F CCV T . The percentage of the relative difference errors and standard deviations are reported in table 2. Note that a positive difference suggests a smaller

Re

Re DMCVT (100) 55 160 -Re eq 100 70 140

U U U ′cus DNS 2 U U U ′cus ROM 2 U U U ′ccvt DNS 2 U U U ′ccvt

Re

Re DMCVT (120) 80 165 -Re eq 90 120 150 error or standard deviation for the Greedy Degenerated CCVT than for CUS. The average reconstruction errors given by ROM prediction are smaller when our sampling method is used. The reduced order model given by Greedy Degenerated CCVT are more robust in all the considered cases, even if in the third case the average error obtained by projection for the constraint uniform sampling is smaller. Indeed, for all the As noted above, in the third case the error of reconstruction computed by the POD model obtained by Degenerated CCVT is smaller than the one computed by the CUS. Finally the ROM gives a good behaviour in terms of reconstruction error also in the case C, starting from constrained points placed on the boundaries of the parameter space, and in presence of a bifurcation.

U U U ′cus DNS 2 U U U ′cus ROM 2 U U U ′ccvt DNS 2 U U U ′ccvt
U U U ′cus DNS 2 U U U ′cus ROM 2 U U U ′ccvt DNS 2 U U U ′ccvt ROM 2 (c) Test case C.
∆ ∆ ∆E E E pro j ∆ ∆ ∆E E E pro j ∆ ∆ ∆E E E pro j ∆ ∆ ∆E E EROM ∆ ∆ ∆E E EROM ∆ ∆ ∆E E EROM ∆ ∆ ∆R R Rpro j ∆ ∆ ∆R R Rpro j ∆ ∆ ∆R R Rpro j ∆ ∆ ∆R R RROM ∆ ∆ ∆R R RROM ∆ ∆ ∆R R RROM
Fig. 8 shows a comparison between projection and prediction residuals after CCVT sampling for all the considered test cases. The models are accurate in terms of residuals, in agreement with the assumption that the residuals can be used as error estimator. Again, in the case C, an error is visible. Indeed, the maximum on the predicted residuals reached across the bifurcation is larger than the maximum residual obtained by projection; this is due to the fact that the dynamics before Re c ≈ 65 are essentially steady solutions. For that, a small perturbation on the steady predicted coefficients given by the low order model is sufficient to give a larger error on the residuals estimation.

Finally, for all the considered cases, those cover an adequate variety of possible situations, the low order models obtained by the proposed sampling method are robust and accurate as in terms of reconstruction error than in terms of residuals estimation. Thus, in a sampling procedure, one can use the Degenerated Greedy CCVT to build robust parameter dependent reduced order model. This avoids huge computational costs by using residuals estimation of the calibrated ROM instead of the approximation error computed by projection. As explained above, this technique can be easily extended to high dimensional parameter space with negligible computational costs compared to those required by a procedure based on the reconstruction error evaluation.

Control based on a linearized model 4.1 POD-based model of the linearized Navier-Stokes equations with control

In this section the linearized reduced order model of the Navier-Stokes equations in presence of control actuation is described. The flow configuration is the same described in fig. 1. We consider a feedback proportional control actuated by the jets sketched in figure 1, using some measurements of vertical velocity given by N v sensors placed at x x x j in the cylinder wake on the centre line. The control law with feedback gains K j is :

c(t) = N v ∑ j=1 K j v(x x x j ,t) (26)
The aim is to find the set of feedback gains K j that stabilizes the vortex shedding in the cylinder wake.

The POD-based linear model is built using the snapshots obtained by a non-linear simulation of the transient flow dynamics, which is started from the steady unstable solution. The starting flow field, which is also the target flow of the controller, is found using the same code, by imposing the velocity field to be symmetric with respect to the symmetry line y = 0 and advancing the simulation in time until a steady state is reached. Indeed, in this particular case, the unstable vortex shedding mode is antisymmetric with respect to y = 0, and the symmetry constraint is a physically-based trick to suppress the instability and, thus, to find the steady unstable solution. More general strategies, which can be straightforwardly applied to an evolutive non-linear code (without the need of deriving linearizations as needed for Newton-like methods) are described in [START_REF] Åkervik | Steady solutions of the navier-stokes equations by selective frequency damping[END_REF], [START_REF] Galletti | Low-order modelling of laminar flow regimes past a confined square cylinder[END_REF].

Snapshots (N t is their number) obtained sampling a part of the transient dynamics, obtained with a particular control law c(t), are used to build a POD model. To this purpose, every snapshot u u u(x x x,t) is decomposed as follows:

w w w(x x x,t) = u u u(x x x,t) -u u u 0 0 0 (x x x) + c(t)u u u c c c (x x x), ( 27
)
where u u u 0 0 0 (x x x) is the unstable steady state and u u u c c c (x x x) is a flow field having a jet velocity equal to 1 and the velocity vanishing on all the other domain boundaries. This is obtained as proposed in [START_REF] Galletti | Accurate model reduction of transient and forced flows[END_REF], i.e. considering the time-averaged flow field (u u u ′ ′ ′ (x x x)) obtained by activating a jet with an intensity c(t) = c * = -0.05 and defining u u u c c c (x x x) as:

u u u c c c (x x x) = 1 c * u u u ′ ′ ′ (x x x) -u u u 0 0 0 (x x x) , (28)
Denoting {φ φ φ n } n=1...N r the N r retained modes obtained by applying the POD to (w w w(x x x,t i )) i=1...N t , the low-dimensional solution is written :

u u u(x x x,t) = u u u 0 0 0 (x x x) + c(t)u u u c c c (x x x) + N r ∑ n=1 a n (t)φ φ φ n (x x x) (29)
The Galerkin projection of the Navier Stokes equations onto the POD modes yields to the low order model [START_REF] Weller | Feedback flow control design by low-order modeling for laminar wake stabilization[END_REF] :

                     ȧr (t) = A r + C kr a k (t) + B ksr a k (t)a s (t) + P r + E r ċ(t) + F r c 2 (t) + G r c(t) + H kr a k (t)c(t) a r (0) = a 0 r 1 ≤ r ≤ N r (30)
INRIA where:

A r = -((u u u 0 • ∇)u u u 0 , φ φ φ r ) + 1 Re (∆u u u 0 , φ φ φ r ) B ksr = -((φ φ φ k • ∇)φ φ φ s , φ φ φ r ) C kr = -((u u u 0 • ∇)φ φ φ k , φ φ φ r ) -((φ φ φ k • ∇)u u u 0 , φ φ φ r ) + 1 Re (∆φ φ φ k , φ φ φ r ) E r = (u c , φ φ φ r ) G r = -((u u u 0 • ∇)u c , φ φ φ r ) -((u c • ∇)u u u 0 , φ φ φ r ) + 1 Re (∆u c , φ φ φ r ) F r = ((u c • ∇)u c , φ φ φ r ) H rk = ((u c • ∇)φ φ φ k , φ φ φ r ) + ((φ φ φ k • ∇)u c , φ φ φ r ) P r = (∇p, φ φ φ r )
We note, as in [START_REF] Weller | Feedback flow control design by low-order modeling for laminar wake stabilization[END_REF], that the modes satisfy the continuity equation. Thus the pressure term P r is equal to to ∂ Ω pφ φ φ r ds. If velocity field is constant at the boundaries, the POD modes are zero there. The pressure term therefore disappears completely.

The POD basis and the resulting model is built using the flow fields w w w(x x x,t), Eq. ( 27), collected using different control laws which derive from different sets of feedback gains. The POD model is calibrated using all the simulations carried out to collect the snapshot database, and the conditioning of the calibration procedure is improved as proposed in [START_REF] Weller | Feedback flow control design by low-order modeling for laminar wake stabilization[END_REF]. In particular, all the terms of the projection matrices are calibrated, with the exception of the convective terms B krs . Moreover, it is imposed that the steady unstable solution u u u 0 0 0 is also a steady solution of the reduced order model and, consequently, the term A r is forced to vanish.

When the feedback control is found using the velocity field of the POD model, Eq. (26) becomes:

c(t) = N v ∑ j=1 K j v(x x x j ,t) = N v ∑ j=1 K j v 0 (x x x j ) + c(t)v c (x x x j ) + N r ∑ r=1 âr (t)φ r v (x x x j ) (31)
where φ r v (x x x j ) are the values of the v-component of the POD modes at the sensors. Note that when steady unstable solution is used as target solution u u u 0 , because of the symmetry, v 0 (x x x j ) = 0 and that c(t) can be found in explicit form from Eq. (31) by trivial manipulation. After algebric manipulation, the low order model in matricial form becomes :

                     ȧ a a(t) = I I I -E E EK K K(I I I -K K Kv c (x x x v )) -1 φ φ φ v (x x x v ) -1 (A A A + C C Ca a a(t) + B B Ba a a(t)a a a(t) + K K K(I I I -K K Kv c (x x x v )) -1 φ φ φ v (x x x v )a a a(t) F F F K K K(I I I -K K Kv c (x x x v )) -1 φ φ φ v (x x x v ) a a a(t) +G G G + H H Ha a a(t)))
a a a(0) = a a a 0 (32)

where x x x v , the vector of the positions of the sensors, and K K K, the set of feedback gains, are used as input parameters.

In order to perform a stability analysis of the target state u u u 0 and to perform an optimisation of the feedback control gains, the POD model is linearized around the equilibrium state a a a * = 0, (which corresponds to the flow field u u u 0 ):

     ȧ a a(t) = L L L(K K K, x x x v )a a a(t) a a a(0) = a a a 0 (33)
where:

L L L(K K K, x x x v ) = I I I -E E EK K K(I I I -K K Kv c (x x x v )) -1 φ φ φ v (x x x v ) -1 C C Ca a a(t) + K K K(I I I -K K Kv c (x x x v )) -1 φ φ φ v (x x x v )a a a(t) G G G
Since the system matrix L L L of the linearized model depends explicitely on the feedback gains and on the position of the sensors, the model is predictive even when those parameters are changed with respect to the reference ones used for calibration. As already stated, the robustness of the model can be increased if, before linearization, a calibration procedure is used including several dynamics chosen by any sampling method, as detailed in [START_REF] Weller | Robust model identification of actuated vortex wakes[END_REF].

The linearized equation ( 33) can be used to perform a classical linear analysis of the dynamical system. Given the position of the sensors and the set of feedback gains K K K, the stable/unstable eigenvalues of the system L L L can be evaluated. For each eigenvalue, the associated eigenvector leads, by means of Eq. ( 27), to an estimation of the corresponding global mode of the linearized Navier-Stokes operator. A good accuracy on the estimation of the unstable modes of the full linearized Navier-Stokes problem allows to use the low order model in a control procedure, as described in the following. Note that the linearized reduced order model is obtained by using a simulation of a non-linear Navier-Stokes code. Moreover, the system matrix L L L depends non-linearly on the feedback gains K K K and on the position of the sensors x x x v , and this does not permit to use classical tools for the linear control. Thus, we propose here an iterative control procedure based on the minimisation of a functional cost, which is described in sec §4.2.

As explained above, the accuracy of the linearized model is an important aspect, and this is briefly investigated in the following. As a first step, it is shown how to reconstruct a global mode associated to an eigenvector of the linearized POD system. The formal solution a a a(t) of the system (33) is:

a a a(t) = R R Re Λ Λ Λt R R R -1 a a a 0 (34)
where Λ Λ Λ is the diagonal matrix of the eigenvalues of L L L, R R R is the matrix whose columns are the corresponding eigenvectors and a a a 0 is the initial condition on a a a(t). When INRIA Eq. ( 34) is substituted in Eq. 27, the fluctuating part of the velocity field u u u ′ (x x x,t) = u u u(x x x,t)u u u 0 0 0 (x x x) is obtained as follows:

u u u ′ (x x x,t) = Q Q QR R Re Λ Λ Λt R R R -1 Q Q Q -1 u u u ′ (x x x, 0) (35)
with

Q Q Q = K K K(I I I -K K Kv c (x x x v )) -1 φ φ φ v (x x x v )
and u u u ′ (x x x, 0) the projection of the initial condition over the POD modes. Thus, assuming that the eigenvalues of the physical system are well approximated by the low order model, we can reconstruct the matrix containing physical eigenmodes :

P P P ≈ P P P = Q Q QR R R (36)
In particular we are interested in the estimation of the unstable modes, which correspond to eigenvalues with positive real part.

In order to asses the accuracy of the feedback linear model described above, we consider a Reynolds number Re = 85, at which the instability is fully developed after a slow transient. In Figure 9 time evolution of the lift coefficient calculated on the cylinder with no control actuation is plotted. We recall that the simulation is carried out by a non-linear Navier-Stokes code. Note the quick growth of the C l after the slow transient regime. In the figure is highlighted, using a continous line, the portion of the transient used to build perform the POD model, which is sampled considering N t = 250 snapshots. This time interval is chosen starting when the lift coefficient reaches a value of C l ≈ 0.001 and including about six quasi-periodic flow oscillations. This choice is motivated by the need of capturing only the most energetic oscillations around the steady state while the system remains in a flow regime for which a linear approximation is still representative. We retain only N r = 6 POD modes to build and calibrate the linearized low order model. This is motivated by the work documented in [START_REF] Galletti | Accurate model reduction of transient and forced flows[END_REF], where it is shown that a model similar to the one built here gives a good approximation of the unstable mode. Thus, the unstable mode estimated by the POD model can be analyzed to explore his observability and to consequently choose the position of the sensors for the feedback control. In particular we used only one sensor of vertical velocity, which is placed in (x = 0.3, y = 0.0), in the area of the first local minimum (maximum in terms of module) of the v-component of the unstable mode. In figure 10 the v-component of the reconstructed mode is plotted; note that very close to the cylinder, before x ≈ 2.5, the value of v is very small; thus, placing feedback sensors near the cylinder is not an optimal choice to control the unstable mode. In order to test the capability of the feedback linear low order model to estimate the physical unstable mode in the presence of an actuation, we performed two numerical simulations of the actuated flow using two different proportional feedback gains for the sensor placed as described above, i.e. k = 0.1 and k = 0.2. In figure 11 simulation are shown together with the two portions of the transient used to build the POD database (solid line). As in the previous case, the two time intervals include six flow oscillations starting from a value of C l ≈ 0.001, with N t = 250 snapshots for each case. Note that the two initial gains are chosen in a random manner and the instability is not stabilized with those proportional paramters, even if the growth is retarded when k = 0.2 is used. The low order model is calibrated over the two dynamics and then the eigenvalues and the estimation of the unstable modes are carried out using the feedback adaptive linearized model for k = 0.1 and k = 0.2. In figure 12 (half of) INRIA the spectrum of L L L(k), for k = 0.1 and k = 0.2, is sketched together the unstable mode estimated by a linearized analsis of the Navier-Stokes operator (denoted in the figure as DNS). As expected, only two unstable conjugate eigenvalues are predicted by the linear low order model. Note that, as the value of the feedback gain is increased, the unstable eigenvalues are displaced closer to the stable region of the complex plane. The estimation of the unstable eigenvalues given by reduced model is very accurate as well as the effect of the increase of the feedback gain on the instability. The percentage error on the estimation of the real and the imaginary part of the unstable eigenvalues are respectively 7.62% and 0.26% when k = 0.1 is used and 0.11% and 0.12% for k = 0.2. Note that the estimation of the frequency of the instability (related to the imaginary part of the unstable eigenvalues) is almost perfect. This is due to the fact that the variation of the frequency is negligible between the two dynamics. In figure 13 the module of the reconstructed unstable mode for the case k = 0.2 and the one found by a linearized analysis of the Navier-Stokes operator are plotted. The prediction of the mode is very accurate in the whole domain; only a little different can be noted at the outflow due to the influence of the imposed boundary conditions in the linear Navier-Stokes code. An analogous result is obtained in the case k = 0.1.

Design of a control strategy based on the linear model

In this section we describe a control optimisation procedure based on the linear feedback low order model. In order to stabilize the steady state the unstable eigenvalues needs to be moved in the stable region of the complex plane. To this aim, while the position of the sensors are kept constant, a function of the gains K K K is proposed, such that its minimisation is equivalent to stabilize the system:

F (K K K) = ∑ N r r=1 tanh(Re(λ r (K K K)) -λ * Re ) + α K ∑ N v j=1 min l=1,...,N v ((K K K j -K K K 0 l ) 2 ) (37)
where λ r are all the N r eigenvalues predicted by the linear feedback model as K K K varies, λ * Re is the stability margin required, K K K 0 0 0 is the set of gains used to build the model the RR n°7092 The minimisation gives an optimal set of parameters K K K * for the present model. This set of gains are tested in a non-linear Navier-Stokes simulation of the transient. If the target state is not stabilized a new reduced order model is built with a database obtained by adding a portion of the transient of the new dynamics to the old POD database. During the optimisation procedure a maximum number of dynamics in the POD database can be fixed a priori, then when the maximum number is reached, a new set of snapshots substitues the one with maximum distance |K K K -K K K * |. Again, a minimisation of the functional (37) is carried out and a new set of parameters are obtained. The procedure is stopped when the steady state is stabilized.

In the test described here, the model built using the databases obtained with k = 0.1 and k = 0.2 is initially used for the optimization. The minimisation of 37 gives a new value of the feedback gain k * = 0.44. A non-linear simulation of the Navier-Stokes equations starting from u u u 0 is carried out, and the resulted flow is completely stabilized, as shwon in figure 14 and15. In figure 14 the lift coefficients obtained with k = 0.1, k = 0.2 and k * are plotted. The use of the optimised feedback gain leads to a stady and vanishing lift coefficient. Thus, the flow is totally controlled as displayed in figure 15, where the vorticity field of the flow obtained with k * at time t = 480 is shown.

Finally, the reduced order model obtained by a non-linear Naiver-Stokes code and then linearized around a steady state, is able to represent, with limited computational costs, the unstable modes of the linearized Navier-Stokes operator, and a control opti- misation based on such a linearized model gives a set of input parameters that stabilizes the actual flow. We recall that the whole procedure can be performed starting from the simulations of a generic non-linear code as those tipically used in enginnering applications.

The obtained results allow us to use the optimisation based on the linear feedback low order model in a control procedure for flow at higher Reynolds numbers and with an higher number of sensors. The main difficulty is to build a POD model which is robust with the parameters variation. Indeed, the most significant computational cost of the procedure is the update of the linear model, which need new DNS simulations. Indeed, a more robust model is characterized by a wider trust region and needs a reduced number of updates to complete the optimization. For this reason, the goal of future works is to couple the two techniques described in this study, i.e. to build a robust low order model to be used in the linearized control design strategy.

  (a) Flow configuration and isocontour of vorticity of a snapshot at Re = 180 (dashed lines represent negative values). (b) Placement of synthetic jet and sensors for control

Figure 1 :

 1 Figure 1: Sketch of the flow configuration with control actuation.

  • case A, correspond to an initial database U [Re 1 ] composed by N t snapshots collected at Re 1 = 120 with Re L = 70 and Re R = 180; • case B, correspond to an initial database U [Re 1 ] composed by N t snapshots collected at Re 1 = 100 with Re L = 40 and Re R = 180; • case C, correspond to an initial database U [Re 1 ;Re 2 ] composed by N t snapshots collected at Re 1 = 40 plus N t snapshots collected at Re 2 = 180 with Re L = 40 and Re R = 180;

65 Figure 2 :

 652 Figure 2: Sketch of the three test cases for sampling.

2 Figure 3 :

 23 Figure 3: Evolution of the error U U U ′ 2 versus the Reynolds number.

2 R

 2 ( U U UDNS) 2 (c) Test case C.

Figure 4 :

 4 Figure 4: Comparison between the mean projection error U U U ′ 2 and the mean residuals R( U U U DNS ) 2 for the three test cases under consideration.

2 R

 2 ( U U UROM) 2 (c) Test case C.

Figure 5 :

 5 Figure 5: Residuals obtained by a POD base built using Re = (40,180). Reynolds number considered between 40 and 180. Both the ROM predicted residuals (circle) and the DNS projection residuals (star) are shown

∆Res

  Non Calibrated ROM 2.1001 • 10 -04 Calibrated ROM 9.6494 • 10 -05

  100, 55, 160 • Case B: I N f = 120, 80, 165 • Case C: I N f = 40, 180, 90, 130 For Constraint Uniform Sampling (CUS), approximated onto the discretized space I h , we have : • Case A: I CUS = 100, 70, 140 • Case B: I CUS = 120, 90, 150 • Case C: I CUS = 40, 180, 85, 135

Fig 7

 7 shows the reconstruction error (21) obtained using the CCVT and the CUS strategies over the Reynolds number subspace I h for the three test cases A, B and C.

  40 180) 90 130 -Re eq 40 180 85 135

Figure 6 :

 6 Figure 6: Comparison of CCVT and CUS reconstruction error for test cases A, B and C. Both the ROM prediction reconstruction (circle) and the DNS projection (star) errors are shown.

Figure 7 :

 7 Figure 7: Histogram of difference between the reconstruction errors and standard deviations using CUS and CCVT sampling.

2 R

 2 ( U U UROM) 2 (c) Test case C.

Figure 8 :

 8 Figure 8: Comparison between projection and prediction residuals after CCVT sampling.

Figure 9 :

 9 Figure 9: Lift coefficient C l time evolution, with no control actuation at Re = 85.

Figure 10 :

 10 Figure 10: Reconstructed component v of the physical unstable mode.

Figure 11 :

 11 Figure 11: Lift coefficient C l time evolution, with proportional gains k = 0.1 and k = 0.2 at Re = 85.

Figure 12 :

 12 Figure12: Spectrum of the eigenvalues of the POD system matrix L vs. unstable physical eigenvalues.

Figure 13 :

 13 Figure 13: Isocontour of the module of the predicted and the physical unstable eigenmode for the case k = 0.2. Plots obtained with the same scale level.

Figure 14 :

 14 Figure 14: Lift coefficient obtained with k = 0.1, k = 0.2 and k * . Sensor position (0.3,0.0) and Re = 85.

Figure 15 :

 15 Figure 15: Vorticity snapshot of controlled flow with k * at time t = 480.

Table 2 :

 2 CCVT sampling efficiency ∆E and robustness ∆R.

		A		B		C
		∆E	∆R	∆E	∆R	∆E	∆R
	DNS Proj 16.150 47.711 11.120 123.373 -4.577 -7.991
	ROM	12.850 41.470 10.854 125.741	9.150 19.804

RR n°7092

In a general way, we have Φ Φ Φ i , Φ Φ Φ j Ω = δ i j , but not φ φ φ i , φ φ φ j Ω = δ i j . So, L i j = δ i j .RR n°7092


Quantum periods -I. Semi-infinite variations of Hodge

structures.


Fabien Lotte
email: fprlotte@i2r.a-star.edu.sg

Cuntai Guan
email: ctguan@i2r.a-star.edu.sg

An Efficient P300-based Brain-Computer Interface with Minimal Calibration Time

Introduction

Brain-Computer Interfaces (BCI) are communication systems that enable subjects to send commands to computers by using only their brain activity [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. Most existing BCI are based on Elec-troEncephaloGraphy (EEG) as the measure of brain activity [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. So far, BCI have been proven to be very promising communication and control tools for disabled people [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. A promising brain signals used in the design of assistive BCI for disabled people is the P300, a positive waveform occuring roughly 300 ms after a rare and relevant stimulus [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF][START_REF] Elshout | Review of brain-computer interfaces based on the P300 evoked potential[END_REF]. In order to use a P300-based BCI, subjects have to focus their attention on a given stimulus randomly appearing among many others, each stimulus corresponding to a given command. The appearance of the desired stimulus being rare and relevant, it is expected to trigger a P300 in the subject's brain activity. As such, detecting the P300 enables the system to identify the desired stimulus and hence the desired command. Interestingly enough, P300-based BCI have been successfully used to control a wheelchair (see, e.g., [START_REF] Rebsamen | Controlling a wheelchair indoors using thought[END_REF]) or to enable severely disabled users to spell words [START_REF] Elshout | Review of brain-computer interfaces based on the P300 evoked potential[END_REF][START_REF] Nijboer | A P300-based brain-computer interface for people with amyotrophic lateral sclerosis[END_REF].

However, current P300-based BCI as well as other BCI systems still suffer from several limitations which prevent them from being widely used [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. One of these limitations is that to use a BCI, many examples of the subject's EEG signals must be recorded in order to calibrate the BCI, which is unconvenient and time consuming. Moreover, this calibration process generally has to be repeated at regular intervals (e.g., from one day to the other) in order to accomodate sources of variations such as changes in electrode positions or changes in the subject's mental state and fatigue level. Therefore, the calibration time should be maintained as brief as possible. Until now, reducing the calibration time of P300-based BCI has been scarcely addressed by the literature. Exceptions are the works of Li et al [START_REF] Li | A self-training semi-supervised svm algorithm and its application in an eeg-based brain computer interface speller system[END_REF] and Lu et al [START_REF] Lu | Unsupervised brain computer interface based on inter-subject information and online adaptation[END_REF]. Li et al suggested to use initially a BCI calibrated with few training samples, and then to incrementally adapt this BCI online, thanks to semi-supervised learning [START_REF] Li | A self-training semi-supervised svm algorithm and its application in an eeg-based brain computer interface speller system[END_REF]. Lu et al proposed to use a subject-independent BCI, previously learnt from the data of many other subjects, also followed by online adaptation [START_REF] Lu | Unsupervised brain computer interface based on inter-subject information and online adaptation[END_REF]. However, the main limitiation of these two approaches is that such BCI would have initially poor detection performances, becoming efficient only after adaptation. An ideal P300-based BCI would have initially high performances, even if trained with very few examples. In this paper, we propose a new P300-based BCI design which can be trained using much fewer examples than current BCI designs, without sacrifying the detection performances.

Method

Preprocessing of EEG signals

In order to classify the EEG signals as P300 or non-P300, we use the EEG epochs from 150 ms to 500 ms after each stimulus, as this epoch should contain the P300, if any. EEG signals from each channels are then low-pass filtered with cut-off frequency at 25 Hz, the P300 being a slow wave [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. Finally, EEG signals are downsampled to 50 Hz in order to reduce the dimensionality of the data. This kind of preprocessing is the one generally used in current P300-based BCI designs [START_REF] Krusienski | A comparison of classification techniques for the P300 speller[END_REF][START_REF] Lu | Unsupervised brain computer interface based on inter-subject information and online adaptation[END_REF][START_REF] Thulasidas | Robust classification of EEG signal for brain-computer interface[END_REF][START_REF] Nijboer | A P300-based brain-computer interface for people with amyotrophic lateral sclerosis[END_REF]. With these designs, the samples from all EEG channels are concatenated into a single feature vector which is then fed into a classifier. In this paper, we term these features the "standard features". However, with our approach, the standard features are not the one fed to the classifier. Rather, we propose to perform a further feature extraction step based on Canonical Correlation Analysis (CCA).

Feature Extraction with regularized Canonical Correlation Analysis

CCA is a machine learning method used to find two directions w X and w Y which maximize the correlation between two random variables X and Y , once projected onto these directions [START_REF] Hardoon | Canonical correlation analysis: An overview with application to learning methods[END_REF]. It can be cast as the following optimization problem:

max w X ,w Y w T X X T Y w Y Xw X Y w Y (1)
where X and Y are represented as data matrices, with examples as rows and features as columns. This problem can be solved as a generalized eigenvalue problem [START_REF] Hardoon | Canonical correlation analysis: An overview with application to learning methods[END_REF]. If the first directions w X and w Y do not capture all the relationship between the two variables, multiple directions can be identified. One advantage of CCA over the classical correlation is its ability to find linear relationships between multivariate variables whereas the classical correlation can only identify linear relationships between univariate variables.

In this paper we propose to use CCA as a linear feature extraction method for P300-based BCI. More precisely, we use CCA to find directions w F and w L which maximize the correlation between the "standard features" F , i.e., the EEG samples, and the class labels L (P300 or non-P300). In other words we use CCA to find linear combinations of the EEG samples which best explain the class labels. Thus, the new features are the standard features F projected onto a number of directions w F identified by CCA. Moreover, by considering only a few directions, CCA gives a new subset of features with a reduced dimensionality, making it more suitable for learning with few examples.

It should be noted that solving equation 1 requires the estimation of the sample covariance matrix of the data. Unfortunately, learning with few examples may cause this sample covariance matrix to be poorly estimated. To overcome this issue, we propose to improve the robustness of CCA by using a regularized estimator for the covariance matrix. More precisely, we used the estimator of Ledoit and Wolf [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] in order to regularized this covariance matrix in a fully automatic way. We denote the resulting CCA as Regularized CCA (RCCA).

Classification using regularized Linear Discriminant Analysis

We used a Linear Discriminant Analysis (LDA) as our classifier. LDA uses a linear hyperplane to separate data from two classes, assuming that data are normally distributed [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces[END_REF]. Despite its simplicity, this classifier is widely used for BCI designs as it has been shown to be one of the most efficient [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces[END_REF], especially for the design of P300-based BCI [START_REF] Krusienski | A comparison of classification techniques for the P300 speller[END_REF].

LDA also requires the estimation of the sample covariance matrix of the data, which may be an issue when learning from few examples. Therefore, we also regularized the LDA by using Ledoit and Wolf's estimator of the covariance matrix [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] instead of the classical sample estimator. We denote the resulting LDA as Regularized LDA (RLDA).

3 Evaluation

EEG Data acquisition

EEG signals were recorded from 10 healthy subjects who participated in 2 sessions of the P300 speller, as described in [START_REF] Thulasidas | Robust classification of EEG signal for brain-computer interface[END_REF][START_REF] Lu | Unsupervised brain computer interface based on inter-subject information and online adaptation[END_REF]. With this speller, subjects can spell characters by focusing their attention on a given character displayed on a computer screen, within a 6 × 6 matrix of characters. The rows and columns of this matrix are successively and randomly flashing. When the row or column flashing is the one that contains the desired character, this triggers a P300 in the subject's brain activity. Indeed, this flash is a rare stimulus, relevant to the subject. The detection of this P300 enables the system to identify which row and which column contain the desired character, and as such, what is this desired character. In each session, subjects had to spell 41 given characters. For each character, ten rounds of flashing were performed, each round corresponding to the flashing of each row and each column once. The first session was used as the training set, whereas the second one was used as the test set. For each subject, both sets were composed of 4920 examples (41 characters × (6 rows + 6 columns) × 10 rounds of flashing), 820 of which corresponding to a P300 signal, i.e., the row or column flashing contained the desired character.

EEG signals were recorded using 64 EEG channels, sampled at 250 Hz. For the purpose of this study, we only considered channels Fz, Cz, Pz, PO7, PO8, OZ, P3 and P4 as they were identified as the ones giving the best classification accuracy for the P300 speller [START_REF] Krusienski | A comparison of classification techniques for the P300 speller[END_REF].

Results

In this study, we did not directly evaluate the classification accuracy of our BCI for the P300 speller. Indeed, in the classical P300 speller, EEG epochs are averaged over the different rounds of flashing, which makes the classification task easier but the application slower [START_REF] Thulasidas | Robust classification of EEG signal for brain-computer interface[END_REF]. Here, we evaluated the performance of our BCI for the more challenging task of detecting the P300 in "single-trial", i.e., without averaging. In other words, we used our BCI to decide whether each EEG epoch following a flash contains a P300. This problem being a detection problem, we classically used the Area Under the ROC (Receiver Operating Characteristic) Curve (AUROCC) as the performance measure.

As an evaluation of our BCI design, we computed the AUROCC obtained on the full test set when using various number of training characters to calibrate the BCI. It should be reminded that a character corresponds to 120 examples, 20 of which belongs to the P300 class. We compared the performances obtained by our design with that of the classical P300-based BCI design, i.e., the design using the standard features as input to an LDA classifier. We also compared it with a BCI design using Principal Component Analysis (PCA) for feature extraction, a technique also popular for the design of P300-based BCI (see, e.g., [START_REF] Thulasidas | Robust classification of EEG signal for brain-computer interface[END_REF]). For both CCA-based and PCA-based feature extraction, we extracted 50 features from the initial 136 EEG samples of each epoch (8 channels × 17 samples per channels). To assess the contribution of the different components of our BCI, we also performed these evaluations without using CCA, and/or without using regularization for CCA and/or LDA. The results obtained are displayed on Figure 1.

Discussion and conclusion

Results suggest that the BCI design we proposed is more efficient than both the standard and the PCA-based BCI design, especially when the number of training characters used is small. Interestingly enough, a paired t-test revealed that the average AUROCC (over all subjects and all number of training characters) obtained by our method was significantly higher than that obtained by the standard design and the PCA-based design (p << 0.001).

When using only 3 characters for training, our design reached an average AUROCC of 0.87, i.e., an AUROCC only 0.02 lower than that obtained by the standard design when using all 41 training characters. To reach an AUROCC of 0.87, the standard P300-based BCI design needs

Figure 1 :

 1 Figure 1: Average AUROCC (over the 10 subjects) versus number of characters used for training, for various methods. For reasons of readability, only values in [0.72, 0.9] are displayed on the y-axis.

  5 times as many training characters as with our design, whereas the PCA-based design needs to use all 41 training characters. This suggests that our design can indeed significantly reduce the calibration time as it needs only a few training characters to reach high performances. For instance, collecting 41 training characters requires approximately 20 minutes. With our design, we just need to collect about 3 or 5 training characters to reach the same level of performance, meaning a time for data collection of roughly 2 minutes. As RCCA and RLDA are linear, fast and fully automatic methods, the whole BCI calibration remains computationally very efficient. Results also suggest that each component of our BCI (i.e., CCA, regularization of CCA and regularization of LDA) contributes to the overall performance improvement when learning from few examples. Indeed, the RCCA+RLDA design outperforms the CCA+RLDA design which in turn outperforms the "standard feature"+RLDA design, itself outperforming the "standard feature"+LDA design.To conclude, in this paper we have proposed a new design for the P300-based BCI, in order to reduce the calibration time of the system. Our BCI is based on Regularized Canonical Correlation Analysis for feature extraction and Regularized Linear Discriminant Analysis for classification. Evaluations have suggested that this design can reach good P300 detection performances while using much less training examples than current approaches, hence effectively reducing the calibration time. Future work will explore the use of this approach in actual applications for disabled users as well as nonlinear feature extraction based on kernel CCA[START_REF] Hardoon | Canonical correlation analysis: An overview with application to learning methods[END_REF].

Acknowledgments

The authors would like to thank Dr. David R. Hardoon for his feedback and help about the CCA method, as well as Dr. Brahim Hamadicharef for his comments about this manuscript.