Wolfgang Steiner
email: steiner@liafa.jussieu.fr

On the structure of (-β)-integers

The (-β)-integers are natural generalisations of the β-integers, and thus of the integers, for negative real bases. When β is the analogue of a Parry number, we describe the structure of the set of (-β)-integers by a fixed point of an anti-morphism.

Introduction

The aim of this paper is to study the structure of the set of real numbers having a digital expansion of the form

n-1 k=0 a k (-β) k ,
where (-β) is a negative real base with β > 1, the digits a k ∈ Z satisfy certain conditions specified below, and n ≥ 0. These numbers are called (-β)-integers, and have been recently studied by Ambrož, Dombek, Masáková and Pelantová [START_REF] Ambrož | Numbers with integer expansion in the numeration system with negative base[END_REF].

Before dealing with these numbers, we recall some facts about β-integers, which are the real numbers of the form

± n-1 k=0 a k β k such that 0 ≤ m-1 k=0 a k β k < β m for all 1 ≤ m ≤ n ,
i.e., n-1 k=0 a k β k is a greedy β-expansion. Equivalently, we can define the set of β-integers as

Z β = Z + β ∪ (-Z + β ) with Z + β = n≥0 β n T -n β (0) ,
where T β is the β-transformation, defined by

T β : [0, 1) → [0, 1) , x → βx -⌊βx⌋ .
This map and the corresponding β-expansions were first studied by Rényi [START_REF] Rényi | Representations for real numbers and their ergodic properties[END_REF].

The notion of β-integers was introduced in the domain of quasicrystallography, see for instance [START_REF] Burdík | Beta-integers as natural counting systems for quasicrystals[END_REF], and the structure of the β-integers is very well understood now. We have β Z β ⊆ Z β , the set of distances between consecutive elements of Z β is

∆ β = {T n β (1 -) | n ≥ 0}
, where T n β (x -) = lim y→x-T n β (y), and the sequence of distances between consecutive elements of Z + β is coded by the fixed point of a substition, see [START_REF] Fabre | Substitutions et β-systèmes de numération[END_REF] for the case when ∆ β is a finite set, that is when β is a Parry number. We give short proofs of these facts in Section 2. More detailed properties of this sequence can be found e.g. in [START_REF] Balková | Asymptotic behavior of beta-integers[END_REF][START_REF] Balková | Sequences with constant number of return words[END_REF][START_REF] Bernat | On a class of infinite words with affine factor complexity[END_REF][START_REF] Frougny | Complexity of infinite words associated with betaexpansions[END_REF][START_REF] Klouda | Factor complexity of infinite words associated with non-simple Parry numbers[END_REF].

Closely related to Z + β are the sets

S β (x) = n≥0 β n T -n β (x) (x ∈ [0, 1)),
which were used by Thurston [START_REF] Thurston | Groups, tilings and finite state automata[END_REF] to define (fractal) tilings of R d-1 when β is a Pisot number of degree d, i.e., a root > 1 of a polynomial

x d + p 1 x d-1 + • • • + p d ∈ Z[x]
such that all other roots have modulus < 1, and an algebraic unit, i.e., p d = ±1. These tilings allow e.g. to determine the k-th digit a k of a number without knowing the other digits, see [START_REF] Kalle | Beta-expansions, natural extensions and multiple tilings associated with Pisot units[END_REF]. It is widely agreed that the greedy β-expansions are the natural representations of real numbers in a real base β > 1. For the case of negative bases, the situation is not so clear. Ito and Sadahiro [START_REF] Ito | Beta-expansions with negative bases[END_REF] proposed recently to use the (-β)-transformation defined by T -β : -β β+1 , 1 β+1 , x → -βx -β β+1βx , see also [START_REF] Frougny | On negative bases[END_REF]. This transformation has the important property that T -β (-x/β) = x for all x ∈ -β β+1 , 1 β+1 . Some instances are depicted in Figures 1,3, 4 and 5. The set of (-β)-integers is therefore defined by

Z -β = n≥0 (-β) n T -n -β (0) .
These are the numbers

n-1 k=0 a k (-β) k such that -β β + 1 ≤ m-1 k=0 a k (-β) k-m < 1 β + 1 for all 1 ≤ m ≤ n .
Note that, in the case of β-integers, we have to add -Z + β to Z + β in order to obtain a set resembling Z. In the case of (-β)-integers, this is not necessary because the (-β)-transformation allows to represent positive and negative numbers.

It is not difficult to see that Z -β = Z = Z β when β ∈ Z, β ≥ 2. Some other properties of Z -β can be found in [START_REF] Ambrož | Numbers with integer expansion in the numeration system with negative base[END_REF], mainly for the case when

T n -β -β β+1 ≤ 0 and T 2n-1 -β -β β+1 ≥ 1-⌊β⌋ β for all n ≥ 1. (Note that T n -β -β β+1 ∈ 1 β+1 -⌊β⌋ β , 1-⌊β⌋ β ∪ -β -1 β+1 , 0 implies T n+1 -β -β β+1 > 0.) The set V β = T n -β -β
β+1 | n ≥ 0 plays a similar role for (-β)-expansions as the set {T n β (1 -) | n ≥ 0} for β-expansions. If V β is a finite set, then we call β > 1 an Yrrap number. Note that these numbers are called Ito-Sadahiro numbers in [START_REF] Masáková | Ito-Sadahiro numbers vs[END_REF], in reference to [START_REF] Ito | Beta-expansions with negative bases[END_REF]. However, as the generalised β-transformations in [START_REF] Góra | Invariant densities for generalized β-maps[END_REF] with E = (1, . . . , 1) are, up to conjugation by the map x → 1 β+1x, the same as our (-β)-transformations, these numbers were already considered by Góra and perhaps by other authors. Therefore, the neutral but intricate name (-β)-numbers was chosen in [START_REF] Liao | Dynamical properties of the negative beta-transformation[END_REF], referring to the original name β-numbers for Parry numbers [START_REF] Parry | On the β-expansions of real numbers[END_REF]. The name Yrrap number, used in the present paper, refers to the connection with Parry numbers and to the fact that T -β is (locally) orientation-reversing.

For any Yrrap number β ≥ (1 + √ 5)/2, we describe the sequence of (-β)-integers in terms of a two-sided infinite word on a finite alphabet which is a fixed point of an antimorphism (Theorem 3). Note that the orientation-reversing property of the map x → -βx imposes the use of an anti-morphism instead of a morphism, and that anti-morphisms were considered in a similar context by Enomoto [START_REF] Enomoto | AH-substitution and Markov partition of a group automorphism on T d[END_REF].

For 1 < β < 1+ √ 5

2 , we have Z -β = {0}, as already proved in [START_REF] Ambrož | Numbers with integer expansion in the numeration system with negative base[END_REF]. However, our study still makes sense for these bases (-β) because we can also describe the sets

S -β (x) = lim n→∞ (-β) n T -n -β (x) x ∈ -β β+1 , 1 β+1 ,
where the limit set consists of the numbers lying in all but finitely many sets (-β) n T -n -β (x), n ≥ 0. Taking the limit instead of the union over all n ≥ 0 implies that every y ∈ R lies in exactly one set

S -β (x), x ∈ -β β+1 , 1 β+1 , see Lemma 2. Note that T 2 -β -β -1 β+1 = -β β+1 when β ∈ Z. Other properties of the (-β)-transformation for 1 < β < 1+ √ 5 2
are exhibited in [START_REF] Liao | Dynamical properties of the negative beta-transformation[END_REF].

β-integers

In this section, we consider the structure of β-integers. The results are not new, but it is useful to state and prove them in order to understand the differences with (-β)-integers.

Recall that ∆ β = {T n β (1 -) | n ≥ 0}, and let ∆ * β be the free monoid generated by ∆ β . Elements of ∆ * β will be considered as words on the alphabet ∆ β , and the operation is the concatenation of words. The β-substitution is the morphism

ϕ β : ∆ * β → ∆ * β , defined by ϕ β (x) = 11 • • • 1 ⌈βx⌉-1 times T β (x -) (x ∈ ∆ β ).
Here, 1 is an element of ∆ β and not the identity element of ∆ * β (which is the empty word). Recall that, as ϕ β is a morphism, we have

ϕ β (uv) = ϕ β (u)ϕ β (v) for all u, v ∈ ∆ * β . Since ϕ n+1 β (1) = ϕ n β (ϕ β (1)
) and ϕ β (1) starts with 1, ϕ n β (1) is a prefix of ϕ n+1 β (1) for every n ≥ 0.

Theorem 1. For any β > 1, the set of non-negative β-integers takes the form

Z + β = {z k | k ≥ 0} with z k = k j=1 u j ,
where u 1 u 2 • • • is the infinite word with letters in ∆ β which has ϕ n β (1) as prefix for all n ≥ 0. The set of differences between consecutive β-integers is ∆ β .

The main observation for the proof of the theorem is the following. We use the notation

|v| = k and L(v) = k j=1 v j for any v = v 1 • • • v k ∈ ∆ k β , k ≥ 0. Lemma 1. For any n ≥ 0, 1 ≤ k ≤ |ϕ n β (1)|, we have T n β z k-1 β n , z k β n = [0, u k ) ,
and

z |ϕ n β (1)| = L(ϕ n β (1)) = β n . Proof. For n = 0, we have |ϕ 0 β (1)| = 1, z 0 = 0, z 1 = 1, u 1 = 1
, thus the statements are true. Suppose that they hold for n, and consider

u 1 u 2 • • • u |ϕ n+1 β (1)| = ϕ n+1 β (1) = ϕ β (ϕ n β (1)) = ϕ β (u 1 ) ϕ β (u 2 ) • • • ϕ β (u |ϕ n β (1)| ) . Let 1 ≤ k ≤ |ϕ n+1 β (1)|, and write u 1 • • • u k = ϕ β (u 1 • • • u j-1 ) v 1 • • • v i with 1 ≤ j ≤ |ϕ n β (1)|, 1 ≤ i ≤ |ϕ β (u j )|, i.e., v 1 • • • v i is a non-empty prefix of ϕ β (u j ).
For any x ∈ (0, 1], we have

T β (x -) = βx -⌈βx⌉ + 1, hence L(ϕ β (x)) = βx for x ∈ ∆ β . This yields that z k = L(u 1 • • • u k ) = β L(u 1 • • • u j-1 ) + L(v 1 • • • v i ) = βz j-1 + i -1 + v i and z k-1 = βz j-1 + i -1, hence z k-1 β , z k β = z j-1 + i -1 β , z j-1 + i -1 + v i β ⊆ [z j-1 , z j-1 + u j ) = [z j-1 , z j ) , T n+1 β z k-1 β n+1 , z k β n+1 = T β i -1 β , i -1 + v i β = [0, v i ) = [0, u k ) .
Moreover, we have L(ϕ n+1 β (1)) = β L(ϕ n β (1)) = β n+1 , thus the statements hold for n+1. Proof of Theorem 1. By Lemma 1, we have

z |ϕ n β (1)| = β n for all n ≥ 0, thus [0, 1) is split into the intervals [z k-1 /β n , z k /β n ), 1 ≤ k ≤ |ϕ n β (1)|. Therefore, Lemma 1 yields that T -n β (0) = {z k-1 /β n | 1 ≤ k ≤ |ϕ n β (1)|} , hence n≥0 β n T -n β (0) = {z k | k ≥ 0} . Since u k ∈ ∆ β for all k ≥ 1 and u |ϕ n (1)| = T n β (1 -) for all n ≥ 0, we have {z k -z k-1 | k ≥ 1} = {u k | k ≥ 1} = ∆ β .
For the sets S β (x), Lemma 1 gives the following corollary.

Corollary 1. For any x ∈ [0, 1), we have

S β (x) = {z k + x | k ≥ 0, u k+1 > x} ⊆ x + S β (0) .
In particular, we have S β (x)x = S β (y)y for all x, y ∈ [0, 1) with (x, y] ∩ ∆ β = ∅. From the definition of S β (x) and since x ∈ β T -1 β (x), we also get that

S β (x) = y∈T -1 β (x)
β S β (y)

x ∈ [0, 1) .

This shows that S β (x) is the solution of a graph-directed iterated function system (GIFS) when β is a Parry number, cf. [15, Section 3.2].

(-β)-integers

We now turn to the discussion of (-β)-integers and sets

S -β (x), x ∈ -β β+1 , 1 β+1 . Lemma 2. For any β > 1, x ∈ -β β+1 , 1 β+1 , we have S -β (x) = n≥0 (-β) n T -n -β (x) \ -β β+1 = y∈T -1 -β (x) (-β) S -β (y) .
For any y ∈ R, there exists a unique x

∈ -β β+1 , 1 β+1 such that y ∈ S -β (x). If T -β (x) = x, then S -β (x) = n≥0 (-β) n T -n -β (x), in particular S -β (0) = Z -β . Proof. If y ∈ S -β (x), then we have y (-β) n ∈ T -n -β (x) for all sufficiently large n, thus y ∈ (-β) n T -n -β (x) \ -β β+1
for some n ≥ 0. On the other hand, y ∈ (-β

) n T -n -β (x) \ -β β+1 for some n ≥ 0 implies that T m -β ( y (-β) m ) = T n -β ( y (-β) n ) = x for all m ≥ n, thus y ∈ S -β (x). This shows the first equation. Since x ∈ -β β+1 , 1 β+1 implies that x ∈ (-β) T -1 -β (x)\ -β β+1 , we obtain that S -β (x) = y∈T -1 -β (x) (-β) S -β (y) for all x ∈ -β β+1 , 1 β+1 . For any y ∈ R, we have y ∈ S -β T n -β y (-β) n for all n ≥ 0 such that y (-β) n ∈ -β β+1 , 1 β+1 , thus y ∈ S -β (x) for some x ∈ -β β+1 , 1 β+1 . To show that this x is unique, let y ∈ S -β (x) and y ∈ S -β (x ′ ) for some x, x ′ ∈ -β β+1 , 1 β+1 . Then we have y ∈ (-β) n T -n -β (x) \ -β β+1 and y ∈ (-β) m T -m -β (x ′ ) \ -β β+1 for some m, n ≥ 0, thus x = T n -β y (-β) n = T m -β y (-β) m = x ′ . If T n -β -β β+1 = x = T -β (x), then T n+2 -β -β -1 β+1 = T n+1 -β -β β+1 = T -β (x) = x yields that (-β) n -β β+1 ∈ S -β (x), which shows that S -β (x) = n≥0 (-β) n T -n -β (x) when T -β (x) = x.
The first two statements of the following proposition can also be found in [START_REF] Ambrož | Numbers with integer expansion in the numeration system with negative base[END_REF].

Proposition 1. For any β > 1, we have (-β) Z -β ⊆ Z -β . If β < (1 + √ 5)/2, then Z -β = {0}. If β ≥ (1 + √ 5)/2, then Z -β ∩ (-β) n [-β, 1] = (-β) n , (-β) n+1 ∪ (-β) n+2 T -n-2 -β (0) ∩ -1 β , 1 β 2
for all n ≥ 0, in particular

Z -β ∩ [-β, 1] = {-β, -β + 1, . . . , -β + ⌊β⌋, 0, 1} if β 2 ≥ ⌊β⌋(β + 1), {-β, -β + 1, . . . , -β + ⌊β⌋ -1, 0, 1} if β 2 < ⌊β⌋(β + 1). Proof. The inclusion (-β) Z -β ⊆ Z -β is a consequence of Lemma 2 and 0 ∈ T -1 -β (0). If β < 1+ √ 5 2 , then -1 β < -β β+1 , hence T -1 -β (0) = {0} and Z -β = {0}, see Figure 1 (right). If β ≥ 1+ √ 5 2 , then -1 β ∈ T -1 -β (0) implies 1 ∈ Z -β , thus (-β) n ∈ Z -β for all n ≥ 0. Clearly, (-β) n+2 T -n-2 -β (0) ∩ -1 β , 1 β 2 ⊆ Z -β ∩ (-β) n (-β, 1) . To show the other inclusion, let z ∈ (-β) m T -m -β (0) ∩ (-β) n (-β, 1) for some m ≥ 0. If z = (-β) m -β β+1 , then z (-β) m ∈ -β β+1 , 1 β+1 and z (-β) n+2 ∈ -1 β , 1 β 2 ⊆ -β β+1 , 1 β+1 imply that T n+2 -β z (-β) n+2 = T m -β z (-β) m = 0. If z = (-β) m -β β+1 , then T n+2 -β z (-β) n+2 = T n+2 -β (-β) m-n-1 β+1 = T m+2 -β -β -1 β+1 = T m+1 -β -β β+1 = T -β (0) = 0 , where we have used that z (-β) n+2 ∈ -β β+1 , 1 β+1 implies m ≤ n. Therefore, we have z ∈ (-β) n+2 T -n-2 -β (0) for all z ∈ Z -β ∩ (-β) n (-β, 1). Consider now n = 0, then Z -β ∩ [-β, 1] = {-β, 1} ∪ {z ∈ (-β, 1) | T 2 -β (z/β 2 ) = 0} . Since -⌊β⌋ β ≥ -β β+1 if and only if β 2 ≥ ⌊β⌋(β + 1), we obtain that (-β) T -1 -β (0) = {0, 1, . . . , ⌊β⌋} if β 2 ≥ ⌊β⌋(β + 1), {0, 1, . . . , ⌊β⌋ -1} if β 2 < ⌊β⌋(β + 1). If T 2 -β (z/β 2 ) = 0, then z = -a 1 β + a 0 with a 0 ∈ (-β) T -1 -β (0), a 1 ∈ {0, 1, . . . , ⌊β⌋}, and
Z -β ∩ [-β, 1] consists of those numbers -a 1 β + a 0 lying in [-β, 1].
Proposition 1 shows that the maximal difference between consecutive (-β)-integers exceeds 1 whenever β 2 < ⌊β⌋(β +1), i.e., T -β -β β+1 < 0. For an example, this was also proved in [START_REF] Ambrož | Numbers with integer expansion in the numeration system with negative base[END_REF]. In Examples 3 and 4, we see that the distance between two consecutive (-β)-integers can be even greater than 2, and the structure of Z -β can be quite complicated. Therefore, we adapt a slightly different strategy as for Z β .

In the following, we always assume that the set

V ′ β = V β ∪ {0} = T n -β -β β+1 | n ≥ 0 ∪ {0} is finite, i.e.
, β is an Yrrap number, and let β be fixed. For

x ∈ V ′ β , let r x = min y ∈ V ′ β ∪ 1 β+1 | y > x , x = x+rx 2 , J x = {x} and J x = (x, r x ) . Then {J a | a ∈ A β } forms a partition of -β β+1 , 1 β+1 , where A β = V ′ β ∪ V ′ β , with V ′ β = { x | x ∈ V ′ β } . We have T -β (J x ) = J T -β (x)
for every x ∈ V ′ β , and the following lemma shows that the image of every J

x , x ∈ V ′ β , is a union of intervals J a , a ∈ A β . Lemma 3. Let x ∈ V ′ β and write J x ∩ T -1 -β (V ′ β ) = {y 1 , . . . , y m } , with x = y 0 < y 1 < • • • < y m < y m+1 = r x . For any 0 ≤ i ≤ m, we have T -β (y i , y i+1 ) = J x i with x i = lim y→y i+1 - T -β (y) , i.e., x i = T -β y i +y i+1 2
, and β(y i+1y i ) = λ(J x i ), where λ denotes the Lebesgue measure.

Proof. Since T -β maps no point in (y i , y i+1 ) to -β β+1 ∈ V ′ β , the map is continuous on this interval and λ(T -β ((y i , y i+1 ))) = β(y i+1y i ). We have

x i ∈ V ′ β since x i = T -β (y i+1 ) in case y i+1 < 1
β+1 , and

x i = -β β+1 in case y i+1 = 1 β+1 . Since y i = max{y ∈ T -1 -β (V ′ β ) | y < y i+1 }, we obtain that r x i = lim y→y i + T -β (y), thus T -β ((y i , y i+1 )) = (x i , r x i ).
In view of Lemma 3, we define an anti-morphism

ψ β : A * β → A * β by ψ β (x) = T -β (x) and ψ β ( x) = x m T -β (y m ) • • • x 1 T -β (y 1 ) x 0 (x ∈ V ′ β
), with m, x i and y i as in Lemma 3. Here, anti-morphism means that

ψ β (uv) = ψ β (v)ψ β (u) for all u, v ∈ A * β . Now, the last letter of ψ β ( 0) is t, with t = max{x ∈ V β | x < 0}, and the first letter of ψ β ( t ) is 0. Therefore, ψ 2n β ( 0) is a prefix of ψ 2n+2 β ( 0) = ψ 2n β (ψ 2 β ( 0)) and ψ 2n+1 β ( 0) is a suffix of ψ 2n+3 β ( 0) for every n ≥ 0.
Theorem 2. For any Yrrap number β ≥ (1 + √ 5)/2, we have

Z -β = {z k | k ∈ Z, u 2k = 0} with z k = k j=1 λ(J u 2j-1 ) if k ≥ 0 , -|k| j=1 λ(J u -2j+1 ) if k ≤ 0 , where • • • u -1 u 0 u 1 • • • is the two-sided infinite word on the finite alphabet A β such that u 0 = 0, ψ 2n β ( 0) is a prefix of u 1 u 2 • • • and ψ 2n+1 β ( 0) is a suffix of • • • u -2 u -1 for all n ≥ 0. Note that • • • u -1 u 0 u 1 • • • is a fixed point of ψ β , with u 0 being mapped to u 0 .
The following lemma is the analogue of Lemma 1. We use the notation

L(v) = k j=1 λ(J v j ) if v = v 1 • • • v k ∈ A k β . Note that u 2k ∈ V ′ β and u 2k+1 ∈ V ′ β for all k ∈ Z, thus λ(J u 2k ) = 0 for all k ∈ Z. Lemma 4. For any n ≥ 0, 0 ≤ k < |ψ n β ( 0)|/2, we have T n -β z (-1) n k (-β) n = u (-1) n 2k , T n -β z (-1) n k (-β) n , z (-1) n (k+1) (-β) n = J u (-1) n (2k+1) ,
and

z (-1) n (|ψ n β ( 0)|+1)/2 = (-1) n L ψ n β ( 0) = λ(J 0 ) (-β) n = r 0 (-β) n .
Proof. The statements are true for n = 0 since |ψ 0 β ( 0)| = 1, z 0 = 0, z 1 = λ(J 0 ) = r 0 . Suppose that they hold for even n, and consider

u -|ψ n+1 β ( 0)| • • • u -2 u -1 = ψ n+1 β ( 0) = ψ β ψ n β ( 0) = ψ β (u |ψ n β ( 0)| ) • • • ψ β (u 2 )ψ β (u 1 ) . Let 0 ≤ k < |ψ n+1 β ( 0)|/2,
and write

u -2k-1 • • • u -1 = v -2i-1 • • • v -1 ψ β (u 1 • • • u 2j ) with 0 ≤ j < |ψ n β ( 0)|/2, 0 ≤ i < |ψ β (u 2j+1 )|/2, i.e., u -2i-1 • • • u -1 is a suffix of ψ β (u 2j+1 ). By Lemma 3, we have L(ψ β ( x)) = β λ(J x ) for any x ∈ V ′ β . This yields that -z -k-1 = β L(u 1 • • • u 2j ) + L(v -2i-1 • • • v -1 ) = β z j + L(v -2i-1 • • • v -1 ) and -z -k = β z j + L(v -2i • • • v -1 )
. By the induction hypothesis, we obtain that

T n+1 -β z -k (-β) n+1 = T n+1 -β z j (-β) n - L(v -2i • • • v -1 ) (-β) n+1 = T -β (u 2j ) = ψ β (u 2j ) = u -2k if i = 0, T -β x + L(v -2i • • • v -1 )/β = T -β (y i ) = v -2i = u -2k if i ≥ 1,
where the y i 's are the numbers from Lemma 3 for x = u 2j+1 , and

T n+1 -β z -k (-β) n+1 , z -k-1 (-β) n+1 = T -β (y i , y i+1 ) = J v -2i-1 = J u -2k-1 .
Moreover, we have L(ψ n+1 β ( 0)) = β L(ψ n β ( 0)) = r 0 β n+1 , thus the statements hold for n + 1. The proof for odd n runs along the same lines and is therefore omitted.

Proof of Theorem 2. By Lemma 4, we have z (-1) n (|ψ n β ( 0)|+1)/2 = r 0 (-β) n for all n ≥ 0, thus [0, r 0 ) splits into the intervals z (-1) n k (-β) -n , z (-1) n (k+1) (-β) -n and points z (-1

) n k (-β) -n , 0 ≤ k < |ψ n β ( 0)|/2, hence T -n -β (0) ∩ [0, r 0 ) = z (-1) n k (-β) -n | 0 ≤ k < |ψ n β ( 0)|/2, u (-1) n 2k = 0 . Let m ≥ 1 be such that β 2m r 0 ≥ 1 β+1 . Then we have -β β+1 , 1 β+1 ⊆ (-β 2m+1 r 0 , β 2m r 0 ), and
T -n -β (0) \ -β β+1 ⊆ (-β) 2m T -n-2m -β (0) ∩ [0, r 0 ) ∪ (-β) 2m+1 T -n-2m-1 -β (0) ∩ [0, r 0 ) , thus n≥0 (-β) n T -n -β (0) \ -β β+1 = n≥0 (-β) n T -n -β (0) ∩ [0, r 0 ) = {z k | k ∈ Z, u 2k = 0} .
Together with Lemma 2, this proves the theorem.

As in the case of positive bases, the word • • • u -1 u 0 u 1 • • • also describes the sets S -β (x). Theorem 2 and Lemma 4 give the following corollary.

Corollary 2. For any

x ∈ V ′ β , y ∈ J x , we have S -β (x) = {z k | k ∈ Z, u 2k = x} and S -β (y) = {z k + y -x | k ∈ Z, u 2k+1 = x} .
Lemma 2 and Corollary 2 imply that S -β (x) is the solution of a GIFS for any Yrrap number β ≥ (1 + √ 5)/2, x ∈ -β β+1 , 1 β+1 , cf. the end of Section 2. Recall that our main goal is to understand the structure of Z -β (in case β ≥ (1+ √ 5)/2), i.e., to describe the occurrences of 0 in the word • • • u -1 u 0 u 1 • • • defined in Theorem 2 and the words between two successive occurrences. Let

R β = {u 2k u 2k+1 • • • u 2s(k)-1 | k ∈ Z, u 2k = 0} with s(k) = min{j ∈ Z | u 2j = 0, j > k} be the set of return words of 0 in • • • u -1 u 0 u 1 • • • . Note that s(k) is defined for all k ∈ Z since (-β) n ∈ Z -β for all n ≥ 0 by Proposition 1.
For any w ∈ R β , the word ψ β (w0) is a factor of • • • u -1 u 0 u 1 • • • starting and ending with 0, thus we can write ψ β (w0) = w 1 • • • w m 0 with w j ∈ R β , 1 ≤ j ≤ m, and set

ϕ -β (w) = w 1 • • • w m .
This defines an anti-morphism ϕ -β : R * β → R * β , which plays the role of the β-substitution.

Since • • • u -1 u 0 u 1 • • • is generated by u 1 = 0, as described in Theorem 2, we consider w β = u 0 u 1 • • • u 2s(0)-1 . We have [0, 1] = 0, 1 β+1 ∪ 1 β+1 , 1 , T -β (-β) -1 1 β+1 , 1 = -β β+1 , 0 , thus L(w β ) = 1 and
w β = 0 0 x 1 x 1 • • • x m x m x -ℓ x -ℓ • • • x -1 x -1 ,
where the x j are defined by

V ′ β = {x -ℓ , . . . , x -1 , 0, x 1 , . . . , x m }, x -ℓ < • • • < x -1 < 0 < x 1 < • • • < x m .
Theorem 3. For any Yrrap number β ≥ (1 + √ 5)/2, we have

Z -β = {z ′ k | k ∈ Z} with z ′ k = k j=1 L(u ′ j ) if k ≥ 0 , -|k| j=1 L(u ′ -j ) if k ≤ 0 , where • • • u ′ -2 u ′ -1 u ′ 1 u ′ 2 • • • is the two-sided infinite word on the finite alphabet R β such that ϕ 2n -β (w β ) is a prefix of u ′ 1 u ′ 2 • • • and ϕ 2n+1 -β (w β ) is a suffix of • • • u ′ -2 u ′ -1 for all n ≥ 0. The set of distances between consecutive (-β)-integers is ∆ -β = {z ′ k+1 -z ′ k | k ∈ Z} = {L(w) | w ∈ R β } . Note that the index 0 is omitted in • • • u ′ -2 u ′ -1 u ′ 1 u ′ 2 • • • for reasons of symmetry. Proof. The definitions of • • • u -1 u 0 u 1 • • • in Theorem 2 and of ϕ -β imply that • • • u ′ -2 u ′ -1 u ′ 1 u ′ 2 • • • is the derived word of • • • u -1 u 0 u 1 • • • with respect to R β , that is u ′ k = u |u ′ 1 •••u ′ k-1 | • • • u |u ′ 1 •••u ′ k |-1 , u ′ -k = u -|u ′ -k •••u ′ -1 | • • • u -|u ′ 1-k •••u ′ -1 |-1 (k ≥ 1) with {|u ′ 1 • • • u ′ k-1 | | k ≥ 1} ∪ {-|u ′ -k • • • u ′ -1 | | k ≥ 1} = {k ∈ Z | u k = 0} .
Here, for any v ∈ R * β , |v| denotes the length of v as a word in

A * β , not in R * β . Since z ′ k = k j=1 L(u ′ j ) = |u ′ 1 •••u ′ k |-1 j=0 λ(J u j ) = |u ′ 1 •••u ′ k | j=1 λ(J u j ) , z ′ -k = - k j=1 L(u ′ -j ) = - |u ′ -k •••u ′ -1 | j=1 λ(J u -j ) for all k ≥ 0, Theorem 2 yields that {z ′ k | k ∈ Z} = Z -β . It follows from the definition of R β that ∆ -β = {L(w) | w ∈ R β }.
It remains to show that R β is a finite set. We first show that the restriction of ψ β to V ′ β is primitive, which means that there exists some m ≥ 1 such that, for every

x ∈ V ′ β , ψ m β ( x) contains all elements of V ′ β .
The proof is taken from [START_REF] Góra | Invariant densities for generalized β-maps[END_REF]Proposition 8]. If β > 2, then the largest connected pieces of images of J x under T -β grow until they cover two consecutive discontinuity points 1 β+1 -a+1 β , 1 β+1 -a β of T -β , and the next image covers all intervals J y ,

y ∈ V ′ β . If 1+ √ 5 2
< β ≤ 2, then β 2 > 2 implies that the largest connected pieces of images of J x under T 2 -β grow until they cover two consecutive discontinuity points of

T 2 -β . Since T 2 -β -β β+1 , β -2 β+1 -1 β = -β 3 +β 2 +β β+1 , 1 β+1 , T 2 -β β -2 β+1 -1 β , -β -1 β+1 = -β β+1 , β 2 -β-1 β+1 , T 2 -β -β -1 β+1 , β -2 β+1 = -β β+1 , 1 β+1 , T 2 -β β -2 β+1 , 1 β+1 = -β β+1 , β 2 -β-1 β+1
, the next image covers the fixed point 0, and further images grow until after a finite number of steps they cover all intervals J y ,

y ∈ V ′ β . The case β = 1+ √ 5 2 is treated in Example 1. If T n -β -β β+1 = 0 for all n ≥ 0, then T n -β is continuous at all points x ∈ -β β+1 , 1 β+1 with T n
-β (x) = 0, thus u 2k = 0 is equivalent to u 2k+1 = 0 (see also Proposition 2 below). Hence we can consider the return words of 0 in • • • u -1 u 0 u 1 • • • instead of the return words of 0. Since ψ m β ( x 0 x 1 x 2 ) has at least two occurrences of 0 for all x 0 , x 1 , x 2 ∈ V ′ β , there are only finitely many such return words. If T n -β -β β+1 = 0, then ψ n β (x 0 x 1 x 2 ) starts and ends with 0 for all x 0 , x 1 , x 2 ∈ V ′ β , hence R β is finite as well. For details on derived words of primitive substitutive words, we refer to [START_REF] Durand | A characterization of substitutive sequences using return words[END_REF]. We remark that, for practical reasons, the set R β can be obtained from the set R = {w β } by adding to R iteratively all return words of 0 which appear in ψ β (w0) for some w ∈ R until R stabilises. The final set R is equal to R β . Now, we apply the theorems in the case of two quadratic examples.

Example 1. Let β = 1+ √ 5 2 , i.e., β 2 = β + 1, and t = -β β+1 = -1 β , then V β = V ′ β = {t, 0}. Since J t = (t, 0) = t, -1 β 3 ∪ -1 β 3 ∪ -1 β 3 , 0 , J 0 = 0, 1 β 2 , see Figure 1 (middle), the anti-morphism ψ β on A * β is defined by ψ β : t → 0 , t → 0 t t , 0 → 0 , 0 → t . Its two-sided fixed point • • • u -1 u 0 u 1 • • • is • • • 0 ψ β (0) 0 t t ψ β ( t ) 0 ψ β (t) t ψ β ( 0) 0 ψ β (0) 0 t t ψ β ( t ) 0 ψ β (t) t ψ β ( 0) 0 ψ β ( 0) 0 t t ψ β ( t ) 0 ψ β (0) 0 t t ψ β ( t ) 0 ψ β (t) t ψ β ( 0) 0 ψ β (0) t ψ β ( 0) 0 ψ β (0) • • • ,
where 0 marks the central letter u 0 . The return word of 0 starting at u 0 is w β = 0 0 t t. The image ψ β (w β 0) = 0 0 t t 0 t 0 contains the return words w β and 0 t. Since ψ β (0 t 0) = 0 0 t t 0, there are no other return words of 0, i.e., R β = {A, B} with A = 0 0 t t, B = 0 t. Therefore,

• • • u ′ -2 u ′ -1 u ′ 1 u ′ 2 • • • is a two-sided fixed point of the anti-morphism ϕ -β : A → AB , B → A , with u ′ -13 • • • u ′ -1 u ′ 1 • • • u ′ 21 = AABAABABAABAB AABAABABAABAABABAABAB . We have λ(J 0 ) = 1 β 2 , λ(J t ) = 1 β , thus L(A) = 1, L(B) = 1 β = β -1
, and some (-β)-integers are shown in Figure 2. Note that (-β) n can also be represented as

(-β) n+2 + (-β) n+1 . -β 3 -β 3 +β 2 -β -β 3 +β 2 -β +1 -β -β +1 0 1 β 2 -β +1 β 2 β 4 -β 3 +β 2 -β β 4 -β 3 +β 2 -β +1 β 4 -β 3 +β 2 β 4 -β +1 β 4 A A B A B A A B A A B A B Figure 2. The (-β)-integers in [-β 3 , β 4 ], β = (1 + √ 5)/2. Example 2. Let β = 3+ √ 5
2 , i.e., β 2 = 3β -1, then the (-β)-transformation is depicted in Figure 3, where

t 0 = -β β+1 , t 1 = T -β (t 0 ) = β 2 β+1 -2 = -β -1 β+1 , T -β (t 1 ) = 1 β+1 -1 = t 0 . Therefore, V ′ β = {t 0 , t 1 , 0} and the anti-morphism ψ β : A * β → A * β is defined by ψ β : t 0 → t 1 , t 0 → t 0 t 1 t 1 0 0 t 0 t 0 , t 1 → t 0 , t 1 → 0 , 0 → 0 , 0 → t 0 t 1 t 1 ,
which has the two-sided fixed point

• • • 0 ψ β (0) 0 ψ β ( t 1 ) t 0 ψ β (t 1 ) t 0 t 1 t 1 0 0 t 0 t 0 ψ β ( t 0 ) t 1 ψ β (t 0 ) t 0 t 1 t 1 ψ β ( 0) 0 ψ β (0) 0 ψ β ( t 1 ) t 0 ψ β (t 1 ) t 0 t 1 t 1 0 0 t 0 t 0 ψ β ( t 0 ) • • • ,
where 0 marks the central letter u 0 . We have w β = 0 0 t 0 t 0 t 1 t 1 and

ψ β : 0 0 t 0 t 0 t 1 t 1 0 → 0 0 t 0 t 0 t 1 t 1 0 0 t 0 t 0 t 1 t 0 t 1 t 1 0 , 0 0 t 0 t 0 t 1 t 0 t 1 t 1 0 → 0 0 t 0 t 0 t 1 t 1 0 0 t 0 t 0 t 0 t 0 t 1 t 1 0 0 t 0 t 0 t 1 t 0 t 1 t 1 0 , 0 0 t 0 t 0 t 0 t 0 t 1 t 1 0 → 0 0 t 0 t 0 t 1 t 1 0 0 t 0 t 0 t 1 t 0 t 1 t 1 0 0 t 0 t 0 t 1 t 0 t 1 t 1 0 .
Note that 0 0 t 0 t 0 t 1 t 0 t 1 t 1 and 0 0 t 0 t 0 t 0 t 0 t 1 t 1 differ only by a letter in V ′ β , and correspond therefore to intervals of the same length. Since the letters t 0 and t 1 are never mapped to 0, we identify these two return words. This means that R β = {A, B} with A = 0 0 t 0 t 0 t 1 t 1 , B = 0 0 t 0 t 0 {t 0 , t 1 } t 0 t 1 t 1 , and

• • • u ′ -2 u ′ -1 u ′ 1 u ′ 2 • • • = • • • ABBABABBABBAB ABBABABBABBAB • • • is a two-sided fixed point of the anti-morphism ϕ -β : A → AB , B → ABB .
We have L(A) = 1, L(B) = β -1 > 1, and some (-β)-integers are shown in Figure 3.

0 0 t 0 t 0 1 β+1 1 β+1 t 1 t 1 J t 0 J t 1 J 0 -β 3 -β 3 +1 -β 3 +β 2 -2β +1 -β 3 +β 2 -β -β 3 +β 2 -β +1 -β 3 +β 2 -β 3 +β 2 +1 -β 3 +2β 2 -2β +1 -β 3 +2β 2 -β -β 3 +2β 2 -β +1 -2β +1 -β -β +1 0 1 β 2 -2β +1 β 2 -β β 2 -β +1 β 2 A B B A B A B B A B B A B A B B A B Figure 3. The (-β)-transformation and Z -β ∩ [-β 3 , β 2 ], β = (3 + √ 5)/2.
We remark that it is sufficient to consider the elements of V ′ β when one is only interested in Z -β . This is made precise in the following proposition.

Proposition 2. Let β and • • • u -1 u 0 u 1 • • • be as in Theorem 2, t = max{x ∈ V β | x < 0}.
For any k ∈ Z, u 2k = 0 is equivalent to u 2k-1 = t or u 2k+1 = 0.

If 0 ∈ V β or the size of V β is even, then u 2k = 0 is equivalent to u 2k-1 = t. If 0 ∈ V β or the size of V β is odd, then u 2k = 0 is equivalent to u 2k+1 = 0.

Proof. Let k ∈ Z and m ≥ 0 such that z k /β 2m ∈ -β β+1 , 1 β+1 . Then we have

• u 2k = 0 if and only if T 2m -β (z k /β 2m ) = 0, • u 2k-1 = t if and only if lim y→z k -T 2m -β (y/β 2m ) = 0, • u 2k+1 = 0 if and only if lim y→z k + T 2m -β (y/β 2m ) = 0. Thus u 2k = 0, u 2k-1 = t and u 2k+1 = 0 are equivalent when T 2m -β is continuous at z k /β 2m . Assume from now on that z k /β 2m is a discontinuity point of T 2m -β . Then T ℓ -β (z k /β 2m ) = -β β+1
for some 1 ≤ ℓ ≤ 2m and, if ℓ is minimal with this property, lim

y→z k - T 2⌊ℓ/2⌋+1 -β (y/β 2m ) = -β β+1
and lim

y→z k + T 2⌈ℓ/2⌉ -β (y/β 2m ) = -β β+1 .
Hence, if 0 ∈ V β , we cannot have u 2k = 0, u 2k-1 = t or u 2k+1 = 0 at a discontinuity point, which proves the proposition in this case. If 0 ∈ V β , then T

#V β -1 -β -β β+1 = 0, thus • T j -β (z k /β 2m ) = 0 if and only if j ≥ ℓ + #V β -1, • lim y→z k -T j -β (y/β 2m ) = 0 if and only if j ≥ 2⌊ℓ/2⌋ + #V β , • lim y→z k + T j
-β (y/β 2m ) = 0 if and only if j ≥ 2⌈ℓ/2⌉ + #V β -1. Since 2⌊ℓ/2⌋ ≥ ℓ -1 and 2⌈ℓ/2⌉ ≥ ℓ, we obtain u 2k = 0 whenever u 2k-1 = t or u 2k+1 = 0. If #V β is even, then u 2k = 0 implies that u 2k-1 = t since 2m ≥ ℓ + #V β -1 implies that 2m ≥ 2⌊ℓ/2⌋ + #V β . If #V β is odd, then u 2k = 0 implies that u 2k+1 = 0 since 2m ≥ ℓ + #V β -1 implies that 2m ≥ 2⌈ℓ/2⌉ + #V β -1. This proves the proposition.

By Proposition 2, it suffices to consider the anti-morphism

ψ β : V ′ β * → V ′ β * defined by ψ β ( x) = x m • • • x 1 x 0 when ψ β ( x) = x m T -β (y m ) • • • x 1 T -β (y 1 ) x 0 (x ∈ V ′ β ).
Then ∆ -β is given by the set R β which consists of the return words of 0 when 0 ∈ V β or the size of V β is odd. When 0 ∈ V β and the size of V β is even, as in Example 1, then R β consists of the words w t such that t w is a return word of t.

t 0 = -β β+1 , t 1 = β 2 β+1 -2 = β -1 -2 β+1 , t 2 = 2β-1 β+1 -1 = β -2 β+1 , t 3 = -β -1 β+1 , t 4 = 1 β+1 -1 = t 0 , see Figure 4. The anti-morphism ψ β : V ′ β * → V ′ β * is therefore defined by ψ β : t 0 → t 2 t 0 , t 1 → t 0 t 1 t 3 0 , t 3 → 0 t 2 , 0 → t 3 , t 2 → t 0 t 1 .
Since 0 ∈ V β , we consider return words of 0 in the ψ β -images of w β = 0 t 2 t 0 t 1 t 3 :

0 t 2 t 0 t 1 t 3 → 0 t 2 t 0 t 1 t 3 0 t 2 t 0 t 0 t 1 t 3 , 0 t 2 t 0 t 0 t 1 t 3 → 0 t 2 t 0 t 1 t 3 0 t 2 t 0 t 2 t 0 t 0 t 1 t 3 , 0 t 2 t 0 t 2 t 0 t 0 t 1 t 3 → 0 t 2 t 0 t 1 t 3 0 t 2 t 0 t 2 t 0 t 0 t 1 t 2 t 0 t 0 t 1 t 3 , 0 t 2 t 0 t 2 t 0 t 0 t 1 t 2 t 0 t 0 t 1 t 3 → 0 t 2 t 0 t 1 t 3 0 t 2 t 0 t 2 t 0 t 0 t 1 t 0 t 1 t 3 0 t 2 t 0 t 2 t 0 t 0 t 1 t 2 t 0 t 0 t 1 t 3 , 0 t 2 t 0 t 2 t 0 t 0 t 1 t 0 t 1 t 3 → 0 t 2 t 0 t 1 t 3 0 t 2 t 0 t 0 t 1 t 3 0 t 2 t 0 t 2 t 0 t 0 t 1 t 2 t 0 t 0 t 1 t 3 . Hence R β = {A, B, C, D, E} with A = 0 t 2 t 0 t 1 t 3 , B = 0 t 2 t 0 t 0 t 1 t 3 , C = 0 t 2 t 0 t 2 t 0 t 0 t 1 t 3 , D = 0 t 2 t 0 t 2 t 0 t 0 t 1 t 2 t 0 t 0 t 1 t 3 , E = 0 t 2 t 0 t 2 t 0 t 0 t 1 t 0 t 1 t 3 , and Z -β is described by the anti-morphism ϕ -β : R * β → R * β given by ϕ -β : A → AB , B → AC , C → AD , D → AED , E → ABD .
The (-β)-integers in [-β 3 , β 4 ] are represented in Figure 4, and we have

L(A) = 1, L(B) = β -1, L(C) = β 2 -β -1, L(D) = β 2 -β ≈ 2.659, L(E) = β.
Note that L(D) > β > 2. Moreover, the cardinality of ∆ -β is larger than that of V β , which in turn is larger than the algebraic degree β+1 . We have t 5 = -1 β+1 = t 6 . The anti-morphism

d of β (#∆ -β = 5, #V β = 4, d = 3). 0 0 t 0 t 0 1 β+1 1 β+1 t 3 t 1 t 1 t 2 t 2 t 3 J t 0 J t 1 J t 3 J 0 J t 2 -β 3 -β 3 +1 -β 3 +β 2 -β -β 3 +β 2 -β +1 -β 3 +β 2 -β 3 +β 2 +1 -β -β +1 0 1 β 2 -β β 2 -β +1 β 2 β 2 +1 β 4 -2β 3 +β 2 +1 β 4 -β 3 -β β 4 -β 3 -β +1 β 4 -β 3 β 4 -β 3 +1 β 4 -β 3 +β 2 -β β 4 -β 3 +β 2 -β +1 β 4 -β 3 +β 2 β 4 -β 3 +β 2 +1 β 4 -β β 4 -β +1 β 4 A C A B A D A B A C A B A E D A B A C A B A D A B
ψ β : V ′ β * → V ′ β * is therefore given by ψ β : t 0 → t 3 t 5 , t 2 → t 4 t 0 t 2 , t 3 → t 5 t 1 0 t 4 t 0 t 2 t 3 t 5 t 1 0 , t 5 → t 2 t 3 , t 1 → 0 t 4 t 0 , 0 → t 5 t 1 , t 4 → t 0 t 2 t 3 .
In order to deal with shorter words, we group the letters forming the words

a = 0 t 4 , b = t 0 t 2 t 3 t 5 t 1 , c = t 0 t 2 t 3 t 5 , d = t 2 t 3 t 5 t 1 , e = t 0 t 2 , f = t 4 , g = t 0 t 2 t 3 , h = t 5 t 1 ,
which correspond to the intervals J a = 0, 1 β+1 , J b = (t 0 , 0), J c = (t 0 , t 1 ), J d = (t 2 , 0), J e = (t 0 , t 3 ), J f = t 4 , 1 β+1 , J g = (t 0 , t 5 ), J h = (t 5 , 0), occurring in iterated images of J a . The anti-morphism ψ β acts on these words by

ψ β : a → b , b → ababac , c → dabac , d → ababae , e → f c , f → g , g → habac , h → ag .
J t 0 J t 2 J t 3 J t 5 J t 1 J 0 J t 4 -β -β +1 -β +2 0 1 2 β 2 -3β +1 β 2 -3β +2 β 2 -2β β 2 -2β +1 β 2 -2β +2 β 2 -β β 2 -β +1 β 2 -β +2 β 2 A A B A A C A B A A B A A B

Conclusions

With every Yrrap number β ≥ (1 + √ 5)/2, we have associated an anti-morphism ϕ -β on a finite alphabet. The distances between consecutive (-β)-integers are described by a fixed point of ϕ -β . In [START_REF] Ambrož | Numbers with integer expansion in the numeration system with negative base[END_REF], the anti-morphism is described explicitely for each β > 1 such that T n -β -β β+1 ≤ 0 and T 2n-1 -β -β β+1 ≥ 1-⌊β⌋ β for all n ≥ 1. Examples 3 and 4 show that the situation can be quite complicated when this condition is not fulfilled. Although ϕ -β can be obtained by a simple algorithm, it seems to be difficult to find a priori bounds for the number of different distances between consecutive (-β)-integers or for their maximal value. Only the case of quadratic Pisot numbers β is completely solved; here, we know from [START_REF] Ito | Beta-expansions with negative bases[END_REF][START_REF] Ambrož | Numbers with integer expansion in the numeration system with negative base[END_REF] that #V β = #∆ -β = 2.

Recall that the maximal distance between consecutive β-integers is 1, and the number of different distances is equal to the cardinality of the set {T n β (1 -) | n ≥ 0}. Example 3 shows that the (-β)-integers do not satisfy similar properties. By generalising Example 4 to β > 1 with β 6 = (m+1)β 5 +mβ 4 +mβ 3 +β 2 -mβ -1, m ≥ 2, one sees that the maximal distance can be arbitrarily close to 4 for algebraic integers of degree 6 and #V β = 6.

In a forthcoming paper, we associate anti-morphisms ϕ -β on infinite alphabets with non-Yrrap numbers β, by considering the intervals occurring in the iterated T -β -images of 0, 1 β+1 , cf. Example 4, and we show that the distances between consecutive (-β)-integers can be unbounded, e.g. for β > 1 satisfying -β β+1 = ∞ k=1 a k (-β) -k where a 1 a 2 • • • = 31232 1 2 31232 2 • • • is a fixed point of the morphism 3 → 31232, 2 → 2, 1 → 1. For Yrrap numbers β, this implies that there is no bound for the distance between consecutive (-β)-integers which is independent of β. However, large distances occur probably only far away from 0 and when #V β is large, and it would be interesting to quantify these relations.

Another topic that is worth investigating is the structure of the sets S -β (x) for x = 0, and of the corresponding tilings when β is a Pisot unit. A related question is whether Z -β can be given by a cut and project scheme, cf. [START_REF] Berthé | Tilings associated with beta-numeration and substitutions[END_REF][START_REF] Gazeau | Geometric study of the beta-integers for a Perron number and mathematical quasicrystals[END_REF].

Figure 1 . 2 ≈

 12 Figure 1. The (-β)-transformation for β = 2 (left), β = 1+ √ 5 2 ≈ 1.618 (middle), and β = 1 β + 1 β 2 ≈ 1.325 (right).

Example 3 .

 3 Let β > 1 with β 3 = 2β 2 + 1, i.e., β ≈ 2.206, and let t n = T n -β -β β+1 for n ≥ 0. Then we have

Figure 4 .Example 4 .

 44 Figure 4. The (-β)-transformation and Z -β ∩ [-β 3 , β 4 ] from Example 3. Example 4. Let β > 1 with β 6 = 3β 5 + 2β 4 + 2β 3 + β 2 -2β -1, i.e., β ≈ 3.695, then the (-β)-transformation is depicted in Figure 5, where t n = T n -β -β

Since 0

 0 only occurs at the beginning of a, the return words of 0 with their ψ β -images are ab → ab ab acb , aed → ab ab aef cb , acb → ab ab acd ab acb , aef cb → ab ab acd ab acgf cb , acd → ab ab aed ab acb , acgf cb → ab ab acd ab acgh = acb ab acd ab acb .Therefore, Z -β is described by the anti-morphism ϕ -β : R * β → R * β which is defined byϕ -β : A → AAB , L(A) = 1 , B → AACAB , L(B) = β -2 ≈ 1.695 , C → AADAB , L(C) = β 2 -3β -1 ≈ 1.569 , D → AAE , L(D) = β 3 -3β 2 -2β -1 ≈ 1.104 , E → AACAF , L(E) = β 4 -3β 3 -2β 2β -2 ≈ 2.081 , F → AACABACAB , L(F ) = β 5 -3β 4 -2β 3 -2β

Figure 5 .

 5 Figure 5. The (-β)-transformation and Z -β ∩ [-β, β 2 ] from Example 4.

  2 + β -2 ≈ 3.12 .

	1
	β+1
	0 t 4
	t 1
	t 5
	t 3
	t 2
	t 0	t 0	t 2	t 3	t 5	t 1	0	t 4	1 β+1


V Laurain

M Gilson

R Tóth
email: r.toth@tudelft.nl

H Garnier

Direct identification of continuous-time LPV models

Controllers in the linear parameter-varying (LPV) framework are commonly designed in continuous-time (CT) requiring accurate and low-order CT models of the system. Nonetheless, most of the methods dedicated to the identification of LPV systems are addressed in discrete-time (DT) settings. In practice when discretizing models which are naturally expressed in CT, the dependency on the scheduling variables becomes non-trivial and over-parameterized. Consequently, direct identification of CT LPV systems in an input-output setting is investigated. To provide consistent model parameter estimates in this setting, a refined instrumental variable (IV) approach is proposed. The statistical properties of this approach is illustrated through a relevant Monte Carlo simulation example.

I. INTRODUCTION

The framework of linear parameter-varying (LPV) systems was introduced in the 1990s with the purpose to handle in a simple but efficient way the often nonlinear or timevarying nature of systems encountered in practice. The LPV system class forms an intermediate step between linear timeinvariant (LTI) systems and nonlinear/time-varying plants as the signal relations in LPV systems are considered to be linear just as in the LTI case, but the parameters are assumed to be functions of a measurable time-varying signal, the socalled scheduling variable p : Z → P. Here the compact set P ⊂ R n P denotes the scheduling space. This LPV modeling concept allows for a wide representation capability of physical processes, but the real practical significance of the LPV framework lays in its well worked out and industrially reputed control synthesis approaches, e.g. [START_REF] Apkarian | A convex characterization of gainscheduled H∞ controllers[END_REF], [START_REF] Scherer | Mixed H 2 /H∞ control for time-varying and linear parametrically-varying systems[END_REF], [START_REF] Wu | Gain-scheduling control of LFT systems using parameter-dependent Lyapunov functions[END_REF], that have led to many successful applications of LPV control in practice [START_REF] Dettori | LPV design for a CD player: An experimental evaluation of performance[END_REF], [START_REF] Marcos | Development of linear-parameter-varying models for aircraft[END_REF], [START_REF] Østergaard | Rate bounded linear parameter varying control of a wind turbine in full load operation[END_REF], [START_REF] Wijnheijmer | Modeling and LPV control of an electro-hydraulic servo system[END_REF].

However a major drawback of the LPV framework today is that, despite the advances of the LPV control field, identification of such systems is not well developed as the current methods are unable to support practical control design. Commonly LPV controllers are synthesized in continuous time (CT) as stability and performance requirements of the closed loop behavior can be more conveniently expressed in CT, like in a mixed-sensitivity setting [START_REF] Zhou | Essentials of Robust Control[END_REF].

However, LPV identification methods are almost exclusively developed for discrete-time (DT) (for a recent survey see [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]), as in this setting it is much easier to handle the estimation of parameter-varying dynamics. Nonetheless, the absence of CT methods represents a gap between the available identification approaches and the needs of LPV control synthesis.

There is therefore a growing need of the LPV framework for efficient identification methods that directly deliver reliable CT models.

In practice, CT systems can only be identified based on sampled measured data records. Thus in general, for delivering a CT model estimate, the available approaches in system identification can be categorized as follows:

• Indirect approaches: These methods involve the identification of a DT model in a completely DT setting which is followed by the transformation of the DT model estimate into a CT form. • Direct approaches: The methods formulate the identification of the CT model directly based on samples of the measured CT signals. Unfortunately, transformation of DT-LPV models to CT-LPV models is more complicated than in the LTI case and despite recent advances in LPV discretization theory (see [START_REF] Tóth | On the discretization of LPV state-space representations[END_REF], [START_REF] Tóth | Discretization of linear fractional representations of LPV systems[END_REF]) the theory of CT realization of DT models is still in an immature state. The discretization of a CT LPV model results in a system order increase and more importantly in complicated dynamic dependency on p (dependency of the model coefficients on time-shifted versions of p with non-trivial rational functions) for which the available LPV identification methods are not well suited.

Even for a very simple CT-LPV model, estimation of a DT model with the purpose of obtaining afterwards a CT realization is a tedious task with many underlaying problems for which there are no general theoretical solutions available.

Unlike an indirect approach, a direct solution offers a way to efficiently overcome these problems but presents intrinsic difficulties mostly linked to the inaccessibility to signal timederivatives from the acquired sampled data. The offered solutions in the LTI case often require the use of signals prefiltering [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF]. These filters depend on hyperparameters input by the user and their efficiency is strongly linked to the adequacy with the considered system. In the CT-LTI case, one of the methods for relaxing the need of prefiltering is the Refined Instrumental Variable for Continuous-time (RIVC) method. This method is attractive in the sense that it provides consistent estimates under the realistic assumption of an unknown noise model and it achieves similar performance as prediction-error-minimization (PEM) methods [START_REF] Rao | Identification of continuous-time systems: direct or indirect ?[END_REF].

Another problem related to CT identification based on sampled data is the mathematical complexity of the CT random process used to describe the noise added to the system. An efficient way to overcome this problem is to consider a discrete-time Box-Jenkins noise model leading to hybrid models (see [START_REF] Pintelon | Box-Jenkins continuoustime modeling[END_REF], [START_REF] Johansson | Identification of Continuous-Time models[END_REF]). In order to avoid the different issues linked to the simultaneous use of CT and DT filters, this paper focuses on the case where the noise added onto the output is a white gaussian noise: the Output Error (OE) case.

Recently, an LPV identification approach has been introduced, which uniquely among other approaches, addresses the minimization of the prediction error for LPV-OE models by using a Multiple Input Single Ouput (MISO)-LTI reformulation of the data-generating LPV system [START_REF] Laurain | Refined instrumental variable methods for identification of LPV output-error and Box-Jenkins models[END_REF]. This paper aims at providing the very first step towards bridging the existing gap between LPV control and identification via the introduction of a direct CT identification approach that benefits from the properties of RIV methods in the continuous-time case and uses the recent advances of the prediction error minimization framework [START_REF] Laurain | Refined instrumental variable methods for identification of LPV output-error and Box-Jenkins models[END_REF], [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF].

The paper is organized as follows: in Section II, the general class of CT-LPV systems in an IO representation form is introduced. In Section III, the proposed CT LPV-IV method is described and analyzed, while its performance is illustrated in Section IV through a representative simulation example.

II. PROBLEM DESCRIPTION

A. System description

Consider the data generating CT LPV system described by the following equations

S o A o (p t , d )χ o (t) = B o (p t , d )u(t) y(t) = χ o (t) + e o (t), (1)
where d denotes the differentiation operator w.r.

A o (p t , d ) = d na + na i=1 a o i (p t )d na-i , (2)
B o (p t , d ) = n b j=0 b o j (p t )d n b -j . (3)
Note that a o i and b o j are functions of p at time t, which is called static dependence. In LPV system theory, a more general p-dependence of coefficients than static is required to establish equivalence of representations. In particular, it can be required that the coefficients a o i and b o j depend also on time derivatives of p, which is called dynamic dependence [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]. In order to simplify the upcoming discussion, we restrict our attention to static dependence. Nevertheless, the established results hold also in the case of dynamic dependence of (1) and of the proposed model structure.

In terms of identification we can assume that sampled measurements of (y, p, u) are available with a sampling period T s > 0. Hence, we will denote the discrete-time samples of these signals as u(t k ) = u(kT s ), where k ∈ Z.

B. Model structure considered

The process model is denoted by G ρ and defined in a form of an LPV-IO representation with a static scheduling dependence:

G ρ : (A(p t , d , ρ), B(p t , d , ρ)) (4)
where the p-dependent polynomials A and B given as

A(p t , d , ρ) = d na + na i=1 a i (p t )d na-i , B(p t , d , ρ) = n b j=0 b j (p t )d n b -j ,
are parameterized as

a i (p t ) = a i,0 + nα l=1 a i,l f l (p t ) i = 1, . . . , n a b j (p t ) = b j,0 + n β l=1 b j,l g l (p t ) j = 0, . . . , n b
In this parametrization, {f l } nα l=1 and {g l } n β l=1 are meromorphic functions of p, with static dependence, allowing the identifiability of the model (they can be chosen for example as linearly independent functions on P). The associated model parameters are stacked columnwise:

ρ = [ a 1 . . . a na b 0 . . . b n b ] ⊤ ∈ R nρ , (5)
where

a i = [ a i,0 a i,1 . . . a i,nα ] ∈ R nα+1 b j = [ b j,0 b j,1 . . . b j,n β ] ∈ R n β +1
and n ρ = n a (n α + 1) + (n b + 1)(n β + 1). Introduce also G = {G ρ | ρ ∈ R nρ }, as the collection of all process models in the form of [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF].

With respect to the considered OE structure, the signal relations of the LPV-BJ model, denoted in the sequel as M ρ , are defined as:

M ρ A(p k , d , ρ)χ(t) = B(p k , d , ρ)u(t) y(t k ) = χ(t k ) + e(t k ) (6)
Based on this model structure, the model set, denoted as

M = {Mρ | ρ ∈ R nρ } = G,

C. Predictors and prediction error

Similar to the LTI case, in the LPV prediction error framework, one is concerned about finding a model in a given LPV model structure M, which minimizes the statistical mean of the squared prediction error based on past samples of (y, u, p). However in the LPV case, no transfer function representation of systems is available. Furthermore, multiplication with d is not commutative over the p-dependent coefficients [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF], meaning that

d (B(p, d )u(t)) = B(d p, d )d u(t) which is not equal to B(p, d )d u(t).
1) System reformulation and prediction error: Following the same idea developed in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV output-error and Box-Jenkins models[END_REF] and if the system belongs to the model set defined with a deterministic p signal, it is possible to express the CT LPV system as a CT MISO LTI system by rewriting the signal relations of (1) as

χ (na) o (t) + na i=1 a o i,0 χ (na-i) o (t) Fo(d )χo(t) + na i=1 nα l=1 a o i,l f l (p(t))χ (na-i) o (t) χ o i,l (t) = n b j=0 n β l=0 b o j,l g l (p(t))u (n b -j) (t) u j,l (t) (7)
where g 0 (t) = 1 and the superscript (n) for a signal, like u (n) , denotes the n th time-derivative of the signal, e.g.

u (n) (t) = d n u(t). Furthermore, F o (d ) = d na + na i=1 a i,0 d na-i while u (n) (t k ) represents the value of the signal u (n) (t) sampled at time instance t k .
Note that in this way, the time variation of the coefficients is transposed onto the signals χ o i,l (t) and u j,l (t):

χ o i,l (t) = f l (p(t))χ (na-i) o (t) {i, l} ∈ {1 . . . n a , 1 . . . n α }, u j,l (t) = g l (p(t))u (n b -j) (t) {j, l} ∈ {1 . . . n b , 1 . . . n α }.
Therefore, the process part of the LPV-BJ model is rewritten as a Multiple-Input Single-Output (MISO) system with (n b + 1)(n β + 1) + n a n α inputs {χ o i,l } na,nα i=1,l=1 and {u j,l } n b ,n β j=0,l=0 . By using ( 7), ( 6) can be rewritten in terms of the sampled output signal y(t k ) as

y(t k ) = - na i=1 nα l=1 a o i,l F o (d ) χ o i,l (t k ) +   n b j=0 n β l=0 b o j,l F o (d ) u k,j   (t k ) + e o (t k ), (9)
which is a sampled LTI representation of the system (1).

2) Prediction Error Model: Similarly to the LTI case, the one-step-ahead prediction error can be expressed and defined as [START_REF] Ljung | System identification : theory for the user -Second Edition[END_REF]:

ε ρ (t k ) = y(t k ) -ŷρ (t k ), (10)
where ŷρ (t k ) is the one step ahead predictor based on the model ( 6) written as a MISO LTI form (9) and defined as:

ŷρ (t k ) = - na i=1 nα l=1 a i,l F (d , ρ) χ i,l (t k ) +   n b j=0 n β l=0 b j,l F (d , ρ) u k,j   (t k ). ( 11
)
3) Prediction error minimization:

Denote D N = {y(t k ), u(t k ), p(t k )} N
k=1 a data sequence of S o . Then to provide an estimate of ρ based on the minimization of ε ρ , an identification criterion W (D N , ρ) can be introduced, like the least square criterion

W (D N , ρ) = 1 N N k=1 ε 2 ρ (t k ), (12)
such that the parameter estimate is ρN = arg min ρ∈R nρ+nη W (D N , ρ).

(13)

4) CT filtering and sampled data:

The CT representation of the model ( 6) consists in a CT filtering operation. In this paper, we considered the practically feasible situation such that only sampled measurements of the CT signals (y, p, u) are available. In order to apply a CT filter on sampled data one can either interpolate the samples to obtain a continuoustime signal and apply the CT filter on this reconstructed signal or use a numerical approximation, i.e. DT approximation of the considered system. This is a common problem for simulation of CT systems. For simulation purposes, DT approximation of the system can efficiently be dealt with, by using powerful numerical algorithms available [START_REF] Atkinson | An Introduction to Numerical Analysis[END_REF].

Note that to derive an accurate DT approximation of the system, it is often sufficient in terms of the classical discretization theory to assume that the sampled free CT signals of the system are restricted to be constant in the sampling period [START_REF] Goodwin | Control System Design[END_REF]. This has been also shown in case of LPV systems with static dependence [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]. This provides the hypothesis, also used in [START_REF] Pintelon | Box-Jenkins continuoustime modeling[END_REF], [START_REF] Johansson | Identification of Continuous-Time models[END_REF], that if CT (p, u) are piecewise constant between two samples, then the trajectory of y is completely determined by its observations at the sample period T s k. Therefore, under these inter-sampling conditions, DT filtering and numerical approximation of CT filtering operations commute [START_REF]Identification of Continuous-time Models from Sampled Data[END_REF]. Nevertheless, it is important to notice that the numerical approximation method used for the evaluation of a CT filter does not have any impact on the coefficients to be estimated which remain, in terms of [START_REF] Ljung | Initialisation aspects for subspace and output-error identification methods[END_REF], the coefficients of the parsimonious CT model.

D. Identification problem statement

Based on the previous considerations, the identification problem addressed in the sequel can now be defined.

Problem 1: Given a CT-LPV data generating system S o defined as in [START_REF] Apkarian | A convex characterization of gainscheduled H∞ controllers[END_REF] Based on the MISO-LTI formulation [START_REF] Ljung | Initialisation aspects for subspace and output-error identification methods[END_REF], it becomes possible in theory to achieve optimal PEM using linear regression [START_REF] Laurain | Refined instrumental variable methods for identification of LPV output-error and Box-Jenkins models[END_REF]. This allows to extend the Refined Instrumental Variable (RIV) approach of the LTI identification framework to provide an efficient way of identifying CT LPV models.

A. Linear Regression for CT LPV-BJ models

Using the LTI model ( 6) reformulated as in [START_REF] Ljung | Initialisation aspects for subspace and output-error identification methods[END_REF], y(t k ) can be written in the regression form:

y (na) (t k ) = ϕ ⊤ (t k )ρ + ṽ(t k ) (14) where, ϕ(t k ) = [-y (na-1) (t k ) . . . -y(t k ) -χ 1,1 (t k ) . . . . . . -χ na,nα (t k ) u 0,0 (t k ) . . . u n b ,n β (t k )] ⊤ ρ = [a 1,0 . . . a na,0 a 1,1 . . . a na,nα b 0,0 . . . b n b ,n β ] ⊤ ṽ(t k ) = F (d , ρ)e(t k ).
The extended regressor in ( 14) contains the noise-free output terms {χ i,k }. Therefore, by momentary assuming that {χ i,l (t k )} na,nα i=1,l=0 are known a priori, the prediction error ε ρ (t k ) for ( 14) is given in terms of [START_REF] Ljung | System identification : theory for the user -Second Edition[END_REF] as:

ε ρ (t k ) = (F (d , ρ)y f ) (t k ) - na i=1 nα l=1 a i,l χ f i,l (t k ) + n b j=0 n β l=0 b j,l u f k,j (t k ) (15)
where y f (t k ), u f j,l (t k ) and χ f i,l (t k ) represent the outputs of an hybrid prefiltering operation, involving the continuous-time filter (see [START_REF] Young | Refined instrumental variable identification of continous-time hybrid Box-Jenkins models[END_REF]):

Q c (d , ρ) = 1 F (d , ρ) , (16)
Based on [START_REF] Pintelon | Box-Jenkins continuoustime modeling[END_REF], the associated linear-in-the-parameters model takes the form [START_REF] Young | Refined instrumental variable identification of continous-time hybrid Box-Jenkins models[END_REF]:

y (na) f (t k ) = ϕ ⊤ f (t k )ρ + ṽf (t k ), (17)
where

ϕ f (t k ) = [-y (na-1) f (t k ) . . . -y f (t k ) -χ f 1,1 (t k ) . . . . . . -χ f na,nα (t k ) u f 0,0 (t k ) . . . u f n b ,n β (t k )] ⊤ ṽf (t k ) = Q c (d , ρ)ṽ(t k ) = e(t k ).

B. The refined instrumental variable approach

Under the assumption that the CT filter Q c (d , ρ) and {χ i,l (t k )} na,nα i=1,l=0 are known a priori, traditional parametric estimation methods from the LTI framework could provide efficient estimates of ρ. However, in a practical situation, Q c (d , ρ) is unknown and only some estimates will be available.

Furthermore, it is important to notice here that the regressors in [START_REF] Rao | Identification of continuous-time systems: direct or indirect ?[END_REF] and ( 14) contain some time-derivatives of y and u which, in the assumed framework considering sampled data, can only be approximated. It is well-known that the approximation of derivatives requires a low pass filtering on the input and output. The most commonly used filters for this purpose are Poisson's filters, or state-variable filters [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF]. The drawback of these filters is that they require the choice of a design variable. Therefore, it is a particular strength of the presented method that the estimated filter F (d , ρ) is not only used for the minimisation of the prediction error but it also provides the filtering for the approximation of the time derivatives. In other words, the regressor ϕ f in ( 17) can be well-approximated numerically whereas the regressor ϕ from ( 14) cannot as it requires prefiltering of the data which must be chosen by the user. In order to estimate the parameter vector in [START_REF] Rao | Identification of continuous-time systems: direct or indirect ?[END_REF] without the prior knowledge of Q c (d , ρ), the RIV method is chosen for the following reasons:

• RIV methods lead to optimal estimates in the LTI case if S o ∈ M (see [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | The refined instrumental variable method: Unified estimation of discrete and continuous-time transfer function models[END_REF], [START_REF] Young | Refined instrumental variable identification of continous-time hybrid Box-Jenkins models[END_REF]). • In a practical situation of identification, G o ∈ G might be fulfilled due to first principle or expert's knowledge. However, it is commonly fair to assume that the model is not OE. In such case, RIV methods have the advantage of providing consistent estimates whereas methods such as extended LS are biased and more advanced PEM methods need robust initialisation [START_REF] Ljung | Experiments with identification of continuous time models[END_REF]. • The RIV algorithm has been successfully used for LTI model with similar CT structure, in the case of linear models [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF], [START_REF] Young | The refined instrumental variable method: Unified estimation of discrete and continuous-time transfer function models[END_REF] and nonlinear ones [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous time Box-Jenkins models[END_REF]. Aiming at the extension of the RIV approach for the estimation of CT LPV models, consider the relationship between the process input and output signals as in [START_REF] Østergaard | Rate bounded linear parameter varying control of a wind turbine in full load operation[END_REF]. Based on this form, the extended-IV estimate is given as [START_REF] Young | The refined instrumental variable method: Unified estimation of discrete and continuous-time transfer function models[END_REF]:

ρXIV (N ) = arg min ρ∈R nρ 1 N N k=1 ζ f (t k )ϕ ⊤ f (t k ) ρ - 1 N N t=1 ζ f (t k )y (na) f (t k ) 2 W , (18)
where ζ(t k ) is the instrument, x 2 W = x T W x, with W a positive definite weighting matrix and the filtered variables ζ f , ϕ f and y f are filtered using a stable prefilter. If G o ∈ G, the extended-IV estimate is consistent under the following two conditions 2

: C1 Ē{ζ f (t k )ϕ ⊤ f (t k )} is full column rank. C2 Ē{ζ f (t k )ṽ f (t k )} = 0.
Moreover it has been shown in [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF] and [START_REF] Young | The refined instrumental variable method: Unified estimation of discrete and continuous-time transfer function models[END_REF] that the minimum variance estimator can be achieved if:

C3 W = I. C4
ζ is chosen as the noise-free version of the extended regressor in [START_REF] Østergaard | Rate bounded linear parameter varying control of a wind turbine in full load operation[END_REF] and is therefore defined in the present LPV case as:

ζ(t k )= -χ (na-1) (t k ) . . . -χ(t k ) -χ 1,1 (t k ) . . . . . . χ na,nα (t k ) u 0,0 (t k ) . . . u n b ,n β (t k ) ⊤ 2
The notation Ē{.} = lim N →∞ 1 N P N t=1 E{.} is adopted from the prediction error framework of [START_REF] Ljung | System identification : theory for the user -Second Edition[END_REF].

C5

G o ∈ G and n ρ is equal to the minimal number of parameters required to represent G o with the considered model structure.

C6

The CT filter used is chosen as the filter [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF]. While conditions C1, C2, C3 and C5 are quite straightforward to fulfill (see [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF]), the obtention of the suitable instrument that fulfills C4 and of the optimal filter fulfilling C6 is not trivial in practical case. The RIV algorithm involves an iterative (or relaxation) algorithm in which, at each iteration, an 'auxiliary model' is used to generate the instrumental variables (which guarantees C2), as well as the associated prefilters. This auxiliary model is based on the parameter estimates obtained at the previous iteration. Consequently, if convergence occurs, C4 and C6 are fulfilled. Thus, the RIV is the most suitable method to simultaneously i) efficiently estimate the parameter vector ρ in [START_REF] Rao | Identification of continuous-time systems: direct or indirect ?[END_REF] , where F (d , ρ(τ) ) is as given in [START_REF] Johansson | Identification of Continuous-Time models[END_REF].

Step 4 Use the CT filter Qc (d , ρ(τ) ) as well as χ(t k ) in order to generate the estimates of the derivatives which are needed: Step 5 Build the filtered estimated regressor φf (t k ) and, in terms of C4, the filtered instrument ζf (t k ).

Step 6 The solution of the IV optimization problem is then

ρ(τ+1) (N ) = N k=1 ζf (t k ) φ⊤ f (t k ) -1 N k=1 ζf (t k )y (na) f (t k ) (19)
Step 7 If ρ (τ +1) has converged or the maximum number of iterations is reached, then stop, else increase τ by 1 and go to Step 2.

IV. SIMULATION EXAMPLE In order to show the relevance of the presented algorithm, the following data-generating system is considered:

S o A o (d , p) = d 2 + a o 1 (p)d + a o 2 (p) B o (d , p) = b o 0 (p)d + b o 1 (p)
where

a o 1 (p) = 1 -0.5p, (20a) a o 2 (p) = 5 + 3p (20b) b o 0 (p) = 2 + p, (20c) b o 1 (p) = 5 -p. ( 20d
)
3000 samples are collected from a 15sec simulation (T s = 0, 005sec). The input signal u is chosen as a uniformly distributed sequence U(-1, 1) while the scheduling variable is chosen as p(t) = sin( 2 3 πt). The following model structure in terms of ( 4) is considered to capture the dynamics of S o :

M A(d , p) = d 2 + a 1 (p)d + a 2 (p) B(d , p) = b 0 (p)d + b 1 (p)
where

a 1 (p) = a 1,0 + a 1,1 p, (21a) a 2 (p) = a 2,0 + a 2,1 p, (21b) b 0 (p) = b 0,0 + b 0,1 p, (21c) b 1 (p) = b 1,0 + b 1,1 p. (21d)
As previously pointed out, the efficiency of the LPV-RIVC estimator can not be yet proven. Nonetheless, in order to analyse its statistical properties on this example the model is estimated using both the LPV-RIVC algorithm and the MATLAB LSQNONLIN method. The LSQNONLIN method is a nonlinear statistically optimal optimization method but such nonlinear method is also sensitive to initialisation [START_REF] Ljung | Initialisation aspects for subspace and output-error identification methods[END_REF]. Therefore, in order to put this latter method at its best, it is initialized on the true parameters.

Monte Carlo simulation results obtained using the different methods are presented in Table I. The statistical properties of each method is evaluated using the estimated parameters mean and standard deviation. These results are based on N run = 100 random realizations under a Signalto-Noise Ratio (SNR) of 10dB with: SNR = 10 log

P χo P eo , (22)
where P x is the power of signal x.

It can be seen from Table I that according to the theory, the estimated parameters using the LPV-RIVC algorithm are unbiased. Moreover, in case of correct parametrization and on this example, the LPV-RIVC method performs equivalently to the optimal LSQNONLIN method which is theoretically optimal but is also known to be sensitive to initialisation in comparison to RIV based algorithms [START_REF] Ljung | Initialisation aspects for subspace and output-error identification methods[END_REF].

Consequently, the presented algorithm constitutes the first direct continuous-time method aiming at minimizing the error prediction. It looks from this example that the empirically accepted properties of RIV based method might apply to the LPV case even though this cannot be yet proven. Moreover, this method does not requires any hyperparameters from the user for the signals time-derivative approximation.

V. CONCLUSION

The proposed approach provides on of the very first direct global LPV identification method that is able to give consistent estimates of LPV-IO models in continuous-time and has a low computational load. The proposed algorithm has been tested on a representative numerical simulation example and it has been shown that the proposed procedure is robust to noise and can compete with the optimal nonlinear optimization method even in the case where the latter is initialized knowing the true parameters. Furthermore, based on previous work on CT-LTI systems operating in closed loop, this methods opens the possibility for closed-loop CT-LPV identification.

  corresponds to the set of candidate models in which we seek the model that explains data gathered from S o the best, under a given identification criterion (cost function).

  and a data set D N collected from S o . Based on the CT LPV model structure M ρ defined by (6), estimate the parameter vector ρ using D N under the

	following assumptions:
	A1 A2	S o ∈ M. In the parametrization A ρ and B ρ , {f l } nα l=1 and {g l } n β l=1 are chosen such that (G o ) is identifiable
		for any trajectory of p.
	A3	u(t k ) is not correlated to e o (t k ).
	A4	D N is informative with respect to M.
	A5	S

o is globally BIBO stable, i.e. for any trajectory of p : R → P and any bounded input signal u, the output of S o is bounded

[START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]

.

III. REFINED INSTRUMENTAL VARIABLE FOR LPV SYSTEMS

  in the case S o ∈ M and ii) consistently estimate ρ in the practical assumption where the noise model is misspecified.C. The LPV-RIVC AlgorithmStep 1 The usual initialisation for CT-RIV algorithm is a DT model estimate issued from an LS method or a DT-RIV algorithm. In the LPV case however, the transformation of a DT model into a CT model is not trivial. Consequently, the initial estimate proposed for the LPV-RIVC algorithm is an LTI-RIVC estimate of M ρ , i.e. ρ(0) is given.

	Set τ = 0.
	Step 2 Compute an estimate of χ(t k ) via numerical ap-
	proximation of
	A(p t , d , ρ(τ) ) χ(t) = B(p t , d , ρ(τ) )u(t),
	where ρ(τ) is estimated in the previous iteration.
	Based on M ρ(τ) , deduce χ(t k ) which is bounded
	according to Assumption A5.
	Step 3 Compute the estimated continuous-time filter 1 Qc (d , ρ(τ) ) = F (d , ρ(τ) )

A function f is called meromorphic if f = g h where g, h are holomorphic (analytic) functions and h is not the zero function.


David Folio

Viviane Cadenat

A sensor-based controller able to treat total image loss and to guarantee non-collision during a vision-based navigation task

This paper deal with the problem of executing a vision-based task in an unknown environment. During such a task, two unexpected events may occur: the image data loss due to a camera occlusion and the robot collision with obstacles. We first propose a method allowing to compute the visual data when they are totally lost, before addressing the obstacle avoidance problem. Then, we design a sensor-based control strategy to perform safely vision-based tasks despite complete loss of the image. Simulation and experimental results validate our work.

I. INTRODUCTION

Visual servoing techniques aim at controlling the robot motion using visual features provided by a camera [START_REF] Chaumette | Visual servo control, part I: Basic approaches[END_REF]. Therefore, these techniques cannot be used anymore if these features are lost during the robotic task. A first way to deal with this kind of problem is to use methods allowing to preserve the visual data visibility [START_REF] Mezouar | Avoiding self-occlusions and preserving visibility by path planning in the image[END_REF][3][4] [START_REF] Remazeilles | Qualitative visual servoing: application to the visibility constraint[END_REF]. However, these works are mostly applied to manipulators. The problem is slightly different in mobile robotics. Indeed, when a robot navigates in an unknown environment using a camera, it is necessary to guarantee both image features visibility and robot safety. Some works in this area propose to fulfill simultaneously these two objectives [START_REF] Folio | A controller to avoid both occlusions and obstacles during a vision-based navigation task in a cluttered environment[END_REF] [START_REF]Using redundancy to avoid simultaneously occlusions and collisions while performing a vision-based task amidst obstacles[END_REF]. However, these methods appear to be limited because they do not allow to perform robotic tasks whose success depends on the ability of tolerating the visual data loss.

In this paper, we consider the problem of executing a vision-based task in an unknown environment. During such a task, obstacles may occlude totally the camera or induce a risk of collision. Our first goal is to build a method allowing to compute the image data when they are entirely unavailable. Then, we address the obstacle avoidance problem using potential field techniques. Finally, thanks to these results, we design a multi-sensor-based control law able to tolerate the image total loss and guarantee non collision.

The article is organized as follows. Section II is dedicated to the problem statement and to our estimation method. Section III addresses the avoidance problem and describes the chosen control strategy. Finally, last section presents our robotic system and some successful experimental results.

II. VISUAL DATA ESTIMATION

In this section, we focus on the visual data estimation problem. Obstacle avoidance will be treated in next section.

Authors are with LAAS-CNRS; Université de Toulouse; 7, avenue du Colonel Roche, F-31077 Toulouse, France. and Université de Toulouse; UPS; Email: {dfolio,cadenat}@laas.fr A. Modelling and Preliminaries 1) System modelling: Let us consider figure 1(a). We define the successive frames: F M (M, -→ x M , -→ y M , -→ z M ) linked to the robot, F P (P, -→ x P , -→ y P , -→ z P ) attached to the panplatform, and F C (C, -→ x c , -→ y c , -→ z c ) linked to the camera. Let ϑ be the direction of the pan-platform wrt. -→ x M , P the pan-platform center of rotation and D x the distance between the robot reference point M and P . The control input is defined by: q = (v, ω, ̟) T , where v and ω are the cart linear and angular velocities, and ̟ is the pan-platform angular velocity wrt. F M . For this specific mechanical system, the kinematic screw v c is related to the control input by the robot jacobian J: v c = J q. As the camera is constrained to move horizontally, it is sufficient to consider a reduced kinematic screw

v c r = (V-→ y C , V-→ zC , Ω-→ xC )
T , and a reduced jacobian matrix J r as follows:

v c r =   -sin(ϑ) Dx cos(ϑ) + Cx Cx cos(ϑ) Dx sin(ϑ) -Cy -Cy 0 -1 -1     v ω ̟   = Jr q (1)
where C x and C y are the coordinates of C along axes -→ x P and -→ y P (see figure 1(a)). 2) Preliminaries: Our goal is to realize a vision-based navigation task with respect to a given static landmark despite possible total occlusions. We assume that this landmark can be characterized by n interest points which can be extracted by our image processing. Therefore, the visual data are represented by a 2n-dimensional vector s made of the coordinates (X, Y ) of each projected point P in the image plane as shown on figure 1(b). For a fixed landmark, the variation of the visual signals ṡ is related to the reduced camera kinematic screw v c r thanks to the interaction matrix L (s) as shown below [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]:

ṡ = L (s) v c r = L (s) J r q (2)
In the case of n points,

L (s) = [L (P 1 ) T , . . . , L (Pn)
T ] T where L (P ) is classically given by [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]:

L (P ) = L (X,z) L (Y,z) = 0 X z XY f -f z Y z f + Y 2 f (3)
where z represents the depth of each projected point P in the image plane, and f is the camera focal (see figure 1(b)).

B. Motivations and hypotheses

Now, we focus on the problem of estimating image data whenever they are lost during the vision-based task. Thus, the key-assumption of this work is that the whole image is considered to be momentarily totally unavailable. Therefore, methods dedicated to partial occlusions such as [START_REF] Comport | Robust model-based tracking for robot vision[END_REF][10] cannot be used here. It is also possible to consider other approaches based on signal processing techniques or tracking methods such as [11][12]. . . However, most of this kind of techniques rely on measures from the image which is considered to be totally unavailable in our case. Moreover, recalling that our estimation method will be used in the control law design, it must provide the estimated visual signals sufficiently rapidly with respect to the control sampling period. Therefore, we need a simple and fast method. We have then thought of using a 3D model of the landmark together with projective geometry. However, this method requires to determine accurately the calibration parameters to reconstruct successfully the lacking data. Thus, we propose another approach consisting in solving (2) using the previous visual measurements and the control inputs q. Therefore, we have made the two following hypotheses: first, we assume that the image features can be measured at the beginning of the robotic task. Second, we suppose that few errors occur on the model and on the measures (we will see in section IV that this hypothesis is reliable).

C. Problem Statement

Our idea is to solve the dynamic system (2) to determine the image data s. However, this last system depends on depth z which must then be determined. As our robot is not equipped with any sensor able to measure this data (see section IV), we have to reconstruct it. To this aim, different approaches may be used [START_REF] Matthies | Kalman filter-based algorithms for estimating depth in image sequences[END_REF] [START_REF] Jerian | Structure from motion: A critical analysis of methods[END_REF][15] [START_REF] Oliensis | Exact two-image structure from controlled motion[END_REF]. However, most of these techniques rely on measures from the image and cannot be applied here. In this paper, we propose another solution consisting in simultaneously estimating depth z and image data s. It is then necessary to express the analytical relation between the variation of z and the camera motion. As the considered landmark is made of n points, it can be easily shown that, for one 3D point p projected into a point P (X, Y ) in the image plane as shown in figure 1(b), the depth variation ż is related to the camera motion according to: ż = L (z) J r q, where L (z) = (0, -1, z f Y ). Thus, the dynamic system to be solved for one point P (X, Y ) expresses as follows:

     Ẋ = X z V-→ z C + XY f Ω-→ x C Ẏ = -f z V-→ y C + Y z V-→ zC + f + Y 2 f Ω-→ xC ż = -V-→ zC -z f Y Ω-→ xC (4)
Now, we address the resolution problem. A first idea is to integrate the above system for any t ∈ [t 0 , t f ] where t 0 and t f are respectively the initial and final instants of the task. However, the computations are quite difficult to carry out. Therefore, in this work, we propose to discretize the problem and to solve system (4) during a time control interval [t k , t k+1 ]. Defining ψ= [X1, Y1, . . . , Xn, Yn, z1, . . . , zn] T and

L (ψ) = L (s) T , L (z)
T T , the differential system to be solved for n points is given by:

ψ(t) = L (ψ k ) J r (ϑ(t)) q(t k ) ψ(t k ) = ψ k = s T k , z T k T , ∀t ∈ [t k , t k+1 ] (5)
where

q(t k )= (v k , ω k , ̟ k ) T is the control input at t k .
Note that, on our robot, the control input is hold during the control law sampling period T s . We will then consider in the sequel that q(t k ) remains constant during this time interval.

D. Our estimation method

A first idea is to solve differential system (5) using numerical schemes as done in [START_REF] Folio | A new controller to perform safe visionbased navigation tasks amidst possibly occluding obstacles[END_REF] [START_REF]Using simple numerical schemes to compute visual features whenever unavailable[END_REF]. The main advantage of such approaches lies in their generality. Indeed, they can be applied to any set of visual features and to any kind of robots, provided that the estimation problem can be expressed as system [START_REF] Remazeilles | Qualitative visual servoing: application to the visibility constraint[END_REF]. However, with such methods, the accuracy of the obtained solution depends on the chosen algorithm and on the integration step size. This is the reason why we propose in this work to solve analytically system [START_REF] Remazeilles | Qualitative visual servoing: application to the visibility constraint[END_REF]. In this way, we obtain its exact solution on [t k , t k+1 ] and the accuracy of the result is significantly improved. After some computations (see [START_REF] Folio | Stratégies de commande référencées multi-capteurs et gestion de la perte du signal visuel pour la navigation d'un robot mobile[END_REF] for a detailed proof), it can be shown that, for any t ∈ [t k , t k+1 ], the evolution of abscissa X and ordinate Y of one point P in the image expresses as:

         X(t)= z k X k z(t) Y (t) =f ż+v k cos(ϑ(t))+Dxω k sin(ϑ(t))-Cy(ω k +̟ k ) (ω k +̟ k )z(t) z(t) =c 1 sin (A 1 (t -t k )) + c 2 cos (A 1 (t -t k )) -D x cos (ϑ(t)) + v k ω k sin (ϑ(t)) -C x (6)
where:

   A1 = (ω k + ̟ k ) c1 = -Y k z k f + Dx sin (ϑ k ) + v k ω k cos (ϑ k ) -Cy c2 = z k + Dx cos (ϑ k ) -v k ω k sin (ϑ k ) + Cx
As one can see, the solution (6) require the determination of ϑ(t). This angle can simply be computed by integrating θ = ̟ between t k and t. Some straightforward calculus leads to ϑ(t) = ̟ k (tt k ) + ϑ k , where ϑ k is the panplatform angular value at t k , which is usually provided by the embedded encoder. The solution for the set of the n points of the landmark is then easily obtained by applying the above solution (6) on each component of the vector ψ.

III. CONTROL STRATEGY

Let us recall that our goal is to perform vision-based navigation tasks in environments cluttered with occluding obstacles, which requires to deal with occlusions and collisions. The first problem will be treated thanks to our estimation method, and the second one using a potential field technique. So, first of all, we will briefly describe how to detect the risks of occlusion and collision. Then, we will detail our control strategy, considering two cases when the occlusion occur in the free space or in the obstacle neighborhood. Finally, we will end this part by some simulation results. In this part, our goal is to detect the risks of collision and occlusion. To this aim, we design two parameters µ coll and µ occ which smoothly increase between 0 (no risk) and 1 (high risk). From the embedded ultrasonic sensors we can compute the shortest distance d coll and the relative orientation α between the robot and the obstacle (see figure 2(a)). Defining by ξ + , ξ 0 , ξ -three envelopes surrounding each obstacle at distance d + > d 0 > d -, it is possible to choose µ coll (d coll ) so that it smoothly increases from 0 when the robot is far from the obstacle (d coll > d0) to 1 when it is close to it (d coll < d-). Assuming that the image processing is able to identify any occluding object O, µ occ can be determined by detecting their left and right borders. Indeed, from these data, we can deduce the shortest distance d occ between s and O (see figure 2(b)). Introducing three envelopes Ξ + , Ξ 0 , Ξ - around the occluding object located at D + > D 0 > D - from it, µ occ (d occ ) can be chosen to smoothly increase from 0 when O is far from the visual features (docc > D0) to 1 when it is close to them (docc < D-). Further details about these modelling, and different µ coll and µ occ expression can be found in [START_REF] Folio | A controller to avoid both occlusions and obstacles during a vision-based navigation task in a cluttered environment[END_REF] and [START_REF]Using redundancy to avoid simultaneously occlusions and collisions while performing a vision-based task amidst obstacles[END_REF].

A. Collision and Occlusion Detection

B. Control in the free space

Our goal is to position the embedded camera with respect to a landmark made of n points despite possible occlusions. To this aim, we have applied the visual servoing technique given in [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] to mobile robots as in [START_REF] Pissard-Gibollet | Applying visual servoing techniques to control a mobile hand-eye system[END_REF]. The proposed approach relies on the task function formalism [START_REF] Samson | Robot Control : The task function approach[END_REF] and consists in expressing the visual servoing task by the following task function to be regulated to zero:

e VS = L (s ⋆ ) + (s -s ⋆ ) (7)
where s ⋆ represents the desired value of the image data. Classically, a controller making e VS vanish can be designed by imposing an exponential decrease, that is: ėVS = -λ VS e VS , where λVS is a positive scalar or a positive definite matrix. Then, the visual servoing controller can be written as follows:

q(s) = J -1 r (-λ VS )L (s ⋆ ) + (s -s ⋆ ) (8)
Now, we suppose that an obstacle enters the camera field of view without inducing a danger of occlusion. In such a case, this controller cannot be computed anymore, leading to a task failure. Here, our idea is to embed our estimation technique in control law [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] to treat successfully this unexpected event. We simply propose to replace in expression [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] the real visual features s by their values s computed by our method. In case of occlusions, the visual servoing controller is then given by: q( s) = J -1 r (-λ VS )L (s ⋆ ) + ( ss ⋆ ). Now, it remains to switch between q(s) and q( s) depending on the landmark visibility. We propose the following final visual servoing controller:

qVS = (1 -µ occ ) q(s) + µ occ q( s) (9)
Remark 1: is no need to smooth controller (9) when the features are lost. Indeed, when the occlusion occurs, as the last provided data are used to feed our reconstruction algorithm, the values of s and s are close and so are q(s) and q( s) . Usually, the same reasoning holds when the occluding object leaves the image. However, some smoothing may be useful if the camera motion has been unexpectedly perturbed during the estimation phase or if the algorithm has been given too inaccurate initial conditions.

C. Control in the obstacle neighborhood

Now, we address the problem of controlling the robot in the obstacle neighborhood. This is the most complex part of our work because, during this phase, both risks of collision and occlusion may occur simultaneously (or not!). Then, we propose to design a multi-sensor-based controller able to guarantee the robot safety and to treat total image loss if needed. To this aim, we propose to extend the controller elaborated in [START_REF] Cadenat | A controller to perform a visually guided tracking task in a cluttered environment[END_REF] which only allows to avoid obstacles while tracking the landmark. We first recall the chosen control strategy and then explain how to embed our estimation technique to deal with the occlusions.

1) Obstacle avoidance strategy:

The chosen strategy consists in defining around each obstacle a rotative potential field so that the repulsive force is orthogonal to the obstacle when the robot is close to it (d coll < d -), parallel to the obstacle when the vehicle is at a distance d 0 from it, and progressively directed towards the obstacle between d 0 and d + (see figure 2(a)). The interest of such a potential is that it can make the robot move around the obstacle without requiring any attractive force, reducing local minima problems. We propose the following potential function U (d coll ) [START_REF] Cadenat | A controller to perform a visually guided tracking task in a cluttered environment[END_REF]:

U (d coll ) = 1 2 k 1 1 d coll -1 d + 2 + 1 2 k 2 (d coll -d + ) 2 if d coll ≤ d + U (d coll ) =0 otherwise
where k 1 and k 2 are positive gains to be chosen. The virtual repulsive force is defined to generate the desired rotative potential field around the obstacle by the couple (F, β), where F = -∂U ∂d coll is the modulus of the virtual repulsive force and β = απ 2d0 d coll + π 2 its direction wrt. F M . The mobile base velocities v coll and ω coll are then given by [START_REF] Cadenat | A controller to perform a visually guided tracking task in a cluttered environment[END_REF]:

qmb = v coll ω coll T = kvF cos β kω Dx F sin β T ( 10
)
where k v and k ω are positive gains to be chosen. qmb drives only the mobile base in the obstacle neighborhood. However, if the pan-platform remains uncontrolled, it will be impossible to switch back to the vision-based task at the end of the avoidance phase. Therefore, our idea is to control the pan-platform to keep the landmark in the camera line of sight during the avoidance motion. To this aim, as the camera is constrained to move within an horizontal plane, it is sufficient to regulate to zero the error e gc = Y gc -Y ⋆ gc where Y gc and Y ⋆ gc are the current and desired ordinates of the target gravity center. Rewriting equation ( 1) as: v c r = J mb qmb + J ̟ ̟ coll , and imposing an exponential decrease to regulate e gc to zero ( ėgc = -λgcegc, λgc > 0), we finally obtain [START_REF] Cadenat | A controller to perform a visually guided tracking task in a cluttered environment[END_REF]:

̟ coll = -1 L (Ygc) J̟ (λgcegc + L (Ygc) J mb qmb ) (11)
where L (Ygc) is the 2 nd row of L (P ) evaluated for Y gc .

2) Extension to deal with occlusions during avoidance: In the obstacle vicinity, the visual data may enter and leave the image plane depending on the relative robot-obstacle configuration and on the obstacle height. We may even imagine that the camera could be occluded by another obstacle than the one which is currently avoided. We need then to design a control law able to deal with these sudden and unexpected modifications. Our idea is then to use our estimation method to provide the necessary visual data to the pan-platform control law whenever needed. Thus, when an occlusion occurs, our algorithm is fed with the last available measures and computes an estimation s of the "real" visual features s from which we can deduce Y gc . The pan-platform controller during an occlusion phase will then be determined by replacing Y gc by its estimate Y gc in [START_REF] Favaro | Seeing beyond occlusions (and other marvels of a finite lens aperture)[END_REF]. We get:

̟ coll = -1 L ( Ygc) J̟ (λgc egc + L ( Ygc) J mb qmb ), (12)
where e gc = Y gc -Y ⋆ gc , and L ( Ygc) is deduced from (3). Now, it remains to select the suitable controller depending on the availability of the image. Recalling that µ occ allows to detect occlusions, we propose the following avoidance controller:

qcoll = v coll , ω coll , (1 -µocc)̟ coll + µocc ̟ coll T (13)
Thanks to this controller, it is possible to avoid the obstacle while treating the image total loss (if needed).

D. Global control strategy

So, we have previously designed two controllers given by ( 9) and ( 13). The first one allows to execute a vision-based task despite occlusions; the second one guarantees obstacle avoidance while treating image loss if necessary. To execute our vision-based navigation task, we have to apply to the robot the right controller depending on the risks of occlusions and collisions. Several approaches allow to sequence tasks [START_REF] Pissard-Gibollet | Applying visual servoing techniques to control a mobile hand-eye system[END_REF][23] [START_REF] Mansard | Tasks sequencing for visual servoing[END_REF]. Here, we have chosen a method relying on convex combinations between the successive controllers [START_REF] Cadenat | A controller to perform a visually guided tracking task in a cluttered environment[END_REF], as it allows to carry out applications more easily. We propose the following global control law:

q = (1 -µ coll ) qVS + µ coll qcoll (14)
where qVS and qcoll = (v coll , ω coll , ̟ coll ) T are respectively given by equations ( 9) and ( 13).

Remark 2: µ coll and µocc are defined to be maintained to 1 once they have reached this value [START_REF] Folio | A controller to avoid both occlusions and obstacles during a vision-based navigation task in a cluttered environment[END_REF]. Moreover, the different envelops are chosen close enough to reduce the transition phase duration. Thus the control strategy is built to insure that the robot will be rapidly controlled by the most relevant controller. In this way, the risks of instability, target loss or collisions during the switch are significantly reduced and the task feasibility can be considered to guaranteed.

E. Simulation results

To validate our work, we have simulated a mission whose objective is to position the camera relatively to a landmark made of n = 9 points. Moreover, the environment has been cluttered with two obstacles which may occlude the camera or represent a danger for the base. Let us notice that the sampling period has been defined to be close to our real robot one, that is Ts = 50ms. As shown on figure 3, the vision-based navigation task is correctly performed despite the presence of occluding obstacles. At the beginning of the task, there is no risk of collision. The robot is only controlled by qVS = q(s) and starts converging towards the target. After 50 steps, we simulate an image processing failure (for example due to a sudden illumination change). Our estimation procedure is then enabled and the robot is driven using qVS = q( s) during 200 steps. Meanwhile, the vehicle enters the vicinity of the wall. Then, the global controller progressively switches to qcoll and the robot is constrained to follow the security envelope ξ 0 while tracking a virtual target until the image data is available anew. This event occurs when the circular obstacle leaves the camera field of view, allowing the panplatform control to smoothly switches from ̟ coll to ̟ coll . When there is no more danger, the control switches back to qVS = q(s) and the robot perfectly realizes the desired task. Figure 3 shows also the values of the estimation errors norms ss and zz which both remain close to zero. This result appears to be consistent, as our estimation method is based on a model which reflects precisely the behaviour of the system. We have implemented the above control law on the mobile robot Super Scout II 1 (see figure 4). This is cylindric cartlike vehicle dedicated to indoor navigation. Incremental encoders provide approximate velocity of the robotic system. The robot is equipped with 16 ultrasonic sensors from which we compute the shortest distance d coll and the relative orientation α between the robot and obstacle. A DFW-VL500 Sony color digital IEEE1394 camera captures in YUV4 : 2 : 2 format with 640 × 480 resolution. The robot is controlled by an on-board laptop computer running Linux on which is installed a specific control architecture called G en oM (Generator of Module) [START_REF] Fleury | G en oM : User Manual[END_REF]. The implementation of our control strategy requires different basic services. They are provided by two already existing modules respectively dedicated to the robotic system control and to data sensory processing. First, we have designed a specific module to manage the camera features. It allows to initialize and configure the camera, acquire pictures and process them to extract the visual data s and the occluding object position from which we compute d occ . The image processing pipeline is shown on figure 5. It relies on the following hypotheses to reduce the computation time: all the occluding objects share the same color, the occluding object and the target colors are different, and finally, an 1 The Super Scout II is provided by the AIP-PRIMECA.

IV. EXPERIMENTAL ASPECTS

A. The robotic platform

occluding object detected in the image can never be positioned behind the target in the 3D scene. We have then embedded our control strategy in a second module which offers the required services to realize vision-based navigation tasks amidst obstacles: robotic system initialization, software emergency stop, parameters settings (choice of the different control gains, of the desired visual features s * . . . ), execution of a vision-based task with or without taking into account the occlusion and/or collision phenomena. We present below the obtained experimentation results.

B. Experimental Results

We consider a vision-based navigation task which consists in positioning the embedded camera in front of a given landmark made of n = 4 points. The environment is cluttered with one obstacle, and envelopes ξ -, ξ 0 and ξ + are respectively located at d -= 40cm, d 0 = 56cm and d + = 70cm from the obstacle. The obtained results are presented on figure 6. Figure 6(a) shows the trajectory performed by the robot. At the beginning of the task, there is no risk of collision, and the image data are available. The robot is controlled using qVS and starts converging towards the landmark. It enters the vicinity of the obstacle, inducing a danger of collision. µ coll increases to reach 1 (see figure 6(c)) and the control law applied to the vehicle smoothly switches to controller qcoll (13). Thus, the linear velocity decreases while the angular velocities ω coll and ̟ coll vary to make the robot avoid the obstacle and track the target (see figure 6(b)). When the obstacle is overcome, µ coll vanishes, and the control progressively switches back to qVS = q(s) . However, at the same time, the landmark is manually occluded and the robot is driven using controller q( s) instead of q(s) , keeping on converging towards the target until s reaches s * . Thus, thanks to our method, it is possible to complete a task which would have been aborted otherwise.

Finally, we compare the obtained simulated and experimental results. To this aim, we analyse the quality of the image features reconstruction which appears to be better in simulation as one could have expected. To exhibit the parameters which influence the estimation error, we first recall that our algorithm is initialized using the last visual data available before their lost and on the base of the robot velocities. Therefore, the estimation errors depends mainly on the feature extraction noise (here about 1 pixel) and on the accuracy of the real robot velocities measurement. As the velocities executed by the robot are provided by the embedded encoders with a small delay and drift, we can observe an increase of the estimation error when the control input q switches between the successive controllers (see figures 6(b) and 6(d)). Furthermore, as our platform runs under Linux which is not a real-time operating system, the sampling period T s is not known accurately. These two aspects reduces the efficiency of reconstruction when validating on our real robot. However, despite these problems, figures 6(d) and 6(e) show that the estimation errors remain small, demonstrating the validity of our approach in an experimental context.

V. CONCLUSIONS AND FUTURE WORKS

In this paper, we have elaborated a multi-sensor-based control law allowing a robot to perform safely a visionbased task in an unknown environment despite occlusions. We have first proposed a technique able to reconstruct the visual features when they are suddenly lost. Then, we have integrated it our control law to realize the desired mission. The method relies on the switch between different controllers depending on the risks of collision and occlusion. We have finally validated our work by experimental results.

However, the developed analytical solutions are restricted to landmarks which can be characterized by points and dedicated to a particular robotic system. Thus a natural extension of these works would be to solve this problem for other kinds of visual features and more complex robotic platform. We have also considered a perfect modelling. Therefore, it would be also interesting to evaluate signal processing approaches such as Kalman filtering to take into account the effects of noise on the reconstruction. Finally, a dynamical sequence of the controllers could also be interesting to provide better theoretic feasibility conditions.

( a )Fig. 1 .

 a1 Fig. 1. The systems modelling.

  (a) Obstacle detection (from ultrasonic sensors). (b) Occlusion detection (provided by image processing).

Fig. 2 .

 2 Fig. 2. Avoidance modeling.

Fig. 3 .

 3 Fig. 3. Simulation results.

Fig. 4 .

 4 Fig. 4. The Nomadic Super Scout II.

Fig. 5 .

 5 Fig. 5. Visual features and occluding object extraction algorithm pipeline.

  (a) Robot trajectory. (b) Robot velocities. (c) Collision risk µ coll . (d) Image data estimation errors: ||s -s|| (pixel).(e) z-depth estimation errors: ||z -z|| (m).

Fig. 6 .

 6 Fig. 6. Visual servoing with obstacle avoidance, experimental results.