OPTIMAL L 2 -CONTROL PROBLEM IN COEFFICIENTS FOR A LINEAR ELLIPTIC EQUATION.

. In this paper we study we study a Dirichlet optimal control problem associated with a linear elliptic equation the coeﬃcients of which we take as controls in the class of integrable functions. The characteristic feature of this control object is the fact that the skew-symmetric part of matrix-valued control A ( x ) belongs to L 2 -space (rather than L ∞ ). In spite of the fact that the equations of this type can exhibit non-uniqueness of weak solutions, the corresponding OCP, under rather general assumptions on the class of admissible controls, is well-posed and admits a nonempty set of solutions [9]. However, the optimal solutions to such problem may have a singular character. We show that some of optimal solutions can be attainable by solutions of special optimal control problems in perforated domains with ﬁctitious boundary controls on the holes.

Minimize I(A, y) = y − y d 2 L 2 (Ω) + Ω (∇y, A sym ∇y) R N dx subject to the constraints − div A sym ∇y + A skew ∇y = f in Ω, where (A sym , A skew ) ∈ L ∞ (Ω; R N ×N ) × L 2 (Ω; R N ×N ) are respectively the symmetric and antisymmetric part of the control A, y d ∈ L 2 (Ω) and f ∈ H −1 (Ω) are given distributions, and A Ad denotes the class of admissible controls which will be precised later.
The characteristic feature of this problem is the fact that the skew-symmetric part of matrix A(x) belongs to L 2 -space (rather than L ∞ ). As a result, the existence and uniqueness of the weak solutions to the corresponding boundary value problem (1) are usually drastically different from the properties of solutions to the elliptic equations with L ∞ -matrices in coefficients. In most of the cases, the situation can deeply change for the matrices A with unremovable singularity. As a rule, some of the weak solutions can be attained by the weak solutions to the similar boundary value problems with L ∞ -approximated matrix A. However, this type does not exhaust all weak solutions to the above problem. There is another type of weak solutions called non-variational [20,22], singular [3,13,14,19], pathological [16,17] and others. As for the optimal control problem (1) we have the following result [9] (see [8] for comparison): for any approximation {A * k } k∈N of the matrix A * ∈ L 2 Ω; S N skew with properties {A * k } k∈N ⊂ L ∞ (Ω; S N skew ) and A * k → A * strongly in L 2 (Ω; S N skew ), optimal solutions to the corresponding regularized OCPs associated with matrices A * k always lead in the limit as k → ∞ to some admissible (but not optimal in general) solution ( A, y ) of the original OCP (1). Moreover, this limit pair can depend on the choice of the approximative sequence {A * k } k∈N . However, as follows from counter-example, given in [9], it is possible a situation when none of optimal solutions to OCP (1) can be attainable in such way. Therefore, the aim of this paper is to discuss a scheme of approximation for OCP (1) in order to attain the other types of optimal solutions, and derive the first order optimality system to this problem.
In order to illustrate the difficulties on the approximations of the OCPs due to the possible existence of variational and non-variational solutions, we present some numerical simulations in section 5. In section 3 we give a precise description of the class of admissible controls A ad ⊂ L 2 Ω; R N ×N which guarantee that non-variational solutions can be attained through the sequence of optimal solutions to OCPs in special perforated domains with fictitious boundary controls on the boundary of holes. Namely, we consider the following family of regularized OCPs subject to the constraints − div A sym ∇y + A skew ∇y = f in Ω ε , y = 0 on ∂Ω, ∂y/∂ν A = v on Γ ε , y ∈ H 1 0 (Ω ε ; ∂Ω),

OCP FOR LINEAR ELLIPTIC EQUATIONS 597
where Ω ε is the subset of Ω such that ∂Ω ⊂ ∂Ω ε , σ > 0, and A(x) S N := max i,j=1,...,N |a ij (x)| ≤ ε −1 a.e. in Ω ε . Here, v stands for the fictitious control. We show that OCP (2) has a nonempty set of solutions (A 0 ε , v 0 ε , y 0 ε ) for every ε > 0. Moreover, as follows from (2) 1 , the cost functional I ε seems to be rather sensitive with respect to the fictitious controls. Due to this fact, we prove that the sequence (A 0 ε , y 0 ε ) ε>0 gives in the limit an optimal solution (A 0 , y 0 ) to the original problem.
The main technical difficulty, which is related with the study of the asymptotic behaviour of OCPs (2) as ε → 0, deals with the identification of the limit of two weakly convergent sequences. Due to the special properties of the skew-symmetric parts of admissible controls A ∈ A ad ⊂ L 2 Ω; S N , we show that this limit can be recovered in an explicit form. We also show in this section that the energy equalities to the regularized boundary value problems can be specified by two extra terms which characterize the presence of the-called hidden singular energy coming from L 2 -properties of skew-symmetric components A skew of admissible controls.
In conclusion, in Section 4, we derive the optimality conditions for regularized OCPs (2) and show that the limit passage in optimality system for the regularized problems (2) as ε → 0 leads to the optimality system for the original OCP (1).
Let us mention briefly that the problem dealt here can be of multiple interests. First, as we explained it in [9] such optimal control problems with singular coefficients arise naturally in the context of porous media. Second, singular matrices in elliptic equations can be viewed as a model of diffusion coefficients that become unbounded in some region of a heated domain and are sometimes called thermical leak. In that spirit, the optimal control problem presented here is the one hand natural if one wants to optimize some heating according to some goals in the presence of such thermical leak while, on the other hand, once you know the singular sets of the coefficients, in order to make some prediction one has to make computations. In some sense this paper indicates some difficulties to handle the good choice of numerical recipes since different types of solution may appear.

Notation and preliminaries.
Let Ω be a bounded open connected subset of R N (N ≥ 2) with Lipschitz boundary ∂Ω. By H 1 0 (Ω) we denote the closure of C ∞ 0 (Ω)-functions in the Sobolev space H 1 (Ω), while H −1 (Ω) denotes the dual of H 1 0 (Ω). Let Γ be a part of the boundary ∂Ω with positive (N − 1)-dimensional measures. We consider C ∞ 0 (R N ; Γ) = ϕ ∈ C ∞ 0 (R N ) : ϕ = 0 on Γ , and denote H 1 0 (Ω; Γ) its closure with respect to the norm y = Ω ∇y 2 Here, S N skew stands for the set of all skew-symmetric matrices C = [c ij ] N i,j=1 , whereas S N sym is the set of all N × N symmetric matrices.
Let L 2 (Ω) N (N −1) 2 = L 2 Ω; S N skew be the normed space of measurable squareintegrable functions whose values are skew-symmetric matrices. By analogy, we can define the space L 2 (Ω) N (N +1) 2 = L 2 Ω; S N sym . Let A(x) and B(x) be given matrices such that A, B ∈ L 2 (Ω; S N skew ). We say that these matrices are related by the binary relation on the set L 2 (Ω; S N skew ) (in symbols, A(x) B(x) a.e. in Ω), if

THIERRY HORSIN, PETER I. KOGUT AND OLIVIER WILK
Here, L N (E) denotes the N -dimensional Lebesgue measure of E ⊂ R N defined on the completed borelian σ-algebra. We define the divergence div A of a matrix A ∈ L 2 Ω; M N as a vector-valued distribution d ∈ H −1 (Ω; R N ) by the following rule where a i stands for the i-th row of the matrix A.
For fixed two constants α and β such that 0 < α ≤ β < +∞, we define M β α (Ω) as a set of all matrices A = [a i j ] in L ∞ (Ω; S N sym ) such that Let A ∈ L 2 Ω; M N be an arbitrary matrix. In view of the representation A = A sym + A skew , we can associate with A the form ϕ(·, ·) A : By analogy with [9], we introduce the following concept.
with some constant c depending only of y and A skew .
As a result, having set we see that the bilinear form [y, ϕ] A can be defined for all ϕ ∈ H 1 0 (Ω) using (6) and the standard rule where {ϕ ε } ε>0 ⊂ C ∞ 0 (Ω) and ϕ ε → ϕ strongly in H 1 0 (Ω). Let ε be a small parameter, I ε : U ε × Y ε → R be a cost functional, Y ε be a space of states, and U ε be a space of controls. Let be a set of all admissible pairs linked by some state equation. We consider the following constrained minimization problem: Since the sequence of constrained minimization problems (8) lives in variable spaces U ε × Y ε , we assume that there exists a Banach space U × Y with respect to which a convergence in the scale of spaces {U ε × Y ε } ε>0 is defined (for the details, we refer to [12,21]). In the sequel, we use the following notation for this convergence (u ε , y ε ) In order to study the asymptotic behavior of a family of (CMP ε ), the passage to the limit in (8) as the small parameter ε tends to zero has to be realized. Following the scheme of the direct variational convergence [12], we adopt the following definition for the convergence of minimization problems in variable spaces. Definition 1.2. A problem inf (u,y)∈Ξ I(u, y) is the variational limit of the sequence (8) as ε → 0 in symbols, inf if and only if the following conditions are satisfied: (dd) For every (u, y) ∈ Ξ ⊂ U × Y there are a constant ε 0 > 0 and a sequence {(u ε , y ε )} ε>0 (called a Γ-realizing sequence) such that 12]). Assume that the constrained minimization problem inf (u,y)∈Ξ0 is the variational limit of sequence (8) in the sense of Definition 1.2 and this problem has a nonempty set of solutions For every ε > 0, let (u 0 ε , y 0 ε ) ∈ Ξ ε be a minimizer of I ε on the corresponding set Ξ ε . If the sequence {(u 0 ε , y 0 ε )} ε>0 is relatively compact with respect to the µ-convergence in variable spaces U ε × Y ε , then there exists a pair (u 0 , y 0 ) ∈ Ξ opt 0 such that inf (u,y)∈ Ξ0 2. Setting of the optimal control problem. Let f ∈ H −1 (Ω) and y d ∈ L 2 (Ω) be given distributions. By choosing an appropriate control A ∈ L 2 (Ω; M N ). More precisely, we are concerned with the following OCP subject to the constraints

THIERRY HORSIN, PETER I. KOGUT AND OLIVIER WILK
To define the class of admissible controls A ad , , we introduce the following sets.
where A * ∈ L 2 (Ω; S N skew ) is a given matrix, c is a positive constant, Q is a nonempty convex compact subset of L 2 (Ω; S N skew ) such that the null matrix A ≡ [0] belongs to Q, and Definition 2.1. We say that a matrix A = A sym + A skew is an admissible control to the Dirichlet boundary value problem (16)- (17) We have the following result.

Proposition 1 ([9]
). The set A ad is nonempty, convex, and sequentially compact with respect to the strong topology of L 2 (Ω; M N ).
The distinguishing feature of optimal control problem (15)- (18) is the fact that the matrix-valued control A ∈ A ad is merely measurable and belongs to the space L 2 Ω; M N (rather than the space of bounded matrices L ∞ Ω; M N ). The unboundedness of the skew-symmetric part of matrix A ∈ A ad can have a reflection in non-uniqueness of weak solutions to the corresponding boundary value problem. It means that there exists a matrix A ∈ L 2 Ω; M N such that the corresponding state y ∈ H 1 0 (Ω) may be not unique. Definition 2.2. We say that (A, y) is an admissible pair to the OCP (15)- (18) if A ∈ A ad ⊂ L 2 Ω; M N , y ∈ D(A) ⊂ H 1 0 (Ω), and the pair (A, y) is related by the integral identity Ω ∇ϕ, A sym ∇y + A skew ∇y We denote by Ξ the set of all admissible pairs for the OCP (15)- (18). Let τ be the topology on the set of admissible pairs Ξ ⊂ L 2 Ω; M N × H 1 0 (Ω) which we define as the product of the strong topology of L 2 Ω; M N and the weak topology of H 1 0 (Ω). We say that a pair (A 0 , y 0 ) ∈ L 2 Ω; M N × D(A 0 ) is optimal for problem (15)- (18) if (A 0 , y 0 ) ∈ Ξ and I(A 0 , y 0 ) = inf (A,y)∈ Ξ I(A, y).
As immediately follows from (7), every weak solution y ∈ D(A) to the problem (16)-(17) satisfies the energy equality Ω A sym ∇y, ∇y where the value [y, y] A may not of constant sign for all y ∈ D(A). Hence, the energy equality (24) does not allow us to derive a reasonable a priory estimate in H 1 0 -norm for the weak solutions (see [9]).
As was shown in [9], OCP (15)-(18) is always regular, i.e. Ξ = ∅, and moreover, for each f ∈ H −1 (Ω) and y d ∈ L 2 (Ω), this problem admits at least one solution. However, the main point is that for any approximation {A * k } k∈N of the matrix A * ∈ L 2 Ω; S N skew with properties {A * k } k∈N ⊂ L ∞ (Ω; S N skew ) and A * k → A * strongly in L 2 (Ω; S N skew ), optimal solutions to the corresponding regularized OCPs associated with matrices A * k always lead in the τ -limit as k → ∞ to some admissible (but not optimal in general) solution ( A, y ) of the original OCP (15)- (18). Moreover, this limit pair can depend on the choice of the approximative sequence {A * k } k∈N . However, as follows from counter-example, given in [9], it is possible a situation when none of optimal solutions to OCP (15)- (18) can be attainable in such way. In particular, the main result of [9] says that if some optimal pair ( A, y ) ∈ L 2 (Ω; M N )× H 1 0 (Ω) to OCP (15)-(18) is attainable through the above L ∞ -approximation of matrix A * , then this pair is related by energy equality Hence, the question is what kind of approximation to OCP (15)- (18) should be applied in order to attain the other types of optimal solutions which do not hold true the energy equality (25).
3. On approximation of non-variational solutions to OCP (15)- (18). We begin this section with some auxiliary results and notions. Let A ∈ A ad be a fixed matrix and let L(A) be a subspace of H 1 0 (Ω) such that i.e., L(A) is the set of all weak solutions of the homogeneous problem − div A∇y = 0 in Ω, y = 0 on ∂Ω.
Let ε be a small parameter. Assume that the parameter ε varies within a strictly decreasing sequence of positive real numbers which converge to 0. Hereinafter in this section, for any subset E ⊂ Ω, we denote by |E| its N -dimensional Lebesgue measure L N (E).
For every ε > 0, let T ε : R → R be the truncation function defined by The following property of T ε is well known (see [10]). Let g ∈ L 2 (Ω) be an arbitrary function. Then we have: Let A * ∈ L 2 Ω; S N skew be a matrix mentioned in the control constraints (21). For a given sequence {ε > 0}, we define the cut-off operators T ε : S N skew → S N skew as follows T ε (A * ) = T ε (a * ij ) N i,j=1 for every ε > 0. We associate with such operators the following set of subdomains {Ω ε } ε>0 of Ω where Definition 3.1. We say that a matrix A * ∈ L 2 Ω; S N skew is of the F-type, if there exists a strictly decreasing sequence of positive real numbers {ε} converging to 0 such that the corresponding collection of sets {Ω ε } ε>0 , defined by (30), possesses the following properties: (i) Ω ε are open connected subsets of Ω with Lipschitz boundaries for which there exists a positive value δ > 0 such that The surface measure of the boundaries of holes Q ε = Ω \ Ω ε is small enough in the following sense: in Ω, and for each element h ∈ D(A), there is a constant c = c(h) depending on h and independent of ε such that . Thus, if A * is of the F-type, each of the sets Ω ε is locally located on one side of its Lipschitz boundary ∂Ω ε . Moreover, in this case the boundary ∂Ω ε can be divided into two parts ∂Ω ε = ∂Ω ∪ Γ ε . Observe also that if A * ∈ L ∞ Ω; S N skew then the estimate (33) is obviously true for all matrices A ∈ L 2 (Ω; M N ) such that A skew A * . Remark 1. As immediately follows from Definition 3.1, the sequence of perforated domains {Ω ε } ε>0 is monotonically expanding, i.e., Ω ε k ⊂ Ω ε k+1 for all ε k > ε k+1 , and perimeters of Q ε tend to zero as ε → 0. Moreover, because of the structure of subdomains Q ε (see (31)) and L 2 -property of the matrix A * , we have This entails the property: |Ω \ Ω ε | = o(ε 2 ) and, hence, lim ε→0 |Ω ε | = |Ω|. Besides, in view of the condition (ii) of Definition 3.1, we have Remark 2. As follows from [4], F-property of the skew-symmetric matrix A * implies the so-called strong connectedness of the sets {Ω ε } ε>0 which means the existence of extension operators P ε from H 1 0 (Ω ε ; ∂Ω) to H 1 0 (Ω) such that, for some positive constant C independent of ε, Remark 3. It is easy to see that in view of the conditions (1)-(ii) of Definition 3.1 and the Sobolev Trace Theorem [1], for all ε > 0 small enough, the inequality holds true with a constant C = C(Ω) independent of ε.
As a direct consequence of Definition 3.1, we have the following obvious result.
(37) (37) and (35) = where Here, y d ∈ L 2 (Ω) and f ε ∈ L 2 (Ω) are given functions, ν is the outward normal unit vector at Γ ε to Ω ε , v ∈ H − 1 2 (Γ ε ) is considered as a fictitious control, and σ is a positive number such that Using the fact that A ∈ L ∞ (Ω ε ; M N ) for every ε > 0 and each A ∈ A ε ad , we arrive at the following obvious result.
In order to study the asymptotic behavior of the sequences of admissible solu- in the scale of variable spaces, we adopt the following concept.
and sup We are now in a position to state the main result of this section.
Theorem 3.5. Assume that the matrix A * ∈ L 2 Ω; S N skew is of the F-type. Let {Ω ε } ε>0 be a sequence of perforated subdomains of Ω associated with matrix A * . Let f ∈ H −1 (Ω) and y d ∈ L 2 (Ω) be given distributions. Then the original optimal control problem inf (A,y)∈Ξ I(A, y) , where the sequence f ε ∈ L 2 (Ω) ε>0 is such that χ Ωε f ε → f strongly in H −1 (Ω), is variational limit of the sequence (38)-(40) as the parameter ε tends to zero.
Proof. Since each of the optimization problems inf (A,v,y)∈Ξε I ε (A, v, y) lives in the , we have to show that in this case all conditions of Definition 1.2 hold true. To do so, we divide this proof into two steps.
Step 1. We show on this step that condition (dd) of Definition 1.2 holds true. Let (A, y) ∈ Ξ be an arbitrary admissible pair to the original OCP (15)- (18). We will indicate two cases. Case 1. The set L(A), defined in (26), is a singleton. It means that h ≡ 0 is a unique solution of homogeneous problem (27); Case 2. The set L(A) is not a singleton. So, we suppose that the set L(A) is a linear subspace of H 1 0 (Ω) and it contains at least one non-trivial element of D(A) ⊂ H 1 0 (Ω). We start with the Case 2. Let h ∈ D(A) be a element of the set L(A) such that h is a non-trivial solution of homogeneous problem (27). In the sequel, the choice of element h ∈ L(A) will be specified (see (65)). Then we construct a (Γ, 0)-realizing sequence {(A ε , v ε , y ε ) ∈ Ξ ε } ε>0 in the following way: (j) A ε = A for all ε > 0. In view of definition of the set A ε ad , we obviously have that A ε ∈ A ε ad ⊂ L 2 (Ω; M N ) ε>0 is a sequence of admissible controls to the problems (38). Note that in this case the properties (43)-(46) are obviously true for the sequence {A ε } ε>0 .
where distributions w ε are such that (jjj) y ε ∈ H 1 0 (Ω ε ; ∂Ω) ε>0 is the sequence of weak solutions to the corresponding boundary value problems Since Hence, due to the Lax-Milgram lemma and the superposition principle, the sequence y ε ∈ H 1 0 (Ω ε ; ∂Ω) ε>0 is defined in a unique way and for every ε > 0 we have the following decomposition y ε = y ε,1 + y ε,2 , where y ε,1 and y ε,2 are elements of H 1 0 (Ω ε ) such that (hereinafter, we suppose that the functions y ε of H 1 0 (Ω ε , ∂Ω) are extended by operators P ε outside of Ω ε ) By the skew-symmetry property of Then (53)-(54) lead us to the energy equalities By the initial assumptions, we have h ∈ L(A). Then the condition (iii) of Definition 3.1 implies that (for the details we refer to [11]) Thus, using the continuity of the embedding H As a result, we arrive at the following the a priori estimates Hence, the sequences y ε,1 ∈ H 1 0 (Ω ε ; ∂Ω) ε>0 and y ε,2 ∈ H 1 0 (Ω ε ; ∂Ω) ε>0 are weakly compact with respect to the weak convergence in variable spaces [21], i.e., we may assume that there exists a couple of functions y 1 and y 2 in Now we can pass to the limit in the integral identities (53)-(54) as ε → 0. Using (50), (62), (57), L 2 -property of A ∈ A ad , and the fact that χ Ωε f ε → f strongly in H −1 (Ω), we finally obtain for every ϕ ∈ C ∞ 0 (Ω). Hence, y 1 and y 2 are weak solutions to the boundary value problem (16)- (17) and (27), respectively. Hence, y 2 ∈ L(A) and y 1 ∈ D(A) (see [9]). As a result, we arrive at the conclusion: the pair (A, y 1 + h) belongs to the set Ξ, for every h ∈ L(A). Since by the initial assumptions (A, y) ∈ Ξ, it follows that having set in (49) we obtain h ∈ L(A) and y ε = y ε,1 + y ε,2 y in H 1 0 (Ω ε ; ∂Ω) as ε → 0.
Therefore, in view of (66), (57), (50), we see that Thus, the property (10) holds true. It is worth to notice that in the Case 1, we can give the same conclusion, because we originally have h ≡ 0. Hence, the solutions to boundary value problems (63)-(63) are unique and, therefore, we can claim that y = y 1 , y 2 = 0, and h = 0. It remains to prove the inequality (11). To do so, it is enough to show that where the sequence {(u ε , v ε , y ε ) ∈ Ξ ε } ε>0 is defined by (49) and (65).
In view of this, we make use the following relations by (37) and (66) In order to obtain the convergence lim ε→0 Ωε we apply the energy equality which comes from the condition (A, y) ∈ Ξ and make use of the following trick. It is easy to see that the integral identity for the weak solutions y ε to boundary value problems (40) can be represented in the so-called extended form

THIERRY HORSIN, PETER I. KOGUT AND OLIVIER WILK
where h * is an arbitrary element of L. Indeed, because of the equality we have an equivalent identity to the classical definition of the weak solutions of boundary value problem (40). As follows from (57), (66), and the Sobolev Trace Theorem, the numerical sequences are bounded. Therefore, we can assume, passing to a subsequence if necessary, that there exists a value ξ 1 ∈ R such that Since y ε y weakly in H 1 0 (Ω ε ; ∂Ω) and y ∈ D(A), it follows that there exists a sequence of smooth functions {ψ ε ∈ C ∞ 0 (Ω)} ε>0 such that ψ ε → y strongly in H 1 0 (Ω). Therefore, following the extension rule (7), we have Because of the initial assumptions, we can assume that the element h * ∈ L(A) is such that So, due to this observation, we specify the choice of element h * ∈ L(A) as follows or, in other words, we aim to ensure the condition ξ 1 − ξ 2 − ξ 3 + [y, y] A = 0. As a result, we have: h * is an element of L(A) such that Having put ϕ = y ε and h * = h * in (71) and using the fact that Ω ∇y ε , A skew ∇y ε R N χ Ωε dx = 0, we arrive at the following energy equality for the boundary value problem (40) As a result, taking into account the properties (37), (66), (75), we can pass to the limit as ε → 0 in (76). This yields Hence, turning back to (67), we see that this relation is a direct consequence of (68) and (77). Thus, the sequence {(u ε , v ε , y ε ) ∈ Ξ ε } ε>0 , which is defined by (49) and (65), is Γ-realizing. The property (dd) is established.
Step 2. We prove the property (d) of Definition 1.2. Let {(A k , v k , y k )} k∈N be a sequence such that (A k , v k , y k ) ∈ Ξ ε k for some ε k → 0 as k → ∞, and the sequence of fictitious controls v k ∈ H − 1 2 (Γ ε k ) k∈N satisfies inequality (48).
In view of Definition 3.4 it means that (A k , v k , y k ) w → (A, y) as k → ∞. Our aim is to show that (A, y) ∈ Ξ and I(A, y) ≤ lim inf It is easy to see that the limit matrix A is an admissible control to OCP (15)-(18), i.e. A ∈ A ad . Since the integral identity holds true for every k ∈ N, we can pass to the limit in (80) as k → ∞ using Definition 3.4 and the estimate v k , ϕ coming from inequality (48). Then proceeding as on the Step 1, it can easily be shown that the limit pair (A, y) is admissible to OCP (15)- (18). Hence, the condition (79) 1 is valid. As for the inequality (79) 2 , we see that As a result, the lower semicontinuity of L 2 -norm with respect to the weak convergence, immediately leads us to the inequality Thus, in order to prove the inequality (79) 2 , it remains to combine relations (81), (82), and take into account the following estimate The proof is complete.
In conclusion of this section, we consider the variational properties of OCPs (38)-(40). To this end, we apply Theorem 1.3. Theorem 3.6. Let A * ∈ L 2 Ω; S N skew be a matrix of the F-type. Let y d ∈ L 2 (Ω) and f ∈ H −1 (Ω) be given distributions. Let (A 0 ε , v 0 ε , y 0 ε ) ∈ Ξ ε ε>0 be a sequence of optimal solutions to regularized problems (38)-(40), where χ Ωε f ε → f strongly in H −1 (Ω). Then there exists an optimal pair (A 0 , y 0 ) ∈ A ad to the original OCP (15)- (18), which is attainable in the following sense Proof. In order to show that this result is a direct consequence of Theorem 1.3, it is enough to establish the compactness property for the sequence of optimal solutions (A 0 ε , v 0 ε , y 0 ε ) ∈ Ξ ε ε>0 in the sense of Definition 3.4. Let h ∈ C ∞ 0 (Ω) be a non-zero function such that div (A sym ∇h + A * ∇h) ∈ L 2 (Ω), where we assume that A = A sym + A * is an admissible control, A ∈ A ad .
We set v ε = ∂h ∂ν A Γε ∈ H − 1 2 (Γ ε ). In view of the initial assumptions and estimate (see [11] for the details) there is a constant C > 0 independent of ε such that be a corresponding solution to boundary value problem (40). Then following (60), we come to the estimate where the constant C is also independent of ε. As a result, we get Since ε −σ H N −1 (Γ ε ) → 0 as ε → 0, it follows that the minimal values of the cost functional (39) bounded above uniformly with respect to ε. Thus, the sequence of optimal solutions (A 0 ε , v 0 ε , y 0 ε ) ε>0 to the problems (38)-(40) uniformly bounded in L 2 (Ω; M N ) × H − 1 2 (Γ ε ) × H 1 0 (Ω ε ) and, hence, in view of Proposition 1 , it is relatively compact with respect to the weak convergence in the sense of Definition 3.4. For the rest of proof, it remains to apply Theorem 1.3.

Remark 5.
We note that variational properties of optimal solutions, given by Theorem 3.6, do not suffice to assert that the convergence of optimal states P ε (y 0 ε ) to y 0 is strong in H 1 0 (Ω). Indeed, the convergence which comes from (84)-(85), does not imply the norm convergence in H 1 0 (Ω). At the same time, combining relation (86) with energy identities and Ω ∇y 0 , A 0 sym ∇y 0 rewritten for optimal solutions of the problems (51)-(52) and (16)- (17), respectively, we get lim It gives us another example of the product of two weakly convergent sequences that can be recovered in the limit in an explicit form. Moreover, this limit does not coincide with the product of their weak limits.
Our next remark deals with a motivation to put forward another concept of the weak solutions to the approximated boundary value problem (40) which can be viewed as a refinement of the integral identity (53). Definition 3.7. Let {Ω ε } ε>0 be a sequence of perforated subdomains of Ω associated with matrix A by the rule (30)-(31). We say that a function y ε = y ε (A, f, v) ∈ H 1 0 (Ω ε ) is a weak solution to the boundary value problem (40) for given holds true for all h ∈ L(A), ϕ ∈ C ∞ 0 (Ω), and ψ ∈ C ∞ 0 (Ω). Since for every A ∈ A ad and h ∈ D(A) the bilinear form [h, ϕ] A can be extended by continuity (see (7)) onto the entire space H 1 0 (Ω), it follows that the integral identity (88) can be rewritten as follows Hence, using the skew-symmetry property of the matrix A skew ∈ L 2 Ω; S N skew and the fact that the set L(A) is closed with respect to the strong topology of H 1 0 (Ω), we conclude: for every ε > 0 there exist an element h ε in L(A) such that the relation (89) can be reduced to the following energy equality . (90) Thus, in contrast to the "typical" energy equality to the boundary value problem (40), relation (90) includes some extra term which coming from the singular energy of the boundary value problem (16)-(17) that was originally hidden in approximated problem (40). However, in contrast to the similar functional effect for Hardy inequalities in bounded domains (see [18]), the term Ω ∇y ε , A sym ∇h ε R N dx + [h ε , y ε ] A is additive to the total energy, and, hence, its influence may correspond to the increasing or decreasing of the total energy and may even constitute the main part of it. 4. Optimality system for regularized OCPs associated with perforated domains Ω ε and its asymptotic analysis. As follows from Theorem 3.3, for each ε > 0 small enough, the optimal control problem inf (A,v,y)∈Ξε I ε (A, v, y) , where the cost functional I ε : Ξ ε → R and its domain Ξ ε ⊂ A ε ad × H − 1 2 (Γ ε ) × H 1 0 (Ω ε ; ∂Ω) are defined by (39)-(40), is a well-posed controllable system. Hence, to deduce an optimality system for this problem, we make use of the following well-know result.
Theorem 4.1 (Ioffe and Tikhomirov [6,5]). Let Y , U , and V be Banach spaces, let J : Y × U → R be a cost functional, let F : Y × U → V be a mapping, and let U ∂ be a convex subset of the space U containing more than one point. Let ( u, y) ∈ U × Y be a solution to the problem J(u, y) → inf, For each u ∈ U ∂ , let the mapping y → J(u, y) and y → F (u, y) be continuously differentiable for y ∈ O( y), where O( y) is some neighbourhood of the point y, and let Im F y ( u, y) be closed and it has a finite codimension in V . In addition, for y ∈ O( y), let the function u → J(u, y) be convex, the functional J is Gâteausdifferentiable with respect to u at the point ( u, y), and the mapping u → F (u, y) is continuous from U to Y and affine, i.e., Then there exists a pair (λ, where the Lagrange functional L is defined by equality If Im F y ( u, y) = V , then it can be assumed that λ = 1 in (91)-(92).
Taking into account the fact that the mapping is an epimorphism (see Theorem 1.1.4 in [5]), from (117) it follows that Thus, in view of (116) and (118), relation (114) takes the form . Applying the same arguments as before, we finally conclude that As a result, having gathered relations (112), (116), and (119), we arrive at the boundary value problem (105)-(106). Moreover, by the regularity of solutions to the problem (105)-(106), we have p ε ∈ H 2 (Ω ε ) ∩ H 1 0 (Ω ε ; ∂Ω) [7]. In order to end of the proof of this theorem, it remains to show the validity of the relations (107)-(108). With that in mind, we note that, in view of the structure (94)-(96), condition (92) takes the form Here, we have used the fact that H 1 2 (Γ ε ) can be reduced to a Hilbert space with respect to an appropriate equivalent norm, and, hence, H − 1 2 (Γ ε ) is a dual Hilbert space as well (for the details we refer to Lions and Magenes [15, p.35]).

Remark 7.
In view of the assumption (102), we make use of the following observation. Let {(A ε , v ε , y ε ) ∈ Ξ ε } ε>0 be a weakly convergent sequence in the sense of Definition 3.4. Since in this case y ε ∈ H 1 0 (Ω ε ; ∂Ω) ε>0 are the solutions to the boundary value problem (99)-(100) with A = A ε , and g = f ε ∈ L 2 (Ω), and w = v ε ∈ H − 1 2 (Γ ε ), it follows that the sequence div A ε ∇y ε χ Ωε ε>0 is obviously bounded in L 2 (Ω). However, because of the non-symmetry of L 2 -matrices {A ε } ε>0 , it does not imply the same property for the sequence div A skew ε ∇y ε χ Ωε ε>0 . In order to guarantee this property, we make use of the notion of divergence div A of a skew-symmetric matrix A ∈ L 2 Ω; S N skew . We define it as a vector-valued distribution d ∈ H −1 (Ω; R N ) following the rule where a i stands for the i-th column of the matrix A. As a result, we can give the following conclusion: if div A skew ε ∈ L ∞ (Ω; R N ) for all ε > 0 and the sequence div A skew ε ε>0 is uniformly bounded in L ∞ (Ω; R N ), then there exists a constant C > 0 independent of ε such that sup ε>0 χ Ωε div A skew ε ∇y ε L 2 (Ω) ≤ C. (123) for any ψ ε , ϕ ∈ C ∞ 0 (Ω ε ) (due to the fact that div A skew ε ∈ L ∞ (Ω; R N ) for all ε > 0), it follows that this relation can be extended by continuity to the following one Hence, × ∇y ε L 2 (Ωε;R N ) < +∞.
In what follows, we divide the proof onto several steps.
Step 2. On this step we study the limit passage in inequality (108) as ε → 0. To this end, we rewrite it as follows where By Theorem 3.6 (see (85)), we have by the compactness of the embedding H 1 0 (Ω) → L 2 (Ω), and lim ε→0 ε −σ v 0 = 0 by Theorem 3.6 (see estimate (83)), it follows from (136) that (138) Step 3. As for the term J ε 3 (A), we see that (138) and (137) and lim ε→0 Ωε as the limit of product of weakly and strongly convergence sequences in L 2 (Ω; R N ). Hence, combining relations (139) and (140), we get Step 4. At this step we study the asymptotic behaviour of the term J ε 1 (A) in (133) as ε → 0. To this end, we note that in view of the property (5), the lower semicontinuity of L 2 -norm with respect to the weak convergence, immediately leads us to the inequality However, because of inequality in (142), we cannot assert that the limit values are related as follows In order to guarantee this relation, we assume the converse, namely, there exists a matrix A ∈ A ad such that J 1 (A ) < J 2 − J 3 (A ). That is, in view of (138),(141), and (142), this leads us to the relation The direct computations show that, in this case, we arrive at the inequality However, this contradicts with the Lagrange principle, and therefore, the inequality (143) remains valid. Thus, following (143), we finally get for all A ∈ A ad . This concludes the proof.
Remark 8. As Theorem 4.3 indicates, the limit passage in optimality system (103)-(108) for the regularized problems (38)-(40) as ε → 0 leads to the optimality system for the original OCP (15)- (18). However, a strict substantiation of this passage requires rather strong assumptions in the form of Hypotheses (H1)-H2). At the same time, the verification of these Hypotheses becomes trivial provided and Indeed, in this case the relation (124) takes the form and it holds obviously true provided y ε y in H 1 0 (Ω), p ε → p in H 1 0 (Ω), and A skew A * ∈ L ∞ (Ω; S N skew ). Hence, Hypothesis (H1) is valid. As for Hypothesis (H2), we see that admissible controls A ∈ A ad with extra property (146) form a close set with respect to the strong convergence in L 2 (Ω; S N skew ). Moreover, in this case we have that the sequence χ Ωε div A 0 ε skew ∇y 0 ε ε>0 is uniformly bounded in L 2 (Ω) (see Remark 7). Hence, the sequence of adjoint states {p ε } ε>0 , given by (105)-(106), is bounded in H 2 (Ω ε ) by the regularity of solutions to the problem (105)-(106). Hence, within a subsequence, we can suppose that the sequence {P ε (p ε )} ε>0 is weakly convergent in H 2 (Ω). This proves Hypothesis (H2).
5. Numerical simulations. The main issue of this section is to present numerical simulations that tend to ascertain our approaches developed above. We restrict ourselves to the case when Ω is the unit ball of R 2 or R 3 . The numerical simulations have been conducted according three guidelines. For this we consider some matrix A d ∈ L 2 (Ω) N ×N and y d in H 1 0 (Ω), and set f = f d := − div (A d ∇y d ).
We focus on the following test case: with the uniform ellipticity condition on A sym given by (5). For this problem under view the algorithm used should allow to recover the pair (A d , y d ), because the minimum of (147) is clearly 0. Once validated, we return to the original OCP (1), for which we consider singular y d and A d in two manners: we still consider A d , y d and f d with A d possibly singular at some point ξ of the unit ball Ω in R 2 or R 3 . We triangulate Ω by a triangulation τ such that no vertices of τ is ξ and such that no edges of τ contains ξ.
We proceed to the classical gradient algorithm.
In this case, we expect, but cannot prove, that the algorithm converges to a variational solution. Indeed, when projecting on the grid, due to our assumption, we cannot distinguish between singular and non singular data. Moreover, for each projected matrix A in the admissible set, the projected matrix gives rise to a unique solution, thus the projected problem changes in its behavior. And of course as already said, due to the non-singular situation, we are led to think that the sequence of approximate solutions constructed will give rise to a variational solution.
In the final simulation procedure, we have punctured our domain and discretized the OCP given in (2). Accordingly, there is now no singularity in the punctured domain. We, afterwards, consider refining the punctured domain by reducing the size of the hole.
In the following sections we describe more precisely each scheme and present some numerical results with some interpretations in each case that, we do think, clarifies the situation.
Of course a specific numerical investigation should be made further on to try to discriminate better between the variationnal and non-variationnal situations. We postpone a precise numerical analysis in that scope to future works and reports.
Let us nevertheless mention that what show the following computations is that the most realistic and at least fast numerical approximation of our optimal control problems is the prevalence of non variationnal solutions and this hypothesis should be, as we just said, later investigated.

5.1.
Validation. Throughout this section and the following ones, we will take A d of the following form: In the 2d-case whilst in the 3d-case For the case of the unpunctured domain, the gradient is obtained by using the adjoint state p (see, for instance, (120) and further).
Let p be the solution of We get where W ∈ L 2 (Ω) N ×N . We adopt a finite element method for y and p such that A is constant for each triangular element of the mesh. In order for the algorithm to be more efficient, we use more data than these discrete components of A. We set n different pairs (u i d , f i d ), i = 1, n . To reduce the value of n, we choose to use a spatial smoothing for each component of G test . In order to do so, several options are possible ( [24], [25]). The new cost functional modified according to these n tests is now (with y i solution of the state equation (148) where The gradient becomes hereafter a mean of terms obtained in (152). For the two-dimensional case, we use 16 pairs (y i d , f i d ) associated to a combination of sinusoidal functions useful to capture sufficient information. Each state y i d verify the state problem with f equal to f i d and A equal to the reference A d (Figure 1). The coefficient ε 0 is equal to 10 6 . The initial matrix A is by its coefficients (A 11 (x), A 12 (x); A 21 (x), A 22 (x)) = (1, 0.2; 0.1, 1.1).
The results ( Figure 2) show a coherent convergence.
For the three-dimensional case, the simulation durations prevent to use the same level of discretization than for the two-dimensional cases. We use 48 pairs (y i d , f i d ) associated to a combination of sinusoidal functions equivalent to the 2D-cases. We use 11929 points and 72946 cells for the mesh (without hole). So we work on 72946 variables for each component of A.
The number of pairs (y i d , f i d ) and the smoothing are useful and allow us to control all theses variables, but with difficulties. We must parallelize our control problem. The n pairs (y i d , f i d ) create n different state problems, each of them can be computed on different core. We use this characteristic to reduce to a few days the simulation duration. We test our three-dimensional program with a singular asymmetric A d . The results are shown on Figures 3 and 4. The results are as consistent as for the 2D-problem.

5.2.
Discretization in the unpunctured domain. We return to the original OCP (1). We use the same pairs (y i d , f i d ), i = 1, n but now the real A d should be considered as unknown that is to say that we now consider a real optimization problem, while the preceding test cases, could be considered as an identification or inverse problem.
The figures 5 and 6 showthe results of the direct simulation in two-and threedimensional cases, respectively (without the trick with puncturing of the singularity region). We use these results to compare with the next results associated to the OCP (2). 5.3. Discretization in the punctured domain. At this step we consider the approximation of the original OCP in the form of (2). In this case, we must add thep to the adjoint p state solution of (151) for each pairs (u i d , f i d ) where Ω is replaced by Ω ε andp = − q ε σ on Γ ε , where q satisfies (denoting B ε :) q − ∆q = 0 in B ε , (here Ω = Ω ε B ε ) ∂q ∂ν = v on Γ ε .     For the two-dimensional case, the pictures 7, 8 show the results. The second case uses a smaller hole. For the three-dimensional case, the figure 9 shows the results. We can note that the values of the functional is always smaller than the cases without hole. For the second 2D case with a smaller hole, the components become more different than these obtained with the OCP (1).
Of course these results do not validate the existence of variational and nonvariational solutions. However, according to Zhikov [23], or if we believe that the uniqueness and regularity results in [2] lead to the absence of non-variational solutions in dimension 2, the numerical simulations above tends to show that arguably this does exist in dimension 2. However, due to computational performance and refinement requirement, it is probably very difficult to ascertain that our numerical simulations do prove the prevalence of non-variational solutions or not to OCP (1) on the class of admissible controls A with unremovable singularity.