Multiobjective Optimal Control Methods for Fluid Flow Using Reduced Order Modeling

# Multiobjective Optimal Control Methods for Fluid Flow Using Reduced Order Modeling

Sebastian Peitz Department of Mathematics, University of Paderborn, Warburger Str. 100, D-33098 Paderborn Sina Ober-Blöbaum Department of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, UK Michael Dellnitz Department of Mathematics, University of Paderborn, Warburger Str. 100, D-33098 Paderborn
###### Abstract

In a wide range of applications it is desirable to optimally control a dynamical system with respect to concurrent, potentially competing goals. This gives rise to a multiobjective optimal control problem where, instead of computing a single optimal solution, the set of optimal compromises, the so-called Pareto set, has to be approximated. When the problem under consideration is described by a partial differential equation (PDE), as is the case for fluid flow, the computational cost rapidly increases and makes its direct treatment infeasible. Reduced order modeling is a very popular method to reduce the computational cost, in particular in a multi query context such as uncertainty quantification, parameter estimation or optimization. In this article, we show how to combine reduced order modeling and multiobjective optimal control techniques in order to efficiently solve multiobjective optimal control problems constrained by PDEs. We consider a global, derivative free optimization method as well as a local, gradient based approach for which the optimality system is derived in two different ways. The methods are compared with regard to the solution quality as well as the computational effort and they are illustrated using the example of the two-dimensional incompressible flow around a cylinder.

## 1 Introduction

In many applications from industry and economy, one is interested in simultaneously optimizing several criteria. For example, in transportation one wants to reach a destination as fast as possible while minimizing the energy consumption. This example illustrates that in general, the different objectives contradict each other. Therefore, the task of computing the set of optimal compromises between the conflicting objectives, the so-called Pareto set, arises. This leads to a multiobjective optimization problem (MOP) or, if the control variable is a function, a multiobjective optimal control problem (MOCP). Based on the knowledge of the Pareto set, a decision maker can use this information either for improved system design or for changing control parameters during operation as a reaction on external influences or changes in the system state itself.

Multiobjective optimization is an active area of research. Different approaches exist to address MOPs, e.g. deterministic approaches [Mie12, Ehr05], where ideas from scalar optimization theory are extended to the multiobjective situation. In many cases, the resulting solution method involves solving multiple scalar optimization problems consecutively. Continuation methods make use of the fact that under certain smoothness assumptions the Pareto set is a manifold that can be approximated by continuation methods known from dynamical systems theory [Hil01]. Another prominent approach is based on evolutionary algorithms [ALG05, CCLVV07], where the underlying idea is to evolve an entire set of solutions (population) during the optimization process. Set oriented methods provide an alternative deterministic approach to the solution of MOPs. Utilizing subdivision techniques, the desired Pareto set is approximated by a nested sequence of increasingly refined box coverings [DSH05, SWOBD13].

When dealing with control functions, a multiobjective optimal control problem (MOCP) needs to be solved. Similar to MOPs, ideas from scalar optimal control theory (see e.g. [HPUU09] for an overview of optimal control methods for PDE-constrained problems) can be extended to take into account multiple objectives. By applying a direct method, the MOCP is transformed into a high-dimensional, nonlinear MOP such that the methods mentioned before can be applied. Another approach is based on the transformation of the MOCP into a sequence of scalar optimal control problems and the use of well established optimal control techniques for their solution. Examples for MOCPs can be found e.g. in [DEF16, LHDvI10, OBRzF12, SWOBD13] with ODE constraints. In [LKBM05] as well as [ARFL09], nonlinear PDE constraints are taken into account but the model is treated as a black box, i.e. no special treatment of the constraints is required. The first articles explicitly taking into account PDE constraints are [IUV13, ITV15], where multiobjective optimal control problems are solved with a weighted sum approach and model order reduction techniques subject to linear and semilinear PDE constraints, respectively. Fluid flow applications have been considered in [OBPG15, PD15].

All approaches to MOP / MOCP have in common that a large number of function evaluations is typically required. Thus, the direct computation of the Pareto set can quickly become numerically infeasible. This is frequently the case for problems described by (nonlinear) partial differential equations such as the Navier-Stokes equations. Standard optimization methods for PDEs [HPUU09] often make use of a discretization by finite elements, finite volumes or finite differences which results in a high-dimensional system of ODEs. In a multi query context (such as optimization, parameter identification, etc.), this approach often exceeds the limits of today’s computing power. Hence, one aims for methods which reduce the computational cost significantly. This can be achieved by approximating the PDE by a reduced order model of low dimension.

In recent years, major achievements were made in the field of reduced order modeling [SvdVR08, BMS05]. In fact, different methods for creating low dimensional models exist for linear systems (e.g. [WP02]) as well as for nonlinear systems (e.g. [BMNP04, GMNP07]), see [ASG01] for a survey. Many researchers focus their attention to the Navier-Stokes equations (e.g. [CBS05, Row05, XFB14]), where Proper Orthogonal Decomposition (POD) [HLB98] has proven to be a powerful tool. In the context of fluids, this technique has first been introduced by Lumley [Lum67] to identify coherent structures. Reduced order models using POD modes and the method of snapshots go back to Sirovich [Sir87].

Due to the wide spectrum of potential applications (optimal mixing, drag reduction, HVAC (Heating, Ventilation and Air Conditioning), etc.) and the progress in computational capabilities, a lot of research is devoted to the control of the Navier-Stokes equations, either directly [GeH89, Lac14] or via reduced order modeling [GPT99, Fah00, IR01, BCB05, BC08, PD15]. In many cases, the energy consumption plays an important role. Thus, ideally, one wants to consider the control cost in addition to the main objective in order to choose a good trade-off between the achievement of the desired objective and the respective control cost. Since this causes a considerable computational effort in case of systems described by the Navier-Stokes equations, we show in this article how reduced order modeling and multiobjective optimal control methods can be combined to solve MOCPs with nonlinear PDE constraints. Using model order reduction via POD and Galerkin projection, we present several methods to compute the Pareto set for the conflicting objectives minimization of flow field fluctuations and control cost for the two-dimensional flow around a cylinder at . We discuss the advantages and disadvantages of the different approaches and comment on the respective computational cost. In fact, it is shown that the different methods strongly differ in their respective efficiency and that gradient based approaches have a better performance, provided that the gradient is sufficiently accurate. We also show that the choice of the algorithm and the corresponding particular numerical realization of the control function lead to different optimal control behavior.

The article is organized as follows. In Section 2 we present the problem setting and formulate the MOCP. After an introduction to multiobjective optimal control in Section 3, the reduced order model and the resulting reduced MOCP are derived in Section 4. We then present our results in Section 5 and draw a conclusion in Section 6.

## 2 Problem formulation Figure 1: Sketch of the domain Ω⊂R2. The length is L=25d, the height H=15d, and the cylinder center is placed at (5d,7.5d).

The two-dimensional, viscous flow around a cylinder is one of the most extensively studied cases for separated flows in general as well as for flow control problems [GPT99, GBZI04, BCB05]. In this paper, we consider the laminar case described by the incompressible Navier-Stokes equations at a Reynolds Number computed with respect to the far field velocity (throughout the paper, we will use bold symbols for vector valued quantities), the kinematic viscosity and the cylinder diameter :

 =−∇p(x,t)ρ+1Re∇2U(x,t), (2.1a) ∇⋅U(x,t) =0, (2.1b) (U(x,0),p(x,0)) =(U0(x),p0(x)), (2.1c) for x∈Ω, t∈[t0,te],

where is the two-dimensional fluid velocity and the pressure. is the standard Sobolev space (cf. e.g. [HPUU09]). The domain (cf. Figure 1) is denoted by . We impose Dirichlet boundary conditions at the inflow as well as the upper and lower walls . At the outflow , we impose a standard no shear stress condition [HRT96]:

 U(x,t) =(U∞,0) for x∈ΓD, (2.1d) p(x,t) n =1Re∂U(x,t)∂n for x∈ΓN, (2.1e)

where is the outward normal vector of the boundary. On the cylinder , we prescribe a time-dependent Dirichlet BC such that it performs a rotation around its center with the angular velocity :

 (2.1f)

with according to Figure 1. The cylinder rotation serves as the control mechanism for the flow. The Hilbert space is equipped with the inner product and the norm .

Following [Fah00, BCB05], we introduce the weak formulation of (2.1a). Consider the divergence free Hilbert space of test functions . Then, a function which satisfies

 (∂U∂t+(U⋅∇)U,ψ)=(pρ,∇⋅ψ)−[pρψ]−1Re(∇ψ,(∇U)⊤)+1Re[(∇U)⊤ψ] (2.2)

for all is called a weak solution of (2.1a). Here, is the boundary integral (e.g. ) and is the inner product for vector-valued quantities (e.g. ). Note that (2.1b) is automatically satisfied by design of the test function space .

### 2.1 Finite volume discretization

The system (2.1a2.1f) is solved by the software package OpenFOAM [JJT07] using a finite volume discretization and the PISO scheme [FP02]. We have chosen OpenFOAM because it contains a variety of efficient solvers for various fluid flow applications. Since we will utilize finite elements for the computation of the reduced order model (evaluation of inner products etc.), we then map our solution to a finite element mesh with degrees of freedom (Figure 2(a)). This is done in the spirit of data driven modeling, where we collect data which does not necessarily have to come from a numerical method. Finally, the velocity field can be written in terms of the FEM basis:

 U(x,t)=⎛⎝∑Nj=1Udj(t)ϕj(x)∑Nj=1Udj+N(t)ϕj(x)⎞⎠, (2.3)

where are the FEM basis functions and are the nodal values of the two velocity components, the superscript denoting that this is a quantity defined on the grid nodes. In the following, the nodal values of all quantities will be denoted by a superscript . All finite element related computations are performed with the FEniCS toolbox [LMW12] using linear basis functions. Figure 2: (a) FEM discretization of the domain Ω by a triangular mesh (N=17.048). (b) A snapshot of the solution to (2.1a) – (2.1f) for a non-rotating cylinder (γ(t)=0), the coloring is according to the velocity magnitude. The pattern is the well-known von Kármán vortex street.

For a non-rotating cylinder, i.e. , the system possesses a periodic solution, the well-known von Kármán vortex street (Figure 2(b)), where vortices separate alternatingly from the upper and lower surface of the cylinder, respectively. The effect is observed frequently in nature and is one of the most studied phenomena in fluid mechanics, also in the context of flow control, where the objective is to stabilize the flow and to reduce the drag.

### 2.2 Multiobjective optimal control problem

In many applications, the control cost is of great interest. This is immediately clear when the goal of the optimization is to save energy such that in this case, the control effort needs to be taken into account. In scalar optimization problems, this is often done by adding an additional term of the form to the cost functional where is a weighting parameter. Here, we want to consider the two objectives flow stabilization, i.e. the minimization of the fluctuations around the mean flow field , and the minimization of the control cost separately which leads to the following multiobjective optimal control problem:

 minU,γˆJ(U,γ)=minU,γ(∫tet0∥u(⋅,t)∥2L2dt∥γ∥2L2),

where and satisfies (2.1a2.1f). In contrast to bounded domains (cf. [FGH98]), the proof of existence of a solution is an open problem for cases with no-shear or do nothing boundary conditions [Ran00]. Nevertheless, based on numerical experiences [Ran00], we will from now on assume that there exists a unique solution for each and hence, we denote by the solution for a fixed and consider the reduced cost functional which leads to the following multiobjective optimal control problem:

 minγJ(γ)=minγ(∫tet0∥u∥2L2dt∥γ∥2L2). (MOCP)

In general, the solution to this problem does not consist of isolated points but a set of optimal compromises between the two objectives. In the following section, we give a short introduction to multiobjective optimal control theory and solution methods.

## 3 Multiobjective Optimal Control

This section is concerned with the treatment of general multiobjective optimal control problems. We will give a short introduction to the general theory before addressing the two algorithms used later on in combination with model order reduction techniques.

### 3.1 Theory of multiobjective optimal control

Consider the general multiobjective optimal control problem:

 minγJ(γ)=minγ⎛⎜ ⎜⎝J1(γ)⋮Jk(γ)⎞⎟ ⎟⎠, (2)

with and , . The space in which the control functions live is denoted as the decision space and the function is a mapping to the -dimensional objective space. The set of feasible functions is the feasible set in decision space. We denote its image as the feasible set in objective space which consists of the feasible points . In contrast to single objective optimization problems, there exists no total order of the objective function values in . Therefore, the comparison of values is defined in the following way [Mie12]:

###### Definition 3.1.

Let . The vector is less than , if for all . The relation is defined in an analogous way.

A consequence of the lack of a total order is that we cannot expect to find isolated optimal points. Instead, the solution to (2) is the set of optimal compromises, the so-called Pareto set named after Vilfredo Pareto:

###### Definition 3.2.
1. A function dominates a function , if and .

2. A feasible function is called (globally) Pareto optimal if there exists no feasible function dominating . The image of a (globally) Pareto optimal function is a (globally) Pareto optimal point.

3. The set of nondominated feasible functions is called the Pareto set, its image the Pareto front.

Consequently, for each function that is contained in the Pareto set (cf. the red line in  Figure 3(a)), one can only improve one objective by accepting a trade-off in at least one other objective. Figuratively speaking, in a two-dimensional problem, we are interested in finding the ”lower left” boundary of the feasible set in objective space (cf. Figure 3(b)). More detailed introductions to multiobjective optimization can be found in [Mie12, Ehr05].

### 3.2 Solution methods

Various methods exist to solve problem (2). In this work, we present two approaches that are fundamentally different. The first one is a reference point method [RBW09] for which the distance between a feasible point (i.e. an objective value that lies in the feasible set in objective space) and an infeasible target in objective space is minimized. The method yields a moderate number of single objective optimization problems that are solved consecutively. The second approach is a global, derivative free subdivision algorithm [DSH05] for which the Pareto set is approximated by a nested sequence of increasingly refined box coverings. However, its applicability depends critically on both the decision space dimension and the numerical effort of function evaluations. The main properties are summarized in Table 3.2.

#### 3.2.1 Reference point method

The reference point method presented here belongs to the category of scalarization techniques for which the solution set of (2) is approximated by a finite set of points, each computed by solving a scalar optimization problem. In the beginning, one Pareto optimal point has to be known. This can be achieved by solving a scalar optimization problem for some weighted sum of all objectives (i.e.  for the case of two objectives), including the scalar optimization with respect to any of the objectives of (MOCP). Then, a so-called target is chosen such that it lies outside the feasible set in objective space, e.g. by shifting the solution of the first Pareto point (). We then solve the scalar optimization problem . As a result, the corresponding optimal point lies on the boundary of the feasible set and it is not possible to further improve all objectives at the same time. Thus, the point is (locally) Pareto optimal. By adjusting the target position based on targets and Pareto points already known, multiple points on the Pareto front (i.e. ) are computed recursively (cf. Figure 4 for an illustration). For being continuous, the change in decision space is small when the target position changes only slightly [Hil01] and hence, the current solution is a good initial guess for the next scalar problem which accelerates the convergence considerably. Figure 4: Reference point method in image space. (a) Determination of the i-th point on the Pareto front by solving a scalar optimization problem. (b) Computation of new target point Ti+1 and predictor step in decision space (γp,i+1, cf. line 16 in Algorithm 1).

Theoretically, the method is not restricted to low objective space dimensions. However, setting the targets properly to get a proper approximation of the Pareto front (i.e. an approximation of the entire front by evenly distributed points) becomes complicated in higher dimensions. In our case, we are dealing with two objectives and the targets can be determined easily using linear extrapolation (cf. Figure 4 and Algorithm 1) as proposed in [RBW09]. This also allows us to compute the whole front in at most two search directions (, line 2 in Algorithm 1). From the initial point, we first proceed in one direction, e.g. decreasing . When, at some point, is increasing again (line 6 in Algorithm 1), we have reached the extremal point of the feasible set (cf. Figure 3(b)) and the scalar optimum of . We then return to the initial point and proceed in the opposite direction (lines 8, 9 in Algorithm 1) until the other extremal point is reached.

###### Remark 3.3.

Using system knowledge, the algorithm can be further simplified. In the case of (MOCP), for example, we know that is Pareto optimal since it is the scalar minimum of and therefore has the lowest possible value of . Knowing this, we do not need to solve the initial optimization problem and moreover, we only need to move along the front in one direction until the other extremal point is reached.

The scalar optimization problems can be solved using any suitable method. Here, we use a line search approach [NW06] and compute the derivatives using an adjoint approach. Note that since the scalar optimization routine often is of local nature, the method overall is also local and depends on the initial guesses as well as the choice of the target points. However, under certain smoothness assumptions, the Pareto set is a -dimensional manifold [Hil01] such that once the first point is computed correctly, the method is promising to find the globally optimal Pareto set for many problems.

###### Remark 3.4.

For the reference point method, it turned out to be numerically beneficial that . Otherwise, the computation of new target points may become sensitive to the step length parameters.

#### 3.2.2 Global subdivision algorithm

When the image space dimension increases or, alternatively, the Pareto front is disconnected (see e.g. [DEF16]), continuation methods like the reference point method may fail. Moreover, the resulting scalar optimization problems are often solved by algorithms of local nature which may also not be sufficient if global optima are desired. In addition to that, derivatives are hard to compute or not available at all in many applications. The subdivision algorithm presented here overcomes these problems. It is described in more detail in [DSH05] including a proof of convergence. The version using gradient information is based on concepts for the computation of attractors of dynamical systems [DH97], it is however not considered here. Instead, we focus on the derivative free approach. It is very robust and since we directly utilize the concept of dominance when comparing solutions, we theoretically do not need to make any assumptions about the problem except that the control is real-valued, i.e. . In the numerical realization however, we approximate sets by sample points and hence, we further require that the objective is a continuous function of .

The subdivision algorithm (cf. Algorithm 2 and Figure 5) is a set oriented method where the Pareto set is approximated by a nested sequence of increasingly refined box coverings (cf. Figure 5). During the algorithm, a sequence of box collections covering the Pareto set is constructed, starting with a sufficiently large initial set . (This results in box constraints for the control .) In each iteration, the algorithm performs the steps subdivision and selection until a given precision criterion is satisfied. This way, a subset of the previous box collection remains that is a closer covering of the desired set. In the selection step, the boxes are compared pairwise and all boxes that are dominated by another box are eliminated from the collection. Numerically, this is realized by a representative set of sampling points in each box (cf. Figure 5(a)), at which the objectives are evaluated. In this case, we say that a box is dominated if all sampling points are dominated by at least one point from another box. The remaining boxes are then subdivided into two boxes of half the size and we proceed with the next selection step. The subdivision step is performed cyclically in the decision space dimensions to . Hence, in order to divide each dimension times, we require subdivision steps in total. Figure 5: Global subdivision algorithm, example with γ∈R2. (a) Sampling. (b) Nondominance test. (c) Elimination of dominated boxes. (d) Subdivision, sampling and nondominance test.

Since the applicability of the subdivision algorithm depends critically on both the decision space dimension and the numerical effort of function evaluations, it is in practice limited to low decision space dimensions. Therefore, we introduce a sinusoidal control, i.e. . The choice is motivated by the fact that the uncontrolled dynamics of the problem at hand are also periodic and we include a phase shift to allow for controls with non-zero initial conditions. This way, we transform (MOCP) into an MOP with :

 minA,ω,τ∈R˜J(A,ω,τ)=minA,ω,τ∈R⎛⎝∫tet0∥u∥2L2 dt∫tet0(Asin(2πωt+τ))2 dt⎞⎠. (MOP-3D)

## 4 Reduced Order Modeling

The problem (MOCP) could now be solved using either of the methods presented in Section 3.2. However, both require a large number of evaluations of the cost functional and consequently, evaluations of the system (2.1a2.1f). A finite element or finite volume discretization yields a large number of degrees of freedom such that solving (MOCP) quickly becomes computationally infeasible. Hence, we need to reduce the cost for solving the dynamical system significantly. This is achieved by means of reduced order modeling where we replace the state equation by a reduced order model of coupled, nonlinear ordinary differential equations that can be solved much faster.

### 4.1 Reduced Order Model

The standard procedure for deriving a reduced order model is by introducing a Galerkin ansatz [HLB98]:

 U(x,t)=l∑j=1αj(t)ψj(x),(x,t)∈Ω×[t0,te], (3)

where are time-dependent coefficients and are basis functions. In contrast to FEM, these basis functions may in general have full support, i.e. they are global. The objective is now to find a basis as small as possible () while allowing for a good approximation of the flow field.

Using (3), we can derive a reduced model of dimension . By selecting a basis that is divergence-free, the mass conservation equation (2.1b) is automatically satisfied. (When computing this basis using data from a divergence-free flow field, the resulting basis elements are also divergence-free [BHL93].) Inserting (3) into the weak formulation of the Navier-Stokes equations (2.2) yields a system of ordinary differential equations of the form

 ˙α(t)=˜Bα(t)+˜C(α(t)),

where the coefficient matrices are computed by evaluating the inner products:

 (˜B)ij =1Re(∇ψi,∇ψj)L2, (˜C(α(t)))j =α(t)⊤Qjα(t), with (Qj)ik =((ψi⋅∇)ψk,ψj)L2.
###### Remark 4.1.

In the above as well as all following formulations, we have neglected the influence of the pressure term. In [NPM05], it was shown that including the pressure term is favorable for the accuracy of the reduced model for open systems. Since we consequently perform a model stabilization in order to increase the accordance with the flow field data, this results in a further modification of the model coefficients. Consequently, we have (following [CEF09]) decided to neglect the influence of the pressure term in this work.

This model, however, is not controllable in the sense that we can no longer influence the solution by choosing . Moreover, the basis functions need to be computed from data with homogeneous boundary conditions ( for ) such that they also possess homogeneous boundary conditions [Sir87]. Otherwise, it cannot be guaranteed that the Galerkin ansatz (3) satisfies the PDE boundary conditions. The idea is therefore to introduce the Galerkin expansion for a modified flow field which satisfies homogeneous boundary conditions for all . In order to realize this, we introduce the following flow field decomposition [Rav00, BCB05]:

 U(x,t) =⟨U(x,t)⟩+γ(t)Uc(x)+u(x,t), (4)

where is a so-called control function which is a solution to (2.1a), (2.1b), (2.1c) with a constant cylinder rotation and zero boundary conditions elsewhere. The choice of the rotating velocity for the computation of is arbitrary but has an influence on the resulting reduced model such that the model performance might be influenced by this choice. Since the main intention of this article is to develop multiobjective optimal control algorithms for PDE-constrained problems, we do not address this problem further and choose such that the cylinder surface rotates with a velocity of . The decomposition is now computed in two steps:

 ˜U(x,t) =U(x,t)∣∣γ(t)=γref(t)−γref(t)γcUc(x), u(x,t) =˜U(x,t)−⟨˜U(x,t)⟩,

where is the reference control for which the data was collected.

In this way, satisfies homogenous boundary conditions for all . (Note that this does not hold exactly on , however, the resulting error is small and hence neglected [BCB05].) Inserting (4) into (2.2) yields an extended reduced order model:

 ˙α(t) =A+Bα(t)+C(α(t))+D˙γ(t)+(E+Fα(t))γ(t)+Gγ2(t), (5) αj(t0) =(u(⋅,t0),ψj)L2=:α0,j, (6)

where . The expressions for the coefficient matrices to are given in Appendix A, their size depends on the dimension of the ROM: with .

###### Remark 4.2.

Since the time derivative of the control occurs in (5), we require higher regularity, i.e. .

Due to the truncation of the POD basis, the higher modes covering the small scale dynamics, i.e. the dynamics responsible for dissipation, are neglected. This can lead to incorrect long term behavior and hence, the model needs to be modified to account for this. Several researchers have addressed this problem and proposed different solutions, see e.g. [Rem00, SK04, NPM05, CEF09, BBI09]. Here, we achieve stabilization of the model (5) by solving an augmented optimization problem [GBZI04]. In this approach the entries of the matrix are manipulated in order to minimize the difference between the ROM solution and the projection of the data onto the reduced basis (). The resulting initial value problem is solved using a Runge-Kutta method of fourth order.

Finally, we replace the objectives in (MOCP) by their equivalents in the reduced order model. Since the reference point method is more stable for objectives of the same order of magnitude (cf. Remark 3.4), we additionally multiply the second component by (i.e. the number of basis functions)

 minγ⎛⎝∫tet0∫Ω∥∑li=1αi(t)ψi(x)∥22dxdtl∫tet0γ2(t)dt⎞⎠. (7)

### 4.2 Proper orthogonal decomposition (POD)

In order to achieve high numerical efficiency, we want to compute a basis as small as possible. To this end, we utilize the well-known Proper Orthogonal Decomposition (POD). Using POD, we can compute the orthonormal basis which, for every , optimally represents a set of flow field realizations [Sir87]. Basically, the computation is realized via a singular value decomposition of the so-called snapshot matrix , where we collect measurements (i.e. the velocity components at nodes) at different time instants. If , as it is often the case for FEM discretizations, it is more efficient to solve the -dimensional eigenvalue problem [KV02, Fah00]:

 S⊤MSvi=σivi,i=1,…,m, (8)

where is the finite element mass matrix. Using the eigenvalues and eigenvectors from (8), the POD modes can be computed by making use of the FEM basis functions :

 ψdi =1√σiSvi, ψi(x) =⎛⎝∑Nj=1ψdi,jϕj(x)∑Nj=1ψdi,j+Nϕj(x)⎞⎠.

The eigenvalues are a measure for the information contained in the respective modes and hence, the error between the Galerkin ansatz and the snapshot data is determined by the ratio of the truncated eigenvalues [Sir87], i.e. . By choosing a value for , typically or , the basis size can be determined. For many applications, the eigenvalues decay fast such that a truncation to a small basis is possible.

The error is only known for the reference control at which the data was collected. If one allows for multiple evaluations of the finite element solution, the reduced model can be updated when necessary. This is realized using, e.g., trust-region methods [Fah00, BC08]. An alternative approach is to use a reference control that yields sufficiently rich dynamics such that the model can be expected to be valid for a larger variety of controls [BCB05]. Since it is computationally cheaper, we follow the second approach and take snapshots in the interval with for a chirping reference control (cf. Figure 6). This leads to a basis of size for . Figure 7 depicts the comparison between the uncontrolled solution (i.e. ) of the reduced model and the projection of the finite volume solution for reduced models computed with a reference control and , respectively. Figure 6: Chirping function (γchirp=−4sin(2πt/120)cos(2πt/3−18sin(2πt/60))) used for the computation of the snapshot matrix [BCB05]. Figure 7: First six modes of (5), (6) (—) and projection of FEM solution (∘) for two reduced order models with different reference controls. (a) γref=0. (b) γref=γchirp.

Figure 8 shows the first four POD modes for the uncontrolled solution. The modes occur in pairs, the second one slightly shifted downstream. This is due to symmetries (see the double eigenvalues in Figure 9(a)) in the problem, namely in the horizontal axis through the cylinder as well as on the upper and lower boundary, respectively. The first four modes already account for of the information, cf. Figure 9, where in (a) the eigenvalues are shown for both the uncontrolled flow and the flow controlled by the chirping function. The error, i.e. the ratio of the truncated eigenvalues, is visualized in (b). Figure 8: The first four POD modes for an uncontrolled solution. Figure 9: Eigenvalues (a) and error at current basis size (b) for an uncontrolled solution and a solution controlled by the chirping function.

Making use of the orthonormality of the POD basis, the first objective in (7) can be simplified:

 ∫Ω∥∥l∑i=1αi(t)ψi(x)∥∥22dx =∫Ω(l∑i=1αi(t)ψi(x))⋅(l∑i=1αi(t)ψi(x))dx=l∑i=1α2i(t) with (ψi,ψj)L2=δi,j,

which allows us to replace (7) by a simpler formulation:

 minγJ(α,γ)=minγ⎛⎝∫tet0∑li=1α2i(t)dtl∫tet0γ2(t)dt⎞⎠. (R-MOCP)

We would like to use a gradient based algorithm for the scalar optimization problems occurring in the reference point method (Section 3.2.1). We realize this with an adjoint approach. The first step is to transform (R-MOCP) into a single objective optimization problem where the euclidean distance to a target is to be minimized:

 minγ¯¯¯¯J(α,γ) =minγ∥T−J(α,γ)∥22, (SOP) with¯¯¯¯J(α,γ) =(T1−J1(α,γ))2+(T2−J2(α,γ))2 =(T1−∫tet0l∑i=1α2i(t)dt)2+(T2−l∫tet0γ2(t)dt)2, (9)

where according to (R-MOCP) and satisfies the reduced order state equation (5), (6).

The task is now to solve a sequence of scalar optimization problems with varying targets . We choose a line search strategy [NW06] to address (SOP), where the direction is computed with a conjugate gradient method and the step length via backtracking and a valid length has to satisfy the Armijo rule. The necessary gradient information is computed with an adjoint approach. To this end, we introduce the adjoint state and the Lagrange functional:

 L(α,λ,γ)=∫tet0¯¯¯¯J(α,γ)−λ⊤(˙α−A−Bα−C(α)−D˙γ−(E+Fα)γ−Gγ2)dt, (10)

which is stationary for optimal values of and the corresponding state and adjoint state . Using integration by parts, this leads to the following system of equations111Note, that since this approach makes use of variational calculus, it is formally only applicable in the case of stronger assumptions, namely for (since we apply integration by parts).:

 ˙α =A+Bα+C(jα)+D˙γ+(E+Fα)γ+Gγ2, (11) ˙λ =−∂¯¯¯¯J∂α−(B+∂C(α)∂α+Fγ)⊤λ, (12) 0 =∂¯¯¯¯J∂γ+(E+Fα+2Gγ)⊤λ−D⊤˙λ=:Dγ¯¯¯¯J, (13)

with the Fréchet derivatives

 ∂¯¯¯¯J/∂αj =4(∫tet0l∑i=1α2i(t)dt−T1)αj(t), ∂¯¯¯¯J/∂γ =4l(l∫tet0γ2(t)dt−T2)γ(t),

and the respective boundary conditions and . A detailed derivation of the optimality system is given in Appendix B.

Equations (11) – (13) form the so-called optimality system which is seldom solved explicitly. Instead, the state equation (11) and the adjoint equation (12) are solved by forward and backward integration, respectively. The optimality condition (13) provides the derivative information of the cost functional with respect to the control . This information can be used to improve an initial guess for the optimal control within an iterative optimization scheme as described in Algorithm 3.

In [BCB05], the condition is the only boundary condition for the adjoint equation and the condition (cf. Appendix B) is neglected. In our computations, using this optimality system caused convergence issues. The reason for this is that the gradient obtained by the adjoint approach appears to be wrong, cf. Figure 10(a) for a comparison to a finite difference approximation. (Since the state equation is very accurate, the gradient can be approximated well this way.) The reason for this might be that we neglect the additional condition for the adjoint equation. However, it has been reported before [Fah00, DV01] that it is favorable to include information of the adjoint state of the PDE in the reduced order model. In fact, the adjoint based gradient can easily become inaccurate or even incorrect if based on information about the state only. Consequently, it is difficult to say whether the convergence problems stem from the missing boundary condition or are a flaw of the model reduction approach itself. To further investigate this, we derive a second optimality system based on an alternative formulation of the reduced model. Following [Rav00], we define an augmented state and introduce a new control . Figure 10: Comparison of the gradient computed using the adjoint approaches with an approximation by finite differences (FD). (a) Optimality system (11) – (13) (FD: Dγ,i¯¯¯¯J=(¯¯¯¯J(γ+hei)−¯¯¯¯J(γ))/h, where hei is the variation by h of the i-th entry of the discretization of γ). (b) Optimality system (14) – (18) (FD: Dv,i¯¯¯¯J=(¯¯¯¯J(v+hei)−¯¯¯¯J(v))/h).

It is known from optimal control theory that if a dynamical system depends linearly on the control, the control is bounded by box constraints and does not appear in the cost functional, the optimal control is of bang bang type and might include singular arcs [Ger12]. This is exactly the situation for the second formulation and since the computation of such solutions is numerically challenging, we avoid this behavior by adding a regularization term to the second objective in (R-MOCP), i.e. , where is a small number (here: ). The optimality system can be derived in an analogous way to the one described above, where the adjoint states are now denoted by and :

 ˙α =A+Bα+C(α)+Dv+(E+Fα)γ+Gγ2, (14) ˙γ =v, (15) ˙λ =−∂¯¯¯¯J∂α−(B+∂C(α)∂α+Fγ)⊤λ, (16) ˙μ =−∂¯¯¯¯J∂γ−(E+Fα+2Gγ)⊤λ, (17) 0 =∂¯¯¯¯J∂v+D⊤λ+μ=:Dv¯¯¯¯J, (18)

with the Fréchet derivatives:

 ∂¯¯¯¯J/∂γ =4l(l∫tet0γ2(t)dt+β∫tet0v2(t)dt−T2)γ(t), ∂¯¯¯¯J/∂v =4β(l∫tet0γ2(t)dt+β∫tet0v2(t)dt−T2)v(t),

as before and the respective boundary conditions

 α(0)=α0, λ(T)=0, μ(0)=0, μ(T)=0. (19)

This yields a boundary value problem for the state equations (14), (15) and the adjoint equations (16), (17) which can, for a given control , be solved by an iterative scheme. Applying a shooting method, the initial value is computed such that is satisfied using Newton’s method (cf. Algorithm 4). This is computationally more expensive than the simple forward-backward integration in Algorithm 3. However, the system (14) – (18) yields strongly improved convergence and decreased sensitivity to the optimization parameters in comparison to the system (11) – (13), which stems form the improved gradient accuracy (cf. Figure 10(b)). This is also evident from Figure 11, where the Pareto front based on (14) – (18) was computed executing Algorithm 4 once, starting in the already known point of zero control cost. In contrast to that, multiple attempts were made with Algorithm 3, starting from different points on the Pareto front that were computed with Algorithm 4. All of them either divert quickly from the front or stop prematurely. The results shown in Section 5 are therefore based on (14) – (18). From our experience, the approach based on the adjoint equation can be very sensitive, also with respect to stabilization methods where the coefficients are manipulated to better match the state equation. Hence, it would be favorable for future work to utilize direct approaches where only accuracy of the state equation is of importance. Figure 11: Comparison of the Pareto fronts computed with the two optimality systems (11) – (13) and (14) – (18).

## 5 Results

In this section we present solutions to (R-MOCP) obtained with the algorithms presented in Section 3. For all cases, we fix the time interval to which corresponds to roughly two vortex shedding cycles. In order to use the subdivision method, we transform (MOCP) into a multiobjective optimization problem with moderate decision space dimension (). One way to do this is by introducing a sinusoidal control (cf. Section 3.2.2) which then yields (MOP-3D). Another way is to approximate by a cubic spline with break points. In this case, the parameters to be optimized are the spline values at these break points. For the reference point method, we discretize with a constant step length of and solve Equations (14) – (17) with a fourth order Runge-Kutta scheme.

Figure 12 shows the box covering of the Pareto set for (MOP-3D), the respective Pareto front is shown in Figure 13. It is noteworthy that, apart from a few spurious points at caused by numerical errors (Figure 5(a)), the largest part of the set is restricted to a small part of the frequency domain. This is comprehensible since the rotation counteracts the natural dynamics of the von Kármán vortex street. The trade-off between control cost and stabilization is almost exclusively realized by changing the amplitude of the rotation. Figure 12: Box covering of the solution of (MOP-3D), obtained with global subdivision algorithm after 27 subdivision steps. There are spurious solutions with almost zero amplitude, the right plot shows only the physically relevant part which is constrained to a small part of the frequency domain. The trade-off between the two objectives is mainly realized by changing the amplitude.

When moving towards the two scalar optimal solutions, i.e. minimal fluctuations and minimal control cost, respectively, further improvements are only possible by accepting large trade-offs in the other objective which is a common phenomenon in multiobjective optimization. When designing a system, one could now accept a small increase in the main objective, i.e. flow stabilization, in order to achieve a large decrease in the control cost. This becomes more clear when looking at Figure 14(b), where different optimal compromises for the control are shown. For a relatively small improvement of from to (), the control cost increases by . Figure 13: (a) Comparison of Pareto fronts for different MOCP solution methods. (b) Comparison of the Pareto optimal controls with equal cost (J2=15) for different MOCP solution methods. Figure 14: Different Pareto points for the reference point method (a) and (MOP-3D) (b).

Figures 13(a) and 14 also show the Pareto fronts and Pareto optimal controls, respectively, computed with the other methods, i.e. two different spline approximations ( and ) computed with the subdivision algorithm on the one hand and arbitrary control functions determined by the reference point method on the other hand. We observe that a spline with 5 breakpoints is too restrictive to find control functions of acceptable quality. When using 10 breakpoints on the other hand, the Pareto front clearly surpasses the solution of (MOP-3D). It is, however, numerically expensive to compute. The best solution is computed with the reference point method. The improvement over the spline based solution is relatively low, the reason being that the controls are almost similar (see Figure 13(b), where solutions of the different methods with the same control cost are compared). When considering scenarios with larger time intervals than 10 seconds, the spline dimension would have to be increased further, leading to again much higher computational cost whereas the cost for the reference point method increases only linearly with time due to the fact that only the number of time steps in the forward and backward integration (14) – (17) is affected. Finally, Figure 15 shows two solutions of (2.1a) – (2.1f) with controls computed by the reference point method, corresponding to different values of the control cost and hence different degrees of stabilization. Due to the relatively short integration time of seconds, changes are only apparent closely past the cylinder. Still, we observe increased stabilization when allowing higher control cost. Figure 15: FV solution at t=10 with two different Pareto optimal solutions computed with the reference point method. (a) Relatively low control cost (J2=1.84), almost no reduction of fluctuations (J1=85). (b) Larger control cost (J2=30.12), stronger reduction of fluctuations (J1=75). Due to the short integration time, the effect is only visible in the vicinity of the cylinder.

In order to evaluate the quality of the solution obtained with the reduced order modeling approach, we have evaluated the cost function using the system state obtained from a PDE evaluation. This is depicted in Figure 16. Since we only introduce an error in the first objective by the reduced order model, the error results in a horizontal shift of the Pareto points. In general, we observe a good agreement between the solutions, the error being less than for all points that were tested. However, especially in regions with a steep gradient