# Optimal Control Theory and the Swimming Mechanism of the Copepod Zooplankton

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

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

To cite this version:

Piernicola Bettiol, Bernard Bonnard, Alice Nolot, Jérémy Rouot. Optimal Control Theory and the Swimming Mechanism of the Copepod Zooplankton. 2016. hal-01387443

HAL Id: hal-01387443

https://inria.hal.science/hal-01387443

Submitted on 25 Oct 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# OPTIMAL CONTROL THEORY AND THE SWIMMING MECHANISM OF THE COPEPOD ZOOPLANKTON

P. BETTIOL, B. BONNARD, A. NOLOT, J. ROUOT

ABSTRACT. In this article, the model of swimming at low Reynolds number introduced by D. Takagi (2015) to analyze the displacement of an abundant variety of zooplankton is used as a testbed to analyze the motion of symmetric microswimmers in the framework of optimal control theory assuming that the motion occurs minimizing the energy dissipated by the fluid drag forces in relation with the concept of efficiency of a stroke. The maximum principle is used to compute periodic controls candidates as minimizing controls and is a decisive tool combined with appropriate numerical simulations using indirect optimal control schemes to determine the most efficient stroke compared with standard computations using Stokes theorem and curvature control. Also the concept of graded approximations in SR-geometry is used to evaluate strokes with small amplitudes providing a fixed displacement and minimizing the dissipated energy.

## CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. A geometric analysis of the problem in the frame of SR-geometry</td><td>4</td></tr><tr><td>2.1. Geodesic computation</td><td>4</td></tr><tr><td>2.2. Computation in the copepod case</td><td>5</td></tr><tr><td>2.3. SR-classification in dimension 3 and strokes with small amplitude for the copepod swimmer</td><td>6</td></tr><tr><td>2.4. Computations of the normal form of order 1</td><td>10</td></tr><tr><td>3. A powerful approach using optimal control theory and numerical simulations</td><td>13</td></tr><tr><td>3.1. The maximum principle with periodic controls</td><td>14</td></tr><tr><td>3.2. Second-order necessary optimality condition</td><td>15</td></tr><tr><td>3.3. Applications: numerical simulations</td><td>16</td></tr><tr><td>4. Geometric efficiency and optimality of the abnormal stroke</td><td>22</td></tr><tr><td>4.1. Some comments on the efficiency term</td><td>22</td></tr><tr><td>4.2. Some comments on normality for the Copepod model</td><td>24</td></tr><tr><td>References</td><td>28</td></tr></table>
# 1. INTRODUCTION

This article is entirely devoted to the analysis combining optimal control theory and sub-Riemannian (SR-) geometry of the swimming process of a variety of zooplankton observed by Takagi (2015) and modeled in the framework of swimming at low Reynolds number. See Fig.1 for the picture of the copepod (left) and the 2-link symmetric micro-robot swimmer to mimic the animal mechanism.

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

In micro-robot modeling, to produce the displacement along the line $Ox_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 constraint $0 \le \theta_1 \le \theta_2 \le \pi$.

Using the swimming model at low Reynolds number, we relate the speed of the displacement to the speed of the shape variable $\theta$ by the equation

$$ (1.1) \qquad \dot{x}_0 = \frac{\sum_{i=1}^{2} \dot{\theta}_i \sin(\theta_i)}{\sum_{i=1}^{2} (1 + \sin^2(\theta_i))}. $$

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} = \sum_{i=1}^{2} u_i F_i(q) $$

with $q = (x_0, \theta)$, $\theta = (\theta_1, \theta_2)$. Moreover we have state constraints given by a triangle $T$ in the shape variables: $\theta_i \in [0, \pi]$, $i = 1, 2$, and $\theta_1 \le \theta_2$. $u_1, u_2$ are periodic controls producing strokes, which are closed curves in the $\theta$-plane, and the reference problem can be analyzed in the framework of optimal control theory introducing a cost function. A choice of particular interest for the cost to minimize, in particular in relation with the concept of efficiency defined by Lightwill (1960), is the mechanical energy dissipated
by the drag forces:

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

where 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}$$

Using (1.1) this amounts to minimize the quadratic form

$$E = \int_{0}^{T} (a_{1}(\theta)u_{1}^{2} + 2b_{1}(\theta)u_{1}u_{2} + c_{1}(\theta)u_{2}^{2})dt$$

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)}.$$

One denotes g the associated Riemannian metric in the θ-space, the optimal control problem is rewritten as a sub-Riemannian problem

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

with appropriate boundary conditions associated with periodic control

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

the triangle inequality constraints $0 \le \theta_1 \le \theta_2 \le \pi$. $x_0(T) - x_0(0)$ represents the displacement of a stroke and one can set $x_0(0) = 0$. According to Maupertuis principle, this is equivalent to minimize the length

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

and using the energy minimization point of view, the period T of a stroke can be fixed to $T = 2\pi$. We emphasize that the problem is equivalent to a time minimal control by fixing the energy level g = 1 of the strokes.

From this point of view, the optimal control problem consists in computing the points of the sub-Riemannian sphere $S_{q_0}(r)$ formed by extremities of the minimizers starting from $q_0$ and with fixed length, and requiring that the optimal control is periodic. This is equivalent to fix the displacement $x_0(2\pi) - x_0(0)$ and to compute strokes minimizing the length.

The concept of geometric efficiency has a clear meaning in the SR-geometry context. Assuming $x_0(0) = 0$, the geometric efficiency is the ratio

$$(1.2) \qquad \mathcal{E} = x_0(T)/L$$
where $L$ is the length of the stroke producing the displacement $x_0(T)$, which is proportional, for a fixed $g$, to $\mathcal{E}' = x_0^2(T)/E$, which is introduced by Takagi (2015).

A further concept of efficiency can be similarly used, see for instance Chambrion (2014). It takes the ratio between the energy used to move the swimmer at constant speed $\bar{v}$ to produce the displacement $x_0(T)$ and the mechanical energy where the shape variables are hold on at $\theta(0)$, that is

$$ (1.3) \qquad Eff = \frac{|\bar{v}|^2 M_{11}(\theta(0))}{1/TE}, \quad M = (M_{ij}), \quad \bar{v} = \frac{1}{T} \int_0^T \dot{x}_0 dt. $$

We have

$$ Eff \sim \frac{x_0^2(T)}{E} M_{11}(\theta(0)) $$

If the concept of geometric efficiency is related to SR-geometry, the problem of maximizing the efficiency can be equally treated by techniques of optimal control. This will be the main achievement of this article.

The paper will be organized in two sections. The first section represents a geometric analysis, in relation with SR-geometry, and is devoted to the problem of computing optimal strokes with small amplitudes. The second section is a direct application of the maximum principle in the frame of periodic optimal control complemented by second order optimality conditions and numerical simulations to compute strokes with the problem of maximizing the different of efficiencies. A final section is devoted to the analysis of the optimality of the triangle abnormal stroke.

## 2. A GEOMETRIC ANALYSIS OF THE PROBLEM IN THE FRAME OF SR-GEOMETRY

### 2.1. Geodesic computation.
Consider the energy minimization problem

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

where the set of admissible controls $\mathcal{U}$ is the set of bounded measurable mapping valued in $\mathbb{R}^2$. 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)$. According to the maximum principle (see Vinter (2000)), minimizers are found among extremals curves, which are solutions 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.

*Normal case.* Assume $p_0 \neq 0$ and it can be normalized to $p_0 = -1/2$. Corresponding controls are given by

$$ u_1 = \frac{cH_1 - bH_2}{ab - b^2}, \quad u_2 = \frac{aH_2 - bH_1}{ab - b^2}. $$
where $H_1, H_2$ are the Hamiltonian lifts of the vector fields $F_1, F_2$ and plugging such controls in $H$ yields the normal Hamiltonian

$$ (2.1) \qquad H_n = \frac{1}{2(ac - b^2)} (aH_2^2 - 2bH_1H_2 + cH_1^2). $$

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

$$ (2.2) \qquad H_1 = H_2 = 0. $$

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) $$

and if $H_F = p \cdot F(q)$, $H_G = p \cdot G(q)$, the Poisson bracket is: $\{H_F, H_G\}(z) = dH_F(H_G) = p \cdot [F, G](q)$. Differentiating twice (2.2) with respect to time, abnormal controls are given by

$$ (2.3) \qquad \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} $$

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

**Definition 2.1.** The exponential mapping is, for fixed $q(0)$ the map: $\exp_{q_0}: (t, p(0)) \mapsto \Pi(\exp t\overrightarrow{H_n}(z(0)))$ where $\Pi$ is the standard projection: $\Pi: (q, p) \mapsto 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$ and $t_{1c}$ denotes the first conjugate time and $q(t_{1c})$ is called the first conjugate point along the reference geodesic $t \mapsto q(t)$.

**Definition 2.2.** 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 in the copepod case.** One has:

$$ F_i = \frac{\sin(\theta_i)}{\Delta} \frac{\partial}{\partial x_0} + u_i \frac{\partial}{\partial \theta_i}, \quad i = 1, 2. $$
with $\Delta = \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) = \frac{2 \sin(\theta_1) \sin(\theta_2) (\cos(\theta_1) - \cos(\theta_2))}{\Delta^2}.$$

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 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$ and $[F_1, F_2]$ are coplanar but at least 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 $T: 0 \le \theta_1 \le \theta_2 \le \pi$ are contact points.

2) The sides of the triangle (vertices excluded) are Martinet point and the triangle is a (piecewise smooth) abnormal stroke.

Geometric comments. Hence the observed stroke by Takagi (2015) corresponds to the policy:

$$\theta_2: 0 \to \pi, \theta_1: 0 \to \pi, \theta_i: \pi \to 0, i = 1, 2 \text{ with } \theta_1 = \theta_2.$$

where the copepod swimmer follows the triangle boundary of the physical domain corresponding to the unique abnormal stroke.

Moreover it has a nice geometric interpretation using Stokes' theorem and curvature control methods.

**Lemma 2.6.** One has:

$$1) \oint \sum_{i=1}^{2} \frac{\sin(\theta_i)}{\Delta} d\theta_i = \int \left[ \frac{\partial}{\partial \theta_2} \left( \frac{\sin(\theta_1)}{\Delta} \right) - \frac{\partial}{\partial \theta_1} \left( \frac{\sin(\theta_2)}{\Delta} \right) \right] d\theta_1 \wedge \theta_2 = \int d\omega.$$

2) The points where $d\omega = 0$ are precisely the abnormal triangle, and $d\omega < 0$ in the interior domain and $d\omega > 0$ in the exterior.

## 2.3. SR-classification in dimension 3 and strokes with small amplitude for the copepod swimmer.
2.3.1. The contact case. The crucial results applicable to our model come from Alaoui et al. (1996). Consider a standard SR-problem $(D, g)$, where $\Delta$ is a distribution and $g$ is a SR-metric, near a point $q_0 \in \mathbb{R}^3$ identified with 0, one has:

* Heisenberg-Brockett nilpotent model. The nilpotent (of order $-1$) is the so-called Heisenberg-Brockett model where $(D, g)$ is defined by the orthonormal frame: $D = \operatorname{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} $$

with $q = (x, y, z)$ and the graduation 1 for $x, y$ and 2 for $z$, forming a set of privileged coordinates.

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.

We introduce $\tilde{Q}(w) = (1 + Q(w))$ and 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. We use 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_1 = p \cdot \frac{\partial}{\partial z} $$

and the normal Hamiltonian is

$$ H_n = \frac{1}{2}(H_1^2 + H_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 constant.

$$ \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, one can introduce the parametrization

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

Denoting by $'$ the derivative with respect to $s$, we get

$$ H'_{1} = H_{2}, \quad H'_{2} = -H_{1}. $$

Hence we deduce

**Lemma 2.7.** In the s-parameter, the normal controls are solutions of the linear pendulum equation $H''_{1} + H_{1} = 0$ and are trigonometric functions.

A more precise analysis requires higher order Lie brackets computations. Note that, in particular, the relations

$$ [[F, G], F] = \frac{\partial f}{\partial x} \frac{\partial}{\partial z}, \quad [[F, G], G] = \frac{\partial f}{\partial y} \frac{\partial}{\partial z}. $$

The remaining equations to be integrated using (2.4) are

$$ (2.5) \qquad 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)}. $$

The solution can be estimated using the following expansion, associated with the gradation, setting

$$ \begin{gathered} x = \varepsilon X, \quad y = \varepsilon Y, \quad z = \varepsilon^2 Z, \\ p_x = p_X, \quad p_y = p_Y, \quad p_z = p_Z / \varepsilon. \end{gathered} $$

so that the Darboux form is homogeneous.

Denoting $Q = (X, Y, Z)$, $P = (P_X, P_Y, P_Z)$ the solution can be obtained in the expansion

$$ \begin{align} (2.6) \qquad & X(t) = X_0(t) + \varepsilon X_1(t) + o(\varepsilon), \nonumber \\ & Y(t) = Y_0(t) + \varepsilon Y_1(t) + o(\varepsilon), \nonumber \\ & Z(t) = Z_0(t) + \varepsilon Z_1(t) + o(\varepsilon). \end{align} $$

Clearly by identification one gets that $t \mapsto (X_0(t), Y_0(t), Z_0(t))$ is the solution obtained by the (nilpotent) Heisenberg-Brockett model and similarly for the higher order expansions.

Heisenberg-Brockett solution. We recall the standard computation Brockett (1982) for the solutions starting from $q(0) = 0$. We have

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

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

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

where $A, \lambda, \varphi$ are parameters.

In particular, in relation with the swimmer problem, in the shape variables identified with $(x, y)$ we get a one-parameter family of circles on each
energy level, each of them deduced by a proper rotation $R_α$ along the z-axis (standard symmetry of revolution of the geodesic equations), see Fig.2.

<figure><img src="image_4.png" /><figcaption>FIGURE 2. One parameter family of circles which are the geodesics of the Heisenberg-Brockett problem.</figcaption></figure>

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.

Clearly due to the symmetries, the model is not stable and higher order terms have to be taken into account to analyze strokes.

Application to the copepod analysis. We choose a point $\theta(0) = (\theta_1(0), \theta_2(0))$ in the interior of the triangle identified with 0 using the translation

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

and using Taylor expansions one can approximate the displacement using

$$\dot{z} = \frac{u_1 \sin(\theta_1) + u_2 \sin(\theta_2)}{2 + \sin^2(\theta_1) + \sin^2(\theta_2)} = c_1 u_1 + c_2 u_2 + o(1)$$

where $c_1, c_2$ are constant. Using $\dot{x} = u_1$ and $\dot{y} = u_2$ we set

$$(2.10) \qquad Z = z - c_1x - c_2y$$

to get an equation of the form

$$\dot{Z} = u_1 d_1(w) + u_2 d_2(w),$$

with $d_i = o(1)$, and $w = (x, y)$. We have

**Lemma 2.8.** $q = (x, y, Z)$ will form near $(0, \theta_1(0), \theta_2(0))$ a graded system of coordinates with respective weight $(1, 1, 2)$.

Geometric comments. The relation (2.10) relates locally the physical displacement variable to the corresponding displacement variable in the Heisenberg-Brockett model. Further computations are necessary to compute the normal form of order 1 and analyze the geodesics using the expansion procedure in (2.6). Next they are given in a special case.
2.4. **Computations of the normal form of order 1.** Setting $x = \theta_1 - \pi/4$, $y = \theta_2 - 3\pi/4$ and $z = x_0$, the Taylor expansion of order 4 of the vector fields $F_1, F_2$ near $(x, y) \sim (0, 0)$ is given by

$$F_1 = \begin{bmatrix} 1 \\ 0 \\ \frac{\sqrt{2}(16x^3 - 15x^2y - 10y^3 - 39x^2 + 6yx + 6y^2 + 36x + 18y + 54)}{324} \end{bmatrix}$$

$$F_2 = \begin{bmatrix} 0 \\ 1 \\ \frac{\sqrt{2}(10x^3 + 15y^2x - 16y^3 + 6x^2 + 6yx - 39y^2 - 18x - 36y + 54)}{324} \end{bmatrix}$$

**Notation 2.9.** Before the change of variables, the old coordinates are $(x, y, z)$ and the new ones are $(X, Y, Z)$. After each change of variables, the resulting vector fields are denoted by $F_1, F_2$ and they are written in the coordinates $(x, y, z)$.

Change of variables and feedback computations.

(1) we remove the constant terms in $F_{13}, F_{23}$; the $x, x^2, x^3$'s terms in $F_{13}$ and the $y, y^2, y^3$'s terms in $F_{23}$.

After the change of variables $[X, Y, Z] = \varphi_1(x, y, z) = [x, y, 162\sqrt{2}z - 4x^4 + 4y^4 + 13x^3 + 13y^3 - 18x^2 + 18y^2 - 54x - 54y]$, the resulting vector fields are

$$F_1 = [1, 0, -15x^2y - 10y^3 + 6yx + 6y^2 + 18y]$$

$$F_2 = [0, 1, 10x^3 + 15y^2x + 6x^2 + 6yx - 18x]$$

(2) we normalize the linear terms of $F_{13}$ et $F_{23}$ in the Heisenberg form and we remove the $y^2, y^3$'s terms in $F_{13}$ and the $x^2$'s term in $F_{23}$.

After the change of variables $[X, Y, Z] = \varphi_2(x, y, z) = [x, y, z/36 - \frac{5yx^3}{18} + \frac{5xy^3}{18} - 1/18x^3 - 1/6x^2y - 1/6xy^2 - 1/18y^3]$, the resulting vector fields are

$$F_1 = [1, 0, -5/4x^2y - 1/6x^2 - 1/6xy + y/2]$$

$$F_2 = [0, 1, 5/4xy^2 - 1/6xy - 1/6y^2 - x/2]$$

(3) we remove the terms of order 0 in $F_{13}$ et $F_{23}$. We consider the change of variables $\varphi_4$ given by

$$\varphi_4^{-1}(P) = \begin{pmatrix} X - \alpha_1 Z + \alpha_1 XY/2 + \alpha_2 Y^2/2 \\ Y - \alpha_2 Z - \alpha_1 X^2/2 - \alpha_2 XY/2 \\ Z + \alpha_1 YZ/2 - \alpha_2 XZ/2 + \text{high order terms} \end{pmatrix}.$$
$$ \text{More precisely, } \varphi_4^{-1}(X, Y, Z) = \begin{pmatrix} x = X + 1/12 XY + 1/12 Y^2 - Z/6 \\ y = Y - 1/12 X^2 - 1/12 XY - Z/6, \\ z = Z - 1/24 X^3 - 1/24 X^2Y \\ -1/24 XY^2 - 1/24 Y^3 - 1/12 XZ \\ + \frac{5X^2Z}{144} + 1/12 YZ + \frac{XYZ}{72} \\ + \frac{5Y^2Z}{144} - \frac{XZ^2}{108} - \frac{YZ^2}{108} \end{pmatrix} $$

Remark 2.10. $\varphi_4^{-1} \neq (X, Y, Z) + \theta(X, Y, Z) + \dots$ so $\varphi_4(x, y, z) \neq (x, y, z) - \theta(x, y, z) + \dots$

The inverse transformation is

$$ \varphi_4(x, y, z) = \begin{pmatrix} X = x - 1/12 xy - \frac{x^2 y}{72} - 1/12 y^2 + \frac{y^3}{72} + z/6 - 1/18 yz \\ Y = y + 1/12 x^2 + \frac{x^3}{72} + 1/12 xy - \frac{xy^2}{72} + z/6 + 1/18 xz \\ Z = z + 1/24 x^3 + 1/24 x^2 y + 1/24 xy^2 + 1/24 y^3 + 1/12 xz - \frac{x^2 z}{144} - 1/12 yz - \frac{xyz}{72} - \frac{y^2 z}{144} \end{pmatrix} $$

Applying the change of variables $\varphi_4$ and the feedback transformation $u_1 \leftarrow u_1 + u_2(\alpha_1 x + \alpha_2 y)$, $u_2 \leftarrow u_2 - u_1(\alpha_1 x + \alpha_2 y)$, the resulting vector fields are

$$ F_1 = [1, 0, y/2 + \frac{y^3}{288} - \frac{x^3}{72} - \frac{359 x^2 y}{288} + 1/48 x y^2], $$

$$ F_2 = [0, 1, -x/2 + \frac{y^3}{72} - \frac{x^3}{288} + \frac{359 x y^2}{288} - 1/48 x^2 y]. $$

(4) we normalize the terms of order 1 in $F_{13}$ and in $F_{23}$ to introduce the quadratic form $Q$. After the change of variables $[X, Y, Z] = \varphi_5(x, y, z) = [x, y, z + \frac{x^4}{288} + \frac{5yx^3}{16} - \frac{5xy^3}{16} - \frac{y^4}{288}]$, the resulting vector fields are

$$ F_1 = [1, 0, -\frac{89 x^2 y}{288} - \frac{89 y^3}{288} + y/2 + 1/48 x y^2] $$

$$ F_2 = [0, 1, \frac{89 x^3}{288} + \frac{89 x y^2}{288} - x/2 - 1/48 x^2 y] $$

in other words,

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

$$ F_2 = F_2^{(-1)} - \frac{x}{2} Q(x, y) \frac{\partial}{\partial z} $$

where $Q(x, y) = -89/144(x^2 + y^2) + 1/24xy$ and $F_1^{(-1)}, F_2^{(-1)}$ are the vector fields of order $-1$.

Remark 2.11. The change of coordinates doesn't let invariant the angular variables $\theta_1, \theta_2$. To relate the nilpotent system to the real one, we shall use the transformation $T(x, y, z) = \varphi_5 \circ \varphi_4 \circ \varphi_3 \circ \varphi_2 \circ \varphi_1(x, y, z)$.
2.4.1. The Martinet case. The analysis at a Martinet point $q_0$ belonging to the vertices of the triangle is more intricate and can be analyzed using the results of Bonnard et al. (2003).

Generic model at a Martinet point $q_0$ identified with 0. There exists local coordinates $q = (x, y, z)$ such that the SR-geometry is given by $(D, g)$ with

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

* The metric $q$ is of the form

$a(q)dx^2 + c(q)dy^2$ and we have

* Nilpotent model (flat Martinet case): $a = c = 1$.

* Generic model:

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

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

where $\alpha, \beta, \gamma$ are parameters.

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 := 1/2(H_1^2 + H_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{array}{l} \dot{x} = \frac{\cos \chi}{\sqrt{a}}, \quad \dot{y} = \frac{\sin \chi}{\sqrt{b}}, \quad \dot{z} = \frac{y^2 \cos \chi}{2\sqrt{a}}, \\ (2.11) \qquad \dot{\chi} = \frac{1}{\sqrt{ac}}(y\lambda - \alpha \cos \chi - \beta \sin \chi) \end{array}$$

We introduce the new parameterization

$$(2.12) \qquad \sqrt{ac} \frac{d}{dt} = \frac{d}{ds}$$

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)$$

and the second order differential equation

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

As a consequence, we obtain the following result.

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

• The flat case is $\alpha = \beta = \gamma = 0$ and corresponds to the standard pendulum.

• In the generic case we have two subcases
  – $\beta = 0$ and (2.13) is integrable using elliptic functions.
  – $\beta \neq 0$, due to dissipation, we are in the non-integrable case.

Application to the copepod. In the integrable case: $\beta = 0$ models of periodic strokes with elliptic functions with modulus $k$ are

• $k = 0$: circles

• $k \approx 0.65$: eight shape (Bernoulli lemniscates). Note that $\beta = 0$ is not a stable case, moreover the triangle constraint $0 \le \theta_1 \le \theta_2 \le \pi$ is not taken into account.

Conclusion about SR-model. Candidates as strokes using the SR-models are represented on Fig.3 with respect to the triangle of constraints: $0 \le \theta_1 \le \theta_2 \le \pi$.

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Simple candidates as strokes for the SR models.</figcaption></figure>

They are:

• simple curves: circles ①,

• limaçons: perturbation of a simple curve by period doubling ②,

• eight shape: Bernoulli lemniscates ③.

Note that the orientation is imposed by Lemma 2.6.

### 3. A POWERFUL APPROACH USING OPTIMAL CONTROL THEORY AND NUMERICAL SIMULATIONS

In this section, the problem is analyzed using the maximum principle applied for optimal control with periodic controls (see Vinter (2000)) and
complemented by necessary second order optimality conditions corresponding to the concept of conjugate point.

**3.1. The maximum principle with periodic controls.** The crucial point is the existence of a maximum principle suitable to analyze the problem of maximizing different concept of efficiencies. The system and the energy are written in the extended state space with $\bar{q} = (q, q^0)$ and the corresponding dynamic

$$ 
\begin{aligned}
\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{aligned}
 $$

and the problem is to minimize a cost of the form

$$ \min_{u(.)} h(\bar{q}(0), \bar{q}(2\pi)) $$

where the end-point conditions are 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, p_0)$ 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_i \cdot F_i(q)$.

From Vinter (2000), an optimal control pair $(q, u)$ is satisfying 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.

$$ (3.1) \quad (\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)) $$

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 mechanical energy

$$ E = \int_{0}^{2\pi} (u_{1}^{2} + u_{2}^{2})dt. $$
In this case we deduce the periodicity condition on $p_{\theta} = (p_{\theta_1}, p_{\theta_2})$ dual to $\theta$:

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

to produce a smooth stroke in the normal case $p_0 \neq 0$.

Moreover $(p_x, p_0)$ at the final point have to be collinear to the gradient of the set $g(x_0, x^0) = c$, see Fig.4.

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Accessibility set and transversality conditions with a cost function $h$.</figcaption></figure>

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

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

where $m$ is a chosen smooth function.

In this case, (3.2) becomes:

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

hence producing a *jump* of the adjoint vector.

## 3.2. Second-order necessary optimality condition.
It is the standard necessary optimality condition related to existence of conjugate point (see Bonnard et al. (2003)).

**Proposition 3.1.** 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[$.

It can be checked numerically using the HamPath code.
**3.3. Applications: numerical simulations.** We present a sequence of simulations using the HamPath software. These are based on our computations for the $\int_0^{2\pi}(u_1^2 + u_2^2)dt$ cost.

* Fig.5: Different kind of normal strokes: simple loop, limaçon and eight and computation of conjugate points. Only the simple loops are candidates to be optimal strokes.

* There is a one-parameter family of simple strokes, each of them associated with a different energy. The corresponding efficiency is represented in Table 1 and compared with the efficiency of the abnormal stroke, producing the maximum efficiency value. We have numerically checked that it corresponds to the stroke with the transversality condition provided by the maximum principle (3.1), see Fig.10.

* As in Chambrion (2014), we consider a cost depending upon $\theta(0)$, namely

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

where $m(\theta_1(0)) = 2 - \cos^2(\theta_1(0))$. In Fig.11 is illustrated the corresponding optimal solution satisfying the transversality conditions (3.1), it is a non smooth stroke and it can be compared to the previous smooth solution of Fig.10.

*Remark 3.2.* From the simulations, we have the following property:

There is a one-parameter family of simple loops whose projection on the $\theta$-plane is symmetric with respect to the straight line ($D$): $\theta_2 = -\theta_1 + \pi$, see Fig.12.

<table><caption>TABLE 1. Value of the geometric efficiency for abnormal solution and different normal strokes for the Copepod swimmer.</caption><thead><tr><th>Types of &gamma;</th><th>x<sub>0</sub>(T)</th><th>L(&gamma;)</th><th>x<sub>0</sub>(T)/L(&gamma;)</th></tr></thead><tbody><tr><td rowspan="4">Simple loops</td><td>1.400 &times; 10<sup>-1</sup></td><td>3.785</td><td>3.698 &times; 10<sup>-2</sup></td></tr><tr><td>1.700 &times; 10<sup>-1</sup></td><td>4.340</td><td>3.917 &times; 10<sup>-2</sup></td></tr><tr><td>2.000 &times; 10<sup>-1</sup></td><td>4.946</td><td>4.043 &times; 10<sup>-2</sup></td></tr><tr><td>2.100 &times; 10<sup>-1</sup></td><td>5.109</td><td>4.110 &times; 10<sup>-2</sup></td></tr><tr><td rowspan="6"><b>Normal stroke</b><br>Fig.10</td><td><b>2.169 &times; 10<sup>-1</sup></b></td><td><b>5.180</b></td><td><b>4.187 &times; 10<sup>-2</sup></b></td></tr><tr><td>2.200 &times; 10<sup>-1</sup></td><td>5.354</td><td>4.109 &times; 10<sup>-2</sup></td></tr><tr><td>2.300 &times; 10<sup>-1</sup></td><td>5.624</td><td>4.089 &times; 10<sup>-2</sup></td></tr><tr><td>2.500 &times; 10<sup>-1</sup></td><td>6.305</td><td>3.965 &times; 10<sup>-2</sup></td></tr><tr><td>2.740 &times; 10<sup>-1</sup></td><td>9.046</td><td>3.028 &times; 10<sup>-2</sup></td></tr><tr><td>Abnormal</td><td>2.742 &times; 10<sup>-1</sup></td><td>10.73</td><td>2.555 &times; 10<sup>-2</sup></td></tr><tr><td>Limaçon</td><td>2.000 &times; 10<sup>-1</sup></td><td>6.147</td><td>3.253 &times; 10<sup>-2</sup></td></tr><tr><td>Eight</td><td>2.000 &times; 10<sup>-1</sup></td><td>6.954</td><td>3.307 &times; 10<sup>-2</sup></td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>FIGURE 5. One parameter family of simple loops, limaçons and Bernoulli lemniscates normal strokes</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Efficiency curve for simple loop normal strokes and efficiency of the abnormal stroke.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Normal stroke of the Copepod swimmer with limaçon shape. The first conjugate point is computed (indicated by the cross)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Normal stroke of the Copepod swimmer with eight shape. The first conjugate point is computed (indicated by the cross)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Normal stroke of the Copepod swimmer with simple loop shape. There is no conjugate point on $[0, 2\pi[$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 10. Normal stroke of the Copepod swimmer for the geometric efficiency, obtained by the transversality conditions of the maximum principle (3.1).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIGURE 11. Normal stroke of the Copepod swimmer for the geometric efficiency depending upon the initial angle $\theta(0)$, obtained by the transversality conditions of the maximum principle (3.1)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 12. One parameter family of simple loops symmetric with respect to the straight line $\mathcal{D}: \theta_2 = -\theta_1 + \pi$ with converges to a point when the displacement tends to 0.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIGURE 13. Normal stroke tangent to the triangle.</figcaption></figure>
4. GEOMETRIC EFFICIENCY AND OPTIMALITY OF THE ABNORMAL STROKE

4.1. **Some comments on the efficiency term.** In this section we shall consider the following optimal control problem for the Copepod model:

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

in which $x_T > 0$ is a given number which represents the Copepod desired displacement, $T = 2\pi$ in the time period, $U = \mathbf{R}^2$, $f: \mathbf{R}^3 \times \mathbf{R}^2 \to \mathbf{R}^3$ is the function

$$
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).
$$

Here, the state variable is $q = (x, \theta = (\theta_1, \theta_2))$ in which $x$ stands for the swimmer displacement, and $\theta = (\theta_1, \theta_2)$ are the link amplitudes.

We shall establish a relationship between minimizers of the optimal control problem (4.1) and minimizers of an optimal control problem in which we minimize the geometric efficiency $\mathcal{E}$ (eventually modified with a suitable penalization term), studying the case in which we minimize a simplified cost which represent the energy of the system (this is the sub-Riemannian problem): $L(q, u) = u_1^2 + u_2^2$. More precisely, suppose that we are given a minimizer $(q^*, u^*)$ for (4.1). Write $E^* := \int_0^T L(q^*(t), u^*(t))dt$ the corresponding energy and consider the following optimal control problem:

$$
(4.2) \quad \left\{
\begin{array}{ll}
\displaystyle \text{Minimize} & J_{\mathcal{E}}((q_0(\cdot), q(\cdot)), u(\cdot)) := \left[-\frac{q_1(T)}{q_0(T)} + K|q_0(T) - E^*|\right] \\
\text{over arcs} & (q_0, q(\cdot)) \in W^{1,1}([0,T]; \mathbb{R}^4) \text{ s.t.} \\
\displaystyle \dot{q}(t) = f(q(t), u(t)), & \text{for a.e. } t \in [0, T] \\
\displaystyle \dot{q}_0(t) = L(q(t), u(t)), & \text{for a.e. } t \in [0, T] \\
u(t) \in U & \text{for a.e. } t \in [0, T] \\
q_0(0) = 0, & q_1(0) = 0 \\
q_2(0) = q_2(T), & q_3(0) = q_3(T),
\end{array}
\right.
$$

**Proposition 4.1.** Let $(q^*(\cdot), u^*(\cdot))$ be a local minimizer for (4.1), then we can find a positive constant $K$ such that $((q_0^*(\cdot), q^*(\cdot)), u^*(\cdot))$ is a local minimizer for (4.2).

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

$$ (4.3) \qquad \int_0^T L(q^*(t), u^*(t))dt \le \int_0^T L(q(t), u(t))dt, $$

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

$$ (4.4) \qquad \|q^*(.) - q(.)\|_{L^\infty} \le \varepsilon_0, \qquad \|u^*(.) - u(.)\|_{L^\infty} \le \varepsilon_0. $$

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

$$ (4.5) \qquad J_{\varepsilon}((q_0^*(.), q^*(.)), u^*(.)) \le J_{\varepsilon}((q_0(.), q(.)), u(.)), $$

for all trajectory-control pair $((q_0(.), q(.), u(.))$ satisfying all the conditions of the control system in (4.2), and with

$$ (4.6) \qquad \|(q_0^*(.), q^*(.)) - (q_0(.), q(.))\|_{L^\infty} \le \varepsilon_1, \qquad \|u^*(.) - u(.)\|_{L^\infty} \le \varepsilon_1. $$

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

$$ (4.7) \qquad J_{\varepsilon}((q_0^\varepsilon(.), q^\varepsilon(.), u^\varepsilon(.)) < J_{\varepsilon}((q_0^*(.), q^*(.), u^*(.)) = -\frac{q_1^*(T)}{q_0^*(T)} $$

and

$$ (4.8) \qquad \|(q_0^*(.), q^*(.)) - (q_0^\varepsilon(.), q^\varepsilon(.))\|_{L^\infty} \le \varepsilon, \qquad \|u^*(.) - u^\varepsilon(.)\|_{L^\infty} \le \varepsilon. $$

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

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

and

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

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

$$ J_{\varepsilon}((\tilde{q}_0(.), \tilde{q}(.)) = -\frac{\tilde{q}_1(T)}{\tilde{q}_0(T)}. $$

Observe that taking $\varepsilon_1 \in (0, \varepsilon_0)$ small enough, then for all $\varepsilon \in (0, \varepsilon_1)$, we obtain that

$$ (4.9) \qquad \|q^*(.) - \tilde{q}(.)\|_{L^\infty} \le \varepsilon_0, \qquad \|u^*(.) - \tilde{u}(.)\|_{L^\infty} \le \varepsilon_0. $$

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

$$ (4.10) \qquad \left| \frac{\tilde{q}_1(T)}{\tilde{q}_0(T)} - \frac{q_1^\varepsilon(T)}{q_0^\varepsilon(T)} \right| \le K|q_0^\varepsilon(T) - \tilde{q}_0(T)| = K|q_0^\varepsilon(T) - E^*|. $$
From (4.7) and (4.10) we would eventually deduce that

$$
\begin{align*}
J_{\mathcal{E}}((\tilde{q}_0(\cdot), \tilde{q}(\cdot)) &= -\frac{\tilde{q}_1(T)}{\tilde{q}_0(T)} \\
&\leq -\frac{q_1^{\varepsilon}(T)}{q_0^{\varepsilon}(T)} + K|q_0^{\varepsilon}(T) - E^*| = J_{\mathcal{E}}((q_0^{\varepsilon}(\cdot), q^{\varepsilon}(\cdot)), u^{\varepsilon}(\cdot)) \\
&< J_{\mathcal{E}}((q_0^{*}(\cdot), q^{*}(\cdot)), u^{*}(\cdot)),
\end{align*}
$$

which contradicts the (local) minimality properties of $(q^*(.), u^*(.))$ (cf. (4.3)-(4.4)).

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

Remark 4.2. A similar result can be obtained considering a slightly different geometric efficiency term, in which $\frac{q_1^2(T)}{q_0^2(T)}$ replaces $\frac{q_1(T)}{q_0(T)}$:

$$
\text{Minimize } \left[ -\frac{q_1^2(T)}{q_0(T)} + K|q_0(T) - E^*| \right].
$$

4.2. Some comments on normality for the Copepod model. In this section to simplify the notation, we write $(q(.), u(.))$ instead of $((q_0(.), q(.)), u(.))$. Consider now the following problem in which $x_T > 0$ is represents a lower bound for the Copepod displacement:

$$
(4.11) \quad \left\{
\begin{array}{ll}
\text{Minimize } \tilde{J}_{\mathcal{E}}(q(.), u(.)) := \left[-\frac{q_1^2(T)}{q_0(T)}\right] & \\
\text{over arcs } q(\cdot) = (q_0(\cdot), q_1(\cdot), q_2(\cdot), q_3(\cdot)) \in W^{1,1}([0,T]; \mathbb{R}^4) \text{ s.t.} & \\
(q_1, q_2, q_3)(t) = f(q(t), u(t)), \quad \text{for a.e. } t \in [0,T] & \\
q_0(t) = L(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) \ge x_T & \\
q_2(0) = q_2(T), \quad q_3(0) = q_3(T). &
\end{array}
\right.
$$

We consider the pseudo Hamiltonian (also referred to as ‘unmaximized’
Hamiltonian) $\mathcal{H}: \mathbb{R}^4 \times \mathbb{R}^4 \times \mathbb{R}^2 \rightarrow \mathbb{R}$ is the function

$$
\mathcal{H}(q, p, u) := \langle (p_1, p_2, p_3), f(q, u) \rangle - p_0 L(q, u).
$$

We are interested in establishing normality properties of the Copepod model when we minimize an efficiency cost like $\tilde{J}_{\mathcal{E}}$. More precisely we shall prove that the optimal control problem (4.11) does not allow abnormal minimizers. This is in particular valid for the case of interest of this report, in which we consider two different Lagrangians:

$$
(\text{I}) \quad L(q, u) = u_1^2 + u_2^2 \quad (\text{this is the sub-Riemannian problem}).
$$
(II) the case in which the Lagrangian represents the mechanical energy of the system.

(In either case, the Lagrangian $L$ 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)$.)

## The Maximum Principle

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

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

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

(iii) $\mathcal{H}(\bar{q}(t), p(t), \bar{u}(t)) = \max_{u \in U} \mathcal{H}(\bar{q}(t), p(t), u) = r$ for a.e. $t$, for some $r \in \mathbf{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^2(T)}{q_0^2(T)}$, and $p_1(T) = \alpha + 2\lambda \frac{q_1(T)}{q_0(T)}$, with $\alpha = 0$ if $\bar{q}_1(T) > x_T$.

Observe that, since $f$ and $L$ 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(.) \equiv \lambda \frac{q_1^2(T)}{q_0^2(T)}$ and $p_1(.) \equiv \alpha + 2\lambda \frac{q_1(T)}{q_0(T)}$) and

$$ (4.12) \qquad -\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}L(\bar{q}(t), \bar{u}(t)) \quad \text{a.e. } t \in [0, T] $$

$$ (4.13) \qquad -\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}L(\bar{q}(t), \bar{u}(t)) \quad \text{a.e. } t \in [0, T], $$

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.3.** We say that the (local) minimizer for (4.11) $(\bar{q}, \bar{u})$ is normal if the Maximum Principle applies with $\lambda \neq 0$, that is the multiplier $p_0$, associated with the state variable $q_0$, is non-zero. The minimizer $(\bar{q}, \bar{u})$ is called abnormal if the necessary conditions are valid with $\lambda = 0$ or, equivalently, $p_0 = 0$.

**Proposition 4.4.** All minimizers for (4.11) are necessarily normal.

### Proof.

**Step 1.** Let $(\bar{q}, \bar{u})$ be a minimizer for (4.11), and assume that the Maximum Principle applies with $\lambda = 0$ (i.e. $(\bar{q}, \bar{u})$ is abnormal). Then $\bar{q}_1(T) = x_T$.
Indeed, supposing by contradiction that $\bar{q}_1(T) > x_T$, and bearing in mind the transversality condition for the adjoint variable $p_1$, we would obtain that $\alpha = 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)).$$

We deduce that $p_2(.)$, $p_3(.)$ are the smooth periodic functions:
(4.14)

$$p_2(t) = -p_1 \frac{\sin(\bar{q}_2(t))}{2 + \sin^2(\bar{q}_2(t)) + \sin^2(\bar{q}_3(t))}, \quad p_3(t) = -p_1 \frac{\sin(\bar{q}_3(t))}{2 + \sin^2(\bar{q}_2(t)) + \sin^2(\bar{q}_3(t))}.$$

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) = x_T$.

**Step 2.** We introduce a new optimal control problem in which we replace the constraint inequality $q_1(T) \ge x_T$ by a penalty term in the cost to minimize.

(4.15)
$$
\left\{
\begin{array}{l}
\text{Minimize } \tilde{J}_{\varepsilon}((q_0(.), q(.)), u(.)) := \left[-\frac{q_1^2(T)}{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) = L(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_2(0) = q_2(T), \quad q_3(0) = q_3(T).
\end{array}
\right.
$$

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^2(T)}{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 (4.11), then it is a minimizer also for (4.15). Indeed, if $(\bar{q}, \bar{u})$ is a minimizer for (4.11), then

(4.16) $$ \tilde{J}_{\varepsilon}(\bar{q}(\cdot), \bar{u}(\cdot)) \le \tilde{J}_{\varepsilon}(q(\cdot), u(\cdot)), $$

for all trajectory-control pair $(q(.), u(.))$.

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

(4.17) $$ J'_{\varepsilon}(\hat{q}(\cdot), \hat{u}(\cdot)) < J'_{\varepsilon}(\bar{q}(\cdot), \bar{u}(\cdot)). $$
From the choice of $K$ we necessarily have that $\bar{q}_1(T) \ge x_T/2$. Consider the trajectory/control pair $(\tilde{q}(\cdot), \tilde{u}(\cdot))$ satisfying the following properties:

$$ (4.18) \qquad \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} $$

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

$$ (4.19) \qquad (\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{u}(2s-T) & \text{if } s \in (T/2, T] \end{cases} $$

(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) = \bar{q}_1(T) \ge x_T$, and a straightforward calculation provides

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

We would deduce that

$$ \begin{aligned} J_{\mathcal{E}}(\tilde{q}(\cdot), \tilde{u}(\cdot)) &= -\frac{\tilde{q}_1^2(T)}{\tilde{q}_0(T)} = -\frac{\hat{q}_1^2(T)}{\hat{q}_0(T)} \le J'_{\mathcal{E}}(\hat{q}(\cdot), \hat{u}(\cdot)) \\ &< 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 (4.11).

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

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

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

(iii) $\mathcal{H}(\bar{q}(t), p(t), \bar{u}(t)) = \max_{u \in U} \mathcal{H}(\bar{q}(t), p(t), u) = r$ for a.e. $t$, for some $r \in \mathbf{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^2(T)}{q_0^2(T)}$, and $p_1(T) \in 2(-\lambda \frac{q_1(T)}{q_0(T)} + [-K, 0])$.

From condition (ii) we know that $p_0(.)$ and $p_1(.)$ are constants, and system (4.12)-(4.13) are valid. Then, from the relations $(iv)'$ above we would deduce that $p_0 \equiv 0$ and $p_1 \equiv 0$, and therefore (4.12)-(4.13) 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 (4.11), then it would be an abnormal minimizer for (4.15) as well, but this is not admissible owing to step 3. $\square$
REFERENCES

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

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

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

T. Chambrion, L. Giraldi, A. Munnier. Optimal strokes for driftless swimmers: A general geometric approach.

M.J. Lighthill. Note on the swimming of slender fish. *J. Fluid Mech.* **9**, 305–317 (1960)

J. Lohéac, J.-F. Scheid. Time optimal swimmers and Brockett integrator. (2015) hal-01164561

D. Takagi. Swimming with stiff legs at low Reynolds number. Phys. Rev. E 92. (2015)

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