# Sub-Riemannian geometry and swimming at low Reynolds number: the Copepod case

Piernicola Bettiol, Bernard Bonnard, Alice Nolot, Jérémy Rouot

► To cite this version:

Piernicola Bettiol, Bernard Bonnard, Alice Nolot, Jérémy Rouot. Sub-Riemannian geometry and swimming at low Reynolds number: the Copepod case. ESAIM: Control, Optimisation and Calculus of Variations, 2019, 25 (9), 10.1051/cocv/2017071. hal-01442880v3

HAL Id: hal-01442880

https://inria.hal.science/hal-01442880v3

Submitted on 23 Nov 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques 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.
# SUB-RIEMANNIAN GEOMETRY AND SWIMMING AT LOW REYNOLDS NUMBER: THE COPEPOD CASE*

P. BETTIOL<fnref n="1" />, B. BONNARD<fnref n="2" />, A. NOLOT<fnref n="3" /> AND J. ROUOT<fnref n="4" />

**Abstract.** In [22], based on copepod observations, Takagi proposed a model to interpret the swimming behaviour of these microorganisms using sinusoidal paddling or sequential paddling followed by a recovery stroke in unison, and compares them invoking the concept of efficiency. Our aim is to provide an interpretation of Takagi's results in the frame of optimal control theory and sub-Riemannian geometry. The Maximum principle is used to select two types of periodic control candidates as minimizers: sinusoidal up to time reparameterization and the sequential paddling, interpreted as an abnormal stroke in sub-Riemannian geometry. Geometric analysis combined with numerical simulations are decisive tools to compute the optimal solutions, refining Takagi computations. A family of simple strokes with small amplitudes emanating from a center is characterized as an invariant of SR-geometry and allows to identify the metric used by the swimmer. The notion of efficiency is discussed in detail and related with normality properties of minimizers.

**Résumé.** Dans [22], à partir de l'observation du mécanisme de nage d'une famille de microorganismes appelés copépodes, Takagi propose un modèle pour interpréter ces nages. Deux types de contrôles périodiques associés sont proposés : le premier correspond à des contrôles sinusoïdaux formant une courbe de nage lisse et simple et le second est un contrôle constant par morceaux produisant un triangle. Notre objectif est d'interpréter cela dans le cadre du contrôle optimal et de la géométrie sous-Riemannienne. Le principe du Maximum est utilisé pour sélectionner des nages géodésiques de deux types : des nages formant des courbes simples associées à des contrôles sinusoïdaux à une reparamétrisation du temps près et le triangle est interprété comme une nage anormale. L'analyse géométrique combinée avec des simulations numériques permet de générer une famille de nages de petites amplitudes ce qui par continuation permet de calculer la nage la plus efficace. La notion d'efficacité est discutée en détails en relation avec le concept de normalité.

**1991 Mathematics Subject Classification.** 70Q05, 93C10, 49K15.

November 22, 2017.

<footnote>Keywords and phrases: Stokes flow, Optimal control theory, sub-Riemannian geometry, Abnormal closed geodesics</footnote>

<footnote>* The second author is supported by the ANR-DFG research program Explosys (Grant No. ANR-14-CE35-0013-01; GL203/9-1) and the fourth author is supported by the European Research Council (ERC) through an ERC-Advanced Grant for the TAMING project.</footnote>

<footnote n="1">¹ Laboratoire de Mathématiques Unité CNRS UMR 6205, Université de Bretagne Occidentale, 6, Avenue Victor Le Gorgeu, 29200 Brest, France</footnote>

<footnote n="2">² Institut de Mathématiques de Bourgogne et Inria Sophia Antipolis, 9 avenue Savary, 21078 Dijon, France</footnote>

<footnote n="3">³ Institut de Mathématiques de Bourgogne, 9 avenue Savary, 21078 Dijon, France</footnote>

<footnote n="4">⁴ LAAS-CNRS, Université de Toulouse, LAAS, 31031 Toulouse cedex 4, France</footnote>
# 1. INTRODUCTION

This paper combines mathematical tools coming from both optimal control theory and sub-Riemannian (SR-) geometry to investigate the copepod swimming model suggested by Takagi in a recent paper [22]. Copepod is a widespread zooplankton species, living in a low Reynolds number environment. (Fig.1 (left) illustrates a copepod picture together with its characteristic dimension). Takagi's swimming model is probably the simplest one based on *slender body theory for Stokes flow*: it has two symmetric links, and the body is reduced to an infinitesimal sphere (see Fig.1 (right)).

<figure><img src="image_1.png" /><figcaption>FIGURE 1. (left) Observation of a zooplankton. (right) Sketch of the 2-link symmetric swimmer.</figcaption></figure>

It represents an excellent case study to understand the swimming mechanism and might be helpful to investigate more complicated models: for instance, Purcell swimmers [17, 19] or swimmer with several pairs of symmetric links.

This article completes along different directions preliminary results obtained on the copepod swimmer in [6], where this model was briefly studied using mainly numerical simulations, the objective being to understand particular properties of optimal strokes carrying out the analysis up to second order sufficient optimality conditions: in that context the focus was mainly on the Purcell swimmer (which has two symmetries), whilst the copepod, being simpler and having just one symmetry degree, represents a sort of intermediate model to compare with the Purcell swimmer.

To produce the displacement along a reference line $0x_0$, we use a pair of two symmetric links with equal length normalized to $l = 1$, $\theta_1, \theta_2$ are the respective angles of the two links and they satisfy the triangle constraint $\mathcal{T}: 0 \le \theta_1 \le \theta_2 \le \pi$.

It will appear that the constraint is not active in our study, since the optimal stroke does not saturate the constraint. Moreover, the model is not valid in the whole triangle, due to the interaction between links and is investigated by the analytic continuation method. A practical problem on the experiments is to check that the optimal stroke stays in the physical domain of validity of the model.

Stokes theory relates the speed of the displacement variable $x_0$ to the speed of the shape vector $\theta = (\theta_1, \theta_2)$ by the equation

$$ \dot{x}_0 = \sum_{i=1}^{2} \frac{\dot{\theta}_i \sin \theta_i}{\Delta(\theta)}, \quad \text{where } \Delta(\theta) = \sum_{i=1}^{2} (1 + \sin^2 \theta_i). \qquad (1) $$

To parameterize the motion as a control system, one introduces the dynamics

$$ \dot{\theta}_1 = u_1, \quad \dot{\theta}_2 = u_2. $$

It provides a control system written as

$$ \dot{q} = u_1 F_1(q) + u_2 F_2(q) \quad q = (x_0, \theta) \qquad (2) $$
with state constraints $\mathcal{T}: 0 \le \theta_1 \le \theta_2 \le \pi$. A stroke $\gamma$ is a closed curve in the $\theta$-plane produced by a periodic control $u(t)$, $u = (u_1, u_2)$.

To compare different strokes or different swimmers one introduces a metric in the shape variable. A choice of interest, in particular in relation with the concept of efficiency defined by [15], is the mechanical energy dissipated by drag forces:

$$E = \int_{0}^{T} \dot{q}^{\top} M(q) \dot{q} dt \qquad (3)$$

where $T$ is the period of the stroke and $M$ is the matrix

$$M = \begin{pmatrix} 2 - 1/2(\cos^2 \theta_1 + \cos^2 \theta_2) & -1/2 \sin \theta_1 & -1/2 \sin \theta_2 \\ -1/2 \sin \theta_1 & 1/3 & 0 \\ -1/2 \sin \theta_2 & 0 & 1/3 \end{pmatrix} \qquad (4)$$

Using (1) this amounts to minimize the quadratic form

$$E = \int_{0}^{T} (a(q)u_{1}^{2} + 2b(q)u_{1}u_{2} + c(q)u_{2}^{2}) dt$$

where

$$a(q) = \frac{1}{3} - \frac{\sin^2 \theta_1}{2\Delta(\theta)}, \quad b(q) = \frac{-\sin \theta_1 \sin \theta_2}{2\Delta(\theta)}, \quad c(q) = \frac{1}{3} - \frac{\sin^2 \theta_2}{2\Delta(\theta)}.$$

This model is similar to the model observed in [4] as a limit case of a symmetric Purcell swimmer but was analyzed with different mathematical tools.

For geometric and numeric computations using continuation methods, a simplified metric is introduced: this is the Euclidean (flat) metric associated with the quadratic cost

$$\int_{0}^{T} (u_{1}^{2} + u_{2}^{2}) dt.$$

We underline the fact that the outcome of the qualitative analysis does not change in its main features, where we vary the parameters of the metric involved in the problem. We also observe that it is consistent with observation of the zooplankton behavior [14].

Denoting by $g$ such a metric, the length of a stroke $\gamma = (\theta_1, \theta_2)$ is

$$l(\gamma) = \int_{0}^{T} g(\theta, \dot{\theta})^{1/2} dt$$

and our optimal control problem is rewritten, owing to Maupertius principle, as a SR-problem:

$$\dot{q} = \sum_{i=1}^{2} u_i F_i(q), \quad \min_{u(.)} \int_{0}^{T} g(\theta, \dot{\theta}) dt$$

with appropriate boundary conditions associated with periodic controls

$$\theta(0) = \theta(T),$$

and with state constraints $\mathcal{T}: 0 \le \theta_1 \le \theta_2 \le \pi$. $x_0(T) - x_0(0)$ stands for the displacement of a stroke and we can set $x_0(0) = 0$. Using the energy minimization point of view, the period of a stroke can be fixed to $T = 2\pi$. Note also that the problem is equivalent to a time minimal control problem by parameterizing the arc-length (which experimentally will allow to compare different strokes or swimmers).
From the general point of view of SR-geometry, the problem consists in computing for each $q(0) = q_0$, the SR-sphere $S_{q_0}(r)$ with radius $r$, formed by extremities of minimizers starting from $q_0$, having fixed length and moreover requiring that the optimal control is periodic. Note that not every radius is producing a minimizer and the main problem is to compute radii (see Fig.2).

This is equivalent to fix the displacement $x_0(2\pi)$ ($x_0(0) = 0$) and to compute strokes minimizing the length. In this geometric frame one can also introduce the concept of geometric efficiency of a stroke $\gamma$ that is the ratio

$$ \mathcal{E} = x_0(2\pi)/l(\gamma) \qquad (5) $$

which turns out to be equivalent to the concept of efficiency in fluid mechanics used in [14, 22] and defined by $\mathcal{E}' = x_0(T)^2/E(\gamma)$.

Note that more general concepts of efficiency can be used as in [11]: Eff $\sim x_0(T)^2 M_{11}(\gamma(0))/E(\gamma)$, with $M = (M_{ij})_{ij}$ defined in (4), which depends on $\theta(0)$ and more generally any function $h(\theta(0), x_0(T), E(\gamma))$.

The main contribution of this article is that in the optimal control theory frame, one can construct a one parameter family of geodesics strokes for the metric defined by $g$, providing a foliation of the interior of the triangle $\mathcal{T}$ whose boundary corresponds to the abnormal geodesic stroke, see Fig.2. A unique stroke of this family is the solution which maximizes the efficiency. A whole section is devoted to the analysis of the optimal status of the abnormal triangle stroke, completing the explicit evaluation of the efficiency, outlined in [6].

<figure><img src="image_2.png" /><figcaption>FIGURE 2. One parameter family of geodesic strokes for the Euclidean cost (left) and for the mechanical cost (3) (right).</figcaption></figure>

This construction is based on the application of the Maximum Principle [18, 23] supplemented by second order optimality conditions [7] and such curves are numerically evaluated using continuation and shooting methods implemented in the HamPath software [12]. The numerical simulations are performed by continuation from strokes with small amplitudes. These differ from those ones appearing in [6] where the continuation is performed using initialization from strokes with large amplitudes near the boundary of the physical domain (and using the Bocop toolbox). With the HamPath code, we can compute second order necessary conditions corresponding to the concept of conjugate points.

Our family of geodesics is obtained starting from strokes of small amplitudes and having a center $C$, and can be evaluated by constructing a graded normal form in SR-geometry [3, 5]. Asymptotics of the corresponding displacement and its length can be also computed. Further asymptotics can be computed for larger amplitudes near the boundary of the triangle formed by the sequential paddling followed by the recovery stroke in unison. An interpretation of this triangle stroke is given in the frame of SR-optimality.
SR-geometry in relation with optimal control was introduced in the seminal article [10], which (fortunately) already contains the two crucial ingredients used in our article: the Heisenberg-Dido model which can be interpreted as a *rough* swimmer model, and normal forms to evaluate the SR-distance and spheres with small radius. This boils down to estimate conjugate and cut loci. Such computations were also detailed in [1] for a generalization of the Dido problem, in connection with the analysis of a similar system (which concerns the estimation of geodesics for a particle in a magnetic field) and are applicable to our study. The relation of this problem with microswimming was already pointed in the earliest reference [16].

The article is organized in three sections. The first section presents a geometric analysis of the copepod swimmer in the frame of sub-Riemannian geometry to compute geodesic stokes with small amplitudes. Near a point interior to the triangle $\mathcal{T}$, they are obtained by perturbation of small circles from the Brockett-Heisenberg nilpotent model and simple loops and limaçons can appear. An algorithm is described to determine the centre of swimming observed in Fig.2, and deduced from the equivalence between the model of order $-1$ (nilpotent model) and the model of order 0. At a point on the edges of the triangle, sub-Riemannian study predicts the existence of additional eight shape strokes. In section 3, the combination of the Maximum Principle, second-order optimality conditions and numerical simulations allow to compute strokes maximizing different concepts of efficiency. A final section is a mathematical analysis relating efficiency and normality for the copepod model. This provides a supplementary analytical argument to confirm the numerical simulations displayed in the previous section, underlying the fact that the triangular abnormal stroke is not efficient.

## 2. SR-GEOMETRY PROPERTIES OF THE COPEPOD MODEL AND STROKES WITH SMALL AMPLITUDES

### 2.1. Geodesics computation

Geodesic equations are presented here for the copepod swimmer, but can be generalized to any (swimmer) model having the following variables decomposition: $q = (q', q'')$, where $q'$ represents the displacement variable, and $q''$ stands for the shape variable. The dimension of the space in which $q''$ is taken provides the dimension of the nonholonomic constraint, and corresponds to the number of the controlled swimmer links: for instance $q'' \in \mathbb{R}^2$ for the copepod model here considered and for the Purcell swimmer (cf [19] and [20]); $q'' \in \mathbb{R}^l$ with $l \ge 3$ for more sophisticated copepod models (cf [22]).

For the copepod control system here investigated, if we consider the energy minimization criterion, then the reference optimal control problem becomes

$$ \min \int_{0}^{T} L(q(t), u(t))dt, $$

over arcs $q(.)$ on $[0, T]$ s.t. $\dot{q}(t) = F(q(t), u(t)) = \sum_{i=1}^{2} u_i F_i(q)$ a.e. . The class of admissible controls is the set $\mathcal{U}$ of bounded measurable mapping taking values in $U = \mathbb{R}^2$. The boundary conditions are:

$$ \begin{aligned} q_1(0) &= 0, & q_1(T) &= x_T > 0 \\ q_2(0) &= q_2(T) & q_3(0) &= q_3(T) \end{aligned} $$

where $T$ can be fixed to $2\pi$. The dynamics are

$$ F(q, u) = \left( u_1 \frac{\sin q_2}{2 + \sin^2 q_2 + \sin^2 q_3} + u_2 \frac{\sin q_3}{2 + \sin^2 q_2 + \sin^2 q_3}, u_1, u_2 \right). \quad (6) $$

Here, the state variable is $q = (x_0, \theta)$ where $\theta = (\theta_1, \theta_2)$ represent the link angles and $x_0$ stands for the displacement. The cost is defined by a Riemannian metric on the shape variables: $L(\theta, u_1, u_2) = a(\theta)u_1^2 + 2b(\theta)u_1u_2 + c(\theta)u_2^2$ which can be either the mechanical energy or the Euclidean metric.
**Maximum principle in SR-geometry** [18]. We introduce the pseudo-Hamiltonian $H(q, p, u) = p\dot{q}+p_0(a(q)u_1^2+2b(q)u_1u_2+c(q)u_2^2)$, where $p_0$ is a constant. According to the Maximum principle, minimizers are found among extremal curves, which are solutions a.e. of the following equations

$$ \dot{q} = \frac{\partial H}{\partial p}, \quad \dot{p} = -\frac{\partial H}{\partial q}, \quad \frac{\partial H}{\partial u} = 0. $$

This leads to consider the following two fundamental cases.

**Normal case.** Assume $p_0 \neq 0$ and it can be normalized to $p_0 = -1/2$. Controls, solving $\partial H/\partial u = 0$, are given by

$$ u_1 = \frac{cH_1 - bH_2}{ac - b^2}, \quad u_2 = \frac{aH_2 - bH_1}{ac - b^2} $$

where $H_i = p_i \cdot F_i$ denotes $F_i$'s Hamiltonian lifts and plugging such controls in $H$ yields the normal Hamiltonian

$$ H_n = \frac{1}{2} \frac{c H_1^2 - 2b H_1 H_2 + a H_2^2}{ac - b^2} \quad (7) $$

The corresponding solutions $z = (q, p)$ are called normal extremals.

A remarkable symmetry property is the following.

**Lemma 2.1.** For the Euclidean case and the mechanical energy case, the geodesic flow defined by (7) is invariant under the transformation $\sigma : (x_0, \theta_1, \theta_2) \mapsto (x_0, \pi - \theta_2, \pi - \theta_1)$.

**Abnormal case.** If $p_0 = 0$, additional extremals $z = (q, p)$ appear and they are called abnormal. They are solutions of the implicit equations

$$ H_1(z(t)) = H_2(z(t)) = 0 \quad (8) $$

and they can be computed by (time) differentiation.

Recall that if $F, G$ are two smooth vector fields, the Lie bracket is computed as:

$$ [F, G](q) = \frac{\partial F}{\partial q}(q)G(q) - \frac{\partial G}{\partial q}(q)F(q). $$

If $H_F(z) = p \cdot F(q)$, $H_G(z) = p \cdot G(q)$, the Poisson bracket is

$$ \{H_F, H_G\}(z) = dH_F(\vec{H}_G) = p \cdot [F, G](q) $$

with $\vec{H}_G = \frac{\partial H_G}{\partial p} \frac{\partial}{\partial q} - \frac{\partial H_G}{\partial q} \frac{\partial}{\partial p}$.

Hence differentiating twice (8) and using $\dot{z} = u_1 \vec{H}_F(z) + u_2 \vec{H}_G(z)$, the abnormal controls are given by

$$ \begin{aligned} H_1 &= H_2 = \{H_1, H_2\} = 0 \\ u_1 \{\{H_1, H_2\}, H_1\} + u_2 \{\{H_1, H_2\}, H_2\} &= 0 \end{aligned} \quad (9) $$

and they can be (generically) computed using (9) provided that one Poisson bracket $\{\{H_1, H_2\}, H_i\}$, $i = 1, 2$, is non zero.

**Definition 2.2.** The exponential mapping is, for fixed $q_0 = q(0)$ the map: $\exp_{q_0}(t, p(0)) \mapsto \Pi(\exp(t \vec{H}_n(z(0)))$ where $\Pi$ is the projection: $(q, p) \to q$. A projection of an extremal is called a geodesic. It is called strictly normal if it is the projection of a normal extremal but not an abnormal one. A time $t_c$ is a conjugate time if the exponential mapping is not of full rank at $t_c$. Denoting by $t_{1c}$ the first conjugate time, $q(t_{1c})$ is called the first conjugate point along the reference geodesic $t \mapsto q(t)$.
We recall also the following facts. Fixing $q_0$, the wave front $W(q_0, r)$ is the set of extremities of geodesics (normal or abnormal) with length $r$ and the sphere $S(q_0, r)$ is the set of extremities of minimizing geodesics. The conjugate locus $C(q_0)$ is the set of first conjugate points of normal geodesics starting from $q_0$ and the cut locus $C_{cut}(q_0)$ is the set of points where geodesics cease to be optimal.

**Definition 2.3.** According to the previous definitions, a stroke is called (strictly) normal if it is a (strictly) normal geodesic with periodic control while an abnormal stroke is a piecewise smooth abnormal geodesic with periodic control.

## 2.2. Computation for the copepod case

One has

$$F_i = \frac{\sin(\theta_i)}{\Delta(\theta)} \frac{\partial}{\partial x_0} + u_i \frac{\partial}{\partial \theta_i}, \quad i = 1, 2.$$

with $\Delta(\theta) = \sum_{i=1}^{2} (1 + \sin^2(\theta_i))$. We get

$$[F_1, F_2](q) = \tilde{f}(\theta_1, \theta_2) \frac{\partial}{\partial x_0}$$

with

$$\tilde{f}(\theta_1, \theta_2) = 2 \frac{\sin(\theta_1) \sin(\theta_2) (\cos(\theta_1) - \cos(\theta_2))}{\Delta^2(\theta)}. \quad (10)$$

Furthermore,

$$[[F_1, F_2], F_i] = \frac{\partial \tilde{f}}{\partial \theta_i} \frac{\partial}{\partial x_0}, \quad i = 1, 2$$

and we have simple formulas to generate all Lie brackets.

**Definition 2.4.** A point $q_0$ is called a Darboux or a contact point if at $q_0$ $F_1, F_2$ and $[F_1, F_2]$ are linearly independent and a Martinet point if $F_1, F_2$ are linearly independent, $[F_1, F_2] \in \text{span}\{F_1, F_2\}$ but at least for one $i = 1, 2$, $[[F_1, F_2], F_i] \notin \text{span}\{F_1, F_2\}$.

According to this terminology and Lie brackets computations, we have

**Proposition 2.5.** (1) All interior points of the triangle $\mathcal{T}$: $0 \le \theta_1 \le \theta_2 \le \pi$ are contact points.

(2) The boundary of the triangle $\mathcal{T}$ represents the abnormal (piecewise smooth) stroke, and each point, vertices excluded, is a Martinet point.

**Geometric comment.** Hence the observed stroke by Takagi [22] of sequential paddling followed by a recovery stroke in unison corresponding to the policy: $\theta_2: 0 \to \pi$, $\theta_1: 0 \to \pi$, $\theta_i: \pi \to 0$, $i = 1, 2$ with $\theta_1 = \theta_2$ where the copepod swimmer followed the triangle $\mathcal{T}$, boundary of the physical domain, is the unique abnormal stroke.

The controlled dynamics (6) allow a suggestive interpretation of the copepod displacement in terms of Stokes' Theorem and curvature control methods.

**Proposition 2.6.** Take the smooth one-form on $\mathbb{R}^2$: $\omega := \sum_{i=1}^{2} \frac{\sin \theta_i}{\Delta(\theta)} d\theta_i$.

Then

(1) $d\omega = -\tilde{f}(\theta_1, \theta_2)d\theta_1 \wedge d\theta_2$ (where $\tilde{f}$ is defined in (10)).

(2) for any bounded Stokes domain $D \subset \mathbb{R}^2$, we have

$$\oint_{\partial D} \sum_{i=1}^{2} \frac{\sin(\theta_i)}{\Delta} d\theta_i = \int_{D} -\frac{\tilde{f}(\theta_1, \theta_2)}{\Delta^2} d\theta_1 \wedge d\theta_2 = \int_{D} d\omega.$$

(3) $d\omega < 0$ in the interior of the triangle $\mathcal{T}$, and $d\omega$ vanishes on the boundary of $\mathcal{T}$.
**Corollary 2.7.** (1) Consider a piecewise smooth stroke $\gamma(\cdot)$, and let $D$ be the bounded Stokes domain such that $\gamma = \partial D$. Then the stroke displacement is

$$x_0(T) = \int_D d\omega.$$

(2) The piecewise smooth stroke confined in $\mathcal{T}$ having the biggest displacement is the abnormal stroke.

## 2.3. SR-classification in dimension 3 and strokes with small amplitudes for the copepod swimmer

In the geometric approach of SR-geometry, the problem of computing strokes with small amplitudes is precisely the problem of computing the points of the SR-sphere with small radius associated with periodic normal strokes. The geometry is related to the action of the group $\mathcal{G}$ preserving the distribution and acting on the metric defined on the shape space. The key notion for computing such spheres is the concept of *privileged coordinates*, see [5, 13], forming near each point a set of graded coordinates. The concept of nilpotent approximation or model of order $-1$ allows to get estimates of the SR-balls with small radii and more precise approximations have to be used, especially to evaluate the conjugate locus. In particular, this notion was developed in [2, 3] and [8, 9] respectively in the contact and Martinet case. This kind of computations are recalled next, and are crucial in our analysis.

### 2.3.1. The contact case (see [3] for concepts and details of computations).

The SR-problem is written as $(D, g)$, where $D = \text{span}\{F_1, F_2\}$ is a two-dimensional distribution and $g$ is a SR-metric. If $D$ is a contact distribution, near a point $q_0 \in \mathbb{R}^3$ identified with 0, one has the following.

* **Heisenberg-Brockett nilpotent model.** The nilpotent model of order $-1$ is the so-called Heisenberg-Brockett model introduced in [10] where the SR-problem is defined by the orthonormal frame $D = \text{span}\{\hat{F}, \hat{G}\}$

$$\hat{F} = \frac{\partial}{\partial x} + y \frac{\partial}{\partial z}, \quad \hat{G} = \frac{\partial}{\partial y} - x \frac{\partial}{\partial z}$$

in which $q = (x, y, z)$ are (local) graded privileged coordinates; the weights of $x, y$ are one and the weight of $z$ is two.

* **Generic model.** Using this gradation, the normal form of order 0 is similar and the generic model is given by the normal form of order 1

$$F = \hat{F} + yQ(w)\frac{\partial}{\partial z}, \quad G = \hat{G} - xQ(w)\frac{\partial}{\partial z},$$

$w = (x, y)$ and $Q$ is quadratic in $w$, $Q = \alpha x^2 + 2\beta xy + \gamma y^2$, where $\alpha, \beta, \gamma$ are parameters.

Introducing $\tilde{Q}(w) = (1 + Q(w))$, one can write

$$F = \frac{\partial}{\partial x} + y\tilde{Q}\frac{\partial}{\partial z}, \quad G = \frac{\partial}{\partial y} - x\tilde{Q}\frac{\partial}{\partial z}.$$

Computing, one has

$$[F, G] = (2\tilde{Q} + x\frac{\partial Q}{\partial x} + y\frac{\partial Q}{\partial y})\frac{\partial}{\partial z}$$

and using Euler formula

$$[F, G] = 2(1 + 2Q)\frac{\partial}{\partial z} = f(w)\frac{\partial}{\partial z}.$$
**Geodesics equations.** They are computed in this approximation, using Poincaré coordinates associated with the frame $(F, G, \frac{\partial}{\partial z})$:

$$H_1 = p \cdot F, \quad H_2 = p \cdot G, \quad H_3 = p \cdot \frac{\partial}{\partial z}.$$

The normal Hamiltonian is

$$H_n = (H_1^2 + H_2^2)/2.$$

One has

$$\dot{H}_1 = dH_1(\vec{H}_n) = \{H_1, H_2\}H_2 = (p \cdot [F, G])H_2 = p_z f(w) H_2$$

with $p_z = H_3$ is a constant (isoperimetric situation).

$$\dot{H}_2 = dH_2(\vec{H}_n) = -\{H_1, H_2\} H_1 = -p_z f(w) H_1.$$

Since $f(w) = 2 + O(|w|^2)$ and for $p_z$ non zero, we can reparameterize using

$$ds = p_z f(w) dt. \qquad (11)$$

Denoting $\phi'$ the derivative of the function $\phi$ with respect to $s$, we obtain

$$H_1' = H_2, \quad H_2' = -H_1$$

**Lemma 2.8.** In the $s$-parameter, the normal controls are solutions of the linear pendulum equation $H_1'' + H_1 = 0$ and are trigonometric functions.

**Geometric comments.** This calculation confirms up to reparameterization the choice of sinusoidal paddling in [22] and more generally validates Fourier expansions to compute strokes.

The remaining equation to be integrated are

$$x' = \frac{H_1}{p_z f(w)}, \quad y' = \frac{H_2}{p_z f(w)}, \quad z' = \frac{(H_1 y - H_2 x) \tilde{Q}}{p_z f(w)}. \qquad (12)$$

The solution can be estimated with micro-local expansions, see [2, 8] for the problem of computing conjugate points. Note that the weights are similar to analyze periodicity since conjugate times correspond to periods.

**Heisenberg-Brockett solution.** For the model of order $-1$, the geodesics equation can be computed, invoking Liouville integrability, and we recall the following [10].

**Lemma 2.9.** The solution starting from $q(0) = 0$ is parameterized by

$$x(t) = \frac{A}{\lambda} (\sin(\lambda t + \varphi) - \sin(\varphi)) \qquad (13)$$

$$y(t) = \frac{A}{\lambda} (\cos(\lambda t + \varphi) - \cos(\varphi)) \qquad (14)$$

$$z(t) = \frac{A^2}{\lambda} t - \frac{A^2}{\lambda^2} \sin(\lambda t) \qquad (15)$$

where $A, \lambda, \varphi$ are parameters defined by $p(0)$.
<figure><img src="image_3.png" /><figcaption>FIGURE 3. Two parameters families of circles (obtained by varying the amplitude and applying the symmetry of revolution) which are projections of geodesics of the Heisenberg-Brockett problem.</figcaption></figure>

**Relation with the copepod problem.** In relation with the swimmer problem, some information can be obtained taking $(x, y)$ as the shape variables. Indeed projecting the geodesics flows, we get a one parameter family of circles on each energy level associated with $A$, each of them deduced by a proper rotation $R_\alpha$ along the $z$-axis, associated with a symmetry of the SR-model, see Fig.3.

The associated displacement is given by

$$z = \oint (ydx - xdy) = \int d\omega$$

where $d\omega = 2(dy \wedge dx)$ is proportional to the standard $\mathbb{R}^2$-volume form.

**Lemma 2.10.** After a period $t = 2\pi/\lambda$, $z(2\pi/\lambda) = A^2 2\pi/\lambda^2$ but one has $\dot{z} > 0$ and hence $z$ is always increasing.

Since for Copepod strokes there always exists a backward motion, to produce a final forward displacement (at the end of each time period), clearly the $z$ variable cannot be identified with the displacement. It is an observed fact (in this zooplankton behaviour) which is mathematically justified by using a small square stroke and Baker-Campbell-Hausdorff formula [19].

The question of identifying the displacement variable for strokes of small amplitude has an answer introducing the privileged coordinates. The analysis goes as follows.

We choose a point $\theta(0) = (\theta_1(0), \theta_2(0))$ in the interior of the triangle which can be set to 0. Invoking the translation

$$x = \theta_1 - \theta_1(0), \quad y = \theta_2 - \theta_2(0)$$

and using Taylor expansion one has the following expansion

$$\dot{z} = \frac{u_1 \sin \theta_1 + u_2 \sin \theta_2}{2 + \sin^2 \theta_1 + \sin^2 \theta_2} = u_1 c_1 + u_2 c_2 + u_1 \alpha_1(x, y) + u_2 \alpha_2(x, y)$$

where $c_1, c_2$ are suitable constants and $\alpha_i = o(1)$, $i = 1, 2$. Using $\dot{x} = u_1, \dot{y} = u_2$ we set

$$Z = z - c_1 x - c_2 y \qquad (16)$$

to get an equation of the form

$$\dot{Z} = u_1 \alpha_1(x, y) + u_2 \alpha_2(x, y).$$

This normalization has been displayed for the Euclidean model. Similar computations can be employed for the mechanical energy.

Hence, we obtain

**Lemma 2.11.** $q = (x, y, Z)$ provide, near $(\theta_1(0), \theta_2(0), 0)$, a graded system of coordinates with respective weights $(1, 1, 2)$, establishing a link between the physical coordinates and the normal coordinates.
Still, the nilpotent model is not enough to make a clear picture, due to the symmetry of revolution described in Fig.3: indeed every projection on the $(x, y)$-plane is periodic, and this property is not stable by perturbation. Hence to go further in the analysis one has to make use of higher order approximations. A natural approach consists in constructing the normal form of order zero, bearing in mind that, in this construction, shape and displacement variables are mixed up, the gauge of the rotation $R_\alpha$ is also fixed.

## 2.4. Computations of the normal form of order 1

### 2.4.1. Contact case.

We consider the following vector fields

$$ 
\begin{aligned}
F_1(x, y, z) &= \frac{\partial}{\partial x} + (a_{30}^0 x^3 + a_{21}^0 x^2 y + a_{12}^0 xy^2 + a_{03}^0 y^3 + a_{20}^0 x^2 + a_{11}^0 xy + a_{02}^0 y^2 + a_{10}^0 x \\
&\qquad + a_{01}^0 y + a_{00}^0) \frac{\partial}{\partial z} \\
F_2(x, y, z) &= \frac{\partial}{\partial y} + (-a_{30}^0 y^3 - a_{21}^0 xy^2 - a_{12}^0 x^2 y - a_{03}^0 x^3 + a_{20}^0 y^2 + a_{11}^0 xy + a_{02}^0 x^2 - a_{10}^0 y \\
&\qquad - a_{01}^0 x + a_{00}^0) \frac{\partial}{\partial z}
\end{aligned}
\quad (17) $$

where $a_{ij}^0 \in \mathbb{R}, \forall i, j = 0, 1, 2, 3$.

We apply at each step, an action on the vector fields as described in the algorithm of [2]. At each step, we write $(\tilde{x}, \tilde{y}, \tilde{z})$ and $(x, y, z)$ respectively the 'new' and the 'old' variables, employing a change of coordinates $(\tilde{x}, \tilde{y}, \tilde{z}) = \varphi(x, y, z)$. The 'new' vector fields are then given by the standard push-forward map: $\tilde{F} = \varphi_* F$ and at the end of each step, we abuse notation and we rewrite $\tilde{F}(\tilde{x}, \tilde{y}, \tilde{z})$ as $F(x, y, z)$.

Step 1. We remove the terms $x^i \partial_i / \partial z$, $i \ge 0$ in $F_1$ and the terms $y^i \partial_i / \partial z$, $i \ge 0$ in $F_2$ by a change of coordinates $\varphi_1$ such that

$$ (\tilde{x}, \tilde{y}, \tilde{z}) = \varphi_1(x, y, z) = (x, y, z - \frac{1}{4} a_{30}^0 x^4 - \frac{1}{3} a_{20}^0 x^3 - \frac{1}{2} a_{10}^0 x^2 - a_{00}^0 x \\ + \frac{1}{4} a_{03}^0 y^4 + \frac{1}{3} a_{02}^0 y^3 + \frac{1}{2} a_{01}^0 y^2 - a_{00}^0 y) $$

and the resulting vector fields are $\tilde{F}_i = \varphi_{1*} F_i$, $i = 1, 2$ which are denoted again by $F_1, F_2$ as

$$ 
\begin{aligned}
F_1 &= \frac{\partial}{\partial x} + (a_{21}^1 x^2 y + a_{12}^1 xy^2 + a_{03}^1 y^3 + a_{11}^1 xy + a_{02}^1 y^2 + y/2) \frac{\partial}{\partial z} \\
F_2 &= \frac{\partial}{\partial y} + (-a_{21}^1 xy^2 - a_{12}^1 x^2 y - a_{03}^1 x^3 + a_{11}^1 xy + a_{02}^1 x^2 - x/2) \frac{\partial}{\partial z}
\end{aligned}
$$

Step 2. To remove the terms of order 0 in $F_1, F_2$, a change of coordinates $\varphi_2$ is performed where

$$ 
\begin{aligned}
\varphi_2(x, y, z) = (x + c_{300} x^3 + c_{210} x^2 y + c_{120} xy^2 + c_{030} y^3 + c_{110} (xy + y^2) + c_{011} yz + c_{001} z, \\
y + c_{030} x^3 + c_{120} x^2 y + c_{210} xy^2 + c_{300} y^3 - c_{110} (xy + x^2) - c_{011} xz + c_{001} z, \\
z + d_{300} x^3 + d_{210} x^2 y + d_{120} xy^2 + d_{030} y^3 + d_{101} xz + d_{011} yz),
\end{aligned}
\quad (18) $$

followed by the feedback

$$ u_1 \leftarrow u_1 + u_2 (2c_{001} x + 2c_{001} y), \quad u_2 \leftarrow u_2 - u_1 (2c_{001} x + 2c_{001} y). $$
The resulting vector fields are written as

$$ 
\begin{aligned}
F_1 &= \frac{\partial}{\partial x} + (a_{30}^2 x^3 + a_{21}^2 x^2 y + a_{12}^2 xy^2 + a_{03}^2 y^3 + y/2) \frac{\partial}{\partial z} \\
F_2 &= \frac{\partial}{\partial y} + (-a_{30}^2 y^3 - a_{12}^2 x^2 y - a_{21}^2 xy^2 - a_{03}^2 x^3 - x/2) \frac{\partial}{\partial z}
\end{aligned}
$$

Step 3. We remove the terms $a_{30}^3 x^3 \partial/\partial z$ and $a_{12}^3 xy^2 \partial/\partial z$ in $F_1$ and the terms $-a_{30}^3 y^3 \partial/\partial z$ and $-a_{12}^3 x^2 y \partial/\partial z$ in $F_2$ by the change of coordinates

$$ \varphi_3(x, y, z) = (1/\sqrt{2} x - 1/\sqrt{2} y, 1/\sqrt{2} x + 1/\sqrt{2} y, z - a_{30}^3/4 x^4 + a_{30}^3/4 y^4) $$

followed by the feedback

$$ u_1 \leftarrow 1/\sqrt{2} u_1 + 1/\sqrt{2} u_2, \quad u_2 \leftarrow -1/\sqrt{2} u_1 + 1/\sqrt{2} u_2. $$

Therefore, we obtain

$$ 
\begin{aligned}
F_1 &= \frac{\partial}{\partial x} + \left( \left(-\frac{a_{21}^2}{2} + \frac{a_{12}^2}{2} + \frac{3}{2} a_{03}^2\right) x^2 y + \left(\frac{a_{21}^2}{2} - \frac{a_{12}^2}{2} + \frac{a_{03}^2}{2}\right) y^3 + \frac{y}{2} \right) \frac{\partial}{\partial z} \\
F_2 &= \frac{\partial}{\partial y} + \left( \left(\frac{a_{21}^2}{2} + \frac{a_{12}^2}{2} - \frac{3}{2} a_{03}^2\right) xy^2 + \left(-\frac{a_{21}^2}{2} - \frac{a_{12}^2}{2} - \frac{a_{03}^2}{2}\right) x^3 - \frac{x}{2} \right) \frac{\partial}{\partial z}
\end{aligned}
$$

Step 4. We normalize the terms of order 1 with a quadratic form $Q$ such that $\partial^2 Q/(\partial x \partial y) = 0$ by the diffeomorphism

$$ \varphi_4(x, y, z) = (x, y, z - 1/4(a_{03}^2 - a_{21}^2) x^3 y - 1/4(a_{21}^2 - a_{03}^2) xy^3) $$

$$ F_1 = \frac{\partial}{\partial x} + \frac{y}{2} (1 + Q(x, y)) \frac{\partial}{\partial z}, \quad F_2 = \frac{\partial}{\partial y} - \frac{x}{2} (1 + Q(x, y)) \frac{\partial}{\partial z} \quad (19) $$

where

$$ Q(x, y) = (1/2 a_{21}^2 + a_{12}^2 + 3/2 a_{03}^2) x^2 + (1/2 a_{21}^2 - a_{12}^2 + 3/2 a_{03}^2) y^2 $$

**Application.** We apply this algorithm for the copepod model where we consider the initial point on the axis of symmetry $\Sigma$ (see Fig.2) $\theta_1(0) = \theta_{10}$, $\theta_2(0) = \pi - \theta_{10}$ and computations are done using the Euclidean model.

One sets: $x = \theta_1 - \theta_{10}$, $y = \theta_2 + \theta_{10} - \pi$ to compute the expression near $(x, y) \sim (0, 0)$ and we denote $z$ the displacement variable. Computing the Taylor expansion of $F_1, F_2$ near 0 up to order 4, the resulting vector fields can be put in the form of (17).

In this step 2, the important point is to keep track of the mixing of the shape variables with the displacement variable. Up to order 2, we have the following.

**Lemma 2.12.** After step 2, up to order 2 with respect to the gradation (1, 1, 2) the relation between the new coordinates and Heisenberg coordinates is given by

$$ \tilde{x} = \varphi_2^1(x, y, z), \quad \tilde{y} = \varphi_2^2(x, y, z) \quad (20) $$

where $\varphi_2^j$ is the j-th component of $\varphi_2$. In particular, this transformation mixes up the shape and the displacement variable.

Lemma 2.12 is the key tool to solve the periodicity equations on the model of order 0 :

$$ \theta(0) = \theta(2\pi), \quad \dot{\theta}(0) = \dot{\theta}(2\pi). $$
**Definition 2.13.** Let $\mathcal{G}'$ the subgroup of $\mathcal{G}$ of local diffeomorphisms $\varphi$ of $\mathbb{R}^3$ defined by $\varphi(\theta, z) = (\varphi_1(\theta), \varphi_2(\theta, z))$ where $\varphi_1 : \mathbb{R}^2 \to \mathbb{R}^2$ and $\varphi_2 : \mathbb{R}^3 \to \mathbb{R}$.

Restricting $\varphi_2$ to $\mathcal{G}'$ imposes $c_{001} = c_{011} = 0$ where $c_{001}, c_{011}$ appear in (18). Symbolic computations yield the following expressions

$$c_{001}(\theta_{10}) = 1/2 \frac{\cos^4(\theta_{10}) + 3 \cos^2(\theta_{10}) - 2}{(\cos^2(\theta_{10}) - 2) \sin(\theta_{10}) \cos(\theta_{10})} \\ c_{011}(\theta_{10}) = 1/2 \frac{\cos^8(\theta_{10}) + 6 \cos^6(\theta_{10}) + 5 \cos^4(\theta_{10}) - 12 \cos^2(\theta_{10}) + 4}{(\cos^6(\theta_{10}) - 5 \cos^4(\theta_{10}) + 8 \cos^2(\theta_{10}) - 4) \cos^2(\theta_{10})} = -2 (c_{001}(\theta_{10}))^2$$

We check numerically that these two coefficients $c_{001}$ and $c_{011}$ are both zero for one value of $\theta_{10} \in [-\pi, \pi]$ which is $\theta_{10}^* \simeq 0.7236888$. Moreover $(\theta_{10}^*, \pi - \theta_{10}^*)$ corresponds to the center of the family of simple strokes presented in Fig.2.

Therefore, we have

**Proposition 2.14.** $P = (\theta_{10}^*, \pi - \theta_{10}^*)$ is the only point of the line $\Sigma : \theta_2 = \pi - \theta_1$ such that the transformation $\varphi_2$ given in (18) is an element of $\mathcal{G}'$.

Finally, for the copepod model, the quadratic form in the normal form (19) is

$$Q(x, y) = 1/8 \frac{(-11 \cos^6(\theta_{10}) + 17 \cos^4(\theta_{10}) - 8 \cos^2(\theta_{10}) + 4) x^2}{(\cos^6(\theta_{10}) - 5 \cos^4(\theta_{10}) + 8 \cos^2(\theta_{10}) - 4) \cos^2(\theta_{10})} \\ + 1/8 \frac{(5 \cos^8(\theta_{10}) - 33 \cos^6(\theta_{10}) + 2 \cos^4(\theta_{10}) + 28 \cos^2(\theta_{10}) - 8) y^2}{(\cos^6(\theta_{10}) - 5 \cos^4(\theta_{10}) + 8 \cos^2(\theta_{10}) - 4) \cos^2(\theta_{10})}$$

and at the center point $(\theta_{10}^*, \pi - \theta_{10}^*)$, we have

$$Q(x, y) = -0.7165898586 x^2 - 0.7379854942 y^2.$$

**Remark 2.15.** Truncating the system to order 0, we obtain geodesics corresponding to periodic controls whose projection on a two-dimensional surface diffeomorphic to a plane is represented by periodic simple loops (see Fig.2). This surface correspond to the $\theta$-plane only in the case where $(\theta_1(0), \theta_2(0))$ is taken as the center of swimming.

**Remark 2.16.** The shape variable and the displacement variable are mixed only in step 2. Therefore, when this algorithm is applied at the center $(\theta_{10}^*, \pi - \theta_{10}^*)$, we use only diffeomorphisms of $\mathcal{G}'$.

**Remark 2.17.** The family of simple strokes and the center of swimming is presented in Fig.2 both for the Euclidean cost and the mechanical cost.

### 2.4.2. The Martinet case

We recall that each point of the boundary of the triangle $\mathcal{T}$, vertices excluded, is a Martinet point. The analysis at a Martinet point is more intricate and as before relies on the construction of a generic model, which is in this case of order 0.

**Martinet model of order 0.** The point $q_0$ is identified to zero and there exist local coordinates $q = (x, y, z)$ such that the SR-geometry is given by $(D, g)$ where:

* $D = \text{span}\{F, G\}$, $F = \frac{\partial}{\partial x} + \frac{y^2}{2} \frac{\partial}{\partial z}$ and $G = \frac{\partial}{\partial z}$ with $q = (x, y, z)$ are graded coordinates with respective weights $(1, 1, 3)$.

* The metric $g$ is of the form $a(q)dx^2 + b(q)dy^2$ where we have
  - Model of order $-1$. (Flat Martinet case): $a = b = 1$
– Generic model:

$$a = (1 + \alpha y)^2 \sim 1 + 2\alpha y \text{ (order zero)}$$

$$b = (1 + \beta x + \gamma y)^2 \sim 1 + 2\beta x + 2\gamma y \text{ (order zero)}$$

where $\alpha, \beta, \gamma$ are parameters and the square in the formula is related to the crucial computation presented next.

**Geodesic equations.** We introduce the orthonormal frame

$$F_1 = \frac{F}{\sqrt{a}}, \quad F_2 = \frac{G}{\sqrt{c}}, \quad F_3 = \frac{\partial}{\partial z}$$

and denoting $H_i = p \cdot F_i$, the normal Hamiltonian is given by $H_n := (H_1^2 + H_2^2)/2$.

We parameterize by arc-length: $H_1^2 + H_2^2 = 1$, $H_1 = \cos \chi$, $H_2 = \sin \chi$ and $H_3 = p_z = \lambda$ constant, assuming $\lambda \neq 0$. Hence the geodesic equations become

$$\begin{aligned}
\dot{x} &= \frac{\cos \chi}{\sqrt{a}}, & \dot{y} &= \frac{\sin \chi}{\sqrt{b}}, & \dot{z} &= \frac{y^2 \cos \chi}{2\sqrt{a}}, \\
\dot{\chi} &= \frac{1}{\sqrt{ac}}(y\lambda - \alpha \cos \chi - \beta \sin \chi).
\end{aligned} \quad (21)$$

We introduce the new parameterization

$$\sqrt{ac} \frac{d}{dt} = \frac{d}{ds} \quad (22)$$

and denoting by $\phi'$ the derivative of a function $\phi$ with respect to $s$, we get the equations

$$y' = \sin \chi(1 + \alpha y), \quad \chi' = (y\lambda - \alpha \cos \chi - \beta \sin \chi)$$

which are equivalent to the second order differential equation

$$\chi'' + \lambda \sin \chi + \alpha^2 \sin \chi \cos \chi - \alpha \beta \sin^2 \chi + \beta \chi' \cos \chi = 0. \quad (23)$$

As a consequence, we obtain the following result.

**Proposition 2.18.** *The generic case projects, up to a time reparameterization, onto a two dimensional equation (23), associated with a generalized dissipative pendulum depending on the parameters $\alpha, \beta$ only.*

**Geometric application.** This pendulum equation provides the dynamics of the controls $u_1 = H_1 = \sin \chi$, $u_2 = H_2 = \cos \chi$ and describes the main properties of the generic model of order 0.

* Flat case: $\alpha = \beta = \gamma = 0$, it corresponds to the standard pendulum which is integrable using elliptic functions.

* Generic case:

- $\beta = 0$: it corresponds to the integrable case with the existence of an additional first integral, linear in the momentum and identified with $p_x$. Moreover using (21) again, the equations can be integrated using elliptic functions.

- $\beta \neq 0$: due to dissipation, we are in the non-integrable case.

An additional property being:

**Lemma 2.19.** *The abnormal line identified with $t \to (t, 0, 0)$ is strictly abnormal if and only if $\alpha \neq 0$.*
**Application to the copepod.** For the copepod case, previous computations tell us that periodic strokes can appear for different values of the modulus *k* of the elliptic functions and they are:

* *k* = 0: circles,

* *k* ≃ 0.65: eight shapes called Bernoulli lemniscates.

Note that in this analysis $\beta = 0$ is not a stable model, and the state constraints represented by the triangle $\mathcal{T}$ are not taken into account.

### 2.4.3. Numerical simulations.

Numerical simulations using the HamPath software are presented on Fig.4 based on the geometric analysis and to achieve the computations of strokes corresponding respectively to simple loops, limaçons and eight shape curves. Note that they are performed using the Euclidean cost.

* *Simple loops.* They are obtained by continuation of small circles.

* *Limaçons.* They are generated by perturbations of simple loops, by period doubling.

* *Eight-shapes.* They appear only on the sides of the triangle.

<figure><img src="image_4.png" /><figcaption>FIGURE 4. One parameter family of simple loops, limaçons and Bernoulli lemniscates normal strokes for the Euclidean metric.</figcaption></figure>

## 3. COPEPOD PROPERTIES COUPLING OPTIMAL CONTROL THEORY AND NUMERICAL SIMULATIONS

### 3.1. Maximum principle with periodic controls

First of all, one needs a neat statement of the *transversality conditions* suitable to analyze the maximizing different concepts of efficiencies with periodic control and we use general tools in optimal control like those ones appearing in [23]. These techniques allow us to establish important properties on the Copepod model where in the problem we insert an efficiency term to compare strokes (for a discussion on the significance of the efficiency term as a physical criterion to select 'optimal' strokes, we refer the reader to [22]). The problem is written in the extended state space $\bar{q} = (q^0, q)$:

$$ 
\left\{ 
\begin{array}{l} 
\min_{u(\cdot) \in \mathcal{U}} h(\bar{q}(0), \bar{q}(2\pi)) \\ 
\dot{q} = \sum_{i=1}^{2} u_i F_i(q) = F(q, u), \quad u \in \mathbb{R}^2 \\ 
\dot{q}^0 = \sum_{i=1}^{2} u_i^2, \quad q^0(0) = 0. 
\end{array} 
\right. 
\qquad (24) 
$$

with end-points conditions of the form $(\bar{q}(0), \bar{q}(2\pi)) \in C$, where $C \subset \mathbb{R}^3 \times \mathbb{R}^3$ is a given closed set.

We denote $\bar{p} = (p_0, p)$ the extended adjoint vector. The pseudo-Hamiltonian takes the form

$$ H(\bar{q}, \bar{p}, u) = \sum_{i=1}^{2} u_i H_i + p_0(u_1^2 + u_2^2) $$
and $H_i = p \cdot F_i(q)$.

From the maximum principle (cf. [23]), an optimal control pair $(q, u)$ satisfies the following necessary conditions that we split into two distinct parts

**Standard conditions.**

$$ \dot{\bar{q}} = \frac{\partial H}{\partial \bar{p}}, \quad \dot{\bar{p}} = -\frac{\partial H}{\partial \bar{q}}, \quad \frac{\partial H}{\partial u} = 0. $$

**Transversality conditions.**

$$ (\bar{p}(0), -\bar{p}(2\pi)) \in \lambda \nabla h(\bar{q}(0), \bar{q}(2\pi)) + N_C(\bar{q}(0), \bar{q}(2\pi)) \quad (25) $$

where $N_C$ is the (limiting) normal cone to the (closed) set $C$, $(\bar{p}, \lambda) \neq 0$, $\lambda \ge 0$.

**Application.** $q = (x_0, \theta_1, \theta_2)$

* Maximizing the geometric efficiency with periodic condition

$$ \theta(0) = \theta(2\pi), \quad h = -\frac{x_0^2(2\pi)}{E} $$

where $E$ is the energy

$$ E = \int_{0}^{2\pi} (u_{1}^{2} + u_{2}^{2})dt. $$

In this case, from (25) we deduce

**Lemma 3.1.** Denoting $p_\theta = (p_{\theta_1}, p_{\theta_2})$ we have: $p_\theta(0) = p_\theta(2\pi)$, to produce a smooth stroke in the normal case $p_0 \neq 0$. Moreover at the final point $(p_0, p_{x_0})$ has to be collinear to the gradient of $h$ at the point $(\bar{q}^0(2\pi), \bar{x}_0(2\pi))$ on the level set $h(q^0, x_0) = c$, where $c$ is the minimum value of the cost.

* Maximizing an efficiency depending on $\theta(0)$, with periodic condition:

$$ \theta(0) = \theta(2\pi), \quad h = -\frac{m(\theta(0)) x_0^2(2\pi)}{E}, \quad (\text{with } E = q^0(2\pi)) $$

where $m$ is a chosen smooth function.

In this case (25) leads to

$$ p_{\theta}(0) - p_{\theta}(2\pi) = \lambda \frac{\partial h}{\partial \theta(0)} \quad (26) $$

hence producing a jump of the adjoint vector at $t = 2\pi$.

## 3.2. Symmetric and non symmetric case

* Using Lemma 2.1, we have: in the $\theta$-plane, the one-parameter family of simple loops, represented in Fig.2, is symmetric with respect to the straight line $\Sigma: \theta_2 = \pi - \theta_1$ and corresponds to a velocity vector $\dot{\theta}(0)$ transverse to $\Sigma$.

* A non symmetric case is presented in Section 3.4 (see Fig.8) obtained by numerical continuation from the symmetric case.

## 3.3. Second-order necessary conditions

Standard second-order necessary optimality conditions are related to the non existence of conjugate points, see [7] and [24].

**Proposition 3.2.** Let $(x_0(t), \theta(t))$, $t \in [0, 2\pi]$ be a strictly normal stroke. Then a necessary optimality condition is the non existence of conjugate time $t_c \in ]0, 2\pi[$.

The existence of conjugate times can be checked numerically using the HamPath code.
### 3.4. Applications: numerical simulations

We present below a sequence of simulations using the HamPath software, considering the Euclidean cost $\int_{0}^{2\pi}(u_{1}^{2}+u_{2}^{2})dt$, and the mechanical energy cost (3). Also we consider a non symmetric case where we minimize the cost $\int_{0}^{2\pi}(u_{1}^{2}-\cos(\theta_{1})/2u_{1}u_{2}+1/3u_{2}^{2})dt$ and the simple loop strokes are not anymore symmetric with respect to the line $\Sigma$ (see Fig.8).

**Remark 3.3.** The existence of minimizers is guaranteed to the fact that efficiency of small strokes tends to zero when the amplitude goes to zero, thanks to Green formulae.

In particular we shall consider:

* Different kinds of normal strokes for the mechanical cost: simple loop, limaçon, eight and computations of conjugate points (see Fig.5). Only simple loops are candidates to be optimal strokes.

* A one-parameter family of simple strokes, each of them associated with a different energy. The corresponding efficiency is represented in Table 1 for the Euclidean case, in Table 2 for the mechanical energy and in Table 3 for the non symmetric case. They are compared with the efficiency of the abnormal stroke and we deduce the normal stroke with the best efficiency.

<figure><img src="image_5.png" /><figcaption>FIGURE 5. Normal strokes: simple loop (left), limaçon with inner loop (right) and eight case (bottom). First conjugate points on $[0, 2\pi]$ appear with a cross except for the simple loop stroke.</figcaption></figure>

Numerical simulations establish also that the abnormal triangle stroke, has a very low efficiency and, in particular, is not optimal. Indeed, we can produce the same displacement by doing twice a normal stroke, but requiring smaller energy.

In section 4, we shall provide an analytical interpretation of this numerical results.

**Remarks about algorithms.** Two algorithms can be used for the problem of efficiency in the smooth case, both implementable in the HamPath code.
* Method 1. The displacement is fixed to $x_0(2\pi) = x_T$ and the one parameter family is computed solving the shooting equation

$$ \theta(0) = \theta(2\pi), \quad p_{\theta}(0) = p_{\theta}(2\pi) $$

to generate the one parameter family $\mathcal{F}_{\lambda}$ of periodic strokes initiated from $\theta_0, p_0$ with $\lambda = 0$ and the one parameter family is generated by continuation. The efficiency maximization is obtained by using the Jacobi equation (also called variational equation) to compute the derivative of the mapping $\lambda \to h(\exp 2\pi \overleftarrow{H}_n(\theta(\lambda), p_{\theta}(\lambda)))$. Note that this algorithm is comparable with the method of computing conjugate points. Table (1), Table (2) and Table (3) provide numerical results which allow to detect optimal strokes for the Euclidean case, the mechanical energy and the non symmetric case respectively.

* Method 2. A second method is based on Maximum principle and the transversality condition (25), solving the associated shooting equation.

Both methods have been used in the smooth case, providing the same solutions: Fig.6 and Fig.7 represent the optimal stroke maximizing the geometric efficiency respectively for the Euclidean case, the mechanical case and the non symmetric case.

**Non smooth case.** As in [11], we consider a cost depending on $\theta(0)$, namely

$$ h = - \frac{x_0^2(2\pi)m(\theta(0))}{E} $$

where $m(\theta(0)) = 2 - \cos^2(\theta_1(0))$ and $E$ is the mechanical energy (3). In Fig.9 is illustrated the corresponding optimal solution satisfying the transversality conditions (25) which leads to (26). We call this case "non-smooth" since the presence of the term $m(\theta(0))$ in the cost produces a jump of the adjoint variable at the end of each stroke, cf (26).

<table><caption>TABLE 1. Geometric efficiency for the abnormal stroke and different normal strokes with the Euclidean cost.</caption><thead><tr><th>Types of strokes</th><th>x<sub>0</sub>(T)</th><th>l(γ)</th><th>x<sub>0</sub>(T)/l(γ)</th></tr></thead><tbody><tr><td rowspan="5">Simple loops</td><td>5.500 × 10<sup>-2</sup></td><td>1.984</td><td>2.520 × 10<sup>-2</sup></td></tr><tr><td>1.400 × 10<sup>-1</sup></td><td>3.785</td><td>3.698 × 10<sup>-2</sup></td></tr><tr><td>1.700 × 10<sup>-1</sup></td><td>4.340</td><td>3.917 × 10<sup>-2</sup></td></tr><tr><td>2.000 × 10<sup>-1</sup></td><td>4.946</td><td>4.043 × 10<sup>-2</sup></td></tr><tr><td>2.100 × 10<sup>-1</sup></td><td>5.109</td><td>4.110 × 10<sup>-2</sup></td></tr><tr><td rowspan="6"><strong>Optimal stroke</strong><br>Fig.6 (left)</td><td><strong>2.169 × 10<sup>-1</sup></strong></td><td><strong>5.180</strong></td><td><strong>4.187 × 10<sup>-2</sup></strong></td></tr><tr><td>2.200 × 10<sup>-1</sup></td><td>5.354</td><td>4.109 × 10<sup>-2</sup></td></tr><tr><td>2.300 × 10<sup>-1</sup></td><td>5.624</td><td>4.089 × 10<sup>-2</sup></td></tr><tr><td>2.500 × 10<sup>-1</sup></td><td>6.305</td><td>3.965 × 10<sup>-2</sup></td></tr><tr><td>2.740 × 10<sup>-1</sup></td><td>9.046</td><td>3.028 × 10<sup>-2</sup></td></tr><tr><td>Abnormal</td><td>2.742 × 10<sup>-1</sup></td><td>10.73</td><td>2.555 × 10<sup>-2</sup></td></tr><tr><td>Limaçon</td><td>2.000 × 10<sup>-1</sup></td><td>6.147</td><td>3.253 × 10<sup>-2</sup></td></tr></tbody></table>
<table><caption>TABLE 2. Geometric efficiency for the abnormal stroke and different normal strokes with the mechanical cost.</caption><thead><tr><th>Types of strokes</th><th>x<sub>0</sub>(T)</th><th>l(γ)</th><th>x<sub>0</sub>(T)/l(γ)</th></tr></thead><tbody><tr><td rowspan="6">Simple loops</td><td>0.500 × 10<sup>-1</sup></td><td>9.935 × 10<sup>-1</sup></td><td>5.03 × 10<sup>-2</sup></td></tr><tr><td>1.500 × 10<sup>-1</sup></td><td>1.860</td><td>8.06 × 10<sup>-2</sup></td></tr><tr><td>1.700 × 10<sup>-1</sup></td><td>2.021</td><td>8.41 × 10<sup>-2</sup></td></tr><tr><td>2.000 × 10<sup>-1</sup></td><td>2.280</td><td>8.77 × 10<sup>-2</sup></td></tr><tr><td>2.100 × 10<sup>-1</sup></td><td>2.501</td><td>8.84 × 10<sup>-2</sup></td></tr><tr><td>2.200 × 10<sup>-1</sup></td><td>2.472</td><td>8.89 × 10<sup>-2</sup></td></tr><tr><td rowspan="4"><strong>Optimal stroke</strong><br>Fig.6 (right)</td><td><strong>2.228 × 10<sup>-1</sup></strong></td><td><strong>2.561</strong></td><td><strong>8.902 × 10<sup>-2</sup></strong></td></tr><tr><td>2.300 × 10<sup>-1</sup></td><td>2.586</td><td>8.895 × 10<sup>-2</sup></td></tr><tr><td>2.500 × 10<sup>-1</sup></td><td>2.854</td><td>8.76 × 10<sup>-2</sup></td></tr><tr><td>2.600 × 10<sup>-1</sup></td><td>3.044</td><td>8.54 × 10<sup>-2</sup></td></tr><tr><td>Abnormal</td><td>2.742 × 10<sup>-1</sup></td><td>4.933</td><td>5.56 × 10<sup>-2</sup></td></tr><tr><td>Limaçon</td><td>2.500 × 10<sup>-1</sup></td><td>3.353</td><td>7.46 × 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>TABLE 3. Geometric efficiency for the abnormal stroke and different normal strokes for the non symmetric case.</caption><thead><tr><th>Types of strokes</th><th>x<sub>0</sub>(T)</th><th>l(γ)</th><th>x<sub>0</sub>(T)/l(γ)</th></tr></thead><tbody><tr><td rowspan="8">Simple loops</td><td>0.1</td><td>2.09</td><td>4.79 × 10<sup>-2</sup></td></tr><tr><td>0.15</td><td>2.77</td><td>5.42 × 10<sup>-2</sup></td></tr><tr><td>0.16</td><td>2.91</td><td>5.5 × 10<sup>-2</sup></td></tr><tr><td>0.18</td><td>3.21</td><td>5.61 × 10<sup>-2</sup></td></tr><tr><td>0.19</td><td>3.37</td><td>5.64 × 10<sup>-2</sup></td></tr><tr><td>0.196</td><td>3.47</td><td>5.65 × 10<sup>-2</sup></td></tr><tr><td>0.198</td><td>3.5</td><td>5.655 × 10<sup>-2</sup></td></tr><tr><td>0.2</td><td>3.54</td><td>5.657 × 10<sup>-2</sup></td></tr><tr><td rowspan="3"><strong>Optimal stroke</strong><br>Fig.7</td><td><strong>0.2033</strong></td><td><strong>3.59</strong></td><td><strong>5.658 × 10<sup>-2</sup></strong></td></tr><tr><td>0.21</td><td>3.71</td><td>5.653 × 10<sup>-2</sup></td></tr><tr><td>0.22</td><td>3.91</td><td>5.627 × 10<sup>-2</sup></td></tr><tr><td>Abnormal</td><td>2.742 × 10<sup>-1</sup></td><td>8.432</td><td>3.25 × 10<sup>-2</sup></td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>FIGURE 6. Optimal stroke of the Copepod swimmer for the Euclidean cost (left) and the mechanical energy (right), obtained by the transversality conditions of the maximum principle (25).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 7. Optimal stroke for the non symmetric case obtained by the transversality condition (25).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIGURE 8. One parameter family of simple loops for the non symmetric case. We perform a continuation considering the Hamiltonian $H_{\lambda} = \lambda H_{I} + (1 - \lambda) H_{F}$, $\lambda \in [0, 1]$ where $H_{I}, H_{F}$ are respectively the Hamiltonians associated with the Euclidean case and the non symmetric case. The path of centers of these Hamiltonians $H_{\lambda}, \lambda \in [0, 1]$ is represented.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 9. Non-smooth normal stroke of the Copepod swimmer for the cost depending on the initial angle $\theta(0)$ and with the Euclidean cost, obtained by the transversality conditions of the maximum principle (25)</figcaption></figure>
## 4. GEOMETRIC EFFICIENCY AND OPTIMALITY OF THE ABNORMAL STROKE

In the previous section, we have displayed numerical simulations indicating that the abnormal stroke is not optimal when the performance criterion involves an efficiency term. This section is devoted to establish an analytical interpretation of this property.

More precisely we shall discuss how the geometric efficiency can be related to a reference optimal control problem for the Copepod model in which we minimize integral cost criteria and how an efficiency term can provide information on minimizers for this model. We shall consider different right end-point conditions: the ‘exact displacement’ case ($x_0(T) = x_T$), the fixed energy case and the case in which we are given a target for the Copepod displacement (for instance $x_0(T) \in [x_T, 2x_T]$). For the latter we shall derive also normality conditions for minimizers.

Standard existence results (cf. [23, Section 2.8]) make sure that in all the problems considered in this section, minimizers exist.

### 4.1. Preliminary comments on the efficiency term

We start considering here the following optimal control problem for the Copepod model:

$$ 
\left\{ 
\begin{array}{l} 
\displaystyle \text{Minimize} \quad \int_0^T \ell(q(t), u(t)) dt \\ 
\text{over arcs } q(\cdot) \in W^{1,1}([0,T]; \mathbb{R}^3) \text{ s.t.} \\ 
\dot{q}(t) = F(q(t), u(t)), \quad \text{for a.e. } t \in [0,T] \\ 
u(t) \in U \quad \text{for a.e. } t \in [0,T] \\ 
q_1(T) = x_T, \quad q_1(0) = 0 \\ 
q_2(0) = q_2(T), \quad q_3(0) = q_3(T), 
\end{array} 
\right. 
\qquad (27) 
$$

in which $x_T > 0$ is a given number which represents the Copepod desired displacement, $T = 2\pi$ in the time period, $U = \mathbb{R}^2$, $\ell: \mathbb{R}^3 \times \mathbb{R}^2 \to \mathbb{R}$, and $F: \mathbb{R}^3 \times \mathbb{R}^2 \to \mathbb{R}^3$ is the function defined in (6).

We provide first a simple relationship between minimizers of the optimal control problem (27) and minimizers of an optimal control problem in which we maximize the geometric efficiency $\mathcal{E} = \frac{(q_1(T))^2}{q^0(T)}$ where $q^0(T) = \int_0^T \ell(q(t), u(t))dt$.

**Remark 4.1.** We observe first the fact that, for any integral cost term $\ell = \ell(q, u)$ (Lipschitz w.r.t. $q$ and Borel measurable w.r.t. $u$), if $(q^*, u^*)$ is a minimizer for (27), then $(q^*, u^*)$ is a minimizer also for the optimal control problem in which we replace the integral term in (27) with $J_{\mathcal{E}}(q(\cdot), u(\cdot)) := -\frac{(q_1(T))^2}{q^0(T)}$ (i.e. a problem in which we want to maximize the geometric efficiency $\mathcal{E}$).

This is immediate to prove. Indeed if $(q^*, u^*)$ is a minimizer for (27), then $q^0^*(T) \le q^0(T)$ for all admissible trajectories for problem (27). Taking the reciprocal, since the right end-point has to be the same $(q_1(T) = q_1^*(T) = x_T)$, we immediately obtain:

$$ -\frac{(q_1^*(T))^2}{q^0^*(T)} \le -\frac{(q_1(T))^2}{q^0(T)}, $$

for all admissible trajectories $q(\cdot)$ for the control system (27).

**Remark 4.2.** Observe that if we take $J_{\varepsilon} = -(q_1(T))^2/q^0(T)$ as a cost to minimize, then we can even leave free the right end-point of the displacement variable $q_1$ (not imposing a fixed or a minimal displacement). Indeed, from Hölder inequality, the term $(q_1(T))^2/q^0(T)$ is bounded (for all admissible trajectories for (28)) and well defined.

As a consequence the fact itself of considering the cost $J_{\varepsilon}$, produces minimizers with non-zero displacement.
We study next the case in which our aim is to maximize the displacement, once we are given the energy of the system in the euclidean case: $\ell(q, u) = u_1^2 + u_2^2$. This is equivalent to minimize the cost $J_{\varepsilon}$ involving the geometric efficiency and the final energy $E^* > 0$ becomes a right-end point condition for the $q^0$-variable. (Due to the symmetry of the Copepod model, we can always suppose that $q_1(T) > 0$.) We consider then:

$$ 
\left\{
\begin{array}{l}
\text{Minimize } J_{\mathcal{E}}((q^0(\cdot), q(\cdot)), u(\cdot)) := -\frac{(q_1(T))^2}{q^0(T)} \\
\text{over arcs } (q^0(\cdot), q(\cdot)) \in W^{1,1}([0, T]; \mathbb{R}^4) \text{ s.t.} \\
\dot{q}(t) = F(q(t), u(t)), \quad \text{for a.e. } t \in [0, T] \\
\dot{q}^0(t) = \ell(q(t), u(t)) = u_1^2(t) + u_2^2(t), \quad \text{for a.e. } t \in [0, T] \\
u(t) \in U \quad \text{for a.e. } t \in [0, T] \\
q_1(0) = 0, \quad q^0(0) = 0, \quad q^0(T) = E^* \\
q_2(0) = q_2(T), \quad q_3(0) = q_3(T),
\end{array}
\right.
\qquad (28) $$

where $E^* > 0$ is fixed. Suppose that we have a minimizer $((q^0^*, q^*), u^*)$ for (28). Necessarily we have $q^0^*(T) \neq 0$ and $q_1^*(T) \neq 0$. We have $E^* = q^0^*(T) = \int_0^T \ell(q^*(t), u^*(t))dt$ the corresponding energy and consider the following optimal control problem (for some constant $K > 0$):

$$ 
\left\{
\begin{array}{l}
\text{Minimize } J_{\mathcal{E},K}((q^0(\cdot), q(\cdot)), u(\cdot)) := \left[-\frac{(q_1(T))^2}{q^0(T)} + K|q^0(T) - E^*|\right] \\
\text{over arcs } (q^0(\cdot), q(\cdot)) \in W^{1,1}([0, T]; \mathbb{R}^4) \text{ s.t.} \\
\dot{q}(t) = F(q(t), u(t)), \quad \text{for a.e. } t \in [0, T] \\
\dot{q}^0(t) = \ell(q(t), u(t)) = u_1^2(t) + u_2^2(t), \quad \text{for a.e. } t \in [0, T] \\
u(t) \in U \quad \text{for a.e. } t \in [0, T] \\
q^0(0) = 0, \quad q_1(0) = 0 \\
q_2(0) = q_2(T), \quad q_3(0) = q_3(T),
\end{array}
\right.
\qquad (29) $$

Observe that we have now an optimal control problem with the free right-end displacement and energy. The geometric efficiency in the cost is now supplemented with a penalization term which forces the control system to stay close to the reference energy level $E^*$.

**Proposition 4.3.** Let $((q^0^*(\cdot), q^*(\cdot)), u^*(\cdot))$ be a weak local minimizer for (28), then we can find a positive constant $K$ such that $((q^0^*(\cdot), q^*(\cdot)), u^*(\cdot))$ is a weak local minimizer for (29).

*Proof.* Assume that $((q^0^*(\cdot), q^*(\cdot)), u^*(\cdot))$ is local a minimizer for (28). Then, there exists $\varepsilon_0 \in (0, x_T/2)$ such that

$$ J_{\mathcal{E}}((q^0^*, q^*)(\cdot), u^*(\cdot)) \leq J_{\mathcal{E}}((q^0, q)(\cdot), u(\cdot)), \qquad (30) $$

for all trajectory-control pair $(q(\cdot), u(\cdot))$ satisfying all the conditions of the control system in (27), and with

$$ \|q^*(\cdot) - q(\cdot)\|_{L^{\infty}} \leq \varepsilon_0, \quad \|u^*(\cdot) - u(\cdot)\|_{L^{\infty}} \leq \varepsilon_0. \qquad (31) $$

We claim that we can find $\varepsilon_1 \in (0, \varepsilon_0)$ small enough such that

$$ J_{\mathcal{E},K}((q^0^*(\cdot), q^*(\cdot)), u^*(\cdot)) \leq J_{\mathcal{E},K}((q^0(\cdot), q(\cdot)), u(\cdot)), \qquad (32) $$

for all trajectory-control pair $((q^0(\cdot), q(\cdot)), u(\cdot))$ satisfying all the conditions of the control system in (29), and with

$$ \|(q^0^*(\cdot), q^*(\cdot)) - (q^0(\cdot), q(\cdot))\|_{L^\infty} \leq \varepsilon_1, \quad \|u^*(\cdot) - u(\cdot)\|_{L^\infty} \leq \varepsilon_1. \qquad (33) $$

Suppose by contradiction that for each $\varepsilon \in (0, \varepsilon_0)$ we can find a trajectory/control pair $((q^0^\varepsilon(\cdot), q^\varepsilon(\cdot)), u^\varepsilon(\cdot))$ satisfying all the requirements of the control system in (29), such that

$$ J_{\mathcal{E},K}((q^0^\varepsilon(\cdot), q^\varepsilon(\cdot)), u^\varepsilon(\cdot)) < J_{\mathcal{E},K}((q^0^*(\cdot), q^*(\cdot)), u^*(\cdot)) = -\frac{(q_1^*(T))^2}{q^0^*(T)} \qquad (34) $$
and

$$ \|(q^{0*}(.), q^{*}(.)) - (q^{0\epsilon}(.), q^{\epsilon}(.))\|_{L^{\infty}} \le \epsilon, \quad \|u^{*}(.) - u^{\epsilon}(.)\|_{L^{\infty}} \le \epsilon. \quad (35) $$

Then, we define $\rho := \sqrt{\frac{q^{0*}(T)}{q^{0\epsilon}(T)}}$ and take the trajectory/control pair $((\tilde{q}^0(.), \tilde{q}(.)), \tilde{u}(.))$ where

$$ \tilde{u}(.) := \rho u^{\epsilon}(.), \quad \tilde{q}^0(0) = 0, \quad \tilde{q}_1(0) = 0 $$

and

$$ \tilde{q}_2(.) = \rho q_2^{\epsilon}(.), \quad \tilde{q}_3(.) = \rho q_3^{\epsilon}(.). $$

It is straightforward to see that $((\tilde{q}^0(.), \tilde{q}(.)), \tilde{u}(.))$ is a solution of the control system in (29), and, since $\tilde{q}^0(T) = q^{0*}(T)$,

$$ J_{\mathcal{E},K}((\tilde{q}^0(.), \tilde{q}(.)) = -\frac{(\tilde{q}_1(T))^2}{\tilde{q}^0(T)} = -\frac{(\tilde{q}_1(T))^2}{q^{0*}(T)}. $$

Observe that taking $\epsilon_1 \in (0, \epsilon_0)$ small enough, then for all $\epsilon \in (0, \epsilon_1)$, we obtain that

$$ \|q^*(.) - \tilde{q}(.)\|_{L^\infty} \le \epsilon_0, \quad \|u^*(.) - \tilde{u}(.)\|_{L^\infty} \le \epsilon_0. \quad (36) $$

And choosing $K > 0$ big enough, we also have that

$$ \left| \frac{(\tilde{q}_1(T))^2}{\tilde{q}^0(T)} - \frac{(q_1^\epsilon(T))^2}{q^{0\epsilon}(T)} \right| \le K |q^{0\epsilon}(T) - \tilde{q}^0(T)| = K |q^{0\epsilon}(T) - E^*| . \quad (37) $$

From (34) and (37) we would eventually deduce that

$$ \begin{aligned} J_{\mathcal{E}}((\tilde{q}^0(.), \tilde{q}(.)), \tilde{u}) &= -\frac{(\tilde{q}_1(T))^2}{\tilde{q}^0(T)} \\ &\le -\frac{(q_1^\epsilon(T))^2}{q^{0\epsilon}(T)} + K |q^{0\epsilon}(T) - E^*| = J_{\mathcal{E},K}((q^{0\epsilon}(.), q^{\epsilon}(.)), u^{\epsilon}(.)) \\ &< J_{\mathcal{E}}((q^{0*}(.), q^{*}(.)), u^{*}(.)), \end{aligned} $$

which contradicts the (local) minimality properties of $(q^*(.), u^*(.))$ (cf. (30)-(31)). □

**Corollary 4.4.** Let $((q^{0*}(.), q^{*}(.)), u^{*}(.))$ be a weak local minimizer for (28), then the Maximum Principle applies in the normal form.

*Proof.* If $((q^{0*}(.), q^{*}(.)), u^{*}(.))$ is a local minimizer for (28), then from Proposition 4.3 we can find a positive constant $K$ such that $((q^{0*}(.), q^{*}(.)), u^{*}(.))$ is a local minimizer also for (29). Invoking a standard Maximum Principle for problem (29), we easily deduce that the Lagrange multiplier $\lambda$ associated with the cost $J_{\mathcal{E}}$ is necessarily non-zero. □

We observe that ‘normal form’ in Corollary 4.4 means that the Lagrange multiplier $\lambda$ associated with the cost $J_{\mathcal{E}}$ is non-zero. This does not necessarily implies that also the Lagrange multiplier $p_0$ associated with $q^0$ is non-zero. (In the next section we will see that if we impose a lower bound for the Copepod displacement, then we will obtain that both $\lambda \ne 0$ and $p_0 \ne 0$.)
## 4.2. Geometric efficiency and normality for the Copepod model

In this section we consider the following optimal control problem in which we are given a target $[x_T, x_T + d]$ with $x_T > 0$, for the Copepod displacement:

$$ 
\left\{
\begin{array}{l}
\text{Minimize } J_{\mathcal{E}}((q^0(\cdot), q(\cdot)), u(\cdot)) := \left[-\frac{(q_1(T))^2}{q^0(T)}\right] \\
\text{over arcs } (q^0(\cdot), q_1(\cdot), q_2(\cdot), q_3(\cdot)) \in W^{1,1}([0, T]; \mathbb{R}^4) \text{ s.t.} \\
(\dot{q}_1, \dot{q}_2, \dot{q}_3)(t) = F(q(t), u(t)), \quad \text{for a.e. } t \in [0, T] \\
\dot{q}^0(t) = \ell(q(t), u(t)), \quad \text{for a.e. } t \in [0, T] \\
u(t) \in U \quad \text{for a.e. } t \in [0, T] \\
q^0(0) = 0, \\
q_1(0) = 0, \quad q_1(T) \in [x_T, x_T + d] \\
q_2(0) = q_2(T), \quad q_3(0) = q_3(T),
\end{array}
\right.
\qquad (38)
$$

where $d \ge x_T(> 0)$. We are interested in establishing normality properties of the Copepod model (38). More precisely we shall prove that the optimal control problem (38) does not allow abnormal minimizers having displacement smaller than $x_T + d$. As a consequence, choosing suitable $x_T$ and $d$, we can restrict attention merely to normal extremals since the only possible candidate to be abnormal extremal (in $\mathcal{T}$) is the triangular stroke which has a specific displacement (cf. Proposition 2.6 and Corollary 2.7). This is valid in particular for the case of interest of this report, in which we consider two different Lagrangians:

(I) $\ell(q, u) = u_1^2 + u_2^2$ (euclidean cost).

(II) the case in which the Lagrangian represents the mechanical energy of the system, cf (3)

(In either case, the Lagrangian $\ell$ does not depend of $q_1$, and can be expressed in terms of regular functions of $\sin(q_2)$, $\sin(q_3)$, $\cos(q_2)$, $\cos(q_3)$.) Therefore, henceforth in this section, we shall assume that either (I) or (II) above are satisfied.

## The Maximum Principle

Consider the pseudo-Hamiltonian (also referred to as ‘unmaximized’ Hamiltonian) $H: \mathbb{R}^4 \times \mathbb{R}^4 \times \mathbb{R}^2 \to \mathbb{R}$:

$$ H(q, p, u) := \langle (p_1, p_2, p_3), F(q, u) \rangle - p_0 \ell(q, u). $$

The Maximum Principle for problem (38) takes the following form. Let $(\bar{q}, \bar{u})$ be a (local) minimizer for (38). Then, there exist a vector-valued function $p = (p_0, p_1, p_2, p_3) \in W^{1,1}([0, T]; \mathbb{R}^4)$, $\alpha_1 \ge 0$, $\alpha_2 \ge 0$, $\beta_1 \in \mathbb{R}$, $\beta_2 \in \mathbb{R}$ and $\lambda \ge 0$ such that

(i) $\lambda + \|p\|_{L^\infty} + |\beta_1| + |\beta_2| + \alpha_1 + \alpha_2 \ne 0$;

(ii) $-\dot{p}(t) = \partial_q H(\bar{q}(t), p(t), \bar{u}(t))$ a.e.;

(iii) $H(\bar{q}(t), p(t), \bar{u}(t)) = \max_{u \in U} H(\bar{q}(t), p(t), u) = r$ for a.e. $t$, for some $r \in \mathbb{R}$;

(iv) $p_2(0) = p_2(T) = \beta_1$, $p_3(0) = p_3(T) = \beta_2$, $p_0(T) = \lambda \frac{(q_1(T))^2}{(q^0(T))^2}$, and $p_1(T) = \alpha_1 - \alpha_2 + 2\lambda \frac{q_1(T)}{q^0(T)}$, with $\alpha_1 = 0$ if $\bar{q}_1(T) > x_T$ and $\alpha_2 = 0$ if $\bar{q}_1(T) < x_T + d$.

Observe that, since $F$ and $\ell$ do not depend on $q^0$ or $q_1$, condition (ii) provides the following relations:

$$ \dot{p}_0(t) = 0, \quad \dot{p}_1(t) = 0 \quad \text{a.e.} $$

(which yield $p_0(\cdot) \equiv \lambda \frac{(q_1(T))^2}{q^0(T)}$ and $p_1(\cdot) \equiv \alpha_1 - \alpha_2 + 2\lambda \frac{q_1(T)}{q^0(T)}$ and

$$ -\dot{p}_2(t) = p_1(\bar{u}_1(t))\partial_{q_2}\varphi_1(\bar{q}(t)) + \bar{u}_2(t)\partial_{q_2}\varphi_2(\bar{q}(t)) - p_0\partial_{q_2}\ell(\bar{q}(t), \bar{u}(t)) \quad \text{a.e. } t \in [0, T] \qquad (39) $$

$$ -\dot{p}_3(t) = p_1(\bar{u}_1(t))\partial_{q_3}\varphi_1(\bar{q}(t)) + \bar{u}_2(t)\partial_{q_3}\varphi_2(\bar{q}(t)) - p_0\partial_{q_3}\ell(\bar{q}(t), \bar{u}(t)) \quad \text{a.e. } t \in [0, T], \qquad (40) $$
in which

$$ \varphi_1(q) = \varphi_1(q_2, q_3) := \frac{\sin(q_2)}{2 + \sin^2(q_2) + \sin^2(q_3)}, \quad \varphi_2(q) = \varphi_2(q_2, q_3) := \frac{\sin(q_3)}{2 + \sin^2(q_2) + \sin^2(q_3)}. $$

**Definition 4.5.** We say that the (local) minimizer for (38) $(\bar{q}, \bar{u})$ is normal if the Maximum Principle applies with $\lambda \neq 0$.

**Remark 4.6.** Notice that for problem (38) we shall prove that if a minimizer $(\bar{q}, \bar{u})$ is such that $\bar{q}_1(T) \in [x_T, x_T + d[$, then both $\lambda \neq 0$ and $p_0 \neq 0$ (the multiplier $p_0$ is associated with the state variable $q^0$). Indeed, since $\lambda$ and $p_0$ turn out to be proportional, the minimizer $(\bar{q}, \bar{u})$ can be considered abnormal if the (first order) necessary conditions for optimality are applicable with $\lambda = 0$ or, equivalently, $p_0 = 0$.

**Proposition 4.7.** If a minimizer $(\bar{q}, \bar{u})$ for (38) is abnormal, then necessarily $\bar{q}_1(T) = x_T + d$.

*Proof. Step 1.* Let $(\bar{q}, \bar{u})$ be a minimizer for (38), and assume that the Maximum Principle applies with $\lambda = 0$ (i.e. $(\bar{q}, \bar{u})$ is abnormal). Then either $\bar{q}_1(T) = x_T$ or $\bar{q}_1(T) = x_T + d$. Indeed, supposing by contradiction that $\bar{q}_1(T) \in [x_T, x_T + d[$, and bearing in mind the transversality condition for the adjoint variable $p_1$, we would obtain that $\alpha_1 = \alpha_2 = 0$ and, therefore $p_1(.) \equiv 0$. On the other hand the Maximality condition (iii) yields:

$$ p_2(t) = -p_1\varphi_1(\bar{q}(t)), \quad p_3(t) = -p_1\varphi_2(\bar{q}(t)). \qquad (41) $$

Then from the formulae derived for $p_2(.)$, $p_3(.)$ we would also have $p(.) \equiv 0$, obtaining that $(p(.)$, $\lambda) = (0, 0)$. But, this would contradict the non-triviality condition (i) of the Maximum Principle. Then, we deduce that $\bar{q}_1(T) \in \{x_T, x_T + d\}$.

**Step 2.** We introduce a new optimal control problem in which we partially replace the right end-point constraint $q_1(T) \in [x_T, x_T + d]$ by a penalty term in the cost to minimize.

$$ \left\{ \begin{array}{l} \text{Minimize } \tilde{J}_{\mathcal{E}}((q^0(.), q(.)), u(.)) := \left[-\frac{(q_1(T))^2}{q^0(T)} + K \max\{-(q_1(T) - x_T); 0\}\right] \\ \text{over arcs } (q^0, q(.)) \in W^{1,1}([0, T]; \mathbb{R}^4) \text{ s.t.} \\ \dot{q}(t) = F(q(t), u(t)), \quad \text{for a.e. } t \in [0, T] \\ \dot{q}^0(t) = \ell(q(t), u(t)), \quad \text{for a.e. } t \in [0, T] \\ u(t) \in U \quad \text{for a.e. } t \in [0, T] \\ q^0(0) = 0, \quad q_1(0) = 0, q_1(T) \le x_T + d, \\ q_2(0) = q_2(T), \quad q_3(0) = q_3(T). \end{array} \right. \qquad (42) $$

Here we take $K > 0$ such that

$$ K > \frac{2C}{x_T}, $$

in which $C > 0$ is a constant such that

$$ \left| \frac{(q_1(T))^2}{q^0(T)} \right| \le C, \quad \text{for all trajectory-control pairs } (q(.), u(.)). $$

(Observe that such a constant $C$ always exists, take for instance $C = \frac{\sqrt{2\pi}}{2}$.)

We claim that if $(\bar{q}, \bar{u})$ is a minimizer for (38), then it is a minimizer also for (42). Indeed, if $(\bar{q}, \bar{u})$ is a minimizer for (38), then

$$ J_{\mathcal{E}}(\bar{q}(.), \bar{u}(.)) \le J_{\mathcal{E}}(q(.), u(.)), \qquad (43) $$
for all trajectory-control pair $(q(\cdot), u(\cdot))$.

Assume, by contradiction, that we can find a trajectory/control pair $(\hat{q}(\cdot), \hat{u}(\cdot))$ satisfying all the requirements of the control system in (42), such that

$$ \tilde{J}_{\mathcal{E}}(\hat{q}(\cdot), \hat{u}(\cdot)) < \tilde{J}_{\mathcal{E}}(\bar{q}(\cdot), \bar{u}(\cdot)). \quad (44) $$

From the choice of $K$ we necessarily have that $\tilde{q}_1(T) \in [x_T/2, x_T]$. Consider the trajectory/control pair $(\tilde{q}(\cdot), \tilde{u}(\cdot))$ satisfying the following properties:

$$ \tilde{u}(s) := \begin{cases} 2\hat{u}(2s) & \text{if } s \in [0, T/2] \\ 2\hat{u}(2s-T) & \text{if } s \in (T/2, T] \end{cases} \quad (45) $$

$\tilde{q}^0(0) = 0$, $\tilde{q}_1(0) = 0$, and

$$ (\tilde{q}_2, \tilde{q}_3)(s) = \begin{cases} (\hat{q}_2, \hat{q}_3)(2s) & \text{if } s \in [0, T/2] \\ 2(\hat{q}_2, \hat{q}_3)(2s - T) & \text{if } s \in (T/2, T]. \end{cases} \quad (46) $$

(Roughly speaking, employing a 'bigger' control, we construct a trajectory which, in the $(q_2, q_3)$-variables, does twice the path of $(\hat{q}(\cdot), \hat{u}(\cdot))$ on the same time interval.) As a consequence, we obtain $\tilde{q}_1(T) \in [x_T, 2x_T] \subset [x_T, x_T + d]$, and a straightforward calculation provides

$$ \frac{(\tilde{q}_1(T))^2}{\tilde{q}^0(T)} = \frac{(\hat{q}_1(T))^2}{\hat{q}^0(T)}. $$

We would deduce that

$$ \begin{aligned} J_{\mathcal{E}}(\tilde{q}(\cdot), \tilde{u}(\cdot)) &= -\frac{(\tilde{q}_1(T))^2}{\tilde{q}^0(T)} = -\frac{(\hat{q}_1(T))^2}{\hat{q}^0(T)} \le \tilde{J}_{\mathcal{E}}(\hat{q}(\cdot), \hat{u}(\cdot)) \\ &< \tilde{J}_{\mathcal{E}}(\bar{q}(\cdot), \bar{u}(\cdot)) = J_{\mathcal{E}}(\bar{q}(\cdot), \bar{u}(\cdot)), \end{aligned} $$

which contradicts the minimality of $(\bar{q}, \bar{u})$ for (38).

**Step 3.** We claim that the optimal control problem (42) has no abnormal minimizers. To see this consider the Maximum Principle for problem (42), which asserts that if $(\bar{q}, \bar{u})$ is an abnormal minimizer for (42), then, there exist a vector-valued function $p = (p_0, p_1, p_2, p_3) \in W^{1,1}([0, T]; \mathbb{R}^4)$, $\beta_1 \in \mathbb{R}$, $\beta_2 \in \mathbb{R}$, $\alpha_2 \ge 0$ and $\lambda \ge 0$ such that

(i)' $\lambda + ||p||_{L^\infty} + |\beta_1| + |\beta_2| \ne 0$;

(ii) $-\dot{p}(t) = \partial_q H(\bar{q}(t), p(t), \bar{u}(t))$ a.e.;

(iii) $H(\bar{q}(t), p(t), \bar{u}(t)) = \max_{u \in U} H(\bar{q}(t), p(t), u) = r$ for a.e. $t$, for some $r \in \mathbb{R}$;

(iv)' $p_2(0) = p_2(T) = \beta_1$, $p_3(0) = p_3(T) = \beta_2$, $p_0(T) = \lambda \frac{(\bar{q}_1(T))^2}{\bar{q}^0(T)}$, and $-p_1(T) \in \lambda(-2\frac{\bar{q}_1(T)}{\bar{q}^0(T)} + [-K, 0]) + \alpha_2$, where $\alpha_2 = 0$ if $\bar{q}_1(T) < x_T + d$.

From condition (ii) we know that $p_0(.)$ and $p_1(.)$ are constants, and system (39)-(40) are valid. Then, from the relations $(iv)'$ above we would deduce that $p_0 \equiv 0$ and $p_1 \equiv 0$, and therefore (39)-(40) yield also $p_2 \equiv 0$ and $p_3 \equiv 0$. This would mean that $(p(.)$, $\lambda) = (0, 0)$, which contradicts the non-triviality condition $(i)'$ of the Maximum Principle.

In conclusion, from step 2 above, if $(\bar{q}, \bar{u})$ is an abnormal minimizer for (38), then it would be an abnormal minimizer for (42) as well, but this is admissible only if $\bar{q}_1(T) = x_T + d$ owing to step 3. $\square$

Write $x_{ab} := \int_T d\omega$ the displacement of the (piecewise smooth) stroke confined in $T$ having the biggest displacement (cf. Proposition 2.6 and Corollary 2.7).
**Corollary 4.8.** Take $x_T \in ]0, x_{ab}[ \text{ and } d \ge x_T \text{ such that } x_T + d < 2x_{ab} \text{ and } x_T + d \ne x_{ab}$. Then each minimizer for (38) is normal.

*Proof.* Suppose, by contradiction, that a minimizer $(\bar{q}, \bar{u})$ is abnormal. Then, from Proposition 4.7 it follows that $\bar{q}_1(T) = x_T + d$. But, being abnormal, from the Maximum Principle (cf. in particular conditions (39)-(41)) we deduce also that $(\bar{q}, \bar{u})$ should satisfy also the equation:

$$ \sin(\bar{q}_2(t)) \sin(\bar{q}_3(t)) [\cos(\bar{q}_2(t)) - \cos(\bar{q}_3(t))] = 0. $$

As a consequence $\bar{q}(\cdot)$ would be the piecewise smooth arc staying on the boundary of the triangle $\mathcal{T}$, which would imply $\bar{q}_1(T) = x_{ab}$. But, $x_{ab} \ne x_T + d$ by assumption. A contradiction. □

## 5. CONCLUSIONS

In this article we have presented the key concepts and results of sub-Riemannian geometry providing a rich setting to study micro-swimmers. In particular the existence of center of strokes and their computations is related to compute an invariant of the problem. This analysis can be generalized to more complicated models e.g. the long-standing classical Purcell swimmer. A preliminary geometric analysis of this swimmer is initialized in [6] and several centers of swimming have been numerically detected. Tough simple, the geometric structure of the Purcell swimmer is very difficult to handle. Its nilpotent model has a growth vector (2, 3, 5) and is called the Cartan flat model. But the microlocal analysis is highly non-trivial, and the effect of perturbations from the nilpotent model to a generic model is multiple. The copepod swimmer is a simpler model which can be interpreted as a limit case of the symmetric Purcell swimmer. It can be seen as an embedding of the Heisenberg and Martinet case in the Cartan flat model. It constitutes a favourable setting to develop geometric and analytical tools to better understand swimmer model at low Reynolds number. The design of micro-swimmers has widespread applications (for instance in medicine), and models can be studied at a large scale using viscous fluids. Experiments on the copepod prototype are currently performed in Pr. D. Takagi's Hawaii laboratory, where observations of the (real) swimming mechanism of this zooplankton are conducted in parallel.

## REFERENCES

[1] A. Agrachev, J.P. Gauthier On the Dido problem and plane isoperimetric problems, Acta Appl. Math. 57, 3 (1999) 287-338

[2] El-H.C. Alaoui, Les métriques sous riemanniennes en dimension 3, Ph.D. thesis, Université de Rouen, France, 1996

[3] El-H.C. Alaoui, J.P. Gauthier, I. Kupka, Small sub-Riemannian balls on $R^3$, J. Dynam. Control Systems. 2, 3 (1996) 359-421

[4] J. E. Avron, and O. Raz, A geometric theory of swimming: Purcell's swimmer and its symmetrized cousin, New Journal of Physics 10, 6 (2008): 063016.

[5] A. Bellaïche, The tangent space in sub-Riemannian geometry, J. Math. Sci. (New York) 83, 4 (1997) 461-476

[6] P. Bettiol, B. Bonnard, J. Rouot, Optimal strokes at low Reynolds number: a geometric and numerical study of Copepod and Purcell swimmers, Submitted.

[7] B. Bonnard, J.-B. Caillaud, E. Trélat, Second order optimality conditions in the smooth case and applications in optimal control, ESAIM Control Optim. Calc. Var. 13 (2007) 207-236

[8] B. Bonnard, G. Charlot, R. Ghezzi, G. Janin, The sphere and the cut locus at a tangency point in two dimensional almost-Riemannian geometry, J. Dynam. Control Systems. 17 (I) (2011) 141-161

[9] B. Bonnard, M. Chyba, Singular trajectories and their role in control theory, Mathématiques & Applications 40, Springer-Verlag, Berlin (2003)

[10] R.W. Brockett, Control theory and singular Riemannian geometry, New directions in applied mathematics (Cleveland, Ohio, 1980) Springer, New York-Berlin (2016) 11-27.

[11] T. Chambrion, L. Giraldi and A. Munnier, Optimal strokes for driftless swimmers: A general geometric approach, Accepted in ESAIM Control Optim. Calc. Var. (2017)

[12] O. Cots, Contrôle optimal géométrique : méthodes homotopiques et applications, PhD thesis, Institut Mathématiques de Bourgogne, Dijon, France (2012)

[13] I. Kupka, Géométrie sous-riemannienne, Astérisque, Séminaire Bourbaki 1995/96, (1997) 351-380

[14] P.H. Lenz, D. Takagi, D.K. Hartline, Choreographed swimming of copepod nauplii, Journal of The Royal Society Interface 12, 112 20150776 (2015)
[15] M.J. Lighthill, Note on the swimming of slender fish, *J. Fluid Mech.* **9** (1960) 305–317

[16] R. Montgomery, Isoholonomic problems and some applications, *Commun. Math. Phys.* **128**, 3 (1990) 565–592.

[17] E. Passov, Y. Or, Supplementary notes to: Dynamics of Purcells three-link microswimmer with a passive elastic tail, *EPJ E* **35** (2012) 1–9

[18] L.S. Pontryagin, V.G. Boltyanskii, R.V. Gamkrelidze and E.F. Mishchenko, The mathematical theory of optimal processes, *Interscience Publishers John Wiley & Sons, Inc.* New York-London (1962)

[19] E.M. Purcell, Life at low Reynolds number, *Am. J. Phys.* **45** (1977) 3–11

[20] J. Rouot, *Méthodes géométriques et numériques en contrôle optimal et applications au transfert orbital poussée faible et à la nage à faible nombre de Reynolds.*, Ph.D. thesis, Inria Sophia Antipolis Méditerranée, France (2016)

[21] J. Rouot, P. Bettiol, B. Bonnard, A. Nolot. Optimal control theory and the efficiency of the swimming mechanism of the Copepod Zooplankton, *Proc. 20th IFAC World Congress*, Toulouse (2017)

[22] D. Takagi, Swimming with stiff legs at low Reynolds number, *Phys. Rev. E* **92**, (2015)

[23] R.B. Vinter, Optimal control, *Systems & Control: Foundations & Applications* (2000) xviii–507

[24] Q. Wang and J.L. Speyer, Necessary and sufficient conditions for local optimality of a periodic process, *SIAM J. Control Optim.* n.2 **28** (1990) 482–497