Pat Fitzsimmons

Yves Le

Jan Jay Rosen
email: jrosen30@optimum.net

Yves Le Jan
email: yves.lejan@math.u-psud.fr

Jay Rosen

Loop measures without transition probabilities

Keywords:

Introduction

To the best of our knowledge, loop measures first appeared in the work of Symanzik on Euclidean quantum field theory, [START_REF] Symanzyk | Euclidean quantum field theory[END_REF], where they are referred to as 'blob measures'. They next appear in the work of Lawler and Werner,[6]. In both works, the loop measure is that associated with Brownian motion. In [7] loop measures associated with a large class of Markov processes are defined and studied. In all these cases it is assumed that the underlying Markov process has transition densities. The goal of this paper is to define and study loop measures for Markov processes without transition densities.

Let X = (Ω, F t , X t , θ t , P x ) be a transient Borel right process [START_REF] Sharpe | General theory of Markov processes[END_REF] with state space S, which we assume to be locally compact with a countable base. We use the canonical representation of X in which Ω is the set of right continuous paths paths ω : [0, ∞) → S ∆ = S ∪ ∆ with ∆ / ∈ S, and is such that ω(t) = ∆ for all t ≥ ζ = inf{t > 0 |ω(t) = ∆}. Set X t (ω) = ω(t).

Let m be a Borel measure on S which is finite on compact sets. We assume that with respect to m, X has strictly positive potential densities u α (x, y), α ≥ 0, which satisfy the resolvent equations. We set u(x, y) = u 0 (x, y), and assume that u(x, y) is excessive in x for each fixed y.

Let h z (x) = u(x, z). If we assume that u is finite, then the h z -transform of X is a right process on S, see [10, Section 62], with probabilities P x/hz . Let 1 Q z,z = u(z, z)P z/hz . We can then define the loop measure as

µ (F ) = Q z,z F ζ dm(z), (1.1)
for any F measurable function F . Loop measures for processes with finite potential densities but without transition densities are discussed in [START_REF] Fitzsimmons | Markovian loop soups: permanental processes and isomorphism theorems[END_REF]. In the present paper we assume that the potential densities u(x, y) are infinite on the diagonal, but finite off the diagonal. In this case the construction of the measures P z/hz given in [START_REF] Sharpe | General theory of Markov processes[END_REF] breaks down. Assuming that all points are polar, we show how to construct a family of measures Q z,z , z ∈ S, which generalize the measures Q z,z = u(z, z)P z/hz in the case of finite u(x, y).

After constructing Q z,z , z ∈ S and defining the loop measure µ using (1.1), we show how to calculate some important moments. We assume that sup x K (u(x, y) + u(y, x)) 2 dm(y) < ∞ (1.2) for any compact K ⊆ S. For exponentially killed Brownian motion in R d this means that d ≤ 3.

Theorem 1.1 For any k ≥ 2, and bounded measurable functions f 1 , . . . , f k with compact support

µ   k j=1 ∞ 0 f j (X t j ) dt j   (1.3) = π∈P k u(x 1 , x 2 )u(x 2 , x 3 ) • • • u(x k , x 1 ) k j=1 f π j (x j ) dm(x j ),
where P k denotes the set of permutations of [1, k] on the circle. (For example, (1, 2, 3), (3, 1, 2) and (2, 3, 1) are considered to be one permutation π ∈ P 3 .)

Our assumption (1.2) will guarantee that the right hand side of (1.3) is finite. Note that if k = 1 our formula would give

µ ∞ 0 f (X t ) dt j = u(x, x)f (x) dm(x) = ∞, (1.4)
for any f ≥ 0, by our assumption that the potentials u(x, y) are infinite on the diagonal.

For f 1 , . . . , f k as above consider more generally the multiple integral

M f 1 ,...,f k t = π∈T k 0≤r 1 ≤•••≤r k ≤t f π(1) (X r 1 ) • • • f π(k) (X r k ) dr 1 • • • dr k , (1.5)
where T k denotes the set of translations π of [1, k] which are cyclic mod k, that is, for some i, π(j) = j + i, mod k, for all j = 1, . . . , k. In the proof of Theorem 1.1 we first show that

µ M f 1 ,...,f k ∞ = u(x 1 , x 2 )u(x 2 , x 3 ) • • • u(x k , x 1 ) k j=1 f j (x j ) dm(x j ). (1.6) (1.3) will then follow since k j=1 ∞ 0 f j (X t j ) dt j = π∈P k M fπ 1 ,...,fπ k ∞ . (1.7)
There is a related measure which we shall use which gives finite values even for k = 1. Set ν(F ) = Q z,z (F ) dm(z).

(1.8)

Assume that for any α > 0, any compact K ⊆ S, and any K which is a compact neighborhood of K

sup z∈ K c , x∈K u α (z, x) < ∞.
(1.9)

Theorem 1.2 For any k ≥ 1, α > 0, and bounded measurable functions

f 1 , . . . , f k with compact support ν   k j=1 ∞ 0 f j (X t j ) dt j e -αζ   (1.10) = π∈P k u α (z, x 1 )u α (x 1 , x 2 ) • • • u α (x k , z) k j=1 f π j (x j ) dm(x j ) dm(z),
where P k denotes the set of permutations of [1, k], and both sides are finite.

We call µ the loop measure of X because, when X has continuous paths, µ is concentrated on the set of continuous loops with a distinguished starting point (since Q x,x is carried by loops starting at x). Moreover, in the next Theorem we show that it is shift invariant. More precisely, let ρ u denote the loop rotation defined by

ρ u ω(s) = ω(s + u mod ζ(ω)), if 0 ≤ s < ζ(ω) ∆, otherwise.
Here, for two positive numbers a, b we define a mod b = a -mb for the unique positive integer m such that 0 ≤ a -mb < b.

For the next Theorem we need an additional assumption: for any δ > 0 and compact K ⊆ S K P δ (z, dx)u(x, z) dm(z) < ∞.

(1.11)

Theorem 1.3 µ is invariant under ρ u , for any u.

Note that if we have transition densities p δ (z, x) then

K P δ (z, dx)u(x, z) dm(z) = K p δ (z, x)u(x, z) dm(x) dm(z) (1.12) = K ∞ δ p t (x, x) dt dm(x).
In our work on processes with transition densities, it was always assumed that sup x ∞ δ p t (x, x) dt < ∞ for any δ > 0, which indeed gives (1.11). For the next Theorem we assume that the measure m is excessive. With this assumption there is always a dual process X (essentially uniquely determined), but in general it is a moderate Markov process. We assume that the measures Û (•, y) are absolutely continuous with respect to m for each y ∈ S.

For CAF's L ν 1 t , . . . L ν k t with Revuz measures ν 1 , . . . , ν k , let

A ν 1 ,...,ν k t = π∈T k 0≤r 1 ≤•••≤r k ≤t dL ν π(1) r 1 • • • dL ν π(k) r k . (1.13)
We refer to A ν 1 ,...,ν k t as a multiple CAF.

Theorem 1.4 For any k ≥ 2, and any CAF's L ν 1 t , . . . L ν k t with Revuz measures ν 1 , . . . , ν k ,

µ (A ν 1 ,...,ν k ∞ ) = u(x 1 , x 2 )u(x 2 , x 3 ) • • • u(x k , x 1 ) k j=1 dν j (x j ).
(1.14)

and µ   k j=1 L ν j ∞   = π∈P k u(x 1 , x 2 )u(x 2 , x 3 ) • • • u(x k , x 1 ) k j=1 dν π j (x j ). (1.15)
The finiteness of the right hand side of (1.15) will depend on the potential densities u(x, y) and the measures ν 1 , . . . , ν k . For a more thorough discussion see [8, (1.5)] and the paragraph there following (1.5).

With the results of this paper, most of the results of [START_REF] Fitzsimmons | Markovian loop soups: permanental processes and isomorphism theorems[END_REF]8,[START_REF] Le Jan | Intersection local times, loop soups and permanental Wick powers[END_REF] on loop measures, loop soups, CAF's and intersection local times will carry over to processes without transition densities.

2 Construction of Q z,z
Let us fix z ∈ S and consider the excessive function h z (x) := u(x, z), finite and strictly positive on the subspace S z := {x ∈ S : x = z}. Doob's h-transform theory yields the existence of laws P x,z , x ∈ S z , on path space under which the coordinate process is Markov with transition semigroup

P z t (x, dy) := P t (x, dy) h z (y) h z (x) . (2.1)
See, for example, [10, pp. 298-299]. Now consider the family of measures η z t (dx) := P t (z, dx)h z (x).

(2.2)

Since we assume that the singleton {z} is polar, the transition semigroup (P t ) will not charge {z}, so these may be viewed as measures on S z or on S.

Adopting the latter point of view, it is immediate that (η z t ) t>0 is an entrance law for (P z t ). There is a general theorem guaranteeing the existence of a right process with one-dimensional distributions (η z t ) and transition semigroup (P z t ); see [START_REF] Getoor | Constructing Markov Processes with Random Times of Birth and Death[END_REF]Proposition (3.5)]. The law of this process is the desired Q z,z . Aside from the entrance law identity η z t P z s = η z t+s , their result only requires that each of the measures η z t be σ-finite, which is clearly the case in the present discussion.

With this we immediately obtain, for 0

< t 1 < • • • < t k , Q z,z   k j=1 f j (X t j )   = η z t 1 (dx 1 ) k j=2 P z t j -t j-1 (x j-1 , dx j ) k j=1 f j (x j ) (2.3) = k j=1 P t j -t j-1 (x j-1 , dx j ) k j=1 f j (x j )u(x k , z)
with t 0 = 0 and x 0 = z. Hence

Q z,z   0<t 1 <•••<t k <∞ k j=1 f j (X t j ) dt j   (2.4) = u(z, x 1 )u(x 1 , x 2 ) • • • u(x k , z) k j=1 f j (x j ) dm(x j ), so that Q z,z   k j=1 ∞ 0 f j (X t j ) dt j   (2.5) = π∈P k u(z, x 1 )u(x 1 , x 2 ) • • • u(x k , z) k j=1 f π j (x j ) dm(x j ). Returning to (2.3) with 0 < t 1 < • • • < t k and using the fact that ζ > t k implies that ζ = t k + ζ • θ t k we have Q z,z   k j=1 f j (X t j )e -αζ   (2.6) = Q z,z   k j=1 f j (X t j )e -αt k e -αζ • θ t k   = k j=1 P α t j -t j-1 (x j-1 , dx j ) k j=1 f j (x j )h z (x k )P x k ,z e -αζ .
Note that by (2.1) and the fact that X has α-potential densities for all α ≥ 0

P x k ,z ∞ 0 e -αt 1 {S} (X t ) dt = ∞ 0 e -αt P x k ,z 1 {S} (X t ) dt (2.7) = ∞ 0 e -αt S P t (x k , dy) h z (y) h z (x k ) dt = u α (x k , y) h z (y) h z (x k ) dm(y).
Combining this with our assumption that the α-potential densities satisfy the resolvent equation we see that

h z (x k )P x k ,z e -αζ = h z (x k )P x k ,z 1 -α ζ 0 e -αt dt (2.8) = h z (x k ) -αh z (x k )P x k ,z ∞ 0 e -αt 1 {S} (X t ) dt = u(x k , z) -α u α (x k , y)u(y, z) dm(y) = u α (x k , z).
Using this in (2.6) we obtain

Q z,z   k j=1 f j (X t j )e -αζ   (2.9) = k j=1 P α t j -t j-1 (x j-1 , dx j )u α (x k , z) k j=1 f j (x j ).
We then have

Q z,z   e -αζ 0<t 1 <•••<t k <∞ k j=1 f j (X t j ) dt j   (2.10) = u α (z, x 1 )u α (x 1 , x 2 ) • • • u α (x k , z) k j=1 f j (x j ) dm(x j ),
and consequently

Q z,z   e -αζ k j=1 ∞ 0 f j (X t j ) dt j   (2.11) = π∈P k u α (z, x 1 )u α (x 1 , x 2 ) • • • u α (x k , z) k j=1 f π j (x j ) dm(x j ).
3 The loop measure and its moments Set

µ(F ) = Q z,z F ζ dm(z). (3.1)
Proof of Theorem 1.1: We use an argument from the proof of [4, Lemma 2.1], which is due to Symanzik, [START_REF] Symanzyk | Euclidean quantum field theory[END_REF].

It follows from the resolvent equation that the potential densities u β (x, y) are continuous and monotone decreasing in β, for x = y. Using this together with the resolvent equation and the monotone convergence theorem we obtain that for

x k = x 1 S u α (x k , z)u α (z, x 1 ) dm(z) = - d dα u α (x k , x 1 ). (3.2)
Hence using (2.10)

Q z,z e -αζ M f 1 ,...f k ∞ dm(z) (3.3) = - π∈T k u α (x 1 , x 2 )u α (x 2 , x 3 ) • • • u α (x k-1 , x k ) d dα u α (x k , x 1 ) k j=1 f π j (x j ) dm(x j ) = - d dα u α (x 1 , x 2 )u α (x 2 , x 3 ) • • • u α (x k-1 , x k )u α (x k , x 1 ) k j=1 f j (x j ) dm(x j ).
For the last step we used the product rule for differentiation and the fact that in the middle line we are summing over all translations mod k. Since, as mentioned, u α (x, y) is monotone decreasing in α for x = y, v(x, y) = lim Proof of Theorem 1.2: The formula (1.10) follows immediately from (2.11). When k ≥ 2, the right hand side of (1.10) can be shown to be finite by repeatedly using the Cauchy-Schwarz inequality, our assumption (1.2) and the fact that u α (x, z) is integrable in z for any α > 0. When k = 1, if K is a compact set containing the support of f 1 and K is a compact neighborhood of K, then

u α (z, x)u α (x, z)f 1 (x) dm(x) dm(z) = K u α (z, x)u α (x, z)f 1 (x) dm(x) dm(z) + K c u α (z, x)u α (x, z)f 1 (x) dm(x) dm(z). Using (1.2) K u α (z, x)u α (x, z)f 1 (x) dm(x) dm(z) (3.6) ≤ m( K) sup z u α (z, x)u α (x, z)f 1 (x) dm(x) < ∞,
and using (1.9)

K c u α (z, x)u α (x, z)f 1 (x) dm(x) dm(z) (3.7) ≤ C u α (x, z) dm(z) f 1 (x) dm(x) < ∞.

Subordination

The basic idea in our proof that the loop measure is shift invariant is to show that the loop measure can be obtained as the 'limit' of loop measures for processes with transition densities. These processes will be obtained from the original process by subordination. We consider a subordinator T t which is a compound Poisson process with Levy measure cψ so that

E x (f (X Tt )) = ∞ j=1 (ct) j j! e -ct ∞ 0 E x (f (X s )) ψ * j (ds). (4.1)
If we take ψ to be exponential with parameter θ, then ψ * j (ds) = s j-1 θ j Γ(j) e -sθ ds so that we have

E x (f (X Tt )) = ∞ j=1 (ct) j j! e -ct ∞ 0 E x (f (X s )) s j-1 θ j Γ(j) e -sθ ds. (4.2)
Hence the subordinated transition semigroup

P t (x, dy) = ∞ j=1 (ct) j j! e -ct ∞ 0 P s (x, dy) s j-1 θ j Γ(j) e -sθ ds. (4.3)
Noting that

∞ 0 P s (x, A)s j-1 e -sθ ds = d j-1 dθ j-1 ∞ 0 P s (x, A)e -sθ ds = d j-1 dθ j-1 U θ (x, A), (4.4)
we see that P t (x, dy) is absolutely continuous with respect to the measure m on S, and we can choose transition densities p t (x, y).

From now on we take θ = c = n, and use (n) as a superscript or subscript to denote objects with respect to the subordinated process, denoted by

X (n) t . Lemma 4.1 u α (n) (x, y) = 1 (1 + α/n) 2 u α/(1+α/n) (x, y). (4.5)
In particular,

u α (n) (x, y) ≤ u (n) (x, y) = u(x, y). (4.6)
Proof:

U α (n) (x, dy) = ∞ 0 e -αt P (n) t (x, dy) dt (4.7) = 1 α + n ∞ j=1 ∞ 0 n α + n j P s (x, dy) s j-1 n j Γ(j) e -sn ds. = n 2 (α + n) 2 ∞ 0 P s (x, dy)   ∞ j=1 s j-1 (n 2 /(α + n)) (j-1) Γ(j) e -sn   ds = n 2 (α + n) 2
∞ 0 e sn 2 /(α+n) e -sn P s (x, dy) ds

= n 2 (α + n) 2 ∞ 0 e -s n-n 2 α+n P s (x, dy) ds = n 2 (α + n) 2 u nα/(α+n) (x, y) dm(y) = 1 (1 + α/n) 2 u α/(1+α/n) (x, y) dm(y).
Lemma 4.2 For any α, α j , j = 1, . . . , k and continuous compactly supported f j , j = 1, . . . , k

lim n→∞ R k + k j=1 e -α j t j Q z,z (n)   k j=1 f j (X (n) t j )e -αζ   k j=1 dt j (4.8) = R k + k j=1 e -α j t j Q z,z   k j=1 f j (X t j )e -αζ   k j=1 dt j .
Proof: Recall (2.9). For 0

< t 1 < • • • < t k Q z,z   k j=1 f j (X t j )e -αζ   (4.9) = k j=1 P α t j -t j-1 (x j-1 , dx j )u α (x k , z) k j=1 f j (x j )
with the corresponding

Q z,z (n)   k j=1 f j (X (n) t j )e -αζ   (4.10) = k j=1 P (n),α t j -t j-1 (x j-1 , dx j )u α (n) (x k , z) k j=1 f j (x j ).
Using (4.10) we see that

R k + k j=1 e -α j t j Q z,z (n)   k j=1 f j (X (n) t j )e -αζ   k j=1 dt j (4.11) = π∈P k {0≤t 1 ≤•••≤t k ≤∞} k j=1 e -α j t j Q z,z (n)
  k j=1 f π(j) (X (n) t j )e -αζ   k j=1 dt j = π∈P k {0≤t 1 ≤•••≤t k ≤∞} k j=1 P (n),α+ k l=j α l t j -t j-1 (x j-1 , dx j ) u α (n) (x k , z) k j=1 f π(j) (x j ) k j=1 dt j = π∈P k k j=1 U α+ k l=j α l (n) (x j-1 , dx j )u α (n) (x k , z) k j=1 f π(j) (x j ) = π∈P k k j=1 u α+ k l=j α l (n) (x j-1 , x j )u α (n) (x k , z) k j=1 f π(j) (x j ) dm(x j ).
(4.8) now follows from (4.6), (1.2) and the dominated convergence theorem.

It follows from (4.8) that for a.e.

t 1 , • • • , t k lim n→∞ Q z,z (n)   k j=1 f j (X (n) t j )e -αζ   = Q z,z   k j=1 f j (X t j )e -αζ   . ( 4
.12) Lemma 4.3 For any α, α j , j = 1, . . . , k and continuous compactly supported f j , j = 1, . . . , k

lim n→∞ R k + k j=1 e -α j t j ν (n)   k j=1 f j (X (n) t j )e -αζ   k j=1 dt j (4.13) = R k + k j=1 e -α j t j ν   k j=1 f j (X t j )e -αζ   k j=1 dt j .
Proof: By (4.11)

R k + k j=1 e -α j t j ν (n)   k j=1 f j (X (n) t j )e -αζ   k j=1 dt j (4.14) = π∈P k k j=1 u α+ k l=j α l (n) (x j-1 , x j )u α (n) (x k , z) k j=1 f π(j) (x j ) dm(x j ) dm(z).
If k ≥ 2, then using the resolvent equation we see that

k l=1 α l R k + k j=1 e -α j t j ν (n)   k j=1 f j (X (n) t j )e -αζ   k j=1 dt j (4.15) = π∈P k k j=2 u α+ k l=j α l (n) (x j-1 , x j )u α (n) (x k , x 1 ) k j=1 f π(j) (x j ) dm(x j ) - π∈P k k j=2 u α+ k l=j α l (n) (x j-1 , x j )u α+ k l=1 α l (n) (x k , x 1 ) k j=1 f π(j) (x j ) dm(x j ),
and (4.13) for k ≥ 2 then follows from (4.6), (1.2) and the dominated convergence theorem. Here we repeatedly use the Cauchy-Schwarz inequality. See the proof of [8, Lemma 3.3].

When k = 1,

R + e -α 1 t 1 ν (n) f 1 (X (n) t 1 )e -αζ dt 1 (4.16) = u α+α 1 (n) (z, x)u α (n) (x, z)f 1 (x) dm(x) dm(z).
Note that by (4.5) we have u α (n) (x, z) ≤ u α/2 (x, z) for n sufficiently large. We can then use the argument from the proof of Theorem 1.2 and the dominated convergence theorem to get (4.13) for k = 1.

It follows from (4.13) that

lim n→∞ R k + k j=1 e -α j t j ν (n)   k j=1 f j (X (n) t j )g(ζ)   k j=1 dt j (4.17) = R k + k j=1 e -α j t j ν   k j=1 f j (X t j )g(ζ)   k j=1 dt j
for all continuous exponentially bounded functions g. We will be particularly interested in g of the form

g(ζ) = k j=1 (1 -e -β j ζ ) k j=1 (1 -e -α j ζ ) e -αζ h s (ζ) (4.18)
where 0

≤ h s (ζ) ≤ 1 is a continuous function with h s (ζ) = 0 for ζ ≤ s.

Invariance under loop rotation

µ(F ) = ν F ζ (5.1)
is equivalent to that of the measure ν. We note that by (1.11)

ν (f (X δ )) = P δ (z, dx)f (x)u(x, z) dm(z) < ∞ (5.2)
for any δ > 0 and bounded measurable f with compact support. We next recall some ideas from [START_REF] Fitzsimmons | Markovian loop soups: permanental processes and isomorphism theorems[END_REF]. Let us define the process X to be the periodic extension of X; that is,

X t = X t-qζ , if qζ ≤ t < (q + 1)ζ, q = 0, 1, 2, . . . X t , if ζ = ∞ (5.
3)

It will be convenient to write

I α (f ) := ∞ 0 e -αt f (X t ) dt, I α (f ) := ∞ 0 e -αt f (X t ) dt. (5.4)
The key observation is that

I α (f ) = I α (f ) 1 -e -αζ , (5.5)
for all α > 0. This follows from

I α (f ) : = ∞ 0 e -αt f (X t ) dt = ∞ q=0 (q+1)ζ qζ e -αt f (X t ) dt = ∞ q=0 e -αqζ ζ 0 e -αt f (X t ) dt = I α (f ) 1 -e -αζ .
The rotation invariance of µ or ν is equivalent to the following Lemma.

Lemma 5.1 ν   k j=1 f j (X t j +r )1 {t k <ζ}   = ν   k j=1 f j (X t j )1 {t k <ζ}   (5.6)
for all 0 < t 1 < • • • < t k and r > 0 and all f j ≥ 0 continuous with compact support.

Let 0 ≤ h s (ζ) ≤ 1 be a continuous function with h s (ζ) = 0 for ζ ≤ s. To prove Lemma 5.1 we first prove the following. Lemma 5.2 For all k ≥ 1, and

0 ≤ t 1 , • • • , t k < s and all f j ≥ 0 continuous with compact support ν   k j=1 f j (X t j +r )h s (ζ)   = ν   k j=1 f j (X t j )h s (ζ)   .
(5.7)

Proof of Lemma 5.2: Using first (5.5) and then (4.17) we have that

lim n→∞ R k + k j=1 e -α j t j ν (n)   k j=1 f j X (n) t j (1 -e -β j ζ )e -αζ h s (ζ)   k j=1 dt j = lim n→∞ R k + k j=1 e -α j t j ν (n)   k j=1 f j X (n) t j (1 -e -β j ζ ) (1 -e -α j ζ ) e -αζ h s (ζ)   k j=1 dt j = R k + k j=1 e -α j t j ν   k j=1 f j (X t j ) (1 -e -β j ζ ) (1 -e -α j ζ ) e -αζ h s (ζ)   k j=1 dt j = R k + k j=1 e -α j t j ν   k j=1 f j (X t j )(1 -e -β j ζ )e -αζ h s (ζ)   k j=1 dt j . (5.8)
It follows from this that for a.e. 0

≤ t 1 , • • • , t k lim n→∞ ν (n)   k j=1 f j X (n) t j (1 -e -β j ζ )e -αζ h s (ζ)   (5.9)
= ν   k j=1 f j (X t j )(1 -e -β j ζ )e -αζ h s (ζ)   .
The same calculations show that for any r > 0,

lim n→∞ ν (n)   k j=1 f j X (n) t j +r (1 -e -β j ζ )e -αζ h s (ζ)   (5.10) = ν   k j=1 f j (X t j +r )(1 -e -β j ζ )e -αζ h s (ζ)   for a.e. 0 ≤ t 1 , • • • , t k . Since ν (n) is invariant under loop rotation, see [4, Lemma 2.4
] for the simple proof, it follows from our last two displays that for a.e. 0

≤ t 1 , • • • , t k ν   k j=1 f j (X t j +r )(1 -e -β j ζ )e -αζ h s (ζ)   (5.11) = ν   k j=1 f j (X t j )(1 -e -β j ζ )e -αζ h s (ζ)   .
We now use an argument from [START_REF] Fitzsimmons | Markovian loop soups: permanental processes and isomorphism theorems[END_REF] (see from (5.31) there until the end of the paragraph). By Fubini we can find a set T ⊆ (0, s) with full measure such that for all t 1 ∈ T we have that (5.11) holds for a.e. t 2 , . . . , t k ∈ (0, s). Using (5.2) with δ = t 1 , the boundedness and continuity of the f j and the right continuity of Xt it follows from the Dominated Convergence Theorem that (5.11) holds for all (t 1 , t 2 , . . . , t k ) ∈ T × [0, s) k-1 . Let now f 1,n be a sequence of continuous functions with compact support with the property that f 1,n ↑ 1. By the above, (5.11) with f 1 replaced by f 1,n holds for all (t 1 , t 2 , . . . , t k ) ∈ T n ×[0, s) k-1 for an appropriate T n ⊆ (0, s) with full measure . In particular T * = ∩ n T n = ∅, and we can apply the Monotone Convergence Theorem with t 1 ∈ T * to conclude that

ν   (1 -e -β 1 ζ ) k j=2 f j (X t j +r )(1 -e -β j ζ )e -αζ h s (ζ)   (5.12) = ν   (1 -e -β 1 ζ ) k j=2 f j (X t j )(1 -e -β j ζ )e -αζ h s (ζ)  
for all t 2 , . . . , t k < s. Applying once again the Monotone Convergence Theorem for β j → ∞, α → 0 we obtain

ν   k j=2 f j (X t j +r )h s (ζ)   = ν   k j=2 f j (X t j )h s (ζ)   (5.13)
for all t 2 , . . . , t k < s. Since k is arbitrary, we obtain our Lemma.

Proof of Lemma 5.1:

Fix 0 < t 1 < • • • < t k . Choose a sequence s n ↓ t k .
It is clear that we can choose h sn so that h sn (ζ) ↑ 1 {t k <ζ} . Lemma 5.1 then follows from Lemma 5.2 by the Monotone Convergence Theorem.

The loop measure and continuous additive functionals

Before proving Theorem 1.4 we will need two facts about continuous additive functionals (CAFs). The first says that to each CAF A of X is associated a measure ν A on S such that for any measurable function f

U A f (x) := E x ∞ 0 f (X t ) dA t = S u(x, y)f (y) ν A (dy), ∀x ∈ S. (6.1)
ν A is referred to as the Revuz measure of A. The second fact we need is that if a CAF has Revuz measure ν with respect to X, it has Revuz measure h • ν with respect to the h-transform of X. Following the proof of Theorem 1.4 we will discuss these facts and provide references. Proof of Theorem 1.4:

To prove (1.15) it is enough to prove the additive functional version of (2.4). We consider first our Borel right process X. These considerations will then be applied to the h-transform of X using h z = u(•, z) for fixed z ∈ S.

Let A j (j = 1, 2, . . .) be CAFs of X with Revuz measures ν j . Using the Markov property, see for example Theorems 28.7 and 22.8 of [START_REF] Sharpe | General theory of Markov processes[END_REF], and (6.1) at the last step

E x {0<t 1 <t 2 <•••<tn<∞} n j=1 dA j t j (6.2) = E x   ∞ 0   {0<t 2 <•••<tn<∞} n j=2 dA j t j   • θ t 1 dA 1 t 1   = E x   ∞ 0 E Xt 1   {0<t 2 <•••<tn<∞} n j=2 dA j t j   dA 1 t 1   = S u(x, x 1 )E x 1   {0<t 2 <•••<tn<∞} n j=2 dA j t j   ν 1 (dx 1 ),
and then by induction

E x {0<t 1 <t 2 <•••<tn<∞} n j=1 dA j t j (6.3) = S n u(x, x 1 )u(x 1 , x 2 ) • • • u(x n-1 , x n ) n j=1 ν j (dx j ).
Notice that by our assumption that u(x, x 1 ) is excessive in x for each x 1 , the expressions in (6.3) are excessive functions of x. Thus if η = (η t ) is an entrance law, then writing E η for the measure under which the one-dimensional distributions are given by the entrance law we have

E η {0<t 1 <t 2 <•••<tn<∞} n j=1 dA j t j (6.4) =↑ lim t↓0 S η t (dx)E x {0<t 1 <t 2 <•••<tn<∞} n j=1 dA j t j =↑ lim t↓0 S η t (dx) S n u(x, x 1 )u(x 1 , x 2 ) • • • u(x n-1 , x n ) n j=1 ν j (dx j ) = S n g(x 1 )u(x 1 , x 2 ) • • • u(x n-1 , x n ) n j=1 ν j (dx j ),
where g(x 1 ) :=↑ lim t↓0 η t (dx)u(x, x 1 ). Now apply the above to the h-transform of the original process X, with h z = u(•, z) for a fixed z ∈ S, as described in Section 2. This process has potential density u hz (x, y) = u(x, y)/h z (x) with respect to the measure h z (y) m(dy). Also, if a CAF has Revuz measure ν with respect to X, it has Revuz measure h z • ν with respect to the h-transform process. Thus by (6.3)

E x,z {0<t 1 <t 2 <•••<tn<∞} n j=1 dA j t j (6.5) = S n u(x, x 1 ) h z (x) u(x 1 , x 2 ) h z (x 1 ) • • • u(x n-1 , x n ) h z (x n-1 ) n j=1 h z (x j ) ν j (dx j ) = 1 h z (x) S n u(x, x 1 )u(x 1 , x 2 ) • • • u(x n-1 , x n ) h z (x n ) n j=1 ν j (dx j ) = 1 h z (x) S n u(x, x 1 )u(x 1 , x 2 ) • • • u(x n-1 , x n ) u(x n , z) n j=1 ν j (dx j ).
When we use the entrance law η z t (dx) = P t (z, dx)h z (x), the function g of the preceding paragraph is

↑ lim t↓0 η z t (dx)u hz (x, x 1 ) =↑ lim t↓0 P t (z, dx)u(x, x 1 ) = u(z, x 1 ). (6.6)
Thus, using the definition of Q z,z from Section 2,

Q z,z 0<t 1 <•••<tn<∞ n j=1 dA t t j = S n u(z, x 1 )u(x 1 , x 2 ) • • • u(x n , z) n j=1
ν j (dx j ).

(6.7) Similar considerations work for the α-potentials, and the argument given in the proof of Theorem 1.1 proves (1.15).

We now discuss the facts mentioned at the beginning of this section. Given a right-continuous strong Markov process X (more precisely, a Borel right Markov process) and an excessive measure m, there is always a dual process X (essentially uniquely determined), but in general it is a moderate Markov process: the Markov property holds only at predictable times.

In what follows f and g are non-negative Borel functions on S. By duality where the kernel Û is the potential kernel of the moderate Markov dual of X.

Under our assumptions it follows from [1, VI, Theorem 1.4] that the potential density u can be chosen so that x → u(x, y) is excessive for each y, and y → u(x, y) is co-excessive (that is, excessive with respect to the moderate Markov dual process X) for each x. (6.8) implies that Û f (y) = S u(x, y)f (x) m(dx), (

for m-a.e. y. Since both sides of (6.9) are co-excessive, they agree for all y. where ν A is the Revuz measure of the CAF A with respect to m. Feeding (6.9) into (6.10) and varying f we find that U A g(x) = S u(x, y)g(y) ν A (dy), (6.11)

first for m-a.e. x, then for all x because both sides of (6.11) are excessive. This proves (6.1). One subtlety: the laws P x of X are only determined modulo a class of sets ("m-exceptional") defined in [3], see (3.4) for the definition of the term, and then Remark (5.14); but that class is not charged by ν, so the exception causes no problem.

To establish the second fact that we needed, let τ t be the right continuous inverse of A t , and let f be a positive measurable function. Using the change of variables formula, [2, Chapter 6, (55.1)] and then Fubini

E x/h ∞ 0 f (X t ) dA t = E x/h ∞ 0 f (X τ (u) ) du (6.12) = ∞ 0 E x/h f (X τ (u) ) du.
Using [10, (62.20)] and then Fubini we have

∞ 0 E x/h f (X τ (u) ) du = 1 h(x) ∞ 0 E x f (X τ (u) )h(X τ (u) ) du (6.13) = 1 h(x) E x ∞ 0 f (X τ (u) )h(X τ (u) ) du .
Using the change of variables formula once again, the last two formulas show that This shows that if a CAF has Revuz measure ν with respect to X, then it has Revuz measure h • ν with respect to the h-transform of X. (Recall that we use u(x,y) h(x) for the potential densities of the h-transform process).

E x/h ∞ 0 f (X t ) dA t = 1 h(x) E x ∞ 0 f (X t )

  α→∞ u α (x, y) (3.4) exists and v(x, y)f (y) dm(y) = lim α→∞ ∞ 0 e -αt P t (x k , dy)f (y) dt = 0. (3.5) Hence v(x, y) = 0 for m-a.e. y. Integrating (3.3) with respect to α from 0 to ∞ and using Fubini's theorem we then obtain (1.6). (1.3) then follows by (1.7). To show that the right hand side of (1.3) is finite we repeatedly use the Cauchy-Schwarz inequality and our assumption (1.2). See the proof of [8, Lemma 3.3].

Proof of Theorem 1 . 3 :

 13 Because the lifetime ζ is rotation invariant (ζ(ρ v ω) = ζ(ω) so long as ζ(ω) < ∞), the rotation invariance of the loop measure

Sf

  (x)U g(x)m(dx) = S Û f (y)g(y)m(dy), (6.8)

  By [3, (5.13)] we have the Revuz formula S f (x)U A g(x) m(dx) = S Û f (y)g(y)ν A (dy), (6.10)

  h(X t ) dA t .

							(6.14)
	Using (6.1) we see that
	E x/h	0	∞	f (X t ) dA t =	1 h(x) S	u(x, y)f (y)h(y) ν A (dy).	(6.15)

* Research of J. Rosen was partially supported by grants from the National Science Foundation and PSC CUNY.


Frédéric Giraud

Member IEEE Betty Semail

Member IEEE Alain Bouscayrol

Modelling and Control of an Ultrasonic Motor: Application to a Mechanical Claw Drive

Properties of the travelling Wave Ultrasonic Motor (TWUM ) -among which we find a high-torque to mass ratio, locking without supply, ultrasonic operation -are very attractive for embedded application, such as human prothesis for example. But unfortunately, its control for high precision or high dynamic motion is not easy to achieve because of many non linearities induced by its particular energy conversion process. This is why this paper deals with the modelling of a TWUM, based on a causal approach. According to the inversion principle, control strategies are deduced as well as a novel torque estimation structure. A mechanical claw control using this torque estimator is then achieved and experimental results check the accuracy of the proposed control scheme.

I. INTRODUCTION

Travelling wave Ultrasonic Motors (TWUM ) exploits a piezoelectrically generated flexural wave which propagates at the surface of a stator. This wave is able to propel by contact a rotor strongly pressed on the stator. The friction caused by the contact mechanism provides us with numerous advantages, including: breaking without supply, a high torque to mass ratio and high torque/low speed characteristics. Thus, while a speed reducer is often needed with an electromagnetic motor, it becomes useless in applications using TWUM, therefore leading to lightweight and compact applications. The last property is particularly important for embedded applications such as medical prothesis ones. Moreover, the working principle doesn't rely on magnetic field interactions which may be awkward for human uses.

Unfortunately, because of the mechanical non linearities, TWUM modelling is not easy to achieve, especially in order to control its behavior, such as torque, speed or position. This is due to the different energy conversions and more specifically to the complications caused by the contact mechanism between the stator and the rotor.

Several references have proposed modelling and control of the TWUM [START_REF] Maas | Model based control for ultrasonic motors[END_REF][2] [START_REF] Ghouty | Hybrid modeling of a traveling wave piezo electric motor[END_REF] The one used in this paper describes the real internal variables behavior. It will be first reminded and explained according to a causal formalism. Then, the control laws will be described following the inversion principles. For the particular application of an active claw control, a torque estimator will be explained. Experimental results will check the validity of the control laws as well as the modelling and estimator accuracy.

II. MODELLING AND CONTROL OF THE TWUM

A. Modeling overview

This section deals with a modelling we proposed for taking into account transitory operations [4][5]. This modelling relies on the assumption of a punctual contact between stator and rotor. Punctual means that the stator is in contact with the ideal rotor at one point only along a wavelength. Moreover, no sliding is considered between the ideal rotor and the stator, so both have the same instantaneous speed. On the figure 1, we have represented the travelling wave which propels the ideal rotor by only one point along a wavelength. In the fame fixed to the contact point ( -→ u r , -→ u θ , -→ u z ) the travelling wave propagates toward the right (along -→ u θ ) although the ideal rotor moves toward the left (--→ u θ ). We also define the zaxis about which the motor and the contact point turn. The location of the contact point is defined by the rotation angle θ c measured from a fixed frame (O, x, y).

1) Causal modelling in (α, β) frame:

The TWUM is composed of a ring shaped stator on which small piezo-electric elements are bonded. These elements have alternate polling directions, and are divided into two groups, named phase α and phase β. These two phases are supplied by two independent voltages, v α and v β . Because of the piezo-electric effect, the first energy conversion stage converts these voltages into forces inside the stator; these forces bend the stator. And if the excitation frequency of the two supply voltages is sufficiently close to the resonant frequency of the stator, two stationary waves appear at its surface: the bending waves, whose height are named w α and w β . This second stage of the energy conversion can be described by the equation ( 1)

m ẅα + d s ẇα + cw α = N v α -F rα m ẅβ + d s ẇβ + cw β = N v β -F rβ (1)
with

• m the vibrating mass of the stator,

• c the vibrating stiffness of the stator,

• d s the damping coefficient,

• N the force factor. In this equation, F rα and F rβ have been added to take into account the effect of external forces on the wave propagation. Their expression which depends on F N (the force pressing the rotor against the stator along the axis of the motor) and F ′ T (the tangential force deriving from the load torque T ) is quite difficult to determine, and strongly depends on the contact condition at the stator-rotor interface [START_REF] Hagood | Modeling of a piezoelectric rotary ultrasonic motor[END_REF] [START_REF] Maas | Simulation model for ultrasonic motors powered by resonant converters[END_REF]. However, in the ideal case of a punctual contact condition, these forces are simply expressed using the location of the contact point θ c , and the well known rotational matrix, usually used in the study of electromagnetic machines. This calculation has been developed in [START_REF] Giraud | Causal modeling and identification of a travelling wave ultrasonic motor[END_REF], and derives from power balance of the system:

F rα F rβ = R(kθ c ) F N F ′ T R(θ) = cos(θ) -sin(θ) sin(θ) cos(θ) (2) with
F ′ T = k h b 2 T (3)

And

• h the thickness of the stator • b the radius of the rotor • k = 2πb λ with λ the wavelength of the excited mode. The contact points are located at the travelling wave's crest. It can be shown that they are situated at an angle given by w α and w β using the expression shown in ( 4)

kθ c = atan w β w α (4)
Moreover, as it can be seen in figure [START_REF] Maas | Model based control for ultrasonic motors[END_REF], the contact point speed which propels the rotor, has two components:

• V N id along the axis of the motor, lifts up the rotor to enable the propagation of the travelling wave. During steady state, this speed is null.

• V ′
T id in the tangential direction , from which the rotational speed of the rotor ω id derives. A kinetic study of the vibrating stator shows then that the expression of those speeds can be written using ẇαβ and R(kθ c ) once again:

V N id V ′ T id = R(-kθ c ) ẇβ ẇα (5)
The rotational speed of the ideal rotor ω id and V ′ T id is given by equation (6) [START_REF] Sashida | An introduction to Ultrasonic Motors[END_REF]:

ω id = k h b 2 V ′ T id (6)
However, the speed ω id does not actually correspond to the rotor speed, because the contact conditions observed at the stator-rotor interface are very different from those of the ideal contact conditions: in fact, because of friction, the actual rotational speed ω of the rotor decreases with the load torque T . Several models ( [START_REF] Maas | Simulation model for ultrasonic motors powered by resonant converters[END_REF][8] [START_REF] Maeno | Finite element analysis of the rotor/stator contact on e ring-type ultrasonic motor[END_REF]) exist to describe this behavior, but they are too complex to be taken into account in a control scheme. In order to achieve a straightforward modelling which may be useful for actuator control, a linear modelling for the global behavior of the motor can be proposed based on experimental results. This is the third energy conversion, given by equations [START_REF] Sashida | An introduction to Ultrasonic Motors[END_REF].

T = f 0 (ω id -ω) F N = K N (V N id -V N ) dt (7)
Of course, this process is very non linear, and the value of the parameters f 0 and K N should vary as a function of the operating point. However, we find here the key variables governing the evolution of the motor torque: the amplitude of the travelling wave, because ω id derives using it, and the rotational speed of the rotor itself. This roughly outlines the torque estimator described in [START_REF] Maas | Optimized drive control for inverterfed ultrasonic motors[END_REF] from the speed and the stator's deformation measurements.

In addition, the dynamic of the rotor, along the two degrees of freedom, leads to lay down two other equations:

m r dV N id dt = F N -F τ J r dω dt = T -T r (8)
Where:

• m r and J r are respectively the mass and rotational inertia of the rotor, • T r is the external load torque applied on the rotor, • F τ is the normal force pressing the rotor against the stator. Finally, all these equations yield a modeling which can be represented by way of an Energetic Macroscopic Representation (EMR) as shown in figure 2. The EMR is useful to highlight the energy conversion process and the electromechanical conversion. Following this formalism, the mechanical transformation processor are depicted with triangles, otherwise, the energy accumulators are described using barred rectangles.

2) Causal modelling in a rotating reference frame: the two rotational matrices, on the stationary waves (eq. 5) and on the forces (eq. 2), lead to setting down two new variables V d and V q . As it is for the electromagnetic machines where current or voltages are expressed in a rotating reference frame fixed to the rotating electromagnetic field, V d and V q represent the supply voltages in a rotating reference frame fixed to the rotating travelling wave. So, this leads us to set down:

v α v β = R(kθ c ) V d V q (9)
Using V d and V q , equations 1, 2, 5, are revised to produce two new equations respectively on d axis and q axis respectively, where the couplings are removed:

m V ′ N id + dsV ′ N id + (c -mk( θc ) 2 ) V ′ N id dt = N V d -F ′ N (10)
2m V ′ T id + dsV ′ T id = N V q -F ′ T ( 11
)
In fact, in equation ( 11), the variables of the tangential axis appear only on the q axis, though equation [START_REF] Maas | Optimized drive control for inverterfed ultrasonic motors[END_REF] shows that the normal variables only appear on the d axis. This equation is of second order type, and completely describes the resonant behavior of the motor.

Upper part of figure 3 shows the EMR of the actuator in the rotating reference frame, reduced to the q axis as it deals with the torque.

B. Control strategies 1) Maximal control structure:

From EMR, it is easy to deduce the maximal structure of control, according to inversion principle. Transformation relations are directly inverted, but energetic storage elements are inverted using regulators. The lower part of figure 3 shows the speed control scheme as it may be deduced following these rules.

2) Non linearities considerations: The torque-speed relation of 7 has been approximated from the curves depicted at figure 4. The direct inversion of the relation is quite easy to achieve, but obviously may lead to an inaccurate controlled torque. To cope with this problem, two approaches are proposed: a torque regulation in a closed loop, or a global Fig. 3. EMR of the TWUM in the rotating reference frame -axis q onlyand the derived control scheme. control including a behavior modelling regulator as detailed in [START_REF] Giraud | Precise position control of a travelling wave ultrasonic motor[END_REF]. Nevertheless, this last solution is available in case of speed or position control. For torque control, a closed torque loop is necessary, as well as the knowledge of the actual output torque. In the specific application of a claw dedicated to a robot arm edge or a hand prothesis, the torque variations have to be monitored so as to indicate an object tightening for instance. The torque estimation may be achieved according to the modelling as explained in II-A.

3) The torque estimator: In order to avoid a speed sensor, the torque estimator will rely on equation 11. However, inverting this equation is not straightforward, because of the derivative term 2m V ′ T id , which should be approximated. A closed loop estimator could then be proposed, whose output would converge toward the solution. But considering a high performing estimator is required especially in the case of steady state operations -or slow variations of V ′ T id , the problematic term 2m V ′ T id is deleted in equation ( 11):

dsV ′ T id = N V q -F ′ T ( 12
)
Finally, using equations 3 in 12 leads to

T = b 2 kh [N V q -d s V ′ T id ] (13)
Achieving the torque estimator.

III. EXPERIMENTAL IMPLANTATION OF THE TORQUE

ESTIMATOR [START_REF] Giraud | A torque estimator for a travelling wave ultrasonic motor -application to an active claw[END_REF] The TWUM used during the test is a Shinsei USR30 [START_REF] Shinsei | Internet web site[END_REF]. This motor has a 30mm diameter and can provide a rated torque of 0.1N m for a rated speed of 200rev/min. On the stator of this motor, an extra electrode is glued to the ring of piezo material. Due to the direct piezo-electrical conversion, the voltage measured on this sensor, named v EA , is directly proportional to the deformation of the stator at this point. Unfortunately, this deformation is a combination of w α and w β : thus it is difficult to build θ c from the equation (4) because w α and w β are not measured independently. In that case, measuring V q from V α , V β and ( 9) is impossible. This leads us to set down the assumption that the amplitude of the two stationary waves are sinusoidal functions of time in quadrature, so:

w α (t) = W cos(2πf t) w β (t) = W sin(2πf t) (14)
With f the frequency of the supply voltages and W the travelling wave's amplitude.

To fulfill this assumption, the construction of the actuator must be sufficiently careful to ensure a good positioning of the stator's electrodes on the piezo ring, thus avoiding dissymmetry between the two phases. However, this condition is often achieved in the Shinsei USR30 when it is supplied by two sinusoidal voltages in quadrature with the same amplitude:

v α (t) = V cos(2πf t + Ψ) v β (t) = V sin(2πf t + Ψ) (15)
with V the amplitude of supply voltages, and Ψ, the phase difference between v α (t) and w α (t). On the experimental test bench, V is kept constant, though f varies.

Then, equations ( 4) and ( 14) yield to:

θ c = atan W sin(2πf t) W cos(2πf t) = 2πf t (16)
To complete the identification of T from equation [START_REF] Shinsei | Internet web site[END_REF], it is necessary to know V ′ T id , but however this is an internal variable. It may be calculated thanks to equation [START_REF] Giraud | Causal modeling and identification of a travelling wave ultrasonic motor[END_REF], and from the measurement of ẇα and ẇβ . Once again, to reach these values, we use the previous assumption which leads to equation [START_REF] Maas | Averaged model of inverter-fed ultrasonic motors[END_REF]. Then, we find the relationship between the wave amplitude and the tangential speed of the rotor, using equations ( 16) and the rotational matrix (5):

V ′ T id = -ẇα sin(2πf t) + ẇβ cos(2πf t) = (2πf )W (17)
Therefore, V ′ T id is proportional to the travelling wave's amplitude. Now, W can be deduced by measuring the voltage amplitude of v EA ; the different state variables of the torque estimator can then be measured. But parameters appearing in the equation ( 13) have also to be identified.

A. Measurement of the stator parameters 1) Identification of d s :

For an unloaded motor which has reached its steady state, the travelling wave's amplitude W derives from the equation 13 where T is set to 0. In these conditions:

d s 2πf W = N V q (18)
Identifying d s can then be achieved by keeping V constant; plotting the evolution of W against V q should result in a straight line, whose slope gives the ratio N ds2πf . However, this method is quite sensitive to non linearities which induce variations of d s with W [START_REF] Maas | Averaged model of inverter-fed ultrasonic motors[END_REF]. This is why, instead of V , W is kept constant. So, for each value of V q , V is tuned in order to have W equal to a value we fixed before.

It appears that d s depends more on the operating wave amplitude W and less on the rotational direction. In order to take into account this non linearity in the torque estimator, we are led to representing in figure 5 the evolution of d s in W , or rather in 1 W because this results in a linear shape. Equation 19 then fits this curve:

d s = -d s0 + ρ W (19)
The values of d s0 and ρ are obtained from figure 5:

d s0 = 19.2N.s.m -1 ρ = 47.8.10 -6 N.s
So, taking into account the non-linearity of d s , the torque estimator can be achieved by merging the equations 13, 17 and 19:

T = b 2 kh (d s0 2πf W + N V q -2πf ρ) (20) 2) Identification of k h b 2 :
The identification of this parameter derives once again from equation [START_REF] Shinsei | Internet web site[END_REF]. But now, T will vary and V q is adjusted to keep W constant. So, the evolution of V q in T should be a straight line whose slope is k h N b 2 . We represents in figure ?? the result curves for several values of W .

The results show that the experimental trials and the modelling are consistent, since the curves are nearly straight. However, the slope of these lines depends on the level of the stator's deformation W and also on the direction of rotation. These variations reveal non-linearities in the modelling; in fact, they express the limits of this modelling, and more specifically the limit imposed by the assumption that there is only one contact point. This is why, an average value of k h b 2 is then calculated from the trials, and the slope of the curves drawn on figure 6 yielding k h b 2 = 67.2m -1 .

B. Experimental runs

Different trials have been carried out to illustrate the performance of this torque estimator. The first measurements shown in figure 7(a) depict the evolution of the estimated load torque and the measured load torque. During this trial, the stator wave's amplitude was kept constant (W ≃ 0.5µm) while the load torque T was time varying.

The second measurement of figure 7(b) shows the evolution of the measured load torque compared to its estimation when W was time varying.

Considering that the torque estimator is quite straightforward, the estimated load torque is found to be consistent with the measurements of the experimental trials of figures 7(b) and 7(a). Estimation errors can however be measured on figure 7(b), revealing that the estimation depends on W . More obviously, errors are amplified for braking operations while they are acceptable for other conditions. This is mainly due to the averaged value of k h b 2 we had to use in equation 20 and which is a tradeoff between all the values identified for all the operating points of the chapter III-A. Performances of the estimator could then be improved by taking into account the variation of k h b 2 in W , but this would be achieved at the expense of complexity. Despite this problem, applications can take advantage of this torque estimator, and this is dealt with in the next section.

IV. APPLICATION: FORCE CONTROL OF TWUM ACTUATED CLAW

A. Working principle of a claw

Because of torque vs speed characteristics, Travelling Wave Ultrasonic Motors are very suitable for applications which need a low speed but a high torque. In this field of applications, the electrically actuated claws are devices which are widely used, for example at the end of a robot arm to allow object prehension.

Different controls of the claw can be achieved. The first control type is a position control of the actuator of the claw between an opened position till the closed position. This is well suited for objects which all have the same dimensions. But if variable objects of random sizes have to be picked up, the larger objects are held extremely tight, while the smallest ones are not even caught (figure 8). This is why, instead of having a claw position control, a force limitation is preferred. In this operating mode, the claw closes up at a constant speed until the force at the claw 's tip reaches a limited value.

This last control method is quite easy to achieve with electromagnetic motors, because the torque is proportional to the supply current which is often controlled, and can thus be limited. But those motors have high speed characteristics, so they need a speed reducer in order to adapt the output speed of the shaft to the desired speed of the claw when closing. This increases the system's bulk and its cost.

In this application, a TWUM is interesting, because it can be directly mounted on the claw since its low-speed characteristic is well adapted to the system. The torque estimator developed in this article is implemented to ensure the limitation of the tightening force.

The control scheme is depicted in figure 9. The torque estimator is depicted and its output is used as an input for a strategy block. This strategy block is inserted to determine the reference wave amplitude: below the limit torque value, the wave amplitude reference is kept constant; when the torque limit is reached, the wave amplitude is set to zero, so as to lock the claw. The next section presents the experimental runs of the claw actuated by a TWUM, with its torque limitation.

B. Experimental runs

The claw fixed on the TWUM's shaft is depicted in figure 10. Several objects with different sizes can be positioned, and for each run, the internal variable of the motor can be recorded. Thus, figure 11 depicts the evolution of V q , the motor rotational angle, the measured and estimated torque and finally the wave amplitude W -which is also controlled. These curves are obtained with T lim = 0.05N m. The cycle begins at t = 0 with the starting signal; a wave amplitude reference of about 1µm is ordered, and the claw closes. The estimated torque as well as the measured one are initially null, but at t = 50ms they both start to increase because the claw is in contact with the object to be grasped. Consequently, V q increases in accordance equation (20), since W is kept constant.

At t = 85ms, the estimated torque reaches the limitation: a reference value of W of 0µm is ordered which is obtained in almost 10ms. During this short time, the torque actually applied to the claw can be a little larger due to the time the wave amplitude needs to become null. Then the torque estimator is switched off because the very low wave amplitude leads to erroneous estimated torque values. So, the claw stops its grasping action when it meets an object. This contact is detected through an increasing in the motor's torque, which is estimated on line. Of course, the force with which the object is taken depends on the value of T lim .

V. CONCLUSION

In this article, we have presented a method to achieve a TWUM control for a mechanical claw function application. According to he causal modelling of the actuator, several control structures have been described. More specifically, we have developed a torque estimator in order to achieve the torque monitoring for a tightening detection. This estimator relies on the causal modelling of the motor which helps to deduce the estimator's equation. Compared to the torque estimators which rely on the torque/speed characteristics of the motor, we emphasize that the estimator presented here doesn't need the measurement of the rotational speed of the motor; this allows to remove the speed sensor, and the measurements have confirmed the global structure of the modelling, while the estimator is successfully tested on the claw application in the final section.

Future work should increase the precision of this torque estimator, more specifically at low rotational speeds, to allow sensorless torque regulators. This implies the control of voltages phases leading to the loss of quadrature. A new modelling should be introduced to take into account this voltage supply.

Fig. 1 .

 1 Fig. 1. View of the stator; under punctual contact assumption, the rotor is propelled by one point only along a wavelength.

Fig. 2 .

 2 Fig. 2. EMR of the TWUM.

Fig. 4 .

 4 Fig. 4. Torque-speed characteristics of the TWUM for given wave amplitude.

Fig. 5 .

 5 Fig. 5. Evolution of ds as a function of 1 W for both rotational direction and fitting.

Fig. 6 .

 6 Fig. 6. Vq as a function of T and comparison with a straight line.

  Fig. 7. Comparison between the actual load torque and the estimator's output. Experimental runs.

Fig. 8 .

 8 Fig. 8. Effect of a controlled position claw when the object is sufficiently big, and in the case where it is too small.

Fig. 9 .

 9 Fig. 9. Control structure for an active claw application.

Fig. 10 .

 10 Fig. 10. The experimental Setup.

Fig. 11 .

 11 Fig. 11. Internal variables of the motor for a closing sequence of a TWUM actuated claw , with limitation of the torque.

  VI. ACKNOWLEDGMENT This work has been carried out within the framework of the INRIA Alcove project and is supported by the IRCICA (Institut de Recherche sur les Composants logiciels et matériels pour l'Information et la Communication Avancée).


W Lhomme

P Delarue

F Giraud

B Lemaire-Semail

A Bouscayrol

-Simulation of a PV system using EMR, synopsis EPE-PEMC 2012 --1

The Energetic Macroscopic Representation (EMR) has been developed in 2000 to develop control of electric systems. Since 2002 this graphical tool has been introduced to teach control of electrical systems in France, then Canada and Switzerland. A Photovoltaic Conversion System is simulated using EMR in order to present an actual and valuable example for students within renewable energy applications..

Introduction

Nowadays, control of electrical systems is of main interest in electrical engineers teaching [START_REF] Blaabjerg | An industry/university collaboration in power electronics and drives[END_REF]- [START_REF] Gambôa | Modern Technologies for Experimental Education in Industrial Electronics and Electric Drives[END_REF]. But it requires skills in various scientific fields such as power electronics, electrical machines, automatic control, mechanics and control electronics.

Several graphical description has been introduced in order to represent in a uniform way multi-physical systems, such as Bond-Graph (BG) [START_REF] Paynter | Analysis and design of engineering systems[END_REF], Power Oriented Graph (POG), Causal Ordering Graph (COG) [START_REF] Hautier | Causal Ordering Graph[END_REF] [START_REF] Hautier | The causal ordering graph -A tool for modelling and control law synthesis[END_REF], or Energetic Macroscopic Representation (EMR) [START_REF] Bouscayrol | Multimachine Multiconverter System: application for electromechanical drives[END_REF]. BG and POG are structural descriptions more dedicated to the system design. COG and EMR are functional descriptions more focused on the control design, because of the use of the exclusive integral causality [START_REF] Iwasaki | Causality and model abstraction[END_REF] [START_REF] Rubin | The development of vehicular powertrain system modeling methodologies: philosophy and implementation[END_REF]. In that case, the control scheme of the system can be directly deduced from an inversion of the model description [START_REF] Barrre | Inversionbased control of electromechanical systems using causal graphical descriptions[END_REF]. These two descriptions have been successfully used for control of various applications [START_REF] Barrre | Inversionbased control of electromechanical systems using causal graphical descriptions[END_REF]- [START_REF] Bouscayrol | Power strategies for Maximum Control Structure of a wind energy conversion system with a synchronous machine[END_REF].

Because of the system-oriented description for modelling and control, EMR is a valuable methodology to teach control of electrical systems. This formalism is now taught in some French universities [14][15], at University of Québec Trois Rivières (Canada), at Ecole Polytechnique Fédérale de Lausanne (Switzerland), and some recent experiences have been realized at Tsinghua University (China), Aalto University (Finland) and Polytechnical University of Catalunya (Barcelona, Spain) [START_REF]EMR[END_REF]. Indeed such a graphical description enables a unified way for causal description of the components of electromechanical. Moreover, summer schools are regularly organized worldwide (Lille, France, in 2006, Hrabin, China, 2008, Trois-Rivières, Canada, 2009, Lausanne, Switzerland, 2011, Mardi, Spain, 2012) [START_REF]EMR[END_REF].

But, because EMR is a new systemic approach, realistic and simple examples have to illustrate the different concepts. Moreover, practice is required in the manipulation of the different pictograms, their connections and the deduced control, in particular during summer schools when the students have only 3 days of lectures on this topic. Simulation of simple systems have thus been developed in order to illustrated the modelling and the control methods using EMR in very short simulation sessions (typically two sessions of four hours each). First the simulation of an electric vehicle using a DC machine has been developed using EMR library and Matlab-Simulink© [START_REF] Bouscayrol | Teaching drive control using Energetic Macroscopic Representation -initiation level[END_REF]. But other systems have to be simulated to demonstrate the flexibility of the method.

In this paper, the simulation of a Photovoltaic (PV) system is studied in order to supply a hydraulic pump. Energy is extracted from the system using a MMPT strategy and is stored in a battery tank. The pump is moved by a DC machine. The EMR of the system is described, and the control scheme is derived from the inversion of this EMR. Finally the entire system is simulated in Matlab Simulink using an EMR library.

I. EMR fundamentals

Interaction principle -The system is decomposed in basic subsystems in interactions (Table 1): energy sources (green ovals), accumulation elements (orange rectangles), conversion element without energy accumulation (various pictograms) and coupling elements for energy distribution (orange overlapped pictograms). All the elements are connected according to the action and reaction principle using exchange variable (arrows). The product of action and reaction variables between two elements leads to the instantaneous power exchanged.

Causality principle -As in COG [START_REF] Hautier | The causal ordering graph -A tool for modelling and control law synthesis[END_REF], only the integral causality is considered in EMR. This property leads to define accumulation elements by a time-dependant relationship between its variables: its outputs are integral functions of its inputs. Other elements are described using relationships without time dependence. In order to respect the integral causality specific association rules are defined.

Inversion principle -

The inversion based control theory has been initiated by COG [START_REF] Hautier | The causal ordering graph -A tool for modelling and control law synthesis[END_REF]. The control structure of a system is obtained from an inversion of the system modelling: the control has to define the inputs to apply to the system given the desired output (Fig. 1). In this method, relationships which are not time-dependant (i.e. the elements without energy accumulation) are directly inverted (with neither control nor measurement). Because the derivative causality is forbidden, a direct inversion of time-dependant relationships is not possible. An indirect inversion is thus made using a controller and measurements. These inversion rules have been extended to EMR (blue pictograms, see Table 1): conversion elements are directly inverted, accumulation elements are inverted using controller and coupling devices are inverted using distribution inputs. It is another way to locate controllers and measurements.

II. Modelling and EMR of the studied PV system

Studied PV system -The studied system aims to supply a pump in desert area. It is composed of a PV panel, a LC filter, a boost chopper, a battery, a buck chopper and a DC machine to activate the pump (Fig. 2).

Modelling of PV system -

The different relationships of each component will be detailed in the final paper. For example, the capacitor voltage of the filter u C is expressed in function of the PV current i PV and the inductor current i L :

C L PV C u R i i u dt d C 1    ( 1 )
where C and R are the capacitance and resistance of the capacitor. As a second example, the boost chopper is described by an average relationship connecting the chopper voltage u ch1 to battery voltage u bat , and the chopper current i ch1 to the inductor current i L :

     L ch ch bat ch ch i m i u m u 1 1 1 1 with 1 1 1 ch ch m    ( 2 )
where  ch1 is the duty cycle of the chopper. EMR of PV system -Each component is translated into EMR elements and their inputs and outputs are defined according to the causality principle. Moreover, connection between elements must respect the interaction principle. For example, the capacitor ( 1) is an accumulation element, where the voltage u C is the output in both directions and the currents i PV and i L are the inputs (see Fig. 3). A crossed orange pictogram is thus used to represent this energy storage. All other EMR elements of the WECS will be described in the final paper. The EMR of the whole system is given in the upper part of Fig. 3.

III. Control of the studied PV system

Inversion-based control of PV system -The control scheme of the PV system is obtained from an inversion of its EMR (see lower part of Fig. 3). The first objective is to extract the maximum of the solar power by acting on chopper 1. The second objective is to impose the right machine velocity to regulate the water flow, by acting on chopper 2. Two control parts are thus organized.

Each relationship is inverted according to inversion rules. For example relationship ( 2) is directly inverted in order to obtain the modulation ratio reference m ch1-ref from the voltage reference u ch1-ref :

ref ch meas bat ref ch u u m     1 1 1 ( 3 )
A blue parallelogram is used as pictogram for (3). Because relationship ( 1) is integral, it can not be directly inverted. A closed-loop control is used for an indirect inversion in order to obtain the current reference i L-ref from the voltage reference u C-ref and the compensation of the PV current:

meas PV meas C ref C ref L i ) u u )( t ( C i        ( 4 )
with C(t) the controller, which could be P, PI or other kind. A crossed-blue parallelogram is used as pictogram for (3). The inversion of others elements will be described in the final paper.

MPPT strategy -A "Pertub & Observ" Maximum Power Point Tracking strategy [START_REF] Femia | Optimization of perturb and observe maximum power point tracking method[END_REF] is used to define the reference voltage to impose on the PV panel in order to get the maximal power in function of the irradiance and the temperature (Fig. 4). This strategy will be detailed in the final paper.

chopper 2 pump u C i PV u ch1-ref MPPT DC machine i L-ref PV i L u C i L u ch1 i ch1 i ch2 u ch2 u bat u bat i dcm i dcm e dcm T dcm T p   battery chopper 1 filter m ch1 u C-ref m ch2 u ch2-ref i dcm-ref T dcm-ref  ref

IV. Simulation of the studied PV System

Simulation requirements -Matlab-Simulink TM is chosen as simulation software. The students have three programs at their disposal: an initialisation program containing the PV system parameters (PV look-up table, filters, battery, and DC machine), a Simulink program with the PV source and the pump, and a Simulink library containing the EMR basic elements.

The objective of the simulation is to study the PV system using a Maximum Power Point Tracking strategy to extract power from the sun, to store this energy when the pump is off, and to use this energy to impose a constant water flow when the pump is on. The students have to build the simulation model of the PV system. They have also to determine all the control structure of both parts, to tune all the controllers and to implement this control in the simulation software. Power electronics is simulated using average values. The MPPT strategy block is provided.

Simulation sessions -Two simulation sessions of four hours each are organized. The first session is focused on the development of the simulation of the system, using the EMR library. The second session is focused on the control part. In order to organise their simulation, the students are invited to describe EMR of the system and the deduced control. Different variables will be plotted in the final paper.

The whole simulink program is achieved by most of the students (Fig. 5). They can study the power extracted from the PV panel under different irradiance (Fig. 6), the battery charge and the rotation speed of the pump. More details will be given in the final paper on that extension.

Conclusion

The use of EMR is an efficient methodology to develop student skills on drive control. A unified and physical modelling is developed, and inversion rules are used to find the control scheme. By using this intermediary step, they can easily connect other scientific fields such as power electronics, electrical machines, mechanics and automatic control.

The simulation of a PV System is proposed to acquire the concepts of EMR and of the drive control. First the description of the system is developed in a causal and systemic way. Secondly, the control scheme is deduced and the students can better understand which variable to measure, how many controllers to use, and where they has to be locate. This example is really motivating for the students because it is a renewable energy application.

Fig. 3 : 2 MPPTFig. 4 :

 324 Fig. 3: EMR and deduced control of the studied WECS

Fig. 5 : 6 :

 56 Fig. 5: Matlab-Simulik TM model of the studied PV system and its control Fig. 6: Extracted power for a variable irradiance

Table 1 :

 1 Elements of EMR and of control

		Source element		Accumulation		Indirect
		(source of			element			inversion
		energy)		(energy storage)		(closed-loop
									control)
		Mono-physical		Mono-physical		Direct inversion
		conversion			coupling			(open-loop
		element (no		element (energy		control)
		energy storage)		distribution)
		Multi-physical		Multi-physical		Coupling
		conversion			coupling			inversion
		element (no		element (energy		(energy criteria)
		energy storage)		distribution)
				PV panels	filter	battery	DC machine	pump
	u(t)	System	y(t)		i L	L	i ch1	i ch2	i dcm
	u tun (t)	C ontrol ? ?	y ref (t)	i pv					T dcm c	J	T p
						u C	u ch1	u bat	u ch2		
						chopper 1	chopper 2
	Fig. 1: Inversion-based control principle				Fig. 2: Studied PV system


L Audibert
email: lorenzo.audibert@edf.fr

A Girard

H Houssem

Imaging of anomalous components in unknown background

Keywords: Inverse scattering, Sampling methods, differential measurement

We introduce a qualitative method capable of imaging defects in an unknown complex environment using differential measurements. The main difficulty is that the background medium is unknown and too complex to obtain a reliable estimate of the associated Green function. To overcome this difficulty our approach exploits two measurements of the farfield operators, one without defects and one with defects. The analysis of our method relies on the recently introduced Generalized Linear Sampling Methods (GLSM) and the link to the solutions of the interior transmission problems. We give numerical examples related to non destructive testing in concrete-like materials, illustrating the performance of our method.

Introduction

We are interested in the imaging of defects inside unknown heterogeneous medium from multistatic measurements of waves at a fixed frenquency. The main difficulty of our setting is that the background medium is unknown and complex, which means that one cannot obtain a good estimate for the background Green function. This discards the possibility of using classical (qualitative) imaging methods [START_REF] Cakoni | Qualitative methods in inverse scattering theory[END_REF]. We assume that two sets of measurements, one for the defect free and one for the defect containing medium, are at our disposal and we design a numerical inversion algorithm that exploits those supplementary data to visualize the defect location. The analysis of our method relies on the GLSM and the link it has with the solution to the interior transmission problems (see [START_REF] Audibert | A generalized formulation of the linear sampling method with exact characterization of targets in terms of farfield measurements[END_REF] and [START_REF] Audibert | Identifying defects in an unknown background using differential measurements[END_REF]). Ones ends up with an indicator function that combines the indicator functions from the GLSM and a filtered difference term computed without additional significant numerical cost. We shall briefly introduce the scattering problem by a heterogeneous medium and the needed notations in Section 2. In Section 3 we review the main theoretical results of [START_REF] Audibert | Identifying defects in an unknown background using differential measurements[END_REF]. We give in Section 4 some numerical experiments on concrete like materials, illustrating the performance of our method.

Scattering by an inhomogeneous medium

We restrict ourselves to the case of scalar time harmonic waves and we focus on full aperture farfield measurements associated to incident plane waves. For a wave number k > 0, the total field solves the Helmholtz equation:

∆u + k 2 nu = 0 in R d
for d = 2 or 3 and n the refractive index, where (n) ≥ 0. We denote by D the support of n -1 and assume that D is a bounded domain with Lipschitz boundary and connected complement. We are interested by the case where u is generated by an incident plane wave, u i (θ, x) := e ikx•θ for x ∈ R d and θ ∈ S d-1 . We also introduce the scattered field u s defined by:

     u s (θ, •) := u -u i (θ, •) in R d , lim r→∞ |x|=r ∂u s ∂r -iku s 2 ds = 0. (1)
We introduce the farfield u ∞ (θ, x) defined through the following expansion:

u s (θ, x) = e ik|x| |x| (d-1)/2 (u ∞ (θ, x) + O(1/|x|)) for |x| → ∞ and for all (θ, x) ∈ S d-1 × S d-1 .
Leading to the farfield operator:

F g(x) := S d-1 u ∞ (θ, x)g(θ)ds(θ).
It is well known that the farfield operator admits two factorisations

F = GH = H * T H. The compact operator H : L 2 (S d-1 ) → L 2 (D)
is defined by :

Hg := S d-1 e ikx•θ g(θ)ds(θ), g ∈ L 2 (S d-1 ), x ∈ D, (2)
and is dense in

{v ∈ L 2 (D) s.t. ∆v + k 2 v = 0 in D}. Its adjoint H * : L 2 (D) → L 2 (S d-1
) is defined by :

H * ϕ(x) := D e -iky.x ϕ(y)dy, ϕ ∈ L 2 (D), x ∈ S d-1 .
We define the compact operator

G : R(H) ⊂ L 2 (D) → L 2 (S d-1 ) defined by: Gv := w ∞ where w ∞ is the farfield of w ∈ H 1 loc (R d ) associated to the incident wave v :      ∆w + nk 2 w = k 2 (1 -n)v in R d , lim r→∞ |x|=r ∂w ∂r -ikw 2 ds = 0. (3)
Finally we define T : L 2 (D) → L 2 (D) by:

T v := -k 2 (1 -n)(v + w). (4)
In the following we will use the operator

F # = | (F )| + | (F )|
, which can be factorised as F # = H * T # H, where the operator T # is a real selfadjoint operator [START_REF] Kirsch | The factorization method for inverse problems[END_REF] that satisfies (under hypothesis 1) :

|(T # h, h)| = (T # ) 1 2 h 2 ≥ µ h 2 ∀ h ∈ R(H), (5)
where µ > 0 is a constant independent of h.

Theoretical Results

The GLSM relies on the solvability of the so-called interior transmission problem defined by

(u, v) ∈ L 2 (D) × L 2 (D) such that u -v ∈ H 2 (D) and ITP(D, f, g, n, n ) =        ∆u + k 2 nu = 0 in D, ∆v + k 2 n v = 0 in D, (u -v) = f on ∂D, ∂ ∂ν (u -v) = g on ∂D, (6)
We denote by σ(D, n, n ) the set of wave number k ∈ R for which the ITP(D, f, g, n, n ) is not well posed for all f ∈ H 1 2 (∂D) and g ∈ H -1 2 (∂D).

Hypothesis 2. We assumed that k 2 ∈ R + and n are such that for all f ∈ H

1 2 (∂D) and g ∈ H -1 2 (∂D) the ITP(D, f, g, n) has a unique solution (u, v) ∈ L 2 (D) × L 2 (D) such that u -v ∈ H 2 (D).
This hypothesis is verified [START_REF] Sylvester | Discreteness of transmission eigenvalues via upper triangular compact operators[END_REF] for all k 2 ∈ R + except a countable set without finite point of accumulation if n verifies 1/(n -1) ∈ L ∞ (D) and (n -1) is either positive or negative definite in the neighbourhood of ∂D. We introduce the farfield pattern of the Green function:

φ z (x) := e -ikx•z
and the key ingredient of the GLSM:

Theorem 1. Assume that k / ∈ σ(D, n, 1)
. Then G is compact, injective with dense range and

φ z ∈ R(G) if and only if z ∈ D. Moreover, if z ∈ D then G(v) = φ z if and only if there exists u ∈ L 2 (D) such that (u, v) is a solution of ITP(D, Φ z , ∂Φz
∂ν , n, 1). We outline the main results of the GLSM in the case of noisy data (see [START_REF] Audibert | Identifying defects in an unknown background using differential measurements[END_REF]). The noisy operators, F δ and F δ # are such that F δ -F ≤ cδ and F δ # -F # ≤ δ where c is a real constant.Let g α,δ z ∈ L 2 (S d-1 ) be the minimizer of

J δ α (φ z ; g) := α(| F δ # g, g | + δα -η g 2 ) + F δ g -φ z 2 , (7)
for α > 0, δ > 0, η ∈]0, 1[ and φ z ∈ L 2 (S d-1 ). The functional

A α,δ (g) := | F δ # g, g | + α -η δ g 2 (8)
gives a characterisation of D through the following result.

Theorem 2. Under hypothesis 2 and 1 we have:

• z ∈ D implies lim sup α→0 lim sup δ→0 A α,δ (g α,δ z ) < ∞, • z / ∈ D implies lim inf α→0 lim inf δ→0 A α,δ (g α,δ z ) = ∞.
When we have two measurements campaigns, the same results applies to D0 = supp(n 0 -1) ⊂ D where A α,δ 0 is defined as above using F δ 0 (the farfield associated with n 0 and D 0 ) and

g α,δ 0,z = arg min g∈L 2 (S d-1 ) α(| F δ 0,# g, g | + δα -η g 2 ) + F δ 0 g -φ z 2 .
However we are interested in supp(n -n 0 ). The filtered difference term defined by:

D δ (g, g 0 ) := | F δ 0,# (g -g 0 ), g -g 0 | + δ g -g 0 2 ,
will be used to image the simply connected part of D 0 that have been modified between the two measurements. We denote this domain by D0 . Let Ω be the part of supp(n -n 0 ) such that Ω ∩ D0 = ∅. We introduce the following indicator function:

I α,δ T (g α,δ z , g α,δ 0,z ) = 1 A α,δ (g α,δ z ) 1 + A α,δ (g α,δ z )D δ (g α,δ z , g α,δ 0,z ) -1 , (9)
which will image a domain larger than the defect as follows:

Theorem 3. If we assume that k / ∈ σ(D, n, 1) ∪ σ(D 0 , n 0 , 1) ∪ σ(D, n, n 0 ) and, n and n 0 verify hypothesis 1, then z ∈ D0 ∪ Ω if and only if lim α→0 lim inf δ→0 I α,δ T (g α,δ z , g α,δ 0,z ) > 0.

Numerical Simulations

In order to fix the ideas, we shall limit ourselves to the two dimensional case and will introduce the algorithms for the discrete setting. We identify S 1 with the interval [0, 2π[. In order to collect the data of the inverse problem we solve numerically (3) for N incident fields u i ( 2πj N , •), j ∈ {0...N-1} using a finite element method implemented with Freefem++ [8]. The discrete version of F is then the matrix F := (u ∞ ( 2πj N , 2πk N )) 0≤j,k≤N-1 . We add some noise to the data to build a noisy far field matrix F δ where (F δ ) j,k = (F ) j,k (1 + σN ij ) for σ > 0 and N ij an uniform complex random variable in [-1, 1] 2 . We similarly generate F δ 0 We denote by φ z ∈ C N , the vector defined by φ z (j) = exp(-ik(z 1 cos( 2πj N )+z 2 sin( 2πj N )) for 0 ≤ j ≤ N-1. The analysis of previous sections suggests to consider

g GLSM z := argmin g∈L 2 (S 1 ) α (F δ ) 1 2 g 2 L 2 (S 1 ) + α 1-η δ g 2 L 2 (S 1 ) + F δ g -φ z 2 L 2 (S 1 ) .
The minimizer is explicitly given by

g GLSM z = (αF δ + α 1-η δId + F δ * F δ ) -1 F δ * φ z .
We similarly construct g GLSM 0,z using F δ 0 . In our numerical simulations we choose η = 0 (which corresponds to the one used in [START_REF] Audibert | A generalized formulation of the linear sampling method with exact characterization of targets in terms of farfield measurements[END_REF]) and set α with the same heuristic as in [START_REF] Audibert | A generalized formulation of the linear sampling method with exact characterization of targets in terms of farfield measurements[END_REF]. We then look at z → I α,δ T (g α,δ z , g α,δ 0,z ) as an indicator function. All our experiments are conducted for the background medium n 0 shown in Figure 1. This background medium is a simplified numerical description of a concrete like material. The wave frequency is 150kHz , the celerity of the medium is 4300m.s -1 (which means a wavelength of 2.87cm) and the celerity inside the inclusion is 5700m.s -1 .

Figure 1. The background medium, n = 0.57 inside the yellow inclusions Our theoretical analysis is only valid for inhomogeneous perturbations n and n 0 . One example of this setting is shown in Figure 2 where we modified the celerity (3 times higher) in two of the inclusions between the two measurements.

The main concern with concrete non destructive testing is the case of cracks. Our analysis do not include this case and its extension to it is the subject of an ongoing work. However the results shown in Figure 3 for a crack being either inside or outside the inhomogeneity n 0 gives promising results. To obtain this results we do not test with φ z but with its normal derivatives as explained in [7]. T (g α,δ z , g α,δ 0,z ),A α,δ 0 (g α,δ 0,z ) -1 and A α,δ (g α,δ z ) -1

Hypothesis 1 .

 1 The index of refraction n and the domain D satisfy n ∈ L ∞ (R d ), supp(n -1) = D, (n) ≥ 0 and there exist a constant n * > 0 such that (n(x) -1) ≥ n * for a.e. x ∈ D.

Figure 2 .

 2 Figure 2. From left to right and up to down : The index n, I α,δT (g α,δ z , g α,δ 0,z ),A α,δ 0 (g α,δ 0,z ) -1 and A α,δ (g α,δ z ) -1

Figure 3 .

 3 Figure 3. Cracks: On the left the medium with the cracks and on the right the corresponding I α,δT (g α,δ z , g α,δ 0,z )


L Audibert
email: lorenzo.audibert@edf.fr

H Haddar

A Girard

Qualitative non-destructive testing of concrete-like materials

In this paper we explore the capacity of Qualitative Inversion Methods to detect macroscopic cracks or a lattice of small cracks in concrete-like materials. These materials are difficult to probe since the heterogeneities size inside the medium and the wavelength of classically used sensors are of the same order of magnitude. We shall demonstrate how this difficulty can be avoided in the case of macroscopic cracks by using so-called differential measurements and application of the Linear Sampling Method. For a lattice of small cracks we rather propose to construct a macroscopic indicator based on the eigenvalues of a suitable transmission problem.

Introduction

We are interested in using elastic waves to perform complete non destructive testing of concrete-like materials. The main difficulty in controlling concrete is its heterogeneous nature. Concrete is made of cement paste, water and aggregates. After a drying period, this mixture results in a heterogeneous material. As far as waves propagation is concerned, the main characteristics is the difference of celerity between aggregates and cement paste, especially as the wavelength and the size of the aggregates are similar, thus, in what follows, we will model concrete as a biphasic material.This concrete-like material has the following properties: the celerity of pressure wave is 5700ms -1 in aggregates and 4300ms -1 in cement paste. Defect First a lattice of small cracks which are located at the interface between aggregates and cement paste. When this lattice is too dense, a macroscopic crack which has a length larger than the aggregates appears and grows until it reaches the surface. Those two types of defects are of interest and we will expose our preliminary results on how to detect them.

Although waves in concrete are elastic, we start with the simpler case of acoustic waves and postpone the treatment of the elastic one. We therefore assume that the pressure field, u is solution to the well-known Helmholtz equations:

     ∆u + k 2 nu = 0, u = u i + u s in R 2 \Γ lim r→∞ √ r( ∂u s ∂r -iku s ) = 0 ∂u ∂ν = 0 on Γ
where Γ is the crack(s) inside the medium, ν is a unit normal vector on Γ, n is the relative index with respect to the celerity in the air and u i is the incident field created by a point source located at the interface air-concrete. For our numerical simulations, in order to simulate the heterogeneities in concrete we used synthetic geometries generated by [START_REF] Escoda | Estimation of local stresses and elastic properties of a mortar sample by FFT computation of fields on a 3D image[END_REF] (See Figure 1, left) and eliminated the aggregates, which have an area smaller than λ 2 10 2 , where λ = 2π k . This results in the medium represented by Figure 1, right. According to the remarks above, the index of the aggregates equals 2, 8.10 -3 and of the cement equals 4, 8.10 -3 .

The inverse crack problem

We here investigate the inverse problem of finding cracks in concrete using the framework of the linear sampling method [START_REF] Cakoni | Qualitative Methods in Inverse Scattering Theory[END_REF]. First we concentrate on the macrocracks and then on the small cracks lattice. We use multistatic measurements with sensor (sources and receivers) located on the interface Σ between air and concrete, namely these measurements are u s (x, and Φ h (right)

Macroscopic Cracks

We introduce the near field operator for differential-measurements (i.e. measurements of u s and u s b ),

[N b g](•) = Σ [u s (•, x 0 ) -u s b (•, x 0 )]g(x 0 )ds(x 0 )
Following the framework of the Linear Sampling Method for cracks developed in [START_REF] Ben Hassen | Application of the linear sampling method to identify cracks with impedance boundary conditions[END_REF], one can find the support of the scattering crack, by looking for every z at the value of g z , where g z is the (regularized

) solution to [N b g z ](•) = ∂Φ b ∂ν(z) (z,
•) (We refer to [START_REF] Ben Hassen | Application of the linear sampling method to identify cracks with impedance boundary conditions[END_REF] for the details on how one copes with the unknown normal ν(z)). Using this algorithm in our configuration yields to a quite good reconstruction as demonstrated by Figure 2 (left). However, although in practice u s b can be known in a differential measurement framework, Φ b will always be unknown. Considering 2, right) that the quality of the reconstruction considerably deteriorates, although one can still distinguish the existence of a defect different from aggregates shape. This can be enough for qualitative inspection. We shall explore in the near future how one can improve the reconstruction by optimizing the choice of the background as in [START_REF] Haddar | Numerical analysis of the Factorization Method for EIT with piece wise constant uncertain background[END_REF]. We shall also numerically analyze the feasibility of control with non-differential measurements. Finally, on the theoretical level, sensors at the interface are a realistic set up which is not yet justified for the LSM and we shall seek for a formal justification.

g z the solution of [N b g z ](•) = ∂Φ h ∂ν(z) (z, •) shows (Figure

The case of lattice of microscopic cracks

When the cracks are small and numerous the LSM algorithm would not be able to locate them even in a differential setting. We shall rather use a macroscopic indicator based on the interior transmission eigenvalues which are computable from the measurements using following near field operator:

[N h g](•) = Σ [u s (•, x 0 ) -u s h (•, x 0 )]g(x 0 )ds(x 0 )
With this operator, since it also contains the contribution of the aggregates in the scattering effect, LSM will produce a support that covers almost all the heterogeneous material. It is well known [START_REF] Cakoni | Qualitative Methods in Inverse Scattering Theory[END_REF] that the LSM fails for some frequencies which are, in our case, the eigenvalues of the following transmission problem:

       ∆v + k 2 v = 0 in D ∆u + nk 2 u = 0 in D\Γ u = v, ∂u ∂ν = ∂v ∂ν
on ∂D ∂u ∂ν = 0 on Γ where D will be the scatterer (including the aggregates convex support and the cracks). Motivated by [START_REF] Cakoni | Transmission eigenvalues for inhomogeneous media containing obstacles[END_REF], we shall investigate the evolution of those frequencies with respect to the presence of cracks. In order to ensure that first transmission eigenvalues stay in the frequency range of the sources, we may introduce an artificial contrast in the equation solved by v, (characterized by an index n a = n) localized in a region of size comparable to the incident wavelength. We finally recall that our interior transmission problem is still open in terms of existence and evolution of the eigenvalues with respect to Γ.

Figure 1 :

 1 Figure 1: Simulated concrete and numerical set-up for direct and inverse problem

  x 0 ) x, x 0 ∈ Σ. In the sequel we introduce the subscript b which indicates the solution of the direct problem without defect and the subscript h which indicates the direct problem without defects and aggregates. The function Φ(z, •) denotes the fundamental outgoing solution with Dirac source at position z.

Figure 2 :

 2 Figure 2: Identification of cracks using Φ b (left) and Φ h (right)