BI-OBJECTIVE OPTIMAL CONTROL OF SOME PDES: NASH EQUILIBRIA AND QUASI-EQUILIBRIA∗

This paper deals with the solution of some multi-objective optimal control problems for several PDEs: linear and semilinear elliptic equations and stationary Navier-Stokes systems. Specifically, we look for Nash equilibria associated with standard cost functionals. For linear and semilinear elliptic equations, we prove the existence of equilibria and we deduce related optimality systems. For stationary Navier-Stokes equations, we prove the existence of Nash quasi-equilibria, i.e. solutions to the optimality system. In all cases, we present some iterative algorithms and, in some of them, we establish convergence results. For the existence and characterization of Nash quasi-equilibria in the Navier-Stokes case, we use the formalism of Dubovitskii and Milyutin. In this context, we also present a finite element approximation and we illustrate the techniques with numerical experiments. Mathematics Subject Classification. 35Q30, 76D05, 76N10. Received April 24, 2020. Accepted May 5, 2021.


Introduction
We consider bi-objective optimal control problems for various PDEs and systems. First, an introductory problem corresponding to a linear elliptic PDE is analyzed with detail. Then, we deal with a similar semilinear elliptic PDE. Finally, we deal with the stationary Navier-Stokes system, that is, the equations satisfied by the time-independent velocity field and pressure of an incompressible viscous fluid flow.
Our aims are to prove existence, characterize efficiently the equilibria and, also, compute numerical solutions to these multi-objective control problems. They are very important from the mathematical viewpoint and appear frequently in the applications; for some previous works on the subject, see for instance [2].
In classical control theory, we usually find a state equation or system and one control, with the task of achieving a predetermined goal. Frequently (but not always), the goal is to minimize a cost functional within a prescribed family of admissible controls; see for instance [10,16]. A different and interesting situation arises when several (in general, conflictive or contradictory) objectives are considered. This may happen, for example, if the cost function is the sum of several terms and it is not clear that an average provides a reasonable criterion. Also, it can be expectable to have more than one control acting on the equation. In these cases, we are led to consider multi-objective control problems. In contrast with the mono-objective case, various strategies for the choice of good or the best controls can appear, depending on many circumstances. Moreover, these strategies can be cooperative or noncooperative (depending on whether or not several controls mutually cooperate in order to achieve prescribed goals).
There exist several concepts of equilibrium for multi-objective problems, with origin in game theory. Each of them can be identified with a strategy. Thus, let us mention the non-cooperative optimization strategy proposed by Nash [18], the Pareto cooperative strategy [19] and the Stackelberg hierarchical-cooperative strategy [24]. In the context of the control of PDEs, up to date, there have been some works on the subject like the seminal papers by Lions [15,17] and other more recent contributions, like [6,21,22].
In this paper, we will be concerned with Nash equilibria and quasi-equilibria associated with standard cost functionals. To be more precise, let us give some details in the case of a linear elliptic equation.
Thus, let the fluid domain be a bounded open set Ω ⊂ R N , with N ≥ 1. Let us introduce four nonempty open subsets, O 1 , O 2 , ω 1 and ω 2 and let us assume that a function u id defined on O i is given for i = 1, 2. We would like to find a couple (f 1 , f 2 ) ∈ L 2 (ω 1 ) × L 2 (ω 2 ) (the control pair) with the following property: the associated state u, that is, the solution to the Dirirchlet problem −ν∆u = f 1 1 ω1 + f 2 1 ω2 in Ω, u = 0 on ∂Ω, (1.1) is such that (f 1 , f 2 ) is a Nash equilibrium for the functionals where a, µ > 0 (see Sect. 2.1). Our first main goal will be to find conditions under which at least one Nash equilibrium exists. The second one will be to characterize these equilibria in terms of first order optimality conditions, i.e. to derive a system that the optimal solution, together with some associated adjoint states must satisfy. The third one will be to show how Nash equilibria can be computed and present related convergent algorithms. These tasks will be successfully accomplished below. Then, we will deal with some more complex equations and systems.
In particular, in what concerns optimality conditions, we will try to extend to this context the techniques usually employed for distributed control problems (see for instance [1,10,16]).
At some point, we will apply the so called formalism of Dubovitskii and Milyutin. This approach was introduced in the context of mathematical programming and has been succesfully applied to the solution of many optimal control problems for ODEs since the 70's. A good presentation of its applications to these areas can be found in Girsanov [11]; see also Flett [9]. Later, these techniques have been applied successfully to some distributed control problems; see [3,4]. In the present context this method is appropriate; note that this would not be so clear if state constraints were imposed.
The plan of this paper is the following. In Section 2, we consider the relatively simple problem given by (1.1)-(1.2). We prove the existence of Nash equilibria, we provide an optimality system and we present some iterative algorithms.
In Section 3, a more complicated problem is analyzed: a semilinear elliptic PDE together with functionals of the same kind. Here, in view of the nonlinearity, we must work with Nash equilibria and Nash quasi-equilibria, that is solutions to the optimality system. Again, existence, characterization and computation-oriented results are established. Section 4 deals with the stationary Navier-Stokes system. New difficulties are found: nonlinearity, lack of uniqueness, lack of regularity of the functionals, etc. We also discuss the existence and optimality characterization of Nash equilibria and quasi-equilibria. Additionally, we present some iterative algorithms and the results of several numerical experiments.
Finally, some additional comments and open questions are presented in Section 5.

Introductory problem: a linear elliptic PDE
In the sequel, we denote by · and (· , ·) the usual L 2 norm and inner product, respectively. The symbol 1 D will be used to denote the characteristic function of the set D. Also, C will stand for a generic positive constant. For simplicity, we will assume that only two controls act on the system and two functionals are taken into account but very similar considerations hold for systems with a larger number of controls and functionals.

Definition of Nash equilibria
Let Ω ⊂ R N be a nonempty bounded connected open set with regular boundary ∂Ω and let us assume that ω 1 and ω 2 are nonempty disjoint open subsets of Ω.
We will consider the problems where the f i ∈ L 2 (ω i ) are the controls and u is the state.
Let O 1 an O 2 be open sets, representing prescribed observation domains and let the J i be given by where the u id ∈ L 2 (O i ) are given functions and a and µ are positive constants. The first bi-objective control problem considered in this paper is the following: Find a Nash equilibria associated with (2.1) and (2.2), that is, a pair ( In this case, since the control-to-state mapping is well-defined, linear and continuous and the cost functionals J i are quadratic, strictly convex and C 1 , it is immediate to deduce that (f 1 ,f 2 ) is a Nash equilibria if and only if

Existence and characterization of Nash equilibria
For future purposes, note that where ϕ i is the i-th adjoint state associated with (f 1 , f 2 ), i.e. the solution to where u is the state associated with (f 1 , f 2 ). We can now present the main result of the section. It deals with the existence and characterization of Nash equilibria.
3. Suppose now that O 1 = O 2 , Then Λ 0 is self-adjoint and positive, that is, Therefore, for any a, µ > 0, (2.7) is uniquely solvable and the solution is moreover the unique minimizer of the functional K : L 2 (ω 1 ) × L 2 (ω 2 ) → R, given by Thus, the proof of part 3 is also done.
Note that, when O 1 = O 2 , the operator Λ 0 is not necessarily self-adjoint. As a consequence, for large a/µ there can be no solution to (2.7). However, since Λ 0 is compact, we can make use of Fredholm's Alternative and affirm that there exist at most β 1 , β 2 , . . . with β n → +∞ such that, if a/µ = β n for all n ≥ 1, a unique Nash equilibrium exists.

Algorithms and convergence
We present in this section several algorithms that can be used for the computation of Nash equilibria. The first one relies on a fixed-point reformulation of (2.6): The following convergence result holds: Theorem 2.2. There exists χ 0 ∈ (0, χ], such that, if a/µ ≤ χ 0 , the controls provided by ALG 1 satisfy , where (f 1 ,f 2 ) is the unique Nash equilibrium associated with J 1 and J 2 . Furthermore, the speed of convergence is at least linear.
The proof is immediate: it suffices to rewrite (2.6) in the form (2.7) and note that ALG 1 is the usual fixed-point iteration method for the mapping Λ : The following two algorithms are inspired by the classical optimal step gradient and conjugate gradient methods.

ALG 2:
(a) Choose f 0 i ∈ L 2 (ω i ), i = 1, 2. (b) Then, for given n ≥ 0 and f n i ∈ L 2 (ω i ), compute the solution u n to (2.9) and the solution ϕ n i to (2.10) and take where g n i = a ϕ n i | ωi + µf n i (2.13) and (2.14) ALG 3: For n = 0, perform one step of ALG 2 and take d 0 Then, for given n ≥ 1 and f n i ∈ L 2 (ω i ), compute the solution u n to to (2.9) and the solution ϕ n i to (2.10) and take Note that, in these iterates, the functions ρ → J 1 (f n 1 − ρg n 1 , f n 2 ), ρ → J 2 (f n 1 , f n 2 − ρg n 2 ), etc. are quadratic and strictly convex. Consequently, the ρ n i can be computed explicitly. To the best of our knowledge, the convergence of ALG 2 and ALG 3 cannot be guaranteed. Note however that we advance in the directions indicated by the partial derivatives of the J i and, consequently, it is reasonable to expect that they do converge. Now, let us assume as before that O 1 = O 2 = O. Then, the Nash equilibrium is the unique minimizer of the functional K in (2.8) and it makes sense to consider the following algorithms:   and take and Step Conjugate Gradient Method for K.
(a) Solve the systems: Then, for given n ≥ 1 and f n i ∈ L 2 (ω i ), compute the solution u n to to (2.9) and the solution ψ n to (2.18) and take Again, the functions ρ → K(f 1 n − ρd 1 n , f 2 n − ρd 2 n ) are quadratic and strictly convex, whence the ρ n can be easily found. The following result holds: is the unique Nash equilibrium associated with J 1 and J 2 .
The proof is immediate in view of the properties of the linear operator Λ 0 . It suffices to directly apply well known classical results; see for instance [7,20]. It is remarkable that, in this particular case, ALG 2 and ALG 2' and ALG 3 and ALG 3' only differ in the definitions of the step parameters ρ n i and ρ n .

Nash equilibria and quasi-equilibria
This section is devoted to introduce Nash optima in the semilinear case. Now, we must distinguish equilibria from quasi-equilibria and take into account the particularities of each of them.
Let us assume that The state equation reads as follows: It is then well known that, for each ( , there exists exactly one (strong) solution u to (3.2). As in Section 2, we will consider the cost functionals From standard optimal control theory, it is known that, under the previous assumptions on φ, the cost functionals J i are again C 1 and satisfy where ϕ i is the unique solution to (that is, ϕ i is the i-th adjoint state corresponding to f 1 and f 2 ) and u is the solution to (3.2); see for instance [10,16].
Definition 3.1. It will be said that (f 1 ,f 2 ) is a Nash quasi-equilibrium for (3.2) and (2.2) if (f 1 ,f 2 ) satisfies (2.4), that is, (f 1 ,f 2 ) solves, together withû and theφ i , the optimality system is also a Nash quasi-equilibrium. However, the converse is not necessarily true.
Proof. We have to prove that, under the previous conditions, any quasi-equilibrium satisfies (2.3). Thus, let us assume that (2.4) holds and let the functionalsJ i be given as follows: Let us first prove that there exists η 0 > 0 and C 0 > 0, only depending on Ω, the u id and M , such that, if a/µ ≤ η 0 , any quasi-equilibrium (f 1 ,f 2 ) satisfies Indeed, if we take the PDEs for u and ϕ i , we respectively multiply by u and ϕ i and we integrate in Ω, we get Now, from the properties of φ and Hölder and Young inequalities, we see that Consequently, if a/µ is sufficiently small one has (3.5) for some C 0 . Observe that, in view of the properties of φ, theJ i are twice continuosly differentiable. Let us see that, if a/µ is sufficiently small,J i (f i ; g i , g i ) > 0 for all f i ∈ L 2 (ω i ), all nonzero g i ∈ L 2 (ω i ) and i = 1, 2.
For instance, let us take i = 1. Then, where ϕ 1 is, together with u, the solution to For any small ε > 0 and any g 1 ∈ L 2 (ω 1 ) N , let us introduce u ε and ϕ ε 1 , with (3.7b) Let us put Note that these limits exist in H 1 0 (Ω). This is easy to see by substracting (3.7a) from (3.7b), dividing by ε and letting ε → 0 (here, we must use that φ ∈ C 2 (R) and φ and φ are uniformly bounded). Furthermore, one has Observe that, from classical elliptic regularity results, one has for some constants C 1 = C 1 (Ω, u 1d , u 2d ) and C 2 = C 2 (Ω). On the other hand, from the PDE satisfied by ψ 1 and z, one also has where C 3 = C 3 (Ω, M ). Therefore, taking into account that, for N ≤ 8, H 2 (Ω) → L N/2 (Ω) with continuous embedding, we see from (3.10) and (3.11) that Clearly, this proves that, if N ≤ 8 and a/µ is sufficiently small,J 1 (f 1 , g 1 , g 1 ) > 0 for all f 1 ∈ L 2 (ω 1 ) and all nonzero g 1 = 0 and, consequently,J i is strictly convex and possesses a unique global minimum atf i . A similar fact holds forJ 2 . That is, (f 1 ,f 2 ) is a Nash equilibrium for (3.2) and (2.2).

Existence of Nash equilibria and quasi-equilibria
We can now prove the existence of Nash equilibria for (3.2) and (2.2): Theorem 3.3. Let the assumptions in Theorem 3.2 be satisfied. The following assertions hold: 1. There exists η 0 ≤ η, only depending on u 1d , u 2d and M such that, if a/µ ≤ η 0 , there exists at least one Nash equilibrium (f 1 ,f 2 ) for (3.2) and (2.2). 2. There exists η 1 ≤ η 0 such that, whenever a/µ ≤ η 1 , the Nash equilibrium is unique.
Proof. To prove the existence of a Nash equilibrium, it suffices to check that, if a/µ is sufficiently small, the optimality system (3.4) possesses at least one solution. To this purpose, we can use Schauder's Fixed-Point Theorem.
Thus, let us consider the mapping Ψ : where ϕ i is the solution to (3.3) and u is the state associated withf 1 andf 2 . It is not difficult to see that Ψ is well-defined, continuous and compact. Furthermore, if a/µ is small enough, Ψ maps the whole space This can be seen from the following estimates: (a) First, from (3.2) and the properties satisfied by φ, one has From (3.12a) and (3.12b), our assertion follows. Hence, we can apply Schauder's Theorem to Ψ and the existence of a Nash equilibrium is ensured. This ends the proof of existence.
In view of what can be said in the linear case, it is reasonable to expect that, when O 1 = O 2 , there exist quasi-equilibria for any a > 0 and any µ > 0. However, to our knowledge, this is unknown.

Algorithms and convergence
In this section, we present some iterative algorithms, similar to those considered in the linear case.
Then, for given n ≥ 0 and (f n 1 , f n 2 ) ∈ L 2 (ω 1 ) × L 2 (ω 2 ), compute the solution u n to −∆u n + φ(u n ) = f n 1 1 ω1 + f n 2 1 ω2 in Ω, u n = 0 on ∂Ω, (3.17) the solution ϕ n i to the system and, finally, take Theorem 3.4. Let the assumptions in Theorem 3.3 be satisfied. There exists η 2 ≤ η 1 such that, if a/µ ≤ η 2 , the couples provided by ALG 4 satisfy (f n This proof is easy. Indeed, note that, under the assumptions in Theorem 3. ALG 5: Then, for given n ≥ 0 and (f n 1 , f n 2 ) ∈ L 2 (ω 1 ) × L 2 (ω 2 ), compute the solution u n to (3.17) and the solution ϕ n i to (3.18) and take (c) Then, for given n ≥ 1, (f n 1 , f n 2 ) ∈ L 2 (ω 1 ) × L 2 (ω 2 ), g n−1 i ∈ L 2 (ω i ) and d n−1 i ∈ L 2 (ω i ), compute the solution u n to (3.17) and the solution ϕ n i to (3.18) and take and Note that in ALG 3 and ALG 6, the coefficient γ n i is given by different expressions. The reason is that, now, the system is nonlinear and, as usual, it seems better to impose Polak condition to ensure convergence.
In the iterates in ALG 5 and ALG 6, we have to compute the steps ρ n 1 and ρ n 2 . Now, the functions to minimize are not quadratic and these numbers cannot be computed exactly. Thus, in practice, we must apply a one-dimensional minimization method to these purposes. More details are given below.

The stationary Navier-Stokes system
This section is devoted to analyze Nash equilibria and quasi-equilibria for the stationary Navier-Stokes equations. In view of the properties of the state system and, in particular, nonlinearity and possible lack of uniqueness, this task will be much more complicated than in Sections 2 and 3. In fact, we will only be able to prove the existence of quasi-equilibria; that Nash equilibria exists is an open question.
The stationary Navier-Stokes equations are the following: in Ω, ∇ · u = 0 in Ω, u = 0 on ∂Ω. (4.1) The state variables are u = (u 1 , . . . , u N ) and p. They can be interpreted as the velocity field and the pressure of a steady viscous Newtonian fluid flow. The controls are f 1 1 ω1 and f 2 1 ω2 and can viewed as external fields of forces respectively applied at the points in ω 1 and ω 2 .

Nash equilibria and quasi-equilibria
The positive constant ν is the kinematic viscosity of the fluid. It must be regarded as a measure of the fluid internal friction.
As in the previous sections, let us introduce the functionals J i with where the u id ∈ L 2 (O i ) N and a, µ > 0. Note that, here, the J i depend not only on the controls f i but also on the associated state u. This is due to the possible non-uniqueness of solution to (4.1), that can be true when ν is not sufficiently large with respect to Ω and the f i L 2 (ωi) . Definition 4.1. It will be said that (f 1 ,f 2 ) ∈ L 2 (ω 1 ) N × L 2 (ω 2 ) N is a Nash equilibrium for (4.1) and (4.2) if there exists an associated state (û,p) satisfying Definition 4.2. It will be said that (f 1 ,f 2 ) ∈ L 2 (ω 1 ) N × L 2 (ω 2 ) N is a Nash quasi-equilibrium for (4.1) and (4.2) if there exists a solution (û,p,φ 1 ,q 1 ,φ 2 ,q 2 ) to the following coupled system

Existence of Nash quasi-equilibria
Let us recall some classical spaces, usual for the analysis of the Navier Stokes equations: They are closed subspaces of L 2 (Ω) N and H 1 0 (Ω) N , respectively; accordingly, they are Hilbert spaces for the inner products (· , ·) and (· , ·) H 1 0 . Also, we have the compact embeddings V → H ≡ H → V where X denotes the dual space of X.
In the sequel, we will have to use the Stokes operator where P : L 2 (Ω) N → H is the usual orthogonal projector. It is known that A can be uniquely extended to a bounded linear operator in L(V ; V ), again denoted by A. Then, A is self-adjoint and one has Let us consider the trilinear continuous forms b(· , · , ·) andb(· , · , ·), with Note that there exist bilinear continuous mappings B andB with The existence of a quasi-equilibrium for (4.1) and (4.2) is guaranteed by the following result: There exists β > 0, only depending on Ω, ν, the ω i , the O i and the u id , such that, if a/µ ≤ β, the optimality system (4.4) possesses at least one solution.
Proof. The proof is not difficult. Thus, let us rewrite (4.4) in the form where Z : V 3 → V 3 is the following continuous and compact mapping: Now, we want to find a fixed-point of Z. It is easy to see that, if a/µ is sufficiently small, the fixed-points of Z are uniformly bounded. Consequently, form the Leray-Schauder Principle, we deduce that, for small a/µ, (4.4) is solvable.
In this case, we cannot ensure the existence of Nash equilibria. This is due to the lack of uniqueness of the control-to-state mapping and the lack of convexity of the related functionals. However, the following holds: Theorem 4.4. Let (f 1 ,f 2 ) be a Nash equilibrium for (4.1) and (4.2). Then, (f 1 ,f 2 ) is a Nash quasi-equilibrium.
We will give a proof of this result relying on the Dubovitsky-Milyutin formalism (see [11]). To this purpose, we must recall some technical results.
Lemma 4.5. Let K 1 , . . . , K n be convex cones in a Banach space X with apex at 0. For each i, we assume that either K i is open or it is a closed subspace. Then the following conditions are equivalent: • There exist linear functionals f i ∈ K * i with i = 1, . . . , n, not all zero, such that f 1 + f 2 + · · · + f n = 0. Here, for any i, we have denoted by K * i the dual cone to K i , that is, K * i := {f ∈ X : f (e) ≥ 0 ∀e ∈ K i }. For the proof, see for instance Lemma 5.11 of [11].
Note that the adjoint system can be equivalently rewritten in the form where Dϕ i = 1 2 (∇ϕ i + (∇ϕ i ) t ) and the pressure has been redefined. This observation will be used in the following lemma and the following section. Lemma 4.6. Let (f 1 , f 2 ) ∈ L 2 (ω 1 ) N × L 2 (ω 2 ) N be given and let u ∈ H 2 (Ω) N ∩ V be (together with p) an associated solution to (4.1). Let R : V → V be the linear mapping defined by Rϕ := (νA) −1 (−B(u, ϕ)), where A is the (extended) Stokes operator in V . Then, for any non-empty open set ω ⊂ Ω, is a norm in N (Id + R).
Proof. From the usual alliptic estimates for the solutions to stationary Navier-Stokes systems, one has u ∈ H 2 (Ω) N , ∇u ∈ L 6 (Ω) N ×N and u ∈ L ∞ (Ω) N , see for instance [23]. Consequently, R is well defined and compact.
Then ϕ n is uniformly bounded in V .
Proof. First, note that, since R is a compact operator, dim(N (Id + R)) < +∞ and R(Id + R) is closed, in view of Fredholm's Alternative Theorem (see for instance [5]). Now, letφ n be, for each n, the unique function in N (Id + R) satisfying Then, ψ n = (ϕ n −φ n ) + R(ϕ n −φ n ). Also, ϕ n −φ n H 1 0 is bounded by a constant C. Indeed, if this were not the case, we could assume that ϕ n −φ n V = dist(ϕ n , N (Id + R)) → +∞.

Let us introduce
Then ζ n + Rζ n = ψ n ϕ n −φ n H 1 0 . (4.6) Since the ζ n V = 1 and ψ n → ψ in V , we see from (4.6) that ζ n → ζ for some ζ ∈ N (Id + R). So, for all n. But this is an absurd, since ζ n → ζ in V . Finally, let us deduce that ϕ n H 1 0 ≤ C. We can write that ϕ n =φ n + η n , with ] ω ( note that theφ n belong to a finite dimensional space) and This ends the proof. Now, we can prove the main result in this section.
Thus, let (f 1 ,f 2 ) be a Nash equilibrium for (4.1) and (4.2). Let us introduce the sets By assumption, there exists a strong solution (û,p) to (4.1) such that (f 1 ,û) and (f 2 ,û) respectively solve the following extremal problems: and Let us introduce some cones associated with (4.7a) and (4.7b). To this end, we set It is then clear that M 1 and M 2 are continuosly differentiable mappings, with We also have Then the cones of descent directions associated with (4.7a) and (4.7b) are the sets and the spaces of tangent directions are the û)).

Let us prove that
Let us introduce the linear mapping S = νA(Id + R), where R is the linear operator in Lemma 4.6. From Fredholm's Alternative, we have that R(S) is closed in V and, from Lemma 4.7, we find that the ϕ n are uniformly bounded in V . As an inmediate consequence, we see that, at least for a subsequence, ϕ n → ϕ weakly in V , with − ϕ| ωi = h and Sϕ = w. In other words, (h, w) ∈ R(M i (f i ,û) * ).
At this moment, we can apply the Dubovitskiy-Milyutin formalism to (4.7a) and (4.7b) and deduce that the descent cone D i and the tangent space T i are disjoint for i = 1, 2 (see [11]): Taking into account the definitions of D i and T i and the fact that R(M i (f i ,û) * ) is closed, we see at once that Hence, there exist ϕ 1 , ϕ 2 ∈ V such that and (4.8) can be rewritten in the form Taking w i = 0, we find that h i = µf i . On the other hand, taking h i = 0 and recalling (4.9), we see that ϕ i | ωi = −µf i for each i = 1, 2 and Therefore, ϕ i solves, together with some q i ∈ L 2 (Ω), the following linear system for each i = 1, 2 and, furthermore,f From (4.1), (4.11a) and (4.11b), we deduce that (f 1 ,f 2 ) is a Nash quasi-equilibrium and the proof is done.

Algorithms
This section is dedicated to the presentation of four iterative algorithms for the computation of a Nash quasi-equilibrium for (4.1) and (4.2). Three of them are similar to those in Sections 2.3 and 3.3; the fourth one is of the Newton kind. for i = 1, 2.
Before presenting the results of some simulations, let us describe another algorithm. It is based on Newton's iterates and, like ALG 7, aims to compute directly a solution to the optimality system (4.4). In practice, this method is much faster but, as is usual for Newton methods and variants, needs a nontrivial starting process (see below).
(a) Choose (f 0 1 , f 0 2 ) ∈ L 2 (ω 1 ) N × L 2 (ω 2 ) N and ν 0 ∈ R + and compute a solution (u 0 , p 0 ) to the Stokes problem  Note that ALG 7 and ALG 10 are designed to compute a solution to the optimality system (4.4) (a quasiequilibrium) that, possibly, is not be a Nash equilibrium. Thus, from the viewpoint of the Calculus of Variations, ALG 7 and ALG 10 can be seen as "indirect methods". Contrarily, ALG 8 and ALG 9 can be viewed as realizations of the "direct method" of Calculus of Variations.

Numerical experiments
In order to illustrate the behavior of the previous algorithms, we will present in this section the results of some numerical experiments. The computations have been performed with the FreeFem++ package (see [14]).
In our experiments, we try to solve numerically the optimality system (4.4), where a is chosen with 0 < a < 2 and µ = 2 − a. In the tests, several values of a and the Reynolds number Re := U L/ν have been chosen. Here, the characteristic speed and length are respectively given by U = max( u 1d , u 2d ) and L = max x,x ∈Ω |x − x |.

Test 1
In this first test, our domain is composed of two rectangles O 1 and O 2 and we assume that the controls act on the narrow bands ω 1 and ω 2 . In order to solve numerically the systems (4.12), (4.13), (4.21), (4.22) and (4.23), we have used meshes like the one in Figure 1 and a mixed finite element formulation with continuous piecewise P 1 -bubble and P 1 functions respectively for the velocity field and the pressure; for details, see [12,13].
The data u id are the following: u 1d = ∇ × ψ 1d , where ψ 1d is the solution to the problem and u 2d ≡ 0. That means that the "desired" (ideal) configurations correspond to a uniformly rotating flow in O 1 and a fluid at rest in O 2 (see Fig. 2). For ALG 8 and ALG 9, the computation of the steps ρ n i has been performed by applying a standard dichotomy process. This provides acceptable values after a reduced amount of iterates. In the case of ALG 7, ALG 8 and ALG 9, the stopping test has been with ε = 10 −6 Re. This choice establishes a criterion with uniform severity. This has also been the stopping criterion for the external iterates in ALG 10. For the internal loops (indexed by k), the stopping test has been u n,k+1 − u n,k L ∞ + ϕ n,k+1 − ϕ n,k L ∞ ≤ ε.  In all cases, the starting control pair has been (f 0 1 , f 0 2 ) = (0, 0). On the other hand, for all choices of a and Re, all the algorithms provide the same solution.
Some results are depicted in Figure 3. There, we can see the computed velocity fields corresponding to some Nash quasi-equilibria, for a = 1.2 and various Re. As expected, the appearance of the velocity field in O 1 (resp. O 2 ) reminds us of u 1d (resp. u 2d ).
We also present in Figure 4 the velocity field computed with ALG 10, a finer mesh and a higher Re.
In order to compare the behavior of the Gradient-like, Conjugate Gradient-like and Newton methods, we present some related numerical values in the Tables in Figure 5. Specifically, we have indicated there the numbers of iterates needed in each case to fulfill the stopping test and produce an approximation. In the case of ALG 10, we have included the total numbers of iterates (that is, the iterates needed for initialization have been taken into account).
Also, we have included in Figure 6 a Table with a comparison of the computation times and required numbers of iterates for each method.
The convergence of the numerical aproximations toward a quasi-equilibrium has been checked in the following way. We have considered several meshes corresponding to triangles of typical sizes 0.2, 0.13, 0.1, . . . , as indicated in Table 1. For h = 0.05, with 23327 vertices and 46012 triangles, we have solved (4.1)-(4.2) for a = 1.8 Figure 6. Test 1 -Computation times (in seconds) and numbers of iterates to reach an error ≤ ε = 10 −6 Re for Re = 141., a = 1.5 and µ = 2 − a = 0.5.  and Re = 141 and we have accepted that the computed solution is, for practical purposes, our Nash quasiequilibrium. Then, we have measured the relative L 2 errors associated with the computations with the other grids. See Table 1 and Figure 7 for quantitative information.

Test 2
We consider the flow of a fluid in a pipe. The fluid enters with a parabolic horizontal profile on the left. It is allowed to exit vertically upwards and downwards before it reaches the observation domains O 1 and O 2 (see Fig. 10). Our objective is to know how we have to act in two specific control regions (ω 1 and ω 2 ) to make the fluid velocity reach a Nash quasi-equilibrium associated with u 1d and u 2d (see Fig. 9). The domain and a mesh are depicted in Figure 8 (see the legend for explanations). What we do now is equivalent to solve a   bi-objective boundary control problem, where the controls act on the parts of ∂ω 1 and ∂ω 2 that "touch" the pipe.
The finite element method has been the same as before. The data u id are similar to those for Test 1: u 1d = ∇ × ψ 1d , where ψ 1d is the solution to (4.25) and u 2d ≡ 0. Thus, the "desired" (ideal) configurations correspond to a uniformly rotating flow in O 1 and a fluid at rest in O 2 (see Fig. 9). The stopping test and the initial controls (f 0 1 , f 0 2 ) have been as in the previous section. Again, all algorithms have provided the same solution in all cases.
We present in Figure 10 the velocity field corresponding to a "free" state, that is, an uncontrolled solution. Some results are depicted in Figure 11 for large a and various Re.
We have also checked the convergence of the numerical aproximations with results similar to those in Section 4.4.1. We have solved (4.1)-(4.2) for a = 1.99, Re = 240 and h = 0.025, which corresponds to a mesh with 30405 vertices and 59448 triangles. We have assumed that this computed solution coincides with the Nash quasi-equilibrium and we have measured the relative L 2 errors associated with the other approximations. For details, see Table 2 and Figure 12.

Conclusions, further comments and open questions
In this paper, we have been able to analyze and sometimes solve several that bi-objective control problems for stationary PDEs. The existence of Nash equilibria (or at least quasi-equilibria) has been established. We  have also deduced related characterizations (optimality systems) and several iterative algorithms have been introduced.
Our numerical experiments, that concern Navier-Stokes control, provide useful illustrations and can be viewed as a confirmation of the interest of the subject. • It is unknown whether convergence results can be established for ALG 2 and ALG 3 (see Sect. 2.3). It is reasonable to expect this at least for small a/µ. The same can be said for ALG 5 and ALG 6 (Sect. 3.3) and ALG 8 and ALG 9 (Sect. 4.3). The reason is that, in these iterates, we respect the gradient or conjugate gradient main idea and, to go from f n i to f n+1 i , we use either the direction ∂J i /∂f i or a suitable approximation. • It is also unknown whether Nash equilibria exist for the Navier-Stokes system (4.1) and the functionals given in (4.2). As can be deduced from the arguments in Section 4.2, the lack of uniqueness of the system and the lack of convexity of the functionals are nontrivial difficulties in this context. Indeed, it does not seem easy to adapt the arguments in the proof of Theorem 3.2. One could try to find an equilibrium directly, forgetting the optimality system and minimizing J 1 and J 2 in an iterative process. But, to our knowledge, no argument has been given up to now. • Finally, the extension of the results in Section 4 to the time-dependent Navier-Stokes system is also a challenge. It is not difficult to prove some related assertions in the 2D case. However, once more, not many things can be said when the spatial dimension is three.