An inverse scattering problem for the time-dependent Maxwell equations : nonlinear optimization and model-order reduction

An inverse scattering problem for the time-dependent Maxwell equations : nonlinear optimization and model-order reduction

Roberta Mancini Stefan Volkwein

Konstanzer Schriften in Mathematik Nr. 297, Februar 2012

ISSN 1430-3558

© Fachbereich Mathematik und Statistik Universität Konstanz

Fach D 197, 78457 Konstanz, Germany

Konstanzer Online-Publikations-System (KOPS) URL: http://nbn-resolving.de/urn:nbn:de:bsz:352-184632

Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nla

An Inverse Scattering Problem for the Time-Dependent Maxwell Equations: Nonlinear Optimization and Model-Order Reduction

R. Mancini^∗, S. Volkwein

University of Constance, Department of Mathematics and Statistics, Germany

SUMMARY

In this paper an inverse scattering problem for the time-dependent Maxwell curl equations is considered.

This problem is formulated as an optimal control problem governed by partial differential equations. First- order necessary optimality conditions are discussed. For the numerical solution a gradient-based algorithm is applied and successfully tested for some numerical examples. Finally, model-order reduction based on proper orthogonal decomposition is utilized to derive a reduced-order model for the time-dependent Maxwell curl equations. Its applicability is tested with respect to different input frequencies. Copyright c2010 John Wiley & Sons, Ltd.

Received . . .

KEY WORDS: Time-dependent Maxwell equations; inverse scattering; optimal control; nonlinear optimization; model reduction; proper orthogonal decomposition.

1. INTRODUCTION

Inverse scattering problems [2,18,19,30,33] are very important in noninvasive imaging. In fact, because electromagnetic (EM) waves can penetrate in various media, where objects may be placed, the solution of inverse scattering problems allows to determine the presence and the properties of these objects hidden in the medium. Hence, one should deal with Maxwell’s equations. They can be considered either in frequency or in time domain. In this paper the time domain version is used. In fact, data in time domain are richer in information because of the large spectrum of possible frequency of the incident field We focus on an unbounded domain, which arises in many applications, e.g., in mine detection [7]. Because of that, in the problem formulation we need to use an infinit radiation condition.

Optimal control problems for partial differential equations are often hard to tackle numerically so that the need for developing novel techniques emerges. One such technique is given by reduced- order methods. For a general overview we refer the reader to [1]. Recently the application of reduced-order models to optimal control problems has received an increasing amount of attention;

see, e.g., [14,25]. The reduced-order approach is based on projecting the dynamical system onto subspaces consisting of basis elements that contain characteristics of the expected solution. This is in contrast to, e.g., finit element techniques, where the elements of the subspaces are uncorrelated to the physical properties of the system that they approximate. In this work we concentrate on the method of proper orthogonal decomposition (POD) [15,35] as one of the reduced-order approaches.

Since the solution of Maxwell equations – especially in three dimensions – has a significantl computational cost, we consider a POD Galerkin scheme for the FDTD scheme. The POD method is used here to reduce the EM fields For other model order reduction technique in the EM context,

∗Correspondence to: Roberta Mancini, University of Constance, Department of Mathematics and Statistics, Universit¨atsstraße 10, D-78457 Konstanz, Germany

Figure 1. A cross section of scenario for an inverse scattering problem in the domain^Ωand scatterer domain D(^Ω. The antennas with the circles are the transmitting antennas in the incidence position, while the

antennas without circles are the receiving ones.

refer to [5,6,10,12]. In [8] POD is applied for a f nite difference discretization using a staggered grid. For an error estimate we refer to [31]. The authors consider parabolic equations and use a f nite difference scheme as a high-dimensional discretization method. An error estimate is derived for the error between the solution to a parabolic equation and the associated POD Galerkin solution.

The paper is organized in the following manner: In Section2we formulate the inverse scattering problem as a nonlinear optimal control problem and present f rst-order necessary optimality conditions. The discretization of the optimality system is introduced in Section 3, where also numerical examples are shown. Section 4 is devoted on the proper orthogonal decomposition technique. Finally, we draw some conclusions in the last section.

2. EM INVERSE SCATTERING PROBLEM IN OPTIMIZATION FRAMEWORK In this section, we formulate a time-domain electromagnetic (EM) inverse scattering problem in the framework of optimization with partial differential equations (PDEs) as constraints; [13,34]. Our objective is to minimize a functional of EM f eld data misf t including a regularization functional for the EM properties that are sought. To characterize the solution of the resulting optimization problem, we derive the corresponding f rst-order optimality system.

2.1. Problem formulation

Consider the scenario illustrated in Figure 1. This scenario consists of an unbounded domain Ω⊂R³. The boundary of Ω– if it exists – is denoted by Γ. Moreover, let T >0 stand for the terminal time andQ= (0, T)×Ω. LetD(Ωdenote thescatterer domainandQD= (0, T)×D. We suppose thatΩ\D has known EM properties, e.g., they can be the free space EM properties with electric permittivityε0>0, the magnetic permeabilityµ0>0and conductivityσ0= 0. In the scatterer domain D we assume that the EM propertiesε,µ, andσdepend on the spatial variable x∈D and that they are different from those of the background medium. Moreover, we consider multiple antennas and receivers surrounding the scatterer domainD. Theinverse scattering problem consists in the problem of reconstructing the EM properties p= (ε, µ, σ)in the scatterer domain without accessing in it, knowing the EM f eld only in some points of the domainΩ\D. Throughout the paper we always assume that the scatterers have non-dispersive and isotropic EM properties (for detailed def nitions, see [11]).

The goal is to estimate the EM properties in the scatterer domain D by utilizing the time- dependent Maxwell curl equations as a mathematical model to predict the electrical and magnetic f eld E,H:Q→R³ atkm specif ed measurement pointsx^m_k ∈Ω\D,1≤k≤km. We suppose that measurement data is available at the pointsx^m_k. Since the EM propertiesp= (ε, µ, σ)enter the Maxwell equations, we can identify these properties in D in such a way that the corresponding f elds E=E(p) and H=H(p) are as close as possible to given data at the measurement points. In particular, we excite the Maxwell curl equations byni different sources in order to get

ni f eld variables (En,Hn), 1≤n≤ni. We start with one incident wave in a certain position xⁱ₁∈Ω\D. Then, we compute the associated EM f elds(E1,H1)in thekmmeasurement points x^m_i ,1≤k≤km over the time interval[0, T]. Then, we repeat the same procedure for the other ni−1incidence points.

We introduce the parameter space

P_=H¹_(D)∩L^∞_(D),

whereH¹(D) =H¹(D)×H¹(D)×H¹(D)andL^∞(D) =L^∞(D)×L^∞(D)×L^∞(D). The set of admissible parameters by

P_ad=

(ε, µ, σ)∈P

_{εa≤ε≤εb, µa ≤µ≤µb, σa ≤σ≤σb}almost everywhere inD with positive lower and upper boundsεa,µa,σa,εb,µb,σb∈L^∞(D)satisfyingεa≤εb,µa≤µb, σa≤σbinDalmost everywhere. Then, for any parameter vectorp∈Padand for thenth incident wave (1≤n≤ni) the EM f elds(En,Hn)satisfy the Maxwell curl equations (see, e.g., [11,28])

∇ ×En+µ0

∂Hn

∂t = 0 inQ\QD,

∇ ×En+µ∂Hn

∂t = 0 inQD,

∇ ×Hn−ε0

∂En

∂t −Jn = 0 inQ\QD,

∇ ×Hn−ε∂En

∂t −σEn= 0 inQD.

(1a)

Together with (1a) we pose the initial conditions

En(0,·) = 0, Hn(0,·) = 0 inΩ (1b) and the Sommerfeld radiation condition at inf nity

r→∞lim r

∇ × En

Hn

+ 1

c0

n_r× ^∂

∂t En

Hn

= 0 inΩr=

x∈R³

_{|x|2< r ,} (1c) where k0=ω√µ0ε0= 2πf√µ0ε0 and n_r is the outward normal vector to ∂Ωr, the surface boundary ofΩr. In (1c) we use thatD(Ωrfor suff ciently large radiusr. Finally, we def ne the f eld excitation. We consider a hard source [32], i.e., we assign a time dependent excitation function to a f eld component while the density currentJnis zero:

Hn(t,xⁱ_n) =f_n(t) in(0, T). (1d) Remark 2.1

Let us mention that (1d) can be replaced byEn(t,xⁱ_n) =f˜_n(t)in(0, T). A third possibility is to apply different current densitiesJnin (1a) forn= 1, . . . , ni. 3 Assumption 1(Existence and uniqueness for the Maxwell curl equations)

For any incident pointxⁱ_n∈Ω\D,1≤n≤ni, and for anyp∈Pad there exist unique f eldsEn

andHnsatisfying

En, Hn ∈C [0, T];H¹(Ω)∩C(Ω)

∩H¹ 0, T;C(Ω)

=:Yi

forn= 1, . . . , ni, whereC(Ω) =C(Ω)×C(Ω)×C(Ω).

Remark 2.2 1) The spaceY_iis a Banach space endowed with the norm kGkY_i= max

t∈[0,T]

kG(t)kH¹(Ω)+kG(t)kC(Ω)

+

Z 1

0 kG(t)k²C(Ω)+kG(t)˙ k²C(Ω)dt 1/2

forG∈Y_i. We def ne the Banach spaces Y=Y_i×. . .×Y_i

| {z }

2n_i−times

, X=Y_×P

endowed with their natural product topology. Moreover letX_ad=Y_×P_ad⊂X.

2) For bounded domains, perfectly conducting boundary conditions and suitable smooth data it is shown in [20] that the Maxwell curl equations admits a unique solution.

3) If Assumption1holds true, we can def ne the solution operator_S:P_ad→Y, where (E,H) =

(En)ⁿ_n=1ⁱ ,(Hn)ⁿ_n=1ⁱ

=S(p), p= (ε, µ, σ)∈P_ad,

denotes the solution pairs to (1) for theni incident waves. Due to the three bilinear terms µ^∂H_∂tⁿ,ε^∂E_∂tⁿ andσEnthe operatorSis nonlinear. 3 To write the inverse scattering problem as an optimization problem, we def ne a cost functional that penalizes the discrepancies between the measured data and the calculated f elds. We aim at minimizing the misf t by varying the unknown EM properties inDwhich def ne the optimization variables. Because of the ill-posedness of the inverse problem, we include a regularization term, R(p). Thus we introduce the cost functionalJ :X_→Rby

J(X) = 1 2

ni

X

n=1 km

X

k=1

Z T

0 |Enk(t)−E^m_nk(t)|²2+η₀²|Hnk(t)−H^m_nk(t)|²2

dt+R(p), (2) whereX = (E,H,p)∈X,E= (E1, . . . ,En_i),H= (H1, . . . ,Hn_i),_{| · |}2denotes the Euclidean norm,η0>0is the characteristic impedance of free space,Enk=En(·,x^m_k)andHnk(·,x^m_k)are the computed evaluated in thekth measurement point for thenth incident wave f elds, whereasE^m_nk andH^m_nkare the associated measurement data. We consider the following regularization term

R(p) =β 2

Z

D

_∇ε(x) ²₊

_∇σ(x) ²₊

_∇µ(x) ²_dx +γ

2 Z

D

_{ε(x)−ε(x)¯} ²₊

_{σ(x)−σ(x)¯} ²₊

_{µ(x)−µ(x)¯} ²_dx

(3)

forp= (ε, µ, σ)∈P. In this regularization functional we distinguish two different terms. The f rst term represents a f rst-order Tikhonov regularization scheme. Minimizing the functionalJ with this type of regularization means that we search the model parameterspin a smooth space with spatially slow varying functions. A result of this choice is that we obtain the edges of the reconstructed objects smeared out by diffusion. The other regularization term represents a zero-order Tikhonov regularization and corresponds to the requirement thatεis as close as possible to a given nominal

¯

ε, and similarlyσis as close as possible toσ¯, andµis as close as possible toµ¯. This regularization term can be used even if we do not have a priori information with knownp^¯= (¯ε,µ,¯ σ)¯ . In this case we choose these values to be zero and the resulting regularization term penalizes large parameters values.

Now, the inverse scattering problem is expressed as the inf nite-dimensional, nonlinear optimization problem

minJ(X) subject to (s.t.) the pair^X= (E,H,p)∈X_adand(E,H) =S(p). (P) If Assumption1is satisf ed, we can introduce thereduced cost functional_J^ˆ:P_ad→Rby

Jˆ(p) =J(S(p),p), p∈P_ad. Then, (P) can be equivalently expressed as thereduced problem

min ˆJ(p) s.t. p∈Pad. (^Pˆ)

Assumption 2(Existence of optimal solutions)

Problem (P) admits at least one local optimal solution X^∗= (E^∗,H^∗,p^∗)∈X_ad. Moreover, p^∗= (ε^∗, µ^∗, σ^∗)is an interior point ofP_ad, i.e., εa< ε^∗< εb, µa< µ^∗< µb and σa < σ^∗< σb

inDalmost everywhere.

Remark 2.3 1) Let Assumption2hold. Then,p^∗solves (^Pˆ).

2) Suppose that Assumption2is satisf ed. Then,p^∗is a local solution to

min ˆJ(p) s.t. p∈P_. (^Pˆeq) In this case we can neglect the inequality constraints for the parameterp. Together with the associated state(E^∗,H^∗) =S(p^∗)the tripleX^∗= (E^∗,H^∗,p^∗)is a local solutionp^∗to the equality constrained problem

minJ(X) s.t. the pairX= (E,H,p)∈Xand(E,H) =S(p). (^Peq) In the numerical experiments carried out in Section3.2we apply a gradient-type algorithm to

(^Pˆeq) respectively (^Peq). 3

2.2. First-order necessary optimality conditions

Problem (^Pˆeq) is an inf nite-dimensional, constrained optimization problem. Thus, we apply techniques from PDE-optimization [13,16,34]. For this purpose we introduce the following Hilbert spaces

Z_i=L²_{(0, T;H}¹_(Ω))∩H¹_{(0, T;L}²_(Ω)), Z₌Z_{i×. . .×}Z_i

| {z }

2n_i-times

with the natural product topology and def ne the Lagrange functionalL:X_×Z_→Rby L(X,Λ) =J(X) +

ni

X

n=1

Z T 0

Z

Ω

λ^E_n ·

∇ ×Hn−ε∂En

∂t −σEn−Jn

dxdt

+

ni

X

n=1

Z T 0

Z

Ω

λ^H_n ·

∇ ×En+µ∂Hn

∂t

dxdt,

(4)

where X= (E,H,p)∈X, ^Λ= (λ^E,λ^H)∈Z with λE= (λ^E₁, . . . ,λ^E_ni), λ^H = (λ^H₁ , . . . ,λ^H_ni) andλ^E_i ,λ^H_i ∈Z_ifor1≤i≤ni.

To solve the constrained minimization problem numerically, we want to make use of the f rst- order optimality conditions. If one can ensure the existence of Lagrange multipliers [22], an optimal solutionX^∗ to (^Peq) can be characterized by a stationary point of the Lagrange functional. Let Assumptions 1and 2hold. By X^∗= (E^∗,H^∗,p^∗)∈X_ad we denote a local optimal solution to (^Peq). Moreover, letY^∗= (E^∗,H^∗) =S(p^∗)∈Y.

Assumption 3(First-order necessary optimality conditions)

The Lagrangian L is continuously Fr´echet-differentiable. Moreover, exists an associated pair of Lagrange multipliers ^Λ∗= (λ^E,∗,λ^H,∗)∈Z satisfying the following f rst-order necessary optimality conditions:

∂L

∂Y(X^∗,Λ^∗)Y = 0 for all^Y = (E,H)∈Y_, (5a)

∂L

∂p(X^∗,Λ^∗)p= 0 for allp∈P_, (5b)

where _∂Y^∂L, ^∂L_∂p denote the Fr´echet derivatives of the Lagrangian _L with respect to Y, and p, respectively.

Remark 2.4

The existence of Lagrange multipliers can be ensured by constraint qualif cation conditions; see [16,22,34], for instance. Due to (1c) the proof is not evident in our case. 3

With Assumption3holding, we proceed formally by computing the directional derivatives of the Lagrange functional. For the details we refer the reader to the thesis [21].

From (5a) we infer that the pair^Λ= (λ^E,∗,λ^H,∗)∈Zsatisf es the following adjoint system

∇ ×λ^E,∗_n −µ0

∂λ^H,∗_n

∂t = 0 inQ\QD,

∇ ×λ^E,∗_n −µ∂λ^H,∗_n

∂t +η²₀

km

X

k=1

(H^∗_n−H^m_n)δx^m_k = 0 inQD,

∇ ×λ^H,∗_n +ε0∂λ^E,∗_n

∂t = 0 inQ\QD,

∇ ×λ^H,∗_n +ε∂λ^E,∗_n

∂t −σλ^E,∗_n +

km

X

k=1

(E^∗_n−E^m_n)δx^m_k = 0 inQD, λ^E,∗_n (T,·) = 0, λ^H,∗_n (T,·) = 0 inΩ,

r→∞lim r ∇ × λ^E,∗_n

λ^H,∗_n

+ 1 c0

n_r× ^∂

∂t λ^E,∗_n

λ^H,∗_n

!

= 0 inΩr=

x∈R³

_{|x|2< r} (6)

for every incidencen∈ {1, . . . , ni}. In (6) we denote byδx^m_k :C(D)→Rthe Dirac delta mapping satisfyingδx^m_k(x) = 1ifx=x^m_k holds. Otherwise, we haveδx^m_k(x) = 0.

Let us have a closer look to system (6). Lagrange multipliers satisfy Maxwell-like equations, in which the signs of the time derivatives are opposite with respect to the standard curl equations.

Moreover, we see that they satisfy terminal conditions and particular boundary conditions, in which there is a change in the time derivative. It can be shown that Lagrange multipliersλ^E,∗_n andλ^H,∗_n are special kind of EM f elds [21,29]. In particular, they run backward in time, i.e., fromt=T to t= 0. We remark also the fact that these particular f elds are excited by the discrepancies between the measurements and the predictions of the model in the km measurement points. If the model could exactly f t the scenario, i.e. the predictions of the model f t the measurements, the adjoint f elds would have zero sources.

From (5b) we deduce the following equations in the dual spaceH¹(D)^′ofH¹(D)

∂L

∂ε(X^∗,Λ^∗) =−β∆ε^∗+γ(ε^∗−ε)¯ −

ni

X

n=1

Z T 0

∂E^∗_n

∂t ·λ^E,∗_n dt= 0,

∂L

∂µ(X^∗,Λ^∗) =−β∆µ^∗+γ(µ^∗−µ) +¯

ni

X

n=1

Z T 0

∂H^∗_n

∂t ·λ^H,∗_n dt= 0,

∂L

∂σ(X^∗,Λ^∗) =−β∆σ^∗+γ(σ^∗−σ)¯ −

n_i

X

n=1

Z T 0

E^∗_n·λ^E,∗_n dt= 0.

It can be shown [13] that the gradientJ^ˆ′(p)of the reduced cost functionalJ^ˆat a given parameter p∈H¹(D)is given by the following functional:

Jˆ^′(p) =

−β∆ε+γ(ε−ε)¯ −

ni

X

n=1

Z T 0

∂En

∂t ·λ^E_ndt,

−β∆µ+γ(µ−µ) +¯

ni

X

n=1

Z T 0

∂Hn

∂t ·λ^H_n dt,

−β∆σ+γ(σ−σ)¯ −

n_i

X

n=1

Z T 0

En·λ^E_n dt

:H¹(D)→R³,

(7)

where(E,H) =S(p)and^Λ∈Zsolve (6) for1≤n≤ni. Hence, we use these three equations in our gradient-based optimization method for updating the EM properties in the scatterer domainD.

r r r r

r r r r

r r r r

r r r r

6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6

- - - - -

- - - - -

- - - - -

- - - - -

r 6

- Ey|i−1/2,j

Ex|i,j−1/2

Hz|i,j

Figure 2. Spatial staggered grid arrangement for the TEzwave.

For the update we have to solve the state equation (1) for the given parameterpand we have to determine the solution^Λto the adjoint equations.

3. DISCRETIZATION OF THE OPTIMALITY SYSTEM

There are several techniques for discretizing Maxwell equations (consider [17] for having a detailed overview of the most used techniques). We considered the f nite difference time domain (FDTD) algorithm. For details about the method, we refer the reader to [9,32,36]. For convergence results we refer the reader to [23], where the authors consider the Maxwell curl equtaions with simpler – compared to the Sommerfeld radiation condition (1c) – boundary conditions.

3.1. The FDTD scheme

The discussion in Section 2is valid in the general case of a three dimensional domain. Now we restrict ourselves to the two dimensional case, i.e., we do not consider variation in thez-direction (see [11,32]). In this context, one can have two different types of waves:transverse electric(TEz) andtransverse magnetic(TMz). We focus on the TEz case, which has the f elds componentsEx, EyandHz:En = (Ex, Ey,0)^⊤andHn = (0,0, Hz)^⊤. For every incident waven, they are coupled accordingly to the following system

∂Ex

∂t = 1 ε0

∂Hz

∂y −Jsx

inQ\QD, ∂Ex

∂t = 1 ε

∂Hz

∂y −

Jsx+σEx

inQD,

∂Ey

∂t = 1 ε0

−∂Hz

∂x −Jsy

inQ\QD, ∂Ey

∂t = 1 ε

−∂Hz

∂x −

Jsy+σEy

inQD,

∂Hz

∂t = 1 µ0

∂Ex

∂y −∂Ey

∂x

inQ\QD, ∂Hz

∂t = 1 µ

∂Ex

∂y −∂Ey

∂x

inQD, where we skipped the indexn.

We consider the FDTD method for discretizing the TEz system. The FDTD scheme is characterized by the fact that the f eld components are considered in a staggered grid in space and in time. In our case, the spatial arrangement of the f elds component is like in Figure2.

This implies that if we consider a space grid of dimension M×N for Hz, there will be M ×(N+ 1) points forEx and(M+ 1)×N forEy. For simplicity, we consider the case of an uniform grid, so we have the step sizes∆xand∆y onxandy directions, respectively. Moreover, magnetic and electric f elds are interleaved in time. Calling∆tthe time step andNtthe number of time-steps, we evaluateExand Ey in the time instants(m+ 1/2)∆tand Hz in the time instants m∆t, withm= 0, . . . , Nt. It is worth noting that time step∆tand step sizes∆xand∆yare related by the Courant condition [9,32]. Considering a square space grid, in whichh= ∆x= ∆y, if one choose the time step ∆t then ∆ should be taken in order to fulf ll c√

2∆t < h, wherec is the speed in the vacuum. On the space-time grid, we considerEx(i∆x,(j−1/2)∆y; (m+ 1/2)∆t)≈ E_x,i,j−1/2^m+1/2 ,Ey((i−1/2)∆x, j∆y; (m+ 1/2)∆t)≈E_y,i−1/2,j^m+1/2 andHz(i∆x, j∆y;m∆t)≈H_z,i,j^m+1,

s 6

- Ey,I,J

Ex,I,J

Hz,I,J

Figure 3. Elementary cell for the TEzcase.

fori= 1, . . . , Mandj= 1, . . . , N,m= 0, . . . , Nt. Identifying the elementary cell by the couple of index(I, J)as in Figure3, the discrete system is given by

E_x,I,J^m+1/2= 2εx,I,J−∆tσx,I,J

2εx,I,J+ ∆tσx,I,J

E_x,I,J^m−1/2+ 2∆t 2εx,I,J+ ∆tσx,I,J

_Hm

z,I,J−H_z,I,J−1^m

∆x −J_sx,I,J^m

, I= 1, . . . , M, J = 2, . . . , N E_y,I,J^m+1/2= 2εy,I,J−∆tσy,I,J

2εy,I,J+ ∆tσy,I,J

E_y,I,J^m−1/2+ 2∆t 2εy,I,J+ ∆tσy,I,J

_Hm

z,I−1,J−H_z,I,J^m

∆x −J_sy,I,Jⁿ

, I= 2, . . . , M, J = 1, . . . , N, H_z,I,J^m+1=H_z,I,J^m + ∆t

µI,J

E_y,I,J^m+1/2−E_y,I+1,J^m+1/2

∆x −E_x,I,J+1^m+1/2 −E_x,I,J^m+1/2

∆y

! ,

I= 1, . . . , M, J = 1, . . . , N m= 0, . . . , Nt.

Clearly, these equations provide an update for the inner grid points. This means that we need other equations for updating the boundary values ofExandEy. Since we consider an unbounded domain, we update the f elds in these points using some absorbing boundary condition (ABC) equations. A very popular ABC in the FDTD community is the perfectly matched layer (PML) [3]. In this approach, one surround the computational domain with a certain number of cells of an anisotropic, lossy medium, in which there are σx, σy, σ^∗_x and σ^∗_y, where the latter 2 are the magnetic conductivities. Because of the anisotropy, in the PML region one need to split the theHz

f eld component, i.e., one need to considerHzxand Hzy. Since next we would like to consider a model reduction technique, we consider another ABC, in order to avoid the reduction of the split f eld components. We have used the second order Mur scheme [32], that follows from the Engquist- Majda one-way wave equations which are

∂Ex

∂t =c∂Ex

∂y for the lower bound,

∂Ex

∂t =−c∂Ex

∂y for the upper bound,

∂Ey

∂t =c∂Ey

∂x for the left bound,

∂Ey

∂t =−c∂Ey

∂x for the right bound,

wherec >0is the speed of the light in the vacuum. Using central differences, we get the second- order Mur absorbing boundary conditions for theExcomponent [24]

E_x,I,1^m+1/2=−E_x,I,2^m−3/2+c∆t−∆y c∆t+ ∆y

E_x,I,2^m+1/2+E_x,I,1^m−3/2

+ 2∆y c∆t+ ∆y

E_x,I,1^m−1/2+E_x,I,2^m−1/2 + (c∆t)²

2∆y(c∆t+ ∆y)

E_x,I+1,1^m−1/2 −2E_x,I,1^m−1/2+E_x,I−1,1^m−1/2 +E_x,I+1,2^m−1/2 −2E_x,I,2^m−1/2+E^m−1/2_x,I−1,2 forI= 2, . . . , M−1andm= 0, . . . , Nt,

E_x,I,N^m+1/2₊₁=−E_x,I,N^m−3/2+c∆t−∆y c∆t+ ∆y

E^m+1/2_x,I,N +E_x,I,N^m−3/2₊₁

+ 2∆y c∆t+ ∆y

E_x,I,N^m−1/2+E_x,I,N^m−1/2₊₁ + (c∆t)²

2∆y(c∆t+ ∆y)

E_x,I+1,N^m−1/2 −2E_x,I,N^m−1/2+E_x,I−1,N^m−1/2 + (c∆t)²

2∆y(c∆t+ ∆y)

E_x,I+1,N+1^m−1/2 −2E^m−1/2_x,I,N+1+E_x,I−1,N+1^m−1/2 fori= 2, . . . , M−1andm= 0, . . . , Nt,

E_x,1,1^m+1/2=E^m−1/2_x,1,2 +c∆t−∆y c∆t+ ∆y

E^m+1/2_x,1,2 −E_x,1,1^m−1/2 , E_x,M,1^m+1/2=E^m−1/2_x,M,2 +c∆t−∆y

c∆t+ ∆y

E^m+1/2_x,M,2 −E_x,M,1^m−1/2 , E_x,1,N^m+1/2₊₁=E_x,1,N^m−1/2+c∆t−∆y

c∆t+ ∆y

E_x,1,N^m+1/2−E_x,1,N^m−1/2₊₁ , E_x,M,N^m+1/2₊₁=E_x,M,N^m−1/2+c∆t−∆y

c∆t+ ∆y

E_x,M,N^m+1/2−E_x,M,N+1^m−1/2 , form= 0, . . . , Nt.

Analogously, we get the equations for the f eld componentEy.

Hence, we have all the equations we need for updating the f elds components in the computational grid: this means that we have a discrete scheme for solving the direct problem. We can combine the update equations in a compact matrix notation. For that purpose we def ne the vectors^Em+1/2x ∈ R^M(N+1),^Em+1/2y ∈R^N(M+1),^Hm+1_{z ∈}R^{N M} form= 0, . . . , Nt

E^m+1/2_x =

E_x,1,1^m+1/2, . . . , E_x,M,1^m+1/2, E_x,1,2^m+1/2, . . . , E_x,M,N^m+1/2, E_x,1,N+1^m+1/2 , . . . , E_x,M,N^m+1/2₊₁⊤

, E^m+1/2_y = E_y,1,1^m+1/2, . . . , E_y,M+1,1^m+1/2 , E_y,1,2^m+1/2, . . . , E_y,M^m+1/2_+1,N−1, E^m+1/2_y,1,N , . . . , E_y,M^m+1/2_+1,N⊤

, H^m+1_z = H_z,1,1^m+1, . . . , H_z,M,1^m+1, H_z,1,2^m+1, . . . , H_z,M,N−1^m+1 , H_z,1,N^m+1, . . . , H_z,M,N^m+1 ⊤

.

Then, the new iterates^Em+1/2x ,^Em+1/2y and^Hm+1z are computed as follows: Solve the linear system M_x 0

0 M_y

E^m+1/2_x E^m+1/2_y

!

= A_E

x 0

0 A_E

y

E^m−1/2_x E^m−1/2_y

! +

A_yH^m_z +J^m_sx A_xH^m

z +J^m

sy

(8a) and set

H^m+1_z =H^m_z +B_yE^m+1/2_x −B_xE^m+1/2_y , (8b) where^Mx,^My,^AEx,^AEy,^Ax,^Ay,^Bx,^Byare matrices that contain the discretization coeff cients.

Next we consider now the discretization of the adjoint problem. Recall that the adjoint f elds λ^E_n and λ^H_n are ingoing waves, with terminal conditions instead of initial conditions, i.e., they backpropagate in time starting from a f nal condition. In principle, we need a discretization scheme that runs with a decreasing time variable, in order to take into account this characteristic. In practice

we can avoid this and using the same FDTD used for solving the direct problem, after performing the transformation^˜t=T−tand thus compute the modif ed Lagrange multipliers that are solution of the following system of equations

∇ טλ^E_n +µ0

∂λ˜^H_n

∂t = 0 inQ\QD,

∇ טλ^E_n +µ∂˜λ^H_n

∂t +η₀²

km

X

k=1

H˜n−H˜^m_n

δx^m_k = 0 inQD

∇ טλ^H_n +ε0

∂λ˜^E_n

∂t = 0 inQ\QD,

∇ טλ^H_n +ε∂λ˜^E_n

∂t −σλ˜^E_n +

km

X

k=1

E˜n−E˜^m_n

δx^m_k = 0 inQD.

In fact, with the change of the time coordinate, the procedure used to compute ^˜λ^E_n andλ^˜H_n goes forward in time. At the end of the FDTD run for the modif ed Lagrange f elds, the adjoint f eldsλ^E_n andλ^H_n can be obtained reading the results from the f nal one to the f rst.

After solving forward and adjoint problem, we can approximate numerically the gradient components of the reduced cost functionalJ^ˆ(p)described in equations (7). The time derivatives are calculated using central differences approximation in time. We computed the time integration of the product terms using the trapezoidal rule. The terms related to the f rst-order regularization require the computation of a Laplace operator∆. This operator is approximated using the f ve-point stencil rule. Having the gradient components, we can use them in an iterative scheme which updates the EM properties.

We have considered a steepest descent method. Hence, in the minimization problem of the functional _J^ˆ, one can f nd a descent directiondi asdi=−Jˆ^′(p_i)∈H¹(D)^′ for every iteration i of the optimization scheme. Then, the update of the optimization variables is performed by p_i+1=p_i+sid˜i, wheresi is a step length and d^˜i= ( ˜d^ε_i,d˜^µ_i,d˜^σ_i)∈H¹(D)is the weak solution to the Neumann boundary value problem

−∆d˜i+d˜i=diinD, ∂d˜i

∂n ^{= 0}on∂D.

We considered also the BFGS for f nding at every step a new descent direction, but we did not observe remarkable improvements in the speed of the reconstruction algorithm.

The line search method gives us the step lengthsi, i.e., it tells how much to move on the descent directiondi. We considered a backtracking procedure that stopped once the Armijo condition was satisf ed (see, e.g., [26]).

3.2. A numerical example

In this subsection we consider a numerical example. We focus on the reconstruction of a dielectric scatterer. The algorithm is valid for multiple scatterers and/or magnetic and conductive scatterers, too.In our experiment, the computational domain is a square whose side-length is three times larger than that of the scatterer square domain D. We assume the side-length of D corresponding to a length of20cells. The number of time steps is set toNt= 200and the time-step size is∆t= 10⁻⁶ seconds. The time domain is(0, T)whereT =Nt∆t seconds. The spatial increment,his equal in thex- andy-direction and it is set equal toh= 2c∆t= 5.99·10²meters, such that it satisf es the Courant limit [32]. The EM properties of the background are those one of the free space i.e., µ0= 4π·10⁻⁷H/m andε0= 8.854212·10⁻¹²F/m. The measurement points and the source points are placed in the region surrounding the domainD. We consider four different incidence points and eight measurement points, i.e.ni= 4andkm= 8; see Figure4. To generate the incidence f elds, we

q q q

q q q

q q

q q

q q

Rx Rx Rx

Rx Rx Rx

Rx Rx

T x T x

T x T x

Figure 4. Picture of the position of transmitters (^{T x}) and receivers (^Rx) in the domain. Two generic scatterers are depicted in the domain^D.

have used hard sources [32]. This means that we set a variation of the magnetic f eld directly where the transmitters are placed. Hence, we need to specifyf_n(t)in equation (1d) that in the TEzcase is Hz=f(t), in the excitation points. As in [29], the excitation of the magnetic f eld was taken equal to the following

Hz(t) = ( ₁

η0

3

10sin ^2πt_T1

+ sin ^2πt_T2

, 0≤t≤T1,

0 T1< t,

for everyni. Hence, we considered a truncated composition of sine wavefunctions with frequency of20kHz and50kHz, respectively.

To evaluate our inversion algorithm, we have considered an hidden dielectric square with side- length of ten cells and relative electrical permittivityεr= 3. This square is placed such that the left lower corner has coordinate(25,25)cell units in the computational domain. This square is depicted in Figure5with a eye-bird view. Hence, we solve the associated direct scattering problem and collect

Figure 5. Simulated scenario (left plot) and reconstructed scatterer after 20 iterations (right plot) the predictions at the receivers for all time-steps in (0, T). We also built synthetic measurements:

we take the f elds obtained using the scatterer domain which we would like to reconstruct and we added noise (thus avoiding inverse crimes). We denote the measurements corresponding to thenth incidence f eld and thekth receiver asE^m_nkandH^m_nk. We have used a white Gaussian noise such that the signal-to-noise ratio is SNR= 10dB. For the two regularization terms we takeβ= 10⁻¹⁰ andγ= 10⁻¹⁰. This choice results from our experience.

To discuss the convergence properties of the algorithm, we introduce a measure of the reconstruction error as follows

ERε= vu ut

PL

l=1 ε^(l)−ε˜^(l)2

PL

l=1 ε0−ε˜^(l)2 ,

whereL= 400is the number of cells of the scatterer domain D and the tilde indicates the true value of the property. Of course, the same formula is valid for the reconstruction of magnetic and conductive objects. In Figure6the decay of the reconstructed error along the iterations is depicted.

Figure 6. Reconstruction error along the iterations.

In the right plot of Figure5, the reconstructed scatterer is presented. After20iterations, we obtain a good imaging of the true object. However, the reconstructed values are not constant in the square.

In fact the edges tend to have greater values than those in the interior of the square. Instead of the true value2.656·10⁻¹¹F/m, we have as maximum value2.695·10⁻¹¹F/m and in the interior we have a minimum value of2.161·10⁻¹¹F/m. This means that we have a maximumεr= 3.044and a minimumεr= 2.440. That is, the optimization algorithm is more accurate in correspondence of the edges.

4. A MODEL REDUCTION APPROACH: PROPER ORTHOGONAL DECOMPOSITION In this section we utilize the POD method to derive a reduced-order model for the FDTD scheme.

The performance of the reduced-order approach is studied for different frequencies.

4.1. The POD Galerkin scheme

First we describe the POD basis construction. Dealing with a TEz wave, we have three f eld components, hence we build three basis set. First, we need to have the so-calledsnapshots matrices.

To do so, we run an EM simulation using the FDTD method and we store the f elds components in three matrices:

E_x= E^1/2_x , . . . ,E^N_x^t+1/2

∈R^{M(N+1)×(Nt+1)}, E_y= E^1/2_y , . . . ,E^N_y^t+1/2

∈R^{N(M+1)×(Nt+1)}, H_z= H¹

z, . . . ,H^Nt+1

z

∈R^{M N×(Nt+1)}.

Because of the staggered grid, the three matrices have different number of rows. Let dx denote the rank of^Ex. Then,dx≤min(M(N+ 1), Nt+ 1). Analogously,dy = rankE_y≤min(N(M+ 1), Nt+ 1)anddz= rankH_z≤min(N M, Nt+ 1).

Let us focus on the f eld component^Ex. The procedure is analogous for^Eyand^Hz. In the POD approach, one would like to represent the columns^Ejx,1≤j≤Nt+ 1, of the matrix^Ex, the so- called snapshotsby ℓx≤dx vectorsψi∈R^M(N+1),i= 1, . . . , ℓx, in an optimal way. The basis

vectorsψiare given by the following optimization problem:

ψ1,...,ψminℓ∈R^dx Nt+1

X

j=1

αj

^Ejx−

ℓx

X

i=1

hE^j_x, ψiiWψi

2 W

s.t.hψi, ψjiW=δijfori, j= 1, . . . , ℓx,

whereh·,·i^W denotes a weighted inner product with a symmetric, positive def nite matrix ^W∈ R^{M(N+1)×M(N+1)}. We need to introduce this new inner product because the EM f elds are functions, hence theL²inner product is the natural product that we need. Introducing an opportune W matrix, we can approximate numerically this product. Moreover, we introduced weights αj, j= 1, . . . , Nt+ 1 for the time dependency of the EM f elds. More precisely, the αj are the trapezoidal weights for the numerical integration over[0, T].

The problem means that one would like to minimize the difference, for every snapshot ^Ejx, between the snapshot itself and its projection in the subspace def ned by the basis vectorsψi, under the constraint of orthonormality property of the basis. Let us def ne^D= diag (α1, . . . , αN_t+1).

Using the Lagrangian framework and considering the ψi vectors as optimization variables, the optimality system is anM(N+ 1)×M(N+ 1)eigenvalue problem [15,35]

YY^⊤ψi=λiψi, i= 1, . . . , ℓx,

with ^Y=W^1/2E_xD^1/2. We assume that the (real and nonnegative) eigenvalues are ordered as follows:λ1≥λ2≥. . .≥λd_x > λd_x+1=. . .=λM(N+1)= 0.

Remark 4.1 1) According to the singular values decomposition properties we can solve the equivalent(Nt+ 1)×(Nt+ 1)eigenvalue problem:

Y^⊤Yvi=λivi, i= 1, . . . , ℓx, Then, theψi’s are given by

ψi= 1

√λi

Yvi, i= 1, . . . , ℓx.

2) If{ψi}^ℓi=1^x denotes the computed POD basis of rankℓx, then we have the error formula

N_t+1

X

j=1

αj

^Ejx−

ℓ_x

X

i=1

hE^j_x, ψiiWψi

2 W=

d_x

X

i=ℓ+1

λi. Thus, if the eigenvalues of^YY⊤decay rapidly, thenPdx

i=ℓ+1λiis small. 3

Summarizing, for every f eld components we get a basis:Ψx∈R^M(N+1)×ℓx,Ψy ∈R^N(M+1)×ℓy, Φz∈R^{N M×ℓz} for the f eldsEx,Ey,Hz, respectively.

The POD bases are used in the full discrete model (8) for a Galerkin ansatz:

E^m+1/2

x ≈ΨxP^m+1/2

x , E^m+1/2

y ≈ΨyP^m+1/2

y , H^m

z ≈ΦzP^m

z,

where^Px∈R^ℓx,^Py∈R^ℓy and^Pz ∈R^ℓz have to be determined. We apply a Galerkin projection projection for the full discrete model. Instead of (8) we arrive at

M^ℓ_x 0 0 M^ℓ_y

P^m+1/2_x P^m+1/2_y

!

= _Aℓ

Ex 0 0 A^ℓ

Ey

P^m−1/2_x P^m−1/2_y

! +

A^ℓ_yP^m_z +W^ℓ_xJ^m_sx A^ℓ_xP^m_z +W_y^ℓJ^m_sy

and set (9a)

P^m+1_z =P^m_z +B^ℓ_yP^m+1/2_x −B^ℓ_xP^m+1/2_y , (9b)

where

M^ℓ_x= Ψ^⊤_xWM_xΨx, A^ℓ_E

x = Ψ^⊤_xWA_E

xΨx, A^ℓ_x= Ψ^⊤_xWA_xΨx, M^ℓ

y= Ψ^⊤_yWM_yΨy, A^ℓ

Ey = Ψ^⊤_yWA_E

yΨy, A^ℓ

y = Ψ^⊤_yWA_yΨy,

B^ℓ_y= Ψ^⊤_zWB_yΨx, B^ℓ_x= Ψ^⊤_zWB_xΨy, W^ℓ_x= Ψ^⊤_xW, W^ℓ_y= Ψ^⊤_yW,

(10)

Notice that all matrices in (10) can be computedoff-line, i.e., before we start the loop overm= 1 toNtin (9). Moreover, we have used in (9b) thatΨ^⊤_zWΨzis the identity matrix.

To measure the error in the reduced-order model we introduce the quantity

Err^Ex(ℓ) =

Nt

P

m=0

αmkE^m+1/2_x −ΨxP^m+1/2_x k^2W

Nt

P

m=0

αmkE^m+1/2_x k^2W ,

whereαm’s are the trapezoidal weights for the interval[0, T]. The quantitiesErr^Ey(ℓ)andErr^Hz(ℓ) are def ned in an analogous way.

Remark 4.2

We mention that one can use a different ansatz for building POD basis as in, e.g. [4]. Hence, one compute the time average of the snapshots of a f eld component and subtract it from the corresponding snapshots matrix. We tested this approach too but its performance are similar to

the approach presented in the section. 3

4.2. Numerical experiments

In this subsection we present numerical results from different test runs, where we study the performance of our POD Galerkin scheme. The amount of energy captured by anℓ-rank basis is [15,27,35]

E(ℓ) = Xℓ

i=1

λi

, _d X

i=1

λi,

where theλi’s are the eigenvalues satisfyingλ1≥λ2≥. . .≥λd>0. According to this, one can choose the percentagepof energy that the basis should keep and then calculate the minimum amount of vector basis required for that, using the formula

ℓ=argminn

E(ℓ) :E(ℓ)≥ p 100

o .

We study now the behavior of the POD scheme using two different kind of source waveform and different frequencies. In the following, the medium used for getting the snapshots is the free space.

First, we consider the reduction of two different sinusoidal waves, with different frequencies.

We choose for the f rst wavef1=5 kHz and for the second onef2=20 kHz. We consider these 2 simulations on the same space-time grid, so we def ne the grid parameter such that we have stability in the case of the highest frequency,f2.

Hence, usingp= 99.99%, the results are plotted in TableI. In both cases, we can compute the f eld variable sinef1= 5kHz sinef2= 20kHz

Ex ℓx= 6 ℓx= 22

Ey ℓy= 6 ℓy = 19

Hz ℓz= 9 ℓz= 24

Table I. Comparison of the number^ℓof POD basis needed in the case of sinusoidal waves with different frequencies (^{p= 99.99%}).

f eld variable sinef1= 5kHz sinef2= 20kHz Ex ErrEx(6) = 3·10⁻³ ErrEx(22) = 9·10⁻³ Ey Err^Ey(6) = 3·10⁻³ Err^Ey(19) = 4·10⁻³ Hz Err^Hz(9) = 1·10⁻² Err^Hz(24) = 1·10⁻²

Table II. Comparison of the POD errors in the case of sinusoidal waves with different frequencies (^p= 99.99%) using the number of basis vectors in TableI.

POD-error in the reconstruction of the three f elds components, see TableII. We can conclude that for keeping the same amount of energy in the higher frequency case, we need more basis vectors than in the lower case. For every f eld component, in the case off2, we need more than three times the basis number used in thef1case.

In the next simulation we considered two Gaussian pulses: the generic expression is stated in the following equation

s(t) = exp

−(t−t0)² τ²

,

wheret0= 4.5τis the time shift used for avoiding starting the waveform int= 0with the maximum value of the pulse andτ determines the width of the pulse (for more details see [9]). In this case, we do not have a monochromatic source, i.e., we do not have a single frequency. Anyway, we can identify the maximum frequency of the pulse. Hence, we are interested in checking what happens if we choose different maximum frequencies. We use f1max=5 kHz andf2max=10 kHz and we usedf2max for having the right time step and space step in order to satisfy stability and sampling conditions. The number of POD basis functions are given in Table III. In TableIVwe report on

f eld variable Gaussf1max= 5kHz Gaussf2max= 10kHz

Ex ℓx= 2 ℓx= 7

Ey ℓy = 2 ℓy = 7

Hz ℓz = 3 ℓz= 8

Table III. Comparison of the number^ℓof POD basis needed in the case of Gaussian pulse waveforms with different maximum frequencies (^{p= 99.99%})

the POD error forp= 99.99%. As before, we see that if we have a higher maximum frequency, we f eld variable Gaussf1max= 5kHz Gaussf2max= 10kHz

Ex ErrEx(2) = 1·10⁻³ ErrEx(7) = 1·10⁻³ Ey ErrEy(2) = 1·10⁻³ ErrEy(7) = 2·10⁻³ Hz ErrHz(3) = 4·10⁻³ ErrHz(8) = 1·10⁻²

Table IV. Comparison of the POD errors in the case of Gaussian pulse waveforms with different maximum frequencies (^{p= 99.99%}) using the number of basis vectors in TableIII.

need more basis functions in order to keep the same amount of energy as the lower frequency case.

Moreover, the number of POD basis is related to the type of waveform used to having the snapshots.

Even having the same frequencies (compare the casesf1for the sinusoidal case andf1maxfor the Gaussian pulse), we need less basis functions in the Gaussian case. Hence, the number of basis functions is inf uenced by the kind of waveform and its frequency.

5. CONCLUSIONS

In this paper the problem of determining the EM properties(ε(x), µ(x), σ(x))of a hidden objects has be considered. This parameter identif cation problem is formulated in terms of a nonlinear