Víctor Hernández-Santamaría

Martin Lazar
email: mlazar@unidu.hr

Enrique Zuazua
email: enrique.zuazua@deusto.es

Greedy optimal control for elliptic problems and its application to turnpike problems

Keywords: parameterized PDEs, optimal control, turnpike property, greedy algorithms, elliptic equations, parabolic equations. MSC2010: 49J20, 49K20, 93C20, 49N05, 65K10

We adapt and apply greedy methods to approximate in an efficient way the optimal controls for parameterized elliptic control problems. Our results yield an optimal approximation procedure that, in particular, performs better than simply sampling the parameter-space to compute controls for each parameter value. The same method can be adapted for parabolic control problems, but this leads to greedy selections of the realizations of the parameters that depend on the initial datum under consideration. The turnpike property (which ensures that parabolic optimal control problems behave nearly in a static manner when the control horizon is long enough) allows using the elliptic greedy choice of the parameters in the parabolic setting too. We present various numerical experiments and an extensive discussion of the efficiency of our methodology for parabolic control and indicate a number of open problems arising when analyzing the convergence of the proposed algorithms.

Introduction

Optimal control problems play a major role in many fields of science and engineering. These problems are commonly subject to partial differential equations (PDE) depending on several parameters. While the PDE describes the underlying system, the parameters are used to identify or specify particular configurations such as material properties, the position of sensors and actuators, initial conditions, among others. In applications, it is of interest to explore within different parameter configurations. However, from the control point of view, this leads to solving a different problem for each new desired configuration which at a computational level may be rather expensive or prohibitive.

To overcome this expensive task, different techniques have been developed in the past years in order to speed up the solution of parameterized optimal control problems. First examples of PDE constrained optimization problems solved by computational reduction techniques have been addressed in [START_REF] Ito | A reduced basis method for control problems governed by PDEs[END_REF] using reduced basis (RB) methods or [START_REF] Ravindran | A reduced-order approach for optimal control of fluids using proper orthogonal decomposition[END_REF] by means of proper orthogonal decomposition (POD).

In the context of parameterized control problems, POD techniques have been successfully applied in [START_REF] Atwell | Proper orthogonal decomposition for reduced basis feedback controllers for parabolic equations[END_REF][START_REF] Kunisch | Control of the Burgers equation by a reduced-order approach using proper orthogonal decomposition[END_REF][START_REF] Kunisch | HJB-POD-based feedback design for the optimal control of evolution problems[END_REF][START_REF] Tröltzsch | POD a-posteriori error estimates for linear-quadratic optimal control problems[END_REF]. In particular, in [START_REF] Tröltzsch | POD a-posteriori error estimates for linear-quadratic optimal control problems[END_REF], the authors estimate the distance between the optimal control and an approximation using perturbation arguments, obtaining an efficient algorithm to solve the optimal control problem. However, the evaluation of the a posteriori error bounds requires a forward-backward solution of the associated state and adjoint equations and, as pointed in [START_REF] Tröltzsch | POD a-posteriori error estimates for linear-quadratic optimal control problems[END_REF], this is computationally expensive.

In the same way, RB techniques have been used in the context of parameter-dependent control problems, see, for instance, [START_REF] Bader | Certified reduced basis methods for parametrized distributed elliptic optimal control problems with control constraints[END_REF][START_REF] Dedè | Reduced basis method and a posteriori error estimation for parametrized linear-quadratic optimal control problems[END_REF][START_REF] Ito | Reduced basis method for optimal control of unsteady viscous flows[END_REF][START_REF] Grepl | Reduced basis a posteriori error bounds for parametrized linearquadratic elliptic optimal control problems[END_REF][START_REF] Kärcher | A Posteriori Error Estimation for Reduced Order Solutions of Parametrized Parabolic Optimal Control Problems[END_REF][START_REF] Negri | Reduced basis method for parametrized elliptic optimal control problems[END_REF]. In such works, reliable error estimates of the reduced order optimal controls are presented.

More recently, in [START_REF] Lazar | Greedy controllability of finite dimensional linear systems[END_REF], the authors develop greedy methods in the context of controllability problems. Their analysis is based on the ideas of reduced modelling and (weak) greedy algorithms for parameter dependent PDE and abstract equations in Banach spaces (see, e.g., [START_REF] Buffa | A priori convergence of the greedy algorithm for the parameterized reduced basis[END_REF][START_REF] Cohen | Kolmogorov widths under holomorphic mappings[END_REF]). The method is shown to be computationally efficient, but it has the drawback that the greedy choice of the parameters depends on the data to be controlled. The problem of developing greedy methods for controllability problems in a way that the greedy choice of the parameters is independent of the data to be controlled is an open problem.

The turnpike theory (see, for instance, [START_REF] Porretta | Long time versus steady state optimal control[END_REF]) roughly states that optimal control problems, under suitable conditions, have the property that in long time-horizons, the optimal trajectories are exponentially close to the steady state ones during most of the time interval. This means that, in some sense, the time-dependent control problem allows to study the corresponding elliptic one. In spite of this simplification, when the system depends on one or more parameters, exploring many configurations becomes prohibitive and this justifies the relevance of applying an efficient approximation tool.

With this in mind, our main contributions can be summarized as follows. Firstly, we present the application of greedy and weak greedy algorithms for elliptic optimal control problems leading to optimal approximations in the sense of the Kolmogorov width. Although the optimal control problem for parameterized elliptic equations is well known and has been studied in numerous papers (see, for instance, [START_REF] Bader | Certified reduced basis methods for parametrized distributed elliptic optimal control problems with control constraints[END_REF][START_REF] Grepl | Reduced basis a posteriori error bounds for parametrized linearquadratic elliptic optimal control problems[END_REF][START_REF] Kärcher | A certified reduced basis method for parametrized elliptic optimal control problems[END_REF][START_REF] Negri | Reduced basis method for parametrized elliptic optimal control problems[END_REF]), as far as the authors knowledge, this greedy and weak greedy algorithms have not been developed in the context of optimal control problems.

On the other hand, it is well known that solutions of optimal control problems for elliptic PDEs can be characterized as the solutions of the corresponding optimality system, which is constituted by the coupling of two elliptic equations: the one fulfilled by the state and the one corresponding to the adjoint state. The second contribution consists in developing cheap surrogates allowing to diminish the computational cost in the offline part of the greedy algorithms.

Lastly, we present a detailed discussion on how the turnpike property may be used for approximating the solution of parabolic optimal control problems. From the turnpike theory we know that under some controllability assumption, the control-state pair (u(t), y(t)) converges to the steady one (ū, ȳ), regardless the stability of the time-dependent problem.

Using our greedy approach, we can compute an optimal approximation (ū , ȳ ) of (ū, ȳ) and apply the obtained time-independent approximate control ū to the corresponding parabolic equation. As observed during experiments, the efficiency of this method will largely depend on the stability of the parabolic equation. In the stable case, by applying the control ū the solution converges to ȳ , no matter what the initial datum is. This is a consequence of the stability of the system and the fact that the applied control is time-independent. On the other hand, the approximated time-independent control is not enough to control the unstable dynamics. The instability of the system imposes the use of time-dependent controls, as the ones the turnpike theory yields, in which an initial rapidly varying arc is followed by a long time period in which the control is nearly steady.

The paper is organized as follows. In section 2, we introduce the optimal control problem for a parabolic equation and then, by means of the turnpike property, we reduce this problem to study a steady state system. This will be the starting point to formulate the greedy approach for elliptic optimal control systems. In Section 3, we present the main results of this paper: firstly, we present a brief summary on greedy algorithms and then we analyze their application in the optimal control of parameterized elliptic equations. Section 4 contains several numerical examples and experiments for the greedy approach of finite-difference discretizations to 2-D elliptic control problems. We devote Section 5 to present the applicability of the (greedy) steady control in time-dependent problems, while in Section 6 we make a detailed analysis of the computational cost of our greedy algorithm. Finally, in Section 7 we make some concluding remarks.

Problem formulation

Let Ω ⊂ R d be an open and bounded Lipschitz domain and consider the parameter dependent parabolic equation with Dirichlet boundary conditions

     yt -div (a(x, ν)∇y) + c(x, ν)y = χωu in Q = Ω × (0, T ), y = 0 on Σ = ∂Ω × (0, T ), y(x, 0) = y0(x) in Ω, (2.1)
where y = y(x, t; ν) is the state, u = u(x, t) is a control function and y0 is given initial datum. Here ω is an open subset of Ω and χω denotes the characteristic function of the set ω where the control is being applied. In (2.1), a(x, ν) is an elliptic, scalar L ∞ (Ω) coefficient depending on some parameter ν ∈ R d and c = c(x, ν) ∈ L ∞ (Ω) is parameter dependent potential. To abridge the notation, we will denote aν = a(x, ν) and cν = c(x, ν)

It is well-known (see, for instance, [START_REF] Evans | Partial differential equations[END_REF]), that for any initial data y0 ∈ L 2 (Ω) and u ∈ L 2 (ω × (0, T )), systems (2.1) admits a unique weak solution y ∈ W (0, T ), where W (0, T ) stands for the Hilbert space

W (0, T ) := y ∈ L 2 (0, T ; H 1 0 (Ω)), yt ∈ L 2 (0, T ; H -1 (Ω)) ,
equipped with the norm

y W (0,T ) = y 2 L 2 (0,T ;H 1 0 (Ω)) + yt 2 L 2 (0,T ;H -1 (Ω)) 1/2 .
Moreover, y satisfies an estimate of the form

y W (0,T ) ≤ C y0 L 2 (Ω) + u L 2 (ω×(0,T )) ,
where C depends on Ω, T , aν and cν ∞ only. Now, consider the following associated optimal control problem min u∈L 2 (0,T ;L 2 (ω))

J T (u) = 1 2 T 0 u(t) 2 L 2 (ω) + β 2 T 0 y(t) -y d 2 L 2 (Ω) , (2.2)
where y is the solution to (2.1), y d ∈ L 2 (Ω) is a desired observation and β > 0 is given. It is classical to prove (see e.g., [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF]) that the minimization problem (2.2) has a unique optimal solution that hereinafter we denote by (u T , y T ). This is due to the fact that the functional J T is strictly convex, continuous and coercive.

Moreover, the optimal control is given by

u T = -χωp T
where p T can be found from the solution (y T , p T ) of the optimality system

         y T t -div (aν ∇y T ) + cν y T = -χωp T , in Q, -p T t -div (aν ∇p T ) + cν p T = β (y T -y d ), in Q, y T = p T = 0, on Σ, y T (x, 0) = y0(x), p T (x, T ) = 0, in Ω. (2.3)
From here, it is clear that the optimal control u T also depends on the parameter ν since the state y equally depends on it. As mentioned before, when studying parameter-dependent problems from the control point of view, this means that one has to solve the minimization problem (2.2) for each new choice of the parameter ν. Although theoretically feasible, the computational effort to compute a control function for every new selection of the parameter is rather expensive and undesirable.

Moreover, the effective computation of optimal controls for problems posed in a long time horizon can be very expensive since it requires iterative methods to solve the coupled optimality system (2.3) combining the forward controlled state equation and the backward adjoint one. Particularly, adjoint equation methods for the solution of control problems posed on long-time intervals may led to the storage of huge quantities of information. Memory saving methods, as for instance discussed in [START_REF] Glowinski | Exact and Approximate Controllability for Distributed Parameter Systems: A Numerical Approach[END_REF] and references within, may help to alleviate this problem.

In this work, we use a combination of (weak) greedy methods and the so-called turnpike property to determine the most relevant values of a parameter-space (to be precised below) and provide the best possible approximation of the set of parameter dependent optimal controls.

Firstly, we will use the turnpike property to reduce the problem of computing the time-dependent optimal controls u T by computing asymptotic simplifications.

To this end, we begin by considering the stationary version of the state equation -div (aν ∇y) + cν y = χωu in Ω, y = 0 on ∂Ω.

(2.4)

Using classical Fredholm theory (see, e.g. [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]), one can prove that for fixed aν , cν ∈ L ∞ (Ω), system (2.4) has a unique weak solution y ∈ H 1 0 (Ω) verifying

α1 y H 1 0 (Ω) ≤ u H -1 (Ω) ≤ α2 y H 1 0 (Ω) . (2.5)
for some constants α1 and α2 (not depending on u), provided the following condition holds

ker(Lν ) = {0}, (2.6)
where Lν is the map defined as

Lν := -div (aν • ) + cI : H 1 0 (Ω) → H -1 (Ω).
Hereinafter, we will briefly discuss the turnpike theory from the point of view where condition (2.6) holds. Observe that when this condition is not fulfilled, the solution to (2.4) is not longer unique and some modifications to the discussion presented below are necessary (see Section 7 and [31,Remark 3.7]), but the results remain valid.

Let us consider the corresponding minimization problem min u∈L 2 (ω)

J(u) = 1 2 u 2 L 2 (ω) + β 2 y -y d 2 L 2 (Ω) . (2.7)
As before, one can prove that, since J is strictly convex, continuous and coercive, the minimization problem (2.7) has a unique optimal solution (ū, ȳ), where the optimal control is characterized by ū = -χω p, and (ȳ, p) solve the optimality system

     -div (aν ∇ȳ) + cν ȳ = -χω p, in Ω, -div (aν ∇p) + cν p = β (ȳ -y d ), in Ω, ȳ = p = 0, on ∂Ω.
(2.8) Also, in this case, it is clear the the optimal control ū depends on ν. Notice that, thanks to condition (2.6), each control ū can be uniquely determined. This is not the case where ker(Lν ) is not trivial, since the solution p to the adjoint equation is defined up to elements of ker(Lν ).

Since each control ū can be uniquely determined, in what follows, we shall write ūν to denote the dependence of the optimal control with respect to the parameter, and similarly for the optimal states yν and pν .

A natural question that arises in this context is to which extent the long horizon optimal controls and states (u T (t), y T (t)) approach the steady ones (ū, ȳ) as T → +∞. According to [START_REF] Porretta | Long time versus steady state optimal control[END_REF], we have the following result: Theorem 2.1 Let us consider the control problem (2.2) and let (u T , y T ) be the optimal control and state. Then, there exists µ > 0 such that

y T (t) -ȳ L 2 (Ω) + u T (t) -ū L 2 (Ω) ≤ K e -µt + e -µ(T -t) , ∀t ∈ [0, T ],
(2.9)

where (ū, ȳ) is the optimal control and state corresponding to (2.8).

This means that we have exponential convergence of the finite horizon control problems to their steady state version as T tends to infinity. In order to prove (2.9), additional assumptions on stabilizability and observability for (2.1) and its adjoint are required. Indeed, although these assumptions are not required when considering solely the optimal control problem (2.2), they are essential for the analysis when the time horizon tends to infinity. Remark 2.2 It is important to mention that for Theorem 2.1 to hold, only the controllability of system (2.1) plays a fundamental role. In other words, positivity or smallness conditions on c(x) guaranteeing the stability of the parabolic problem are not necessary to prove the exponential convergence of the finite time horizon control problem to the steady version.

As pointed out in [START_REF] Porretta | Long time versus steady state optimal control[END_REF], this kind of stationary behavior in the transient time for long horizon control problems is commonly referred in the literature as turnpike property. Such property has been mostly investigated in the finite-dimensional case, as well as in connection to Calculus of Variations. We refer the interested reader to the survey [START_REF] Zaslavski | Turnpike properties in the calculus of variations and optimal control[END_REF] and references therein and to [START_REF] Carlson | Existence of overtaking solutions to infinite dimensional control problems on unbounded time intervals[END_REF][START_REF] Zaslavski | Existence and structure of optimal solutions of infinite-dimensional control problems[END_REF] for some results on control problems in the infinite-dimensional case. We also refer to [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF] for a more recent and systemic discussion on this topic.

Thanks to Theorem 2.1, we now focus our attention on the steady system (2.4) and the approximation of the family of parameterized controls. We begin by assuming the following hypotheses: H1 The parameter ν ranges within a compact set K ⊆ R d .

H2

The coefficient functions aν and cν depend on ν in an analytic manner.

H3 Condition (2.5) holds uniformly for all ν ∈ K.

Then, the main idea is to propose a methodology to determine an optimal selection of a finite number of realizations of the parameter ν so that all controls, for all possible values of ν, are optimally approximated. More precisely, the problem can be formulated as follows:

Problem 2.3 Let us consider the set of controls verifying (2.7) for each possible value ν ∈ K. That is,

U = {ūν : ν ∈ K} .
(2.10)

This control set is compact in L 2 (ω).
Given ε > 0, we seek to determine a family of parameters {ν1, . . . , νn} in K, whose cardinality n depends on ε, so that the corresponding controls, denoted by uν 1 , . . . , uν n are such that for every ν ∈ K, there exists u ν ∈ span{uν 1 , . . . , uν n } such that

u ν -ūν L 2 (ω) ≤ ε.
This problem is motivated by the practical issue of avoiding the computation of a control function uν for each new parameter value ν. By developing suitable greedy algorithms to solve Problem 2.3, we look for the most representative values of ν providing a fast and reliable way to compute approximated optimal controls for any other value of ν.

As noted in [START_REF] Lazar | Greedy controllability of finite dimensional linear systems[END_REF], Hypotheses H1 and H2 makes this goal feasible and, in particular, allows to implement a naive approach where the parameter set K is uniformly sampled in a very fine mesh. Nevertheless, the objective is to minimize the number of spanning controls n and derive the most efficient approximation. As we will see below, H3 enable us to implement a (weak) greedy approach for solving Problem 2.3 in an optimal way. Remark 2. [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF] We can readily distinguish two classes of coefficients verifying Hypothesis H3. These will be of interest in the discussion of Sections 4 and 5.

1. The case where a1 ≤ aν ≤ a2 and cν ≥ 0, for some positive constants a1, a2 and all ν ∈ K. Under these assumptions we can recover (2.5) by Lax-Milgram theorem and classical energy estimates.

3 Greedy optimal control for elliptic problems

Preliminaries on (weak) greedy algorithms

In this section, we present a short introduction about linear approximation theory of parametric problems based on (weak) greedy algorithms. For a more detailed read, we refer, for instance, to [START_REF] Cohen | Approximation of high-dimensional parametric PDEs[END_REF][START_REF] Devore | The Theoretical Foundation of Reduced Basis Methods[END_REF]. We will apply this theory systematically to deal with the optimal control problem of parameterized elliptic equations.

In general, the goal is to approximate a compact set K in a Banach space X by a sequence of finite dimensional subspaces Vn of dimension n. Increasing n improves the accuracy of the approximation.

Although constructing this subspace within a given precision error commonly implies a high computational effort, this process is performed only once as an offline procedure. Then, with a good approximated subspace at hand, one can easily compute online approximations for every vector from K.

Vectors xi, i = 1, ..., n spanning the space Vn are called snapshots of K.

The goal of (weak) greedy algorithms is to construct a family of finite dimensional spaces Vn ≤ X approximating the set K in the best possible way. The algorithm reads as follows.

Algorithm 1: Weak greedy algorithm initialize: fix γ ∈ (0, 1] and a tolerance parameter ε > 0; In the first step, choose x 1 ∈ K such that

x 1 X ≥ γ max x∈K x X ;
At the general step, having found x 1 , . . . , x n denote

V n = span{x 1 , . . . , x n } and σ n (K) := max x∈K dist(x, V n ) ; (3.1) repeat choose x n+1 such that dist(x n+1 , V n ) ≥ γσ n (K); until σ n (K) < ε;
The algorithm produces a finite dimensional space Vn that approximates the set K within precision ε. As mentioned in [START_REF] Devore | The Theoretical Foundation of Reduced Basis Methods[END_REF], it is important to notice that the weak greedy algorithm does not give unique sequences (xn) n≥1 and (σn(K)) n≥1 . However, every chosen sequence decays at the same rate, which under certain assumptions (see Theorem 3.1), is close to the optimal one. Therefore, the algorithm optimizes the number of steps required to satisfy the given tolerance as well as the dimension of the approximated space Vn.

The pure greedy algorithm corresponds to the case γ = 1. As noted in [START_REF] Devore | The Theoretical Foundation of Reduced Basis Methods[END_REF], the relaxation of the pure greedy method (γ = 1) to a weak greedy one (γ ∈ (0, 1]) will not significantly reduce the efficiency of the algorithm, making it, by the contrary, much easier for implementation.

When using the (weak) greedy algorithm one has to choose the elements of the approximation space by exploring the distance defined in (3.1) for all possible values x ∈ K. Two main difficulties arises:

1. the set K in general consists of infinitely many vectors, 2. in practical implementations, the set K is often unknown, e.g., it represents the family of solutions to parameter dependent problems.

The first problem can be avoided by searching over some finite discrete subset of K. Here we use the fact that K, being a compact set, can be covered by a finite number of balls of an arbitrary small radius. To circumvent the second one, instead of considering the distance appearing in (3.1), one uses some surrogate, which is easier to compute.

In order to estimate the efficiency of the weak greedy algorithm we compare its approximation rates σn(K) with the best possible ones.

The best choice of an approximating space Vn is the one producing the smallest approximation error. This smallest error for a compact set K is called the Kolmogorov n-width of K, and is defined as

dn(K) := inf dim Y =n sup x∈K inf y∈Y x -y X .
It measures how well K can be approximated by a subspace in X of a fixed dimension n.

In the sequel we want to compare σn(K) with the Kolmogorov width dn(K), which represents the best possible approximation of K by a n dimensional subspace of the referent Banach space X. A precise estimate in that direction was provided by [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF] in the Hilbert space setting, and later improved and extended to the case of a general Banach space X in [START_REF] Devore | Greedy Algorithms for Reduced Bases in Banach Spaces[END_REF].

Theorem 3.1 (Corollary 3.3 of [START_REF] Devore | Greedy Algorithms for Reduced Bases in Banach Spaces[END_REF]) For the weak greedy algorithm with constant γ in a Hilbert space X we have the following: If the compact set K is such that, for some α > 0 and C0 > 0

dn(K) ≤ C0n -α , n ∈ N, then σn(K) ≤ C1n -α , n ∈ N,
where

C1 := γ -2 2 5α+1 C0.
This theorem implies that the weak greedy algorithms preserve the polynomial decay rates of the approximation errors. A similar estimate also holds for exponential decays (cf. [START_REF] Devore | Greedy Algorithms for Reduced Bases in Banach Spaces[END_REF]).

Definition of the residual

As mentioned before, one of the main goals of this paper is to apply the greedy approach described above to the family of parameter dependent steady state optimal control problems (P) min

u∈L 2 (ω) {Jν (u)} , (3.2)
where J is the cost functional given by

Jν (u) = 1 2 u 2 L 2 (ω) + 1 2 yν (u) -y d 2 L 2 (Ω) , while yν (u) is the solution to (2.4).
As stated in Problem 2.3, the aim is to choose the most representative set of parameter values νi, whose associated controls ui will provide a good approximation of the control manifold (2.10). Essential for an effective application of a greedy algorithm is the construction of a residual by which one can estimate the distance between two (possible unknowns) controls by some easily computable quantity.

To construct an appropriate residual, we begin by computing the optimality condition to the minimization problem (3.2). More precisely, for given ν ∈ K, the optimality condition read as follows

ūν + S * ν (Sν ūν -y d ) = 0 (3.3)
where Sν : L 2 (ω) → L 2 (Ω) is an operator that assigns to the control u the solution yν (u) of the problem (2.4), while S * ν : L 2 (Ω) → L 2 (ω) is its adjoint operator defined by

S * ν y = χωp
where p is the solution to the adjoint problem with the right-hand side equal to βy. As mentioned in the previous section, this leads to the following optimality system ūν = -χω pν ,

-div (aν ∇ȳν ) + cν ȳν = χω ūν , -div (aν ∇pν ) + cν pν = β (ȳν -y d ), (3.4)
where ȳν stands for the optimal state, while pν for the optimal dual variable. As before, H3 plays a key role to define uniquely thesolution to the optimality system (3.4).

To this effect one can consider ∇Jν as a candidate for the residual. Indeed, one can readily verify that ūν

-u L 2 (ω) ≈ ∇Jν (u) L 2 (ω) .
This allows to estimate the distance between the given control u from the (unknown) optimal one ū by checking the minimization performance of the first , i.e, by calculating ∇Jν (u). However, as ∇Jν is expressed in terms of S and S * , that calculation requires solving both the primal and dual problems, which, in general, is expensive (cf. [START_REF] Tröltzsch | POD a-posteriori error estimates for linear-quadratic optimal control problems[END_REF]).

In order to reduce the cost of computing the surrogate, we can rewrite the cost functional in terms of the dual variable p and state variable y. Namely, the gradient of the cost functional is expressed as

∇Jν (u) = u + χωp(u),
where the dual variable p is related to the control by the relation 1ωp(u) = S * ν (Sν u -y d ), i.e.

-div (aν ∇y(u))

+ cν y(u) = χωu -div (aν ∇p(u)) + cν p(u) = β y(u) -y d .
In what follows, to abridge the notation, for any z ∈ H 1 0 (Ω) and any parameter ν ∈ K, we will use the notation

Lν z := -div (aν ∇z) + cν z

In this way, we introduce the residual operator as

Rν (p, y) := Lν y + χωp Lν p -β y -y d . (3.5)
Of course, for optimal variables (pν , ȳν ) the value of the residual equals zero, and the residual can be equivalently written as

Rν (p, y) := Lν (y -ȳν ) + χω(p -pν ) Lν (p -pν ) -β (y -ȳν ) . (3.6)
By means of the introduced operator we shall try to perform greedy approximation of the manifold

P × Ȳ = {(pν , ȳν ) : ν ∈ K} ⊂ H 1 0 (Ω) × H 1 0 (Ω).
We will show that for an arbitrary (p, y) ∈ H 1 0 (Ω) × H 1 0 (Ω) the residual Rν (p, y) provides a good measure of the distance of (p, y) from the optimal pair (pν , ȳν ). To this end, from definition (3.6), we easily obtain that

H -1 (Ω) Rν (p, y), p -pν -(y -ȳν ) H 1 (Ω) = χω(p -pν ) 2 L 2 (Ω) + β y -ȳν 2 L 2 (Ω) ≥ 2 β χω(p -pν ), y -ȳν L 2 (Ω) . (3.7)
Similarly, one gets

H -1 (Ω)
Rν (p, y), y -ȳν 0

H 1 (Ω) ≥ α1 y -ȳν 2 H 1 (Ω) -χω(p -pν ), y -ȳν L 2 (Ω) (3.8)
where α1 is the uniform bound from (2.5).

Combining the last inequality with (3.7) we obtain

c1 y -ȳν 2 H 1 (Ω) ≤ Rν (p, y) H -1 (Ω) ( p -pν H 1 (Ω) + y -ȳν H 1 (Ω) ),
where c1 depends on α1 and β only.

In order to obtain a similar bound for the dual variable, we consider the following product

H -1 (Ω) Rν (p, y), y -ȳν p -pν H 1 (Ω) = H -1 (Ω)
Lν (y -ȳν ), (y -ȳν )

H 1 (Ω) + χω(p -pν ), y -ȳν L 2 (Ω) + H -1 (Ω)
Lν (p -pν ), (p -pν )

H 1 (Ω) -β y -ȳν , p -pν L 2 (Ω) ≥ α1 y -ȳ 2 H 1 (Ω) -χω(p -pν ), y -ȳν L 2 (Ω) + α1 p -pν 2 H 1 (Ω) -β| y -ȳν , p -pν L 2 (Ω) |.
In order to estimate the last term, note that for an arbitrary ε > 0 we have

| y -ȳν , p -pν L 2 (Ω) | ≤ 1 2ε y -ȳν 2 L 2 (Ω) + ε 2 p -pν 2 L 2 (Ω)
.

By taking ε = α 1 β from (3.8), we obtain α1 2 ( y -ȳν 2 H 1 (Ω) + p -pν 2 H 1 (Ω) ) ≤ Rν (p, y) H -1 (Ω) ( p -pν H 1 (Ω) + y -ȳν H 1 (Ω) ) + χω(p -pν ), y -ȳν L 2 (Ω) + 1 2α1 y -ȳν 2 H 1 (Ω) .
Finally, combining the last estimate with (3.7) and (3.8) we obtain

c1( y -ȳν H 1 (Ω) + p -pν H 1 (Ω) ) ≤ Rν (p, y) H -1 (Ω) , (3.9)
for some constant c1 > 0 only depending on α1 and β. This shows, in particular, that the introduced residual provides a good estimate of the distance of (p, y) from the optimal pair (pν , ȳν ).

In order to obtain a reverse type inequality, we explore the regularity estimates (2.5). Combining it with the inclusion of

H 1 0 (Ω) into H -1 (Ω), we obtain Rν (p, y) H -1 (Ω) ≤ c2( y -ȳν H 1 (Ω) + p -pν H 1 (Ω) ), (3.10)
with the constant c2 > 0 depending on α2 and β only.

The last two inequalities enable us to apply a weak greedy procedure for the approximation of the manifold P × Ȳ (see Algorithm 2 below).

The greedy algorithm

The residual Rν (p, y) introduced in the previous subsection enable us to construct a weak greedy algorithm for an effective construction of an approximating linear space of the control manifold U.

The precise description of the offline part of the algorithm is given below in Algorithm 2. This algorithm results in the approximating space (Pn, Yn) = span{(p1, ȳ1) . . . , (pn, ȳn)}, where n is a number of chosen snapshots (specially (Pn, Yn) = {0} for n = 0). In particular, the obtained bounds on the residual (3.9)-(3.10) ensure that Algorithm 2 satisfies the requirements of the weak greedy procedure. More precisely, the following result holds. Theorem 3.2 For a given ε > 0 take the discretisation constant δ such that δ ≤ ε/(2CLc1).

(3.11)

Then the Algorithm 2 provides a weak greedy approximation of the manifold (P, Y) with the constant

γ = c1 2c2 , (3.12)
and the approximation error less than ε/c1. Here CL is the Lipschitz constant of the mapping ν → (pν , ȳν ), while c1 and c2 are constants appearing in estimates (3.9) and (3.10), respectively.

Remark 3.3 Note that the infimum appearing in (3.17) can be expressed as a distance from a suitable space determined by the residual Rν . More precisely, denote by Gν the linear part of the residual Rν , i.e.

Gν (p, y) := Lν y + χωp Lν p -βy , (3.13)
implying

Rν (p, y) = Gν (p, y) + 0 βy d .
Denoting by Gν (Pj, Y j ) the space spanned by Gν (p i , y i ), i = 1, . . . , j, we obtain

|| inf (p,y)∈(P j ,Y j ) Rν (p, y)|| H -1 (Ω) = dist 0 -βy d , Gν (Pj, Y j ) . (3.14)
The last relation enables practical computation of the infimum appearing in (3.17) which can be performed by projecting vector 0 -βy d to the space Gν (Pj, Y j ).

Algorithm 2: Greedy control algorithm -offline part

Initialize: Fix the approximation error ε > 0. STEP 1: (discretisation) Choose a finite subset K such that

(∀ ν ∈ K) dist(ν, K) < δ,
where δ > 0 is a constant determined by (3.11). STEP 2: (Choosing ν 1 ) Check the inequality max

ν∈ K βy d ν H -1 (Ω) < ε 2 . (3.15)
If it is satisfied, stop the algorithm. Otherwise, choose the first parameter value as

ν 1 ∈ argmax ν∈ K { y d ν H -1 (Ω) }. (3.16)
and find corresponding optimal primal and dual states ȳ1 , p1 ; STEP 3: (Choosing ν j+1 ) Having chosen ν 1 , . . . , ν j calculate R ν (p j , 0) and R ν (0, ȳj ) for each ν ∈ K. 10 Check the approximation criteria

max ν∈ K || inf (p,y)∈(Pj ,Y j ) R ν (p, y)|| H -1 (Ω) < ε 2 . (3.17)
11 If the inequality is satisfied, stop the algorithm. Otherwise, determine the next parameter value as

ν j+1 ∈ argmax ν∈ K || inf (p,y)∈(Pj ,Y j ) R ν (p, y)|| H -1 (Ω) , (3.18)
12 Find the corresponding optimal primal and dual states ȳj+1 , pj+1 and repeat Step 3.

Remark 3.4 Observe that the above algorithm definitely stops after finite, n ∈ N0 number of iterations due to the compactness assumption on the parameter set K. In the sequel we exclude the case n = 0, occurring when inequality (3.15) holds, which results in a null approximating space for which the statement trivially holds

Proof of Theorem 3.2. We will adapt the weak greedy procedure presented in Algorithm 1. For this particular problem, we have to show that i) the selected optimal states (pi, ȳi) associated to the parameter values determined by (3.16) and (3.18) satisfy

(pν 1 , ȳν 1 ) H 1 (Ω) ≥ γ max ν∈K (pν , ȳν ) H 1 (Ω) dist (pν j+1 , ȳν j+1 ), ( Pi, Ȳi) ≥ γ max ν∈K dist (pν , ȳν ), ( Pi, Ȳi) , j = 1, . . . , n -1
with the constant γ given by (3.12), ii) the approximation error, defined as

σ (p, ȳ)(K) := max ν dist (ȳν , pν ), ( Pi, Ȳi)
obtained at the end of the algorithm is less than ε/c1.

In order to prove i), we begin by noting that from estimate (3.9) and the definition of the residual Rν , the following inequality holds

(pν , ȳν ) H 1 (Ω) = dist (pν , ȳν ), (0, 0) ≤ 1 c1 Rν (0, 0) H -1 (Ω) = 1 c1 0 βy d ν H -1 (Ω)
, for all ν ∈ K. From the choice of the first snapshot (3.16) for every parameter ν from the discretised set K the above inequality readily implies

(pν , ȳν ) H 1 (Ω) ≤ 1 c1 0 βy d ν 1 , H -1 (Ω) , = 1 c1 Rν 1 (0, 0) H -1 (Ω) ≤ c2 c1 (pν 1 , ȳν 1 ) H 1 (Ω) , (3.19)
where in the last line we have used the estimate (3.10).

In order to obtain a similar inequality for every ν ∈ K, by taking ν ∈ K such that |ν -ν| < δ, by means of (3.11) and (3.19) we have

(pν , ȳν ) H 1 (Ω) ≤ (pν -pν , ȳν -ȳν ) H 1 (Ω) + (pν , ȳν ) H 1 (Ω) ≤ ε 2c1 + c2 c1 (pν 1 , ȳν 1 ) H 1 (Ω) .
Having excluded the case (3.15), we have

ε 2 ≤ Rν 1 (0, 0) H -1 (Ω) ≤ c2 (pν 1 , ȳν 1 ) H 1 (Ω) ,
implying the first inequality in i).

The same arguments can be employed for the general j-th iteration. Recall that the stopping criterion (3.17) is not fulfilled for iterations j < n. Suppose that we have selected ν1, . . . , νi different parameters. From inequalities (3.9), (3.10) and the definition of the next snapshot (see eq. (3.18)), for each ν ∈ K we readily obtain dist (pν , ȳν ), ( Pi, Ȳi)

≤ || 1 c1 inf (p,y)∈(P i ,Y i ) Rν (p, y)|| H -1 (Ω) ≤ || 1 c1 inf (p,y)∈(P i ,Y i ) Rν j+1 (p, y)|| H -1 (Ω) ≤ c2 c1
dist (pν j+1 , ȳν j+1 ), ( Pi, Ȳi) .

As before, the above inequality can be generalized to all ν ∈ K by employing Lipschitz continuity of the mapping ν → (p, ȳ), particularly implying the second inequality in i).

To prove ii), it is straightforward to see from equation (3.9) that after n iterations and for any ν ∈ K by selecting ν ∈ K such that |ν -ν| < δ, the following holds dist (pν , ȳν ), ( Pn, Ȳn) ≤ (pν -pν , ȳν -ȳν ) H 1 (Ω) + dist (pν , ȳν ), ( Pn, Ȳn)

≤ ε 2c1 + 1 c1 Rν ( Pn, Ȳn) H -1 (Ω) < ε c1 ,
since the stopping condition (3.17) has been reached. This ends the proof of Theorem 3.2. It remains to check the effect of the approximation of the manifold P × Ȳ on the control manifold Ū and, in the end, the approximation of the cost functional Jν .

In the next step we want to explore the approximation space constructed in the offline part of the algorithm with the aim of an effective construction of an approximating states and controls for an arbitrary given parameter value. This is the step of the greedy method usually referred to as the online part.

To this effect we propose the following algorithm for the online part of the greedy procedure.

Using the above proposed definition of the approximative optimal control (3.20) and the approximation performance of the offline part of the algorithm we immediately obtain

u ν -ūν L 2 (ω) ≤ p ν -pν H 1 (Ω) ≤ ε c1 , (3.22)
Algorithm 3: Greedy control algorithm -online part Initialize: A parameter value ν ∈ K is given. STEP 1:

Project the vector 0 -βy d to G ν (P n , Y n ), the space spanned by G ν (p i , y i ), i = 1, . . . , n.

Here {(p i , y i ), i = 1, . . . , n} are snapshots selected in the offline part of the algorithm, while G ν is the linear part of the residual given by (3.13).

Denote the projection by P n 0 -βy d . STEP 2: Solve the system

P n 0 -βy d = α i G ν (p i , y i ).

STEP 3:

Define the approximating dual state as

p * ν = α i p i .
STEP 4: 10 Define the approximating optimal control as

u * ν = -χ ω p * ν . (3.20)
STEP 5: 11 Define the approximating optimal state y ν as the solution to

L ν y ν = χ ω u ν . (3.21)
which provides the required estimate proposed by Problem 1 (up to a scaling factor c1).

Similarly, from the regularity bound (2.5) it follows

y ν -ȳν H 1 (Ω) ≤ 1 α1 u ν -ūν L 2 (ω) ≤ ε α1c1 , (3.23)
which provides the estimate on the approximative optimal state obtained by the greedy procedure. Finally, we obtain the estimate for the error in the cost functional. Note that

|J(u ν ) -J(ūν )| = 1 2 ūν 2 L 2 (Ω) -u ν 2 L 2 (Ω) + ȳν -y d 2 L 2 (Ω) -y ν -y d 2 L 2 (Ω) . (3.24)
The first two terms on the right hand side can be bounded as follows:

Ω |ūν | 2 -|u ν | 2 dx = Ω (|ūν | -|u ν |)(|ūν | + |u ν |)dx ≤ Ω |ūν -u ν |(|ūν -u ν | + 2|ūν |)dx ≤ ūν -u ν L 2 (Ω) ūν -u ν L 2 (Ω) + 2 ūν L 2 (Ω) .
(3.25)

By using optimality condition (3.3) one easily obtains the bound

ūν L 2 (Ω) ≤ S * y d L 2 (Ω) , (3.26)
which together with (3.22) enables the first difference in (3.25) to be estimated in terms of the approximation error ε and the given problem data.

Similarly, for the last two terms in (3.24) one gets ȳν -

y d 2 L 2 (Ω) -y ν -y d 2 L 2 (Ω) ≤ ȳν -y ν L 2 (Ω) ȳν -y ν L 2 (Ω) + 2 ȳν -y d L 2 (Ω) ≤ ȳν -y ν L 2 (Ω) ȳν -y ν L 2 (Ω) + 2 α1 S * y d L 2 (Ω) + 2 y d L 2 (Ω)
where the last inequality follows from the regularity assumptions (2.5), the definition of the optimal state (3.4)2 and the estimate (3.26). Thus, combining estimates on the error on approximate control (3.22) and the approximate state (3.23), one obtains bound on the approximation of the minimal cost functional value in terms of the approximation error ε and the desired trajectory y d .

We can summarise the obtained results in the following theorem Theorem 3.5 The proposed greedy algorithm provides the following estimates on the approximate control, approximate optimal state and approximate minimal value of the cost functional

• u ν -ūν L 2 (Ω) ≤ ε, • y ν -ȳν H 1 (Ω) ≤ ε α1 , • |J(u ν ) -J(ūν )| ≤ ε c 1 ε c 1 1 + 1 α 2 1 + 2 1 + 1 α 2 1 S * y d L 2 (Ω) + 2 α 1 y d L 2 (Ω) ,
where u ν and y ν are given by (3.20) and (3.21), respectively, while α1 is the elliptic bound appearing in (2.5).

Optimal approximation rates

One of main advantages of the greedy procedure is that it provides optimal approximation rates in terms of Kolmogorov widths. However, although the Kolmogorov widths of a set of admissible parameters (set K in our problem) is usually easy to estimate, this is not the case for the corresponding set of controls ( Ū in our setting). Fortunately, a result in that direction has been provided recently for holomorphic mappings ( [START_REF] Cohen | Kolmogorov widths under holomorphic mappings[END_REF]) under the assumption of a polynomial decay of Kolmogorov widths. Theorem 3.6 For a pair of complex Banach spaces X and V assume that u is a holomorphic map from an open set O ⊂ X into V with uniform bound. If K ⊂ O is a compact subset of X then for any α > 1 and C0 > 0

dn(K) ≤ C0n -α =⇒ dn(u(K)) ≤ c1n -β , n ∈ N,
for any β < α -1 and the constant C1 depending on C0, α and the mapping u.

Remark 3.7 The proof of the theorem provides an explicit estimate of the constant C1 in dependence on C0, α and the mapping u. However, due to its rather complicated form we do not expose it here.

Going back to our problem, Theorem 3.6 can be applied if we proved that the mapping ν → ūν is analytic, which in terms of the optimality equations (3.4) provides that the mapping K to P × Ȳ is analytic as well.

To this effect we explore analytic version of the implicit function theorem and apply it to the mapping

DuJν (u) = u + S * ν (Sν u -y d ) .
Of course, for an optimal control it holds DuJν (ūν ) = 0, and in order to employ the implicit function theorem we have to check if the operator Du(DuJν )(ūν ) is an isomorphism. As Du(DuJν )(ūν ) = I + S * ν Sν and using that S * ν Sν is positive semi-definite it follows that Du(DuJν )(ūν ) is coercive, bounded operator, thus allowing for the inverse.

Furthermore, the operator Sν inhertis regulrity properties of the coefficients aν , cν in its dependence with respect on ν, and this on all levels. Specially, assuming the mappings ν → aν and ν → cν are analytic the same holds for the map ν → Sν , and similarly for S * ν . Thus the implicit function theorem implies that the unique mapping ν → ūν is analytic.

Combining this result with Theorem 3.6 we obtain the following one. The last result can be extend to the case of infinite number of parameters as well, provided the Kolmogorov widths of the set K decay polynomially.

Numerical experiments

We devote this section to present two numerical examples where the greedy approach is applied. To illustrate the procedure, we consider the two dimensional domain Ω = (0, 1) 2 . In order to approximate the solution of the elliptic problem (2.4) we use the elementary finite difference (FD) method. We shall use uniform meshes, i.e., meshes with constant discretization steps in each direction. That is, for given Nx 1 , Nx 2 ∈ N * , we set

x1,i = ihx 1 , i ∈ [0, Nx 1 + 1], hx 1 = 1 Nx 1 + 1 , x2,j = jhx 2 j ∈ [0, Nx 2 + 1], hx 2 = 1 Nx 2 + 1 .
Thus, by abuse of notation, we denote by yi,j an approximate value of the solution to (2.4) at the grid points (x1,i, x2,j).

We will approximate the elliptic operator Aν = -div(aν ∇ •) with homogeneous Dirichlet boundary conditions, by using the standard 5-point discretization given by

(A ν h y)i,j = -a ν i+ 1 2 ,j yi+1,j -yi,j h 2 x 1 + a ν i-1 2 ,j yi,j -yi-1,j h 2 x 1 -a ν i,j+ 1 2 yi,j+1 -yi,j h 2 x 2 + a ν i,j-1 2 yi,j -yi,j-1 h 2 x 2 (4.1)
where we set a ν i± 1 2 ,j = aν (xi ± hx 1 /2, yj) and so on. Observe that the boundary conditions are taken into account in those formulas by imposing that y0,j = yN x 1 +1,j = 0, 1 ≤ j ≤ Nx 2 and yi,0 = yi,N y +1 = 0, 1 ≤ i ≤ Nx 1 .

The linear system of algebraic equations derived from the FD discretization (4.1) in N ×N grid points for an elliptic PDE has O(N 2 ) equations, so the coefficient matrix has O(N 2 × N 2 ) entries. Even for a small number N = 200, the resulting matrix cannot be stored and handled in most computers. However, the matrix for a self adjoint elliptic operator as in (2.4) obtained from (4.1) is sparse with number of non-zero entries ∼ 5N 2 . Thus, sparse matrix techniques will be used in what follows. We refer to Section 6 for a more detailed description of the computational cost of implementing the greedy approach.

Greedy test # 1

Let us consider a coefficient aν of the form

aν (x) = 1 + ν(x 2 1 + x 2 2 )
where we assume that

ν ∈ [1, 10] = K. (4.2)
In this way, we can readily verify that aν satisfies H1 and H2.

On the other hand, we take the coefficient c as c(x) = sin(πx1) sin(πx2).

which clearly fulfills the condition c ≥ 0. For the discretization of the elliptic problem, we choose Nx 1 = Nx 2 = 400 to have a uniform, equally sized mesh in both directions. Finally, for the optimal control problem, we set the parameter β = 10 4 while the desired target is the x2-independent function y d (x) = sin(2πx1). We take ω as the region (0.2, 0.5) × (0, 5, 0.9) ∪ (0.5, 0.9) × (0.1, 0.9) in the x1x2-plane (see Figure 2b). Since the functional Jν is quadratic and coercive, a standard conjugate gradient (CG) algorithm is a quite natural and simple choice to solve the minimization problem (2.7). In fact, for a given ν ∈ K one can directly solve the minimization problem (3.2). The average time for computing the corresponding control up to a given tolerance of 10 -8 using the CG is around six seconds.

As mentioned in [START_REF] Lazar | Greedy controllability of finite dimensional linear systems[END_REF], Hypotheses H1 and H2 allow us to implement a naive approach for approximating the parameterized control set Ū(K). This approach consists in discretizing the parameter set in a very fine mesh, that we denote it by K, and then computing the corresponding control for each parameter in this finite-dimensional set. If the number of parameters in K is rather small, or the desired precision ε is large enough, one can directly use the naive approach.

For instance, in this particular test, let us take K as the uniform discretization of the interval (4.2) in k = 100 values. Then, the naive approach amounts to solve 100 different times the minimization problem (3.2). The whole process takes around 600 seconds and requires to store all the computed controls. Therefore, one expects that the greedy approach will help us to reduce the computing effort.

We apply the greedy procedure described in Algorithm 2 over the set K and choosing a precision of ε = 0.005. The algorithm stops after seven iterations selecting 7 (out of 100) parameter values. The corresponding optimal states (pν i , ȳν i ) associated to each selected parameter are stored in memory.

The time needed to finish the offline part of the greedy algorithm takes 476 seconds. Compared to the time that takes the naive approach, this might seem as a marginal gain. As shown in Section 6, the offline part is costly. However, this test being an academic example, does not exhibit all of the potential of the greedy theory. In fact, in real life applications where simulation time may take hours, it is indeed much cheaper to compute a residual as (3.5) than solving a large number of control problems associated to that many parameters.

The way the parameters are chosen for this test is illustrated in Figure 1a. One can see that the selected values are taken from different parts of the parameter set in a zig-zag manner, leading to the best possible approximation. Indeed, in Figure 1b, we plot the approximation rate of the greedy algorithm corresponding to σn( K) = max ν∈ K dist (pν , ȳν ), ( Pn, Ȳn) . Such plot suggests an exponential decay of the approximation rate σn, which is in accordance with Corollary 3.8 that provides a polynomial decay at any rate.

Following the discussion in Section 3, once the offline part is completed, we can construct (approximate) optimal controls by choosing a suitable combination of the optimal states pν i . The methodology to construct such approximations is described in Algorithm 3. In Figure 2, we plot the approximation of the real control with the greedy algorithm for ν = √ 2. The approximated control is nearly identical to the obtained, for instance, by minimizing directly functional (2.7) for the associated value ν = √ 2. In fact, one can obtain for this particular experiment that

ū√ 2 -u √ 2 L 2 (Ω) ≈ 4.61 × 10 -5 .
In spite of obtaining almost the same solution, the convergence of conjugate gradient method takes 5.4 s compared to 0.488 s in the online greedy part. This fact also shows the computational efficiency of the proposed algorithm.

In Figure 3a, we plot the solution to (2.4) using the approximated control u √ 2 . As for the control, we can compute the difference between the real optimal state against the approximate optimal state ȳ √ 2 . In this particular test, we obtain the following estimate ȳ√ 2 -y √ 2 L 2 (Ω) ≈ 7.14 × 10 -7 .

Note, however, that the approximation of the state to the desired target (see figure 3b) is quite different from what one may expect. This fact is not related with the greedy procedure. Actually, it is closely related to the geometry of the domain ω and the penalization parameter β. One can obtain a better approximation of the desired target y d by taking a larger value of β, but at the price of increasing the L 2 -norm of the control, or by taking a bigger control set ω, leading to an increasing number of degrees of freedom.

Greedy test # 2

Here, we present a series of experiments for the case that fit the setting of example 2 in Remark 2.4. As mentioned before, this will be of interest when discussing the applicability of the method for the turnpike problem.

To this end, let us consider the case where a(x) ≡ 1 for every parameter ν. It is well-known that the eigenvalues of the Laplacian on a square are given by λnm = π 2 (n 2 + m 2 ), n, m ∈ N. Hence, a straightforward example for testing our greedy methodology is the case where c(x, ν) = ν with ν ∈ [-5π 2 + , -2π 2 -] and > 0. We see that in this case, H3 holds as long as we are sufficiently far from the eigenvalues. We will use as a control region the shape depicted in Figure 4. For this particular test, the average time for computing the optimal control for different values of ν is around 10 seconds. Implementing the naive approach for the refined mesh K implies that at least 1000 seconds are needed to finish this process. We will see that by means of the greedy approach we can reduce the computational effort.

Using our computational tool, we choose again the approximation tolerance ε = 0.005. The greedy algorithm stops after five iterations. We present in figure 5a the selected parameters in the order they were chosen by the program. As in the previous test, the selection of the parameters leads to an exponential decay of the approximation rates σn (see figure 5b). The elapsed time for completing the offline process is 389 seconds. This shows a clear improvement of almost 60% less time with respect to the time consumed in naive approach.

In Figure 6, we plot the approximation of the optimal control obtained by the greedy method for a chosen value of ν = -5π 2 /2, while in Figure 7 we show corresponding controlled state and a comparison with the target y d . For this particular test, the elapsed time to compute the online control is 0.406 seconds while computing the optimal control directly by means of the CG method takes 15.9 seconds. The approximation error with respect to the real control is

ū-5π 2 /2 -u -5π 2 /2 L 2 (Ω) ≈ 5.81 × 10 -4 ,
while the approximation error for the state using the real and approximated control is ȳ-5π 2 /2 -y -5π 2 /2 L 2 (Ω) ≈ 6.69 × 10 -5 .

In Figure 7, we observe that the approximation to the desired target y d is (to a certain extent) better than the one we obtain in greedy test #1. This is because the chosen y d does not change sign in the domain Ω (cf. Figure 3b). This also translates into a less effort by the control (see Figure 6), even if the control set ω is smaller than in the previous case (cf. Figure 2).

Connection with turnpike problems

We devote this section to discuss the utility and limitations of the greedy approach for elliptic equations in the control of evolution problems.

For the moment, let us turn our attention to the time dependent optimization problem min

u∈L 2 (0,T ;L 2 (ω)) J T (u) = 1 2 T 0 u(t) 2 L 2 (ω) + β 2 T 0 y(t) -y d 2 L 2 (Ω) , (5.1)
where y is solution to the parabolic problem

    
yt -div (aν ∇y) + cν y = 1ωu, in Ω × (0, T ), y = 0, on ∂Ω × (0, T ), y(x, 0) = y0(x), in Ω.

(

As stated in Section 2, the turnpike property through Theorem 2.1 ensures that the optimal state and control (u T (t), y T (t)) solution to (5.1) both simplify exponentially in time large to their steady counterpart (ū, ȳ). From this fact, a first intuitive experiment is to put the control ū in (5.2) and test its capability to control the system.

Here, we present a series of experiments related to the (approximated) optimal steady controls in Sections 4.1 and 4.2. We use them to control the corresponding evolution equation and test their efficiency. We will differentiate two main cases to be studied.

The case c(x) ≥ 0

In addition to the parameters already set in the greedy test #1 (see Section 4.1), let us take

• T = 3, • y0(x) ≡ 0.
In this stage, we are going to use the optimal control ū √ 2 as a time-independent control for system (5.2). More precisely, we take u(x, t) = ū √ 2 (x) for t ∈ (0, T ).

(

We recall that this control is shown in Figure 2. Now, by plugging such control into equation (5.2), we obtain the controlled solution displayed in Figure 8. We see that this control steers y away from the initial condition at time t = 0 and reaches in a short amount of time a region near the optimal steady state ȳ √ 2 (cf. Figure 8c). For this experiment, the efficiency of the time-independent control is closely related to the fact that c(x) ≥ 0. In particular, the conditions on the coefficients aν and c allow to prove that the uncontrolled system (5.2) is exponentially stable regardless of the initial datum. Figure 8: Evolution in time for a stable system controlled with the steady control approximated by greedy procedure.

In Figure 9, we illustrate the efficiency of the steady control by plotting different curves representing the time evolution of the L 2 -norm of y solution to (5.2) with different initial datums and taking (5.3) as a control. To compare, we have computed the time-dependent solution y T (t) associated to the optimal control u T (t), obtained by minimizing (5.1), for the given parameter ν = √ 2 and y0(x) = 0. Note that we have split the time horizon in two different intervals for the sake of clarity.

The solid line corresponds to the trajectory obtained by using the turnpike control u T (t) and exhibits the prototypical behavior expected from the turnpike theory. For most of the time, except at the beginning and end, the solution remains exponentially close to corresponding steady state ȳ.

On the other hand, the discontinuous lines represent the trajectories of the system controlled with the steady control (5.3) and differential initial datums. We see that this control drives the solution y to the same region (up to an approximation error inherited from the greedy procedure) and maintains the solution there during the whole time interval independently of the initial data. For more that 95% percent of the time, both solutions are almost indistinguishable (up to the end of the interval, where they are superimposed).

To some extent, the turnpike property, together with the greedy procedure for elliptic optimal control problems, gives a partial answer on how to propose a greedy methodology that is independent of the initial datum for the system under consideration (cf. [START_REF] Lazar | Greedy controllability of finite dimensional linear systems[END_REF]).

The case c(x) < 0

In spite of the above conclusions, it is important to mention that the same behavior cannot be expected for all coefficients c(x) < 0. In fact, if c(x) < -λ1 where λ1 if the first eigenvalue of the operator -div (aν ∇ •), then system (5.2) is shown to be unstable.

As an example, let us take the parameters and results in the greedy test #2 shown in Section 4.2 together with y0(x) = 0 and T = 3. As before, we can put the (approximated) steady control ū -5π 2 /2 (x) in system (5.2) and test its performance. Recall that we have chosen cν ≡ -5π 2 /2, which clearly satisfies cν < -λ1. Thus, in this case, the uncontrolled evolution equation is shown to be unstable.

We show in Figure 10 the solution to the evolution problem using the steady control ū -5π 2 /2 for t ∈ [0, 0, 3]. We see that in this case, the steady control lacks to stabilize the system around the steady state shown in figure 7a and, moreover, continues growing exponentially during the remainder of the time interval.

y(t) -ȳ L 2 (Ω)
2.9 3 time (sec) Figure 9: Difference between the steady state ȳ and the solution to system (5.2) with different controls: turnpike control ( ), steady control with y 0 (x) = 0 ( ) and ỹ0 (x) = χ (0.4,0.9) 2 (x) ( ). In Figure 11, we present further experiments showing that for some initial data the controlled solution y grows exponentially in time. In fact, only for y0(x) = ȳ -5π 2 /2 and a sufficiently small neighborhood of this initial datum, the steady control is effective to control the underlying system.

Note however that this behavior is compatible with the fulfillment of the turnpike property which relies on the controllability of the system and not on its the stability. In Figure 12, we plot the time evolution of the optimal controlled state and control (obtained by minimizing (5.1)) and, as expected, we can see the asymptotic simplification towards the state and control for most of the time horizon. However, we can also see that the control makes a large effort at the beginning of the temporal interval to move the system to this steady state. Thus we conclude that unlike the case c(x) ≥ 0, the steady turnpike controls, and in particular the approximated controls derived from the greedy approach, are in general not enough to control an unstable evolution system (5.2). (Ω)

Figure 11: Time evolution comparison for an unstable system controlled with steady control with different initial data: y 0 (x) = ȳ -5π 2 /2 (x) ( ), y 0 (x) = 0 ( ) and ỹ0 (x) = χ (0.4,0.9) 2 (x) ( ). Quasi-optimal trajectories Real turnpike

Figure 13: Quasi-optimal trajectories (for different T 0 ) passing through the steady state ȳ to emulate the turnpike property.

Remark 5.1 The piecewise control strategy bellow allows to build quasi-turnpike control strategies and trajectories, for which the initial arc is particularly relevant when the system under consideration is unstable. More precisely, given the steady solution ȳ to (2.7) and using the property of controllability to trajectories of (5.2) (see, for instance, [START_REF] Fernández-Cara | Null and approximate controllability for weakly blowing up semilinear heat equations[END_REF],[17, Remark 6.1]), one can construct, for any given T0 > 0, a control v0 steering the system to y(T0) = ȳ. Then, at t = T0, one can switch to the steady control to remain there. The same idea can be used to exit the steady state. In this way, we can approximate the real turnpike behavior with a quasi-optimal trajectory (see Figure 13). Note however that in this strategy T0 is a new design parameter, since the null controllability for parabolic problems holds for any time. The cost of controllability depends on T0 in a non-trivial and very sensitive way, so that it diverges exponentially when T0 → 0 (see, for instance, [START_REF] Fursikov | Controllability of evolution equations[END_REF]). Thus the optimal way to choose T0 such that this strategy emulates as close as possible the turnpike property is an open and interesting problem.

Computational cost

In this section, we make a more detailed analysis on the computational cost of the greedy algorithm versus the naive approach. We rely on classic numerical linear algebra results (see, for instance, [START_REF] Boyd | Convex optimization[END_REF][START_REF] Quarteroni | Numerical mathematics[END_REF]). The analysis is restricted to 2-D problems for which numerical results are presented.

During numerical implementation of Algorithm 2, two main basic problems arise systematically: matrix-vector multiplication and solving the set of linear equations

Ax = b, A ∈ R M ×M . (6.1)
For a given matrix A and any vector x, the cost of computing the product between them is 2M 2 operations and standard methods for solving the linear system (6.1) have a computational cost that grows as O(M 3 ). As we mentioned before, the finite difference discretization of the operator A = -div (aν ∇•) on a N × N grid can be written as a system of O(N 2 ) equations. Arranging them in the form (6.1) yields a coefficient matrix with O(N 2 × N 2 ) entries. In particular, since M = N 2 for our case, even for a small number N of grid points the standard methods become unsuitable. By arranging equations (4.1) in a natural row order, one can see that the non zero entries of the resulting matrix are distributed along 5 diagonals of the matrix. In figure 14, we see an example of the distribution pattern for the nonzero entries of A h for N = 10 grid points.

These kind of matrices are called sparse matrices and have the feature that their number of nonzero entries is of order N instead of N 2 . The matrix-vector multiplication for sparse matrices can be carried out in 2s operations where s is the number of nonzero entries of the matrix. For this particular problem, matrix A h has a maximum number of 5N 2 nonzero entries, thus the computational cost is 10N 2 .

Moreover, for our problem, A h is an special case of sparse matrices with a defined structure called banded matrices. The bandwidth µ is the number such that ai,j = 0 if |i -j| > µ. It can be readily seen that for the general case, the FD matrix associated for a given problem, the bandwidth of A h is µ = N . This special structure will allow us to make a precise (but conservative) estimation of the computational cost for the greedy procedure presented here.

Cost of the offline part

The offline part of the greedy algorithm consists of two main ingredients. On one hand, the search for distinguised parameter values νj by examining the residual (3.5) over the set K and, on the other, the computation of the corresponding snapshots (pν j , ȳν j ). We will estimate the computational cost by differentiating three main steps.

Step 1. Choosing ν1. The first parameter is distinguished by maximizing over all of the possible targets y d ν . To do this, one has to the compute k = card( K) times the Euclidian norm in R N 2 and look for the maximum value. The implementation cost is 2kN 2 .

Step 2. Computing (pν 1 , ȳν 1 ). In order to determine the first snapshot one has to solve the minimization problem (3.2) with the selected ν1 in the previous step. As we mentioned before, we are using a standard CG method to compute the optimal control.

The dominating computations during an iteration of the CG are matrix-vector products and the computation of the operator S * S (see eq. (3.3)). The computation of S * Sv for any v ∈ L 2 (ω) amounts first to solve the system of cascade equations

     Ay + c y = χωv, in Ω A * p + c p = β(y -y d ), in Ω y = p = 0 on ∂Ω.
and then we can take S * Sv = p|ω. Since A is self-adjoint both equations can be solved at the numerical level by taking L h = A h + cI h and then finding for y and p. Since this has to be done during many iterations of the CG, we can exploit the nature of the matrix L h (observe that the term cI only contributes to the diagonal and the sparsity pattern of A h is preserved) and take the LU factorization of L h once at the beginning of the process.

For the matrix L h , which has a bandwidth µ = N , the LU factorization gives matrices L and U with the same lower and upper bandwidth, respectively (see [START_REF] Quarteroni | Numerical mathematics[END_REF]Property 3.4]). The computational cost is 2N 4 .

Once the LU factors are known, we can solve for either y or p by forward-backward substitution in 4N 2 µ = 4N 3 operations.

Finally, the maximum number of iterations of the CG to achieve a given tolerance δ (see, for instance, [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]) for this particular case is

≤ 1 2 √ κ ln 2 δ where κ is the condition number of S * S.
Gathering all together, we see that the cost for computing an optimal control ū for a given parameter ν is

Copt = 2N 4 + 4N 3 + O(N 2 ) (6.
2) where the last term involves all related matrix-vector operations of one iteration of the CG. Remark 6.1 The cost provided here is rather conservative and far from being optimal. Observe that, even if A h has a banded pattern, the band itself is sparse (see Figure 14). LU factorization in this case gives matrices L and U which are filled in the band and have more nonzero elements than A h itself. Different order schemes or matrix compacting techniques can lead to a better performance than computing the corresponding LU matrices. In general, the computational cost for solving large sparse system depends on a complicated way of the A h , the number of nonzero elements, its sparsity pattern and the particular algorithm used. We refer the interested reader to [START_REF] George | Computer solution of large sparse positive definite[END_REF] for a more detailed discussion.

Step 3. Choosing νj+1. Suppose that we have determined the first j snapshots and have constructed the approximating space (Pj, Y j ). The basis of the approximating space has been determined gradually through the previous iterations up to the last vector (pν j , ȳν j ). As explained in the previous step, each vector requires Copt operations.

According to Algorithm 2, the next parameter value is chosen by computing

|| inf (p,y)∈(P j ,Y j ) Rν (p, y)|| H -1 (Ω) (6.3)
For any given parameter ν, equations (3.13) and (3.14) enable the practical implementation of this infimum. First, the residual Gν can be build just by vector and scalar multiplication of sparse matrices, see (3.13). The computational cost of this part can be estimated as 45N 2 .

As stated in Remark 2, the infimum appearing in (6.3) can be computed by projecting the vector 0 -βy d to the space Gν (Pj, Y j ). Indeed, one can compute by means of the Gram-Schmidt procedure an orthonormal basis for the space (Pj, Y j ) and then finding the projection just by vector multiplications. To improve the efficiency, this basis is orthonormalized iteratively just by adding the last computed vector to the already orthonormal base. This has a cost of 4N 2 j. Finding the actual projection can be carried out by merely matrix-vector multiplication, which also has a total cost of 4N 2 j.

Adding up, the total cost of this part of the algorithm equals to (k -j)(45N 2 + 8N 2 j) Remark 6.2 Observe that at the numerical level, the residual Rν can be implemented at the very low cost of O(N 2 ) operations, since it is only composed by matrix-vector multiplications. In other works (see, e.g., [START_REF] Lazar | Greedy controllability of finite dimensional linear systems[END_REF]), where the surrogate is obtained by solving coupled forward-backward systems, the cost is dominated by LU-factorization which (for banded sparse matrices) has a cost of 2N 4 (see Step 2 above). The idea of obtaining cheaper surrogates is, of course, not new and has been successfully explored within the framework of optimal control for elliptic and parabolic equations, see, for instance, [START_REF] Grepl | Reduced basis a posteriori error bounds for parametrized linearquadratic elliptic optimal control problems[END_REF][START_REF] Kärcher | A certified reduced basis method for parametrized elliptic optimal control problems[END_REF][START_REF] Kärcher | A Posteriori Error Estimation for Reduced Order Solutions of Parametrized Parabolic Optimal Control Problems[END_REF]. The difference between these works and the approach presented in this mansucript is in the way the greedy procedure is implemented. In the cited papers the authors employ Galerkin projections in order to determine the next snapshot, while we use orthogonal projection to the constructed space. However, the cost of both approaches is similar.

Step 4. Calculating (pj+1, ȳj+1). To compute the new snapshot it suffices to solve the optimal control problem for νj+1. As before, the cost is given in (6.2).

Additional comments and open problems

We devote this section to present some concluding remarks and state some open problems regarding the turnpike property as well as the greedy approach.

1. Robust control of the evolution equation. The application of the turnpike property and the reduction of the optimal control problem for parabolic equations to the steady state one allows us to employ greedy procedure and construct an approximating space for parameterized controls which is independent of the initial datum. This allows for a robust approach which is not datum sensitive. Thus we provide partial answers to the questions arising in [START_REF] Lazar | Greedy controllability of finite dimensional linear systems[END_REF] in the context of finite-dimensional systems, as well as in [START_REF] Choulli | Lipschitz dependence of the coefficients on the resolvent and greedy approximation for scalar elliptic problems[END_REF], where greedy approach for elliptic problems is explored in a robust manner independently of the given source terms.

2. Stable vs unstable. From the results presented in [START_REF] Porretta | Long time versus steady state optimal control[END_REF], it is clear that we can expect the asymptotic simplification of the time-dependent controls towards their steady version if time large. However, from the practical point of view, two main cases have to be differentiated. As discussed in Section 5, the possibility to use the steady control in the evolution problem largely depends on its stability.

If the system is stable, one can effectively use the time independent control to steer the system near the turnpike and remain there indefinitely. When the system is unstable, one can use the approach discussed in Section 5.2 to first control the system to the steady state and then use the steady control to remain there. In this way, the original problem of computing an optimal control in a large time horizon is reduced to an exact controllability problem in a small time interval plus a steady optimal control one.

3. On hypothesis H3. In Section 2, we assumed that the mapping Lν is uniformly elliptic with respect to ν. hold for all ν ∈ K. This condition allowed us to define the residual Rν (and also to employ systematically (2.5)), which in turn is at the heart of the proof of Theorem 3.2 and the greedy procedure described in Algorithm 2.

In the general case, when hypothesis H3 is not fulfilled, there might be some ν ∈ K such that ker(Lν ) = {0}. For that parameter value the solutions to system (2.4) are not uniquely determined and can be written as y = z + e with z ∈ ker(Lν ) ⊥ , e ∈ ker(Lν ). In this case, the functional J appearing in (7.2) is again strictly convex, continuous and coercive and therefore the minimization problem is well posed and the optimal control ū and the optimal state ȳ can be identified. However, from the definition of the residual (3.5), we see that the operator Lν is involved, and if H3 is not fulfiled the residual is not longer a reliable way to measure the distance between two possiblly unknown controls. Thus for the general case in which H3 is not satisfied a more careful analysis is needed and it remains as an open problem.

4. Wave equations. When the free dynamics of the system under consideration enjoys some stability property, it is expected for the turnpike property to hold. For wave-like models, where solutions of the free dynamics are of oscillatory nature and do not enjoy the property of asymptotic simplification, is less clear. However, as shown in [START_REF] Porretta | Long time versus steady state optimal control[END_REF][START_REF] Zuazua | Large time control and turnpike properties for wave equations[END_REF], the turnpike property still holds for wave-like models under suitable controllability assumptions. More precisely, except for an initial time layer [0, τ ] and a final one [0, T -τ ] during the rest of the time interval the solutions are exponentially close to the steady state ones. Therefore, the greedy methodology and the robust control approach discussed above applies for such models without major changes. In the same line as the unstable parabolic case, it remains to show what is the best strategy to control the system in the initial and final layer.

5. Finite element methods. The numerical simulations presented here were implemented by using finite differences methods. As usual, these methods are suitable for analyzing examples posed in simple cartesian geometries with regular enough meshes and are known to be easily programmable. One can use finite element methods to study more complex problems arising in applications and the same results presented in this work are expected. In particular, for problems where the cost of computing one solution to the optimal control problem is very large, the greedy approach and, in general, reduced bases techniques, are reliable tools to speed up the simulations, allowing to explore different parameter configuration within a reasonable amount of time.

6. Greedy controllability for parabolic equations. In this paper we generalize the greedy approach for controllability problems presented in [START_REF] Lazar | Greedy controllability of finite dimensional linear systems[END_REF] for ODEs to the PDE setting. Nevertheless, to quantify the convergence rate of greedy algorithms, the analytic dependence of the controls with respect to the parameters plays a key role. For heat and parabolic equations, when the unknown parameters enter in the principal part of the diffusion operator, the analytic dependence can be expected (this is not the case for wave equations, see, e.g. [START_REF] Allaire | Optimal design of low-contrast two-phase structures for the wave equation[END_REF]). However, the functional setting in which the controls are analytically dependent of the parameters is not clear. Moreover, classical observability inequalities fail to provide good upper and lower bounds for a residual (surrogate) that allows to measure the distance between to possible unknown controls. By exploring the turnpike property, we can bypass the difficulties of analytical dependence for the controls of parabolic problems and reduce the problem to the elliptic setting where we have an appropriate functional setting (cf. subsection 3.4), which provides optimal approximation rates in terms on Kolmogorov widths. Furthermore, unlike other papers in which the affine parameter dependence is required for efficient implementation of a greedy sampling procedure ( [START_REF] Cohen | Approximation of high-dimensional parametric PDEs[END_REF][START_REF] Grepl | Reduced basis a posteriori error bounds for parametrized linearquadratic elliptic optimal control problems[END_REF][START_REF] Kärcher | A certified reduced basis method for parametrized elliptic optimal control problems[END_REF][START_REF] Kärcher | A Posteriori Error Estimation for Reduced Order Solutions of Parametrized Parabolic Optimal Control Problems[END_REF]), here no special dependence of the coefficients with respect to the parameter is required.

Corollary 3 . 8

 38 Let O ∈ R d be an open domain containing the parameter set K and let the mappings aν , cν : O → L ∞ (Ω) be analytic. Then the greedy control algorithm ensures a polynomial decay of arbitrary order of the approximation rates.

Figure 1 :

 1 Figure 1: Numerical results of the greedy test #1.

  The control set ω

Figure 2 :

 2 Figure 2: Approximated control by greedy methods.

2 andFigure 3 : 2 Figure 4 :

 2324 Figure 3: Controlled solution.

Figure 5 :

 5 Figure 5: Numerical results of the greedy test #2.

2 Figure 6 :

 26 Figure 6: The approximated control x → ū -5π 2 /2 (x).

Figure 7 :

 7 Figure 7: Controlled solution.

  t ∈ (0.1, 3)

3 Figure 10 :

 310 Figure 10: Evolution in time for an unstable system controlled with the steady control approximated by means of the greedy approach.

Figure 12 :

 12 Figure 12: Turnpike property for an unstable system.

Figure 14 :

 14 Figure 14: Sparsity pattern of matrix A h with n = 10. The total number of nonzero elements is 460 (out of 10000).

(7. 1 )2 u 2 L 2 (ω) + β 2 y -y d 2 L 2 (

 122222 and the minimization problem (2.7) has to be replaced by the more general one min J(u) =1 Ω) , u ∈ L 2 (ω), y ∈ H 1 0 (Ω) verifying (7.1) . (7.2)

The case where aν ≡ const and cν ∈ [-λi+1 + , -λi -], where λi, λi+1 are two adjacent eigenvalues of the constant Laplacian. In this case, (2.6) holds and since we are away from the eigenvalues, Lν is invertible with bounded inverse.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 694126-DyCon). Part of this research was done while the second author visited DeustoTech and Univesity of Deusto with the support of the DyCon project. The second author was also partially supported by Croatian Science Foundation under ConDyS Project, IP-2016-06-2468. The work of the third author was partially supported by the Grants MTM2014-52347, MTM2017-92996 of MINECO (Spain) and ICON of the French ANR.

Total cost

Summing the above cost for j = 1, . . . , n, the tocal cost of the algorithm results in

As the cost Copt of solving the optimal control problem is O(N 4 ), this part contributes with most of the computational cost of the greedy algorithm. As the number of chosen snapshots n approaches the number of eligible parameters k = card( K), this part converges, up to lower order terms, to the cost of computing ūν for all the possible values k. This proves that the application of the greedy control algorithm is always cheaper than a naive approach that consists of computing controls for all values of the parameter from a very refined uniform mesh on K.

Cost of the online part

Here, we will estimate the computational cost of the online part of the greedy procedure, which is described in Algorithm 3. Observe that only the first four steps of this algorithm contribute to the total cost of computing the approximated control ū with the greedy procedure.

• Step 1. In the first step, we have to project the vector

where Gν (Pn, Y n ) denotes the space spanned by Gν (p i , y i ), i = 1, . . . , n and {(p i , y i ), i = 1, . . . , n} are the snapshots selected in the offline part. To this effect we first build the linear part of the residual, Gν , for each pair (pi, ȳi), i = 1, . . . , n, with a computational cost of 45N 2 for each one. Then we orthonormalize the set {Gν (p i , y i ), i = 1, . . . , n} with the Gram-Schmidt procedure, which can be achieved in 2N 2 n 2 operations. Finally, the projection can be carried out by by adding and multiplying vectors in this new orthonormal basis with a total cost of 4N 2 n.

• Step 2. Denoting the projection by Pn 0 -βy d , we solve for the system Pn 0 -βy d = αiGν (p i , y i ).

Solving this system with QR decomposition has a cots of 8N 2 n 2 -16 3 n 3 operations. • Step 3. This step can be achieved by scalar-vector multiplications with a total cost of N 2 n.

• Step 4. The computational cost of this part is negligible with respect to the rest of the algorithm.

Adding up these contributions, the total cost for obtaining an approximative control ū ν for some parameter value equals to

Recall that the the computational cost of obtaining the optimal control by the conjugate gradient method for some parameter ν is Copt (see eq. (6.2)) is of order N 4 . The cost for computing the ū ν only depends on n and, in particular, since no iterations are required the total cost of obtaining an approximation of the control by the greedy procedure is less than the cost of obtaining it, for instance, by the CG method.

Finally, Step 5 of the online algorithm amounts to solve the corresponding system with control ū ν and, as mentioned before, we can solve it by computing the LU factors and then by forward-backward subsstitution in 2N 4 + 4N 3 operations.