Jean-Pierre Demailly

Hermitian-Yang-Mills approach to the conjecture of Griffiths on the positivity of ample vector bundles

Keywords: Ample vector bundle, Hermitian metric, Griffths positivity, Nakano positivity, Hermitian-Yang-Mills equation, Monge-Ampère equation, elliptic operator. MSC Classification 2020. 32J25, 53C07

Given a vector bundle of arbitrary rank with ample determinant line bundle on a projective manifold, we propose a new elliptic system of differential equations of Hermitian-Yang-Mills type for the curvature tensor. The system is designed so that solutions provide Hermitian metrics with positive curvature in the sense of Griffithsand even in the dual Nakano sense. As a consequence, if an existence result could be obtained for every ample vector bundle, the Griffiths conjecture on the equivalence between ampleness and positivity of vector bundles would be settled.

Introduction

Let X be a projective n-dimensional manifold. A conjecture due to Griffiths [START_REF] Griffiths | A: Hermitian differential geometry, Chern classes and positive vector bundles, Global Analysis[END_REF] stipulates that a holomorphic vector bundle E → X is ample in the sense of Hartshorne, meaning that the associated line bundle O P(E) (1) is ample, if and only if E possesses a Hermitian metric h such that the Chern curvature tensor Θ E,h = i∇ 2 E,h is Griffiths positive. In other words, if we let rank E = r and (1.1) Θ E,h = i 1≤j,k≤n, 1≤λ,µ≤r c jkλµ dz j ∧ dz k ⊗ e * λ ⊗ e µ in terms of holomorphic coordinates (z 1 , . . . , z n ) on X and of an orthonormal frame (e λ ) 1≤λ≤r of E, the associated quadratic form

(1.2) Θ E,h (ξ ⊗ v) := Θ E,h (ξ, ξ) • v, v h = 1≤j,k≤n, 1≤λ,µ≤r c jkλµ ξ j ξ k v λ v µ
should take positive values on non zero tensors ξ ⊗v ∈ T X ⊗E. A stronger concept is Nakano positivity (cf. [START_REF] Nakano | On complex analytic vector bundles[END_REF]), asserting that (1.3) Θ E,h (τ ) := 1≤j,k≤n, 1≤λ,µ≤r c jkλµ τ jλ τ kµ > 0 for all non zero tensors τ = j,λ τ jλ ∂ ∂z j ⊗ e λ ∈ T X ⊗ E. It is in fact interesting to consider also the curvature tensor of the dual bundle E * , which happens to be given by the opposite of the transpose of Θ E,h , i.e.

(1.4) Θ E * ,h * = -T Θ E,h = -1≤j,k≤n, 1≤λ,µ≤r c jkµλ dz j ∧ dz k ⊗ (e * λ ) * ⊗ e * µ .

This leads to the concept of dual Nakano positivity, stipulating that (1.5) -Θ E * ,h * (τ ) = 1≤j,k≤n, 1≤λ,µ≤r c jkµλ τ jλ τ kµ > 0 for all non zero tensors τ = j,λ τ jλ ∂ ∂z j ⊗ e * λ ∈ T X ⊗ E * . On the other hand, Griffiths positivity of Θ E,h is equivalent to Griffiths negativity of Θ E * ,h * , and implies the positivity of the induced metric on the tautological line bundle O P(E)(1) . By the Kodaira embedding theorem [START_REF] Kodaira | On Kähler varieties of restricted type[END_REF], the positivity of O P(E)(1) is equivalent to its ampleness, hence we see immediately from the definitions that (1.6) Θ E,h dual Nakano positive ⇒ Θ E,h Griffiths positive ⇒ E ample.

In this short note, we consider the following converse problem:

1.7. Basic question. Does it hold that E ample ⇒ Θ E,h dual Nakano positive ?

A positive answer would clearly settle the Griffiths conjecture, in an even stronger form. One should observe that Nakano positivity implies Griffiths positivity, but is in general a more restrictive condition. As a consequence, one cannot expect ampleness to imply Nakano positivity. For instance, T P n is easily shown to be ample (and Nakano semi-positive for the Fubini-Study metric), but it is not Nakano positive, as the Nakano vanishing theorem [START_REF] Nakano | On complex analytic vector bundles[END_REF] would then yield (1.8) H n-1,n-1 (P n , C) = H n-1 (P n , Ω n-1

P n ) = H n-1 (P n , K P n ⊗ T P n ) = 0.
On the other hand, it does not seem that there are any examples of ample vector bundles that are not dual Nakano positive, thus the above basic question is still legitimate, even though it might look very optimistic. We should mention here that subtle relations between ampleness, Griffiths and Nakano positivity are known to hold -for instance, B. Berndtsson [START_REF] Berndtsson | Curvature of vector bundles associated to holomorphic fibrations[END_REF] has proved that the ampleness of E implies the Nakano positivity of S m E ⊗ det E for every m ∈ N. See also [START_REF] Demailly | Relations entre les notions de positivité de P[END_REF] for an earlier direct and elementary proof of the much weaker result that the Griffiths positivity of E implies the Nakano positivity of E ⊗ det E, and [START_REF] Mourougane | Hodge metrics and positivity of direct images[END_REF] for further results analogue to those of [START_REF] Berndtsson | Curvature of vector bundles associated to holomorphic fibrations[END_REF].

So far, the Griffiths conjecture is known to hold when n = dim X = 1 or r = rank E = 1 (in which cases, Nakano and dual Nakano positivity coincide with Griffiths positivity). Proofs can be found in [Ume73, Theorem 2.6] and [START_REF] Campana | A characterization of ample vector bundles on a curve[END_REF]. In both cases, the proof is based on the existence of Harder-Narasimhan filtrations and on the Narasimhan-Seshadri theorem [START_REF] Narasimhan | Stable and unitary vector bundles on a compact Riemann surface[END_REF] for stable vector bundles -the 1-dimensional case of the Donaldson-Uhlenbeck-Yau theorem [START_REF] Donaldson | Anti-self-dual Yang-Mills connections over complex algebraic surfaces and stable vector bundles[END_REF], [START_REF] Uhlenbeck | On the existence of Hermitian-Yang-Mills connections in stable vector bundles[END_REF]. It is tempting to investigate whether techniques of gauge theory could be used to approach the Griffiths conjecture. In this direction, P. Naumann [Nau17] proposed a Kähler-Ricci flow method that starts with a given Finsler metric of positive curvature, and converges to a Hermitian metric. It is however unclear whether the flow introduced in [Nau17] preserves positivity, so it might very well produce in the limit a Hermitian metric that does not have positive curvature. Another related suggestion is V. Pingali's proposal made in [START_REF] Pingali | A vector bundle version of the Monge-Ampère equation[END_REF] to study the vector bundle Monge-Ampère equation (Θ E,h ) n = η Id E , where η is a positive volume form on X. Solving such an equation requires polystability in dimension n = 1, and, in general, a positivity property of (E, h) that is even stronger than Nakano positivity (and thus much stronger than ampleness).

In section 2, we describe a more flexible differential system based on a combination of a huge determinantal equation and a trace free Hermite-Einstein condition. It relies on the well known continuity method, and is designed to enforce positivity of the curvature, actually in the dual Nakano sense -a condition that could eventually still be equivalent to ampleness. We show that it is possible to design a non linear differential system that is elliptic and invertible, at least near the origin of time. It would however remain to check whether one can obtain long time existence of the solution for the said equation or one of its variants. Section 3 is devoted to the discussion of a related extremal problem, and a concept of volume for vector bundles.

Approach via a combination of Monge-Ampère and

Hermitian-Yang-Mills equations

Let E → X be a holomorphic vector bundle equipped with a smooth Hermitian metric h. If the Chern curvature tensor Θ E,h is dual Nakano positive, then the 1 r -power of the (n × r)dimensional determinant of the corresponding Hermitian quadratic form on T X ⊗ E * can be seen as a positive (n, n)-form

(2.1) det T X ⊗E * ( T Θ E,h ) 1/r := det(c jkµλ )
1/r (j,λ),(k,µ) idz 1 ∧ dz 1 ∧ . . . ∧ idz n ∧ dz n . Moreover, this (n, n)-form does not depend on the choice of coordinates (z j ) on X, nor on the choice of the orthonormal frame (e λ ) on E (but the orthonormality of (e λ ) is required). Conversely, given a Kähler metric ω 0 on X, the basic idea is that assigning a "matrix Monge-Ampère equation"

(2.2) det T X ⊗E * ( T Θ E,h ) 1/r = f ω n 0 ,
where f is a smooth positive function, may enforce the dual Nakano positivity of Θ E,h if that assignment is combined with a continuity technique from an initial starting point where positivity is known. For r = 1, we have T Θ E,h = Θ E,h = -i∂∂ log h, and equation (2.2) is a standard Monge-Ampère equation. If f is given and independent of h, Yau's theorem [START_REF] Yau | On the Ricci curvature of a complex Kähler manifold and the complex Monge-Ampère equation I[END_REF] guarantees the existence of a unique solution θ = Θ E,h > 0, provided E is an ample line bundle and X f ω n 0 = c 1 (E) n . One then gets a smoothly varying solution θ t = Θ E,h t > 0 when the right hand side f t of (2.2) varies smoothly with respect to some parameter t. Now, assuming E to be ample of rank r > 1, equation (2.2) becomes underdetermined, since the real rank of the space of hermitian matrices h on E is equal to r 2 , while (2.2) provides only one scalar equation. If E = 1≤j≤r E j splits as a direct sum of ample line bundles and we take a diagonal Hermitian structure h = h j on E, the nr ×nr determinant splits as a product of blocks, and equation (2.2) reduces to

(2.2 s ) 1≤j≤r Θ n E j ,h j 1/r = f ω n 0 .
This "split equation" can be solved for any f = f 1/r j with X f j ω n 0 = c 1 (E j ) n , just by solving the individual equations Θ n E j ,h j = f j ω n 0 , f j > 0, but the decomposition need not be unique. In this case, the Hölder inequality requires X f ω n 0 ≤ ( c 1 (E j ) n ) 1/r , and the equality can be reached by taking all f j 's to be proportional to f . In general, solutions might still exist, but the lack of uniqueness prevents us from getting a priori bounds. In order to recover a well determined system of equations, one needs to introduce (r 2 -1) additional scalar equations, or rather a matrix equation of real rank (r 2 -1). If E is ample, the determinant line bundle det E is also ample. By the Kodaira embedding theorem, we can find a smooth Hermitian metric η 0 on det E so that ω 0 := Θ det E,η 0 > 0 is a Kähler metric on X. In case E is ω 0 -stable or ω 0 -polystable, we know by the Donaldson-Uhlenbeck-Yau theorem that there exists a Hermitian metric h on E satisfying the Hermite-Einstein condition

(2.3) ω n-1 0 ∧ Θ E,h = 1 r ω n 0 ⊗ Id E ,
since the slope of E with respect to ω 0 ∈ c 1 (E) is equal to 1 r . In general, one cannot expect E to be ω 0 -polystable, but Uhlenbeck-Yau have shown that there always exist smooth solutions to a certain "cushioned" Hermite-Einstein equation. To make things more precise, let Herm(E) be the space of Hermitian (non necessarily positive) forms on E, and given a Hermitian metric h > 0, let Herm h (E, E) be the space of h-Hermitian endomorphisms u ∈ Hom(E, E); we denote by

(2.4) Herm(E) → Herm h (E, E), q → q such that q(v, w) = v, w q = q (v), w h
the natural isomorphism between Hermitian quadratic forms and Hermitian endomorphisms, which depends of course on h. We also let

(2.5) Herm • h (E, E) = u ∈ Herm h (E, E) ; tr u = 0
be the subspace of "trace free" Hermitian endomorphisms. In the sequel, we fix a reference Hermitian metric H 0 on E such that det H 0 = η 0 , so that Θ det E,det H 0 = ω 0 > 0. By [UhY86, Theorem 3.1], for every ε > 0, there exists a smooth Hermitian metric q ε on E such that

(2.6) ω n-1 0 ∧ Θ E,q ε = ω n 0 ⊗ 1 r Id E -ε log q ε ,
where q ε is computed with respect to H 0 , and log u denotes the logarithm of a positive Hermitian endomorphism u. The intuitive reason is that the term log q ε introduces sufficient "friction" to avoid any explosion of approximating solutions when using a standard continuity method (see sections 2,3 in [START_REF] Uhlenbeck | On the existence of Hermitian-Yang-Mills connections in stable vector bundles[END_REF]). On the other hand, when ε → 0, the metrics q ε become "more and more distorted" and yield asymptotically a splitting of E in weakly holomorphic subbundles corresponding to the Harder-Narasimhan filtration of E with respect to ω 0 . If we write det q ε = e -ϕ det H 0 and take the trace in (2.6), we find

ω n-1 0 ∧(ω 0 +i∂∂ϕ) = ω n 0 (1+εϕ), hence ω n-1 0 ∧ i∂∂ϕ -εϕ ω n 0 = 0.
A standard application of the maximum principle shows that ϕ = 0, thus (2.6) implies det q ε = det H 0 and log q ε ∈ Herm • H 0 (E, E). In general, for an arbitrary Hermitian metric h, we let

(2.7) Θ • E,h = Θ E,h - 1 r Θ det E,det h ⊗ Id E ∈ C ∞ (X, Λ 1,1 R T * X ⊗ Herm • h (E, E))
be the curvature tensor of E ⊗ (det E) -1/r with respect to the trivial determinant metric h • := h ⊗ (det h) -1/r . Equation (2.6) is equivalent to prescribing det q ε = det H 0 and

(2.8)

ω n-1 0 ∧ Θ • E,q ε = -ε ω n 0 ⊗ log q ε .
This is a matrix equation of rank (r 2 -1) that involves only q • ε and does not depend on det q ε . Notice that we have here log q ε ∈ Herm • H 0 (E, E), but also log q ε ∈ Herm • q ε (E, E). In this context, given α > 0 large enough, it seems natural to search for a time dependent family of metrics h t (z) on the fibers E z of E, t ∈ [0, 1], satisfying a generalized Monge-Ampère equation

(2.9) det T X ⊗E * T Θ E,h t + (1 -t)α ω 0 ⊗ Id E * 1/r = f t ω n 0 , f t > 0,
and trace free Hermite-Einstein conditions

(2.9 • ) ω n-1 t ∧ Θ • E,h t = g t ,
with smoothly varying families of functions

f t ∈ C ∞ (X, R), Hermitian metrics ω t > 0 on X and sections g t ∈ C ∞ (X, Λ n,n R T * X ⊗ Herm • h t (E, E)), t ∈ [0, 1].
Here, we start e.g. with the Yau-Uhlenbeck solution h 0 = q ε of (2.6) (so that det h 0 = det H 0 ), and take α > 0 so large that T Θ E,h 0 + α ω 0 ⊗ Id E * > 0 in the sense of Nakano. If these conditions can be met for all t ∈ [0, 1] without any explosion of the solutions h t , we infer from (2.9) that (2.9

+ ) T Θ E,h t + (1 -t)α ω 0 ⊗ Id E * > 0 in the sense of Nakano for all t ∈ [0, 1]
. At time t = 1, we will then get a Hermitian metric h 1 on E such that Θ E,h 1 is dual Nakano positive. We still have the freedom of adjusting f t , ω t and g t in equations (2.9) and (2.9 • ). We have a system of differential equations of order 2, and any choice of the right hand sides of the form

f t (z) = F (t, z, h t (z), D z h t (z), D 2 z h t (z)) > 0, (2.10) g t (z) = G(t, z, h t (z), D z h t (z), D 2 z h t (z)) ∈ C ∞ (X, Λ 1,1 R T * X ⊗ Herm • (E, E)) (2.10 • )
is a priori acceptable for the sake of enforcing the positivity condition (2.9 + ), although the presence of second order terms D 2 z h t (z)) might affect the principal symbol of the equations. In equation (2.9 • ), the metrics ω t could possibly be taken to depend on t, but unless some commodity reason would appear in next stages of the analysis, it seems simpler to set ω t = ω 0 independent of t. At this stage, we have the following 2.11. Theorem. Let (E, H 0 ) be a smooth Hermitian holomorphic vector bundle such that E is ample and ω t = ω 0 = Θ det E,det H 0 > 0. Then the system of equations (2.9, 2.9 • ) is a well determined (essentially non linear) elliptic system of equations for all choices of smooth right hand sides

f t = F (t, z, h t , D z h t ) > 0, g t = G(t, z, h t , D z h t , D 2 z h t ) ∈ Herm • (E, E),
provided that the symbol η h of the linearized operator u → DG D 2 h (t, z, h, Dh, D 2 h) • D 2 u has an Hilbert-Schmidt norm sup ξ∈T * X ,|ξ| ω 0 =1 η h (ξ) h ≤ (r 2 + 1) -1/2 n -1 for any of the metrics h = h t involved. If a smooth solution h t exists on the whole time interval [0, 1], then E is dual Nakano positive.

Proof. If we write a hermitian metric h on E under the form h(v, w) = h(v), w H 0 with h ∈ Herm h 0 (E, E), we have h = H 0 h in terms of matrices. The curvature tensor is given by the usual formula Θ

E,h = i∂(h -1 ∂h) = i∂( h -1 ∂ H 0 h), where ∂ H 0 s = H -1 0 ∂(H 0 s
) is the (1, 0)component of the Chern connection associated with H 0 on E. For simplicity of notation, we put M := Herm(E),

M h = Herm h (E, E), M • h = Herm • h (E, E).
The system of equations (2.9, 2.9 • ) is associated with the non linear differential operator

P : C ∞ (X, M ) → C ∞ (X, R ⊕ M • h ), h → P (h)
defined by

P (h) = ω -n 0 det T X ⊗E * T Θ E,h +(1 -t)α ω 0 ⊗Id E * 1/r , ω n-1 0 ∧Θ E • ,h -G(t, z, h, Dh, D 2 h) .
It is by definition elliptic at h if its linearization u → (dP ) h (u) is an elliptic linear operator, a crucial fact being that M and R ⊕ M • h have the same rank r 2 over the field R. Our goal is to compute the symbol

σ dP ∈ C ∞ (X, S 2 T R X ⊗ Hom(M, R ⊕ M • h ))
of dP , and to check that u → σ dP (ξ) • u is invertible for every non zero vector ξ ∈ T * X . We pick an infinitesimal variation δh of h in C ∞ (X, M ), and represent it as

δh = u •, • h with u ∈ M h = Herm h (E, E).
In terms of matrices, we have δh = hu, i.e. u = (u λµ ) = h -1 δh is the "logarithmic variation of h". In this setting, we evaluate (dP ) h (u) in orthonormal coordinates (z j ) 1≤j≤n on X relatively to ω 0 . We have h

+ δh = h(Id + u) and (h + δh) -1 = (Id -u)h -1 modulo O(u 2 ), thus dΘ E,h (u) = i∂(h -1 ∂(hu)) -i∂(uh -1 ∂h) = i∂∂u + i∂(h -1 ∂h u) -i∂(uh -1 ∂h) = i∂∂ h * ⊗h u = -i∂ h * ⊗h ∂u, (2.12)
where ∂ h * ⊗h denotes here the (1, 0)-component of the Chern connection on the holomorphic vector bundle Hom(E, E) = E * ⊗ E induced by the metric h * ⊗ h. As a consequence, the order 2 term of the linearized operator is just dΘ E,h (u) [2] = -i∂∂u, and the logarithmic differential of the first scalar component P R (h) of P (h) has order 2 terms given by (2.13)

P R (h) -1 dP R,h (u) [2] = 1 r tr(-θ -1 • T i∂∂u) = - 1 r (det θ) -1 j,k,λ,µ θ jkλµ ∂ 2 u λµ ∂z j ∂z k ,
where θ is the (n × r)-matrix of θ = θ(t, h) = T Θ E,h + (1 -t)αω 0 ⊗ Id E * > 0, θ its co-adjoint and θ -1 = (det θ) -1 T θ, so that P R (h) = ω -n 0 (det θ) 1/r . We also have to compute the order 2 terms in the differential of the second component

h → P • (h) = ω -n 0 ω n-1 0 ∧ Θ • E,h -G(t, z, h, Dh, D 2 h) . Let us set u = 1 r tr u ⊗ Id E + u • , u • ∈ M • , and tr u = λ u λλ ∈ R. Putting τ = 1 r tr u, this actually gives an isomorphism M h → R ⊕ M • h , u → (τ, u • ). Since u • is the logarithmic variation of h • = h(det h) -1/r , we get (2.14) (dP • ) h (u) [2] = ω -n 0 -ω n-1 0 ∧ i∂∂u • -DG D 2 h • D 2 u .
If we fix a Hermitian metric h and take a non zero cotangent vector 0 = ξ ∈ T * X , the symbol σ dP is given by an expression of the form (2.15)

σ (dP ) h (ξ) • u = - (det θ) -1+1/r r ω n 0 j,k,λ,µ θjkλ mu ξ j ξ k u λµ , 1 n |ξ| 2 u • + σ G (ξ) • u
where σ G is the principal symbol of the operator DG D 2 h • D 2 . If g t = G(t, z, h t , Dh t ) is independent of D 2 h t , the latter symbol σ G is equal to 0 and it is easy to see from (2.13) that u → σ (dP ) h (ξ) • u is an isomorphism in Hom(M h , R ⊕ M • h ). In fact, the first summation yields

j,k,λ,µ θjkλµ ξ j ξ k u λµ = j,k,λ,µ θjkλµ ξ j ξ k u • λµ + 1 r j,k,λ θjkλλ ξ j ξ k tr u.
By an easy calculation, we get an inverse operator

R ⊕ M • h → M h , (τ, v) → u where -r ω n 0 (det θ) 1-1/r τ = j,k,λ,µ θjkλµ ξ j ξ k u • λµ + 1 r j,k,λ θjkλλ ξ j ξ k tr u, -v = 1 n |ξ| 2 u • , hence u • = -n |ξ| 2 v and σ (dP ) h (ξ) -1 • (τ, v) = n |ξ| 2 j,k,λ,µ θjkλµ ξ j ξ k v λµ -r ω n 0 (det θ) 1-1/r τ j,k,λ θjkλλ ξ j ξ k Id E - n |ξ| 2 v.
Let us take the Hilbert-Schmidt norms |u| 2 = λ,µ |u λµ | 2 on M h = Herm h (E, E), and c|τ | 2 +|v| 2 on R ⊕M • h (h being the reference metric, and C > 0 a constant). By homogeneity, we can also assume |ξ| = |ξ| ω 0 = 1. Since ( j,k θjkλµ ξ j ξ k ) 1≤λ,µ≤r is a positive Hermitian matrix by the Nakano positivity property, its trace is a strict upper bound for the largest eigenvalue, and we get

j,k,λ θjkλµ ξ j ξ k v λµ 2 ≤ (1 -δ) j,k,λ θjkλλ ξ j ξ k 2 λ |v λµ | 2 .
The Cauchy-Schwarz inequality yields

j,k,λ,µ θjkλµ ξ j ξ k v λµ 2 ≤ r(1 -δ) j,k,λ θjkλλ ξ j ξ k 2 λ,µ |v λµ | 2 . For |ξ| = 1, as Id E ⊥ M • and |Id E | 2 = r, this implies σ (dP ) h (ξ) -1 • (τ, v) 2 ≤ nr 1/2 (1 -δ) 1/2 |v| + r ω n 0 (det θ) 1-1/r j,k,λ θjkλλ ξ j ξ k |τ | 2 r + n 2 |v| 2 < (n 2 r 2 + n 2 )(C|τ | 2 + |v| 2 )
for C large enough. By a standard pertubation argument, (2.13) remains bijective if | σ G (ξ)| h is less than the inverse of the norm of σ (dP ) h (ξ) -1 , i.e. (r 2 + 1) -1/2 n -1 . Similarly, one could also allow the scalar right hand side F to have a "small dependence" on D 2 h t , but this seems less useful.

Our next concern is to ensure that the existence of solutions holds on an open interval of time [0, t 0 [ (and hopefully on the whole interval [0, 1]). In the case of a rank one metric h = e -ϕ , it is well-known that the Kähler-Einstein equation (ω 0 +i∂∂ϕ t ) n = e tf +λϕ t ω n 0 more easily results in getting openness and closedness of solutions when applying the continuity method for λ > 0, as the linearized operator ψ → ∆ ω ϕ t ψ -λψ is always invertible. One way to generalize the Kähler-Einstein condition to the case of higher ranks r ≥ 1 is to take

(2.16) f t (z) = (det H 0 (z)/ det h t (z)) λ a 0 (z), λ ≥ 0,
where a 0 (z) = ω -n 0 det( T Θ E,h 0 + αω 0 ⊗ Id E * ) 1/r > 0 is chosen so that the equation is satisfied by h 0 at t = 0 (the choice λ > 0 has the interest that f t gets automatically rescaled by multiplying h t by a constant, thus ensuring strict invertibility). For the trace free part, what is needed is to introduce a friction term g t that helps again in getting invertibility of the linearized operator, and could possibly avoid the explosion of solutions when t grows to 1. A choice compatible with the Yau-Uhlenbeck solution (2.8) at t = 0 is to take

(2.16 • ) g t = -ε (det H 0 (z)/ det h t (z)) µ ω n 0 ⊗ log h • t , ε > 0, µ ∈ R,
if one remembers that det h 0 = det H 0 . These right hand sides do not depend on higher derivatives of h t , so Theorem 2.11 ensures the ellipticity of the differential system. Moreover:

2.17. Theorem. For ε ≥ ε 0 (h t ) and λ ≥ λ 0 (h t )(1 + µ 2 ) with ε 0 (h t ) and λ 0 (h t ) large enough, the elliptic differential system defined by (2.9, 2.9 • ) and (2.16, 2.16 • ), namely

ω -n 0 det T X ⊗E * T Θ E,h t + (1 -t)α ω 0 ⊗ Id E * 1/r = det H 0 (z) det h t (z) λ a 0 (z) ω -n 0 ω n-1 0 ∧ Θ • E,h t = -ε det H 0 (z) det h t (z) µ log h • t ,
possesses an invertible elliptic linearization. As a consequence, for such values of ε and λ, there exists an open interval [0, t 0 ) ⊂ [0, 1] on which the solution h t exists.

Proof. We replace the operator P :

C ∞ (X, M ) → C ∞ (X, R ⊕ M • h
) used in the proof of Theorem 2.9 by P = ( P R , P • ) defined by

P R (h) = ω -n 0 (det h(z)/ det H 0 (z)) λ det T X ⊗E * T Θ E,h + (1 -t)α ω 0 ⊗ Id E * 1/r , P • (h) = ω -n 0 ω n-1 0 ∧ Θ • E,h + ε (det h(z)/ det H 0 (z)) -µ log h • .
Here, we have to care about the linearized operator dP itself, and not just with its principal symbol. We let again u = h -1 δh ∈ Herm h (E, E) and use formula (2.12) for dΘ E,h (u). This implies

P R (h) -1 d P R,h (u) = λ tr u - 1 r tr T X ⊗E * θ -1 • T i∂ h * ⊗h ∂u .
We need the fact that h • = h • (det h) -1/r possesses, when viewed as a Hermitian endomorphism, a logarithmic variation

( h • ) -1 δ h • = u • = u - 1 r tr u • Id E .
By the classical formula expressing the differential of the logarithm of a matrix, we have

d log g(δg) = 1 0 (1 -t)Id + tg -1 δg (1 -t)Id + tg -1 dt, which implies d log h • (u) = 1 0 (1 -t)Id + t h • -1 h • u • (1 -t)Id + t h • -1 dt.
In the end, we obtain

(d P • ) h (u) = -ω -n 0 ω n-1 0 ∧ i∂ h * ⊗h ∂u • + ε det h(z) det H 0 (z) -µ 1 0 (1 -t)Id + t h • -1 h • u • (1 -t)Id + t h • -1 dt -µ tr u log h • .
In order to check the invertibility, we use the norm |τ | 2 + C|v| 2 on R ⊕ M • h and compute the L 2 inner product (d P ) h (u), (τ, u • ) over X, where τ = 1 r tr u. The ellipticity of operators -i∂ H ∂ implies that it has a discrete sequence of eigenvalues converging to +∞, and that we get Gårding type inequalities of the form

-i∂ H ∂v, v H ≥ c 1 ∇v 2 H -c 2 v 2 H where c 1 , c 2 > 0 depend on H. We apply such inequalities to v = τ , H = 1, and v = u • , H = h * ⊗ h, replacing u with u = τ Id + u • .
From this, we infer

(d P ) h (u), (τ, u • ) ≥ c 1 dτ 2 -c 2 τ 2 + λr τ 2 - 1 r tr T X ⊗E * θ -1 • T i∂ h * ⊗h ∂u • , τ + C c • 1 ∇u • 2 -c • 2 u • 2 + c 3 ε u • 2 -c 4 ε |µ| τ u •
where all constants c j may possibly depend on h. An integration by parts yields

1 r tr T X ⊗E * θ -1 • T i∂ h * ⊗h ∂u • , τ ≤ c 5 ∇u • ( dτ + τ ) ≤ 1 2 c 1 dτ 2 + τ 2 + c 6 ∇u • 2
and we have

c 4 ε |µ| τ u • ≤ 1 2 c 3 ε u • 2 + c 7 εµ 2 τ 2 . If we choose ε ≥ 2c • 2 /c 3 + 1, C ≥ c 6 /c • 1 + 1 and λr ≥ c 2 + 1 2 c 1 + Cc 7 εµ 2 + 1, we finally get (d P ) h (u), (τ, u • ) ≥ 1 2 c 1 dτ 2 + τ 2 + c • 1 ∇u • 2 + 1 2 Cc 3 ε u • 2
and conclude that (d P ) h is an invertible elliptic operator. The openness property at t = 0 then follows from standard results on elliptic PDE's.

2.18. Remarks. (a) Theorem 2.17 is not very satisfactory since the constants ε 0 (h t ) and λ 0 (h t ) depend on the solution h t . It would be important to know if one can get sufficiently uniform estimates to make these constants independent of h t , thereby guaranteeing the long time existence of solutions. This might require modifying somewhat the right hand side of our equations, especially the trace free part, while taking a similar determinantal Monge-Ampère equation that still enforces the dual Nakano positivity of the curvature tensor. The Yau iteration technique used in [START_REF] Yau | On the Ricci curvature of a complex Kähler manifold and the complex Monge-Ampère equation I[END_REF] to get 0 order estimates for Monge-Ampère equations will probably have to be adapted to this situation.

(b) The non explosion of solutions when t → 1 does not come for free, since this property cannot hold when det E is ample, but E is not. One possibility would be to show that an explosion at time t 0 < 1 produces a "destabilizing subsheaf" S contradicting the ampleness of E/S, similarly to what was done in [START_REF] Uhlenbeck | On the existence of Hermitian-Yang-Mills connections in stable vector bundles[END_REF] to contradict the stability hypothesis.

2.19. Variants. (a) The determinantal equation always yields a Kähler metric

β t := tr E Θ E,h t + (1 -t)α ω 0 ⊗ Id E = Θ det E,det h t + r(1 -t)α ω 0 > 0.
An interesting variant of the trace free equation is

( * ) ω -n 0 ω n-1 t ∧ Θ • E,h t = -ε det H 0 (z) det h t (z) µ log h • t
with ω t = 1 rα+1 β t (notice that β 0 = (rα + 1)ω 0 ). It is then important to know whether the corresponding differential system is still elliptic with an invertible linearization. According to equation ( * ), the Herm(E, E) • part of the differential system depends on the functional

P • (h) = ω -n 0 ω n-1 t ∧ Θ • E,h + ε (det h(z)/ det H 0 (z)) -µ log h • ,
and, with respect to the functional used in Theorem 2.17, the differential d P • h (u) acquires one additional term coming from the variation of ω n-1 t . With the same notation as in our previous calculations, we have Θ det E,det h t = -i∂∂ log det(h t ) and δ(β t ) h (u) = -i∂∂ tr u, hence

(d P • ) h (u) = -ω -n 0 ω n-1 t ∧ i∂ h * ⊗h ∂u • + n-1 rα+1 ω n-2 t ∧ i∂∂ tr u ∧ Θ • E,h + ε det h(z) det H 0 (z) -µ 1 0 (1 -t)Id + t h • -1 h • u • (1 -t)Id + t h • -1 dt -µ tr u log h • .
Putting again τ = tr u, this requires to estimate one extra term appearing in the L 2 inner product (d P ) h (u), (τ, u • ) , namely

(ω n 0 ) -1 ω n-2 t ∧ i∂∂τ ∧ Θ • E,h , u • .
We can apply the same integration by part argument as before to conclude that (d P ) h is again invertible, under a similar hypothesis λ ≥ λ 0 (h t )(1 + µ 2 ), at least for t small. A very recent note posted by V.P. Pingali [START_REF] Pingali | A note on Demailly's approach towards a conjecture of Griffiths[END_REF] shows that when E is ω 0 -stable and h 0 is taken to be the Hermite-Einstein metric, the trace free part of the differential system used in Theorem 2.17 has a solution of the form h t = h 0 e -ψ t , thus always "conformal" to h 0 . There are cases where the dual Nakano positivity of h 0 is doubtful. As a consequence, even in that favorable case, it is unclear whether a long time existence result can hold for the total system, unless stronger restrictions on the Chern classes are made. Equation ( * ) does not seem to entail such contraints, and may thus be better suited to the investigated problem.

(b) In a first step towards solving (2.6), [START_REF] Uhlenbeck | On the existence of Hermitian-Yang-Mills connections in stable vector bundles[END_REF] consider equations that have even stronger friction terms, taking the right hand side to be of the form

ω n-1 0 ∧ Θ E,h = ω n 0 ⊗ -ε log h + σ h -1/2 Γ 0 h 1/2 -Γ 0 ), σ > 0,
and letting σ → 0 at the end of the analysis. Here we can do just the same, for instance by adding a term equal to a multiple of ( h • t ) -1/2 Γ t ( h • t ) 1/2 -Γ t in the trace free equation, as such terms are precisely trace free for any Γ t ∈ C ∞ (X, Hom(E, E)).

A concept of Monge-Ampère volume for vector bundles

If E → X is an ample vector bundle of rank r, the associated line bundle

O P(E) (1) → Y = P(E)
is ample, and one can consider its volume c 1 (O P(E) (1)) n+r-1 . It is well known that this number (which is an integer) coincides with the Segre number X (-1) n s n (E), where (-1) n s n (E) is the n-th Segre class of E. Let us assume further that E is dual Nakano positive (if the solution of the Hermitian-Yangs-Mills differential system of §2 is unobstructed, this would follow from the ampleness of E). One can then introduce the following more involved concept of volume, which we will call the Monge-Ampère volume of E :

(3.1) MAVol(E) = sup h X det T X ⊗E * (2π) -1 T Θ E,h 1/r ,
where the supremum is taken over all smooth metrics h on E such that T Θ E,h is Nakano positive. This supremum is always finite, and in fact we have 3.2. Proposition. For any dual Nakano positive vector bundle E, one has

MAVol(E) ≤ r -n c 1 (E) n .
Proof. We take h to be a hermitian metric on E such that T Θ E,h is Nakano positive, and consider the Kähler metric

ω = (2π) -1 Θ det E,det h = (2π) -1 tr E * T Θ E,h ∈ c 1 (E).
If (λ j ) 1≤j≤nr are the eigenvalues of the associated hermitian form (2π) -1T Θ E,h with respect to ω ⊗ h, we have

det T X ⊗E * (2π) -1 T Θ E,h 1/r = j λ j 1/r
ω n and j λ j 1/nr ≤ 1 nr j λ j by the inequality between geometric and arithmetic means. Since

j λ j = tr ω tr E * (2π) -1 T Θ E,h = tr ω ω = n, we conclude X det T X ⊗E * (2π) -1 T Θ E,h 1/r ≤ X 1 nr j λ j n ω n = r -n X ω n = r -n c 1 (E) n .
The proposition follows.

Remarks. (a)

In case E = 1≤j≤r E j and h = 1≤j≤r h j are split, with all metrics h j normalized to have proportional volume forms ((2π) -1 Θ E j ,h j ) n = β j ω n with suitable constants β j > 0, we get β j = c 1 (E j ) n /c 1 (E) n , and the inequality reads

1≤j≤r c 1 (E j ) n 1/r ≤ r -n c 1 (E) n .
It is an equality when E 1 = • • • = E r , thus Proposition 3.2 is optimal as far as the constant r -n is concerned. For E = 1≤j≤r E j split with distinct ample factors, it seems natural to conjecture that

MAVol(E) = 1≤j≤r c 1 (E j ) n 1/r , i.e. that the supremum is reached for split metrics h = h j . In case E is a non split extension 0 → A → E → A → 0 with A an ample line bundle -this is possible if H 1 (X, O X ) = 0, e.g. on an abelian variety -we strongly suspect that MAVol(E) = c 1 (A) n but that the supremum is not reached by any smooth metric, as we have E semi-stable but not polystable.

(b) It would be interesting to characterize the "extremal metrics" h achieving the supremum in (3.1) when they exist. The calculations made in §2 show that they satisfy some Euler-Lagrange equation

X (det θ) 1/r • tr T X ⊗E * θ -1 • T i∂ h * ⊗h ∂u = 0 ∀u ∈ C ∞ (X, Herm(E)),
where θ is the (n × r)-matrix representing T Θ E,h . After performing two integration by parts that free u from any differentiation, we get a fourth order non linear differential system that h has to satisfy. Remark 3.3 (a) leads us to suspect that this system is not always solvable, but the addition of adequate lower order "friction terms" might make it universally solvable. This could possibly yield a better alternative to the more naive order 2 differential system we proposed in §2 to study the Griffiths conjecture.

(c) When r > 1, one may wonder what is the infimum inf h X det T X ⊗E * (2π) -1 T Θ E,h 1/r .

In the split case (E, h) = (E j , h j ), we can normalize Θ E j ,h j to satisfy Θ n E j ,h j = f j ω n with X f j ω n = c 1 (E j ) n , f j > 0. Then

X det T X ⊗E * (2π) -1 T Θ E,h 1/r = X (f 1 • • • f r ) 1/r ω n
and this integral becomes arbitray small if we take the f j 's to be large on disjoint open sets, and very small elsewhere. This example leads us to suspect that one always has inf h X det T X ⊗E * (2π) -1 T Θ E,h 1/r = 0 for r > 1. The "friction terms" used in our differential systems should be chosen so as to prevent any such shrinking of the volume.

* This work is supported by the European Research Council project "Algebraic and Kähler Geometry" (ERC-ALKAGE, grant No. 670846 from September 2015)


Jurado Realpe
email: juradoreal@lirmm.fr

Salih Abdelaziz
email: abdelaziz@lirmm.fr

Philippe Poignet
email: poignet]@lirmm.fr

Model Predictive Controller for a Planar Tensegrity Mechanism with decoupled position and stiffness control

Keywords: Dynamic modeling, MPC controller, Stiffness modulation, Tensegrity mechanism, Trajectory tracking

Precise trajectory tracking and stiffness modulation for tensegrity mechanisms are a challenging topic that can open new horizon of applications for this type of systems. This paper presents a new control strategy of tensegrity mechanisms using a model predictive controller (MPC). Based on a dynamic model, the proposed approach allows to track trajectories with low and relatively high dynamics as well as to modulate the mechanism stiffness by changing only the controller's parameters. Trajectories of 30s, 5s and 1s are performed showing a trajectory tracking improvement of up to 64% in the root mean square error when compared to literature results.

Introduction

The word tensegrity comes from the contraction of tensional and integrity [START_REF] Fuller | Tensile-integrity structures[END_REF]. A tensegrity system is composed by compressive and tensile elements that interacts to define a stable volume [START_REF] Pugh | An introduction to tensegrity[END_REF]. The integrity of the structure is then ensured by only the internal forces produced by the tensile elements (cables, springs). Their large workspace to size ratio, their light weight and their ability to modulate stiffness, make tensegrity systems an interesting solution for structural applications [START_REF] Snelson | Space Frame Structure Made By 3-D Weaving of Rod Members[END_REF] [4], deployable systems [START_REF] Skelton | Tensegrity Systems[END_REF] [START_REF] Sychterz | Deployment and Shape Change of a Tensegrity Structure Using Path-Planning and Feedback Control[END_REF], mobile robots [START_REF] Paul | Design and control of tensegrity robots for locomotion[END_REF] [START_REF] Sunspiral | Super Ball Bot -Structures for Planetary Landing and Exploration[END_REF] or manipulators [START_REF] Furet | Kinetostatic Analysis and Actuation Strategy of a Planar Tensegrity 2-X Manipulator[END_REF] [START_REF] Boehler | From modeling to control of a variable stiffness device based on a cable-driven tensegrity mechanism[END_REF].

Different control approaches have been proposed to control the motion of tensegrity mechanisms. The authors are commonly focused in the generation of gait for locomotion using open loop controllers [START_REF] Masic | Open-loop control of class-2 tensegrity towers[END_REF] [START_REF] Iscen | Learning Tensegrity Locomotion Using Open-Loop Control Signals and Coevolutionary Algorithms[END_REF], neuronal networks [START_REF] Kanchanasaratool | Modelling and control of class NSP tensegrity structures[END_REF] [START_REF] Hustig-Schultz | Morphological design for controlled tensegrity quadruped locomotion[END_REF], kinematic models [START_REF] Boehler | From modeling to control of a variable stiffness device based on a cable-driven tensegrity mechanism[END_REF] or Lyapunov function [START_REF] Wroldsen | Modelling and control of non-minimal non-linear realisations of tensegrity systems[END_REF]. Nevertheless, the trajectory tracking has not been a priority. To our knowledge, only two works [START_REF] Boehler | From modeling to control of a variable stiffness device based on a cable-driven tensegrity mechanism[END_REF] [START_REF] Sabelhaus | Model-Predictive Control with Reference Input Tracking for Tensegrity Spine Robots[END_REF] propose strategies to control finely trajectory tracking. On one hand, a PID controller and a tension distribution algorithm [START_REF] Boehler | From modeling to control of a variable stiffness device based on a cable-driven tensegrity mechanism[END_REF], inspired from cable-driven robots controllers, has been proposed to control the position and the stiffness of a 1 DOF planar tensegrity mechanism. This controller presents good position tracking for trajectories with low dynamic but loose performance when tracking relatively high dynamics trajectories. On the other hand, a Model Predictive Controller (MPC) with inverse static optimization has been introduced to realize the trajectory tracking of a tensegrity spine robot. Nevertheless, stiffness modulation is not addressed.

Our interest here is to fully exploit a tensegrity system by controlling its position and stiffness. A dynamic-based control strategy would be pertinent to deliver good trajectory tracking performance. A MPC using a dynamic model is here proposed to modulate the stiffness and trajectory tracking of the mechanism in an asynchronous way by only modifying the control parameters of the controller. The proposed strategy is experimentally validated on a 1 DOF planar cable-driven tensegrity mechanism.

This paper is organized as follow. The static, stiffness and dynamic models of the planar tensegrity mechanism is developed in section 2. Then a MPC strategy is proposed in section 3. Position tracking and stiffness modulation results are discussed in section 4. Finally conclusions and future work are presented in section 5.

System modelling

The system studied here is a 1 DOF planar tensegrity mechanism, whose end effector is considered as the bar bc (Fig 1). The end effector orientation is denoted θ . The structure is composed by an articulated equilateral parallelogram, of length L, and is driven by two cables. Each cable is attached to one corner of the parallelogram at one end and to a linear spring, of stiffness k, on the other end. The mechanism is actuated by two actuators. The angular position for each motor i is denoted α i . Pulleys, of radius R, are mounted on the actuators so that ρ i = Rα i describes the displacement of the cable i. This redundancy of actuation is performed so as to allow stiffness modulation of the system.

Static and stiffness models

The static model of the system is computed using the virtual work approach. It is described as:

Wτ τ τ = Γ (1)
where W = [Lsin(θ /2), -Lcos(θ /2)]. Γ denotes the torque generated by the two cables on the end effector. The components of the vector τ τ τ = [τ 1 , τ 2 ] T describe the cables tensions, computed as:

τ 1 = k∆ l 1 = k ρ 1 + 2Lcos(θ /2) -L √ 2 τ 2 = k∆ l 2 = k ρ 2 + 2Lsin(θ /2) -L √ 2 (2)
where ∆ l i denotes the elongation of the spring i. The potential energy U of the system is assumed to be only the elastic energy stored in the springs:

U = 1 2 k(∆ l 2 1 + ∆ l 2 2 ) (3)
The static equilibrium of the platform is obtained by solving ∂U ∂ θ = 0:

θ = 2 arctan L √ 2 -ρ 2 L √ 2 -ρ 1 (4)
The angular stiffness K θ at the equilibrium configuration is obtained by:

K θ = ∂ 2 U ∂ θ 2 = - ∂Γ ∂ θ = - L 2 cos(θ /2)τ 1 + sin(θ /2)τ 2 + kL 2 (5)

Dynamic model

The dynamic model of the structure is calculated using Euler-Lagrange formulation. The total kinetic energy E k of the system is defined by:

E k = 5 6 mL 2 θ 2 (6)
where m denotes the mass of one bar of the parallelogram. The equation of motion is therefore given by:

f d + f q + Γ ext = d dt ∂ L ∂ θ - ∂ L ∂ θ (7)
where Γ ext denotes the external torque applied on the end-effector and L = E T -U.

In this study, the friction and the gravity forces, respectively, f d and f q , are neglected. The structure dynamics is therefore described as:

Γ ext = 5 3 mL 2 θ + ∂U ∂ θ (8)
The dynamic of the actuator i is defined using Newton's law:

C i -C v αi -Rτ i = I m d αi dt (9)
with C i and C v being respectively the torque and viscous coefficient of the actuator. Finally, the complete dynamic model of the tensegrity mechanism is expressed in a matrix formulation as:

  RC 1 RC 2 Γ ext   =   I m 0 0 0 I m 0 0 0 5mL 2 3     ρ1 ρ2 θ   +   C v ρ1 C v ρ2 0   +   R 2 τ 1 R 2 τ 2 ∂U ∂ θ   (10)

Proposed control strategy

A control approach that allows to asynchronously track low/high dynamic trajectories and modulate the stiffness of the tensegrity mechanism is hereafter presented. It consists on a linear model predictive controller based on a discrete-time state-space model.

Dynamic model linearization

The linearization of ( 10) is performed around an equilibrium point (ρ 10 , ρ 20 , θ 0 ).

Only ∂U/∂ θ , τ 1 and τ 2 present nonlinear terms. The linearized terms are computed using the first-order Taylor expansion as:

τ 1 τ 2 = kρ 1 -f 1 θ + f 3 kρ 2 + f 2 θ + f 4 ( 11
) ∂U ∂ θ = -f 1 ρ 1 + f 2 ρ 2 -f 5 θ -f 6 (12)
such that:

f 1 = kLsin(θ 0 /2) f 2 = kLcos(θ 0 /2) f 3 = 2 f 2 + f 1 θ 0 -kL √ 2 f 4 = 2 f 1 -f 2 θ 0 -kL √ 2 f 5 = -L √ 2/2( f 1 + f 2 ) + f 1 ρ 20 /2 + f 2 ρ 10 /2 f 6 = L √ 2( f 2 -f 1 ) + {L √ 2/2( f 1 + f 2) -( f 1 ρ 20 /2 + f 2 ρ 10 /2)}θ 0
By defining the state vector x = [ρ 1 , ρ 2 , θ , ρ1 , ρ2 , θ ] T and the control vector as the torque of the motors u = [C 1 , C 2 ] T and by using the linearized expressions, the continuous-time state-space is expressed as:

ẋ = A c x + B c u + d c (13)
with:

A c =          0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 -kR 2 I m 0 R 2 f 1 I m -C v I m 0 0 0 -kR 2 I m -R 2 f 2 I m 0 -C v I m 0 3 f 1 5mL 2 -3 f 2 5mL 2 3 f 5 5mL 2 0 0 0          B c =         0 0 0 0 0 0 R I m 0 0 R I m 0 0         d c =          0 0 0 R 2 f 3 I m R 2 f 4 I m 3 f 6 5mL 2          (14)
The linear discrete-time model is defined as:

x k+1 = A k x k + B k u k + d k (15)

MPC formulation

The objective of the MPC controller is to minimize the error E k over an horizon H p between the desired trajectory T k = [x d k+1 , x d k+2 , . . . , x d k+H p ] T and the predicted output Y k = [x k+1 , x k+2 , . . . , x k+H p ] T . The cost function is defined as:

V k = E T k QE k + ∆ ∆ ∆ U T k R∆ ∆ ∆ U k (16)
with

E k = T k -Y k . Q and R are diagonal weight matrices.
In order to solve the predictive control problem, it is necessary to predict the future states x k+i , i = 1, ..., H p . To do so, we consider that over a horizon H p the linear time varying system (15) can be approximated by a linear invariant system. This yields to consider the matrices A k , B k and d k as constants, which mean that

A k+i = A k = A, B k+i = B k = B and d k+i = d k = d for i = 1, ..., H p .
With this as-sumption the future states can be computed by applying H p times the equation [START_REF] Wroldsen | Modelling and control of non-minimal non-linear realisations of tensegrity systems[END_REF]:

Y k =      x k+1 x k+2
. . .

x k+H p      =      A A 2 . . . A H p      Ψ x k +      B AB + B . . . ∑ H p -1 i=0 A i B      γ u k-1 +      d Ad + d . . . ∑ H p -1 i=0 A i d      D + . . . • • • +      B 0 . . . 0 AB + B B 0 . . . 0 . . . . . . ∑ H p -1 i=0 A i B ∑ H p -2 i=0 A i B . . . ∑ H p -H u i=0 A i B      Θ      ∆ ∆ ∆ u k ∆ ∆ ∆ u k+1 . . . ∆ ∆ ∆ u k+H u -1      ∆ ∆ ∆ U (17)
To take into account the unilateral tensions in the cables, a constrained optimization problem has to be considered. As the control variables of the MPC controller are the torques of the actuators, the tension limits, τ τ τ min = [τ min , τ min ] T and τ τ τ max = [τ max , τ max ] T , have to be converted into torques limits u min = [C min , C min ] T and u max = [C max , C max ] T . The solution of ( 16) can be obtained by solving the constrained quadratic programming (QP) problem:

min ∆ ∆ ∆ U ∆ ∆ ∆ U T H∆ ∆ ∆ U + G T ∆ ∆ ∆ U (18)
subject to :

u min ≤ u ≤ u max (19)
with

H = Θ T QΘ + R and G = 2Θ T Q(T k -Ψ x k -γu k-1 -D)
. This problem has to be solved at each step k. Once the optimal solution ∆ ∆ ∆ U is found, the first control output is applied and kept until the next time step.

Trajectory tracking

The overall control approach is composed of an external MPC controller and an internal tension control loop (cf. Fig2). The control output u k of the MPC is converted into a desired tension τ τ τ d . This vector of tension is then compared to τ τ τ m that contains the real tensions in the cables. These tensions are measured using force sensors. A tension controller C τ is considered to control these tensions in the inner loop. The output of this tension controller is sent to the actuators as desired velocities α α α d = [ αd 1 , αd 2 ] T . During the trajectory tracking, the variables θ and θ in the state vector x are to be controlled. To do so, the weight matrix Q has to be defined with high values for the components Q θ and Q θ and values near to zero for the other components

Q ρ 1 , Q ρ 2 , Q ρ1 . and Q ρ2 : Ψ m C Force sensors d - d α . τ τ Δ✉ τ MPC Optimizer z z-1 ❳ k ( k ❚ k k z -1 ✉ k k-1 γ 1 R
Q = diag(Q ρ 1 , Q ρ 2 , Q θ , Q ρ1 , Q ρ2 , Q θ ) (20)

Stiffness modulation

The stiffness modulation is possible by changing the control parameter Q. During the modulation of the angular stiffness K θ , the mechanism has to conserve a static position θ d . As K θ depends on the cables tensions ( 5), an initial and a final desired tensions, respectively τ τ τ

i = [τ i 1 , τ i 2 ] T and τ τ τ f = [τ f 1 , τ f 2 ]
T , have to be defined. These tensions are computed by solving [START_REF] Fuller | Tensile-integrity structures[END_REF]. The solution has the form [START_REF] Gouttefarde | A Versatile Tension Distribution Algorithm for n-DOF Parallel Robots Driven by n+2 Cables[END_REF]:

τ τ τ = W + Γ + Nλ (21)
where W + denotes the Moore-Penrose pseudoinverse of W, N = [-Lcos(θ /2), -Lsin(θ /2)] T a basis of W null space and λ is a scalar to be defined. Since the mechanisms is in an equilibrium position, the generated torque produced by the cables is equal to zero Γ = 0. The bounds of λ , defined as [λ min , λ max ], can therefore be obtained by solving the inequality:

τ τ τ min ≤ Nλ ≤ τ τ τ max (22)
Once the bounds of λ are obtained, the tensions τ τ τ i and τ τ τ f could be chosen for given values of λ that satisfy (22). Then, using (2), the tensions τ τ τ i and τ τ τ f are converted into initial and final desired cables displacement, respectively ρ ρ ρ

i = [ρ i 1 , ρ i 2 ] T and ρ ρ ρ f = [ρ f 1 , ρ f 2 ] T .
A desired trajectory, from ρ ρ ρ i to ρ ρ ρ f , can therefore be defined in order to modulate K θ .

Since the trajectory of interest is now in terms of the cables displacement, the weight matrix Q has to be selected with high values for Q ρ 1 , Q ρ 2 and with small values for the rest of the components. A special attention has to be made with the component Q θ , because its weight allows to ensure that the orientation of the mechanism will not change during the stiffness modulation.

Results

The control approach is evaluated on the developed prototype shown in Fig 1 . The mechanism is actuated using 2 DC motors. Polyethylene fibre cables are deployed for motion transmission. 2 force sensors are used to measure the cables tensions. Each sensor is placed between one motor and a spring. Moreover, an optical encoder, with a resolution of 0.225 • , is used to measure the angle θ .

Trajectory tracking with low and relatively-high dynamic is first assessed. Three trajectories, going from 90 • to 45 • , are created using fifth-order polynomial equation. The duration of these trajectories are respectively 30s, 5s and 1s. The cables tensions limits are set to τ min = 4N and τ max = 10N and the control parameters of Q are defined in function of the trajectory. For the first trajectory, Q θ = 8e4, Q θ = 50. For relatively-high dynamic trajectories, Q θ = 2.1e5, Q θ = 8. The control parameters literature ones [START_REF] Boehler | From modeling to control of a variable stiffness device based on a cable-driven tensegrity mechanism[END_REF], where similar experiments have been conducted using a PID controller. Reductions of 51% and 64% in the RMSE for the two trajectories, 30s and 5s, are observed (cf. Table 1). Besides, the RMSE increases when the duration of trajectories decreases. This is mainly due the bandwidth of the tension inner loop, limited actually to 1s of time response. A solution could be to reformulate the MPC problem by changing the control vector u from motor's torques to motor's angular velocities. Thus, this reformulation eliminates the tension inner loop, and would improve the time response of the system. A final experiment is performed in order to validate the stiffness modulation of the system. This experiment is divided into two phases. In the first phase, trajectory tracking is performed during 30s, going from 90 • to 70 • . Once the trajectory tracking phase is finished, the stiffness modulation is assessed during 10s, while the position of the mechanism is kept at θ f = 70 • . τ τ τ i is defined as desired tensions at the end of the trajectory tracking. τ τ τ f is defined by setting λ = λ max . A linear displacement of the cables, going from ρ ρ ρ i to ρ ρ ρ f , is therefore defined and sent to the MPC controller. The weights variables of Q are defined as Q the stiffness modulation, K θ varies linearly from 0.98Nm/rad to 0.89Nm/rad. This variation is related to the cables tensions variation, as expected in [START_REF] Skelton | Tensegrity Systems[END_REF] and observed in Fig 4 (b andc). The evolution of the angular stiffness K θ is inversely proportional to the variation of τ τ τ.

Q ρ 1 = Q ρ 2 = Q ρ1 = Q ρ2 = 10
ρ 1 = Q ρ 2 = 7.2e5, Q θ = Q ρ1 = Q ρ2 = 1, Q θ = 10

Conclusion

This work presents a linear MPC control strategy for tensegrity mechanisms. The approach allows the position tracking and the stiffness modulation by only chang-ing the values of the control parameters Q and R. A planar cable-driven tensegrity mechanism was used to validate the performance of the proposed controller. The results have shown a static error of less than 0.225 • and good performance in terms of tracking relatively-high and low dynamics trajectories. A reduction of up to 64% on the RMSE during the trajectory tracking, with respect to literature results, was obtained. Furthermore, the MPC proposed shows a decoupled control of the position and the stiffness. Future work will be the extension of this approach for spatial tensegrity mechanisms. An interesting perspective would be the proposal of a MPC with the motor's angular velocity as the output control vector, in order to avoid the use of an inner tension control loop.

Fig. 1

 1 Fig. 1 planar mechanimt.

Fig. 2

 2 Fig. 2 Strategy of control proposed based on a MPC controller

  -3 and R = diag(9e8, 9e8) are set equal for all the trajectories. H p and H u are set equal to 40 with a sampling time of 7ms. The constrained MPC optimization problem (19) is solved using the Interior point method implemented in the EIgen quadprog solver. The results can be observed in Fig3.

Fig. 3

 3 Fig. 3 Trajectory tracking with low and relatively high dynamics.

  -3 and R = diag(7e7, 7e7). The results are depicted in Fig 4. During the trajectory tracking, the angular stiffness is quite constant. During

Fig. 4

 4 Fig. 4 Results of the trajectory tracking and stiffness modulation experiment

Table 1

 1 Tracking trajectory results for PID controller[START_REF] Boehler | From modeling to control of a variable stiffness device based on a cable-driven tensegrity mechanism[END_REF] and the MPC approach.

	Trajectory PID MPC
	30s 5s	0.97 • 0.47 • 3.6 • 1.28 •

ARK2020, 028, v1: 'Model Predictive for a Planar Tensegrity Mechanism with . . .

JR Jurado Realpe, Salih Abdelaziz, and Philippe Poignet.


Jurado Realpe

G Aiche

S Abdelaziz

P Poignet

Asynchronous and decoupled control of the position and the stiffness of a spatial RCM tensegrity mechanism for needle manipulation.*

This paper introduces a 2-DOF spatial remote center of motion (RCM) tensegrity mechanism, based on a double parallelogram system, dedicated for percutaneous needle insertion. The originality of this mechanism is its ability to be reconfigured and its capacity to perform a decoupled modulation of its stiffness in an asynchronous way. To do so, an analytical stiffness model of the robot is established, and a control methodology is proposed. A prototype of the robot is developed and assessed experimentally. The position tracking is evaluated using a 6-DOF magnetic tracker sensor showing a root mean square error less than 0.8 • in both directions of the needle guide.

I. INTRODUCTION

Tensegrity systems are increasingly attracting the interest of roboticians for their large workspace and their ability to vary its stiffness. They were introduced for the first time in 1962 by Fuller [START_REF] Fuller | Tensile-integrity structures[END_REF]. Tensegrity systems can be defined as an assembly of a set of elements subjected to compression forces and others to tensile forces, and this for all configuration of the system. The set of the compressive elements interacts with the set of tensile components to define a stable structure [START_REF] Pugh | An introduction to tensegrity[END_REF]. The fact of using tensile elements allows the introduction of cables and springs to lighten the structure. In the literature, tensegrity systems have been considered as structures [START_REF] Snelson | Continuous tension, discontinuous compression structures[END_REF], [START_REF] Motro | Tensegrity: Structural Systems for the Future[END_REF], mobile robots [START_REF] Paul | Design and control of tensegrity robots for locomotion[END_REF] or deployable systems [START_REF] Pinaud | Deployment of a Class 2 Tensegrity Boom[END_REF], [START_REF] Masic | Open-loop control of class-2 tensegrity towers[END_REF].

Research on control for tensegrity systems was initiated in the mid-1990s, with efforts dedicated to establish dynamic model of these structures [START_REF] Skelton | Tensegrity Systems[END_REF]. Today, most of the control contributions of these systems are related to disturbance rejection and motion control. Disturbance rejection refers to structures subject to deformation and vibration due in particular to thermal expansion or atmospheric fluctuations such as the wind or seismic disturbances. Motion control addresses the robots whose architecture is based on tensegrity, in particular mobile robots and deployable mechanisms.

For mobile robot and deployable mechanisms, trajectory tracking is not considered as an absolute priority. Indeed, the principal focus of these systems is the production of gait, in the case of mobile robot [START_REF] Paul | Design and control of tensegrity robots for locomotion[END_REF] and the deployment ability for deployable systems. For these latter, several open-loop control approaches have been introduced [START_REF] Pinaud | Deployment of a Class 2 Tensegrity Boom[END_REF], [START_REF] Masic | Open-loop control of class-2 tensegrity towers[END_REF] under the assumption that reconfiguration is performed in quasi-static *This work was supported by French state funds managed by the ANR within the Investissements d'Avenir programme (Robotex ANR-10-EQPX-44, Labex CAMI -ANR-11-LABX-0004) 1 The authors are with the LIRMM, University of Montpellier, CNRS, France. [juradareal, abdelaziz, poignet] @lirmm.fr.

while neglecting dynamic effects. Besides, dynamic has also been taken into account in the development of open-loop control laws, considering in particular the dynamic model as a constraint in the optimization procedure to minimize the execution time [START_REF] Sultan | Modeling, design, and control of tensegrity structures with applications[END_REF] and the energy required for deployment [START_REF] Aldrich | Time-energy optimal control of hyper-actuated mechanical systems with geometric path constraints[END_REF].The problem of rapid tensegrity control has been also addressed in [START_REF] Sultan | Tensegrity Deployment Using Infinitesimal Mechanisms[END_REF] by using infinitesimal mechanisms and non-linear feedback control.

For the locomotion control of mobile tensegrity robots, open-loop control strategies based on genetic and evolutionary algorithms [START_REF] Paul | Design and control of tensegrity robots for locomotion[END_REF], [START_REF] Fujiia | Movement control of tensegrity robot[END_REF], [START_REF] Sunspiral | Super Ball Bot -Structures for Planetary Landing and Exploration[END_REF] and neural networks [START_REF] Kanchanasaratool | Modelling and control of class NSP tensegrity structures[END_REF] have been developed to overcome complex dynamic modeling and precise knowledge of the system's physical environment. Other strategies, this time in a closed loop, have been proposed in order to minimize vibrations during the movement of the structure [START_REF] Van De Wijdeven | Shape change of tensegrity structures : design and control[END_REF] by using on the one hand a H ∞ controller based on a linearization of the dynamic model and on the other hand to compensate for inaccuracies of the dynamic model [START_REF] Aldrich | Control synthesis for a class of light and agile robotic tensegrity structures[END_REF] by using a non-linear decoupling control.

Although these control strategies are of interest to improve the performance of tensegrity systems in terms of disturbance rejection, locomotion and deployment capacity, we are interested here in controlling these systems from a very different perspective. Indeed, our motivation is to control the reconfiguration and trajectory tracking of tensegrity systems, considered in this work as manipulators, while explicitly managing their stiffness. This paper is the continuation of the work of [START_REF] Boehler | From modeling to control of a variable stiffness device based on a cabledriven tensegrity mechanism[END_REF] where a proof of concept of this motivation has been demonstrated through a planar 1-degree of freedom (DOF) mechanism. In this work, this proof of feasibility is evaluated through a new spatial 2-DOF tensegrity mechanism, where analytical stiffness model of the mechanism is derived. Decoupling to control the position or the stiffness along one of the mobilities without affecting the other mobility is also addressed. This paper is organized as follow. In the second section, the requirements and an overview of the mechanism is presented. Section 3 introduces the direct and inverse kinematic models. In section 4, the static and the stiffness models are derived. Section 5 focuses on the device control of the prototype where results on position tracking and stiffness modulation are discussed. Conclusions and future work are developed in section 6. A. Requirements Robotic assistants developed for percutaneous needle insertions [START_REF] Melzer | INNOMOTION for Percutaneous Image-Guided Interventions[END_REF], [START_REF] Bricault | Light Puncture Robot for CT and MRI Interventions[END_REF], [START_REF] Arnaud | Toward Unibody Robotic Structures with Integrated Functions using Multimaterial Additive Manufacturing: Case Study of an MRI-compatible Interventional Device[END_REF], [START_REF] Abdelaziz | Controlling the Stormram 2: An MRI-compatible robotic system for breast biopsy[END_REF] often use rigid structures and target static tumors (prostate, breast, neurological targets, ,...). These mechanisms are not adapted for interventions on organs subjected to disturbances related to physiological motions, such as the liver. To perform a liver biopsy for instance, a needle is inserted in the abdomen. This procedure requires to define a fixed needle insertion point as well as a correct orientation of the needle before its insertion. A maximum needle orientation of 45 • around the insertion point could be necessary to reach the target [START_REF] Kaufman | Vascular and Interventional Radiology[END_REF], [START_REF] Walsh | A Patient-Mounted, Telerobotic Tool for CT-Guided Percutaneous Interventions[END_REF]. Moreover, the procedure requires to hold stiffly the needle during the insertion (patient in apnea) so as to ensure precision. It requires also to release the needle, which moves freely with the organ in movement when the patient breathes. This allows to avoid organ laceration. To perform this gesture using a robotic assistance, a mechanism with a minimum 2 rotational DOF is required to perform the needle orientation. A system with RCM corresponding to the insertion point is also necessary. Lastly, the mechanism should allow stiffness modulation in order to avoid organ laceration.

B. System description

The system, as illustrated in Fig 1, is a variable stiffness 2-DOF tensegrity mechanism with a RCM. It is is composed of two equilateral parallelograms placed at 90 • to each other. The upper bars of the two parallelograms are linked to the needle guide using a sliding pivot. Moreover, the two parallelograms are elevated with an angle α = 15 • with respect to the horizontal plane so as to allow the definition of a common point between the RCM and the needle insertion point. The way of performing the RCM is inspired from [START_REF] Baumann | The PantoScope: a spherical remote-center-of-motion parallel manipulator for force reflection[END_REF] with a modification of the way of actuation. Instead of actuating each parallelogram by one actuator, a redundancy of actuation is performed so as to allow stiffness modulation. Each parallelogram is therefore actuated through a pair of cables. Each cable is attached to a corner of the parallelogram at one end and to a spring before rolling up on the actuated pulley on the other end. The mechanism is hence actuated by 4 rotary actuators and has 2 mobilities which are the two orientations. The 2 degrees of kinematic redundancy exploit the prestress property of the tensegrity structures in order to modulate the stiffness of the mechanism [START_REF] Azadi | Variable antagonistic stiffness element using tensegrity mechanism[END_REF]. The orientation of the parallelograms are measured using an optical encoder, as shown in Fig 1 . A magnetic tracker is also placed inside the needle guide in order to validate the needle guide orientation.

III. KINEMATIC MODELS

The center of the reference frame R 0 = (O, x 0 , y 0 , z 0 ) is defined as the mechanism RCM.

R 1 = (O, x 1 , y 1 , z 1 ) is a fixed frame defined such that 0 R 1 = rot(y 0 , -α) (cf. Fig IIIa). The fixed frame R 1 = (O, x 1 , y 1 , z 1
) is defined such that x 1 and z 1 axes are in the plane (y 0 z 0 ) and that the angle between x 1 and y 0 is equal to α. The fixed frames

R 2 = (O, x 2 , y 2 , z 2 ) and R 2 = (O, x 2 , y 2 , z 2 ) are defined respectively such that 1 R 2 = rot(x 1 , β 0 ) and 1 R 2 = rot(x 1 , -β 0 ).

A. Direct kinematic model

The direct kinematic model allows to express the orientation of the needle guide, defined by X = [η µ 0] T , according to the joint variables β = [β 1 β 2 ] T , as illustrated in Fig III (b) and (c). To do this, the orientation of the frame R n is expressed in the reference frame R 0 :

0 R n = 0 R 1 1 R 3 3 R n = rot(y 0 , -α)rot(x 1 , β 1 )rot(y 3 , φ 1 ) (1)
The angle φ 1 can be computed by solving the closure equation [START_REF] Baumann | The PantoScope: a spherical remote-center-of-motion parallel manipulator for force reflection[END_REF]:

2 R 1 1 R 3 3 R n z = 2 R 2 2 R 1 1 R 3 3 R n z rot(x 1 , β 1 )rot(y 3 , φ 1 )z = rot(z 2 , δ)rot(x 1 , β 2 )rot(y 3 , φ 2 )z (2)
The expression of φ 1 is therefore:

φ1 = atan2(-sin(β 1 )cos(δ) + tan(β 2 )cos(β 1 ), sin(δ)) (3) with: β 1 = β 1 -β 0 and β 2 = β 2 + β 0 .
The result of ( 1) is expressed in the form of: The orientation of the needle guide is finally computed as:

0 R n =   n x o x a x n y o y a y n z o z a z   (4)
η = atan2(-a y , a z ) µ = atan2(a x , a z ) (5)

B. Inverse kinematic model

In the inverse kinematic model, the orientation of the needle guide, defined by (4), is supposed to be known. The objective is to determine the joint variables β 1 and β 2 . To do so, the equation ( 1) is pre-multiplied by rot(x 1 , -β 1 )rot(y 0 , α). From the resulting second row, third column and by replacing -a y /a z = tan(η) and a x /a z = tan(µ), we can compute β 1 :

β 1 = atan2(tan(η), -sin(α)tan(µ) + cos(α)) (6)
The angle φ 1 is computed by considering the first row, third column and third row, third column elements:

φ 1 = atan2 cos(β 1 )(cos(α)tan(µ) + sin(α)), -sin(α)tan(µ) + cos(α) (7)
The angles β 2 and φ 2 are determined by considering the closure equation ( 2):

β 2 = atan2(sin(δ)tan(φ 1 ) + cos(δ)sin(β 1 ), cos(β 1 )) (8) φ 2 = atan2 cos(δ)sin(φ 1 ) -sin(δ)sin(β 1 )cos(φ 1 ), cos(β 1 )cos(φ 1 )/cos(β 2 ) (9)

C. Workspace and singularities points

This mechanism is designed to reach all the points of a cone of 45 • around the RCM, so as to satisfy the needle orientation requirements for a liver biopsy procedure. All the geometric parameters are chosen so as to satisfy this workspace, as shown in the Fig 3.

This mechanism has three different regions of singularities. The first one is produced when θ 1 = 0. In this case, the parallelogram 1 is completely leaning forward and lines up with x 1 -axis, which means the losing of one DOF in the mechanism. The second region of singularity is defined by θ 2 = 0, that produces a similar effect. The third region is determined when the needle guide lies up in the plane (x 0 y 0 ). All these regions of singularities are outside of the reachable workspace.

IV. STATIC AND STIFFNESS MODELS

A. Static model

As aforementioned, the actuation of each parallelogram i is performed using two cables attached to the parallelogram corners D i and C i (cf. Fig 4). These cables are attached at the other end to springs of stiffness k before rolling up on actuators. The vector angular position of these actuators is denoted α = [α 11 α 21 α 12 α 22 ] T . The tensions in the cables applied to the parallelogram i are denoted τ i = [τ 1i τ 2i ] T . The static of the parallelogram i can be described as:

W i τ i = Γ i ( 10
)
where W i = [lsin(θ i /2) lcos(θ i /2)] and Γ i being the torque generated by the cables on the bar A i D i . The angle θ i can

This system of equations can be solved by multiplying the first equation by cos(φ 1 ) and the second equation by sin(φ 1 ). Subtracting the obtained equations which leads to:

d cos(φ1) -le -l cos(φ1)cos(θ1) -l sin(φ1)sin(θ1) = 0 (13)
Using an analogical reasoning although by projecting [START_REF] Sultan | Tensegrity Deployment Using Infinitesimal Mechanisms[END_REF] along x 3 and z 3 axes, θ 2 can be obtained. The angle θ i for the parallelogram i can thus be computed as:

θ i = atan2 (xz + y x 2 + y 2 -z 2 )/l 2 , (yz -x x 2 + y 2 -z 2 )/l 2 (14)
with: x = l cos(φ i ), y = l sin(φ i ) and z = d cos(φ i ) -l e

B. Stiffness model

To elaborate the stiffness model, one should compute the force contribution, generated by the tensions in the cables, of each parallelogram to the needle guide. These forces are represented in the needle guide frame R n and computed by applying the principle of virtual works in each bar of the parallelogram. The contact point between the needle guide and the parallelogram i is denoted O i and the force applied at this point is defined as f i (cf. Fig 4). These forces are expressed first in the frame attached to each parallelogram. f i ,expressed in R 3 , gives :

             3 f 1x = -τ 11 sin(θ 1 /2) + τ 21 cos(θ 1 /2) sin(θ 1 ) -tan(φ 1 )cos(θ 1 ) 3 f 1y = 0 3 f 1z = -τ 11 sin(θ 1 /2) + τ 21 cos(θ 1 /2) cos(θ 1 ) -cot(φ 1 )sin(θ 1 ) (15)
This force 3 f 1 can be expressed in the frame R n as:

n f 1 = n R 3 3 f 1 = rot(y 3 , -φ 1 ) 3 f 1 ( 16
)
After simplification by using trigonometric relations, the expression of n f 1 is defined as:

         n f 1x = -τ 11 sin(θ 1 /2) + τ 21 cos(θ 1 /2) -sin(φ 1 )cos(θ 1 ) + sin(θ 1 )cos(φ 1 ) n f 1y = 0 n f 1z = 0 (17)
Using a similar reasoning, the expression of n f 2 = [ n f 2x n f 2y n f 2y ] could be obtained. At the equilibrium configuration, the sum of the torques applied to the needle guide must be equal to zero:

Γ e + l 1 × n f 1 + l 2 × n f 2 = 0 ( 18
)
with Γ e being the external torque, that is equal to:

Γ e = n f 2y • l 2 -n f 1x • l 1 0 T ( 19
)
with

l i = d 2 + l 2 + 2 • d • l • cos(θ i ) -l 2
e being the distance between O i and the RCM (cf. Fig 4). The stiffness of the needle guide is then obtained by:

K = - dΓ e dX T = k η k µ 0 T ( 20
)
where k η and k µ represent the angular stiffnesses of the needle guide around the mobilities of the mechanism.

For ease of notation, let's define Γ ey = n f 1x • l 1 . The stiffness k µ is computed as:

k µ = dΓ ey dµ = ∂Γ ey ∂θ 1 ∂θ 1 ∂µ + ∂Γ ey ∂φ 1 ∂φ 1 ∂µ (21)
with ∂θ1 ∂µ = ∂θ1 ∂φ1 ∂φ1 ∂µ . This can be calculated using ( 7) and [START_REF] Kanchanasaratool | Modelling and control of class NSP tensegrity structures[END_REF]. The expression of k µ is reformulated as:

k µ = ( ∂Γ ey ∂θ 1 ∂θ 1 ∂φ 1 + ∂Γ ey ∂φ 1 ) ∂φ 1 ∂µ (22)
The computation of k η can be obtained using the same reasoning.

V. CONTROL A. Control Principle

The main objective is to control asynchronously the position and the stiffness of the mechanism. The control strategy (cf. Fig 5) is based on a tension distribution algorithm, inspired from [START_REF] Abdelaziz | Control of cabledriven manipulators in the presence of friction[END_REF]. This algorithm, as shown below, allows the stiffness modulation. The control strategy has an external loop to control the configuration of the parallelograms defined by θ = [θ 1 θ 2 ] T with respect to the reference

θ d = [θ d 1 θ d 2 ]
T . This reference can be computed from the desired orientation X d = [η d µ d 0] T of the needle guide using the equations of the inverse kinematic model and [START_REF] Kanchanasaratool | Modelling and control of class NSP tensegrity structures[END_REF].

The error of the external loop is multiplied by a proportional controller C θ that generates a virtual torque

Γ v = [Γ v 1 Γ v 2 ]
T . Γ v and θ are used as inputs to the tension distribution algorithm.

The output of this algorithm is a reference vector tension

τ d = [τ d 1 τ d 2 ] T . The vector τ d i = [τ d 1i τ d 2i ] T contains IK d d β d * θ d * ( *
) Equations ( 7),( 9) and ( 14). the tensions in the cables that move the parallelogram i. These tensions are supposed to generate the virtual torque Γ v i that should be applied to the parallelogram i. Besides, the algorithm has to ensure that the vector τ d i remains always positive and bounded in a defined interval [τ min , τ max ] for all the possible configurations of the mechanism. These bounds are defined as τ min = [τ min τ min ] T and τ max = [τ max τ max ] T . The tension distribution algorithm solves:

W i τ d i = Γ v i with τ min ≤ τ d i ≤ τ max (23)
The solution to (23) has the form [START_REF] Gouttefarde | A Versatile Tension Distribution Algorithm for n-DOF Parallel Robots Driven by n+2 Cables[END_REF]:

τ d i = W + i Γ v i + N i λ i (24)
W + i denotes the Moore-Penrose pseudoinverse of W i , N i = [-l cos(θ i /2) -l sin(θ i /2)
] T a basis of W i null space. λ i is a scalar to be defined that satisfies [START_REF] Walsh | A Patient-Mounted, Telerobotic Tool for CT-Guided Percutaneous Interventions[END_REF]. The bounds of λ i can be obtained by solving the inequality system:

τ min -W + i Γ v i ≤ N i λ i ≤ τ max -W + i Γ v i (25)
These bounds are defined as [λ i min , λ i max ] for each parallelogram. They are variable and depend on the configuration of the mechanism as well the tensions limit. λ i must be chosen inside this interval so as to satisfy [START_REF] Walsh | A Patient-Mounted, Telerobotic Tool for CT-Guided Percutaneous Interventions[END_REF].

Once τ d is computed using [START_REF] Baumann | The PantoScope: a spherical remote-center-of-motion parallel manipulator for force reflection[END_REF], it is compared to τ that contains the real tensions in the cables. These tensions are measured using force sensors. A tension controller C τ is considered to control these tensions. The output of this tension controller is sent to the actuators as desired velocities αd = [ αd

B. Stiffness Modulation

To simplify the notations in [START_REF] Kaufman | Vascular and Interventional Radiology[END_REF] that expresses the angular stiffness k µ produced by the first parallelogram, let's define:

a 1 = ∂Γey ∂θ1 , a 2 = ∂θ1
∂φ1 , a 3 = ∂Γey ∂φ1 and a 4 = ∂φ1 ∂µ . which yields:

k µ = (a 1 a 2 + a 3 )a 4 (26)
Assuming that the system is in an equilibrium configuration Γ v = 0 and that the desired tensions are reached (τ = τ d ), the equation ( 24) for the parallelogram 1 becomes:

τ 1 = N 1 λ 1 (27)
It can be shown that, at the equilibrium configuration, a 1 = -l1 2sin(θ1+φ1) τ T 1 N 1 and a 3 = 0. Hence:

k µ = -a 2 a 4 l 1 2sin(θ 1 + φ 1 ) τ T 1 N 1 (28)
Pre-multiplying ( 27) by N T 1 yields:

N T 1 τ 1 = N T 1 N 1 λ 1 (29)
Substituting ( 29) in (28) leads to:

k µ = - l 1 a 2 a 4 l 2 λ 1 2sin(θ 1 + φ 1 ) (30)
The relationship between k η and λ 2 can be established by applying similar reasoning. Through this result, it is therefore possible to modulate explicitly the stiffness of the mechanism by varying λ i . It is also important to highlight that it is possible to modulate stiffness only along one mobility of the needle guide without affecting the other mobility.

C. Setup and experimentation

In order to validate the control strategy proposed in the previous subsection, a prototype of the mechanism is developed, as shown in Fig 1, and an experiment is conducted. The mechanism is actuated using 4 DC motors, two for each parallelogram. Polyethylene fibre cables are deployed for motion transmission. 4 force sensors are used to measure the cable tensions. Each sensor is placed between one motor and a spring. Moreover, two optical encoders are used to measure the angles β 1 and β 2 . Finally, a magnetic tracker with a resolution of 0.5 • is placed inside the needle guide in order to validate the needle guide orientation obtained using the developed kinematic models.

The experiment is divided into two phases. In the first phase, t ∈ [0, 20]s, the position tracking performance is evaluated. In the second phase, t ∈ [START_REF] Arnaud | Toward Unibody Robotic Structures with Integrated Functions using Multimaterial Additive Manufacturing: Case Study of an MRI-compatible Interventional Device[END_REF]40]s, the stiffness modulation is assessed. During the first phase, a fifth-order polynomial equation is used to generate the cartesian trajectory along the two mobilities: η d from 0 • to -20 • and µ d from 0 • to 20 • . The admissible tensions in the cables are set between τ min = 4N and τ max = 20N . Finally, the values of λ 1 and λ 2 are chosen respectively at their lower admissible values λ 1 min and λ 2 min so as to reduce the transversal forces, f 1 and f 2 , applied to the needle guide during the positioning. During the second phase, a linear variation of λ 1 and λ 2 is performed, passing from λ i min to λ i max for each parallelogram i. This variation allows to evaluate the capacity to modulate the angular stiffness of the needle guide.

D. Results

Three different curves can be observed in Fig 6 (ab). These curves correspond to the desired angles ( d ), the estimated ( e ) ones using the optical encoders as well as the direct kinematic model, and finally the measured ones ( m ) using the magnetic tracker. During the position tracking (Fig 6, a-b), a dead-zone of about 6s, caused by static frictions, is observed. After this zone, the tracking improves with a root mean square error, respectively, of 1.41 • between η d and η e , and of 1.47 • between µ d and µ e . The static error in both trajectories is less than 2%. Besides, a good tension tracking is observed during the two phases, as seen in Fig 6 (e-f). The robustness of the control strategy with respect to the uncertainty of the kinematic parameters is also analyzed. A variation of + -10% of these parameters, does not present changes in the position static error. Moreover, the root mean square error of the tracking is less than 0.1 • compared to the results obtained without parameters variation.

In order to validate the models of the mechanism, a comparison between the estimated results and the measured ones is realized. A root mean square error of about 0.8 • along µ and η directions is observed, knowing that the precision of the tracker is 0.5 • . This error could be explained by the use of a rigid model in the estimation of the needle guide orientation, whereas the prototype is 3D printed, therefore presents some flexibilities. The mechanical backlashes have also an impact on this error. Besides and as it can be observed in the results, a stick and slip phenomenon occurs during the position tracking. This phenomenon is mainly caused by static frictions. A mechanical improvement of the sliding joint would reduce this frictions and attenuate the dead zone present at the beginning of the trajectories.

Regarding the stiffness modulation, by varying λ 1 and λ 2 the angular stiffnesses, k η and k µ , vary gradually as seen in (Fig 6, c-d). The stiffnesses variation is consistent with the cable tensions variation, as it can be observed in Fig 6 (e-f). A maximal angular stiffness of 1.06Nm/rad is reached along the two mobilities. Theoretically, no motions of the mechanism should occur during the stiffness modulation. However, a 10% position error is observed. This can be explained by the increasing forces f 1 and f 2 , that are due to the cable tensions augmentation.Indeed, these forces produce a small deflection in the sliding joint, given its flexibility. A rigidity improvement of the sliding joint should reduce this orientation error.

VI. CONCLUSION

In this paper, a 2-DOF spatial remote center of motion cable-driven mechanism is proposed. Based on a tensegrity system, this mechanism allows to control asynchronously the reconfiguration while modulating explicitly its stiffness. The stiffness model is derived analytically and linked to the tension distribution algorithm for a decoupled stiffness modulation.

Trajectory tracking is assessed on the developed prototype obtaining a static error less than 2% between the estimated and reference orientations. An evaluation of the position tracking has been performed using a magnetic tracker. The results show that the root mean square error, between the estimated orientation and the measured one is about 0.8 • , knowing that the resolution of the magnetic tracker is 0.5 • .

The experiment has shown also that it is possible to explicitly modify the stiffness of the mechanism by a factor of 3.5. A position static error of 10% has been observed during the stiffness modulation. To our sense, these results are satisfactory, promising and show that it is possible to vary the stiffness of the needle guide without affecting roughly its position.

Future work will be focused in the mechanical improvement of the mechanism, particularly in the re-design of the sliding joint. New control strategies based on friction compensation would be of interest to reduce the stick and slip phenomenon and the dead zone observed in the trajectory tracking. An Interesting perspective would be the implementation of a set up, with sensors placed in the parallelograms, in order to validate the stiffness modulation of the system.

Fig. 1 .

 1 Fig. 1. System overview.

  z 2 and z 2 axes are coincident and are defined as the common normal to x 1 and x 1 axes. The angle between x 1 and x 1 is denoted δ. R 3 = (O, x 3 , y 3 , z 3 ) is a mobile frame attached to the first parallelogram and is defined such that 1 R 3 = rot(x 1 , β 1 ). The same definition applies for R 3 = (O, x 3 , y 3 , z 3 ), a mobile frame attached to the second parallelogram:1 R 3 = rot(x 1 , β 2 ).Finally, the mobile frame R n = (O, x n , y n , z n ) is attached to the needle guide. z n axis is defined along the needle axis. The orientation of R n with respect to R 3 is 3 R n = rot(y 3 , φ 1 ), with φ 1 the angle between z 3 and z n axes, as illustrated in Fig IIIb.The same definition applies for R 3 such that 3 R n = rot(y 3 , φ 2 ).

Fig. 2 .

 2 Fig. 2. (a) Schematic of the robot with the frames placement and principals angles ; (b) Planar view of the first parallelogram of the mechanism ; (c) Definition of the orientation of the needle guide.

ig. 3 .

 3 Reachable workspace of the mechanism.

Fig. 4 .

 4 Fig. 4. Overview of the forces applied to the first parallelogram.

Fig. 5 .

 5 Fig. 5. Position/stiffness control strategy.

11 αd 21 αd 12 αd 22 ]

 11211222 T . The tension inner loop must have a wider bandwidth compared to the position outer loop in order to ensure stability.

Fig. 6 .

 6 Fig. 6. Tracking trajectory and stiffness modulation: Left η variation and first parallelogram results and Right µ variation and second parallelogram results