Audric Drogoul

Gilles Aubert

Didier Auroux

TOPOLOGICAL GRADIENT FOR THE DETECTION OF FINE STRUCTURES IN 2D AND 3D IMAGES

Keywords: Object detection, Fine structures, Image segmentation, Calculus of variations, Topological Gradient

Topological gradient for a fourth order PDE and application to the detection of fine structures in 2D and

3D images

INTRODUCTION

In image processing, segmentation / restoration or detection of fine structures are challenging problems with many applications (in satellite, medical, biological imaging, ...). In this work we give and experiment a variational model to detect fine structures (filaments and points in 2D, surfaces, filaments and points in 3D) by using the topological gradient method. First introduced by Sokolowski [START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF] and Masmoudi [START_REF] Masmoudi | The topological asymptotic[END_REF], this notion consists in the study of the variations of a cost function j(Ω) = J Ω (u Ω ) with respect to a topological variation, where J Ω (u) is of the form J Ω (u) = Ω F (u, ∇u, ∇ 2 u, . . . ) ; u Ω is a solution of a PDE defined on the image domain Ω. In order to calculate the topological gradient, we remove to Ω a small object ω ǫ of size ǫ → 0 centered at x 0 (generally a ball or a crack), and we compute the limit I(x 0 ) = lim ǫ→0 j(Ω\ωǫ)-j(Ω) ǫ d

where d is the dimension of the ambiant space. I(x 0 ) is called the topological gradient at x 0 . A particularity of this method is that the computation of the topological gradient only needs the direct state u Ω and an adjoint state v Ω solution of a similar PDE depending on u Ω . This makes the topological gradient computation easy and very fast. Initially applied in structural mechanics, this method has been used in image processing by several authors, for example by Belaid et al. [START_REF] Auroux | Image restoration and classification by topological asymptotic expansion[END_REF] in restoration / segmentation problems. In this case F (∇u) = |∇u| 2 and u Ω is the solution Contact: drogoula@unice.fr, gaubert@unice.fr, auroux@unice.fr of a Laplace PDE. The focus of the method is to find the most energetic points, associated to high topological gradients. If the use of the gradient operator for the cost function is classical in edge detection, it is known [START_REF] Steger | An unbiased detector of curvilinear structures[END_REF] that this choice is not adapted for filament detection: the "gradient does not see" these structures. To illustrate this fact, let us consider in 1D the function f (x) = 0 if x = 0 and f (0) = 1. This function can be approximated by the function

f η (x) = 0 if |x| ≥ η and f η (x) = 2 η 3 |x| 3 -3 η 2 |x| 2 + 1 if |x| ≤ η. We have f ′ η (0) = 0 but f ′′ η (0) = -6 η 2 , thus f ′ η "does not see" 0 but f ′′
η becomes singular at 0. Other variational models have been proposed in the literature according to applications, see [START_REF] Jacob | 3D shape estimation of DNA molecules from stereo cryo-electron micro-graphs using a projection-steerable snake[END_REF] for the detection of biological filaments or [START_REF] Rochery | New higherorder active contour energies for network extraction[END_REF] for road network detection. In [START_REF] Baudour | Detection and Completion of Filaments: A Vector Field and PDE Approach[END_REF] authors propose a model for detecting objects of codimension two and one in 2D images. Their method is inspired by Ginzburg-Landau models. There exists of course other approaches that are not based on variational calculus. In [START_REF] Zana | Segmentation of vessel-like patterns using mathematical morphology and curvature evaluation[END_REF] a morphological method is presented ; authors use morphological filters and a curvature evaluation to detect vessel-like patterns. In [START_REF] Stoica | A Gibbs Point Process for Road Extraction from Remotely Sensed Images[END_REF] a thin network is simulated by a point process which penalizes unconnected segments and favors aligned segments. The estimate of the network is obtained by minimizing an energy function. Finally let us mention the wavelet approach [START_REF] Liew | Application of wavelet theory for crack identification in structures[END_REF]. In this paper we present a new variational model based on the topological gradient method, inspired by the static PDE of a deflection of a thin plate subject to transverse force and bending moment. In this case F (∇ 2 u) = ∇ 2 u 2 and u Ω is the solution of a PDE based on the Bilaplacian operator. The model presented here allows us to detect filaments and points in 2D blurred and / or noisy images (surfaces and filaments in 3D images). The paper is organized as follows. In section 2 we present the problem and give the main idea to compute the topological gradient (TG) in the case of a crack. In section 3, we develop the numerical algorithm for computing the TG and give some numerical results in 2D and 3D.

PROBLEM STATEMENTS AND COMPUTATION

OF THE TOPOLOGICAL GRADIENT

Problem statement

We suppose that the observed image f writes as f = Ku + b where K is a blurring operator, b a Gaussian noise and u the image to recover. We denote by Ω ǫ the domain Ω\ω ǫ where ω ǫ = x ǫ , x ∈ ω with ω a crack or a ball. We introduce the cost function and the PDE proposed in [START_REF] Aubert | A topological gradient based model for the detection of fine structures in 2D images[END_REF]. The model is inspired by the Kirchhoff thin static plate model subject to pure bending (see [START_REF] Zhou | Boundary Element Methods with Applications to Nonlinear Problems[END_REF]) with a Poisson ratio ν = 0). We denote by J ǫ (u) = J Ωǫ (u) the cost function defined by:

J ǫ (u) = Ωǫ ∇ 2 u 2 , (1)
where for a matrix M , M 2 = tr(M T M ). Let u ǫ = u Ωǫ be a regularization of the observed image f ∈ L 2 (Ω ǫ ), solution of the following minimization problem :

min u∈H 2 (Ωǫ) αJ ǫ (u) + Ku -f 2 L 2 (Ωǫ) (P ǫ ) (2)
where α > 0 is a parameter that we have to tune.

H 2 (Ω ǫ ) = u ∈ L 2 (Ω ǫ ), ∇u ∈ L 2 (Ω ǫ ), ∇ 2 u ∈ L 2 (Ω ǫ ) and K : L 2 (Ω ǫ ) → L 2 (Ω ǫ ) is a linear operator.
The variational formulation of (P ǫ ) writes as:

find u ǫ ∈ H 2 (Ω ǫ ) : a ǫ (u ǫ , v) = l ǫ (v), ∀v ∈ H 2 (Ω ǫ ) (3)
with the following bilinear a ǫ (u, v) and linear l ǫ (v) forms:

a ǫ (u, v) = Ωǫ α 1≤i,j≤2 ∂ 2 u ∂x i ∂x j ∂ 2 v ∂x i ∂x j + KuKv l ǫ (v) = Ωǫ K ⋆ f v. (4)
where K ⋆ denotes the adjoint operator of K. The Euler equations associated with (P ǫ ) are

(P ǫ ) α∆ 2 u ǫ + K ⋆ Ku ǫ = K ⋆ f, on Ω ǫ B 1 u ǫ = B 2 u ǫ = 0, on ∂Ω ǫ (5)
where

B 1 u = ∂ n (∆u) -∂ σ n 1 n 2 ∂ 2 u ∂x 2 1 - ∂ 2 u ∂x 2 2 -(n 2 1 -n 2 2 ) ∂ 2 u ∂x 1 ∂x 2 B 2 u = n 2 1 ∂ 2 u ∂x 2 1 + n 2 2 ∂ 2 u ∂x 2 2 + 2n 1 n 2 ∂ 2 u ∂x 1 ∂x 2
setting n = (n 1 , n 2 ) the outer normal to the domain, and σ = (σ 1 , σ 2 ) the tangent vector such that ( n, σ) forms an orthonormal basis.

Computation of the topological gradient in 2D

The calculus of the topological gradient is very technical. In this section we just give the important steps and the main ideas in the case of the crack with α = 1 and K = I d . Notations. To simplify, we suppose that the crack σ is centered at 0 and writes as σ = {(s, 0), -1 < s < 1}. We denote by σ a fixed smooth closed curve that contains σ and by ω the set such that ∂ ω = σ. For x ∈ σ, u + and u -denote the limit values lim y→x,y∈ ω c u(y) and lim y→x,y∈ ω u(y) and [u] = u + -u -the jump across σ. We set σ ǫ = x, x ǫ ∈ σ , Λ the exterior space R 2 \σ and Ω ǫ the cracked domain Ω\σ ǫ .

Main ideas. To compute the topological gradient we have to evaluate the leading term in the difference J ǫ (u ǫ ) -J 0 (u 0 ) where u 0 is the solution of (3) with ǫ = 0. We introduce for that an adjoint problem v ǫ ∈ H 2 (Ω ǫ ) solution of :

a ǫ (u, v ǫ ) = -L ǫ (u), ∀u ∈ H 2 (Ω ǫ ) (6)
with

L ǫ (u) = Ωǫ K ⋆ (f -2Ku 0 )u.
Then by integration by parts we get

J ǫ (u ǫ ) -J 0 (u 0 ) = σǫ B 1 u 0 [w ǫ ] -B 2 u 0 [∂ n w ǫ ] -A ǫ (7)
with w ǫ = v ǫ -v 0 (v 0 is the solution of ( 6) with ǫ = 0),

A ǫ = u ǫ -u 0 2 L 2 (Ωǫ)
and where [w ǫ ] and [∂ n w ǫ ] denote the jumps across σ ǫ of w ǫ and ∂ n w ǫ . Next we approximate w ǫ by w ǫ = ǫ 2 P x ǫ + e ǫ where e ǫ ∈ H 2 (Ω ǫ ) is such that e ǫ H 2 (Ωǫ) = O(ǫ 2 log(ǫ)) and where P ∈ W 2 (Λ)/P 1 is the solution of the exterior problem defined by (P ext ) ∆ 2 P = 0, on Λ B 1 P = 0, and

B 2 P = V 2 on σ (8)
where

V 2 = -∂ 2 v0 ∂x 2 2 (0) is such that B 2 (v 0 )(x) = V 2 + O(ǫ).
We recall that W 2 (Λ)/P 1 is the weighted Sobolev space of functions defined up to the polynomial functions of degree less than or equal to one. We know [START_REF] Aubert | A topological gradient based model for the detection of fine structures in 2D images[END_REF] that P writes as a multilayer potential :

P (x) = σ λ 1 (y)B 1,y E(x-y)dσ y + σ λ 2 (y)B 2,y E(x-y)dσ y
where denotes the principal Cauchy value and E(x) is the fundamental solution associated to the Bilaplacian operator. We can show that

[P ] = λ 1 = 0 and [∂ n P (x)] = λ 2 (x) = 4 3 V 2 1 -x 2 1
. Thus (7) expresses as

J ǫ (u ǫ ) -J 0 (u 0 ) = - 2π 3 ǫ 2 ∂ 2 u 0 ∂x 2 2 (0) ∂ 2 v 0 ∂x 2 2 (0) + o(ǫ 2 )

Expression of the topological gradient in the general case in 2D

The topological gradient of the function ǫ → J ǫ (u ǫ ) with u ǫ given by ( 3) and for a crack centered at x 0 and of normal n is (see [START_REF] Aubert | A topological gradient based model for the detection of fine structures in 2D images[END_REF]) :

I(x 0 , n) = - 2π 3 ∇ 2 u 0 (x 0 )( n, n)∇ 2 v 0 (x 0 )( n, n) (9)
where u 0 and v 0 are respectively given by ( 5) and ( 6) with ǫ = 0. We deduce from ( 9) an indicator based on the topological gradient by taking the maximum of its absolute value over n :

I(x 0 ) = max n =1 |I(x 0 , n)| ( 10
)
Remark 2.1 In 3D we can model a planar circular crack by σ = {(rcos(θ), rsin(θ), 0), 0 ≤ r < 1, 0 ≤ θ < 2π} written in the orthonormal basis ( u, v, n); by similar calculus we have the following TG expression :

I(x 0 , n) = - 4π 3 ∇ 2 u 0 (x 0 )( n, n)∇ 2 v 0 (x 0 )( n, n) (11)
Remark 2.2 We will check numerically that the maximum magnitude of I(x 0 ) is reached when n is perpendicular to the filament direction.

NUMERICAL ALGORITHM AND EXPERIMENTAL RESULTS

In this section we describe the numerical method to perform the TG [START_REF] Liew | Application of wavelet theory for crack identification in structures[END_REF] and give some experiment results.

Algorithm

To compute the TG given in [START_REF] Liew | Application of wavelet theory for crack identification in structures[END_REF] we need to compute the solution u 0 and v 0 given by ( 3) and ( 6) for ǫ = 0. To do that, we consider that Ω is the unit square and we extend f by symmetry with respect to the boundary of Ω and by periodicity.

Then we can compute u 0 and v 0 by DCT (Discrete Cosine Transform). The computation time is O(N log(N )) where N denotes the number of pixels (or voxels). If we represent the blur K by a convolution operator we get the solutions in the Fourier domain :

u 0 = K ⋆ f α|ν| 4 + K ⋆ K v 0 = 2 K ⋆ K u 0 -K ⋆ f α|ν| 4 + K ⋆ K ( 12
)
where ν is the discrete Fourier variable and K, K ⋆ , f , u 0 and v 0 the array of Fourier coefficients. The discrete solutions u 0 and v 0 are then obtained from (12) by using the IDCT (Inverse DCT). We compute the Hessian matrices ∇ 2 u 0 and ∇ 2 v 0 at each point of the meshgrid by convolution with derivative filters; finally we approximate (10) by writting n in polar or spherical coordinates and by taking the discrete maximum. All calculus are implemented in Matlab 7.5.0 and the experiments are performed on a computer equipped with a processor Intel Core 1.9 GHz.

Experimental results in 2D

In this section we will compare the TG with an indicator based on the Hessian matrix eigenvalues of a Gaussian convolution of the image and used in [START_REF] Steger | An unbiased detector of curvilinear structures[END_REF]. We will denote by I Hes this indicator. The TG given in [START_REF] Auroux | Image restoration and classification by topological asymptotic expansion[END_REF] associated to F (u) = |∇u| 2 is denoted by I Lap , the TG performed in (10) by I Bilap . Figure 1 displays an image containing both an edge and a filament. First we see that I Lap detects more the edge than the filament while I Bilap (and also the Hessian indicator but less clearly) detects more the filament than the edge. In Figure 2 we display both I Hes and I Bilap for a very noisy image containing a filament, and we see that I Hes is less robust with respect to noise than the TG I Bilap . On Figure 3 we compare I Bilap , I Lap and I Hes on a satellite image of a road network. We can see that I Bilap detects the center of the roads while I Lap is high on the edges of the image and on the border of the roads. I Hes is sensitive to texture and is high both on the center of the roads and on edges.

Experimental results in 3D

Here we illustrate the fact that the method can easily be extended to 3D imaging with the TG [START_REF] Aubert | A topological gradient based model for the detection of fine structures in 2D images[END_REF]. The DCT is again

CONCLUSION

In this paper we have proposed a new detector (the topological gradient) of fine structures in 2D or 3D imaging. The main qualities of this detector are its simplicity and rapidity. Moreover it is quite efficient and robust both in the case of noisy and blurred images. More extensive simulations and comparisons on real images as well as results for perforated domains can be found in [START_REF] Drogoul | Numerical analysis of the topological gradient method for fourth order models and applications to the detection of fine structures in imaging[END_REF][START_REF] Aubert | A topological gradient based model for the detection of fine structures in 2D images[END_REF].

Fig. 1 .

 1 Fig. 1. Comparison of I Bilap (α = 10 -4 ) with the indicators I Lap (α = 10 2 ) and I Hes (σ = 5/4) on a simple non noisy image

Figure 4

 4 displays the TG I Bilap and the indicator I Hes for a blurred and noisy image. The result is clear : I Hes is inefficient in this case while the TG always detects the filament. The dimension of images are : 49×49 for Figure 1, 132×112 for Figure 2, 493 × 351 for Figure 3 and 193 × 165 for Figure 4. The computation times for Figures 1, 2, 3 and 4 are about 0.2 sec.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Comparison of I Bilap (α = 0.5) with the indicator I Hes (σ = 7/4) on a noisy image (PSNR=14dB)

Fig. 4 .Fig. 5 .

 45 Fig. 4. Comparison of I Bilap (α = 0.7) with the indicator I Hes (σ = 9/4) on a motion blurred (θ = 90 o , len=7) and Gaussian noisy image (PSNR=16dB)


Nicolas Courty

Rémi Flamary

Devis Tuia

Domain adaptation with regularized optimal transport

We present a new and original method to solve the domain adaptation problem using optimal transport. By searching for the best transportation plan between the probability distribution functions of a source and a target domain, a non-linear and invertible transformation of the learning samples can be estimated. Any standard machine learning method can then be applied on the transformed set, which makes our method very generic. We propose a new optimal transport algorithm that incorporates label information in the optimization: this is achieved by combining an efficient matrix scaling technique together with a majoration of a non-convex regularization term. By using the proposed optimal transport with label regularization, we obtain significant increase in performance compared to the original transport solution. The proposed algorithm is computationally efficient and effective, as illustrated by its evaluation on a toy example and a challenging real life vision dataset, against which it achieves competitive results with respect to state-of-theart methods.

Introduction

While most learning methods assume that the test data X t = (x t i ) i=1,...,Nt , x i ∈ R d and the training data X s = (x s i ) i=1,...,Ns are generated from the same distributions µ t = P(X t ) and µ s = P(X s ), real life data often exhibit different behaviors. Many works study the generalization capabilities of a classifier allowing to transfer knowledge from a labeled source domain to an unlabeled target domain: this situation is referred to as transductive transfer learning [1]. In our work, we assume that the source and target domains are by nature different, which is usually referred to as domain adaptation. In the classification problem, the training data are usually associated with labels corresponding to C different classes. We consider the case where only the training data are associated with a label Y s = (y i ) i=1,...,Ns , y i ∈ {1, . . . , C}, yielding an unsupervised domain adaptation problem, since no labelled data is available in the target domain. In this acceptation, the training (resp. testing) domain is usually referred to as source (resp. target) distribution.

Domain adaptation methods seek to compensate for inter domain differences by exploiting the similarities between the two distributions. This compensation is usually performed by reweighing the contribution of each samples in the learning process (e.g. [START_REF] Sugiyama | Direct importance estimation with model selection and its application to covariate shift adaptation[END_REF]) or by means of a global data transformation that aligns the two distributions in some common feature space (e.g. [START_REF] Zhang | Covariate shift in Hilbert space: A solution via surrogate kernels[END_REF]). Our work departs from these previous works by assuming that there exists a non-rigid transformation of the distribution that can account for the non-linear transformations occurring between the source and target domains. This transformation is conveniently expressed as a transportation of the underlying probability distribution functions thanks to optimal transport (OT). The OT problem has first been introduced by the French mathematician Gaspard Monge in the middle of the 19th century as the way to find a minimal effort solution to the transport of a given mass of dirt into a given hole. The problem reappeared later in the work of Kantorovitch [START_REF] Kantorovich | On the translocation of masses[END_REF], and found recently surprising new developments as a polyvalent tool for several fundamental problems [START_REF] Villani | Optimal transport: old and new[END_REF]. In the domain of machine learning, OT has been recently used for computing distances between histograms [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transportation[END_REF] or label propagation in graphs [START_REF] Solomon | Wasserstein propagation for semi-supervised learning[END_REF].

Contributions Our contributions are twofold: i) First, we show how to transpose the optimal transport problem to the domain adaptation problem, and we propose experimental validations of this idea. To the best of our knowledge, this is the first time that optimal transport is considered in the domain adaptation setting. ii) Second, we propose an elegant group-based regularization for integrating label information, which has the effect of regularizing the transport by adding inter-class penalties. The resulting algorithm exploits a proven efficient optimization approaches and will benefit from any advances in this domain. The proposed optimal transport with label regularization (OT-reglab) allows to achieve competitive state-of-the-art results on challenging datasets.

Related Work

Two main strategies have been considered to tackle the domain adaptation problem: on the one hand, there are approaches considering the transfer of instances, mostly via sample re-weighting schemes based on density ratios between the source and target domains [START_REF] Sugiyama | Direct importance estimation with model selection and its application to covariate shift adaptation[END_REF][START_REF] Bickel | Discriminative learning for differing training and test distributions[END_REF]. By doing so, authors compare the data distributions in the input space and try to make them more similar by weighting the samples in the source domain.

On the other hand, many works have considered finding a common feature representation for the two (or more) domains, or a latent space, where a classifier using only the labeled samples from the source domain generalize well on the target domains [START_REF] Daumé | Frustratingly easy domain adaptation[END_REF][START_REF] Pan | Domain adaptation via transfer component analysis[END_REF]. The representation transfer can be performed by matching the means of the domains in the feature space [START_REF] Pan | Domain adaptation via transfer component analysis[END_REF], aligning the domains by their correlations [START_REF] Kumar | Generalized multiview analysis: A discriminative latent space[END_REF] or by using pairwise constraints [START_REF] Wang | Manifold alignment without correspondence[END_REF]. In most of these works, the common latent space is found via feature extraction, where the dimensions retained summarize the information common to the domains. In computer vision, methods exploiting a gradual alignment of sets of eigenvectors have been proposed: in [START_REF] Gopalan | Domain adaptation for object recognition: An unsupervised approach[END_REF], authors start from the hypothesis that domain adaptation can be better approached if comparing gradual distortions and therefore use intermediary projections of both domains along the Grassmannian geodesic connecting the source and target observed eigenvectors. In [START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF][START_REF] Zheng | A grassmann manifold-based domain adaptation approach[END_REF], authors propose to obtain all sets of transformed intermediary domains by using a geodesic-flow kernel instead of sampling a fixed number of projections along the geodesic. While these methods have the advantage of providing easily computable out-of-sample extensions (by projecting unseen samples onto the latent space eigenvectors), the transformation defined is global and applied the same way to the whole target domain.

An approach combining the two logics is found in [START_REF] Zhang | Covariate shift in Hilbert space: A solution via surrogate kernels[END_REF], where authors extend the sample re-weighing reasoning to similarity of the distributions in the feature space by the use of surrogate kernels. By doing so, a linear transformation of the domains is found, but, as for the feature representation approaches above, it is the same for all samples transferred.

Our proposition strongly differs from those reviewed above, as it defines a local transportation plan for each sample in the source domain. In this sense, the domain adaptation problem can be seen as a graph matching problem for all samples to be transported, where their final coordinates are found by mapping the source samples to coordinates matching the marginal distribution of the target domain. In the authors knowledge, this is the first attempt to use optimal transportation theory in domain adaptation problem

Optimal transportation

In this Section, we introduce the original formulation of optimal transport through the Monge-Kantorovitch problem and its discrete formulation. Then, regularized versions of the optimal transport are exposed.

The Monge-Kantorovitch problem and Wasserstein space

Let us first consider two domains Ω 1 and Ω 2 (in the following, we will assume without further indication that Ω 1 = Ω 2 = R d ). Let P(Ω i ) be the set of all the probability measures over Ω i . Let µ ∈ P(Ω 1 ), and T be an application from Ω 1 → Ω 2 . The image measure of µ by T, noted T#µ, is a probability measure over Ω 2 which verifies:

T#µ(y) = µ(T -1 (y)), ∀y ∈ Ω 2 . (1)
Let µ s = P(Ω 1 ) and µ t = P(Ω 2 ) be two probability measures from the two domains. T is said to be a transport if T#µ s = µ t . The cost associated to this transport is

C(T) = Ω1 c(x, T(x))dµ(x), (2)
where the cost function c : Ω 1 × Ω 2 → R + can be understood as a regular distance function, but also as the energy required to move a mass µ(x) from x to y. It is now possible to define the optimal transport T 0 as the solution of the following minimization problem:

T 0 = arg min T Ω1 c(x, T(x))dµ(x), s.t. T#µ s = µ t (3)
which is the original Monge transportation problem. The equivalent Kantorovitch formulation of the optimal transport [START_REF] Kantorovich | On the translocation of masses[END_REF] seeks for a probabilistic coupling γ ∈ P(Ω 1 × Ω 2 ) between Ω 1 and Ω 2 :

γ 0 = arg min γ Ω1×Ω2 c(x, y)dγ(x, y), s.t. P Ω1 #γ = µ s , P Ω2 #γ = µ t , (4)
where P Ωi is the projection over Ω i . In this formulation, γ can be understood as a joint probability measure with marginals µ s and µ t . γ 0 is the unique solution to the optimal transport problem. It allows to define the Wasserstein distance between µ s and µ t as:

W 2 (µ s , µ t ) = inf γ Ω1×Ω2 c(x, y)dγ(x, y), s.t. P Ω1 #γ = µ s , P Ω2 #γ = µ t , (5)
This distance, also known as the Earth Mover Distance in computer vision community [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF], defines a metric over the space of integrable squared probability measure.

Optimal transport of discrete distributions

Usually one does not have a direct access to µ s or µ t but rather to collections of samples from those distributions. It is then straightforward to adapt the optimal transport problem to the discrete case. The two distributions can be written as

µ s = ns i=1 p s i δ x s i , µ t = nt i=1 p t i δ x t i ( 6
)
where δ xi is the Dirac at location x i ∈ R d . p s i and p t i are probability masses associated to the i-th sample, and belong to the probability simplex, i.e. ns i=1 p s i = nt i=1 p t i = 1. The set of probabilistic coupling between those two distributions is then the set of doubly stochastic matrices P defined as

P = γ ∈ (R + ) ns×nt | γ1 nt = µ s , γ T 1 ns = µ t (7)
where 1 d is a d-dimensional vector of ones. The Kantorovitch formulation of the optimal transport [START_REF] Kantorovich | On the translocation of masses[END_REF] reads:

γ 0 = arg min γ∈P γ, C F (8)
where ., . F is the Frobenius dot product and C ≥ 0 is the cost function matrix of term C(i, j) related to the energy needed to move a probability mass from x s i to x t j . This cost can be chosen for instance as the Euclidian distance between the two locations, i.e. C(i, j) = ||x s ix t j || 2 , but other types of metric could be considered, such as Riemannian distances over a manifold [START_REF] Villani | Optimal transport: old and new[END_REF].

Remark 1. When n s = n t = n and when ∀i, j p s i = p t j = 1/n, the γ 0 is simply a permutation matrix Remark 2. In the general case, it can be shown that γ 0 is a sparse matrix with at most n s + n t -1 non zero entries (rank of constraints matrix). This problem can be solved by linear programming, with combinatorial algorithms such as the simplex methods and its network variants (transport simplex, network simplex, etc.). Yet, the computational complexity was shown to be O(n 2 ) in practical situations [START_REF] Bonneel | Displacement interpolation using lagrangian mass transport[END_REF] for the network simplex (while being O(n 3 ) in theory) which leverages the utility of the method to handle big data. However, the recent regularization of Cuturi [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transportation[END_REF] allows a very fast transport computation as discussed in the next Section.

Regularized optimal transport

When the target and source distributions are high-dimensional, or even in presence of numerous outliers, the optimal transportation plan may exhibit some irregularities, and lead to incorrect transport of points. While it is always possible to enforce a posteriori a given regularity in the transport result, a more theoretically convincing solution is to regularize the transport by relaxing some of the constraints in the problem formulation of Eq.( 8). This possibility has been explored in recent papers [START_REF] Ferradans | Regularized discrete optimal transport[END_REF][START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transportation[END_REF].

In [START_REF] Ferradans | Regularized discrete optimal transport[END_REF], Ferradans and colleagues have explored the possibility of relaxing the mass conservation constraints of the transport, i.e. slightly distorting the marginals of the coupling γ 0 . Technically, this boils down to solving the same minimization problem but with inequality constraints on the marginals in Eq.( 7). As a result, elements of the source and target distributions can remain still. Yet, one major problem of this approach is that it converts the original linear program into more computationally demanding optimizations impractical for large sets.

In a recent paper [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transportation[END_REF], Cuturi proposes to regularize the expression of the transport by the entropy of the probabilistic coupling. The regularized version of the transport γ λ 0 is then the solution of the following minimization problem:

γ λ 0 = arg min γ∈P γ, C F - 1 λ h(γ), (9)
where h(γ) =i,j γ(i, j) log γ(i, j) computes the entropy of γ. The intuition behind this form of regularization is the following: since most of the elements of γ 0 should be zero with high probability, one can look for a smoother version of the transport by relaxing this sparsity through an entropy term. As a result, and contrary to the previous approach, more couplings with non-nul weights are allowed, leading to a denser coupling between the distributions. An appealing result of this formulation is the possibility to derive a computationally very efficient algorithm, which uses the scaling matrix approach of Sinkhorn-Knopp [START_REF] Knight | The sinkhorn-knopp algorithm: Convergence and applications[END_REF]. The optimal regularized transportation plan is found by iteratively computing two scaling vectors u and v such that:

γ λ 0 = diag(u) exp(-λC)diag(v), (10)
where the exponential exp(.) operator should be understood element-wise. Note that while these regularizations allow the inclusion of additional priors in the optimization problem, they do not take into account the fact that the elements of the source distribution belong to different classes. This idea is the core of our regularization strategy.

Domain Adaptation with Label Regularized Optimal Transport

From the definitions above, the use of optimal transport for domain adaptation is rather straightforward: by computing the optimal transport from the source distribution µ s to the target distribution µ t , one defines a transformation of the source domain to the target domain. This transformation can be used to adapt the training distribution by means of a simple interpolation. Once the source labeled samples have been transported, any classifier can be used to predict in the target domain. In this section, we present our optimal transport with label regularization algorithm (OT-labreg) and derive a new efficient algorithm to solve the problem. We finally discuss how to interpolate the training set from this regularized transport.

Regularizing the transport with class labels

Optimal transport aims at minimizing a transport cost linked to a metric between distributions. It does not include any information about the particular nature of the elements of the source domain (i.e. the fact that those samples belong to different classes). However, this information is generally available, as labeled samples are used in the classification step following adaptation. Our proposition to take advantage of label information is to penalize couplings that match together samples with different labels. This is illustrated in Figure 1.c, where one can see that samples belonging to the same classes are only associated to points associated to the same class, contrarily to the standard and regularized versions of the transport (Figures 1.a and 1.b ).

Principles of the label regularization Over each column of γ, we want to concentrate the transport information on elements of the same class c. This is usually done by using ℓ p -ℓ q mixed-norm regularization, among which the ℓ 1 -ℓ 2 known as as group-lasso is a favorite. The main idea is that, even if we do not know the class of the target distribution, we can promote group sparsity in the columns of γ such that a given target point will be associated with only one of the classes. Promoting group sparsity leads to a new term in the cost function [START_REF] Daumé | Frustratingly easy domain adaptation[END_REF], which now reads: where I c contains the index of the lines such that the class of the element is c, γ(I c , j) is a vector containing coefficients of the jth column of γ associated to class c and || • || p q denotes the ℓ q norm to the power of p. η is a regularization parameter that weights the impact of the supervised regularization.

γ 0 = arg min γ∈P γ, C F - 1 λ h(γ) + η j c ||γ(I c , j)|| p q , (11)
The choice of the p, q parameters is particularly sensitive. For p ≥ 1 and q ≥ 1 the regularization term is convex. The parameters p = 1, q = 2 lead to the classical group-lasso that is used, for instance, for joint features selection in multitask learning. The main problem of using the group-lasso in this case is that it makes the optimization problem much more difficult. Indeed, when using an ℓ 2 norm in the objective function, the efficient optimization procedure proposed in [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transportation[END_REF] cannot be used anymore. Moreover there is no particular reason to choose the ℓ 2 norm for regularizing coefficients of a transport matrix. Those coefficients being all positive and associated to probabilities, we propose to use q = 1 that will basically sum the probabilities in the groups. When q = 1, one needs to carefully chose the p coefficient in order to promote group sparsity. In this work we propose to use p = 1/2 < 1. This parameter is a common choice for promoting sparsity, as the square root is non-differentiable in zero and has been used recently for promoting non-grouped sparsity in compressed sensing [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF].An additional advantage of our proposal is that, despite the fact that the proposed regularization is non-convex, a simple approach known as reweighted ℓ 1 can be performed for its optimization, as detailed below.

Majoration Minimization strategy

The optimization problem with a ℓ pℓ 1 regularization boils down to optimizing

γ 0 = arg min γ∈P J(γ) + ηΩ(γ), (12)
with

J(γ) = γ, C F -1 λ h(γ)
and Ω(γ) = j c ||γ(I c , j)|| p 1 . We want to be able to use the optimization in [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transportation[END_REF] to solve the left term, as it is very efficient. First, note that the regularization term can be reformulated as

Ω(γ) = j c g(||γ(I c , j)|| 1 ) (13)
where g(•) = (•) p is a concave function of a positive variable (∀γ ≥ 0). A classical approach to address this problem is to perform what is called Majorizationminimization [START_REF] Hunter | A Tutorial on MM Algorithms[END_REF]. This can be done because the ℓ pℓ 1 regularization term is concave in the positive orthant as illustrated in the left part of figure 2. It is clear from this Figure that the surface can be majored by an hyperplane. For a given group of variable, one can use the concavity of g to majorize it around a given vector ŵ > 0 g(w) ≤ g( ŵ 1 ) + ∇g( ŵ 1 ) ⊤ (wŵ)

with ∇g( ŵ 1 ) = p( ŵ 1 ) p-1 for ŵ > 0. An illustration of the majoration of g(•) can be seen in the right part of Figure 2. For each group, the regularization term can be majorized by a linear approximation. In other words, for a fixed γ

Ω(γ) ≤ Ω(γ) = γ, G F + cst (15)
where the matrix G has components

G(I c , j) = p( γ(I c , j) + ǫ) p-1 , ∀c, j (16)
Note that we added a small ǫ > 0 that helps avoiding numerical instabilities, as discussed in [START_REF] Candes | Enhancing sparsity by reweighted l1 minimization[END_REF]. Finally, solving problem (11) can be performed by iterating the two steps illustrated in Algorithm 1. This iterative algorithm is of particular interest in our case as it consists in iteratively using an efficient Sinkhorn-Knopp matrix scaling approach. Moreover this kind of MM algorithm is known to converge in a small number of iterations.

Algorithm 1 Majoration Minimization for ℓ p -ℓ 1 regularized Optimal Transport

Initialize G = 0 Initialize C0 as in Equation ( 8) [START_REF] Daumé | Frustratingly easy domain adaptation[END_REF] with C G ← Update G with Equation ( 16) until Convergence

repeat C ← C0 + G γ ← Solve problem

Interpolation of the source domain

Once the transport γ 0 has been defined using either Equations ( 8), ( 9) or [START_REF] Kumar | Generalized multiview analysis: A discriminative latent space[END_REF], the source samples must be transported in the target domain using their transportation plan. One can seek the interpolation of the two distributions by following the geodesics of the Wasserstein metric [START_REF] Villani | Optimal transport: old and new[END_REF] (parameterized by t). This allows to define a new distribution µ t such that:

µ t = arg min µ (1 -t)W 2 (µ s , µ) 2 + tW 2 (µ t , µ) 2 . ( 17
)
One can show that this distribution is:

µ t = i,j γ 0 (i, j)δ (1-t)x s i +tx t j . (18)
In our approach, we suggest to compute directly the image of the source samples as the result of this transport, i.e. for t = 1. Those images can be expressed through γ 0 as barycenters of the target samples. Let T γ 0 : R d → R d be the mapping induced by the optimal transport coupling. This map transforms the source elements X s in a target domain dependent version Xs . The mapping T γ 0 can be conveniently expressed as:

Xs = T γ 0 (X s ) = diag((γ 0 1 nt ) -1 )γ 0 X t . (19)
We note that T γ 0 is fully invertible and can be also used to compute an adaptation from the target domain to the source domain by observing that T -1 γ 0 = T γ T 0 . Let us finally remark that similar interpolation methods were used in the domain of color transfer [START_REF] Ferradans | Regularized discrete optimal transport[END_REF].

Experimental validation

In this Section, we validate the proposed algorithm in two domain adaptation examples. On the first one, we study the behavior of our approach on a simple toy dataset. The second one considers a challenging computer vision dataset, used for a comparison with state-of-the-art methods. In every experiment, the original optimal transport (OT-ori) is computed with a network simplex approach [START_REF] Bonneel | Displacement interpolation using lagrangian mass transport[END_REF].

The Sinkhorn transport, which corresponds to the regularized version of the optimal transport (OT-reg) described in Section 3.3, was implemented following the algorithm proposed in [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transportation[END_REF]. Our approach, OT-reglab, follows the Algorithm 1. As expected, these last two methods are generally one order of magnitude faster than the network simplex approach. As for the choice of the weights of Eq. ( 6), the problem can be cast as an estimation of a probability mass function of a discrete variable on the sample space of the source and target distributions. A direct and reasonable choice is to take an uniform distribution, i.e. p s i = 1 ns and p s t = 1 nt . This choice gives the same value for every samples in the two discrete distributions. Alternatively, one can seek to strengthen the weights of samples that are in a high density region, and lower weights for samples in low density regions. This way, outliers should be associated with lower masses. A possible solution relies on a discrete variant of the Nadaraya-Watson estimator [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF] where one enforces the sum-to-1 property:

p s i = ns j=1 k σ (x s i , x s j ) ns j=1 ns i=1 k σ (x s i , x s j ) (20)
where k σ (•, •) is a gaussian kernel of bandwidth σ. The drawback of such an estimator is that it adds an hyper parameter to the method. Yet, while standard approaches [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF] can be used to estimate this parameter, we observed in our experiments, and for large number of samples, that this parameter exerts little influence over the final result (less than a standard deviation) for a large range of values.

Toy dataset

In this first experiment, the behavior of the optimal transport is examined on a simple two-dimensional dataset. We consider a two-class distribution by sampling independently for each class c1 and c2 following the normal distributions N s 1 and N s 2 . The set of all those samples constitute the source domain. The target domain samples are then obtained by sampling the mixture N t 1 + N t 2 . The target distributions N t i , (i = 1, 2) are deduced from N s i , (i = 1, 2) by changing both the scale and translating the distribution mean. The produced domain transformation is thus non-linear and cannot be expressed by a simple 2D transformation of the input space. This makes the problem particularly interesting with respect to our initial assumptions on the nature of the domain change. We then sample randomly from these distributions n s 1 ,n s 2 ,n t 1 and n t 2 samples from N s 1 , N s 2 , N t 1 and N t 2 to form the corresponding learning and test sets. An illustration of this toy dataset is given in Figure 3.a for n s 1 + n s 2 = 100 samples in the source distribution (red and white circles) and n t 1 + n t 2 = 200 samples in the target one (blue crosses). Note that the size of the points in the Figure is proportional to its weight p i and reflects the density of the distribution.

Figure 3.b presents the result of the optimal transport OT-ori coupling as a set of non-nul connections (red and black arcs) between the source and the target distributions. The color of those connections is related to the magnitude of the coupling (up to a global scaling factor). As expected, the coupling matrix γ 0 contains less than 100 + 200 -1 = 299 non-nul entries. One can see that some white and red elements are clearly misled by the transport, but the overall adaptation remains coherent with the test distribution. Figure 4 illustrates the results obtained on this dataset with the regularized versions of the transport OT-reg and OT-reglab for a regularization parameter value of λ = 1. The γ 0 matrix of OT-reg, on the left of the first row of Figure 4, is indeed sparse, but much less than the corresponding one in OT-ori. This can be assessed by comparing the denser connections issued from OT-reg (left panel of the second row of Figure 4) with respect to those observed for OTori (right panel of Figure 3). In the proposed OT-reglab (right column of Figure 4), the sparsity is clearly enforced per class (the rows of the coupling matrix are sorted by class), which yields a sparser coupling matrix with block structure. In the last row of Figure 4 we show the result of the adaptation of the source distribution following the procedure described in Section 4.3. Two additional interesting behaviors are observed in the regions highlighted by red squares, where some of the incoherencies observed in OT-ori and OT-reg of the transport are resolved by the label regularization proposed with OT-reglab. Classification measures. We now consider performances of a classifier trained on the source samples adapted to the target distribution. In those experiments, we use a SVM classifier with a Gaussian kernel. The hyperparameters of the classifier are computed for each trial by a 2-fold cross validation over a grid of potential values. For every setting considered, the data generation / adaptation / classification was conducted 20 times to leverage the importance of the sampling. When informative, we provide the standard deviation of the result. In the first experiment, we examine the importance of the regularization parameter λ over the overall classification accuracy (Figure 5.a). In this case, we set n s 1 = n s 2 = n t 1 = n t 2 = 100. We confirm that the use of the transport for domain adaptation increases the performances significantly (by 8%) over a classification conducted directly with the source distribution as learning set. When varying the λ regularization parameter and using OT-reglab, another very significant increase is achieved (up to 25% for λ = 0.04), which demonstrates the relevance of our transport regularization. In the second experiment, we set n s 1 = n s 2 = 100 and we increase the number of elements in the source target n t 1 and n t 2 equivalently. For this experiment and for the next one, λ is set by a standard cross-validation method. In this case, the standard deviation is omitted as it is constant over the experiments and no informative. One can observe that the performances of the classification are i) consistent with the first experiment and ii) constant over the volume of samples in the target domain as long as the proportions are conserved. In the third experiment, we set n t 1 and n t 2 to the value of 100 samples each and we vary the proportion of the classes through a parameter p ∈ [0, 1] with n s 1 = p * 100 and n s 2 = (1p) * 100. This parameter allows to control the proportion of elements in class 1 and in class 2 in the source distribution. As shown in Figure 5.c, the best result is achieved when the proportion of each class samples is similar in the source and target distributions (at 50%). This somehow highlights one limit of the method: the mass equivalent to each class should match in proportions for both distributions to get the best adaption result. Nevertheless, we can see from Figure 5.c that a variation of ±15% between the source and target distribution still leads to significant performance improvements.

Visual adaptation dataset

We now evaluate our method on a challenging real world dataset coming from the computer vision community. The objective is now a visual recognition task of several categories of objects, studied in the following papers [START_REF] Saenko | Adapting visual category models to new domains[END_REF][START_REF] Gopalan | Domain adaptation for object recognition: An unsupervised approach[END_REF][START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF][START_REF] Zheng | A grassmann manifold-based domain adaptation approach[END_REF]. The dataset contains images coming from four different domains: Amazon (online merchant), the Caltech-256 image collection [START_REF] Griffin | Caltech-256 Object Category Dataset[END_REF], Webcam (images taken from a webcam) and DSLR (images taken from a high resolution digital SLR camera). Those domains are respectively noted in the remainder as A, C, W and D. A feature extraction method is used to preprocess those images; it namely consists in computing SURF descriptors [START_REF] Saenko | Adapting visual category models to new domains[END_REF], which allows to transform each image into a 800 bins histogram, which are then subsequently normalized and reduced to standard scores. We followed the experimental protocol exposed in [START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF]: each dataset is considered in turn as the source domain and used to predict the others. Within those datasets, 10 classes of interest are extracted. The source domain are formed by picking 20 elements per class for domains A,C and W, and 8 for D. The training set is then formed by adapting these samples to the target domain. The latter is composed of all the elements in the test domain. The classification is conducted using a 1-Nearest Neighbor classifier, which avoids cross-validation of hyper-parameters. As for the toy example above, we repeat each experiment 20 times and report the overall classification accuracy and the associated standard deviation. We compare the results of the three transport models (OT-ori, OT-reg and OT-reglab) against both a classification conducted without adaptation (no adapt.) and 3 state-of-the-art methods: 1) the surrogate kernel approach (SuK), which in [START_REF] Zhang | Covariate shift in Hilbert space: A solution via surrogate kernels[END_REF] was shown to outperform both the Transfer Component Analysis method [START_REF] Pan | Domain adaptation via transfer component analysis[END_REF] and the reweighing scheme of [START_REF] Sugiyama | Direct importance estimation with model selection and its application to covariate shift adaptation[END_REF]; 2) the (SGF) method proposed in [START_REF] Gopalan | Domain adaptation for object recognition: An unsupervised approach[END_REF] and 3) the Geodesic Flow Kernel (GFK) approach proposed in [START_REF] Gong | Geodesic flow kernel for unsupervised domain adaptation[END_REF]. Note that this last method can also efficiently incorporate label information: therefore we make a distinctions between methods, which do not incorporate label information (no adapt, SuK, SGF, GFK, OT-ori and OT-reg) and those that do (GFK-lab and OT-reglab). For each setting we used the recommended parameters to tune the competing methods. Results are reported in Table . 1.

When no label information is used, (OT-reg) usually performs best. In some cases (notably when considering the adaptation from (W→A or D→W), it can even surpass the (GFK-lab) method, which uses labels information.OT-ori usually enhances the result obtained without adaptation, but remains less efficient than the competing methods (except in the case of W→A where it surpasses SGF and SuK. Among all the methods, OT-reglab usually performs best, and with a significant increase in the classification performances for some cases (W→C or D→W). Yet, our method does not reach state-of-the-art performance in two cases: A→C and D→A. Finally, the overall mean value (last line of the table) shows a consistent increase of the performances with the proposed OT-reglab, which outperforms in average GFK-lab by 2%. Also note that the regularized unsupervised version OT-reg outperforms all the competing methods by at least 3%.

Conclusion and discussion

We have presented in this paper a new method for unsupervised domain adaptation based on the optimal transport of discrete distributions from a source to a target domain. While the classical optimal transport provide satisfying re-sults, it fails in some cases to provide state-of-the-art performances in the tested classification approaches. We proposed to regularize the transport by relaxing some sparsity constraints in the probabilistic coupling of the source and target distributions, and to incorporate the label information by penalizing couplings that mix samples issued from different classes. This was made possible by a Majoration Minimization strategy that exploits a ℓ pℓ 1 norm, which promotes sparsity among the different classes. The corresponding algorithm is fast, and allows to work efficiently with sets of several thousand samples. With this regularization, competitive results were achieved on challenging domain adaptation datasets thanks to the ability of our approach to express both class relationship and non-linear transformations of the domains.

Possible improvements of our work are numerous, and include: i) extension to a multi-domain setting, by finding simultaneously the best minimal transport among several domains, ii) extension to semi-supervised problems, where several unlabeled samples in the source domain, or labelled samples in the target domain are also available. In this last case, the group sparsity constraint should not only operate over the columns but also the lines of the coupling matrix, which makes the underlying optimization problem challenging. iii) Definition of the transport in a RKHS, in order to exploit the manifold structure of the data.

Fig. 1 .

 1 Fig. 1. Illustration of the transport for two simple distributions depicted in the image. The colored disks represent 3 different classes. The transport solution is depicted as blue lines whose thickness relate to the strength of the coupling. (a) Solution of the original optimal transport solution (OT-ori); (b) using the Sinkhorn transport (OTreg [6]); (c) using our class-wise regularization term (OT-reglab).

1Fig. 2 .

 2 Fig. 2. Illustration of the regularization term loss for a 2D group (left). Illustration of the convexity of g(•) and its linear majoration (right)

Fig. 3 .

 3 Fig.3. Illustration of the transport OT-ori on a simple toy dataset. The initial distributions are depicted in the right image (a) The source distribution is depicted in white and red for respectively class 1 and 2, the target distributions are in blue. In image (b), we show the optimal transport couplings, depicted as links colored with respect to the source class label.

Fig. 4 .

 4 Fig. 4. Comparisons of two versions of the regularized transport: Sinkhorn transport (OT-reg, left column) and Sinkhorn transport with the label regularization (OTreglab, right column). The first row shows the transport coupling matrices γ 0 , the second row their equivalent graphical representations, with connections colored by the source node label. The third row is the adaptation of the source samples induced by γ 0 using Equation (19).

Fig. 5 .

 5 Fig. 5. Classification results for the toy dataset example: (a) influence of the regularization parameter λ; (b) influence of the proportions of samples between the source and the target distributions; (c) influence of the balance of classes on the overall performance of the adaptation.

Acknowledgements We thank the anonymous reviewers for their critics and suggestions. This work has been partially funded by a visiting professor grant from EPFL, by the French ANR under reference ANR-13-JS02-0005-01 (Asterix project). and by the Swiss National Science Foundation (grant 136827, http://p3.snf.ch/project-136827).


Christophe Chalons

Frédéric Coquel

Samuel Kokh

Nicole Spillane

Large Time-Step Numerical Scheme for the Seven-Equation Model of Compressible Two-Phase Flows

 for Euler equations. The main objective is to get a Courant-Friedrichs-Lewy (CFL) condition driven by (slow) contact waves instead of (fast) acoustic waves.

Introduction

We are interested in the computation of compressible two-phase flows with the so-called two-fluid two-pressure or seven-equation model. It was first proposed in Baer & Nunziato [START_REF] Baer | A two phase mixture theory for the deflagration to detonation transition in reactive granular materials[END_REF] and has since aroused more and more interest, see for instance Embid & Baer [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF], Stewart & Wendroff [START_REF] Stewart | Two-phase flow : models and methods[END_REF], Abgrall & Saurel [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF], Gallouët, Hérard & Seguin [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid twopressure approach[END_REF], Andrianov & Warnecke [START_REF] Andrianov | The Riemann problem for the Baer-Nunziato two-phase flow model[END_REF], Karni et al. [START_REF] Karni | Compressible two-phase flows by central and upwind schemes[END_REF] Schwendeman, Wahle & Kapila [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF], Munkejord [START_REF] Munkejord | Comparison of Roe-type methods for solving the two-fluid model with and without pressure relaxation[END_REF], Tokareva & Toro [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF], Ambroso, Chalons, Coquel & Galié [START_REF] Ambroso | Relaxation and numerical approximation of a two fluid two pressure diphasic model[END_REF], Ambroso, Chalons & Raviart [START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF], and the references therein. One of the main features of this model is that it is hyperbolic, at least in the context of subsonic flows. In particular, an interesting property is that the seven-equation model possesses seven real eigenvalues given by λ ± k (u) = u k ± c k , λ 0 k (u) = u k and λI (u) = uI , where u k denote the velocities of both phases k = 1, 2, c k the sound speeds, uI an interfacial velocity and u the vector of unknowns. However from a numerical point of view, the seven-equation model raises some issues. The first difficulty is related to the large size of the model and as a consequence to the Riemann problem that is difficult to solve, even approximately. The second difficulty comes from the presence of nonconservative products and more precisely the fact that the model cannot be equivalently recast in full conservative form. However, the nonconservative products naturally vanish when the void fractions α k are locally constant in space, and the model coincides in that case with two (decoupled) gas dynamics systems. This property will be used in the numerical strategy. Numerous papers are devoted to the numerical study of two-fluid two-pressure models, see again for instance [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid twopressure approach[END_REF], [START_REF] Andrianov | The Riemann problem for the Baer-Nunziato two-phase flow model[END_REF], [START_REF] Karni | Compressible two-phase flows by central and upwind schemes[END_REF], [START_REF] Schwendeman | The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow[END_REF], [START_REF] Munkejord | Comparison of Roe-type methods for solving the two-fluid model with and without pressure relaxation[END_REF], [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF], [START_REF] Ambroso | Relaxation and numerical approximation of a two fluid two pressure diphasic model[END_REF], [START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF] and the references therein. Many of the proposed methods are based on time-explicit, exact or approximate, Godunov-type methods (Roe or Roe-like scheme, HLL or HLLC scheme...). For stability reasons, the time steps ∆t involved in such methods are subject to a usual Courant-Friedrichs-Lewy (CFL) condition that reads max

k,u (|λ ± k (u)|, |λ 0 k (u)|, |λI (u)|) ∆t ≤ 0.5∆x,
where ∆x represents the space step. It is then clear that the definition of ∆t is driven by the fastest eigenvalues λ ± k (u), associated with the so-called acoustic waves. In many applications, like for instance in two-phase flows involved in nuclear reactors, the acoustic waves are not predominant physical phenomena. A CFL condition based on the most influent waves, the so-called contact waves associated with eigenvalues λ 0 k (u) and λI (u) would be more adapted. The idea is then to propose a timeimplicit treatment of the (fast) acoustic waves, in order to get rid of a too restrictive CFL condition, together with an explicit treatment of the (slow) contact waves in order to preserve accuracy. This was recently proposed in Coquel et al. [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] in the context of Euler equations, using a Lagrange-Projection approach. This approach is well-adapted as it naturally decouples the fast and slow waves in the Lagrange and Projection steps respectively. In this paper, we propose a first attempt to extend this approach to the seven-equation model. The idea is to operate a relevant operator splitting between the conservative and nonconservative parts of the original model, in order to make Euler systems for each phase appear. The latter parts are treated as in [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]. The nonconservative products are then discretized so as to maintain conservativity properties of the model on each partial mass, on the total momentum and total energy. Numerical results are proposed. We underline that this work is still in progress.

Governing equations

The model under consideration in this contribution reads as follows:

       ∂tα k + uI ∂xα k = 0, t > 0, x ∈ R, ∂tα k ̺ k + ∂xα k ̺ k u k = 0, ∂tα k ̺ k u k + ∂xα k (̺ k u 2 k + p k ) -pI∂xα k = 0, ∂tα k ̺ k e k + ∂xα k (̺ k e k + p k )u k -pI uI ∂xα k = 0, (1)
with k = 1, 2. Here, α k , ̺ k , u k , e k and p k denote the volume fraction, density, velocity, specific total energy and pressure of the phase k = 1, 2. The two phases are assumed to be non-miscible that is α1 + α2 = 1. The structure of ( 1) is the one of two gas dynamics systems for each phase, coupled with a transport equation on the void fraction α k at speed uI . We note that nonconservative products involving the interfacial pressure pI and velocity uI (to be precised later on) and the space derivative of the void fractions α k are present in the momentum and energy equations. These terms act as coupling terms in the evolution of the two phases. Source terms like external forces, pressure and velocity relaxations, dissipation, heat conduction, phase changes and heat exchanges between the two phases are not taken into account. Each phase is provided with an equation of state

p k = p k (̺ k , ε k ), where ε k = e k -u 2
k /2 is the specific internal energy. So far as the definitions of uI and pI are concerned, we follow [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid twopressure approach[END_REF] and first observe that the characteristic speeds of (1) are always real and given by uI , u k , u k ± c k , k = 1, 2, where c k denotes the speed of sound in phase k. System (1) turns out to be only weakly hyperbolic since there are not enough eigenvectors to span the entire space when uI = u k ± c k for some index k (resonance occurs). When ( 1) is hyperbolic, one can easily check that similarly to the classical gas dynamics equations, the characteristic fields associated with the eigenvalues u k ± c k are nonlinear while the one associated with u k is linearly degenerate. Regarding the characteristic field associated with uI , it is generally required to be linearly degenerate in practice. This property holds as soon as

uI = βu1 + (1 -β)u2, β = χα1̺1 χα1̺1 + (1 -χ)α2̺2 (2)
where χ ∈ [0, 1] is a constant (we refer to [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid twopressure approach[END_REF] for the details), which gives a natural definition for the interfacial velocity uI . The usual choices for χ are 0, 1/2 and 1. Regarding the interfacial pressure pI , we set pI = µp1 + (1 -µ)p2, µ ∈ [0, 1].The choice of the coefficient µ is not detailed here (see again [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid twopressure approach[END_REF]) but is related to the ability to provide the system with an entropy balance equation. Indeed, it can be proved that for a specific choice of µ, smooth solutions of (1) verify the conservation law ∂tη + ∂xq = 0, where (η, q) plays the role of a mathematical entropy pair.

A natural operator splitting

The starting point is to propose an equivalent form of (1) where the space derivatives of α k p k and α k p k u k are decomposed using a chain rule:

       ∂tα k + uI ∂xα k = 0, ∂tα k ̺ k + ∂xα k ̺ k u k = 0, ∂tα k ̺ k u k + ∂xα k ̺ k u 2 k + α k ∂xp k + (p k -pI )∂xα k = 0, ∂tα k ̺ k e k + ∂xα k ̺ k e k u k + α k ∂xp k u k + (p k u k -pI uI )∂xα k = 0. (3)
We then suggest to split (3) into two independent and quasi-classical gas dynamics equations (their Lagrangian forms will be seen to be classical), namely

       ∂tα k = 0, ∂tα k ̺ k + ∂xα k ̺ k u k = 0, ∂tα k ̺ k u k + ∂xα k ̺ k u 2 k + α k ∂xp k = 0, ∂tα k ̺ k e k + ∂xα k ̺ k e k u k + α k ∂xp k u k = 0, (4)
and into the following genuinely nonconservative system:

       ∂tα k + uI ∂xα k = 0, ∂tα k ̺ k = 0, ∂tα k ̺ k u k + (p k -pI )∂xα k = 0, ∂tα k ̺ k e k + (p k u k -pI uI )∂xα k = 0. (5)
This transformation aims at proposing in the next section an implicit-explicit Lagrange-Projection scheme similar to [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF], and at treating separately the nonconservative products. Note from now on that the overall algorithm will be conservative on the partial mass α k ̺ k , total momentum α1̺1u1 + α2̺2u2 and on the total energy α1̺1e1 + α2̺2e2, as it is expected from the original form (1) of the model.

Numerical approximation

This section is devoted to the discretization of (1), using ( 4) and ( 5). Let us introduce a time step ∆t > 0 and a space step ∆x > 0 that we assume to be constant for simplicity. We set λ = ∆t/∆x and define the mesh interfaces x j+1/2 = j∆x for j ∈ Z, and the intermediate times t n = n∆t for n ∈ N. In the sequel, u n j = (α1, u1, u2)

n j where (u k ) n j = (α k ̺ k , α k ̺ k u k , α k ̺ k e k
) n j denotes the approximate value of the unknowns at time t n and on the cell Cj =]x j-1/2 , x j+1/2 [. Implicit-explicit discretization of (4). We first recall that ( 4) is made of two independent quasi-classical gas dynamics systems, whose eigenvalues are given by u k ± c k , u k and 0. As already stated, our aim is to propose an implicit treatment of the fast waves u k ± c k , and an explicit treatment of u k . With this in mind, we follow [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] and adopt a Lagrange-Projection scheme, coupled with a pressure relaxation strategy that is well adapted to this purpose. A Lagrange-Projection splitting strategy applied to (4) amounts to introducing the Lagrangian system

       ∂tα k = 0, ∂tα k ̺ k + α k ̺ k ∂xu k = 0, ∂tα k ̺ k u k + α k ̺ k u k ∂xu k + α k ∂xp k = 0, ∂tα k ̺ k e k + α k ̺ k e k ∂xu k + α k ∂xp k u k = 0, or equivalently        ∂tα k = 0, ∂tτ k -τ k ∂xu k = 0, ∂tu k + τ k ∂xp k = 0, ∂te k + τ k ∂xp k u k = 0, (6)
with τ k = 1/ρ k , and the transport (or projection) system

       ∂tα k = 0, ∂tα k ̺ k + u k ∂xα k ̺ k = 0, ∂tα k ̺ k u k + u k ∂xα k ̺ k u k = 0, ∂tα k ̺ k e k + u k ∂xα k ̺ k e k = 0. (7)
We note that (6) coincides with two classical gas dynamics systems written in Lagrangian coordinates, the eigenvalues of which are given by ±c k and 0. This system is treated using a pressure relaxation approach that consists in introducing a linearized pressure π k (see for instance [START_REF] Chalons | Relaxation approximation of the Euler equations[END_REF] and especially the references therein), such that (π k ) n j = (p k ) n j , and in solving the partial differential system

           ∂tα k = 0, ∂tτ k -τ k ∂xu k = 0, ∂tu k + τ k ∂xπ k = 0, ∂tπ k + a 2 k τ k ∂xu k = 0, ∂te k + τ k ∂xπ k u k = 0, or equivalently            ∂tα k = 0, ∂tI k = 0, ∂tw + k + a k τ k ∂xw + k = 0, ∂tw - k -a k τ k ∂xw - k = 0, ∂te k + τ k ∂xπ k u k = 0, (8)
where the second one (top right) and the third one (bottom) involve several distinct waves. The scheme we propose here is denoted LP implicit and is compared with its explicit version (which amounts to replacing (9) by its time-explicit version) and the well-known Rusanov scheme (see [START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid twopressure approach[END_REF]). We observe that our approach is clearly less diffusive around the contact discontinuities since the CFL condition is well-adapted to the corresponding speed of propagation, but more diffusive around the acoustic waves since it is implicit. Table 1 gives for each test case the number of iterations needed to perform the computations. As expected, the gain is important when using the proposed implicit-explicit algorithm and the corresponding CFL restriction based on the material waves (instead of the acoustic waves as for the explicit scheme). A careful evaluation of the CPU cost necessitates an additional programming effort that has not been implemented yet.

w ± k = π k ± a k u k , I k = π k + a 2 k τ k , and a k is a constant satisfying the subcharacteristic condition a k > ρ k c k . A natural time-implicit discretization of (8) is                  (α k ) n+1= j = (α k ) n j , (I k ) n+1= j = (I k ) n j , (w + k ) n+1= j = (w + k ) n j -λ(τ k ) n j a k (w + k ) n+1= j -(w + k ) n+1= j-1 , (w - k ) n+1= j = (w - k ) n j + λ(τ k ) n j a k (w + k ) n+1= j+1 -(w + k ) n+1= j , (e k ) n+1= j = (e k ) n j -λ(τ k ) n j (π k u k ) n+1= j+1/2 -(π k u k ) n+1= j-1/2 , (9) with
(π k u k ) n+1= j+1/2 = (π k ) n+1= j+1/2 (u k ) n+1= j+1/2 and (π k ) n+1= j+1/2 = 1 2 (w + k ) n+1= j + (w - k ) n+1= j , (u k ) n+1= j+1/2 = 1 2a k (w + k ) n+1= j -(w - k ) n+1= j . 0

Table 1 :

 1 Number of time-iterations for each test case

	1			Rusanov		0.265						Rusanov
	0.9			LP explicit LP implicit		0.26					LP explicit LP implicit Rusanov(reference)
	0.8
						0.255
	0.7
	0.6					0.25
	0.5					0.245
	0.4
						0.24
	0.3
	0.2					0.235
	.1					0.23
	0	0.25	0.5	0.75	1	0	0.25	0.5	0.75	1	1.25	1.5	1.75	2
			1				Rusanov
							LP explicit
			0.9				LP implicit
						Rusanov(reference)
			0.8
			0.7
			0.6
			0.5
			0.4
			0.3
			0.2
			0	0.25	0.5		0.75		1
		Figure 1: Comparison of several schemes with a reference solution (density profile)
					Test 1 Test 2 Test 3
				Rusanov	4231	550	2630
				LP explicit	4297	551	2631
				LP implicit	63	41	151

* christophe.chalons@cea.fr † frederic.coquel@cmap.polytechnique.fr ‡ samuel.kokh@cea.fr § spillane@ann.jussieu.fr

The updated values of u k , τ k and ρ k are recovered from the formulas

k and ρ k = 1/τ k . The computation of (w ± k ) n+1= j is cheap and amounts to solving a tridiagonal system of linear equations, while the time-implicit definition of (e k ) n+1= j explicitly follows. Then, the transport equations involved in [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF] are associated with the following classical time-explicit update formula

and of course (α k ) n+1- j = (α k ) n+1= j .

Discretization of (5). Our objective is to propose a consistent approximation of ( 5) such that the overall algorithm is conservative for each partial mass, for the total momentum and for the total energy, as already motivated. First of all and similarly to [START_REF] Munkejord | Comparison of Roe-type methods for solving the two-fluid model with and without pressure relaxation[END_REF], the transport equation associated with α k is treated as follows:

for the partial mass, so that only the treatments of the momentum and total energy of each phase are now left. We propose

In order to get the expected overall conservativity properties, we pay a particular attention to the definitions of (p k )j, (pI )j, (p k u k )j and (pI uI )j . For any consistent definition of the flux (α k ) n j+1/2 , we set with κj ∈ [0, 1]

(pI uI )j = (pI )j(uI )j.

We choosed in practice (α k ) n j+1/2 = (α k ) n j or equivalently κj = 1.

With such definitions, it can be proved that under a suitable CFL condition based on the velocities u k and uI only, and not on the acoustic waves u k ± c k , the void fractions (α k ) n+1 j belong to (0, 1) if (α k ) n j do. We can also prove that under the same restriction on the time step (ρ k ) n+1 j is positive, as well as (ε k ) n+1- j and (p k ) n+1- j . Unfortunately, the positivity of (ε k ) n+1 j and (p k ) n+1 j is not proved at the moment.

Numerical experiments

For the sake of illustration, we present in this section the results given by our algorithm on three Riemann problems. They are all taken from [START_REF] Ambroso | A Godunov-type method for the seven-equation model of compressible two-phase flow[END_REF] and are fully described therein. Space and time orders of accuracy are one. The first one (top left) corresponds to an isolated contact discontinuity propagating with a positive velocity, while


Quentin Boehler

Marc Vedrines

Salih Abdelaziz

Philippe Poignet

Pierre Renaud

Influence of spring characteristics on the behavior of tensegrity mechanisms

Keywords: tensegrity mechanisms, stiffness computation, workspace estimation

There is today a growing interest for tensegrity mechanisms. Their analysis is however challenging because of their self-stress state. The most popular tensegrity mechanisms use linear springs as tensioned elements. Their synthesis for given user requirements is an open issue. In this article, we propose as a first step to better understand the influence of the spring characteristics, that constitute important design parameters. The influence of spring free length is in particular assessed, considering two planar tensegrity mechanisms. Impact of the spring selection on the workspace, the stiffness and the actuation requirements is observed. The simulation results outline that using nonzero free length springs can be of interest, and conclusions are given on further steps towards a synthesis method.

Introduction and scope of the study

Tensegrity systems can be defined [START_REF] Motro | Tensegrity: Structural Systems for the Future[END_REF] as systems in stable self-equilibrated state comprising a discontinuous set of compressed components, i.e. struts, inside a continuum of tensioned components, either cables or springs. Thanks to the self-stress state, prestress can be imposed for a same topology in order to modify the level of internal forces. A high stiffness-to-mass ratio can thus be reached. Moreover, since the components are axially loaded, they can be very light and so does the system [START_REF] Skelton | Tensegrity Systems[END_REF]. Tensegrity mechanisms use actuators to modify their configuration [START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-DOF tensegrity mechanism[END_REF][START_REF] Shekarforoush | Kinematic and static analyses of statically balanced spatial tensegrity mechanism with active compliant components[END_REF]. The tensioned components are then generally linear springs, in order to systematically reach Quentin Boehler (quentin.boehler@icube.unistra.fr), Marc Vedrines, Pierre Renaud ICube, France Salih Abdelaziz, Philippe Poignet LIRMM, France 1 equilibrated configurations without tension loss. They have recently received a lot of attention [START_REF] Skelton | Controllable tensegrity: a new class of smart structures[END_REF][START_REF] Sultan | Tensegrity flight simulator[END_REF][START_REF] Tran | The reverse displacement analysis of a tensegrity based parallel mechanism[END_REF], and could be interesting solutions for instance in surgical contexts, where lightweight devices with large workspaces are required [START_REF] Abdelaziz | Design considerations for a novel MRI compatible manipulator for prostate cryoablation[END_REF].

The analysis of tensegrity mechanisms is delicate compared to rigid-body mechanisms. The determination of a mechanism configuration requires for instance to take into account joint variables and external loads while respecting self-equilibrium state conditions. The corresponding static model is obtained using so-called form finding methods [START_REF] Tibert | Review of form-finding methods for tensegrity structures[END_REF]. Up to now, related works have been focused on the analysis of tensegrity mechanisms, but not yet on their synthesis. Even the relationship between the design parameters and the mechanism behavior, and the estimation of actuator performances for given specifications remains to be investigated. In this paper, we wish therefore to analyze the sensitivity of the tensegrity characteristics to their design parameters.

The spring characteristics strongly influence the mechanism behavior. We therefore concentrate our efforts on the parameters related to the linear springs used as tensioned components. The influence of spring stiffness is here considered, and more importantly the free length of the springs. To our knowledge, estimating the influence of nonzero free lengths on mechanism performances has not yet been covered in the literature.

As a result, we introduce in this paper an analysis of spring stiffness and free length on the behavior of tensegrity mechanisms. Two planar tensegrity mechanisms of Snelson-cross type [START_REF] Arsenault | Kinematic, static, and dynamic analysis of a planar one-degreeof-freedom tensegrity mechanism[END_REF][START_REF] Arsenault | Kinematic, static and dynamic analysis of a planar 2-DOF tensegrity mechanism[END_REF] are considered, that respectively exhibit 1 and 2-DOF. In section 2, the 1-DOF mechanism is considered. The mechanism is simple enough to allow a detailed analysis of free length influence on mechanism workspace, stiffness and criteria related to actuator requirements, namely actuator stroke, force and energy consumption. In section 3, a 2-DOF mechanism is analyzed. This mechanism can be used as a tool holder in a surgical task [START_REF] Abdelaziz | Design considerations for a novel MRI compatible manipulator for prostate cryoablation[END_REF]. Conclusions are finally given with emphasis on further steps towards a generalization for other tensegrity mechanisms.

2 Sensitivity analysis of a 1-DOF tensegrity mechanism

Mechanism description

The 1-DOF tensegrity mechanism, based on the Snelson-cross planar tensegrity, is represented in Figure 1. It consists in four nodes {A, B,C, D}. Two rigid struts of length L connect respectively A with C and B with D. Three linear springs are installed on the outside of the mechanism, with same stiffness k and free length l 0 . An actuator is used to modify the distance AB = ρ, with the node A anchored to the base and B constrained to move in a horizontal direction. In such a situation, the mechanism exhibits 1 DOF, with nodes C and D that have equal vertical displacements. The output variable is defined as y, the Y -coordinate of C and D, and ρ is the input variable. For sake of simplicity, a so-called neutral configuration is defined as

||AD|| = ||BC|| = ||CD|| = ρ =100

Modeling and analysis criteria

The mechanism modeling consists in deriving expressions of the direct and inverse static models. The output variable y is obtained for a given value of ρ by solving the direct static problem. An energy method, as described in [START_REF] Tibert | Review of form-finding methods for tensegrity structures[END_REF], is here considered. The mechanism configuration is then obtained by determining the local minimum of the potential energy U stored in the mechanism springs. In our context, U can be easily computed as

U = 1 2 k(2(l 1 -l 0 ) 2 + (l 2 -l 0 ) 2 ) (1)
with l 1 = ||AD|| = ||BC|| and l 2 = ||CD||. The lengths l 1 and l 2 can be expressed using simple geometrical relationships as functions of the variables ρ and y and the length L. In the following, we are only interested in positive solutions of y. The analytic expression of y = f (ρ) is obtained by solving ∂U ∂ y = 0 for a given ρ, using a computer algebra system (Mupad, The Mathworks Inc.). Only one admissible solution is found with the simulated conditions where ρ ∈ [70, 130] mm and l 0 ∈ [20, 80] mm. The inverse static model that expresses ρ as a function of y is difficult to obtain if nonzero free length of the springs are considered. As noticed previously in [START_REF] Crane | Kinematic analysis of a planar tensegrity mechanism with pre-stressed springs[END_REF], the use of nonzero free length springs in a tensegrity mechanism significantly increases its analysis complexity. In our situation, no analytical expression could be found. A numerical approach is then adopted, with a Levenberg-Marquardt algorithm to compute the ρ value for a given y as the solution of (yf (ρ)) 2 = 0. The tensegrity mechanism behavior is assessed by using two types of criteria. The first type corresponds to usual mechanism properties, namely the workspace and the stiffness. The workspace is the range of y values for a given range of ρ. The stiffness of the mechanism is defined as K y = δ F y δ u y with δ u y the infinitesimal displacement of nodes C and D along the Y -axis when an infinitesimal external vertical force δ F y is evenly applied on these nodes, while the actuator is locked. As stated in [START_REF] Arsenault | Kinematic, static, and dynamic analysis of a planar one-degreeof-freedom tensegrity mechanism[END_REF], the stiffness K y can be estimated with ∂ 2 U ∂ y 2 . The other type of criteria characterizes the actuation requirements of the mechanism. For a predefined path followed by the mechanism in quasistatic conditions, we compute • ∆ ρ: the corresponding required actuator stroke • F and F max : respectively the mean and maximum forces delivered by the actuator • ∆U: the variation of potential energy.

Meanwhile, the mean value K y and the variation ∆ K y of the stiffness are also estimated. In the presented results, the path is defined as symmetric with respect to the neutral configuration, with a 20-mm displacement. It is discretized in 11 steps to compute the analysis criteria.

Results and discussion

As it can be observed in Eq. ( 1), the configuration corresponding to the local minimum of the potential energy U does not depend on the springs stiffness k. The equilibrium configuration is therefore invariant with respect to k. Consequently, modifying the spring stiffness does not modify the mechanism workspace. The mechanism stiffness K y is proportional to k as it can be observed in Eq. ( 1). The energy variation ∆U and the actuator force F are also proportional to k. The energy variation corresponds to the energy to be delivered by the actuator. The actuation requirements are therefore linearly dependent on the spring stiffness for a given workspace.

The spring free length l 0 has on the contrary a non-linear influence on the analysis criteria. Figure 2 represents the computation of y and K y according to l 0 when ρ varies between 70 and 130 mm in 25 steps, for k = 4 N/mm. In this example, one may observe that the neutral configuration is not dependent on the spring free length. For a given y, increasing l 0 leads to an increase of K y while the force delivered by the actuator decreases. This may appear unintuitive, but it is explained by the reconfiguration of the mechanism when l 0 is increased (bottom of Figure 2). The workspace is on the contrary reduced for the same actuation stroke. For given mechanism specifications, interesting compromises therefore exist for nonzero free lengths.

Figure 3 shows the evolution of the criteria related to the actuator requirements. Increasing l 0 leads to the increase of the actuation stroke ∆ ρ while F, F max are reduced, as well as the potential energy variation ∆U. An optimal value of l 0 could be searched for in order to adapt the mechanism and the actuator force/displacement characteristics. The influence of the free length on the stiffness is complex, with an increase of the mean value K y but a reduction of the variation ∆ K y . As a conclusion of this case, the analysis suggests first that nonzero free lengths can be of interest for the mechanism performance, and that the choice of l 0 is a compromise between mechanism properties K y , ∆ K y and the actuator requirements ∆ ρ, F, F max and ∆U. 3 Analysis of a 2-DOF tensegrity mechanism

Mechanism description

A second actuator is integrated between nodes B and C in the previous tensegrity mechanism (Figure 1). The mechanism becomes a 2-DOF manipulator with two joint variables (ρ 1 , ρ 2 ) and its end-effector in D.

Modeling and analysis criteria

For this 2-DOF mechanism, no analytical direct static model could be obtained if springs have a nonzero free length. Hence, the coordinates of D are determined by minimizing the potential energy U using a Nelder-Mead algorithm. The inverse static model is not implemented since it is not necessary for the evaluation. Contrary to a rigid-body mechanism, the mechanism configuration does not depend only on (ρ 1 , ρ 2 ), but also on one parameter that can be either the angle ABD, the variable x, or y. Therefore U = U 1 (ρ 1 , ρ 2 , x) = U 2 (ρ 1 , ρ 2 , y), and the end-effector stiffnesses K x and K y in X-and Y -direction can be computed as follows

K x = ∂ 2 U 1 (ρ 1 , ρ 2 , x) ∂ x 2 (2) K y = ∂ 2 U 2 (ρ 1 , ρ 2 , y) ∂ y 2 (3)
For this mechanism, analysis criteria are related to the mechanism workspace, stiffness and the actuation requirements:

• A is the area of the reachable workspace • K x and K y are the mean values of respectively K x and K y over the workspace • ∆ K x and ∆ K y designate the variations of respectively K x and K y over the workspace • F and F max are computed as the mean and maximum forces in each actuator.

The workspace is explored with joint ranges that avoid mechanism singularity and loss of tensegrity configuration. A singularity occurs when ρ 1 = ρ 2 = L/2, i.e. BD is vertical. The lower bounds of ρ 1 and ρ 2 are therefore chosen equal to 75 mm. The strut BD must remain in compression to keep the mechanism in a tensegrity configuration. As a result, the upper bound of the joint variables is such that ρ 2 1 + ρ 2 2 < 5L 2 /2. A 5-mm step for each joint variable is considered for the computation. The area A is computed by extracting the edge of the discrete set of positions obtained after the exploration of the joint space. With such joint ranges, and thanks to the mechanism symmetry, F and F max are the same for both actuators.

Results and discussion

Similarly to section 2, the influence of the spring stiffness k on the analysis criteria is simple. It does not affect the equilibrium configurations, and linearly impacts the required actuators forces and the stiffnesses.

The variation of K y is represented in Figure 4, for a spring stiffness k = 3 N/mm. For such a mechanism, an end-effector position can sometimes be obtained with two different configurations and hence two different stiffness values. The map is therefore plotted in the joint space. Figure 5 depicts the corresponding reachable workspace. The l 0 value modifies the shape of the workspace boundaries, and the workspace area. In addition, tension losses can occur that limit the workspace, as it can be observed between the middle and right plots of Figure 5. Further analysis shows that the first tension loss takes place when D is at the middle of AC, and

l 0 (mm) l 0 (mm) l 0 (mm) l 0 (mm) K x (N/mm) ∆ K x (N/mm) K y (N/mm) ∆ K y (N/mm) F(N) F max (N) A(N/mm 2 )
Fig. 6 Evolution of the analysis criteria according to l 0 .

that l 0 = l 0 loss = L/2 ≃ 70 mm. The l 0 value also affects the mean value and the variation of the stiffness over the workspace (left and middle plots of Figure 4). Figure 6 shows the evolution of the analysis criteria according to l 0 . This evolution of the criteria is notably affected by choosing l 0 > l 0 loss , with in particular a fast reduction of the workspace area A. If l 0 < l 0 loss , K x , K y , F and F max decrease when l 0 increases. Although ∆ K x is also decreasing monotonously with respect to l 0 , the evolution of ∆ K y is more complex (as it can be seen on Figure 4). The lowest value of ∆ K x is indeed obtained for a nonzero value of l 0 . The evolution of the workspace area A is also non uniform. Increasing l 0 does not necessarily decrease the size of the workspace.

The criteria sensitivity to the free length value is different from the one observed for the 1-DOF mechanism. Increasing l 0 can lead to tension losses, which is quite logical, and affects stiffness and workspace properties. However, the free length can be used to lower the required actuator forces and an adequate choice of l 0 can help increasing the workspace size.

Conclusions

In this paper, influence of the spring characteristics on tensegrity mechanisms was considered. Using analysis criteria related to the mechanism properties and the actuator requirements, two planar mechanisms were analyzed. The first major remark is that increasing the free length of the springs can have a beneficial impact on the mechanism or the selection of an adequate mechanism. Whereas the spring stiffness influence can be easily assessed, relationships between the selected criteria and the free length is much more complex. The determination of the spring free length during a mechanism synthesis will need to be performed simultaneously with other design parameters related to the mechanism geometry. The observed significant nonlinearities of the sensitivity curves outline also that the selection of the optimization strategy will be delicate, and will be investigated as the next step of this work as well as the generalization for other tensegrity mechanisms.

2 = 2 Fig. 1

 221 Fig. 1 Kinematic schemes of the 1-DOF (left) and 2-DOF (right) tensegrity mechanisms. In black plain lines, rigid struts and in red dotted lines, linear springs.

Fig. 2 Fig. 3

 23 Fig.2Relationships between the spring free length l 0 and the stiffness K y and the workspace. Cases (a), (b) and (c) depict the mechanism at y=95 mm for l 0 =20, 50 and 80 mm respectively.

Fig. 4 Fig. 5

 45 Fig.4Variation of the stiffness K y in N/mm in the joint space for l 0 =20, 50 and 80 mm. In red circles: unreachable configurations because of tension loss.

Acknowledgment

This work was supported by French state funds managed by the ANR within the Investissements d'Avenir programme (Labex CAMI -ANR-11-LABX-0004) and by the Région Alsace.