Alexandre Goldsztejn
email: alexandre.goldsztejn@univ-nantes.fr

Luc Jaulin
email: luc.jaulin@ensieta.fr

INNER APPROXIMATION OF THE RANGE OF VECTOR-VALUED FUNCTIONS

No method for the computation of a reliable subset of the range of vector-valued functions is available today. A method for computing such inner approximations is proposed in the specific case where both domain and co-domain have the same dimension. A general sufficient condition for the inclusion of a box inside the image of a box by a continuously differentiable vector-valued is first provided. This sufficient condition is specialized to a more efficient one, which is used in a specific bisection algorithm that computes reliable inner and outer approximations of the image of a domain defined by constraints. Some experimentations are presented.

Introduction

An image set is a subset of R m that is defined as the image of a subset D of R n by a nonlinear vector-valued function f : R n → R m . This paper deals with the image set problem which consists in computing a reliable description of the image set. More precisely, we want to find two subpavings (i.e. union of boxes) Y -⊆ R m and Y + ⊆ R m such that

Y -⊆ Y ⊆ Y + ,
where Y = {f (x) : x ∈ D} is the image set of D by f . The main contribution of this paper is to present the first method able to solve the image set problem for continuously differentiable nonlinear vector-valued functions f : R n -→ R n (both domain and co-domain have the same dimension).

Interval analysis provides a large number of methods able to compute efficiently an outer approximation Y + of the image set (see [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF][START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF][START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF]). As a consequence, the outer approximation problem can be considered as solved and this is why this paper is mainly devoted to the inner characterization problem. Existing methods are unable to compute an inner approximation Y -of an image set, except for some particular cases :

First, when f maps R n into R, methods based on Kaucher arithmetic (see [START_REF] Herrero | Quantified Set Inversion Algorithm with Applications to Control[END_REF] and extensive references) or on twin arithmetic (see [START_REF] Nesterov | Interval and Twin Arithmetic[END_REF]) have been shown to be able to solve the image set problem. Second, when f is vector-valued and linear the Kaucher arithmetic has been shown to be able to provide an accurate inner approximation of the image set (see [START_REF] Shary | A new technique in systems analysis under interval uncertainty and ambiguity[END_REF][START_REF] Goldsztejn | A Right-Preconditioning Process for the Formal-Algebraic Approach to Inner and Outer Estimation of AE-solution Sets[END_REF]). Finally, some advances [START_REF] Goldsztejn | A Branch and Prune Algorithm for the Approximation of Non-Linear AE-Solution Sets[END_REF][START_REF] Goldsztejn | Inner and Outer Approximations of Existentially Quantified Equality Constraints[END_REF] have been carried out in the general case of constraints ∃y ∈ y, f (x, y) = 0, but none of these methods is able to deal with the problems tackled in the present paper.

Although the method proposed in the present paper is restricted to functions where both domain and co-domain have the same dimension, this class of image set problems is of particular importance for many applications in robotics and control. Let us present two of them.

First, within a set theoretical framework state estimation alternates a prediction step and a correction step (see e.g. [START_REF] Schweppe | Recursive state estimation: unknown but bounded errors and system inputs[END_REF][START_REF] Bertsekas | Infinite-time reachability of state-space regions by using feedback control[END_REF]). The correction step amounts to solving a set inversion problem whereas the prediction step requires the characterization of the image of a subpaving. Existing algorithms based on interval analysis (see e.g. [START_REF] Kieffer | Robust autonomous robot localization using interval analysis[END_REF]) are unable to compute an inner approximation of the feasible set for the state vector at time k. An efficient solution for the image set problem could be really appreciated during the prediction step. Second, the workspace of a serial robot represents the set of all configurations that could be reached by the tool of the robot (see e.g. [START_REF] Khalil | Modlisation, identification et commande des robots[END_REF]). This set can be defined as an image set (as illustrated in the application section of this paper). A guaranteed inner approximation of the workspace is needed during the conception phase of the robot : we need to conceive the robot so that we are certain that it will be able to reach the whole region of interest.

The paper is composed as follows. Section 2 recalls the main notions of interval analysis needed for the proposed developments. A test that makes possible to prove that a box is included in an image set is presented in Section 3. This test is made more efficient using a preconditioning step as detailed in Section 4. Section 5 provides a bisection algorithm that computes an inner and an outer approximation of the image set. The efficiency of the algorithm is demonstrated on two testcases in Section 6. Finally, Section 7 presents some related work.

Interval analysis

The interval theory was born in the 60's aiming at rigorous computations using finite precision computers (see [START_REF] Moore | Interval Analysis[END_REF]). Since its birth, it has been developed and it proposes today original algorithms for solving problems independently of the finite precision of computers computations, although reliable computations using finite precision remains one important advantage of the interval based algorithms (see e.g. [START_REF] Kearfott | Interval Computations: Introduction, Uses, and Resources[END_REF][START_REF] Hayes | A Lucid Interval[END_REF]).

An interval is a set of real numbers

[a, b] = {x ∈ R | a ≤ x ≤ b}
= (a 1 , • • • , a n ) T and b = (b 1 , • • • , b n ) T , is identified with the vector of intervals ([a 1 , b 1 ], • • • , [a n , b n ]) T . The membership is defined componentwise, e.g. x ∈ [a, b] is equivalent to x i ∈ [a i , b i ] for all i ∈ [1.
.n]. The interval vectors are also called boxes. Interval matrices are defined in the same way.

The boundary of an interval x is denoted by ∂x and its interior by int x (these notations also stand for general subsets of R n ). Also, the midpoint of an interval [a, b] is denoted by mid [a, b] = (a + b)/2, its radius is denoted by rad [a, b] = (ba)/2. The same definitions hold for interval vectors and interval matrices (note that rad x ∈ R n while rad A ∈ R m×n ). The magnitude (resp. mignitude

) of an interval [a, b] is |[a, b]| = max{|a|, |b|} (resp. [a, b] = min{|a|, |b|}). For an interval vector [a, b], its magnitude |[a, b]| = max i |[a i , b i ]|.
Given an interval vector x ∈ IR n (resp. an interval matrix A ∈ IR n×m ), we define

Ix := {x ∈ IR n : x ⊆ x} (resp. IA := {A ∈ IR n×m : A ⊆ A}).
Throughout the paper, we use the infinite norms ||x|| = max i {|x i |} and

||A|| = max i j |A ij |.
The core of the interval theory is the extension of real functions to intervals. These extensions allow computing outer approximations of the image of intervals by real functions. Formally, given a real function f : R n -→ R m and an interval function f : IR n -→ IR m , the interval function f is an interval extension of f if and only if, for all x ∈ IR n , range (f, x) ⊆ f (x), where range (f, x) is the image of x by f . The smallest box that contains range (f, x) is called the interval hull of range (f, x) and is denoted by range (f, x). It is the most accurate interval extension that one may hope to construct. Now that the definition of interval extensions is stated, it remains to compute such interval extensions. Computing the interval hull of range (f, x) with arbitrary precision is a NP-hard problem (see [START_REF] Kreinovich | Computational Complexity and Feasibility of Data Processing and Interval Computations[END_REF]). Therefore lots of work has been done to provide interval extensions that are both precise and cheap. Based on the interval arithmetic (that computes exact ranges for elementary functions like +, -, /, ×, exp, cos, etc...) two interval extensions are widely used (see [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF]):

• The natural extension simply consists in replacing the real operation by their interval counterparts in the expression of the function. For example, f (x) = x 2x is an interval extension of f (x) = x 2x. • The mean-value extension consists in rigorously linearizing the function before computing the natural extension of the linearization. For example, the mean-value extension of the function

f (x) = x 2 -x is f (x) = f (mid x)+∆(x-mid x)
where ∆ is an interval that contains {f (x) | x ∈ x} (the interval ∆ can be computed using the natural extension of f ). The quality of an interval extension is usually quantified using its order of convergence: Definition 2.1. Let f : R n -→ R m be a real function. One of its interval extension f : R n -→ R m is said to have a linear order of convergence if for all x ref ∈ IR n there exists γ > 0 such that for all x ∈ Ix ref

|| wid f (x) -wid range (f, x) || ≤ γ|| wid x||.
An interval extension which is locally lipschitz continuous has at least a linear order of convergence. Under mild hypothesis, the natural interval extension has a linear order of convergence. Note that an interval extension which has a linear order of convergence satisfies for all x ref ∈ IR n there exists γ > 0 such that for all x ∈ Ix ref , || wid f (x)|| ≤ γ|| wid x||.

Interval matrices naturally arise when the mean-value extension is used for vector-valued functions. In this case, the derivative of f is the matrix

f (x) =       ∂f 1 ∂x 1 (x) ∂f 1 ∂x 2 (x) • • • ∂f 1 ∂xn (x) ∂f 2 ∂x 1 (x) ∂f 2 ∂x 2 (x) • • • ∂f 2 ∂xn (x) . . . . . . . . . . . . ∂fm ∂x 1 (x) ∂fm ∂x 2 (x) • • • ∂fm ∂xn (x)      
.

The expression of the mean-value extension is then

f (mid x) + J(x -mid x)
where J is an interval matrix that contains {f (x) | x ∈ x}. The following specific definitions about interval matrices will be used in the sequel: consider A ∈ IR n×n , then

• Diag A is the diagonal interval matrix whose diagonal entries are (Diag A) ii = A ii . • Diag -1 A is the diagonal interval matrix whose diagonal entries are (Diag A -1 ) ii = 1/A ii . • OffDiag A is the interval matrix whose diagonal is null and whose off-diagonal entries are (OffDiag A) ij = A ij .

A test for the inclusion in the range of a function

This section presents Theorem 3.1 which provides a sufficient condition for a box y ⊆ R n to be included inside the image of a box x ⊆ R n by a continuously differentiable function f : R n -→ R n . Its corollary 3.1 will be used in the next section as the basis of a specific bisection algorithm that computes an inner approximation of range (f, x). Theorem 3.1. Let x, y ∈ IR n and f : x ⊆ R n -→ R n be a continuous function continuously differentiable in int x. Suppose that the three following conditions are satisfied:

(i) y ∩ range (f, ∂x) = ∅; (ii) y ∩ range (f, Σ) = ∅ where Σ = {x ∈ int x | det f (x) = 0}; (iii) f (x) ∈ y for some x ∈ x.
Then y ⊆ range (f, x).

Proof. x is compact and we can apply Lemma A.2 of Appendix A. Therefore ∂ range (f, x) ⊆ range (f, ∂x) range (f, Σ) .

Using the conditions (i) and (ii), we have

y ∩ ∂ range (f, x) = ∅. (1)
Denote range (f, x) by E (which is compact because x is compact and f is continuous) and f (x) by y. Notice that y ∈ int E because of (1). Consider any z ∈ y and suppose that z / ∈ E. There exists a path included in y between y and z. By Lemma A.1 of Appendix A this path intersects ∂E which is absurd by [START_REF] Bertsekas | Infinite-time reachability of state-space regions by using feedback control[END_REF]. Therefore z ∈ E which concludes the proof. Theorem 3.1 can be used in many ways to construct some inner approximations of the range of a function. The bisection algorithm to be proposed in the next section is based on the following corollary of Theorem 3.1.

Corollary 3.1. Let x ∈ IR n and f : x -→ R n continuous and continuously differentiable in int x. Consider y ∈ IR n and x ∈ x such that f (x) ∈ y. Consider also an interval matrix J ∈ IR n×n that contains all f (x) for x ∈ x. Suppose that 0 / ∈ J ii for all i ∈ [1.

.n]. Then

x + Γ J, (x -x), y -f (x) ⊆ int x (2) =⇒ y ⊆ range (f, x) , where Γ(A, u, b) := (Diag -1 A) b -(OffDiag J)u . Furthermore, (2)
also implies that the matrix J is an H-matrix.

Proof. It is sufficient to prove that the condition (2) implies the three conditions of Theorem 3.1.

(i) Consider an arbitrary x ∈ ∂x. As x ∈ x the mean-value theorem shows that f (x) ∈ f (x) + J(xx) (this is the well-known argument that leads to the interval mean-value extension). It is therefore sufficient to prove that f (x) + J(xx) y = ∅ contradicts (2). So suppose that there exist J ∈ J and y ∈ y such that f (x)+J(x-x) = y. Splitting J to (Diag J)+(OffDiag J) leads to

x = x + (Diag -1 J) y -f (x) -(OffDiag J)(x -x) .
As x ∈ x, y ∈ y, (Diag -1 J) ∈ (Diag -1 J) and (OffDiag J) ∈ (OffDiag J), using the interval arithmetic leads to

x ∈ x + (Diag -1 J) y -f (x) -(OffDiag J)(x -x) .
As x ∈ ∂x this latter appurtenance contradicts (2). (ii) It is sufficient to prove that (2) implies that any real matrix J ∈ J is regular (so Σ = ∅ and (ii) holds). The condition (2) implies

x + (Diag -1 J) y -f (x) -(OffDiag J)(x -x) ⊆ int x for every J ∈ J. As f (x) ∈ y (so 0 ∈ y -f (x) ) the previous inclusion implies x -(Diag -1 J) (OffDiag J)(x -x) ⊆ int x.

This inclusion proves that rad x-(Diag

-1 J) (OffDiag J)(x-x) < rad x.
Using the radius rules proved page 84 of [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF], one concludes that

|(Diag -1 J)(OffDiag J)|(rad x) < rad x. As (Diag -1 J) is diagonal, |(Diag -1 J)(OffDiag J)| = |(Diag -1 J)||(OffDiag J)|
and the i th line of the latter inequality can be written

i =j |J ij | rad x j < |J ii | rad x i As |J ij | rad x j = |J ij (rad x j )| (because rad x j ≥ 0)
, the latter inequality means that J can be scaled to a strictly diagonally dominant matrix. As a consequence, J is an H-matrix and is eventually regular. This also proves that the interval matrix J is an H-matrix (because it contains only real Hmatrices).

(iii) By hypothesis, f (x) ∈ y.

Remark. Corollary 3.1 is closely related to the existence test provided by the Hansen-Sengupta interval operator (see [START_REF] Hansen | Bounding solutions of systems of equations using interval analysis[END_REF]). Indeed the latter is retrieved if y is set to [0, 0] in the expression (2). As a consequence, Theorem 3.1 provides an elementary proof of the existence test related to the Hansen-Sengupta operator.

A preconditioning process

In order to succeed, the test provided by Corollary 3.1 needs the interval evaluation of the derivative of f to be an H-matrix. However, this is not the case in general. The following simple example illustrates this situation.

-→ R 2 defined by f (x) = M x with M = 1 1 -1 1 . Also consider x = ([-2, 2], [-2, 2]) T and y = ([-, ], [-, ]) T .
Obviously, the image of x is obtained from x by a π/4 rotation followed by a √ 2 inflation. Therefore y is a subset of the image of x by f if and only if ≤ 2. The function f has a constant derivative so M is an interval evaluation of its derivative. Now apply Corollary 3.1: on one hand choose x = mid x so x = 0 and f (x) = 0. On the other hand (Diag J) = I and (OffDiag J) = M -I. Therefore Corollary 3.1 leads to the following implication:

[-, ] [-, ] - [-2, 2] -[-2, 2] ⊆ int [-2, 2] [-2, 2] =⇒ y ⊆ range (f, x) .
The sufficient condition is never satisfied whatever is > 0. This is due to the fact that J = M is not an H-matrix.

In such situations, the right thing to do is to precondition the problem in such a way that the interval evaluation of the derivative becomes close to the identity matrix (and therefore possibly an H-matrix). In the present situation, preconditioning the problem consists in replacing the inclusion to be checked, i.e. y ⊆ range (f, x), by the more restrictive one

Cy ⊆ range (Cf, x) , (3)
where C ∈ R n×n is any regular real matrix. The latter is indeed more restrictive than the former: (3) implies Cy ∈ range (Cf, x) for every y ∈ y, that is ∃x ∈ x Cy = Cf (x) . As C is regular, this finally implies y ∈ range (f, x) for every y ∈ y, i.e. y ⊆ range (f, x). As usual, CJ is an interval evaluation of the derivative of Cf provided that J is an interval evaluation of the derivative of f . Therefore, once J is computed using the expressions of the derivative of f , one can choose C ≈ (mid J) -1 in order to obtain a near identity interval evaluation of the derivative of Cf (so that CJ can be an H-matrix). The preconditioned version of Corollary 3.1 leads to the following condition:

x + Γ CJ, (x -x), Cy -Cf (x) ⊆ int x (4) =⇒ y ⊆ range (f, x) .
When the midpoint inverse preconditioning is used (i.e. C ≈ (mid J) -1 ), it is proved in [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF] that CJ is an H-matrix if and only if J is strongly regular. Therefore the strong regularity of the interval evaluation of the derivative of f is a necessary condition for the success of the preconditioned version of Corollary 3.

.5 × [-, ] -0.5 × [-, ] 0.5 × [-, ] + 0.5 × [-, ] - 0 × [-2, 2] 0 × [-2, 2] ⊆ int [-2, 2] [-2, 2] ⇓ y ⊆ range (f, x) .
Therefore y ⊆ range (f, x) is now proved for < 2.

A specific bisection algorithm

This section presents a bisection algorithm that computes a paving included inside the image of a box x by a function f . A bisection of both the domain and co-domain aiming at a direct use of the tests provided in the previous section presents several drawbacks: on one hand, this would result in the bisection of a 2n dimensional space. On the other hand, the bisection would have to provide couples of boxes x and ỹ that satisfy Corollary 3.1. This situation is not easy to reach during a bisection process. Some experimentations have shown that acting in such a way leads to instability of the bisection algorithm. The specific bisection algorithm described in Subsection 5.1 overcomes these difficulties by bisecting only the box x. It is based on a function Inner() which is described in Subsection 5.2. Subsection 5.3 presents a preliminary study of the convergence of the algorithm. Finally Subsection 5.4 provides the modifications needed for the bisection algorithm to deal with domains defined by constraints.

5.1. The bisection algorithm. Given the initial box x ⊆ R n and a function f : x -→ R n , the bisection algorithm bisects only x producing a list of boxes L Domain = {x (1) , • • • , x(N) }. Given x ∈ L Domain , the basic idea is to compute some outer approximation ỹ of range (f, x) using an interval extension of f and to test if ỹ is a subset of range (f, x). This results in a list of boxes L Inside that are proved to be subsets of range (f, x). This is illustrated by Figure 1 for n = 1: the initial interval x has been bisected to twelve intervals (one x being displayed). For each of them, an outer approximation has been computed. The resulting boxes (x, ỹ) T are displayed in white if ỹ ⊆ range (f, x) or in gray if ỹ ⊆ range (f, x). The figure displays seven intervals that are subsets of range (f, x). One can notice that these intervals actually overlap.

The core of this bisection algorithm is therefore a function Inner() that tests if a box ỹ is a subset of the image of x. The information that ỹ is an outer approximation of range (f, x) is of central importance for an efficient implementation of the function Inner() (see Subsection 5.2). Therefore, the function Inner() will have four arguments f , x, x and ỹ. It satisfies

Inner(f, x, x, ỹ) = 1 =⇒ ỹ ⊆ range (f, x) .
The efficient implementation of this function is proposed in Subsection 5.2.

Finally, a simple but efficient bisection algorithm is then easily constructed. It is summarized in Algorithm 1. The algorithm computes both an inner and an outer approximation of range (f, x): on one hand, obviously

y∈L Inside y ⊆ range (f, x) .
On the other hand, the initial box x has been paved into a set of boxes x and for each x an outer approximation of range (f, x) is stored either in the list Inside or in the list Boundary. Therefore, x ← Extract(L Domain ); The next section presents the domain inflation which plays a key role in the implementation of the function Inner(f ,x,x, ỹ). The bisection algorithm has been actually made stable thanks to the domain inflation (experimentations have shown that constructing some bisection algorithms without domain inflation leads to unstable algorithms). where φ is a relation (in the case of the approximation of the range of a function φ(x, y) ≡ f (x) = y). After the bisection of the box x, he may happen to have to prove the quantified proposition

range (f, x) ⊆ y∈ L Inside L Boundary
5 ỹ ← f (x) ∩ f (mid x) + f (x)(x -mid x) ; 6 if || rad x|| ≥ then
∀y ∈ ỹ ∃x ∈ x φ(x, y) , (5)
where x ⊆ x. At this point, it is allowed and certainly useful to inflate x to a bigger box x * that satisfies x ⊆ x * ⊆ x. Then the quantified proposition

∀y ∈ ỹ ∃x ∈ x * φ(x, y) (6)
may be easier to validate. The trivial choice x * = x cannot be done in general because sufficient conditions for the validity of the quantified proposition (6) certainly work only on small enough intervals (e.g. it involves some interval extensions that are precise only for small intervals) while the initial interval x can be very large. Therefore, instead of trying to prove the quantified proposition (5), the task becomes to construct a box x * such x ⊆ x * ⊆ x and such that one can easily prove that the quantified proposition ( 6) is true.

In Algorithm 1, the domain inflation is obligatory as the quantified proposition ( 5) is generally false: indeed the box ỹ is an outer approximation of range (f, x) so the quantified proposition (5) cannot be true. Therefore, the box x has to be inflated before applying Corollary 3.1. This inflation process is illustrated by Figure 2: a box (x, ỹ) T is displayed in the left hand side graphic. Of course ỹ ⊆ range (f, x) because ỹ is an outer approximation of range (f, x) by construction. So x is enlarged to x * in the right hand side graphic in such a way that ỹ ⊆ range (f, x * ) holds. Corollary 3.1 offers an easy and efficient way to enlarge x to x * : we compute the fixed point iteration

x (0) = x, u (k+1) = Γ Cf (x (k) , x (k) -x, Cy -Cf (x) , x (k+1) = x + τ u (k+1) , (7)
where τ > 1 (experimentations have shown that τ = 1.01 is an efficient choice), x := mid x and C = f (x) -1 (note that x and C are fixed during the iteration). We expect this iteration to be contracting, therefore it is stopped as soon as the distance between consecutive iterated does not decrease of at least a given ratio µ (µ = 0.9 has been chosen based on experimentations). In this case, the process fails proving y ⊆ range (f, x). At each step of the iteration, we check if Corollary 3.1 can prove y ⊆ range f, x (k) .

This implementation of the function Inner(f ,x,x,ỹ) is described in Algorithm 2. It is correct as it returns true only if x ⊆ x (Line 6) and if x + Γ(Cf (x), xx, C ỹ -Cf (x)) ⊆ x (Line 8). The latter inclusion allows applying Corollary 3.1 to prove that y ⊆ range (f, x) holds. Finally, termination of Algorithm 2 is obvious: as computations are performed with floating point numbers, the fixed point iteration cannot be µ-contracting forever, and hence the condition d k ≤ µd k-1 will be false after a finite number of steps. 1 cannot succeed if the interval evaluation of the derivative of f is not regular. Therefore, if det f (x) = 0 for some x ∈ x then it will not be possible to prove that f (x) is inside range (f, x) using this x. However, there may exist another x ∈ x such that both det f (x ) = 0 and f (x ) = f (x) and f (x) may eventually be proved to belong to range (f, x). To formalize this idea, the following subset of range (f, x) is defined:

Data: f , x, x ,ỹ Result: Boolean τ ← 1.01; µ ← 0.9; 1 x ← mid x; 2 C ← f (x) -1 ; 3 b ← C ỹ -Cf (x);
range * (f, x) def = y ∈ R n | ∃x ∈ int x y = f (x) ∧ det f (x) = 0 .
Then, the following theorem shows that Algorithm 1 converges asymptotically to range * (f, x). That is, ignoring the finite precision of the computations, any y ∈ range * (f, x) is finally proved to be inside range * (f, x) provided that a small enough is chosen.

Theorem 5.1. Suppose that the interval extensions f and f used in algorithms 1 and 2 locally Lipschitz continuous1 (and hence have a linear order of convergence). Then for all y ∈ range * (f, x) there exists > 0 such that Algorithm 1 succeeds proving that y ∈ range (f, x).

Proof. Let us consider an arbitrary y * ∈ range * (f, x). So, there exists x * ∈ int x such that y * = f (x * ) and det f (x * ) = 0. The theorem is proved by contradiction: we suppose that, no matter how small is chosen, y * is not proved to be inside range (f, x) by the algorithm. Therefore, we can pick up a sequence of boxes x (k) such that x * ∈ x (k) , || wid x (k) || ≤ 1/k and Inner(f, x, x (k) , y (k) ) fails for all k, where y (k) is defined as in Line 6 of Algorithm 1, i.e.

y (k) := f (x (k) ) ∩ f (x (k) ) + f (x (k) )(x (k) -x(k) ) ,
with x(k) := mid x (k) . To build a contradiction, we just have to find a l ≥ 0 for which Inner(f, x, x (l) , y (l) ) succeeds.

For δ ≥ 0 let us define x δ := x * ± (δe), where e i = 1 for all i ∈ {1, . . . , n}. By Proposition A.1 (cases (2) and (3)), there exists δ * , κ * , c 1 > 0 such that, for all x ∈ x δ * , f (x) is nonsingular and ||f (x) -1 || ≤ κ * , and for all δ ≤ δ * and all x ∈ x δ , f (x) -1 f (x δ ) ⊆ I ± (c 1 δ E), where E ij = 1 for all i, j ∈ {1, . . . , n}. Define x * := x δ * and b

* := {f (x) -1 f (x) -f (x) -1 f (x) + f (x) -1 f (x * )(x * -x) : x ∈ x * } (which is well defined because ||f (x) -1 || ≤ κ * ). Note that for all x ∈ x * we have b * ⊇ f (x) -1 f (x * )(x * -x).
For x ∈ x * , consider the function

A x : Ix * -→ IR n×n defined by A x(x) := f (x) -1 f (x). As f is locally Lipschitz continuous, it is Lipschitz continuous inside Ix * , for some Lipschitz constant λ * . As ||f (x) -1 || ≤ κ * , A x is (κ * λ * )- Lipschitz continuous inside Ix * . For some arbitrary δ ≤ min{1/2, δ * }, x ∈ x δ , b ⊆ b * , define g δ,b,x : Ix δ -→ IR n x -→ Γ(A x(x), x -x, b).
Proposition A.2 shows that there exists c 2 > 0 such that the functions g δ,b,x are c 1 c 2 δ-Lipschitz. Furthermore, using Proposition A.1 (case (1)) and Lemma A.4, we show that there exists c 3 > 0 such that

∀x ∈ Ix * , ∀x, ỹ ∈ x * , d g δ,b,x (x), g δ,b,ỹ (x) ≤ c 3 ||x -ỹ||. ( 8
)
Then, define

δ + = min{1/2, δ * , (c 1 (n + 1)) -1 , µ(c 1 c 2 τ ) -1 }
so that both A x(x δ + ) is strictly diagonally dominant and δ ≤ δ + implies that the functions τ g δ,b,x are µ-Lipschitz. Define x + := x δ + ⊆ x * . Lemma A.5 shows that there exist b

+ ∈ IR n such that rad b + > 0, 0 ∈ int b + , b + ⊆ b * and τ g δ + ,b + ,x * (x + ) = x + -x * . ( 9
) Now define b (k) := f (x (k) ) -1 y (k) -f (x (k) ) -1 f (x (k) ) which satisfies b (k) ⊆ b * (see the definition of b * and note that y (k) ⊆ f (x (k) ) + f (x (k) )(x (k) - x(k) )).
As f is a first order interval extension, there exists c > 0 such that || wid y

(k) || ≤ c|| wid x (k) || = c/k, which implies || wid b (k) || ≤ cκ * /k. Finally, define l := cκ * + c 3 + 1 m + and l = max{l , 1/δ + } ( 10
)
where m + := min i |b + i | (which satisfies m + > 0 as 0 ∈ int b + and ±(m + e) ⊆ b + ), so that both x (l) ⊆ x + (because wid x (l) ≤ 1/l ≤ δ + = rad x + while x * ∈ x (l) and x + = x * ± (δ + e)) and l ≥ l .

We are now in position to prove that Inner(f, x, x (l) , y (l) ) succeeds. Fix x = x(l) and define h :

Ix + -→ IR n by h(x) := x + τ g δ + ,b (l) ,x (x) = x + τ Γ(f (x) -1 f (x), x -x, b (l) ).

Note that x computed at Line 11 of Algorithm 2 is actually x ← h(x).

We have proved that h is µ-Lipschitz continuous inside Ix + and we now prove that for all x ∈ Ix + , h(x) ⊆ x + . First, note that by [START_REF] Kearfott | Interval Computations: Introduction, Uses, and Resources[END_REF] we have l ≥ cκ * /m + , which implies that α := m +cκ * /l ≥ 0. For all i ∈ {1, . . . , n} we have

α + wid b (l) i ≤ α + cκ * /l = m + . As 0 ∈ b (l) , we conclude that b (l) i ± α ⊆ ±m + ⊆ b + i , and hence b (l) ± (αe) ⊆ b + . Now, h(x) = x + τ g δ + ,b (l) ,x (x) (using (8)) ⊆ x + τ g δ + ,b (l) ,x * (x) ± (c 3 /l e) ⊆ x * + τ g δ + ,b (l) ,x * (x) ± ((1 + c 3 )/l e) (using (10)) 2 ⊆ x * + τ g δ + ,b (l) ,x * (x) ± (αe) (1 ∈ (f (x) -1 f (x)) ii ) ⊆ x * + τ g δ + ,b (l) ,x * (x) + (±(αe))/(f (x) -1 f (x)) ii (basic properties of IA) = x * + τ g δ + ,b (l) ±(αe),x * (x) (inclusion isotonicity) ⊆ x * + τ g δ + ,b + ,x * (x) (using (9)) = x + .
Therefore, we can apply the Banach fixed point theorem which proves that the fixed point iteration x ← h(x) converges to the unique solution of x = h(x). From the proof of the Banach fixed point theorem, we know that the distance of consecutive iterated decreases of at least a factor µ. Therefore, the condition at Line 6 always succeeds. It just remains to prove that the condition at Line 8 eventually succeeds. This condition can be written x + g(x) ⊆ x with g(x) := g δ + ,b (l) ,x (x).

Let us explicitly enumerate the iterations computed in Algorithm 2 by x(i) , and denote by x * the unique interval vector which satisfies x * = h(x * ) (which hence satisfies τ g(x * ) = x *x). By the Banach fixed point theorem, we know that d(

x (i) , x * ) ≤ µ i /(1 -µ)d(x (0) , x(1) ) =: β (i) and thus d(g(x (i) ), g(x * )) ≤ µτ -1 β (i) (because τ g is µ-Lipschitz continuous). Then, x + g(x (i) ) ⊆ x + g(x * ) ± (µτ -1 β (i) e) = x + τ -1 τ g(x * ) ± (µτ -1 β (i) e) = x + τ -1 (x * -x) ± (µτ -1 β (i) e) ⊆ x + τ -1 (x (i) ± (β (i) e) -x) ± (µτ -1 β (i) e) = x + τ -1 (x (i) -x) ± (γ (i) e) with γ (i) := τ -1 (1 + µ)β (i) . Now note that [-m + , m + ] ⊆ b (l)
j for all j ∈ {1, . . . , n}. Using direct computations, we show that x ∈ x(i) implies both x ∈ xi+1 and b (l) ⊆ xi+1 . By induction, we thus prove that for all i ≥ 1 both x ∈ x(i) and b (l) ⊆ xi hold. Therefore, we can write

x(i) = x+[a (i) , b (i) ] with a (i) j ≤ -m + and m + ≤ b (i) j for all j ∈ {1, . . . , n}. We obtain x + g(x (i) ) ⊆ x + [a (i) /τ -γ (i) , b (i) /τ + γ (i) ]. Now because τ ≥ 1, multiplying m + ≤ b (i) j by (1 -1/τ ) gives rise to b (i) /τ ≤ b (i) -(1 -1/τ )m + , while multiplying m + ≥ a (i) by (1 -1/τ ) gives rise to a (i) /τ ≥ a (i) + (1 -1/τ )m + . Therefore, we have proved that x + g(x (i) ) is a subset of x + [a (i) + (1 -1/τ )m + -γ (i) , b (i) -(1 -1/τ )m + + γ (i) ].
Finally, as γ (i) can be made arbitrary small and (1 -1/τ )m + > 0, we can choose i such that γ (i) ≤ (1 -1/τ )m + . For this i, we have x + g(x (i) ) ⊆ x(i) , and the condition in Line 8 succeeds.

The experimentations presented in Section 6 illustrate this convergence. 5.4. Domains described by constraints. The presented bisection algorithm can be easily extended in order to deal with domains that are defined by constraints. Constrained domains that are considered are of the following form:

D = x E,
where x ∈ IR n and E ⊆ R n . The constraint will be actually described by the following three-valued interval constraint:

c : IR n -→ {0, 1, {0, 1}}.
This function has the following semantic:

• c(x) = 1 =⇒ x ⊆ D; • c(x) = 0 =⇒ x D = ∅; • c(x) = {0,
1} is interpreted as "don't known" so x will have to be bisected.

Example 5.1. Suppose E = {x ∈ R 2 |g(x) ≤ 0}. Then an interval extension g of g can be used in order to obtain the corresponding function c:

• c(x) = 1 if g(x) ≤ 0 and x ⊆ x; • c(x) = 0 if g(x) > 0 or x ∩ x = ∅; • c(x) = {0, 1} otherwise.
Finally, the bisection algorithm is easily updated in order to deal with the constrained domain:

(1) The constrained domain has to be taken into account in the function Inner(f, x, x, ỹ). This is done changing the condition x ⊆ x in Line 6 of Algorithm 2 to

x ⊆ x ∧ c(x) = 1.
As a consequence, the function Inner(f, x, x, ỹ) returns true only if ỹ is an inner approximation of range (f, x) and x is inside D. (2) Line 12 of Algorithm 1 is modified so that only boxes satisfying c(x ) = 0 (c(x ) = 0 respectively) are stored in the list. As a consequence, boxes x (x respectively) such that c(x) = 0 (c(x ) = 0 respectively) are not used anymore. These boxes do not belong to the domain D. Therefore the lists L Inside and L Outside keep their respective semantics.

Note that the asymptotic convergence proved by Theorem 5.1 also holds in the case of constrained domains provided that the interval extensions used to define c are convergent.

Experimentations

This section presents some applications of Algorithm 1. The algorithm has been implemented in C/C++ using the PROFIL/BIAS [START_REF] Knueppel | PROFIL/BIAS -A Fast Interval Library[END_REF] interval library and executed on a PentiumM 1.4Ghz processor. In order to provide clearer results, a regularisation of the overlapping pavings computed by Algorithm 1 has been performed leading to regular pavings. Indicated timings concern only the execution of Algorithm 1. The quality of the couple inner/outer approximations is defined by qual(Inner, Outer) := n vol(Inner) vol(Outer) ,

if vol(Outer) = 0. The closer to one the quality is, the closer inner and outer approximations are. Example 6.1 shows the approximation of the range of a simple function with a constrained domain. It was proposed in [START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF] where a non-reliable approximation was obtained.

= (x, y) T ∈ R 2 | x 2 + y 2 ∈ [1, 2]
, which is represented in the left hand side graph of Figure 3. The following table displays the time needed to compute the approximations and their quality measure for some different values of epsilon. 0.1 0.05 0.025 0.0125 0.00625 qual 0.37 0.63 0.80 0.89 0.94 time 0.09 0.28 0.69 1.9 5.1

In these experimentations, the quality of the approximation obviously improves linearly as decreases, leading to 1 when goes to zero. This is an experimental evidence of the truthfulness of Theorem 5.1. On the other hand, the computation time seems to increase quadratically with : t( ) ≈ 0.000097 -2 +0.017 -1 -0.10. One may be surprised as the boundary of the set to be approximated has dimension 1. However, two pieces of this boundary are built folding the domain and are therefore made of singularities. On these singular boundaries, the algorithm accumulates on a thin surface around the boundary, hence the quadratic dependence w.r.t. -1 . The right hand side graphic of Figure 3 shows the approximation of range (f, D) obtained with Algorithm 1 for = 0.0125.

A more realistic situation is now investigated in Example 6.2. The constraint c(a, θ) is true if and only if the whole robot arm does not meet the obstacle for the controls a and θ. To obtain an expression of this constraint, we can note that d 2 (M, A, B), the square of the distance between a point M to a segment [A, B], is given by

M -A|M -A if B -A|M -A < 0 M -B|M -B if A -B|M -B < 0 M -A|M -A -B-A|M -A 2 B-A|B-A otherwise,
where . | . is the scalar product (writing explicitly these constraints w.r.t. the variables of the problem allows some useful formal simplification). An interval expression of this constraint is easily obtained. The function that computes the working point coordinates from the commands is

f (a, θ) = a + L cos(θ) a + L sin(θ) .
The image of D by f is the workspace of the robot. For our numerical experimentation, we choose C = (3, 1), r = 1 and L = 2. The following table displays the time needed to compute the approximations and their quality measure for some different values of epsilon.

0.1 0.05 0.025 0.0125 0.00625 qual 0.80 0.9 0.95 0.975 0.987 time 0.95 2.3 5.9 15.3 40.2 As previously, the quality of the approximation obviously improves linearly as decreases, leading to 1 when goes to zero. This is again an experimental evidence of the truthfulness of Theorem 5.1. The dependence of the computational time w.r.t. -1 is once more quadratic (t( ) ≈ 0.00072 -2 + 0.14 -1 -0.69). Here again a piece of the boundary is made of singularities and the algorithm accumulates on a thin surface in this area.

The results are plotted on the right hand side of Figure 4 and have been computed using = 0.025.

Related work

In the context of the study of serial robots, some work has been done for providing some description of the workspace of the robot. These works naturally concern the image set problem as workspaces of serial robots consists in the computation of the image set of a vector valued function.

One of the main technique proposed in this framework is to compute the boundary of the workspace (i.e. the boundary of the image set) using some continuation method (see [START_REF] Haug | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF][START_REF] Khalil | Modlisation, identification et commande des robots[END_REF]). The method proposed in the present paper has several advantages: on one hand, the continuation method does not provide a reliable approximation while Algorithm 1 does. On the other hand, Algorithm 1 can deal with constrained domains while the continuation method proposed in [START_REF] Haug | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF] cannot.

Conclusion and perspectives

This paper has presented a method able to compute an inner approximation of an image set in the case where the involved function is vector-valued (with domain and co-domain of the same dimension) nonlinear and differentiable. Although this problem appears in many applications, to our knowledge, it has never been solved before. A general reliable sufficient condition for a box to be a subset of the image of a given box has been provided (Theorem 3.1). It has been specialized to an efficient test (Corollary 3.1) that turned out to be closely related to the existence test associated to the Hansen-Sengupta interval operator. Our algorithm has been made efficient and stable thanks to a preconditioning process and an original domaininflation phase based on a fixed-point iteration. The good behavior of our algorithm has been demonstrated on two academic examples.

The main forthcoming work will be to deal with domain and co-domain that do not have the same dimension. When the domain has a dimension smaller than the one of the co-domain, the image set has a zero volume and an inner approximation is not needed. When the domain has a greater dimension than the co-domain, the image set has an interior and the presented method is not able to compute an inner approximation of the image set without serious adaptations. Now, there exist some applications where the domain has a greater dimension than co-domain: workspaces of overactuated robots (i.e., its number of motors in greater than its number of degrees of freedom) are image sets of vector-valued functions with greater domain dimension than co-domain dimension.

Lemma A.2. Let E ⊆ R n be compact and f : E -→ R n be a continuous function continuously differentiable in int E. Denote {x ∈ int E | det f (x) = 0} by Σ. Then we have ∂ f (E) ⊆ f (∂E) ∪ f (Σ).
Proof. Consider any y ∈ ∂ f (E) . As E is compact and f is continuous, f (E) is also compact. Therefore ∂ f (E) ⊆ f (E) and there exists x ∈ E such that y = f (x). Now suppose that x ∈ int E and x / ∈ Σ. We now prove that it is absurd. As x ∈ int E we have ∃α > 0 B(x, α) ⊆ E . As x / ∈ Σ we can apply the inverse function theorem and we conclude that there exists β such that f is a diffeomorphism in B(x, β). Define γ := min{α, β} and F := f B(x, γ) . Notice that y ∈ F . On one hand F ⊆ f B(x, α) ⊆ f (E). On the other hand, F is the preimage of B(x, γ) by (f

| B(x,β) ) -1 . This latter being continuous, F is open in R n . As a conclusion y ∈ int f (E) which contradicts y ∈ ∂ f (E) . As a conclusion we have x ∈ ∂E or x ∈ Σ. And eventually ∂ f (E) ⊆ f (∂E) ∪ f (Σ).
A.2. Technical results for Section 5. In the rest of the subsection, we define x δ := x * ± (δe) (as previously, e ∈ R n and E ∈ R n×n satisfy respectively e i = 1 for all i ∈ {1, . . . , n} and E ij = 1 for all i, j ∈ {1, . . . , n}). Basic properties of interval arithmetic can be found in [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF] and won't be explicitly mentioned each time they are used.

Proposition A.1. Let f : R n -→ R n be continuously differentiable, and x * ∈ R n such that det f (x * ) = 0. Consider a first order interval extension f of f , e.g. its natural extension. Then, there exist δ * > 0, κ * > 0 and c > 0 such that:

(

1) f (x) -1 is Lipschitz continuous inside x δ * , ( 2
) for all x ∈ x δ * , ||f (x) -1 || ≤ κ * , ( 3
) for all δ ≤ δ * , ∀x ∈ x δ , f (x) -1 f (x δ ) ⊆ I ± (cδE).
Proof.

(1) It is well known that the function M -→ M -1 is continuously differentiable where it is defined. Therefore, f (x) -1 , which is the composition of two continuously differentiable functions, is continuously differentiable where it is defined. As f (x * ) is nonsingular, there exists a neighborhood

x δ * of x * inside which f (x) -1 is continuously differentiable. As x δ * is com- pact, f (x) -1 is Lipschitz continuous inside x * .
(2) ||f (x) -1 || is the composition of two continuous functions, and is therefore continuous. As a consequence, it is bounded above by some κ * > 0 inside the compact set x * .

(3) As f is a first order extension, there exists c such that for all δ ≤ δ * we have || rad f (x δ )|| ≤ c δ. Consider some arbitrary δ ≤ δ * and x ∈ x δ . Then we have

|| rad f (x) -1 f (x δ )|| ≤ ||f (x) -1 || || rad f (x δ )|| ≤ c κ * δ.
Because f (x) ∈ f (x δ ), we have I ∈ f (x) -1 f (x δ ). From the last two results, we obtain f (x) -1 f (x δ ) ⊆ I ± (2c κ * δE). As the latter holds for all x ∈ x δ , this concludes the proof with c = 2c κ * . ii (b * i -

j =i A ij x j ) , Ã-1 ii (b * i - j =i Ãij x j ) .
Now, by Lemma A.3 the latter is less than Thus, the proposition is proved with λ = λ λ + λ .

max i |A * -1 ii | d b * i - j =i A ij x j , b * i - j =i Ãij x j +|b * i - j =i à * ij x j | d(A -1 ii , Ã-1 ii ). Note that |A * -1 ii | = A * ii -1 = (1 -δ) -1 while |b * i -j =i à * ij x j | ≤ ||b * || + (n -1)δ||x * ||. Now d b * i - j =i A ij x j , b * i - j =i Ãij x j ≤ j =i d(A ij x j , Ãij
The proof of the following technical lemma involves generalized intervals and the Kaucher arithmetic (see [START_REF] Shary | A new technique in systems analysis under interval uncertainty and ambiguity[END_REF] and references therein). Details on this generalization of the interval arithmetic are not given here as it is used only in the following short proof. Proof. Using the Kaucher arithmetic group properties, Γ(A, x, b) = x is equivalent to, for all i ∈ {1, . . . , n}

b i = (dual A ii )x i + j =i (dual A ij )(dual x j ) = ± A ii δ + j =i ∓|A ij | δ,
where ∓ := dual ± . Now, rad b i = δ(A iij =i |A ij |), which is strictly positive as A is strictly diagonally dominant, so b is proper. Furthermore mid b i = 0 as it is the sum of zero-centered intervals. Finally, by inclusion monotonicity we have b i = A ii x i + j =i A ij x j = Ax.

  where a, b ∈ R such that a ≤ b are called respectively the lower and upper bounds of [a, b]. Intervals and related objects (e.g. interval functions and interval vectors and matrices) are denoted by boldface letters. The set of intervals is denoted by IR. Interval vectors and interval matrices are defined like in the context of real numbers. They are identified with vectors of intervals and matrices of intervals. For example, the interval vector [a, b], where a

Example 4 . 1 .

 41 Consider the linear function f : R 2

1 .Example 4 . 2 .

 142 Example 4.1 is now revisited thanks to the presented preconditioning: Consider f , x and y as in Example 4.1. Choose C = M -1 so that Cf (x) = x. Therefore the interval evaluation of the derivative of Cf is I. The preconditioned version of Corollary 3.1 leads to the following condition: 0

Figure 1 .

 1 Figure 1. Bisection of the initial domain x leading to a list of boxes x and to some outer approximations ỹ of their image.

1 L

 1 y. Data: f , x , Result: L Inside (list of boxes), L Boundary (list of boxes) L Inside : empty list of boxes; Domain : empty list of boxes (sorted by decreasing radius); 2 Store the box x in L Domain ; 3 while L Domain not empty do

4

 4

7 if 9 else 10 Bisect 17 return 18 Algorithm 1 :

 791017181 Inner(f ,x,x,ỹ) then 8 Store the box ỹ in L Inside ; the box x to obtain x and x ; 11 Store x and x in L Domain ; (L Inside ,L Boundary ); The bisection algorithm

Figure 2 .

 2 Figure 2. Enlargement of the box x to x * aiming at a successful application of Corollary 3.1 applied to x * and ỹ.

4 d 10 x 13 Algorithm 2 : 5 . 3 .

 41013253 k ← +∞; d k-1 ← +∞; /* The while loop is run at least5 twice */ while d k ≤ µd k-1 ∧ x ⊆ x do 6 u ← Γ(Cf (x), xx, b); 7 if x + u ⊆ x then return True; 8 d k-1 ← d k ; 9 d k ← d(x, x + τ u); ← x + τ u; FunctionInner() Convergence of the algorithm. The convergence of the algorithm is now investigated. The test based on the preconditioned version of Corollary 3.

Example 6 . 1 .

 61 Consider the function f (x) = xy x + y and the domain D

Figure 3 .

 3 Figure 3. Domain D (left hand side graph) and the reliable approximation of its image (right hand side graph).

Figure 4 .

 4 Figure 4. A serial robot with an obstacle (left hand side graphic) and its workspace reliable approximation (right hand side graph). The dotted circle in the workspace approximation has been added to represent the obstacle.

Lemma A. 3 .

 3 Let a, b, a , b , a , b ∈ IR be interval such that a ⊆ a, a ⊆ a, b ⊆ b and b ⊆ b. Then d(a b , a b ) ≤ |a|d(b , b ) + |b|d(a , a ). Proof. We have d(a b , a b ) ≤ d(a b , a b ) + d(a b , a b ) by the triangular inequality. Then d(a b , a b ) ≤ |a |d(b , b ) ≤ |a|d(b , b ), and d(a b , a b ) ≤ |b |d(a , a ) ≤ |b|d(a , a ). Lemma A.4. Let A * = I ± (δE) with δ ≤ 1/2, and x * , b * ∈ IR n . Then, there exist λ , λ > 0 such that for all δ ≤ 1/2, Γ( . , . , b) : IA * ×Ix * -→ IR n is λ δ-Lipschitz w.r.t. A and λ δ-Lipschitz w.r.t. x, where λ and λ depend only on x * and b * . Proof. Let us compute the Lipschitz constant w.r.t. A, the Lipschitz constant w.r.t. x being similar and simpler to compute. Let A, Ã ⊆ A * and x ⊆ x * . Then d(Γ(A, x, b * ), Γ( Ã, x, b * )) is equal to max i d A -1

  x j ) ≤ j =i |x * j | d(A ij , Ãij ) ≤ j =i ||x * || 2δ = 2(n -1)||x * ||δ, while d(A -1 ii , Ã-1 ii ) ≤ |A * -2 ii |d(A ii , Ãii ) ≤ (1δ) -2 2δ. So far, we have proved that Γ( . , . , b) is Lipschitz continuous of constant 2(n -1)||x * ||(1δ) -1 δ + 2(||b * || + (n -1)δ||x * ||)(1δ) -2 δ w.r.t. A. The proof is concluded noting that δ ≤ 1/2 which implies (1δ) -1 ≤ 2.Proposition A.2. Let x * ∈ IR * and A : Ix * -→ I I ± (δE) , with δ ≤ 1/2, be λ -Lipschitz. Let x ∈ R n and b * ∈ IR n and an arbitrary b ⊆ b * . Define T : Ix * -→ IR n by T(x) := Γ A(x), xx, b . Then, there exists λ > 0 such that for all δ ≤ 1/2, T is λδ-Lipschitz continuous, where λ depends only on x * and b * . Proof. The proposition is proved using Lemma A.4. Fix an arbitrary δ ≤ 1/2. Then, d( Γ(A(x), xx, b) , Γ(A(y), yx, b) ≤ d( Γ(A(x), xx, b) , Γ(A(y), xx, b) + d( Γ(A(y), xx, b) , Γ(A(y), yx, b) ≤ λ δ d(A(x), A(y)) + λ δ d(xx, yx) ≤ λ λ δ d(x, y) + λ δ d(x, y).

Lemma A. 5 .

 5 Let A ∈ IR n×n be a strictly diagonally dominant interval matrix and 0 < δ ∈ R. Define x ∈ IR n by x i = ±δ. Then there exists an unique b ∈ IR n such that rad b i > 0 and Γ(A, x, b) = x. Furthermore, mid b = 0 and b ⊆ Ax.

The natural interval extensions can be used under mild hypothesis like abs(x) is not evaluated for intervals which contain 0.

Note that[START_REF] Kearfott | Interval Computations: Introduction, Uses, and Resources[END_REF] implies 1 + c3 + cκ * ≤ lm + which implies (1 + c3)/l ≤ α.

Acknowledgments. The authors would like to thanks Nacim Ramdani for his helpful comments and Gilles Chabert and Nicolas Delanoue for their careful reading of the paper and their helpful comments.

Appendix A. Technical results

A.1. Technical results for Section 3.

Lemma A.1. Let E be compact in R n and y ∈ int E and z / ∈ E. Any continuous path connecting y to z intersects ∂E.

We will prove that if the continuous path w : [0, 1] -→ R n such that w(0) = y and w(1) = z does not intersect ∂E then it is disconnected. This fact being absurd because the image of the connected set [0, 1] by the continuous function w must be connected. On one hand

). Also, z ∈ w([0, 1]) ∩ CE and therefore this set is nonempty. Therefore, w([0, 1]) ⊆ (int E ∪ CE) implies w([0, 1]) is disconnected, which is absurd. Therefore, w([0, 1]) ∩ ∂E = ∅.


Hanh Nguyen-Ngoc

Canh Tran-Minh

Sol-gel process for vegetal-cells encapsulation

Keywords: sol-gel, silica, cell entrapment, vegetal cells, Chlorella vulgaris, fluorescence detection ___________________________________________________________________________

A sol-gel process was carried out at room temperature to obtain porous silica matrixes for entrapment of vegetal cells. The process uses sodium silicate as a precursor rather than alkoxides to avoid generation of by-products that is detrimental to Chlorella vulgaris cells.

The influence of several factors influencing the gelification was explored to optimize the encapsulation process. The technique produces a transparent matrix and is then very suitable for vegetal cells entrapment since they require light for their photosynthetic activity. Activity of the algal cells was determined by measurement of the intensity of fluorescence emission at 682 nm. This measurement also allowed to evaluate the efficiency of the immobilization technique and assess the long term activity of the encapsulated whole-cells.

Introduction

Conventional glass-making processes require reaction at elevated temperatures which are not compatible with entrapment of fragile biomolecules. The sol-gel process involves the transition of a system from a molecular precursor state through the formation of nano-sized bricks during the sol colloidal phase into a solid gel phase and finally transitions into a dried ceramic material. Silica matrixes are relatively inexpensive to synthetize and have interesting properties including optical transparency, biocompatibility and chemical inertness. The design flexibility of sol-gel technique and ease of fabrication can fulfill to create the surfaces with structural and chemical features that could be compatible with biomaterials.

Immobilization of chemicals in such inorganic structures is not restricted to organic molecules. Since the sol-gel process can be carried out in moderate conditions, the structural elements of the nano-sized sol bricks are preserved and enable the encapsulation of biomolecules [1][2] and even bacteria [3]. Biomolecules entrapped in sol-gel matrixes typically exhibit improved resistance to thermal and chemical denaturation, and increased storage and operational stability. This technique has been applied to immobilize various biological molecules (proteins, enzymes, antibodies) using the alkoxide route with tetraethyl orthosilicate (TEOS) or tetramethyl orthosilicate (TEMOS) as precursors [4], but they are not really cell compatible because the by-products ethanol or methanol are harmful for the biological species [5]. These alcohols can be removed from the acidic aqueous sol prior gelation [6] or the alkoxides can be replaced by aqueous precursors [7][8]. In this work, micro-algae cells of Chlorella vulgaris as an example of vegetal cells are entrapped in sol-gel matrixes using sodium silicate as a precursor to preserve the biological activity.

Experimental

Chemicals

Silica sources were sodium silicate solutions (purchased from Riedel-de-Haen) and colloidal silica LUDOX HS-40 (from Aldrich). Glycerol was purchased from Aldrich. All solutions were prepared with purified water (Rios Millipore system).

Cells culture

The Chlorella vulgaris strain (CCAP 211/12) was purchased from The Culture Collection of Algae and Protozoa at Cumbria, United Kingdom. The axenic algal strain was grown in the culture medium and under conditions described by the International Organization for Standardization (ISO 8692) [START_REF]Water quality : fresh water growth inhibition test with Scenedesmus subspicatus and Selenastrum capricornutum[END_REF]. The concentration of algae in suspension was determined by optical density (OD).

Cell entrapment

Sodium silicate (0.4M, 4ml) and LUDOX (8.5M, 4ml) were thoroughly mixed to obtain an homogeneous silica solution. A HCl, 4M solution was then added drop by drop until an appropriate pH is reached to induce the gelation process. Immediately, an algal solution containing 1.3x10 8 cells/ml and 10% (w/w) glycerol was introduced under stirring and the resulting solution is cast in moulds or in cuvettes to produce silica matrixes containing algal cells. Gelation was performed at ambient temperature.

Microscopy of entrapped cells

Photomicrographs of silica membranes containing algal cells were obtained from an Olympus BH2 upright microscope. With an objective lens (x100) and an eyepiece (x10), the overall magnification is x1000.

Fluorescence measurements

Fluorescence measurements were achieved with a Spex Fluorolog 2 fluorometer from Jobin-Yvon equipped with a microcomputer for data recording. The algal matrixes were illuminated with a 469-nm excitation light and the fluorescence was collected at 682 nm.

Fluorescence assays were carried out in cuvettes with algal cells in solution or entrapped in a silica matrix. For long-term activity, the algal cells are stored at 4°C between measurements.

Results and Discussion

Sol-gel process to obtain active algal cells in silica matrixes

Whole-cells (bacteria, yeasts, fungi, plant tissues, mammalian tissues) are much larger than enzymes and can thus be immobilized intact using physical entrapment [START_REF] Fukui | [END_REF]. In the past, they have been trapped in various supports: agar or polyacrylamide gels [11], dialysis membranes [12], collagen suspensions [13], glass microfiber filter or alginate gel [15]. The susceptibility of alginate beads to phosphate chelation involve a complex interaction of cation type, concentration, and pH of phosphate solution. Gel swelling and shrinking along with deterioration of alginate structure in the presence of phosphate ions limit the use of this technique. Vegetal cells require light for their photosynthetic activity and detection, so they have been immobilized on the surface of the support [16]. Basically, immobilization of live cells is very similar to the enzyme counterpart. Any entrapment technique in polymeric structure for enzyme immobilization [17] could be applied to whole cells. One of the problems is the mass transfer resistance imposed by the fact that the substrate has to diffuse to the reactive site and inhibitory or toxic products must be removed to the environment. Solgel silica matrixes enable those possibilities because they are mesoporous materials with regularly arranged, uniform mesopores that facilitate chemical transfer between the immobilized cells and the bulk solution.

To entrap the whole cells, the usual sol-gel process is based on acid hydrolysis and condensation reactions in solution followed by immobilization of the biocomponent at neutral pH to preserve the cell activity. Polyvinylalcohol and polyvinylpyrrolidone were used in some cases to minimize cracks in the gel microstructure [18], but glycerol is preferred since they also improve the mechanical properties of the polymer [19]. Unlike methanol or ethanol, glycerol appears to have no inhibition effect on the algal cells, in terms of fluorescence emission. Glycerol at 10% (w/w) was added to the algal solution at the beginning of the immobilization process and the algal activity in the glycerol containing matrix was regularly tested by fluorescence measurements during 4 weeks. No inhibition effect was observed compared to the reference matrix without glycerol. A characteristic feature of the sol-gel silica matrix is the relative optical transparency of the resulting insoluble support. It makes the technique very suitable for vegetal cells immobilization since they require light for their photosynthesis and activity. It also enables the algal cells entrapped inside the support to be seen with an upright microscope and compared to the algal-free matrix (Fig. 1A), particularly when the vegetal cells are illuminated and photosynthesis is triggered. Using the same algal concentration (ie. 6x10 7 cells/ml) and pH 7, the fluorescence emission is 25-30x10 5 a.u. with silica matrix and 10-15x10 5 a.u. in BSA (bovine serum albumin)+glutaraldehyde. On the photograph, some algal cells exhibit higher light intensity under illumination probably due to their proximity to the surface (Fig. 1B). Of course, no algal cells are detected on the reference support .

Microscope images of the silica support and the algae encapsulated matrix

The BET method used gave an average pore size of 10 nm in diameter which corresponds to a pore volume of 0.40 cm 3 /g. Since the surface area also depends on gelation pH (178 m 2 /g for pH 7), the gelation time will affect the porosity of the matrix. The fact that the cells are entrapped in the inner part of the matrix prevent them from microbial attack and leaching out.

On the other hand, all the chemical and biological transformations can be continuously observed and monitored by means of optical techniques: fluorescence, luminescence or light absorption. An optimized polymerization pH should allow an efficient entrapment of the algal cells while preventing loss of cells activity under polymerization conditions. The effect of cells encapsulation pH on the gelation time in the silica matrix was studied for pHs ranging from 6 to 8.5 (Fig. 2). A minimum of gelation time at pH 7-7.5 is noted for different cell concentrations in silica gel. Gelation time becomes higher with increasing pH and no gel is formed beyond pH 8.5. The particles are so charged [START_REF] Iler | Polymerisation of silica[END_REF] that they cannot agglomerate to undergo the polymerization process. The amounts of immobilized algal cells also affect the gelation time. High algal concentration tends to increase the gelation time particularly for pH>7.

Effects of algal concentration and pH on gelation time

Viability of algal cells after encapsulation

The algal cells viability is determined by measuring the chlorophyll fluorescence which reflects their photosynthetic activity. The fluorescence signal also provides information on the fluorescence induction and the activity of the Photosystem II (PSII), which is the target of many herbicides. Fig. 3 Viability of algal cells in terms of fluorescence emission at 682 nm when algal cells are in solution (A and B) or entrapped (C and D) in a silica matrix, compared to a reference (E) silica support without any cells (algal concentration C 1 =13x10 7 cells/ml, C 2 =6x10 7 cells/ml ) Fig. 3 shows that the presence of algal cells can be detected by fluorescence emission at 682 nm. The emission intensity is lower when the cells are entrapped in the silica matrix than in solution. This results from the partial absorption of light by the insoluble support. No fluorescence at 682 nm was emitted with the algal-free matrix. In any case, emission intensity is reduced with decreasing cell concentrations.

Effect of pH on algal cell long-term activity

The activity of algal cells depends on the working pH of the sample solution. The effect of pH on the long-term behaviour of the Chlorella vulgaris immobilized in a silica matrix was studied with 0.2 M phosphate buffer solution. The pH of the solution ranges from 6 to 8.

The cells activity expressed in terms of fluorescence emission is shown in Fig. 4. Although there are fluctuations in the values, a general tendency can be observed from those various curves. Fluorescence intensity increases with high pH values particularly at pH 8.

Fluorescence intensity is higher at pH 8 because the activity of PSII system is affected at this pH, which induces a greater fluorescence. The same effect is observed with anti PSII herbicides. High cells concentrations slightly give rise to a high fluorescence emission for fresh cells but a decrease in fluorescence can be noted after 4 weeks: it could be attributed to the fluorescence quenching by a higher number of aging cells and their metabolites. On the whole, the activity of Chlorella vulgaris remains stable during the 4 week period of assays.

Application to herbicides detection

The principle of herbicides detection is based on the effect of herbicides on algal fluorescence emission. The presence of anti PSII herbicides (diuron, atrazine, simazine and isoproturon) increases the fluorescence emission. Optical biosensors based on fluorescence detection require a translucent matrix in which a fluorescent bioreceptor can be immobilized and placed in front of the tip of an optical fiber. The optical fiber is used to send the excitation radiation to the algal cells and convey the fluorescence radiation up to a fluorometer [14], to generate a significant signal. The optical transparency of the algal sol-gel silica matrix would allow its application to algal biosensors.

Conclusion

In this study, the possibility of using a sol-gel method to immobilize Chlorella vulgaris has been investigated. The use of this vegetal cells entrapment technique is motivated by the optical transparency and chemical inertness of the support. We have shown that aqueous precursors prevents the loss of algal activity during the polymerization process. This technique also enables the use of this algae containing matrix in phosphate environment, which is damaging for alginate gel [START_REF] Dainty | [END_REF] . The measurement of fluorescence emission at 682 nm provides an easy method to control the activity of algal cells and optimize the various immobilization parameters. Association of algal silica matrixes with optical fibers is under way to construct a optical biosensor for determination of pollutants. This sol-gel technique can be extended to other micro-organisms immobilizations for the design and construction of whole-cell biosensors.

Fig. 1 :

 1 Fig. 1 : Silica support without (A) and with (B) entrapped algal cells

Fig. 2 :

 2 Fig. 2 : Effect of pH and algal concentration on the gelation time

Fig. 4

 4 Fig. 4 Long term activity of encapsulated algal cells as a function of pH for various cells concentrations

University of Technology HCM,

Ly Thuong Kiet, Ho Chi Minh city Vietnam * Corresponding author. Tel: +33 (0) 477 420033; fax: +33 (0) 477 420000, E-mail address: tran@emse.fr,

Cours Fauriel, 42023 Saint Etienne cedex 2, France

Acknowledgments

This work has been financially supported by Région Rhône-Alpes (MIRA fellowships and Thématiques program) and by CNRS (ACI program). The authors appreciate Dr. Claude Durrieu, Dr. Martine Naessens and Prof. Jean-Marc Chovelon for fruitful discussions and Anne-Marie Danna for experimental assistance.