Regularization method for calibrated POD reduced-order models ˆ

In this work we present a regularization method to improve the accuracy of reduced-order models based on Proper Orthogonal Decomposition. The bench mark configuration retained corresponds to a case of relatively simple dynamics: a two-dimensional flow around a cylinder for a Reynolds number of 200. Finally, we show for this flow configuration that this procedure is efficient in term of reduction of errors.


Introduction
A reduced-order modeling strategy is crucial to achieve optimization and model-based control in a wide class of complex flow configurations [2].Indeed, if these tasks can always be reduced to constrained optimization problems, their numerical resolution is generally so expensive in terms of CPU and memory storage that strongly limits the applications to simplified configurations [4].The development of flow control for three-dimensional turbulent flows encountered in many engineering applications or real-life systems is thus conditioned by the use of Reduced-Order Models (ROMs).The objective of these models is to capture the essence of the physics of the system while reducing the costs associated to the solution of non-linear state equations.In fluid mechanics, ROMs are mostly derived by Galerkin projection of first-principles equations, the Navier-Stokes equations, onto the Proper Orthogonal Decomposition (POD) modes.These POD ROMs are known to be relatively fragile when used for control design.This is partly due to the lack of adaptation of the POD modes when the flow conditions are changing but this is also strongly linked to the intrinsic truncation effect introduced by keeping only the most energetic modes in the model [9].For these reasons, it is necessary to develop specific numerical methods for improving the accuracy of these reduced-order models, the goal being to determine part or all of their coefficients.
This paper is a sequel of [8] in which we have compared different calibration methods found in the literature and treated the ill-conditioned system by developing a Tikhonov-based regularization technique that was extending the ideas of [1].Here, we introduce a new regularized optimization problem aiming at minimizing a weighted average of a normalized error, and a term measuring the variation of the coefficients from their value obtained by Galerkin projection.Moreover, by analogy of the L-curve method [17], we calculate the optimal value of the regularization parameter by determining an expression for the inflection point of the curve.Finally, we bench mark our method with respect to the Tikhonov-based calibration technique for a two-dimensional wake flow around a cylinder.
This manuscript is organized as follows.Section 2 presents the Proper Orthogonal Decomposition and describes how to derive a POD-based reducedorder model for the Navier-Stokes equations.In Section 3.1, we define the calibration as a minimization problem and introduce a regularized framework in order to reduce the ill-conditioning of the linear system.Then, in Section 3.2, we introduce the concept of filter factors for explaining the ill-conditioning and detail how to calculate optimally the regularization parameter by determining the inflection point of the L-curve.In Section 4, we then describe the bench mark flow configuration used in this study and compare our approach to the Tikhonov-based regularization method introduced in [8].Finally, we conclude and mention some perspectives (Section 5).

POD-based Reduced-Order Model
Reduced-order modeling is a way to remove un-necessary redundancies present in a physical system to describe with a sufficient level of accuracy the dynamics of a large-scale dynamical system with a low-order dynamical system.There is a large variety of ROMs in the literature depending on the targeted applications and on the flow physics.Here, we are focusing on those based on the Proper Orthogonal Decomposition.

Proper Orthogonal Decomposition
Proper Orthogonal Decomposition (POD) [18] is a powerful and efficient method introduced in Turbulence for extracting from a high-dimensional data set obtained numerically or experimentally a low-dimensional approximation that captures much of the phenomena of interest.The mathematical framework at the heart of POD is developed in detail in [6] and two typical applications of POD to fluid flows is described in [7].Here, for self-consistency of the paper, we summarize the most important properties of POD for reduced-order modeling.
Let U = {u(x, t m ) = u m } m=1,...,Ns be a set of N s snapshots taken over a time interval [0, T ], with x ∈ Ω, the spatial domain of interest.The main idea of POD is to determine a subspace S of dimension N POD N s , such that the error E ( u − P S u H ) is minimized.Here, || • || H denotes a norm induced by the inner product (•, •) H defined on a Hilbert space H, P S is the orthogonal projection onto the subspace S, and E(•) is an average operator over the snapshots, for instance an ensemble average (E (u) = 1 Ns Ns m=1 u m ).Note that minimizing E u − P S u 2 H is equivalent to maximizing E P S u 2 H , since u 2 H = u − P S u 2 H + P S u 2 H for any orthogonal projection P S .For our purposes, the space H corresponds to L 2 (Ω), the space of square-integrable functions on Ω.The POD procedure is then equivalent to minimize the expression: , where {Φ j } j=1,...,NPOD is a basis for the subspace S and {a P j } j=1,...,NPOD refer to the temporal coefficients corresponding to the POD expansion (as indicated by the superscript P ).
Solving this optimization problem leads [6] to the eigenvalue problem: with R = E (u ⊗ u * ), the two-point spatial correlation tensor.Here, ⊗ denotes the dyadic product bewteen two vectors u and u * where the superscript * indicates complex conjugate.It can be shown [6] that the operator R is linear, self-adjoint and positive semi-definite on H.If we further assume that R is compact, then there exists a countable set of decreasing non-negative eigenvalues λ j (j = 1, • • • , N POD ) associated to Φ j .These eigenfuntions may be chosen to be orthonormal for the inner product defined on L 2 (Ω), i.e.
with • the notation of the standard Euclidean inner product and δ jk the Kronecker delta.The POD modes are ranked according to the magnitude of their eigenvalue, with λ 1 equal to the largest eigenvalue and λ NPOD equal to the smallest one.
When the input data come from numerical simulations, it is in general much more efficient in terms of computational cost to use the method of snapshots introduced by [23] than the direct method defined in (2.1).The method of snapshots consists of writing the POD eigenfunctions Φ j as linear combinations of the snapshots.After inserting this expression in (2.1), a new eigenvalue problem for a P j = a P j (t 1 ), • • • , a P j (t Ns ) T is obtained: where C is the two-point temporal correlation tensor.This tensor is defined as

Galerkin projection
For deriving the reduced-order model, the next step is to project the governing physical equations onto the POD basis.For incompressible flows, the motion of the fluid is described by the incompressibility condition (∇ • u = 0) and the Navier-Stokes equations are given by: In these equations, all the variables (the velocity vector u and the pressure p) are assumed to be non-dimensional and Re is the Reynolds number.The expansion of u over the POD modes Φ j writes: where u m = E(u).By substituting (2.3) into the Navier-Stokes equations (2.2), and using the eigenfunctions Φ i (i = 1, . . ., N Gal N POD ) as test functions for a variational formulation, we obtain: N Gal corresponds to the number of Galerkin modes retained in the reducedorder model.This number is determined based on the energetic content of the first N Gal POD modes as defined by the ratio GP ijk and the pressure term P i depend explicitly on Φ and u m and are given in Appendix.If the snapshots u m contained in the database are equal to zero on the boundary, then by linear combinations Φ i = 0 on ∂Ω and the pressure term P i vanishes in (2.4).In most of the applications [3,10,14], the contribution of P i is simply neglected as a first approximation.In the traditional POD Galerkin approach, the coefficients and C GP ijk are then determined and the POD ROM (2.4) integrated in time from (2.5).A set of predicted time histories for a R i (t) is obtained and compared to a P i (t).In general, the original dynamics is not reproduced perfectly well.This can be explained [8] by the structural instability of the Galerkin projection, the truncation effect introduced in the Galerkin projection and to a lower level to the neglect of the pressure term.For this reason, it is then necessary to identify, or calibrate, whole or part of the coefficients.
To simplify the presentation, the POD Galerkin system (2.4) is written as: The POD Galerkin system (2.6) can be further simplified as follows: where 3 Calibration method

General formulation
The aim of the POD-based reduced-order model is to reproduce accurately the original dynamics given by the POD temporal eigenfunctions.It can be shown [8] that this is equivalent of solving the minimization problem given by: where the error e(y, t) is defined as follows: + e(0, t) .
In this paper, • To corresponds to the arithmetic time-average on N t equally spaced elements of the interval [0, T o ]: For the norm, we introduce Λ ∈ R N Gal ×N Gal the symmetric definite positive matrix associated to • Λ and define for any e ∈ R N Gal : Λ acts as a weight function by allowing to balance the importance of specific POD modes.When Λ = I N Gal , all the POD modes have the same importance in terms of error.
If Λ is a symmetric matrix, then the minimization problem (3.1) gives rise [8] to the linear system Ay = b, where In practice, the matrix A is often ill-conditioned leading to solutions of (3.2) that are very sensitive to perturbations [8].To overcome this difficulty, an idea is to determine the coefficients of calibration y as solutions of a regularized optimization problem aiming at minimizing a weighted average of the normalized error, and a term measuring the variation of the coefficients from their value obtained by POD Galerkin.Consequently, we introduce a new optimization problem defined as: where α ∈ [0, 1] is a weighting parameter.For α = 0, the calibrated model is fully optimized whereas for α = 1, the coefficients from the Galerkin projection are recovered.In (3.3), Π ∈ R Ny×Ny is a semi-norm on the polynomial vector space.For any y ∈ R Ny , we have where Π ∈ R Ny×Ny is a non-negative symmetric matrix.For Π = I Ny , all the coefficients are calibrated and have the same weight in the calibration.Finally, it can be shown [8] that the minimization of J α amounts to solve the linear system: where The parameter α which sets the level of regularization must be chosen with care.In the next section, the optimal value of α is determined by adapting the L-curve method proposed by [17] to our case.

Filter factors
To emphasize the role of the regularization parameter α in the calibration procedure, we introduce in this section the concept of filter factors through the use of the Singular Value Decomposition [15].Hereafter, we consider to simplify the particular case where1 Λ = I N Gal and Π = I Ny .Equivalently, the solution y α of (3.4) is also solution of The coefficient α is equal to αc with c 2 = J (y GP ) y GP 2 .Furthermore, the matrix A can be written as: where We now apply the SVD to the matrix B. By definition, it comes: where Ny are orthogonal matrices containing the left u j and right v j singular vectors.Here Σ = diag σ 1 , • • • , σ Ny is a diagonal matrix with the singular values σ j arranged in non-increasing order (σ By inserting the SVD of B into the expression of A α , we obtain: Finally, the regularized solution is given by: where are the filter factors.For α = 0 (no regularization), the solution (3.5) becomes This expression clearly illustrates what happens numerically if the linear system is not regularized.Indeed, if the Fourier coefficients v T j b , corresponding to the smallest singular values σ j do not decrease sufficiently fast compared to the singular values, then the solution y α is dominated by the terms in the sum corresponding to the smallest σ j .This behavior can be assessed by inspecting the discrete Picard condition plotted in Figure 1 and discussed in [16].The singular values decay faster than the Fourier coefficients and, as a result, the solution presents many oscillations around zero, and thus appears to be completely random.To regularize the problem, we then filter out the contributions to the solution corresponding to the smallest singular values by using the filter factor ϕ j .Obviously, as σ j decreases, ϕ j tends to zero, meaning that the contributions |v T j b| σ 2 j to the solution y α from the smallest σ j are damped.In addition, as σ j decreases, ψ j tends to one, so that each solution component filtered out by ϕ j is replaced by (v T j y GP )v j .As a consequence, if the majority of the singular values are approximately equal to zero (ϕ j ≈ 0, ψ j ≈ 1), then and we recover the coefficients determined by Galerkin projection.

Optimal value of α
A critical aspect of any regularization method is how to determine optimally the weighting parameter α or in other words how to find a fair balance between the minimization of the residual norm of (3.4) and the norm of y α .By analogy of the so-called L-curve method introduced by [17], a way is to display in log-log scale the variation of I(y) with respect to E(y).As long as the uncorrelated noise in b α dominates the more highly correlated noise in A α , this plot exhibits a typical "L" shape, and the optimal value of the regularization parameter α is considered intuitively to be at the corner of the L wherein both the residual and the norm simultaneously attain low values (see Figure 2).In [5], different iterative methods were used for the determination of a suitable value of the regularization parameter.Hereafter, we determine directly the curvature of the L-curve.The corner may be defined as a point where the curve (log(I(y)), log(E(y))) has its maximum curvature.Let us introduce some notations: The signed curvature of the L-curve is given by: where ρ , η , ρ , η are the first and second derivatives of ρ and η with respect to α.The goal is now to derive the expressions of these different derivatives for estimating κ for our linear problem.The first and second logarithmic derivatives of η and ρ are given by: The next step involves the computation of the first and second derivatives of ρ and η with respect to α.For ρ, we obtain: whereas for η, we find: Finally, the computation of the first and second derivative of y α with respect to α is straightforward since in (3.5),only the filter factors ϕ j and ψ j depend on α.

Numerical results
The calibration approach presented in Section 3 is now applied to a twodimensional incompressible cylinder wake flow at Re = 200.The database was computed using a finite-element code (DNS code Icare, see [13] for details) and consists of N s = 200 snapshots of the flow velocity taken evenly over a time horizon T s = 6 i.e. over more than one period of vortex shedding (T vs = 5).Typical iso-values of the longitudinal velocity are shown in Figure 3.The method of snapshots [23] is then applied to the velocity fluctuation.The relative kinetic energy is plotted in logarithmic scale for the first 40 POD modes in Figure 4. Arbitrarily, 99.9% of the flow energy contained in the cylinder wake is judged sufficient for describing correctly the original dynamics with a reduced-order model.This figure indicates clearly that the first six POD modes meet this criterion.For this reason, we consider N Gal = 6 to derive the POD ROM.After determination of the Galerkin coefficients of (2.4), the POD ROM is integrated in time with a classical fourth-order Runge-Kutta scheme and a time step of 10 −3 T s .A set of predicted time histories a R i (t) is obtained for i = 1, • • • , N Gal , and compared to the POD temporal eigenfunctions a P i (t).As shown in Figure 5, the original dynamics is globally well reproduced but the accuracy is not perfect.Depending on the targeted application of the reduced-order model, these representation errors may have no consequence or may lead to the failure of the application as it is the case for a model-based control approach.
To assess the numerical efficiency of the calibration procedure presented in Section 3, we now use the coefficients obtained as solution of the minimization problem (3.3)where α is determined as in Section 3.2.2. with the determination of all the coefficients (constant, linear and quadratic).The regularization parameter α is determined at the inflection point of the L-curve.
In Figure 6, the temporal evolutions of the POD modes are compared with those predicted by the calibrated reduced-order model.Contrary to the results presented in Figure 5 there is no visible difference in the dynamics.The immediate consequence is that the modal energy distribution associated to the calibrated model now corresponds perfectly to the POD energy (see Figure 7).
The difference between different methods of calibration can be analyzed precisely by introducing the modal errors J i defined for Λ = I N Gal as: In Figure 8, we compare the modal errors obtained with the proposed calibration method to those determined with the calibration method based on Tikhonov regularization developed in [8].We obtain that for the two first POD modes (the most energetic) the Tikhonov-based regularization method performs better than the proposed method, whereas for the higher modes the opposite behavior is found.
These slight reductions of the modal errors for the higher POD modes are an important feature of our method.Indeed, in the practical configurations (3-D turbulent flow obtained by numerical simulations or challenging experimental data) many more POD modes are active and then necessary to derive a realistic POD ROM.It is then critical to have a calibration method that does not only perform well for the higher energetic modes but also for the lower ones.

Conclusions
We have introduced a new regularized calibration method by defining a cost functional aiming at minimizing a weighted combination of normalized error, and a term measuring the variation of the coefficients from their value obtained classically by Galerkin projection.The efficiency of our approach was demonstrated for a two-dimensional cylinder wake flow by comparing it with the Tikhonov-based calibration method introduced in [8].An important feature of our approach is that it performs better for the higher POD modes (less energetic) that are known to be calibrated with more difficulties.The success of our approach is strongly related to the optimal determination of the regularization parameter inspired by the L-curve method.
As perspectives, we propose to reformulate the proposed method as a multicriteria optimization problem and to use a Nash equilibrium concept to solve it.The main idea is to search for the solution as the equilibrium point of a simulated dynamic game in which the set of parameters is split into subsets, each subset being considered as the strategy or territory of a given functional.This framework was already used in [11,20,22] to optimize a business-jet wingshape.We aim also to extend the concept of parametrization adaptivity, introduced and demonstrated in [12,19,21,24,25], in this context of reduced-order models.

2 Figure 1 .
Figure 1.Discrete Picard condition corresponding to the minimization of J (y) with determination of all the coefficients (constant, linear and quadratic) for the cylinder wake configuration presented in Section 4.

Figure 2 .
Figure 2. L-curve corresponding to the minimization of Jα(y) with determination of all the coefficients (constant, linear and quadratic) for the cylinder wake configuration presented in Section 4.

Figure 3 .
Figure 3. Cylinder wake flow at Re = 200.Iso-values of the longitudinal velocity fluctuation where the Bénard von-Kármán vortex shedding is clearly visible.

Figure 4 .
Figure 4. POD eigenvalues in logarithmic scale.E Ns corresponds to twice the energy contained in the database (E Ns = Ns j=1 λ j ).

Figure 5 .
Figure 5.Comparison between the temporal evolutions of the projected (POD) and predicted (POD ROM) mode amplitudes when the Galerkin coefficients are used for integrating the POD ROM.

Figure 6 .
Figure 6.Comparison between the temporal evolutions of the projected (POD) and predicted (POD ROM) mode amplitudes.The POD ROM is calibrated by minimizing Jα, with the determination of all the coefficients (constant, linear and quadratic).The regularization parameter α is determined at the inflection point of the L-curve.

Figure 7 .
Figure 7.Comparison between the modal energetic contents obtained before and after calibration.The POD eigenvalues are reported for reference.The POD ROM is calibrated by minimizing Jα, with the determination of all the coefficients (constant, linear and quadratic).The regularization parameter α is determined at the inflection point of the L-curve.

Figure 8 .
Figure 8. Modal errors J i .Comparison between the results obtained by the proposed approach and those determined with the calibration technique based on Tikhonov regularization developed in[8].