Bernard Bonnard
email: bernard.bonnard@u-bourgogne.fr

J.-B Caillau

Metrics with equatorial singularities on the sphere *

Keywords: two-sphere of revolution, almost and sub-Riemannian metrics, cut and conjugate locus MSC classification. 53C17, 49K15

published or not. The documents may come

Introduction

In the papers [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF][START_REF] Sinclair | The cut locus of a two-sphere of revolution and Toponogov's comparison theorem[END_REF] the authors investigate the cut and conjugate loci of a Riemannian metric on the two-sphere of revolution that can be written in the normal form m(ϕ)dθ 2 + dϕ 2 -where ϕ is the angle along the meridian and θ the angle of revolution-, under the assumption that m(π -ϕ) = m(ϕ) (symmetry with respect to the equator). Motivated by examples in optimal control, our aim is to extend these results to metrics of the form XR(X)dθ 2 + dϕ 2 , X = sin 2 ϕ,

where R is a rational fraction with a single pole at X = 1. The metric is Riemannian on both open hemispheres but has a singularity on the equator. Such metrics arise when one considers the time minimum control of a system ẋ(t) = u 1 (t)F 1 (x(t)) + u 2 (t)F 2 (x(t)), u 2 1 (t) + u 2 2 (t) ≤ 1, with fixed boundary conditions, x(0) = x 0 , x(t f ) = x f , and it is well known that minimizing time amounts to minimizing the length of the curve t → x(t) provided the metric is defined by assuming F 1 and F 2 orthonormal. For such problems, singularities occur when the two vector fields become collinear. This leads to the concept of almost-Riemannian metrics [START_REF] Agrachev | Twodimensional almost-Riemannian structures with tangency points[END_REF][START_REF] Agrachev | A Gauß-Bonnet like formula on two-dimensional almost-Riemannian manifolds[END_REF][START_REF] Bonnard | The sphere and the cut locus at a tangency point in two-dimensional almost-Riemannian geometry[END_REF][START_REF] Boscain | Lipschitz classification of almost-Riemannian distances on compact oriented surfaces[END_REF], and the analysis of such metrics is related to sub-Riemannian geometry [START_REF] Bellaïche | Sub-Riemannian geometry[END_REF][START_REF] Montgomery | A tour of subriemannian geometries, their geodesics and applications[END_REF]. If the distribution {F 1 , F 2 } is bracket generating, every pair of points can be joined by a length minimizing curve. Moreover, if there exists no abnormal trajectory, each minimizing curve is a geodesic, projection on the x-space of the Hamiltonian flow exp t -→ h where

h(z) = 1 2 H 2 1 (z) + H 2 2 (z) , z = (x, p),
and where H i (z) = p, F i (x) are the Hamiltonian lifts of the F i 's. Denoting exp x0 the exponential mapping

exp x0 (t, p 0 ) = Π exp t -→ h (x 0 , p 0 )
where Π is the projection (x, p) → x, the cut and conjugate loci are defined as in the Riemannian setting; the cut locus is the set of points where geodesic curves fail to be minimizing, while the conjugate locus is the image of the critical points of the exponential mapping.

In this paper, we generalize to the singular case [START_REF] Agrachev | A Gauß-Bonnet formula for contact sub-Riemannian manifolds[END_REF] the results in [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF] relating the structure of the cut and conjugate loci to the convexity of the quasi-period of the θ-coordinate: Under appropriate assumptions, the cut locus of a point is reduced to a single segment, and the conjugate locus has at most four cusps. Although similar to the Riemannian case, the singularity of the metric on the equator has two consequences; the injectivity radius is zero and the singularities of the conjugate locus of an equatorial point differ. In Section 2, two examples motivating this study are presented. One stems from quantum mechanics, the other from space mechanics. Both are limit cases of optimal control problems not linear but affine in the control. Section 3 is devoted to the integrability properties of the geodesics of (1), and preliminary computations of the period and quasi-period of coordinates ϕ and θ are made. The optimality status of the geodesics is studied in Section 4 where the main results on the structure of cut and conjugate loci are established. The second example of §2 is given a full treatment in Section 5. In order to be able to apply §4 results, a detailed analysis using a parameterization of geodesics by elliptic curves is presented. The last section accounts for concentration of curvature and collapsing phenomena that allow to reinterpretate the effect of the equatorial singularity on the metric.

Preliminaries

We consider the metric (1) on S 2 where (θ, ϕ) are the coordinates induced by the covering of the sphere minus poles R × (0, π) (θ, ϕ) → (sin ϕ cos θ, sin ϕ sin θ, cos ϕ).

The function R is a rational fraction with a single pole of order p at X = 1, R(X) = p k=0 a k (1 -X) k , a 0 , . . . , a p-1 ≥ 0, a p > 0 (p ≥ 1).

(

) 2
We assume that the normalization condition a 0 + • • • + a p = 1 holds. The metric is Riemannian on hemispheres with a singularity on the equator when ϕ = π/2 since R(X) → ∞ when X → 1. It is a standard fact that Riemannian geometry can be recast in the Hamiltonian setting of optimal control: On each open hemisphere, finding a length minimizing curve connecting two points x 0 , x f , is equivalent to finding a measurable control function u : [0, t f ] → R 2 such that almost everywhere

ẋ(t) = u 1 (t)F 1 (x(t)) + u 2 (t)F 2 (x(t)), u 2 1 (t) + u 2 2 (t) ≤ 1, (3)
x(0) = x 0 , x(t f ) = x f , and such that the final time t f is minimized. In the previous coordinates, the vector fields are (because of the topology of the 2-sphere, one cannot provide an orthonormal frame of vector fields globally defined on S 2 )

F 1 (θ, ϕ) := 1 sin ϕ 1/R(X) ∂ ∂θ , F 2 (θ, ϕ) := ∂ ∂ϕ (4)
as Lemma 1. 1/R(X) has a smooth square root.

Proof. One has 1 R(X) = cos 2p ϕ a p + • • • + a 0 cos 2p ϕ and a p + • • • + a 0 cos 2p ϕ ≥ a p > 0 by virtue of (2).
This optimal control formulation makes sense on the whole sphere: This is how (1) will be understood as a metric defined on all S 2 in this paper. Geometrically, the singularity R(1) = ∞ forces length minimizing curves to be vertical (directed along meridians) when crossing the equator. Remark 1. The singularity here is not a degeneracy of the Riemannian tensor on the tangent space as in, e.g., [START_REF] Pelletier | Sur le théorème de Gauß-Bonnet pour les pseudo-métriques singulières[END_REF] but a degeneracy of a positive tensor on the cotangent space.

An appropriate algebraic setting for such a metric is the notion of generalized sub-Riemannian metrics (see [START_REF] Bellaïche | Sub-Riemannian geometry[END_REF][START_REF] Jean | Sub-Riemannian Geometry[END_REF]) or almost-Riemannian metrics (see [START_REF] Agrachev | A Gauß-Bonnet formula for contact sub-Riemannian manifolds[END_REF][START_REF] Agrachev | Twodimensional almost-Riemannian structures with tangency points[END_REF][START_REF] Agrachev | A Gauß-Bonnet like formula on two-dimensional almost-Riemannian manifolds[END_REF][START_REF] Bonnard | The sphere and the cut locus at a tangency point in two-dimensional almost-Riemannian geometry[END_REF]). Let g 0 be the round metric on the sphere; to (1) is associated the following morphism of fibre bundle:

1 f : (T S 2 , g 0 ) → T S 2
1 That is id-morphism, according to [START_REF] Bourbaki | Variétés différentielles et analytiques[END_REF] terminology: Base points are unchanged by the morphism, and fibres are linearly sent to fibres. such that f is the identity on the fibres above the poles and

f (θ,ϕ) : ∂ ∂θ → 1/R(X) ∂ ∂θ , ∂ ∂ϕ → ∂ ∂ϕ
otherwise. This morphism induces an application between vector fields on the sphere,

f * : Γ(T S 2 ) → Γ(T S 2 ); if F ∈ Γ(T S 2 ), f * (F )(x) := f x (F (x)), x ∈ S 2 .
Although the direct image f (T S 2 ) is not a fibre bundle,

∆ := f * (Γ(T S 2 )) ⊂ Γ(T S 2 )
is a well defined submodule over smooth functions. For x ∈ S 2 and v ∈ ∆ x := {F (x), F ∈ ∆}, define

g x (v) := inf{g 0x (u) | f x (u) = v}.
Outside poles,

g x (v) = inf{u 2 1 + u 2 2 | u 1 F 1 (x) + u 2 F 2 (x) = v}.
Then, given two points x 0 , x f on the sphere, set

d(x 0 , x f ) := inf t f 0 g x(t) ( ẋ(t)) dt
where the infimum is taken over all Lipschitz trajectories x such that ẋ(t) ∈ ∆ x(t) , t ∈ [0, t f ] (a.e.)

x(0) = x 0 , x(t f ) = x f .
The set of such horizontal trajectories is not empty as is clear from the proof of Proposition 1. d defines a complete distance on S 2 that induces the usual topology on the sphere.

To prove this fact, one defines recursively the flag associated with ∆ by means of the Lie bracket of vector fields,

∆ 1 := ∆, ∆ k+1 := ∆ k + [∆, ∆ k ] with [∆, ∆ k ] := {[F, G] | F ∈ ∆, G ∈ ∆ k }, k ≥ 1. Lemma 2. For all x ∈ S 2 , ∆ p+1 x = T x S 2 .
Proof. Outside the equator, F 1 and F 2 have rank two, so the verification is restricted to ϕ = π/2 where F 1 vanishes. Since R has an order p pole,

(ad p F 2 )F 1 (θ, ϕ) = d p dϕ p [(1/ sin ϕ) 1/R(X)] =0 at ϕ=π/2
∂ ∂θ so brackets of length at most p+1 span everywhere the whole tangent space.

Proof of Proposition 1. According to lemma 2, ∆ is bracket generating so there exist, for any pair of points on the sphere, Lipschitz curves almost everywhere tangent to ∆ connecting them by Chow-Rashevsky theorem. On the compact manifold S 2 , Filippov theorem [START_REF] Agrachev | Control Theory from the Geometric Viewpoint[END_REF] then asserts existence of time minimizing curves among those. That the metric induces the canonical manifold topology is another consequence of Chow-Rashevsky [START_REF] Jean | Sub-Riemannian Geometry[END_REF].

Remark 2. As explained in [START_REF] Agrachev | Twodimensional almost-Riemannian structures with tangency points[END_REF], generic (in the sense of Whitney) almost-Riemannian metrics on the sphere are such that, for all x ∈ S 2 , either ∆ x (regular point), ∆ 2

x (Grušin point) or ∆ 3

x (tangency point) span the tangent space. According to Lemma 2, the situation we consider is not generic for p ≥ 3.

Motivating examples

The following Hamiltonian on S 2 is considered in [START_REF] Bonnard | Energy minimization in two-level dissipative quantum control: The integrable case[END_REF]:

H(θ, ϕ, p θ , p ϕ ) = -δ cos ϕ sin ϕ p ϕ + 1 2 p 2 θ tan 2 ϕ + p 2 ϕ . (5)
It originates in quantum mechanics and partly describes the energy minimum control of a spin 1/2 particle in a magnetic field. When the parameter δ is zero, this Hamiltonian corresponds exactly to the metric (1) obtained for R(X) = 1/(1 -X) (geodesics are integral curves of an appropriate quadratic Hamiltonian-see the beginning of §3). As explained in §4, a local model near the singularity ϕ = π/2 is the so-called Grušin metric on R 2 [14, 22]

dx 2 + dy 2 x 2 •
For this reason, the metric defined by (5) when δ = 0 is called Grušin metric on the sphere. In contrast to the analogous metric on the plane, it has peculiarities (e.g., meridional cusps of conjugate loci, see §4) due to the topology of the sphere. Like the second example, it is also connected with space mechanics (see [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF]). Consider a controlled dynamical system of the form

dx dl (l) = m i=1 u i (l)F i (l, x(l)) (6)
on a smooth n-dimensional manifold X where the F i : R × X → X are vector fields parameterized by l and periodic, F (l, x) = F (l + 2π, x) for all l ∈ R and x ∈ X. This dynamics is actually a particular case of an (autonomous) affine in the control system, as we may define the vector fields

F0 (x) := ω(x) ∂ ∂l , Fi (x) := ω(x) n j=1 dx j , F i ∂ ∂x j , i = 1, . . . , m, with x := (l, x) ∈ X := R × X and write dx dt (t) = F0 (x(t)) + m i=1 u i (t) Fi (x(t)).
The pulsation ω is any 2π-periodic in l positive function on X relating the two times of the system, t and l (the latter being understood as an angular length). The periodic vector fields F i induce vector fields on S 1 × X and it is appealing to use averaging [START_REF] Moser | Regularization of Kepler's problem and the averaging method on a manifold[END_REF] to define

F i (x) := 1 2π 2π 0 F i (l, x) dl.
We consider instead the Hamiltonian associated with the L 2 dt-minimization of the control,

H(l, x, p) = ω(l, x) 2 m i=1 p, F i (l, x) 2 .
Then,

H(x, p) := 1 2π
2π 0 H(l, x, p) dl remains a positive quadratic form in the adjoint variable, p. Under the standard assumption that the Fi , i = 0, . . . , m, are bracket generating, one expects the rank of this quadratic form to be maximum and equal to n. The case being, one can find n independent vector fields F i such that

H(x, p) = 1 2 n i=1 p, F i (x) 2 .
Singularities of two types may exist though. First, such implicitly defined vector fields need not be smooth on X, even in the analytic situation. 2 In addition, the rank of the form may not be constant on X and drop at some points. This phenomenon accounts for the existence of an equatorial singularity in the following case. On X = R * + × D, D being the open unit Poincaré disk, set

F 1 := - 3(1 -e 2 )w n 1/3 (1 + e cos v) 2 ∂ ∂n + 2(1 -e 2 ) 2 n 4/3 (1 + e cos v) 2 w (e + cos v) ∂ ∂e + sin v e ∂ ∂θ with v := l -θ, w := 1 + 2e cos(l -θ) + e 2 ,
and

ω(l, x) := n(1 + e cos(l -θ)) 2 (1 -e 2 ) 3/2 •
Here, x = (n, e, θ) ∈ R * + × R * + × R are coordinates on the positive line times the pointed disk. They are used in space mechanics to represent the geometry of plane elliptic trajectories in the controlled two-body problem [START_REF] Bonnard | Energy minimization of single-input orbit transfer by averaging and continuation[END_REF]: n is the mean motion (that is n = a -3/2 where a is the semi-major axis), e is the eccentricity, and θ is the argument of the pericenter. The averaging procedure just described leads to [START_REF] Bonnard | Energy minimization of single-input orbit transfer by averaging and continuation[END_REF] H(n, e, θ, p n , p e , p θ ) = 9n The integration of H may so be performed by integrating h. This two-dimensional Hamiltonian (which retains the linear first integral p θ ) is lifted from the Poincaré disk to S 2 using the following compactification: In the covering (θ, ϕ) ∈ R × (0, π) of the sphere minus poles where

e = sin ϕ 1 + cos 2 ϕ , one has h(θ, ϕ, p θ , p ϕ ) = 1 2 4 cos 4 ϕ sin 2 ϕ(2 -sin 2 ϕ) 2 p 2 θ + p 2 ϕ .
The rank of this quadratic form in (p θ , p ϕ ) drops from 2 to 1 when ϕ = π/2. Actually, h is exactly associated to the metric (1) with

R(X) = 1 -X/2 1 -X 2 = 1 4 1 + 2 1 -X + 1 (1 -X) 2 •
This metric with an order two equatorial singularity will be referred to as the (1, 2, 1) case (according to coefficients involved in the series) and studied in Section 5.

Remark 3. The Hamiltonian (5) in the first example can either be interpretated as stemming from an affine controlled system, or as defining a pseudo-Riemannian metric (see [START_REF] Bonnard | Energy minimization in two-level dissipative quantum control: The integrable case[END_REF]). In constrast with the sub-Riemannian case that is characterized by linearity in the control, time minimization (with a bounded control) and minimization of the L 2 -norm (with a prescribed final time) of affine systems cease to be equivalent.

Integrability properties

As time minimizing curves of the control system (3), geodesics satisfy Pontrjagin maximum principle [START_REF] Agrachev | Control Theory from the Geometric Viewpoint[END_REF][START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]: If x is a shortest time trajectory generated by the optimal control u : [0, t f ] → R 2 , there exist a nonpositive constant p 0 ≤ 0 and a Lipschitzian lift (x, p) : [0, t f ] → T * S 2 of the trajectory to the cotangent bundle such that (p 0 , p) = (0, 0) and

ẋ(t) = ∂H ∂p (x(t), p(t), u(t)), ṗ(t) = - ∂H ∂x (x(t), p(t), u(t)), t ∈ [0, t f ] (a.e.),
where

H : T * S 2 × R 2 → R, H(x, p, u) := p 0 + u 1 H 1 (x, p) + u 2 H 2 (x, p).
The H i : T * S 2 → R are the Hamiltonian lifts of the two vector fields (4),

H i (x, p) := p, F i (x) , i = 1, 2.
Besides, the following maximization condition holds,

H(x(t), p(t), u(t)) = max |v|≤1 H(x(t), p(t), v), t ∈ [0, t f ] (a.e.)
As a result, the Hamiltonian evaluated along the extremal (x, p, u) is almost everywhere equal to a constant; since the final time is free, this constant is zero. The previous relations are homogeneous in (p 0 , p) and there are two cases: Either p 0 = 0 (abnormal case), or p 0 < 0 (normal case).

Lemma 3. Abnormal trajectories are stationary equatorial curves.

Proof. Assume p 0 = 0. Then H = u 1 H 1 + u 2 H 2 has to be zero and maximized along the extremal; necessarily, the two Lipschitz functions H 1 and H 2 must vanish identically on [0, t f ] (the maximum would otherwise be positive). As a result, for any time t, p(t) is orthogonal to the span at x(t) of F 1 and F 2 . If there exists t ∈ [0, t f ] such that ϕ(t) = π/2, the span at such an x(t) is the whole tangent space, so p(t) = 0. As p is solution of the linear differential equation

ṗ(t) = - ∂H ∂x (x(t), p(t), u(t)) = -p(t) [F 1 (x(t)) + F 2 (x(t))] ,
p is identically zero, which contradicts (p 0 , p) = (0, 0). So ϕ ≡ π/2; then

θ(t) = u 1 (t) 1/R(X) | X=1 = 0
so θ is also constant.

We disregard theses trivial curves and normalize p 0 to -1.

Scholium. Geodesics of (1) are integral curves of the quadratic Hamiltonian

h(θ, ϕ, p θ , p ϕ ) := 1 2 p 2 θ XR(X) + p 2 ϕ ( 7
)
restricted to the level set {h = 1/2}.

Proof. Let u be a minimum time control, and let (x, p, u) be the associated extremal. The normal Hamiltonian -1+u 1 H 1 +u 2 H 2 is zero along the extremal, so (H 1 , H 2 ) does not vanish. Because of the maximization condition,

u(t) = (H 1 , H 2 ) H 2 1 + H 2 2 (x(t), p(t)) a.e.
Then,

ẋ(t) = u 1 (t)F 1 (x(t)) + u 2 (t)F 2 (x(t)) (8) = H 1 (x(t), p(t))F 1 (x(t)) + H 2 (x(t), p(t))F 2 (x(t)) since H 2 1 + H 2 2 = 1 on {H = 0} (that is {h = 1/2}), so ẋ(t) = ∂h ∂p (x(t), p(t)) a.e.
One similarly verifies that

ṗ(t) = - ∂h ∂x (x(t), p(t)) a.e.
Remark 4. Such geodesics are arc length parameterized as [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF] implies that

g x(t) ( ẋ(t)) = u 2 1 (t) + u 2 2 (t) = 1 a.e.
As θ is a cyclic variable (symmetry of revolution-see [START_REF] Besson | Géodésiques des surfaces de révolution[END_REF] for a general reference) of h, p θ is a linear first integral and Proposition 3. h is Liouville integrable. The coordinate ϕ is parameterized by a hyperelliptic curve of genus at most p.

Proof.

X = sin 2 ϕ, so Ẋ2 = 4X(1 -X) φ2 . On {h = 1/2}, Ẋ2 = 4X(1 -X) 1 - p 2 θ XR(X) , so Y 2 = 4(1-X)(a p +• • •+a 0 (1-X) p )[X(a p +• • •+a 0 (1-X) p )-p 2 θ (1-X) p ] (9) with Y = (a p + • • • + a 0 (1 -X) p ) Ẋ.
The right hand side of ( 9) has degree at most 2(p + 1), so the genus of the complex curve is at most p.

Set Γ(ϕ) := 1 XR(X) , ϕ ∈ (0, π).
Consider the extremal departing from ϕ 0 = 0 (π) (not a pole), θ 0 being normalized to 0 and defined by a positive p θ (the degenerate case p θ = 0 corresponding to meridians-which are the only extremals passing through the poles) and nonnegative

p ϕ 0 = 1 -Γ(ϕ 0 )p 2 θ . Along the extremal, φ first vanishes when ϕ is equal to ϕ 1 := π -Γ -1 (p -2 θ ) since Lemma 4. Γ is a strictly decreasing one-to-one mapping between (0, π/2] and R + such that Γ(π -ϕ) = Γ(ϕ). Proof. One has dΓ dX = - R(X) + XR (X) (XR(X)) 2 with R (X) = p k=1 ka k (1 -X) k+1 •
Since dX/dϕ = 2 sin ϕ cos ϕ is positive on (0, π/2), the conclusion follows.

As Γ does not vanish on (0, π/2),

1 -Γ(ϕ)p 2 θ = O(ϕ 1 -ϕ)
in the neighbourhood of π -Γ -1 (p -2 θ ), and the following integral (depending on p θ and ϕ 0 ) is well-defined,

t 1 := π-Γ -1 (p -2 θ ) ϕ0 dϕ 1 -Γ(ϕ)p 2 θ • Lemma 5.
The axial symmetry σ 1 with respect to θ(t 1 ) is an inner symmetry of the extremal.

Proof. Set θ(t) := 2θ(t 1 ) -θ(2t 1 -t), pθ (t) := p θ (t), φ(t) := ϕ(2t 1 -t), pϕ (t) := -p ϕ (2t 1 -t),
and check that new curve is still an extremal, passing through the same point of the cotangent bundle at t 1 since p ϕ (t 1 ) = 0.

Necessarily, π -Γ -1 (p -2 θ ) ≥ π -ϕ 0 , so there also exists t 1 ≤ t 1 such that ϕ(t 1 ) = π -ϕ 0 . Using the previous axial symmetry, we deduce the existence of

t 2 := 2t 1 -t 1 ≥ t 1 such that, again, ϕ(t 2 ) = π -ϕ 0 . Using now the equatorial symmetry of Γ, Γ(π -ϕ) = Γ(ϕ),
the following is clear.

Lemma 6. The central symmetry s 2 with respect to (θ(t 2 )/2, π/2) defines another extremal with the same initial condition.

Proof. Set

θ(t) := θ(t 2 ) -θ(t 2 -t), pθ (t) := p θ (t), φ(t) := π -ϕ(t 2 -t), pϕ (t) := p ϕ (t 2 -t).
The new curve is still an extremal since

θ(t) = Γ(π -φ(t))p θ = Γ( φ(t))p θ , ṗϕ (t) = 1 2 Γ (π -φ(t))p 2 θ = - 1 2 Γ ( φ(t))p 2 θ ,
and θ0 = 0 = θ 0 , φ0 = π -(π -ϕ 0 ) = ϕ 0 .
Finally denote t 3 the point such that ϕ(t 3 ) = π/2 ≤ π -ϕ 0 , and remark that the central symmetry σ 2 with respect to (θ(t 3 ), π/2) leaves the extremal invariant.

Since the axial symmetry s 1 with respect to θ = 0 obviously defines another extremal originating from the same point, we conclude that the group generated by s 1 and s 2 acts on the set of extremals with same initial condition, while the group generated by σ 1 and σ 2 defines inner symmetries of each extremal. The composition rules indicate in both cases that the underlying group is the fourorder abelian Klein group,

V = Z/2Z × Z/2Z {id, s 1 , s 2 , s 1 s 2 } {id, σ 1 , σ 2 , σ 1 σ 2 }.
Proposition 4. Given any initial condition, the Klein group acts on the set of extremals issuing from the point. It also defines inner symmetries of any extremal.

An extremal is said to be a pseudo-equator whenever φ(0) = p ϕ (0) is equal to zero, whereas the equator itself cannot be an extremal because of the singularity. Lemma 7. Every extremal that is not a meridian is a pseudo-equator.

Proof. For p θ positive and p ϕ 0 nonnegative (the other cases are deduced by symmetry), there exists ϕ 0 = Γ -1 (p -2 θ ) such that, up to a time shift, the extremal is the pseudo-equator with initial condition ϕ 0 .

Conversely, any pseudo-equator meets ϕ = π/2 as one understands from the analysis of symmetries. Taking ϕ 0 = π/2 as new initial condition and retaining the same value for p θ provides the same geodesic, up again to a time shift. As a result, rather than parameterizing extremals using both ϕ 0 -we set θ 0 = 0 thanks to the symmetry of revolution-and p θ , one may either parameterize by ϕ 0 ∈ (0, π/2) alone using the fact all geodesics (with the exception of meridians, p θ = 0) are pseudo-equators (then, implicitly, p 2 θ = 1/Γ(ϕ 0 ), ϕ 0 = π/2 since the equator is not a geodesic), or parameterize by their Clairaut constant p θ ∈ R * + , considering only the initial condition at singularity, ϕ 0 = π/2. The second point of view reduces the study of geodesics to those starting at singularity. Proposition 5. On every extremal, the coordinate ϕ is periodic with period

T (p θ ) = 4 π/2 Γ -1 (p -2 θ ) dϕ 1 -Γ(ϕ)p 2 θ ,
and θ(t + T ) = θ(t) ± ∆θ (sign depending on the sign of p θ ) with quasi-period

∆θ(p θ ) = 4 π/2 Γ -1 (p -2 θ ) Γ(ϕ)p θ dϕ 1 -Γ(ϕ)p 2 θ •
Proof. According to the previous analysis, it is enough to check the result on pseudo-equators. But then, t 1 = t 2 = t 3 = 2t 4 , so setting T := 2t 1 and using the axial symmetry with respect to θ(t 1 ) gives the result since ϕ(T ) = ϕ(0), p ϕ (T ) = -p ϕ (0) = 0 = p ϕ (0). So θ = Γ(ϕ)p θ is also periodic, which concludes the proof.

As functions of ϕ 0 ,

T (ϕ 0 ) = 4 π/2 ϕ0 dϕ 1 -Γ(ϕ)/Γ(ϕ 0 ) , (10)
and

∆θ(ϕ 0 ) = 4 π/2 ϕ0 Γ(ϕ)dϕ Γ(ϕ 0 ) -Γ(ϕ) • (11)
These relations actually cover the case of meridians p θ = 0 (i.e. ϕ 0 = 0) for which T = 2π and ∆θ = 2π (two instantaneous rotations of angle π when crossing poles at t = π and t = 2π). We end the section with the following result that will be used for §5 computations.

Proposition 6. T = p θ ∆θ .

Proof. Write as in [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF] T (p θ ) = p θ ∆θ(p θ ) + 4

π/2 Γ -1 (p -2 θ ) 1 -Γ(ϕ)p 2 θ dϕ, so T (p θ ) = ∆θ(p θ ) + p θ ∆θ (p θ ) + 4 π/2 Γ -1 (p -2 θ ) ∂ ∂p θ 1 -Γ(ϕ)p 2 θ dϕ -8(Γ -1 ) (p -2 θ )p θ 1 -Γ(ϕ)p 2 θ | ϕ=Γ -1 (p -2 θ ) 0 = p θ ∆θ (p θ ).

Cut and conjugate loci

The cut time along a geodesic x, that is along an extremal of the minimum time problem (3), is the supremum of times t such that the curve restricted to [0, t] is a shortest time trajectory between x(0) and x(t f ):

t cut := sup{t ≥ 0 | x is minimizing on [0, t]}.
When t cut < ∞, x(t cut ) is called a cut point. The set of all cut points on geodesics departing from a given intial point x 0 is the cut locus of x 0 . A separating point along the geodesic is a point x(t M ), t M > 0, such that there exists a different geodesic, y, reaching the point at the same time:

x(t M ) = y(t M ).
The exponential mapping of a fixed point

x 0 ∈ S 2 is exp x0 : R * + × T * x0 S 2 ∩ {h = 1/2} → S 2 (12)
(t, p 0 ) → x(t, x 0 , p 0 ) := Π • exp t -→ h (x 0 , p 0 )
where Π : T * S 2 → S 2 is the canonical projection and exp t -→ h the one-parameter global subgroup generated by the symplectic gradient of the quadratic Hamiltonian [START_REF] Besson | Géodésiques des surfaces de révolution[END_REF],

-→ h (x, p) = ∂h ∂p (x, p) ∂ ∂x - ∂h ∂x (x, p) ∂ ∂p •
The intersection of the fibre T * x0 S 2 with {h = 1/2} is a compact oval diffeomorphic to S 1 outside the equator, or the union two lines {p ϕ 0 = ±1} for an equatorial point. That the subgroup is globally defined in the second case comes from the analysis of Section 3: For any p θ ∈ R and p ϕ 0 = ±1, the coordinates ϕ and θ are periodic or quasi-periodic, respectively. A conjugate point along the geodesic x is a critical value of exp x(0) ; if (t c , p 0 ) is the corresponding critical point, t c is the conjugate time. When t c > 0 is the first conjugate time along the geodesic, x(t c ) is called a first conjugate point. The set of first conjugate points on geodesics departing from x 0 is the conjugate locus of x 0 . Results of optimal control ensure that local optimality is lost after conjugate points [START_REF] Sarychev | The index of second variation of a control system[END_REF]. Besides, extremals of such problems have to be smooth so that broken curves that are concatenations of minimizing geodesics cannot be minimizing, entailing that optimality cannot be preserved after separating points, such as those generated by the symmetries between extremals described in the previous section. Separating and conjugate times are so upper bounds of cut times. The following standard result remains valid in our setting with singularities. Proof. Directly follows from estimates of integrals (10-11) using the fact that Γ does not vanish identically at ϕ = π/2.

Proof of Proposition 7. Let γ(t cut ) be the cut point along the geodesic γ starting from x 0 = (θ 0 , ϕ 0 ) and generated by the adjoint vector p 0 . If the cut point is not a conjugate point, the exponential mapping is a diffeomorphism in a small enough neighbourhood V 0 of (t cut , p 0 ). Since the metric is complete, there are minimizing extremals

γ n joining x 0 to γ n (t n ) = γ(t cut + 1/n), t n < t cut + 1/n, for n ≥ 1. First assume that ϕ 0 = π/2. Then the oval h -1 (θ 0 , ϕ 0 , •)({1/2}
) is compact, and one can extract a converging subsequence of the (p 0n ) n generating the extremals γ n , and thus conclude classically (see, e.g., [START_REF] Sakai | Riemannian Geometry[END_REF]): p 0n → p 0 and γ n (t n ) → γ(t cut ); assuming by contradiction that the point is not a separating one, that is assuming that p 0 = p 0 , implies that for n large enough γ n (t n ) belongs to exp x0 (V 0 ), so that t n = t cut + 1/n, whence the contradiction. Let now ϕ 0 = π/2. Though h -1 (θ 0 , ϕ 0 , •)({1/2}) = {p ϕ 0 = ±1} is not compact anymore, the associated sequence (p θ n ) n still has to be bounded otherwise there would exist a subsequence such that |p θ n | → ∞. Assume by contradiction that this is the case. Because of the central symmetry s 2 , (∆θ(p θ n )/2, π/2) is a separating point on γ n = (θ n , ϕ n ), so t n < T (p θ n )/2 and θ n (t n ) < ∆θ(p θ n )/2; as

| φn | = |p ϕ n | ≤ 1, |ϕ n (t n ) -π/2| ≤ T (p θ n )/2 and (θ n (t n ), ϕ n (t n )) → (0, π/2)
according to Lemma 8. Since γ(t cut + 1/n) = γ n (t n ), this implies γ(t cut ) = x 0 , which is contradictory. The sequence (p θ n ) n is hence bounded, and we can conclude as previously.

Cut loci for an analytical Riemannian metric on the sphere are known to be finite trees whose extremities are conjugate points after the work of Poincaré [START_REF] Myers | Connections between geometry and topology I[END_REF][START_REF] Myers | Connections between geometry and topology II[END_REF][START_REF] Poincaré | Sur les lignes géodésiques des surfaces convexes[END_REF]. In our case, cut loci have peculiarities due to the symmetry of revolution and the equatorial singularity (see also Corollary 3 in §6; more generally on surface of revolutions, see [START_REF] Shiohama | The geometry of total curvature on complete open surfaces[END_REF]).

Theorem 1. Under the assumption that ∆θ is strictly decreasing for p θ > 0, cut loci of the metric (1) are antipodal subarcs. The cut locus of a pole is reduced to the opposite pole, is equal to the equator minus the point for an equatorial point, and to a proper closed subarc of the antipodal parallel otherwise.

Proof. The case of poles is obvious and does not depend on any assumption on ∆θ since the only extremals through them are meridians. Consider now the situation ϕ 0 = π/2, and show that the exponential mapping is injective on the quadrant

p θ >0 (0, T (p θ )/2] × {(p θ , 1)},
that is show that subarcs of extremals defined by t ∈ [0, T (p θ )/2], positive p θ and p ϕ = +1 do not intersect. If p θ > p θ , the arc associated with p θ is strictly below the one associated with p θ . Indeed, note that on the first half of such an arc t ∈ [0, T (p θ )/4) and φ does not vanish so that the curve can be parameterized by ϕ instead of time. There,

f (ϕ, p θ ) := dθ dϕ = Γ(ϕ)p θ 1 -Γ(ϕ)p 2 θ is an increasing function of p θ since ∂f ∂p θ (ϕ, p θ ) = Γ(ϕ) (1 -Γ(ϕ)p 2 θ ) 3/2 > 0.
As geodesics starting from ϕ 0 = π/2 cross again the equator at ∆θ(p θ )/2, the assumption ensures that the aforementioned subarcs do not intersect. We conclude by remarking that the full set of extremals is obtained by considering the action of the Klein group on geodesics with same initial condition (see §3).

First, the central symmetry s 2 that generates intersections at t = T (p θ )/2, then the axial symmetry s 1 with respect to θ = 0 that generates intersections at θ = π, thus not prior to the previous ones since θ(T (p θ )/2) = ∆θ(p θ )/2, and since ∆θ(p θ ) < 2π for p θ > 0 (by assumption, ∆θ is decreasing, and equal to 2π on meridians, i.e. when p θ = 0). So extremals are optimal up to t = T (p θ )/2, and the corresponding point is a separating point. Since the metric is complete, each point of the equator is reached by such an extremal and the set of separating points, hence the cut locus, is the equator minus the initial point itself. Consider finally the case when the initial point is neither polar nor equatorial. Then p 2 θ belongs to (0, 1/Γ(ϕ 0 )), and extremals are again optimal up to t = T (p θ )/2. Indeed, there would otherwise exist shorter extremals, which would lead to the existence of shorter extremals for the initial condition ϕ 0 = π/2 too, contradicting the previous fact. The central symmetry s 2 still generates an intersection at t = T (p θ )/2, and ϕ(T (p θ )/2) = π -ϕ 0 so the corresponding separating point belongs to the antipodal parallel of the starting point. Since ∆θ is decreasing, the extremities of the cut are obtained letting p θ tend to ±(Γ(ϕ 0 )) -1/2 (now finite, since ϕ 0 = π/2), and the subarc is closed.

To study the conjugate loci, we start with some properties of the local model at singularity. Setting x := π/2 -ϕ and y := θ, since 1 -sin 2 ϕ ∼ (π/2 -ϕ) 2 when ϕ tends to π/2, a local model for the metric (1) is

ds 2 = dx 2 + dy 2 x 2p ( 13
)
where p is the order of the pole. The equatorial symmetry of Γ translates into (-x) 2p = x 2p , so the discrete symmetry group is preserved. Such almost-Riemannian metrics are related to sub-Riemannian distributions. For p = 1, the local model is the Grušin metric ds 2 = dx 2 + dy 2 /x 2 , which is actually obtained by projecting the Heisenberg sub-Riemannian distribution [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]. This distribution is indeed defined, up to a renormalization, by the following two vector fields on R 3 ,

F 1 (x, y, z) := ∂ ∂x -y ∂ ∂z , F 2 (x, y, z) := ∂ ∂x + x ∂ ∂z ,
and the corresponding sub-Riemannian Hamiltonian is

H(x, y, z, p x , p y , p z ) := 1 2 (p 2 x + p 2 y ) + 2p z (xp y -yp x ) + (x 2 + y 2 )p 2 z ,
which suggests to use cylindrical coordinates. In these variables

H(r, θ, z, p r , p θ , p z ) = 1 2 p 2 r + (p θ /r + rp z ) 2 .
As θ and z are cyclic, the system is integrable in dimension three, and projects onto a Hamiltonian in the (r, z)-space with the desired singularity,

h(r, z, p r , p z ) := 1 2 (p 2 r + r 2 p 2 z ),
when restricting to p θ = 0. For p = 2, the local model is ds 2 = dx 2 + dy 2 /x 4 , which is connected to the flat Martinet sub-Riemannian distribution. Consider indeed the two vector fields on R 3 (see [START_REF] Agrachev | Sub-Riemannian sphere in the Martinet flat case[END_REF])

F 1 (x, y, z) := ∂ ∂x + y 2 ∂ ∂z , F 2 (x, y, z) := ∂ ∂y , so the sub-Riemannian Hamiltonian is H(x, y, z, p x , p y , p z ) := 1 2 p 2 y + (p x + y 2 p z ) 2 .
The two coordinates x and z are cyclic, and the Hamiltonian projects onto h(y, z, p y , p z ) := (1/2)(p 2 y + y 4 p 2 z ) in the (y, z)-space when restricting to p x = 0, providing the higher order singularity. Going back to the general case, we compute geodesics issuing from the origin on the level set {h = 1/2} of the Hamiltonian h(x, y, p x , p y ) := 1 2 (p 2 x + x 2p p 2 y ) so that the initial adjoint state belongs to the union of the two lines, {p x = ±1}. We set λ := p y and restrict to positive λ by symmetry (the trivial geodesics (±t, 0) being obtained for λ = 0). The coordinate x is then

x(t) = 1 p √ λ q(t p √ λ) , ( 14
)
where q is the solution of

q 2 + q 2p = 1, q(0) = 0, q (0) = 1. ( 15
)
Equivalently,

q -1 (u) = u 0 dv p √ 1 -v 2p , u ∈ [-1, 1].
For p = 1, q is harmonic, elliptic for p = 2, hyperelliptic and reciprocal to a hypergeometric function in general. More precisely,

q -1 (u) = 2 F 1 (1/2, 1/(2p); 1 + 1/(2p); u 2p ) • u ( 16
)
where 2 F 1 (a, b; c; z) is the hypergeometric series

2 F 1 (a, b; c; z) = n≥0 α n z n n! , and
α 0 := 1, α n+1 /α n := (n + a)(n + b)/(n + c).
The reciprocal of q is hence equal to Here,

q -1 (u) =
α 0 = 1, α 1 = 1 2 • 1 2p + 1 , α 2 = 3 4 • 2p + 1 8p 2 + 6p + 1 • • •
which gives the usual Taylor series of the reciprocal of the sine function for p = 1, arcsin u = u + u 3 /6 + 3u 5 /40 + 5u 7 /112... Eventually, ẏ = λx 2p , so

y(t) = 1 ( p √ λ) p+1 r(t p √ λ), (17)
where r is defined by a second quadrature,

r(s) := s 0 q 2p . Lemma 9. r = 1 p + 1 (s -qq ).
Proof. Differentiating and using [START_REF] Boscain | Lipschitz classification of almost-Riemannian distances on compact oriented surfaces[END_REF],

1 p + 1 (s -qq ) = 1 p + 1 (1 -q 2 -qq ) = q 2p ,
hence the result since r(0) = 0.

For any p ≥ 1, symmetry reasons imply that non-trivial minimizing geodesics emanating from the origin first intersect on the y-axis (see, e.g., Fig. 1). As a consequence, the cut locus at the origin of the local model is the axis minus the origin itself (compare with the metric on the sphere, Theorem 1). The conjugate locus of the origin is obtained from the set of critical values of the exponential mapping exp (0,0) (t, λ) = (x(t, λ), y(t, λ)), (t, λ) ∈ (R * + ) 2 . Let t 1c (λ) denote the first conjugate time along the geodesic defined by λ > 0. + dy 2 /x 2 (p = 1); sphere, wavefront and conjugate locus of the origin. The wavefront (in blue) is the image of the exponential mapping for a fixed time, t; it is made of endpoints at time t of geodesics. The subset obtained by ruling out the part after the two first self-intersection points of the wavefront on the y-axis is the sphere; it is made of points at distance (minimum time) t of the origin. The conjugate locus (in red) is partly drawn (here, it is the set y = ±C 1 x 2 minus the origin); it is made of critical values of the exponential (it contains so the first singularity of the wavefront portrayed). Compare with the Heisenberg metric in [START_REF] El Alaoui | Small sub-Riemannian balls on R 3[END_REF].

Lemma 10. t 1c (λ) = s p / p √ λ where s p is the first positive root of the envelope equation q = sq .

Proof. Using ( 14) and ( 17), a pair (t, λ) is a critical point of the exponential if and only if t p √ λ is a solution to qr -(p + 1)q r = 0.

Expressing r according to Lemma 9, one gets the result.

Proposition 8. The conjugate locus at the origin of ds 2 = dx 2 + dy 2 /x 2p is the set y = ±C p x p+1 minus the origin where

C p ∼ 2 √ 2 p , p → ∞.
Proof. According to Lemma 10,

x(t 1c (λ)) = 1 p √ λ q(s p ), y(t 1c (λ)) = 1 ( p √ λ) p+1 r(s p ),
so points in the conjugate locus lie on the curve y = ±C p x p+1 where C p = r(s p )/q p+1 (s p ). The whole curve minus (0, 0) is obtained because of symmetries, and conjugate points accumulate towards the origin when |λ| → ∞. It is clear from ( 15) that, as p → ∞, q pointwisely converges towards the Lipschitz function equal to s → s on [0, 1], s → 2 -s on [START_REF] Agrachev | A Gauß-Bonnet formula for contact sub-Riemannian manifolds[END_REF][START_REF] Agrachev | Twodimensional almost-Riemannian structures with tangency points[END_REF]. In particular, the solution s p to the envelope equation is such that s p → 3-; so q (s p ) → -1/3, q 2p (s p ) = 1 -q 2 (s p ) → 8/9, and q p-1 (s p ) → 2 √ 2/3. Now, by virtue of Lemmas 9 and 10,

C p = 1 p + 1 s p -q(s p )q (s p ) q p+1 (s p ) = s p p + 1 q p-1 (s p ) ∼ 2 √ 2 p , p → ∞.
Remark 5. For p = 1, q(s) = sin s, r(s) = (1/2)(s -sin t cos t); t 1c (λ) = s 1 /λ where s 1 is the first positive root of sin s = s cos s.

For p = 2, one obtains

q(s) = √ 2 2 sn dn (s √ 2), r(s) = 1 3 s + √ 2 2 cn sn dn 3 (s √ 2) ,
in terms of Jacobi functions of modulus k = √ 2/2; t 1c (λ) = s 2 / √ λ where τ = s 2 √ 2 (with s 2 < 3) is the first root solution of sn τ dn τ = τ cn τ.

(Compare with the flat Martinet case in [START_REF] Agrachev | Sub-Riemannian sphere in the Martinet flat case[END_REF].)

Before stating the structure result on the conjugate locus, we finally recall the following.

Scholium. Along a Jacobi field tangent to the level set of a Hamiltonian quadratic in the momentum, the Liouville form is constant.

Proof. In coordinates z = (x, p) ∈ R 2n , let H(x, p) = (1/2)(A(x)p|p) with A(x) symmetric; let γ(t) = (z(t), δz(t)) be a Jacobi field, that is ż(t) = -→ H (z(t)) and δ ż(t) = -→ H (z(t))δz(t).
Along γ, the time derivative of the Liouville form p dx is

d dt (p|δx) = ( ṗ|δx) + (p|δ ẋ) = -(∇ x H|δx) + (p|∇ 2 xp H δx) 2(∇xH|δx) + (p|∇ 2 pp H δp) (∇pH|δp)
= (∇ x H|δx) + (∇ p H|δp).

Let now z 0 (σ) be a local parameterization of the level set {H = c}; for a Hamiltonian curve z(t, z 0 (σ)) with initial condition z 0 (σ), H(z(t, z 0 (σ))) = c, so H (z(t, z 0 (σ))) ∂z ∂z 0 (t, z 0 (σ))z 0 (σ)

=:δz(t) = 0.
Along the curve (z(t, σ), δz(t)), the Liouville form is hence constant.

Here, the level set {h = 1/2} of the quadratic Hamiltonian (7) on S 2 is locally parameterized by p θ and the exponential [START_REF] Bonnard | The sphere and the cut locus at a tangency point in two-dimensional almost-Riemannian geometry[END_REF] writes

exp ϕ0 (t, p θ ) = (θ(t, p θ ), ϕ(t, p θ )), so p θ ∂θ ∂p θ (t, p θ ) + p ϕ (t, p θ ) ∂ϕ ∂p θ (t, p θ ) = 0 (18)
since ∂θ(0, p θ )/∂p θ = ∂ϕ(0, p θ )/∂p θ = 0.

Lemma 11. Critical points (t, p θ ) of the exponential are characterized either by ∂θ(t, p θ )/∂p θ = 0 when φ = 0, or by ∂ϕ(t, p θ )/∂p θ = 0 when p θ = 0.

Proof. The point (t, p θ ) is critical if and only if

θ(t, p θ ) ∂ϕ ∂p θ (t, p θ ) -φ(t, p θ ) ∂θ ∂p θ (t, p θ ) = 0.
When p ϕ (t, p θ ) = φ(t, p θ ) = 0, one can multiply both sides by p ϕ and use [START_REF] Caillau | On some Riemannian aspects of two and three-body controlled problems[END_REF] to get ∂θ ∂p θ (t, p θ ) p θ θ(t, p θ ) + p ϕ (t, p θ ) φ(t, p θ ) =h=1/2 = 0 whence the result. Similar computation when p θ = 0.

Theorem 2. Under the assumption that ∆θ is strictly decreasing and convex for p θ > 0, conjugate loci of the metric (1) are reduced to the opposite pole for poles, have four cusps otherwise.

Proof. Let ϕ 0 = π/2. Consider an extremal defined by a positive p θ and p ϕ 0 = +1. For t in (T (p θ )/4, 3T (p θ )/4), φ = 0 and the extremal can be parameterized by ϕ according to

θ(ϕ, p θ ) = ∆θ(p θ ) 2 + π/2 ϕ f (φ, p θ ) dφ,
where, as before,

f (ϕ, p θ ) = dθ dϕ = Γ(ϕ)p θ 1 -Γ(ϕ)p 2 θ •
The conjugacy condition is ∂θ/∂p θ = 0 (Lemma 11) so the coordinate ϕ 1c (p θ ) of the first conjugate point is solution of

π/2 ϕ1c(p θ ) ∂f ∂p θ (ϕ, p θ ) dϕ = - ∆θ (p θ ) 2 > 0,
in order that ϕ 1c (p θ ) < π/2 (since ∂f /∂p θ > 0). By differentiating the previous equality, one gets

ϕ 1c (p θ ) = ∂f ∂p θ (ϕ 1c (p θ ), p θ ) -1 ∆θ 2 (p θ ) + π/2 ϕ1c(p θ ) ∂ 2 f ∂p 2 θ (ϕ, p θ ) dϕ ,
which is positive first because ϕ 1c (p θ ) < π/2, then because

∂ 2 f ∂p 2 θ = 3Γ 2 (ϕ)p θ (1 -Γ(ϕ)p 2 θ ) 5/2 > 0,
and by virtue of the nonnegativeness of ∆θ . In particular, the parameterization p θ → (θ(ϕ 1c (p θ ), p θ ), ϕ 1c (p θ )) of the conjugate locus is regular, and the set has no cusp for positive p θ . The tangent vector along the conjugate locus is

∂θ ∂ϕ ϕ 1c + ∂θ ∂p θ 0 ∂ ∂θ + ϕ 1c ∂ ∂ϕ ,
which is proportional to f (ϕ 1c (p θ ), p θ )∂/∂θ + ∂/∂ϕ. On the one hand, as in the regular case discussed in [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF], f (ϕ 1c (p θ ), p θ ) tends to 0 when p θ → 0+, and the locus has a first meridional cusp because of the axial symmetry s 1 .

On the other hand, when p θ goes to +∞, the analysis on the local model (Proposition 8) shows that, contrary to the regular case, the conjugate locus has two tangential contacts with the meridian that, combined again with symmetry s 1 , form a second meridional cusp at the initial point where conjugate points accumulate. The central symmetry s 2 gives the symmetric part, whence the result. When ϕ 0 = π/2, the same reasoning shows that there is no cusp for p θ ∈ (0, 1/ Γ(ϕ 0 )); there are two meridional cusps due again to symmetry s 1 , and two cusps tangent to the antipodal parallel containing the cut locus (see Theorem 3.6 of [START_REF] Bonnard | Conjugate and cut loci of a two-sphere of revolution with application to optimal control[END_REF]).

5 The (1, 2, 1) case

The results of Section 4 are applied to the second example discussed in §3 (the application to the first example, the Grušin metric on the sphere, is straightforward; see [START_REF] Caillau | On some Riemannian aspects of two and three-body controlled problems[END_REF]). Consider the metric (1) with p = 2 and (a 0 , a 1 , a 2 ) = (1/4)(1, 2, 1), that is

R(X) = 1 -X/2 1 -X 2 = 1 4 1 + 2 1 -X + 1 (1 -X) 2 • Lemma 12.
The coordinate ϕ is parameterized by the family of elliptic curves

C p θ : Y 2 = 4(1 -X)[X(2 -X) 2 -4p 2 θ (1 -X) 2 ], p θ ∈ R. (19)
Proof. Set Y := (2 -X) Ẋ on the level set 1/2 of the quadratic Hamiltonian [START_REF] Besson | Géodésiques des surfaces de révolution[END_REF].

In the sequel, ℘ denotes the Weierstraß function with invariants g 2 , g 3 (see [START_REF] Jones | Complex Functions. An Algebraic and Geometric Viewpoint[END_REF]).

Proposition 9.

sin 2 ϕ = ℘(z) -4/3 ℘(z) -1/3 , z ∈ R, t = z + 1 ℘ (a) 2ζ(a)z + ln σ(z -a) σ(z + a) z 0
with a such that ℘(a) = 1/3 and invariants

g 2 = 16 3 + 16p 2 θ , g 3 = 64 27 - 16 3 p 2 θ . (20)
Proof. The rational transform

u = 1 3 + 1 1 -X , v = u 2 Y,
sends to infinity the fixed root X = 1 in the right hand side of [START_REF] Faridi | Geodesics and deformed spheres[END_REF], and allows to recast the equation of C p θ under the canonical form v 2 = 4u 3 -g 2 u -g 3 with invariants [START_REF] Figalli | Nearly round spheres look convex[END_REF]. When parameterizing the elliptic curve C p θ by the Weierstraß function, (u, v) = (℘(z), ℘ (z)), only the unbounded component of the real cubic has to be used since X = sin 2 ϕ ∈ (0, 1] (that is ℘(z) > 4/3), so z ∈ R. In this parameterization, the change of time from z to t verifies For a geodesic originating from the singularity, ϕ = π/2 at t = 0 which corresponds to z = 0.

dt dz = 1 + 1 ℘(z) -1/3 > 0. (21)
As a function of z, the coordinate ϕ is a doubly periodic meromorphic function. Its lattice of periods 2ωZ + 2ω Z depends on p θ (the Weierstraß half-periods ω, ω are functions of p θ ) and is real rectangular: ω ∈ R, ω ∈ iR, and the periods can be choosen so that their ratio τ := ω /ω belongs to the Poincaré upper half-plane H = {x + iy ∈ C, y > 0}. Lattices are classified up to conformal transformations; these transformations are Möbius transforms in the Fuschian modular subgroup PSL(2, Z) = SL(2, Z)/ ± id of automorphisms of H, so the moduli space of congruences of lattices is H/PSL(2, Z). The modular function j establishes a one-to-one correspondence between these moduli and the complex plane; in terms of the invariants of the elliptic curve,

j(τ ) = g 3 2 ∆
where ∆ = g 3 2 -27g 2 3 is the discriminant of the elliptic curve with ratio of periods τ . For the family of elliptic curves C p θ (see [START_REF] Figalli | Nearly round spheres look convex[END_REF]),

j(τ (p θ )) = 16(1 + 3p 2 θ ) 3 27p 2 θ (8 + 13p 2 θ + 16p 4 θ ) • (22)
Corollary 1. The number of conformal classes associated with C p θ is -equal to 1 for p θ ∈ (0, 1/2) ∪ {2/3}, -equal to 2 for p θ ∈ {1/2, √ 2},

-equal to 3 for p θ ∈ (1/2, 2/3) ∪ (2/3, √ 2) ∪ ( √ 2, ∞).
The only square lattice is obtained for p θ = 2/3.

Proof. The rational fraction [START_REF] Grušin | A certain class of elliptic pseudodifferential operators that are degenerate on a submanifold (Russian)[END_REF] has exactly one global minimum at p θ = 2/3 and one local maximum at p θ = √ 2. Moreover,

j(τ (p θ )) = 28 27 - 16(p 2 θ -1/4)(p 2 θ -2) 2 27p 2 θ (8 + 13p 2 θ + 16p 4 θ )
showing that the value of the local maximum is also attained when p θ = 1/2. = ln e -4ηa = -4ηa mod 2iπ.

The term in brackets in [START_REF] Halphen | Traité des fonctions elliptiques et de leurs applications[END_REF] has to be positive because of [START_REF] El Alaoui | Small sub-Riemannian balls on R 3[END_REF], so ℘ (a) = +2i by the previous lemma (℘(a) = 1/3) and

T (p θ ) 4 = ω + ωζ(a) -ηa i mod π 2 •
The bounded component of the imaginary cubic is parameterized by z ∈ ω + iR so a = ω + ib and ζ(a) = η + iβ for some real b and β.

In order to verify the assumptions of Theorems 1 and 2, one has to differentiate the quasi-period ∆θ with respect to p θ . By means of Proposition 6, it suffices to compute derivatives of the period T . Define

∂ := 1 δ d dp θ , δ := ∆ 256p θ = p θ (8 + 13p 2 θ + 16p 4 θ ). Lemma 14. ∂ω = -Aω -Bη, ∂b = -Ab -Bβ + D, ∂η = Cω + Aη, ∂β = Cb + Aβ + E, with A, B, C, D, E in R[p θ ], A = 1 3 (4 + 13p 2 θ + 24p 4 θ ), B = -2 + p 2 θ , C = 4 9 (-2 -5p 2 θ + 3p 4 θ ), D = 2 + 6p 2 θ , E = 1 3 (-4 + 9p 2 θ ).
Proof. The derivatives of the (half) period and quasi-period ω and η = ζ(ω) of ℘ and ζ, respectively, with respect to the invariants g 2 , g 3 are known (see [23, p. 307]), and

∂ = 32 3 3 ∂ ∂g 2 - ∂ ∂g 3
according to [START_REF] Figalli | Nearly round spheres look convex[END_REF]. Moreover, Proof. Applying Lemma 14 rules to T , which is bilinear in (ω, η, b, β), one observes that the coefficients of b and β cancel.

℘(a) = ℘(a(p θ ), p θ ) = 1 3 implies ℘ (a)∂a + ∂℘(a) = 0, so ∂a = - 1 2i ∂℘ (a
Proposition 11.

∆θ (p

θ ) = - 4 8 + 13p 2 θ + 16p 4 θ 2 3 (11 + 12p 2 θ )ω + 7η , ∆θ (p θ ) = 4 p θ (8 + 13p 2 θ + 16p 4 θ ) 2 × 2 3 (24 + 181p 2 θ + 830p 4 θ + 480p 6 θ )ω + (-24 + 143p 2 θ + 400p 4 θ )η .
Proof. Lemma 15 combined with Proposition 6 gives the first order derivative of the quasi-period; the second order one is obtained by applying the previous rules anew.

The penultimate result of the section ensures that the cut and conjugate loci in the (1, 2, 1) case have the structure described in Theorem 1 and Theorem 2.

Theorem 3. ∆θ < 0 ≤ ∆θ in the (1, 2, 1) case.

Proof. The first derivative ∆θ , expressed as linear combination of the nonnegative (quasi-)periods ω and η, is obviously negative. To obtain nonnegativity of ∆θ , one has to check the sign of 6 Deformation of the round sphere

We consider the following deformation of the metric (1):

XR(λX)dθ 2 + dϕ 2 (X = sin 2 ϕ), λ ∈ [0, 1]. (24)
For λ = 0, one gets the canonical metric on the sphere as R(0) is normalized to one, so we have an homotopy connecting the round metric to the singular one. For λ < 1, the metric is Riemannian (no equatorial singularity). In the case R(X) = 1/(1 -X), such a deformation of the round metric appears in [START_REF] Bonnard | Geodesic flow of the averaged controlled Kepler equation[END_REF][START_REF] Faridi | Geodesics and deformed spheres[END_REF].

See also [START_REF] Figalli | Nearly round spheres look convex[END_REF] for perturbation results of the round sphere.

Proposition 13. The curvature of the metric ( 24) is a rational fraction homogeneous in R; for λ ∈ [0, 1),

K λ (X) = 1 + λ(4X -3) R (λX) R(λX) + λ 2 X(1 -X) R 2 (λX) R 2 (λX) - 2R (λX) R(λX) .
Proof. The result follows from the fact that, given a metric Γ(ϕ)dθ 2 + dϕ 2 on S 2 , the Gaussian curvature is K

= -( √ Γ) / √ Γ.
Corollary 2 (Concentration of curvature). The following curvature estimates hold for the metric (24):

K λ (X) = 1 + λR (0)(4X -3) + O(λ 2 ) when λ → 0, K 1 (X) ∼ - p(p + 1) 1 -X when X → 1, K λ (1) ∼ p 1 -λ when λ → 1,
where p is the order of the pole of R.

Corollary 3. Cut loci of the metric ( 24) are closed antipodal subarcs for λ close enough to zero.

Proof. Since K λ (X) ∼ 4λ(R /R)(0)+O(λ 2 ) with (R /R)(0) > 0 by virtue of (2), the curvature is monotone non-decreasing along half-meridians and the result follows from [START_REF] Sinclair | The cut locus of a two-sphere of revolution and Toponogov's comparison theorem[END_REF] main theorem.

The existence of conjugate points for metrics with singularities is typical of the almost-Riemannian setting where conjugate points may exist although curvature remains negative whenever defined (see, e.g., [START_REF] Agrachev | A Gauß-Bonnet like formula on two-dimensional almost-Riemannian manifolds[END_REF]). In the (1, 2, 1) case of §5 for instance,

K 1 (X) = - (1 + X)(4 -X) (2 -X)(1 -X) < 0 (X = sin 2 ϕ).
Corollary 2 establishes that although curvature is negative in the neighbourhood of the singularity and even tends to -∞ when ϕ → π/2, there is a concentration of positive curvature on the singularity itself as K λ (1) → +∞ when λ → 1-, responsible for the existence of conjugate points after crossing the singularity (cut, and thus conjugate points being located after the equator by antipodality). An alternative interpretation of the singularity comes from the following fact.

Proposition 14. For R(X) = 1/(1 -X) and λ ∈ [0, 1), the metric ( 24) is conformal to the canonical metric on a oblate ellipsoid of revolution of unit semi-major axis and √ 1 -λ semi-minor axis.

Proof. A parameterization of such an ellipsoid of semi-minor axis µ is x = sin ϕ cos θ, y = sin ϕ sin θ, z = µ cos ϕ.

In these coordinates, the restriction of the flat R 3 metric reads

sin 2 ϕ dθ 2 + (1 -(1 -µ 2 ) sin 2 ϕ)dϕ 2 = (1 -λX)[XR(λX)dθ 2 + dϕ 2 ] with λ = 1 -µ 2 .
When λ tends to 1, µ = √ 1 -λ tends to zero so the oblate ellipsoid collapses onto a two-sided Poincaré disk, each face being endowed with the flat metric.

Corollary 4. For R(X) = 1/(1 -X), the metric (1) is conformal to the following constant curvature metrics on the Poincaré disk: (i) the flat metric dρ 2 + ρ 2 dθ 2 (K = 0), (ii) the canonical Poincaré metric (dρ 2 + ρ 2 dθ 2 )/(1 -ρ 2 ) 2 (K = -1), where (ρ, θ) are polar coordinates on D.

Proof. Setting ρ = sin ϕ, one retrieves the standard polar coordinates on the disk. In these coordinates, dρ 2 = (1 -X)dϕ 2 , so

X 1 -X dθ 2 + dϕ 2 = dρ 2 + ρ 2 dθ 2 1 -ρ 2 •
In this case, crossing the equatorial singularity can thus be interpretated as crossing the boundary of the disk to go from one side of D to the other. This can also be seen, like in the flat case, as generating reflections of geodesics with the boundary. As for the canonical Poincaré metric, those reflections turn out to be orthogonal in general.

Proposition 15. For the metric (1), crossing the equatorial singularity is interpretated on the Poincaré disk as reflecting on the boundary. Reflections on the boundary of the metric [START_REF] Jean | Sub-Riemannian Geometry[END_REF] are tangential (except for meridians) when λ < 1, and orthogonal when λ = 1.

Proof. In polar coordinates on D, the Hamiltonian (7) reads on {h = 1/2}. Reparameterizing time according to dτ = dt √ 1 -X we get dθ dτ ∼ p θ a p (1 -X) p-1/2 → 0 and dρ dτ ∼ 1, X → 1, so contacts with ∂D (i.e. for ρ 2 = X = 1) are orthogonal reflections. Besides, using homotopy [START_REF] Jean | Sub-Riemannian Geometry[END_REF] to replace the singular metric by a Riemannian one when λ < 1 changes the contact for p θ = 0 (meridians, obtained for p θ = 0, obviously remain perpendicular to the boundary). Indeed, the deformed Hamiltonian is

h λ (ρ, θ, p ρ , p θ ) := 1 2 (1 -X)p 2 ρ + p 2 θ XR(λX) ,
and, when X → 1, θ ∼ p θ R(λ) and ρ = p ϕ √ 1 -X since p ϕ = p ρ √ 1 -X, which remains finite. So θ = 0 if p θ = 0 while ρ = 0 at X = 1, and contacts with the boundary are tangential outside meridians. Figure 3: (1, 2, 1) case; geodesics, cut and conjugate loci on the disk. On the left, metric obtained by homotopy from the round one for λ < 1; the initial condition is on the boundary, tangential contacts of the geodesics (in blue) with ∂D are observed. On the right, singular metric (λ = 1) for the same initial condition (compare with Fig. 2); contacts with ∂D are orthogonal. In both cases, the conjugate locus (red dots) is the catacaustic generated by reflections of geodesics on the boundary [START_REF] Bruce | [END_REF]. (In the case of the flat metric on the disk, the catacaustic of geodesics-that is of straight lines-originating from a point on the boundary, i.e. at the singularity, would be a cardioid.) Since reflections are orthogonal (or specular ), the figure can also be interpretated as a billiard on the disk endowed with a particular Riemannian metric.

A Asymptotics in the (1, 2, 1) case Let α, β, γ and δ = 1 be the roots of the degree four polynomial P (X, p θ ) involved in the computation (compare with (19)), P (X, p θ ) := (1 -X)(X(2 -X) 2 -4p 2 θ (1 -X) 2 ).

An alternative quadrature for the period is

T (p θ ) = 4 √ A 2 B 1 Π(ν, k) + 2 -p p -q K(k) (25)
where K and Π are respectively complete elliptic integrals of first and third kind, Π(ν, k) :=

1 0 dv (1 -νv 2 ) √ 1 -v 2 √ 1 -k 2 v 2 , K (

Proposition 7 .Lemma 8 .

 78 Cut points of the metric (1) are either conjugate or separating points. Both T and ∆θ vanish when |p θ | → ∞.

  α n = (a) n (b) n /(c) n for a = 1/2, b = 1/(2p) and c = 1 + 1/(2p), the notation (a) n standing for the Pochhammer symbol (a) n := a(a + 1) • • • (a + n -1).

Figure 1 :

 1 Figure1: Grušin metric, dx 2 + dy 2 /x 2 (p = 1); sphere, wavefront and conjugate locus of the origin. The wavefront (in blue) is the image of the exponential mapping for a fixed time, t; it is made of endpoints at time t of geodesics. The subset obtained by ruling out the part after the two first self-intersection points of the wavefront on the y-axis is the sphere; it is made of points at distance (minimum time) t of the origin. The conjugate locus (in red) is partly drawn (here, it is the set y = ±C 1 x 2 minus the origin); it is made of critical values of the exponential (it contains so the first singularity of the wavefront portrayed). Compare with the Heisenberg metric in[START_REF] El Alaoui | Small sub-Riemannian balls on R 3[END_REF].

Introducing

  Weierstraß functions ζ and σ (℘ = -ζ and ζ = σ /σ), one has (see [23]) ℘ (a) dz ℘(z) -℘(a) = 2ζ(a)z + ln σ(z -a) σ(z + a) •

2 .

 2 Proposition 10. T (p θ ) 4 = (1 + β)ω -bη mod π 2 where η := ζ(ω), b := Im a and β := Im ζ(a). Lemma 13. For all p θ > 0, (1/3, ±2i) belongs to the bounded component of the imaginary cubic C p θ .Proof.4(1/3) 3 -(1/3)g 2 -g 3Proof of Proposition 10. As X = sin 2 ϕ, the period of X is half the period of ϕ so T

)

  (because ℘ (a) = 2i), and one also knows the derivatives of ℘ (and ζ) with respect to g 2 , g 3 (see[23, p. 298]). Then ∂b = Im ∂a. Similarly, ∂[ζ(a(p θ ), p θ )] = -℘(a)∂a + ∂ζ(a) = -1 3 ∂a + ∂ζ(a) and ∂β is obtained taking the imaginary part. Lemma 15. ∂T is R[p θ ]-linear in (ω, η), ∂T (p θ ) 4 = -(A + E)ω -(B + D)η.

Figure 2 :

 2 Figure 2: (1, 2, 1) case; geodesics, cut and conjugate loci on the sphere. Symmetric geodesics (in blue) from opposite hemispheres intersect on the equator that, deprived of the origin, forms the cut locus (in black; Theorem 1). The envelope of geodesics generate the conjugate locus (red dots) which has for meridional cusps (Theorem 2).

= sin 2 ϕ = ρ 2 .( 1 -

 21 X) p and | ρ| = (1 -X)|p ρ |

  k) := Π(0, k), and∆ := 4(β -α)(β -δ)(γ -α)(γ -δ), σ := (α + δ)(β + γ) -2(αδ + βγ), l 1 := σ -√ ∆ (β -γ) 2 , l 2 := σ + √ ∆ (β -γ) 2 ,

  2

Supported by ANR Geometric Control Methods (project no. NT09 504490). Also supported by Conseil Régional de Bourgogne (contract no. 2009-160E-160-CE-160T) and by SADCO Initial Training Network (FP7 grant no. 264735-SADCO).

With the same notation as before,

where the elliptic integral of third kind is now incomplete, and where

Similarly,

In order to compute the asymptotics of Proposition 12, we need expansions of these roots in the neighbourhood of p θ = 0 and p θ = ∞. Such expansions, which are propagated to T using [START_REF] Jones | Complex Functions. An Algebraic and Geometric Viewpoint[END_REF], are available in √ ε-scale as both


Thomas Moulard
email: thomas.moulard@gmail.com

Florent Lamiraux

Karim Bouyarmane

Eiichi Yoshida

RobOptim: an Optimization Framework for Robotics

Keywords: Numerical Optimization, Software, Humanoid Robotics

Numerical optimization is useful for various areas of robotics. However tackling optimization problems properly requires the use of non-trivial algorithms whose tuning is challenging. RobOptim aims at providing a unified framework for different categories of optimization problems while relying on strong C++ typing to ensure efficient and correct computations. This paper presents this software, demonstrates its genericity and illustrates current use by two full scale robotics examples.

I. Introduction

Over the past years, numerical optimization proved itself particulaly suited for various robotics applications such as posture or trajectory optimization [START_REF] Miossec | Development of a software for motion optimization of robots -application to the kick motion of the hrp-2 robot[END_REF], [START_REF] Saab | Dynamic whole-body motion generation under rigid contacts and other unilateral constraints[END_REF], robot control [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF] and more. These applications yield both linear and non-linear optimization problems with equalities and inequalities constraints. Robot control also relies on other types of optimization such as quadratic programming. As robot control algorithms run in real-time, it leads to strong constraints on the implementation efficiency. The design and implementation of a solver is tedious and error-prone. Avoiding numerical precision issues, ensuring that the algorithm behaves properly in all cases even with illconditioned problems is challenging, in particular for roboticists who are not necesseraly experts in optimization techniques. Among available optimization toolboxes, the Matlab Optimization Toolbox [START_REF] Mathworks | Matlab Optimization Toolbox User's Guide[END_REF], the Open Optimization library, OPT++, IPOPT [START_REF] Biegler | Large-scale nonlinear programming using ipopt: An integrating framework for enterprise-wide dynamic optimization[END_REF], SciPy [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF] and the GSL (Gnu Scientific Library) [START_REF] Galassi | Gnu Scientific Library: Reference Manual[END_REF] all provides several optimization algorithms. Unfortunately, these libraries suffer from several drawbacks: some are difficult to use, some do no provide support for advanced algorithms such as support for constrained optimization or have efficiency issues, etc. They also all lack a unified model expressing optimization problems.

RobOptim solves these limitations by introducing a framework to model any continuous optimization problem, constrained or not. The design is focused on providing a easy-touse C++ set of safe and efficient libraries which can be used to prototype robotics applications. Let us emphasize here that our framework does not implement any numerical optimization algorithm. Instead, it provides a model for functions, problems and solvers on the one hand and plug-ins wrapping existing solvers as CFSQP or IPOPT on the other hand. The RobOptim computational model will be first introduced in section II and different applications will be detailed in section III. In particular, an extension of RobOptim for a particular category of problem has been implemented recently, demonstrating the ability of RobOptim to support a large variety of problems. The conclusion will detail advantages and limitations of the current approach and summarize the roadmap for the next developments. The RobOptim framework is divided into three parts: the solver, core and application layers Figure 2. The core layer provides a computational model expressive enough to formulate different types of optimization problems. The solver layer gathers various optimization algorithms. The application layer consists in several application-dependent packages bundling reusable mathematical functions.

A. Mathematical function representation

Continuous optimization problems can be defined as follows: where 𝑓 : R 𝑛 → R is the cost function and X ⊂ R 𝑛 is the subspace of the admissible solutions. This space is usually defined by a set of inequality and equality constraints:

min x∈R 𝑛 𝑓 (x) under the constraint x ∈ X (1)
X { 𝑐𝑖(𝑥) = 0 𝑖 ∈ 𝜉 𝑐𝑗(𝑥) ≤ 0 𝑗 ∈ 𝜈 (2)
𝑐𝑖 : R 𝑛 → R, 𝑐𝑗 : R 𝑛 → R are respectively the set of equality and inequality constraints. 𝑖 and 𝑗 are the indices identifying the constraints.

The mathematical functions 𝑓 and 𝑐 𝑘 (𝑘 ∈ 𝜉∪𝜈) must provide a way to evaluate their result at any point they are defined. Their associated gradient and hessian can also be provided. Finally, each of these function may be categorized as belonging to a set of functions matching a particular structure which may help the resolution. One example is linear functions. The goal of the RobOptim core layer is to express these features through the C++ typing rules.

Depending on the solving algorithm, it may be necessary to obtain the Jacobian and in some cases even the functions Hessian. Hence, functions providing maximum information about themselves will be compatible with a larger proportion of solvers.

In RobOptim, each mathematical function is represented by a different C++ class. All of these types then inherit the abstract kind of mathematical functions they represent. The following kind of mathematical functions are bundled with RobOptim core: function that can evaluate itself (Function type), function that can evaluate itself and its gradient (DifferentiableFunction), function that can evaluate itself, its gradient and its hessian (TwiceDifferentiableFunction), etc. Let 'A <: B' be the relationship "type B is a subtype (subclass) of A" or "B inherits from A". Then a partial order can be defined where: Function <: DifferentiableFunction <: TwiceDifferentiableFunction.

B. Optimization problem definition and problem resolution

Once the cost function and all constraint functions have been implemented, an optimization problem has to be built. RobOptim core provides a meta-class Problem parametrized by two parameters 𝐹 and 𝐶𝐿. 𝐹 the cost function type and 𝐶𝐿 is the list of the constraints types. A non-linear problem with constraints then has the following type:

Problem<DerivableFunction, vector<LinearFunction, DerivableFunction > >

The constraints can be either linear or non-linear. With this type definition, the constraints will be divided into two categories which will help the solver perform efficiently.

Additionally, bounds can be set on the optimization variables and a starting point for the optimization process may be specified. When the constraints are added to the problem, each one of them is associated with a validity interval. If this interval is reduced to a point, the constraint is an equality constraint. At both compile time and run-time, RobOptim checks that only valid problems are built. For instance if one adds a linear constraint then RobOptim checks at compile time that the constraint is a subtype of the LinearFunction type.

A solver that will solve the problem needs to be instantiated. Each solver is parametrized by the same types as optimization problems. Therefore the solver 𝑆 < 𝑃1, 𝐶𝐿1 > can solve the problem 𝑃 < 𝑃2, 𝐶𝐿2 > if the following relation is true:

𝑃1 <: 𝑃2 ∧ ∀𝑖, 𝐶𝐿1(𝑖) <: 𝐶𝐿2(𝑖) (3)
Basically, the problem can be solved if all types provide enough or more information than necessary. For instance if gradients are required, the function may also provide hessian computation but if gradient is lacking the compile time assertions will fail and prevent the user from building an invalid optimization problem.

By separating problem expression from solver, dynamic changes of the solving algorithm are possible. Each solver is bundled as a plug-in which is loaded at run-time. The interest is to let the user freely change the problem complexity during the design process. Other frameworks would require a different API depending on the kind of optimization at hand, where through RobOptim changes are minimal. One may choose to use a more powerful solver than required at first and then refine the choice or implement later a new plug-in providing the best algorithm for one particular application. These features are provided through meta-programming techniques and come with a near zero cost1 at runtime and are unique to RobOptim.

C. Cost and constraint toolbox

Unlike others frameworks where computations are tightly linked to one problem and one solver, RobOptim abstraction mechanism allows user to implement toolboxes of reusable functions. This part of RobOptim is dedicated to robotics. The "trajectory" toolbox from the RobOptim application layer is currently providing trajectores defined as B-Splines and associated mathematical functions to realize minimal time optimization for instance.

D. Toward easy and safe problem design

RobOptim is heavily relying on templated metaprogramming, a C++ language feature allowing in particular to define parametrized types and execute algorithms at compile time [START_REF]ISO 14882 Programming Languages -C++[END_REF]. By defining problems and solvers through parametrized types, the compiler is able to check that the functions used to instantiate the optimization problem are compatible with the type of problem under construction. Thus, only solvers supporting this kind of problem can be used. Through meta-programming, these safety checks are evaluated at compile-time and thus do not impact final performances. Regarding ease of use, having a unified representation of all models which matches closely the underlying theory simplifies the understanding of the implementation process. Finally, RobOptim relies on modern tools such as Boost and Eigen to support natural ways of implementing algebra computation. Additional tools such as gradient checks through finite differentiation is also provided to help ensuring functions correctness.

III. Applications and case study

RobOptim has been used to solve several different types of robotics problems. The two scenarii that will be detailed here are footsteps optimization and posture optimization for a humanoid robot. Another important point is the extensibility of the framework. To demonstrate RobOptim capacities to adapt to new types of problems an example of such extension will be given. These experiments aim at generating motion for the humanoid robot HRP-2 [START_REF] Kaneko | Humanoid robot hrp-2[END_REF].

A. Step planning for humanoid robots

Generating a walking motion in an environment cluttered with obstacles is challenging. One commonly used approach is the Rapidly-exploring Random Trees (RRT) method applied to the robot bounding box or to the robot itself sliding on the ground [START_REF] Dalibard | Small-space controllability of a walking humanoid robot[END_REF]. This probabilistic algorithm will try to create a path between the starting point and the goal point by sampling configurations randomly and is able to find solutions for highly dimensioned problems on a reasonable time. However, the probabilistic nature of these algorithms leads to paths which seem unnatural. One solution to improve these paths is numerical optimization. In this application, RobOptim has been used to optimize a biped robot walking trajectory determined beforehand by a motion planning algorithm.

Let 𝛾 be the initial robot waist trajectory defined as a B-Spline from 𝑡𝑚𝑖𝑛 to 𝑡𝑚𝑎𝑥. A free time trajectory Γ is defined from 𝑇𝑚𝑖𝑛 to 𝑇𝑚𝑎𝑥 = 𝑇𝑚𝑖𝑛 + 𝜆(𝑡𝑚𝑎𝑥 -𝑡𝑚𝑖𝑛) as

Γ 𝜆 (𝑇 ) = 𝛾(𝑡𝑚𝑖𝑛 + 1 𝜆 (𝑇 -𝑇𝑚𝑖𝑛)) (4) Time control point 1 . . . control point 𝑛 𝜆 𝑥 0 𝑦 0 𝜃 0 . . .

𝑥𝑛 𝑦𝑛 𝜃𝑛

Table I Optimization variables for walking motion problem

The use of a free time trajectory allows the solver to optimize both trajectory shape and speed by respectively changing the control points and the scale 𝜆. Optimization variables are illustrated by Table I. The free time trajectory derived from the trajectory 𝛾 is used as the state of the solver.

Cost function 𝑓 (x) = 𝜆

Speed constraint ∀𝑇, foot, (

𝑣 𝑥 foot 𝑣 𝑥 max ) 2 + ( 𝑣 𝑦 foot 𝑣 𝑦 max ) 2 -1 ≤ 0 Distance constraint ∀𝑇, distance(obstacle, Γ(𝑇 )) > 0 Table II Walking optimization problem formulation
The problem cost function is defined as 𝜆. It corresponds to minimizing the time by accelerating the trajectory as much as possible. To preserve a feasible final trajectory while encouraging forward motion, a speed constraint is added which takes separately into account the front speed 𝑣 𝑥 and the lateral speed 𝑣 𝑦 of each foot of the robot. Another constraint is preventing the robot from colliding with obstacles. The problem definition is detailed in Table II. This non-linear problem has been successfully solved by RobOptim while relying on the proprietary solver CFSQP [START_REF] Lawrence | A C code for solving (large scale) constrained nonlinear (minimax) optimization problems, generating iterates satisfying all inequality constraints[END_REF] internally and result is shown on Figure 3.

B. Posture optimization for humanoid robots

Another problem is posture optimization, which consists in choosing the best configuration such that a system will accomplish the objectives it has been assigned (see 1). In this example, the goal is to find a configuration as close as possible to a goal posture while taking into account various constraints. In this case, robots and environment objects are considered as elements of the optimization problem. Constraints include: robot and objects static equilibrium, Newton's third law, Coulomb friction model, fixed grasp model for bilateral contact forces, joint limits, torque limits and FEM-modeled deformation of the environment [START_REF] Bouyarmane | FEM-based static posture planning for a humanoid robot on deformable contact support[END_REF]. This problem is described extensively in [START_REF]Multi-contact stances planning for multiple agents[END_REF], [START_REF] Bouyarmane | Static multi-contact inverse problem for multiple humanoid robots and manipulated objects[END_REF] and demonstrated the ability of RobOptim to solve complex large scale non-linear problems.

C. Extending the framework: least square optimization

Initially, non-linear optimization problems with constraints have been solved with RobOptim. By providing a generic computational framework, it is possible to extend RobOptim to support other types of problems such as least square optimization. A least square optimization problem is an unconstrained problem, the cost function of which is defined by

𝑓 (x) = 𝑚 ∑ 𝑖=1 𝑓𝑖(x) 2 (5)
where 𝑓1, • • • , 𝑓𝑚 are 𝑚 non-linear differentiable functions from R 𝑛 to R.

To build functions matching this definition, another function type called SumOfC1Squares has been derived from the DifferentiableFunction class. The constructor of this new type takes as input the function from R 𝑛 to R 𝑚 the coordinates of which are 𝑓1, • • • , 𝑓𝑚. SumOfC1Squares then redefines evaluation and gradient of ∑ 𝑚 𝑖=1 𝑓𝑖(x) 2 . The CMinPack solver has then been wrapped into a RobOptim solver class which takes as input problem of the type: Problem<SumOfC1Squares, vector<> >.

The introduction of this new type of solver did not require any change in the RobOptim infrastructure and proved the genericity of the approach.

IV. Conclusion

A new approach to optimization problem representation has been exposed in this paper. By expressing numerical optimization problems through C++ typing, the RobOptim optimization framework provides a unified computational model. Moreover, advanced C++ template meta-programming techniques preserve efficiency while giving the capacity to express high-level mathematical objects such as cost functions and constraints. These features proved useful to solve different types of robotics problems.

However though being generic, the RobOptim core layer suffers from limitations. It still lacks support for optimization in non-scalar spaces such as SO [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF]. By importing knowledge about the structure of the optimization variables, solvers could realize more efficient computations. For instance, 3D rotations can be represented by homogeneous matrices, quaternions, vectors and angles, etc. Switching from one representation to another may lead to better convergence and a decrease in the number of necessary mathematical operations. One goal would be to both be able to express information regarding the optimization variables structure and to find a method to help solvers rely on these additional information. To achieve this level of expressiveness, modern C++ features are of great help and no optimization framework is providing yet such features when it comes to design an optimization problem.

Therefore, RobOptim is a step forward toward making optimization techniques available for non-expert by providing easyto-understand and generic model, strong C++ object model ensuring safety and toolboxes bundling robotics oriented functions that may be reused when one builds its own optimization problem.

Figure 1 .

 1 Figure 1. Posture optimization while considering dynamic constraints with RobOptim.

Figure 3 .

 3 Figure 3. Walking trajectory optimization

RobOptim core does not realize copy so the additional runtime cost is only due to calls to virtual functions.