S. IMPERIALE Geoffrey Beck
email: geoffrey.beck@ensta-paristech.fr

Sebastien Imperiale
email: sebastien.imperiale@inria.fr

Patrick Joly
email: patrick.joly@inria.fr

Patrick Joly Poems

MATHEMATICAL MODELLING OF MULTI CONDUCTOR CABLES

Keywords: 2010 Mathematics Subject Classification. 35L05, 35A35, 73R05, 35A40 telegrapher models, Maxwell equations, Asymptotic analysis

1. Introduction. The present paper is the continuation of the article [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF]. We aim at proposing a rigorous justification of simplified 1D models for the propagation of electromagnetic waves in thin conductor cables (also called transmission lines in the literature [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF]). Our approach consists in deriving such models from an asymptotic analysis of 3D Maxwell's equations, considering the transverse dimensions of the cable as the small parameter in the analysis. In this sense, our approach is similar to what has been done in mechanics for deriving the theory of beams from the equations of 3D elasticity [START_REF] Aganovic | A justification of the one-dimensional model of elastic beam[END_REF], [START_REF] Veiga | Asymptotic method applied to a beam with a variable cross section[END_REF]. It differs by the domain of application (electromagnetism) and the fact that we are considering evolution problems. Doing so, we justify and extend the telegrapher's models classically used by electrical engineers [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF], [START_REF] Weeks | Calculation of Coefficients of Capacitance of Multiconductor Transmission Lines in the Presence of a Dielectric Interface[END_REF], [START_REF] Weeks | Multiconductor transmission line theory in the TEM approximation[END_REF].

In [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF], we treated the case of a co-axial cable containing a single metallic wire. In the present paper, we essentially extend and complete the results of [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF] to the multiwires case. Note that, as in [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF], we shall restrict ourselves to the formal derivation of the model. The corresponding error analysis is postponed to a future work. Note however that a first step has been done in this direction in [START_REF] Imperiale | Error estimates for 1D asymptotic models in coaxial cables with non-homogeneous cross-section[END_REF].

We must emphasize that the situation we consider in this paper is rather general: the cable has a variable cross section of arbitrary geometry and the dielectric medium in which the electromagnetic field is propagating is heterogeneous and lossy. This degree of generality is mandatory in the application context that we have in mind, namely the non destructive testing of electric networks: the presence of defects or localized damages modifies the nice structure of a perfect cable. From the mathematical point of view, the absence of any cylindrical structure prevents us from using some tools such has Fourier transform in space or modal decomposition (separation of variables), which is done in many textbooks [6], [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF], [START_REF] Stratton | Electromagnetic Theory[END_REF].

The outline of this article is as follows. Section 2 is devoted to preliminary material about 2D vector fields. The problem under consideration is presented in detail in section 3. In section 4, we perform the formal asymptotic analysis of this problem and establish, using the results of section 2, our generalized telegrapher's model (90). Finally, in section 5, various properties of this model are established.

We intend this paper to put the theory of transmission lines into a fresh perspective that will be of interest to both physicists and mathematicians. Future work will be devoted to completing the mathematical analysis of our models (error estimates) and developing numerical methods for their resolution, as well as the corresponding numerical analysis.

2. Spaces of quasi-static vector fields and their properties.

2.1. Preliminaries and notation. In this section we establish or recall useful results about 2D vector fields. In particular, these results extend, in another context, those of [START_REF] Bermúdez | Mathematical Models and Numerical Simulation in Electromagnetism[END_REF], chapters 7 and 9. In what follows, all manifolds that will be introduced will be supposed at least Lipschitz continuous. The proof of some results sometimes requires additional regularity that will not be explicitly mentioned. We are however convinced, that our results are still valid after removing these regularity assumptions.

We define a 2D domain S in the plane x = (x 1 , x 2 ), with canonical basis (e 1 , e 2 ), that has the following structure

S = O \ N j O j , O k ∩ O = ∅ for k = , (1)
where O is a bounded simply connected open set and the O j 's are simply connected open subsets of O (holes) whose closures do not intersect each other. As a consequence, the boundary of S is the union of an exterior boundary and N interior boundaries, each of them being closed curves:

∂S = ∂S e N j=1
Σ j with ∂S e = ∂O and Σ j = ∂O j .

(

) 2
Along ∂S we shall denote n = (n 1 , n 2 ) the unit normal vector which is outgoing with respect to S. Finally, we assume that there exists N cuts Γ j such that (i) each Γ j is a (possibly curved) line joining Σ j to ∂S e , (ii

) denoting Γ = N j=1 Γ j , S Γ := S \ Γ is simply connected MATHEMATICAL MODELLING OF MULTI CONDUCTOR CABLES 3
In the sequel, we shall say that:

A 2D domain S satisfying (1), ( 2), (i) and (ii) belongs to the class C N .

In what follows, each Γ j will be parameterized by a curvilinear abscissa τ ∈ [0, L j ] (L j is the length of Γ j ) that orients (when τ increases) Γ j from ∂S e towards Σ j . Denoting also τ the corresponding unit tangent vector along Γ j , we shall define n along Γ j as the unit normal vector to Γ j such that (τ, n) is a direct basis of the plane. All these assumptions and notation are summarized by figure 1. In the sequel 2D vector fields will be implicitly identified to 3D vector fields whose third component vanish and scalar fields to 3D vector fields which are normal to the plane (x 1 , x 2 ). Therefore, for any

v = (v 1 , v 2 ) ∈ C 2 , one can write n × v = -v × n = n 1 v 2 -n 2 v 1 , e 3 × v = -v × e 3 = (-v 2 , v 1 ), (4)
so that (e 3 × v)

• n = v × n.
(5) The divergence and scalar rotational of a 2D vector field v = (v 1 , v 2 ) are defined by (using an obvious notation for partial derivatives)

div v = ∂ 1 v 1 + ∂ 2 v 2 , rot v = ∂ 1 v 2 -∂ 2 v 1 . (6)
and one has the obvious relationships

div v = rot (e 3 × v), rot v = div (v × e 3 ). (7)
Analogously, one defines the gradient and vector rotational of a scalar field ϕ as

∇ϕ = (∂ 1 ϕ, ∂ 2 ϕ), rot ϕ = (∂ 2 ϕ, -∂ 1 ϕ), (8)
in such a way that

∇ϕ = -rot ϕ × e 3 , rot ϕ = -e 3 × ∇ϕ, (9)
as well as rot (∇ϕ) = 0, div (rot ϕ) = 0.

Given ϕ ∈ C 0 S Γ , we shall denote the jump of ϕ, denoted [ϕ], along Γ as:

∀ x ∈ Γ j , [ϕ](x) := lim η 0 u(x + η n) -u(x -η n) (11)
We shall use the following Green's formula, valid for (ϕ, v) .

∈ C 1 S Γ × C 1 S Γ 2 : SΓ div v ϕ = - SΓ v • ∇ϕ + ∂S v • n ϕ + Γ [ϕ v] • n (12) SΓ rot v ϕ = SΓ v • rot ϕ - ∂S v × n ϕ - Γ [ϕ v] × n (13)
Remark 1. Both formulas [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF][START_REF] Stratton | Electromagnetic Theory[END_REF] are in fact deduced the one from the other thanks to [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF]. Moreover, it is well known that formulas [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF][START_REF] Stratton | Electromagnetic Theory[END_REF] can appropriately be extended to H(div, S) × H 1 (Γ) and H(rot, S) × H 1 (Γ) respectively, integrals over Γ being changed into duality brackets.

2.2. Quasi-static vector fields. We shall say that ρ(x) : S → C is in the class N (S) if it is bounded and its real part is strictly positive:

a. e. x ∈ S, 0 < ρ -≤ Re ρ(x). (14)
Note that in ρ ∈ N (S), ρ -1 ( ≡ 1/ρ ) ∈ N (S). Denoting

V(S) = v ∈ L 2 (S) 2 / rot v = 0}, (15)
we introduce the two spaces of quasi-static vector fields respectively of electric and magnetic type (which are known respectively as the harmonic Dirichlet and Neumann fields when ρ is constant):

E(ρ; S) = v ∈ V(S) / div(ρ v) = 0, v × n = 0 on ∂S}, (16)
H(ρ; S) = v ∈ V(S) / div(ρ v) = 0, v • n = 0 on ∂S}. (17)
By adapting, to the dimension 2 and to the case where the function ρ is not constant, the proof of the propositions 3.14 and 3.18 of [START_REF] Amrouche | Vector potentials in three-dimensional nonsmooth domains[END_REF] (see also [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations, Theory and Algorithms[END_REF], section 3.1, theorem 3.1 for an analogous result), one shows that E(ρ; S) and H(ρ; S) are N -dimensional spaces and more precisely that

E(ρ; S) = span ∇ϕ j (ρ, S), 1 ≤ j ≤ N , H(ρ; S) = span ∇ψ j (ρ, S), 1 ≤ j ≤ N ,
where the notation ∇ is defined in Remark 2 and the ϕ j (ρ, S) ∈ H 1 (S) are the solutions of the boundary value problems:

         div ρ ∇ϕ j (ρ, S) = 0, in S, ϕ j (ρ, S) = 0, on ∂S e , ϕ j (ρ, S) = δ jk , on Σ k , k ≤ N, (18)
while the ψ j (ρ, S) ∈ H 1 (S Γ ) are the solutions of

         div ρ ∇ψ j (ρ, S) = 0, in S Γ , ∇ψ j (ρ, S)•n = 0, on ∂S, ψ j (ρ, S) = δ jk , ρ ∂ n ψ j (ρ, S) = 0 on Γ k , k ≤ N. (19)
The well-posedness of (18) and (19) follows from Lax-Milgram's lemma. The solution of ( 19) is unique only up to an additive constant (its gradient is thus unique).

Remark 2. For any ψ ∈ H 1 (S Γ ) , the field ∇ψ is defined by ∇ψ ∈ L 2 (S) 2 and ∇ψ = ∇ψ in D (S Γ ).

One defines analogously rot ψ.

Note that if ψ j (ρ, S) does depend on the position of the cuts Γ k , ∇ψ j (ρ, S) does not. This is not so obvious but will become clear later thanks to Lemma 2.4, (27).

Remark 3. Referring to the weak formulation of (18) and ( 19) respectively, one sees that:

       S ρ ∇ϕ j (ρ, S) • ∇ϕ * = 0, ∀ ϕ * ∈ H 1 0 (S), S ρ ∇ ψ j (ρ, S) • ∇ψ * = 0, ∀ ψ * ∈ H 1 (S). ( 20
)
It is easy to check that the fields ∇ϕ j (ρ, S) on one hand, and similarly ∇ψ j (ρ, S) on the other hand, are linearly independent:

• If α j ∇ϕ j (ρ, S) = 0, by connectedness α j ϕ j (ρ, S) is constant in S and this constant is 0 since each ϕ j (ρ, S) vanishes on ∂S e . Then writing α j ϕ j (ρ, S) = 0 on Σ k yields α k = 0.

• If α j ∇ψ j (ρ, S) = 0, by connectedness α j ψ j (ρ, S) is constant in S. In particular α j [ψ j (ρ, S)] = 0 on Γ k , which yields α k = 0.

Moreover, these bases satisfy particular generalized bi-orthogonality relations:

Lemma 2.1. Let (ρ 1 , ρ 2 ) ∈ N (S) 2 , for any (k, ) ∈ {1, ..., N } 2 , S ∇ϕ (ρ 1 , S) • rot ψ k (ρ 2 , S) ≡ - S rotϕ (ρ 1 , S) • ∇ψ k (ρ 2 , S) = δ k . (21)
Proof. We drop the argument S for the sake of simplicity. Let us apply Green's formula [START_REF] Stratton | Electromagnetic Theory[END_REF] with ϕ = ψ k (ρ 2 ) and v = ∇ϕ (ρ 1 ). Using [START_REF] Imperiale | Error estimates for 1D asymptotic models in coaxial cables with non-homogeneous cross-section[END_REF] and the fact that, since ϕ (ρ 1 ) is constant along ∂S e and each Σ j , v × n vanishes on ∂S, we get

S ∇ϕ (ρ 1 ) • rot ψ k (ρ 2 ) ≡ SΓ ∇ϕ (ρ 1 ) • rot ψ k (ρ 2 ) = - N m=1 Γm ψ k (ρ 2 ) ∇ϕ (ρ 1 ) × n = - Γ k ∇ϕ (ρ 1 ) × n.
We conclude after noticing that -

Γ k ∇ϕ (ρ 1 ) × n = ϕ (ρ 1 )| Γ k -ϕ (ρ 1 )| ∂Se .
We introduce the "generalized Gramm N × N matrices"

Y(ρ, S) = Y ij (ρ, S) ij , Y ij (ρ, S) = S ρ ∇ϕ i (ρ, S) • ∇ϕ j (ρ, S), (22)
Z(ρ, S) = Z ij (ρ, S) ij , Z ij (ρ, S) = S ρ ∇ψ i (ρ, S) • ∇ψ j (ρ, S). ( 23
)
Lemma 2.2. The matrices Y(ρ, S) and Z(ρ, S) are invertible and symmetric. They are real symmetric positive definite if ρ is real.

Proof. We give the proof for Y(ρ, S). The proof for Z(ρ, S) is similar.

To show the symmetry, we apply (20) with ϕ * = ϕ j (ρ, S)ϕ j (ρ, S) to observe that

Y ij (ρ, S) = S ρ ∇ϕ i (ρ, S) • ∇ϕ j (ρ, S). ( 24
)
For the invertibility, let u ∈ C N and ϕ ρ (u) := N j=1 u j ϕ j (ρ, S)

∀ (u, v) ∈ C N × C N , Z(ρ, S) u, v C N = S ρ ∇ϕ ρ (u) • ∇ϕ ρ (v). (25)
In particular

Z(ρ, S) u, u C N = 0 =⇒ S (Re ρ) |∇ϕ ρ (u)| 2 = 0 =⇒ ϕ ρ (u) = 0 =⇒ u = 0.
This also shows that when Im ρ = 0, Z(ρ) is real symmetric positive definite.

Let us consider the linear map in L 2 (S) 2

J ρ : v -→ ρ e 3 × v (26)
which is an isomorphism with inverse J -1 ρ = -J 1/ρ .

Lemma 2.3. J ρ is an isomorphism from E(ρ ; S) into H(ρ -1 ; S) and from H(ρ ; S) into E(ρ -1 ; S).

Proof. If v ∈ E(ρ; S) then, using ( 7)

div ρ -1 J ρ v = div (e 3 × v) = -rot (v) = 0, rot J ρ v = rot (e 3 × ρv) = div (ρv) = 0, J ρ v • n = ρ (e 3 × v) • n = ρ v × n = 0.
which means that J ρ v ∈ H(ρ -1 ; S). Since the two spaces have dimension N , J ρ is an isomorphism. A similar argument stands for the second statement of the lemma.

Lemma 2.4. One has the formulas

∇ψ j (ρ, S) = -ρ -1 N =1 Z j (ρ, S) rot ϕ (ρ -1 , S), (27)
∇ϕ j (ρ, S) = ρ -1 N =1 Y j (ρ, S) rot ψ (ρ -1 , S), (28)
and moreover Y(ρ, S) -1 = Z(ρ -1 , S), and thus Z(ρ, S) -1 = Y(ρ -1 , S).

Proof. For simplicity, we note E(ρ) for E(ρ; S), ϕ j (ρ) for ϕ j (ρ, S) and so on.

According to lemma 2.3, J ρ ∇ψ j (ρ) ∈ E(ρ -1 ) which implies, using ( 9)

ρ ∇ψ j (ρ) = N =1 α j rot ϕ (ρ -1 ). Therefore Z jk (ρ) = S ρ ∇ψ j (ρ) • ∇ψ k (ρ) = N =1 α j S rot ϕ (ρ -1 ) • ∇ψ k (ρ).
Using (21) with ρ 1 = ρ -1 and ρ 2 = ρ shows that α jk = -Z jk (ρ) and thus (27). ( 28) is shown analogously.

Taking ρ -1 e 3 × (27), once (27) has been written for ρ -1 instead of ρ, leads to

ρ -1 rot ψ (ρ -1 ) = N k=1 Z k (ρ -1 ) ∇ϕ k (ρ). (29)
Substituting (29) into (28) leads to

∇ϕ j (ρ) = N k=1 N =1 Y j (ρ) Z k (ρ -1 ) ∇ϕ k (ρ). ( 30
)
By linear independence of the ∇ϕ j (ρ), this shows Y(ρ) Z(ρ -1 ) = Id.

Our next properties concern some invariance properties of the matrices Y(ρ, S) and Z(ρ, S). In the sequel T c denotes a conformal mapping, defined here as a diffeomorphism of R 2 whose Jacobian matrix DT c (x) is proportional to a unitary matrix at each point:

T c ∈ C 1 (R 2 , R 2 ) / ∀ x ∈ R 2 , DT c (x) t DT c (x) = α 2 (x) Id. (31)
Note that conformal mappings include linear transformations such as rotations, symmetries, homotheties and products of them.

Given ρ ∈ N (S), one defines ρ Tc ∈ N (S Tc ), where S Tc := T c (S) as

ρ Tc = ρ • T -1 c . (32)
Lemma 2.5. For any conformal mapping T , one has

Y(ρ Tc , S Tc ) = Y(ρ, S), Z(ρ Tc , S Tc ) = Z(ρ, S). (33)
Proof. Let us give the proof of the first equality, the second being deduced by lemma 2.4. From ϕ j (ρ, S) ∈ H 1 (S) solution of (18), we define ϕ Tc,j (ρ, S) ∈ H 1 (S Tc ) by

ϕ j,Tc (ρ, S) = ϕ j (ρ, S) • T -1 c
. Using (31), one sees that ϕ j,Tc (ρ, S) solves (the details are left to the reader)

         div ρ Tc ∇ϕ j,Tc (ρ, S) = 0, in S Tc , ϕ j,Tc (ρ, S) = 0, on ∂S e Tc := T c (∂S e ), ϕ j,Tc (ρ, S) = δ jk , on Σ k,Tc := T c (Σ k ), k ≤ N,
which means that ϕ j,Tc (ρ, S) = ϕ j (ρ Tc , S Tc ) namely

y = T c (x) ⇒ ϕ j (ρ Tc , S Tc ) (y) = ϕ j (ρ, S) (x).
This yields

y = T c (x) ⇒ ∇ϕ j (ρ Tc , S Tc ) (y) = DT c (x) -1 ∇ϕ j (ρ, S) (x).
Therefore, using the change of variable y = T c (x), with Jacobian α 2 (x) because we are in R 2 , we have

Y ij (ρ Tc , S Tc ) = S Tc ρ τ ∇ϕ i (ρ Tc , S Tc ) • ∇ϕ j (ρ Tc , S Tc ) dy = S Tc ρ DT t c DT c -1 ∇ϕ i (ρ, S) • ∇ϕ j (ρ, S) α 2 dx, which yields Y ij (ρ Tc , S Tc ) = Y ij (ρ, S) since α 2 DT t c DT c -1 = Id.
Finally, let us give a property of the matrices Y(ρ, S) and Z(ρ, S) when ρ is real.

Lemma 2.6. If ρ is real valued, then the matrix Y(ρ, S) is a strict M-matrix, more precisely

Y ij (ρ, S) < 0 for i = j, Y ii (ρ, S) > j =i |Y ij (ρ, S)|, ∀ 1 ≤ i ≤ N, (34)
and the matrix Z(ρ, S) satisfies:

Z ij (ρ, S) > 0, ∀ 1 ≤ i, j ≤ N. (35)
Proof. (35) is obtained from (34), lemma 2.4 (last line) and well-known properties of inverses of M-matrices. (34) appears as a consequence of the maximum principle applied to (18), which implies that 0 ≤ ϕ j (S, ρ) ≤ 1. In particular, ϕ j (S, ρ) reaches its minimum on S along each Σ i , for i = j. Then, ∂ n ϕ j (S, ρ) ≤ 0 along Σ i . Using Green's formula we deduce

Y ij (ρ, S) = Σi ρ ∂ n ϕ j (S, ρ) ≤ 0. ( 36
)
In fact, the above inequality is strict since Y ij (ρ, S) = 0 would imply ∂ n ϕ j (S, ρ) = 0 along Σ i . As ϕ j (S, ρ) = 0 along Σ i , by unique continuation, one would deduce that ϕ j (S, ρ) vanishes on S which contradicts the fact that ϕ j (S, ρ) = 1 along Σ j .

Then, the second inequality of (36) amounts to

j Y ij (ρ, S) > 0. We have j Y ij (ρ, S) = S ρ ∇ϕ i (S, ρ)
• ∇ϕ e (S, ρ) where ϕ e (S, ρ) := j ϕ j (S, ρ)

satisfies the boundary value problem

         div ρ ∇ϕ e (ρ, S) = 0, in S, ϕ e (ρ, S) = 0, on ∂S e , ϕ e (ρ, S) = 1, on Σ k , k ≤ N. (37)
Again, by the maximum principle, ϕ e (ρ, S) is maximal on each Σ k , thus ∂ n ϕ e (S, ρ) ≥ 0 on Σ i , while Green's formula leads to the identity

j Y ij (ρ, S) = Σi ρ ∂ n ϕ e (S, ρ) ≥ 0.
The proof that this inequality is strict is similar to the proof that inequality (36) is strict.

3. A mathematical model for multi-coaxial cables. A cable, with axis x 3 , will be defined as the union of its cross sections (in the (x 1 , x 2 ) plane)

Ω = ∪ x3∈R S x3 × {x 3 }, S x3 ∈ C N cf. (3) , for each x 3 . ( 38
)
We assume that, in addition, ∂Ω is a Lipschitz continuous manifold such that

∂Ω = ∪ x3∈R ∂S x3 × {x 3 }.
Along ∂Ω, we define the 2D vector field of unit normal vectors n : ∂Ω → R 2 such that, along ∂S x3 , n is the unit normal vector to ∂S x3 , outgoing with respect to S x3 (this is coherent with the notation of section 2). On ∂Ω, a field of (non unitary) normal vectors to ∂Ω, outgoing with respect to Ω, will be defined by:

n := (n, g) : ∂Ω -→ R 3 , n : ∂Ω -→ R 2 , g : ∂Ω -→ R. (39)
In other words g denotes the longitudinal component of the normal vector to ∂Ω whose projection on the (x 1 , x 2 ) plane is n. According to the notation of section 2,

S x3 = O x3 \ N j=1 O x3 j ,
where the holes O x3 j in S x3 , with boundaries Σ x3 j , define the metallic wires (filled by perfectly conducting metal)

W j = ∪ x3∈R O x3 j × {x 3 }. ( 40
)
All this notation is summarized in figure 2 in the case N = 3.

x 1

x 2

x 3 Example : A fundamental example is the case where the cross sections are obtained by smooth deformations of a reference section, the domain S of section 2. More precisely, let {T x3 , x 3 ∈ R} be a family of diffeomorphisms of R 2 from which we define T : R 3 -→ R 3 such that

T (x, x 3 ) := T x3 (x), (41)
we define the cross sections as

∀ x 3 ∈ R, S x3 = T x3 S, (42)
so that S x3 has the same topological properties and regularity as S (see [START_REF] Delfour | Zolésio Shapes and geometries. Analysis, differential calculus, and optimization, Advances in Design and Control SIAM[END_REF]). If we assume that T ∈ C 1 (R 3 , R 3 ), we can define the vector field V ∈ C 0 (R 3 , R 3 ) of Eulerian velocities (we use this vocabulary by analogy with fluid mechanics, assuming that the variable x 3 plays the role of the time), such that

y = T (x, x 3 ) ≡ T x3 (x) =⇒ V(y, x 3 ) = ∂ 3 T (x, x 3 ). ( 43
)
It is then an exercise in differential geometry to check that ∂Ω is a Lipschitz manifold and that (39) is satisfied with

g = -V • n. (44)
We shall consider a lossy dielectric material whose characteristics coefficients, electric permittivity, magnetic permeability, electric and magnetic conductivities (responsible for losses), satisfy the usual assumptions:

0 < ε -≤ ε(x) ≤ ε + , 0 < µ -≤ µ(x) ≤ µ + , x ∈ Ω, 0 ≤ σ e,-≤ σ e (x) ≤ σ e,+ , 0 ≤ σ m,-≤ σ m (x) ≤ σ m,+ , x ∈ Ω. ( 45
)
We consider a family of (thin) domains parametrized by a small but strictly positive scalar δ (the diameter of the cross section in the plane

x 3 = 0) Ω δ = G δ Ω ,
where G δ is the scaling transformation G δ : (x 1 , x 2 , x 3 ) -→ (δx 1 , δx 2 , x 3 ). We assume that the material properties of the thin cables are defined from the ones of the reference cable Ω according to this scaling:

ε δ = ε • G -1 δ , µ δ = µ • G -1 δ , σ δ e = σ e • G -1 δ , σ δ m = σ m • G -1 δ . Remark 4. Abuse of notation.
In what follows, we use the passage from Ω to Ω δ via the change of variable G δ . To make thing clearer, we use the same letters (x, ...) for the coordinates in Ω or Ω δ .

The equations that govern the propagation of electromagnetic waves in the cable Ω δ are Maxwell's equations, whose unknowns are the two time dependent vector fields E δ (the electric field) and H δ (the magnetic field) ( ∇× is the 3D curl operator)

   ε δ ∂ t E δ + σ δ e E δ -∇ × H δ = j δ , x ∈ Ω δ , t > 0, µ δ ∂ t H δ + σ δ m H δ + ∇ × E δ = 0, x ∈ Ω δ , t > 0, (46)
completed with perfectly conducting boundary conditions

E δ × n δ = 0 x ∈ ∂Ω δ , t > 0, (47)
and zero initial conditions:

E δ (x, 0) = 0, H δ (x, 0) = 0, x ∈ Ω δ . ( 48
)
The source current j δ is also defined by rescaling,

j δ = j • G -1 δ
where, for simplicity:

j = ( j T , 0 ) t , j T = ( j 1 , j 2 ) t , div j T = ∂ 1 j 1 + ∂ 2 j 2 = 0, x ∈ Ω. ( 49
)
For the forthcoming asymptotic analysis, it is useful to work in the reference geometry Ω with the rescaled electromagnetic field:

E δ = E δ • G δ , E δ = E δ • G δ (50)
The longitudinal variable x 3 , being invariant by G δ , plays a different role compared to the transverse variables x := (x 1 , x 2 ), it is therefore natural to distinguish the transverse components of the electromagnetic fields from their longitudinal ones:

E δ T = ( E δ 1 , E δ 2 ) t , H δ T = ( H δ 1 , H δ 2 ) t from E δ 3 , H δ 3 .
We can rewrite (46) with these new unknowns and the notation of section 2 (note that a δ -1 factor appears where a transverse derivative, in x 1 or x 2 , is involved):

                   ε ∂ t + σ e E δ T -∂ 3 e 3 × H δ T -δ -1 rot H δ 3 = j T , x ∈ Ω, t > 0, ε ∂ t + σ e E δ 3 -δ -1 rot H δ T = 0, x ∈ Ω, t > 0, µ ∂ t + σ m H δ T + ∂ 3 e 3 × E δ T + δ -1 rot E 3 = 0, x ∈ Ω, t > 0, µ ∂ t + σ m H δ 3 + δ -1 rot E δ T = 0, x ∈ Ω, t > 0. ( 51
)
as well as the boundary conditions (47) (multiplied by δ -1 for convenience)

δ -1 E δ T × n = 0, δ -1 E δ 3 n -g E δ T = 0. ( 52
)
In the next section, our goal is to analyze the behavior of the solution of (51, 52) when δ tends to 0.

4. Asymptotic analysis.

As in section 2, our (formal) developments in this section are fully justified provided that the functions we manipulate are smooth enough . This requires implicitly smoothness assumptions of the geometry of the cable and of the coefficients of the problem. However, we are convinced that these assumptions could be removed.

4.1.

A formal asymptotic expansion. The idea is to conjecture a formal power series expansion in δ of the solution of (51, 52). It is useful to pass to the frequency domain via a time Fourier transform:

E δ (x, t) -→ E δ (x, ω), H δ (x, t) -→ H δ (x, ω). ( 53
)
Taking into account zero initial data, we obtain for each value of the frequency ω, the following problem for the Fourier transform of the electromagnetic field (we use the notation of section 2.1, indicate by the index T (resp. 3) the transverse (and resp. longitudinal) of a 3D vector field):

                   iω ε + σ e E δ T -∂ 3 e 3 × H δ T -δ -1 rot H δ 3 = j T , x ∈ Ω, (a) iω ε + σ e E δ 3 -δ -1 rot H δ T = 0, x ∈ Ω, (b) iω µ + σ m H δ T + ∂ 3 e 3 × E δ T + δ -1 rot E δ 3 = 0, x ∈ Ω, (c) iω µ + σ m H δ 3 + δ -1 rot E δ T = 0, x ∈ Ω, (d) (54)
with the boundary conditions ( deduced from (52) )

δ -1 E δ T × n = 0, (a) δ -1 E δ 3 n -g E δ T = 0. (b) (55)
We shall also use the "hidden" divergence free equations, that are easily deduced from (54) (see for instance [START_REF] Monk | Finite element methods for Maxwell's equations[END_REF])

   δ -1 div (iω ε + σ e ) E δ T + ∂ x3 (iω ε + σ e ) E δ 3 = 0, x ∈ Ω, (a) δ -1 div (iω µ + σ m ) H δ T + ∂ x3 (iω µ + σ m ) H δ 3 = 0, x ∈ Ω, (b) (56)
and the "hidden" boundary condition for the magnetic field (see again [START_REF] Monk | Finite element methods for Maxwell's equations[END_REF])

δ -1 H δ T • n -(V • n) H δ 3 = 0, x ∈ ∂Ω. (57)
Our approach consists in postulating a priori an asymptotic expansion in power series expansion of the form

E δ = E 0 + δ E 1 + δ 2 E 2 + • • • , H δ = H 0 + δ H 1 + δ 2 H 2 + • • • , (58)
and to identify and characterize the formal limit field ( E 0 , H 0 ). For this, we substitute the expansions (58) into equations ( 54) to (57) and identify each power of δ in the resulting series, which begins with the δ -1 terms ( because of the δ -1 factors in (54) to (57) ).

Limit longitudinal fields

Concerning, the fields E 0 3 and H 0 3 , it suffices to proceed as in [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF] (we omit here the details) to show that E 0 3 = H 0 3 = 0 (59) which means that the limit field is transversely polarized. This is considered as an assumption in many textbooks on transmission lines [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF].

Identification of the δ -1 terms.

Equations (54-(d)), (56-(a)) and (55-(b)) lead to rot E 0 T = 0, div (ε + σ e /iω) E 0 T = 0 in Ω, E 0 T × n = 0, on ∂Ω. (60)
Writing (60) for each x 3 , considered as a parameter, one realizes that, according to section 2.2, we have

∀ x 3 ∈ R, E 0 T (., x 3 ) ∈ E ε(., x 3 , ω), S x3 , ε(., x 3 , ω) := ε(., x 3 ) + σ e (., x 3 ) iω . (61)
Therefore, setting

ϕ e,j (•, x 3 , ω) := ϕ j ε(., x 3 , ω), S x3 ∈ H 1 (S x3 ), (62)
which defines N scalar fields ϕ e,j (x 3 , ω) : Ω → C, we know that there exist N 1D (in space) complex valued functions V j (x 3 , ω) (electric potentials), such that:

∀ x 3 ∈ R, E 0 T (., x 3 , ω) = N j=1 V j (x 3 , ω) ∇ ϕ e,j (•, x 3 , ω). (63)
In the same way, equations (54-(b)), (56-(b)) and (57) lead to

rot H 0 T = 0, div (µ + σ m /iω) H 0 T = 0 in Ω, H 0 T • n = 0 on ∂Ω. (64)
Writing (60) for each x 3 , one sees that, according to section ,

∀ x 3 ∈ R, H 0 T (., x 3 ) ∈ H µ(., x 3 , ω), S x3 , µ(., x 3 , ω) := µ(., x 3 ) + σ m (., x 3 ) iω . (65)
Therefore, setting

ψ m,j (•, x 3 , ω) := ψ j µ(., x 3 , ω), S x3 ∈ H 1 (S Γ x3 ), (66)
which defines N scalar fields ψ m,j (x 3 , ω) : Ω → C, we know that there exist N 1D (in space) complex valued functions I j (x 3 , ω) (electric currents), such that:

∀ x 3 ∈ R, H 0 T (., x 3 , ω) = N j=1 I j (x 3 , ω) ∇ ψ m,j (•, x 3 , ω). (67)
Formulas ( 63) and (67) express a generalized separation of variables, valid at the limit δ → 0. The fields ∇ϕ e,j (•, x 3 , ω) and ∇ψ m,j (•, x 3 , ω) can be precomputed by solving 2D boundary value problems in the cross sections S x3 . It remains to determine the equations that determine the functions V j (x 3 , ω) and I j (x 3 , ω), which will provide us our effective model, the generalized telegrapher's model. For this we need to proceed to the next step of the identification process.

Identification of the δ 0 terms.

Using ( 63) and (67), the δ 0 terms in equations (54-(a)) and (54-(c)) are

           iω ε + σ e N j=1 V j ∇ ϕ e,j -∂ 3 e 3 × N j=1 I j ∇ ψ m,j -rot H 1 3 = j T , (a) iω µ + σ m N j=1 I j ∇ ψ m,j + ∂ 3 e 3 × N j=1 V j ∇ ϕ e,j + rot E 1 3 = 0, (b) (68)
while ( 55) gives

E 1 3 n -g E 0 T = 0 (69)
To get rid of the (unknown) terms H 1 3 and E 1 3 , we take the L 2 scalar product, in each cross section S x3 , of (68-(a)) with each ∇ ϕ e,i and of (68-(b)) with each ∇ ψ m,i . This is motivated by the following identities (70) and (71).

First, using Green's formula, we have, in each section

S x3 rot H 1 3 , ∇ ϕ e,i L 2 (Sx 3 ) = ∂Sx 3 H 1 3 ∇ ϕ e,i × n = 0, (70)
since ∇ ϕ e,i × n = 0 along ∂S x3 and rot ∇ ϕ e,i = 0.

On the other hand, thanks to the jump relations satisfies by ψ m,i , the reader will easily verify that rot ∇ ψ m,i = 0 in S x3 . Therefore, by Green's formula

rot E 1 3 , ∇ ψ m,i L 2 (Sx 3 ) = ∂Sx 3 E 1 3 ∇ ψ m,i × n.
Then, using (69), we have, along ∂S x3 :

E 1 3 ∇ ψ m,i × n ≡ -E 1 3 rot ψ m,i • n = -g E 0 T • n rot ψ m,i • n = -g E 0 T • rot ψ m,i
, where the last equality results from the fact that ∂ n ψ m,i = 0 means that rot ψ m,i is normal to ∂S x3 . Consequently, using Green's formula and ( 63)

rot E 1 3 , ∇ ψ m,i L 2 (Sx 3 ) = - Sx 3 g E 0 T • rot ψ m,i = - N j=1 ∂Sx 3 g ∇ ϕ e,j • rot ψ m,i V j . (71)
For the second terms in (68-(a) and (b)), we shall use the following lemma, whose proof will be omitted since it is a straightforward adaptation of lemma 3.4 in [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF] Lemma 4.1.

(i) ∂ 3 rot ϕ e,j , ∇ ψ m,i L 2 (Sx 3 ) + ∂Sx 3 g ∇ ϕ e,j • rot ψ m,i = 0. (ii) ∂ 3 rot ψ m,j , ∇ϕ e,i L 2 (S Γ x 3 ) = 0.
Using e 3 × ∇ = -rot, the bi-orthogonality relations (21) and lemma 4.1, we have

-∂ 3 e 3 × N j=1 I j ∇ ψ m,j , ∇ ϕ e,i L 2 (Sx 3 ) = ∂ 3 I i (72) ∂ 3 e 3 × N j=1 V j ∇ ϕ e,j , ∇ ψ m,i L 2 (Sx 3 ) = ∂ 3 V i + N j=1 ∂Sx 3 g ∇ ϕ e,j • rot ψ m,i V j . (73)
Finally, taking the L 2 scalar product in S x3 of (68-(a)) with ∇ ϕ e,i and (68-(b)) with ∇ ψ m,i , we obtain, using (70, 71, 72, 73), the equations for the potentials V j and the currents I j .

To write this model in a compact way, we introduce the vectors

V(x 3 , ω) := V j (x 3 , ω) 1≤j≤N , I(x 3 , ω) := I j (x 3 , ω) 1≤j≤N (74)
and define the matrices (with the notation ( 22) and (23) of section 2)

Y e (x 3 , ω) := Y( ε(•, x 3 , ω), S x3 ), Z m (x 3 , ω) := Z( µ(•, x 3 , ω), S x3 ), (75)
where ε and µ have been defined in (61) and (65). One then gets

   iω Y e (x 3 , ω) V(x 3 , ω) + ∂ 3 I(x 3 , ω) = J T (x 3 , ω), (a) iω Z m (x 3 , ω) I(x 3 , ω) + ∂ 3 V(x 3 , ω) = 0, (b) (76)
where the vector of sources is defined by

J T (x 3 , ω) = J T,i (x 3 , ω) 1≤i≤N , J T,i (x 3 , ω) = Sx 3 j T • ∇ ϕ e,i (ω). (77)
4.2. Derivation of the generalized telegrapher's model in time domain.

Our objective is to write the time domain formulation of equation (76). This is done formally by an inverse Fourier transform that can be applied on (76). We shall do it in a way that emphasizes energy preservation (or dissipation) relations and related stability properties of the resulting evolution problem. Towards this goal, we introduce the high frequency behavior of the functions ϕ e,j and ψ m,j .

These high frequency limits are naturally defined after noticing that ε → ε and µ → µ when ω → ±∞.

This justifies to introduce

ϕ ∞ e,j (•, x 3 ) := ϕ j ε(., x 3 ), S x3 , ψ ∞ e,j (•, x 3 ) := ψ j µ(., x 3 ), S x3 , (78)
as well as the (frequency dependent) residual terms ϕ r e,j (•, x 3 , ω) := ϕ e,j (•,

x 3 , ω) -ϕ ∞ e,j (•, x 3 ), ψ r m,j (•, x 3 , ω) := ψ m,j (•, x 3 , ω) -ψ ∞ m,j (•, x 3 ). ( 79
)
By definition of ϕ ∞ e,j , we remark that ϕ r e,j ∈ H 1 (S x3 ) is the unique solution of

   div ε ∇ ϕ r e,j = - 1 iω div σ e ∇ϕ ∞ e,j , in S,
ϕ r e,j = 0, on ∂S.

(

) 80
Multiplying this equation by ϕ r e,j , we get the relation

Sx 3 (iωε + σ e ) |∇ ϕ r e,j | 2 = - Sx 3 σ e ∇ ϕ ∞ e,j • ∇ ϕ r e,j ,
which enables us to deduce the following estimates (by simply taking the imaginary and the real part of the previous equality)

∇ ϕ r e,j L 2 (Sx 3 ) ≤ 1 ε -|ω| σ e ∇ ϕ ∞ e,j L 2 (Sx 3 ) , ( 81
)
σ 1/2 e ∇ ϕ r e,j L 2 (Sx 3 ) ≤ σ 1/2 e ∇ ϕ ∞ e,j L 2 (Sx 3 ) . ( 82
)
Similarly, the residual function ψ r m,j ∈ H 1 (S x3 ) is the unique solution of

     div µ ∇ ψ r m,j = - 1 iω div σ m ∇ψ ∞ m,j , in S, ∇ ψ r m,j • n = 0, on ∂S, (83)
and satisfies

∇ ψ r m,j L 2 (Sx 3 ) ≤ 1 µ -|ω| σ m ∇ ψ ∞ m,j L 2 (Sx 3 ) , σ 1/2 m ∇ ψ r m,j L 2 (Sx 3 ) ≤ σ 1/2 m ∇ ψ ∞ m,j L 2 (Sx 3 )
.

The high frequency decomposition (79) yields the following proposition Proposition 1. The matrices Y e (x 3 , ω) and Z m (x 3 , ω) satisfy

iω Y e (x 3 , ω) = iω C ∞ (x 3 ) + G ∞ (x 3 ) + K e (x 3 , ω), ( 84
) iω Z m (x 3 , ω) = iω L ∞ (x 3 ) + R ∞ (x 3 ) + K m (x 3 , ω), (85)
where, the N × N real symmetric non negative matrices

(C ∞ , G ∞ , L ∞ , R ∞ ) are defined by (C ∞ ) ij := Sx 3 ε ∇ϕ ∞ e,i • ∇ϕ ∞ e,j dx, (G ∞ ) ij := Sx 3 σ e ∇ϕ ∞ e,i • ∇ϕ ∞ e,j dx, (L ∞ ) ij := Sx 3 µ ∇ψ ∞ m,i • ∇ψ ∞ m,j dx, (R ∞ ) ij := Sx 3 σ m ∇ψ ∞ m,i • ∇ψ ∞ m,j dx.
(86) The matrices C ∞ and L ∞ are positive definite, C ∞ is a strict M-matrix and L ∞ has strictly positive entries. Finally, the frequency dependent symmetric matrices ( K e , K m ) are defined by

( K e ) ij = Sx 3 σ e ∇ ϕ r e,i • ∇ϕ ∞ e,j dx, ( K m ) ij = Sx 3 σ m ∇ ψ r m,i • ∇ψ ∞ m,j dx. ( 87
)
Proof. We prove only (84), the argument being similar to prove (85). From the definition (22, 75) of Y e (x 3 , ω) and the intermediate result ( 24) we have

Y e = Sx 3 ε ∇ ϕ e,i • ∇ ϕ e,j dx = Sx 3 ε ∇ ϕ e,i • ∇ϕ ∞ e,j dx,
where the last equality is obtained using Remark 3 with ρ = ε and ϕ * = ϕ r e,j . Using again Remark 3 with ρ = ε and ϕ * = ϕ r e,i , we can simply further the previous expression:

Y e = Sx 3 ε ∇ϕ ∞ e,i • ∇ϕ ∞ e,j dx + Sx 3 σ e iω ∇ ϕ r e,i • ∇ϕ ∞ e,j dx,
which conclude the proof of the decomposition. The symmetry property of ( K e , K m ) is a direct consequence of the obvious symmetry of (C ∞ , G ∞ , L ∞ , R ∞ ) and of ( Y e , Z m ) (as proven in lemma 2.2 ). The properties of the matrices C ∞ and L ∞ are straightforward application of lemma 2.6.

Lemma 4.2. The function of the frequency

( K e ) ij (x 3 , •), ( K m ) ij (x 3 , •) satisfy K e (x 3 , •) ∈ L 2 R, M N (C) , K m (x 3 , •) ∈ L 2 R, M N (C) ,
and, for all

u ∈ C N , inf ω∈R Re G ∞ (x 3 ) + K e (x 3 , •) u, u C N ≥ 0, inf ω∈R Re R ∞ (x 3 ) + K m (x 3 , •) u, u C N ≥ 0.
Proof. Using (87), Cauchy-Schwartz inequality and estimate (81), we get

|( K e ) ij | ≤ ∇ ϕ r e,j L 2 (Sx 3 ) σ e ∇ ϕ ∞ e,j L 2 (Sx 3 ) ≤ 1 µ -|ω| σ e ∇ ϕ ∞ e,j 2
L 2 (Sx 3 ) .
In a similar way, using (82) we can deduce an uniform, frequency independent, bound on ( K e ) ij (x 3 , •):

|( K e ) ij | ≤ σ 1/2 e ∇ ϕ r e,j L 2 (Sx 3 ) σ 1/2 e ∇ ϕ ∞ e,j L 2 (Sx 3 ) ≤ σ 1/2 e ∇ ϕ ∞ e,j 2
L 2 (Sx 3 )
.

Combining the two previous estimates we deduce the first result of the lemma which can be proven for K m in a very similar way. Next the positivity properties follow from the simple observation that

Re[G ∞ (x 3 ) + K e (x 3 , •)] ij = Re[iω Y e ] ij = Sx 3 σ e ∇ ϕ e,i • ∇ ϕ e,j dx,
which define a positive matrix since σ e ≥ 0 (the argument uses the expression (25) of the scalar product ).

Looking at (80) suggests to introduce the causal function ϕ r e,j (x, x 3 , t) defined by ϕ r e,j (x 3 , x, t) = 0 for t < 0 and is solution, for t ≥ 0, of the Cauchy problem

           div ε ∂ ∂t ∇ϕ r e,j + σ e ∇ϕ r e,j = 0, in S x3 , ϕ r e,j = 0, on ∂S x3 , div ε ∇ϕ r e,j = -div σ e ∇ϕ ∞ e,j , in S x3 at t = 0. (88)
The well-posedness of ( 88) is proved as in [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF], Section 5, with

ϕ r e,j (•, x 3 , •) ∈ C ∞ (R + ; H 1 0 (S x3 )
). These results are due to the fact that (88) can be rewritten as an abstract evolution problem in H 1 0 (S x3 ) of the form (see [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF])

d dt ϕ r e,j (•, x 3 , •) + A e (x 3 ) ϕ r e,j (•, x 3 , •) = 0, A e (x 3 ) ∈ L H 1 0 (S x3 ) .
With the same arguments the function ψ r m,j (x 3 , x, t) is the causal function uniquely defined for t > 0 as the solution of the evolution problem:

           div µ ∂ ∂t ∇ψ r m,j + σ m ∇ψ r m,j = 0, in S x3 , ∇ψ r m,j • n = 0, on ∂S x3 , div µ ∇ψ r m,j = -div σ m ∇ψ ∞ m,j , in S x3 at t = 0. ( 89
) which satisfies ψ r m,j (•, x 3 , •) ∈ C ∞ (R + ; H 1 (S x3
)). We can now deduce a convenient time domain formulation of the generalized telegrapher's model from the equations (76) with the coefficients decomposed following proposition 1. Using the standard properties of the Fourier transform with respect to convolution, we obtain

         C ∞ (x 3 ) ∂ t V + G ∞ (x 3 ) V + t 0 K e (x 3 , t -s) V(x 3 , s) ds + ∂ 3 I = J T , L ∞ (x 3 ) ∂ t I + R ∞ (x 3 ) I + t 0 K m (x 3 , t -s) I(x 3 , s) ds + ∂ 3 V = 0, (90)
this set of equations being completed with zero initial conditions

V(x 3 , 0) = I(x 3 , 0) = 0. (91)
In (90) the convolution kernels

(K e (x 3 , •), K m (x 3 , •)) ∈ L 2 R + ; M N (R)
are the inverse Fourier transform in time of ( K e (x 3 , •), K m (x 3 , •)). According to (87), it is easy to see that they are given by the following explicit expressions

(K e ) ij (x 3 , t) = Sx 3 σ e (x, x 3 ) ∇ϕ r e,i (x, x 3 , t) • ∇ϕ ∞ e,j (x, x 3 ) dx, (K m ) ij (x 3 , t) = Sx 3 σ m (x, x 3 ) ∇ψ r m,i (x, x 3 , t) • ∇ψ ∞ m,j (x, x 3 ) dx. (92)
The model ( 90) is our effective 1D model in which:

• the effective coefficients C ∞ (x 3 ), G ∞ (x 3 ), L ∞ (x 3 ), R ∞ (x 3 ) are deduced, through (86)
, from the solution of 2D elliptic problems in S x3 , see (18, 19, 78), • the kernels (K e (x 3 , •), K m (x 3 , •)) are deduced, through formulas (92), from the solution of 2D evolution problems in S x3 , namely (88) and (89). Finally it is worthwhile mentioning that from the solution (V, I), one can reconstruct (approximately) the 3D electromagnetic field through the following formulas, deduced from formulas (63, 67) and the decompositions (79):

E δ (x, x 3 , t) ∼ N j=1 ∇ϕ ∞ e,j x δ , x 3 , t V j (x 3 , t) + N j=1 t 0 ∇ϕ r e,j x δ , x 3 , t -s V j (x 3 , s) ds, H δ (x, x 3 , t) ∼ N j=1 ∇ψ ∞ m,j x δ , x 3 , t I j (x 3 , t) + N j=1 t 0 ∇ψ r m,j x δ , x 3 , t -s I j (x 3 , s) ds. (93)
Remark 5. We call a generalized telegrapher's model because the unknowns I and V are vector valued, and not scalar, and because of the time convolution terms that represent memory effects. Let us emphasize that the presence of these memory terms is due to the conjugated presence of losses in the medium and heterogeneity of the cross section. More precisely, as in [START_REF] Imperiale | Mathematical modeling of electromagnetic wave propagation in heterogeneous lossy coaxial cables with variable cross section[END_REF], it is easy to show that:

K e (x 3 , •) = 0 ⇐⇒ (σ e /ε (x 3 , •) is constant in S x3 , K m (x 3 , •) = 0 ⇐⇒ (σ m /µ (x 3 , •) is constant in S x3 . (94)
5. Properties of the generalized telegrapher's model.

5.1.

Well-posedness of the limit model.

Theorem 5.1. If the source term as the regularity

J T ∈ L 1 (R + , L 2 (R 3 ))
and there exist strictly positive constants C -and L -such that, for all

u ∈ R inf x3∈R (C ∞ (x 3 )u, u) R N ≥ C -|u| 2 R N , inf x3∈R (L ∞ (x 3 )u, u) R N ≥ L -|u| 2 R N , (95)
then, there exists a unique weak solution (V, I) ∈ C 0 R + , L 2 (R) 2 of (90, 91) and

V(•, t) 2 L 2 (R) + I(•, t) 2 L 2 (R) ≤ 1 min(C -, L -) t 0 J T (•, s) L 2 (R) ds. ( 96
)
The proof of the above result relies on the introduction of the following two bilinear forms in L 2 0, t; L 2 (R 3 ) :

a e t; U, U = t 0 R G ∞ (x 3 ) U(x 3 , s), U(x 3 , s) R N dx 3 ds + t 0 s 0 R K e (x 3 , s -r) U(x 3 , r), U(x 3 , s) R N dx 3 dr ds, a m t; U, U = t 0 R R ∞ (x 3 ) U(x 3 , s), U(x 3 , s) R N dx 3 ds + t 0 s 0 R K m (x 3 , s -r) U(x 3 ,
∀ U ∈ L 2 0, t; L 2 (R 3 ), a t; U, U ≥ 0, = e, m. (97)
Proof. We give the proof for a e . We extend U outside [0, t] into a function U such that U (•, s) = U(•, s) for s ∈ [0, t] and U (•, s) = 0 otherwise, then we have

a e t; U, U = R R (G ∞ (x 3 ) U (x 3 , s), U (x 3 , s)) R N dx 3 ds + R R R
(K e (x 3 , sr) U (x 3 , r), U (x 3 , s)) R N dx 3 dr ds.

Denoting U the time Fourier transform of U , we get, by Parseval identity

a e t; U, U = R R Re [G ∞ (x 3 ) + K e (x 3 , ω)] U (x 3 , ω), U (x 3 , ω)
C N dω dx 3 . We then conclude thanks to lemma 4.2.

Proof of theorem 5.1. We first assume that j T is regular enough and that a corresponding regular and integrable solution of (90, 91) exists. We multiply the two equations of (90) by V and I respectively, integrate in space over R and in time between 0 and t and add the two resulting equalities. We use (95) then to obtain the following inequality:

C - 2 V(•, t) 2 L 2 (R) + L - 2 I(•, t) 2 L 2 (R) + a e t; V, V + a m t; I, I ≤ t 0 R (J T , V) R N dx 3 ds, (98)
Lemma 5.2 then gives

C - 2 V 2 L 2 (R) + L - 2 I 2 L 2 (R) ≤ t 0 J T (•, t) L 2 (R) V L 2 (R) ds,
and Gronwall's lemma enables us to derive the energy estimate (96). This estimate is the key tool to deduce existence and uniqueness of the solution in the right functional space. We shall not reproduce the (rather standard) proof but refer the reader to [START_REF] Cockburn | Maxwell equations in polarizable media[END_REF].

Theorem 5.1 guarantees the well-posedness of the generalized telegrapher's equation under the assumptions (95). In the next section, we investigate the assumptions to be done on the original 3D problem (these are not straightforward) in order that the effective coefficients given by (86) do satisfy this assumption.

then (C ∞ (x 3 )u, u) C N ≥ ε -λ D (B * ) m * |u| 2 C N , (L ∞ (x 3 )u, u) C N ≤ µ + λ D (B * ) -1 m -1 * |u| 2 C N , (100)
where λ D (B * ) denotes the lowest eigenvalue for the Dirichlet problem associated to the Laplace operator in the domain D.

(ii) If one makes the stronger assumption (the reader will easily check that (101) implies (99) ) that there exists

S * = O * \ N j=1 O * j ∈ C N
and for each x 3 ∈ R, a conformal mapping T c,x3 such that (see figure 3)

T c,x3 (S x3 ) ⊂ S * , (101)
then there exists a positive constant A(S * ), depending only on S * , such that (a characterization of A(S * ) is given in the proof below) Proof. We simply need to show (100, 102) when T c,x3 = Id. The proof can then be naturally extended to any conformal mapping using lemma 2.5.

(C ∞ (x 3 )u, u) C N ≤ ε + A(S * ) |u| 2 C N , (L ∞ (x 3 )u, u) C N ≥ µ -A(S * ) -1 |u| 2 C N . (102
(i) Setting ϕ(u) := N j=1 u j ϕ ∞ e,j (•, x 3 ), we have

(C ∞ (x 3 )u, u) C N = Sx 3 ε(•, x 3 ) ∇ϕ(u) 2 ≥ ε -(x 3 ) Sx 3 ∇ϕ(u) 2 .
We extend ϕ(u) as a function ϕ(u) ∈ H 1 0 (B * ), piecewise constant outside S x3 :

ϕ(u) = ϕ(u) in S x3 , ϕ(u) = 0 in B * \ O x3 , ϕ(u) = u i in O x3 j , j = 1, ..., N,
in such a way that

Sx 3 ∇ϕ(u) 2 = B * ∇ ϕ(u) 2 . Therefore, since O x3 j ⊂ B * , (C ∞ (x 3 )u, u) C N ≥ ε -(x 3 ) λ D (B * ) N j=1 O x 3 j ϕ(u) 2 = ε -(x 3 ) λ D (B * ) N j=1 |u j | 2 meas(O x3 j ).
It is then easy to obtain the first inequality of (100). To obtain the second inequality, in the above reasoning, it suffices to replace ε by µ

-1 , ϕ ∞ e,j (•, x 3 ) by ϕ ∞ m,j (•, x 3 ) := ϕ(µ -1 (•, x 3 ), S x3 ) (103)
and to use lemma 2.4 that says that L

∞ (x 3 ) -1 = Y(µ -1 (•, x 3 ), S x3
). The details are left to the reader. Note that one obtains a lower bound for L ∞ (x 3 ) -1 , which results into an upper bond for L ∞ (x 3 ).

(ii) By the Dirichlet principle, we have

(C ∞ (x 3 )u, u) C N = inf ϕ∈V (x3,u) Sx 3 ε(•, x 3 ) ∇ϕ 2 , where V (x 3 , u) = ϕ ∈ H 1 (S x3 ) / ϕ| Σ x 3 j = u j , ϕ| ∂O x 3 = 0 . Let us set ϕ * j = ϕ j (S * , 1) ∈ H 1 (S * ) and ϕ * (u) = N j=1 u j ϕ * j ∈ H 1 (S * ).
We can extend ϕ * (u) in S x3 as a function ϕ * (u) in H 1 (S x3 ) by

ϕ * (u) = ϕ * (u) in S * , ϕ * (u) = 0 in O x3 \ O * , ϕ * (u) = u i in O x3 j \ O * j , 1 ≤ j ≤ N. By construction ϕ * (u) ∈ V (x 3 , u). Thus (C ∞ (x 3 )u, u) C N ≤ Sx 3 ε(•, x 3 ) ∇ ϕ * (u) 2 = S * ε(•, x 3 ) ∇ϕ * (u) 2 .
It is then easy to obtain the first inequality of (102) with

A(S * ) = sup u∈C N , |u|=1 S * ∇ϕ * (u) 2 .
To prove the second inequality, we use again that L ∞ (x 3 ) -1 = Y(µ -1 (•, x 3 ), S x3 ) and apply the above reasoning to obtain un upper bound for L ∞ (x 3 ) -1 . This leads to the announced result.

Propagation velocities.

In the absence of source, the generalized telegrapher's model (90) is a zero order perturbation of the symmetric hyperbolic system (as it is the case for instance in [START_REF] Canadas | Speed of Propagation of Solutions of a Linear Integro-differential Equation with Nonconstant Coefficients[END_REF])

   C ∞ (x 3 ) ∂ t V + ∂ 3 I = 0, L ∞ (x 3 ) ∂ t I + ∂ 3 V = 0. (104)
In the homogeneous case (C ∞ (x 3 ) = C ∞ , L ∞ (x 3 ) = L ∞ ), looking for traveling wave solutions of the form

V(x 3 , t) = V a f (x 3 -λt), I(x 3 , t) = I a f (x 3 -λt), λ ∈ R
leads to characterize the propagation velocities λ via the eigenvalue problem

C -1 ∞ L -1 ∞ I a = λ 2 I a ⇐⇒ C -1 ∞ L -1 ∞ V a = λ 2 V a , (V a = L ∞ I a )
. In other words, ± λ is an admissible velocity associated to (104) if and only if

λ 2 ∈ σ C -1 ∞ L -1 ∞ , where σ C -1 ∞ L -1 ∞ ⊂ R + * since C ∞ and L ∞ are real symmetric positive definite.
By extension, the (variable in space) propagation velocities associated to the system (104) at point x 3 are the reals ± λ(x 3 ) such that

λ(x 3 ) 2 ∈ σ C ∞ (x 3 ) L ∞ (x 3 ) . ( 105
)
Our goal is to provide lower and upper bounds for these velocities. Not surprisingly, these bounds will be related to the function

c(x, x 3 ) = ε(x, x 3 ) µ(x, x 3 ) -1 2 , (106)
which represents the velocity of electromagnetic waves at point (x, x 3 ). The first result in this direction is:

Lemma 5.4. Assume that the velocity of electromagnetic waves is constant in the cross section

S x3 : c(x, x 3 ) = c(x 3 ) (107) then C -1 ∞ (x 3 ) L -1 ∞ (x 3 ) = c(x 3 ) 2 Id. Proof. Let us remind that C ∞ (x 3 ) = Y ε(•, x 3 ), S x3 = c(x 3 ) -2 Y µ(•, x 3 ) -1 , S x3 ,
where we have used (106) and (107). According to lemma 2.4,

C ∞ (x 3 ) = c(x 3 ) -2 Z µ(•, x 3 ), S x3 -1 = c(x 3 ) -2 L ∞ (x 3 ) -1 .
The above result is not true in general as it can be illustrated by the following numerical results. We consider two cross sections with the same geometry and three holes. These two configurations are described in figure 4.

In one case, the function c(x) is constant and on the other case it is not. The results that we obtain via a finite element calculation are the following: for the first situation c(x) = 1, and the computed matrices (C ∞ , L ∞ ) are given by It is easy to see that saying λ(x 3 ) 2 ∈ σ C ∞ (x 3 ) -1 L ∞ (x 3 ) -1 is equivalent to

C ∞ = L ∞ (x 3 ) -1  
|λ(x 3 )| ∈ σ C ∞ (x 3 ) -1 2 L ∞ (x 3 ) -1 2 ,
thus, there exists w ∈ C N , with |w| C N = 1, such that

L ∞ (x 3 ) -1 2 w = |λ(x 3 )| C ∞ (x 3 ) 1 2 w.
Therefore, choosing u = C ∞ (x 3 ) -1 2 w and v = L ∞ (x 3 ) -1 2 w, we have We then get

(u, v) C N = C ∞ (x 3 ) -1 2 w, L ∞ (x 3 ) -
(u, v) C N ≤ c -(x 3 ) -1 Sx 3 ε(•, x 3 ) -1 ∇ϕ(u) 2 1 2 Sx 3 µ(•, x 3 ) -1 ∇ψ(v)

A remark on undetectable defects.

We finish this section by a remark linked to the use of our models in the context of non destructive testing. Then as a consequence of lemma 2.5, the two effective 1D telegrapher's models associated to the two reference cables are the same.

In the context of non destructive testing, let us first consider an ideal perfectly cylindrical cable corresponding to: 31)), then it is not possible to detect the defect with the use of electromagnetic waves by simply using our effective model.

Ω * = S ×
In such a situation, either the original 3D model, or hopefully a higher order asymptotic model, is needed.

An example: the Joukowsky's transformation. To illustrate our discussion we consider the specific case of two homogeneous straight coaxial cables whose sections are obtained one form the other using the Joukoswky's conformal mapping: J (x 1 , x 2 ) = x 1 + x 1 /(x 2 1 + x 2 2 ) , x 2x 2 /(x 2 1 + x 2 2 ) T , which, in particular maps circles into ellipses. In the numerical results represented in Figure 5, the reference sections S considered are annuli of inner radius r -and outer radius r + . In that specific case (one hole, homogenous straight cables) the (scalar) capacitances for the section S and J (S) are the same and are given by C ∞ = ε ∇ϕ 1 (1, S) 2 S = ε ∇ϕ 1 1, J (S) 2 S . Figure 5 we plot the annulus and its transform by J in the case (r -, r + ) = (1.2, 2). For the annulus it is is know that the (scalar) capacitance is given by C ∞ = 2πε ln(r + /r -)

12.30005899 ε.

For the transformed domain J (S), we have used a Q5-finite element method and the computational mesh Figure 5 to compute the capacitance numerically, which gave C ∞ 12.3001 ε.

Figure 1 .

 1 Figure 1. The geometry of the domain S.

Figure 2 .

 2 Figure 2. The geometry of the domain Ω.

5. 2 .

 2 Bounds for the effective coefficients. Lemma 5.3. (i) Assume that there exists a simple connected bounded domain B * such that, for each x 3 ∈ R, there exists a conformal mapping T c,x3 such that T c,x3 (S x3 ) ⊂ B * and m * := inf

Figure 3 .

 3 Figure 3. Illustration of the assumption (101).

1 2 w

 2 C N = |λ(x 3 )|, while (C ∞ (x 3 )u, u) C N = |w| 2 C N = 1 and (L ∞ (x 3 )v, v) C N = |w| 2 C N = 1. Then, (109) leads to |λ(x 3 )| ≤ c + (x 3 ).For the lower bound, we repeat the same reasoning as above with the functions defined by (103), ψ e,j (•,x 3 ) := ψ(ε -1 (•, x 3 ), S x3 ) and ϕ(u) := N j=1 u j ϕ ∞ m,j (•, x 3 ) and ψ(v) := N j=1v j ψ ∞ e,j (•, x 3 ).

2

 2

1 2 that 2 (L ∞ (x 3 ) 1 2 1 2

 122311 is to say, instead of (109),(u, v) C N ≤ c -(x 3 ) -1 (C ∞ (x 3 ) -1 u, u) C N 1 -1 v, v) u = C ∞ (x 3 ) w and v = L ∞ (x 3 ) w, leads to |λ(x 3 )| ≥ c -(x 3 ).

First, let us

  consider two reference cables (Ω, ε, µ, σ e , σ m ) and ( Ω, ε, µ, σ e , σ m ), such that for any x 3 ∈ R, there exists a conformal mappingT c,x3 (cf. (31)) such that ∀ x 3 ∈ R, S x3 = T c,x3 (S x3 )(111)as well asε(•, x 3 ) = ε(•, x 3 ) • T c,x3 , µ(•, x 3 ) = µ(•, x 3 ) • T c,x3 , σ e (•, x 3 ) = σ e (•, x 3 ) • T c,x3 , σ m (•, x 3 ) = σ m (•, x 3 ) • T c,x3 .

R

  and coefficients (ε * , µ * , σ * e , σ * m ) that are only functions of the transverse variable x. Let us consider a cable with a local defect defined, according to our example, byS x3 = T c,x3 (S), where each T c,x3 is a diffeomorphism of R 2 such that T c,x3 ≡ Id for x 3 / ∈ [a, b], T c,x3 = Id for x 3 ∈ [a, b],which corresponds to a geometrical defect localized in the portion [a, b] of the cable. Assume that the coefficients for this cable are given byε(•, x 3 ) = ε * • T c,x3 , µ(•, x 3 ) = µ * • T c,x3 , σ e (•, x 3 ) = σ * e • T c,x3 , σ m (•, x 3 ) = σ * m • T c,x3. Then, if one assumes that, for x 3 ∈ [a, b], T c,x3 is a conformal mapping (cf. (

Figure 5 .

 5 Figure 5. Left: Computational mesh. Right: level lines of ϕ 1 1, S) and ϕ 1 1, J (S) .

  r), U(x 3 , s) R N dx 3 dr ds, that have the important positivity property: Lemma 5.2. The bilinear forms a e t; •, • and a m t; •, • satisfy:

Acknowledgments. The authors would like to thank Michel Sorine (Inria Rocquencourt) for helpful discussion and suggestions. They also acknowledge the support of the ANR project SODDA-2011-BS-013-01.

c(x, x 3 ), we have:

Using the biorthogonality formulas (21) with

which can be rewritten, using the definition of C ∞ (x 3 ) and L ∞ (x 3 ), (u, v)

(109)


Dang The

Huynh

Chung Shue Chen

Siu-Wai Ho
email: siuwai.ho@unisa.edu.au

Exploiting User Movement for Position Detection

Keywords: Positioning system, localization algorithm, user movement, mobile device, smart applications

The major issue of indoor localization system is the trade-off between implementation cost and accuracy. A low-cost system which demands only few hardware devices could save the cost but often it turns out to be less reliable. Aiming at improving classical triangulation method that requires several reference points, this paper proposes a new method, called Two-Step Movement (2SM), which requires only one reference point (RP) by exploiting useful information given by the position change of a mobile terminal (MT), or the user movement. This method can minimize the number of reference points required in a localization system or navigation service and reduce system implementation cost. Analytical result shows that the user position can be thus derived and given in simple closed-form expression. Finally, simulation is conducted to demonstrate its effectiveness under noisy environment.

I. INTRODUCTION

Positioning systems are crucial to today's digital society. They provide geographic information about devices that facilitates many human activities. For instance, vehicle navigation systems are indispensable for drivers in big cities. Some location-based services are deployed in commercial malls so that customers can get navigation while walking in complex environment and can receive promotion advertisement from shops. The market of indoor and outdoor location-based services has grown rapidly in the last decade.

Global positioning system (GPS) is very popular and widely used for user localization. When line-of-sight to at least four GPS satellites is available, location (latitude, longitude, and elevation) and timing information can be obtained. Although GPS is very convenient in outdoors, its quality is susceptible to weather conditions, for example when sky view is poor due to fog, rain, cloud, etc., or being blocked by tall buildings in urban areas. These issues can significantly degrade the accuracy. As expected, GPS is not for indoors due to the lack of line-of-sight. There also exist cellular-based positioning systems [START_REF] Gartner | Location Based Services and TeleCartography II: From Sensor Fusion to Context Models[END_REF] which are built on measuring signal strength from three or more base stations for tracking mobile user's location. However, these solutions also do not work well for indoors.

Various indoor positioning systems have been developed, see e.g., [START_REF] Liu | Survey of wireless indoor positioning techniques and systems[END_REF]- [START_REF] Kushki | WLAN Positioning Systems: Principles and Applications in Location-Based Services[END_REF]. They can be categorized into network-based or non-network-based solutions. The network-based approach, which takes advantages of existing network infrastructure such as wireless local area networks (WLANs), without demanding new infrastructure, can maintain low deployment cost. The non-network-based approach uses dedicated positioning infras-tructure and often can provide higher reliability but with extra cost. For example, solutions based on ultrasound and infrared have high deployment cost. One may also consider simple proximity-based solution like iBeacon [START_REF] Cavallini | iBeacons Bible[END_REF] which however is only able to offer an approximate location. Some systems consider using visible light to construct an indoor positioning system with high accuracy [START_REF] Yasir | Indoor positioning system using visible light and accelerometer[END_REF]. A good positioning system should be cost-effective and also be able to offer high accuracy.

Constructing an efficient and simple positioning system is always challenging. Technically, it would depend on: 1) the number of available reference points (RPs); 2) the technologies used (e.g., RF-based, ultrasound, infrared, etc.) and; 3) the characteristics of the environment. In this study, we propose a geometry-based positioning method which can determine user position by only using one RP and exploiting user's simple movement, for instance walking or waving user's hand-held device, and some simple information. As the solution requires only one RP and can provide either exact result in noiseless environment or accurate positioning in noisy condition, our approach brings competitive advantages compared to other methods, thanks to its simplicity and effectiveness. Meanwhile, the method is interesting and may have a high potential to improve today's technology or existing solutions.

II. RELATED WORK

Indoor positioning problem has attracted a lot of interest over years [START_REF] Gu | A survey of indoor positioning systems for wireless personal networks[END_REF], [START_REF] Harle | A survey of indoor inertial positioning systems for pedestrians[END_REF]. Studies have been done extensively and many possible solutions have been proposed. There are mainly four major approaches to solve this problem: triangulation, fingerprinting, scene analysis, and proximity.

Triangulation is used to estimate the position of a user or mobile terminal (MT) if the geographical coordinates of the RPs are known and assume that the MT is capable of measuring the distance between itself and the RPs. A priori, this method requires three RPs to construct a distinct geometric intersection of three circles, which indicates the position of the MT. The principle is illustrated in Fig. 1(a) (or see e.g., [START_REF] Liu | Survey of wireless indoor positioning techniques and systems[END_REF]). Note that not all schemes based on triangulation requires three circles (see e.g., Fig. 1(b) and (c)). For instance, given angleof-arrival (AoA) information, using only one RP is sufficient to locate the MT.

Fingerprinting [START_REF] Mirowski | Probabilistic radio-frequency fingerprinting and localization on the run[END_REF] is to estimate device position by using pre-measured location-related data. This method consists of two phases: an offline training phase and an online position estimation phase. In the offline phase, location-related data is collected at different positions in the area. During the online position determination phase, real-time location-related data is measured and then matched with the set of data gathered during the offline phase to estimate the device's location.

Scene analysis [START_REF] Gu | A survey of indoor positioning systems for wireless personal networks[END_REF] is a localization method based on a set of images or scenes received by one or multiple cameras. This approach in principle does not require user (to be tracked) to carry any extra device. However, the solution is usually expensive because it requires one or many cameras to perform tracking and may prone to a high computation cost due to image or video processing.

Proximity detects if a MT is nearby or for example in the coverage area of a RP. However, it is hard to provide accurate position with high reliability.

Each of the above method also has some variants or hybrid schemes. Our proposed geometry-based solution is built on triangulation. We will explain and discuss in comparison other methods stemmed from this branch. The cost and accuracy of triangulation method primarily rely on the number of RPs required. Traditionally, one would need at least three RPs to determine the position of the MT.

Figure 1(b) shows a variant of traditional triangulation method, which requires two RPs and the last estimated previous position of the mobile terminal so as to eliminate one of the two intersection points of the two circles constructed by the two RPs. In such case, the location closer to the last estimated position would be selected. Or, the system has to be able to predict user mobility pattern in order to select one. Note that this method still requires more than one RP. A variant of the above triangulation method is to use only one RP but requires the information of angle-of-arrival (AoA) provided by an array of antennas either implemented in the user terminal (MT) or at the RP [START_REF] Peng | Angle of arrival localization for wireless sensor networks[END_REF], see Fig. 1(c). However, such an array of antennas is often costly and cumbersome.

III. SYSTEM DESIGN

Here, we propose a new method called "Two-Step Movement (2SM)". It aims to improve the classical triangulation approach and requires only one RP. It makes use of the changes in the position of the MT relative to the RP. The changes is caused by either active movement (e.g., a user may wave his/her MT to assist) or natural movement (e.g., the user is walking or moving). Therefore, 2SM turns out to have a competitively low deployment cost and without extra or expensive tracking hardware such as antenna array and is able to determine user position in exact closed-form solution. The simplicity and effectiveness would highly facilitate practical indoor positioning systems. Table I gives a comparison of the above methods and outlines their key difference. In our proposed 2SM method, the MT is assumed be capable to measure its movement using its embedded sensors and software.

A. One-Step Movement (1SM)

Our method exploits useful information generated by user movement. For the sake of simplicity, the 2SM is presented as a combination of two One-Step Movements.

One-Step Movement (1SM) relies on one position change (one move) to identify the two possible locations (position candidates) of the MT with the following assumptions:

• The position of the RP is known.

• The MT is capable of measuring the distance between itself and the RP. • The MT is capable of measuring the distance and the angle (direction) of its movement. Figure 2 illustrates the system design, where

• A is the RP with a known position (x A , y A ).

• B is the initial position of the MT that is unknown and yet to be found. It is denoted by coordinates (x B , y B ). • C is the position of MT right after the first movement, (x C , y C ), which is also unknown. • MT is capable of measuring the distance between itself and the RP. That is, the distances AB and AC are given for example by measuring the received signal strength or standard techniques. • MT is capable of measuring the distance and the angle of its movement, thus BC and the angle α ∈ (0, 2π] (with respect to the positive x-axis) are also measurable.

Theorem 1 Suppose that A(x A , y A ), AB, BC, AC, and α are known, the One-Step Movement (1SM) will give two estimated locations, denoted by generic point B(x B , y B ), whose x and y coordinates satisfy:

x B cos α + y B sin α = x A cos α + y A sin α - (AB 2 + BC 2 -AC 2 ) 2BC . ( 1
)
Proof: Using Fig. 2, from the two measured distances AB and AC, the equations of the two circles centered at A(x A , y A ) on which the MT probably lies can be expressed as

(x B -x A ) 2 + (y B -y A ) 2 = AB 2 (x C -x A ) 2 + (y C -y A ) 2 = AC 2 (2)
where

x C = x B + BC cos α, y C = y B + BC sin α. (3)
From (2), we have

AB 2 -AC 2 = (x B -x C )(x B + x C -2x A ) +(y B -y C )(y B + y C -2y A ). (4)
A(x A , y A ) B(x B , y B ) C(x C , y C ) (a) Three RPs A(x A , y A ) B(x B , y B )
The last estimated position of the MT Substitute x C and y C in (3) to (4), we can have

AB 2 -AC 2 = -BC cos α(2x B + BC cos α -2x A )
-BC sin α(2y B + BC sin α -2y A ) which can be re-written as

AB 2 + BC 2 -AC 2 = -2BC(x B cos α -x A cos α +y B sin α -y A sin α). Hence, x B cos α + y B sin α = x A cos α + y A sin α - (AB 2 + BC 2 -AC 2 ) 2BC .
We can solve (1) as follows • If sin α = 0, thus cos α = ±1, (1) becomes:

x B = x A ± (AB 2 + BC 2 -AC 2 ) 2BC .
It is then straightforward to compute the values of x B and y B , by substituting the value of x B to (2).

• If sin α = 0, by dividing (1) by sin α, we have:

y B = -cot αx B + x A cot α + y A - AB 2 + BC 2 -AC 2 2BC sin α . Let b = x A cot α + y A -(AB 2 + BC 2 -AC 2 )/(2BC sin α)
and a = -cot α. We see that now y B can be expressed as a function of x B such that y B = ax B + b. Substituting y B to the first equation of (2), we have

(x B -x A ) 2 + (ax B + b -y A ) 2 = AB 2 .
Then

(1+a 2 )x 2 B -2x B (x A -a(b-y A ))+x 2 A +(b-y A ) 2 -AB 2 = 0.
(5) The above quadratic equation ( 5) can be solved easily. Algorithm 1 shows in detail how to perform 1SM. It outputs two points B1(x B1 , y B1 ) and B2(x B2 , y B2 ), which are the possible solution of B.

Remark 1

It is clear that one of the two points, B1(x B1 , y B1 ) and B2(x B2 , y B2 ), must be the position of the MT (or both of them are, if B1 and B2 are identical).

B. Two-Step Movement (2SM)

After the first movement, we have two possible locations of the MT given by 1SM using Algorithm 1, but cannot determine which one is the true location. We need to resolve this ambiguity. It is natural to think about performing an additional movement. The basic idea is simple: a Two-Step Movement (2SM) is a combination of two consecutive 1SM's where each move gives two possible positions (in which one of these two positions must be the true position). It is clear that by comparing the results of two 1SM's, we can determine the location of the MT, given that the results of the two 1SM's are not redundant.

Fig. 3 depicts how 2SM works. The MT makes the second movement from C to D in the direction of angle β, which

Algorithm 1 One-Step Movement algorithm

Require: A(x A , y A ), AB, AC, BC, α;

1: function ONESTEP(A(x A , y A ), AB, AC, BC, α)

2:
if sin α == 0 then 3:

if cos α == 1 then 4:

x B = x A -(AB 2 + BC 2 -AC 2 )/(2BC); 5:
else 6:

x B = x A + (AB 2 + BC 2 -AC 2 )/(2BC); 7:
end if 8:

y B1 = y A + AB 2 -(x B -x A ) 2 ;
9: a = -cot α;

y B2 = y A -AB 2 -(x B -x A ) 2 ;
14: b = x A cot α + y A -(AB 2 + BC 2 - AC 2 )/(2BC sin α); 15: ⊲ Compute x B , y B ; 16: ∆ = (x A -a(b -y A )) 2 -(1 + a 2 )(x 2 A + (b - y A ) 2 -AB 2 ); 17: x B1 = (x A -a(b -y A ) + √ ∆)/(1 + a 2 );
18:

y B1 = ax B1 + b; 19: x B2 = (x A -a(b -y A ) - √ ∆)/(1 + a 2 ); 20: y B2 = ax B2 + b; 21: return {B1(x B1 , y B1 ), B2(x B2 , y B2 )}; 22:
end if 23: end function is measured from the positive x-axis counter-clockwise. The distance CD and β are known by the MT, whereas the distance AD from the MT to the RP is measured from the received signal strength by standard techniques. The underlying idea is that, we now consider the movement of 2SM case similarly as that of 1SM case in which the starting point is now B and the ending point is D. We can compute the distance BD and the angle γ analytically (see Algorithm 2: line 5-10) and then use the method of Algorithm 1 to determine B. Algorithm 2 details how 2SM works. By comparing the results from the two 1SM's computation, we determine the location of the MT.

Remark 2 Note that the directions of the two movements should not be in parallel, i.e., β = α and β = α±π, otherwise the ambiguity cannot be resolved since the system of equations generated by the second movement would be equivalent to that of the first one.

In practice with estimation error or system imperfection, say the existence of noise, we may not have a common solution from the two 1SM's computation. In this case, the first movement may give us two possible solutions denoted by B1(x B1 , y B1 ) and B2(x B2 , y B2 ), but the second movement may give us another two possible solutions denoted by B3(x B3 , y B3 ) and B4(x B4 , y B4 ). However {B1, B2} and {B3, B4} may not have a common point as shown in Fig. 4. To solve this problem, we can choose the pair of points that have the smallest distance, i.e., solving min{d(P 1, P 2)|P 1 = P 2}, for P 1, P 2 ∈ {B1, B2, B3, B4} where d(P 1, P 2) denotes the Euclidean distance of points P 1 and P 2. After that, we take their mean (e.g., the mid-point of B1 and B3 in Fig. 4) as the estimate of the MT's position for minimizing the error. In general, one can formulate it as an optimization problem and find the optimal result.

B2(x B2 , y B2 ) B4(x B4 , y B4 ) B1(x B1 , y B1 ) B3(x B3 , y B3 )
Fig. 4: Ambiguity elimination in case of noise.

IV. SIMULATION

Simulation is performed to investigate the performance of the proposed scheme (2SM) under noisy environment. The RP is placed at the center of a room, i.e., A = (0, 0). We are going to determine the MT's location, denoted by B(x B , y B ), which is randomly distributed in the room. In the following analysis, we consider three distances (1, 5, and 10 meters) between the MT (denoted by B, in Fig. 3) and the RP. Also, we assume that the direction from B to A is uniformly distributed in (0, 2π].

For a given AB, the movement from B to C or from C to D is equal to 0.1, 0.2, and 0.5 times of AB. The directions of the movement, i.e., α and β, are uniformly distributed in Algorithm 2 Two-Step Movement algorithm Require: A(x A , y A );

1: function TWOSTEP(A(x A , y A ))

2:

MT makes the first movement from B to C; measure AB, AC, BC, α;

3:

{B1(x B1 , y B1 ), B2(x B2 , y B2 )} = ONESTEP(A(x A , y A ), AB, AC, BC, α); ⊲ Obtain B1 and B2

4:

MT makes the second movement from C to D; measure CD, AD, β; make sure that β = α and β = α ± π ;

BD = √ X 2 + Y 2 ; 8:
cos γ = X/BD;

9:
sin γ = Y /BD; 10:

Compute γ ∈ (0; 2π] from cos γ and sin γ;

11: {B3(x B3 , y B3 ), B4(x B4 , y B4 )} = ONESTEP(A(x A , y A ), AB, AD, BD, γ); ⊲ Obtain B3 and B4 12: B(x B , y B ) = {B1(x B1 , y B1 ), B2(x B2 , y B2 )} ∩ {B3(x B3 , y B3 ), B4(x B4 , y B4 )};
⊲ Determine MT location B(x B , y B ) from the set of B1, B2, B3 and B4; 13: return B(x B , y B ); 14: end function (0, 2π]. Estimation error to the measurement of distances AB, AC, AD, and BC, is considered to be bounded in [-1%, 1%], [-2%, 2%], and [-5%, 5%], for comparison. We use e d to denote the above bound such that e d = 1%, 2%, and 5%, respectively. Estimation error to the measurement of angles α and β is considered to be bounded in

[-1 • , 1 • ], [-2 • , 2 • ],
and [-5 • , 5 • ]. The bound on the angle measurement error is denoted by e a such that e a = 1, 2, and 5 degrees. For each (e d , e a ) setup, the errors are randomly generated to corrupt the proposed algorithm in determining B(x B , y B ). Results are shown in Figs. 5,6, and 7 for BC = CD is equal to 0.1×AB, 0.2 × AB and 0.5 × AB, respectively. Note that each curve in the figures is obtained by 10,000 runs. During these runs, we observe that about 10% of the time the system fails to find the MT position (i.e., the quadratic equation ( 5) has no solution since ∆ in Algorithm 1 is negative) due to the noise (which would be accumulated to ∆). We find that when ∆ < 0, the system is indeed heavily corrupted. In this case, the movement is considered as bad and is not used to locate the MT. Note that it would be interesting to derive the position of the MT even when ∆ < 0 and see how to extract useful information to optimize results. This is left as future work.

As shown in Figs. 567, the estimation error in determining the position of the MT increases as e d increases. Note that the estimation error is defined by the distance between the real position of the MT and the result given by Algorithm 2. Clearly, e d = 1% (curves in "red") results in smaller estimation error than that e d = 3% or e d = 5% (curves in "blue" and "black", respectively) makes, given that e a is the same.

As expected, the estimation error in determining the position of the MT also increases as e a increases. However, when e d is relatively large (5%), the impact of the considered e a is relatively less significant. This can be clearly shown by Figs. 6 and7. Roughly speaking, e d is more dominating.

In comparing Figs. 567, we observe that when increasing BC and CD from 0.1 × AB to 0.5 × AB, the error in estimating the position of the MT decreases quite substantially. In Fig. 6 and7, the curves shift to the left. The distance of the movement is a significant factor. We can improve the system performance by requiring a longer movement distance. However, a longer distance may be less favorable in some scenarios. Furthermore, from the obtained simulation results which are not shown in this paper, we see that the improvement is indeed decreasing and starts to get flat at 0.5 × AB.

Table II shows the average error in determining the position of the MT under different AB (at 1, 5, and 10 meters) and various BC, CD, and noise levels (e d , e a ). For AB = 5, the results are plotted in Figs. 567. For AB = 1 and 10, the results have characteristics very similar to those in Figs. 5-7 so that they are not shown in this paper. Comparing the results at AB = 1, 5, and 10, we see that the magnitude of the error increases roughly proportional to AB. It is clear that the estimation error is minimized when (e d , e a ) are small and the movement distance is relatively large. Roughly speaking, if BC = CD = 0.5 × AB, the performance is quite desirable for e d ≤ 2% and e a ≤ 5%. When the movement distance is at the level of 0.2 × AB, the same performance can be achieved for a smaller e d ≤ 1%. The average error can be limited to within about 10% of AB. In the best case, the average error can be less than 5% of AB.

V. CONCLUSION

In this paper, we have proposed a new method called Two-Step Movement (2SM) to estimate the position of MT. It requires only one reference point (RP) by exploiting useful information given by the position change of the MT or user movement. One can therefore reduce the number of RPs required and also the system cost. Analytical result shows that the user position can be derived and given in simple closed-form expression with low complexity. Simulation is conducted to study its performance under noisy environment. It is possible to achieve average error within about 10% of the distance between the RP and MT, or even less. Note that further analysis of noise impact and issues related to reflection and refraction of signals are important to improve the proposed method. Our method, thanks to the reliance on a single reference point, makes a lot of sense in the context of Internet of Things (IoT) such as home or business office area. It should be also noted that our method can be easily extended to localization in 3D coordinates and to device-todevice (D2D) applications in which both devices could be mobile. The practical implementation is left as future work.

Fig. 1 :

 1 Fig. 1: Positioning techniques using different number of reference points (RPs).

Fig. 2 :

 2 Fig. 2: One-Step Movement (1SM).

10 :⊲

 10 return {B1(x B , y B1 ), B2(x B , y B2 )}; Pre-compute a, b such that y B = ax B + b; 13:

Fig. 3 :

 3 Fig. 3: Two-Step Movement (2SM).

5 :X

 5 = BC cos α + CD cos β; ⊲ The change in x-coordinate after the second move 6: Y = BC sin α + CD sin β; ⊲ The change in y-coordinate after the second move 7:

  e d = 1% e d = 1% e d = 1% e d = 2% e d = 2% e d = 2% e d = 5% e d = 5% e d = 5% ea = 1°ea = 2°ea = 5°ea = 1°ea = 2°ea = 5°ea = 1°ea = 2°ea = 5°A B = 1 (BC=CD= 0

Fig. 5 :

 5 Fig. 5: Resulting error when AB = 5 meters, BC = CD = 0.1 × AB under various (e d , e a ).

Fig. 6 :

 6 Fig. 6: Resulting error when AB = 5 meters, BC = CD = 0.2 × AB under various (e d , e a ).

Fig. 7 :

 7 Fig. 7: Resulting error when AB = 5 meters, BC = CD = 0.5 × AB under various (e d , e a ).

TABLE I :

 I Requirement comparison between triangulation and proposed 2-Step Movement (2SM) method.

TABLE II :

 II Average error (in meter) under various AB, BC, CD, and noise levels (e d , e a ).

ACKNOWLEDGMENT

The work presented in this paper has been carried out at LINCS (www.lincs.fr). The authors would like to thank Fabien Mathieu and Philippe Jacquet for their discussion and valuable comments.