Optimization and Discretization in 2D Problems of Electromagnetic Invisible Cloaking Gennady Alekseev1,2 , Aleksey Lobanov3 , and Yuliya Spivak1 1 Far Eastern Federal University, Vladivostok, Russia 2 Far Eastern State Fisheries University, Vladivostok, Russia 3 Institute of Applied Mathematics FEB RAS, Vladivostok, Russia alekseev@iam.dvo.ru Abstract. We study control problems for the 2D electromagnetic field model describing scattering TM-polarized electromagnetic waves in unbounded ho- mogeneous medium containing a penetrable inhomogeneous dielectric obstacle with the boundary partially coated for masking. These problems arise when developing the design technologies of electromagnetic cloaking devices using optimization method. Two constitutive parameters: variable refraction index of the obstacle and surface conductivity of the coated part of the boundary play the role of controls. Solvability of control problems is proved, the opti- mality system which describes the necessary conditions of extremum is derived, uniqueness and stability of optimal solutions are established. Two numerical algorithms are proposed and discussed. The first of them is based on strategy “optimize-then-discretize” and the second one is based on opposite strategy “discretize-then-optimize”. Keywords: TM-polarization · invisibility cloaking · design · regularization · discretization · control problem · optimality system · uniqueness · stability 1 Introduction. Statement of Direct Problem In recent years significant research has focused on design of devices cloaking material objects from detection of radar system. Beginning with papers [8, 11, 13] the large number of papers (see, e.g., [5, 7, 10, 12, 15]) was devoted to developing different schemes of cloaking. These schemes include metamaterial cloaking based on transformation optics (TO) proposed by Pendry et al. [13], conformal method proposed by Leonhardt [11], plasmonic cloaking method based on scattering cancellation proposed by Alú and Engheta [5], mantle cloaking, impedance cloaking, etc. Development of the above-mentioned approaches have opened up the opportunities for creation the invisibility cloaking design strategies. They obtained the name of di- rect design strategies as they were based on solving the forward electromagnetic (or acoustic) problems. It should be noted that the invisibility devices (hereafter, cloaks) designed on the basis of direct strategies possess serious drawbacks. The main one is the Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A. Kononov et al. (eds.): DOOR 2016, Vladivostok, Russia, published at http://ceur-ws.org 126 G. Alekseev et al. difficulty of their technical realization. For example, the design of the TO-based cloaks involves extreme values of constitutive parameters and spatially varying distributions of the permittivity and permeability tensors which are very difficult to implement [18]. That is why the another cloak design strategy began develop recently. It obtained the name of inverse design as it is related with solving inverse electromagnetic (or acoustic) problems. The optimization method forms the core of this inverse design methodology. This enables us to solve some substantial limitations of previous cloak- ing solutions. A growing number of papers is devoted to applying the inverse design methodology in various cloaking problems. Among them we mention [14, 16, 17] where numerical optimization algorithms are applied for finding the unknown material pa- rameters of TO-based cloak. It was shown there that the optimized multi-layer cloak essentially outperforms the similarly sized metamaterial cloak designed by using the TO approach. In [18] the authors review the invisibility cloak design methodologies and discuss the recent transition from forward design to inverse design. We also men- tion papers [1–4] where the mathematical apparatus for solving impedance cloaking problem on the basis of optimization approach is developed. This paper is devoted to theoretical analysis of control problems for 2D electromag- netic wave scattering model. These problems arise when optimization method is applied for solving cloaking problems for respective 2D electromagnetic scattering model. We begin with formulation of the direct scattering problem. Let Ω be a bounded domain in R2 with a connected complement Ω c := R2 \Ω and Lipschitz boundary Γ consisting of two parts ΓD and ΓI . We consider the scattering problem for TM-polarized electromagnetic waves in homogeneous medium containing penetrable inhomogeneous dielectric obstacle Ω with coated partially (for masking) boundary. Mathematically this problem is reduced to finding functions w in Ω and u = uinc + us in Ω c satisfying equations ∆w + k 2 n(x)w = 0 in Ω, 4u + k 2 u = 0 in Ω c , (1) mixed transmission conditions on the boundary Γ ∂w ∂u ∂w ∂u w − u|Γ = 0, − |ΓD = 0, − |Γ = iη(x)u (2) ∂ν ∂ν ∂ν ∂ν I and the Sommerfeld radiation condition in R2 √ ∂us lim r( − ikus ) = 0, r = |x|. (3) r→∞ ∂r Here uinc is the incident wave, us is the scattered wave, k is the wave number, k 2 = ε0 µ0 ω 2 where ω is an angular frequency, ε0 and µ0 are constant electric permittivity and magnetic permeability, n(x) > 0 is a variable index of refraction of dielectric obstacle Ω, η is the surface conductivity of the coated part ΓI of Γ , i is an imaginary unit, ν is the outward (relative to Ω) unit normal on Γ . One can find the formulation and brief analysis of problem (1)–(3) in [6]. Besides, in [6] the inverse scattering problem of recovering the shape and surface conductivity of a partially coated dielectric infinite cylinder from the far field data was also stud- ied. The control problems considered in our paper consist of minimization of certain cost functionals dependent on the state (electromagnetic field) and unknown functions Optimization and Discretization in 2D Problems of Invisible Cloaking 127 Fig. 1. Geometry of cloaking problem (controls or design parameters) satisfying equations (1)–(3). As the cost functional we choose one of the following: Z Z d 2 I1 (U ) = |U − u | dx, I2 (U ) = |U − ud |2 dσ. (4) Q Γr Here U is the function equalled to w in Ω and to u in Ω c , function ud ∈ L2 (Q) (or ud ∈ L2 (Γr )) describes the field measured in some subdomain Q ⊂ Ω c or on the boundary Γr of the disc Br of the radius r containing Ω inside. In the case when ud = uinc the functional I1 (U ) (or I2 (U )) has the sense of squared mean-square integral norm of the scattered field us over Q (or over Γr ). As controls we choose index of refraction n and surface conductivity η. We assume that n and η are elements of Sobolev spaces H r (Ω) and H s (ΓI ) and define the following regularized functional: α0 α1 α2 Jj (U, n, η) = Ij (U ) + knk2H r (Ω) + kηk2H s (ΓI ) . (5) 2 2 2 Here j = 1 or 2, α0 , α1 and α2 are nonnegative parameters specifying the relative importance of each term in (5). We want to find controls n, η and the associated state – electromagnetic field U = (w, u), such that the functional Jj (U, n, η) defined in (5) is minimized subject to state equations (1)–(3). The rest of the paper is organised as follows. In Section 2 we reduce unbounded problem (1)–(3) to an equivalent problem considered in bounded domain and prove the correct solvability of the latter problem. In Section 3 we prove the existence of a solution of control problem and derive an optimality system. Based on analysis of the optimality system we prove further in Section 4 uniqueness and stability of optimal solutions. Then in Section 5 we propose and discuss two numerical algorithms for solving our control problem. 2 Functional Spaces. Solvability of Direct Problem Let us introduce the function spaces to be used in the subsequent analysis. Let BR be the disc of radius R containing Ω, Ωe := Ω c ∩ BR . We will use the spaces H s (Ω), 128 G. Alekseev et al. H 1 (Ωe ), H 1 (BR ), L2 (Q), L2 (ΓI ), H 1/2 (ΓR ), H −1/2 (ΓR ), L∞ (ΓI ), H s (ΓI ) with norms k · ks,Ω , k · k1,Ωe , k · k1,BR , k · kQ , k · kΓI , k · k1/2,ΓR , k · k−1/2,ΓR , k · kL∞ (ΓI ) and k · ks,ΓI , respectively. The scalar products and norms in H r (Ω), L2 (Ω), H s (ΓI ) and L2 (ΓI ) will be denoted by (·, ·)r,Ω , k · kr,Ω , (·, ·)Ω , k · kΩ , (·, ·)s,ΓI , k · ks,ΓI and (·, ·)ΓI , k · kΓI , respectively. Along with the space H 1 (Ω) we will consider it’s subspace H 1 (∆, Ω) := {w : w ∈ H (Ω), ∆ w ∈ L2 (Ω)} equipped with the norm kwk2H 1 (∆,Ω) = kwk21,Ω + k∆wk2Ω . It 1 is well known (see [9, p. 28]) that any function w ∈ H 1 (∆, Ω) has the trace γ1 w ≡ ∂w/∂ν|Γ ∈ H −1/2 (Γ ) and the following Green formula holds: Z ∂w (∆w, w)Ω = −(∇w, ∇w)Ω + wdσ ∀w ∈ H 1 (Ω). (6) Γ ∂ν Here and below integral over Γ (or over ΓR ) denotes the duality pairing h·, ·iΓ between H 1/2 (Γ ) and H −1/2 (Γ ) (or between H 1/2 (ΓR ) and H −1/2 (ΓR )). Similar formula holds and for the domain Ωe . We also need the space X = H 1 (BR ) with the norm k · kX := k · k1,BR and the space inc H ≡ H inc (Ωe ) = {u ∈ H 1 (Ωe ) : ∆u + k 2 u = 0 in Ωe } with the norm kuk1,Ωe . These spaces will be used for describing properties of the weak solutions of problem (1)–(3) and for describing incident waves uinc , respectively. By X ∗ we denote dual of space X. Let L∞ ∞ r r n0 (Ω) = {n ∈ L (Ω) : n(x) ≥ n0 }, Hn0 (Ω) = {n ∈ H (Ω) : n(x) ≥ n0 }, ∞ ∞ s s Lη0 (ΓI ) = {η ∈ L (ΓI ) : η(x) ≥ η0 } and Hη0 (ΓI ) = {η ∈ H (ΓI ) : η(x) ≥ η0 } where n0 = const > 0, η0 = const > 0, r > 0, s > 0. These sets will serve for describing properties of index of refraction n and conductivity η. We note that continuous compact embeddings H r (Ω) ⊂ L∞ (Ω) at r > 1 and H s (ΓI ) ⊂ L∞ (ΓI ) at s > 1/2 take place (if ΓI ∈ C 1,1 ) and the following estimates hold: knkL∞ (Ω) ≤ Cr0 knkr,Ω ∀n ∈ H r (Ω), kηkL∞ (ΓI ) ≤ Cs kηks,ΓI ∀η ∈ H s (ΓI ). (7) Here Cr0 (or Cs ) is the constant dependent on r and Ω (or on s and ΓI ). Now we are in position to study problem (1)–(3). We begin with reducing problem (1)–(3) to an equivalent problem considered in the disc BR . For this purpose we define the Dirichlet-to-Neumann (DtN) operator T : H 1/2 (ΓR ) → H −1/2 (ΓR ) that maps every function g ∈ H 1/2 (ΓR ) to a function ∂ ũ/∂ν ∈ H −1/2 (ΓR ) where ũ is a solution of the exterior Dirichlet problem for the Helmholtz equation ∆ũ + k 2 ũ = 0 in Ω c \B R with condition ũ|ΓR = g. It is well known that problem (1)–(3) is equivalent to problem (1), (2) considered in the disc BR under the following boundary condition for scattered field us on ΓR : ∂us /∂ν = T us on ΓR . (8) We will refer to (1) in Ω ∪ Ωe , (2) and (8) as problem 1. Let us multiply the first equation in (1) by Φ|Ω where Φ ∈ X is a test function, integrate over Ω and apply the Green formula (6). We obtain Z Z ∂w ∇Φ · ∇w − k 2 nΦw dx =  Φ dσ. (9) Ω Γ ∂ν Here and below Φ denotes the complex conjugate of Φ. In a similar manner, we multiply the second equation in (1) by Φ|Ωe , integrate over Ωe , apply the Green formula for the Optimization and Discretization in 2D Problems of Invisible Cloaking 129 domain Ωe and add with (9). Using the boundary conditions in (2) and condition (8) for us we arrive at the following identity with respect to function U := (w, u) ∈ X: aλ (U, Φ) := a0 (U, Φ) − an (U, Φ) − aη (U, Φ) = hf, Φi ∀Φ ∈ X. (10) Here and below index λ denotes the pair (n, η), a0 , an , aη and f are sesquilinear and linear forms defined by Z Z Z 2 a0 (U, Φ) := ∇Φ · ∇U dx + (∇Φ · ∇U − k ΦU )dx − ΦT U dσ, (11) Ω Ωe ΓR Z Z an (U, Φ) = k 2 (nU, Φ)Ω := k 2 nΦU dx, aη (U, Φ) := i(ηU, Φ)ΓI := i ηΦU dσ, (12) Ω ΓI ∂uinc Z Z inc hf, Φi := − ΦT u dσ + Φ dσ, λ = (n, η). (13) ΓR ΓR ∂ν The solution U ∈ X of problem (10) is called a weak solution of problem 1. Using the embedding theorems, trace theorem and the properties of DtN operator T it is easy to derive the following estimates for forms a0 , an , aη , f : ka0 k ≤ C1 , kan k ≤ C1 knkL∞ (Ω) , kaη k ≤ C1 kηkL∞ (ΓI ) , kf kX ∗ ≤ C1 kuinc k1,Ωe . (14) Here C1 is a constant dependent on Ω, k and R. We note that the sesquilinear form aλ introduced in (10) defines operator Aλ : X→X ∗ by hAλ U, Φi = aλ (U, Φ) for all U ∈ X, Φ ∈ X and problem (10) for U ∈ X is equivalent to equation Aλ U = f. (15) Using properties of forms a0 , an , aη and operator T one can prove arguing as in [2] that the operator Aλ is an isomorphism. Denote by A−1 ∗ λ : X → X the inverse of the operator Aλ . Let C˜λ = kA−1 λ k. It follows from the results above that for any element f ∈ X ∗ equation (15) has a unique solution Uλ ∈ X which satisfies the estimate kUλ kX ≤ C̃λ kf kX ∗ . Besides, in the case when index of refraction n and conductivity η belong to nonempty bounded sets K1 ⊂ Hnr 0 (Ω), r > 1 and K2 ⊂ Hηs0 (ΓI ), s > 1/2 one can show proceeding as in [2] that the solution Uλ of (15) satisfies the estimate kUλ kX ≤ C̃0 kf kX ∗ where constant C̃0 is independent of λ. Using estimate kf kX ∗ ≤ C1 kuinc k1,Ωe from (14) and setting C0 = C̃0 C1 we rewrite this estimate as kUλ kX ≤ C0 kuinc k1,Ωe ∀λ = (n, η) ∈ K1 × K2 . (16) Let us formulate the result obtained as the next theorem. Theorem 1. Let Γ ∈ C 0,1 , ΓI ∈ C 1,1 and let K1 ⊂ Hnr 0 (Ω) and K2 ⊂ Hηs0 (ΓI ) be nonempty bounded sets where r > 1, s > 1/2. Let uinc ∈ H inc be an incident field. Then for any pair (n, η) ∈ K1 × K2 problem (10) has a unique solution Uλ ∈ X which satisfies estimate (16) with constant C0 independent of λ. 130 G. Alekseev et al. 3 Solvability of Control Problem. Optimality System In this Section we formulate and study our control problem. We assume that controls n and η can change in certain sets K1 and K2 . More precisely, the following conditions are assumed to hold: (j) Γ ∈ C 0,1 ΓI ∈ C 1,1 ; α0 > 0; K1 ⊂ Hnr 0 (Ω) and K2 ⊂ Hηs0 (ΓI ) are nonempty convex closed sets where n0 =const > 0, η0 =const > 0, r > 1, s > 1/2. Let K = K1 × K2 , λ = (n, η). Defining the operator G : X × K × H inc → X ∗ by hG(U, λ, uinc ), Φi = a0 (U, Φ) − k 2 (nU, Φ)Ω − i(ηU, Φ)ΓI − hf, Φi for all Φ ∈ X consider the constrained minimization problem α0 α1 α2 J(U, λ) := I(U ) + knk2r,Ω + kηk2s,ΓI → inf, 2 2 2 G(U, λ, uinc ) = 0, (U, λ) ∈ X × K, λ = (n, η). (17) Here I : X → R is a weakly lower semicontinuous cost functional. Denote by Zad = Zad (uinc ) := {(U, λ) ∈ X × K : G(U, λ, uinc ) = 0, J(U, λ) < ∞} the set of admissible pairs for problem (17). Theorem 2. Let under conditions (j) I : X → R be a weakly lower semicontinuous functional, uinc ∈ H inc and Zad be nonempty set. Suppose that α1 ≥ 0, α2 ≥ 0 and K is bounded set or α1 > 0, α2 > 0 and functional I is bounded below. Then problem (17) has at least one solution (U, λ) ∈ X × K. Proof. Let (Um , λm ) ∈ Zad where λm := (ηm , nm ) be minimizing sequence for which lim J(Um , λm ) = inf J(U, λ) := J ∗ . m→∞ (U,λ)∈Zad From conditions (j) and Theorem 1 the following estimates follow: knm kr,Ω ≤ c1 , kηm ks,ΓI ≤ c2 , kUm kX ≤ c3 . Here c1 , c2 , c3 are some constants which are independent of m ∈ N = {1, 2, ...}. By definition of Zad the pair (Um , λm ) satisfies the identity a0 (Um , Φ) − k 2 (nm Um , Φ)Ω − i(ηm Um , Φ)ΓI = hf, Φi ∀Φ ∈ X, m ∈ N. (18) It follows from the estimates above that there exist weak limits n∗ ∈ K1 ⊂ Hns 0 (Ω), η∗ ∈ K2 ⊂ Hηs0 (ΓI ) and U∗ ∈ X of some subsequences of the sequences {nm }, {ηm } and {Um }. Using this fact and compactness of embeddings H r (Ω) ⊂ L∞ (Ω) at r > 1, H s (ΓI ) ⊂ L∞ (ΓI ) at s>1/2 we conclude (passing if necessary to subsequen- cies) that Um → U∗ ∈ X weakly in X and Um |Ω → U∗ |Ω weakly in L2 (Ω), Um |ΓI → U∗ |ΓI weakly in L2 (ΓI ), nm → n∗ ∈ K1 strongly in L∞ (Ω), ηm → η∗ ∈ K2 strongly in L∞ (ΓI ). Let us pass to the limit in (18) when m → 0. Using (14) we obtain that the pair (U∗ , λ∗ ) where λ∗ := (n∗ , η∗ ) satisfies a0 (U∗ , Φ) − k 2 (nU∗ , Φ)Ω − i(η∗ U∗ , Φ)ΓI = hf, Φi ∀Φ ∈ X. (19) Optimization and Discretization in 2D Problems of Invisible Cloaking 131 This means that G(U∗ , λ∗ , uinc ) = 0. Since J is weakly lower semicontinuous on X × K, we have J(U∗ , λ∗ ) = J ∗ which proves the theorem. We note that the assertion of Theorem 2 is valid for both functionals I1 (U ) and I2 (U ) since they are nonnegative and are weakly lower semicontinuous. The next stage in the study of control problem (17) is to establish sufficient con- ditions on the input data under which its solution is unique and stable for particular cost functionals. For this purpose we make use the approach developed in [2, 3]. It is based on the derivation and analysis of an optimality system describing the first-order necessary conditions for an extremum in problem (17). Arguing as in [3] one can prove the following result. Theorem 3. Let under conditions (j) the triple (Û , n̂, η̂) ∈ X × K be a solution of problem (17) where the functional I(U ) is continuously differentiable with respect to U in the point Û . Then there exists a unique nonzero Lagrange multiplier P ∈ X which satisfies the Euler-Lagrange equation a0 (Ψ, P ) − k 2 (n̂Ψ, P )Ω − i(η̂Ψ, P )ΓI = −(α0 /2)hIU0 (Û ), Ψ i ∀Ψ ∈ X (20) and the minimum principle holds which is equivalent to inequalities α1 (n̂, n − n̂)r,Ω − k 2 Re((n − n̂)Û , P )Ω ≥ 0 ∀n ∈ K1 , (21) α2 (η̂, η − η̂)s,ΓI − Re[i((η − η̂)Û , P )ΓI ] ≥ 0 ∀η ∈ K2 . (22) Direct problem (10), the Euler-Lagrange equation (20) which has the sense of ad- joint problem for the adjoint state P ∈ X and variational inequalities (21), (22) con- stitute the optimality system for control problem (17). The optimality system plays an important role in studying the properties of solutions of the control problem. On its basis, efficient numerical algorithms of solving problem (17) can be developed. In addition, using analysis of the optimality system one can establish the sufficient condi- tions on the initial data providing the uniqueness and stability of solutions of particular extremal problems. 4 Uniqueness and Stability of Optimal Solutions We assume that the incident field uinc can change in a bounded set K inc ⊂ H inc . Denote by (U1 , n1 , η1 ) ∈ X × K a solution of (17) corresponding to given field uinc = uinc 1 ∈ K inc . By (U2 , n2 , η2 ) ∈ X × K we denote a solution of problem ˜ λ) = α0 I(U J(U, ˜ ) + α1 knk2r,Ω + α2 kηk2s,Γ → inf, 2 2 2 I G(U, λ, ũinc ) = 0, (U, λ) ∈ X × K, λ = (n, η). (23) It is obtained from (17) by replacing functional I(U ) by another functional I(U˜ ) and replacing incident field uinc by another incident field ũinc = uinc 2 ∈ K inc . We assume 132 G. Alekseev et al. that the set K is bounded and derive one important inequality for the difference of solutions of problems (17) and (23). We note firstly that by Theorem 1 the following estimates hold for Ul , l = 1, 2: kUl kX ≤ MU = C0 sup kuinc k1,Ωe , uinc ∈ K inc . (24) Denote by Pl ∈ X , l = 1, 2 Lagrange multipliers corresponding to solutions (Ul , nl , ηl ). By Theorem 3 Pl , l = 1, 2 satisfy identity a0 (Ψ, Pl ) − k 2 (nl Ψ, Pl )Ω − i(ηl Ψ, Pl )ΓI = −(α0 /2)hIU0 (Ul ), Ψ i ∀Ψ ∈ X. (25) Set uinc = uinc inc 1 − u2 , n = n1 − n2 , η = η1 − η2 , U = U1 − U2 , P = P1 − P2 , f = f1 − f2 . (26) We subtract the identity (10) written for U2 , n2 , η2 , uinc 2 from (10) for U1 , n1 , η1 , uinc 1 to obtain the following equation for the difference U = U1 − U2 : a0 (U, Φ) − k 2 (n2 U, Φ)Ω − i(η2 U, Φ)ΓI = = k 2 (nU1 , Φ)Ω + i(ηU1 , Φ)ΓI + hf, Φi ∀Φ ∈ X. (27) Using estimates (7), (14) and (24) we deduce that |k 2 (nU1 , Φ)Ω | ≤ C1 Cr0 knkr,Ω MU kΦkX , |(ηU1 , Φ)ΓI | ≤ C1 Cs kηks,ΓI MU kΦkX . (28) It follows from (28) and Theorem 1 applied to problem (27) for U that kU kX ≤ C0 (Cr0 MU knkr,Ω + Cs MU kηks,ΓI + kuinc k1,Ωe ), C0 := C1 C̃0 . (29) We set n = n1 in (21) written at n̂ = n2 , Û = U2 , P = P2 and then we set n = n2 in (21) written at n̂ = n1 , Û = U1 , P = P1 . Using (26) we obtain α1 (n2 , n)r,Ω − k 2 Re[(n U2 , P2 )Ω ] ≥ 0, −α1 (n1 , n)r,Ω + k 2 Re[(n U1 , P1 )Ω ] ≥ 0. Adding these inequalities we arrive at the following inequality for n, U and P : −k 2 Re[(n U, P1 )Ω + (n U2 , P )Ω ] ≤ −α1 knk2r,Ω . (30) In a similar manner we obtain the following inequality for η, U and P : −Re[i(η U, P1 )ΓI + i(η U2 , P )ΓI ] ≤ −α1 kηk2s,ΓI . (31) Let us subtract (25) at l = 2 from (25) at l = 1. Setting Ψ = U we obtain a0 (U, P ) − k 2 (n2 U, P )Ω − k 2 (nU, P1 )Ω − i(η2 U, P )ΓI − −i(ηU, P1 )ΓI = −(α0 /2)hIU0 (U1 ) − I˜U0 (U2 ), U i. (32) We set Φ = P in (27), subtract from (32) and add the real part of obtained result with (30) and (31). Using relation (nU, P1 )Ω + (nU2 , P )Ω − (nU1 , P )Ω = (nU, P1 )Ω − (nU, P )Ω = (nU, P2 )Ω and analogous one for terms in (32) containing η we arrive at the inequality (α0 /2)Re[hIU0 (U1 ) − I˜U0 (U2 ), U i] ≤ −α1 knk2r,Ω − α2 kηk2s,ΓI + +Re[k 2 (nU, P1 + P2 )Ω + i(ηU, P1 + P2 )ΓI − hf, P i]. (33) Let us formulate obtained results as the Lemma. Optimization and Discretization in 2D Problems of Invisible Cloaking 133 Lemma 1. Let in addition to conditions (j) K and K inc ⊂ H inc be bounded sets and let the triples (U1 , n1 , η1 ) and (U2 , n2 , η2 ) be solutions of problems (17) at uinc = uinc 1 ∈ K inc and (23) at ũinc = uinc 2 ∈ K inc , respectively. Let functionals I(U ) and ˜ ) be continuously differentiable and let Pl be Lagrange multipliers corresponding I(U to (Ul , nl , ηl ), l = 1, 2. Then the estimate (29) for the difference U := U1 − U2 and inequality (33) for differences defined in (26) hold. Based on Lemma 1 we are able now to study uniqueness and stability of solutions of control problem Jj (U, n, η)→ inf, G(U, λ, uinc ) = 0, (U, λ)∈X×K, λ = (n, η), j = 1, 2 (34) for particular cost functionals defined in (4), (14). We begin with the case j = 1 corresponding to functional I1 (U ) := kU − ud k2Q . Denote by (U1 , n1 , η1 ) the solution of (34) at j = 1 corresponding to given functions ud = ud1 ∈ L2 (Q) and uinc = uinc 1 ∈ K inc ⊂ H inc . By (U2 , n2 , η2 ) we denote the solution of (34) at j = 1 corresponding to perturbed functions ũd = ud2 ∈ L2 (Q) and ũinc = uinc 2 ∈ K inc . Setting ud = ud1 − ud2 in addition to relation (26) we have hI10 (Ul ), Ψ i = 2(Ul − udl , Ψ )Q , hI10 (U1 ) − I˜10 (U2 ), U )i = 2(kU k2Q − (U, ud )Q ). Then the identity (25) for Lagrange multiplier Pl ∈ X and inequality (33) for differences (26) take the form a0 (Ψ, Pl ) − k 2 (nl Ψ, Pl )Ω − i(ηl Ψ, Pl )ΓI = −α0 (Ψ, Ul − udl )Q ∀Ψ ∈ X, (35) α0 (kU k2Q − Re(U, ud )Q ) ≤ Re[k 2 (nU, P1 + P2 )Ω + +i(ηU, P1 + P2 )ΓI − hf, P i] − α1 knk2r,Ω − α2 kηk2s,ΓI . (36) Firstly we estimate multipliers P1 , P2 and the term hfi , P i entering into the right- hand side of (36). To this end we consider the problem (35) for Lagrange multiplier Pl which is equivalent to the following equation: A∗λl Pl = α0 fl ∈ X ∗ , hfl , Ψ i = −(Ψ, Ul − udl )Q , l = 1, 2. Here A∗λl is an adjoint operator of Aλl . It is defined by hA∗λl P, Φi = aλl (Φ, P ) = hAλl Φ, P i for all P ∈ X, Φ ∈ X. Since |(Ψ, Ul − udl )Q | ≤ [kUl kX + max(kud1 kQ , kud2 kQ )]kΨ kX then using the properties of adjoint operators and Theorem 1 we derive the estimate kPl kX ≤ C̃0 α0 MU0 , l = 1, 2, MU0 = MU + max(kud1 kQ , kud2 kQ ). (37) Taking into account (14), (26) and (37) we deduce that |hf, P i| ≤ C1 kuinc kX kP1 + P2 kX ≤ α0 akuinc kX , a := 2C0 MU0 . (38) Using estimates (7), (14), (29), (37) and Young’s inequality 2cd ≤ c2 + (1/)d2 for all c ≥ 0, d ≥ 0,  > 0 at  = 1 we have |k 2 (nU, P1 + P2 )Ω | ≤ C1 knkL∞ (Ω) kU kX kP1 + P2 kX ≤ ≤ 2C1 C0 C˜0 α0 MU0 Cr0 knkr,Ω (Cr0 MU knkr,Ω + Cs MU kηks,ΓI 134 G. Alekseev et al. +kuinc k1,Ωe ) ≤ α0 b(4Cr02 MU2 knk2r,Ω + Cs2 MU2 kηk2s,ΓI + kuinc k21,Ωe ), (39) |i(ηU, P1 + P2 )ΓI | ≤ C1 kηkL∞ (ΓI ) kU kX kP1 + P2 kX ≤ 0 ≤ α0 b(Cr2 MU2 knk2r,Ω + 4Cs2 MU2 kηk2s,ΓI + kuinc k21,Ωe ), b := C02 MU0 MU−1 . (40) We assume that the following conditions take place: α1 (1 − ε) > 5α0 bCr02 MU2 , α2 (1 − ε) > 5α0 bCs2 MU2 (41) where ε ∈ (0, 1) is an arbitrary constant. Using (41) we derive from (39), (40) Re|k 2 (nU, P1 + P2 )Ω + i(ηU, P1 + P2 )ΓI | ≤ ≤ α1 (1 − ε)knk2r,Ω + α2 (1 − ε)kηk2s,ΓI + 2α0 bkuinc k21,Ωe . (42) Taking into account (38) and (42) from (36) we obtain α0 kU k2Q ≤ α0 Re(U, ud )Q − εα1 knk2r,Ω − εα2 kηk2s,ΓI + α0 ϕ(kuinc k1,Ωe ). (43) Here function ϕ(·) is defined by 1/2 ϕ(kuinc k1,Ωe ) = akuinc k1,Ωe + 2bkuinc k21,Ωe (44) where constants a and b are defined in (38) and (39). Omitting nonpositive term −εα1 knk2r,Ω − εα2 kηk2s,ΓI in the right-hand side of (43) we have kU k2Q ≤ kU kQ kud kQ + akuinc k1,Ωe + 2bkuinc k21,Ωe . This is a quadratic inequality for kU kQ . Solving it for kU kQ = kU1 − U2 kQ we obtain the estimate kU1 − U2 kQ ≤ kud1 − ud2 kQ + ϕ(kuinc inc 1 − u2 k1,Ωe ). (45) If uinc 1 = uinc 2 (45) transforms to the estimate kU1 − U2 kQ ≤ kud1 − ud2 kQ . Using estimate (45) and inequality kU kQ kud kQ ≤ kU k2Q +(1/4)kud k2Q which follows from Young’s inequality we obtain from (43) that εα1 knk2r,Ω + εα2 kηk2s,ΓI ≤ α0 [(1/2)kud kQ + ϕ(kuinc k1,Ωe )]2 . (46) From (46) and (29) we deduce the estimates: p p kn1 − n2 kr,Ω ≤ α0 /εα1 ∆, kη1 − η2 ks,ΓI ≤ α0 /εα2 ∆, p p kU1 − U2 kX ≤ C0 (Cr0 MU α0 /εα1 ∆ + Cs MU α0 /εα2 ∆ + kuinc inc 1 − u2 k1,Ωe ) (47) where ∆ = (1/2)kud1 − ud2 kQ + ϕ(kuinc inc 1 − u2 k1,Ωe ). (48) The estimates (47) have the sense of stability estimates of the solution (Û , n̂, η̂) of problem (34) at j = 1 with respect to small perturbations of functions ud ∈ L2 (Q) and uinc ∈ H inc . We formulate the obtained result as Optimization and Discretization in 2D Problems of Invisible Cloaking 135 Theorem 4. Let in addition to conditions (j) K := K1 × K2 and K inc ⊂ H inc be bounded sets and let the triple (Ul , nl , ηl ) ∈ X × K be a solution of problem (34) at j = 1 corresponding to given functions udl ∈ L2 (Q) and uinc l ∈ K inc , l = 1, 2, where Q ⊂ Ωe is a nonempty open subset. Suppose that conditions (41) take place. Then the stability estimates (45) and (47) hold where ∆ is given by (48). Similar result holds and for problem (34) at j = 2 corresponding to I2 (U ). 5 Numerical Algorithms Optimality system (10), (20), (21), (22) derived above can be used to design efficient numerical algorithms for solving control problem (34). The simplest one (Algorithm 1) for I1 (U ) can be obtained by applying simple iteration method for solving the optimality system. The m-th iteration of this algorithm consists of finding unknown values U m , P m , nm+1 and η m+1 for given nm and η m by sequentially solving following problems: a0 (U m , Φ) − k 2 (nm U m , Φ)Ω − i(η m U m , Φ)ΓI = hf inc , Φi ∀Φ ∈ X, (49) a0 (Ψ, P m ) − k 2 (nm Ψ, P m )Ω − i(η m Ψ, P m )ΓI = −α0 (Ψ, U m − ud )Q ∀Ψ ∈ X, (50) m+1 m 2 m+1 m m α1 (n , n − n )r,Ω − k Re((n − n )U , P )Ω ≥ 0 ∀n ∈ K1 , (51) m+1 m m+1 m m α2 (η , η − η )s,ΓI − Re[i((η − η )U , P )ΓI ] ≥ 0 ∀η ∈ K2 . (52) For discretization and solving problems (49), (50) one can use open source software free FEM++ (www.freefem.org) based on using finite element method. For discretization of (51), (52) it is convenient to look for solutions n and η as N X M X n(x) = nj ϕj (x), x ∈ Ω, η(x) = ηk ψk (x), x ∈ Γ1 . (53) j=1 k=1 Here ϕj ∈ H r (Ω) and ψk ∈ H s (ΓI ) are nonnegative basis functions in H+ r (Ω) and s H+ (ΓI ) and nj ≥ 0 and ηk ≥ 0 are unknown coefficients. Similar algorithm which is based on the strategy “optimize-then-discretize” can be used and for functional I2 (U ). Now we discuss another algorithm (Algorithm 2) which is based on the opposite strategy: “discretize-then-optimize”. The idea of this algorithm consists of seeking un- known controls – refraction index n and surface conductivity η in the form (53). Here nj and ηk are unknown coefficients which one should define from the condition of minimum of the discrete analogue of functional I1 (U ) (or I2 (U )) in (4) which has the form Z I1 (n1 , . . . , nN , η1 , . . . , ηM ) = |U (n1 , . . . , nN , η1 , . . . , ηM ) − ud |2 dx. (54) Q Here U (n1 , . . . , nN , η1 , . . . , ηM ) is a solution of the direct problem (1)–(3) for the case when parameters n(x) and η(x) have the form (53). In such a case the discrete analogue of problem (34) takes the form N M α0 α1 X 2 α2 X 2 I1 (n1 , . . . , nN , η1 , . . . , ηM ) + nj + ηk → inf, (55) 2 2 j=1 2 k=1 136 G. Alekseev et al. where nj ≥ 0, ηk ≥ 0, j = 1, . . . , N, k = 1, . . . , M . Problem (55) represents the finite-dimensional problem of conditional minimization which can be solved numerically using known methods of solution of discrete extremum problems. The formal comparison of both algorithms shows that Algorithm 1 is more complicated and more expensive in terms of CPU time and memory space. This is due to the fact that Algorithm 1 is based on solving the optimality system (10), (20), (21), (22) which involves a coupled system of state and adjoint equation together with two variational inequalities for sought for controls. 6 Concluding Remarks We studied control problems for the 2D electromagnetic field model describing scat- tering TM-polarized electromagnetic waves by a penetrable inhomogeneous dielectric obstacle. These problems arise when optimization method is applied for solving cloaking problems for respective scattering model. The refraction index n(x) of the inhomoge- neous medium filling the obstacle and the boundary conductivity η(x) of the coated part of the boundary play the role of controls. We studied some new properties of solutions of the direct problem, proved the solvability of control problems and de- rived the optimality systems describing the necessary conditions of extremum. Based on analysis of the optimality system we established the uniqueness and stability esti- mates of optimal solutions. Besides, we proposed two numerical algorithms for solving our cloaking problems. Separate paper by the authors will be devoted to comparative study of properties of these algorithms and to detailed analysis of results of numerical experiments. Acknowledgments. The first and the third authors were supported by Russian Foun- dation for Basic Research (project no. 16-01-00365-a), the second author was supported by the Russian Science Foundation (project no. 14-11-00079). References 1. Alekseev, G.V.: Optimization in problems of material-body cloaking using the wave-flow method. Dokl. Phys. 58(4), 147–151 (2013) 2. Alekseev, G.V.: Control of boundary impedance in two-dimensional material-body cloak- ing by the wave flow method. Comp. Math. Mathem. Phys. 53(12), 1853–1869 (2013) 3. Alekseev, G.V.: Cloaking via impedance boundary condition for 2-D Helmholtz equation. Appl. Anal. 93(2), 254–268 (2014) 4. Alekseev, G.V.: Analysis and Optimization in Problems of Cloaking of Material Bodies for the Maxwell Equations. Differential Equations 52, 366–377 (2016) 5. Alú, A., Engheta, N.: Achieving transparency with plasmonic and metamaterial coatings. Phys. Rev. E. 72, 016623 (2005) 6. Cakoni, F., Colton, D., Monk, P.: The inverse electromagnetic scattering problem for a partially coated dielectric. J. Comp. Appl. Math. 204(2), 256–267 (2007) 7. Cummer, S.A., Popa, B.-I., Schurig, D., Smith, D.R., Pendry, J., Rahm, M., Starr, A.: Scattering theory derivation of a 3D acoustic cloaking shell. Phys. Rev. Lett. 100, 024301 (2008) Optimization and Discretization in 2D Problems of Invisible Cloaking 137 8. Dolin, L.S.: On a possibility of comparison of three-dimensional electromagnetic systems with nonuniform anisotropic filling. Izv. Vuzov Radiofizika 4, 964–967 (1961) 9. Girault, V., Raviart, P.A.: Finite element method for Navier–Stokes equations. Theory and algorithms. Springer-Verlag, Berlin (1986) 10. Greenleaf, A., Kurylev, Y., Lassas, M., Uhlmann, G.: Invisibility and inverse prolems. Bull. Amer. Math. Soc. 46, 55–97 (2009) 11. Leonhardt, U.: Optical conformal mapping. Science 312, 1777–1780 (2006) 12. Norris, A.N.: Acoustic cloaking theory. Proc. R. Soc. Lond. A. 464, 2411–2434 (2008) 13. Pendry, J.B., Shurig, D., Smith, D.R.: Controlling electromagnetic fields. Science 312, 1780–1782 (2006) 14. Popa, B.-I., Cummer, S.A.: Cloaking with optimized anisotropic layers. Phys. Rev. A. 79, 023806 (2009) 15. Romanov, V.G., Chirkunov, Y.A.: Nonscattering acoustic objects in an anisotropic medium of special kind. Dokl. Math. 87, 73–75 (2013) 16. Urzhumov, Y., Landy, N., Driscoll, T., Basov, D., Smith, D.R.: Thin low-loss dielectric coating for free-space cloaking. Opt. Lett. 38, 1606–1608 (2013) 17. Wang, X., Semouchkinaa, E.: A route for efficient non-resonance cloaking by using mul- tilayer dielectric coating. Appl. Phys. Lett. 102, 113506 (2013) 18. Xu, S., Wang, Y., Zhang, B., Chen, H.: Invisibility cloaks from forward design to inverse design. Science China 56, 120408 (2013)