Model Reduction For Parametrized Optimal Control Problems in Environmental Marine Sciences and Engineering

# Model Reduction For Parametrized Optimal Control Problems in Environmental Marine Sciences and Engineering

Maria Strazzullo mathLab, Mathematics Area, SISSA, Via Bonomea 265, I-34136 Trieste, Italy Francesco Ballarin mathLab, Mathematics Area, SISSA, Via Bonomea 265, I-34136 Trieste, Italy Renzo Mosetti National Institute of Oceanography and Experimental Geophysics, OGS, Via Beirut 2, I-34151 Trieste, Italy Gianluigi Rozza mathLab, Mathematics Area, SISSA, Via Bonomea 265, I-34136 Trieste, Italy
###### Abstract

In this work we propose reduced order methods as a suitable approach to face parametrized optimal control problems governed by partial differential equations, with applications in environmental marine sciences and engineering. Environmental parametrized optimal control problems are usually studied for different configurations described by several physical and/or geometrical parameters representing different phenomena and structures. The solution of parametrized problems requires a demanding computational effort. In order to save computational time, we rely on reduced basis techniques as a reliable and rapid tool to solve parametrized problems. We introduce general parametrized linear quadratic optimal control problems, and the saddle-point structure of their optimality system. Then, we propose a POD-Galerkin reduction of the optimality system. Finally, we test the resulting method on two environmental applications: a pollutant control in the Gulf of Trieste, Italy and a solution tracking governed by quasi-geostrophic equations describing North Atlantic Ocean dynamic. The two experiments underline how reduced order methods are a reliable and convenient tool to manage several environmental optimal control problems, for different mathematical models, geographical scale as well as physical meaning.

Keywords: reduced order methods, proper orthogonal decomposition, parametrized optimal control problems, PDEs state equations, environmental marine applications, quasi-geostrophic equation.

AMS: 49J20, 76N25, 35Q35

## 1 Introduction

Parametrized optimal control problems (OCP()s) governed by parametrized partial differential equations (PDE()s) are usually complex and demanding, computationally speaking. In this case the parameter could represent physical or geometrical features. In order to study different configurations, a rapid and suitable approach based on reduced order models could allow to face OCP()s in a low dimensional framework. Computational methods for OCP()s are a quite widespread tool in many contexts and fields: in shape optimization (see e.g. [16, 20, 27]), in flow control (see e.g. [12, 13, 30, 32]), in environmental applications (see e.g. [15, 33, 34]). Reduced order methods are a strategy that decreases the computational costs and simplifies the resolution of simulations governed by partial differential equations (see e.g. [21, 35, 36, 26]).
In this work, we focus on reduced order modelling for optimal control problems with quadratic cost functional constrained to linear PDE()s dealing with applications in environmental marine sciences. Optimal control theory fits well in these fields, as it is at the basis of forecasting models and it could be used to make previsions on several natural phenomena (see e.g. [18, 23, 44]). In order to describe different configurations, simulations have to be run for several values of . For this reason reduced order methods are presented as an important resource to manage these problems. The main novelty of this work deals with the use of consolidated reduced order modelling techniques in geographically realistic experiments involved in marine sciences with environmental purposes. Two examples will be discussed:

1. a first example dealing with a pollutant control in the Gulf of Trieste, Italy. In this case we aimed at underlining how parametrized optimal control problems could be useful in order to study several potential configurations describing different phenomena in this specific geographic area. These experiments could improve the monitoring strategies (presented, for example, in [28, 41]) of the marine environment of the Gulf of Trieste;

2. an Oceanographic solution tracking governed by quasi-geostrophic equations (see e.g. [10, 24, 38]), describing the North Atlantic Ocean dynamics. This problem could be linked to data assimilation based forecasting models (see e.g. [4, 8, 9, 18, 23, 44]) which are used to simulate possible scenarios and to analyse and predict climatological phenomena.

One of the main purposes of this work is to present reduced order methods as a reliable and useful tool to manage environmental simulations, both for problems characterized by large scales as well as small ones. Indeed, the two experiments that we have analysed are different under many aspects: equations used, scale analysis, geographical regions considered, dynamics involved. Despite that, they have in common a physical parametrized setting describing several configurations and modelling different natural phenomena. Many computational resources are needed to solve an optimal control problem in environmental sciences, most of all when parameters-dependent simulations have to be run many times in order to study and analyse several results, representing very different physical and natural aspects. Reduced order methods allow to recast a computational demanding problem, the “truth” problem, into a new fast and reliable low-dimensional formulation, the reduced problem. The latter is formulated as a Galerkin projection into reduced spaces, generated by basis functions chosen through a proper orthogonal decomposition sampling algorithm, as presented in [1, 21].
This work will show how convenient the reduced formulation is, since it is capable to give real-time results, while environmental simulations based on classical approximation methods may take a very long time. The computational time saved for the reduced simulations could be invested in the study of many scenarios in order to achieve a deeper knowledge of ecological and climatological phenomena, that are very thorny to be analyzed and understood, since they are linked to several natural aspects, as well as anthropic features.
The paper is outlined as follows. In section 2, first of all the saddle-point structure for linear quadratic parametrized optimal control problems is briefly discussed, then its Finite Element (FE) approximation is presented (see [17, 19]). Section 3 aims at introducing the reduced order approximation (following [21, 22]) and POD sampling algorithm for OCP()s with a brief mention of aggregated reduced space strategy (used in [14, 30, 31]) and affine decomposition (see e.g. ). Finally, in section 4, numerical results dealing with reduced order methods applied to environmental marine sciences are detailed. Conclusions follow in section 5.

## 2 Problem Formulation and Finite Element Approximation

This section aims at describing linear quadratic OCP()s exploiting their saddle-point formulation. This strategy is illustrated in many works and in many applications (see [30, 31, 37, 39]). Saddle-point formulation is advantageous, since the results about its well-posedness are well known in literature (see [6, 7]). Then, we will briefly focus on the Finite Element (FE) “truth” approximation of OCP()s.

In the treatment of linear quadratic OCP()s and their saddle-point structures we will essentially follow [29, 30, 31]. Let us consider an open and bounded domain with Lipschitz boundary . Let and be Hilbert spaces for state and control, respectively. Moreover, let be an Hilbert space, in which the observation is taken. Furthermore, let us define a compact set of parameters , for , where is considered. Taking into account the adjoint Hilbert space , the linear constraint equation is defined by:

 a(y,q;μ)=c(u,q;μ)+⟨G(μ),q⟩∀q∈Q, (2.1.1)

where represents the state operator, describes the role of the control in the formulation of the problem and gives information about forcing and boundary terms. Given a constant , the quadratic objective functional is given by:

 J(y,u;μ)=12m(y−yd(μ),y−yd(μ);μ)+α2n(u,u;μ), (2.1.2)

where , and are bilinear forms representing the objective on the state variable and a penalization for the control variable, respectively.
An OCP() problem can be formalized as follows: given , solve

 min (y(μ),u(μ))∈Y×UJ(y(μ),u(μ);μ)such that (y(μ),u(μ))∈Y×U % satisfies (???). (2.1.3)

In order to recast the problem (2.1.3) in a saddle-point formulation, let us define . Being and two elements of , we can endow with the scalar product and with the norm . Now let us consider

 A:X×X→R A(x,w;μ)=m(y,z;μ)+αn(u,v;μ) ∀x,w∈X, B:X×Q→R B(w,q;μ)=a(z,q;μ)−c(v,q;μ) ∀w∈X,∀q∈Q, F(μ)∈X∗ ⟨F(μ),w⟩=m(yd(μ),z;μ) ∀w∈X.

Thanks to these quantities, the following functional can be defined:

 J(x;μ)=12A(x,x;μ)−⟨F(μ),x⟩.

In [5, 31], it is shown that minimizing with respect to is equivalent to minimize with respect . So problem (2.1.3) is equivalent to

 (2.1.4)

The constrained optimization problem (2.1.4) can be recast into an unconstrained optimization problem by defining the Lagrangian functional as

 L(x,p;μ)=J(x;μ)+B(x,p;μ)−⟨G(μ),p⟩, (2.1.5)

where is the adjoint variable. It is well known that the minimization problem (2.1.4) is equivalent to find the saddle-point of the Lagrangian Functional (2.1.5) and this leads to the saddle-point formulation (see [5, 39] as references).
The new formulation of the problem is: given , find such that

 {A(x(μ),w;μ)+B(w,p(μ);μ)=⟨F(μ),w⟩∀w∈X,B(x(μ),q;μ)=⟨G(μ),q⟩∀q∈Q. (2.1.6)

If , the existence and the uniqueness of the solution could be provided (see [5, 29, 31]) by the fulfillment of the hypotheses of the Brezzi’s theorem reported in .

### 2.2 Finite Element “Truth” Approximation of OCP(μ)s

Let be a triangulation over . In a Finite Element approximation and , where

 XNr={vN∈C0(¯¯¯¯Ω):vN|K∈Pr,∀K∈TN}.

and represents the space of polynomials of degree at most equal to and a triangle of . To have Brezzi’s hypotheses guaranteed, we assume also that (see [11, 29, 31]). Let us consider the discrete product space . The Galerkin Finite Element discretization of the saddle-point problem (2.1.6) is: given , find such that

 {A(xN(μ),wN;μ)+B(wN,p(μ)N;μ)=⟨F(μ),wN⟩∀vN∈XN,B(xN(μ),qN;μ)=⟨G(μ),qN⟩∀qN∈QN. (2.2.1)

Let us focus on the algebraic structure of the system associated to (2.2.1). The dimension of and are respectively indicated with and . Let us define the basis of the finite dimensional spaces and respectively as:

 {φi}NXi=1{ψj}NQj=1.

We now can rewrite the solution as:

 (xN(μ)=NX∑i=1xμiφi,pN(μ)=NQ∑j=1pμjψj).

Let us define and as follows:

 Aij(μ)=A(φi,φj;μ),Bml(μ)=B(φl,ψm;μ),Fk(μ)=⟨F(μ),φk⟩,Gs(μ)=⟨G(μ),ψs⟩.

From those quantities, we can build the following linear system, with a block structure:

 (A(μ)BT(μ)B(μ)0)(xμpμ)=(F(μ)G(μ)), (2.2.2)

where and . The approach proposed is known as optimize-then-discretize (see [17, 19]). In all the applications presented, the linear system is solved through one-shot method (see [40, 42]).

## 3 Reduced Basis Methods for Parametrized Optimal Control Problems

In this section reduced basis methods for OCP()s are described. First of all, the general idea of reduced order approximation for OCP()s is given as proposed in [15, 30, 31]. Then, proper orthogonal decomposition (POD, see [1, 21] as references) and the theory of aggregated spaces will be introduced with some considerations about the efficiency of the method through affinity assumption, as presented in .

### 3.1 Problem Formulation and Solution Manifold

In section 2.1, we have already affirmed that a linear quadratic OCP()s could be formulated as a saddle-point problem of the form (2.1.6). Let us recall that . In several cases, the association defines a smooth solution manifold of the form:

 M={(x(μ),p(μ))|μ∈P}.

When the full order problem (2.2.1) is solved, one finds the approximated solution manifold:

 MN={(xN(μ),pN(μ))|μ∈P}.

Reduced basis methods aim at building a good approximation of through linear combination of properly chosen snapshots and , assuming that the approximated manifold has a smooth dependence from . In other words, the reduced spaces are built with full order solutions computed for suitable parameters in . Let us suppose to have already built and as reduced product space and reduced adjoint space, respectively (the reduced spaces will be specified in section 3.3). Than, the reduced problem is formulated as follows: given , find such that

 {A(xN(μ),wN;μ)+B(wN,pN(μ);μ)=⟨F(μ),wN⟩∀wN∈XN,B(xN(μ),qN;μ)=⟨G(μ),qN⟩∀qN∈QN. (3.1.1)

### 3.2 POD Algorithm for OCP(μ)s

Let us focus our attention on POD algorithm used as sampling procedure for the construction of the reduced bases. In order to apply it, a discrete and finite dimensional subset is needed. For this specific set of parameters, the discrete solution manifold is defined as:

 MN(Ph)={(xN(μ),pN(μ))|μ∈Ph}.

The cardinality of is . Naturally it holds since . If is fine enough, is a good approximation of the discrete manifold . From now on, we will refer to the set of the linear combinations of elements of as . The algorithm of POD is based on two processes:

1. sampling the parameter space in order to compute the full order solutions at selected parameters,

2. a compression phase, where one discards the redundant information, respectively for state, control and adjoint variables.

The -spaces resulting from the POD algorithm minimize the following quantities, respectively:

 √1M∑μ∈Phinf zN∈YN∥yN(μ)−zN∥2Y,√1M∑μ∈Phinf vN∈UN∥uN(μ)−vN∥2U,√1M∑μ∈Phinf qN∈QN∥pN(μ)−qN∥2Q. (3.2.1)

Let us introduce an ordering on the parameters . This induces an ordering on the full order solutions , and . In order to build the POD-spaces, we define the symmetric and linear operators:

 Cy:M→M Cy(zN)=1MM∑m=1(zN,yN(μm))yN(μm), zN∈M, Cu:M→M Cu(vN)=1MM∑m=1(vN,uN(μm))uN(μm), vN∈M, Cp:M→M Cp(qN)=1MM∑m=1(qN,pN(μm))pN(μm), qN∈M.

Let us consider their eigenvalues and the corresponding eigenfunctions , with , verifying:

Let us assume that the eigenvalues satisfy , for . The orthogonal POD basis functions are given by , and and they span . We can take into consideration the first eigenfunctions for the sake of reduction, respectively for state, control and adjoint space: the reduced spaces and will be defined by them.

###### Remark 3.2.1

The POD algorithm can be also seen under an algebraic point of view. For example, let us consider the control variables111The concept proposed is easily extended to state and adjoint variables, analogously. for . Let be the correlation matrix of the control snapshots, that is:

 Cumq=1M(uN(μm),uN(μq))U,1≤m,q≤M.

Then, the -largest eigenvalue-eigenvector pairs solve the problem

 Cuvn=λnvn,1≤n≤N,

with . Giving a descending order to the eigenvalues , the orthogonal basis functions satisfy . The basis is given by:

 ξun=1√MM∑m=1(vn)muN(μm),1≤n≤N,

where is m-th component of the control eigenvector .

### 3.3 Aggregated Reduced Spaces and Affinity Assumption

We now focus on the conditions needed to guarantee stability and efficiency of the proposed reduced order method. In order to prove the well-posedness of the problem (3.1.1), the reduced inf-sup condition of the bilinear form must be fulfilled, in other words, it must exists a constant such that

 βN(μ)=infqN∈QNsupwN∈XNB(wN,qN;μ)∥wN∥X∥qN∥Q>βN0∀μ∈P. (3.3.1)

In order to ensure the hypothesis (3.3.1), must be assumed. As underlined in , simply building the reduced spaces as linear combinations of snapshots may not lead to the fulfillment of the reduced inf-sup condition. Then, we adopted the solution of aggregated spaces, used in [30, 31], already presented in . This technique is based on the definition of an enriched space

 ZN=span {yN(μn),pN(μn),n=1,…,N}.

Now, let us define the reduced spaces and , where

 UN=span {uN(μn),n=1,…,N}.

Thanks to this choice, the hypothesis is recovered and so the saddle-point problem (3.1.1) verifies the reduced inf-sup condition.
Let us briefly introduce the affinity assumption222If the problem does not fulfill the affinity assumption, it can be recovered thanks to the empirical interpolation method (see e.g.  and [21, Chapter 5]). that guarantees efficiency of reduced order methods. The problem (2.1.6) admits affine decomposition if we can rewrite the bilinear forms and the functionals involved as:

 A(x,w;μ)=QA∑q=1ΘqA(μ)Aq(x,w),B(w,p;μ)=QB∑q=1ΘqB(μ)Bq(w,p),⟨G(μ),s⟩=QG∑q=1ΘqG(μ)⟨Gq,s⟩,⟨F(μ),w⟩=QF∑q=1ΘqF(μ)⟨Fq,w⟩, (3.3.2)

for some finite , where are dependent smooth functions, whereas ,
are independent bilinear forms and functionals. This hypothesis allows us to divide the resolution of the reduced order approximation of (2.1.6) in two stages:

1. offline: in this stage the reduced spaces are built and all the independent quantities are assembled. It is performed only once and it may be very costly;

2. online: in this phase the dependent quantities are assembled and the reduced system is solved. This stage is performed every time we want the model to be simulated at a new value of , representing a new configuration for our system.

## 4 Applications in Environmental Marine Sciences and Engineering

This section aims at applying the proposed reduced order method (ROM) to parametrized optimal control problems involved in environmental marine sciences and engineering problems. One of the purpose is to show the computational savings enabled by the use of a ROM in place of the usual FE approximation strategies. Two specific examples are proposed:

1. A Pollutant Control on the Gulf of Trieste
This first example involves an advection-diffusion pollutant control problem set in the Gulf of Trieste, Italy. The latter is a physical basin particularly windy and it has very peculiar flora and fauna population (as underlined in [28, 41]). Moreover its analysis is important from a social point of view since it has a great impact on the local community: the city of Trieste overlooks the sea and depends on the Gulf and on its structures from harbours as well as from tourist infrastructures. For these reasons it needs to be monitored and kept under control.

2. A Solution Tracking of the Large Scale Ocean Circulation Model
The Solution tracking is an optimal control problem that aims at making a solution the most similar to a given observation. As a second application, we propose a solution tracking problem of the large scale Ocean circulation model, governed by quasi-geostrophic equations. This OCP() example fits in the framework of a data assimilation technique (see [4, 8, 9, 18] as references) that allows the model to be modified adding information from experimental data. The importance of studying Ocean Circulations Models lies in forecasting analysis of future meteorological and climatological scenarios in order to prevent catastrophic events, as presented in .

Both the presented applications are characterized by several physical parameters, and so reduced order methods could be an useful tool to decrease the time required by numerical simulations.
In order to reach more realistic results, we have used meshes derived from satellite images representing the geographic area to be studied, as shown in Figure 4.0.1. Figure 4.0.2 gives an idea about the work needed in order to build realistic meshes: some details of the meshes overlapped to the satellite images are shown. For our analysis, having these specific meshes was very important to give physical meaning to our problems and to their formulations, in order to achieve reliable results that could be potentially compared to real experimental data, which are strictly linked to the geographical area where they are collected. Figure 4.0.1: Satellite images. (a) North Atlantic Ocean, West coast; (b) North Atlantic Ocean, East coast; (c) Gulf of Trieste, Italy. Figure 4.0.2: Meshes overlapping satellite images. (a) North Atlantic Ocean, particular of the Florida peninsula; (b) North Atlantic Ocean, particular of Portugal and North Africa; (c) Gulf of Trieste, particular of the city of Trieste.

The simulations have been run using FEniCS  for the full order solutions and RBniCS for the reduced order ones [21, 2]. The machine used for the simulations has a processor AMD A8-6410 APU with 8 GB of RAM.

### 4.1 Reduced Basis Applied to a Pollutant Control on the Gulf of Trieste

The proposed problem aims at limiting the impact of a pollutant tracer on touristic and natural areas of the Gulf of Trieste. The OCP() is governed by advection-diffusion state equation (see [15, 33, 34]). Let be an open, bounded and regular domain representing the Gulf of Trieste (see Figure 4.1.1 (b)), with boundary and , where homogeneous Dirichlet and Neumann boundary conditions are imposed on and , respectively. The coasts are considered in , while the open sea represents (see Figure 4.1.1 (c)).
Let us define the state and the control spaces and the adjoint space .
The non-dimensional OCP() reads: given , find such that:

 min (y,u)∈Y×UJ(y,u)=min (y,u)∈Y×U12∫ΩOBS(y−yd)2dΩOBSsuch that a(y,q;μ)=c(u,q),∀q∈Q. (4.1.1)

where the state is the pollutant concentration and represents the safety threshold of the pollutant tracer. (a) Boundaries: light blue: open sea; brown: coasts. Subdomains: red: monitored natural area; green: source of the pollutant tracer.

The bilinear forms and are defined as:

 a(y,q,μ)=∫Ω(ν(μ)∇y⋅∇q+β(μ)⋅∇yq)dΩ, c(u,q)=Lu∫ΩuqdΩu.

where represents the diffusivity action of the state equation, while is the transport field acting on the Gulf. Then, the parameter influences the circulation of the current in the Gulf. In this case is . The constant makes the system non-dimensional. For the transport field we decided to take into consideration in proximity of the observation domain333Constant transport field will be sufficient to simulate the most interesting configurations for the circulation dynamic of the Gulf of Trieste.

 β1(μ2)≡μ2,β2(μ3)≡μ3.

The parameter space considered is . The plot of in Figure 4.1.1 (c) shows the considered subdomains: in green , where the pollutant loss is, in red , where we want to monitor the pollutant concentration: it represents the swimming touristic area of the city and Miramare natural area (red subdomain in Figure 4.1.1). Two argumentations drove us in the choice of (as underlined in ):

1. its unique ecological flora and fauna marine population,

2. it is an area crowded by Trieste citizens and by many tourists.

The parameter is a physical parameter that describes the dynamic of the currents deriving from the specific winds blowing on this geographical area (the winds acting on the Gulf will be introduced later). Varying the parameter, it is possible to simulate several configurations in order to study how the wind could affect the diffusion of a dangerous pollutant in the natural area of Miramare. Thanks to reduced order modelling, many scenarios could be analysed, with a great saving of computational time resources. Several possible results must be taken into consideration in this context and, then, a deeper analysis can be made in order to better understand possible scenarios. The monitoring of the diffusion of a pollutant is necessary in order to preserve natural areas or to safeguard an ecological polluted habitat in case of ecological accident, and an accurate and fast model is helpful in planning a program of action.
In order to recast the problem in the framework (2.1.6), let be the product space of state and control spaces. Let and be elements of , whereas an element of . Moreover, we define the bilinear forms and as follows:

 m(y,z)=∫ΩOBSyzdΩOBS,n(u,v)≡0.

Furthermore, we define the forms and as follows:

 A:X×X→R A(x,w)=m(y,z), B:X×Q→R B(w,q;μ)=a(z,q;μ)−c(v,q), F:X→R ⟨F,w⟩=yd∫ΩOBSzdΩOBS.

To build the aggregated reduced spaces for state and adjoint we used the POD-Galerkin algorithm introduced in section 3. In this specific example the reduced space for the control does not need to be reduced and thus we set . For this problem the affinity assumption is guaranteed: with , and the affine decomposition of the problem is given by

 Θ1A=1 A1(x,w)=A(x,w), Θ1B=μ1 B1(x,q)=∫Ω∇y⋅∇qdΩ, Θ2B=μ2 B2(x,q)=∫Ω∂y∂x1qdΩ, Θ3B=μ3 B3(x,q)=∫Ω∂y∂x2qdΩ, Θ1F=1 ⟨F1,w⟩=⟨F,w⟩.

In the following, some numerical results linked to two different configurations are shown: the data of the experiments are reported in Table 4.1. We analysed how the wind action could influence pollutant diffusion. We simulated the net water transport due to Bora, a wind blowing from East to North-West, described by , and Scirocco, a wind blowing from South-East, corresponding to . From the lower value of the cost functional reported in Table 2, one can deduce that Bora makes the polluted water removed from , while Scirocco acts is the opposite way, pushing it towards the coast. Furthermore, Table 2 shows the differences between FE and ROM performances, in terms of dimension of the systems, time of resolution and cost functional. In Table 3 the speed up index behaviour with respect to the increasing of the basis numbers is presented. The spped up index represents the number of reduced problems solved in the time needed for a full order simulation. With the terminology basis number we refer to as the number such that (and , being ), where are the number of online degrees of freedom for state, control and adjoint variable, respectively. In other words if the basis number is , we are solving a system of dimension .
The top left plot in Figure 4.1.2 shows an “uncontrolled” Bora configuration (corresponding to ), resulting from the simulation of the state equation (only) with a value as a forcing term. The top right and the bottom left plots in 4.1.2 represent the optimal control problem with FE discretization and ROM, respectively. As one can see, the FE approximation and the ROM one match. Another proof of the reliability of reduced basis POD-Galerkin method is the pointwise error shown in the bottom right of the same Figure: the maximum value reached is for .

In Figure 4.1.3 errors444The errors considered for state, control and adjoint are, respectively: between FE solution and ROM solution with respect to the basis number are presented. Some plots are cut before , since the error would be below machine epsilon for a larger value of . The obtained results show that the ROM allows to get fast and accurate simulations, since very few basis functions are required to have very small errors.
The bottom right plot illustrates a comparison between the POD approach presented in section 3 (labeled as partitioned) and a different reduction in which only one POD is carried out for the variables at the same time on the space (labeled as monolithic). The results show that it is preferable to use a partitioned POD approach rather than a monolithic POD algorithm. The improvement is significant, since the latter approach gave an error of the order of with , while, for the same value of and with the partitioned option, the sum of the state, control and adjoint errors reaches the value of . Since the partitioned approach gave better values of the errors, we decided to exploit it for the Oceanographic application which we will present in the next section.

### 4.2 Reduced Basis Applied to an Ocean Circulation Solution Tracking

The general Ocean circulation model describes large scale flow dynamics. It is strictly linked to wind action: it can be considered as a coupled system of Ocean and Atmosphere. The theory associated to this topic is deeply analysed in [10, Chapter 3]. The model is governed by the following non-dimensional PDE(), known as quasi-geostrophic equation:

 μ3F(ψ,Δψ)+∂ψ∂x=f−μ1Δψ+μ2Δ2ψ,

where, given a suitable spaces , the non-linearity of the expression is given by defined as:

 F(ψ,q)=∂ψ∂x∂q∂y−∂ψ∂y∂q∂x.

The physical parameter is considered in the parametrized space . The parameter represents how much the non-linear term affects the flow dynamics, while and stand for diffusive action, respectively.
The parameter describes the North Atlantic Ocean dynamics completely, since it gives information about how the large scale Ocean circulation is affected by different phenomena, such as location (typically described by ) and intensity variations of the gyres and the currents in the Ocean. Let us recall that Ocean dynamic is strongly linked to wind stress and atmospheric behaviour. Then, the parameter describes the dynamic of a very complex physical system, taking into account several natural factors and phenomena. It is very important to run many simulations for different values of the parameter , in order to achieve a full knowledge of this system describing Oceanic and Atmospheric dynamics, linked to climatological forecasting issues.
As specified in [10, Section 3.2], the forcing term depends on wind stress by the following relation:

 f=^k⋅rot τ,

where is the third reference spatial unit vector. In our application we considered a bounded and regular bi-dimensional domain555Experiments showed that the analysis of the underwater depth did not affect the dynamics of the equation, so we decided to exploit a simpler bi-dimensional model. Although, adding the bathimetric effect is quite simple: let us suppose that the ocean floor is described by a smooth function , then one can consider in order to treat a more complete model. As one can see, the bathimetry affects only the non-linear term. representing the North Atlantic Ocean (see Figure 4.1.1 (a)).

Furthermore, we focused on the linear version of the quasi-geostrophic equation (). The solution tracking problem constrained to this particular state equation is:

 min (ψ,u)∈Y×UJ(ψ,u)=min (ψ,u)∈Y×U12∫Ω(ψ−ψd)2dΩ+α2∫Ωu2dΩ (4.2.1) such that ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂ψ∂x=u−μ1Δψ+μ2Δ2ψin Ω,ψ=0on ∂Ω,Δψ=0on ∂Ω,

where is our state variable, is the forcing term to be controlled, where and are two suitable functions spaces and is the penalization term. The problem aims at making the solution the most similar to .
In some applications can represents experimental data. In this sense, our experiment could be seen as a prototype of a data assimilation model with forecasting purposes (see [4, 8, 9, 18, 23, 44] as references). This particular technique changes the model in order to achieve a solution comparable with real experimental data. Data assimilation techniques are very costly, and reduced order methods fit very well in this context. The proposed experiment helps to understand how reduced order techniques could be exploited in order to simulate several climatological scenarios in a low dimensional and accurate framework. The opportunity of running the reduced model many times allow us to have a deeper comprehension of the Ocean dynamic, and of atmospheric phenomena and climatological scenarios.
In order to manage a handier problem666This version of the problem does not ensure the coercivity of the state equation, that is proved in  for the state equation of the system (4.2.1), we rewrite the previous system as it follows:

 min ((ψ,q),u)∈Y×UJ((ψ,q),u)=min ((ψ,q),u)∈Y×U12∫Ω(ψ−ψd)2dΩ+α2∫Ωu2dΩ (4.2.2) such that ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩q=Δψin Ω,∂ψ∂x=u−μ1q+μ2Δqin Ω,ψ=0on ∂Ω,q=0on ∂Ω,

where the spaces are defined as and . The weak formulation reads as:

 a((ψ,q),(ϕ,r);μ)=c(u,(ϕ,r))∀ϕ,r∈H10(Ω), (4.2.3)

where and are given by:

 a((ψ,q),(ϕ,r);μ) =∫Ω∂ψ∂xrdΩ+μ2∫Ω∇q⋅∇rdΩ+ +μ1∫ΩqrdΩ+∫ΩqϕdΩ+∫Ω∇ψ⋅∇ϕdΩ, (4.2.4) c(u,(ϕ,r)) =∫ΩurdΩ.

In this case is .
Since we are facing a linear quadratic optimal control problem, it can be recast in saddle-point formulation (2.1.6). Let us define the product space and let and be two elements of , whereas let be an element of the adjoint space (). Furthermore, one has to specify:

 A:X×X→R A(x,w)=m((ψ,q),(χ,t))+αn(u,v), B:X×Q→R B(w,s;μ)=a((χ,t),(ϕ,r),μ)−c(v,(ϕ,r)), F:X→R ⟨F,w⟩=∫ΩψdχdΩ,

where and are defined as

 m((ψ,q),(χ,t))=∫ΩψχdΩ,n(u,v)=∫ΩuvdΩ.

In order to build the aggregated reduced spaces of the type proposed in section 3.3 a POD-Galerkin algorithm has been exploited. As we underlined in section 3.3, the affinity assumption must be guaranteed for the efficiency of the reduced problem. Indeed, with , and the affine decomposition of the problem is given by

 Θ1A=1 A1(x,w)=A(x,w), Θ1B=μ1 B1(x,s)=∫ΩqrdΩ, Θ2B=μ2 B2(x,s)=∫Ω∇q⋅∇rdΩ, Θ2B=1 B3(x,s)=∫Ω∂ψ∂xrdΩ+∫ΩqϕdΩ+∫Ω∇ψ⋅∇ϕdΩ−∫ΩurdΩ, Θ1F=1 ⟨F1,w⟩=⟨F,w⟩.

In Table 4 the data of the experiment are shown. In Figure 4.2.1 the desired value to be reached is presented (top left plot) with the FE and ROM solutions (top right and bottom left plot, respectively). The approximated solutions match. In the bottom right of the figure, the pointwise error is shown: the maximum value reached is .

The following tables represent the comparison between ROM and FE performances, in terms of: system dimension, cost functional optimal value, time of resolution (Table 5) and speed up index with respect to the basis number such that and (Table 6), which results in a solution of a linear system. We can deduce how ROM method is a suitable and convenient approach to study large scale phenomena in oceanography, a field dealing with parametrized simulations that require days of CPU times for complex configurations.

In Figure 4.2.2 the error norm777The errors considered for state, control and adjoint are, respectively: between the FE approximation and the ROM discretization is shown for all the variables, state and , control , and adjoint and , respectively. For the state and the adjoint , we plot the error until , the value above which the error goes below machine precision.

## 5 Conclusions and Perspectives

In this work we have exploited reduced order methods in environmental parametrized optimal control problems dealing with marine sciences and engineering. We proposed two specific examples: one representing the ecological issue of pollutant control in a specific naturalist area, the Gulf of Trieste, Italy, the other one consisting in a large scale solution tracking OCP() governed by quasi-geostrophic equations. We showed how reduced order methods could be a very useful tool in environmental sciences, like oceanography and ecology, where parametrized simulations are usually very demanding and costly. Reduced order methods are a suitable approach to face these issues. We have used a POD-Galerkin method on a saddle-point formulation for sampling and for the projection, by exploiting an aggregated space strategy. Reduced order methods performances have been compared to FE approximation, classically used to study these phenomena, in order to prove how convenient the reduced order approach could be in this particular field of applications: in the North Atlantic problem, the reduced time of resolution decreases of one order of magnitude with respect to the full order one, and the error of the state variable is negligible. To the best of our knowledge, the main novelty of this work is in the POD-Galerkin reduction of a solution tracking optimal control problem governed by quasi-geostrophic equations.
Let us expose some improvements of this work, focusing on the optimal control problem governed by quasi-geostrophic equation. A possible development would involve a time dependent optimal control problem considering also the non-linear case. This kind of formulation is of the utmost importance in climatological applications, in order to forecast and predict possible scenarios in a reliable way. Time dependency and non-linearity will make the problems more and more realistic and suited to actual ecological and climatological challenges, as well as more and more computational demanding. In this sense, reduced order modelling appears, again, to be a suitable and versatile approach to be used. Time dependent non-linear optimal control problems insert themselves in the framework of data assimilation techniques, that, as briefly introduced in section 4.2, allow to modify the model in order to reach more reliable results in the forecasting applications, thanks to a solution tracking where the solution desired to be reached represents real experimental data.
For both the examples presented, a further step could be the development of three-dimensional marine model that could take into consideration bathimetry effect. Finally, the problems could be inserted in a reduced order uncertainty quantification context (see e.g. ), when it is not possible to assign specific values for the parameters by classical statistical methods.

## Acknowledgements

We acknowledge the support by European Union Funding for Research and Innovation – Horizon 2020 Program – in the framework of European Research Council Executive Agency: Consolidator Grant H2020 ERC CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics”. We also acknowledge the INDAM-GNCS project “Metodi numerici avanzati combinati con tecniche di riduzione computazionale per PDEs parametrizzate e applicazioni”.

## References

•  F. Ballarin, A. Manzoni, A. Quarteroni, and G. Rozza. Supremizer stabilization of POD–Galerkin approximation of parametrized steady incompressible Navier–Stokes equations. International Journal for Numerical Methods in Engineering, 102(5):1136–1161, 2015.
•  F. Ballarin, A. Sartori, and G. Rozza. RBniCS - reduced order modelling in FEniCS.
•  M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera. An Empirical Interpolation Method: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathematique, 339(9):667–672, 2004.
•  D. W. Behringer, M. Ji, and A. Leetmaa. An improved coupled model for ENSO prediction and implications for Ocean initialization. Part I: The Ocean data assimilation system. Monthly Weather Review, 126(4):1013–1021, 1998.
•  P. B. Bochev and M. D. Gunzburger. Least-squares finite element methods, volume 166. Springer-Verlag, New York, 2009.
•  D. Boffi, F. Brezzi, and M. Fortin. Mixed finite element methods and applications, volume 44. Springer-Verlag, Berlin and Heidelberg, 2013.
•  F. Brezzi. On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers. Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique, 8(2):129–151, 1974.
•  J. A. Carton, G. Chepurin, X. Cao, and B. Giese. A simple Ocean data assimilation analysis of the global upper Ocean 1950–95. Part I: Methodology. Journal of Physical Oceanography, 30(2):294–309, 2000.
•  J. A. Carton and B. S. Giese. A reanalysis of ocean climate using Simple Ocean Data Assimilation (SODA). Monthly Weather Review, 136(8):2999–3017, 2008.
•  F. Cavallini and F. Crisciani. Quasi-geostrophic theory of Oceans and atmosphere: topics in the dynamics and thermodynamics of the Fluid Earth, volume 45. Springer Science & Business Media, New York, 2013.
•  P. Chen, A. Quarteroni, and G. Rozza. Multilevel and weighted reduced basis method for stochastic optimal control problems constrained by Stokes equations. Numerische Mathematik, 133(1):67–102, 2015.
•  J. C. de los Reyes and F. Tröltzsch. Optimal control of the stationary Navier-Stokes equations with mixed control-state constraints. SIAM Journal on Control and Optimization, 46(2):604–629, 2007.
•  L. Dedè. Optimal flow control for Navier-Stokes equations: Drag minimization. International Journal for Numerical Methods in Fluids, 55(4):347–366, 2007.
•  L. Dedè. Reduced basis method and a posteriori error estimation for parametrized linear-quadratic optimal control problems. SIAM Journal on Scientific Computing, 32(2):997–1019, 2010.
•  L. Dedè. Adaptive and reduced basis method for optimal control problems in environmental applications. PhD thesis, Politecnico di Milano, 2008. Available at http://mox.polimi.it.
•  M. C. Delfour and J. Zolésio. Shapes and geometries: metrics, analysis, differential calculus, and optimization, volume 22. SIAM, Philadelphia, 2011.
•  E. Fernández Cara and E. Zuazua Iriondo. Control theory: History, mathematical achievements and perspectives. Boletín de la Sociedad Española de Matemática Aplicada, 26, 79-140., 2003.
•  M. Ghil and P. Malanotte-Rizzoli. Data assimilation in meteorology and oceanography. Advances in geophysics, 33:141–266, 1991.
•  M. D. Gunzburger. Perspectives in flow control and optimization, volume 5. SIAM, Philadelphia, 2003.
•  J. Haslinger and R. A. E. Mäkinen. Introduction to shape optimization: theory, approximation, and computation. SIAM, Philadelphia, 2003.
•  J. S. Hesthaven, G. Rozza, and B. Stamm. Certified reduced basis methods for parametrized partial differential equations. SpringerBriefs in Mathematics, 2015, Springer, Milano.
•  K. Ito and SS. Ravindran. A reduced-order method for simulation and control of fluid flows. Journal of computational physics, 143(2):403–425, 1998.
•  E. Kalnay. Atmospheric modeling, data assimilation and predictability. Cambridge university press, 2003, Cambridge.
•  T.Y. Kim, T. Iliescu, and E. Fried. B-spline based finite-element method for the stationary quasi-geostrophic equations of the Ocean. Computer Methods in Applied Mechanics and Engineering, 286:168–191, 2015.
•  A. Logg, K.A. Mardal, and G. Wells. Automated Solution of Differential Equations by the Finite Element Method. Springer-Verlag, Berlin, 2012.
•  Y Maday, AT Patera, and G Turinici. Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods. J. Fluids Eng, 124(1):70–80, 2002.
•  B. Mohammadi and O. Pironneau. Applied shape optimization for fluids. Oxford University Press, New York, 2010.
•  R. Mosetti, C. Fanara, M. Spoto, and E. Vinzi. Innovative strategies for marine protected areas monitoring: the experience of the Istituto Nazionale di Oceanografia e di Geofisica Sperimentale in the Natural Marine Reserve of Miramare, Trieste-Italy. In OCEANS, 2005. Proceedings of MTS/IEEE, pages 92–97. IEEE, 2005.
•  F. Negri. Reduced basis method for parametrized optimal control problems governed by PDEs. Master thesis, Politecnico di Milano, 2011.
•  F. Negri, A. Manzoni, and G. Rozza. Reduced basis approximation of parametrized optimal flow control problems for the Stokes equations. Computers & Mathematics with Applications, 69(4):319–336, 2015.
•  F. Negri, G. Rozza, A. Manzoni, and A. Quarteroni. Reduced basis method for parametrized elliptic optimal control problems. SIAM Journal on Scientific Computing, 35(5):A2316–A2340, 2013.
•  M. Pošta and T. Roubíček. Optimal control of Navier–Stokes equations by Oseen approximation. Computers & Mathematics With Applications, 53(3):569–581, 2007.
•  A. Quarteroni, G. Rozza, L. Dedè, and A. Quaini. Numerical approximation of a control problem for advection-diffusion processes. In IFIP Conference on System Modeling and Optimization, pages 261–273, Ceragioli F., Dontchev A., Futura H., Marti K., Pandolfi L. (eds) System Modeling and Optimization. CSMO 2005. vol 199. Springer, Boston, 2005.
•  A. Quarteroni, G. Rozza, and A. Quaini. Reduced basis methods for optimal control of advection-diffusion problems. In Advances in Numerical Mathematics, pages 193–216. RAS and University of Houston, Moscow, 2007.
•  G. Rozza, D.B.P. Huynh, and A. Manzoni. Reduced basis approximation and a posteriori error estimation for Stokes flows in parametrized geometries: Roles of the inf-sup stability constants. Numerische Mathematik, 125(1):115–152, 2013.
•  G. Rozza, D.B.P. Huynh, and A.T. Patera. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics. Archives of Computational Methods in Engineering, 15(3):229–275, 2008.
•  G. Rozza, A. Manzoni, and F. Negri. Reduction strategies for PDE-constrained oprimization problems in Haemodynamics. pages 1749–1768, ECCOMAS, Congress Proceedings, Vienna, Austria, September 2012.
•  O. San and T. Iliescu. A stabilized proper orthogonal decomposition reduced-order model for large scale quasigeostrophic Ocean circulation. Advances in Computational Mathematics, 41(5):1289–1319, 2015.
•  J. Schöberl and W. Zulehner. Symmetric indefinite preconditioners for saddle point problems with applications to PDE-constrained optimisation problems. SIAM Journal on Matrix Analysis and Applications, 29(3):752–773, 2007.
•  V. Schulz and I. Gherman. One-shot methods for aerodynamic shape optimization. In MEGADESIGN and MegaOpt-German Initiatives for Aerodynamic Simulation and Optimization in Aircraft Design, pages 207–220. Springer, 2009.
•  T. Shiganova and A. Malej. Native and non-native ctenophores in the Gulf of Trieste, Northern Adriatic Sea. Journal of Plankton Research, 31(1):61–71, 2009.
•  S. Taasan. One shot methods for optimal control of distributed parameter systems 1: Finite dimensional control. Tech. report 91-2, Institute for Computer Applications in Science and Engineering, Hampton, VA, 1991.
•  D. Torlo, F. Ballarin, and G. Rozza. Weighted stabilized reduced basis methods for parametrized advection dominated problems with random inputs. Submitted, 2017.
•  E. Tziperman and W. C. Thacker. An optimal-control/adjoint-equations approach to studying the Oceanic general circulation. Journal of Physical Oceanography, 19(10):1471–1485, 1989.
•  H. Yang, G. Lohmann, W. Wei, M.I. Dima, M. Ionita, and J. Liu. Intensification and poleward shift of subtropical western boundary currents in a warming climate. Journal of Geophysical Research: Oceans, 121(7):4928–4945, 2016. 