Control Problem Related to 2D Stokes Equations with Variable Density and Viscosity

: We study an optimal control problem for the stationary Stokes equations with variable density and viscosity in a 2D bounded domain under mixed boundary conditions. On in-ﬂow and out-ﬂow parts of the boundary, nonhomogeneous Dirichlet boundary conditions are used, while on the solid walls of the ﬂow domain, the impermeability condition and the Navier slip condition are provided. We control the system by the external forces (distributed control) as well as the velocity boundary control acting on a ﬁxed part of the boundary. We prove the existence of weak solutions of the state equations, by ﬁrstly expressing the ﬂuid density in terms of the stream function (Frolov formulation). Then, we analyze the control problem and prove the existence of global optimal solutions. Using a Lagrange multipliers theorem in Banach spaces, we derive an optimality system. We also establish a second-order sufﬁcient optimality condition and show that the marginal function of this control system is lower semi-continuous.


Introduction
In this work, we study an optimal control problem for the stationary Stokes equations with variable density and viscosity in a domain Ω ⊂ R 2 , which is assumed to be a bounded and connected set with boundary Γ def = ∂Ω of class C 1 . Specifically, we consider the following system of partial differential equations (PDEs): Here, the unknown functions are the velocity field u : Ω → R 2 , the hydrostatic pressure π : Ω → R, and the density ρ : Ω → R + for a nonhomogeneous fluid that flows through the domain Ω. The vector function f : Ω → R 2 represents a distributed control that acts on the motion equation and belongs to a closed and convex set F ; the function µ : Ω → R + describes the dynamic viscosity of the fluid; D(u) is the symmetric part of the velocity gradient (the deformation rate tensor), i.e., D(u) def = 1 2 ∇u + (∇u) .
In order to describe boundary conditions for the velocity field and the density function in (1), we divide the boundary in three open disjoint parts: Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 , where Γ 1 and Γ 2 are in-flow/out-flow parts and Γ 3 = Γ \ (Γ 1 ∪ Γ 2 ) is the lateral part (solid walls) of the boundary. Then, we impose the following boundary conditions: The functions ρ 0 and u 0 are defined on Γ 1 and describe a Dirichlet boundary condition. The function g describes a Dirichlet boundary control for u on Γ 2 and belongs to a closed convex set G u 0 .
Additionally, we assume the following: Γ 1 is an arcwise connected and closed set on Γ, with |Γ 1 | > 0; (4) either u 0 · n > 0 on Γ 1 (out-flow) or u 0 · n < 0 on Γ 1 (in-flow), where |Γ 1 | denotes the length of the curve Γ 1 and n is the outward unit normal vector to the curve Γ. Observe that condition (5) implies that Γ 1 is either a part of the boundary where the fluid flows outward or a part where the fluid flows inward. An example of the flow domain Ω is given in Figure 1.
For the sake of brevity, the nonhomogeneous Stokes system (1) with mixed boundary conditions (2), (3), (6) will be referred to as problem NS-MBC in what follows.
The Navier slip boundary condition was proposed by Navier [1]; it is based on a balance between the fluid velocity tangent to the surface and the rate of strain at the boundary, i.e., the tangential component of the viscous stress at the boundary should be proportional to the tangential velocity. The components of the normal velocity to the surface is naturally zero, as mass cannot penetrate an impermeable solid surface. This type of boundary condition is involved when one studies boundary layer problems, such as in channels or Couette flows and is well justified by Jäger-Mikelić [2,3]. In addition, the real number α ≥ 0 is the friction coefficient, which measures the tendency of the fluid to slip on the part Γ 3 . When α = 0, the fluid slips on Γ 3 without friction and there are no boundary layers, while if α goes to +∞, then the friction is so intense that the fluid is almost at rest near the boundary [4] (see also [5]). Thus, from the physical point of view, boundary conditions (6) make more sense than the classical Dirichlet boundary conditions. Moreover, in order to define a correct variational formulation (weak solutions concept) of problem NS-MBC, the classical Green identities are not applicable in this case. The above leads us to study other results of integration by parts (see Lemma 1,below). Additionally, in order to correctly control the H σ -norm in terms of the L 2 -norm of the deformation rate tensor, i.e., the norm D(u) 2 , we must employ the Korn inequality (see [6], p. 52), which is not usual.
Some examples of flow phenomena that might require the introduction of a Navier slip boundary condition have been presented by Fujita [7,8], namely: a drainage or canal in which its bottom is covered with a layer of mud and pebbles; the blood flow in the vein of a patient suffering from arterial sclerosis; an avalanche of water and rocks; a flow of iron coming out of a smelting. Some other applications can be consulted in [9][10][11].
It is well known that fluids in which the density and viscosity are variable have several important applications, both in natural phenomena as well as in industry. These fluids present a wider range of interesting phenomena, compared to fluids with constant density and viscosity (also called homogeneous fluids), and there are several mathematical models that describe different physical situations. In particular, some examples of the physical modeling of mixtures of incompressible fluids can be found in [12] and the respective mathematical analysis of these models were developed in [13][14][15]. We emphasize that all these works only consider fluids with constant viscosity and Dirichlet boundary conditions. Furthermore, as far as we know, there are no previous studies on the existence of weak solutions of problem NS-MBC, and analyses of optimal control problems have not been carried out either.
The above motivates the analysis of problem NS-MBC. The main mathematical difficulty that arises in the study of this problem in relation to the Stokes system with variable density and viscosity is due to several nonlinearities in the problem (−div [µD(u)], ρf in Equation (1) 1 and u · ∇ρ in Equation (1) 2 ), which, due to the product with the unknown density and viscosity, are much harder to deal with than the corresponding terms ones appearing in the Stokes equations in which the density and viscosity are known constants. Moreover, since the nonlinearity of the term u · ∇ρ is not related to a monotone operator, we must work with arguments that are different than what is typical when dealing with elliptic or parabolic problems and nonlinearities and restrictions related to monotone operators. To overcome these difficulties, we will use the stream formulation presented by Frolov [16] that is, the fluid density is represented in terms of its stream function through another function determined by the boundary conditions. This representation allows us to drop the continuity equation, which facilitates the introduction of the weak solution concept of problem NS-MBC. The stream formulation was used in [17][18][19][20][21].
The purpose of this work is to study an optimal control problem related to weak solutions of problem NS-MBC. We would like to point out that the literature to the mathematical analysis of optimal control problems associated with systems of partial differential equations modeling the motion of viscous incompressible fluids with variable density and viscosity is scarce. The main difficulty lies in the mathematical handling. In fact, the systems that model the behavior of viscous and incompressible fluids with constant density can have two types of character: either elliptical, for the steady state, or parabolic, for the non-stationary state. In both cases, when an approximation argument is used, such as the Galerkin method, for instance, higher-order estimates are usually obtained for the approximate solutions, which facilitates the passage to the limit in the nonlinear terms and, consequently, it is more simple to achieve the required results. On the other hand, for systems that model the motion of fluids with variable density, the equations that compose the system present nonlinearities that are more difficult to handle since they makes it impossible to obtain higher-order estimates and make it difficult to pass to the limit in the nonlinear terms. The system of equations that we have considered in our work (problem NS-MBC), the first equation of system (1) has an elliptic character but is coupled with a first-order transport equation for the density ρ (see Equation (1) 1 ), which gives a hyperbolic character to the model. These drawbacks make it difficult to obtain estimates with high regularity and increase the difficulties in passing to the limit.
Let us mention the available literature on the analysis of optimal control problems for nonhomogeneous fluid flows. The beginning of the study of such problems dates back to the paper of Illarionov [17]. He investigated optimal boundary control for a model of 2D steady-state flows of a nonhomogeneous incompressible fluid under the assumption that the viscosity µ is constant. Using the external forces field as a control function, Mallea-Zepeda et al. [19] proved the existence of optimal solutions to this model and obtained the first-order necessary optimality conditions. Boundary control for a 2D stationary system of micropolar fluid with variable density was studied in [20]. Certain classes of optimal control problems for the 2D Boussinesq equations with variable density and constant viscosity are analyzed in [21].
We also mention that there is a large number of mathematical works devoted to the study of optimal control problems for PDEs describing flows of a homogeneous fluid (see, for example, refs. [22][23][24][25][26] and the numerous references therein). Such problems are currently fairly well understood, while control and optimization problems for nonhomogeneous fluid flows remain a serious challenge. Motivated by this fact, we performed our study, which can be considered a starting point for future investigations of non-Newtonian fluid models, where more general inhomogeneities can be analyzed.
The outline of the present paper is as follows: In Section 2, we fix the notation, introduce the function spaces to be used and establish the concept of weak solutions of problem NS-MBC, using the stream formulation. In Section 3, we prove the existence and uniqueness of weak solutions of system NS-MBC, in order to establish that the admissible solutions set is nonempty (see (39) below). In Section 4, we introduce the optimal control problem to this system and prove the existence of global optimal solutions. Moreover, using a Lagrange multipliers theorem in Banach spaces, we derive an optimality system, and, establishing a coercivity condition for the Lagrangian function, we obtain a second-order sufficient optimality condition. Finally, we show that the marginal function related to this control system is lower semi-continuous with respect to the one-sided Hausdorff distance.
For the reader's convenience, in Appendix A (see Table A1), we collect the main symbols used in this paper and explain their meaning.

Preliminaries
In this section, we introduce the notation, function spaces and principal results that we will use throughout the work.

Notation and Function Spaces
Let E be a Banach space. By 2 E we denote the collection of all subsets of E (the power set of E). By definition, put Let U and W be subsets of E. By d E (U , W ) we denote the one-sided Hausdorff distance from the set U to the set W, that is, For m ∈ N and 1 ≤ p ≤ +∞, we will use the Sobolev space H m (Ω) def = W m,2 (Ω) and the Lebesgue space L p (Ω) with norms denoted by · H m and · L p , respectively. In particular, when p = 2, L 2 (Ω) is a Hilbert space; in this case, the L 2 -norm and L 2 -scalar product are represented by · and (·, ·), respectively. In addition, the L p (Γ)-norm is denoted by · L p (Γ) . The corresponding spaces of vector-valued functions are denoted by bold script; for instance, H 1 (Ω), L p (Ω), and so on. Recall that the restriction of a function is the trace operator (see [27], § 2.4.2).
We also use the following divergence-free Banach spaces (solenoidal spaces): in Ω and v · n = 0 on Γ 3 } equipped with the usual norm and scalar product of H 1 (Ω); By C l (Ω) (resp. C l (R)) we denote the Banach space of functions that are l times continuously differentiable on Ω (resp. on R) and by C l b (Ω) (resp. C l b (R)) the space of l times continuously differentiable and bounded functions on Ω (resp. on R).
Moreover, if X is a general Banach space, its dual topological is denoted by X and the pairing duality by ·, · X or ·, · , when this does not lead to ambiguities. In particular, H σ denotes the dual of H σ , and V σ denotes the dual of V σ .
Furthermore, the space of the traces of vector functions from H 1 (Ω) is defined by is the lifting (trace extension) operator such that γ 0 • L Ω = Id (see, for example, ref. [28], Chapter III, Theorem III.2.22).
For a subset S of Γ, we introduce the subspace H 1/2 0 (S) as follows: Note that the injection H 1/2 0 (S) → L 2 (S) is continuous; that is, there exists a positive constant C which depends only on Ω such that ϕ L 2 (S) ≤ C ϕ H 1/2 0 (S) (see [27], § 2.4.2). The letter C denotes a positive constant, independent of the state (u, π, ρ) and the control (f, g), but its value may change from line to line.

Weak Solutions
Let us assume that the following conditions are fulfilled: there exists constants µ * and µ * such that 0 < µ * ≤ µ(x) ≤ µ * , for any x ∈ Ω, To define the concept of weak solutions of problem NS-MBC, following the ideas of Frolov [16] (see, also [17][18][19][20][21], for more details), we express the fluid density in terms of the stream function. Namely, assuming conditions (4) and (5), we can express the fluid density as follows: where η is a function from the space C 0 b (R), and N : H σ → H 2 (Ω) is the linear and continuous operator (see [17], Lemma 2.1) that takes each vector function u ∈ H σ to a function ψ (stream function) such that where e 3 def = (0, 0, 1), u def = (u 1 , u 2 , 0), x 0 is the initial point of the curve Γ 1 , and Γ 1 (x 0 , x) is the portion of Γ 1 lying between the points x 0 and x. In the Cartesian coordinate system, the first equality of (11) is equivalent to For any η ∈ C 1 b (R), the velocity field u and the density ρ, defined as in (10), satisfy relation (1) 2 . Indeed, using the chain rule and (12), we obtain where η denotes the first derivative of η. Therefore, we obviously have Note also that if η ∈ C 0 b (R), then condition (1) 2 is satisfied at least in the weak sense, that is, This identity is not difficult to establish by applying the procedure of regularization for the function η and passage to the limit. Now it remains to choose a function η so that the density function ρ satisfies boundary condition (3). To this end, we introduce the function ψ 0 : Γ 1 → R by the following formula: Clearly, ψ 0 ∈ C 0 (Γ 1 ). Moreover, from conditions (4) and (5) it follows that ψ 0 is strictly monotone on the curve Γ 1 , and hence, there exists the inverse function ψ −1 (4) and (5) hold, In the sequel, we assume that the function η satisfying (13) is fixed.
Let us show that boundary condition (3) is satisfied. Indeed, taking into account (10) and (13), we obtain for any u ∈ H σ such that u = u 0 on Γ 1 and for all x ∈ Γ 1 .

Remark 1.
If the curve Γ 1 is of class C k and the following inclusions hold, with k ≥ 1; then there exists a function η ∈ C k (R) satisfying relations (13) (see [17]).
To obtain the weak formulation of problem NS-MBC, as usual, we take the L 2 -scalar product of equality (1) 1 with a test function v ∈ V σ . Next, we use the relation (∇π, v) = 0 and the following result on integration by parts: Lemma 1 (see [29]). Suppose that a vector function u ∈ H 2 (Ω) satisfies boundary conditions (6), a vector function v ∈ H 1 (Ω) satisfies the condition v · n = 0 on Γ, and µ ∈ C 1 (Ω). Then, we have where the symbol: denotes the component-wise product of matrices.
Thus, we arrive at the following definition of weak solutions to system NS-MBC.

Definition 1 (Weak solution)
. We say that a pair (u, ρ) is a weak solution of problem NS-MBC if u ∈ H σ , ρ = η(Nu) and the following weak formulation holds: where the operator F :

Existence and Uniqueness of Weak Solutions
In this section, using the Schauder fixed-point theorem, we prove the existence and uniqueness of weak solutions of system NS-MBC; as far as we know, there are no previous results on the existence of weak solutions of this problem.
In order to prove the existence of a solution of (14), we reduce this problem to a new problem with homogeneous boundary conditions in terms of new unknown u. The following result allows us to make the desired reduction.
Now, rewriting u ∈ H σ in the form u = u e + u, where u ∈ V σ is a new unknown function, we obtain that problem (14) is equivalent to the following problem: Find u ∈ V σ such that We define the operator T : where the vector function u satisfies the following relation: We observe that a fixed-point of the operator T is a solution of problem (16). In the next lemma, we prove that the operator T satisfies the hypotheses of the Schauder fixed-point theorem.
Lemma 3. The operator T : V σ → V σ is well-defined and completely (weak-to-strong) continuous. Moreover, the following estimate holds: where C > 0 is a constant, which depends only on the domain Ω and the constants µ * and µ * .
Proof. We introduce the bilinear symmetric form a : V σ × V σ → R and the linear functional L : V σ → R given by Using (19) and (20), we rewrite problem (17) as follows: Find u ∈ V σ such that Clearly, the form a(·, ·) is continuous on V σ × V σ . Moreover, from (8) and (19) it follows that a( u, u) ≥ µ * D( u) 2 for all u ∈ V σ . Thus, we deduce that the form a(·, ·) is coercive.
On the other hand, using the Hölder and Poincaré inequalities, we find where the constant C satisfies the inequality D(u) ≤ C ∇u , for all u ∈ H σ . Hence, the linear functional L is continuous.
Therefore, by the Lax-Milgram lemma (see, for example, ref. [30], Theorem 9.14), we deduce that there exists a unique u ∈ V σ satisfying (21), which implies that the operator T is well-defined. Now, from (21)-(23), taking into account (15), we find which implies estimate (18). Finally, let us prove that the operator T is completely continuous. We consider a sequence { u m } m≥1 ⊂ V σ and an element u belonging to V σ such that u m → u weakly in V σ , as m → +∞.
Then, since the injection of H 2 (Ω) in C 0 (Ω) is compact and the operator N is continuous, we deduce that for some subsequence of { u m } m≥1 , still denoted by { u m } m≥1 , the following convergence holds: In view of u m = T( u m ), we have Then, taking the difference between (17) and (25), we obtain Setting v = u m − u into the last equality, by the Hölder and Poincaré inequalities we derive the estimate as follows: Therefore, passing to the limit in (26) as m goes to +∞, and taking into account the strong convergence (24), we obtain This means that T( u m ) → T( u) strongly in V σ , as m → +∞.
Thus, the proof is complete.
From Lemma 3, we have the following result on the existence of solutions to system (14). (4), (5) and (7)-(9) hold; then problem NS-MBC has at least one weak solution in the sense of Definition 1. Moreover, if a pair (u, ρ) is a weak solution of problem NS-MBC, then the function u satisfies the following inequality

Theorem 1 (Existence). Suppose conditions
where C is a positive constant, which depends only on the domain Ω and the constants µ * and µ * .
Proof. From Lemma 3, it follows that the operator T, defined in (17), satisfies the conditions of the Schauder fixed-point theorem. Thus, we deduce that there exists u ∈ V σ such that T u = u, which is a solution of system (16). Therefore, we conclude the existence of u ∈ H σ , with u = u + u e , satisfying system (14). Finally, using inequalities (15) and (18), we obtain which implies (27).

Theorem 2 (Uniqueness).
In addition to the assumptions of Theorem 1, we suppose that the function η belongs to the space C 1 b (R) and the constant µ * is sufficiently large such that where C is a fixed positive constant, which depends only on Ω. Then, the solution u ∈ H σ of system (14), provided by Theorem 1, is unique.

Optimal Control Problem
In this section, we give the statement of the control problem under study. Then, we prove the existence of at least one global optimal solution. Using a Lagrange multipliers result on the existence of Lagrange multipliers theorem in Banach spaces, we derive an optimality system and a second-order sufficient optimality condition. In addition, we show that the marginal function related to this control system (see problem (38) below) is lower semi-continuous with respect to the one-sided Hausdorff distance.
Let us consider sets F and G such that Moreover, assume the following: Note that from the inclusion G ∈ P c, cl [H 1/2 0 (Γ 2 )] it follows that the set G u 0 is convex and closed in the space H 1/2 0 (Γ 2 ). We consider a vector function f ∈ F describing a distributed control for the motion equations in Ω and a vector function g ∈ G u 0 describing a boundary control for u on Γ 2 .
For simplicity, we use the product X def = H σ × F × G u 0 , and we introduce a cost functional J : X → R by the following formula: where α f and α g are nonnegative real constants, and J : H σ → R is a weakly lower semi-continuous functional with inf u∈H σ J(u) > −∞. (36) Examples of weakly lower semi-continuous functionals, interesting from the physical point of view, are as follows: The functional J 1 describes the deviation of the flow velocity u from a given desired vector field u d . The functional J 2 measures the vorticity of the velocity field u. The functional J 3 describes the total resistance in a fluid due to viscous friction (see [17,31]).
The constants α f and α g , given in (35), measure the cost of the control. Assume that at least one of the following two conditions holds: Thus, we define the following constrained minimization problem for system (14): Find a triplet (u, f, g) ∈ X that minimizes the functional , subject to (u, f, g) satisfies system (14).
The set of admissible solutions of control problem (38) is defined by Related to problem (38), we have the following definitions.
Definition 2 (Global optimal solution). A triplet ( u, f, g) ∈ Z ad (F , G u 0 ) is called a global optimal solution of control problem (38) if By Z opt (F , G u 0 ) we denote the set of all global optimal solutions to control problem (38).

Existence of Global Optimal Solutions
In this subsection, we will prove the solvability of optimization problem (38).

Theorem 3 (Existence of optimal solutions).
Under the assumptions of Theorem 1, suppose that conditions (33), (34) and (36) hold. Moreover, suppose that at least one of conditions (i) and (ii) given in (37) is satisfied. Then, optimal control problem (38) has at least one global optimal solution ( u, f, g) ∈ Z ad (F , G u 0 ) in the sense of Definition 2.

Proof. From Theorem 1 it follows that
Since the functional J is bounded from below, we see that there exists a minimizing Moreover, if at least one of the conditions given in (37) is satisfied, then there exists a positive constant C such that On the other hand, by definition of the set Z ad (F , G u 0 ), for each m ∈ N, the sequence {(u m , f m , g m )} m≥1 satisfies system (14); thus, from estimate (27) we have that there exists a constant C > 0, independent of m, such that Therefore, using estimates (41) and (42) and the fact that the admissible controls set F × G u 0 ⊂ L 2 (Ω) × H 1/2 0 (Γ 2 ) is closed and convex (in particular, is weakly closed in the space L 2 (Ω) × H 1/2 0 (Γ 2 )), we deduce that there exist a triplet ( u, f, g) ∈ X and a subsequence of {(u m , f m , g m )} m≥1 , still denoted by {(u m , f m , g m )} m≥1 , such that the following convergences hold as m → +∞: Since u m = u 0 on Γ 1 and u m = g m on Γ 2 , from (43) it follows that u = u 0 on Γ 1 , u = g on Γ 2 and η(N u) = ρ 0 on Γ 1 . Thus, we deduce that the triplet ( u, f, g) satisfies the boundary conditions given in (14).
Next, we observe that for all v ∈ V σ the following estimate holds: where Since the operator N : Using (43) 2 and (45), we derive from (44) the following: Then, the convergences (43) and (46) allow us pass to the limit in system (14) written by (u m , f m , g m ), as m goes to +∞; thus we obtain that ( u, f, g) is a solution of (14). Consequently, the triplet ( u, f, g) belongs to the set Z ad (F , G u 0 ) and On the other hand, since the functional J is weakly lower semi-continuous, we have which, together with (47), imply (40). Thus, control problem (38) has at least one global optimal solution.

Optimality System and Second-Order Sufficient Optimality Condition
In this subsection, we derive an optimality system for optimal control problem (38) and establish a second-order sufficient optimality condition.
In order to obtain first-order necessary optimality conditions and derive an optimality system for local optimal solutions of control problem (38), we reformulate this problem in the abstract context given by Zowe and Kurcyusz [32]. The method for obtaining first-order necessary optimality conditions, provided by [32], was also previously used in [33,34], among others.
Notice that the set of admissible solutions of problem (49) is

Remark 2.
For simplicity, in what follows we carry out our study taking the following: Then, for problem (38), the cost functional is Thus, we can easily deduce the following results concerning to the differentiability of the functional J and the operator S. Lemma 4. The functional J : X → R is Fréchet differentiable and the Fréchet derivative of J at the point r = ( u, f, g) ∈ X in the direction s = (w, y, z) ∈ X is given by J ( r)[s] = ( u − u d , w) + (µD( u), D(w)) + α f ( f, y) + α g g, z Γ 2 .

Lemma 5.
Let η ∈ C 1 b (R). The operator S : X → W is continuously Fréchet-differentiable and the Fréchet derivative of S at the point r = ( u, f, g) ∈ X in the direction s = (w, y, z) ∈ X is the linear operator S ( r)[s] = (S 1 ( r)[s], S 2 ( r)[s]), where and the linear operator Q : Following [32], we say that r = ( u, f, g) is a regular point for problem (49) if for each pair (a, b) ∈ W there exists an element s = (w, y, z) ∈ H σ × C( f) × C( g) such that is the conical hull of ( f, g) in F × G u 0 . Lemma 6. Suppose that r = ( u, f, g) ∈ Z ad (F , G u 0 ); then r is a regular point.
Proof. Since the pair (0, 0) belongs to C( f) × C( g), we deduce that it is sufficient to prove the existence of w ∈ H σ satisfying the following system The existence of w satisfying (51) can be obtained by arguing similarly as in the proof of Theorem 1. Now, we prove the existence of Lagrange multipliers for problem (49); thereafter, we derive an optimality system. Theorem 4. Let r = ( u, f, g) ∈ Z ad (F , G u 0 ) be a local optimal solution of problem (49) and assume that η ∈ C 1 b (R). Then, there exist Lagrange multipliers (λ,ϕ) ∈ V σ × H −1/2 0 (Γ \ Γ 3 ) such that the following variational inequality holds: for all s = (w, y, z) ∈ H σ × C( f) × C( g).

Remark 3.
Suppose α f and α g are strictly positive. Then, taking into account that F × G u 0 is a closed convex set, from optimality conditions (54) and (55) and a well-known theorem on the projection onto a closed convex set (see, for example, ref. [35], Chap. 5, Theorem 5.2), we derive the following: and g = Proj Finally, we establish a second-order sufficient optimality condition for control problem (49) via the derivation of an X-coercivity condition on the second derivative of the Lagrangian function, which assures that an admissible solution r = ( u, f, g) is a local optimal solution.
We observe that the Lagrangian L related to optimal control problem (49) is given by where the operators S 1 and S 2 are defined in (48) and the functional J is given in (50).

Lemma 7.
Let r = ( u, f, g) ∈ Z ad (F , G u 0 ) be a local optimal solution of problem (49) and assume that η ∈ C 1 b (R). If the constant µ * is sufficiently large such that where C is a positive constant depends only on Ω. Then, the Lagrange multiplier λ, provided by Theorem 4, satisfies the following inequality Proof. Taking w = λ ∈ V σ in (53), we have Applying the Hölder and Young inequalities and Sobolev embeddings and taking into account that the operator N is continuous, we obtain Then, using (60)-(62), we derive from (59) the following estimate: From hypothesis (57), it follows that µ * > 2C η L ∞ f . Thus, the result follows from (63).
In the following result, we establish a second-order sufficient condition for optimal control problem (49).
Consequently, r is a local optimal solution of control problem (49).
Proof. Let s = (w, y, z) ∈ X. Notice that the Lagrangian L related to control problem (49) is twice Fréchet-differentiable, and the second derivative of L at the admissible solution r = ( u, f, g) in the directions (s, s) ∈ X × X is given by Using the Hölder inequality, Sobolev embeddings and the continuity of the operator N, we derive Furthermore, from the Korn inequality (see [6], p. 52) it follows that there exists a constant C K = C K (Ω) > 0 such that Taking into account (65) and (66), we deduce from (64) the following estimate: Next, from (58) we derive which, together with (67), imply .
In particular, if s = (w, y, z) ∈ ker(S ( r)) with (y, z) ∈ C( f) × C( g), then from [36] we conclude that the triplet r is a local optimal solution of control problem (49).

Marginal Function
In the studying of optimal solutions, it is important to investigate the case when the collection of all admissible controls (in our problem, the set F × G u 0 ) can be expanded/reduced. Following the ideas developed in [37][38][39], we introduce the concept of the marginal function, which shows how the minimal value of the cost functional J changes under a variation of the set F × G u 0 .

Conclusions
In this article, we have studied an optimal control problem for the 2D stationary Stokes equations with variable density and viscosity, using nonhomogeneous Dirichlet boundary conditions on one part of the boundary of the flow domain and the Navier slip boundary conditions on the other part. Expressing the fluid density in terms of the stream function, we proved the existence and uniqueness of weak solutions of the dynamical equations with the same regularity (H 1 -regularity) as weak solutions of the classical Stokes system with constant density and viscosity under Dirichlet boundary conditions acting on the whole boundary. For the optimal control problem, we controlled the system, applying a control on a part of the boundary and another control acting as an external force on the domain, i.e., a distributed type control. We proved the existence of at least one global optimal solution and derived an optimality system. Establishing a coercivity property for the Lagrangian function, we obtained a second-order sufficient optimality condition for a local optimal solution. Furthermore, we introduced the concept of the marginal function, which shows how the minimal value of the cost functional changes under a variation of the set of admissible controls. Then, we proved that the marginal function is lower semi-continuous. Therefore, one can deduce that the control system is stable in the sense that it is not possible to achieve a significant improvement in the optimal value of the cost functional without an essential extension of the admissible controls set. Finally, considering that the physical phenomena of fluid flow are intrinsically transient, for future research, we will consider the analysis of an unsteady-state situation with time-dependent control.

Conflicts of Interest:
The authors declare no conflict of interest.