Camilla Fiorini

Régis Duvigneau

Jean-Antoine Désidéri

Optimization of an Unsteady System Governed by PDEs using a Multi-Objective Descent Method

Keywords: optimization, unsteady PDE, multi-objective * Opale Project-Team optimisation, EDP instationnaire, multiobjectif

The aim of this work is to develop an approach to solve optimization problems in which the functional that has to be minimized is time dependent. In the literature, the most common approach when dealing with unsteady problems, is to consider a time-average criterion. However, this approach is limited since the dynamical nature of the state is neglected. These considerations lead to the alternative idea of building a set of cost functionals by evaluating a single criterion at different sampling times. In this way, the optimization of the unsteady system is defined as a multi-objective optimization problem, that will be solved using an appropriate descent algorithm. Moreover, we also consider a hybrid approach, for which the set of cost functionals is built by doing a time-average operation over multiple intervals. These strategies are illustrated and applied to a non-linear unsteady system governed by a one-dimensional convection-diffusionreaction partial differential equation.

Introduction

Over the years, optimization has gained an increasing interest in the scientific world, due to the fact that it can be highly valuable in a wide range of applications.

The simplest optimization problem consists in individuating the best point, according to a certain criterion, in a set of admissible points: this can easily be reconducted to finding the minimum (or the maximum) of a given function. However, this model is often too simple to be used in real applications, because of the presence of many different quantities of different nature that interact among themselves, sometimes in a competitive way, sometimes in a synergistic way: this is the framework of multi-objective optimization.

Multi-objective optimization problems arise from many different fields: economics is the one in which, for the first time, has been introduced the concept of optimum in a multi-objective framework as we intend it in this work; the idea was introduced by Pareto in [START_REF] Pareto | Cours d'Economie Politique[END_REF]. Other fields in which multi-objective optimization can be applied are, for instance: game theory problems (for the relationship with the Pareto-concepts, see [START_REF] Rao | Game theory approach for multiobjective structural optimization[END_REF]); operations research, that often deals with problems in which there are many competitive quantities to be minimized; finally, an important field is engineering design optimization (see, for instance, [START_REF] Abou | Multilevel strategies for parametric shape optimization in aerodynamics[END_REF]).

In this work, we propose an algorithm to identify a region of Pareto-optimal points. The starting point is the Multiple Gradient Descent Algorithm, that is a generalization of the classical Steepest-Descent Method. It has been introduced in its first version, MGDA I, in [3]; some first improvements to the original version gave birth to MGDA II (see [4]). Finally, a third version, MGDA III, has been proposed in [START_REF] Désidéri | Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF] from a theoretical point of view. In this work, we implemented MGDA III and we tested it on some simple cases: the testing process highlighted some limitations of the algorithm, therefore we proposed further modifications, generating in this way MGDA III b.

The main topic of this work is to use MGDA to optimize systems governed by unsteady partial differential equations. More precisely, we considered a model problem defined by a onedimensional nonlinear advection-diffusion equation with an oscillatory boundary condition, for which a periodic source term is added locally to reduce the variations of the solution: this could mimic a fluid system including an active flow control device, for instance, and was inspired by [START_REF] Duvigneau | Optimization of a synthetic jet actuator for aerodynamic stall control[END_REF], in which, however, single-objective optimization techniques were used. The objective of this work is to apply MGDA to optimize the actuation parameters (frequency and amplitude of the source term), in order to minimize a cost functional evaluated at a set of times.

In the literature, the most common approach when dealing with unsteady problems, is to consider time-average quantities (for instance, in [START_REF] Duvigneau | Optimization of a synthetic jet actuator for aerodynamic stall control[END_REF] and [1] the cost functional considered is the time-averaged drag force): however, this approach is limited since the dynamical nature of the state is neglected. For this reason, in this work we propose a new approach: the main idea is to built a set of cost functionals by evaluating a single cost functional at different sampling times, and, consequently, applying multi-objective optimization techniques to this set. On one hand, this method is computationally more expensive than a single-objective optimization procedure; on the other hand, it guarantees that the cost functional decrease at every considered instant. Finally, an hybrid approach is proposed, which will be referred to as windows approach: in this case, the set of cost functionals is built by doing a time-average operation over multiple intervals.

The work is organized as follows:

Section 1: in the first part of the section, starting from the definition of a single-objective optimization problem, we define rigorous formulation of a multi-objective optimization one, we recall the Pareto-concepts such as Pareto-stationarity and Pareto-optimality and we analyse the relation between them. Finally we briefly illustrate the classical Steepest-Descent Method in the single-objective optimization framework. In the second part of the section, we recall some basic concepts of differential calculus in Banach spaces and we illustrate a general optimization problem governed by PDEs. Finally, we explain how all these concepts will be used in this work.

Section 2: initially, we briefly explain the first two versions of MGDA; then follows an accurate analysis of MGDA III: we explain the algorithm in details with the help of some flow charts, we introduce and prove some properties of it, and we illustrate the simple test cases and the problems that they put in light. Finally, we propose some modification and we illustrate the new resulting algorithm.

Section 3: first, we illustrate the general solving procedure that has been applied in order to obtain the numerical results: in particular the Continuous Sensitivity Equation method (CSE) is introduced, which has been used instead of the more classical adjoint approach; then, the specific PDEs considered are introduced, both linear and nonlinear, along with the numerical schemes adopted to solve them: in particular, two schemes have been implemented, a first order and a second order one. Finally, the first numerical results are shown: the order of convergence of the schemes is verified in some simple cases whose exact solution is known, and the CSE method is validated. Moreover, in this occasion, the source term and the boundary condition used in most of the test cases are introduced.

Section 4: here, the main numerical results are presented, i.e. the results obtained by applying MGDA to the PDEs systems described in section 3. First MGDA III b is applied to the linear case: these results show how this problem is too simple and, therefore, not interesting. For this reason, we focused on the nonlinear PDEs: to them, we applied MGDA III and III b. The results obtained underlined how the modifications introduced allow a better identification of the Pareto-optimal zone. Then, we applied MGDA III b in a space of higher dimension (i.e. we increased the number of parameters). Finally, we tested the windows approach.

Inria 1 Multi-objective optimization for PDEs constrained problems

In this section, first we recall some basic concepts, like Pareto-stationarity and Pareto-optimality, in order to formulate rigorously a multi-objective optimization problem; then we recall some basic concepts of differential calculus in Banach spaces to illustrate a general optimization problem governed by PDEs.

From optimization to multi-objective optimization

A generic single-objective optimization problem consists in finding the parameters that minimize (or maximize) a given criterion, in a set of admissible parameters. Formally, one can write:

min p∈P J(p) (1)
where J is the criterion, also called cost functional, p is the vector of parameters and P is the set of admissible parameters. Since every maximization problem can easily be reconduced to a minimization one, we will focus of the latter. Under certain conditions (for instance, P convex closed and J convex) one can guarantee existence and uniqueness of the solution of (1), that we will denote as p opt . For more details about the well posedness of this kind of problems, see [START_REF] Nocedal | Numerical Optimization[END_REF], [2]. This format can describe simple problems, in which only one quantity has to be minimized. However, it is often difficult to reduce a real situation to a problem of the kind (1), because of the presence of many quantities of different natures. This can lead to a model in which there is not only one, but a set of criteria {J i } n i=0 to minimize simultaneously, and this is the framework of multi-objective optimization. In order to solve this kind of problems, we need to formalize them, since (1) does not make sense for a set of cost functionals. In fact, in this context it is not even clear what means for a configuration p 0 to be "better" than another one p 1 . First, we give some important definitions that should clarify these concepts.

Basic concepts of multi-objective optimization

The first basic concept regards the relation between two points. Definition 1 (Dominance). The design point p 0 is said to dominate p 1 if:

J i (p 0 ) ≤ J i (p 1 ) ∀i and ∃ k : J k (p 0 ) < J k (p 1 ).
From now on, we consider regular cost functionals such as J i ∈ C 1 (P) ∀i so that it makes sense to write ∇J i (p), where we denoted with ∇ the gradient with respect to the parameters p. Now, we want to introduce the concepts of stationarity and optimality in the multi-objective optimization framework.

Definition 2 (Pareto-stationarity). The design point p 0 is Pareto-stationary if there exist a convex combination of the gradients ∇J i (p 0 ) that is zero, i.e.: In two or three dimension it can be useful to visualize the gradients of the cost functionals as vectors, all applied in the same point. In fact this is a quick and easy way to understand whether or not a point is Pareto-stationary (see Figure 1). Definition 3 (Pareto-optimality). The design point p is Pareto-optimal if it is impossible to reduce the value of any cost functional without increasing at least one of the others. This means that p is Pareto-optimal if and only if p ∈ P that dominates it. The relation between Pareto-optimality and stationarity is explained by the following theorem: Theorem 1. p 0 Pareto-optimal ⇒ p 0 Pareto-stationary.

∃ α = {α i } n i=0 , α i ≥ 0 ∀i, n i=0 α i = 1 : n i=0 α i ∇J i (p) = 0. RR n°8603 ∇J 1 ∇J 2 ∇J 3 ∇J 4 ∇J 5 ∇J 6 ∇J 7 ∇J 1 ∇J 2 ∇J 3 ∇J 4 ∇J 5
Proof. Let g i = ∇J i (p 0 ), let r be the rank of the set {g i } n i=1 and let N be the dimension of the space (g i ∈ R N ). This means that r ≤ min{N, n}, n ≥ 2 since this is multi-objective optimization and N ≥ 1.

If r = 0 the result is trivial. If r = 1, there exist a unit vector u and coefficients β i such that g i = β i u. Let us consider a small increment ε > 0 in the direction u, the Taylor expansion will be:

J i (p 0 + εu) = J i (p 0 ) + ε ∇J i (p 0 ), u + O(ε 2 ) = J i (p 0 ) + εβ i + O(ε 2 ) ⇒ J i (p 0 + εu) -J i (p 0 ) = εβ i + O(ε 2 ) If β i ≤ 0 ∀ i, J i (p 0 + εu) ≤ J i (p 0 ) ∀ i but this is impossible since p 0 is Pareto-optimal.
For the same reason it can not be β i ≥ 0 ∀i. This means that ∃ j, k such that β j β k < 0. Let us assume that β k < 0 and β j > 0. The following coefficients satisfy the Definition 2:

α i =      -β k βj -β k i = k βj βj -β k i = j 0 otherwise then p 0 is Pareto-stationary.
If r ≥ 2 and the gradients are linearly dependent (i.e r < n, with r ≤ N ), up to a permutation of the indexes, ∃ β i such that:

g r+1 = r i=1 β i g i ⇒ g r+1 + r i=1 µ i g i = 0 (2)

Inria

where µ i = -β i . We want to prove that µ i ≥ 0 ∀ i. Let us assume, instead, that µ r < 0 and define the following space:

V := span{g i } r-1 i=1 . Let us now consider a generic vector v ∈ V ⊥ , v = 0. It is possible, since V ⊥ = {0}, in fact: dimV ≤ r -1 ≤ N -1 ⇒ dimV ⊥ = N -dimV ≥ 1.
By doing the scalar product of v with (2) and using the orthogonality, one obtains:

0 = v, g r+1 + r i=1 µ i g i = (v, g r+1 ) + µ r (v, g r ) , then follow that: v, ∇J r+1 (p 0 ) = -µ r v, ∇J r (p 0 ) ,
and since µ r < 0, v, ∇J r+1 (p 0 ) and v, ∇J r (p 0 ) have the same sign. They cannot be both zeros, in fact if they were, they would belong to V , and the whole family {g i } n i=1 would belong to V , but this is not possible, since:

dimV ≤ r -1 < r = rank{g i } n i=1 . Therefore, ∃ v ∈ V ⊥ , such that v, ∇J r+1 (p 0 ) = -µ r v, ∇J r (p 0
) is not trivial. Let us say that v, ∇J r+1 (p 0 ) > 0, repeating the argument used for r = 1 one can show that -v is a descent direction for both J r and J r+1 , whereas leaves the other functionals unchanged due to the orthogonality, but this is a contradiction with the hypothesis of Pareto-optimality. Therefore µ k ≥ 0 ∀ k and, up to a normalization constant, they are the coefficients α k of the definition of Pareto-stationarity. Finally, this leave us with the case of linearly independent gradients, i.e. r = n ≤ N . However, this is incompatible with the hypothesis of Pareto optimality. In fact, one can rewrite the multi-objective optimization problem as follows:

min p∈P J i (p) subject to J k (p) ≤ J k (p 0 ) ∀ k = i. (3)
The Lagrangian formulation of the problem (3) is:

L i (p, λ) = J i (p) + k =i λ k J k (p) -J k (p 0 ) ,
where λ k ≥ 0 ∀k = i. Since p 0 is Pareto optimal, it stands:

0 = ∇ p L i | (p 0 ,λ) = ∇ p J i (p 0 ) + k =i λ k ∇ p J k (p 0 ),
but this is impossible since the gradients are linearly independent. In conclusion, it cannot be r = n, and we always fall in one of the cases previously examinated, which imply the Paretostationarity.

Unless the problem is trivial and there exist a point p that minimizes simultaneously all the cost functionals, the Pareto-optimal design point is not unique. Therefore, we need to introduce the concept of Pareto-front. A Pareto-front represents a compromise between the criteria. An example of Pareto-front in a case with two cost functionals is given in Figure 2: each point corresponds to a different value of p.

Given these definitions, the multi-objective optimization problem will be the following:

given P, p start ∈ P not Pareto-stationary, and

{J i } n i=0 find p ∈ P : p is Pareto-optimal. ( 4
)
To solve the problem (4), different strategies are possible, for instance one can build a new cost functional:

J = n i=0 α i J i α i > 0 ∀i,
and apply single-objective optimization techniques to this new J. However, this approach presents heavy limitations due to the arbitrariness of the weights. See [START_REF] Giacomini | Multiple-gradient descent algorithm for isogeometric shape optimization[END_REF] for more details. In this work we will consider an algorithm that is a generalization of the steepest descent method: MGDA, that stands for Multiple Gradient Descent Algorithm, introduced in its first version in [3] and then developed in [4] and [START_REF] Désidéri | Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF]. In particular we will focus on the third version of MGDA, of which we will propose some improvements, generating in this way a new version of the algorithm: MGDA-III b.

Steepest Descent Method

In this section, we recall briefly the steepest descent method, also known as gradient method. We will not focus on the well posedness of the problem, neither on the convergence of the algorithm, since this is not the subject of this work. See [START_REF] Quarteroni | Numerical mathematics[END_REF] for more details.

Inria

Let F (x) : Ω → R be a function that is differentiable in Ω and let us consider the following problem:

min x∈Ω F (x).
One can observe that:

∃ρ > 0 : F (a) ≥ F (a -ρ∇F (a)) ∀a ∈ Ω.
Therefore, given an initial guess x 0 ∈ Ω, it is possible to build a sequence {x 0 , x 1 , x 2 , . . .} in the following way:

x n+1 = x n -ρ∇J(x n ) such that F (x 0 ) ≥ F (x 1 ) ≥ F (x 2 ) ≥ . . ..
Under certain conditions on F , Ω and x 0 (for instance, Ω convex closed, F ∈ C(Ω) and x 0 "close enough" to the minimum), it will be:

lim n→+∞ x n = x = argmin x∈Ω F (x).
1.2 Optimization governed by PDEs

Basic concepts

We start this section by introducing briefly the notion of a control problem governed by PDEs, of which the optimization is a particular case. For a more detailed analysis, see [START_REF] Lions | Optimal Control of Systems governed by Partial Differential Equations[END_REF], [START_REF] Tröltzsch | Optimal Control of Partial Differential Equations: Theory, Methods, and Applications[END_REF].

In general, a control problem consists in modifying a system (S), called state, in our case governed by PDEs, in order to obtain a suitable behaviour of the solution of (S). Therefore, the control problem can be expressed as follows: finding the control variable η such that the solution u of the state problem is the one desired. This can be formalized as a minimization problem of a cost functional J = J(η) subject to some constraints, described by the system (S).

If the system (S) is described by a PDE, or by a system of PDEs, the control variable η normally influences the state through the initial or boundary condition or through the source term.

We will focus on problems of parametric optimization governed by PDEs: in this case, the control variable η is a vector of parameters p = (p 1 , . . . , p N ) ∈ R N . Moreover, in this work we will consider only unsteady equations and the optimization parameters will be only in the source term. Therefore, the state problem will be:

∂u ∂t + L(u) = s(x, t; p) x ∈ Ω, 0 < t ≤ T + b.c. and i.c., (5)
where L is an advection-diffusion differential operator, and the optimization problem will be:

min p∈P J(u(p)) subject to (5), (6)
where P ⊆ R N is the space of admissible parameters. Note that J depends on p only through u.

Theoretical results

In this subsection, we will present some theoretical results about control problems. For a more detailed analysis, see [START_REF] Salsa | Note su ottimizzazione in spazi di Banach e Hilbert e sul controllo per EDP[END_REF]. First, we recall some basic concepts of differential calculus in Banach spaces. Let X and Y be two Banach spaces, and let L(X, Y ) be the space of linear functionals from X to Y .

Definition 5. F : U ⊆ X → Y , with U open, is Fréchet differentiable in x 0 ∈ U , if there exists an application L ∈ L(X, Y ), called Fréchet differential, such that, if x 0 + h ∈ U , F (x 0 + h) -F (x 0 ) = Lh + o(h),
that is equivalent to:

lim h→0 F (x 0 + h) -F (x 0 ) -Lh Y h X = 0.
The Fréchet differential is unique and we denote it with: dF (x 0 ). If F is differentiable in every point of U , it is possible to define an application dF : U → L(X, Y ), the Fréchet derivative:

x → dF (x).

We say that F ∈ C 1 (U ) if dF is continuous in U, that is:

dF (x) -dF (x 0 ) L(X,Y ) → 0 if x -x 0 X → 0.
We now enunciate a useful result, known as chain rule.

Theorem 2. Let X, Y and Z be Banach spaces.

F : U ⊆ X → V ⊆ Y and G : V ⊆ Y → Z, where U and V are open sets. If F is Fréchet differentiable in x 0 and G is Fréchet differentiable in y 0 , where y 0 = F (x 0 ), then G • F is Fréchet differentiable in x 0 and: d(G • F )(x 0 ) = dG(y 0 ) • dF (x 0 ).
We now enunciate a theorem that guarantees existence and uniqueness of the solution of the control problem, when the control variable is a function η living in a Banach space H (more precisely η ∈ H ad ⊆ H, H ad being the space of admissible controls). Theorem 3. Let H be a reflexive Banach space, H ad ⊆ H convex closed and J : H → R. Under the following hypotheses:

1. if H ad is unbounded, then J(η) → +∞ when η H → +∞, 2. inferior semicontinuity of J with respect to the weak convergence, i.e.

η j η ⇒ J(η) ≤ lim inf J(η j ) j → +∞.
There exists η ∈ H ad such as J(η) = min η∈H ad J(η). Moreover, if J is strictly convex than η is unique.

Let us observe that the Theorem 3 guarantees the well posedness of a control problem. The well posedness of a parametric optimization problem as [START_REF] Duvigneau | Optimization of a synthetic jet actuator for aerodynamic stall control[END_REF] is not a straightforward consequence: it is necessary to add the request of differentiability of u with respect to the vector of parameters p.

Objective of this work

In this work, as we already said, we will consider only unsteady equation: this will lead to time dependent cost functional J(t). A common approach in the literature is to optimize with respect to a time-average over a period, if the phenomena is periodic, or, if it is not, over an interesting Inria window of time. Therefore, chosen a time interval (t, t + T ), one can build the following cost functional:

J = t+T t J(t)dt, (7)
and apply single-objective optimization techniques to it. However this approach, although straightforward, is limited since the dynamical nature of the state is neglected. Considering only a time-averaged quantity as optimization criterion may yield undesirable effects at some times and does not allow to control unsteady variations, for instance. The objective of this work is thus to study alternative strategies, based on the use of several optimization criteria representing the cost functional evaluated at some sampling times:

J i = J(t i ) for i = 1, . . . , n. (8)
To this set of cost functionals {J i } n i=1 we will apply the Multiple Gradient Descent Algorithm and Pareto-stationarity concept, already mentioned in Subsection 1.1.1, that will be described in details in the next section. Moreover, an other approach will be investigated in this work that will be referred to as windows approach. Somehow, it is an hybrid between considering a time average quantity like [START_REF] Abou | Multilevel strategies for parametric shape optimization in aerodynamics[END_REF] and a set of instantaneous quantities like [START_REF] Giacomini | Multiple-gradient descent algorithm for isogeometric shape optimization[END_REF]. The set of cost functional which the MGDA will be applied to in this case is:

J i = ti+1 ti J(t)dt for i = 1, . . . , n, (9)
i.e the average operation is split on many subintervals. Since the idea of coupling MGDA to the same quantity evaluated at different times is new, only scalar one-dimensional equation will be considered, in order to focus more on the multiobjective optimization details, that are the subject of this thesis, and less on the numerical solution of PDEs.

Descent algorithms for multi-objective problems

The Multiple Gradient Descent Algorithm (MGDA) is a generalization of the steepest descent method in order to approach problems of the kind (4): according to a set of vectors {∇J(p)} n i=0 , a search direction ω is defined such as -ω is a descent direction ∀J i (p), until a Pareto-stationary point is reached. Definition 6 (descent direction). A direction v is said to be a descent one for a cost functional J in a point p if ∃ ρ > 0 such as:

J(p + ρv) < J(p).
If the cost functional is sufficiently regular with respect to the parameters (for instance if

J ∈ C 2 ), v is a descent direction if (v, ∇J(p)) < 0.
In this section we will focus on how to compute the research direction ω. That given, the optimization algorithm will be the one illustrated in Figure 3.

p 0 , {J i } k = 0 compute ∇J i (p k ) ∀i MGDA ω ? = 0 p k+1 = p k -ρω k = k + 1 p opt = p k = 0 = 0

MGDA I and II: original method and first improvements

The original version of MGDA defines ω as the minimum-norm element in the set U of the linear convex combinations of the gradients:

U = {u ∈ R N : u = n i=1 α i ∇J i (p), with α i ≥ 0 ∀i, n i=1 α i = 1}. (10)
In order to understand intuitively the idea behind this choice, let us analyse some two and three dimensional examples, shown in Figure 4 and Figure 5. The gradients can be seen as vectors of R N all applied in the origin. Each vector identifies a point x i ∈ R N , and d is the point identified by ω. In this framework, U is the convex hull of the points x i , highlighted in green in the Figures. The convex hull lies in an affine subspace of dimension at most n -1, denoted with A n-1 and drawn in red. Note that if n > N , A n-1 ≡ R N . We briefly recall the definition of projection on a set. Definition 7. Let x be a point in a vector space V and let C ⊂ V be a convex closed set. The projection of x on C is the point z ∈ C such as:

Inria O x 1 x 2 O ⊥ ≡ d ∇J 2 ∇J 1 ω U A 1 O x 1 x 2 ≡ d O ⊥ ∇J 2 ≡ ω ∇J 1 U A 1 Figure 4: 2D examples. O x 1 ∇J 1 x 2 ∇J 2 x 3 ∇J 3 d ≡ O ⊥ ω U O x 1 ∇J 1 x 2 ∇J 2 x 3 ∇J 3 O ⊥ d ω U A 2
z = argmin y∈C x -y V .
Given this definition, it is clear that finding the minimum-norm element in a set is equivalent to finding the projection of the origin on this set. In the trivial case O ∈ U (condition equivalent to the Pareto stationarity), the algorithm gives ω = 0.

Otherwise, if O / ∈ U , we split the projection procedure in two steps: first we project the origin O on the affine subspace A n-1 , obtaining O ⊥ ; then, two scenarios are possible: if O ⊥ ∈ U then d ≡ O ⊥ (as in the left subfigures) and we have ω, otherwise it is necessary to project O ⊥ on U finding d and, consequently, ω (as shown in the right subfigures).

Let us observe that in the special case in which O ⊥ ∈ U , we are not "favouring" any gradients, fact explained more rigorously in the following Proposition.

(ω, ∇J i (p)) = ω 2 ∀i = 1, . . . , n.
Proof. It is possible to write every gradient as the sum of two orthogonal vectors:

∇J i (p) = ω + v i , where v i is the vector connecting O ⊥ and x i . (ω, ∇J i (p)) = (ω, ω + v i ) = = (ω, ω) + (ω, v i ) = (ω, ω) = ω 2 .
The main limitation of this first version of the algorithm is the practical determination of ω when n > 2. For more details on MGDA I, see [3].

In the second version of MGDA, ω is computed with a direct procedure, based on the Gram-Schmidt orthogonalization process (see [START_REF] Strang | Linear algebra ant its applications[END_REF]), applied to rescaled gradients:

J i = ∇J i (p) S i .
However, this version has two main problems: the choice of the scaling factors S i is very important but arbitrary, and the gradients must be linearly independent to apply Gram-Schmidt, that is impossible, for instance, in a situation in which the number of gradients is greater than the dimension of the space, i.e. n > N . Moreover, let us observe that in general the new ω is different from the one computed with the first version of MGDA. For more details on this version of MGDA, see [4]. For these reasons, a third version has been introduced in [START_REF] Désidéri | Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF].

MGDA III: algorithm for large number of criteria

Description

The main idea of this version of MGDA is that in the cases in which there are many gradients, a common descent direction ω could be based on only a subset of gradients. Therefore, one can consider only the most significant gradients, and perform an incomplete Gram-Schmidt process on them. For instance, in Figure 6 there are seven different gradients, but two would be enough to compute a descent direction for all of them. In this contest, the ordering of the gradients is very important: in fact, only the first m out of n will be considered to compute ω, and hopefully

Inria 0 g 1 g 2 g 3 g 4 g 5 g 6 g 7 0 g 1 g 2 g 3 g 4 g 5 g 6 g 7
Figure 6: 2D example of trends among the gradients. Two good vectors to compute a common descent direction could be the ones in red.

it will be m n. Let us observe that if the gradients live in a space of dimension N , it will be, for sure, m ≤ N . Moreover, this eliminates the problem of linear dependency of the gradients. In this subsection we will indicate with {g i } n i=1 the set of gradients evaluated at the current design point p, corresponding respectively to the criteria {J i } n i=1 , and with {g (i) } n (i)=1 the reordered gradients, evaluated at p, too. The ordering criterion proposed in [START_REF] Désidéri | Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF] selects the first vector as follows:

g (1) = g k where k = argmax i=1,...,n min j=1,...,n (g i , g j ) (g i , g i ) , (11)
while the criterion for the following vectors will be shown afterwards. The meaning of this choice will be one of the subject of the next subsection.

Being the ordering criterion given, the algorithm briefly consists in: building the orthogonal set of vectors {u i } m i=1 by adding one u i at a time and using the Gram-Schmidt process; deciding whether or not to add another vector u i+1 (i.e. stopping criterion); and, finally, computing ω as the minimum-norm element in the convex hull of {u 1 , u 2 , . . . , u m }:

ω = m k=1 α k u k , where α k = 1 1 + m j=1 j =k u k 2 u j 2 . ( 12
)
The choice of the coefficients α k will be justified afterwards. Now we analyse in detail some of the crucial part of the algorithm, with the support of the flow chart in Figure 7. The algorithm is initialized as follows: first g (1) is selected (according to ( 11)), then one impose u 1 = g (1) and the loop starts with i = 2. In order to define the stopping criterion of this loop (TEST in Figure 7), we introduce an auxiliary structure: the lower-triangular matrix C. Let c k, be the element of the matrix C in the k th row, th column.

c k, =        (g (k) ,u ) (u ,u ) < k k-1 =1 c k, = k (13)
This matrix is set to 0 at the beginning of the algorithm, then is filled one column at a time (see Figure 7). The main diagonal of the matrix C contains the sum of the elements in that row and is updated every time an element in that row changes. Given a threshold value a ∈ (0, 1), chosen Fiorini, Duvigneau & Désidéri

{g j } n j=1 select g (1) u 1 = g (1) i = 2 fill(i -1) th column of C TEST pos. neg. compute ω ω i = i + 1 u i = 0? = 0 = 0 compute u i select g (i)
P-S TEST pos.

?

ω = 0 ω go to MGDA-I Figure 7: Flow chart of third version of MGDA
by the user and whose role will be studied in subsection 2.3.1, the test consists in checking if c , ≥ a ∀ > i: if it is true for all the coefficient the test is positive, the main loop is interrupted and ω is computed as in [START_REF] Nocedal | Numerical Optimization[END_REF] with m being the index of the current iteration, otherwise one continues adding a new u i computed as follows: first there is the identification of the index such as: = argmin j=i,...,n

c j,j , (14)
then it is set g (i) = g (these two steps correspond to "select g (i) " in Figure 7), and finally the Gram-Schmidt process is applied:

u i = 1 A i g (i) - k<i c i,k u k , (15)
where

A i := 1 - k<i c i,k = 1 -c i,i .
The meaning of the stop criterion is the following:

c , ≥ a > 0 ⇒ ∃ i : c ,i > 0 ⇒ g ( ) u i ∈ - π 2 , π 2
This means that a descent direction for J ( ) can be defined from the {u i } already computed, and adding a new u computed on the base of g ( ) would not add much information. Once computed a new u i , before adding it to the set, we should check if it is nonzero, as it is shown in the flow chart. If it is zero, we know from the (15) that is:

g (i) = k<i c i,k u k ( 16
)

Inria

It is possible to substitute in cascade ( 15) in ( 16), obtaining:

g (i) = k<i c i,k 1 A k g (k) - <k c k, u = = k<i c i,k   1 A k   g (k) - <k c k,   1 A   g ( ) - j< c ,j u j         = = . . . (17)
Therefore, since

u 1 = g (1)
, one can write:

g (i) = k<i c i,k g (k) (18)
Computing the c i,k is equivalent to solving the following linear system:

      g (1) g (2) • • • g (i-1)       =A         c 1,i c 2,i . . . . . . c i-1,i         =x =       g (i)       =b
Since N is the dimension of the space in which the gradients live, i.e.

g i ∈ R N ∀i = 1, . . . , n, the matrix A is [N × (i -1)],
and for sure is i -1 ≤ N , because i -1 is the number of linearly independent gradients we already found. This means that the system is overdetermined, but we know from ( 17)-( 18) that the solution exists. In conclusion, to find the c i,k we can solve the following system:

A T Ax = A T b ( 19
)
where A T A is a square matrix, and it is invertible because the columns of A are linearly independent (for further details on linear systems, see [START_REF] Strang | Linear algebra ant its applications[END_REF]). Once found the c i,k , it is possible to apply the Pareto-stationarity test (P-S TEST in Figure 7), based on the following Proposition.

Proposition 5. c i,k ≤ 0 ∀k = 1, . . . , i -1 ⇒ Pareto-stationarity.
Proof. Starting from the (18), we can write:

g (i) - k<i c i,k g (k) = 0 ⇒ n k=1 γ k g (k) = 0
where:

γ k =      -c i,k k < i 1 k = i 0 k > i γ k ≥ 0 ∀k. If we define Γ = n k=1 γ k , then the coefficients { γ k Γ } n k=1 satisfy the Definition 2.
Note that the condition given in Proposition 5 is sufficient, but not necessary. Therefore, if it is not satisfied we fall in an ambiguous case: the solution suggested in [START_REF] Désidéri | Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF] is, in this case, to compute ω according to the original version of MGDA, as shown in Figure 7.

Properties

We now present and prove some properties of MGDA III. Note that these properties are not valid if we fall in the ambiguous case. Proposition 6. α i , i = 1, . . . , n defined as in [START_REF] Nocedal | Numerical Optimization[END_REF] are such that:

m k=1 α k = 1 and (ω, u k ) = ω 2 ∀k = 1, . . . , m.
Proof. First, let us observe that:

1 α k = 1 + m j=1 j =k u k 2 u j 2 = u k 2 u k 2 + m j=1 j =k u k 2 u j 2 = m j=1 u k 2 u j 2 = u k 2 m j=1 1 u j 2
Therefore, it follows that:

m k=1 α k = m k=1 1 u k 2 m j=1 1 u j 2 = m k=1 1/ u k 2 m j=1 1 u j 2 = 1 m j=1 1 u j 2 m k=1 1 u k 2 = 1.
Let us now prove the second part:

(ω, u i ) = m k=1 α k u k , u i = α i (u i , u i ) = α i u i 2 = 1 m k=1 1 u k 2 =: D
where we used the definition of ω, the linearity of the scalar product, the orthogonality of the u i and, finally, the definition of α i . Let us observe that D does not depend on i. Moreover,

ω 2 = (ω, ω) = m k=1 α k u k , m =1 α u = m j=1 α 2 j (u j , u j ) = = m j=1 α j D = D m j=1 α j = D = (ω, u i ) ∀i = 1, . . . , m. Proposition 7. ∀i = 1, . . . , m (ω, g (i) ) = (ω, u i ).
Proof. From the definition of u i in [START_REF] Quarteroni | Numerical mathematics[END_REF], it is possible to write g (i) in the following way:

g (i) = A i u i + k<i c i,k u k = 1 - k<i c i,k u i + k<i c i,k u k = = u i - k<i c i,k u i + k<i c i,k u k

Inria

and doing the scalar product with ω one obtains:

(g (i) , ω) = (u i , ω) - k<i c i,k (u i , ω) + k<i c i,k u k , ω =( )
We want to prove that ( ) = 0. As we have shown in the proof of Proposition 6:

α k = 1 u k 2 m j=1 1 u j 2 ⇒ α i u i 2 = 1 m j=1 1 u j 2
Using the definition of ω, the linearity of the scalar product and the orthogonality of the u i one obtains:

( ) = - k<i c i,k α i u i 2 + k<i c i,k α k u k 2 = = - k<i c i,k 1 m j=1 1 uj 2 + k<i c i,k 1 m j=1 1 uj 2 = 0. Proposition 8. ∀i = m + 1, . . . , n (g (i) , ω) ≥ a ω 2 .
Proof. It is possible to write g (i) as follows:

m k=1 c i,k u k + v i
where v i ⊥ {u 1 , . . . , u m }. Since i > m, we are in the case in which the test c i,i ≥ a is positive. Consequently:

(g (i) , ω) = m k=1 c i,k (u k , ω) = m k=1 c i,k ω 2 = c i,i ω 2 ≥ a ω 2 .
In conclusion, as a straightforward consequence of the Propositions above, we have:

(g (i) , ω) = ω 2 ∀i = 1, . . . , m (g (i) , ω) ≥ a ω 2 ∀i = m + 1, . . . , n,
which confirms that the method provides a direction that is a descent one even for the gradients not used in the construction of the Gram-Schmidt basis and, consequently, not directly used in the computation of ω.

MGDA III b: Improvements related to the ambiguous case

In this work we propose some strategies to avoid going back to MGDA-version I, when we fall in the ambiguous case, since it is very unefficient when a large number of objective functional are considered. In order to do this, we should investigate which are the reasons that lead to the ambiguous case. Let us observe that in the algorithm there are two arbitrary choices: the value of the threshold a for the stopping test and the algorithm to order the gradients.

0 g 1 g 2 g 3 R 2 0 g 1 g 2 ≡ g (1) ≡ u 1 g 3 u 2 R 2

Choice of the threshold

Let us start with analysing how a too demanding choice of a (i.e. of a too large value) can affect the algorithm. We will do this with the support of the example in Figure 8, following the algorithm as illustrated in the precedent subsection. For example, assume that the gradients are the following:

g 1 = 5 2 g 2 = 0 5 g 3 = -5 2 . ( 20
)
Using the criterion (11) the gradient selected as first is g 2 . At the first iteration of the loop the first column of the matrix C is filled in as in [START_REF] Pareto | Cours d'Economie Politique[END_REF], with the following result:

C =   1 0.4 0.4 0.4 0 0.4  
If a is set, for instance, at 0.5 (or any value greater than 0.4) the stop test fails. However, looking at the picture, it is intuitively clear that ω = u 1 is the best descent direction one can get (to be more rigorous, it is the descent direction that gives the largest estimate possible from Proposition 8). If the test fails, the algorithm continues and another vector u i is computed. This particular case is symmetric, therefore the choice of g (2) is not very interesting. Let us say that g (2) = g 3 , than u 2 is the one shown in the right side of Figure 20 and, since it is nonzero, the algorithm continues filling another column of C:

C =   1 0.4 0.4 0.4 -0.6 -0.2  
Even in this case the stopping test fails (it would fail with any value of a due to the negative term in position (3,3)). Since we reached the dimension of the space, the next u i computed is zero: this leads to the Pareto-stationarity test, that cannot give a positive answer, because we are not at a Pareto-stationary point. Therefore, we fall in the ambiguous case.

To understand why this happens, let us consider the relation between a and ω, given by Proposition 8 that we briefly recall:

∀i = m + 1, . . . , n (g (i) , ω) ≥ a ω 2 .
This means that, for the gradients that haven't been used to compute ω, we have a lower bound for the directional derivative along ω: they are, at least, a fraction a of ω 2 . The problem is Inria that it is impossible to know a priori how big a can be. For this reason, the parameter a has been removed. The stopping test remains the same as before, with only one difference: instead of checking c , ≥ a one checks if c , > 0.

However, removing a is not sufficient to remove completely the ambiguous case, as it is shown in the next subsubsection.

Ordering criteria

In this subsection we will analyse the meaning of the ordering criterion used in MGDA III and we will present another one, providing some two-dimensional examples in order to have a concrete idea of different situations.

The ordering problem can be divided in two steps:

(i) choice of g (1) (ii) choice of g (i) , given g (k) ∀k = 1, . . . i -1
The naïve idea is to order the gradients in the way that makes the MGDA stop as soon as possible. Moreover, we would like to find an ordering criterion that minimizes the times in which we fall in the ambiguous case (i.e. when we find a u i = 0 and the P-S test is not satisfied). Definition 8. Let us say that an ordering is acceptable if it does not lead to the ambiguous case, otherwise we say that the ordering is not acceptable. First, we discuss point (ii), since it is simpler. If we are in the case in which we have to add another vector, it is because the stopping test gave a negative result, therefore:

∃ , i ≤ ≤ n : c , < 0.
We chose as g (i) the one that violates the stopping test the most. This explains the [START_REF] Quarteroni | Numerical models for differential problems[END_REF], that we recall here:

g (i) = g k : k = argmin c , i ≤ ≤ n.
The point (i) is more delicate, since the choice of the first vector influences strongly all the following choices. Intuitively, in a situation like the one shown in Figure 9, we would pick as first vector g 5 and we would have the stopping test satisfied at the first iteration, in fact ω = g 5 is a descent direction for all the gradients. Formally, this could be achieved using the criterion 11,

0 g 1 g 2 g 3 g 4 g 5 g 6 g 7 g 8 g 9 0 g 1 g 2 g 3 g 4 g 5 g 6 g 7 g 8 g 9 Figure 9: Example 1. 0 g 1 g 2 g 3 g 4 g 5 0 g 1 g 2 ≡ ω g 3 g 4 g 5
Figure 10: Choice of the first gradient with the criterion ( 11)- [START_REF] Quarteroni | Numerical models for differential problems[END_REF].

0 g 1 g 2 g 3 0 g 1 g 2 ≡ u 1 g 3 u 2 ω
Figure 11: Criterion ( 11)-( 14) leads to a not acceptable ordering.

that we recall:

g (1) = g k , where k = argmax i min j (g j , g i ) (g i , g i ) (21)
However, it is clear that this is not the best choice possible in cases like the one shown in Figure 10, in which none of the gradients is close to the bisecting line (or, in three dimensions, to the axis of the cone generated by the vectors). In this particular case, it would be better to pick two external vectors (g 1 and g 5 ) to compute ω. Anyway, with the criterion ( 11)-( 14) we still obtain an acceptable ordering, and the descent direction is computed with just one iteration. However, there are cases in which the criterion ( 11)-( 14) leads to a not acceptable ordering. Let us consider, for example, a simple situation like the one in Figure 11, with just three vectors, and let us follow the algorithm step by step. Let the gradients be the following:

g 1 = 1 0 g 2 = 0.7 0.7 g 3 = -1 0.1 (22)
Given the vectors (22), we can build a matrix T such as:

[T ] i,j = (g j , g i ) (g i , g i ) (23)
obtaining:

T =   1 0.7 -1 0.714 1 -0.643 -0.99 -0.624 1   min = -1 min = -0.643 min = -0.99 Inria {g j } n j=1 select g (1) u 1 = g (1) i = 2 fill(i -1) th column of C TEST pos. neg. compute ω ω i = i + 1 u i = 0? = 0 = 0 compute u i select g (i)
P-S TEST pos.

? The internal minimum of [START_REF] Lions | Optimal Control of Systems governed by Partial Differential Equations[END_REF] is the minimum of each row. Among this, we should select the maximum: the index of this identifies the index k such as g (1) = g k . In this case, k = 2. We can now fill in the first column of the triangular matrix C introduced in (13), obtaining:

ω = 0 ω compute ω using u k , k < i compute d j = (ω, g j ) ∀j ω d j > 0 ∀j ∀j : d j < 0 g new j = γg j ∃j : d j < 0 {g new j } N j=1
C =   1 0.714 0.714 -0.642 0 -0.642  
Since the diagonal terms are not all positive, we continue adding g (2) = g 3 , according to [START_REF] Quarteroni | Numerical models for differential problems[END_REF]. Note that, due to this choice, the second line of the matrix C has to be switched with the third. The next vector of the orthogonal base u 2 is computed as in [START_REF] Quarteroni | Numerical mathematics[END_REF], and in Figure 11, right subfigure, it is shown the result. Having a new u 2 = 0, we can now fill the second column of C as follows:

C =   1 -0.642 -0.642 0.714 -1.4935 -0.77922  
This is a situation analogue to the one examinated in the example of subsection 2.3.1, therefore we fall in the ambiguous case.

A possible solution to this problem is the one proposed in the flow chart in Figure 12: compute ω with the {u i } that we already have and check for which gradients it is not a descent direction (in this case for g 1 , as shown in the right side of Figure 11). We want this gradient to be "more considered": in order to do that we double their length and we restart the algorithm from the

0 g 1 ≡ u 1 g 2 g 3 u 2 ω
Figure 13: Criterion ( 11)-( 14) after doubling.

beginning. This leads to a different matrix T :

T =   1 0.35 -0.5 1.428 1 -0.643 -1.98 -0.624 1   min = -0.5 min = -0.643 min = -1.98
With the same process described above, one obtains g (1) = u 1 = g 1 . The final output of MGDA is shown in Figure 13, with ω = [0.00222, 0.0666] T . However, rerunning MGDA from the beginning is not an efficient solution, in particular in the cases in which we have to do it many times before reaching a good ordering of the gradients. For this reason, a different ordering criterion is proposed: we chose as a first vector the "most external", or the "most isolated" one. Formally, this means:

g (1) = g k , where k = argmin i min j (g j , g i ) (g i , g i ) (24)
On one hand with this criterion it is very unlikely to stop at the first iteration, but on the other hand the number of situations in which is necessary to rerun the whole algorithm is reduced. Note that, in this case, instead of doubling the length of the non-considered gradients, we must make the new gradients smaller.

These solutions are not the only ones possible, moreover there is no proof that the ambiguous case is completely avoided. For instance, another idea could be the following: keeping the threshold parameter a for the stopping test and, every time the ambiguous case occurs, rerunning the algorithm with a reduced value of a.

Inria 3 Numerical methods for state and sensitivity analyses

In this section, first we will illustrate the general solving procedure that has been applied to solve multi-objective optimization problems with PDE constraints, then we will introduce the specific PDEs used as test cases in this work and finally we will present the first numerical results, related to the convergence of the schemes and the validation of the methods used.

Solving method

In this subsection we want to explain the method used in this work to solve problem of the kind [START_REF] Duvigneau | Optimization of a synthetic jet actuator for aerodynamic stall control[END_REF]. First, let us recall all the ingredients of the problem:

• state equation:          ∂u ∂t + L(u) = s(x, t; p) x ∈ Ω, 0 < t ≤ T ∇u • n = 0 Γ N ⊆ ∂Ω, 0 < t ≤ T u = f (t) Γ D ⊆ ∂Ω, 0 < t ≤ T u = g(x)
x ∈ Ω, t = 0

• cost functional J that we assume quadratic in u:

J(u(p)) = 1 2 q(u, u), (25)
where q(•, •) is a bilinear, symmetric and continuous form. For instance, a very common choice for J in problems of fluid dynamics is J = 1 2 ∇u 2 L 2 (Ω) . Let us observe that every classical optimization method requires the computation of the gradient of the cost functional with respect to the parameters p, that we will indicate with ∇J. In order to compute it, we can apply the chain rule, obtaining:

∇J(p 0 ) = d(J • u)(p 0 ) = dJ(u 0 ) • du(p 0 ),
being u 0 = u(p 0 ). Therefore, ∇J is composed by two pieces. Let us start analysing the first one. Proposition 9. If J is like in (25), than its Fréchet differential computed in u 0 and applied to the generic increment h is: dJ(u 0 )h = q(u 0 , h)

Proof. Using the definition of Fréchet derivative one obtains:

J(u 0 + h) -J(u 0 ) = 1 2 q(u 0 + h, u 0 + h) - 1 2 q(u 0 , u 0 ) = = q(u 0 , h) + 1 2 q(h, h) = q(u 0 , h) + O( h 2 ),
where the last passage is justified by the continuity of q(•, •).

The second piece, i.e. du(p 0 ), is just the derivative of u with respect to the parameters. This means that the i -th component of the gradient of J is the following:

[∇J( p 0 )] i = q u 0 , u pi (p 0 ) , (26)
where

u pi (p 0 ) = ∂u ∂p i (p 0 ),
and it is called sensitivity (with respect to the parameter p i ). We now apply this procedure to the example introduced above, J = 1 2 ∇u 2 L 2 (Ω) . We observe that the form q is:

q(u, v) = (∇u, ∇v) L 2 (Ω) ,
that is bilinear, symmetric and continuous. Therefore we can apply the result of Proposition 9, obtaining:

[∇J( p 0 )] i = ∇u 0 , ∇u pi (p 0 ) L 2 (Ω) .
Summing up: in order to solve a minimization problem like (6), we need to compute the gradient of J with respect to the parameters. To do that, we see from (26) that it is necessary to know u 0 = u(p 0 ) (therefore to solve the state equation with p = p 0 ) and to know u pi . The last step we have to do is finding the equations governing the behaviour of the sensitivities u pi . To do that, we apply the Continuous Sensitivity Equation method, shorten in CSE method from now on. See [START_REF] Kozub | Interactive and immersive simulation using the continuous sensitivity equation method with application to thermal transport[END_REF]. It consists in formally differentiating the state equation, the boundary and initial condition, with respect to p i :

           ∂ ∂pi ∂u ∂t + L(u) = ∂ ∂pi s(x, t; p) x ∈ Ω, 0 < t ≤ T ∂ ∂pi (∇u • n) = 0 Γ N ⊆ ∂Ω, 0 < t ≤ T ∂ ∂pi u = 0 Γ D ⊆ ∂Ω, 0 < t ≤ T ∂ ∂pi u = 0 x ∈ Ω, t = 0.
After a permutation of the derivatives, if L is linear and Ω does not depend on p, one obtains the sensitivities equations:

         ∂up i ∂t + L(u pi ) = s pi (x, t; p) x ∈ Ω, 0 < t ≤ T ∇u pi • n = 0 Γ N ⊆ ∂Ω, 0 < t ≤ T u pi = 0 Γ D ⊆ ∂Ω, 0 < t ≤ T u pi = 0 x ∈ Ω, t = 0. (27)
For every i = 1, . . . , N , this is a problem of the same kind of the state equation and independent from it, therefore from a computational point of view it is possible to use the same solver and to solve the N + 1 problems in parallel (N sensitivity problems plus the state). However, if L is nonlinear,

L(u pi ) = ∂ ∂p i L(u),
and we need to introduce a new operator L = L(u, u pi ) = ∂ ∂pi L(u). Therefore, the sensitivities equations in this case will be:

         ∂up i ∂t + L(u, u pi ) = s pi (x, t; p) x ∈ Ω, 0 < t ≤ T ∇u pi • n = 0 Γ N ⊆ ∂Ω, 0 < t ≤ T u pi = 0 Γ D ⊆ ∂Ω, 0 < t ≤ T u pi = 0 x ∈ Ω, t = 0. (28)
They are no more independent from the state equation, but they remain independent ones from the others. Let us observe that in this case it is necessary to implement a different (but similar) solver than the one used for the state.

Inria

Test cases

In this work we will consider only unsteady one-dimensional scalar advection-diffusion PDEs, both linear and nonlinear. Let Ω = (x a , x b ) be the domain. The linear equation will be the following:

           ∂ t u -b∂ 2 x u + c∂ x u = s(x, t; p) in (x a , x b ), 0 < t ≤ T u(x a , t) = f D (t) 0 < t ≤ T ∂ x u(x b , t) = f N (t) 0 < t ≤ T u(x, 0) = g(x) in (x a , x b ) (29)
where b and c are constant coefficients b, c > 0, c b (i.e. dominant advection), s(x, t), f D (t), f N (t) and g(x) are given functions and p is the vector of control parameters. Being c > 0 it makes sense to impose Dirichlet boundary condition in x a , since it is the inlet, and Neumann boundary condition in x b , since it is the outlet. We observe that the problem (29) is well posed, i.e. the solution exists and is unique. For more details on the well posedness of this kind of problems see [START_REF] Salsa | Equazioni a derivate parziali: Metodi, modelli e applicazioni[END_REF], [START_REF] Quarteroni | Numerical models for differential problems[END_REF].

Applying the CSE method, introduced in subsection 3.1, one obtains the following sensitivity equations, for every i = 1, . . . , N :

           ∂ t u pi -b∂ 2 x u pi + c∂ x u pi = s pi (x, t; p) in (x a , x b ), 0 < t ≤ T u pi (x a , t) = 0 0 < t ≤ T ∂ x u pi (x b , t) = 0 0 < t ≤ T u pi (x, 0) = 0 in (x a , x b ) (30)
Let us observe that the sensitivity equations (30) are independent one from the other and all from the state. Moreover, the same solver can be used for all of them. On the other hand, for the nonlinear case we consider the following equation:

           ∂ t u -b∂ 2 x u + cu∂ x u = s(x, t; p) in (x a , x b ), 0 < t ≤ T u(x a , t) = f D (t) 0 < t ≤ T ∂ x u(x b , t) = f N (t) 0 < t ≤ T u(x, 0) = g(x) in (x a , x b ). (31)
We will impose source term, initial and boundary condition such as u(x, t) > 0 ∀x ∈ (x a , x b ), ∀t ∈ (0, T ). In this way the advection field cu is such that x a is always the inlet and x b the outlet, and we can keep the same boundary condition and, most important thing, the same upwind scheme for the convective term, introduced in the next subsection. However, in this case the sensitivity equations are:

           ∂ t u pi -b∂ 2 x u pi + cu∂ x u pi + cu pi ∂ x u = s pi (x, t; p) in (x a , x b ), 0 < t ≤ T u pi (x a , t) = 0 0 < t ≤ T ∂ x u pi (x b , t) = 0 0 < t ≤ T u pi (x, 0) = 0 in (x a , x b ). (32)
These are linear PDEs of reaction-advection-diffusion with the advection coefficient cu(x, t) and the reaction coefficient c∂ x u(x, t) depending on space and time. In this case it is not possible to use the same solver for state and sensitivity, moreover the sensitivity equations depend on the solution of the state.

Numerical schemes

We now introduce the numerical schemes used to solve the problems (29)-( 30)-( 31)-(32). In every case, the space has been discretized with a uniform grid {x i } N i=0 where x i = x a + i∆x and ∆x is the spatial step, that can be chosen by the user.

x 0 x 1 x 2 x 3 x N -1 x N . . . ∆x
Also in time, we considered a uniform discretization: {t n } K n=0 where t n = n∆t and ∆t is chosen such as two stability conditions are satisfied:

∆t < 1 2 ∆x 2 b ∧ ∆t < ∆x c . (33)
These conditions (33) are necessary, since we use explicit methods in time, as shown in the next subsubsections.

A first order and a second order scheme have been implemented, by using finite differences in space for both orders, while in time explicit Euler for first order, and a Runge-Kutta scheme for second order. See [START_REF] Quarteroni | Numerical models for differential problems[END_REF] and [START_REF] Leveque | Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems[END_REF] for more details on these schemes. The resulting schemes are the following: all of them are initialized using the initial condition and the Dirichlet boundary condition:

u 0 i = g(x i ) ∀i = 0, . . . , N u n 0 = f D (t n ) ∀n = 0, . . . , K.
• First order scheme for the linear case ∀ n = 1, . . . , K, ∀ i = 1, . . . , N -1:

u n+1 i = u n i -∆t c u n i -u n i-1 ∆x + ∆t b u n i+1 -2u n i + u n i-1 ∆x 2 + ∆t s(x i , t n ) (34)
Neumann boundary condition i = N :

u n+1 N = u n+1 N -1 + ∆xf N (t n+1 )
• Second order scheme for the linear case ∀ n = 1, . . . , K, ∀ i = 2, . . . , N -1:

             u n+ 1 2 i = u n i -∆t 2 c 3u n i -4u n i-1 +u n i-2 2∆x + ∆t 2 b u n i+1 -2u n i +u n i-1 ∆x 2 + ∆t 2 s(x i , t n ) u n+1 i = u n i -∆t c 3u n+ 1 2 i -4u n+ 1 2 i-1 +u n+ 1 2 i-2 2∆x + +∆t b u n+ 1 2 i+1 -2u n+ 1 2 i +u n+ 1 2 i-1 ∆x 2 + ∆t s(x i , t n+ 1 2 ) (35)
Neumann boundary condition i = N :

u n+1 N = 1 3 (4u n+1 N -1 -u n+1 N -2 + 2∆x f N (t n+1 ))
Since the numerical dependency stencil for the internal nodes (∀i = 1, . . . , N -1) is the following:

Inria u n+ 1 2 i-2 u n+ 1 2 i-1 u n+ 1 2 i u n+ 1 2 i+1 u n+1 i ∆t 2
∆x for the computation of u n+1 1 it is necessary to know the value of the solution in a point outside the domain. Different choices for this node lead to different schemes. Unless otherwise indicated, in this work we decided to impose u n -1 = u n 0 .

With these two schemes (34)-( 35) it is possible to solve the linear state equation ( 29) and the linear sensitivity equations (30).

For the nonlinear state equation (31) the scheme is very similar to the one shown above, the only difference is in the advection term.

• First order scheme for the nonlinear case ∀ n = 1, . . . , K, ∀ i = 1, . . . , N -1:

u n+1 i = u n i -∆t cu n i u n i -u n i-1 ∆x + ∆t b u n i+1 -2u n i + u n i-1 ∆x 2 + ∆t s(x i , t n ) (36)
• Second order scheme for the nonlinear case ∀ n = 1, . . . , K, ∀ i = 2, . . . , N -1:

             u n+ 1 2 i = u n i -∆t 2 cu n i 3u n i -4u n i-1 +u n i-2 2∆x + ∆t 2 b u n i+1 -2u n i +u n i-1 ∆x 2 + ∆t 2 s(x i , t n ) u n+1 i = u n i -∆t cu n+ 1 2 i 3u n+ 1 2 i -4u n+ 1 2 i-1 +u n+ 1 2 i-2 2∆x + +∆t b u n+ 1 2 i+1 -2u n+ 1 2 i +u n+ 1 2 i-1 ∆x 2 + ∆t s(x i , t n+ 1 2 ) (37)
Note that, since the numerical dependency stencil is the same as in the linear case, we have the same problem for computation of the solution in the first node x 1 .

Finally, to solve the equation (32), the linear schemes can be used, with two simple modification. The first one is that the coefficient c, constant in (34)-( 35), has to be substituted with c n i in (34) and in the first equation of (35), and with c

n+ 1 2 i
in the second equation of (37), being c n i cu(x i , t n ). Note that the equation is linear since we are solving the sensitivity equations in u p k and u(x, t) is given. The second one is that the reaction term c(∂ x u) n i (u p k ) n i must be added. Let us observe that in this case, for the second order scheme, it is necessary to know the state u on a temporal grid twice as fine as the one used for the sensitivity equations.

Verification of the state analysis

In this subsection we will present some numerical convergence results in order to verify the actual order of the finite differences schemes in some simple cases. This subsection is not meant to be exhaustive on the convergence topic, but only a brief check in some simple cases in which it is easy to compute the analytical solution. First let us observe that the discretization of the diffusion term does not change from the first to the second order scheme, therefore we focused on only advection problems.

Advection problem with a steady solution

The first simple case is the following:

     ∂ t u + c ∂ x u = x in (0, 1), 0 < t ≤ T u(0, t) = 0 0 < t ≤ T u(x, 0) = 0 in (0, 1). ( 38
)
Since both the source term and the Dirichlet boundary condition do not depend on time, there will be a steady solution, that is u ex = x 2 2c . If the final time T is big enough, also the numerical solution will be steady. The error is computed as follows:

e(T ) = u ex (x, T ) -u h (x, T ) L 2 (0,1) ,
where u ex is the exact solution and u h the numerical solution. In Figure 14 is shown the error with respect to the grid: the theoretical slope is in red. The first order scheme has the expected convergence. For the second order, two different results are plotted: the green one correspond to the choice of

✶ ✵ ✶ ✁ ✶ ✷ ✶ ✸ ✶ ✹ ✶ ✲✹ ✶ ✲✸ ✶ ✲✷ ✶ ✲✁ ✶ ✵ •✂✄☎✆✝ ✞ ✟ ✠ ✡ ☛ ☞ ☞ ✌ ❋✍✎✏✑ ✂✎✒✓✎ |✓✏✔•✑✏ ❚✕✓✂✎✍✖✗• ✂✎✒✓✎ ✶ ✵ ✶ ✁ ✶ ✷ ✶ ✸ ✶ ✹ ✶ ✲✘ ✶ ✲✙ ✶ ✲✚ ✶ ✲✹ ✶ ✲✸ ✶ ✲✷ ✶ ✲✁ ✶ ✵ ❙✓✖✂✛✒ ✂✎✒✓✎ •✂✄☎✆✝ ✞ ✟ ✠ ✡ ☛ ☞ ☞ ✌ |✓✏✔•✑ ✶ |✓✏✔•✑ ✜ ❚✕✓✂✎✍✖✗• ✂✎✒✓✎
u n -1 = u ex (-∆x, t n ), the blue one to u n -1 = f D (t n ).
As one can notice, the difference is not significant. However, in both cases the order is not the one expected. This can be explained by the fact that numerically we are imposing two boundary condition to a problem that requires only one. The results obtained removing the Neumann boundary condition and imposing in the last node the same scheme of the internal nodes are shown in Figure 15. These results show that the order is the one expected if we use u n -1 = f D (t n ), otherwise, using the exact solution for the external node leads to a problem that is too simple: the error is comparable to the machine error for all the grids. In fact, we are imposing the exact solution in two nodes (x -1 and x 0 ). Moreover, developing by hand the scheme for the computation of u n 1 and imposing u n i = u n+1 i (because of the stationarity), one finds that the scheme yields u 1 = ∆x 2 2c , that is the exact solution in x 1 . Since u ex is a parabola, two nodes inside the domain (plus the external one) are enough to represent the exact solution.

Inria ✶ ✵ ✶ ✁ ✶ ✷ ✶ ✸ ✶ ✹ ✶ ✲✁✂ ✶ ✲✁✄ ✶ ✲✁✹ ✶ ✲✁✷ ✶ ✲✁✵ ✶ ✲✂ ✶ ✲✄ ✶ ✲✹ ✶ ✲✷ ✶ ✵ •☎✆✝✞✟ ✠ ✡ ☛ ☞ ✌ ✍ ✍ ✎ |✏✑✒•✓ ✶ |✏✑✒•✓ ✔ ❚✕✏☎✖✗✘✙•

Advection unsteady problem

The second simple case considered is the following:

     ∂ t u + c ∂ x u = 0 in (0, 1), 0 < t ≤ T u(0, t) = 1 + α sin(f t) 0 < t ≤ T u(x, 0) = 0 in (0, 1). (39)
Let us observe that the function u(x, t) = 1 + α sin(f (ct -x)) satisfies the equation and the boundary condition. Moreover, if the time t is big enough (let us say t > t), the influence of the initial condition is negligible. The results shown in this subsection are obtained with c = 1, f = 3 2 π, t = 3 and T = 6. The error is computed as follows:

e(t) = u ex (x, t) -u h (x, t) L 2 (0,1) err = e(t) L ∞ (t,T ) ,
and it is shown in Figure 16 with respect to the grid. On the left there is the first order scheme, on the right the second order. The theoretical slope is in red: when the grid becomes fine enough, both schemes have the order expected, for every time considered.

Consistency of the sensitivity analysis

The purpose of this subsection is to validate the CSE method introduced in subsection 3.1. To do that, we will consider the problems ( 29) and ( 31) with specific boundary and initial conditions:

f D (t) = k + α sin(2πf t) f N (t) ≡ 0 g(x) ≡ k (40)
where α and f are given parameters, k is big enough to guarantee u(x, t) ≥ 0. We consider the following source term:

s(x, t) = √ A sin(2πf t + ϕ) sin 2 ( x-xc L π -π 2 ), if x ∈ (x c -L 2 , x c + L 2 ) 0, otherwise ( 41
)
where L is the width of the spatial support, x c is the centre of the support and A and ϕ are the design parameters, therefore we have p = (A, ϕ). The source term at some sample times is shown below.

x c 1 L x Source term t 1 t 2 t 3 t 4

Inria

These being the conditions of the state equation, we now introduce the source terms of the sensitivity equations ( 30) and (32):

s A (x, t) = 1 2 √ A sin(2πf t + ϕ) sin 2 ( x-xc L π -π 2 ), if x ∈ (x c -L 2 , x c + L 2 ) 0 otherwise ( 42
)
and

s ϕ (x, t) = √ A cos(2πf t + ϕ) sin 2 ( x-xc L π -π 2 ), if x ∈ (x c -L 2 , x c + L 2 ) 0 otherwise. ( 43
)
Note that one must solve each problem ((30) and ( 32)) twice: once with (42) as source term, and once with (43). To validate the method, we observe that a first order Taylor expansion of the solution with respect to the parameters gives:

u(A + dA, ϕ) = u(A, ϕ) + dA u A (A, ϕ) + O(dA 2 ), (44a)
u(A, ϕ + dϕ) = u(A, ϕ) + dϕ u ϕ (A, ϕ) + O(dϕ 2 ). (44b)
We now define u F D as the solution computed with the finite differences scheme, and u CSE as an approximation of the solution obtained as follows:

u CSE (x, t; p i + dp i ) = u(x, t; p i ) + dp i u pi (x, t; p i ),
where p i can be either A or ϕ, and u(x, t; p i ) is called reference solution. We want to measure the following quantity:

diff(T ) := u F D (x, T ) -u CSE (x, T ) L 2 (0,1) ,
and we expect, due to the Taylor expansion (44), that it will be diff(T ) O(dp 2 i ). To do that, we varied the design parameters one at a time, following the procedure illustrated in Algorithm 1. Figure 17 shows in a logarithmic scale diff A (on the left) and diff ϕ (on the right), computed as in Algorithm 1 Validation of the CSE method solve state equation with (A, ϕ + dϕ) -→ store solution in u ϕ+dϕ 9:

1: A ← A refer , ϕ ← ϕ refer 2: r ← r start r = ratio = dA A = dϕ
diff A [k] ← u A+dA (x, T ) -u refer (x, T ) -dAu A (x, T ) L 2 (0,1) 10: diff ϕ [k] ← u ϕ+dϕ (x, T ) -u refer (x, T ) -dϕu ϕ (x, T ) L 2 (0 ,1) 11:
r * = 2 12: end for line 9-10, Algorithm 1, for the linear case. The red line is the theoretical slope. Figure 18 shows the same results for the nonlinear case. One can notice that, in the nonlinear case, the slope is the one expected only starting from a certain threshold. This is due to the approach chosen in this work, "differentiate, then discretize", opposed to the other possible approach "discretize, then

✶ ✲✁ ✶ ✲✂ ✶ ✲✄ ✶ ✲☎ ✶ ✲✆ ✶ ✵ ✶ ✲✆☎ ✶ ✲✆✵ ✶ ✲✝ ✶ ✲✞ ✶ ✲✂ ✶ ✲☎ ✶ ✵ ❞✟✠✟ ✡ ☛ ☞ ☞ ✌ ✍ ✎ ❞✏✑✑ ❆ ❈ ❞✟ ✷ ✶ ✲✁ ✶ ✲✂ ✶ ✲✄ ✶ ✲☎ ✶ ✲✆ ✶ ✵ ✶ ✲✆✵ ✶ ✲✒ ✶ ✲✝ ✶ ✲✓ ✶ ✲✞ ✶ ✲✁ ✶ ✲✂ ✶ ✲✄ ✶ ✲☎ ✶ ✲✆ ✶ ✵ ❞ ❢ ✠ ❢ ✡ ☛ ☞ ☞ ✌ ✍ ✎ ❞✏✑✑ ✔ ❈ ❞❢ ✷
Figure 17: Validation of the CSE method, linear case. On the left diff A , on the right diff ϕ . In red the theoretical slope.

✶ ✲✁ ✶ ✲✂ ✶ ✲✄ ✶ ✲☎ ✶ ✲✆ ✶ ✵ ✶ ✲✆✵ ✶ ✲✝ ✶ ✲✞ ✶ ✲✟ ✶ ✲✠ ✶ ✲✁ ✶ ✲✂ ✶ ✲✄ ✶ ✲☎ ✶ ✲✆ ✶ ✵ ❞✡☛✡ ☞ ✌ ✍ ✍ ✎ ✏ ✑ ✶ ✲✁ ✶ ✲✂ ✶ ✲✄ ✶ ✲☎ ✶ ✲✆ ✶ ✵ ✶ ✲✆✵ ✶ ✲✝ ✶ ✲✞ ✶ ✲✟ ✶ ✲✠ ✶ ✲✁ ✶ ✲✂ ✶ ✲✄ ✶ ✲☎ ✶ ✲✆ ✶ ✵ ❞ ❢ ☛ ❢ ☞ ✌ ✍ ✍ ✎ ✏ ✑ Figure 18:
Validation of the CSE method, nonlinear case. On the left diff A , on the right diff ϕ .

In red the theoretical slope.

Inria differentiate": in fact Taylor expansion (44) is carried out on the exact solution. By discretizing the equations, one adds discretization error, so:

diff(T ) = O(dp 2 i ) + O(h r ),
where r is the order of the scheme used. The change of slope of the numerical results in Figure 18 corresponds to the point in which O(h r ) gets comparable with O(dp 2 i ). Note that in this case a spatial step ∆x = 0.05 has been used.

On the contrary, in the "discretize, then differentiate" approach, the Taylor expansion is consistent with the numerical solution computed, whatever grid is used. In the linear case, one can notice that both approaches yield to the same solution, which explains the results.

Numerical results

In this section we present the numerical results obtained by applying MGDA III b, to problems governed by PDEs like (29) and (31), by following the procedure described in subsection 3.1.

First of all, we need to introduce the time dependent cost functional J. The same J has been used in all the test cases, and it is the one introduced, as an example, in subsection 3.1, i.e.:

J(p) = 1 2 ∇u(p) 2 L 2 (Ω) = 1 2 (∇u(p), ∇u(p)) L 2 (Ω) . (45)
We also recall that in this case the gradient of J with respect to the parameters is:

∇J(p) = (∇u(p), ∇u pi (p)) L 2 (Ω) , (46)
whose computation requires to solve the state and the sensitivity equations. First we will present some results obtained with the instantaneous approach, that means, we define the set of objective functionals as:

J k (p) = J(u(t k ; p)) ∀k = 1, . . . , n,
later, we will present a result obtained with the windows approach, where:

J k (p) = t k+1 t k J(u(t; p))dt ∀k = 1, . . . , n.
The solving procedure is schematically recalled in the Algorithm 2, in the case of instantaneous approach. Note that the algorithm is general, the only significant difference for the windows approach is in the evaluation of the cost functionals J k (line 8).

In subsection 3.5 we introduced some specific source terms, boundary and initial conditions: unless otherwise indicated, all the results of this subsection are obtained with the same data, with fixed values of the non-optimized parameters, reported in Table 1. These choices for the values of the parameters are related to the choice of the time and grid steps: in order to have a spatial grid fine enough to represent accurately the source term we chose

∆x = L 20 ,
that is, being L = 0.2, ∆x = 0.01. For the temporal step, we recall that it has to satisfy the stability condition (33).

As we said in the first subsection, unless all the cost functionals have their minimum in the same point, there is not a unique optimal solution: this means that running the multi-objective optimization algorithm from different starting points p start we will obtain different optimal points. solve sensitivity equations with p -→ store solution in u pi 8:

J k ← J(u(t k )) 9:
compute ∇u(t k ) and ∇u(t k ) pi 10:

compute ∇ p J k as in (46) 11: solve state equation with p temp -→ store solution in u temp 17: end if 24: end while For this reason, first we find the optimal points of the single cost functionals p opt k (using, for instance, a steepest descent method), and in a second time we run MGDA from different starting points around the zone identified by the p opt k . Let us observe that, being the vector of the optimization parameters composed by an amplitude and a phase, i.e. p = (A, ϕ), the set of admissible parameters will be:

∇ p J k -→ MGDA -→ ω 12: ρ ← ρ start
J temp k ← J(u temp (t k )) 18: if J temp k < J k ∀k then
P = (0, +∞) × [0, 2π], (47)
that is an open set, therefore there is no theoretical proof of existence and uniqueness of the minimum of the single cost functional J. However, we verified numerically case by case at least the existence. Moreover, every time the parameters are updated it is controlled whether or not they are in the admissible set (see Algorithm 2, line 14): if the phase falls outside [0, 2π] it is set to a value on the boundary (for instance, if ϕ < 0 ⇒ ϕ new = 2π), while if the amplitude is negative it is set to a threshold, chosen by the user. Finally, while running the single objective optimization algorithm, we noticed an anisotropy problem: the functionals are much flatter in A then in ϕ, i.e. to have a significant difference in the value of the cost functionals, big variations of A are necessary (big if compared with the variation necessary for ϕ to obtain the same difference for the cost functionals), and this can affect significantly the numerical optimization method. To solve this problem quickly, we multiplied by a factor µ > 1 the source term (41) and, consequently, the source terms of the sensitivities. This is equivalent to a changement of variable, therefore we are now optimizing with respect to a new parameter A new , related to the old one as follows:

µ √ A new = √ A old . (48)
From (48), it is clear that, if µ is big enough, big variations of A old correspond to small variations of A new .

Instantaneous approach

We start by considering the instantaneous approach: first, it requires to chose suitable (and interesting) times for the evaluation of the cost functionals. Let us observe that, since in our case the Dirichlet boundary condition and the source term are both periodic, with the same frequency f = 4, we expect also the solution u, and therefore the cost functional J, to be periodic with frequency f = 4. For this reason it makes sense to fix all the instant t k considered inside the same period. Moreover, we would like to skip the transient part of the phenomena. This given, the choice made for the instant is the one in Table 2. Let us observe that it would be better if the instant considered t k were actual points of the temporal grid, i.e. if ∃m ∈ N : t k = m∆t ∀k. The values of the parameters corresponding to the minimum of the single objective functionals are shown in Table 3. As one can notice, all the values are very close to each other, therefore it is not an interesting case from a multi-objective optimization point of view. Anyway, in Figure 19 it is shown the final result: the red points correspond to the optimal value of the parameters for each functional; each blue point, instead, correspond to one iteration of the while loop of Algorithm 2, this means that for each blue point at least three equations have been solved (in fact in the cases in which the control of line 18 does not succeed it is necessary to solve the state equation more than once). From Figure 19 we can observe that there is a region strictly contained in the convex hull of the "red points" in which all the points are Pareto-stationary. Moreover, for smaller amplitudes (0.875 ∼ 0.885), it is possible to identify quite clearly the border of the set of Pareto-stationary points. However, being the region of the Pareto-stationary parameters so small, it is not possible to appreciate a significant decrement in the cost functional. In Figure 20 it is shown the cost functional J with respect to time, over a period, for three different choices of parameters: one outside the Pareto-stationarity zone, one inside and one on the border. As one can notice, it is almost impossible to distinguish one from the other: to do that it is necessary to zoom in, as done for Figure 21. However, it is clear that this is not a particularly interesting case, therefore from now on we will consider nonlinear equation.

Nonlinear case

For the nonlinear case (31)-(32) we chose the following values for advection and diffusion coefficients, and for µ: c = 1 b = 0.1 µ = 50.

We followed the same procedure as in the linear case, therefore at first we computed the optimal points of the single cost functionals, reported in Table 4 and drawn in red in Figure 22. In this case, the optimal points found are much more varied, compared to the linear case, therefore from a multi-objective optimization point of view it is a more interesting case. For this reason we decided to run MGDA with both the ordering criteria proposed in subsection 2.3.2 to see the differences. The results are shown in Figure 22. In order to analyse the differences, we briefly recall the ordering criteria:

• Criterion (11) selects the first gradient as:

g (1) = g k where k = argmax i=1,...,n min j=1,...,n (g i , g j ) (g i , g i ) , (49)
Inria ✸ ✸ ✁✁✁✂ ✸ ✁✁✁✄ ✸ ✁✁✁✸ ✸ ✁✁✁☎ ✸ ✁✁✁✆ ✸ ✁✁✁✝ ✸ ✁✁✁✞ ✸ ✁✁✁✟ ✸ ✁✁✁✠ ✸ ✁✁✂ ✟ ✞☎ ✟ ✞✝ ✟ ✞✟ ✟ ✟ ✟ ✟✄ ✟ ✟☎ ✟ ✟✝ ✟ ✟✟ ✟ ✠ t✡☛☞ ❏ ❖✌t☛✍✎✏ ✑✒ ✓✔✕✏ ■✕☛✍✎✏ ✑✒ ✓✔✕✏ ❖✕ t✖✏ ✗✔✘✎✏✘ ✔✙ ✑✒ ✓✔✕✏
(g i , g j ) (g i , g i ) , (50)
• The next gradients g (i) are chosen in the same way from both the criteria, that is:

g (i)
= g where = argmin j=i,...,n c j,j , First, let us observe that the shape of the Pareto-stationary set is qualitatively the same in the two cases, even if in case (22b) is much more defined. This means that, globally, the choice of the first gradient does not affect too much the final result. However, it is clear from the Figures that starting from the same point the two criteria leads to different final points. Moreover, in the result shown in Figure 22a there seems to be a tendency to converge at the closest singleobjective optimal, if there is one close: this can be explained by the fact that the criterion (49) tends to select the "middle vector" and therefore is very likely to stop the MGDA loop at the first iteration and to compute ω on the base only of one gradient; on the other hand, the criterion (50) prefers as first vector an external one and therefore is very unlikely to stop MGDA at the first iteration, this implies that in most of the cases the direction ω is computed on the base of two gradients (obviously, in this case the maximum number of gradients considered from MGDA is two, since it is the dimension of the space).

In Figure 23 it is shown the cost functional J with respect to time, for different choices of parameters: the red ones are obtained evaluating the cost functional in some parameters outside the Pareto-stationarity zone, while the green ones correspond to parameters inside the region (see the exact value in Table 5). From this Figure it is clear that the functional is periodic and that, sampling the time from only one period, we minimize the cost functional also for all the other times, as expected. In Figure 24 there is the same image, zoomed over a period. The vertical lines denote the sampling times: at a first view can seem strange that, even at some sampling times, the green functionals are greater than the red ones. However, this can happen because the notion of Pareto-stationarity is somehow local: in fact it is based on the evaluation of the gradients in that point. What we can assure is the monotonicity of all the cost functionals during the MGDA loop: for instance, considering Figure 22, we have monotonicity "along every blue line", but we cannot say anything if we skip from one to the other. In Figures 2526-27-28 it is shown the monotonicity "along the blue lines" in some particular cases: in each subfigure, the red line represents the cost functional evaluated in the starting parameters, while the purple and the green lines are the cost functional evaluated in the convergence points of MGDA obtained using, respectively, criterion (49) and criterion (50).

✥ ✁ ✥ ✂ ✄ ✄ ✄ ✄ ☎ ✆ ✄ ✆ ☎ ✆ ✁ ✆ ✂ ✥ ✥ ✄ ♣ ✝ ✞ ✟ ✠ ❆ (a) criterion (49). ✥ ✁ ✥ ✂ ✄ ✄ ✄ ✄ ☎ ✆ ✄ ✆ ☎ ✆ ✁ ✆ ✂ ✥ ✥ ✄ ♣ ✝ ✞ ✟ ✠ ❆ (b) criterion (50).
✵ ✶ ✷ ✸ ✹ ✺ ✻ ✵ ✵ ✺ ✶ ✶ ✺ ✷ ✷ ✺ ✸ ✸ ✺ ✹ ✹ ✺ t ❏ ❖✁t✂✄☎✆ ✝✞ ✟✠✡✆ ■✡✂✄☎✆ ✝✞ ✟✠✡✆
✸ ✸ ✁✂ ✸ ✄ ✸ ✄✂ ✸ ☎ ✸ ☎✂ ✄ ✄ ✂ ☎ ☎ ✂ ✸ ✸ ✂ ✹ ✹ ✂ t ❏ ❖✆t✝✞✟✠ ✡☛ ☞✌✍✠ ■✍✝✞✟✠ ✡☛ ☞✌✍✠
Inria ✸ ✸ ✁✂ ✸ ✄ ✸ ✄✂ ✸ ☎ ✸ ☎✂ ✄ ✄ ✂ ☎ ☎ ✂ ✸ ✸ ✂ ✹ ✹ ✂ t ❏ Figure 26: Monotonicity of MGDA -starting parameters: A = 0.2, ϕ = 1.6. ✸ ✸ ✁✂ ✸ ✄ ✸ ✄✂ ✸ ☎ ✸ ☎✂ ✄ ✄ ✂ ☎ ☎ ✂ ✸ ✸ ✂ ✹ ✹ ✂ t ❏ Figure 27: Monotonicity of MGDA -starting parameters: A = 1, ϕ = 2.3. RR n°8603 ✸ ✸ ✁✂ ✸ ✄ ✸ ✄✂ ✸ ☎ ✸ ☎✂ ✄ ✄ ✂ ☎ ☎ ✂ ✸ ✸ ✂ ✹ ✹ ✂ t ❏

Three sources case

In order to test MGDA in dimension higher than 2, we increased the number of parameters adding two source terms, of the same type of (41), with different centres x c . A possible result of the complete source term is shown below.

x c1 x c2 x c3 1 L x Three sources t 1 t 2 t 3
As one can notice, the width L of the three sources is the moreover, also the frequency f is the same (and for this simulation we kept the same values in the previous subsections). In this way there are 6 parameters to optimize (three amplitudes and three phases), with 7 different cost functionals. It is clear that it makes no sense to find the single-objective optimal points and then run MGDA systematically in a grid around these values, because there is no way to visualize the region. Therefore, in this case we run MGDA starting from some random points. For the first starting point we set all the amplitudes at 0.5 and all the phases at π. Figure 29 shows the different cost functionals with respect to the iterations in a logarithmic scale: one can observe that at each iteration none of the functionals increases and at least one of them decreases. In Figure 30 are compared the cost functional pre and post optimization (respectively in red and blue) over a period. In Table 6 we report the values of the parameters during the main MGDA loop: one can notice that the second source's amplitude tends to zero.

We here present the same results described above for a different starting point. One can notice that, even in this case, the second source's amplitude tends to zero.

Windows approach

The final result presented in this work is obtained with the windows approach applied to the nonlinear case. In Table 8 are reported the intervals (t k , t k+1 ) considered in this case. All the non-optimized parameters have the same values used for the nonlinear case, with the simple source, and the criterion (50) has been used. In Figure 33 are shown in red the optimal points of the single objective optimization at in blue the MGDA iterations, as in the previous cases. The shape of the Pareto-stationary zone is significantly different from the one obtained in the same case with the instantaneous approach, however, one can say that approximatively in both cases the zone is centred around the same value (A c , ϕ c ) ∼ (0.6, 2.0).

To see how the window approach and the instantaneous one work, and if there is any difference between them, we selected two starting points (in particular (A 1 , ϕ 1 ) = (1, 1.8) and (A 2 , ϕ 2 ) = (0.6, 2.3)) that are outside the Pareto-stationary zone in both cases, and we plotted J(t) with respect to time, computed in the optimal parameters obtained with the two approaches. The results are shown in Figure 34: what appears from this images is that none of the approaches is better that the other one, it depends on the final result one is interested in. Moreover, it is not fair to compare two methods whose aim is to identify a region of points, only on the base of two or three points in the identified region. The results shown in Figure 34 are meant as an example of how MGDA, coupled with different approaches, works. In Tables 9-10-11-12 the precise values are reported: J pre correspond to the cost functional pre optimization, J win correspond to the cost functional optimized with the window approach and J inst is the one optimized with the instantaneous approach. The values shown in the Tables are a further confirmation of the monotonicity of MGDA. Moreover, we recall that the approaches (windows and instantaneous) are not alternative: in fact, it is possible to consider a set of cost functionals of different nature, some of them generated by instantaneous evaluation of the cost functional, and some generated by an average operation over a window of time.

✸ ✸ ✁✂ ✸ ✄ ✸ ✄✂ ✸ ☎ ✸ ☎✂ ✄ ☎ ✸ ✹ ✂ ✻ ✼ ✽ ✾ ✄✁ ✄✄ t ❏ ✆ ✝✞✟ ✠✝t ✆ ✝✠✡t ✠✝t
Finally, we report that, in this case (i.e. with only one source term and a nonlinear state equation), the optimal parameters obtained by minimizing with a single objective technique the time-average over the whole period of the cost functional are ( Â, φ) = (0.528949, 1.96501). We observe that these values are in the middle of both the Pareto-stationarity zones, the one obtained with the instantaneous approach and the one with the windows approach. Moreover, we observe that this single objective approach is equivalent to minimizing the sum of the cost functionals used for the windows approach. In Table 13 we report the values (instantaneous and time-averaged over the intervals) of the cost functional computed in the optimal parameters ( Â, φ): as expected, with this approach there are some instant at which the cost functional increases, and also some intervals over which the average increases; this is balanced by other instant, or intervals, in which the reduction of the cost functional value is bigger than the one obtained with the multi-objective approach.

Conclusion

In this work we have presented a procedure that identifies a zone of parameters that are Paretostationary, i.e. a region of interesting values for the parameters of a complex problem. This identification has a very high computational cost that, however, can be reduced on three levels:

(i) it is possible to parallelize the whole process of identification of the Pareto-stationary zone at the highest level, since the application of MGDA from one starting point is independent from its application to another starting point;

(ii) it is also possible to parallelize the solution of the sensitivities, since they are independent one from the other, and, in the linear case, also the solution of the state; this can be useful as the number of parameters increases;

(iii) finally, it is possible to use reduced methods to solve the PDEs more efficiently.

A new ordering criterion of the gradients has been proposed in this work, to perform an incomplete Gram-Schmidt process: the numerical results have shown how this criterion leads to a better identification of the Pareto-stationary region, which is more defined.

Moreover, we would like to underline the generality of this multi-objective optimization method: in fact, even if this work was focused on systems governed by unsteady PDEs, the method can be applied to any kind of problem, it is sufficient to provide a set of cost functionals.

One can apply this method also to problems where there are different temporal phases and one wants to minimize different cost functionals for each phase, as shown in Figure 35. In fact the quantities that have to be minimized during a transient phase are usually different from the ones to minimize when the steady state has been reached. For instance, one can think to the optimization of an airfoil of an aircraft: the interesting quantities during the take-off are different from the one interesting during the steady state or during the landing. The classical approach in these cases is to minimize a convex combination of the different cost functionals, however the multi-objective optimization approach guarantees that all the cost functionals considered decrease (or, in the worst case, remain at the same value) at every iteration, which instead is not guaranteed using the classical approach.

Inria

Figure 1 :

 1 Figure 1: On the left a Pareto-stationary point, on the right a non Pareto-stationary point.

Figure 2 :

 2 Figure 2: Example of Pareto-Front

Figure 3 :

 3 Figure 3: Flow chart of the optimization algorithm. The optimization loop is highlighted in green.

Figure 5 :

 5 Figure 5: 3D examples.

Proposition 4 .

 4 If O ⊥ ∈ U , all the directional derivatives are equal along ω:

Figure 8 :

 8 Figure 8: Example to illustrate a wrong choice for a.

Figure 12 :

 12 Figure 12: Flow chart of MGDA-III b

Figure 14 :

 14 Figure 14: Convergence advection steady problem

Figure 15 :

 15 Figure 15: Convergence advection steady problem, without Neumann boundary condition

Figure 16 :

 16 Figure 16: Convergence advection unsteady problem.

ϕ 3 : 7 :

 37 solve state equation with (A, ϕ) -→ store solution in u refer 4: solve sensitivity equations with (A, ϕ) -→ store solutions in u A , u ϕ 5: for k = 0 . . . N do 6: dA ← rA, dϕ ← rϕ solve state equation with (A + dA, ϕ) -→ store solution in u A+dA8:

Inria Algorithm 2 4 :Nonlinear then 5 :

 245 Solving procedure 1: p ← p start 2: while convergence do 3: solve state equation with p -→ store solution in u if set u as advection-reaction field in sensitivity equations

4. 1 . 1

 11 Linear caseFor the linear case (29)-(30) we chose the following values for the advection and diffusion coefficients, and for µ: c = 1 b = 0.01 µ = 10.

Figure 20 :

 20 Figure 20: Cost functional for different parameters with respect to time, over a period.

Figure 21 :

 21 Figure 21: Cost functional for different parameters with respect to time -zoom.

Figure 22 :

 22 Figure 22: Result of MGDA applied to the nonlinear case from different starting points.Inria

Figure 23 :

 23 Figure 23: Cost functional with respect to time for different parameters values -nonlinear case.

Figure 24 :Figure 25 :

 2425 Figure 24: Cost functional with respect to time for different parameters values over a periodnonlinear case.

Figure 28 :

 28 Figure 28: Monotonicity of MGDA -starting parameters: A = 1.1, ϕ = 1.6.

Inria

Figure 29 :

 29 Figure 29: Cost functionals with respect to iteration -three sources -symmetric starting point.

Figure 30 :

 30 Figure 30: Cost functionals with respect to time -three sources -symmetric starting point.

Inria

Figure 31 :

 31 Figure 31: Cost functionals with respect to iteration -three sources -asymmetric starting point.

Figure 32 :Figure 33 :

 3233 Figure 32: Cost functionals with respect to time -three sources -asymmetric starting point.

( 3 .

 3 00, 3.04) (3.04, 3.08) (3.08, 3.12) (3.12, 3.16) (3.16, 3.20)(3.20, 3.25)

  (a) Starting parameters (A, ϕ) = (1, 1.8).

  Figure 34: Comparing the windows and the instantaneous approach.

J S-O I 1 =

 1 (3.00, 3.04) 0.0736 I 2 = (3.04, 3.08) 0.0568 I 3 = (3.08, 3.12) 0.0631 I 4 = (3.12, 3.16) 0.0744 I 5 = (3.16, 3.20) 0.0669 I 6 = (3.20, 3.25) 0.1076

Figure 35 :

 35 Figure 35: Example of different cost functionals for different phases.

Table 1 :

 1 Values of non-optimized parameters

	Parameter	Symbol Value
	frequency	f	4
	source centre	x c	0.2
	source width	L	0.2
	b.c. average	k	5
	b.c amplitude	α	1

Table 2 :

 2 Instant considered

	t 1	t 2	t 3	t 4	t 5	t 6	t 7
	3.00 3.04 3.08 3.12 3.16 3.20 3.24

Table 3

 3

			: Single-objective optimization results -linear case
	time t k [s]	3.00	3.04	3.08	3.12	3.16	3.20	3.24
	A	0.907108 0.882811 0.876467 0.906197 0.887208 0.873657 0.904139
	ϕ	5.12258	5.11216	5.12904	5.1204	5.1116	5.12704	5.12732
								Inria

Table 4 :

 4 Single-objective optimization results -nonlinear case

	time t k [s]	3.00	3.04	3.08	3.12	3.16	3.20	3.24
	A	0.34512 1.00652 0.644476 0.47365 1.15389 0.663846 0.216462
	ϕ	2.36373 2.04267 1.76051 2.29071 2.1838	1.65968	2.22826

Table 5 :

 5 Values of the parameters chosen for Figure23-24

	A	0.15	0.2	1.0	1.1	0.45	0.4	0.6
	ϕ	2.35	1.6	2.3	1.6	1.85	2.15	2.0
	collocation outside outside outside outside inside inside inside

Table 6 :

 6 Parameters during MGDA main loop -three sources -symmetric starting point.

	A 1	A 2	A 3	ϕ 1	ϕ 2	ϕ 3
	0.5	0.5	0.5	3.14159 3.14159 3.14159
	0.326183	0.01	0.497761 2.54601	3.031	3.55003
	0.422592	0.0155267	0.276467	2.285	3.02927 3.51984
	0.43974	0.0033347	0.248952 2.26425 3.02819 3.51827
	0.434359 0.00126382	0.242303 2.26661 3.02778 3.52063
	0.366366	0.01	0.126234 2.29823 3.02349 3.55691
	0.348423	0.01	0.01	2.13807 3.01764 3.59255
	0.349058 0.00062438 0.0102434	2.137	3.01759 3.59259
	0.349064 0.000237505 0.0102312 2.13699 3.01759 3.59259
	0.349065 9.08649e-05 0.0102276 2.13698 3.01759 3.59259
	0.349065	3.4122e-05	0.0102267 2.13698 3.01759 3.59259
	0.349065 1.16204e-05 0.0102264 2.13698 3.01759 3.59259

Table 7 :

 7 Parameters during MGDA main loop -three sources -asymmetric starting point. .000162565 0.024526 2.17177 2.9569 3.75495 0.370198 2.57116e-05 0.0245225 2.17177 2.9569 3.75495 0.370198 5.76837e-06 0.0245223 2.17177 2.9569 3.75495

	A 1	A 2	A 3	ϕ 1	ϕ 2	ϕ 3
	0.25	0.5	0.75	3.14159 3.14159 3.14159
	0.209686	0.01	0.671903 2.70767 2.9716 3.70705
	0.28742	0.01	0.516188 2.56182 2.96944 3.7173
	0.32335	0.01	0.385196 2.46421 2.96649 3.7297
	0.352622	0.0151327	0.325739 2.40728 2.96619 3.72979
	0.416572	0.01	0.209501 2.29981 2.96446 3.73147
	0.432156 0.00440427	0.179948	2.2774 2.96396 3.73208
	0.367466	0.01	0.0574286 2.25788 2.9599 3.74702
	0.369833	0.01	0.0247545 2.1731 2.95696 3.75488
	0.37018	0.00183642 0.0246056 2.17192 2.9569 3.75495
	0.370197 0

Table 8 :

 8 Intervals considered -window approach.

	I 1	I 2	I 3	I 4	I 5	I 6

Table 9 :

 9 Comparing the windows and the instantaneous approach -starting parameters (A, ϕ) = (1, 1.8) -instantaneous values.

	J pre	J win	J inst
	t 1 = 3.00 3.4651 3.0819 2.7550
	t 2 = 3.04 1.4981 1.3320 1.4196
	t 3 = 3.08 1.4826 1.4694 1.3487
	t 4 = 3.12 2.6260 2.2950 2.1037
	t 5 = 3.16 2.0030 1.7714 1.8955
	t 6 = 3.20 1.6249 1.6738 1.5410
	t 7 = 3.24 3.5074 3.1961 2.7891
			Inria

Table 10 :

 10 Comparing the windows and the instantaneous approach -starting parameters (A, ϕ) = (1, 1.8) -averaged values. = (3.00, 3.04) 0.1015 0.0883 0.0847 I 2 = (3.04, 3.08) 0.0509 0.0495 0.0505 I 3 = (3.08, 3.12) 0.0838 0.0770 0.0687 I 4 = (3.12, 3.16) 0.1006 0.0871 0.0861 I 5 = (3.16, 3.20) 0.0618 0.0597 0.0620 I 6 = (3.20, 3.25) 0.1415 0.1339 0.1167

	J pre	J win	J inst
	I 1

Table 11 :

 11 Comparing the windows and the instantaneous approach -starting parameters (A, ϕ) = (0.6, 2.3) -instantaneous values. = 3.00 2.1546 2.2237 2.1386 t 2 = 3.04 1.7273 1.5738 1.6303 t 3 = 3.08 2.2081 2.0684 1.9736 t 4 = 3.12 1.6072 1.6375 1.5941 t 5 = 3.16 1.7478 1.6514 1.7390 t 6 = 3.20 2.5253 2.3834 2.2855 t 7 = 3.24 2.4536 2.5198 2.3795

	J pre	J win	J inst
	t 1

Table 12 :

 12 Comparing the windows and the instantaneous approach -starting parameters (A, ϕ) = (0.6, 2.3) -averaged values.

	J pre	J win	J inst
	I 1 = (3.00, 3.04) 0.0702 0.0691 0.0696
	I 2 = (3.04, 3.08) 0.0833 0.0759 0.0753
	I 3 = (3.08, 3.12) 0.0763 0.0747 0.0709
	I 4 = (3.12, 3.16) 0.0641 0.0636 0.0652
	I 5 = (3.16, 3.20) 0.0850 0.0789 0.0795
	I 6 = (3.20, 3.25) 0.1293 0.1288 0.1213

Table 13 :

 13 Results of the single-objective optimization.

	J S-O
	t 1 = 3.00 2.2961
	t 2 = 3.04 1.4410
	t 3 = 3.08 1.4575
	t 4 = 3.12 1.7497
	t 5 = 3.16 1.8121
	t 6 = 3.20 1.7042
	t 7 = 3.24 2.3827

RR n°8603

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399