# Optimal Strokes: a Geometric and Numerical Study of the Copepod Swimmer

Bernard Bonnard, Monique Chyba, Daisuke Takagi, Jeremy Rouot, Rong Zou

<figure><img src="image_1.png" /></figure>

To cite this version:

Bernard Bonnard, Monique Chyba, Daisuke Takagi, Jeremy Rouot, Rong Zou. Optimal Strokes: a Geometric and Numerical Study of the Copepod Swimmer. 2015. hal-01162407v3

HAL Id: hal-01162407

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

Submitted on 8 Jan 2016

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.

Public Domain
OPTIMAL STROKES : A GEOMETRIC AND NUMERICAL
STUDY OF THE COPEPOD SWIMMER

BERNARD BONNARD, MONIQUE CHYBA, JÉRÉMY ROUOT, DAISUKE TAKAGI,
RONG ZOU

BERNARD BONNARD

Inria Sophia Antipolis et Institut de Mathématiques de Bourgogne
9 avenue Savary
21078 Dijon, France

MONIQUE CHYBA, DAISUKE TAKAGI

2565 McCarthy the Mall
Department of Mathematics
University of Hawaii
Honolulu, HI 96822, USA

JÉRÉMY ROUOT

Inria Sophia Antipolis
2004 route des lucioles
06902 Sophia Antipolis, France

RONG ZOU

Department of Systems Innovation,
Graduate School of Engineering Science
Osaka, Japan

ABSTRACT. The objective of this article is to make a geometric and numerical analysis concerning the optimal displacements of a larval copepod swimming at low Reynolds number. A simplified model of locomotion is analyzed in the framework of Sub-Riemannian geometry. In particular, the role of both normal and abnormal geodesics is related to observed geometric motions in relation with the mechanical power dissipated by the swimmer.

1. **Introduction.** Swimming microorganisms employ a variety of mechanisms of propulsion, and they have inspired numerous models starting with undulating sheets and filaments introduced in the fifties [24, 15]. Different types of strokes are observed depending upon the circumstances which can be used to design micro-robots in relation with the purpose, e.g. military or medical applications.

<footnote n="2010">2010 Mathematics Subject Classification. 49K15, 93C10, 70Q05.</footnote>

<footnote>*Key words and phrases.* Low Reynolds number, Copepod swimmer, SR-geometry, periodic optimal control.</footnote>

<footnote>B. Bonnard and M. Chyba are partially supported by the National Science Foundation (NSF) Division of Mathematical Sciences, award #1109937, J. Rouot is supported by the French Space Agency CNES, R&T action R-S13/BS-005-012 and by the région Provence-Alpes-Côte d'Azur, and R. Zou is partially supported by the National Science Foundation (NSF) Division of Division Of Graduate Education, award #0841223.</footnote>
Recent studies have explored optimal strategies for swimming with minimal amount of mechanical work, an important criterion for assessing the fitness of different organisms and for designing efficient robotic swimmers [17]. Previous studies have computed optimal solutions in the framework of variational analysis or optimal control [3],[5],[23],[6],[12].

A well-studied model is the Purcell swimmer [21], which consists of three rigid and slender rods representing respectively the leg, the body, and the arm. The configuration of the swimmer is described by two angles $\theta = (\theta_1, \theta_2)$ with three other variables $q = (x, y, \Phi)$ representing respectively the position and the orientation of the body. The system can be written as

$$ \dot{q} = D(\Phi)G(\theta)\dot{\theta}, $$

$$ \dot{\theta} = H(\theta)\tau, $$

where $D(\Phi)$ is the rotation matrix

$$ D(\Phi) = \begin{pmatrix} \cos(\Phi) & -\sin(\Phi) & 0 \\ \sin(\Phi) & \cos(\Phi) & 0 \\ 0 & 0 & 1 \end{pmatrix}. $$

By denoting $u = \dot{\theta}$, the mechanical power is

$$ \tau u = uH^{-1}u $$

and the energy minimization problem becomes

$$ \int_{0}^{T} (uH^{-1}u)dt $$

where $u$ is taken as the control variable (see [19] for a complete description of the matrices $G$ and $H$). This control problem falls into the Sub-Riemannian (SR) framework [18]. It is complex even locally and is related to the Cartan flat SR-model [6]. A simplified cost is to replace $\tau$ by $u$ (neglecting the fluid interaction) and leads to minimize the integral:

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

The simplest SR-case called the Heisenberg problem provides us with some informations in the swimming problem. The Heisenberg control system is given by

$$ \dot{x} = u_1\theta_2 - u_2\theta_1 $$

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

while minimizing the integral

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

The integration of the non planar geodesics starting from the origin gives

$$ \theta_{1}(t) = \frac{A}{\lambda}(\sin(\lambda t + \varphi) - \sin(\varphi)) $$

$$ \theta_{2}(t) = \frac{A}{\lambda}(\cos(\lambda t + \varphi) - \cos(\varphi)) $$

$$ x(t) = \frac{A^2}{\lambda}t - \frac{A^2}{\lambda^2}\sin(\lambda t) $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Sketch of a symmetric swimmer consisting of two pairs of legs.</figcaption></figure>

where A, λ and ϕ are parameters related to the initial velocity. Both angular variables are periodic motion corresponding to strokes to provide the displacement x(2π/λ), whose average is given by A²/λ², the period of the stroke being 2π/λ. A standard computation of conjugate points shows that such geodesic is optimal up to t = 2π/λ (included). In this example, the optimal displacement is shown to correspond to a simple stroke: indeed the first conjugate time is at t = 2π/λ and moreover it corresponds to the cut time.

While this model can provide some insights on optimal locomotion, it is too primitive because:

1. The geodesic flow is integrable due to a symmetry of revolution along Ox and every θ-motion is periodic.

2. The model is quasi-homogeneous, θ₁ and θ₂ are of weight 1 and x is of weight 2, and invariant in the Heisenberg group.

Recently a new model was developed to mimic the locomotion of larval copepods, an abundant type of zooplankton thriving in the ocean [22, 20]. The simplest form of the model, hereafter referred to as the copepod swimmer, is a symmetric body consisting of two pairs of legs, with the first pair making an angle θ₁ and the second pair making an angle θ₂ with respect to the displacement direction Ox (Fig.1).

The swimming velocity at x₀ is given by

$$ \dot{x}_0 = \frac{\dot{\theta}_1 \sin(\theta_1) + \dot{\theta}_2 \sin(\theta_2)}{2 + \sin^2(\theta_1) + \sin^2(\theta_2)} \quad (1) $$

and the controls are the angular velocities

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

We also have the state constraints $\theta_i \in [0, \pi]$, $i = 1, 2$, $\theta_1 \le \theta_2$.

A simplified cost neglecting the fluid interaction can be identified as

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

but the true cost corresponding to the mechanical energy of the system is given by the quadratic form

$$ \dot{q}^T M \dot{q}, \quad 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} \quad (2) $$
where $q = (x_0, \theta_1, \theta_2)$.

Using (1), this amounts to minimize the quadratic cost

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

with

$$ a = \frac{1}{3} - \frac{\sin^2 \theta_1}{2(2 + \sin^2 \theta_1 + \sin^2 \theta_2)}, $$

$$ b = -\frac{\sin \theta_1 \sin \theta_2}{2(2 + \sin^2 \theta_1 + \sin^2 \theta_2)}, $$

$$ c = \frac{1}{3} - \frac{\sin^2 \theta_2}{2(2 + \sin^2 \theta_1 + \sin^2 \theta_2)}. $$

This copepod swimmer serves as a suitable model for computing optimal controls in the framework of SR-geometry. The system is three-dimensional (two controls and one variable), which is arguably simpler than the five-dimensional system (two controls and three variables) of the previously studied Purcell swimmer where moreover the expression of the control fields are complicated see [19]. It is a global model of SR-geometry which can be analyzed in detail, showing in particular the role of normal and abnormal geodesics in the motion. In addition, the optimal controls could be compared with observations of copepods to determine whether they are optimizing their strokes to minimize energy. Copepods must swim in order to find food and escape from predators, and they have had a chance to adapt and evolve over millions of years, but it remains unknown to what extent they have adapted their strokes to maximize their swimming efficiency. Thus the model optimization could offer new insight into biological behaviour.

This article is organized in two sections. In section 2, we recall some properties of the copepod swimmer [22] and the mathematical tools from geometric optimal control (see [8] for a general reference). Section 3 contains the contribution of this article based on a geometric analysis and numerical simulations to describe the normal strokes in relation with the classification of periodic planar curves [4]. Finally, the optimal strokes satisfying the constraints are numerically computed using the two softwares: Bocop(www.bocop.org, [7]) and HamPath(http://cots.perso.enseeiht.fr/hampath/, [13]).

## 2. Preliminary results.

### 2.1. Geometric analysis of a copepod swimmer.

A (general) stroke of period $T$ consists in a periodic motion in the shape variables $(\theta_1, \theta_2)$. Assuming $x_0(0) = 0$, the corresponding displacement is $x_0(T)$.

In [22], two types of geometric motions are described:

**First case:** (Fig.2) The two legs are assumed to oscillate sinusoidally with period $2\pi$ according to

$$ \theta_1 = \Phi_1 + a \cos(t), \quad \theta_2 = \Phi_2 + a \cos(t + k_2) $$

with $a = \pi/4$, $\Phi_1 = \pi/4$, $\Phi_2 = 3\pi/4$ and $k_2 = \pi/2$. This produces a displacement $x_0(2\pi) = 0.2$.

Parameters $a, \Phi_1, \Phi_2$ and $k$ are designed to maximize the efficiency.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Two legs oscillating sinusoidally according to $\theta_1 = \pi/4 + a \cos t$ and $\theta_2 = 3\pi/4 + a \cos(t + \pi/2)$, where $a = \pi/4$ is the amplitude. The second leg (blue) oscillates about $\Phi_2 = 3\pi/4$, while the first leg (red) oscillates about $\Phi_1 = \pi/4$ with a phase lag of $\pi/2$. The swimmer position $x_0$ translates about a fifth of the leg length after one cycle.</figcaption></figure>

**Second case:** (Fig. 3) The two legs are paddling in sequence followed by a recovery stroke performed in unison. In this case the controls $u_1 = \dot{\theta}_1$, $u_2 = \dot{\theta}_2$ produce bang arcs to steer the angles between from the boundary $0$ of the domain to the boundary $\pi$, while the unison sequence corresponds to a displacement from $\pi$ to $0$ with the constraint $\theta_1 = \theta_2$.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Two legs paddling in sequence. The legs perform power strokes in sequence and then a recovery stroke in unison, each stroke sweeping an angle $\pi$.</figcaption></figure>
Our main objective is to relate these policies to geometric optimal control.

## 2.2. Abnormal curves in the copepod swimmer.

We introduce $q = (x_0, \theta_1, \theta_2)$, then the system is written as a driftless affine control system

$$ \dot{q}(t) = \sum_{i=1}^{2} u_i(t) F_i(q(t)) $$

where the control vector fields are given by

$$ F_i = \frac{\sin(\theta_i)}{\Delta} \frac{\partial}{\partial x_0} + \frac{\partial}{\partial \theta_i} $$

with $\Delta = 2 + \sin^2(\theta_1) + \sin^2(\theta_2)$. We denote by $D$ the distribution generated by the two vector fields: $D = \operatorname{span}\{F_1, F_2\}$.

The Lie bracket of two vector fields $F, G$ is computed with the convention

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

Finally, we denote by $p = (p_1, p_2, p_3)$ the adjoint vector associated with $q$.

We first recall basic facts concerning the local classification of two-dimensional distributions, in relation with abnormal curves.

### 2.2.1. Local classification of two-dimensional distributions in dimension three and abnormal curves.

Let $D = \operatorname{span}\{G_1, G_2\}$ be the distribution generated by two vectors fields $G_1, G_2$ in $\mathbb{R}^3$. Let $z = (q, p)$ and denote $H_i(z) = \langle p, G_i(q) \rangle$, $i = 1, 2$ the Hamiltonian lifts. The Poisson bracket is given by

$$ \{H_1, H_2\}(z) = dH_1(\vec{H}_2)(z) = \langle p, [G_1, G_2](q) \rangle. $$

Abnormal curves are defined by

$$ H_1(z) = H_2(z) = 0, $$

and differentiating using the dynamics

$$ \frac{dz}{dt} = \sum_{i=1}^{2} u_i \vec{H}_i(z) $$

we obtain the relations

$$ \{H_1, H_2\}(z) = 0 \\ u_1\{\{H_1, H_2\}, H_1\}(z) + u_2\{\{H_1, H_2\}, H_2\}(z) = 0 $$

defining the corresponding abnormal controls.

Tools from singularity theory can be used to classify the distributions, see [25]. Here we present only the two (stable) models related to our study.

**Contact case.** We say that $q_0$ is a contact point if $\{G_1, G_2, [G_1, G_2]\}$ is of dimension 3 at $q_0$. At a contact point, identified to 0, there exists a system of local coordinates $q = (x, y, z)$ such that

$$ D = \ker(\alpha), \quad \alpha = ydx + dz. $$

Observe that $d\alpha = dy \wedge dx$ (Darboux form) and that $\frac{\partial}{\partial z}$ is the characteristic direction of $d\alpha$. This form is equivalent to

$$ D = \ker(\alpha'), \quad \alpha' = dz + (xdy - ydx). $$
with

$$ 
\begin{align} 
D &= \operatorname{span}\{G_1, G_2\}, & G_1 &= \frac{\partial}{\partial x} + y \frac{\partial}{\partial z}, \tag{4} \\ 
G_2 &= \frac{\partial}{\partial y} - x \frac{\partial}{\partial z}, & G_3 &= [G_1, G_2] = 2 \frac{\partial}{\partial z}. \nonumber
\end{align} 
$$

**The Martinet case.** A point $q_0$ is a Martinet point if at $q_0$, $[G_1, G_2] \in \operatorname{span}\{G_1, G_2\}$ and at least one Lie bracket $[[G_1, G_2], G_1]$ or $[[G_1, G_2], G_2]$ does not belong to $D$. Then, there exist local coordinates $q = (x, y, z)$ near $q_0$ identified to 0 such that

$$ D = \ker \omega, \quad \omega = dz - \frac{y^2}{2} dx $$

where

$$ 
\begin{align} 
G_1 &= \frac{\partial}{\partial x} + \frac{y^2}{2} \frac{\partial}{\partial z}, & G_2 &= \frac{\partial}{\partial y}, & G_3 &= [G_1, G_2] = y \frac{\partial}{\partial z} \tag{5} \\ 
[[G_1, G_2], G_1] &= 0, & [[G_1, G_2], G_2] &= \frac{\partial}{\partial z}. \nonumber 
\end{align} 
$$

The surface $\Sigma: y = 0$ is called the *Martinet surface* and is foliated by abnormal curves, solutions of $\frac{\partial}{\partial x}$. In particular, through the origin it corresponds to the curve $t \to (t, 0, 0)$.

### 2.2.2. Computations in the copepod case.

We have

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

with

$$ 
\begin{align} 
f(\theta_1, \theta_2) &= \frac{2 \sin(\theta_1) \sin(\theta_2) (\cos(\theta_1) - \cos(\theta_2))}{\Delta^2}, \nonumber \\ 
[[F_1, F_2], F_1] &= \frac{\partial f}{\partial \theta_1}(\theta_1, \theta_2) \frac{\partial}{\partial x_0}, \quad [[F_1, F_2], F_2] = \frac{\partial f}{\partial \theta_2}(\theta_1, \theta_2) \frac{\partial}{\partial x_0}. \nonumber 
\end{align} 
$$

In particular we deduce the following lemma.

**Lemma 2.1.** The singular set $\Sigma: \{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 collinear, is given by $2 \sin(\theta_1) \sin(\theta_2)(\cos(\theta_1) - \cos(\theta_2)) = 0$ which corresponds to

* $\theta_1 = 0$ or $\pi$,

* $\theta_2 = 0$ or $\pi$,

* $\theta_1 = \theta_2$.

It is formed by the boundary of the physical domain: $\theta_i \in [0, \pi]$, $\theta_1 \le \theta_2$, with respective controls $u_1 = 0$, $u_2 = 0$ or $u_1 = u_2$.

**Remark 1.** The previous lemma provides the interpretation of the policy represented in Fig.3. In the shape space $(\theta_1, \theta_2)$ it corresponds to a *triangle*. The edges of the triangle are abnormal curves (where by definition the linearized system is not controllable).

**Remark 2.** A recent contribution [16] proves that a minimizer with such a corner cannot be optimal.

**Remark 3.** Using the results from [1], [9] allow to study the integrability properties and the optimality issues for the strokes with small amplitude.

To analyse the first situation of Fig.2, the mechanical energy has to be used in relation with SR-geometry.
## 2.3. Sub-Riemannian geometry.

The problem is written

$$ \dot{q} = \sum_{i=1}^{2} u_i G_i(q), \quad \min_{u(\cdot)} \int_{0}^{T} (u_1^2 + u_2^2) dt, $$

where the cost is defined for a fixed final time $T$ and corresponds to the energy. In this representation, we assume that the vector fields $G_1, G_2$ are orthonormal.

**Remark 4.** For the swimming copepod swimmer an orthonormal frame can be computed as followed. Using the following feedback transformation

$$ \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} = \begin{pmatrix} \cos(\alpha) & \sin(\alpha) \\ -\sin(\alpha) & \cos(\alpha) \end{pmatrix} \begin{pmatrix} v_1 \\ v_2 \end{pmatrix} $$

where $\alpha = \arctan \begin{pmatrix} \sin(\theta_1) \\ \sin(\theta_2) \end{pmatrix}$, the mechanical cost takes the form

$$ \frac{1}{3} v_1^2 + \frac{1}{6} \frac{2+\cos^2(\theta_2)+\cos^2(\theta_1)}{4-\cos^2(\theta_1)-\cos^2(\theta_2)} v_2^2. $$

Introducing $w_1 = \frac{1}{\sqrt{3}}v_1$, $w_2 = \sqrt{\frac{1}{6} \frac{2+\cos^2(\theta_2)+\cos^2(\theta_1)}{4-\cos^2(\theta_1)-\cos^2(\theta_2)}} v_2$ one gets the metric $w_1^2 + w_2^2$.

The admissible controls are bounded measurable mappings. According to Pontryagin maximum principle, we introduce the pseudo-Hamiltonian in the normal case

$$ H(z, u) = \sum_{i=1}^{2} u_i H_i(z) - \frac{1}{2} \sum_{i=1}^{2} u_i^2, $$

where $H_i$ are the Hamiltonian lifts $\langle p, G_i(q) \rangle$. The maximization condition is equivalent to $\frac{\partial H}{\partial u_i} = 0, i = 1, 2$. It follows that $u_i = H_i$ and plugging this expression for $u_i$ into $H$ produces the true Hamiltonian in the normal case

$$ H_n = \frac{1}{2} (H_1^2 + H_2^2). $$

In the contact situation, it corresponds to the Heisenberg case while in the Martinet situation, it corresponds to the flat Martinet case. In both cases it amounts to impose that $G_1, G_2$ are orthonormal and that the associated distribution is nilpotent.

**Definition 2.2.** A normal stroke is a solution of $\overrightarrow{H}_n$ such that $\theta_1$ and $\theta_2$ are periodic with period $T$.

According to the transversality conditions of the maximum principle the dual variables $p_2$ and $p_3$ are such that $p_2$ and $p_3$ are both periodic of period $T$ (to produce a smooth solution).

**Second order optimality condition.** In the normal case, the first conjugate point corresponds to the first point where a normal geodesic ceases to be minimizing which respect to the $C^1$-topology on the set of curves and they can be computed using the HamPath software [13].

This leads to the following definition.

**Definition 2.3.** A normal stroke is called $C^1$-optimal on $[0, T]$ if there exists no conjugate point on the interval $[0, T]$.

## 3. Computations and analysis in the copepod swimmer.
### 3.1. Simplified cost.

We start by considering the simplified cost

$$ \min_{u(\cdot)} \int_{0}^{T} (u_{1}^{2} + u_{2}^{2}) dt $$

in relation with the contact case.

Outside the singular set $\Sigma$, we have only contact points. Introducing the Hamiltonian lifts: $H_i = \langle p, F_i(q) \rangle$ for $i = 1, 2$, and $H_3 = \langle p, [F_1, F_2](q) \rangle$, the set $\{q, H_1, H_2, H_3\}$ are coordinates and we have

$$ H_1 = \frac{p_1 \sin \theta_1}{\Delta} + p_2, \quad H_2 = \frac{p_1 \sin \theta_2}{\Delta} + p_3, \qquad (6) $$

$$ H_3 = \frac{2p_1 \sin \theta_1 \sin \theta_2 (\cos \theta_1 - \cos \theta_2)}{\Delta^2}. $$

Moreover the problem is isoperimetric since $x_0$ is a cyclic coordinate, i.e. $p_1$ is a first integral: $\dot{p}_1 = 0$.

Straightforward computations lead to the expressions

$$ \dot{H}_1 = dH_1(\vec{H}_n) = dH_1 \left( \frac{1}{2} \overrightarrow{(H_1^2 + H_2^2)} \right) = \{H_1, H_2\} H_2, $$

$$ \dot{H}_2 = dH_2(\vec{H}_n) = dH_1 \left( \frac{1}{2} \overrightarrow{(H_1^2 + H_2^2)} \right) = \{H_2, H_1\} H_1. $$

This can be expressed in the following condensed way:

$$ \dot{H}_1 = H_2 H_3, \quad \dot{H}_2 = -H_1 H_3. $$

Moreover,

$$ \dot{H}_3 = dH_3(\vec{H}_n) = dH_3(\frac{1}{2}(\overrightarrow{H_1^2} + \overrightarrow{H_2^2})) = \{H_3, H_1\} H_1 + \{H_3, H_2\} H_2 $$

with

$$ \{H_3, H_1\}(z) = \langle p, [[F_1, F_2], F_1](q) \rangle, \quad \{H_3, H_2\}(z) = \langle p, [[F_1, F_2], F_2](q) \rangle. $$

At a contact point $\{F_1, F_2, F_3\}$ forms a frame, therefore

$$ [[F_1, F_2], F_1](q) = \sum_{i=1}^{3} \lambda_i(q) F_i(q) $$

and computing one gets,

$$ \lambda_1 = \lambda_2 = 0, \quad \frac{\partial f}{\partial \theta_1} = \lambda_3 f. $$

Similarly,

$$ [[F_1, F_2], F_2](q) = \sum_{i=1}^{3} \lambda'_{i}(q) F_{i}(q), $$

with

$$ \lambda'_{1} = \lambda'_{2} = 0, \quad \frac{\partial f}{\partial \theta_{2}} = \lambda'_{3} f. $$

We conclude that

$$ \begin{aligned} \dot{H}_1 &= H_2 H_3, \quad \dot{H}_2 = -H_1 H_3, \\ \dot{H}_3 &= H_3 (\lambda_3 H_1 + \lambda'_3 H_2). \end{aligned} \qquad (7) $$
We introduce a new time reparameterization with $ds = H_3dt$, and we obtain

$$ \frac{dH_1}{ds} = H_2, \quad \frac{dH_2}{ds} = -H_1, \quad \frac{dH_3}{ds} = \lambda_3 H_1 + \lambda'_3 H_2. $$

Hence we have the harmonic oscillator since $H_1'' + H_1 = 0$ when differentiating with respect to the new time $s$.

Furthermore $H_3$ can be computed using the remaining equation (7). Observe that with the approximation $\lambda_3, \lambda'_3$ constant, the equation is

$$ \frac{dH_3}{ds} = A \cos(s + \rho). $$

In this formalism, the Heisenberg case corresponds to $\lambda_3 = \lambda'_3 = 0$.

## 3.2. True Cost.

In this case the direct computation is more intricate since $F_1, F_2$ are not orthonormal.

The pseudo-Hamiltonian in the normal case becomes

$$ H(q,p) = u_1 H_1(q,p) + u_2 H_2(q,p) - \frac{1}{2} \left( a(q)u_1^2 + 2b(q)u_1u_2 + c(q)u_2^2 \right). $$

and the normal controls are computed solving the equations

$$ \frac{\partial H}{\partial u_1} = 0, \quad \frac{\partial H}{\partial u_2} = 0. $$

From an easy calculation, we obtain

$$ u_1 = - \frac{3(4H_1 + 2H_1 \sin^2 \theta_1 + 3H_2 \sin \theta_1 \sin \theta_2 - H_1 \sin^2 \theta_2)}{\sin^2 \theta_1 + \sin^2 \theta_2 - 4}, $$

$$ u_2 = - \frac{9H_1 \sin \theta_1 \sin \theta_2 + 6H_2(2 + \sin^2 \theta_2) - 3H_2 \sin^2 \theta_1}{\sin^2 \theta_1 + \sin^2 \theta_2 - 4}. $$

Plugging such $u$ into the pseudo-Hamiltonian gives the true Hamiltonian

$$ \begin{aligned} H_n &= u_1 H_1 + u_2 H_2 - \frac{1}{2} \left( a(q)u_1^2 + 2b(q)u_1u_2 + c(q)u_2^2 \right) \\ &= \frac{1}{2} \left( a(q)u_1^2 + 2b(q)u_1u_2 + c(q)u_2^2 \right) \\ &= -\frac{3}{2} \left( 2(\sin^2 \theta_1 + \sin^2 \theta_2)p_1^2 + (4 \sin^3 \theta_1 + 4 \sin^2 \theta_2 \sin \theta_1 + 8 \sin \theta_1)p_1 p_2 \right. \\ &\quad \left. + (2 \sin^4 \theta_1 + \sin^2 \theta_1 \sin^2 \theta_2 - \sin^4 \theta_2 + 8 \sin^2 \theta_1 + 2 \sin^2 \theta_2 + 8)p_2^2 \right. \\ &\quad \left. + (-\sin^4 \theta_1 + \sin^2 \theta_1 \sin^2 \theta_2 + 2 \sin^4 \theta_2 + 2 \sin^2 \theta_1 + 8 \sin^2 \theta_2 + 8)p_3^2 \right. \\ &\quad \left. + (4 \sin^2 \theta_1 \sin \theta_2 + 4 \sin^3 \theta_2 + 8 \sin \theta_2)p_1 p_3 + (6 \sin \theta_2 \sin^3 \theta_1 + 6 \sin^3 \theta_2 \sin \theta_1 \right. \\ &\quad \left. + 12 \sin \theta_1 \sin \theta_2)p_2 p_3 \right) / \left( (\sin^2 \theta_1 + \sin^2 \theta_2 - 4)(\sin^2 \theta_1 + \sin^2 \theta_2 + 2) \right) \end{aligned} \quad (8) $$

Again the problem is isoperimetric and $p_1$ is a first integral.
### 3.3. Numerical computations.

The period T is fixed to $2\pi$ in our simulations. We use the HamPath software [13] at two levels:

1. The shooting equations associated with the problem are

$$ 
\begin{aligned}
x_0(0) &= 0, & x_0(2\pi) &= x_f, \\
\theta_{1|2}(0) &= \theta_{1|2}(2\pi), & p_{2|3}(0) &= p_{2|3}(2\pi).
\end{aligned}
$$

2. The Hampath code is also used to show that the normal stroke is optimal testing the nonexistence of conjugate points using the variational equation to compute Jacobi fields. Recall that according to [8], given a reference curve $(q(t), p(t))$ solution of $\vec{H}_n$, a time $t_c \in [0, 2\pi]$ is a conjugate time if there exists a Jacobi field $\delta z = (\delta q, \delta p)$, that is a non-zero solution of the variational equation

$$ \dot{\delta z}(t) = \frac{\partial \vec{H}_n}{\partial z}(q(t), p(t))\delta z(t) \qquad (9) $$

such that $\delta q(0) = \delta q(t_c) = 0$. We denote $\delta z_i = (\delta q_i, \delta p_i)$, $i = 1...n$, $n$-independent solutions of (9) with initial condition $\delta q(0) = 0$. At time $t_c$ we have the following rank condition

$$ \mathrm{rank}\{\delta q_1(t_c), \dots, \delta q_n(t_c)\} < n. \qquad (10) $$

#### 3.3.1. Commented numerical results.
We present a sequence of numerical simulations for both costs, not taking into account the angular constraints. In the two cases we obtain similar results.

* A sequence of three identical strokes is represented on Fig. 4-5 for the two costs where $x_0(2\pi) = 0.2$. Numerical simulations show the non existence of conjugate points for both cases.

* Fig. 6-7-8 illustrate three different strokes confirming the complexity of the model and are related to the generic classification of periodic planar curves [4]. More complicated choreographies can be obtained, see Fig. 9. Conjugate points are also computed to check the second order optimality conditions. There is no conjugate points on $[0, 2\pi]$ in the case of the simple loop whereas they appear for the limaçon case, the eight case and more complicated cases. Hence, the only candidates for optimality are the simple loops.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. A sequence of three identical strokes for the simplified cost $\int_0^{2\pi} (u_1^2 + u_2^2)dt$ (state, adjoint and control variables).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 5. A sequence of three identical strokes for the mechanical cost (state, adjoint and control variables).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Normal stroke for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost (left) and the mechanical cost (right): simple loop.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Normal stroke for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost (left) and the mechanical cost (right): limaçon with inner loop.</figcaption></figure>

### 3.3.2. Optimal curves circumscribed in the triangle of constraints.
We use a combination of the Bocop and HamPath softwares.

**Bocop** software: This software is suitable to take into account de state constraints on the shape variables. Fig. 10 gives numerical simulations with this software, describing a creeping normal stroke in accordance with the abnormal triangle policy.
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Normal stroke for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ (left) and the mechanical cost (right): eight case.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. More complex normal strokes for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost.</figcaption></figure>

Fig.11 describes a single loop tangent to the boundary which is used to initialize the shooting algorithm of the HamPath software.

HamPath software: This software cannot be directly applied to compute the optimal solution using the Maximum Principle with state constraints, due to the complexity of the different principles [11].

Fig.12 describes a normal stroke tangent to the boundary.

### 3.3.3. Comparisons of the geometric efficiency of the strokes.
To compare the different normal and abnormal solutions corresponding to different displacements and in relation with the SR-interpretation we represent the ratio $L/x_0$ where L is the length of the stroke and $x_0$ is the corresponding displacement (such a quantity is not depending upon the parameterization).

For the triangle, a displacement along the vertical or horizontal edge gives $x_0 = \frac{2\sqrt{3}}{3} \text{arctanh}(\frac{\sqrt{3}}{3})$ and along the hypotenuse $x_0 = -\sqrt{2} \text{arctanh}(\frac{\sqrt{2}}{2})$ and the total displacement is $2.742 \cdot 10^{-1}$.

The length of a normal stroke $\gamma$ is $L(\gamma) = \int_0^{2\pi} \sqrt{\langle \dot{q}, \dot{q} \rangle} dt$ and easily computed using the energy level $H_n = \frac{1}{2}\langle \dot{q}, \dot{q} \rangle = c$ and is $2\pi\sqrt{2c}$. Comparisons of the ratio are presented into the tables 1-2 associated with the two costs: simplified vs mechanical energy. Note that the ratio $L/x_0$ is different from the concept of efficiency of the literature which takes into account the parametrization and the initial shape of the
stroke [12].

In table 2 we represent this ratio for simple loops with different amplitudes. This is in accordance with results of [2].

Finally, taking into account the constraints, the optimal numerical solution is the simple loop tangent to the triangle.

<figure><img src="image_11.png" /><figcaption>FIGURE 10. Creeping normal stroke for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost obtained by the Bocop software.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Tangential normal stroke for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost obtained by the Bocop software.</figcaption></figure>

<table><caption>TABLE 1. Ratio $L/x_0$ for abnormal solution and different normal strokes with the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost.</caption><thead><tr><th>Types of γ</th><th>x<sub>0</sub></th><th>L(γ)</th><th>L(γ)/x<sub>0</sub></th></tr></thead><tbody><tr><td>Abnormal</td><td>2.742·10<sup>-1</sup></td><td>10.73</td><td>39.13</td></tr><tr><td>Simple loop (Fig.10)</td><td>2.700·10<sup>-1</sup></td><td>7.689</td><td>28.48</td></tr><tr><td>Simple loop (Fig.12, left)</td><td>2.490·10<sup>-1</sup></td><td>6.314</td><td><b>25.36</b></td></tr><tr><td>Limaçon (Fig.14, left)</td><td>2.220·10<sup>-1</sup></td><td>6.602</td><td>29.74</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Normal stroke for the $\int_0^{2\pi}(u_1^2+u_2^2)dt$ cost (left) and the mechanical cost (right) where the constraints are satisfied (simple loop).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13. Second order sufficient condition checked on a stroke for the simplified cost $\int_0^{2\pi}(u_1^2+u_2^2)dt$ (left) and the mechanical cost (right). The smallest singular value associated with the rank condition (10) doesn't vanish on $[0, 2\pi]$, and there is no conjugate time.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 14. Normal stroke for the $\int_0^{2\pi}(u_1^2+u_2^2)dt$ cost (left) and the mechanical cost (right) where the constraints are satisfied (limaçon with inner loop).</figcaption></figure>

<table><caption>TABLE 2. Ratio $L/x_0$ for abnormal solution and different normal strokes with the mechanical cost.</caption><thead><tr><th>Types of γ</th><th>x<sub>0</sub></th><th>L(γ)</th><th>L(γ)/x<sub>0</sub></th></tr></thead><tbody><tr><td>Abnormal</td><td>2.742·10<sup>-1</sup></td><td>4.933</td><td>17.99</td></tr><tr><td>Simple loop (Fig.12, right)</td><td>2.600·10<sup>-1</sup></td><td>3.046</td><td><b>11.71</b></td></tr><tr><td>Limaçon (Fig.14, right)</td><td>2.500·10<sup>-1</sup></td><td>3.353</td><td>13.41</td></tr><tr><td>Small Amplitudes (Fig.15, right)</td><td>1.500·10<sup>-1</sup></td><td>1.860</td><td>12.40</td></tr><tr><td>Small Amplitudes (Fig.15, right)</td><td>0.500·10<sup>-1</sup></td><td>9.935·10<sup>-1</sup></td><td>19.87</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>FIGURE 15. Simple loops with different amplitudes for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost (left) and the mechanical cost (right).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 16. Illustrations of several simple loops with different amplitudes in the $(\theta_1, \theta_2)$ plane for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost (left) and the mechanical cost (right) (and corresponding displacements).</figcaption></figure>

4. **Conclusion.** The main contributions of this article are

* Complex strokes are simulated and their respective optimality are compared using the concept of conjugate points showing that the simple loops are the only candidates.

* The abnormal triangle is not optimal due to the existence of corners.

* The efficiency of the simple loops with respect to the amplitude is numerically investigated.

* The geometric optimal control relates the abnormal and normal strokes to the observation of [22].
Additional work is to complete our analysis with the Maximum Principle with state constraints and to use LMI techniques implemented in the *GloptiPoly* software to evaluate the global optimum [10].

Finally this article is a first step to analyze the more complex Purcell swimmer. It is not a symmetric swimmer but the displacements are invariant for the group of Euclidean motions. An important question is the existence of abnormal smooth strokes.

## REFERENCES

[1] Agrachev, A., Alaoui, El-H.C., Gauthier, J.P., Kupka, I.: Generic singularities of sub-Riemannian metrics on $\mathbb{R}^3$. C. R. Acad. Sci. Paris Sér. I Math. **322**, 4 377–384 (1996)

[2] Alouges, F., DeSimone, A., Giraldi, L., Zoppello, M.: Self-propulsion of slender microswimmers by curvature control: N-link swimmers. Int. J. Nonlinear Mech. **56** (2013)

[3] Alouges, F., DeSimone, A., Lefebvre, A.: Optimal strokes for low Reynolds number swimmers: an example. J. Nonlinear Sci. **18**, 277–302 (2008)

[4] Berger, M.: La taxonomie des courbes. Pour la science, 297 56–63 (2002)

[5] Becker, L.E., Koehler, S.A., Stone, H.A.: On self-propulsion of micro-machines at low Reynolds number: Purcell’s three-link swimmer. J. Fluid Mech. **490**, 15–35 (2003)

[6] Bettiol, P., Bonnard, B., Giraldi, L., Martinon, P., Rouot, J.: The three links Purcell swimmer and some geometric problems related to periodic optimal controls. Rad. Ser. Comp. App. **18**, Variational Methods, Ed. by M. Bergounioux et al. (2016)

[7] Bonnans, F., Giorgi, D., Maindrault, S., Martinon, P., Grélard, V.: Bocop - A collection of examples, Inria Research Report, Project-Team Commands. **8053** (2014)

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

[9] Bonnard, B., Chyba, M.: Méthodes géométriques et analytiques pour étudier l'application exponentielle, la sphère et le front d'onde en géométrie sous-riemannienne dans le cas Martinet. ESAIM Control Optim. Calc. Var. **4**, Springer-Verlag, Berlin 245–334 (1999)

[10] Bonnard, B., Claeys, M., Cots, O., Martinon, P.: Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance. Acta Appl. Math. **135**, 5–45 (2015)

[11] Bonnard, B., Faubourg, L., Trélat, E.: Mécanique céleste et contrôle des véhicules spatiaux. Mathématiques & Applications **51**, Springer-Verlag Berlin (2006)

[12] Chambrion, T., Giraldi, L., Munnier, A.: Optimal strokes for driftless swimmers: A general geometric approach. Submitted (2014)

[13] Cots, O.: Contrôle optimal géométrique: méthodes homotopiques et applications. PhD thesis, Université de Bourgogne (2012)

[14] Giraldi, J., Martinon, P., Zoppello, M.: Optimal design of the three-link Purcell swimmer. Phys. Rev. E **91** (2015)

[15] Gray, J., Hancock, G. J.: The propulsion of sea-urchin spermatozoa. J. Exp. Biol. **32**, 802–814 (1955)

[16] Hakavuori, E., Le Donne, E.: Non-minimality of corners in sub-Riemannian geometry. Preprint (2015)

[17] Lauga, E., Powers, T.R.: The hydrodynamics of swimming microorganisms. Rep. Progr. Phys. **72**, 9 (2009)

[18] Montgomery, R.: A tour of subriemannian geometries, their geodesics and applications. American Mathematical Society, Providence, RI. **91** (2002)

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

[20] Lenz, P. H., Takagi, D., Hartline, D. K.: Choreographed swimming of copepod nauplii. Submitted (2015)

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

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

[23] Tam, D., Hosoi, A.E.: Optimal Stroke Patterns for Purcell Three-Link Swimmer. Phys. Rev. Lett. **98** (2007)

[24] Taylor, G.I.: Analysis of the swimming of microscopic organisms. Proc. Roy. Soc. London. Ser. A. **209**, 447–461 (1951)
[25] Zhitomirskiǐ, M.: Typical singularities of differential 1-forms and Pfaffian equations. American Mathematical Society, Providence, RI. **113** (1992)

E-mail address: bernard.bonnard@u-bourgogne.fr

E-mail address: mchyba@math.hawaii.edu

E-mail address: jeremy.rouot@inria.fr

E-mail address: dtakagi@hawaii.edu

E-mail address: zourong@sigmath.es.osaka-u.ac.jp