Bernard Bonnard

Monique Chyba

Jérémy Rouot

Geometric and Numerical Optimal Control -Application to Swimming at Low Reynolds Number and Magnetic Resonance Imaging

Keywords: spins with B 1 -inhomogeneity

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Preface

The motivation for the notes presented in this volume of BCAM SpringerBriefs comes from a multidisciplinary graduate course offered to students in Mathematics, Physics or Control Engineering (at the University of Burgundy, France and at the Institute of Mathematics for Industry Fukuoka, Japan). The content is based on two real applications, which are the subject of current academic research programs and are motivated by industrial uses. The objective of these notes is to introduce the reader to techniques of geometric optimal control as well as to provide an exposure to the applicability of numerical schemes implemented in HamPath [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF],

Bocop [START_REF] Bonnans | Bocop -A collection of examples[END_REF] and GloptiPoly [START_REF] Henrion | GloptiPoly: global optimization over polynomials with Matlab and SeDuMi[END_REF] software.

To highlight the main ideas and concepts, the presentation is restricted to the fundamental techniques and results. Moreover the selected applications drive the exposition of the different methodologies. They have received significant attention recently and are promising, paving the way for further research by our potential readers. The applications have been chosen based on the existence of accurate mathematical models to describe them, models that are suitable for a geometric analysis, and the possibility of implementing results from the analysis in a practical manner.

The notes are self-contained, moreover, the simpler geometric computations can be reproduced by the reader using our presentation of the maximum principle. The weak maximum principle covers the case of an open control domain which is the standard situation encountered in the classical calculus of variations, and is suitable for analysis of the first application, motility at low Reynolds number, although a good understanding of the so-called transversality conditions is necessary. For the second application, control of the spin dynamics by magnetic fields in nuclear magnetic resonance, the use of the general maximum principle is required since the control domain is bounded. At a more advanced level, the reader has to be familiar with the numerical techniques implemented in the software used for our calculations. In addition, symbolic methods have to be used to handle the more complex computations.

The first application is the swimming problem at low Reynolds number describing the swimming techniques of microorganisms. It can be easily observed in nature, but also mechanically reproduced using robotic devices, and it is linked to v vi Preface medical applications. This example serves as an introduction to geometric optimal control applied to sub-Riemannian geometry, a non-trivial extension of Riemannian geometry and a 1980's tribute of control theory to geometry under the influence of R. Brockett [START_REF] Brockett | Control theory and singular Riemannian geometry[END_REF]. We consider the Purcell swimmer [START_REF] Purcell | Life at low Reynolds number[END_REF], a three-link model where the shape variables are the two links at the extremities and the displacement is modeled by both the position and the orientation of the central link representing the body of the swimmer. To make a more complete analysis in the framework of geometric control, we use a simplified model from D. Takagi called the copepod swimmer [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF], where only line displacement is considered using symmetric shape links, and which is also the swimming model for an abundant variety of zooplankton (copepods). This fact is particularly interesting with respect to validating the correlation between the observed and predicted displacement using the mathematical model. Also from the mathematical point of view, the copepod model is the simplest slender body model and can be obtained as a limit case of more complex systems e.g. the Purcell swimmer.

For this problem, only the weak maximum principle is necessary and thus will be presented first, with its simple proof nevertheless containing all the geometric ingredients of the general maximum principle (see the historical paper by R.V. Gamkrelidze about the discovery of the maximum principle [START_REF] Gamkrelidze | Discovery of the maximum principle[END_REF]). Moreover, in this case, under proper regularity assumptions, the second-order conditions can be easily explained and numerically implemented using the concepts of conjugate points and the Jacobi equation. More specifically, using the optimal control framework, the sub-Riemannian problem is expressed as: dx dt

(t) = ∑ i=1,••• ,m u i (t)F i (x(t)), min u(.) T 0 ∑ i=1,••• ,m u 2 i (t)dt,
where x ∈ M, M is a smooth manifold, and the sub-Riemannian metric is defined by the orthonormal sub-frame {F 1 , • • • , F m } that determines the so-called nonholonomic constraints on the set of curves: ẋ(t) ∈ D(x(t)), where D is the distribution Span{F 1 , • • • , F m }. The relation to the swimming problem, modeled by some of the earliest prominent scientists (e.g. Newton, Euler and more recently Stokes, Reynolds, and Purcell), is straightforward in the framework of control theory. The state variable x of the system decomposes into (x 1 , x 2 ) where x 1 represents the displacement of the swimmer and x 2 is the shape variable representing the periodic shape deformation of the swimmer's body (called stroke) necessary to produce a net displacement for a given stroke. The mathematical model relates the speed of the displacement of x 1 to the speed of the shape deformation x 2 , thus characterizing the sub-Riemannian constraints, while the expended mechanical energy defines the metric. The model comes from hydrodynamics and is subject to vital approximations. First, at the scale of the micro-swimmer's life, it implies that inertia can be neglected [START_REF] Happel | Low Reynolds number hydrodynamics with special applications to particulate media[END_REF]. Second, according to the resistive force theory [START_REF] Hancock | The self-propulsion of microscopic organisms through liquids[END_REF] the interaction of the swimmer with the fluid is reduced to a drag force depending linearly upon the velocity. Finally, each of our swimmers is approximated by a slender body. Due to these approximations, experiments are crucial to validate the models. This theo-retical research also prompted experimentation using mechanical prototype devices (see for instance [START_REF] Or | Dynamics and stability of low-Reynolds-number swimming near a wall[END_REF]).

Our objective is to provide a self-contained presentation of the model, of the underlying concepts of sub-Riemannian geometry and of the techniques needed to conduct the mathematical analysis. The application of optimization techniques to the problem is recent. Our contribution's goal is to present a complete analysis using geometric and numerical techniques in the case of the copepod swimmer. It provides an excellent introduction to these methods, which have to be developed in the case of the Purcell swimmer based on our numerical results.

The second example concerns the optimal control of systems in nuclear magnetic resonance (NMR) and magnetic resonance imaging (MRI). The dynamics is modeled using the Bloch equation (1946), which describes at the macroscopic scale the evolution of the magnetization vector of a spin 1/2 particle depending on two relaxation parameters T 1 and T 2 , which are the chemical signatures of the chemical species (e.g. water, fat) and controlled by an Rf-magnetic pulse perpendicular to the strong polarizing field applied in the z-axis direction (see Bloch equations [START_REF] Bloch | Nuclear induction[END_REF]). At the experimental level, optimal control was introduced in the early 2000 in the dynamics of such systems, the objective being the replacement of the heuristic MRI pulse sequences used in hospital settings (in vivo), which means replacing in near future the standard industrial software by a new generation of software, producing a double gain: a better image in a shorter time.

Since the Bloch equations describe the evolution of the dynamics of the process with great accuracy and the computed control strategies can be implemented easily, this application provides an ideal platform to test the geometric optimal control framework presented in this volume. Clearly, the theory has to be developed to handle the mathematical problems and the analysis has to be supplemented by the use of a new generation of specific software dedicated to optimal control (HamPath [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF],

Bocop [START_REF] Bonnans | Bocop -A collection of examples[END_REF], GloptiPoly [START_REF] Henrion | GloptiPoly: global optimization over polynomials with Matlab and SeDuMi[END_REF]). With this in mind, the reader is introduced to two important problems in NMR and MRI. The first one is simply to saturate in minimum time the magnetization vector, which corresponds to driving its amplitude to zero. For this problem, we must first introduce the most general maximum principle since the applied Rf-magnetic field is of bounded amplitude. The second step is to compute, using geometric techniques, the structure of the optimal law as a closed loop function. This is the so-called concept of optimal synthesis. The second problem is the contrast problem in MRI where we must distinguish within a given picture between two heterogeneously distributed species, e.g. healthy versus cancer cells, that are characterized thanks to the Bloch equation by different responses to the same Rf-field due to different relaxation parameters. The actual MRI technology enables the transformation of this observation problem into an optimal control problem of the Mayer form: dx/dt(t) = f (x(t), u(t)) with |u(t)| ≤ M, min u(.) c(x(t f )), where t f is a fixed transfer time directly related to the image processing time and the cost function measures the contrast. The dynamics represents the coupling of the two Bloch equations controlled by the same Rf-field including the respective viii Preface parameters associated to the two species to be distinguished, parameters which can be easily computed experimentally.

We use three numerical software based on different approaches:

• Bocop. The so-called direct approach transforms the infinite dimensional optimal control problem into a finite dimensional optimization problem. This is done by a discretization in time applied to the state and control variables, as well as to the dynamics equation. These methods are usually less precise than indirect methods based on the maximum principle, but more robust with respect to the initialization. • Hampath. The HamPath software is based upon indirect methods: simple and multiple shooting; differential continuation (or homotopy) methods; and computation of the solutions of the variational equations needed to check the secondorder conditions of local optimality. Shooting methods consist in finding a zero of a specific function and use Newton-like algorithms. While simple shooting leads to solution of a two-point boundary value problem, multiple shooting takes into account intermediate conditions and the structure of the optimal solution has to be determined. This can be done using the Bocop software, which also allows initialization of the shooting equation. The Jacobian of the homotopic function is computed using variational equations to calculate the Jacobi fields. Moreover the Jacobi fields are used to check the necessary second-order optimality conditions. • LMI. The direct and indirect methods provide local optimal solutions. By comparing the different paths of zeros, one obtains a strong candidate solution whose global optimality must be analyzed. This can be done by the moment approach.

The moment approach is a global optimization technique that transforms a nonlinear, possibly infinite-dimensional optimization problem into convex and finitedimensional relaxations in the form of linear matrix inequalities (LMI). The first step consists in linearizing the problem by formulating it as a linear program on a measure space, a problem often referred to as a generalized moment problem. This can be performed by the use of so-called occupation measures, encoding admissible trajectories. The second step is to exploit the problem's structure, here given by its polynomial data, to manipulate the measures by their moment sequences. This leads to a semi-definite program (SDP) with countably many decision variables, one for each moment. The third and last step is to truncate this last problem to a finite set of those moments, leading to a relaxation in the form of LMI. The method is used through the GloptiPoly software. This approach is developed in the MRI problem thanks to the algebraic structure of Bloch equations and is crucial in this problem to discriminate the global optimum from the multiple local optimum solutions.

Numerical methods are supplemented by symbolic computations to handle or to check more complicated calculations. The combination of geometric, numerical and symbolic computations represents the main originality of the book and leads to the development of a modern and non-trivial computational framework. Another originality of the work presented here is its connection to real experiments. For the swimming problem, the copepod represents a variety of zooplankton observed at the University of Hawaii in Prof. Takagi's laboratory and is a model for the design of swimming robots. We represent in Fig. 0.1 the copepod observed by Takagi and the associated micro-robot model. For the MRI problem the numerical computations were implemented by Prof. Glaser at UTM in in vitro experiments and more recently in vivo experiments were performed at Creatis (INSA Lyon) by the group of Prof. Ratiney. On Figs. 0.2 and 0.3 are represented the in vivo and in vitro experiments realized in the project. Note that the numerical computations were performed using the Grape algorithm [START_REF] Lapert | Exploring the physical limits of saturation contrast in magnetic resonance imaging[END_REF].

Historical part -Calculus of variations

The calculus of variations is an old mathematical discipline and historically finds its origins in the introduction of the brachistochrone problem at the end of the 17th century by Johann Bernoulli to challenge his contemporaries to solve it. Here, we briefly introduce the reader to the main results. First, we introduce the fundamental formula of the classical calculus of variations following the presentation by Gelfand Fomin in [START_REF] Gelfand | Calculus of Variations[END_REF]. The originality of this presentation lies in the fact that it provides a general formula rather than starting with the standard Euler-Lagrange formula derivation and dealing with general variations. The fundamental formula leads to a derivation of the standard first order necessary conditions: Euler-Lagrange equation, tranversality conditions, Erdmann-Weierstrass conditions for a broken extremal and the Hamilton Jacobi equation. Second, we present a derivation of the second order necessary conditions in relation with the concept of conjugate points and the Jacobi equation. The main idea is to introduce the so-called accessory problem replacing the positivity test of the second order derivative by a minimization problem of the associated quadratic form [START_REF] Gelfand | Calculus of Variations[END_REF]. The modern interpretation in terms of the spectral theory of the associated self-adjoint operator (Morse theory) is also stated.

Statement of the Problem in the Holonomic Case

We consider the set C of all curves x : [t 0 ,t 1 ] → R n of class C 2 , where the initial and final times t 0 ,t 1 are not fixed, and the problem of minimizing a functional over C :

C(x) = t 1 t 0 L(t, x(t), ẋ(t)) dt
where L is C 2 . Moreover, we can impose extremity constraints: x(t 0 ) ∈ M 0 , x(t 1 ) ∈ M 1 where M 0 , M 1 are C1 -submanifolds of R n . The distance between two curves x(.), x * (.) is given by 1 ρ(x, x * ) = max t∈J∩J * x(t) -x * (t) + max t∈J∩J * ẋ(t) -ẋ * (t) + d(P 0 , P * 0 ) + d(P 1 , P * 1 )

where P 0 = (t 0 , x 0 ), P 1 = (t 1 , x 1 ), J, J * are the domains of x, x * and • is any norm on R n while d is the usual distance mapping on R n+1 . Note that a curve is interpreted in the time-extended space (t, x). If the two curves x(•), x * (•) are closed, they are by convention C 2 -extended on J ∪ J * .

Proposition 1 (Fundamental formula of the classical calculus of variations) We adopt the standard notation of classical calculus of variations, see [START_REF] Gelfand | Calculus of Variations[END_REF]. Let γ(•) be a reference curve with extremities (t 0 , x 0 ), (t 1 , x 1 ) and let γ(•) be any curve with extremities (t 0 + δt 0 , x 0 + δ x 0 ), (t 1 + δt 1 , x 1 + δ x 1 ). We denote by h(•) the variation:

h(t) = γ(t) -γ(t). Then, if we set ∆C = C( γ) -C(γ), we have ∆C = t 1 t 0 ∂ L ∂ x - d dt ∂ L ∂ ẋ | γ •h(t) dt + ∂ L ∂ ẋ | γ • δ x t 1 t 0 + L - ∂ L ∂ ẋ • ẋ | γ δt t 1 t 0 + o(ρ( γ, γ)) (1.1)
where • denotes the scalar product in R n .

Proof. We write

∆C = t 1 +δt 1 t 0 +δt 0 L(t, γ(t) + h(t), γ(t) + ḣ(t)) dt - t 1 t 0 L(t, γ(t), γ(t)) dt = t 1 t 0 L(t, γ(t) + h(t), γ(t) + ḣ(t)) dt - t 1 t 0 L(t, γ(t), γ(t)) dt + t 1 +δt 1 t 1 L(t, γ(t)+h(t), γ(t)+ ḣ(t)) dt - t 0 +δt 0 t 0 L(t, γ(t)+h(t), γ(t)+ ḣ(t)) dt.
We develop this expression using Taylor expansions keeping only the linear terms in h, ḣ, δ x, δt. We get

∆C = t 1 t 0 ∂ L ∂ x | γ • h(t) + ∂ L ∂ ẋ | γ • ḣ(t) | γ dt + [L(t, γ, γ)δt] t 1 t 0 + o(h, ḣ, δt).
The derivative of the variation ḣ is depending on h, and integrating by parts we obtain

∆C ∼ t 1 t 0 ∂ L ∂ x - d dt ∂ L ∂ ẋ | γ • h(t) dt + ∂ L ∂ ẋ | γ • h(t) t 1 t 0 + L | γ δt t 1 t 0 .
We observe that h, δ x, δt are not independent at the extremities and we have for t = t 0 or t = t 1 the relation

h(t + δt) = h(t) + o( ḣ, δt). So h(t) ∼ δ x -ẋδt.
Hence, we obtain the following approximation:

∆C ∼ t 1 t 0 ∂ L ∂ x - d dt ∂ L ∂ ẋ | γ • h(t) dt + ∂ L ∂ ẋ | γ • δ x t 1 t 0 + L - ∂ L ∂ ẋ ẋ |γ δt t 1 t 0
where all the quantities are evaluated along the reference trajectory γ(•). In this formula h, δ x, δt can be taken independent because in the integral the values h(t 0 ), h(t 1 ) do not play any special role.

From 1.1, we deduce that the standard first-order necessary conditions of the calculus of variations.

Corollary 1 Let us consider the minimization problem where the extremities (t 0 , x 0 ), (t 1 , x 1 ) are fixed. Then, a minimizer γ(•) must satisfy the Euler-Lagrange equation

∂ L ∂ x - d dt ∂ L ∂ ẋ | γ = 0. (1.2)
Proof. Since the extremities are fixed we set in (1.1) δ x = 0 and δt = 0 at t = t 0 and t = t 1 . Hence

∆C = t 1 t 0 ∂ L ∂ x - d dt ∂ L ∂ ẋ | γ • h(t) dt + o(h, ḣ)
for each variation h(•) defined on [t 0 ,t 1 ] such that h(t 0 ) = h(t 1 ) = 0. If γ(•) is a minimizer, we must have ∆C ≥ 0 for each h(•) and clearly by linearity, we have

t 1 t 0 ∂ L ∂ x - d dt ∂ L ∂ ẋ | γ • h(t) dt = 0 for each h(•). Since the mapping t → ( ∂ L ∂ x -d dt ∂ L ∂ ẋ ) | γ
is continuous, it must be identically zero along γ(•) and the Euler-Lagrange equation 1.2 is satisfied.

Hamiltonian Equations

The Hamiltonian formalism, which is the natural formalism to use for the maximum principle, appears in the classical calculus of variations via the Legendre transformation.

Definition 1

The Legendre transformation is defined by

p = ∂ L ∂ ẋ (t, x, ẋ) (1.3)
and if the mapping ϕ : (x, ẋ) → (x, p) is a diffeomorphism we can introduce the Hamiltonian:

H : (t, x, p) → p• ẋ -L(t, x, ẋ). (1.4)
Remark 1 In mechanics, the Lagrangian L is of the form V (x) -T (x, ẋ) where V is the potential and T is the kinetic energy and T is strictly convex with respect to ẋ.

Proposition 2

The formula (1.1) takes the form

∆C ∼ t 1 t 0 ∂ L ∂ x - d dt ∂ L ∂ ẋ | γ • h(t) dt + pδ x -Hδt t 1 t 0 (1.5)
and if γ(•) is a minimizer it satisfies the Euler-Lagrange equation in the Hamiltonian form

ẋ(t) = ∂ H ∂ p (t, x(t), p(t)), ṗ(t) = - ∂ H ∂ x (t, x(t), p(t)). (1.6)
Proof. Computing, one has

dH = ∂ H ∂t dt + ∂ H ∂ p dp + ∂ H ∂ x dx = (p - ∂ L ∂ ẋ )d ẋ + ẋdp - ∂ L ∂ x dx - ∂ L ∂t dt.

Hamilton-Jacobi-Bellman Equation

Definition 2 A solution of the Euler-Lagrange equation is called an extremal. Let P 0 = (t 0 , x 0 ) and P 1 = (t 1 , x 1 ). The Hamilton-Jacobi-Bellman (HJB) function is the multivalued function defined by S(P 0 , P 1 ) =

t 1 t 0 L(t, γ(t), γ(t)) dt
where γ(•) is any extremal with fixed extremities x 0 , x 1 . If γ(•) is a minimizer, it is called the value function.

Proposition 3 Assume that for each P 0 , P 1 there exists a unique extremal joining P 0 to P 1 and suppose that the HJB function is C 1 . Let P 0 be fixed and let S : P → S(P 0 , P). Then, S is a solution of the Hamilton-Jacobi-Bellman equation

∂ S ∂t (P 0 , P) + H(t, x, ∂ S ∂ x ) = 0. (1.7)
Proof. Let P = (t, x) and P + δ P = (t + δt, x + δ x). Denote by γ(•) the extremal joining P 0 to P and by γ(•) the extremal joining P 0 to P + δ P. We have

∆ S = S(t + dt, x + dx) -S(t, x) = C( γ) -C(γ)
and from (1.5) it follows that:

∆ S = ∆C ∼ t t 0 ∂ L ∂ x - d dt ∂ L ∂ ẋ | γ • h(t) dt + pδ x -Hδt t t 0 , where h(•) = γ(•) -γ(•). Since γ(•)
is a solution of the Euler-Lagrange equation, the integral is zero and

∆ S = ∆C ∼ pδ x -Hδt t t 0 .
In other words, we have

d S = pdx -H dt.
Identifying, we obtain

∂ S ∂t = -H, ∂ S ∂ x = p. (1.8)
Hence we get the HJB equation. Moreover p is the gradient to the level sets {x ∈

R n ; S(t, x) = c}.
Other consequences of the general formula are the so-called transversality and Erdmann Weierstrass (1877) conditions. They are presented in the exercises below.

Exercise 1.1. Consider the following problem: among all smooth curves t → x(t) whose extremity point P 1 = (t 1 , x 1 ) lies on a curve y(t) = Ψ (t), find the curve for which the functional t 1 t 0 L(t, x, ẋ) dt has an extremum. Deduce from the general formula the transversality conditions

L + L ẋ( Ψ -ẋ) = 0 at t = t 1 .
Exercise 1.2. Let t → x(t) be a minimizing solution of t 1 t 0 L(t, x, ẋ) dt with fixed extremities. Assume that t → x(t) is a broken curve with a corner at t = c ∈]t 0 ,t 1 [. Prove the Erdmann Weierstrass condition

L ẋ(c-) = L ẋ(c+), [L -L ẋ ẋ] (c-) = [L -L ẋ ẋ] (c+).
Give an interpretation using Hamiltonian formalism.

Second Order Conditions

The Euler-Lagrange equation has been derived using the linear terms in the Taylor expansion of ∆C. Using the quadratic terms we can get necessary and sufficient second order conditions. For the sake of simplicity, from now on we assume that the curves t → x(t) belong to R, and we consider the problem with fixed extremities: x(t 0 ) = x 0 , x(t 1 ) = x 1 . If the map L is taken C 3 , the second derivative is computed as follows:

∆C = t 1 t 0 L(t, γ(t) + h(t), γ(t) + ḣ(t)) -L(t, γ(t), γ(t)) dt = t 1 t 0 ∂ L ∂ x - d dt ∂ L ∂ ẋ | γ • h(t) dt + 1 2 t 1 t 0 ∂ 2 L ∂ x 2 | γ h 2 (t) + 2 ∂ 2 L ∂ x∂ ẋ | γ h(t) ḣ(t) + ∂ 2 L ∂ ẋ2 | γ ḣ2 (t) dt + o(h, ḣ) 2 .
If γ(t) is an extremal, the first integral is zero and the second integral corresponds to the intrinsic second-order derivative δ 2 C, that is:

δ 2 C = 1 2 t 1 t 0 ∂ 2 L ∂ x 2 | γ h 2 (t) + 2 ∂ 2 L ∂ x∂ ẋ | γ h(t) ḣ(t) + ∂ 2 L ∂ ẋ2 | γ ḣ2 (t) dt. (1.9)
Using h(t 0 ) = h(t 1 ) = 0, it can be written after an integration by parts as

δ 2 C = t 1 t 0 P(t) ḣ2 (t) + Q(t)h 2 (t) dt (1. 10
)
where

P = 1 2 ∂ 2 L ∂ ẋ2 | γ , Q = 1 2 ∂ 2 L ∂ x 2 - d dt ∂ 2 L ∂ x∂ ẋ | γ .
Using the fact that in the integral (1.10) the term P ḣ2 is dominating [START_REF] Gelfand | Calculus of Variations[END_REF], we get the following proposition.

Proposition 4 If γ(•) is a minimizing curve for the fixed extremities problem then it must satisfy the Legendre condition:

∂ 2 L ∂ ẋ2 | γ ≥ 0. (1.11)

The Accessory Problem and the Jacobi Equation

The intrinsic second-order derivative is given by

δ 2 C = t 1 t 0 P(t) ḣ2 (t) + Q(t)h 2 (t) dt, h(t 0 ) = h(t 1 ) = 0,
where P, Q are as above. It can be written as

δ 2 C = t 1 t 0 (P(t) ḣ(t)) ḣ(t) + (Q(t)h(t))h(t) dt
and integrating by parts using h(t 0 ) = h(t 1 ) = 0, we obtain

δ 2 C = t 1 t 0 Q(t)h(t) - d dt (P(t) ḣ(t)) h(t) dt.
Let us introduce the linear operator D : h → Qhd dt (P ḣ). Hence, we can write

δ 2 C = (D h, h) (1.12)
where (, ) is the usual scalar product on L 2 ([t 0 ,t 1 ]). The linear operator D is called the Euler-Lagrange operator.

Definition 3 From (1.12), δ 2 C is a quadratic operator on the set C 0 of C 2 -curves h : [t 0 ,t 1 ] → R satisfying h(t 0 ) = h(t 1 ) = 0, h = 0.
The so-called accessory problem is:

min h∈C 0 δ 2 C.

Definition 4

The Euler-Lagrange equation corresponding to the accessory problem is called the Jacobi equation and is given by

D h = 0 (1.13)
where D is the Euler-Lagrange operator: D h = Qhd dt (P ḣ). It is a second-order linear differential operator.

Definition 5

The strong Legendre condition is P > 0, where P =

1 2 ∂ 2 L ∂ ẋ2 | γ
. If this condition is satisfied, the operator D is said to be nonsingular.

Conjugate Point and Local Morse Theory

In this section, we present some results from [START_REF] Gregory | Quadratic form theory and differential equations[END_REF] and [START_REF] Milnor | Morse theory[END_REF].

Definition 6 Let γ(•) be an extremal. A solution J(•) ∈ C 0 of D J = 0 on [t 0 ,t 1 ] is called a Jacobi curve. If there exists a Jacobi curve along γ(•) on [t 0 ,t 1 ] the point γ(t 1
) is said to be conjugate to γ(t 0 ).

Theorem 1 (Local Morse theory [START_REF] Milnor | Morse theory[END_REF]). Let t 0 be fixed and let us consider the Euler-Lagrange operator (indexed by t > t 0 ) D t defined on the set C t 0 of curves on [t 0 ,t] satisfying h(t 0 ) = h(t) = 0. By definition, a Jacobi curve on [t 0 ,t] corresponds to an eigenvector J t associated to an eigenvalue λ t = 0 of D t . If the strong Legendre condition is satisfied along an extremal γ : [t 0 , T ] → R n , we have a precise description of the spectrum of D t as follows. There exists t 0 < t 1 < • • • < t s < T such that each γ(t i ) is conjugate to γ(t 0 ). If n i corresponds to the dimension of the set of the Jacobi curves J t i associated to the conjugate point γ(t i ), then for any T such that t

0 < t 1 < • • • < t k < T < t k+1 < • • • < T we have the identity n - T = k ∑ i=1 n i (1.14)
where n - T = dim{linear space of eigenvectors of D T corresponding to strictly negative eigenvalues }. In particular if T > t 1 we have

min h∈C 0 T t 0 (Q(t)h 2 (t) + P(t) ḣ2 (t)) dt = -∞.
(1.15)

1.7 From calculus of variations to optimal control theory and Hamiltonian dynamics

An important and difficult problem is to generalize the first and second order conditions from classical calculus of variations to optimal control theory (OCT).

In OCT, the problem is stated as dq dt = F(q, u) min

u(•) T 0 L(q, u) dt
with smooth data but the set of admissible controls is the set of bounded measurable mappings valued in a control domain U, thus the set of admissible trajectories is the set of absolutely continuous curves t → q(•). Minimizers are found among the set of extremals (q, p, u) solutions of the Hamilton equations

dq dt = ∂ H ∂ p , dp dt = - ∂ H ∂ q (1.16)
where H is the so-called unmaximized Hamiltonian

H(q, p, u) = p• F(q, u)
where the controls have to satisfy the maximization condition H(q, p, u) = max v∈U H(q, p, v).

(1.17)

1.7 From calculus of variations to optimal control theory and Hamiltonian dynamics Solving this equation leads in general to several true Hamiltonian H i (q, p), i = 1, . . . , k and the optimal solution is found by concatenation of trajectories of the vector fields -→ H i 's.

Remark 1.1. OCT is a non trivial extension to the so-called Lagrange problem in calculus of variations since there exists no restriction of the control domain.

The concept of conjugate points can be extended in optimal control and is related to losing optimality for some prescribed topology on the set of controls but practical computation is intricate.

A major problem in the analysis is due to bad controllability properties of the so-called abnormal trajectories. This problem stopped further developments of calculus of variations in the forties [START_REF] Bliss | Lectures on the Calculus of Variations[END_REF]. It was revived recently in optimal control theory when dealing with SR-geometry and more geometrically investigated, see for instance [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF].

Also in the frame of Hamiltonian formulation of the Maximum Principle defined by (1.16), (1.17) a bridge is open between Hamiltonian dynamics and variational principles. Indeed Hamiltonian and Lagrangian can be related with some regularity assumptions using the Legendre-Fenchel transform

H(q, p) = max v (p• v -L(x, v))
and interaction between Hamiltonian dynamics and optimal control is a rich and active domain, see for instance [START_REF] Agrachev | On curvatures and focal points of dynamical Lagrangian distributions and their reductions by first integrals[END_REF], [START_REF] Arcostanzo | Tonelli Hamiltonians without conjugate points and C 0 -integrability[END_REF].

Chapter 2

Weak Maximum Principle and Application to Swimming at low Reynolds Number

Pre-requisite of Differential and Symplectic Geometry

We refer to [START_REF] Helgason | Differential geometry, Lie groups, and symmetric spaces[END_REF][START_REF] Arnol˜d | Mathematical methods of classical mechanics[END_REF][START_REF] Godbillon | Geométrie différentielle et mécanique analytique[END_REF] for more details about the general concepts and notations introduced in this section.

Notations. Let M be a smooth (C ∞ or C ω ) connected and second-countable manifold of dimension n. We denote by T M the fiber bundle and by T * M the cotangent bundle. Let V (M) be the set of smooth vector fields on M and Di f f (M) the set of smooth diffeomorphisms.

Definition 7 Let X ∈ V (M) and let f be a smooth function on M. The Lie derivative is defined as:

L X f = d f (X). If X,Y ∈ V (M), the Lie bracket is given by ad X(Y ) = [X,Y ] = L Y • L X -L X • L Y . If x = (x 1 , • • • , x n ) is a local system of coordinates we have: X(x) = n ∑ i=1 X i (x) ∂ ∂ x i L X f (x) = ∂ f ∂ x X(x) [X,Y ](x) = ∂ X ∂ x (x)Y (x) - ∂Y ∂ x (x)X(x).
The mapping (X,Y ) → [X,Y ] is R-linear and skew-symmetric. Moreover, the Jacobi identity holds:

[X, [Y, Z]] + [Y, [Z, X]] + [Z, [X,Y ]] = 0.
Definition 8 Let X ∈ V (M). We denote by x(t, x 0 ) the maximal solution of the Cauchy problem ẋ(t) = X(x(t)), x(0) = x 0 . This solution is defined on a maximal open interval J containing 0. We denote by exptX the local one parameter group associated to X, that is: exptX(x 0 ) = x(t, x 0 ). The vector field X is said to be complete if the trajectories can be extended over R.

Definition 9 Let X ∈ V (M) and ϕ ∈ Di f f (M). The image of X by ϕ is ϕ * X = dϕ(X • ϕ -1 ).
We recall the following results.

Proposition 5 Let X,Y ∈ V (M) and ϕ ∈ Di f f (M). We have:

1. The one parameter local group of Z = ϕ * X is given by:

exptZ = ϕ • exptX • ϕ -1 . 2. ϕ * [X,Y ] = [ϕ * X, ϕ * Y ].
3. The Baker-Campbell-Hausdorff (BCH) formula is:

exp sX exptY = exp ζ (X,Y )
where ζ (X,Y ) belongs to the Lie algebra generated by [X,Y ] with:

ζ (X,Y ) = sX + tY + st 2 [X,Y ] + st 2 12 [[X,Y ],Y ] - s 2 t 12 [[X,Y ], X] - s 2 t 2 24 [X, [Y, [X,Y ]]] + • • • .
The series is converging for s,t small enough in the analytic case. 4. We have

exptX exp εY exp -tX = exp η(X,Y ) with η(X,Y ) = ε ∑ k≥0 t k k! ad k X(Y )
and the series converging for ε,t small enough in the analytic case. 5. The ad-formula is:

exptX * Y = ∑ k≥0 t k k! ad k X(Y )
where the series is converging for t small enough.

Definition 10 Let V be a R-linear space of dimension 2n. This space is said to be symplectic if there exists an application ω : V × V → R which is bilinear, skewsymmetric and nondegenerate, that is: if ω(x, y) = 0 for all x ∈ V , then y = 0. Let W be a linear subspace of V . We denote by W ⊥ the set

W ⊥ = {x ∈ V ; ω(x, y) = 0 ∀y ∈ W }. The space W is isotropic if W ⊂ W ⊥ . An isotropic space is said to be Lagrangian if dimW = dim V 2 . Let (V 1 , ω 1 ), (V 2 , ω 2 ) be two symplectic linear spaces. A linear mapping f : V 1 → V 2 is symplectic if ω 1 (x, y) = ω 2 ( f (x), f (y)) for each x, y ∈ V 1 .
Proposition 6 Let (V, ω) be a linear symplectic space. Then there exists a basis {e 1 , ..., e n , f 1 , ..., f n } called canonical defined by ω(e i , e j ) = ω( f i , f j ) = 0 f or 1 ≤ i, j ≤ n and ω(e i , f j ) = δ i j (Kronecker symbol). If J is the matrix 0 I -I 0 where I is the identity matrix of order n, then we can write ω(x, y) = Jx, y where , is the scalar product (in the basis (e i , f j )). In the canonical basis, the set of all linear symplectic transformations is represented as the symplectic group defined by Sp(n, R) = {S ∈ GL(2n, R); S JS = J}.

Definition 11 Let M be a C ∞ -manifold of dimension 2n. A symplectic structure on M is defined by a 2-form ω such that dω = 0 and such that ω is regular, that is: ∀x ∈ M, ω x is nondegenerate.

Proposition 7 For any C ∞ -manifold of dimension n, the cotangent bundle T * M admits a canonical symplectic structure defined by ω = dα where α is the Liouville form. If x = (x 1 , ..., x n ) is a coordinate system on M and (x, p) with (p 1 , ..., p n ) the associated coordinates on T * M, the Liouville form is written locally as α = ∑ n i=1 p i dx i and ω = dα = ∑ n i=1 dp i ∧ dx i . Proposition 8 (Darboux) Let (M, ω) be a symplectic manifold. Then given any point in M, there exists a local system of coordinates called Darboux coordinates, (x 1 , ..., x n , p 1 , ..., p n ) such that ω is given locally by ∑ n i=1 dp i ∧ dx i . (Hence the symplectic geometry is a geometry with no local invariant).

Definition 12 Let (M, ω) be a symplectic manifold and let X be a vector field on M. We note i X ω the interior product defined by i X ω(Y ) = ω(X,Y ) for any vector field Y on M. Let H : M → R a real-valued function. The vector field denoted by -→ H and defined by i-→ H (ω) = -dH is the Hamiltonian vector field associated to H. If (x, p) is a Darboux coordinate system, then the Hamiltonian vector field is expressed in these coordinates as:

- → H = n ∑ i=1 ∂ H ∂ p i ∂ ∂ x i - ∂ H ∂ x i ∂ ∂ p i .
Definition 13 Let F, G : M → R be two mappings. We denote by {F, G} the Poissonbracket of F and G defined by {F, G} = dF( -→ G ).

Proposition 9 (Properties of the Poisson-bracket)

1. The mapping (F, G) → {F, G} is bilinear and skew-symmetric.

The Leibniz identity holds:

{FG, H} = G{F, H} + F{G, H}.

3. In a Darboux coordinate system, we have

{F, G} = n ∑ i=1 ∂ G ∂ p i ∂ F ∂ x i - ∂ G ∂ x i ∂ F ∂ p i . 4. If the Lie bracket is defined by [ - → F , - → G ] = - → G • - → F - - → F • - →
G , then its relation with the Poisson bracket is given by: [

- → F , - → G ] = { --→ F, G}.
→ H X = X ∂ ∂ x -∂ X ∂ x p ∂ ∂ p .
Each diffeomorphism ϕ on M can be lifted into a symplectic diffeomorphism -→ ϕ on T * M defined in a local system of coordinates as follows. If

x = ϕ(y), then - → ϕ : (y, q) → (x, p) = ϕ(y), ∂ ϕ -1 ∂ y q .
Theorem 2 (Noether) Let (x, p) be Darboux coordinates on T * M, X a vector field on M and -→ H X its Hamiltonian lift. We assume -→ H X to be a complete vector field and we denote by ϕ t the associated one parameter group. Let F : T * M → R and let us assume that F • ϕ t = F for all t ∈ R. Then H X is a first integral for -→ F .

Definition 17 Let M be a manifold of dimension 2n + 1 and let ω be a 2-form on M. Then for all x ∈ M, ω x is bilinear , skew-symmetric and its rank is ≤ 2n. If for each x, the rank is 2n, we say that ω is regular. In this case ker ω is of rank one and is generated by an unique vector field X up to a scalar. If α is a 1-form such that dα is of rank 2n, the vector field associated with dα is called the characteristic vector field of α and the trajectories of X are called the characteristics.

Proposition 10 On the space T * M × R with coordinates (x, p,t) the characteristics of the 1-form ∑ n i=1 (p i dx i -H dt) project onto solutions of the Hamilton equations:

ẋ(t) = ∂ H ∂ p (x(t), p(t),t), ṗ(t) = - ∂ H ∂ x (x(t), p(t),t).
Definition 18 Let ϕ : (x, p,t) → (X, P, T ) be a change of coordinates on T * M × R.

If there exist two functions K(X, P, T ) and S(X, P, T ) such that p dx -H dt = P dX -K dT + dS, then the mapping ϕ is a canonical transformation and S is called the generating function of ϕ.

Proposition 11 For a canonical transformation the equations

ẋ(t) = ∂ H ∂ p (x(t), p(t),t), ṗ(t) = - ∂ H ∂ x (x(t), p(t),t) transform onto dX dT (T ) = ∂ K ∂ P (X(T ), P(T ), T ), dP dT (T ) = - ∂ K ∂ X (X(T ), P(T ), T ). If T = t,
H(x, ∂ S ∂ x ) = K(X).
Since X(t) = (X 1 (0), . . . , X n (0)) is fixed, if we can integrate the previous equation we get solutions to the Hamilton equations. A standard method is by separating the variables. This is called the Jacobi method to integrate the Hamilton equations. In particular, this leads to a classification of integrable mechanical systems in small dimension, see [START_REF] Landau | Physique théorique[END_REF].

Definition 19 A polysystem D is a family {V i ; i ∈ I} of vector fields. We denote by the same letter the associated distribution, that is the mapping x → span{V (x);V ∈ D}. The distribution D is said to be involutive if

[V i ,V j ] ⊂ D, for all V i ,V j ∈ D.
Definition 20 Let D be a polysystem. We design by D L.A. the Lie algebra generated by D, it is constructed recursively as follows:

D 1 = span{D}, D 2 = span{D 1 + [D 1 , D 1 ]}, . . . , D k = span{D k-1 + [D 1 , D k-1 ]} and D L.A. = ∪ k≥1 D k .
By construction the associated distribution D L.A. is involutive. If x ∈ M, we associate the following sequence of integers:

n k (x) = dim D k (x).
Definition 21 Consider a control system ẋ = f (x, u) on M with u ∈ U. We can associate to this system the polysystem D = { f (•, u); u constant, u ∈ U}. We denote by S T (D) the set

S T (D) = {expt 1 V 1 • • • expt k V k ; k ∈ N, t i ≥ 0 and k ∑ i=1 t i = T, V i ∈ D}
and by S(D) the local semi-group: ∪ T ≥0 S T (D). We denote by G(D) the local group generated by S(D), that is

G(D) = {expt 1 V 1 • • • expt k V k ; k ∈ N, t i ∈ R, V i ∈ D}.
Properties.

1. The accessibility set from x 0 in time T is:

A(x 0 , T ) = S T (D)(x 0 ).
2. The accessibility set from x 0 is the orbit of the local semi-group:

A(x 0 ) = S(D)(x 0 ).

Definition 22

We call the orbit of x 0 the set O(x 0 ) = G(D)(x 0 ). The system is said to be weakly controllable if for every x 0 ∈ M, O(x 0 ) = M.

Controllability Results

Sussmann-Nagano Theorem

When the rank condition is satisfied (rank ∆ = constant, ∆ : x → D L.A. (x)) we get from the Frobenius theorem a description of all the integral manifolds near x 0 . If we only need to construct the leaf passing through x 0 the rank condition is clearly too strong. Indeed, if D = {X} is generated by a single vector field X, there exists an integral curve through x 0 . For a family of vector fields this result is still true if the vector fields are analytic.

Theorem 3 (Nagano-Sussman Theorem [START_REF] Sussmann | Orbits of families of vector fields and integrability of distributions[END_REF]) Let D be a family of analytic vector fields near x 0 ∈ M and let p be the rank of ∆ : x → D L.A. (x) at x 0 . Then through x 0 there exists locally an integral manifold of dimension p.

Proof. Let p be the rank of ∆ at x 0 . Then there exists p vector fields of D L.A. :

X 1 , • • • , X p such that span{X 1 (x 0 ), • • • , X p (x 0 )} = ∆ (x 0 ). Consider the map α : (t 1 , • • • ,t p ) → expt 1 X 1 • • • expt p X p (x 0 ).

It is an immersion for

(t 1 , • • • ,t p ) = (0, • • • , 0)
. Hence the image denoted by N is locally a submanifold of dimension p. To prove that N is an integral manifold we must check that for each y ∈ N near x 0 , we have T y N = ∆ (y). This result is a direct consequence of the equalities

D L.A. (exptX i (x)) = d exptX i (D L.A. (x)), i = 1, • • • , p
for x near x 0 , and t small enough. To show that the previous equalities hold, let

V (x) ∈ D L.A. (x) such that V (x) = Y (x) with Y ∈ D L.A.
. By analycity and the adformula for t small enough we have

(d exptX i )(Y (x)) = ∑ k≥0 t k k! ad k X i (Y )(exptX i (x)).
Hence for t small enough, we have

(d exptX i )(D L.A. (x)) ⊂ D L.A. (exptX i (x)).
Changing t to -t we show the second inclusion.

C ∞ -Counter Example

To prove the previous theorem we use the following geometric property. Let X,Y be two analytic vector fields and assume X(x 0 ) = 0. From the ad-formula, if all the vector fields ad k X(Y ), k ≥ 0 are collinear to X at x 0 , then for t small enough the vector field Y is tangent to the integral curve exptX(x 0 ). Hence is is easy to construct a C ∞ -counter example using flat C ∞ -mappings. Indeed, let us take f : R → R a smooth map such that f (x) = 0 for x ≤ 0 and f (x) = 0 for x > 0. Consider the two vector fields on R 2 :

X = ∂ ∂ x and Y = f (x) ∂ ∂ y . At 0, D L.A. is of rank 1. Indeed, we have [X,Y ](x) = -f (x) ∂
∂ y = 0 at 0 and hence [X,Y ](0) = 0. The same is true for all high order Lie brackets. In this example the rank D L.A. is not constant along exptX(0), indeed for x > 0, the vector field Y is transverse to this vector field.

Chow-Rashevskii Theorem

Theorem 4 ( [START_REF] Chow | Über Systeme von linearen partiellen Differentialgleichungen erster Ordnung[END_REF][START_REF] Rashevskii | About connecting two points of complete non-holonomic space by admissible curve[END_REF]) Let D be a C ∞ -polysystem on M. We assume that for each x ∈ M, D L.A. (x) = T x M. Then we have

G(D)(x) = G(D L.A. (x)) = M, for each x ∈ M.
Proof. Since M is connected it is sufficient to prove the result locally. The proof is based on the BCH-formula. We assume M = R 3 and D = {X,Y } with rank {X,Y, [X,Y ]} = 3 at x 0 ; the generalization is straightforward. Let λ be a real number and consider the map

ϕ λ : (t 1 ,t 2 ,t 3 ) → exp λ X expt 3 Y exp -λ X expt 2 Y expt 1 X(x 0 ).
We prove that for small but nonzero λ , ϕ λ is an immersion. Indeed, using the BCH formula we have

ϕ λ (t 1 ,t 2 ,t 3 ) = exp(t 1 X + (t 2 + t 3 )Y + λt 3 2 [X,Y ] + • • • )(x 0 ), hence ∂ ϕ λ ∂t 1 (0, 0, 0) = X(x 0 ), ∂ ϕ λ ∂t 2 (0, 0, 0) = Y (x 0 ), ∂ ϕ λ ∂t 3 (0, 0, 0) = Y (x 0 ) + λ 2 [X,Y ](x 0 ) + o(λ ).
Since X,Y, [X,Y ] are linearly independent at x 0 , the rank of ϕ λ at 0 is 3 for λ = 0 small enough.

Weak Maximum Principle

We consider the autonomous control system

ẋ(t) = f (x(t), u(t)), x(t) ∈ R n , u(t) ∈ Ω (2.1)
where f is a C 1 -mapping. The initial and target sets M 0 , M 1 are given and we assume they are C 1 -submanifolds of R n . The control domain is a given subset Ω ⊂ R m . The class of admissible controls U is the set of bounded measurable maps u : [0, T (u)] → Ω . Let u(•) ∈ U and x 0 ∈ R n be fixed. Then, by the Caratheodory theorem [START_REF] Lee | Foundations of optimal control theory[END_REF], there exists a unique trajectory of (2.1) denoted x(•, x 0 , u) such that x(0) = x 0 . This trajectory is defined on a nonempty subinterval J of [0, T (u)] and t → x(t, x 0 , u) is an absolutely continuous function solution of (2.1) almost everywhere.

To each u(•) ∈ U defined on [0, T ] with corresponding trajectory x(•, x 0 , u) issued from x(0) = x 0 ∈ M 0 defined on [0, T ], we assign a cost

C(u) = T 0 f 0 (x(t), u(t)) dt (2.2)
where f 0 is a C 1 -mapping. An admissible control u * (•) with corresponding trajectory x * (•, x 0 , u) and defined on [0, T * ] such that x * (0) ∈ M 0 and x * (T * ) ∈ M 1 is optimal if for each admissible control u(•) with corresponding trajectory x(•, x 0 , u) on [0, T ], x(0) ∈ M 0 and x(T ) ∈ M 1 , then

C(u * ) ≤ C(u).

The Augmented System

The following remark is straightforward but is geometrically very important to understand the maximum principle. Let us consider f = ( f , f 0 ) and the corresponding system on R n+1 defined by the equations ẋ = f ( x(t), u(t)), i.e.:

ẋ(t) = f (x(t), u(t)), (2.3)
ẋ0 (t) = f 0 (x(t), u(t)). (2.4)
This system is called the augmented system. Since f is C 1 , according to the Caratheodory theorem, to each admissible control u(•) ∈ U there exists an admissible trajectory x(t, x0 , u) such that x0 = (x 0 , x 0 (0)), x 0 (0) = 0 where the added coordinate x 0 (•) satisfies x 0 (T ) = T 0 f 0 (x(t), u(t)) dt. Let us denote by ÂM 0 the accessibility set ∪ u(•)∈U x(T, x0 , u) from M0 = (M 0 , 0) and let M1 = M 1 × R. Then, we observe that an optimal control u * (•) corresponds to a trajectory x * (•) such that x * ∈ M0 and intersecting M1 at a point x * (T * ) where x 0 is minimal. In particular x * (T ) belongs to the boundary of the Accessibility set ÂM 0 .

Related Problems

Our framework is a general setting to deal with a large class of problems. Examples are the following:

1. Nonautonomous systems:

ẋ(t) = f (t, x(t), u(t)).
We add the variable t to the state space by setting dt ds = 1,t(s 0 ) = s 0 . 2. Fixed time problem. If the time domain [0, T (u)] is fixed (T (u) = T for all u(•))

we add the variable t to the state space by setting dt ds = 1,t(s 0 ) = s 0 and we impose the following state constraints on t : t = 0 at s = 0 and t = T at the free terminal time s. Some specific problems important for applications are the following.

1. If f 0 ≡ 1, then min T 0 f 0 (x(t), u(t)) dt = min T and we minimize the time of global transfer. 2. If the system is of the form: ẋ(t) = A(t)x(t) + B(t)u(t), where A(t), B(t) are matrices and C(u) = T 0 L(t, x(t), u(t)) dt where L(•, x, u) is a quadratic form for each t, T being fixed, the problem is called a linear quadratic problem (LQ-problem).

Singular Trajectories and the Weak Maximum Principle

Definition 23 Consider a system of R n : ẋ(t) = f (x(t), u(t)) where f is a C ∞ -map from R n × R m into R n . Fix x 0 ∈ R n and T > 0. The end-point map (for fixed x 0 , T ) is the map E x 0 ,T : u(•) ∈ U → x(T, x 0 , u). If u(•) is a control defined on [0, T ] such that the corresponding trajectory x(•, x 0 , u) is defined on [0, T ], then E x 0 ,T is defined on a neighborhood V of u(•) for the L ∞ ([0, T ]) norm.

First and Second Variations of E x 0 ,T It is a standard result, see for instance [START_REF] Sontag | Mathematical control theory. Deterministic finite-dimensional systems[END_REF], that the end-point map is a C ∞ -map defined on a domain of the Banach space L ∞ ([0, T ]). The formal computation of the successive derivatives uses the concept of Gâteaux derivative. Let us explain in details the process to compute the first and second variations.

Let v(•) ∈ L ∞ ([0, T ]) be a variation of the reference control u(•) and let us denote by x(•) + ξ (•) the trajectory issued from x 0 and corresponding to the control u(•) + v(•). Since f is C ∞ , it admits a Taylor expansion for each fixed t:

f (x + ξ , u + v) = f (x, u) + ∂ f ∂ x (x, u)ξ + ∂ f ∂ u (x, u)v + ∂ 2 f ∂ x∂ u (x, u)(ξ , v) + 1 2 ∂ 2 f ∂ x 2 (x, u)(ξ , ξ ) + 1 2 ∂ 2 f ∂ u (x, u)(v, v) + • • •
Using the differential equation we get

ẋ(t) + ξ (t) = f (x(t) + ξ (t), u(t) + v(t)).
Hence we can write ξ as:

δ 1 x + δ 2 x + • • • where δ 1 x is linear in v, δ 2
x is quadratic, etc. and are solutions of the following differential equations:

δ 1 x = ∂ f ∂ x (x, u)δ 1 x + ∂ f ∂ u (x, u)v (2.5) δ 2 x = ∂ f ∂ x (x, u)δ 2 x + ∂ 2 f ∂ x∂ u (x, u)(δ 1 x, v) + 1 2 ∂ 2 f ∂ x 2 (x, u)(δ 1 x, δ 2 x) + 1 2 ∂ 2 f ∂ u 2 (x, u)(v, v). (2.6)
Using ξ (0) = 0, these differential equations have to be integrated with the initial conditions

δ 1 x(0) = δ 2 x(0) = 0. (2.7)
Let us introduce the following notations:

A(t) = ∂ f ∂ x (x(t), u(t)), B(t) = ∂ f ∂ u (x(t), u(t)).

Definition 24

The system

δ x(t) = A(t)δ x(t) + B(t)δ u(t)
is called the linearized system along (x(•), u(•)).

Let M(t) be the fundamental matrix on [0, T ] solution almost everywhere of

Ṁ(t) = A(t)M(t), M(0) = Id.
Integrating (2.5) with δ 1 x(0) = 0 we get the following expression for δ 1 x:

δ 1 x(T ) = M(T ) T 0 M -1 (t) B(t) v(t) dt. (2.8)
This implies the following lemma.

Lemma 1 The Fréchet derivative of E x 0 ,T at u(•) is given by

E x 0 ,T (v) = δ 1 x(T ) = M(T ) T 0 M -1 (t)B(t)v(t) dt.

Definition 25

The admissible control u(•) and its corresponding trajectory x(•, x 0 , u) both defined on [0, T ] are said to be regular if the Fréchet derivative E x 0 ,T is surjective. Otherwise they are called singular.

Proposition 12 Let A(x 0 , T ) = ∪ u(•)∈U x(T, x 0 , u) be the accessibility set at time T from x 0 . If u(•) is a regular control on [0, T ], then there exists a neighborhood U of the end-point x(T, x 0 , u) contained in A(x 0 , T ).

Proof. Since E x 0 ,T is surjective at u(•), we have using the open mapping theorem that E x 0 ,T is an open map.

Theorem 5 Assume that the admissible control u(•) and its corresponding trajectory x(•) are singular on [0, T ]. Then there exists a vector p(•) ∈ R n \ {0} absolutely continuous on [0, T ] such that (x, p, u) are solutions almost everywhere on [0, T ] of the following equations:

dx dt (t) = ∂ H ∂ p (x(t), p(t), u(t)), dp dt (t) = - ∂ H ∂ x (x(t), p(t), u(t)) (2.9) ∂ H ∂ u (x(t), p(t), u(t)) = 0 (2.10)
where H(x, p, u) = p, f (x, u) is the pseudo-Hamiltonian, , being the standard inner product.

Proof. We observe that the Fréchet derivative is a solution of the linear system

δ x(t) = A(t)δ 1 x(t) + B(t)v(t).
Hence, if the pair (x(•), u(•)) is singular this system is not controllable on [0, T ].

We use an earlier proof on controllability to get a geometric characterization of this property. The proof which is the heuristic basis of the maximum principle is given in detail. By definition, since u(•) is a singular control on [0, T ] the dimension of the linear space

T 0 M(T )M -1 (t) B(t) v(t) dt; v(•) ∈ L ∞ ([0, T ])
is less than n. Therefore there exists a row vector p ∈ R n \ {0} such that

pM(T )M -1 (t)B(t) = 0
for almost everywhere t ∈ [0, T ]. We set

p(t) = pM(T )M -1 (t).
By construction p(•) is a solution of the adjoint system

ṗ(t) = -p(t) ∂ f ∂ x (x(t), u(t)).
Moreover, it satisfies almost everywhere the following equality:

p(t) ∂ f ∂ u (x(t), u(t)) = 0.
Hence we get the equations (2.9) and (2.10) if H(x, p, u) denotes the scalar product p, f (x, u) .

Geometric interpretation of the Adjoint Vector

In the proof of Theorem 5 we introduced a vector p(•). This vector is called an adjoint vector. We observe that if u(•) is singular on [0, T ], then for each 0

< t ≤ T, u | [0,T ] is singular and p(t) is orthogonal to the image denoted K(t) of E x 0 ,T evaluated at u | [0,t] . If for each t, K(t)
is a linear space of codimension one then p(t) is unique up to a factor.

The Weak Maximum Principle

Theorem 6 Let u(•) be a control and x(•, x 0 , u) the corresponding trajectory, both defined on [0, T ]. If x(T, x 0 , u) belongs to the boundary of the accessibility set A(x 0 , T ), then the control u(•) and the trajectory x(•, x 0 , u) are singular.

Proof. According to Proposition 12, if u(•) is a regular control on [0, T ] then x(T ) belongs to the interior of the accessibility set.

Corollary 2 Consider the problem of maximizing the transfer time for system

ẋ(t) = f (x(t), u(t)), u(•) ∈ U = L ∞ , with fixed extremities x 0 , x 1 . If u * (•)
and the corresponding trajectory are optimal on [0, τ * ], then u * (•) is singular.

Proof. If u * (•) is maximizing then x * (T ) must belong to the boundary of the accessibility set A(x 0 , T ) otherwise there exists ε > 0 such that x * (T -ε) ∈ A(x 0 , T ) and hence can be reached by a solution x(•) in time T :

x * (T -ε) = x(T ).
It follows that the point x * (T ) can be joined in a time T > T . This contradicts the maximality assumption.

Corollary 3 Consider the system ẋ(t) = f (x(t), u(t)) where u(•) ∈ U = L ∞ ([0, T ])
and the minimization problem: min

u(•)∈U T 0 L(x(t), u(t))
dt, where the extremities x 0 , x 1 are fixed as well as the transfer time T . If u * (•) and its corresponding trajectory are optimal on [0, T ], then u * (•) is singular on [0, T ] for the augmented system:

ẋ(t) = f (x(t), u(t)), ẋ0 (t) = L(x(t), u(t)). Therefore there exists p * (t) = (p(t), p 0 ) ∈ R n+1 \ {0} such that ( x * , p * , u * ) satisfies ẋ(t) = ∂ Ĥ ∂ p ( x(t), p(t), u(t)), ṗ(t) = - ∂ Ĥ ∂ x ( x(t), p(t), u(t)) ∂ Ĥ ∂ u ( x(t), p(t), u(t)) = 0 (2.11)
where x = (x, x 0 ) and Ĥ( x, p, u) = p, f (x, u) + p 0 L(x, u). Moreover p 0 is a nonpositive constant.

Proof. We have that x * (T ) belongs to the boundary of the accessibility set Â( x0 , T ). Applying (2.9), (2.10) we get the equations (2.11) where ṗ0

= -∂ Ĥ ∂ x 0 = 0 since Ĥ is independent of x 0 .

Abnormality

In the previous corollary, p * (•) is defined up to a factor. Hence we can normalize p 0 to 0 or -1 and we have two cases: Case 1: u(•) is regular for the system ẋ(t) = f (x(t), u(t)). Then p 0 = 0 and can be normalized to -1. This is called the normal case (in calculus of variations), see [START_REF] Bliss | Lectures on the Calculus of Variations[END_REF].

Case 2: u(•) is singular for the system ẋ(t) = f (x(t), u(t)). Then we can choose p 0 = 0 and the Hamiltonian Ĥ evaluated along (x(•), p(•), u(•)) doesn't depend on the cost L(x, u). This case is called the abnormal case.

Second order conditions and conjugate points

In this section we make a brief introduction to the concept of conjugate point in optimal control, in relation with second order conditions, generalizing the similar concepts in calculus of variations presented in section 2.5.7.

The underlying geometric framework is elegant and corresponds to the concept of Lagrangian manifold [START_REF] Mishchenko | Lagrangian manifolds and the Maslov operator[END_REF] and singularity of projection of Lagrangian manifold [START_REF] Arnol˜d | Singularities of differentiable maps[END_REF][START_REF] Zakaljukin | Lagrangian and Legendre singularities[END_REF]. They can be numerically computed using rank tests on Jacobi fields which is one of the key components of the HamPath code [START_REF] Cots | Contrôle optimal géométrique: méthodes homotopiques et applications[END_REF]. Also this concept is well known to be related to the zero eigenvalue of self-adjoint operators associated to the intrinsic second order derivative [START_REF] Klingenberg | Riemannian geometry[END_REF].

Lagrangian manifold and Jacobi equation

Definition 26 Let (M, ω) be a (smooth) symplectic manifold of dimension 2n. A regular submanifold L of M of dimension n is called Lagrangian if the restriction of ω to T x L × T x L is zero.
Definition 27 Let L be a Lagrangian submanifold of T * M and let Π : z = (x, p) → x be the canonical projection. A tangent non zero vector v of L is called vertical if dΠ (v) = 0. We call caustic the set of points x of L such that there exists at least one vertical field.

Definition 28 Let

-→ H be a (smooth) Hamiltonian vector field on T * M, ϕ t = expt -→ H the associated one parameter group, L 0 the fiber T x M and L t = ϕ t (L 0 ). The set of caustics is called the set of conjugate loci of L.

Definition 29 Let

-→ H be a (smooth) Hamiltonian vector field on T * M and let z(t) = (x(t), p(t)) be a reference trajectory of

- → H defined on [0, T ]. The variational equation δ z(t) = ∂ - → H ∂ z (z(t))δ z(t)
is called Jacobi equation. We called Jacobi field J(t) = (δ x(t), δ p(t)) a non trivial solution of Jacobi equation. It is said to be vertical at time t if δ x(t) = 0. A time t c is called conjugated if there exists a Jacobi field vertical at times 0 and t c and the point x(t c ) is called geometrically conjugate to x(0).

Numerical computation of the conjugate loci along a reference trajectory

Verticality test

Let z(t) = (x(t), p(t)) be a reference trajectory of -→ H and x 0 = x(0). The set of Jacobi fields forms an n-dimensional linear subspace. Let (e 1 , ..., e n ) be a basis of T *

x 0 M and let J i (t) = (δ x i (t), δ p i (t)), i = 1, ..., n the set of Jacobi fields (vertical at t = 0), such that δ x i (0) = 0, δ p i (0) = e i . Therefore the time t c is geometrically conjugate if and only if the rank of dΠ

z(t c ) (J 1 (t c ), ..., J n (t c ))
is strictly less than n.

Sub-Riemannian Geometry

In this section a quick introduction to sub-Riemannian (SR-geometry) is presented which is the proper geometry framework for the swimming problem at low Reynolds number.

Sub-Riemannian manifold

Definition 30 A sub-Riemannian manifold is a triple (M, D, g) where M is a smooth connected manifold, D is a smooth distribution of rank m on M and g is a riemannian metric on M.

An horizontal curve is an absolutely continuous curve t → γ(t), t ∈ I such that γ(t) ∈ D(γ(t)). The length of a curve γ is defined by l(γ) = I g( γ(t)) 1/2 dt and its energy is given by E(γ) = 1/2 T 0 g( γ(t)) dt where the final time T can be fixed at 1.

Controllability

Let D 1 = D, D k = D 1 + [D 1 , D k-1 ].
We assume that there exists for each x ∈ M an integer r(x), called the degree of non holonomy, such that D r(x) = T x M. Moreover at a point x ∈ M, the distribution D is characterized by the growth vector (n 1 , n 2 , ..., n r )

where n k = dim D k (x).

Distance

According to Chow's theorem, for each pair (x, y) ∈ M, there exists an horizontal curve γ : [0, 1] → M such that γ(0) = x, γ(1) = y. We denote by d the sub-Riemannian distance (SR-distance):

d(x, y) = inf γ {l(γ)
; γ is an horizontal curve joining x to y}.

Geodesics equations

According to Maupertuis principle the length minimization problem is equivalent to the energy minimization problem. Additionally if we parametrize the curves by arclength, then the length minimization problem is equivalent to the time minimization problem.

To compute the geodesics equations it is convenient to minimize the energy E. We proceed for the calculations as follows. We choose a local orthonormal frame {F 1 , ..., F m } of D, and we consider the minimization problem:

dx dt (t) = m ∑ i=1 u i (t)F i (x(t)), min u(.) 1 2 1 0 ∑ i u 2 i (t) dt.
According to the weak maximum principle (corresponding to a control domain U = R m ) we introduce the pseudo-Hamiltonian:

H(x, p, u) = m ∑ i=1 u i H i (x, p) + p 0 m ∑ i=1 u 2 i
where H i (x, p) = p, F i (x) is the Hamiltonian lift of F i . By homogeneity p 0 can be normalized to 0 or -1 2 .

Normal case:

p 0 = -1/2.
According to the maximum principle the condition ∂ H ∂ u = 0 leads to u i = H i . Plugging this last expression for u i into H leads to the true Hamiltonian in the normal case:

H n (z) = 1 2 m ∑ i=1 H 2 i (z)
where z = (x, p). A normal extremal is a solution of the Hamiltonian system associated to the true Hamiltonian, and its projection on the state space is called a normal geodesic.

Abnormal case: p 0 = 0.

In this case, the maximum principle leads to the conditions H i = 0, i = 1, ..., m, thus defining implicitly the abnormal curves related to the structure of the distribution D.

The solutions are called abnormal extremals, and their projections on the state space are the abnormal geodesics.

Next we introduce the basic definitions related to the analysis of the geodesics equations and generalizing the Riemannian concepts. Definition 32 Let x ∈ M be fixed. The set of points at a SR-distance less or equal to r from x form the ball of radius r centered at x and the sphere S(x, r) is formed by the set of points at a distance r from x.

Evaluation of the Sub-Riemannian ball

The computation of the Sub-Riemannian ball (SR-ball), even with small radius is a very complicated task. One of the most important result in SR-geometry is an approximation result about balls of small radius, in relation with the structure of the distribution.

Definition 33 Let x ∈ M, and let f be a germ of a smooth function at x. The multiplicity of f at x is the number µ( f ) defined by:

• µ( f ) = min{n; there exist X 1 , ..., X n ∈ D(x) such that: (L X 1 • .... • L X n f )(x) = 0}, • if f (x) = 0 then µ( f ) = 0, and µ(0) = +∞. Definition 34 Let f be a germ of a smooth function at x, f is called privileged at x if we have that µ( f ) is equivalent to min{k; d f x (D k (x)) = 0}. A coordinate system {x 1 , ..., x n } : V → R defined on an open subset V of x is called privileged if all the coordinates functions x i , 1 ≤ i ≤ n are privileged at x.

Nilpotent Approximation

Let us fix a privileged coordinate system at x = (x 1 , ..., x n ), where the weight of x i is given by µ(x i ). Each smooth vector field V at x has a formal expansion V ∼ ∑ j≥-1 V j , where each V j = ∑ n i=1 P j i (x 1 , ..., x n ) ∂ ∂ x i is homogeneous of degree j for the weights associated with the coordinate system, and the weight of ∂ ∂ x i is -µ(x i ). P j i (x 1 , ..., x n ) is a homogenous polynomial of degree j.

Proposition 13 Let {F 1 , ..., F m } be the orthonormal subframe of the distribution D and set Fi = F -1 i , i = 1, ..., m in the formal expansion. Then, the family Fi is a first order approximation of {F 1 , ..., F m } at x since they generate a nilpotent Lie algebra with a similar growth vector. Moreover, for small x it gives the following estimate of the SR-

norm |x| = d(0, x) |x 1 | 1/w 1 + ...|x n | 1/w n .
See [START_REF] Bellaïche | The tangent space in sub-Riemannian geometry[END_REF], [START_REF] Kupka | Géométrie sous-riemannienne[END_REF] and [START_REF] Jean | Control of Nonholonomic Systems: from Sub-Riemannian Geometry to Motion Planning[END_REF] for the details of the construction of privileged coordinates. In addition, note that [START_REF] Maciejewski | The nilpotent tangent 3-dimensional sub-Riemannian problem is nonintegrable[END_REF] contains also the relation of the integrability issues which is important for the practical implementation.

Conjugate and cut loci in SR-geometry

The standard concepts of conjugate and cut point from Riemannian geometry can be generalized in optimal control and thus in SR-geometry. Consider the SR-problem:

ẋ(t) = m ∑ i=1 u i (t)F i (x(t)), min u(.) T 0 m ∑ i=1 u 2 i (t) 1/2 dt.
Definition 35 Let x(.) be a reference (normal or abnormal) geodesic defined on [0, T ]. The time t c is called the cut time if the reference geodesic stops to be optimal at t = t c , i.e. is no longer optimal for t > t c , and x(t c ) is called the cut point. Taking all geodesics starting from x 0 = x(0), their cut points will form the cut locus C cut (x 0 ). The time t 1c is called the first conjugate time if it is the first time such that the reference geodesic is no longer optimal for t > t 1c for the C 1 -topology on the set of curves, and the point x(t 1c ) is called the first conjugate point. Calculated over all geodesics, the set of first conjugate points will form the (first) conjugate locus C(x 0 ).

An important step is to relate the computation of the geometric conjugate locus (using a test based on Jacobi fields) to the computation of the conjugate locus associated to optimality. It can be done under suitable assumptions in both the normal and the abnormal case [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] but for simplicity we shall restrict ourselves to the normal case.

Conjugate locus computation

Using Maupertuis principle, the SR-problem is equivalent to the (parametrized) energy minimization problem:

min u(.) T 0 m ∑ i=1 u 2 i (t) dt
where T is fixed, and we can choose

T = 1. Let H i (z) = p, F i (x) and let H n (z) = 1 2 ∑ m i=1 H 2 i (z)
be the Hamiltonian in the normal case. Take a reference normal geodesic x(.) defined on [0, 1] and let z(.) = (x(.), p(.)) be a symplectic lift solution of -→ H n . Moreover assume that x(.) is strict, which means that it is not a projection of an abnormal curve. Then the following proposition holds.

Proposition 14

The first conjugate time t 1c along x(.) corresponds to the first geometric conjugate point and can be computed numerically using the test of Section 2.4.

Integrable case

If the geodesic flow is Liouville integrable, then the Jacobi equation is integrable and the conjugate points can be computed using the parametrization of the geodesic curve. This result is a consequence of the following standard lemma from differential geometry.

Lemma 2 Let J(t) = (δ x(t), δ p(t)) be a Jacobi curve along z(t) = (x(t), p(t)), t ∈ [0, 1] and vertical at time t = 0, i.e. δ x(0) = 0. Let α(ε) be any curve in T *

x 0 M defined by p(0) + εδ p(0) + o(ε). Then:

J(t) = d dε |ε=0 exptH n (x(0), α(ε)).

Nilpotent models in relation with the swimming problem

The models in dimension 3 are related to the classification of stable 2-dimensional distributions, see [START_REF] Zhitomirskiȋ | Typical singularities of differential 1-forms and Pfaffian equations[END_REF], and will be used for the copepod swimmer. See also [START_REF] Brockett | Control theory and singular Riemannian geometry[END_REF] for the analysis of the Heisenberg case.

Contact case. A point x 0 ∈ R 3 is a contact point of the distribution D = span{F 1 , F 2 } if [F 1 , F 2 ](x 0 ) / ∈ D(x 0 )
and the growth vector is (2, 3). A normal form at x 0 ∼ 0 is given by:

x = (x 1 , x 2 , x 3 ), D = ker α, α = x 2 dx 1 + dx 3 .
Observe that

• dα = dx 2 ∧ dx 1 : Darboux form, • ∂ ∂ x 3 is equal to the Lie bracket [F 1 , F 2 ]
and is the characteristic direction of dα.

This form is equivalent to the so-called Dido representation:

D = ker α , α = dx 3 + (x 1 dx 2 -x 2 dx 1 ) with D = span{F 1 , F 2 }, F 1 = ∂ ∂ x 1 + x 2 ∂ ∂ x 3 , F 2 = ∂ ∂ x 2 -x 1 ∂ ∂ x 3 . If we set F 3 = ∂ ∂ x 3
, we have that [F 1 , F 2 ] = 2 F 3 and the corresponding so-called Heisenberg SR-case is given by:

ẋ(t) = 2 ∑ i=1 u i (t)F i (x(t)), min u(.) T 0 (u 2 1 (t) + u 2 2 (t)) dt.
It corresponds to minimizing the Euclidean length of the projection of the curve t → x(t) on the (x 1 , x 2 )-plane. Starting from the origin (0, 0, 0), we observe that

x 3 (T ) = T 0 ( ẋ1 (t)x 2 (t) -ẋ2 (t)x 1 (t)) dt
is proportional to the area swept by the curve t → (x 1 (t), x 2 (t)). The Heisenberg SR-case is therefore dual to the Dido problem: among the closed curves in the plane with fixed length, find those for which the enclosed area is maximal. The solutions are well known and they are arcs of circles. They can be easily obtained using simple computations as follows. The geodesic equations written in the (x, H) coordinates where

H = (H 1 , H 2 , H 3 ), H i = p, F i , i = 1, 2, 3 are given by: ẋ1 = H 1 , ẋ2 = H 2 , ẋ3 = H 1 x 2 -H 2 x 1 , Ḣ1 = 2H 2 H 3 , Ḣ2 = -2H 1 H 3 , Ḣ3 = 0.
Since H 3 is constant we can introduce H 3 = λ /2 with λ ∈ R, and we obtain the equation of a linear pendulum Ḧ1 + λ 2 H 1 = 0. The integration can be done directly since we can observe that:

ẍ3 - λ 2 d dt (x 2 1 + x 2 2 ) = 0.
Since λ = 0, which can be assumed positive, we obtain the well known parametrization for the geodesics:

x 1 (t) = A λ (sin(λt + ϕ) -sin(ϕ)) x 2 (t) = A λ (cos(λt + ϕ) -cos(ϕ)) x 3 (t) = A 2 λ t - A 2 λ 2 sin(λt) with A = H 2 1 + H 2 2
and ϕ is the angle of the vector ( ẋ1 , -ẋ2 ). If λ = 0, the geodesics are straight lines.

Conjugate points. Computations of first conjugate points are straightforward using the parameterization above for the normal geodesics. Only geodesics whose projections are circles have a first conjugate point given by t c = 2π/λ which corresponds to the first intersection of the geodesic with the axis Ox 3 . Geometrically, it is due to the symmetry of revolution along this axis which produces a one-parameter family of geodesics starting from 0 and intersecting at such point. This point is also a cut point and a geodesic is optimal up to this point (included).

Note that the SR-Heisenberg case will lead to interesting geometric consequences in the swimming problem: the circles projections correspond to the concept of stroke. But while this model can provide some insights on optimal swimming, it is too primitive because:

1. The geodesic flow is integrable due to the symmetries and every (x 1 , x 2 ) motion is periodic; 2. The model is quasi-homogeneous where x 1 , x 2 are of weight 1 and x 3 is of weight 2.

Martinet case. A point x 0 is a Martinet point if at x 0 , [F 1 , F 2 ] ∈ span{F 1 , F 2 } and at least one Lie bracket [[F 1 , F 2 ], F 1 ] or [[F 1 , F 2 ], F 2 ]
does not belong to D . Hence the growth vector is (2, 2, 3). Then, there exist local coordinates near x 0 identified to the origin such that:

D = ker ω, ω = dx 3 - x 2 2 2 dx 1
where

F 1 = ∂ ∂ x 1 + x 2 2 2 ∂ ∂ x 3 , F 2 = ∂ ∂ x 2 , F 3 = [F 1 , F 2 ] = x 2 ∂ ∂ x 3 .
The surface Σ : det(F 1 , F 2 , [F 1 , F 2 ]) = 0 is identified to x 2 = 0 and is called the Martinet surface. This surface is foliated by abnormal curves which are integral curves of ∂ ∂ x 1 . In particular abnormal curves passing through the origin and parameterized by arc-length corresponds to the curve t → (t, 0, 0).

Those two cases are nilpotent Lie algebras associated to nilpotent approximations of the SR-metric in the copepod swimmer and are respectively the Heisenberg case and the Martinet flat case. Also it can be easily checked that this second case leads to integrable geodesic flow using elliptic functions.

Swimming problems at low Reynolds number

2.6.1 Purcell's 3-link swimmer.

The 3-link swimmer is modeled by the position of the center of the second stick x = (x, y) as well as the angle α between the x-axis and the second stick (the orientation of the swimmer). The shape of the swimmer is modeled by the two relative angles θ 1 and θ 2 (see Fig 2.1). We also denote respectively by L and L 2 the length of the two external arms and central link. In what follows, x (resp. x ) corresponds to (x, y) (resp. to (α, θ 1 , θ 2 )).

Dynamics via Resistive Force Theory.

We approximate the non local hydrodynamic forces exerted by the fluid on the swimmer with local drag forces depending linearly on the velocity. For each i ∈ {1, 2, 3}, we denote by e i and e ⊥ i the unit vectors parallel and perpendicular to the i-th link, and we also introduce v i (s) the velocity of the point at distance s from the extremity of the i-th link, that is:

v 1 (s) = ẋ - L 2 2 αe ⊥ 2 -s( α -θ1 )e ⊥ 1 , s ∈ [0, L], v 2 (s) = ẋ -(s - L 2 2 ) αe ⊥ 2 , s ∈ [0, L 2 ], v 3 (s) = ẋ + L 2 2 αe ⊥ 2 + s( α -θ2 )e ⊥ 3 , s ∈ [0, L].
The force f i acting on the i-th segment is taken as:

f i (s) := -c t v i (s) • e i e i -c n v i (s) • e ⊥ i e ⊥ i
where c t and c n are respectively the drag coefficients in the directions of e i and e ⊥ i . Neglecting inertia forces, Newton laws are written as:

f = 0 , e z • T x = 0 (2.12)
where f is the total force exerted on the swimmer by the fluid and e z = e x ∧ e y ,

f = L 0 f 1 (s) ds + L 2 0 f 2 (s) ds + L 0 f 3 (s) ds
and T x is the corresponding total torque computed with respect to the central point x,

T x = L 0 (x 1 (s) -x 1 ) × f 1 (s) ds + L 2 0 (x 2 (s) -x 1 ) × f 2 (s) ds + L 0 (x 3 (s) -x 1 ) × f 3 (s) ds
where x i = (x i , y i ), for i = 1, 2, 3, corresponds to the left-end point of the i-th link, and x i (s) = x i + se i .

Since the f i (s) are linear in ẋ, α, θ1 , θ2 , the system (2.12) can be rewritten as

A(q) • ẋ α -B(q) • θ1 θ2 = 0
where q(t) = (θ 1 , θ 2 , x, y, α)(t). The matrix A(q) is invertible (see [START_REF] Alouges | Self-propulsion of slender microswimmers by curvature control: N-link swimmers[END_REF]). Then, the dynamics of the swimmer is finally expressed as the system

q(t) = f (q, θ1 , θ2 ) = θ1 (t) F 1 (q(t)) + θ2 (t) F 2 (q(t))
where

F 1 (q) F 2 (q) := I 2 A -1 (q)B(q)
with I 2 the 2 × 2 identity matrix. The equations of the dynamics take the form

  ẋ ẏ α  = 1 G R α   g 11 g 12 g 21 g 22 g 31 g 32   θ1 θ2 , θ = u = S(θ )τ (2.13)
where τ is the torque, R α is the rotation matrix R α =   cos(α) -sin(α) 0 sin(α) cos(α) 0 0 0 1   and g i j , G and S are functions depending only on (θ 1 , θ 2 ) which have long expressions (cf. [START_REF] Passov | Supplementary notes to: Dynamics of Purcells three-link microswimmer with a passive elastic tail[END_REF] for a details).

The cost function u is minimizing the expanded mechanical power

T 0 τ • u dt (2.14)
where τu = uH -1 u and H -1 (θ ) is the symmetric matrix described in [START_REF] Passov | Supplementary notes to: Dynamics of Purcells three-link microswimmer with a passive elastic tail[END_REF]. It can be computed as

T 0 L 0 f 1 • v 1 + L 2 0 f 2 • v 2 + L 0 f 3 • v 3 .
Expressions of the controlled vector fields and the mechanical energy.

Normalizing L = 2, L 2 = 1, c t = 1, c n = 2,
we write the swimming control system (2.13) as

q(t) = 2 ∑ i=1 u i (t)F i (q(t)), (2.15)
and we obtain the following expressions of the vector fields F 1 , F 2 :

F 1 = ∂ ∂ θ 1 + ∂ ∂ x f 13 + ∂ ∂ y f 14 + ∂ ∂ α f 15 , (2.16
)

F 2 = ∂ ∂ θ 2 + ∂ ∂ x f 23 + ∂ ∂ y f 24 + ∂ ∂ α f 25 (2.17) where δ = 1692 + 336 cos (θ 1 -θ 2 ) + 84 cos (2 θ 1 ) -24 cos (θ 1 + 2 θ 2 ) -48 cos (θ 1 + θ 2 ) + 816 cos (θ 2 ) + 72 cos (-2 θ 2 + θ 1 ) + 816 cos (θ 1 ) -6 cos (2 θ 1 + 2 θ 2 ) + 18 cos (-2 θ 2 + 2 θ 1 ) + 84 cos (2 θ 2 ) -24 cos (2 θ 1 + θ 2 ) + 72 cos (-θ 2 + 2 θ 1 ) in . f 13 = 1/δ 4 sin (α -2 θ 2 ) -sin (α + 2 θ 2 -θ 1 ) + 18 sin (α -θ 1 -θ 2 ) + 3 sin (α -θ 1 -2 θ 2 ) + 2 sin (α -2 θ 1 + 2 θ 2 ) -9 sin (α + θ 1 -2 θ 2 ) -21 sin (α + θ 1 + 2 θ 2 ) -126 sin (α + θ 1 + θ 2 ) -30 sin (α -θ 1 + θ 2 ) -2 sin (α + 2 θ 1 -2 θ 2 ) + 2 sin (α -2 θ 1 ) -78 sin (α + θ 1 -θ 2 ) + 16 sin (α -θ 2 ) -104 sin (α + θ 2 ) -8 sin (α + 2 θ 1 -θ 2 ) -24 sin (α + 2 θ 2 ) -18 sin (α + 2 θ 1 ) -36 sin (α) -262 sin (α + θ 1 ) + 26 sin (α -θ 1 ) , . f 14 = 1/δ 18 cos (α + 2 θ 1 ) + 24 cos (α + 2 θ 2 ) + 30 cos (α -θ 1 + θ 2 ) -3 cos (α -θ 1 -2 θ 2 ) + 126 cos (α + θ 1 + θ 2 ) + 78 cos (α + θ 1 -θ 2 ) -18 cos (α -θ 1 -θ 2 ) + 21 cos (α + θ 1 + 2 θ 2 ) + 9 cos (α + θ 1 -2 θ 2 ) -26 cos (α -θ 1 ) + 104 cos (α + θ 2 ) -16 cos (α -θ 2 ) + 8 cos (α + 2 θ 1 -θ 2 ) -4 cos (α -2 θ 2 ) + 36 cos (α) + 262 cos (α + θ 1 ) + cos (α + 2 θ 2 -θ 1 ) -2 cos (α -2 θ 1 ) -2 cos (α -2 θ 1 + 2 θ 2 ) + 2 cos (α + 2 θ 1 -2 θ 2 ) , . f 15 = 1/δ -216 -4 cos (2 θ 1 ) + 6 cos (θ 1 + 2 θ 2 ) + 12 cos (θ 1 + θ 2 ) -204 cos (θ 1 ) -18 cos (-2 θ 2 + θ 1 ) -84 cos (θ 1 -θ 2 ) -4 cos (-2 θ 2 + 2 θ 1 ) + 8 cos (2 θ 2 ) , . f 23 = 1/δ 21 sin (α + θ 2 + 2 θ 1 ) -2 sin (α + 2 θ 1 -2 θ 2 ) -2 sin (α -2 θ 2 ) + 9 sin (α + θ 2 -2 θ 1 ) + 2 sin (α -2 θ 1 + 2 θ 2 ) + 30 sin (α + θ 1 -θ 2 ) + 8 sin (α + 2 θ 2 -θ 1 ) -3 sin (α -θ 2 -2 θ 1 ) -18 sin (α -θ 1 -θ 2 ) + 126 sin (α + θ 1 + θ 2 ) + 78 sin (α -θ 1 + θ 2 ) + sin (α + 2 θ 1 -θ 2 ) + 262 sin (α + θ 2 ) + 104 sin (α + θ 1 ) -4 sin (α -2 θ 1 ) -16 sin (α -θ 1 ) -26 sin (α -θ 2 ) + 24 sin (α + 2 θ 1 ) + 18 sin (α + 2 θ 2 ) + 36 sin (α) , . f 24 = 1/δ 4 cos (α -2 θ 1 ) -2 cos (α -2 θ 1 + 2 θ 2 ) -8 cos (α + 2 θ 2 -θ 1 ) + 2 cos (α -2 θ 2 ) -18 cos (α + 2 θ 2 ) + 26 cos (α -θ 2 ) -24 cos (α + 2 θ 1 ) -cos (α + 2 θ 1 -θ 2 ) + 2 cos (α + 2 θ 1 -2 θ 2 ) -30 cos (α + θ 1 -θ 2 ) -21 cos (α + θ 2 + 2 θ 1 ) -126 cos (α + θ 1 + θ 2 ) -78 cos (α -θ 1 + θ 2 ) + 3 cos (α -θ 2 -2 θ 1 ) -9 cos (α + θ 2 -2 θ 1 ) + 18 cos (α -θ 1 -θ 2 ) + 16 cos (α -θ 1 ) -104 cos (α + θ 1 ) -262 cos (α + θ 2 ) -36 cos (α) , . f 25 = 1/δ -2168 cos (2 θ 1 ) + 12 cos (θ 1 + θ 2 ) + 6 cos (2 θ 1 + θ 2 ) -4 cos (2 θ 2 ) -18 cos (2 θ 1 -θ 2 ) -204 cos (θ 2 ) -4 cos (-2 θ 2 + 2 θ 1 ) -84 cos (θ 1 -θ 2 ) .
Moreover, writing the integrand of the cost function (2.14) as au 2 1 + 2bu 1 u 2 + cu 2 2 , the coefficients a, b, c are given by -a(q) = 1/κ 3 cos

(2 θ 1 + 2 θ 2 ) -6 cos (-2 θ 2 + 2 θ 1 ) -12 cos (2 θ 1 -θ 2 ) + 24 cos (2 θ 1 + θ 2 ) + 72 cos (2 θ 1 ) -84 cos (2 θ 2 ) -492 cos (θ 2 ) -1233 , -b(q) = 1/κ cos (2 θ 1 + 2 θ 2 )-246 cos (θ 1 )-246 cos (θ 2 )+12 cos (2 θ 1 + θ 2 )- 6 cos (2 θ 1 -θ 2 ) + 12 cos (θ 1 + 2 θ 2 ) + 84 cos (θ 1 + θ 2 ) -276 cos (θ 1 -θ 2 ) -6 cos (-2 θ 2 + θ 1 ) -4 cos (2 θ 2 ) -4 cos (2 θ 1 ) -153 , -c(q) = 1/κ 3 cos (2 θ 1 + 2 θ 2 ) -492 cos (θ 1 ) -6 cos (-2 θ 2 + 2 θ 1 ) +24 cos (θ 1 + 2 θ 2 )-12 cos (-2 θ 2 + θ 1)+72 cos (2 θ 2 )-84 cos (2 θ 1 )-1233 .
where κ = 36 cos

(θ 1 -2 θ 2 ) -222 cos (2 θ 1 ) -1116 cos (θ 2 ) -222 cos (2 θ 2 ) + 18 cos (-2 θ 2 + 2 θ 1 ) -72 cos (2 θ 1 + θ 2 ) -72 cos (θ 1 + 2 θ 2 ) -180 cos (θ 1 + θ 2 ) +36 cos (2 θ 1 -θ 2 )-1116 cos (θ 1 )+36 cos (θ 1 -θ 2 )-12 cos (2 θ 1 + 2 θ 2 )-3258.

Copepod swimmer

It is a simplified model proposed by [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF] of a symmetric swimming where only line displacement is authorized, see also [START_REF] Avron | A geometric theory of swimming: Purcell's swimmer and its symmetrized cousin[END_REF]. It consists in two pairs of symmetric links of equal lengths with respective angles θ 1 , θ 2 with respect to the displacement directions Ox while the body is assumed to be an infinitesimal sphere, see Fig. 2.2. The swimming velocity at x 0 is given by

x 0 (t) x ✓ 1 (t) ✓ 2 (t) ✓ = ⇡ ✓ = 0 ✓ 2 (t) ✓ 1 (t)
ẋ0 = θ1 sin θ 1 + θ2 sin θ 2 2 + sin 2 θ 1 + sin 2 θ 2 (2.18) and θ1 = u 1 , θ2 = u 2 .
The mechanical energy is the quadratic form q M qt where q = (x 0 , θ 1 , θ 2 ) is the state variable and M is the symmetric matrix

M =   2 -1/2 (cos 2 θ 1 + cos 2 θ 2 ) -1/2 sin θ 1 -1/2 sin θ 2 -1/2 sin θ 1 1/3 0 -1/2 sin θ 2 0 1/3   .
The corresponding Riemannian metric defines the associated SR-metric thanks to the relation between ẋ0 and θ1 , θ2 .

Some geometric remarks

In order to analyze the swimming problem one must introduce the concept of stroke.

Definition 36 A stroke is a periodic motion of the shape variables associated with a periodic control producing a net displacement of the displacement variable after one period. Observe that due to the SR-structure one can fix the period of the stroke to 2π.

A first geometric analysis is to consider bang-bang controls and the associated strokes. For a single link one gets the famous scallop theorem.

Theorem 7 A scallop cannot swim.

Proof. The relation between the displacement and angular velocity is given by the relation

ẋ0 = sin(θ ) θ 2 -cos 2 (θ ) , θ = u
where θ is the angle of the symmetric link with respect to the axis. Let γ be the angle with respect to the vertical and a stroke is given by

u = 1 : θ : π/2 -γ → π/2 u = -1 : θ : π/2 → π/2 -γ
and the control u = 1 produces a displacement: x 0 → x 1 while the control u = -1 reverses the motion: x 1 → x 0 . The net displacement of the stroke is zero and clearly is related to the reversibility of the SR-model.

A similar computation can be done on the Purcell swimmer using a square stroke like in the original paper ( [START_REF] Purcell | Life at low Reynolds number[END_REF]). Considering the controlled system (2.15), the displacement associated with the sequence stroke described in Fig. 2.3 is given by

β (t) = (exptF 2 exp -tF 1 exp -tF 2 exptF 1 ) (q(0)), q = (θ 1 , θ 2 , x, y, α),
and using Baker-Campbell-Hausdorff formula one has

β (t) = exp(t 2 [F 1 , F 2 ] + o(t 2 ))(q(0))
which gives for small stroke t a displacement of

β (t) ∼ q(0) + t 2 [F 1 , F 2 ](q(0)).
This shall be compared with [START_REF] Becker | On self-propulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer[END_REF]. Hence for a small square stroke the displacement can be evaluated using (2.16), (2.17).

In the case of the copepod swimmer, due to the constraints θ i ∈ [0, π], θ 1 ≤ θ 2 on the shape variable, a geometric stroke corresponds to a triangle in the shape variable and is defined by θ 2 : 0 → π; θ 1 : 0 → π and θ 1 = θ 2 : π → 0. See in the specific analysis of the copepod swimmer the interpretation of this stroke (see Fig. 2.20 (right)).

Purcell swimmer

Due to the mathematical complexity of the expressions of the vector fields F 1 and F 2 (cf. Section 2.6.1) employed in this model, the nilpotent approximation will play a crucial role in our analysis. First, as a consequence of the integrability of the associated normal extremals in the class of elliptic functions, the nilpotent approximation will allow us to make a micro-local analysis of the different kinds of strokes and to establish the existence of conjugate points using a suitable time rescaling. Second, the abnormal extremals forming piecewise smooth strokes can be easily computed in this approximation and, then, the optimality of these strokes can be studied using the concept of the (corresponding) conjugate point.

The flat nilpotent model

The Purcell system (2.13) can be written as a control system of the form q = F(q)u = ∑ 2 i=1 u i F i (q), where q = (θ 1 , θ 2 , x, y, α) ∈ R 5 . Even though the vectors fields F 1 , F 2 have a complicated expression, they provide a 2-distribution with growth (2, 3, 5) (see [START_REF] Bettiol | The three links Purcell swimmer and some geometric problems related to periodic optimal controls[END_REF]). There exists a unique nilpotent model associated with a 2-dimensional distribution in dimension 5 with growth vector (2, 3, 5), see [START_REF] Cartan | Les systèmes de Pfaff a cinq variables et les équations aux derivées partielles du second ordre[END_REF][START_REF] Sachkov | Symmetries of flat rank two distributions and sub-Riemannian structures[END_REF].

Definition 2.1. We call the flat Cartan model the 2-dimensional distribution in dimension five defined by the two vector fields:

F1 ( x) = ∂ ∂ x1 , F2 ( x) = ∂ ∂ x2 + x1 ∂ ∂ x3 + x3 ∂ ∂ x4 + x2 1 ∂ ∂ x5 (2.19) where x = (x 1 , x2 , x3 , x4 , x5
) are the privileged coordinates with the following weights: 1 for x1 and x2 , 2 for x3 , and 3 for x4 and x5 .

Computations of the nilpotent approximation

The nilpotent approximation of the Purcell model is computed at the origin. It provides a nilpotent approximation for the SR-problem with the simplified cost 2π 0 (u 2 1 (t) + u 2 2 (t)) dt. The two-jets of F 1 and F 2 at q = (0, 0, 0, 0, 0) are given by:

F 1 (q) = ∂ ∂ q 1 + - 1 6 q 5 - 4 27 q 1 - 2 27 q 2 ∂ ∂ q 3 + 1 6 - 1 12 q 5 2 - 2 27 q 5 q 2 - 4 27 q 5 q 1 - 1 27 q 1 2 - 1 27 q 1 q 2 - 1 36 q 2 2 ∂ ∂ q 4 + - 7 27 + 2 81 q 1 2 - 2 81 q 1 q 2 - 5 162 q 2 2 ∂ ∂ q 5 + O(|q| 3 ) F 2 (q) = ∂ ∂ q 2 + 1 6 q 5 + 4 27 q 2 + 2 27 q 1 ∂ ∂ q 3 + - 1 6 + 1 12 q 5 2 + 4 27 q 5 q 2 + 2 27 q 5 q 1 + 1 36 q 1 2 + 1 27 q 1 q 2 + 1 27 q 2 2 ∂ ∂ q 4 + - 7 27 - 5 162 q 1 2 - 2 81 q 1 q 2 + 2 81 q 2 2 ∂ ∂ q 5 + O(|q| 3 ) .
The local diffeomorphism ϕ, which transforms F 1 , F 2 into the nilpotent approximation F1 , F2 , can be explicitly written using a sequence ϕ = ϕ N o ... o ϕ 1 : R 5 → R 5 , where N = 13 (see [START_REF] Bettiol | The three links Purcell swimmer and some geometric problems related to periodic optimal controls[END_REF]). This leads to a complicated transformation whose role is to relate the privileged coordinates to the physical coordinates (θ 1 , θ 2 , x, y, α) in particular we have a 'stability' property for the shape variables as stated in the next proposition.

Proposition 2.1. The shape variables θ = (θ 1 , θ 2 ) corresponds to the ( x1 , x2 ) coordinates.

Integration of normal extremal trajectories

Computing with (2.19), we obtain:

F1 ( x) = ∂ ∂ x1 , F2 ( x) = ∂ ∂ x2 + x1 ∂ ∂ x3 + x3 ∂ ∂ x4 + x2 1 ∂ ∂ x5 , [ F1 , F2 ]( x) = - ∂ ∂ x3 -2 x1 ∂ ∂ x5 , [[ F1 , F2 ], F1 ]( x) = -2 ∂ ∂ x5 , [[ F1 , F2 ], F2 ]( x) = ∂ ∂ x4 .
All brackets of length greater than 3 are zero. Let us introduce ẑ = ( x, p). Employing the corresponding Hamiltonian lifts, we have:

H 1 (ẑ) = p, F1 ( x) = p1 , H 2 (ẑ) = p, F2 ( x) = p2 + p3 x1 + p4 x3 + p5 x2 1 , H 3 (ẑ) = p, [ F1 , F2 ]( x) = -p3 -2 x1 p5 , H 4 (ẑ) = p, [[ F1 , F2 ], F1 ]( x) = -2 p5 , H 5 (ẑ) = p, [[ F1 , F2 ], F2 ]( x) = p4 . The SR-Cartan flat case is ẋ(t) = 2 ∑ i=1 u i (t) Fi ( x(t)), min u 2π 0 (u 2 1 (t) + u 2 2 (t)) dt.
and the normal Hamiltonian takes the form

H n = 1/2 (H 2 1 + H 2 2 ). (2.20)
More precisely, using the Poincaré coordinates, the control system can be written as:

ẋ1 = H 1 , ẋ2 = H 2 , ẋ3 = H 2 x1 , ẋ4 = H 2 x3 , ẋ5 = H 2 x2 1 .
(2.21)

By differentiating with respect to the time variable, we obtain:

Ḣ1 = dH 1 ( -→ H n ) = {H 1 , H 2 }H 2 = p, [ F1 , F2 ]( x) H 2 = H 3 H 2 , Ḣ2 = -H 3 H 1 , Ḣ3 = H 1 H 4 + H 2 H 5 , Ḣ4 = 0 hence H 4 = c 4 , Ḣ5 = 0 hence H 5 = c 5 .
We fix the energy level H 2 1 + H 2 2 to 1, and we introduce H 1 = cos ϑ and H 2 = sin ϑ which implies:

Ḣ1 = -sin ϑ θ = H 2 H 3 = sin ϑ H 3 . It follows that θ = -H 3 and θ = -(H 1 c 4 + H 2 c 5 ) = -c 4 cos ϑ -c 5 sin ϑ = -ω 2 sin(ϑ + φ ) (2.22)
where ω and φ are constant. More precisely, we have:

ω = p2 40 + 4 p2 50 1/4 , φ = arctan (-2 p50 / p40 ) .
• First, we consider the degenerate case which corresponds to θ = 0. Therefore, ϑ (t) = p30 t + ϑ 0 where ϑ 0 , p30 are constant and for p30 = 0, the solutions x1 , x2 of (2.21) are expressed as:

x1 (t) = x10 + 1/ p30 sin( p30 t + ϑ 0 ), x2 (t) = x20 -1/ p30 cos( p30 t + ϑ 0 ) (2.23)
where x10 , x20 are constant. • Second, the case corresponding to θ = 0 leads to a pendulum equation. Indeed, by introducing ψ = ϑ + φ , (2.22) becomes:

1/2 ψ2 -ω 2 cos(ψ) = B, (2.24)
where B is the constant

B = 1/2 ( p30 + 2 x10 p50 ) 2 -p10 p40 -2 p50 p20 -2 p50 p40 x30 .
We have the following two possible cases.

-Oscillating case. We introduce

k 2 = 1/2 + B/(2 ω 2 ) with 0 < k < 1 so that (2.24) becomes ψ2 = 4ω 2 k 2 -sin 2 (ψ/2)
and, using standard relations on elliptic functions (cf. [START_REF] Lawden | Elliptic functions and applications[END_REF]), we obtain

sin(ψ/2) = k sn(u, k), cos(ψ/2) = dn(u, k)
where u = ωt + ϕ 0 . cn and dn are elliptic functions of the first kind and the solutions of (2.21), x1 , x2 , are expressed as

ω x1 (u) = ω x10 + -2 k sin (φ ) cn (u ) + (-u + 2 E (u )) cos (φ ) ω x2 (u) = ω x20 + -2 k cos (φ ) cn (u ) + (u -2 E (u )) sin (φ ) (2.25)
where x10 and x20 are constant, and E(.) is the elliptic integral of the second kind.

-Rotating case. We introduce k 2 = 2 ω 2 /(B+ω 2 ) with 0 < k < 1 so that (2.24) becomes ψ2 = 4ω 2 /k 2 1 -k 2 sin 2 (ψ/2) .

Invoking again elliptic functions properties ( [START_REF] Lawden | Elliptic functions and applications[END_REF]) we have

sin(ψ/2) = sn(u/k, k), cos(ψ/2) = cn(u/k, k)
where u = ωt + ϕ 0 . Still sn and cn are elliptic functions of the first kind. The solutions of (2.21), x1 , x2 , satisfy the relations

ω x1 (u) = ω x10 + 1 - 2 k 2 + 2 E (k) k 2 K (k) cos (φ ) u + 2 k cos (φ ) Z u k -sin (φ ) dn u k ω x2 (u) = ω x20 + 2 k 2 -1 -2 E (k) k 2 K (k) sin (φ ) u - 2 k sin (φ ) Z u k + cos (φ ) dn u k (2.26)
where x10 and x20 are constant, K(k), E(k) are respectively the complete elliptic integrals of the first and second kind, Z(.) is the Jacobi's Zeta function.

Computations of strokes with small amplitudes using the nilpotent approximation

We recall that the physical variables q are related to x using the transformation ϕ.

The adjoint variables p are obtained by a Mathieu transformation associated with ϕ. More precisely, according to Proposition 2.1, recall that the shape variables θ = (θ 1 , θ 2 ) correspond to the ( x1 , x2 ) coordinates.

Strokes with small amplitudes such that q(0) = 0 are computed from the nilpotent approximation in the following way:

• Degenerate case: The corresponding solutions xi (.), i = 1, 2 of (2.23) yield the periodic shape variables θ i (t) = xi (t), i = 1, 2 of period 2π/ p30 . Moreover, the constants x10 , x20 , ϑ 0 may be chosen so that q(0) = (θ 1 (0), θ 2 (0), x(0)) = 0. • Oscillating case:

The modulus k can be expressed as (2.28)

k( p(0)) = 1
We have θ

1 (t) = x1 (ωt + ϕ 0 ), θ 2 (t) = x2 (ωt + ϕ 0 )
where x1 , x2 are explicitly written in (2.26). We choose p(0) so that H 1 ( x(0), p(0)) 2 + H 2 ( x(0), p(0)) 2 = 1, k( p(0)) ∈ (0, 1) and such that the denominator of k( p(0)) is nonzero. As k( p(0)) tends to 0, the linear terms of x1 (u), x2 (u) of (2.26) tend to 0. This is the case when p40 → 0 and p50 → 0, and at the limit, equation (2.22) reduces to the equation of the degenerate case: θ = 0.

Abnormal case

We can reduce the problem by considering the minimal time problem for the singleinput affine system (cf. [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]):

ẋ(t) = F1 ( x(t)) + u(t) F2 ( x(t))
where u(.) is now a scalar control. We denote by x(.) a reference minimum time trajectory, and since we consider abnormal extremals it follows from the Pontryagin maximum principle that along the extremal lift of x(.), the identity H 2 ( x, p) = 0 must hold and, differentiating with respect to t, it implies that {H 1 , H 2 }( x, p) = 0 must hold too. Differentiating once more time, the extremals associated with the controls:

u a ( x, p) = -{H 1 , {H 2 , H 1 }}( x, p) / {H 2 , {H 1 , H 2 }}( x, p) = 2 p5 / p4
satisfy the relation H 2 = {H 1 , H 2 } = 0 along ( x(.), p(.)) and are solutions of:

ẋ(t) = ∂ H a ∂ p ( x(t), p(t)), ṗ(t) = - ∂ H a ∂ x ( x(t), p(t)),
where H a is the true Hamiltonian:

H a ( x, p) = H 1 ( x, p) + u a H 2 ( x, p) = p1 + 2 p5 p2 + p3 x1 + p4 x3 + p5 x2 1 / p4 .
From the Pontryagin maximum principle, we also have that H 1 ( x(.), p(.)) = 0. The extremal system subject to the constraints H 1 = H 2 = {H 1 , H 2 } = 0 is integrable and the corresponding solutions can be written as:

x1 (t) = t + x10 , x2 (t) = 2 p50 / p40 t + x20 , x3 (t) = p50 / p40 t 2 + 2 p50 x10 / p40 t + x30 , x4 (t) = 2/3 p2 50 / p2 40 t 3 -2 p50 / p2 40 p50 x2 10 + p30 x10 + p20 t -p50 p30 / p2 40 t 2 + x40 , x5 (t) = 2/3 p50 / p40 t 3 + (4 p50 x10 + p30 ) / p40 t 2 + 2 2 p50 x2 10 + p30 x10 + x30 p40 + p20 / p40 t + x50 , p1 (t) = -2 p50 p30 -4 p2 50 x10 / p40 t + p10 , p2 (t) = p20 , p3 (t) = -2 p50 t + p30 , p4 (t) = p40 , p5 (t) = p50 with ( x10 , x20 , x30 , x40 , x50 , p10 , p20 , p30 , p40 , p50 ) are constant satisfying p10 = 0, p20 = p50 x2 10 -p40 x30 , p30 = -2 p50 x10 .
Remark 2.2. The θ -projection of abnormals are straight lines and form triangular strokes.

Numerical results

This section presents the numerical simulations performed on the Purcell swimmer problem. Simulations are performed using both direct and indirect methods, respectively with the solvers Bocop and HamPath . We use the multipliers from the solutions of the direct method to initialize the adjoint variables in the indirect approach. We display the optimal trajectories obtained for both the nilpotent approximation as well as for the true mechanical system.

BOCOP.

Bocop (www.bocop.org, [START_REF] Bonnans | Bocop -A collection of examples[END_REF]) implements a so-called direct transcription method. More precisely, a time discretization is used to rewrite the optimal control problem as a finite dimensional optimization problem (i.e nonlinear programming), solved by an interior point method (IPOPT). We recall below the optimal control problem, formulated with the state q = (θ 1 , θ 2 , x, y, α) and control u = ( θ1 , θ2 ):

         min u T 0 E(u(t)) dt q(t) = F 1 (q(t)) u 1 (t) + F 2 (q(t)) u 2 (t) x(0) = y(0) = 0, x(T ) = x f y(T ) = y f , α(T ) = α(0), θ i (T ) = θ i (0), i = 1, 2.
(2.29)

HamPath .

The HamPath software (http://www.hampath.org/, [START_REF] Cots | Contrôle optimal géométrique: méthodes homotopiques et applications[END_REF]) is based on indirect methods to solve optimal control problems using simple shooting methods and testing the local optimality of the solutions. More precisely two purposes are achieved with HamPath :

• Shooting equations: to compute periodic trajectories for the Purcell swimmer, we consider the true Hamiltonian H given by the Pontryagin maximum principle and the associated transversality conditions associated. The normal and regular minimizing curves are the projection of extremals solutions of the following boundary value problem:

       q = ∂ H ∂ p , ṗ = -∂ H ∂ q , x(0) = x 0 , x(T ) = x f , y(0) = y 0 , y(T ) = y f θ i (T ) = θ i (0), i = 1, 2 α(T ) = α(0), p θ i (T ) = p θ i (0), i = 1, 2 p α (T ) = p α (0) (2.30)
where q = (θ 1 , θ 2 , x, y, α), p = (p θ 1 , p θ 2 , p x , p y , p θ ) and the final time T > 0 is fixed. Due to the sensitivity of the initialization of the shooting algorithm, the latter is initialized with direct methods namely the Bocop toolbox.

• Local optimality: to show that the calculated normal stroke is optimal, we perform a rank test on the subspaces spanned by the solutions of the variational equation with suitable initial conditions [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF].

Using proposition 14, in the normal case it allows us to check the necessary optimality condition related to the concept of conjugate point. The same holds in the abnormal case using [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF].

Nilpotent approximation

Notations. The state variables are given by x = ( x1 , x2 , x3 , x4 , x5 ), the adjoint by p = ( p1 , p2 , p3 , p4 , p5 ), and F1 , F2 are the vector fields of the normal form given by (2.19). The Hamiltonian lifts are respectively denoted H 1 and H 2 .

Normal case

In the normal case, we consider the extremal system given by the true Hamiltonian described in (2.20). We compute the optimal trajectories with HamPath , and we display on Fig. 2.4 the state and adjoint variables as functions of time. We also illustrate the conjugate points computed according to the algorithm in [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF], as well as the smallest singular value for the rank test. Property on the first conjugate point. Let us consider the fixed energy level

(H 2 1 + H 2
2 ) |t=0 = 1 along the extremals and the initial state x(0) = 0. We take a large number of random initial adjoint vectors p(0) and numerically integrate the extremal system. For each normal extremal, we compute the first conjugate time t 1c , the pulsation ω = (p Let γ(.) be a normal extremal starting at t = 0 from the origin and defined on [0, +∞[. As illustrated on Fig. 2.5, there exists a first conjugate point along γ corresponding to a conjugate time t 1c satisfying the inequality: 0.3ωt 1c -0.4 < K(k) < 0.5ωt 1c -0.8. Remark 2.3. In section 2.6.4 u = ωt + ϕ 0 is the normalized parametrization of the solutions. 1 + H 2 2 = 1. We observe: 0.3ωt 1c -0.4 < K(k) < 0.5ωt 1c -0.8.

Abnormal case

Fig. 2.6 illustrates the time evolution of the state variables for an abnormal extremal. We check the second order optimality conditions with the algorithm described in [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]. The determinant test and the smallest singular value for the rank condition both indicate that there is no conjugate time for abnormal extremals (Fig. 2.7).

True mechanical system

We now consider the optimal control problem (2.29) consisting in minimizing either the mechanical energy (2.14) or the criterion |u| 2 .

Direct method. In the first set of simulations performed by Bocop , we set T = 10, x f = 0.5, and the bounds a = 3 large enough so that the solution is actually unconstrained. The state and the control variables for the optimal trajectory are shown on Fig. 2.8, 2.9 and 2.10, and we observe that the trajectory is actually a sequence of identical strokes. Fig. 2.11 displays the phase portrait for the shape angles θ 1 , θ 2 , which is an ellipse. The constant energy level satisfied by the optimal trajectory implies that the phase portrait of the controls is a circle for the |u| 2 criterion, but not for the energy criterion. The adjoint variables (or more accurately in this case, the multipliers associated to the discretized dynamics) are shown on Fig. 2.12-2.13. Indirect method. Now we use the multipliers from the Bocop solutions to initialize the shooting algorithm of HamPath . Fig. 2.14-2.15 and Fig. 2.16 represent respectively non intersecting strokes and an eight shape stroke. We check the second order optimality conditions according to [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF] and observe that there is no conjugate point on [0, 2π] for the non intersecting case while a conjugate point is found on [0, 2π] for the eight shape stroke.

Continuation method

Finally, we construct for the Purcell swimmer, a one parameter family of simple loops strokes using continuation methods. For the Purcell swimmer, the two families presented in Fig. 2.18 are compared in Fig. 2.19 using the efficiency concept defined as

E (γ(•)) = x(T ) 2 + y(T ) 2 /l(γ(•))
where l(γ(•)) is the length of the stroke.

Copepod swimmer

Geometric analysis of a copepod swimmer

In [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF], two types of geometric motions are described.

First case: (Fig. 2.20 (left) ) The two legs are assumed to oscillate sinusoidally according to

θ 1 = Φ 1 + a cos(t), θ 2 = Φ 2 + a cos(t + k 2 )
with a = π/4, Φ 1 = π/4, Φ 2 = 3π/4 and k 2 = π/2. This produces a displacement x 0 (2π) = 0.2. Second case: (Fig. 2.20 (right)) The two legs are paddling in sequence followed by a recovery stroke performed in unison. In this case the controls u 1 = θ1 , u 2 = θ2 produce bang arcs to steer the angles between the boundary θ i = 0 of the domain to the boundary θ i = π, while the unison sequence corresponds to a displacement from π to 0 with the constraint θ 1 = θ 2 .

Our first objective is to relate these properties to geometric optimal control.

Abnormal curves in the copepod swimmer

Let q = (x 0 , θ 1 , θ 2 ), then the system takes the form:

q(t) = 2 ∑ i=1 u i (t)F i (q(t))
where the control vector fields are given by:

F i = sin θ i ∆ ∂ ∂ x 0 + ∂ ∂ θ i , ∆ = 2 + sin 2 θ 1 + sin 2 θ 2 .
The Lie brackets in the copepod case are easily calculated and are given by:

F 3 = [F 1 , F 2 ] = f (θ 1 , θ 2 ) ∂ ∂ x 0 with f (θ 1 , θ 2 ) = 2 sin θ 1 sin θ 2 (cos θ 1 -cos θ 2 ) ∆ 2 , [[F 1 , F 2 ], F 1 ] = ∂ f ∂ θ 1 (θ 1 , θ 2 ) ∂ ∂ x 0 , [[F 1 , F 2 ], F 2 ] = ∂ f ∂ θ 2 (θ 1 , θ 2 ) ∂ ∂ x 0 .
Lemma 3 The singular set Σ : {q; det(F 1 (q), F 2 (q), [F 1 , F 2 ](q)) = 0}, where the vector fields F 1 , F 2 , [F 1 , F 2 ] are coplanar, is given by 2 sin θ 1 sin θ 2 (cos θ 1 -cos θ 2 ) = 0 which is equivalent to:

• θ i = 0 or π i = 1, 2, • θ 1 = θ 2
and corresponds to the boundary of the physical domain: θ i ∈ [0, π], θ 1 ≤ θ 2 , with respective controls u 1 = 0, u 2 = 0 or u 1 = u 2 forming a stroke of triangular shape in the phase portait of the variables θ 1 , θ 2 .

Remark 2.4. Each point of the boundary is a Martinet point except at the non smooth points (vertices).

The previous lemma provides the interpretation of the triangle shape stroke in terms of abnormal curves.

To understand smooth stroke strategies via optimal control we must introduce the cost function related to the mechanical energy. Recall that according to [START_REF] Passov | Supplementary notes to: Dynamics of Purcells three-link microswimmer with a passive elastic tail[END_REF] the mechanical energy of the copepod swimmer is given by: T 0 qt M q dt where q = (x 0 , θ 1 , θ 2 ) and M is the symmetric matrix:

M =   2 -1/2(cos 2 (θ 1 ) + cos 2 (θ 2 )) -1/2 sin(θ 1 ) -1/2 sin(θ 2 ) -1/2 sin(θ 1 ) 1/3 0 -1/2 sin(θ 2 ) 0 1/3   . (2.31)
Taking into account the constraints on the velocities, the integrand can be written as:

a(q)u 2 1 + 2b(q)u 1 u 2 + c(q)u 2 2 where a = 1 3 - sin 2 θ 1 2(2 + sin 2 θ 1 + sin 2 θ 2 ) , b = - sin θ 1 sin θ 2 2(2 + sin 2 θ 1 + sin 2 θ 2 ) , c = 1 3 - sin 2 θ 2 2(2 + sin 2 θ 1 + sin 2 θ 2 )
.

The pseudo-Hamiltonian is then expressed as:

H(q, p, p 0 ) = u 1 H 1 (q, p) + u 2 H 2 (q, p) + p 0 a(q)u 2 1 + 2b(q)u 1 u 2 + c(q)u 2 2 .
Taking p 0 = -1/2, the normal controls are computed by solving the equations:

∂ H ∂ u 1 = 0, ∂ H ∂ u 2 = 0.
We obtain:

u 1 = - 3(4H 1 + 2H 1 sin 2 θ 1 + 3H 2 sin θ 1 sin θ 2 -H 1 sin 2 θ 2 ) sin 2 θ 1 + sin 2 θ 2 -4 , u 2 = - 9H 1 sin θ 1 sin θ 2 + 6H 2 (2 + sin 2 θ 2 ) -3H 2 sin 2 θ 1 sin 2 θ 1 + sin 2 θ 2 -4 .
and plugging this control u back into the pseudo-Hamiltonian provides the true Hamiltonian which we denote by H n .

Note also that H n can also be obtained by constructing an orthonormal basis of the metric using a feedback transformation u = β (q)v to transform the problem into:

q = (Fβ (q))(v), min v(.) T 0 (v 2 1 (t) + v 2 2 (t)) dt
where F is the matrix (F 1 , F 2 ). Writing Fβ = (F 1 , F 2 ), F 1 , F 2 will form an orthonormal frame. The computation is straightforward and the normal Hamiltonian H n takes the form

H n = 1 2 (H 1 2 + H 2 2
) where H i is the Hamiltonian lift of F i .

The concept of efficiency

To compare strokes with different amplitudes we introduce the following definition of efficiency [START_REF] Lighthill | Note on the swimming of slender fish[END_REF].

Definition 37

The efficiency of a stroke γ(•) is defined by:

E(γ(•)) = x 0 (T )/L(γ(•))
where x 0 is the displacement of the swimmer and L is the length of the curve γ(.).

The transversality condition given in Exercise 1.1 can be generalized, see [START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF]. For instance, for the copepod swimmer, considering the augmented adjoint vector (p, p 0 ), the transversality condition implies that:

(p(T ), p 0 (T )) is collinear to the gradient of the set E = c, where c is a constant.

Geometric classification of smooth strokes

The expected strokes are related to the classification of smooth periodic curves in the plane up to a diffeomorphism, assuming that in our discussion we relax the state constraints on the shape variable. This problem was studied by Whitney (1937) and Arnold (1994) , see [START_REF] Berger | La taxonomie des courbes[END_REF]. In this classification we have in particular the three cases of Fig. 2.21. Each of this curve has a specific physical interpretation for the swimmer problem.

Numerical computations

• Micro-local analysis. First, we compute the normal strokes using the Maximum Principle to recover the strokes displayed in Fig. 2.21. Below, we present the numerical calculations of these strokes using the weak Maximum Principle. An important point is to account for the transversality conditions associated with the periodicity requirement θ i (0) = θ i (2π), i = 1, 2 which are given by:

p θ i (0) = p θ i (2π), i = 1, 2.
The solutions are computed via a shooting method using the HamPath code. Finally, we evaluate numerically the value function which reduces to 2πH n the given reference geodesic, since H n is constant. • Second order optimality. Conjugate points are computed for each type of stroke which leads to select simple loops as candidates for minimizers, see Fig. 2.22. • Abnormal triangle. To deal with the global optimality problem we use the geometric efficiency E = x 0 /L for single loops constrained in the triangle (see Fig. 2.24 and Table 2.1). From our analysis we deduce that the (triangle) abnormal stroke is not optimal. Indeed, one can choose a normal stroke (inside the triangle) such that the displacement is x0 /2 with x0 = 2.742 and length < L/2 where L =length of the triangle. Applying twice the normal stroke, we obtain the same displacement x0 than with the abnormal stroke but with a length < L. Therefore, we proved the following theorem.

Theorem 8 The abnormal triangle is not optimal for both costs: minimizing length with fixed displacement or maximizing the efficiency. 2.500e-1 3.353 7.456e-2 Simple loop with small amplitude 0.500e-1 9.935e-1 5.033e-2 Table 2.1 Ratio x 0 /L for the abnormal stroke and different normal strokes corresponding to the mechanical cost.

So far, the copepod microswimmer was analyzed using mainly simulations but a complete analysis can be obtained combining mathematical analysis based on numerical evidence. We proceed as follows.

First, to simplify the computations and to have a clear interpretation of the pictures in the Euclidean frame, we replace the mechanical energy by the Euclidean cost T 0 |u(t)| 2 dt. Note that the true cost case can be analyzed using a numerical continuation between the two costs (HamPath software).

Using the nilpotent approximation and Lemma 3, one must consider two cases with respect to the triangle T associated with the state constraints: 0 ≤ θ 1 ≤ θ 2 ≤ π.

Point interior to the triangle. Take such a point q = (x 0 , θ 1 , θ 2 ). Then near the chosen point, there exists privileged coordinates x = ( x1 , x2 , x3 ) such that the nilpotent SRmodel is given by the Dido model:

ẋ = u 1 F1 ( x) + u 2 F2 ( x), min u T 0 (u 1 (t) 2 + u 2 2 (t)) dt with F1 = ∂ ∂ x1 + x2 ∂ ∂ x3 , F2 = ∂ ∂ x2 -x1 ∂ ∂ x3 .
This model implies that starting from each q we have a one parameter family of symmetric simple strokes (see Fig. Points on the sides of the triangle but different of the vertices. Take such a point q = (x 0 , θ 1 , θ 2 ). Then the SR-nilpotent model is the Martinet flat case. Thus, one can find privileged coordinates x = ( x1 , x2 , x3 ) such that the model is:

ẋ = u 1 F1 ( x) + u 2 F2 ( x), min u T 0 (u 2 1 (t) + u 2 2 (t)) dt where F1 = ∂ ∂ x1 + x2 2 2 ∂ ∂ x3 , F2 = ∂ ∂ x2 .
This model leads to the calculation of eight strokes parameterized by elliptic functions which correspond to lemniscates of Bernoulli.

All these models are not stable models and higher order approximations can be used to generate strokes with small amplitudes. Also by perturbation at a interior point of the triangle, we can obtain limac ¸on's strokes by doubling the period. This is indeed confirmed by numerical simulations using the true model and represented on Fig. 2.26. Moreover for the true system with the Euclidean cost, the numerical simulations show the existence of a one-parameter family of simple strokes symmetric with respect to the axis D : θ 2 = -θ 1 + π. They are obtained by integrating from D identified to a cross-section and with a tangent vector taken normal to D, each stroke being associated with a different energy level, see Fig.

It leads to the following proposition.

Proposition 15 There exists a one parameter family of simple strokes, symmetric with respect to the D-axis and foliating the interior of the triangle T , each associated to a different energy level.

The final result of our analysis is captured in the following theorem.

Theorem 9 Among this family of strokes, there exists a unique stroke with maximal efficiency among all the strokes of the copepod swimmer for the

T 0 (u 2 1 (t)+u 2 2 (t)) dt cost.

Sketch of the proof

First we have the following lemma.

Lemma 2.1. For the Euclidean case (or the mechanical energy case) the geodesic flow is invariant under the transformation δ :

(θ 1 , θ 2 , x 0 ) → (π -θ 2 , π -θ 1 , x 0 ).
From this, we deduce that the one parameter family of simple loops represented on Fig. 2.27 is symmetric with respect to the straight line D. The center of this family can be calculated as follows. We choose a point θ (0) = (θ 1 (0), θ 2 (0)) on the line D which can be identified to (0, 0) if we introduce the new coordinates x = θ 1 -θ 1 (0), y = θ 2 -θ 2 (0). Using a transformation of the form Z = zc 1 xc 2 y we get a graded set of coordinates (x, y, Z) with weights (1, 1, 2) establishing a link between the physical coordinates (θ 1 , θ 2 , x 0 ) and the privileged coordinates identified to (x, y, Z). Using this gradation, the nilpotent (order -1) SR-model is given by the Dido model. This model is not stable under perturbation and higher order terms have to be taken into account. In particular, using the weights (1, 1, 2) the model of order 0 can be computed. Using the analysis of [START_REF] Brockett | Control theory and singular Riemannian geometry[END_REF], the model of order zero can be identified with the model of order -1 using diffeomorphism and feedback preserving the Euclidean energy. A precise computation detailed in [START_REF] Bettiol | Sub-Riemannian geometry and swimming at low Reynolds number: the copepod case[END_REF] shows that the only point θ (0) such that the diffeomorphism is not mixing the shape variable θ with the displacement variable x 0 corresponds to the center θ (0) (0.72, 2.41) of Fig. 2.27. Hence, we proved that there exists only one point to generate such a family of simple loops (compare with Fig. 2.18 in the Purcell case). Now, we must prove that the only strokes candidates as minimizers in the interior of the triangle are simple strokes. This can be proved using the Stokes theorem and the following lemma. Lemma 2.2. Consider the smooth one-form on R 2 : ω :

= 2 ∑ i=1 sin θ i ∆ (θ ) dθ i with ∆ (θ ) = 2 + sin 2 θ 1 + sin 2 θ 2 and introduce f (θ ) = 2 sin θ 1 sin θ 2 (cos θ 1 -cos θ 2 )/∆ (θ ) 2 .
Then,

1. dω = -f (θ 1 , θ 2 )dθ 1 ∧ dθ 2 .
2. For any bounded Stokes domain D ⊂ R 2 , we have

∂ D ω = D dω and if γ is a piecewise smooth stroke with γ = ∂ D the associated displacement is x 0 (T ) = D dω.
3. dω < 0 in the interior of the triangle T : 0 ≤ θ 1 ≤ θ 2 ≤ π, and dω vanishes on the boundary of T formed by the abnormal stroke.

In particular this lemma allows to compare efficiency of simple loops versus limac ¸ons and eight shape strokes in the interior of the triangle.

Another method from optimal control theory is to compute conjugate points. This can be performed by numerical computations but more theoretical computations are related to conjugate loci computations on the SR-sphere. In particular, for limac ¸ons with small amplitudes, conjugate points can be estimated as follows. According to the Dido model, the only strokes with small amplitudes can be either simple loops or limac ¸ons, obtained by perturbation of a simple loop followed twice. For the Dido model, using the explicit computation, the first conjugate point appears on a simple loop after exactly one period. By perturbation, for a simple stroke with small amplitude, the first conjugate time corresponds approximately to the period. Hence a limac ¸on of small amplitude produced by period doubling has necessarily a conjugate point. This gives a rigorous proof of the existence of conjugate point for limac ¸ons with small amplitude.

Conclusion and bibliographic remarks

We made a short presentation of the problem of microswimming using the Purcell and the copepod case in the frame of SR-geometry, combining analytic and numeric methods in optimal control based on the analysis of the geodesic flow to determine the most efficient stroke. A different approach combining Stokes theorem to determine the shape of optimal strokes and direct numeric methods using Fourier analysis were used earlier in a series of articles, see for instance [START_REF] Avron | A geometric theory of swimming: Purcell's swimmer and its symmetrized cousin[END_REF].

Note also that the copepod case is the analog of a limit of symmetric Purcell swimmer described and analyzed in [START_REF] Avron | A geometric theory of swimming: Purcell's swimmer and its symmetrized cousin[END_REF].

The approaches are complementary. The main result of this theory is the existence of center of swimmings from which are emanating a one parameter family of simple strokes to compute the most efficient stroke. See [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF] for an earlier computation using a shooting method.

Note also the (geometric) link of microswimmers in SR-geometry with the geodesic motion of a 2D-particle in a magnetic field very well presented in [START_REF] Montgomery | Isoholonomic problems and some applications[END_REF]. This leads to a fine and technical study in [START_REF] Agrachev | On the Dido problem and plane isoperimetric problems[END_REF] as a generalization of the Dido problem, to compute conjugate and cut loci for small lengths. Such results being applicable to generate in general conjugate and cut loci, using numeric continuation methods.

Chapter 3

Maximum Principle and Application to Nuclear Magnetic Resonance and Magnetic Resonance Imaging

Maximum Principle

In this section we state the Pontryagin maximum principle and we outline the proof. We adopt the presentation from Lee and Markus [START_REF] Lee | Foundations of optimal control theory[END_REF] where the result is presented into two theorems. The complete proof is complicated but rather standard, see the original book from the authors [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF].

Theorem 10 We consider a system of R n : ẋ(t) = f (x(t), u(t)), where f : R n+m → R n is a C 1 -mapping. The family U of admissible controls is the set of bounded measurable mappings u(•), defined on [0, T ] with values in a control domain Ω ⊂ R m such that the response x(•, x 0 , u) is defined on [0, T ]. Let ū(•) ∈ U be a control and let x(•) be the associated trajectory such that x(T ) belongs to the boundary of the accessibility set A(x 0 , T ). Then there exists p(•) ∈ R n \ {0}, an absolutely continuous function defined on [0, T ] solution almost everywhere of the adjoint system:

ṗ(t) = -p(t) ∂ f ∂ x ( x(t), ū(t)) (3.1)
such that for almost every t ∈ [0, T ] we have

H( x(t), p(t), ū(t)) = M( x(t), p(t)) (3.2) where H(x, p, u) = p, f (x, u) and M(x, p) = max u∈Ω H(x, p, u). Moreover t → M( x(t), p(t)) is constant on [0, T ].
Proof. The accessibility set is not in general convex and it must be approximated along the reference trajectory x(•) by a convex cone. The approximation is obtained by using needle type variations of the control ū(•) which are closed for the L 1 -topology. (We do not use L ∞ perturbations and the Fréchet derivative of the endpoint mapping computed in this Banach space.)

Needle type approximation

We say that 0 ≤ t 1 ≤ T is a regular time for the reference trajectory if

d dt | t=t 1 t 0 f ( x(τ), ū(τ))dτ = f ( x(t 1 ), ū(t 1 ))
and from measure theory we have that almost every point of [0, T ] is regular. At a regular time t 1 , we define the following L 1 -perturbation ūε (•) of the reference control: we fix l, ε ≥ 0 small enough and we set

ūε (t) = u 1 ∈ Ω constant on [t 1 -lε,t 1 ] ū(t) otherwise on [0, T ] .
We denote by xε (•) the associated trajectory starting at xε (0) = x 0 . We denote by ε → α t (ε) the curve defined by α t (ε) = xε (t) for t ≥ t 1 . We have

xε (t 1 ) = x(t 1 -lε) + t 1 t 1 -lε f ( xε (t), ūε (t))dt where ūε = u 1 on [t 1 -lε,t 1 ], Moreover x(t 1 ) = x(t 1 -lε) + t 1 t 1 -lε f ( x(t), ū(t))dt
and since t 1 is a regular time for x(•) we have xε (t 1 ) -x(t 1 ) = lε( f ( x(t 1 ), u 1 ) -f ( x(t 1 ), ū(t 1 )) + o(ε).

In particular if we consider the curve ε → α t 1 (ε), it is a curve with origin x(t 1 ) and whose tangent vector is given by

v = l( f ( x(t 1 ), u 1 ) -f ( x(t 1 ), ū(t 1 ))). (3.3)
For t ≥ t 1 , consider the local diffeomorphism: φ t (y) = x(t,t 1 , y, ū) where x(•,t 1 , y, ū) is the solution corresponding to ū(•) and starting at t = t 1 from y. By construction we have α t (ε) = φ t (α t (ε)) for ε small enough and moreover for t

≥ t 1 , v t = d dε | ε=0 α t (ε)
is the image of v by the Jacobian ∂ φ t ∂ x . In other words v t is the solution at time t of the variated equation

dV dt = ∂ f ∂ x ( x(t), ū(t))V (3.4)
with condition v t = v for t = t 1 . We can extend v t on the whole interval [0, T ]. The construction can be done for an arbitrary choice of t 1 , l and u 1 . Let Π = {t , l, u 1 } be fixed, we denote by v Π (t) the corresponding vector v t .

Additivity property

Let t 1 ,t 2 be two regular points of ū(•) with t 1 < t 2 and l 1 , l 2 small enough. We define the following perturbation

ūε (t) =    u 1 on [t 1 -l 1 ε,t 1 ] u 2 on [t 2 -l 2 ε,t 2 ] ū(t); otherwise on [0, T ]
where u 1 , u 2 are constant values of Ω and let xε (•) be the corresponding trajectory. Using the composition of the two elementary perturbations

Π 1 = {t 1 , l 1 , u 1 } and Π 2 = {t 2 , l 2 , u 2 } we define a new perturbation Π : {t 1 ,t 2 , l 1 , l 2 , u 1 , u 2 }. If we denote by v Π 1 (t), v Π 2 (t)
and v Π (t) the respective tangent vectors, a computation similar to the previous one gives us:

v Π (t) = v Π 1 (t) + v Π 2 (t), for t ≥ t 2 .
We can deduce the following lemma.

Lemma 4 Let Π = {t 1 , • • • ,t s , λ 1 l 1 , • • • , λ s l s , u 1 , • • • , u s } be a perturbation at regular times t i ,t 1 < • • • < t s , l i ≥ 0, λ i ≥ 0, ∑ s i=1 λ i = 1
and corresponding to elementary perturbations Π i = {t i , l i , u i } with tangent vectors v Π i (t). Let xε (•) be the associated response with perturbation Π . Then we have

xε (t) = x(t) + s ∑ i=1 ελ i v Π i (t) + o(ε) (3.5)
where o(ε) ε → 0, uniformly for 0 ≤ t ≤ T and 0 ≤ λ i ≤ 1. Definition 38 Let ū(•) be an admissible control and x(•) its associated trajectory defined for 0 ≤ t ≤ T . The first Pontryagin's cone K(t), 0 < t ≤ T is the smallest convex cone at x(t) containing all elementary perturbation vectors for all regular times t i .

Definition 39 Let v 1 , • • • , v n be linearly independent vectors of K(t), each v i being formed as convex combinations of elementary perturbation vectors at distinct times. An elementary simplex cone C is the convex hull of the vectors v i .

Lemma 5 Let v be a vector interior to K(t). Then there exists an elementary simplex cone C containing v in its interior.

Proof. In the construction of the interior of K(t), we use the convex combination of elementary perturbation vectors at regular times not necessarily distinct. Clearly by continuity we can replace such a combination with n distinct times.

Approximation lemma

An important technical lemma is the following topological result whose proof uses the Brouwer fixed point theorem.

Lemma 6 Let v be a nonzero vector interior to K(t), then there exists λ > 0 and a conic neighborhood N of λ v such that N is contained in the accessibility set A(x 0 , T ).

Proof. See [START_REF] Lee | Foundations of optimal control theory[END_REF].

Separation step

To finish the proof, we use the geometric Hahn-Banach theorem. Indeed if x(T ) ∈ ∂ A(x 0 , T ) there exists a sequence x n / ∈ A(x 0 , T ) such that x n → x(T ) when n → +∞ and the unit vectors x n -x(T )

|x n -x(T )| have a limit ω when n → ∞. The vector ω is not interior to K(T ) otherwise from Lemma 6 there would exist λ > 0 and a conic neighborhood of λ ω in A(x 0 , T ) and this contradicts the fact that x n / ∈ A(x 0 , T ) for any n. Let π be any hyperplane at x(T ) separating K(T ) from ω and let p be the exterior unit normal to π at x(T ). Let us define p(•) as the solution of the adjoint equation

ṗ(t) = -p(t) ∂ f ∂ x ( x(t), ū(t))
satisfying p(T ) = p. By construction we have p(T )v(T ) ≤ 0 for each elementary perturbation vector v(T ) ∈ K(T ) and since for t ∈ [0, T ] the following equations hold:

ṗ(t) = -p(t) ∂ f ∂ x ( x, ū), v(t) = ∂ f ∂ x ( x, ū)v we have d dt p(t)v(t) = 0. Hence p(t)v(t) = p(T )v(T ) ≤ 0, ∀t. Assume that the maximization condition (3.2)
is not satisfied on some subset S of 0 ≤ t ≤ T with positive measure. Let t 1 ∈ S be a regular time, then there exists

u 1 ∈ Ω such that p(t 1 ) f ( x(t 1 ), ū(t 1 )) < p(t 1 ) f ( x(t 1 ), u 1 ).
Let us consider the elementary perturbation Π 1 = {t 1 , l, u 1 } and its tangent vector

v Π 1 (t 1 ) = l [ f ( x(t 1 ), u 1 ) -f ( x(t 1 ), ū(t 1 ))] .
Then using the above inequality we have that p(t 1 )v Π 1 (t 1 ) > 0 which contradicts p(t 1 )v Π 1 (t 1 ) ≤ 0, for all t. Therefore the inequality

H( x(t), p(t), ū(t)) = M( x(t), p(t))
is satisfied almost everywhere on 0 ≤ t ≤ T . Using a standard reasoning we can prove that t → M( x(t), p(t)) is absolutely continuous and has zero derivative almost everywhere on 0 ≤ t ≤ T , see [START_REF] Lee | Foundations of optimal control theory[END_REF].

Theorem 11 Let us consider a general control system:

ẋ(t) = f (x(t), u(t))
where f is a continuously differentiable function and let M 0 , M 1 be two C 1 submanifolds of R n . We assume the set U of admissible controls to be the set of bounded measurable mappings u

: [0, T (u)] → Ω ∈ R m ,
where Ω is a given subset of R m . Consider the following minimization problem:

min u∈U C(u), C(u) = T 0 f 0 (x(t), u(t))dt where f 0 ∈ C 1 , x(0) ∈ M 0 , x(T ) ∈ M 1
and T is not fixed. We introduce the augmented system:

ẋ0 (t) = f 0 (x(t), u(t)), x 0 (0) = 0, (3.6)
ẋ(t) = f (x(t), u(t)), (3.7)
x(t) = (x 0 (t), x(t)) ∈ R n+1 , f = ( f 0 , f ). If (x * (•), u * (•) is optimal on [0, T * ], then there exists p * (•) = (p 0 , p(•)) : [0, T * ] → R n+1 \ {0} absolutely continuous, such that ( x * (•), p * (•), u * (•))
satisfies the following equations almost everywhere on 0 ≤ t ≤ T * :

ẋ(t) = ∂ Ĥ ∂ p (x(t), p(t), u(t)), ṗ(t) = - ∂ Ĥ ∂ x (x(t), p(t), u(t)) (3.8)
Ĥ(x(t), p(t), u(t)) = M(x(t), p(t)) (

where Ĥ(x(t), p(t), u(t)) = p, f (x, u) , M( x, p) = max u∈Ω Ĥ( x, p, u).

Moreover, we have M(x(t), p(t)) = 0, ∀t, p 0 ≤ 0 (3.10)

and the boundary conditions (transversality conditions):

x * (0) ∈ M 0 , x * (T * ) ∈ M 1 , (3.11)
p * (0) ⊥ T x * (0) M 0 , p * (T * ) ⊥ T x * (T ) M 1 . (3.12)
Proof. (For the complete proof, see [START_REF] Lee | Foundations of optimal control theory[END_REF] or [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF].) Since (x * (•), u * (•)) is optimal on [0, T * ], the augmented trajectory t → x * (t) is such that x * (T ) belongs to the boundary of the accessibility set Â(x * (0), T * ). Hence by applying Theorem 10 to the augmented system, one gets the conditions (3.8), (3.9) and M constant. To show that M ≡ 0, we construct an approximated cone K (T ) containing K(T ) but also the two vectors ± f (x * (T ), u * (T )) using time variations (the transfer time is not fixed).

To prove the transversality conditions, we use a standard separation lemma as in the proof of Theorem 10.

Definition 40 A triplet (x(•), p(•), u(•)) solution of the maximum principle is called an extremal.

Special cases Minimal Time

Consider the time minimum case: f 0 = 1. In this case, an optimal control u * on [0,t * ] is such that the corresponding trajectory x * (.) is such that for each t > 0, x * (t) belongs to the boundary of the accessibility set A(x * (0),t). The pseudo-Hamiltonian of the augmented system is written:

Ĥ( x, p, u) = H(x, p, u) + p 0
with H(x, p, u) is the reduced pseudo-Hamiltonian and since p 0 ≤ 0, conditions 3.9, 3.10 become

H(x * (t), p * (t), u * (t)) = M(x * (t), p * (t)) a.e.
with M(x, p) = Max u∈Ω H(x, p, u) and M(x * (t), p * (t)) ≥ 0 is constant everywhere.

Mayer Problem

A Mayer problem is an optimal control problem for a system dx dt = f (x, u), u ∈ Ω , x(0) = x 0 , where the cost to be minimized is of the form:

Min u∈Ω c(x(t f ))
where c : R n -→ R is smooth the transfer time t f is fixed and the final boundary conditions are of the form: g(x(t f )), with g : R n → R k is smooth. In this case the maximum principle and the geometric interpretation of this principle lead to the following:

• Each optimal control u * on [0,t f ] with response x * (.) is such that x * (t f ) belongs to the boundary of the accessibility set A(x 0 ,t f ) and at the final point the adjoint vector p * (t f ) is orthogonal to the manifold defined by the boundary conditions and the cost function:

M : { x; g(x) = 0, c(x) = m }
where m is the minimum.

Introducing the pseudo-Hamiltonian

H(x, p, u) = p, f (x, u)
the necessary optimality conditions are:

dx * dt = ∂ H ∂ p (x * , p * , u * ), dp * dt = - ∂ H ∂ x (x * , p * , u * ), (3.13)
H(x * , p * , u * ) = max u∈Ω H(x * , p * , u) (3.14)
and the following boundary conditions

f (x * (t f )) = 0, (3.15)
p * (t f ) = p 0 . ∂ c ∂ x (x * (t f )) + δ . ∂ g ∂ x (x * (t f )),
p 0 ≤ 0 (transversality conditions).

(3.16)

Exercise 3.1. Write the necessary optimality conditions for a Bolza problem where the cost problem is of the form:

C(u) = g(x(t f )) + t f 0 f 0 (x(t), u(t))dt.

Application to NMR and MRI

Optimal control was very recently applied very successively to a general research project initiated by S. Glaser: the control of spins systems with applications to Nuclear Magnetic Resonance (NMR) and Magnetic Resonance Imaging (MRI). Such success is partially explained by an accurate representation of the control problem by the Bloch equations introduced in 1946 and F. Bloch and E.M. Purcell were awarded the 1952 Nobel Prizes for Physics for "their development of new ways and method for NMR", Purcell providing a nice link between our two working examples.

Next, we make a mathematical introduction of Bloch equations and the concept of resonance, in order to model the class of associate problems objects of our research program.

Model

The Bloch equations are a set of macroscopic equations which accurately describe the experimental model in NMR and MRI [START_REF] Levitt | Spin dynamics: basics of nuclear magnetic resonance[END_REF] based on the concept of the dynamics of a spin-1/2 particle. At this level it is represented by a magnetization vector M = (M x , M y , M z ) in the laboratory reference frame which evolves according to

dM dτ = γM ∧ B + R(M) (3.17)
where γ is the gyromagnetic ratio, B = (B x , B y , B z ) is the applied magnetic field which decomposes into a strong polarizing field B z = B 0 in the z-direction, while B x , B y are the components of a Rf-magnetic field in the transverse direction and corresponds to the control field and R(M) is the dissipation of the form:

- M x T 2 , - M y T 2 , - (M z -M 0 ) T 1
where T 1 , T 2 are the longitudinal and transverse relaxation parameters characteristic of the chemical species, e.g. water, fat, blood.

The parameter M 0 is normalized to 1 up to a rescaling of M. We denote ω 0 = -γB 0 the resonant frequency and let introduce the control components: u(τ) = -γB y and v(τ) = -γB x . The Bloch equations in the stationary frame can be written in the matrix form:

d dτ   M x M y M z   =   0 -ω 0 u(τ) ω 0 0 -v(τ) -u(τ) v(τ) 0     M x M y M z   -    M x T 2 M y T 2 M z -1 T 1    . (3.18)
The Bloch equations can be represented in a rotating frame of reference: S(τ) = exp(τωΩ z ), M = S(τ)q, q = (x, y, z),

Ω z =   0 -1 0 1 0 0 0 0 0  
and introducing the control representation:

u 1 = u cos ωτ -v sin ωτ u 2 = u sin ωτ + v cos ωτ,
one gets the Bloch equations in the moving frame:

d dτ   x y z   =   0 -∆ ω u 2 ∆ ω 0 -u 1 -u 2 u 1 0     x y z   -    x T 2 y T 2 z-1 T 1    (3.19)
where ∆ ω = ω 0 -ω is the resonance offset.

The control is bounded by m, m = 2π × 32.3Hz being the experimental intensity of the experiments. Assuming ∆ ω = 0 (resonance), and using the normalized time t = τm, denoting Γ = 1/mT 2 , γ = 1/mT 1 and the physical parameters satisfying the constraint: 2Γ ≥ γ ≥ 0, the system is normalized to:

dx dt = -Γ x + u 2 z dy dt = -Γ y -u 1 z dz dt = γ(1 -z) + u 1 y -u 2 x, (3.20)
where |u| ≤ 1. Moreover since 2Γ ≥ γ ≥ 0, one has that the Bloch ball |q| ≤ 1 is invariant for the dynamics. This equation describes the evolution of the magnetization vector in NMR. The choice of ω = ω 0 corresponding to resonance neutralized the existence of the strong polarizing field B 0 , except the side effect of a stable unique equilibrium, corresponding to the North pole N = (0, 0, 1) of the Bloch equation for the uncontrolled motion.

In MRI, the situation is more complex due to the spatial position of the spin in the image and one must control an ensemble of spins corresponding to each pixel. Due to this localization they are some distortions corresponding to B 0 and B 1 inhomogeneities. The variation of B 0 producing a resonance offset and ∆ ω belongs to some intervals, while B 1 -inhomogeneity introduces a scaling factor a i ≥ 0 depending on the spatial position of the spin in the image modeling the distortion of the amplitude of the control field and the equation transforms into

dx dt = -Γ x + a i u 2 z dy dt = -Γ y -a i u 1 z dz dt = γ(1 -z) + a i (u 1 y -u 2 x). (3.21)
In the general case one must consider both inhomogeneities producing a detuning and amplitude alteration. Note that such distortions cannot be modelized and have to be experimentally determined.

To relate Bloch equation to imaging we associate to the amplitude |q| of the normalized magnetization vector a grey level where |q| = 1 corresponds to white while the center of the Bloch ball defined by |q| = 0 corresponds to black.

The problems

Having introduced the control systems in relation with Bloch equations taking into account B 0 and B 1 inhomogeneities one can present the fundamental problems studied in NMR and MRI.

Saturation problem

The objective of the saturation problem for a single spin is to bring the magnetization vector q from the North pole N = (0, 0, 1) of the Bloch ball (which is the equilibrium of the free system) to the center O = (0, 0, 0) of the Bloch ball, recalling that |q| corresponds to a grey level where the sphere |q| = 1 corresponds to white and |q| = 0 to black.

A direct generalization being to consider an ensemble of spin particles corresponding to the same chemical species and to bring each spin of this ensemble from the North pole to the center, corresponding to the multisaturation problem.

The contrast problem

In the contrast problem in NMR called ideal contrast problem we consider two pairs of (uncoupled) spin-1/2 systems corresponding to different chemical species, each of them solutions of the Bloch equations (3.21) with respective parameters (γ 1 ,Γ 1 ) and (γ 2 ,Γ 2 ) controlled by the same magnetic field. Denoting each system by dq i dt = f (q i ,Λ i , u), Λ i = (γ i ,Γ i ) and q i = (x i , y i , z i ) is the magnetization vector representating each spin particle, i = 1, 2. This leads to the consideration of the system abbreviated as: dq dt = f (q, u), where q = (q 1 , q 2 ) is the state variable. The constrast problem by saturation is the following optimal control problem: starting from the equilibrium point q 0 = ((0, 0, 1), (0, 0, 1)) where both chemical species are white and hence indistinguishable, reach in a given transfer time t f the final state q 1 (t f ) corresponding to saturation of the first spin while maximizing |q 2 (t f )| 2 , the final observed contrast being |q 2 (t f )|.

Obvious generalization of the problems in MRI, taking into account B 0 and B 1 inhomogeneities, is to consider in the image an ensemble of N pairs of chemical species, e.g. water or fat, and distributed in the image and the objective is to provide multisaturation of the ensemble of spins of the first species and to reach for the second species a small ball centered at |q 2 (t f )| where |q 2 (t f ) corresponds to the contrast calculated in NMR.

The objective in MRI is to produce a robust control taking into account the B 0 and B 1 inhomogeneities. In the sequel and in order to present the concepts and the theoretical tools, we shall restrict to the saturation problem of a single spin and the contrast problem by saturation in NMR. It is the preliminary step to the analysis of an ensemble of spins which is in the applications treated numerically using adapted software e.g. Bocop and

HamPath representative respectively of direct and indirect methods in numeric optimal control.

The saturation problem in minimum time for a single spin

The saturation problem in minimum time was first analyzed in [START_REF] Lapert | Développement de nouvelles techniques de contrôle optimal en dynamique quantique : de la Résonance Magnétique Nucléaire à la physique moléculaire[END_REF] and was an important step to the applications of geometric optimal control to the dynamics of spins particles.

Preliminaries

First of all, since the transfer is from the North pole N = (0, 0, 1) to the center of the Bloch ball O = (0, 0, 0) which belongs to the z-axis of revolution of the system corresponding to polarization the system can be restricted to the two-dimensional plane of the Bloch ball and the control u = (u 1 , u 2 ) reduces to the u 1 component. The system is compactly written as: dq dt = F(q) + u 1 G(q), while the control is bounded by |u| ≤ 1 and q = (y, z). We have

F = -Γ y ∂ ∂ y + γ(1 -z) ∂ ∂ z G = -z ∂ ∂ y + y ∂ ∂ z . (3.22)
According to the maximum principle an optimal trajectory is a concatenation of bang arcs where u(t) = sign p(t), G(q(t)) and singular arcs where p(t), G(q(t)) = 0. The following Lie brackets are relevant in our analysis. Denoting δ = γ -Γ , we have

[G, F] = (-γ + δ z) ∂ ∂ y + δ y ∂ ∂ z [[G, F], F] = (γ(γ -2Γ ) -δ 2 z) ∂ ∂ y + δ 2 y ∂ ∂ z [[G, F], G] = 2δ y ∂ ∂ y + (γ -2δ z) ∂ ∂ z .

Singular trajectories and optimality

The singular trajectories are located on the set S : det( G, [G, F] ) = 0 which is given by y(-2δ z + γ) = 0. Hence it is formed by

• the z-axis of revolutions y = 0,

• the horizontal line z = γ/(2δ ). This line intersects the Bloch ball |q| < 1 when 2Γ > 3γ and moreover z is negative. The singular control is given by D + u s D = 0, where

D = det( G, [[G, F], G] ) and D = det( G, [[G, F], F] ).
• for y = 0, D = -z(γ -2δ z) and D = 0. The singular control is zero and a singular trajectory is a solution of ẏ = -y, ż = γ(1 -z) where the equilibrium point (0, 1) is a stable node if γ = 0. • for z = γ/(2δ ), D = -2δ y 2 , D = yγ(2Γ -γ) and u s = γ(2Γ -γ)/(2δ y), 2Γγ ≥ 0. Hence along the horizontal direction, the flow: ẏ = -Γ yγ 2 2Γ -γ 4δ 2 y and |u s | → ∞ when y → 0.

An easy computation gives the following proposition.

Proposition 16 If γ = 0, the singular control along the singular line is L 1 but not L 2 .

The maximum principle selects the singular line but the high order maximum principle and the so-called generalized Legendre-Clebsch condition [START_REF] Krener | The high order maximum principle and its application to singular extremals[END_REF] has to used to distinguish between small time minimum and maximum solution. It can be easily understood using the two seminal examples:

ẋ = 1 -u 2 , ẋ = 1 + u 2 , ẏ = u, |u| ≤ 1, ẏ = u, |u| ≤ 1
where in both case the x-axis is the singular line and is time minimizing in the first case and time maximizing in the second case. The optimality condition takes the following form in our case. Let D = det(G, F) = γz(z -1) + Γ γ 2 . The set C : D = 0 is the collinear set. If γ = 0, this set forms an oval joining the North pole to the center of the Bloch ball and the intersection with the singular line is empty.

Denoting D = det(G, [[G, F], G]
) the singular lines are fast displacement direction if DD > 0 and slow if DD < 0. From this condition, one deduces that the z-axis of revolution is fast if 1 > z > z = γ/(2δ ) and slow if z = γ/(2δ ) > z > -1, while the horizontal singular line is fast.

From the analysis we deduce first Lemma 7 If the condition 2Γ > 3γ is not satisfied the horizontal singular line doesn't intersect the Bloch ball |q| < 1 and the optimal solution is the standard inversion sequence used in practices: apply a bang control to steer (0, 1) to (0, - * ).

Followed by u = 0 to relax the system to (0, 0) along the z-axis.

If 2Γ > 3γ the existence of the fast displacement horizontal line will determine the optimal policy. First of all observe that since u S → ∞, when q → 0 along this line, it is saturating the constraint |u| < 1 at a point of this line. Hence this line has to be quitted before this point. The exact exit point is determined by the maximum principle because such point has to be a switching point at both extremities for the corresponding bang arc. Such an arc is called a bridge.

Note that in this analysis we assume that the applied RF-field is large enough, which correspond to the experimental situation. We deduce the following theorem, see [START_REF] Lapert | Développement de nouvelles techniques de contrôle optimal en dynamique quantique : de la Résonance Magnétique Nucléaire à la physique moléculaire[END_REF] for further details.

Theorem 12 If 2Γ > 3γ, in the time minimal saturation problem is of the form: δ + δ h s δ + δ v s , concatenating the bang arc to quit the North pole to the horizontal singular line, followed by the bridge and relaxation to 0 along the z-axis of revolution. Remark 3.1. The bridge can be empty and in this case the optimal policy is δ + δ v s .

y z ® δ + ® δ v s ® δ + ® δ h s ® δ + ® δ v s N • O • • • • • Fig. 3
This gives a complete solution to the saturation problem using a careful geometric analysis to understand the interaction between the two singular lines. Moreover a similar analysis leads to a complete understanding of the time minimum synthesis to transfer any point of the Bloch ball to the center.

Extension of this type of results to an ensemble of two or more spins is an important issue. The complexity is related to the analysis of singular extremals at two levels. First of all, in general the symmetry of revolution due to z-polarization cannot be invoked to reduced the bi-inputs case to the single single-input case. Secondly, in dimension ≥ 3, the analysis of the singular flow even in the single-input case is a complicated task. Next, we shall present this complexity in the contrast problem and present some achievements.

The maximum principle in the contrast problem by saturation

The system is written as: q = F 0 (q) + u 1 F 1 (q) + u 2 F 2 (q), |u| ≤ 1 where q = (q 1 , q 2 ) ∈ { |q 1 | ≤ 1, |q 2 | ≤ 1 } and q 1 , q 2 represents the normalized magnetization vector of the first and second spin, q i = (x i , y i , z i ), i = 1, 2. Using the notation of the section 3.2 for a Mayer problem, the cost function is c(q(t f )) = -|q 2 (t f )| 2 and the final boundary condition is F(q(t f )) = q 1 (t f ) = 0. Splitting the adjoint vector into p = (p 1 , p 2 ) ∈ R 3 × R 3 , the transversality condition is:

p 2 (t f ) = -2p 0 q 2 (t f ), p 0 ≤ 0 and if p 0 = 0 it can be normalized to p 0 = -1/2. We denote z = (q, p), H i = p, F i (q) , i = 0, 1, 2, the Hamiltonian lift of the sys-

tem ż = - → H 0 + ∑ 2 i=1 - → H i (z). If (H 1 , H 2 ) = 0
, the maximization condition of the maximum principle leads to the following parametrization of the controls

u 1 = H 1 H 2 1 + H 2 , u 2 = H 2 H 2 1 + H 2 2 .
Define the switching surface:

Σ : H 1 = H 2 = 0.
Plugging such u into the pseudo-Hamiltonian gives the true Hamiltonian:

H n = H 0 + H 2 1 + H 2 2 .
The corresponding extremal solutions are called zero. Besides those generic extremals, additional extremals are related to Lie algebraic properties of the system and a careful analysis is the key factor to determine the properties of the optimal solutions.

Lie bracket computations

Due to the bilinear structure of the Bloch equations, Lie brackets can be easily computed, which is crucial in our analysis.

Recall that the Lie bracket of two vectors fields F, G is computed with the convention

[F, G](q) = ∂ F ∂ q (q)G(q) - ∂ G ∂ q (q)F(q)
and if H F , H G are the Hamiltonian lifts, recall that the Poisson bracket is

{H F , H G }(z) = dH F ( - → G )(z) = H [F,G] (z).
To simplify the computation, each spin system is lifted on the semi-direct Lie product GL(3, R) × R 3 acting on the q-space using the action (A, a).q = Aq + a. The Lie bracket computation rule is

((A, a), B, b) = ([A, B], Ab-Ba) where [A, B] = AB-BA. Introducing F 0 = (A 0 , a 0 ), with A 0 = diag(-Γ 1 , -Γ 1 , -γ 1 , -Γ 2 , -Γ 2 , -γ 2
) and a 0 = (0, 0, γ 1 , 0, 0, γ 2 ) whereas the control fields (F 1 , F 2 ) are identified to B 1 = diag(C 1 ,C 1 ) and B 2 = diag(C 2 ,C 2 ) where C 1 ,C 2 are the antisymmetric matrices

C 1 = E 32 -E 23 , C 2 = (E 13 -E 31
) with E i j = (δ i j ) (Kronecker symbol). See [START_REF] Bonnard | Singular Trajectories and the Contrast Imaging Problem in Nuclear Magnetic resonance[END_REF] for more details.

Next, we present in details the Lie brackets needed in our computations, each entry form by a couple (v 1 , v 2 ) and we use the notation omitting the indices. We set δ = γ -Γ .

• Length 1:

F 0 = (-Γ x, -Γ y, γ(1 -z)) F 1 = (0, -z, y) F 2 = (z, 0, -x).
• Length 2: Let z = (q, p) be a curve solution of

[F 0 , F 1 ] = (0, γ -δ z, -δ y) [F 0 , F 2 ] = (-γ + δ z, 0, δ x) [F 1 , F 2 ] = (-y, x, 0). • Length 3: [[F 1 , F 2 ], F 0 ] = 0 [[F 1 , F 2 ], F 1 ] = F 2 [[F 1 , F 2 ], F 2 ] = -F 1 [[F 0 , F 1 ], F 1 ] = (0, -2δ y, -γ + 2δ z) [[F 0 , F 1 ], F 2 ] = (δ y, δ x, 0) = [[F 0 , F 2 ], F 1 ] [[F 0 , F 2 ], F 2 ] = (-2δ x, 0, 2δ z -γ) [[F 0 , F 1 ], F 0 ] = (0, -γ(γ -2Γ ) + δ 2 z, -δ 2 y) [[F 0 , F 2 ], F 0 ] = (γ(γ -2Γ ) -δ 2 z, 0, δ 2 x).
- → H 0 + u 1 - → H 1 + u 2 - → H 2 .
Differentiating H 1 and H 2 along such a solution, one gets:

Ḣ1 = {H 1 , H 0 } + u 2 {H 1 , H 2 } Ḣ2 = {H 2 , H 0 } + u 1 {H 2 , H 1 }. (3.23)
Hence we have:

Proposition 17 Let z 0 ⊂ Σ 1 = Σ \ {H 1 , H 2 } = 0
and define the control u s by:

u s (z) = (-{H 0 , H 2 }(z), {H 0 , H 1 }(z)) {H 1 , H 2 }(z) , (3.24)
and plugging such u s into H defines the true Hamiltonian

H s (z) = H 0 (z) + u s,1 (z)H 1 (z) + u s,2 (z)H 2 (z)
which parameterizes the singular solutions of the bi-input system contained in Σ 1 .

This gives the first stratum of the surface Σ . Moreover, the behaviors of the extremals of order zero near a point z 0 of Σ 1 can be easily analyzed using (3.23) and a nilpotent model where all Lie brackets at z 0 ∈ Σ 1 of length ≥ 3 are zero. Denoting:

{H 1 , H 0 }(z 0 ) = a 1 , {H 2 , H 0 }(z 0 ) = a 2 , {H 2 , H 1 }(z 0 ) = b
and using polar coordinates H 1 = r cos θ , H 2 = r sin θ , then (3.23) becomes:

ṙ = a 1 cos θ + a 2 sin θ θ = 1 r (b -a 1 sin θ + a 2 cos θ ). (3.25)
To analyze this equation, we write: At such a point z 0 of Σ 1 , the singular control (3.24) is such that

a 1 sin θ -a 2 cos θ = A sin(θ + φ ) with A tan φ = -a 2 /a 1 , A = a 2 1 + a 2 2 .
u 2 s,1 + u 2 s,2 = a 2 1 + a 2 2 b 2 > 1
and hence is not admissible.

Next we analyze more degenerated situations and one needs the following concept.

Goh condition

Higher order necessary optimality conditions along singular extremals in the biinput case are related to finiteness of the index of the quadratic forms associated with the second order derivative [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF] known as Goh condition which is the relation:

{H 1 , H 2 } = 0. (3.26)
Using H 1 = H 2 = {H 1 , H 2 } = 0 and (3.23), one gets the additional relations:

{H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = 0. (3.27)
Then differentiating again along a solution leads to the relations:

{{H 1 , H 2 }, H 0 } + u 1 {{H 1 , H 2 }, H 1 } + u 2 {{H 1 , H 2 }, H 2 } = 0 (3.28) {{H 0 , H 1 }, H 0 } + u 1 {{H 0 , H 1 }, H 1 } + u 2 {{H 0 , H 1 }, H 2 } = 0 {{H 0 , H 2 }, H 0 } + u 1 {{H 0 , H 2 }, H 1 } + u 2 {{H 0 , H 2 }, H 2 } = 0 (3.29)
This leads in general to three relations to compute two control components, and for a generic system such conditions are not satisfied [START_REF] Chitour | Genericity results for singular curves[END_REF], but in our case, according to Lie brackets computations, we have:

Lemma 9 If H 1 = H 2 = 0, one has {{H 1 , H 2 }, H 0 } = {{H 1 , H 2 }, H 1 } = {{H 1 , H 2 }, H 2 } = 0
and (3.28) is satisfied.

The equation (3.29) are then written: Ã + Bu and if det( B) = 0, the corresponding singular control is given by:

u s (z) = -B-1 (z) Ã(z) (3.30)
Using the relations:

H 1 = H 2 = {H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = 0, the vector p is orthogonal to F 1 , F 2 , [F 1 , F 2 ], [F 0 , F 1 ], [F 0 , F 2 ]. Introducing: A = A 1 A 2 , B = B 1 B 3 B 2 B 4 , C = (F 1 , F 2 , [F 1 , F 2 ], [F 0 , F 1 ], [F 0 , F 2 ]), with A 1 = det(C, [[F 0 , F 1 ], F 0 ]), A 2 = det(C, [[F 0 , F 2 ], F 0 ]),
and

B 1 = det(C, [[F 0 , F 1 ], F 1 ]), B 2 = det(C, [[F 0 , F 2 ], F 1 ]), B 3 = det(C, [[F 0 , F 1 ], F 2 ]), B 4 = det(C, [[F 0 , F 2 ], F 2 ]),
the relation (3.29) leads to: A + Bu = 0, and if det B = 0, one gets the singular control given by the feedback: u s (q) = -B -1 (q)A(q) (3.31)

and the associated vector field:

Q s = F 0 + u s,1 F 1 + u s,2 F 2 .
Moreover, the singular control has to be admissible: |u s | ≤ 1. We introduce the stratum:

Σ 2 : H 1 = H 2 = {H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } \ det B = 0.
Hence we have:

Lemma 10 1. On the stratum Σ 2 , there exists singular extremals satisfying Goh condition where the singular control is given by the feedback (3.30). 2. For the contrast problem:

det B =(x 1 y 2 -x 2 y 1 ) 4 (δ 1 -δ 2 )(2δ 1 z 1 -γ 1 )(2δ 2 z 2 -γ 2 ) 2(δ 2 1 γ 2 z 1 -δ 2 2 γ 1 z 2 ) -γ 1 γ 2 (δ 1 -δ 2 ) -2δ 1 δ 2 (γ 1 z 2 -γ 2 z 1 ) , (3.32)
The behaviors of the extremals of order zero near a point z 0 ∈ Σ 2 is a complicated problem. Additional singular extremals can be contained in the surface:

Σ 3 : H 1 = H 2 = {H 1 , H 2 } = {H 0 , H 1 } = {H 0 , H 2 } = det B = 0,
and they can be computed using the property that the corresponding control has to force the surface det B = 0 to be invariant. Some have an important meaning, due to the symmetry of revolution of the Bloch equations. They correspond to control the system, imposing u 2 = 0. In this case, one can restrict the system to

Q = {q = (q 1 , q 2 ) ∈ R n : |q 1 | ≤ 1, |q 2 | ≤ 1, x 1 = x 2 = 0}.
The computations of the corresponding extremals amount to replace in the relations; H 2 by εH 2 and to impose ε = 0. The remaining relations are then:

H 1 = {H 0 , H 1 } = 0
and from (3.29) one gets the relations:

{{H 0 , H 1 }, H 0 } + u 1,s {{H 0 , H 1 }, H 1 } = 0, (3.33)
and thus, this defines the singular control:

u 1,s = - {{H 0 , H 1 }, H 0 } {{H 0 , H 1 }, H 1 } (3.34)
and the associated Hamiltonian H 1,s = H 0 + u 1,s H 1 . We have the following result:

Proposition 18 The singular extremals of the single-input case with u 2 ≡ 0 are extremals of the bi-input case with the additional condition:

x 1 = p x 1 = x 2 = p x 2 = 0.
Moreover from the geometric interpretation of the maximum principle for a Mayer problem, in order to be optimal the generalized Legendre-Clebsch condition has to be satisfied:

∂ ∂ u 1 d 2 dt 2 ∂ H ∂ u 1 = {H 1 , {H 1 , H 0 }}(z) ≤ 0. (3.35)
Observe that if we impose u 2 = 0, the classification of the extremals near the switching surface, which reduces to H 1 = 0, is a standard problem [START_REF] Kupka | Geometric theory of extremals in optimal control problems. i. the fold and Maxwell case[END_REF].

Finally, another important property of the extremal flow, again a consequence of the symmetry of revolution is given next, in relation with Goh condition. It is a consequence of Noether integrability theorem.

Proposition 19

In the contrast problem, for the Hamiltonian vector field -→ H n whose solutions are extremals of order zero, the Hamiltonian lift H

(z) = {H 1 , H 2 }(z) = (p x 1 y 1 -p y 1 x 1 ) + (p x 2 y 2 -p y 2 x 2 ) is a first integral.
Exercise 3.2 (Generalization to the case of B 1 and B 0 inhomogeneities). It is interesting to compare to the case of an ensemble of two spins of the same spin particle with B 0 and B 1 inhomogeneities which is left to the reader. More precisely:

• B 1 -inhomogeneities.
In this case, the control directions of the second spin are relaxed by a factor and the Lie brackets computations can be used to stratified. It can be applied to the multisaturation problem. • B 0 -inhomogeneities.

In this case the vector field F 0 of the second spin contains a non-zero detuning.

Clearly this introduces modifications in the Lie brackets computations. Again it can be applied to multisaturation problem. It explains the following phenomenon: in the precense of detuning both controls (u 1 , u 2 ) have to be used.

Next, motivating by the fact that due to the symmetry of revolution and the observed numerical experiments, we shall restrict our study to the single-input case. It is an important theoretical step since we can reduce the analysis of the singular flow for a 4-dimensional system with one input vs a 6-dimensional system. This complexity will be illustrated by the computations presented next.

The classification of the singular extremals and the action of the feedback group

Preliminairies

Restricting to the single input case, the research program concerning the contrast problem or the multisaturation problem for an ensemble of two spins is clear.

Saturation problem for a single spin and bridge phenomenon

In the case of a single spin the complete geometric analysis requires the computations of the two singular line and the understanding of the singularity associated with their intersection, which causes the saturation of the singular control and the occurrence of a bang arc called a bridge to connect both singular arcs. This phenomenon generalizes to higher dimension and it tells you that the analysis of the singular flow codes all the information of the optimal solution which is a sequence of arcs of the form δ ± δ S δ ± ...δ S , where δ ± denotes bang arcs with u = ±1, while δ S are singular arcs. This will be presented in details next, in relation with the action of the feedback group.

Computations of the singular flow

Consider a control system of the form: dq dt = F(q) + uG(q), q ∈ R n and relaxing the control constraints: u ∈ R. Denoting H F and H G the Hamiltonian lifts of F and G, if the denominator is not vanishing, a singular control is given by:

u S (z) = - {{H G , H F }, H F }(z) {{H G , H F }, H G }(z) . (3.36)
Plugging such u S into the pseudo-Hamiltonian one gets the true Hamiltonian: H S = H F + u S H G and the singular extremals are solutions of the constrained Hamiltonian equation:

dz dt = - → H S (z), z ∈ Σ : H G = {H G , H F } = 0.
This set of equations defines a Hamiltonian vector field on the surface

Σ \ {{H G , H F } , H G } = 0
, restricting the standard symplectic from ω = dp ∧ dq.

We We get the system:

dq ds -DF -D F, dp ds = -p D ∂ F ∂ q -D ∂ G ∂ q
restricted to the surface Σ .

In the contrast problem, since the state space is of dimension four, using the two constraints H G = {H G , H F } = 0 and the homogeneity with respect to p, equation (3.36) can be reduced to the explicit form:

dq dt = F(q) - D (q, λ ) D(q, λ ) G(q)
where λ is a one-dimension time-dependant parameter whose dynamics is deduced from the adjoint equation.

Using the previous remark, the optimal problem can be analyzed by understanding the behavior of the corresponding trajectories and the singularities of the flow near the set D = 0, which codes the switching sequence. This is a very complicated task, in particular because the system is depending upon four parameters and simplifications have to be introduced to simplify this task. Two simplifications can be introduced. First, we can restrict to some specific parameters corresponding to some experimental cases. For instance, in the water case, saturation of a single spin amounts to the standard inversion sequence. Second, a projection of the singular flow which is physically relevant can be introduced. A natural choice is to consider the case where the transfer time t is not fixed. Then according to the maximum principle this leads to the additional constraint: M = Max u(.) H F + uH G = 0, which gives in the singular case the additional constraint: H F (z) = 0. This case is called the exceptional case using the terminology of [START_REF] Bonnard | Théorie des singularités de l'application entrée/sortie et optimalité des trajectoires singulières dans le problème du temps minimal[END_REF].

With this constraint, the adjoint vector can be eliminated and the singular control in this exceptional case is the feedback:

u e S = - D (q) D(q) , where D = det(F, G, [G, F], [[G, F], G]), D = det(F, G, [G, F], [[G, F], F]
) with the corresponding vector field X e defined by dq dt = F(q) -D (q) D(q) G(q) which can again be desingularized using the reparametrization ds = dt/D(q(t)) and this gives the smooth vector field X e r = DF -D G.

Feedback classification

Definition 41 Let E and F be two R-vector spaces and let G be a group acting linearly on E and F. A homomorphism

X : G → R \ {0} is called a character. A semi-invariant of weight X is a map λ : E → R such that for all g ∈ G and all x ∈ E, λ (g, x) = X(g)λ (x); it is an invariant if X = 1. A map λ : E → F is a semi-covariant of weight X if for all g ∈ G and for all x ∈ E, λ (g.x) = X(g)g.λ (x); it is called a covariant if X = 1.
More about invariant theory can be found in [START_REF] Dieudonné | Invariant theory, old and new[END_REF].

The key concept in analyzing the role of relaxation parameters in the control problem is the action of the feedback group G on the set of systems. We shall restrict our presentation to the single-input case and we denote C = {F, G} the set of such (smooth) systems on the state space Q R n , see [START_REF] Bonnard | Feedback equivalence for nonlinear systems and the time optimal control problem[END_REF] for the details.

Definition 42 Let (F, G), (F , G ) be two elements of C . They are called feedback equivalent if there exist a smooth diffeomorphism ϕ of R n and a feedback u = α(q)+ β (q)v, β invertible such that:

• F = ϕ * F + ϕ * (Gα), G = ϕ * (Gβ ).
where ϕ * z denotes the image of the vector field.

Definition 43 Let (F, G) ∈ C and let λ 1 be the map which associated the constrained Hamiltonian vector field ( -→ H S , Σ ) (see (3.36)) to (F, G). We define the action of (ϕ, α, β ) of G on ( -→ H S , Σ ) to be the action of the symplectic change of coordinates:

-

→ ϕ : q = ϕ(Q), p = P ∂ ϕ ∂ x -1
in particular the feedback acts trivially.

Theorem 13 ([20])

The mapping λ 1 is a covariant.

Next, we detail the induced action restricting to exceptional singular trajectories when dim Q = 4.

The exceptional singular trajectories and the feedback classification

Notation. Let ϕ be a diffeomorphism of Q. Then ϕ acts on the mapping F : Q → R according to ϕ.F = Foϕ and on vector fields X as ϕ.X = ϕ * X ( image of X): this corresponds to the action on tensors.

The feedback group acts on the vector field X e by change of coordinates only and this can be checked as a consequence of the following lemma.

Lemma 11 • D F+αG,β G = β 4 D F,G . • D F+αG,β G = β 3 D F,G + αD F,G . • D ϕ * F,ϕ * G (q) = det ∂ ϕ ∂ q -1 D F,G (ϕ(q)). • D ϕ * F,ϕ * G (q) = det ∂ ϕ ∂ q -1 D F,G (ϕ(q)).
From which we deduce the following crucial result in our analysis.

Theorem 14 We have the following:

• λ 2 : (F, G) → X e is a covariant. • λ 3 : (F, G) → D is a semi-covariant. • λ 4 : (F, G) → X e r = DF -D G is a semi-covariant.
The classification program. Having introduced the concepts and results, the contrast problem is related to the following classification program (up to change of coordinates)

• Classification of the vector fields X e r = DF -DF and the surfaces:

D = 0, D = D = 0.
Interpretation.

• The singular control is u e S = -D D and will explode at D = 0 except if D = 0, taking into account the (non isolated) singularities of X e r (if D = D = 0, X e r = 0). Collinear set. The collinear set of F, G is a feedback invariant which has also an important meaning in our classification. Remark 3.2. In our classification program we use semi-covariants and in the set of parameters Λ = (γ 1 ,Γ 1 , γ 2 ,Γ 2 ) it amounts to work in the projective space. It is also clear from our reparametrization of time. Now, the problem is to test the computational limits of our program which is clearly:

• Compute the surfaces D = 0, D = D = 0, • Compute the equilibrium points of X e r . Clearly, in the framework of computational methods in real algebraic geometry it is a complicated task which has been achieved in two cases.

• The multisaturation problem of two spins taking into account B 1 -inhomogeneity.

• The contrast problem when the first spin system corresponds to water (γ 1 = Γ 1 ).

The second problem has application in in vivo, where the parameters are varying, in particular in the brain.

We shall present the results in details in the first case.

Algebraic classification in the multisaturation of two spins with B 1 -inhomogeneity

The point N = ((0, 1), (0, 1)) is a singular point of X r e and under a translation N is taken as the origin of the coordinates. We have:

F 0 = (-Γ y 1 , -γz 1 , -Γ y 2 , -γz 2 ), F 1 = ((-(z 1 + 1), y 1 ), (1 -ε)(-(z 2 + 1), y 2 ))
where (1 -ε) denotes the control rescaling of the second spin.

We have D = (1 -ε) D, where D is a quadric which decomposes into h 2 + h 3 + h 4 where h i are the homogeneous part of degree i:

h 2 = (2Γ -γ) h2 h2 = Γ (2Γ -γ) ((ε -1) y 1 + y 2 ) 2 + γ 2 (ε -1) 2 z 2 1 -γ 2 2 -2 ε + ε 2 z 2 z 1 + γ 2 z 2 2 h 3 = 2(γ -Γ ) h3 h3 = (γ -2Γ ) γ + 2Γ (ε -1) 2 z 2 y 2 1 + (γ -2Γ ) (γ + 2Γ ) (ε -1) (y 2 z 1 + z 2 y 2 ) y 1 - γ 2 ε (ε -2) z 2 z 2 1 + (γ -2Γ ) 2Γ + (ε -1) 2 γ y 2 2 + γ 2 ε (ε -2) z 2 2 z 1 h 4 = 4(γ -Γ ) 2 h4 h4 = γ + (ε -1) 2 Γ z 2 2 y 2 1 + 2 (ε -1) (γ + Γ ) z 2 y 2 z 1 y 1 + Γ + (ε -1) 2 γ y 2 2 z 2 1 D = 2γ 2 (Γ -γ)(2Γ -γ)(1 -ε)(z 1 -z 2 )((ε -1) z 1 y 2 + z 2 y 1 ).
In particular we deduce (compare with [START_REF] Bonnard | Algebraic geometric classification of the singular flow in the contrast imaging problem in nuclear magnetic resonance[END_REF] in the contrast problem):

Proposition 20 The quadric D reduces to a cubic form which is factorized into a linear and a quadratic (homogeneous) forms.

Singular analysis

We assume γ > 0 and 2Γ > 3γ. It implies γ = Γ and γ = 2Γ . The main result is the following:

Theorem 15 Provided ε = 1 the equilibrium points of X r e = DF 0 -D F 1 are all contained in {D = D = 0}.
A simple proof exists, but we present a method based on symbolic computation and Gröbner basis.

Proof. Obviously, every point of {D = 0} ∩ {D = 0} is a singularity of X r e . Conversely, let us assume ε = 1. We first divide X r e by 1 -ε. We still assume that Γ = 0. We consider the equations {(X r e ) y 1 = 0, (X r e ) z 1 = 0, (X r e ) y 2 = 0, (X r e ) z 2 = 0} and remark that the last third are dividable by γ. By homogeneity, changing γ into γΓ , we get rid of Γ . So we may assume Γ = 1. The resulting system is denoted Σ r . We add the two polynomials ((ε -1) z 1 y 2 + z 2 y 1 ) a 1 -1 and (z 1 -z 2 ) a 2 -1, and the polynomials γg -1, (γ -1)g 1 -1, (γ -2)g 2 -1. We denote Σr this new system, involving four new variables g 1 , g 2 , a 1 , a 2 . We compute a Gröbner basis with total degree with reverse lexicographic order on (y 1 , y 2 , z 1 , z 2 , ε, g, g 1 , g 2 , a 1 , a 2 ) and get {1}. Hence, provided γ is different from 0, 1, 2, there is no singular point of X r e outside of {D = 0} ∩ {D = 0}.

The remaining of the section is devoted to the singularity resolution. From the factorized form of D (Proposition 20) we get:

Proposition 21 {D = 0} ∩ {D = 0} is an algebraic variety of algebraic dimension 2 whose components are located in the hyperplane z 1 = z 2 and in the hypersurface (ε -1) z 1 y 2 + z 2 y 1 = 0. These components are studied in the following analysis, and explicitly expressed in Lemmas 12, 13, 14, 15.

• Case

A: components of {D = 0} ∩ {D = 0} in z 1 = z 2 .
Under the constraint z 1 = z 2 , we have a factorization D = p 1 p 2 with:

p 1 = 2 (γ -Γ ) z 1 + γ -2Γ
and:

p 2 = 2 (γ -Γ ) γ + (ε -1) 2 Γ z 1 + Γ (ε -1) 2 (γ -2Γ ) y 2 1 + (4 (γ -Γ ) (γ + Γ ) (ε -1) z 1 + 2Γ (ε -1) (γ -2Γ )) y 2 y 1 + 2 (γ -Γ ) Γ + (ε -1) 2 γ z 1 + Γ (γ -2Γ ) y 2 2 .
The first polynomial has one root

z 1 = z γ,Γ z γ,Γ = 1 2 2Γ -γ γ -Γ which corresponds to the plane-solution {(y 1 , z γ,Γ , y 2 , z γ,Γ ), (y 1 , y 2 ) ∈ R 2 }.
We put:

d 2 (y 1 , y 2 ) = γ + (ε -1) 2 Γ y 2 1 + 2 (ε -1) (γ + Γ ) y 2 y 1 + Γ + (ε -1) 2 γ y 2 2 .
The discriminant of d 2 with respect to y 1 is -4 (ε -2) 2 γ Γ ε 2 y 2 2 which is strictly negative provided ε = 0. So d 2 is non-zero outside y 1 = y 2 = 0. So, provided y 2 1 + y 2 2 = 0, d 2 = 0, and p 2 = 0 is solved with respect to z 1 . We get

z 1 = r 2 (y 1 , y 2 ) with r 2 (y 1 , y 2 ) = Γ (2Γ -γ) ((ε -1) y 1 + y 2 ) 2 2 (γ -Γ ) d 2 (y 1 , y 2 ) -D = 2γ 2 (Γ -γ)(2Γ -γ)(1 -ε)(z 2 -z1 )[(-1 + z1 )y 2 (ε -1) + (-1 + z2 )y 1 ].
Conclusion: these computations allow to evaluate the equilibrium points and the behaviors of the solutions near such point, using linearization methods. A first step towards the global behavior is the following result.

Lemma 16

The surface y 1 = y 2 = 0 is foliated by lines solutions connecting O to the north pole N, the singular control being zero.

Proposition 22 Singular points on y 1 = y 2 = 0, z 1 = z 2 = z are such that: in the coordinates q = (y 1 , y 2 , z 1 -z 2 , z 1 ) the system takes the form

q = A q + R( q) where A =     0 0 0 0 0 0 0 0 0 0 0 0 0 -γ 3 ε z2 (ε -2) (2 δ + γ + 2 zδ ) 0 0     . • At the North Pole, A = 0, R q = O(| q| 3 ). • At the point S = (0, z s , 0, z s ) where z s = γ-2Γ 2(Γ -γ) , A = 0, R( q) = O(| q| 2 ).
Locally the trajectories can be computed using a blowing-up.

Numerical simulations, the ideal contrast problem

This section is devoted to numerical simulation in the ideal control problem using three complementary softwares:

• Bocop : direct method,

• HamPath : indirect method,

• GloptiPoly : Lmi technique to estimate the global optimum.

The algorithms based on the softwares are presented in details in [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF].

The ideal contrast problem by saturation in the single-input case, can be summarized this way:

       c(q(t f )) = -|q 2 (t f )| 2 -→ min u(•) , fixed t f q = F 0 (q) + u 1 F 1 (q), q(0) = q 0 q 1 (t f ) = 0 (ICPS) where q = (q 1 , q 2 ), q i = (y i , z i ) ∈ R 2 , |q i | ≤ 1, i = 1, 2.
The initial condition for each spin is q i (0) = (0, 1). The vector fields F 0 and F 1 are given by:

F 0 (q) = ∑ i=1,2 (-Γ i y i ) ∂ ∂ y i + (γ i (1 -z i )) ∂ ∂ z i , F 1 (q) = ∑ i=1,2 -z i ∂ ∂ y i + y i ∂ ∂ z i ,
where Λ i = (γ i ,Γ i ) are the physical parameters representing each spin. We present the simulations using the numerical methods (see [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF] for a complete description of the algorithms).

The simulations correspond to the two following sets of experimental data, with the relaxation times in seconds and T min the solution of the time minimal saturation problem for a single spin, from section 3.3.3. Optimal solutions of the contrast problem are concatenations of bang and singular extremals. For the following sections, we introduce some notations. We note BS the sequence composed by one bang arc (δ + or δ -) followed by one singular arc (δ s ), and nBS, n > 1, the concatenation of n BS-sequences.

First results with fixed final time

The first difficulty comes from the discontinuities of the optimal control structure. We need to know the control structure (meaning the number of Bang-Singular sequences) before calling the multiple shooting method. The indirect method also typically requires a reasonable estimate for the control switching times, as well as the states and costates values at the initial and switching times. We use the Bocop software based upon direct methods to obtain approximate optimal solutions in order to initialize the indirect shooting, within the HamPath code. We recall that the costate (or adjoint state) for Pontryagin's Principle corresponds to the Lagrange multipliers for the dynamics constraints in the discretized problem, and can therefore be extracted from the solution of the direct method.

The only a priori information is the value of the minimum time transfer T min , used to set the final time t f in the [T min , 3T min ] range. We note t f = λ T min with λ in [START_REF] Agrachev | On curvatures and focal points of dynamical Lagrangian distributions and their reductions by first integrals[END_REF][START_REF] Agrachev | Abnormal sub-Riemannian geodesics: Morse index and rigidity[END_REF]. The state variables are initialized as constant functions equal to the initial state, i.e. y 1 (

•) = y 2 (•) = 0, z 1 (•) = z 2 (•) = 1.
For the control variables we use the three constant initializations u 1 (•) ∈ {0.1, 0.25, 0.5}. The discretization method used is implicit midpoint (2nd order) with a number of time steps set to λ × 100. In order to improve convergence, we add a small regularization term to the objective to be minimized, ε reg t f 0 |u(t)| 2 dt, with ε reg = 10 -3 .

We repeat the optimizations for λ in {1.1, 1.5, 1.8, 2.0, 3.0} with the three control initializations, see Table . 3.1. The solutions from Bocop are used to initialize the continuations in HamPath , and we discuss in the following sections the results obtained with the indirect method. Both methods confirm the existence of many local solutions, as illustrated on Fig. 3.2 for λ = 1.5, due in particular to symmetry reasons.

Second order conditions

According to proposition 3.2 from [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF], the non-existence of conjugate points on each singular arc of a candidate solution is a necessary condition of local optimality. See [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF] for details about conjugate points in the contrast problem. Here, we compute for each singular arc of all the solutions from subsection 3.3.8, the first conjugate point along the arc, applying the algorithm presented in Sect. 4.3 from [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF]. None of the solutions has a conjugate point on a singular arc. Hence all the solutions satisfy the second order necessary condition of local optimality. Fig. 3.3 represents the computations of the two conjugate points (since the structure is 2BS) of the best solution with λ = 2.0 from subsection 3.3.8. [START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF] is evaluated along the two singular arcs. See [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] for details on the concept of conjugate times. On the left subplot, for each singular arc, the curve is reparameterized so that the final time corresponds to the abscissa 1 (vertical blue dashed line); the determinant associated with the rank condition is plotted (top subgraph), so there is a conjugate time whenever it vanishes (vertical red dashed lines). One observes that conjugate times on each arc are located after the (normalized to 1) final time, satisfying necessary condition of local optimality of the trajectory. At the bottom, the smallest singular value of the matrix whose rank we test is plotted, extracting only the relevant information to detect the rank drops. On the right subplot is presented a zoom of top-left subgraph near the two conjugate times.

Influence of the final time

Given that the initial point (the North pole) is a stationary point, the constrast is an increasing function of t f acting as a parameter. Indeed, applying a zero control at t = 0 leaves the system in its initial state so there is an inclusion of admissible controls between problems when the final time is increased (and the bigger the set of controls, the larger the maximum contrast). Having increasing bounded (by one, which is the maximum possible contrast given the final condition on spin no. 1) functions, it is natural to expect asymptotes on each branch.

In both cases P 1 and P 2 , the contrast problem has many local solutions, possibly with different control structures. Besides, the structure of the best policy can change depending on the final time. The possible change of structure along a single path of zeros is emphasized in Fig. 3.4. In this figure, the branch made of 2BS solutions is represented in blue, whereas the 3BS branch is the dashed red line. We also show a crossing between two value functions of two different paths of zeros in Fig. 3.5.

Then for each solution of each branch the second order necessary condition is checked as in subsection 3.3.8: the first conjugate point of each singular extremal is computed. There is no failure in this test, hence all the solutions satisfy the necessary second order condition of local optimality. Sub-optimal syntheses in fluid and blood cases

We give the syntheses of locally optimal solutions obtained in the blood and fluid cases. Note that in the special case t f = T min , for both cases the solution is 2BS and of the form δ + δ s δ + δ s . For the fluid case, the left subplot of Fig. 3.7 represents all the different branches we obtained by homotopy on λ . The greatest two value functions intersect around t f = 1.048T min . The right subplot shows the sub-optimal synthesis. The best policy is: δ + δ s δ + δ s for λ ∈ [1.000, 1.006], δ + δ s δ + δ s δ -δ s for λ ∈ For the blood case, the results are excerpted from [START_REF] Cots | Contrôle optimal géométrique: méthodes homotopiques et applications[END_REF]. The left subplot of Fig. 3.8 shows the contrast for five different components of {h = 0}, for final times t f ∈ Fig. 3.8 Blood case, sub-optimal synthesis. Illustration on the left subplot, of local solutions (each branch corresponds to a control structure). Best policy as solid lines, local solutions as dashed lines. The suboptimal synthesis is plotted on right subplot. The colors are black for BS structure, blue for 2BS and red for 3BS. The best policy is BS for t f ∈ (1, 1.294)T min and 3BS of the form δ + δ s δ -δ s δ -δ s for t f ∈ (1.294, 2]T min . In the special case t f = T min , the solution is 2BS and of the form δ + δ s δ + δ s

Sub-optimal syntheses compared to global results

We now apply the lmi method to the contrast problem, described in [START_REF] Bonnard | Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance[END_REF], in order to obtain upper bounds on the true contrast. Comparing these bounds to the contrast of our solutions then gives an insight about their global optimality. Table 3.2 shows the evolution of the upper bound on the contrast in function of lmi relaxation order, for the fluid case with t f = T min . As expected, the method yields a monotonically non-increasing sequence of sharper bounds. Relaxations of orders 4 and 5 yield very similar bounds, but this should not be interpreted as a termination criterion for the lmi method. Figs. 3.9 and 3.10 compare the tightest upper bounds found by the lmi method with the best candidate solutions found by Bocop and HamPath , in both the blood and fluid cases. The figures also represent the relative gap between the methods defined as (C LMI -C H )/C H , where C LMI is the lmi upper bound and C H is the contrast found with HamPath . As such, this measure characterizes the optimality gap between the methods. It does not, however, specify which of the method(s) could be further improved. At the fifth relaxation, the average gap is around 11% in the blood case, which, given the application, is satisfactory on the experimental level. For the fluid case, the average gap on the contrast is about 2% at the fifth relaxation, which strongly suggest that the solution is actually a global optimum. The gap is even below the 1% mark for t f ≤ 2 T min .

3.3.9 Numerical simulations, the multisaturation of two spins with B 1 -inhomogeneity.

In this section we give an illustration of our techniques applied to the saturation of two spins combining geometric analysis and numerical simulations to deduce the solution. We proceed in two steps.

• Step 1: Time minimal saturation of a single spin In the single-spin case the time minimal solution is described in Fig. 3.1 leads to construct the optimal solution for a continuation on the set of parameters where λ = 0 corresponds to the case of deoxygenated blood, λ = 1 corresponds to the case 2Γ = 3γ and λ = λ f is the water case: Γ = γ. According to Figs.3.11-3.12-3.13, due to the control bound, the bifurcation occurs not exactly at λ = 1 when the horizontal singular line z = γ/2/δ leaves the Bloch ball but at λ 0.99, since for λ > λ this line is no more accessible from the north pole.

• Step 2: We describe in Figs.3.14-3.15-3.16-3.17 the BC-extremal for the multisaturation problem with B 1 -inhomogeneity using the same continuation on the set of parameters. The control is computed using HamPath software in combination with Bocop in order to determine the structure of the extremal trajectory for λ = 0. Figs.3.14-3.15 show a control with the same structure δ + δ s δ + δ s δ + δ s , that is a sequence of three bang-singular arcs. A bifurcation occurs at λ 0.94 where the first singular arc disappears. Figs.3.16-3.17 show a control with a different structure δ + δ s δ + δ s . In each picture, we have represented the critical alti- tude z = γ/2/δ (on horizontal dotted line). At λ = λ f , the extremal is simply δ + δ 0 δ + δ 0 : singular arcs are obtained by applying a zero control.

Conclusion

The two cases studied in this book show the practical interest of combining geometric optimal control with numeric computations using the developed software to solve industrial type problems. The application to microswimmers is very recent and validate results obtained from fluid mechanics practitioners based on curvature control and Fourier analysis. The SR-geometry framework allows to compare different strokes and different swimmers, using the mechanical energy cost. The copepod mathematical swimmer is the simplest slender body model. Normal and abnormal strokes have interpretation in terms of sinusoidal and sequential paddlings. This leads to design a simple macroscopic copepod robot to validate the theoretical computations of the most efficient stroke. Another validation of the mathematical model using Resistive Force Theory for Stokes' flow is coming from the observations [START_REF] Lenz | Choreographed swimming of copepod nauplii[END_REF] showing the agreement between observed and predicted displacements. The mathematical developments lead to solve the inverse problem of identifying the cost used for the copepod nauplii displacement. The developments motivated by MRI are more profound and lead to intricate numerical investigations to deal with an highly complex optimal control problem with many local optimal solutions. Nevertheless we believe that the techniques validate by in vitro and in vivo experiments realized under the auspices of the ANR project DFG Explosy will find in a very near future applications in MRI diagnosis.

FIG. 1 . 2 Fig. 0. 1

 121 FIG.1. Sketch of the upper half of a swimmer paddling along the x axis, the line of symmetry.I. INTRODUCTION

Fig. 0. 2 Fig. 0. 3

 23 Fig. 0.2 Experimental results: The inner circle shape sample mimics the deoxygenated blood, the outside moon shape sample corresponds to the oxygenated blood. Left: Without control, Right: Optimized contrast.

Definition 31

 31 Parametrizing the normal geodesics solutions of → H n (z) and fixing x ∈ M, the exponential map is defined by exp x : (p,t) → Π (expt → H n (z)) where z = (x, p) and Π is the projection (x, p) → x.

2 LFig. 2 . 1

 221 Fig. 2.1 Purcell's 3-link swimmer.

FIG. 1 .

 1 FIG.1. Sketch of the upper half of a swimmer paddling along the x axis, the line of symmetry.I. INTRODUCTION

Fig. 2 . 3

 23 Fig. 2.3 Purcell stroke.

+ 4 p2 50 1/ 4 .

 504 k( p(0)) such that the linear terms of θ 1 (t) = x1 (ωt +ϕ 0 ), θ 2 (t) = x2 (ωt + ϕ 0 ) of (2.25) vanish, leads to periodic strokes with eight shapes of periodT = 4 K(k)/ p2 40The constants x10 , x20 are chosen such that ϑ (0) = 0. The initial adjoint vector p(0) has to verify the conditions H 1 ( x(0), p(0)) 2 + H 2 ( x(0), p(0)) 2 = 1, k( p(0)) ∈ (0, 1) and p2 40 + 4 p2 50 = 0. We integrate numerically the stroke in the physical variables starting from (q(0) = 0, p(0)) and show that the stroke has a conjugate point on [0, T ].• Rotating case: The modulus k can be expressed as k( p(0))

5 Fig. 2 . 4

 524 Fig. 2.4 Nilpotent approximation (normal case): state, adjoint variables and first conjugate point (blue cross), with the smallest singular value of the rank test.

Fig. 2 . 5

 25 Fig. 2.5 Computations of the complete elliptic integral K(k, ωt c ) and of the first conjugate point t 1c for normal strokes on the energy level H 21 + H 2 2 = 1. We observe: 0.3ωt 1c -0.4 < K(k) < 0.5ωt 1c -0.8.

2 Fig. 2 . 6

 226 Fig.2.6 Abnormal case: state variables for x(0) = (1, 0, 1, 0, 0), p(0) = (0, 0, -2, 1, 1).

Fig. 2 . 8

 28 Fig. 2.8 Optimal trajectory for |u| 2 (left) and the energy criterion (right) -displayed are the state variables x, y, α.

2 θFig. 2 . 9 2 Fig. 2 . 2 Fig. 2 .

 2292222 Fig. 2.9 Optimal trajectory for |u| 2 (left) and the energy criterion (right) -displayed are the state variables θ 1 , θ 2 .

Fig. 2 .

 2 Fig. 2.12 Optimal trajectory for |u| 2 (left) and the energy criterion (right) -displayed are the adjoint variables p x , p y and p α .

Fig. 2 .

 2 Fig. 2.13 Optimal trajectory for |u| 2 (left) and the energy criterion (right) -displayed are the adjoint variables p θ 1 , p θ 2 .

Fig. 2 .

 2 Fig. 2.14 (Left) State and adjoint variables for the Purcell swimmer minimizing the mechanical cost. (Right) Test of conjugate points (no conjugate point on [0, 2π]).

Fig. 2 .

 2 Fig. 2.15 (Left) State and adjoint variables for the Purcell swimmer minimizing the mechanical cost. (Right) Test of conjugate points (no conjugate point on [0, 2π]).

Fig. 2 .

 2 Fig. 2.16 (Left) State and adjoint variables for the Purcell swimmer minimizing the mechanical cost. (Right) Test of conjugate points. The cross on the trajectories on the left indicates the location of the first conjugate point.

Fig. 2 .+ u 2 2 ) 2 Fig. 2 .0 (u 2 1 + u 2 2 )

 222222 Fig.2.17 Continuation on the amplitude: x(T ) 2 + y(T ) 2 = c 1 for the T 0 (u 2 1 + u 2 2 ) dt cost.

1 Fig. 2 .

 12 Fig. 2.19 Efficiency curves for the two families of strokes presented in Fig.2.18.

Fig. 2 .

 2 Fig. 2.21 Closed periodic planar curves: non intersecting curve, eight curve and limac ¸on curve.

Fig. 2 .Fig. 2 .Fig. 2 .

 222 Fig. 2.22 (Left) Normal stroke where the constraints are satisfied: simple loop with no conjugate point on [0, T ]. (Right) Limac ¸on with inner loop with one conjugate point on [0, T ].

Fig. 2 .

 2 Fig. 2.25 One parameter family of circles which are the geodesics of the Heisenberg-Brockett problem.

Fig. 2 .

 2 Fig. 2.26 One parameter family of simple loops, limac ¸ons and Bernoulli lemniscates normal strokes.

Fig. 2 .

 2 Fig. 2.27 One parameter family of simple loops symmetric with respect to the straight line D : θ 2 = -θ 1 + π with converges to a point when the displacement tends to 0.

. 1 (

 1 Fig. 3.1 (left) Time minimal solution compared with (right) inversion sequence .

3. 3 . 5

 35 Stratification of the surface Σ : H 1 = H 2 = 0 and partial classification of the extremal flow near Σ

Lemma 8

 8 Hence the equation θ = 0 leads to the relation A sin(θ + φ ) = b, which has two distinct solutions on [0, 2π[ denoted θ 0 , θ 1 if and only if A > |b|, one solution if A = |b| and zero solution if |A| < |b|. Moreover θ 1 -θ 0 = π if and only if b = 0. Plugging θ 0 , θ 1 in (3.25), one gets two solutions of (3.25). Hence we deduce: If a 2 1 + a 2 2 > |b| and b = 0, we have a broken extremal formed by concatenating two extremals of order zero at each point z 0 of Σ 1 .

  use the notation D = {{H G , H F }, H G } and D = {{H G , H F }, H F }. The differential equation (3.36) can be desingularized using the time reparametrization ds = dt/D(z(t)) which amounts to analyze the one dimensional foliation.

P 1 :

 1 Fluid case. Spin 1: Cerebrospinal fluid: T 1 = 2, T 2 = 0.3; Spin 2: Water: T 1 = 2.5 = T 2 . T min = 26.17040. P 2 : Blood case. Spin 1: Deoxygenated blood: T 1 = 1.35, T 2 = 0.05; Spin 2: Oxygenated blood: T 1 = 1.35, T 2 = 0.2.T min = 6.7981.

1 Fig. 3 . 2

 132 Fig. 3.2 Fluid case: Two local solutions for λ = 2.0. Trajectories for spin 1 and 2 in the (y,z)plane are portrayed in the first two subgraphs of each subplot. The corresponding control is drawn in the bottom subgraph. The two bang arcs have the same sign for the left solution, whereas for the right solution, the two bang arcs are of opposite sign.

Fig. 3 . 3

 33 Fig.3.3 Fluid case: second order conditions. Second order necessary condition checked on the best solution with λ = 2.0 from subsection 3.3.8. The rank condition from the algorithm presented in subsection 4.3 from[START_REF] Bonnard | Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance[END_REF] is evaluated along the two singular arcs. See[START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] for details on the concept of conjugate times. On the left subplot, for each singular arc, the curve is reparameterized so that the final time corresponds to the abscissa 1 (vertical blue dashed line); the determinant associated with the rank condition is plotted (top subgraph), so there is a conjugate time whenever it vanishes (vertical red dashed lines). One observes that conjugate times on each arc are located after the (normalized to 1) final time, satisfying necessary condition of local optimality of the trajectory. At the bottom, the smallest singular value of the matrix whose rank we test is plotted, extracting only the relevant information to detect the rank drops. On the right subplot is presented a zoom of top-left subgraph near the two conjugate times.

Fig. 3 .Fig. 3 . 4

 334 Fig.3.4 Fluid case: influence of the final time. On the left subgraph are shown the control laws of solutions at λ = 2 and λ = 1.32 from path from the right subplot. For λ = 1.32, we can see the saturating singular arc around the normalized time τ = 0.92 (the time is normalized to be between 0 and 1 for each solution). The 2BS solution at λ = 1.32 is used to initialize a multiple shooting with a 3BS structure and then to perform a new homotopy from λ = 1.32 to λ = 1. On the right subgraph is portrayed the two homotopies: the first from λ = 2 to λ = 1.32 and the second to λ = 1, with one more BS sequence. The value function, the norm of the initial adjoint vector, the norm of the shooting function and the switching times along the path are given. The blue color represents 2BS solutions while the red color is for 3BS structures. The dashed red lines come from the extended path after the change of structure detected around λ = 1.32.

1 Fig. 3 . 5

 135 Fig. 3.5 Fluid case: influence of the final time.Crossing between two branches with 3BS solutions. The crossing is around λ = 1.0484, see top subgraph. Thus for λ ≤ 1.0484, the best solution, locally, has a 3BS structure of the form δ + δ s δ + δ s δ -δ s (bottom-left subgraph) while for λ ∈ [1.0484, 1.351] the best solution is of the form δ + δ s δ -δ s δ -δ s (bottom-right subgraph). On the two bottom subgraphs, the trajectories for spin 1 and 2 in the (y,z)-plane are portrayed with the corresponding control, both for λ = 1.0484.

Fig. 3 . 6

 36 Fig. 3.6 Fluid case: influence of the final time. Second order necessary condition checked along the extended path from Fig. 3.4. For all solutions from λ = 1 to λ = 3 are computed the first conjugate times along each singular arc. For λ ∈ [1, 1.32], the structure is 3BS and there are 3 singular arcs. For λ ∈ [1.32, 3], there are 2 singular arcs. Each singular interval is normalized in such a way the initial time is 0 and the final time is 1. The lower dashed horizontal line represents the final time 1.There is no conjugate time before the normalized final time 1 which means that all solutions satisfy the second order necessary condition of local optimality. Note that at a magenta cross, around (1.32, 1), the control of the first singular arc saturates the constraint |u| = 1, and so no conjugate time is computed after this time.

Fig. 3 . 9

 39 Fig.3.9 Fluid case. Best upper bounds (dashed line) by the lmi method compared with best solutions by HamPath (solid line), and relative gap between the two.

Fig. 3 .

 3 Fig. 3.10 Blood case. Best upper bounds (dashed line) by the lmi method compared with best solutions by HamPath (solid line), and relative gap between the two.

Fig. 3 .

 3 Fig. 3.11 Time-minimal saturation of a single spin (λ = 0).

Fig. 3 .

 3 Fig. 3.12 Time-minimal saturation of a single spin (λ = 0.9941 λ ).

  Fig. 3.13 Time-minimal saturation of a single spin (λ = λ f ).

Fig. 3 .Fig. 3 .

 33 Fig. 3.14 BC-extremal for the multisaturation problem with λ = 0.

3 ×Fig. 3 .

 33 Fig. 3.16 BC-extremal for the multisaturation problem with λ = 0.948 > λ .

Fig. 3 .

 3 Fig. 3.17 BC-extremal for the multisaturation problem with λ = λ f .

  101

  5. The Jacobi identity is satisfied:{{F, G}, H} + {{G, H}, F} + {{H, F}, G} = 0.

	Definition 14 Let say that F is a first integral for -→ H be a Hamiltonian vector field on (M, ω) and F : M → R. We -→ H if F is constant along any trajectory of -→ H , that is dF( -→ H ) = {F, H} = 0. Definition 15 Let (x, p) be a Darboux coordinate system and H : M → R. The co-ordinate x 1 is said to be cyclic if ∂ H ∂ x 1 = 0. Hence F : (x, p) → p 1 is a first integral.
	Definition 16 Let M be a n-dimensional manifold and let (x, p) be Darboux coordi-nates on T * M. For any vector field X on M we can define a Hamiltonian vector field -→ H X on T * M by H(x, p) = p, X(x) ; -→ H X is called the Hamiltonian lift of X and -

  2 40 + 4 p 2 50 ) 1/4 , and the complete elliptic integral K(k), where k is the amplitude

	k =	1 2	2 p2 40 + 4 p2 50 + p2 30 -2 p10 p40 -4 p50 p20 -4 p50 p40 x30 p2 40 + 4 p2

50

.

  ).

	x 10 4					det
	10
	8
	6
	4
	2
	0 0	2	4	6	8	10	12	14	16	18	20
						t
						svd
	1
	0.8
	0.6
	0.4
	0.2
	0 0	2	4	6	8	10	12	14	16	18	20
						t
	Fig. 2.7 Abnormal case: the second order sufficient condition indicates there is no con-jugate point.

Table 3 .

 3 1 Fluid case: Batch optimizations (Direct method). For each value of λ we test the three initializations for the control u, and record the value of the objective (i.e. the contrast), as well as the control structure (i.e. the signs of bang arcs). CPU times for a single optimization are less than one minute on a Intel Xeon 3.2GHz.

  Fluid case, sub-optimal synthesis. Illustration on the left subplot, of local solutions (each branch corresponds to a control structure). The suboptimal synthesis is plotted on right subplot. The colors are blue for 2BS structure, red for 3BS and green for 4BS. The best policy is δ + δ s δ + δ s δ -δ s for λ ≤ 1.0484, and δ + δ s δ -δ s δ -δ s for λ ∈ [1.0484, 1.351]. Then, for λ ∈ [1.351, 3], the best policy is 2BS and of the form δ + δ s δ -δ s .[1, 2]T min . The three black branches are made only of BS solutions whereas the two others are made of 2BS and 3BS solutions. To maximize the contrast, the best policy, drawn as solid lines, is:δ + δ s δ + δ s for λ ∈ [1.000, 1 + ε], ε > 0 small δ + δ s for λ ∈ [1 + ε, 1.294], δ + δ s δ -δ s δ -δ s for λ ∈ [1.294, 2.000].

		0.7							0.7
		0.68							0.68
		0.66							0.66
	contrast	0.64						contrast	0.64
		0.62							0.62
		0.6							0.6
		0.58							0.58
		1	1.5		2	2.5	3		1	1.5	2	2.5	3
					λ						λ
	Fig. 3.7 (3.38)
		0.5						0.5
		0.48						0.48
								0.46
	Contrast	0.44 0.46					Contrast	0.44
								0.42
		0.42						0.4
		1 0.4	1.2	1.4	1.6	1.8	2	1 0.38	1.5	2		2.5	3
				λ						λ

Table 3 .

 3 2 Fluid case, t f = T min : upper bounds on contrast -J r M , numbers of moments N r and CPU times t r in function of relaxation order r.

	r 1 0.8474 63 0.7 -J r M N r t r 2 0.7552 378 3 3 0.6226 1386 14 4 0.6069 3861 332 5 0.6040 9009 8400

Historical part -Calculus of variations

Weak Maximum Principle and Application to Swimming at low Reynolds Number

Acknowledgements

supported by the ANR Project -DFG Explosys (Grant No. ANR-14-CE35-0013-01;GL203/9-1). M. Chyba is partially supported by the Simons Foundation with Collaboration Grant for Mathematicians: Geometric Optimal Control and Its Application. J. Rouot is supported by the European Research Council (ERC) through an ERC-Advanced Grant for the TAMING project.

• This book benefited from the support of the FMJH PGMO and from the support of EDF, Thales and Orange.

and (y 1 , r 2 (y 1 , y 2 ), y 2 , r 2 (y 1 , y 2 )) (defined for (y 1 , y 2 ) = (0, 0)) vanishes both D and D . Finally, if y 1 = y 2 = 0, we have the solution (0, z, 0, z), z ∈ R. We summarize the case z 1 = z 2 in:

Lemma 12 {D = 0} ∩ {D = 0} ∩ {z 1 = z 2 } is the union of an affine plane z 1 = z 2 = z γ,Γ , a rational surface z 1 = z 2 = r 2 (y 1 , y 2 ) (defined for (y 1 , y 2 ) = (0, 0)), and the line {(0, z, 0, z), z ∈ R}.

-Assume first that y 1 = 0 and z 1 = z 2 . We have z 1 y 2 = 0.

• If y 1 = z 1 = 0, then:

We have:

Observe that the polynomial d 1 vanishes if and only if z 1 equals zγ,Γ with

and in this case, there is no solution such that z 2 = z 1 . Provided d 1 (z 1 ) = 0, one gets z 2 = r 1 (z 1 ):

which is a rational function of z 1 . And the intersection with

is the union of two lines of {y 1 = z 1 = 0} intersecting at N and a rational curve

We can eliminate z 2 using:

and, substituting in y 2 1 D we get the factorization y 2 1 D = q 1 q 2 , with:

Provided d 1 = 0 (that is z 1 = zγ,Γ ), we solve q 2 = 0 with respect to y 2 , and then we get the value of (y 2 , z 2 ):

.

We put d 3

Its discriminant with respect to y 1 is:

2 and its sign changes when z 1 reaches zγ,Γ . Provided d 3 (y 1 , z 1 ) = 0, we solve q 1 with respect to y 2 , and then we get the value of (y 2 , z 2 ):

.

• Analysis of the behaviors of the solutions of X r e near O. We set zi = 1 + z i and we have the following approximations: