Pei Li
email: pei.li@telecom-sudparis.eu

Elisabeth Brunet
email: elisabeth.brunet@telecom-sudparis.eu

Raymond Namyst
email: raymond.namyst@inria.fr

High Performance Code Generation for Stencil Computation on Heterogeneous Multi-device Architectures

Keywords: GPGPUs, OpenCL, Stencil computations, Multidevice, Code generation, Heterogeneous architectures

Heterogeneous architectures have been widely used in the domain of high performance computing. On one hand, it allows a designer to use multiple types of computing units and each able to execute the tasks that it is best suited for to increase performance; on the other hand, it brings many challenges in programming for novice users, especially for heterogeneous systems with multi-devices.

In this paper, we propose the code generator STEPOCL that generates OpenCL host program for heterogeneous multidevice architecture. In order to simplify the analyzing process, we ask user to provide the description of input and kernel parameters in an XML file, then our generator analyzes the description and generates automatically the host program. Due to the data partition and data exchange strategies, the generated host program can be executed on multi-devices without changing any kernel code. The experiment of iterative stencil loop code (ISL) shows that our tool is efficient. It guarantees the minimum data exchanges and achieves high performance on heterogeneous multi-device architecture.

I. INTRODUCTION

High performance computing (HPC) is closely tied to scientific computing and industries. Because of the increasing complexity and growing amount of data for practical problem, we always demand better performance to achieve faster time to solution. Hence, there are two paths to combine: the enhancement of the algorithms and better material architectures. Recently HPC system architectures are shifting from the traditional homogeneous multi-core systems to heterogeneous systems such as GPGPU. Compared to the standard multi-core CPUs, GPGPUs offer a significantly higher floating point peak and a better power efficiency.

However, this novel architecture presents new challenges at the application developing level: time and effort are needed to exploit such kind of material. Moreover, many heterogeneous systems start having multiple computing devices. It can be more difficult and error-prone since developing programs that make best use of the characteristics of different computing devices increases the programmer's burden. Balancing the workload between several available computing devices can be also complicated, especially given that they have different performance characteristics. Besides, the communication and the exchange of intermediary results between several devices should also be considered. It is costly and difficult to design applications for heterogeneous multi-device systems. Thus, there is a huge demand for programming tools that help the novices designing applications for heterogeneous multi-device systems.

In this paper, we propose the code generator STEPOCL that automatically generates the parallel host OpenCL code for heterogeneous multi-device systems. It enables OpenCL programs written for a single compute device to run on systems with multiple devices without any modification. The architecture of the system is completely transparent to user. The information of available devices is obtained at run time and the workload is distributed to each device with optimal strategy. Thus, OpenCL kernels are executed in parallel. The host program manages automatically the communication and data exchanges between devices and results are retrieved from each device at the end of execution.

The rest of paper is organized as follows. In Section II, we present our contribution in detail. Section III discusses the evaluation of our generated code. Section IV presents related works. Finally, Section V concludes the paper.

II. STEPOCL STEPOCL aims to facilitate programming on heterogeneous multi-GPUs systems through the open standard

OpenCL [START_REF] Munshi | OpenCL Programming Guide, ser. OpenGL[END_REF], a language especially designed to address heterogeneous platforms consisting of multi-core CPUs, GPUs and other modern processors. Instead of writing the errorprone code for heterogeneous multi-device system, user only needs to provide the basic description of kernel argument, space information and the kernel function for one device. From these description, STEPOCL automatically generates an entire OpenCL source code for multiple devices architectures dealing with all the necessary technical aspects including not only the basic ones with the tuning of the initialization phase (library, devices declaration, etc.), or kernels launching and retrieving of the results, but also trickier aspects in order to determine the best data and computation distribution or the exchange of intermediary results.

f o r ( i n t t = 0 ; t <T ; ++ t ) { f o r ( i n t i = 1 0 ; i <N-10;++ i ) { f o r ( i n t j = 1 0 ; j <N-10;++ j ) { A[ i ] [ j ] =CNST * ( B[ i ] [ j +1]+B[ i ] [ j -1] +B[ i -10][ j ] +B[ i + 1 0 ] [ j ] ) ; } } swap (A, B) ; } Listing 1. Stencil Loop Example
Indeed, if we consider an iterative stencil computationa case widely used in many scientific domains as depicted in Listing 1, the computation of each element needs the access to a set of neighboring elements according to a fixed pattern. Thus, in a multi-device version, the distribution may imply the allocation of some of those neighboring elements in the memory of another device. Then, some data need to be shared in some way -by exchange or by replicationand furthermore, need to be updated after each iteration in order to maintain the data coherency. All these operations make programming iterative stencil code on heterogeneous architectures more difficult and error-prone.

Thus, in order to generate a complete multi-device code, we need to face the following challenges:

• Workload and data partitioning: the kernel space and data space should be characterized in order to be efficiently partitioned by rows, columns or grids. The shared data regions [START_REF] Meng | A performance study for iterative stencil loops on gpus with ghost zone optimizations[END_REF] -called here ghost zones -should be distinguished from useful data.

• Managing the data transmission between devices: in order to maximize the reutilization of data transferred in devices memory, intermediary results need to be exchanged in a way that guarantees the minimum data transfer, that is to say only the effective data necessary to pursue the computation on the other devices.

• Code generation: it should be transparent to users. Users should not care about the number or the type of available devices. The OpenCL host code will be automatically generated from the description of kernel information for single device. The rest of this section presents more details about each of mentioned points and how we consider to implement them in STEPOCL.

A. Workload and Data partitioning 1) Background -OpenCL kernel space: An OpenCL application consists of two distinct parts: the host program and a collection of computation kernels addressing computation potentially heterogeneous devices such as CPUs or GPUs. Kernels are typically simple functions that transform input memory objects into output memory objects. Nevertheless, the host part is the driver of the execution: it is in charge of the data transfers to and from the devices memory and of kernels

G x =20 W x =4 G y =15 W y =3
Fig. 1. The kernel index space launching on the devices. When the host submits a kernel, the OpenCL runtime system creates an integer index space based on two arguments: the number of the work-item and the size of work-groups. An instance of the kernel is executed for each point in this index space.

Each instance of an executing kernel is called a work-item. Work-item are organized into work-groups. The work-groups provide a more coarse-grained decomposition of the index space. All the work-groups are the same size in corresponding dimensions, and this size evenly divides the global size in each dimension. A unique ID is assigned to each work-group following the same dimensionality as the index space used for the work-items. Each work-item is referenced by a unique local ID within a work-group, so that a single work-item can be uniquely identified by its global ID or by a combination of its local ID and work-group ID. For example, Figure 1 shows a two-dimensional index space with a global size of index space of (20 × 15) and a size of work-group of (5 × 5). Hence the number of work-groups in the whole index space is (4 × 3)

During the execution, each work-item uses the same sequence of instructions defined by a single kernel. In order to execute the kernel on multiple devices, the number of work-item should be adjusted according to the number of device. The corresponding data space should also be split and distributed to the memory of each device.

2) Determining the kernel and data space: As we said in introduction, for the first prototype of STEPOCL, we solicit the assistance of the user to describe how the space can be split. Therefore, the kernel and data space is described in XML as a property of arguments as depicted in Listing 2. Still on the same example, the argument A is defined as a 2 dimension memory object with a size of 4096 × 1026. The number of work-item in kernel space is 4096×256. After analyzing these description with a XML parser, we can easily determine the size of kernel and data space. Nevertheless, the objective of our further work is to automatically analyze the kernel and data space from the real kernel program thanks to efficient compilers as PIPS [START_REF] Irigoin | Semantical interprocedural parallelization: an overview of the pips project[END_REF] or Insieme [START_REF]Insieme compiler project[END_REF]. 3) Distinguishing the useful data from ghost zones [START_REF] Meng | Performance modeling and automatic ghost zone optimization for iterative stencil loops on gpus[END_REF] and keeping them up-to-date : If we manage an ISL distribution, the ghost zones need to be exchanged among different processing elements at the beginning of the execution and after each iteration involving significant overhead in terms of communication and synchronization. Hence, larger ghost zones may be created to replicate stencil operations, reducing communication and synchronization costs at the expense of redundantly computing some values on multiple processing elements [START_REF] Meng | A performance study for iterative stencil loops on gpus with ghost zone optimizations[END_REF]. The optimal size of ghost zone can improve the performance for ISL on GPUs. Nevertheless, the objective of our work is to maximize the memory utilization to allow the kernels scaling. In this way, in our model, we always consider the minimal ghost zone size -even if it may be the case for the data space, computation space will not be replicated-. Just as in the previous paragraph, in the current STEPOCL prototype, we ask the programmer to give the relative information but the goal is to take advantage of the data dependency analysis from compilers to do it automatically.

We extract the domain of ghost zone and useful data region from the data description in the XML file given by the programmer. Then we split the useful data region into subregion, and allocate relative ghost sub-regions for each data sub-regions. Still on the same example, ghost zone is required and may be described in the file as depicted below in Listing 3. It indicates that the offset of data region from ghost zone in X dimension is (-10,10) and in Y dimension is (-1,1). The data region and ghost zone are presented in Figure 2.

Meanwhile, by analyzing the shape of ghost zone, we propose a simple partitioning strategy that guarantees a minimum For example, the Listing 1 above presents a part of typical stencil loop code. In this case, updating the matrix A depends more on the data which are relatively allocated on X axes. If we partitioning the kernel space by row (Fig. 3) instead of partitioning by column (Fig. 4), we can reduce significantly the communication and data transfer.

4) Partitioning with abstract number of available devices:

The number of available devices is unknown until the runtime in order to be able to exploit correctly the actual devices involved in the execution. Thus at compilation time, the number can only be presented with variable num device. At first, we analyze the kernel space to evaluate the maximum parallel capacity. In Figure 1, the number of work-group in this kernel space is 12. If we want to parallelize stencil code without modifying kernel code, we need at most 12 devices. And the number of regions in dimension X should not exceed 4 (num group x), equally the number of regions in dimension Y should not exceed 3 (num group y). Afterwards with the dependent points which are provided by user, we can determine a priority list. The list presents the splitting order of dimensions. For instance, the dependent points in Listing 1 are (0,1), (0,-1), (10,0) and (-10,0). The projection of ghost region size on X is dx = (|10| + | -10|) = 20 and the projection on Y is dy = (1 + | -1|) = 2. By comparing the size of projection, we can decide which axis should be split first. In this example, we should split the data in row first, which means splitting dimension Y . And, if there are still available devices, we can split the data on dimension X afterwards. So, the list of priorities is list = {y, x}.

Following the list of priority, we calculate the greatest common divisor (gcd) of num workgroup i and the number of available devices num device for each dimension. num device is updated with num device/ = gcd(num workgroup i , num device) after each gcd operation. The list of gcd presents the data distribution, and the result of N i=0 gcd(i) is the number of devices that will be used for computing. In our application, we may decrease the number of total devices in a certain range (default value is 100) to achieve the maximum parallelization. In the example of Listing 1, with a priority list list = {y, x}, if OpenCL detected 2 available devices, the list of gcd should be {1,2}. This means that the data in dimension X remain in one part, but the data in dimension Y will be split into two parts, as shown in Figure 3. After the data partition, each segment (we call it local zone) should keep the following information: the global zone ID, the relative ID in each dimension and the range of indexes in each dimension. Each local zone is composed of written data region (we call it the write-zone) and ghost zone (corresponding to shared data). Since we know the size of local region and all the dependence points, the index domain of write-zone can also be determined. Thereby, local regions represent EXACT READ data regions and write-regions, EXACT W RIT E ones. As ghost zones are partial write-zone projection of neigbors, after each iteration of stencil loops, they need to be updated to prepare the next iteration. Thus, we need to determine the neighbors global ID list for each local zone.

The global ID of zone region can be represented with relative ID (re ID x, re ID y) and device partitioning information (the number of devices in the first dimension numdev x, the number of devices in the second dimension numdev y).

global ID = re ID y × numdev y+ re ID x(2D : global ID) (1)
global ID = re ID z × numdev y × numdev x+ re ID y × numdev y+ re ID x(3D : global ID) (2) A dependent point indicates a direction of dependent zone. For instance, on Figure 5: the dependent point of device 3 is (-1,-1). We suppose that the length and the height of zone is M and N which are larger than 1. So the offset of relative ID is still (-1,-1). The ID of dependent zone is (1, 1) + (-1, -1) = (0, 0). So the global ID can be calculated as follows:

global ID = 0 × 2 + 0 × 2 = 0; (3)
In this way, we can find all the dependent zones with the list of dependent points. Then, the list of dependent zones (we call it neighbor list) will be used during the data transmission or communication process.

B. Managing the data transmission between devices

1) Read/write memory objects in OpenCL: OpenCL lets users create three kinds of memory objects: buffers, 2D images and 3D images. These memory objects are stored in the host memory (typically, in RAM) or in the device memory (typically, in GRAM directly on the graphic card). There are several functions that can be used to read and write memory object. The Table I presents five functions that read and write buffer object [START_REF] Munshi | OpenCL Programming Guide[END_REF]. 2) Data transmission between multiple devices: After each iteration, the data that needs to be transferred from device A to device B can be calculated in the following way:

region aT ob = exact write A ∩ exact read B (4)
The data region is an EXACT WRITE region; the zone region is an EXACT READ region. Thus, the formula [START_REF]Insieme compiler project[END_REF] can be redefined as:

region aT ob = data region A ∩ zone region B (5)
Then all the data can be transferred with the neighbor list following the process presented in Listing 4. OpenCL does not assume that data can be transferred directly between devices, so commands only exist to move from a host to device, or from a device to host. Copying data from one device to another requires an intermediate transfer to the host. The pointer of host memory cannot be simultaneously possessed by several devices, one device cannot communicate with host until the data transmissions of other devices have finished. This sequential process degrades the overall performance. Increasing additional ghost zone may reduce the frequency of communication, and this option is also provided in our generator.

C. Code generation

The general procedure is presented in Algorithm 1. In fact, the inputs of Algorithm 1 are described in an XML file. STEPOCL reads the script, and generates the host code without changing kernel code.

The generation procedure is composed of 3 parts: initialization of OpenCL environment, data partition and data transmission. OpenCL initialization identifies all available devices by creating corresponding Command Queues.

In the second part, subroutine DataP artition() split the data space into several regions with the partitioning strategy which are mentioned in Section II. According to the information of ghost zone and relative ID, we can predict with which region the data exchange will take place. Thus, subroutine f indN eighborRegion() creates a neighbor list for each region. After allocating the memory objects for each region CopyDeviceT oHost(P host , dev i , RegionList(i)); on corresponding devices and adjusting the kernel arguments, the kernels are launched simultaneously.

In the last part, the data transmission happens after each iteration. By following the neighbor list of each region, subroutine transf erData() transfers and updates the data for preparing the next execution. If all data is updated, the next execution of kernel is launched. If there are only one device, the neighbor list will be empty, and the data transmission is not permitted.

Two synchronizations are used in this algorithm. The first ensures that the data will not be transferred during writing process while the second ensures that all data are ready for the next iteration.

The generated host code can directly be executed with the original kernel.

III. EVALUATION

In order to evaluate our implementation, we generated a 4points Jacobi 2D stencil with 50 iterations. We executed it on a machine with 4 NVIDIA GTX-460 devices. The system is Red Hat Enterprise Linux Server release 5.5. We varied the size of the input data to highlight the new scalibity offered by our generated program.

P erf ormance = N um iter × Size date × F lt Opt/T n (6)
The N um iter presents the number of iterations, while F lt Opt means the number of operations in the inner loop and T n means the total time used for executing stencil code which also includes the time of communication between each computing device. On one hand, the results show that single GPU gains better performance with small data sizes as the number of work-groups in the kernel is not big enough to make all devices busy. In this case, the overhead of data transmission is significant in the whole execution time. Figure 8 helps to picture exactly the impact of the communication on the overall time execution in percentage (PTE): PTE decreases as the size of data growth. Though we use our partition strategy to avoid unnecessary extra data transmission, the PTE value is still very high. Thus, even if we will investigate in the short term how to take care of the communication -by overlapping them in order to avoid occupying CPU communication by several GPUs at the same time for example -, we above all want to take in consideration this saturation threshold as a parameter of a dynamic scheduling strategy with the purpose of determining the best compromise between data size and number of devices to use. The goal is to use only the necessary resources and not to just occupy them and avoid their exploitation for another computation. On another hand, as expected, the performance curves show aggressive growth with the increase of data. With four GPU devices, we achieved 61 GFLOPS which is 3.5 times faster than using only one GPU device and almost 1.8 times faster than using two GPU devices. The generated program even achieves to process a data input size which is impossible to treat on a single device program.

IV. RELATED WORK

A number of recent studies have focused on the parallelization of OpenCL code for multi-devices GPGPU systems. The Amdahl Software provides similar application -OpenCL CodeBench [START_REF]OpenCL Codebench by Amdahl Software[END_REF] . It enables developers to rapidly generate and optimize OpenCL code. The main difference between STEPOCL and OpenCL CodeBench is that the host code generated by OpenCL CodeBench is more general and users need to define their own ways of communication between several accelerators. Jungwon Kim and Honggyu Kim [START_REF] Kim | Achieving a single compute device image in opencl for multiple gpus[END_REF] propose an OpenCL framework that treats multiple GPUs as a single compute device. This framework analyzes the OpenCL kernel index space at run time and it performs a sampling run just before the kernel is executed. The sampling run obtains buffer access ranges of each affine array references for different GPUs. Using this information, the runtime distributes the kernel work-group index space efficiently. Sylvain Henry provides an OpenCL implementation which is called SOCL [START_REF] Henry | Programmation unifiée multiaccélérateur OpenCL[END_REF]. It is based on StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF]. It gives a unified access to every available OpenCL device: applications can now share entities such as Events, Contexts or Command Queues between several OpenCL implementations. In addition, the Command Queues that are created without specifying a device provide automatic scheduling of the submitted commands on OpenCL devices contained in the context to which the command queue is attached. On the other hand, this implementation use dynamic analysis and scheduling the available devices at runtime. Considering our data partitioning also happens at runtime, it will be very promising to combine our research achievement.

Data transfer between the CPU and GPUs can degrade the performance. Overlapping [START_REF] Kim | Achieving a single compute device image in opencl for multiple gpus[END_REF], [START_REF] Gelado | Cuba: an architecture for efficient cpu/co-processor data communication[END_REF] the data transfer and GPU computation is a solution to reduce the overhead of data transmission.

V. CONCLUSION

In this paper, we introduced the design and implementation of the new OpenCL code generator STEPOCL which provides the facilities for developing code on heterogeneous multi-device systems. Instead of developing a tedious code implementation, user only needs to provide the basic description of kernel argument, space information and the kernel function for one device. Then, STEPOCL automatically generates OpenCL code for multi-device without changing kernel functions. STEPOCL builds model of kernel and data space from the description then partition the workload with best strategy. Communication management and exchange of intermediary results between several devices are also generated. Preliminary experiments on an iterative stencil show that the generated code achieved high performance on multi-device architectures.

Further works are planned at different levels. First, we will massively enhance our tool by relaxing the information demanded to the user thanks to a cooperation with static compilation techniques: kernel and data space and ghost zones will be automatically detected when possible. Next, as STEPOCL is able to generate parametric kernel, it is able to generate non uniform distribution. We plan to collect information on available heterogeneous devices at runtime with the purpose of applying dynamic scheduling strategies and performing an accurate partitioning.

  <a r g u m e n t> <name>A</ name> <p r o p e r t y>o u t p u t</ p r o p e r t y> <d a t a T y p e> f l o a t</ d a t a T y p e> <cl mem>1</ cl mem> <a r g s i z e> <d i m s i z e>4098</ d i m s i z e> <d i m s i z e>1026</ d i m s i z e> </ a r g s i z e> </ a r g u m e n t> <k e r n e l> <name> s t e n c i l . c l</ name> <g l o b a l> <g s i z e>4096</ g s i z e> <g s i z e>256</ g s i z e> </ g l o b a l> < l o c a l> < l s i z e>64</ l s i z e> < l s i z e>16</ l s i z e> </ l o c a l> </ k e r n e l> Listing 2. Kernel arguments information

  <shadow> <p o i n t>{ -10 ,0}</ p o i n t> <p o i n t>{10 ,0}</ p o i n t> <p o i n t>{0 ,1}</ p o i n t> <p o i n t>{0 , -1}</ p o i n t> </ shadow> Listing 3. Ghost zone information

Fig. 2 .

 2 Fig. 2. Data space

Fig. 3 .Fig. 4 .

 34 Fig. 3. Partition the data by row

Fig. 5 .

 5 Fig. 5. The dependent points of device 3

  f o r e a c h d e v i c e A l o a d ( d a t a r e g i o n a ) from d e v i c e A f o r e a c h d e v i c e B i n n e i g h b o r l i s t A l o a d ( z o n e r e g i o n b ) from d e v i c e B t r a n s f e r D a t a ( r e g i o n a T o b ) w a i t ( ) ; Listing 4. The data transmission between all devices

Fig. 6 .

 6 Fig. 6. Data transmission between two devices

Algorithm 1 : 1 N 2 RegionList ← N ewRegionList; 3 RegionList← 5 foreach dev i in deviceList do 6 dev i ← 7

 1123567 Host code generation Input: N W I: Number of Work Item, W S: Work group size, DS: Data Size, DP L: Dependent Points List, DRL: Dependent Regions List, T I: Total Iteration Output: P host :Host Program P host ← N ewHostP rogram(); um Device ← clGetDeviceIDs; DataP artition(W S, DS, DL); 4 DRL ← f indN eighborRegion(RegionList, DP L); AllocateDeviceBuf f er i withSizeOf (RegionList(i)); P host ← 8 CopyHostT oDevice(P host , dev i , RegionList(i)); foreach i in T I do 9 N W I ← 10 KernelSpaceP artition(N W I, RegionList); P host ← 11 InvokeKernel(N W I, W S, OtherArguments); synchronization(); 12 foreach region i in RegionList do 13 foreach region b in DRL i do 14 P host ← transf erData(Region iT ob); 15 ∀device, swap(inputBuf f er, outputBuf f er); 16 synchronization(); 17 foreach dev i in deviceList do 18 P host ← 19

Fig. 8 .

 8 Fig. 7. Performance of stencils across architectures

TABLE I READ

 I

AND WRITE BUFFER OBJECTS


Kaushik Kumar Nagarajan
email: kaushik@ctfd.cmmacs.ernet.in

Laurent Cordier
email: laurent.cordier@univ-poitiers.fr

Christophe Airiau
email: christophe.airiau@imft.fr

Development and Application of a Reduced Order Model for the Control of Self-Sustained Instabilities in Cavity Flows

Keywords: AMS subject classifications: 93C05, 93C10, 76N25 Reduced order modelling, proper orthogonal decomposition, cavity flows, feedback control

Flow around a cavity is characterized by a self-sustained mechanism in which the shear layer impinges on the downstream edge of the cavity resulting in a feedback mechanism. Direct Numerical Simulations of the flow at low Reynolds number has been carried out to get pressure and velocity fluctuations, for the case of un-actuated and multi frequency actuation. A Reduced Order Model for the isentropic compressible equations based on the method of Proper Orthogonal Decomposition has been constructed. The model has been extended to include the effect of control. The Reduced Order dynamical system shows a divergence in time integration. A method of calibration based on the minimization of a linear functional of error, to the sensitivity of the modes, is proposed. The calibrated low order model is used to design a feedback control of cavity flows based on an observer design. For the experimental implementation of the controller, a state estimate based on the observed pressure measurements is obtained through a linear stochastic estimation. Finally the obtained control is introduced into the Direct Numerical Simulation to obtain a decrease in spectra of the cavity acoustic mode.

Introduction

Self sustained instabilities in compressible cavities arise due to the impingement of the shear layer on the downstream edge of the cavity, and is caused by the energetic coherent structure, resulting in a pressure feedback mechanism which induces the far field noise in the cavity. The flow is characterized by a strong coupling between the hydrodynamic instability and the acoustic propagation which can be sometimes intense. Applications of such flows arise in airframe fuel vents, landing gears, weapon bays. Accurate modelling of the cavity noise is important from the point of view of a search for quieter aircrafts. The traditional tools for the modelling of the cavity like Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) pose difficulties in terms of computational resources, also the constraint increases when the ultimate aim of study is in the application of control. The self sustained instabilities in cavity has been studied in the past by many authors [START_REF] Rowley | On self-sustained oscillations in two-dimensional compressible flow over rectangular cavities[END_REF][START_REF] Gloerfelt | Direct computation of the noise radiated by a subsonic cavity flow and application of integral methods[END_REF][START_REF] Bres | Three-dimensional instabilities in compressible flow over open cavities[END_REF][START_REF] Larcheveque | Large-eddy simulation of a subsonic cavity flow including asymmetric three-dimensional effects[END_REF][START_REF] Rowley | Dynamics and control of high-reynolds number flow over cavities[END_REF]. The principles of Reduced Order Modelling (ROM) based on Proper Orthogonal Decomposition (POD) and a Galerkin projection, has been in the past applied to model and control to flows with self-sustained instabilities, like wake behind cylinder [START_REF] Bergmann | Optimal rotary control of the cylinder wake using POD Reduced Order Model[END_REF] and high lift configurations [START_REF] Luchtenburg | A generalized mean-field model of the natural and high-frequency actuated flow around a high-lift configuration[END_REF]. Flow past an open cavity has already been studied using ROM by [START_REF] Rowley | Model reduction for compressible flows using POD and Galerkin projection[END_REF] and [START_REF] Gloerfelt | Compressible Proper Orthogonal Decomposition/Galerkin reduced order model of self sustained oscillations in a cavity[END_REF] but without any application to flow control. More recently, ROM for controlled configurations has been proposed by [START_REF] Weller | Robust model identification of actuated vortex wakes[END_REF][START_REF] Kasnakoglu | Reduced order modeling, nonlinear analysis and control methods for flow control problems[END_REF]. In [START_REF] Samimy | Feedback Control of Subsonic Cavity Flows Using Reduced-order Models[END_REF] the ROM for flows issued from an experiment has been used to design a controller. The major hurdle in using the ROM for control applications is the accuracy of the model in predicting the dynamics of the system even for short periods, also difficulty arises when the control parameters are changed as in a real time simulation. Various numerical strategies termed as calibration techniques has been developed in the recent past based on the solution of an optimization problem and can be found in the works of [START_REF] Bergmann | Optimal rotary control of the cylinder wake using POD Reduced Order Model[END_REF][START_REF] Galletti | Low-order modelling of laminar flow regimes past a confined square cylinder[END_REF][START_REF] Couplet | Calibrated Reduced-Order POD-Galerkin system for fluid flow modelling[END_REF][START_REF] Perret | Polynomial identification of POD based low-order dynamical system[END_REF][START_REF] Kalb | An intrinsic stabilization scheme for proper orthogonal decomposition based low-dimensional models[END_REF]. A summary of the application of various method has been discussed in detail [START_REF] Cordier | Calibration of POD Reduced-Order models using Tikhonov regularization[END_REF] where a method of calibration based on the Tikhonov regularization has been proposed. The main contribution of this paper is then to complete the full development as applied to cavities, like building up the ROM, including the effect of control, calibrating the model and finally performing control studies. We give further developments to the calibration technique by performing a sensitivity analysis to weight the importance of the errors in the different modes of the ROM, and accurately predict the dynamics of the system. The ROM thus obtained has been applied to perform the feedback control of cavity flows to obtain a reduction in the noise level of the cavity.

The paper is organized as follows. In Section 2 the derivation of the plant model used for control studies has been discussed. The general principles of model reduction is discussed in Subsection 2.1 where the usual ROM based on the principle of Galerkin projection is obtained as a special case. In Subsection 2.2 we discuss the general principles POD, followed by a discussion of the ROM for the isentropic equations in Subsection 2.2.1. Extension of ROM to include the effect of actuation is presented in Subsection 2.3. Section 3 concerns the calibration of the ROM to accurately represent the temporal dynamics, the various errors that define the dynamics of the ROM are presented in Subsection 3.1. A new approach to Tikhonov based regularization based on the sensitivity analysis is presented in Subsection 3.2 to obtain an accurate representation of the dynamics. The ROM has been applied to the cavity flow configuration in Section 4. Finally feedback control law based on the estimation of the observer dynamics has been presented in Section 5.

Derivation of the plant model 2.1 Model reduction principles

To present the general ideas consider a physical system of the form S : ȧ(t) = f (t,a(t),γ(t)), b(t) = g(t,a(t),γ(t)).

(2.1)

The RHS can be considered as any evolutionary model representing the dynamics, for example the Navier-Stokes equations. Although the ideas are easily extended to infinite dimension we present them for a finite dimensional case which is justified considering the usual approach of numerically resolving the equations by discretization in n space dimension. a ∈ R n represents the state variables, γ ∈ R m represents the inputs or the control applied to the flow, and b ∈ R p represents the output of observables. The RHS of (2.1) can be linearized to obtain a state-space form as

S LT I : E ȧ(t) = Aa(t)+ Bγ(t), b(t) = Ca(t)+ Dγ(t), (2.2)
where E ∈ R n×n , A ∈ R n×n , B ∈ R n×m , C ∈ R p×n and D ∈ R p×m . E need not necessarily be invertible. The dimension of the model usually corresponds to the dimension of the spatial discretization which can be very large, hence the principle of model reduction seeks a subspace V of dimension r ≪ n. The non-linear reduced order model corresponding to (2.1) can be written as

S : ˙â(t) = f (t, â(t),γ(t)), â ∈ R r , b(t) = ĝ(t, â(t),γ(t)), b ∈ R p , with r ≪ n, ( 2.3)
and the corresponding linear reduced model as

S LT I : Ê ȧ(t) = Â â(t)+ B γ(t), b(t) = Ĉ â(t)+ D γ(t). (2.4)
The reduced subspace has to have the following desirable properties when used as an approximation to the high fidelity model (2.1):

1. b-b < ǫ× γ ∀γ where ǫ is the tolerance.

2. The model need to preserve the stability of the attractor, and be passive, i.e. entropy conserving approximation to the high fidelity model S.

3. Must be efficient to reproduce the most important features of the high fidelity model.

To present the idea of projection, from now on we just consider the linear model S. We seek biorthogonal matrices V and W of size R n×r such that W T QV = I r where Q ∈ R n×n is a weight matrix. The Reduced Order Model is obtained by a projection of the state variable a on the matrix V such that a = V â and ŷ ≃ y. One defines the residue R to measure the accuracy of the projection, a ≃ V â as

R = EV ȧ(t)-AV â(t)-Bγ(t), (2.5a) b(t) = CV â(t)+ Dγ(t).
(2.5b)

The Petrov-Galerkin projection is obtained by requiring the residue orthogonal to the approximated space, i.e. W T QR = 0 r . The projection matrices of the system ŜLTI can now be obtained after a simple calculation as

 = W T QAV, B = W T QB, Ĉ = CV, D = D, Ê = W T QEV.
The well known Galerkin projection is obtained by taking V ≡ W. Several projection methods exist for the linear system (2.2). Examples are the Krylov subspace method for the controllability matrix which relies on the identification of the moments of transfer function [START_REF] Gugercin | A survey of model reduction by balanced truncation and some new results[END_REF], projection on the dominant modes of controllability and observability matrix [START_REF] Moore | Principal component analysis in linear systems: Controllability, observability, and model reduction[END_REF][START_REF] Rowley | Model reduction for fluids using balanced proper orthogonal decomposition[END_REF], projection of global stability modes [START_REF] Barbagallo | Closed-loop control of an open cavity flow using reduced order models[END_REF], POD based projection method on the subspace determined from the snapshots [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF] in which the flow fields are decomposed into energy ranked coherent structures. This is the method used in this work and is briefly described in the next section.

Proper orthogonal decomposition

The basic idea of a POD is to decompose the flow field into energy ranked coherent structures represented by mathematical modes the details of which can be found in [START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF][START_REF] Cordier | Proper Orthogonal Decomposition: an overview[END_REF]. Let H denote any Hilbert space with an inner product denoted by (•,•) and the norm induced by the inner product denoted by • . In the fluid dynamic context the space H consist of functions defined on some spatial domain Ω on which the fluid evolves, usually H = L 2 (Ω). Let q(x,t) be any vector of flow variable, where x and t ∈ [0,T] denote respectively the spatial variables and the time of the numerical simulation or experimental measurements. We also assume a discrete sampling of our data in time and space. We seek an expansion for q of the form

q(x,t) = q(x)+ n ∑ i=1 a i (t)φ i (x), (2.6)
where q denotes the average of the n flow snapshots. The spatial modes φ i (x) and the temporal coefficients a i (t) are determined by solving an eigenvalue problem involving the time correlation tensor that minimizes the average projection error:

E( q-P S q ). (2.7)
E denotes the averaging operator, for instance the ensemble average (E = 1 n ∑ n i=1 ) and P S denote the projection operator over the space S ⊂ H and dimension m. Note that the problem of minimizing E( q-P S q ) is equivalent to maximizing E( P S q 2 ), the "energy" of the orthogonal projection since by Pythogoraus theorem we have q 2 = q-P S q 2 + Pq 2 . Solving the optimization problem leads to an eigenvalue problem for the spatial modes φ given by Rφ = λφ, (

where R : H → H is the linear operator defined as

R = E(q k ⊗q * k ).
(2.9)

Here q * ∈ H * denotes the dual of q, given by q * (•) = (•,q) and ⊗ is the usual tensor product † . One can easily verify that R is self adjoint i.e. (R x,y)=(x,R y) and hence by the spectral theorem the eigenfunctions φ can be chosen to be orthonormal, The eigenvalues λ are then determined by taking the inner product of Eq. (2.8) with φ to obtain

λ = E |(q k ,φ)| 2 .
(2.10)

The eigenvalues λ in the above equation represents the average energy in the projection of the ensemble onto φ where energy is defined in sense of the induced norm. We also conclude that the R is positive semi-definite, and the eigenfunctions φ j which maximize E( P S q k 2 ) are the eigenfunctions corresponding to the largest eigenvalues. Also the eigenfunctions reproduce almost every member of the ensemble (except on a set of measure zero). One can also verify that the range of the operator R is contained in the span of the ensemble q k and that any eigenfunction φ can be written as linear combination of snapshots:

φ = n ∑ k=1 c k q k , (2.11)
where the coefficients c k ∈ R. If the average operator E over the snapshots q k is given as

E( f (q)) = n ∑ k=1 α k f (q k ), (2.12)
where α k > 0 satisfies ∑ n k=1 α k = 1 (usually giving equally weights to snapshots i.e. α k = 1/m). The eigenvalue problem (2.8) may be written in terms of the coefficients c k as

R c = λc, ( 2.13)
where c = (c 1 ,••• ,c m ) and R is an m×m correlation matrix with R ij = α i q j ,q i . This is the usual method of snapshots as described [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF]. Here we note that the direct problem involves solving an eigenvalue problem possibly on an infinite dimensional space H which may be very large, on the other hand the method of snapshots involve solving only an m dimensional eigenvalue problem and the method proves computationally efficient if m is small compared with the dimension of H. Having described the POD procedure we now proceed to build a ROM as explained in the next section.

ROM for cavity flows

The Reduced Order Model is obtained by performing a Galerkin projection of the governing dynamics (for instance the Navier-Stokes equations) onto the spatial modes determined from the POD procedure. Two ROM models for cavities has been proposed in [START_REF] Gloerfelt | Compressible Proper Orthogonal Decomposition/Galerkin reduced order model of self sustained oscillations in a cavity[END_REF], one using the full compressible equations and the other with an isentropic approximation. In this work we use the isentropic equations for compressible flows [START_REF] Rowley | Model reduction for compressible flows using POD and Galerkin projection[END_REF], which is a valid physical models for flows at moderate Mach number. Scaling the stream-wise and span-wise component of velocities u, v by the free stream velocity U ∞ , the local sound speed c by the ambient sound speed c ∞ , the lengths by the cavity depth D (see Fig. 2), and time by D/U ∞ , the equations are given by

u t +uu x +vu y + 1 M 2 2 γ-1 cc x = 1 Re (u xx +u yy ), v t +uv x +vv y + 1 M 2 2 γ-1 cc y = 1 Re (v xx +v yy ), c t +uc x +vc y + γ-1 2 c(u x +v y ) = 0,
where M = U ∞ /a ∞ is the Mach number and Re = U ∞ D/ν is the Reynolds number. If we denote q = (u,v,c) the vector of flow variables, the above equations can be recasted as

q = 1 Re L(q)+ 1 M 2 Q 1 (q,q)+Q 2 (q,q), (2.14)
with

L(q) =   u xx +u yy v xx +v yy 0   , Q 1 (q 1 ,q 2 ) = - 2 γ-1   c 1 c 2 x c 1 c 2 y 0   , Q 2 (q 1 ,q 2 ) = -      u 1 u 2 x +v 1 u 2 y u 1 v 2 x +v 1 v 2 y u 1 c 2 x +v 1 c 2 y + γ-1 2 c 1 (u 2 x +v 2 y )      .
To obtain the reduced order model by means of a Galerkin projection we define an inner product on the state space as

(q 1 ,q 2 ) Ω = Ω u 1 u 2 +v 1 v 2 + 2α γ-1 c 1 c 2 dΩ,
where α is a constant and γ is the ratio of specific heats. In this work we choose the value of α = 1, which gives the definition of stagnation enthalpy while calculating the norm.

The above definition ensures the stability of the origin of the attractor as demonstrated by [START_REF] Rowley | Model reduction for compressible flows using POD and Galerkin projection[END_REF]. We use expansion (2.6) and the definition of the inner product given above to perform the Galerkin projection of the isentropic equations onto the first n≪N POD spatial eigenfunctions, to obtain the Reduced Order Model as

ȧR i (t) = 1 Re C 1 i + 1 M 2 C 2 i +C 3 i + n ∑ j=1 1 Re L 1 ij + 1 M 2 L 2 ij + L 3 ij a R j (t) + n ∑ j,k=1 1 M 2 Q 1 ijk +Q 2 ijk a R j (t)a R k (t) = C i + n ∑ j=1 L ij a R j (t)+ n ∑ j,k=1 Q ijk a R j (t)a R k (t) = f i (C i , L i ,Q i y ,a R (t)) = f (y,a R (t)), (2.15)
where f i is a polynomial of degree 2 in a R and the coefficients of the projection are given by

C 1 i = (φ i , L(q)) Ω , L 1 ij = φ i , L(φ j ) Ω , Q 1 ijk = φ i ,Q 1 (φ j ,φ k ) Ω , C 2 i = (φ i ,Q 1 (q,q)) Ω , L 2 ij = φ i ,Q 1 (q,φ j )+Q 1 (φ j ,q) Ω , Q 2 ijk = φ i ,Q 2 (φ j ,φ k ) Ω , C 3 i = (φ i ,Q 2 (q,q)) Ω , L 3 ij = φ i ,Q 2 (q,φ j )+Q 2 (φ j ,q) Ω .

ROM for the actuated flow

Having discussed the ROM for the un-actuated case we now move on to introduce the effect of actuation. The advantage of the reduced order models can be fully exploited when they are capable of being used in control studies. We would also like to have a dynamical system where the actuation effect is naturally embedded. There has been recent attempts to take care of the introduction of actuation in [START_REF] Weller | Robust model identification of actuated vortex wakes[END_REF][START_REF] Caraballo | Control input separation methods for reduced-order model-based feedback flow control[END_REF] for a feedback control where snapshots from the controlled case are augmented with the un-actuated modes (baseline flow as it is usually called) to obtain a global representation for the controlled case. The regions where the control is introduced is explicitly identified while performing the projections, hence are termed as input separation methods. To explain further the flow domain Ω into two sub-regions, such that Ω = Ω ac +Ω unac where Ω ac , Ω unac represents the actuated and the un-actuated part of the domain respectively. The inner product for the corresponding domains is defined as

(•,•) = (•,•) Ω ac +(•,•) Ω unac ,
with the corresponding boundary condi- tion imposed for Ω ac . This procedure yields an autonomous system of ODE's of the form, on using the summation convention as.

ȧk = C k + L ik a k +Q ijk a i a j +h 1k γ+h 2ik a i γ+h 3k γ 2 , (2.16)
where γ is the actuation signal. Another procedure of introducing the actuation in the ROM is to identify the system coefficients of the actuated model based on a stochastic estimation technique, to give a ROM of the form (2.16) as in [START_REF] Caraballo | Control input separation methods for reduced-order model-based feedback flow control[END_REF]. The main disadvantages of the above methods is that they are quite cumbersome, as one has to explicitly identify the region of separation where the actuation is introduced, where there is a possibility of mismatch with the base line modes. Also the information of the input actuation is hidden in the model. One then looks for an alternative separation method to include the effect of actuation on the base line model as discussed in [START_REF] Kasnakoglu | Reduced order modeling, nonlinear analysis and control methods for flow control problems[END_REF], i.e. one seeks an expansion of the form

q(x,t) = q(x)+ n ∑ i=1 a i (t)φ i (x)+γ(t)ψ(x), (2.17)
where γ is the actuation and ψ is the spatial actuated coefficient and which should ideally satisfy the following conditions:

1. The actuation effect should be explicitly available in our ROM and should be able to reproduce the un-actuated dynamics in case the value of actuation tends to zero.

2. Most of the actuated dynamics are well captured by the un-actuated POD modes and hence the completion of the full actuated dynamics is obtained by an extension of the bases of the un-actuated case. This condition is an additional constraint to condition 1 in that we require a constraint on the dimension of our actuation, by specifying that the base-line modes capture most of the dynamics. This is the crux of the whole problem.

3. We further strengthen the above condition by imposing an optimality condition by requiring that the energy not captured by (2.17) is minimal.

Fig. 1 summaries the above condition. The conditions above leads to the following algorithm due to [START_REF] Kasnakoglu | Reduced order modeling, nonlinear analysis and control methods for flow control problems[END_REF].

Unactuated Space

∑ N i=1 a i φ i γψ Figure 1
: Schematic representation of the actuated expansion, the actuated mode is orthogonal to the unactuated space.

Algorithm 2.1

1. To start with, let the actuated snapshot sets be denoted as q ac k ,γ k m k=1 , where γ k = γ(t k ) is the value of the actuation, q ac k = q ac (x,t k ) and m is the number of actuated snapshots.

2. We subtract the mean q of the un-actuated base flow from the snapshot set. We define a new set of realizations by an innovation operator given as

qk = q ac k -P S q ac k = q ac k - n ∑ i=1 (q ac k ,φ i ) Ω φ i
to care of the part of the actuated mode which can be captured well by the un-actuated subspace.

3. We then wish to construct an orthogonal subspace to the un-actuated space to capture the effect of actuation. This is done by solving an L 2 minimization problem for the functional given by

J (ψ) = E qk -γ k ψ 2 ,
where E is any averaging operator.

4. The solution of the above minimization problem is given by

ψ = E[γ k qk ] E[γ 2 k ]
.

5.

The expansion for the flow field can now be written for the actuated case as

q ac (x,t) = q ac (x)+ n ∑ i=1 a ac i (t)φ i (x)+γ(t)ψ(x).
With this expansion the Galerkin projection of (2.14) gives

ȧac i (t) = C i + n ∑ j=1 L ij a ac j (t)+ n ∑ j,k=1
Q ijk a ac j (t)a ac k (t)+h 1i γ(t)

+ n ∑ j=1 h 2ij a ac j (t)γ(t)+h 3i γ 2 (t), (2.18)
where

h 1i = 1 Re (φ i , L(ψ)) Ω +(φ i ,Q(q,ψ)) Ω +(φ i ,Q(ψ,q)) Ω , h 2ij = φ i ,Q(φ j ,ψ) Ω + φ i ,Q(ψ),φ j ) Ω , h 3i = (φ i ,Q(ψ,ψ)) Ω , with Q(•,•) = 1 M 2 Q 1 (•,•)+Q 2 (•,•).
To study the effect of actuation on the unactuated baseline modes we observe from (2.18) the actuated temporal modes are given by

a ac i (t) = (φ i ,q ac -q ac -γ(t)ψ) Ω .
Here we make the assumption that the average of the mean flow in un-actuated and actuated cases are equal (q ac = q), since the value of actuation introduced is small. The term q ac -q ac can now be interpreted as a translation of the reduced order subspace. On noting this and the fact that the actuated mode is orthogonal to the un-actuated subspace i.e. ψ ⊥φ the term a ac i (t)-→0. This fact is important in that while calibrating the ROM, it is sufficient to just calibrate the un-actuated parts, as will be explained in the next section. Another point worth mentioning is that the time derivative of control vanishes due to the orthogonal property of the actuated modes. We also introduce, an error to take care of the difference in the average values while performing the Galerkin projections as

ε i = (δq ac ,φ i ) (φ i ,q ac ) . (2.19)

Calibration of the ROM

The reduced order model (2.18) has an intrinsic problem of converging to the wrong attractor when used to simulate the flow for long time periods. Hence there is a need to identify the coefficients of the dynamical system so as to minimis the error between the actual time coefficients a P i (t) from the POD and that obtained from the ROM a R i (t) using a suitable norm for the error. The reasons for the inaccurate behavior of the ROM can be attributed to the truncation of the POD bases where the dissipative scales of the higher POD modes are neglected. An analogous problem occurs in the Large Eddy Simulation (LES) of flows when there is lack of dissipation due to the smaller scales. Even including all the modes in Galerkin projection may still lead to the wrong attractor due to structural instability as has been demonstrated in [START_REF] Rempfer | On low-dimensional Galerkin models for fluid flow[END_REF]. Other problems may arise due to the contribution of pressure at the boundaries of the domain, which is usually neglected [START_REF] Noack | The need for a pressure-term representation in empirical Galerkin models of incompressible shear flows[END_REF]. The stability properties of the compressible POD-Galerkin approximation has been studied by [START_REF] Iollo | Stability properties of POD-Galerkin approximations for the compressible Navier-Stokes equations[END_REF].

Definition of errors

Three definitions of error have been defined in [START_REF] Cordier | Calibration of POD Reduced-Order models using Tikhonov regularization[END_REF]. The natural choice is in defining e 1 ( f ,t)=a P (t)-a R (t) but this leads to a nonlinear constrained problem (see Section 3.1.1). The other two choices being state calibration error e 2 without the dynamical constraint given by e 2 (y,t) = a P (t)-a P (0)-T 0 f(y,a P (τ))dτ and the flow calibration e 3 ( f ,t) = ȧP (t)f (a P (t)) where f is characterized by the coefficients of (2.15).

Nonlinear constrained problem

In this method we have to find the coefficients of the dynamical system such that the error e 1 is minimized under the constraints that the coefficients (2.15). We rather seek to minimize I 1 ( f ) = e 1 ( f ,t) 2 Λ T where • T is a time averaging operator over [0,T], for N t equally spaced elements on [0,T] we have

C i , L i ,Q i (i = 1,••• ,n) satisfy
f (t) T = 1 N t N t ∑ i=1 f (t i ), with t i = (i-1)∆t and ∆t = T N t -1 ,
. Λ is a norm on R N . For any positive definite matrix Λ ∈ R N×N the norm of any vector e ∈ R N is given by e 2 Λ = e T Λe.

The matrix Λ acts as a weight function giving importance to the specific POD modes, when Λ = I N it means that all the POD modes have the same importance in terms of the error. This leads us to a constrained optimization problem as in [START_REF] Bergmann | Optimal rotary control of the cylinder wake using POD Reduced Order Model[END_REF], in which we solve backward in time, an adjoint equation given by

ξi (t) = - n ∑ j=1 L ji ξ j (t)- n ∑ j,k=1 ξ j (t)(Q jik +Q jki )a R i (t)-2 a R i (t)-a P i (t) s = g i (C i , L i ,Q i y ,a R (t),ξ R (t))+s = g(y,a R ,ξ R (t))+s, (3.1)
where ξ i is the adjoint variable and s the source term. The terminal condition of the above equation is given by ξ i (T) = 0. The gradients of the cost functional I (1) with respect to the coefficients C i and L i of the ROM are

∇I 1 C i = T 0 ξ i (t)dt, ∇I 1 L i = T 0 ξ i (t)a R j (t)dt
and represent the sensitivity of the functional to the coefficients. The corresponding optimality system (state equations, adjoint equations, optimality conditions) is solved by a gradient based method to obtain the coefficients of the reduced order model.

A weighted approach to Tikhonov regularization

In this section we introduce a new method of calibration. As demonstrated in [START_REF] Couplet | Calibrated Reduced-Order POD-Galerkin system for fluid flow modelling[END_REF] min-

imization of I 3 ( f ) = e 3 ( f ,t) 2
T leads to the solution of a linear system for the coefficients y defining f . Unfortunately, this linear system is not well conditioned and lead to numerical divergence when the calibrated coefficients are used to integrate in time (2.15). For that purpose, [START_REF] Couplet | Calibrated Reduced-Order POD-Galerkin system for fluid flow modelling[END_REF] have introduced a new cost functional defined as a weighted sum of a term measuring the normalized error between the behavior of the model (2.15) with f and with the coefficients determined directly by Galerkin projection f GP and another term linked to the distance between f and f GP . The value of the weighting factor which represents the cost of calibration was chosen rather arbitrarily and hence was user dependent. In [START_REF] Cordier | Calibration of POD Reduced-Order models using Tikhonov regularization[END_REF], a Tikhonov regularization method was suggested to improve the conditioning of the linear system. The idea of this method is to seek a regularized solution y ρ as the minimizer of the following weighted functional

Φ ρ (y) = Ay-b 2 2 +ρ L(y-y 0 ) 2 2 ,
where the first term corresponds to the residual norm, and the second to a side constraint imposed on the solution. L represents the discrete approximation matrix of a differential operator of order d and y 0 a reference solution. The value ρ is chosen so as to compromise between the minimization of the norm of the residual for the linear system and the seminorm of the solution. Here, the regularization parameter ρ is determined by the classical L-curve method, as described in [START_REF] Cordier | Calibration of POD Reduced-Order models using Tikhonov regularization[END_REF]. We have the error defined for the minimization of the functional I i as

I i (y) = e i (y,t) 2 Λ T .
Usually the matrix Λ is chosen as Identity, which means that we give equal weights to all the modes in the definition of calibration. However this matrix can be chosen in suitable way so as to include the effect of mode selection in the definition. Two ways of defining the weights can be proposed:

1. We consider that the main interest is in modelling the effect of the energetic structures and hence the eigen-spectra themselves serve as a measure of the relative importance of the modes, which is the most natural choice of the weights for the definition of error.

2. The error can be based on an overall sensitivity of the model with respect to a cost functional.

The weight matrix Λ for the definition of error for the first case can be simply written as a diagonal matrix:

Λ ii = σ i max σ i for i = {1,••• , N},
where σ is the singular value obtained as a solution of the eigenvalue problem of the time correlation matrix. For the second case consider the state equations

ȧR = f(y,a R ). (3.2)
Variation of any convex cost functional J with respect to the state variables

a R = {a R i } N i=1
gives the adjoint equation of (3.2) as

ξR = g(y,a R ,ξ R ), (3.3)
where ξ R (t)={ξ R i } N i=1 is the adjoint variable. The overall sensitivity of the cost functional J with respect to a R is obtained as

S = dJ da R = a R (t)ξ R (t) , (3.4)
where • is any time averaging operator. We can then define the weight matrix Λ with respect to the sensitivity as

Λ ii = S i max S i for i = {1,••• , N}.
In this study we have taken the cost functional J based on the energy of the temporal modes as

J = 1 2 T 0 N ∑ i=1 (a R i (t)) 2 dt. (3.5)
The above functional is minimized subject to the constraint (3.2), by the method of Lagrange multipliers with the adjoint state equation of (3.3) given by

ξ(t) = -a i (t)- N ∑ i=1 L ij + N ∑ i=1 Q jik +Q jki a k (t) ξ j (t). (3.6)

Application of POD ROM for cavity flows

The typical cavity flow configuration used in this study is shown in Fig. 2. The cavity is of an L e /D ratio of 2. The flow is initialized by a laminar boundary layer so as to have a thickness of δ/D = 0.28 at the leading edge of the cavity. The Reynolds number of the flow based on the cavity depth is 1500 and the flow Mach number is 0.6. The equations of Navier-Stokes are discretize using a 4 th order accurate scheme in time and space. Snapshots are taken once the flow has stabilized and 56 snapshots are uniformly sampled which corresponds to 1 period of the flow oscillation (2.8 in non dimensional time) corresponding to the first Rossiter mode [START_REF] Delprat | Rossiter's formula: A simple spectral model for a complex amplitude modulation process?[END_REF]. Fig. 3 demonstrates a degenerate eigen-spectrum showing eigenvalues which occur in pairs. Also the first 4 eigenmodes capture around 98.5% of the total fluctuation energy as shown by the Relative Information Content (RIC) defined as RIC(i) = ∑ i j=1 σ j /∑ M j=1 σ j . The temporal coefficients shown in Fig. 4 and demonstrates a phase difference of π/2 between the paired modes. The representation of the first 4 spatial POD modes for the vorticity is shown in Fig. 5.

Although the modes occur in pairs and their values are distinct the representation is topologically equivalent. We note that the cascade of the energy in the POD representation in terms of the size of the eddies represented. The vorticity being a hydrodynamic phenomenon represents the low frequency dynamics of the flow. In Fig. 6 a comparison of the two weight matrix used for the purpose of calibration is given. We note a similar behavior between the two weights as an energy based cost functional is used in both the cases † . The sensitivity based weight is different for mode 3, and this allows us to speculate that in the case where the ROM is not built using the POD bases, but on a different criteria, such as the observability as described in [START_REF] Jordan | Identifying noisy and quiet modes in a jet[END_REF], the above analysis may yield a sensitivity weight for the lower modes. The time evolution and phase plot for the case of the calibrated model is shown in † Although the sense of energy is different in two cases, it is a energy of projection in the case of POD and the energy contained of the flow dynamics, for the sensitivity based analysis. a small errors for the most energetic modes. We are unable to explain the behavior of the error at the certain higher modes which seem pretty large. Since we are dealing with periodic flow it is natural to expect the temporal dynamics be valid for a time longer than the period of snapshot acquisition. As shown by [START_REF] Sirisup | A spectral viscosity method for correcting the long-term behavior of POD model[END_REF] although the system is initialized with correct state the solution may drift away for a long period of integration even after calibration. Fig. 12 shows that the calibrated model predicts the dynamics for a nondimensional time 11, which corresponds to around 4 cycle of the flow period, after which it diverges rapidly. This is different from the results where the calibrated model is shown to work for a much longer period. The reason is that the neglected modes which contribute to the regularization of the system is not modelled during calibration. The validity of ROM for a longer time of integration is still open, where the coefficients need to be given suitable weights to model the neglected term. While performing control, since the time period where the control law is applicable is much larger than the period of validity of the model, it is imperative that we calibrate more than one period of the flow as demonstrated in the next section.

a 1 (t) a 2 (t) a 3 (t) a 4 (t) t a k (t)

Feedback design

In this section we propose a feedback law based on an observer design, obtained from the pressure measurements inside the cavity. The plant model utilized in this study is shown in Fig. 13. It consists of the system i.e. linearized ROM and a feedback controller. When all the states of the system are difficult to measure we build an observer to estimate the state while only measuring the output. Regarding the use of ROM for optimal control studies one can refer [START_REF] Bergmann | Optimal rotary control of the cylinder wake using POD Reduced Order Model[END_REF][START_REF] Graham | Optimal Control of Vortex Shedding Using Low Order Models. Part 1. Open-Loop Model Development[END_REF][START_REF] Graham | Optimal Control of Vortex Shedding Using Low Order Models. Part 2: Model-based control[END_REF] where the optimal control has been applied for the control of wake flows behind a cylinder. Since the ROM is designed for a particular flow condition, [START_REF] Fahl | Trust-region methods for flow control based on Reduced Order Modeling[END_REF] propose a method where the model is refined during the optimization process, to take care of the validity of the model to the variation of control parameters. This method has been applied for the optimal control of wake flows by [START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by Trust-Region methods and POD Reduced Order Models[END_REF]. Control of cavity flows, using experimental data to build the ROM has been performed in [START_REF] Samimy | Feedback Control of Subsonic Cavity Flows Using Reduced-order Models[END_REF], but without any calibration of the temporal dynamics. In this work, the feedback law is determined using the linear part of the dynamical system, neglecting the influence of the quadratic terms. When information of the full state a is not available, an observer C is construct from the measured output y. This is true in case of experimental data, where wall sensors to measure the instantaneous pressure are used. Prediction of flow field from the observed states has been studied in the past by performing a Linear Stochastic Estimation (LSE) by [START_REF] Adrian | Stochastic estimation of organized turbulent structure: homogeneous shear flow[END_REF][START_REF] Bonnet | Stochastic estimation and proper orthogonal decomposition: Complementary techniques for identifying structures[END_REF], where sensor signals are correlated with flow data base information. A higher order correlation is also demonstrated in [START_REF] Buffoni | A non-linear observer for unsteady three-dimensional flows[END_REF], where a Quadratic Stochastic Estimation (QSE) is performed. The main drawback of the LSE and QSE is that they take into account the correlations at a particular time only without considering the time history of the previous measurements, hence they are termed as 'static' estimators.

In an experimental setup this poses a problem due to the measurement noise and also on the limitations of the number of sensors that can be placed. A combination of the isentropic model with the data obtained from experiments has been studied in [START_REF] Samimy | Feedback Control of Subsonic Cavity Flows Using Reduced-order Models[END_REF]. The technique of LSE which is used in this work for observer design is briefly explained in the next section.

Linear Stochastic Estimation (LSE)

The LSE technique enables us to determine the observer matrix C, as shown in Fig. 13, given a limited number of pressure measurements obtained from sensors along the walls of the cavity. For each sensor i ∈ [1 N s ], we can estimate the pressure pi (t) coefficients by a series expansion of the discrete temporal POD coefficients a j (t) as

pi (t) = n ∑ j=1 Cij a j (t)+ n ∑ j=1 n ∑ k=1 Dijk a j (t)a k (t)+••• . (5.1)
The expansion can be truncated above the linear term or the quadratic term corresponding to the LSE or QSE problem. The coefficients Cij and Dijk are determined so as to minimize the mean square error between the estimated pressure signal pi (t) and the one obtained from measurement p i (t)

e i = T 0 [ pi (t)-p i (t)] 2 dt = [ pi (t)-p i (t)] 2 . ( 5.2)
Since the number of snapshots is much larger than pressure sensors, the above system is overdetermined and the coefficients are determined using the method of least squares by solving an overdetermined linear system given by

de i d Cij = 0 and de i d Dijk = 0. (5.3)
On noting the orthogonality of the temporal modes a j ,a k = σ i δ jk , the pressure field can be simply expanded in terms of the temporal modes as

pi (t) = N ∑ i Cij a j (t), where Cij = p i ,a j a j ,a j . ( 5.4)
In this study N s = 4 sensors along the wall are used to measure the pressure. It is also worthwhile mentioning that solving the system of equation given by (5.3) leads to an ill conditioned problem. A luxury afforded by numerical simulations is that we can modify the condition number by changing the location of the pressure sensors, however for data obtained from experimental simulations where the sensor positions are fixed we can use the method of Tikhonov regularization mentioned in Section 3.2. The observer dynamics for the construction of the pressure signals, issued from (5.4) is shown in Fig. 14, and provides a good estimation of the wall pressure, due to the orthogonality property of the eigenmodes.

Linearization of the plant

Once through the observer design, the plant (2.18) is built by linearizing the ROM around a suitable equilibrium point a 0 . To determine the influence of various terms of the ROM while linearizing, we perform a sensitivity analysis. Sensitivity studies help in analyzing the variation of a given quantity (generally given by a functional) with respect to a parameter of the model. In this study a sensitivity analysis is performed to identify the terms to be neglected in the linearization of the model. We seek the sensitivity of the time-independent coefficients h 1i , h 2ij and h 3i of the actuated model, with respect to an energy based functional given by (3.5). A variation of the augmented Lagrangian functional based on an adjoint vector ξ and the state equation (2.18) yield the sensitivities:

∂J ∂h 1i = ξ i (t),γ(t) , ∂J ∂h 2ij = a j (t),ξ i (t)γ(t) , ∂J ∂h 3i = ξ i (t),γ 2 (t) , (5.5)
and

∂J ∂γ = n ∑ i h 1i ξ i (t)+ n ∑ i n ∑ j h 2ij ξ i (t)a j (t)+2 n ∑ i h 3i ξ i (t)γ(t),
where the adjoint equation is given by:

ξi (t) = g i (C i , L i ,Q i ,a(t),ξ(t))+ n ∑ j h 2ji ξ j (t)γ(t)-a i (t). (5.6)
The sensitivities are plotted in Fig. 15 and shows that the constant term of the actuated model h 1i has a higher sensitivity when compared to the quadratic terms h 3i by an order of magnitude and hence plays an important role in linearization.

It is also important to note that the first two modes which drives the actuation dynamics do not contribute significantly to the sensitivity of the terms h 1i and h 3i as compared to h 2ij , this signifies that more care has to be taken while linearizing these terms.

To determine the equilibrium point a 0 to be used for linearization we solve a nonlinear algebraic matrix equation for the unactuated case γ = 0. Eq. (2.18) can be recasted for the

j S ij i = 1 i = 2 i = 3 i = 4 i = 5 i = 6 (b)
f (a) := C+ La+    a T Q 1 a . . . a T Q N a   . (5.7)
In this study a Newton's iterative method is used to calculate the equilibrium point a 0 using the following algorithm:

a k+1 = a k -J -1 (a k ) f (a k ), (5.8)
where J(a) is the Jacobian matrix: J = ∂ f (a) ∂a . Also note that the solution of the Newton method (Eq. (5.8)) depends on initial conditions and is not unique. To find the equilibrium point, we have chosen the initial conditions of the dynamical system at t = 0. On using a linearization at the origin, the plant model can be written in the form (refer to appendix for details) ȧ = Ã ã+ Bγ, p = C ã.

(5.9)

Controller

After linearizing the plant, one can perform the feedback control. The full plant, control and estimator model are shown in Fig. 13. The model for the controller is given by the equations for a LQG design as ȧ = Ãa+ Bγ+w 1 , y = Ca+w 2 .

(5.10)

where for the sake of clarity we denote a = ã and y = p. The controller noise w 1 and the measurement noise w 2 are assumed to be uncorrelated, zero mean, white Gaussian processes with modelled spectral densities Q 1 and Q 2 . In this work the matrices Q 1 and Q 2 matrices are set to identity. We also introduce the noise perturbation vector w = [w 1 w 2 ].

The response of the state can be separated into three components, the natural response a 0 ( Ã,a(0);t) which depends on the eigenvalues of the main plant matrix à and on the initial condition a(0), the response to a given control γ(t) denoted by a γ ( Ã, B,γ;t) and the response to the noise w 1 given by a w ( Ã,w 1 ;t Assuming a full information state, the feedback law is given by γ = -K c a, (5.12) where the matrix K c is chosen so as to stabilize the plant and minimize the quadratic cost functional given by:

J = T 0 a T Wa+ℓ 2 γ 2 dt, (5.13)
where W is a weight matrix (here it is identity) and ℓ represents the cost of the control. By a simple variational approach, it can be demonstrated that the gain vector K c is given by

K c = 1 ℓ 2 BT X, (5.14)
where X is the solution of the steady Riccati equation to get infinite horizon stabilization:

ÃT X +X Ã-1 ℓ 2 X B BT X +W = 0.

(5.15)

Observer

In many cases the state a is known only from a limited number of output measurements y, and must be built from an observation ŷ which is an approximation of the real value of the measurement y. Denoting â the approximation of the state, the observer equations can be written as: ȧ = Ã â+ Bγ+K o (yŷ), ŷ = C â.

(5.16)

The observer gain matrix K o is obtained by solving the Ricatti equation for Y and the observer matrix C and is given by

ỸA T + ÃY -Y CT CY +Q 2 = 0.
(5.17)

As described in [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF] the observer problem is dual of the controller problem described in (5.15). The estimation gain K o is obtained finally as

K o = Y CT .
(5.18)

Stabilization of the full system

The full system for the controller and observer, defined in terms of the state and the estimator error e a = a-â, with the state denoted by x = [a, e a ], can be written as:

ẋ = A t x+ D t w, ( 5.19)
where

A t = Ã-BK c BK c (0) dim à Ã-K o C , D t = I d 0 I d -K o .
To get the system the feedback is applied on the observer state, â, as γ = -K c â. The eigenvalue of the full matrix A t are the same as the eigenvalue of the matrices Ã-BK c and Ã-K o C, and are chosen stable by the Riccati equations. In the absence of any noise, the state errors goes to zero in the infinite horizon, hence the estimation state â goes to the real state a, and in the same way, the observer error e y = y-ŷ = Ce a goes to zero. In presence of Gaussian noise w, errors e a and e y remains under a given threshold, and are functions of the correlation matrices Q 1 and Q 2 as stated in [START_REF] Bagheri | Input-output analysis and control design applied to a linear model of spatially developing flows[END_REF]. Following the notation of (5.11), the solution can written as:

x(t) = x 0 (A t ,x(0);t)+ x w (A t ,D t w;t),
where ||x 0 || → 0 when t → ∞.

Application to cavity flow

The control has been designed for different values of the cost ℓ ∈ [0.01,5000], and the result was found to be roughly independent of the value of ℓ. Finally the value of ℓ = 1 was chosen for this study. In Fig. 16 the eigenvalues of the open loop matrix à and the closed loop matrix Ã-BK c are represented. The feedback control clearly mirrors the unstable eigenvalues, with dominant eigenvalues of the closed loop λ max =-0.038±i2.26. For values of ℓ>5000, it seems numerically difficult to use the Matlab system toolbox for the solution of the Ricatti equation, to find an accurate and stable eigenvalue of Ã-BK c .

The value of the feedback gain ||K c ||≈950 seems quite large comparing to some usual experimental or numerical data. When ℓ → ∞ the cost of the control goes to zero and the system can not be stabilized. In the same manner the norm of the gain goes to zero, similarly when ℓ is small, the efficiency increases due to the large amplitude of the control. However in practical cases the value of the control that can be fed into the actuator is limited by a range of the input signal as a large gain may damage the actuator [START_REF] Samimy | Feedback Control of Subsonic Cavity Flows Using Reduced-order Models[END_REF].

To scale down the large gain and to decrease the amplitude of the control, a new parameter α is introduced such that γ α (t) = αK c with α ∈ [0,1.1]. In Fig. 17 the real part of the largest eigenvalue of Ã-α BK c is plotted with respect to α. The value of α = 1 corresponds to the optimal control, while α = 0 corresponds to the uncontrolled case. The value of α = 0.5 corresponds to the stability limit of the controlled system, and is the result of mirroring the unstable eigenvalues about the imaginary axis. It means that to stabilize the plant, the norm of the gain control vector should be at least equal to 475, which is a large value for physical application.

The observer response is directly related to the largest eigenvalue of the matrix Ã-K o C. For our case we found a real value of -0.45 which is sufficient to ensure the decrease of the estimation errors e a with time. To determine the optimal control input γ(t) opt , we solve the full system (Eq. (5.19)) and we use the feedback law γ(t) opt = -K c (a-e a ). To remind of the overall procedure, a simulation without control is performed, and the unactuated ROM model is constructed using a Galerkin projection. The model is further calibrated as explained in Section 3.2 to avoid the divergence of solution when integrated in time. To introduce the effect of actuation, the high fidelity model is forced with a chirp excitation of the form

γ g (t) = A 1 sin(2πSt 1 t)×sin(2πSt 2 t-A 2 sin(2πSt 3 t)) (5.20)
at the upstream corner of the cavity where

A 1 =0.1, A 2 =27, f 1 =1/60, f 2 =2/3, f 3 =1/30.
The form of the excitation is shown in Fig. 18.

The spectra of the chirp is represented Fig. 19 depicts a large spectral band in the range 0.2-1.5 with a concentration around the value of 1.5. 300 snapshots are sampled in the time interval t ∈ [0,T] with T = 15, corresponding to half the period of excitation of the full chirp. The actuated mode and the plant model are constructed with the help of these realizations. The plant model with the input γ g (t) is simulated. The feedback law γ(t) = -K c a g is built from the response of the plant a g (t):

γ(t) = -K c e Ãt a(0)+ t 0 e Ã(t-τ) B γ g (τ) dτ .
(5.21)

The new updated control is found to be of large magnitude max γ(t) = G max ≈ 250. The evolution of the final optimal control in comparison with the initial solution is as shown in Fig. 20.

A curve fitting approach demonstrates that γ(t) is given by

γ(t) = G max sin(ωt+ ϕ), ( 5.22)
where ω = 2.27, corresponding to a non-dimensional frequency f = 0.3617 and a phase ϕ = -1.406. The frequency is in good approximation to the imaginary part of the largest eigenvalue of the controlled plant (matrix Ã-BK c ). As is seen from (5.21), the gain vector K c acts as a filter and signals out the frequency which correspond to the shear layer instability of the second Rossiter mode.

To avoid numerical divergence of the numerical solution of the Navier-Stokes equations with the introduction of control, the amplitude (gain) G 0 of any control is set to a maximum value of G 0 = 0.08. So the updated of control is of the same form as γ(t) but with an amplitude G 0 . The power spectral density for the vertical component of velocity at a point in the shear layer is compared with the simulation from the unactuated case as shown in Fig. 21 and shows a decrease in amplitude of the second Rossiter mode.

However there is an increase in amplitude at lower frequency due to the effect of actuation as we add energy to the flow. The additional energy goes in furthering the impingement length of the shear layer on the downstream cavity edge, thereby reducing the strength of the pressure signals propagating upstream. The increase in the lower frequency does not correspond to any acoustic modes. The full system response corresponding to (5.19) is shown in Fig. 22 showing an asymptotic stability.

The errors are plotted in Fig. 23, also the response of the observer model is shown in Fig. 24. The response of the observer shows that the state for the pressure goes to zero asymptotically, but we find large values of gain for the short period of the DNS simulation, where we implement the control, which necissates the scaling of the final control introduced.

Conclusion

In this work different tools of reduced order modelling for the simulation and control fluid systems with instabilities has been presented with an application to cavity flow. The methods presented in this work for the case of a cavity are as well applicable to other systems with instabilities. Numerical simulations were performed to obtain the variables of the isentropic compressible model, which was used to construct the POD Galerkin based reduced order model. Results for the cavity flow configuration has been presented. The extension of the ROM to the actuated case is obtained by solving an optimal problem for capturing the dynamics not captured by the unactuated modes. A new error to measure the effect of actuation on the mean flow has been proposed.

Calibration of the ROM is performed using a Tikhonov based regularization obtained from a sensitivity analysis to obtain an accurate representation of the dynamics. It is also shown that the model diverges for a long time of integration even after calibration, and it is necessary to calibrate over a long period of time, especially when performing control studies.

A feedback control law based on the estimation of the observer dynamics has been presented. The observer matrix is constructed using a linear stochastic estimation. A sensitivity study of the actuated dynamics has been performed to determine the relevant terms in the linearization of the model. Finally an LQG based feedback law is obtained to find the optimal solution. The optimal control law is shown to have a single frequency, corresponding to the second Rossiter mode of the cavity flow. The LQG control provides stabilization over infinite time horizon with a huge control amplitude, unrealistic to be utilized in the Direct Numerical Simulation which are solved over a short time interval. To resolve this problem, the control amplitude has been scaled down before introducing in Direct Numerical Simulation. The control provides a decrease in spectra corresponding to the second Rossiter mode.

As a future perspective it is observed that the POD modes which is an optimal representation of the energy fields is not suitable to resolve the far field acoustic, one of the extension of POD to include the effect of far field acoustic noise is the Most Observable Decomposition, as demonstrated in [START_REF] Jordan | Identifying noisy and quiet modes in a jet[END_REF] which represent a good choice of bases to perform control studies for acoustics.

( ã+a 0 ) T Q 1 ( ã+a 0 )

. . . where

( ã+a 0 ) T Q N ( ã+a 0 )   +h1γ+    (h 1 21j γ) T ( ã+a 0 ) . . . (h N 2Nj γ) T ( ã+a 0 )    = C+ La 0 +    a T 0 Q 1 a 0 . . . a T 0 Q N a 0   + L ã+    ãT Q 1 a 0 +a T 0 Q 1 ã . . . ãT Q N a 0 +a T 0 Q N ã  +    ãT Q 1 ã . . . ãT Q 1 ã  
L = L+    ãT 0 (Q 1 +Q 1 ) T ) . . . ãT 0 (Q N +Q N ) T )   , h1 = h 1 +    (h 1 21j ) T . . . (h N 2Nj ) T   a0.
The observer model can be written with a perturbation p = p-p 0 as: p = p-p 0 = Ca-Ca 0 = C ã.

(A.4)

The plant model (5.9) is obtained by taking à = L and B = h1 .

Figure 2 :

 2 Figure 2: Cavity configuration and computational domain.

Figure 3 :

 3 Figure 3: Real part of the eigenvalues and Relative Information Content (RIC).

6 Figure 4 : 4 Figure 5 :

 6445 Figure 4: First 6 POD temporal coefficients.

Figure 6 :Figure 7 :

 67 Figure 6: Comparison of the weight matrices used in definition of errors. (a) Weight matrix based on the eigen spectra of the POD, i.e. σ k maxσ k ; (b) Weight matrix based on the sensitivity analysis for the cost functional based on energy, i.e. S k maxS k .

Figure 8 :

 8 Figure8: Phase portrait comparison of the first 6 modes, with the POD coefficients:ROM calibrated (solid line), reference POD dynamics (o). POD ROM is calibrated using the I(3) minimization for all coefficients. The Tikhonov regularization is applied for (L = 0, y 0 = 0), and the weight matrix Λ based on the overall sensitivity.

Fig. 7 and 1 -

 71 Fig.7and Fig.8respectively. For the case of the introduction of the actuation mode, snapshots are obtained from the DNS by introducing an actuation of the form Asin(ωt) just before the leading edge of the cavity (x ∈ [-0.15;-0.05] and y = 0), where the flow is more sensitive to actuation. The spatial modes for the components of velocity are shown in Fig.9and exhibit a local behavior capturing the effect of actuation which is the salient feature of the L 2 optimization method. The difference in average values δq ac for the stream-wise and span-wise component of velocity u, v is shown in Fig.10. The value of the error in the projection ε i as described in Eq. (2.19) is shown in Fig.11, showing

Figure 9 :

 9 Figure 9: u and v velocity components of the actuation mode ψ corresponding to an actuation defined by v wall = 0.2 sin(0.4t).

Figure 10 :

 10 Figure 10: u and v of the difference in the average value δq ac between the actuated and the unactuated case.

Figure 11 :

 11 Figure 11: Average projection error ε k .

Figure 12 :

 12 Figure 12: Evolution of first 4 temporal modes of the state of the ROM a k (t) for a time integration greater than the period of validity of the model. The model diverges after a time t = 11.

Figure 13 :

 13 Figure 13: Plant and observer for feedback control. C: observer, L: Estimator, K: Feedback.

Figure 14 :

 14 Figure 14: LSE reconstruction of the pressure signals. The pressure value from the sensor are denoted by line, and the coefficients obtained from LSE denoted by • (one over 3• are plotted). The sensors are located at positions (0,D/4), (0,3D/4), (L,D/4), (L,3D/4), (D/2,0), (3D/4,0) along the cavity walls and in the shear layer. Unactuated case.

Figure 15 :

 15 Figure 15: Sensitivity of the actuated coefficients. Constant: S i = ∂J ∂h 1i ; quadratic: S i = ∂J ∂h 3i (rescaled by 10);

Figure 16 :

 16 Figure 16: Circle: eigenvalues of the plant; (o): eigenvalues of Ã; (+): eigenvalues of Ã-BK c .

Figure 17 :

 17 Figure 17: Re(λ max ) with respect to α.

1 tFigure 18 :Figure 19 :

 11819 Figure 18: Chirp excitation γ e : forcing the DNS.

Figure 20 :

 20 Figure 20: Comparison of the initial forcing γ g (t) and of the final control law γ(t).

Figure 21 :

 21 Figure 21: Comparison of Power Spectral Density (PSD) of the original solution and the optimal solution of the normal component of velocity sampled at y = 0 and x = 1.8D: Actuated (dashed line); Unactuated (continuous line).

Figure 23 :

 23 Figure 23: Error on the components of the state e x .

6 Figure 24 :

 624 Figure 24: Output prediction y, the pressure fluctuations for the 6 modes. Mode 1 and 2 have been rescaled by 100.

2 )

 2 Since a 0 is the equilibrium point of the unactuated model it satisfies (5.7). The terms in (A.2) can be rearranged as the function of the perturbation of the state ã as ã

  Figure22: System response of the full system for the state components a 1 and a 2 .

		0,16
										e x 1 (t)
		0,12								e x2 (t)
										e x3 (t)
	x i	0,08								e x 4 (t) e x5 (t)
	e									e x6 (t)
		0,04
		0
	-0,04	0	1	2	3	t	4	5	6	7

a 1 (t)

† We have (u⊗v * )(ψ) = u(ψ,v), ∀u,v,ψ ∈ H.

Acknowledgments

This research project has been supported by a Marie Curie Early Stage Research Training Fellowship of the European Community's Sixth Framework Programme under contract number MEST CT 2005 020301 [41]. We would also like to thank CALMIP center in Toulouse for providing the computational resources. This article has been written in memory of Professor Pierre Comte from PPrime lab, Poitiers, who has made available his DNS code and who tragically passed away in September 2011.

Appendix

To obtain the linearization equation (2.18) can be written as

In this system the quadratic term γ 2 is neglected. As the equilibrium point is unaltered in the presence of a feedback, we introduce a perturbation ã = a-a 0 to write the above system with the perturbation and the equilibrium state: