Xuan Vu

Caroline Chaux
email: caroline.chaux@univ-amu.fr

Nadège Thirion-Moreau

Sylvain Maire
email: maire@univ-tln.frthi-thanh-xuan.vu@lsis.org

Elfrida Mihaela Carstea

A new penalized nonnegative third order tensor decomposition using a block coordinate proximal gradient approach: application to 3D fluorescence spectroscopy

In this article, we address the problem of tensor factorization subject to certain constraints. We focus on the Canonical Polyadic Decomposition (CPD) also known as Parafac. The interest of this multi-linear decomposition coupled with 3D fluorescence spectroscopy is now well established in the fields of environmental data analysis, biochemistry and chemistry. When real experimental data (possibly corrupted by noise) are processed, the actual rank of the "observed" tensor is generally unknown. Moreover, when the amount of data is very large, this inverse problem may become numerically ill-posed and consequently hard to solve. The use of proper constraints reflecting some a priori knowledge about the latent (or hidden) tracked variables and/or additional information through the addition of penalty functions can prove very helpful in estimating more relevant components rather than totally arbitrary ones. The counterpart is that the cost functions that have to be considered can be non convex and sometimes even non differentiable making their optimization more difficult, leading to a higher Caroline Chaux is the corresponding author.

computing time and a slower convergence speed. Block alternating proximal approaches offer a rigorous and flexible framework to properly address that problem since they are applicable to a large class of cost functions while remaining quite easy to implement. Here, we suggest a new block coordinate variable metric forwardbackward method which can be seen as a special case of Majorize-Minimize (MM) approaches to derive a new penalized nonnegative third order CPD algorithm. Its interest, efficiency, robustness and flexibility are illustrated thanks to computer simulations carried out on both simulated and real experimental 3D fluorescence spectroscopy data.

Index Terms

Constrained optimization -Proximal approaches -Block alternating minimization -Nonnegative tensor factorization (NTF) -3D fluorescence spectroscopy I. INTRODUCTION Acquired data sets in numerous modern applications are now often organized into multi-way arrays of numerical values because they are obtained thanks to ever more performant acquisition, data transmission and massive data storage systems. This is typically the case for the spectrofluorimetry [START_REF] Lacowicz | Principle of fluorescence spectroscopy[END_REF] which is targeted in this article. This analysis technique is generally used to study the composition of solutions in order to detect the fluorescent chemical species (also called fluorophores) that are present. Once a coordinate basis is fixed, a tensor of order N can always be represented as an organized N -way array of numerical values. Since tensors and multi-linear algebra constitute a rigorous and natural mathematical framework for the formulation of many models and for the resolution of the problems that come along, they have been the subject of a growing interest in recent years.

The order of a tensor corresponds to the dimensionality of the array needed to represent it (i.e. the number of indices required to reference the elements of the array). In this article, we will focus on third order tensors, yet, the approach suggested here can be generalized to higher tensor orders. Then, the tensors can be processed which means either directly compressed or analyzed using a great variety of tensor decompositions.

The most popular one certainly remains the low tensor rank decomposition also known by a number of other names among which are Canonical Polyadic Decomposition (CPD 1 ), Candecomp [START_REF] Carroll | Analysis of individual differences in multi-dimensional scaling via n-way generalization of eckart-young decomposition[END_REF], CanD or Parafac (for PARAllel FACtor analysis) [START_REF] Harshman | Foundation of the Parafac procedure: models and conditions for an explanatory multimodal factor analysis[END_REF]. This factorization method constitutes an informative and compact model which has proven to be relevant in many application fields, including those in which we are most interested here, namely chemistry and chemometry [START_REF] Bro | Parafac: tutorial and applications[END_REF] [START_REF] Bro | Multi-way analysis in the food industry: models, algorithms and applications[END_REF], process analysis/monitoring [START_REF] Meng | On-line monitoring of batch processes using a parafac representation[END_REF] and environmental data mining [START_REF] Stedmon | Tracing dissolved organic matter in aquatic environments using a new approach to fluorescence spectroscopy[END_REF][8] [START_REF] Stedmon | Resolving the variability in dissolved organic matter fluorescence in a temperate estuary and its catchment using parafac analysis[END_REF]. The Polyadic Decomposition consists of decomposing a tensor into a sum of R rank-1 tensors (with R suitably large). When this number R of rank-1 terms is minimal, the decomposition is referred to as the Canonical Polyadic Decomposition and R is 1 It is the acronym that we will use in the rest of this article October 25, 2016 DRAFT called the rank of the tensor. However one of the main difficulties faced by practical fluorescence spectroscopy applications, is that the number of fluorophores is generally unknown and thus the actual rank of the "observed" tensor (i.e. the tensor of the acquired Fluorescence Emission Excitation Matrices, FEEM) is mostly unknown making the estimation of reliable fluorescent chemical compounds more difficult. Several works have been dedicated to the tensor rank estimation problem (see COre CONsistency DIAgnostic or CORCONDIA [START_REF] Bro | A new efficient method for determining the number of components in parafac models[END_REF] and Threshold-CORCONDIA [START_REF] Costa | Robust methods based on the hosvd for estimating the model order in parafac models[END_REF] to cite a few), however this problem is still hard and open. Consequently, we are interested here in new penalized methods ensuring the sparsity of the estimated loading matrices and thus insensitive to overestimation problems.

The CPD can also be regarded as a generalization of the matrix Singular Value Decomposition (SVD) to tensors or as a special case of another tensor decomposition known as the Tucker decomposition [START_REF] Tucker | Some mathematical notes on three-mode factor analysis[END_REF] by restricting the core tensor to be "diagonal". Another interesting advantage of this decomposition is its uniqueness under mild conditions [START_REF]Three way arrays: Rank and uniqueness of trilinear decompositions with apllication to arithmetic complexity and statistics[END_REF][14] [START_REF] Sidiropoulos | On the uniqueness of multilinear decomposition of n-way arrays[END_REF], involving that on perfectly multilinear, noiseless models of known rank, the use of constraints is not necessary.

When real experimental data are processed, there might be noises, the constituent vectors of the loading matrices (also known as the loading factors or vectors) can be close and the observed system can be dynamic (appearance or disappearance of compounds). The actual rank of the observed tensor is unknown and may be difficult to estimate. Moreover, when the amount of data becomes very large, this inverse problem may become numerically ill-posed [START_REF] Paatero | Construction and analysis of degenerate Parafac models[END_REF] and very difficult to solve [17][18]. The use of proper constraints reflecting some a priori knowledge about the latent or hidden variables that are tracked can be very helpful for the "unsupervised" estimation of reliable components. The whole processing chain can be rendered more robust and easier to use which is interesting when real-time monitoring or automated control systems are considered. Problems such as the automatic detection of water pollution could be tackled by this kind of approach.

Therefore, in the context of the 3D fluorescence spectroscopy analysis [START_REF] Coble | Characterization of marine and terrestrial dom in seawater using excitation-emission matrix spectroscopy[END_REF][20] [START_REF] Stedmon | Tracing dissolved organic matter in aquatic environments using a new approach to fluorescence spectroscopy[END_REF] (assuming the absence of errors coming either from the pre-processing used to remove Raman and Rayleigh scattering or possible bad settings of the devices) nonnegativity constraints should be considered given the physical nature of the hidden variables. In fact, in this particular application (see [START_REF] Royer | A regularized nonnegative canonical polyadic decomposition algorithm with preprocessing for 3D fluorescence spectroscopy[END_REF] for a reminder of the links that exist between 3D fluorescence spectroscopy and CPD), the loading vectors stand for physical quantities intrinsically nonnegative since they are related to emission and excitation spectra and concentrations through the samples acquired at different times in monitoring applications (or locations or pH for other kind of applications). The rank of the CP model that will approximate the "observed" tensor is closely linked to the number of fluorescent chemical compounds that are present in the studied samples.

The nonnegativity of the considered datasets as well as the one of the quantities that have to be estimated is also crucial for numerous other leading applications of CPD, especially those encountered in the area of October 25, 2016 DRAFT image processing (for example hyperspectral imaging [START_REF] Karami | Noise reduction of hyperspectral images using kernel non-negative tucker decomposition[END_REF][23], computer vision [START_REF] Wang | Compact representation of multidimensional data using tensor rank-one decomposition[END_REF][25], biomedical image processing [START_REF] Kopriva | Blind decomposition of low-dimensional multi-spectral image by sparse component analysis[END_REF][27] or functional magnetic resonance imaging (fMRI) for brain mapping [START_REF] Andersen | Structure-seeking multilinear methods for the analysis of fmri data[END_REF][29], etc.). It is the reason why it has given rise to numerous nonnegative CPD algorithms (see [START_REF] Cichocki | Non negative matrix and tensor factorizations: Application to exploratory multi-way data analysis and blind separation[END_REF] or [START_REF] Royer | A regularized nonnegative canonical polyadic decomposition algorithm with preprocessing for 3D fluorescence spectroscopy[END_REF] for example for an overview of those NTF algorithms).

If nonnegativity constraints are sufficient on known simulated models, sometimes they may not be for the correct estimation of latent variables in case of noise, model errors or real experimental data (i.e. unknown mixtures of an unknown number of compounds possibly corrupted by noise). When dealing with complicated scenarios, most algorithms are unable to identify the relevant components, leaving the end-user to decide which components have a chemical meaning. To automate decisions, algorithms can be helped to recover more reliable components. This is achieved by the use of additional information or constraints. These are the reasons why we are concerned with possibly sparse loading matrices and/or much more continuous than discontinuous loading factors. We will address these issues with the help of the addition of penalty/regularization functions. The counterpart is that the cost functions that have to be considered are not necessarily convex and sometimes are even non differentiable making their optimization more difficult. The introduction of constraints may also involve a higher computing time and a reduced convergence speed.

To properly address the problem of tensor factorizations subject to certain constraints, we suggest, here, to consider a Block Coordinate Variable Metric Forward-Backward (BC-VMFB) approach [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF] where the forward stage consists of a gradient step and the backward stage consists of a proximal step. This algorithm could thus also be called a Block Coordinate Proximal Gradient algorithm. The term "variable metric" means that a preconditioning is used. The aim of preconditioning is to increase the convergence rate. Such an approach has already been successfully used in the case of sparse non Negative Matrix Factorization (NMF) problems i)

without preconditioning and the resulting algorithm is denoted by PALM (for Proximal Alternating Linearized Minimization) [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] and ii) with preconditioning [START_REF] Repetti | A preconditioned forward-backward approach with application to large-scale nonconvex spectral unmixing problems[END_REF] (but we notice that in this case the authors were using an alphabet of known spectra for the endmembers). An unpreconditioned version has been proposed for third order tensor decompositions in [START_REF] Xu | A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion[END_REF], where both CPD and Tucker decomposition have been studied. NTF problem has also been solved recently in [START_REF] Zhang | Fast nonnegative tensor factorization based on accelerated proximal gradient and low-rank approxiamation[END_REF] using an accelerated projection gradient based algorithm. Finally, we recall that BC-VMFB approaches can be seen as a special case of block alternating Majorize-Minimize approaches.

MM approaches [START_REF] Hunter | A tutorial on MM algorithms[END_REF] offer a rigorous and flexible framework to properly address the problem of penalized non-negative third order CPD since they are applicable to a large class of cost functions while remaining quite easy to implement. This article is organized as follows: after a brief introduction, the notations are given and some recalls about constrained optimization are performed in Section II. We focus more specifically on proximity operators and block-coordinate forward-backward methods. This general framework is then used to tackle the problem of the penalized nonnegative canonical polyadic decomposition of third order tensors in Section III. To that aim, the model and the objective functions that are considered in this case are first introduced. Then, it is shown how the CPD problem can be integrated into the variational approaches framework. The gradient matrices are recalled and the preconditioning matrices that will be used are provided.

II. THEORETICAL BACKGROUND

A. Context

A way to address inverse problems is to formulate them under a variational approach. This means that we seek a solution to an optimization problem and more precisely to a minimization problem in which the involved functional to be minimized often consists of two terms: one linked to the noise properties, named "data fidelity term" and another one linked to a priori information on the target solution, named "regularization" [START_REF] Chaux | A variational formulation for frame based inverse problems[END_REF]. The minimization problem can thus be expressed as minimize

x∈R L F(x) Fidelity + R(x) Regularization . (1)
Concerning the data fidelity term (here denoted by F), one may choose e.g., in relation to the noise statistics that may corrupt the data, a quadratic term (linked to a Gaussian noise) or a Kullback-Leibler divergence (linked to a Poisson noise). The regularization R can aim at enforcing the solution sparsity possibly in a transformed domain (e.g. wavelet transform, Fourier Transform, time-frequency domain, etc.) [START_REF] Repetti | Euclid in a taxicab: Sparse blind deconvolution with smoothed l1/l2 regularization[END_REF] or at favoring a piecewise smooth behavior (Total-Variation). It can also model some hard constraints such as, for example, a range constraint or more precisely, a nonnegativity constraint. In this article, F and R are assumed to satisfy assumptions H 1 and H 2 stated as follows.

H 1 : F and R are proper 2 lower semi-continuous functions 3 .

H 2 : F is differentiable with a β-Lipschitz gradient 4 on the domain domR and its gradient is denoted by ∇F(x). In other words it means that F is sufficiently regular.

These technical assumptions play a prominent role in the algorithm's convergence proof.

In Problem (1), the regularization term R(x) may be split in a sum of J terms, leading to

R(x) = J j=1 R j (x) (2)
where for all j = 1, . . . , J, R j : R L → ]-∞, +∞] satisfies H 1 and H 3 .

H 3 : R j is assumed to be bounded from below by an affine function, and its restriction to its domain is continuous.

This allows to consider simultaneously various a priori on the solution (e.g. sparsity, nonnegativity, regularity, and so on).

Instead of performing the optimization on the whole set of unknowns at once, it is sometimes fruitful (either because of the intrinsic structure of the data, or performance issues or high size of the dataset) to compute updates and consider regularization on smaller sets of unknowns (called blocks in which the unknowns are stored 5 ).

The optimization is then performed on one block at a time. The different blocks are inspected according to a certain scanning rule, but after a given number of iterations all the blocs must have been swept. Such an idea is also used in randomization or stochastic optimization [START_REF] Mahoney | Randomized algorithm for matrices and data[END_REF][40] to be able to tackle the computational problems that occur when large scale datasets are processed. It can also be interesting if local instead of global constraints have to be enforced. If data blocks are denoted by x (j) , we thus have x = (x (j) ) 1≤j≤J and each block belongs to R Lj where J j=1 L j = L. The regularization term thus becomes

R(x) = J j=1 R j (x (j) ) (3)
where for all j = 1, . . . , J, R j : R Lj → ]-∞, +∞] satisfies H 1 and H 3 .

2 which means that domF and domR are non empty.

3 5 In the case of the CPD problem, we will see in the next section, that one "natural" partitioning of the unknowns is the one provided by the factor matrices themselves (such a partitioning is used in all popular alternate algorithms (ALS, HALS, fast HALS, etc.)), but other subsets of unknowns could be used October 25, 2016 DRAFT At this stage, it should be mentioned that:

F (resp. R) is lower-semicontinuous at x ∈ R L if, for every sequence (x i ) i∈N ∈ R L , lim i→+∞ x i -x = 0 ⇒ F (x) ≤ lim inf F (x i ) (resp. R(x) ≤ lim inf R(x i )). 4 which means that ∀(x, y) ∈ (domF ) 2 , ∇F (x) -∇F (y) ≤ β F (x) -F (y) where β ∈ ]0, +∞[ (β is called the Lipschitz constant)
1) the resulting functional to be minimized (Eqs. ( 1)-( 3)) is often non differentiable and is sometimes not convex which makes the minimization problem difficult to solve,

2) regarding the blocks, a strategy to scan them has to be adopted. Indeed, they can be checked either i) randomly, ii) cyclically or iii) quasi-cyclically, making sure to check each block at least once regularly.

These different scanning strategies play a role in the algorithm convergence proof but also in the algorithm convergence rate.

In this work, we will use a block alternating minimization method based on proximal gradient steps as described here after. Indeed, the Block Coordinate Variable Metric Forward-Backward algorithm (BC-VMFB) [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF], [START_REF] Chouzenoux | Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF] allows to tackle the general minimization problem described by Eqs. ( 1)-( 3) and additionally the convergence to a critical point is guaranteed. Furthermore, the introduction of preconditioners can help to accelerate the convergence. Before describing the algorithm, we now briefly review the notion of proximity operators.

B. An introduction to proximity operators

We mainly work, here, with functions ϕ satisfying H 1 (also denoted by ϕ ∈ Γ 0 (R)). Their proximity operator [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] is defined as

prox ϕ : R → R : v → arg min u∈R 1 2 u -v 2 + ϕ(u), (4)
prox ϕ ξ = sign(ξ) max{|ξ| -α, 0} with sign(ξ) =            +1 if ξ > 0 -1 if ξ < 0 0 if ξ = 0 (5)
where max{•, •} means the maximum of the two values given in argument. The proximity operator of the ℓ 1norm function reduces to a soft thresholding operation. Another example of interest is the proximity operator of the (squared) ℓ 2 -norm function. More precisely, let α > 0, and set ϕ : R → R : ξ → α|ξ| 2 . Then, for every

ξ ∈ R, prox ϕ ξ = ξ 2α + 1 . (6)
More recently, some authors [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF] introduced the notion of proximity operators associated with a Symmetric

Positive Definite (SPD) matrix P that is considering the metric induced by P in the definition of the proximity operator. First, we recall that ∀x ∈ R L , x 2 P = x, Px , where •, • is the inner product. The proximity operator of a function ϕ ∈ Γ 0 (R L ) associated with a SPD matrix P is defined as [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF] prox P,ϕ :

R L → R L : v → arg min u∈R L 1 2 u -v 2 P + ϕ(u). (7)
Note that if P reduces to the identity matrix, then the definition [START_REF] Stedmon | Tracing dissolved organic matter in aquatic environments using a new approach to fluorescence spectroscopy[END_REF] of the proximity operator associated with a SPD matrix reduces to the classical definition of the proximity operator given in (4).

C. Block-Coordinate Variable Metric Forward-Backward (BC-VMFB) algorithm

A proximal algorithm based on proximity operators associated with a SPD matrix P [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF] is presented here.

The considered minimization problem reads minimize

x∈R L F(x) + J j=1 R j (x (j) ). (8)
where, as aforementioned, F is smooth, satisfies H 1 , as well as each function R j : R Lj → ]-∞, +∞] for all j = 1, . . . , J which additionally satisfy H 3 .

The optimization problem defined in ( 8) can be solved using a Block-Coordinate Variable Metric Forward-Backward (BC-VMFB) algorithm whose principle is summarized below in Algorithm 1.

Algorithm 1 Block-Coordinate Variable Metric Forward-Backward (BC-VMFB) algorithm. Let P j k (x k ) be a SPD matrix // Construction of the preconditioner P j k (x k )

5:

Let ∇ j k F(x k ) be the Gradient // Calculation of Gradient 6:

x(j k ) k = x (j k ) k -γ k P j k (x k ) -1 ∇ j k F(x k )
// Updating of block j k according to a Gradient step 7:

x

(j k ) k+1 ∈ prox γ -1 k Pj k (x k ),Rj k x(j k ) k
// Updating of block j k according to a Proximal step 8:

x jk k+1 = x jk k where j = {1, ..., J} \ {j} // Other blocks are kept unchanged 9: end for

We can see that it consists mainly of two steps:

➊ the computation of the gradient of F (more precisely a partial gradient w.r.t. the chosen block j (or j k if it is considered at the k -th iteration of the algorithm (see Algorithm 1)),

➋ the computation of the proximity operator of each R j (j ∈ {1, . . . , J}) associated to the metric P.

It is important to notice that if the proximal stage ➋ (ensuring regularization and/or hard constraints) is removed, the algorithm goes back to an alternating preconditioned gradient method. The preconditioning matrix enables to accelerate the algorithm.

Finally, it has been proven that the convergence of this algorithm is guaranteed under assumptions [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF] that are recalled in Appendix.

III. APPLICATION TO THE PENALIZED NONNEGATIVE TENSOR FACTORIZATION PROBLEM

Due to our application in 3D fluorescence spectroscopy, we now focus on the problem of the penalized nonnegative CP decomposition. We explain how this problem can be naturally expressed in the form of ( 8), justifying the use of the BC-VMFB approach to tackle the CPD problem.

A. Nonnegative Canonical Polyadic Decomposition of third order tensors

The (Canonical) polyadic decomposition consists of decomposing an original tensor into a (minimal) sum of rank-1 6 terms. In the case of third order tensor, it reads:

T = R r=1 ā(1) r • ā(2) r • ā(3) r = [[ Ā(1) , Ā(2) , Ā(3) ]], (9)
where

T ∈ R I1×I2×I3 + is a nonnegative third order tensor, • is the outer product of vectors ā(n) r (n = 1, . . . , 3)
which are called the loading factors, and

ā(n) r = (a (n)
in,r ) in ∈ R In with i n = 1, . . . , I n for all n ∈ {1, 2, 3} and for all r ∈ {1, . . . , R}. The matrix

Ā(n) = [ā (n) 1 , ā(n) 2 , . . . , ā(n) R ] ∈ R In×R + is called n-th loading matrix.

Denote by T (n)

In,I-n ∈ R

In×I-n +

the matrix obtained by unfolding the tensor T in the n-th mode where the size

I -n is equal to I 1 I 2 I 3 /I n .
The previous model can be written in the matrix form as follows

T (n) In,I-n = Ā(n) (Z (-n) ) ⊤ , n ∈ {1, 2, 3} (10)
where

Z (-1) = Ā(3) ⊙ Ā(2) ∈ R I-1×R + , Z (-2) = Ā(3) ⊙ Ā(1) ∈ R I-2×R + , Z (-3) = Ā(2) ⊙ Ā(1) ∈ R I-3×R + ,
with ⊙ standing for the Khatri-Rao product [START_REF] Khatri | Solutions to some functional equations and their applications to characterization of probability distributions[END_REF] and (.) ⊤ denoting the transpose of a matrix. Solely from an observation T of T (e.g. T can be a noisy version of T ), our aim is to estimate the hidden variables i.e. the loading matrices Ā(n) for all n = 1, . . . 3. To reach that goal, we are going to express this problem under a variational approach. It will lead us to a minimization problem of the form given by ( 8), whose solution will constitute an estimate T of T .

B. Proposed approach

The tensor structure naturally leads, in the case of third order (resp. N -th order) tensor, to consider 3 (resp. N ) blocks corresponding to the loading matrices A (1) , A (2) and A (3) (resp. A (1) , . . ., A (N ) ). Thus, we propose to solve the following optimization problem with regularization for A (1) , A (2) and A (3) minimize

A (n) ∈R In×R ,n∈{1,2,3}
F(A (1) , A (2) , A (3) )+R 1 (A (1) ) + R 2 (A (2) ) + R 3 (A (3) ).

We opt for a quadratic data fidelity term F(A (1) , A (2) , 1) , A (2) ,

A (3) ) = 1 2 T -[[A (
A (3) ]] 2 F = 1 2 T (n) In,I-n -A (n) Z (-n) ⊤ 2 F ( 12
)
and the penalization terms R n (A (n) ) n∈{1,2,3} are defined by

R n (A (n) ) = In in=1 R r=1 ρ n (a (n) inr ) ∀n ∈ {1, 2, 3} (13)
where the loading matrices are defined element wise as

A (n) = (a (n)
inr ) (in,r)∈{1,...,In}×{1,...,R} and where

ρ n (ω) =      α (n) |ω| π (n) if η (n) min ≤ ω ≤ η (n) max +∞ otherwise (14)
and

α (n) ∈]0, +∞[, π (n) ∈ N * , η (n) min ∈ [-∞, +∞[ and η (n) max ∈ [η (n) min , +∞]
. This shows that, in the example developed in this article, the regularization parameters are chosen block dependent but are constant within a block (in our case a block will correspond to a loading matrix but other splittings could have been chosen).

Then, the minimization can be performed using Algorithm 1 with J = 3 (i.e. considering 3 blocks namely A (1) , A (2) and A (3) ). At the convergence (i.e. after k max iterations), the algorithm will provide estimations

A (n) of loading matrices Ā(n) for all n ∈ {1, 2, 3}.
As mentioned in the previous Section, algorithm 1 requires to compute the partial gradient matrices of F with respect to A (n) for all n = 1, . . . , 3 but also the proximity operators of each regularization term R n (A (n) ) n∈{1,2,3} associated to the metric P (n) .

We recall that the gradient matrices of F with respect to A (n) for all n = 1, . . . , 3, are defined as

∇ n F(A (1) , A (2) , A (3) ) = -(T (n) In,I-n -A (n) Z (-n) ⊤ )Z (-n) . (15)
In Algorithm 1, the gradient matrix ∇ n F(A (1) [k], A (2) [k], A (3) [k]) is replaced by a more compact notation

∇ n [k]
where k stands for the iteration step.

Another key point is to derive good preconditioners. Inspired by [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] and [START_REF] Repetti | A preconditioned forward-backward approach with application to large-scale nonconvex spectral unmixing problems[END_REF], we generalized the preconditioners proposed in [START_REF] Repetti | A preconditioned forward-backward approach with application to large-scale nonconvex spectral unmixing problems[END_REF] (dedicated to NMF) to the case of NTF problems. To this end, we adopted a similar approach using the n-th mode unfolding of the tensor given in (10) (i.e. a matrix form of the NTF problem).

The majorant function Q involved in the algorithm's convergence proof (see [START_REF] Zhang | Tensors methods for hyperspectral data processing: a space object identification study[END_REF] of Appendix) writes (1) [k], A (2) [k],A (3) [k]) = F(A (1) [k], A (2) [k], A (3) [k]) + trace((

Q n (A (n) |A
A (n) -A (n) [k]) ⊤ ∇ n [k]) (16)
+ 1 2 trace((A (n) -A (n) [k]) ⊤ (P (n) [k] (A (n) -A (n) [k]))), ∀n ∈ {1, 2, 3},
where denotes the Hadamard product between two matrices and P (n) [k] is a compact notation standing for 1) [k], A (2) [k], A (3) [k]).

P (n) (A (
The matrix P for the n-th block arising in ( 16) can be defined as follows

P (n) (A (1) , A (2) , A (3) ) = A (n) (Z (-n) ⊤ Z (-n) ) ⊘ A (n) , ∀n ∈ {1, 2, 3}, (17)
where ⊘ denotes the Hadamard division between two matrices.

Finally, the proximity operator associated with each (R n ) n∈{1,2,3} and corresponding to the computation of

prox γ[k] -1 P (n) [k],Rn ( A (n) [k]
) still has to be delineated. To rigorously define it, we need to vectorize the data 7 , that is

• Let define each element of the preconditioner matrix

P (n) [k] in (17) as (p (n)
i [k]) i∈{1,...,RIn} and let define the diagonal matrix

P (n) [k] as P (n) [k] = Diag(p (n) 1 [k], ..., p (n) RIn [k]), • Let A (n) [k] be stored in a vector a (n) [k] ∈ R RIn .
Then, the computation of prox γ

[k] -1 P (n) [k],Rn ( A (n) [k]) consists in computing prox γ[k] -1 P (n) [k],Rn ( a (n) [k])
whose definition is given in [START_REF] Stedmon | Tracing dissolved organic matter in aquatic environments using a new approach to fluorescence spectroscopy[END_REF].

The regularization being separable and given the shape of the preconditioning matrices P (n) [k], proximity operator can be expressed as

(∀y = (y (i) ) i∈{1,...,RIn} ∈ R RIn ) prox γ[k] -1 P (n) [k],Rn (y) = prox γ[k] -1 p (n) i [k],ρn (y (i) ) i∈{1,...,RIn} . (18)
For all i ∈ {1, ..., RI n }, we have [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] (∀υ ∈ R)

prox γ[k] -1 p (n) i ,ρn (υ) = min η (n) max , max η (n) min , prox γ[k]α (n) (p (n) i [k]) -1 | . | π (n) (υ) (19)
[Fig. 1

about here.]

The proximity operator of ρ n is illustrated in Fig. 1 for fixed parameter values α (n) = 2 and [η

(n) min , η (n) max ] = [0, 4].
We see that the positivity constraint is applied (negative values are projected to 0) but also that the maximum value η

(n) max is respected. Furthermore, when the ℓ 1 -norm is considered, we recognize the thresholding rule (up to 2 which is the value of the regularization parameter α (n) chosen here) which enforces data sparsity. Note that when the ℓ 2 -norm is considered, this is no longer the case.

For all the simulations that will be performed in the next section, we have considered π (n) ∈ {1, 2} (the proximity operator when π (n) = 1 (resp. when π (n) = 2) is recalled in Eq. ( 5) (resp. Eq. ( 6)). Regularization parameter α (n) and range constraint parameters being simulation dependent, they will be defined in the numerical simulation section.

C. Derived algorithm

At this stage, all the elements involved in the BC-VMFB algorithm suggested, here, to solve the penalized nonnegative CPD problem have been defined. The general scheme of the resulting penalized nonnegative CPD algorithm is depicted in Fig. 2.

[Fig. 2 about here.]

Finally, we notice that we have chosen to update the blocks according to a quasi cyclic rule. It means that the block that is updated is chosen randomly according to a uniform law (see [START_REF] Vervliet | A randomized block sampling approach to canonical polyadic decomposition of large-scale tensors[END_REF] for the importance of randomization in the choice of blocks). Moreover, we make sure that each block is chosen at least once every K iterations (with K = 100).

IV. NUMERICAL SIMULATIONS

In this section dedicated to computer simulations, we are going to show the interest of the BC-VMFB algorithm in the context of environmental data analysis. The method is applied to both synthetic (i.e. numerically simulated)

3D fluorescence spectroscopy data and raw data coming from an experiment of water quality monitoring. In order to better assess the performance of the suggested algorithm and to be able to compare it with other algorithms of the literature, two error indices are first introduced.

A. Error measures

When simulated data are used, the true tensor rank R is known, while the tensor rank that will be used for the model and thus for the decomposition is denoted by R. In the case of simulated data, we can consider the two error indices that have already been used in [START_REF] Vu | Study of different strategies for the canonical polyadic decomposition of nonnegative third order tensors with application to the separation of spectra in 3D fluorescence spectroscopy[END_REF] instead of the reconstruction error 8 T -T 2 F which is classically used with real experimental data. The first error index, denoted by E 1 , measures the error of estimation but discarding the over-factoring part. The second one, denoted by E 2 measures the error induced by the over-factoring part only. The purpose is to evaluate more accurately the quality of the results. In this case, 8 or the normalized reconstruction error

T -T 2 F T 2 F
, where T stands for the reconstructed (estimated) tensor.

the inherent indetermination of the CPD problem (i.e. permutation and scaling ambiguities) also have to be taken into account: if we apply to a given solution a permutation of their loading factors or an appropriate scaling of their loading matrices, the result is also a solution and defines the same decomposition. Therefore, permutation and scaling-independent measurements of the errors are necessary. For 3-way tensors, the loading matrices A (n) , for all n ∈ {1, . . . , 3}, are first normalized in such a way that each column of A (n) for n ∈ {1, ..., 2} is normalized in l 1 and each column of A (3) carries the weight.

The estimated normalized solutions of the CPD algorithm are denoted by A (n) . Then, the R column vectors of A (n) are permuted such that its Euclidean distance to A (n) is minimized. The permuted normalized estimates are denoted by

A (n) σ while A (n)
σ (1 : R) means that only its R first columns are considered and σ is the considered permutation. Thus, the estimation error E 1 is defined by

E 1 (σ) = 3 n=1 A (n) σ (1 : R) -A (n) 1 3 n=1 A (n) 1 (20) ⇒      E 1 = min σ E 1 (σ) or E 1dB = 10 log 10 (E 1 ) σ opt = arg min σ E 1 (σ) (21)
The over-factoring error E 2 concerns the remaining components R + 1, . . . , R in [START_REF] Stedmon | Resolving the variability in dissolved organic matter fluorescence in a temperate estuary and its catchment using parafac analysis[END_REF]. It is computed as follows

E 2 = R r=R+1
a (1) σopt,r • a (2) σopt,r • a (3) σopt,r 1

or E 2dB = 10 log 10 (E 2 ).

B. Synthetic case

In this section, we illustrate the behavior of the BC-VMFB algorithm on two synthetic data sets: the first one, is noiseless whereas the second one, is corrupted by a Gaussian noise (with a SNR = 17.6 dB). Here, the Signal to Noise Ratio (SNR) is defined as SNR = 20 log 10

T F T -T F
and the noisy data are truncated: all negative values are set to 0. In both cases, the tensor rank, R, is equal to 5, yet, it will be overestimated and the tensor decomposition will be performed considering that it is equal to 6 (i.e. R = 6). The purpose is to evaluate the robustness versus noise of the penalized nonnegative CPD algorithm suggested here but also its ability to overcome model errors. The performance of the proposed algorithm will also be compared with other algorithms of the literature namely the Bro's N -way algorithm [START_REF] Bro | Parafac: tutorial and applications[END_REF] and the fast HALS algorithm described in [START_REF] Phan | Fast alternating ls algorithms for high order candecomp/parafac tensor factorizations[END_REF].

The studied data sets are built as follows. The tensors T and T are 100×100×100. The excitation and emission spectra are created using sums of (shifted) density function drawn from a generalized normal distribution and can be either monomodal or bimodal. They are truncated to ensure the non-negativity of the spectra. The concentrations are generated according to a uniform law between 0 and 10.

For our simulations, regarding the BC-VMFB algorithm, we have used the following parameters: η (n) min ≡ η min = 2e -16 and η (n) max ≡ η max = 1000. The value of the exponent π (n) is chosen based on the maximum-likelihood method [52, p. 225] for each block and is equal to 1 for n = 1, 2 and to 2 for n = 3 in these tests. The stepsize γ[k] ≡ γ has been fixed to 0.99. Each time, we consider the two following cases: (i) with regularization parameter α (n) ≡ α = 0.05 for all blocks (constant within blocks) and (ii) without regularization i.e. with α = 0 for all blocks.

The estimated Fluorescence Emission Excitation Matrices (FEEM) are displayed in Fig. 3 (resp. Fig. 5) in the noiseless case (resp. in the noisy case for a SNR = 17.6 dB). The estimated emission and excitation spectra as well as the estimated concentrations are provided in Fig. 4 (resp. Fig. 6) in the noise-free case (resp. noisy case for a SNR = 17.6 dB).

[Fig. 3 The estimation error E 1 is found to be equal to -15dB (resp. -12.4dB, -11.2dB, -12.5dB) and the overfactoring error E 2 is found to be equal to -409dB (resp. 25.6dB, -409dB, 30.6dB) in the penalized noiseless case (resp. non penalized noiseless case, penalized noisy case and non penalized noisy case). In the noisy case, the final SNR is 32.3dB (resp. 31.3dB) in the penalized case (resp. non penalized case). We clearly observe that the penalized version of the BC-VMFB algorithm accurately estimates the loading matrices (even in the presence of noise) but is also able to overcome the problem of overfactoring by detecting the absence of a sixth component. In both noise-free and noisy cases, the regularized version of the BC-VMFB algorithm outperforms its non regularized version.

Finally, in the noisy, overestimated case, we compare the BC-VMFB algorithm ((i) with regularization α = 0.05 and (ii) without regularization) with two other classical NTF algorithms of the literature ((iii) the Bro's N -way algorithm under non-negative constraints and (iv) the fast HALS algorithm under non-negative constraints) both in terms of performance and computational cost per iteration. To that aim, we draw 100 initial values randomly from a uniform distribution U (0, 1). The following stopping conditions were used: (a) the number of iterations is either equal to k max = 10 5 or (b) the relative diminishing rate of the quadratic criterion

F [•+l]-F [•] F [•]
is smaller than a given tolerance tol = 10 -8 . Here l is equal to 1 for (iii) the N -way algorithm and (iv) the fast HALS algorithm, while for versions (i) and (ii) of the BC-VMFB algorithm, l corresponds to a block of 500 iterations. The obtained performance versus the different initializations are given in Fig. 7 (they are sorted in ascending order). The error index E 1 is displayed at the top of this figure while the (over-factoring) error index E 2 is displayed at bottom of this figure. The computation time per iteration for each algorithm is given in Table I.

[Fig. 7

about here.]

[TABLE 1

about here.]

We observe that in most cases (98% of the initializations that were tested here), the penalized version of the BC-VMFB algorithm outperforms all the other algorithms (no fake compound is estimated as it is clearly emphasized by the over-factoring error index E 2 ), yet the counterpart is that, in the noisy case, the performance measured thanks to the error index E 1 is slightly below to that obtained with the other algorithms among which is the non penalized version of the BC-VMFB algorithm (the observed small difference is due to the fact that the noise is only distributed on existing compounds since non existing compounds are enforced to be equal to 0). In the remaining cases (2% of the tested initializations), the algorithm exhibits the same behavior as the other algorithms: the chosen initializations seem to be too bad to be able to recover the true solution. We also notice in Table I that the BC-VMFB algorithm remains competitive regarding the computation time (since the fast-HALS and this implementation of the Bro's N -way algorithm can be regarded as the fastest algorithms of the literature).

C. Experimental case: a water monitoring campaign

In this section, we now test our penalized BC-VMFB method on a real experimental data set [START_REF] Carstea | On-line parafac analysis of fluorescence spectra[END_REF]. In this experiment, the data were acquired automatically every 3 minutes, during a 10 days monitoring campaign performed on water extracted from an urban river. The size of the baseline data set is: 36 × 111 × 2594. The excitation wavelengths range from 225nm to 400nm with a 5nm bandwidth, whereas the emission wavelengths range from 280nm to 500nm with a 2nm bandwidth. The FEEM have been pre-processed using the Zepp's method [START_REF] Zepp | Dissolved organic fluorophores in southeastern us coastal waters: correction method for eliminating rayleigh and raman scattering peaks in excitation-emission matrices[END_REF] implying that all remaining negative values were set to 0 (but it should have negligible impact on the overall analysis, since on random chosen FEEMs, less than 10 negative points out of 3996 were found (i.e. ≃ 0.25%). One example of the FEEM before and after pre-processing is shown in Fig. 8.

[Fig. 8

about here.]

Additional pre-processing were applied: the data corresponding to the first 6 emission slits were removed. The first 1200 FEEMs (data acquired during the first six days) were discarded too and 2 other FEEMs (corresponding to the 1737th and 1738th acquisition time) were suppressed due to noticeable acquisition problems. Finally, the size of data that are processed in this example is 36 × 105 × 1392. Different tensor ranks (i.e. different estimates of the number of chemical compounds) were tested: R = 4/5/6 but the obtained results are presented in two cases only ( R = 4 and 6). Concerning the penalized BC-VMFB method, a regularization π (n) ≡ π = 1 was applied on all loading matrices. When searching for R = 4 components (resp. R = 6 components), the following value α (n) ≡ α = 2 × 10 4 (resp. α = 8 × 10 3 ) were chosen for the regularization parameters for each loading factors. The reconstructed Fluorescence Excitation-Emission Matrices are displayed on the top of Fig. 9 for R = 4 compounds and at its bottom for R = 6 compounds. Results obtained with our method are given on the left part of this figure whereas there are given on its right part for the Bro's N -way algorithm.

The (scaled) reconstructed concentrations are given on Fig. 10 (on its left with our method and on its right with the Bro's N -way algorithm). The (scaled) reconstructed concentrations for R = 4 compounds are given on the top of Fig. 10 and on its bottom for R = 6 compounds.

[Fig. 9 Thanks to the obtained results, we are now able to conclude that only four fluorescent chemical compounds were present in the studied data set (the FEEM of the first and fourth estimated compounds in the case R = 6 are nearly null as well as their concentrations, see the left bottom part of the Fig. 9 and the left bottom part of the Fig. 10). During this experiment, a contamination with diesel oil [START_REF] Fan | Analysis of diesel fuels/biodiesel blends and identification of biodiesel using time-resolved laser-induced fluorescence spectroscopy (trlfs)[END_REF] appeared 7 days after the beginning of the monitoring campaign: it can be clearly observed on the Fig. 9. In fact, before the 1480th sample, one single fluorescent chemical compound was mainly present (compound labeled (1) in the case R = 4, or compound labeled [START_REF] Meng | On-line monitoring of batch processes using a parafac representation[END_REF] in the case R = 6) whereas two others occur in trace amounts (compounds labeled [START_REF] Harshman | Foundation of the Parafac procedure: models and conditions for an explanatory multimodal factor analysis[END_REF] and (4) in the case R = 4, or compounds labeled (2) and (3) in the case R = 6). After this time, a fourth compound occurs (compound labeled (2) in the case R = 4, or compound labeled (5) in the case R = 6) and we also observe an important increase of the concentrations of the two aforementioned compounds. Moreover, the concentration curves of those three compounds seem to exhibit a similar behavior. Even if our estimated spectra, concentrations and FEEM are sometimes close to those estimated with the Bro's N -way algorithm, they remain indeed different (it is particularly true for the concentrations and the FEEM). The main advantage of the BC-VMFB algorithm is that thanks to penalization the estimated spectra and concentrations are stable with respect to the tested ranks. It is not the case with the N -way algorithm. With our method, we were also able to decide that only four components were effectively present. We also observe the influence of the regularization parameters on the obtained results: the smallest regularization does not seem to lead to satisfying results in the case of 6 sought compounds, whereas the highest one is not adapted in the case of 4 sought compounds.

V. CONCLUSION

In this article, we addressed the problem of tensor factorizations subject to certain constraints (nonnegativity, sparsity, regularity, etc.). We tackled this problem within the broader framework of Block Coordinate Variable Metric Forward-Backward (BC-VMFB) approaches. The main interest of BC-VMFB approaches is to offer a clear theoretical and mathematical framework, since the conditions under which the sequence generated by this family of algorithms converges to a critical point of the objective function have been established in previous works of other authors. Through this general framework, we were able to derive a new penalized nonnegative third order CPD algorithm. In our case, the forward stage consists of a gradient step and the backward stage consists of a proximal step. Moreover, a preconditioning is also introduced in order to increase the convergence rate. Attention must still be drawn to the fact that some regularization parameters have to be set and that they may have an impact on the obtained results: that is why different strategies regarding the choice of the regularization terms that are added have been investigated. Computer simulations have been provided in order to enlighten the effectiveness and the robustness of the proposed approach in the applicative context of 3D fluorescence spectroscopy. Both simulated and real experimental data have been considered. Even if we only took advantage of a very small part of the enormous potential of the BC-VMFB approaches on those examples, we were able to illustrate some of their very interesting properties: reliability, robustness versus noise, good performance despite model errors and relative quickness. On real experimental data, identifying relevant components with traditional CPD algorithms is not always so straightforward, leaving the end-user to decide which components have a chemical meaning. The family of algorithms presented here can help to automate decisions. We focused on third order tensors but we have already extended these approaches to tensors of higher orders. The problem of possible missing data under the BC-VMFB framework will be addressed in future works.

APPENDIX

The BC-VMFB algorithm's convergence theorem reads Theorem A.1. [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF]Theorem 3.1] Assume that Assumptions A.1 hold and that F + R satisfies the Kurdyka-Łojasiewicz inequality9 [START_REF] Łojasiewicz | Une propriété topologique des sous-ensembles analytiques réels[END_REF], [START_REF] Kurdyka | w f -stratification of subanalytic functions and the Łojasiewicz inequality[END_REF]. Then the sequence (x k ) k∈N converges to a critical point x of F + R.

Moreover, (F(x k ) + R(x k )) k∈N is a non increasing sequence converging to F(x) + R(x).

It is guaranteed under the following assumptions [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF] Assumption A.1. 1) Let k ∈ N and let j k ∈ {1, ..., J}. The quadratic function defined as:

For every x (j k ) ∈ R j k , Q j k (x (j k ) |x k ) = F(x k ) + x (j k ) -x (j k ) k , ∇ j k F(x k ) + 1 2 x (j k ) -x (j k ) k 2 Pj k (x k ) (23)
is a majorant function of the restriction of F to its j k -th block on domR j k , i.e., for every x (j k ) ∈ domR j k ,

F(x (1)
k , . . . , x

(j k -1) k , x (j k ) k , x (j k +1) k , . . . , x (J) k ) ≤ Q j k (x (j k ) |x k ) (24)
Moreover, the eigenvalues of P j k (x k ) are lower and upper bounded by positive values.

2) Blocks (j k ) k∈N are updated according to an essentially cyclic rule, i.e., there exists K ≥ J such that, for every k ∈ N, {1, . . . , J} ⊂ {j k , . . . , j k + K -1}.

3) One of the following statements holds:

a) There exists

(γ 1 , γ2 ) ∈]0, +∞[ 2 such that, for every k ∈ N, γ1 ≤ γ k ≤ 1 -γ2 .
b) For every j ∈ {1, . . . , J}, R j is a convex function and there exists (γ 1 , γ2 ) ∈]0, +∞[ 2 such that, for every k ∈ N, γ1 ≤ γ k ≤ 2 -γ2 .

Gradient step

Proximal step Is stopping criterion reached ?

prox γ -1 P (n) [k],Rn A (n) [k] = A (n) [k] -γ∇ n [k] ⊘ P (n) [k] ( A (n) [k])
Other blocks unchanged at k matrices A (n)

Estimated loading

Outputs:

Eqs. ( 18)-( 19)

No k ← k + 1

where

  arg min means the position at which the function is minimum. The proximity operator can be viewed as a generalization of the projection operator. Indeed, when ϕ is the indicator function ι C of a nonempty closed convex subset C of R, i.e. it takes on the value 0 in C and +∞ in R \ C, prox ιC reduces to the projector Π C onto C. Explicit forms of this operator are known for numerous functions ϕ ∈ Γ 0 (R) [43][37][44]. Of particular interest in this work is the proximity operator of the ℓ 1 -norm function. More precisely, let α > 0, and set ϕ : R → R : ξ → α|ξ| (where | • | means the absolute value). Then, for every ξ ∈ R,

  about here.] [Fig. 4 about here.] [Fig. 5 about here.] [Fig. 6 about here.]

  about here.] [Fig. 10 about here.] Finally on Fig. 11, we compare the (scaled) estimated excitation and emission spectra obtained when we assume that R = 4 compounds are present, using either our BC-VMFB algorithm with different values of the regularization parameters ((b) α = 3 × 10 4 , (c) α = 2 × 10 4 and (d) α = 10 4 ) or the Bro's N -way algorithm with nonnegativity constraints ((a)). Whereas, on Fig. 12, the same study is performed considering that R = 6, (for the BC-VMFB algorithm: (b) α = 8 × 10 3 , (c) α = 7 × 10 3 and (d) α = 6 × 10 3 ). [Fig. 11 about here.] [Fig. 12 about here.]

2 Fig. 1 .

 21 Fig. 1. Plots of prox ρn (υ) in function of υ ∈ [-2, 22] when α (n) = 2 and [η (n) min , η (n) max ] = [0, 4].

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .

 23456 Fig.2. BC-VMFB algorithm to minimize[START_REF] Costa | Robust methods based on the hosvd for estimating the model order in parafac models[END_REF].

Fig. 7 .

 7 Fig. 7. Performance versus the different initializations sorted in ascending order in the noisy, overestimated case: error index E 1 (top), overfactoring error index E 2 (bottom).

Fig. 8 .

 8 Fig. 8. One of the 2594 original FEEM acquired during the monitoring campaign and the same FEEM after pre-processing thanks to the Zepp's method

Fig. 9 .Fig. 11 .Fig. 12 .

 91112 Fig. 9. Estimated FEEM using the penalized BC-VMFB algorithm (left), and Bro's N -way algorithm (right). Case R = 4 (top), case R = 6 (below)

  R + stands for the set of real nonnegative numbers and N for the set of natural numbers. The domain of R is denoted by domR: it is the set of values for which the function R is defined. • , • F and • 1 are respectively the norm, the Frobenius norm and the ℓ 1 -norm. Π C is the projector onto a closed subset C, whereas •, • is the inner product. The inverse of a matrix is denoted by (•) -1 , its transpose by (.) ⊤ . The outer product of vectors is denoted by •, ⊙ stands for the Khatri-Rao product, ⊘ for the Hadamard division between two matrices and for the Hadamard product between two matrices. The matrix trace is denoted by trace(.) and Diag(., . . . , .) is a diagonal matrix whose diagonal elements are those given in arguments.

	Finally, algorithms devoted to the problem
	of CPD under different constraints can be derived. In Section IV, computer simulations are performed in the
	specific context of 3D fluorescence spectroscopy. Both simulated and real experimental data are considered. The
	obtained results illustrate the interest and the good behavior of the proposed algorithm for such applications.
	The suggested algorithm is also compared to other existing approaches. Finally, a conclusion is drawn and
	perspectives are delineated.
	Notations: in this article, scalars, vectors, matrices and tensors are denoted by lower case (e.g. x), bold lower
	case (e.g. a), bold upper case (e.g. A) and calligraphic (e.g. T ) letters, respectively. R stands for the set of
	real numbers,

TABLE I COMPUTATION

 I TIME COMPARISON OF BC-VMFB IN TWO CASES: WITH OR WITHOUT PENALTY, WITH N-WAY AND FAST HALS USING THE SAME INITIAL VALUE IN NOISY DATA AND IN NOISELESS DATA.

		Elapsed time (s)	BC-VMFB without penalty	BC-VMFB with penalty	N-way	fast HALS
		For 50 iterations	0.2	0.2	11	0.5
	Noisy case	To reach stopping conditions	102	75	8	8
		(actual number of iterations)	(48500)	(36500)	(43)	(1856)
		(SNR,E 1 , E 2 ) dB	(31.3, -12.5, 30.6)	(32.7, -11.2, -409)	(31.3, -12.5, 30.6)	(31.3, -12.5, 30.6)
	Noiseless case	To reach stopping conditions	202	74	80	3.7
		(actual number of iterations)	(100000)	(36500)	(838)	(308)
		(RRE,E 1 , E 2 ) dB	(-75.1,-12.4,25.6)	(-44.7, -15, -409)	(-127.9,-8.7, 31.7)	(-63.9, -6.1, 31.7)
	October 25, 2016					DRAFT

October 25, 2016 DRAFT

A rank-1 tensor of order N (resp. 3) is defined as the outer product of N (resp. 3) vectors.

Note that this data vectorization is not applied in practice: for the implementation elementwise operations are performed thus avoiding memory issues.

Note that the Kurdyka-Łojasiewicz inequality holds for a wide class of functions[START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF].October 25, 2016 DRAFT

This work has been supported by Labex Archimède (ANR-11-LABX-0033).

The R = 6 emission spectra, excitation spectra and concentrations estimated with the BC-VMFB algorithm in two cases: without regularization and with regularization α = 0.05 in the case of noiseless data. In solid red line: the reference spectra, in dashed blue line: BC-VMFB without penalty and in dash-dot green line: BC-VMFB with penalty. . . . . . . . . . . . . . . . . . . . . 26 5 The FEEM of reference (left column) and the estimated FEEM using the BC-VMFB algorithm, in two cases: without regularization (middle column) and with regularization α = 0.05


P Bettiol

B Bonnard

A Nolot

J Rouot

OPTIMAL CONTROL THEORY AND THE SWIMMING MECHANISM OF THE COPEPOD ZOOPLANKTON

 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.

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

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 [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF] 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.

x 0 (t)

x ✓ 1 (t)

✓ 2 (t) 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, θ 1 , θ 2 are the respective angles of the two links, and they satisfy the constraint 0 ≤ θ 1 ≤ θ 2 ≤ π.

✓ = ⇡ ✓ = 0 ✓ 2 (t) ✓ 1 (t)
Using the swimming model at low Reynolds number, we relate the speed of the displacement to the speed of the shape variable θ by the equation

(1.1) ẋ0 = 2 i=1 θi sin(θ i ) 2 i=1 1 + sin 2 (θ i )
.

To parameterize the motion as a control system, one introduces the dynamics: θ1 = u 1 , θ2 = u 2 .

It provides a control system written as q = 2 i=1 u i F i (q) with q = (x 0 , θ), θ = (θ 1 , θ 2 ). Moreover we have state constraints given by a triangle T in the shape variables: θ i ∈ [0, π], i = 1, 2, and θ 1 ≤ θ 2 . u 1 , u 2 are periodic controls producing strokes, which are closed curves in the θ-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 = T 0 ( q M (θ) q)dt
where M is the matrix

M =   2 -1/2(cos 2 θ 1 + cos 2 θ 2 ) -1/2 sin θ 1 -1/2 sin θ 2 -1/2 sin θ 1 1/3 0 -1/2 sin θ 2 0 1/3  
Using (1.1) this amounts to minimize the quadratic form

E = T 0 (a ( θ)u 2 1 + 2b(θ)u 1 u 2 + c(θ)u 2 2 )dt with a = 1 3 - sin 2 θ 1 2(2 + sin 2 θ 1 + sin 2 θ 2 ) , b = - sin θ 1 sin θ 2 2(2 + sin 2 θ 1 + sin 2 θ 2 ) , c = 1 3 - sin 2 θ 2 2(2 + sin 2 θ 1 + sin 2 θ 2 )
.

One denotes g the associated Riemannian metric in the θ-space, the optimal control problem is rewritten as a sub-Riemannian problem q = 2 i=1 u i F i (q), min u(.) T 0 g(θ, θ)dt with appropriate boundary conditions associated with periodic control θ(0) = θ(T ), the triangle inequality constraints 0 ≤ θ 1 ≤ θ 2 ≤ π. 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 = T 0 g(θ, θ) 1/2 dt
and using the energy minimization point of view, the period T of a stroke can be fixed to T = 2π. 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π) -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) 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 E = x 2 0 (T )/E, which is introduced by [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF].

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 v to produce the displacement x 0 (T ) and the mechanical energy where the shape variables are hold on at θ(0), that is

(1.3) Ef f = |v| 2 M 11 (θ(0)) 1/T E , M = (M ij ), v = 1 T T 0 ẋ0 dt.
We have

Ef f ∼ x 2 0 (T ) E M 11 (θ(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.

A geometric analysis of the problem in the frame of SR-geometry

2.1. Geodesic computation. Consider the energy minimization problem

q = 2 i=1 u i F i (q), min u(.) T 0 (a(q)u 2 1 + 2b(q)u 1 u 2 + c(q)u 2 2 )dt
where the set of admissible controls U is the set of bounded measurable mapping valued in R 2 . We introduce the pseudo-Hamiltonian H(q, p, u) = p q + p 0 (a(q)u 2 1 + 2b(q)u 1 u 2 + c(q)u 2 ). According to the maximum principle (see [START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF]), minimizers are found among extremals curves, which are solutions of the following equations

q = ∂H ∂p , ṗ = - ∂H ∂q , ∂H ∂u = 0.
This leads to consider the following. : Normal case. Assume p 0 = 0 and it can be normalized to p 0 = -1/2. Corresponding controls are given by

u 1 = cH 1 -bH 2 ab -b 2 , u 2 = 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)

H n = 1 2(ac -b 2 ) (aH 2 2 -2bH 1 H 2 + cH 2 1 ).
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)

H 1 = H 2 = 0.
If F, G are two smooth vector fields, the Lie bracket is computed as

[F, G](q) = ∂F ∂q (q)G(q) - ∂G ∂q (q)F (q)
and if

H F = p • F (q), H G = p • G(q), the Poisson bracket is: {H F , H G }(z) = dH F ( -→ H G ) = p • [F, G](q)
. Differentiating twice (2.2) with respect to time, abnormal controls are given by (2.3)

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
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)) → Π(exp t -→ H n (z(0))) where Π is the standard projection: Π : (q, p) → 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 → 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 = sin(θ i ) ∆ ∂ ∂x 0 + u i ∂ ∂θ i , i = 1, 2. with ∆ = 2 i=1 (1 + sin 2 (θ i )). We get [F 1 , F 2 ](q) = f (θ 1 , θ 2 ) ∂ ∂x 0 with f (θ 1 , θ 2 ) = 2 sin(θ 1 ) sin(θ 2 )(cos(θ 1 ) -cos(θ 2 )) ∆ 2 .
Furthermore,

[[F 1 , F 2 ], F i ] = ∂ f ∂θ i ∂ ∂x 0 , 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 ] / ∈ 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 ≤ θ 1 ≤ θ 2 ≤ π 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 [START_REF] Takagi | Swimming with stiff legs at low Reynolds number[END_REF] corresponds to the policy:

θ 2 : 0 → π, θ 1 : 0 → π, θ i : π → 0, i = 1, 2 with θ 1 = θ 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)

2 i=1 sin(θ i ) ∆ dθ i = ∂ ∂θ 2 sin(θ 1 ) ∆ - ∂ ∂θ 1 sin(θ 2 ) ∆ dθ 1 ∧ θ 2 = dω.
2) The points where dω = 0 are precisely the abnormal triangle, and dω < 0 in the interior domain and dω > 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 [START_REF] El | Small sub-Riemannian balls on R 3[END_REF]. Consider a standard SR-problem (D, g), where ∆ is a distribution and g is a SR-metric, near a point q 0 ∈ 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 = span{ F , Ĝ} F = ∂ ∂x + y ∂ ∂z , Ĝ = ∂ ∂y -x ∂ ∂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 = F + yQ(w) ∂ ∂z , G = Ĝ -xQ(w) ∂ ∂z , w = (x, y) and Q is quadratic in w, Q = αx 2 + 2βxy + γy 2
, where α, β, γ are parameters. We introduce Q(w) = (1 + Q(w)) and one can write

F = ∂ ∂x + y Q ∂ ∂z , G = ∂ ∂y -x Q ∂ ∂z .
Computing, one has

[F, G] = (2 Q + x ∂Q ∂x + y ∂Q ∂y ) ∂ ∂z and using Euler formula [F, G] = 2(1 + 2Q) ∂ ∂z = f (w) ∂ ∂z .
Geodesics equations. We use Poincaré coordinates associated with the frame (F, G, ∂ ∂z ):

H 1 = p • F, H 2 = p • G, H 1 = p • ∂ ∂z and the normal Hamiltonian is H n = 1 2 (H 2 1 + H 2 2 ). One has Ḣ1 = dH 1 ( -→ H n ) = {H 1 , H 2 }H 2 = (p • [F, G])H 2 = p z f (w)H 2 with p z = H 3 is constant. Ḣ2 = dH 2 ( -→ 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) ds = p z f (w)dt.
Denoting by the derivative with respect to s, we get

H 1 = H 2 , 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 ] = ∂f ∂x ∂ ∂z , [[F, G], G] = ∂f ∂y ∂ ∂z .
The remaining equations to be integrated using (2.4) are

(2.5) x = H 1 p z f (w) , y = H 2 p z f (w) , z = (H 1 y -H 2 x) Q p z f (w) .
The solution can be estimated using the following expansion, associated with the gradation, setting

x = εX, y = εY, z = ε 2 Z, p x = p X , p y = p Y , p z = p Z /ε.
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 (2.6)

X(t) = X 0 (t) + εX 1 (t) + o(ε), Y (t) = Y 0 (t) + εY 1 (t) + o(ε), Z(t) = Z 0 (t) + εZ 1 (t) + o(ε).
Clearly by identification one gets that t → (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

x(t) = A λ (sin(λt + ϕ) -sin(ϕ)) (2.7) y(t) = A λ (cos(λt + ϕ) -cos(ϕ)) (2.8) z(t) = A 2 λ t - A 2 λ 2 sin(λt) (2.9)
where A, λ, ϕ 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. 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 θ(0) = (θ 1 (0), θ 2 (0)) in the interior of the triangle identified with 0 using the translation

x = θ 1 -θ 1 (0), y = θ 2 -θ 2 (0)
and using Taylor expansions one can approximate the displacement using

ż = u 1 sin(θ 1 ) + u 2 sin(θ 2 ) 2 + sin 2 (θ 1 ) + sin 2 (θ 2 ) = c 1 u 1 + c 2 u 2 + o(1)
where c 1 , c 2 are constant. Using ẋ = u 1 and ẏ = u 2 we set (2.10)

Z = z -c 1 x -c 2 y
to get an equation of the form

Ż = 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, θ 1 (0), θ 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 = θ 1π/4, y = θ 2 -3π/4 and z = x 0 , the Taylor expansion of order 4 of the vector fields F 1 , F 2 near (x, y) ∼ (0, 0) is given by

F 1 =     1 0 √ 2(16 x 3 -15 x 2 y-10 y 3 -39 x 2 +6 yx+6 y 2 +36 x+18 y+54) 324     F 2 =     0 1 √ 2(10 x 3 +15 y 2 x-16 y 3 +6 x 2 +6 yx-39 y 2 -18 x-36 y+54) 324    
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] = ϕ 1 (x, y, z) = [x, y, 162 √ 2 z-4 x 4 + 4 y 4 + 13 x 3 + 13 y 3 -18 x 2 + 18 y 2 -54 x -54 y], the resulting vector fields are

F 1 = [1, 0, -15 x 2 y -10 y 3 + 6 yx + 6 y 2 + 18 y] F 2 = [0, 1, 10 x 3 + 15 y 2 x + 6 x 2 + 6 yx -18 x]
(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] = ϕ 2 (x, y, z) = [x, y, z/36 - 5 yx 3
18 + 5 xy 3 18 -1/18 x 3 -1/6 x 2 y -1/6 xy 2 -1/18 y 3 ], the resulting vector fields are

F 1 = [1, 0, -5/4 x 2 y -1/6 x 2 -1/6 xy + y/2] F 2 = [0, 1, 5/4 xy 2 -1/6 xy -1/6 y 2 -x/2]
(3) we remove the terms of order 0 in F 13 et F 23 .

We consider the change of variables ϕ 4 given by

ϕ -1 4 (P ) =   X -α 1 Z + α 1 XY /2 + α 2 Y 2 /2 Y -α 2 Z -α 1 X 2 /2 -α 2 XY /2 Z + α 1 Y Z/2 -α 2 XZ/2 + high order terms   . More precisely, ϕ -1 4 (X, Y, Z) =         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 2 Y -1/24 XY 2 -1/24 Y 3 -1/12 XZ + 5 X 2 Z 144 + 1/12 Y Z + XY Z 72 + 5 Y 2 Z 144 -XZ 2 108 -Y Z 2 108         Remark 2.10. ϕ -1 4 = (X, Y, Z) + θ(X, Y, Z) + . . . so ϕ 4 (x, y, z) = (x, y, z) -θ(x, y, z) + . . .

The inverse transformation is

ϕ 4 (x, y, z) =      X = x -1/12 xy -x 2 y 72 -1/12 y 2 + y 3 72 + z/6 -1/18 yz Y = y + 1/12 x 2 + x 3 72 + 1/12 xy -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 -x 2 z 144 -1/12 yz -xyz 72 -y 2 z 144     
.

Applying the change of variables ϕ 4 and the feedback transformation

u 1 ← u 1 + u 2 (α 1 x + α 2 y), u 2 ← u 2 -u 1 (α 1 x + α 2
y), the resulting vector fields are

F 1 = [1, 0, y/2 + y 3 288 - x 3 72 - 359 x 2 y 288 + 1/48 xy 2 ], F 2 = [0, 1, -x/2 + y 3 72 - x 3 288 + 359 xy 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] = ϕ 5 (x, y, z) = [x, y, z + x 4 288 + 5 yx 3 16 -5 xy 3 16 -y 4 288 ], the resulting vector fields are

F 1 = [1, 0, - 89 x 2 y 288 - 89 y 3 288 + y/2 + 1/48 xy 2 ] F 2 = [0, 1, 89 x 3 288 + 89 xy 2 288 -x/2 -1/48 x 2 y]
in other words,

F 1 = F (-1) 1 + y 2 Q(x, y) ∂ ∂z F 2 = F (-1) 2 - x 2 Q(x, y) ∂ ∂z
where Q(x, y) = -89/144(x 2 + y 2 ) + 1/24xy and F (-1) 1

, F

(-1) 2 are the vector fields of order -1.

Remark 2.11. The change of coordinates doesn't let invariant the angular variables θ 1 , θ 2 . To relate the nilpotent system to the real one, we shall use the transformation

T (x, y, z) = ϕ 5 • ϕ 4 • ϕ 3 • ϕ 2 • ϕ 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 [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]. 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 = span{F, G} and F = ∂ ∂x + y 2 2 ∂ ∂z , G = ∂
∂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 + αy) 2 ∼ 1 + 2αy (order zero) c = (1 + βx + γy) 2 ∼ 1 + 2βx + 2γy (order zero)
where α, β, γ are parameters. Geodesic equations. We introduce the orthonormal frame

F 1 = F √ a , F 2 = G √ c , F 3 = ∂ ∂z
and denoting H i = p•F i , the normal Hamiltonian is given by

H n := 1/2(H 2 1 + H 2
2 ). We parameterize by arc-length: H 2 1 + H 2 2 = 1, H 1 = cos χ, H 2 = sin χ and H 3 = p z = λ constant, assuming λ = 0. Hence the geodesic equations become

(2.11) ẋ = cos χ √ a , ẏ = sin χ √ b , ż = y 2 cos χ 2 √ a , χ = 1 √ ac (yλ -α cos χ -β sin χ)
We introduce the new parameterization 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 α, β only.

Geometric application.

• The flat case is α = β = γ = 0 and corresponds to the standard pendulum. • In the generic case we have two subcases -β = 0 and (2.13) is integrable using elliptic functions.

-β = 0, due to dissipation, we are in the non-integrable case. Application to the copepod. In the integrable case: β = 0 models of periodic strokes with elliptic functions with modulus k are

• k = 0: circles • k 0.65: eight shape (Bernoulli lemniscates). Note that β = 0 is not a stable case, moreover the triangle constraint 0 ≤ θ 1 ≤ θ 2 ≤ π 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 ≤ θ 1 ≤ θ 2 ≤ π. They are:

• simple curves: circles x,

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

• eight shape: Bernoulli lemniscates z. Note that the orientation is imposed by Lemma 2.6.

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 [START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF]) 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 q = (q, q 0 ) and the corresponding dynamic

q = 2 i=1 u i F i (q) = F (q, u), u ∈ R 2 q0 = 2 i=1 u 2 i , q 0 (0) = 0.
and the problem is to minimize a cost of the form min u(.)

h(q(0), q(2π))

where the end-point conditions are of the form (q(0), q(2π)) ∈ C, where

C ⊂ R 3 × R 3 is a given closed set.
We denote p = (p, p 0 ) the extended adjoint vector. The pseudo-Hamiltonian takes the form

H(q, p, u) = 2 i=1 u i H i + p 0 (u 2 1 + u 2 2 )
and

H i = p • F i (q).
From [START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF], an optimal control pair (q, u) is satisfying the following necessary conditions that we split into two distinct parts Standard conditions.

q = ∂H ∂ p , ṗ = - ∂H ∂ q , ∂H ∂u = 0.
Transversality conditions.

(3.1) (p(0), -p(2π)) ∈ λ∇h(q(0), q(2π)) + N C (q(0), q(2π))

where N C is the (limiting) normal cone to the (closed) set C, (p, λ) = 0, λ ≥ 0.

Application. q = (x 0 , θ 1 , θ 2 )

• Maximizing the geometric efficiency with periodic condition

θ(0) = θ(2π), h = - x 2 0 (2π) E where E is the mechanical energy E = 2π 0 (u 2 1 + u 2 2 )dt.
In this case we deduce the periodicity condition on p θ = (p θ 1 , p θ 2 ) dual to θ:

(3.2) p θ (0) = p θ (2π),
to produce a smooth stroke in the normal case p 0 = 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.

h = c h = m, m = min h x 0 x 0 (t) = t 0 (u 2 1 + u 2 2 )dt
A c c e s s i b i l i t y s e t a t t i m e 2 π • Maximizing an efficiency depending on θ(0), with periodic conditions:

θ(0) = θ(2π), h = - m(θ(0))x 2 0 (2π) E
where m is a chosen smooth function. In this case, (3.2) becomes:

(3.3) p θ (0) -p θ (2π) = λ ∂h ∂θ(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 [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]).

Proposition 3.1. Let (x 0 (t), θ(t)), t ∈ [0, 2π] be a strictly normal stroke. Then a necessary optimality condition is the non existence of conjugate time t c ∈]0, 2π[. 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 2π 0 (u 2 1 + 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 θ(0), namely h = -x 2 0 (2π)m(θ( 0)) E where m(θ 1 (0)) = 2 -cos 2 (θ 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 θ-plane is symmetric with respect to the straight line (D) : θ 2 = -θ 1 + π, see Fig. 12.

Types of γ

x 0 (T ) L(γ)

x 0 (T )/L(γ)

Simple loops 1.400 × 10 -1 3.785 3.698 × 10 -2 1.700 × 10 -1 4.340 3.917 × 10 -2 2.000 × 10 -1 4.946 4.043 × 10 -2 2.100 × 10 -1 5.109 4.110 × 10 -2 Normal stroke 2.169 × 10 -1 5.180 4.187 × 10 -2 Fig. 10 2.200 × 10 -1 5.354 4.109 × 10 -2 2.300 × 10 -1 5.624 4.089 × 10 -2 2.500 × 10 -1 6.305 3.965 × 10 -2 2.740 × 10 -1 9.046 3.028 × 10 -2 Abnormal 2.742 × 10 -1 10.73 2.555 × 10 -2 Limaçon 2.000 × 10 -1 6.147 3.253 × 10 -2 Eight 2.000 × 10 -1 6.954 3.307 × 10 -2 Table 1. Value of the geometric efficiency for abnormal solution and different normal strokes for the Copepod swimmer. 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)

              
Minimize T 0 L(q(t), u(t))dt over arcs q(.) ∈ W 1,1 ([0, T ]; R 3 ) s.t. q(t) = f (q(t), u(t)), for a.e. t ∈ [0, T ] u(t) ∈ U for a.e. t ∈ [0, T ] q 1 (T ) = x T , q 1 (0) = 0 q 2 (0) = q 2 (T ), q 3 (0) = q 3 (T ) , in which x T > 0 is a given number which represents the Copepod desired displacement, T = 2π in the time period, U = R 2 , f : R 3 × R 2 → R 3 is the function f (q, u) := u 1 sin(q 2 ) 2 + sin 2 (q 2 ) + sin 2 (q 3 ) + u 2 sin(q 3 ) 2 + sin 2 (q 2 ) + sin 2 (q 3 ) , u 1 , u 2 .

Here, the state variable is q = (x, θ = (θ 1 , θ 2 )) in which x stands for the swimmer displacement, and θ = (θ 1 , θ 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 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 2 1 + u 2 2 . More precisely, suppose that we are given a minimizer (q * , u * ) for (4.1). Write E * := T 0 L(q * (t), u * (t))dt the corresponding energy and consider the following optimal control problem:

(4.2)                     
Minimize J E ((q 0 (.), q(.)), u(.)) := -q 1 (T ) q 0 (T ) + K|q 0 (T ) -E * | over arcs (q 0 , q(.)) ∈ W 1,1 ([0, T ]; R 4 ) s.t. q(t) = f (q(t), u(t)), for a.e. t ∈ [0, T ] q0 (t) = L(q(t), u(t)), for a.e. t ∈ [0, T ] u(t) ∈ U for a.e. t ∈ [0, T ] q 0 (0) = 0, q 1 (0) = 0 q 2 (0) = q 2 (T ), q 3 (0) = q 3 (T ) , Proposition 4.1. Let (q * (.), u * (.)) be a local minimizer for (4.1), then we can find a positive constant K such that ((q * 0 (.), q * (.)), u * (.)) is a local minimizer for (4.2).

Proof.

Assume that (q * (.), u * (.)) is local a minimizer for (4.1). Then, there exists ε 0 ∈ (0, x T /2) such that (4.3)

T 0 L(q * (t), u * (t))dt ≤ T 0 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) q * (.) -q(.)

L ∞ ≤ ε 0 , u * (.) -u(.) L ∞ ≤ ε 0 .
We claim that we can find ε 1 ∈ (0, ε 0 ) small enough such that (4.5) J E ((q * 0 (.), q * (.)), u * (.)) ≤ J E ((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) (q * 0 (.), q * (.)) -(q 0 (.), q(.))

L ∞ ≤ ε 1 , u * (.) -u(.) L ∞ ≤ ε 1 .
Suppose by contradiction that for each ε ∈ (0, ε 0 ) we can find a trajectory/control pair ((q ε 0 (.), q ε (.)), u ε (.)) satisfying all the requirements of the control system in (4.2), such that (4.7)

J E ((q ε 0 (.), q ε (.)), u ε (.)) < J E ((q * 0 (.), q * (.)), u * (.)) = -q * 1 (T ) q * 0 (T ) and (4.8) (q * 0 (.), q * (.)) -(q ε 0 (.), q ε (.))

L ∞ ≤ ε, u * (.) -u ε (.) L ∞ ≤ ε .
Then, we define ρ := q * 0 (T ) q ε 0 (T ) and take the trajectory/control pair ((q 0 (.), q(.)), ũ(.)) where ũ(.) := ρu ε (.), q0 (0) = 0, q1 (0) = 0 and q2 (.) = ρq ε 2 (.), q3 (.) = ρq ε 3 (.) . It is straightforward to see that ((q 0 (.), q(.)), ũ(.)) is a solution of the control system in (4.2), and, since q0 (T ) = q * 0 (T ), J E ((q 0 (.), q(.)) = -q1 (T ) q0 (T ) .

Observe that taking ε 1 ∈ (0, ε 0 ) small enough, then for all ε ∈ (0, ε 1 ), we obtain that (4.9) q * (.) -q(.)

L ∞ ≤ ε 0 , u * (.) -ũ(.) L ∞ ≤ ε 0 .
And choosing K > 0 big enough, we also have that

(4.10) q1 (T ) q0 (T ) - q ε 1 (T ) q ε 0 (T ) ≤ K|q ε 0 (T ) -q0 (T )| = K|q ε 0 (T ) -E * | .
(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 (q, ū) be a (local) minimizer for (4.11). Then, there exist a vector-valued function

p = (p 0 , p 1 , p 2 , p 3 ) ∈ W 1,1 ([0, T ]; R 4 ), α ≥ 0, β 1 ∈ R, β 2 ∈ R and λ ≥ 0 such that (i) λ + ||p|| L ∞ + |β 1 | + |β 2 | + α = 0; (ii) -ṗ(t) = ∂ q H(q(t), p(t), ū(t)) a.e.;
(iii) H(q(t), p(t), ū(t)) = max u∈U H(q(t), p(t), u) = r for a.e. t, for some r ∈ R;

(iv) p 2 (0) = p 2 (T ) = β 1 , p 3 (0) = p 3 (T ) = β 2 , p 0 (T ) = λ q 2 1 (T ) q 2
0 (T ) , and p 1 (T ) = α + 2λ q 1 (T ) q 0 (T ) , with α = 0 if q1 (T ) > x T . Observe that, since f and L do not depend on q 0 or q 1 , condition (ii) provides the following relations: ṗ0 (t) = 0, ṗ1 (t) = 0 a.e.

(which yield p 0 (.) ≡ λ q 2 1 (T ) q 2 0 (T ) and p 1 (.) ≡ α + 2λ q 1 (T ) q 0 (T ) ) and (4.12) -ṗ2 (t) = p 1 (ū 1 (t)∂ q 2 ϕ 1 (q(t))+ū 2 (t)∂ q 2 ϕ 2 (q(t)))-p 0 ∂ q 2 L(q(t), ū(t)) a.e. t ∈ [0, T ] (4.13) -ṗ3 (t) = p 1 (ū 1 (t)∂ q 3 ϕ 1 (q(t))+ū 2 (t)∂ q 3 ϕ 2 (q(t)))-p 0 ∂ q 3 L(q(t), ū(t)) a.e. t ∈ [0, T ] , in which ϕ 1 (q) = ϕ 1 (q 2 , q 3 ) := sin(q 2 ) 2 + sin 2 (q 2 ) + sin 2 (q 3 ) , ϕ 2 (q) = ϕ 2 (q 2 , q 3 ) := sin(q 3 ) 2 + sin 2 (q 2 ) + sin 2 (q 3 ) .

Definition 4.3. We say that the (local) minimizer for (4.11) (q, ū) is normal if the Maximum Principle applies with λ = 0, that is the multiplier p 0 , associated with the state variable q 0 , is non-zero. The minimizer (q, ū) is called abnormal if the necessary conditions are valid with λ = 0 or, equivalently, p 0 = 0.

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

Proof.

Step 1. Let (q, ū) be a minimizer for (4.11), and assume that the Maximum Principle applies with λ = 0 (i.e. (q, ū) is abnormal). Then q1 (T ) = x T .

Indeed, supposing by contradiction that q1 (T ) > x T , and bearing in mind the transversality condition for the adjoint variable p 1 , we would obtain that α = 0 and, therefore p 1 (.) ≡ 0. On the other hand the Maximality condition (iii) yields:

p 2 (t) = -p 1 ϕ 1 (q(t)), p 3 (t) = -p 1 ϕ 2 (q(t)) .
We deduce that p 2 (.), p 3 (.) are the smooth periodic functions: (4.14)

p 2 (t) = -p 1 sin(q 2 (t)) 2 + sin 2 (q 2 (t)) + sin 2 (q 3 (t))

, p 3 (t) = -p 1 sin(q 3 (t)) 2 + sin 2 (q 2 (t)) + sin 2 (q 3 (t))

.

Then from the formulae derived for p 2 (.), p 3 (.) we would also have p(.) ≡ 0, obtaining that (p(.), λ) = (0, 0). But, this would contradict the non-triviality condition (i) of the Maximum Principle. Then, we deduce that q1 (T ) = x T .

Step 2. We introduce a new optimal control problem in which we replace the constraint inequality q 1 (T ) ≥ x T by a penalty term in the cost to minimize.

(4.15)

                    
Minimize J E ((q 0 (.), q(.)), u(.)) := -q 2 1 (T ) q 0 (T ) + K max{-(q 1 (T ) -x T ); 0} over arcs (q 0 , q(.)) ∈ W 1,1 ([0, T ]; R 4 ) s.t. q(t) = f (q(t), u(t)), for a.e. t ∈ [0, T ] q0 (t) = L(q(t), u(t)), for a.e. t ∈ [0, T ] u(t) ∈ U for a.e. t ∈ [0, T ] q 0 (0) = 0, q 1 (0) = 0, q 2 (0) = q 2 (T ), q 3 (0) = q 3 (T ) .

Here we take K > 0 such that

K > 2C x T ,
in which C > 0 is a constant such that q 2 1 (T ) q 0 (T ) ≤ C, for all trajectory-control pairs (q(.), u(.)) .

(Observe that such a constant C always exists, take for instance C = √ 2π 2 .) We claim that if (q, ū) is a minimizer for (4.11), then it is a minimizer also for (4.15). Indeed, if (q, ū) is a minimizer for (4.11), then (4.16) JE (q(.), ū(.)) ≤ JE (q(.), u(.)) , for all trajectory-control pair (q(.), u(.)). Assume, by contradiction, that we can find a trajectory/control pair (q(.), û(.)) satisfying all the requirements of the control system in (4.15), such that (4.17) J E (q(.), û(.)) < J E (q(.), ū(.)) .

From the choice of K we necessarily have that q1 (T ) ≥ x T /2. Consider the trajectory/control pair (q(.), ũ(.)) satisfying the following properties: (Roughly speaking, employing a 'bigger' control, we construct a trajectory which, in the (q 2 , q 3 )-variables, does twice the path of (q(.), û(.)) on the same time interval.) As a consequence, we obtain q1 (T ) = q1 (T ) ≥ x T , and a straightforward calculation provides q2 1 (T ) q0 (T ) = q2 1 (T ) q0 (T ) .

We would deduce that J E (q(.), ũ(.)) = -q2 1 (T ) q0 (T ) = -q2 1 (T ) q0 (T ) ≤ J E (q(.), û(.))

< J E (q(.), ū(.)) = J E (q(.), ū(.)) , which contradicts the minimality of (q, ū) 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 (q, ū) is an abnormal minimizer for (4.15), then, there exist a vector-valued function p = (p 0 , p 1 , p 2 , p 3 ) ∈ W 1,1 ([0, T ]; R 4 ),

β 1 ∈ R, β 2 ∈ R and λ ≥ 0 such that (i)' λ + ||p|| L ∞ + |β 1 | + |β 2 | = 0;
(ii) -ṗ(t) = ∂ q H(q(t), p(t), ū(t)) a.e.;

(iii) H(q(t), p(t), ū(t)) = max u∈U H(q(t), p(t), u) = r for a.e. t, for some r ∈ R; (iv)' p 2 (0) = p 2 (T ) = β 1 , p 3 (0) = p 3 (T ) = β 2 , p 0 (T ) = λ q 2 1 (T ) q 2 0 (T ) , and p 1 (T ) ∈ 2(-λ 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 ≡ 0 and p 1 ≡ 0, and therefore (4.12)-(4.13) yield also p 2 ≡ 0 and p 3 ≡ 0. This would mean that (p(.), λ) = (0, 0), which contradicts the non-triviality condition (i) of the Maximum Principle. In conclusion, from step 2 above, if (q, ū) 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.

FIG. 1 .Figure 1 .

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

Figure 2 .

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

  φ the derivative of a function φ with respect to s, we get the equationsy = sin χ(1 + αy), χ = (yλα cos χβ sin χ)and the second order differential equation (2.13) χ + λ sin χ + α 2 sin χ cos χαβ sin 2 χ + βχ cos χ = 0.

Figure 3 .

 3 Figure 3. Simple candidates as strokes for the SR models.

Figure 4 .

 4 Figure 4. Accessibility set and transversality conditions with a cost function h.

Figure 5 .

 5 Figure 5. One parameter family of simple loops, limaçons and Bernoulli lemniscates normal strokes

Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .

 6789101112 Figure 6. Efficiency curve for simple loop normal strokes and efficiency of the abnormal stroke.

Figure 13 .

 13 Figure 13. Normal stroke tangent to the triangle.

  , q3 )(2s) if s ∈ [0, T /2] 2û(2s -T ) if s ∈ (T /2, T ] .

From (4.7) and (4.10) we would eventually deduce that J E ((q 0 (.), q(.)) = -q1 (T ) q0 (T )

≤ -q ε 1 (T ) q ε 0 (T ) + K|q ε 0 (T ) -E * | = J E ((q ε 0 (.), q ε (.)), u ε (.))

< J E ((q * 0 (.), q * (.)), u * (.)) , which contradicts the (local) minimality properties of (q * (.), u * (.)) (cf. (4.3)-(4.4)).

Remark 4.2. A similar result can be obtained considering a slightly different geometric efficiency term, in which

q 0 (T ) replaces q 1 (T ) q 0 (T ) :

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)

Minimize JE (q(.), u(.)) := -q 2 1 (T ) q 0 (T )

over arcs q(.) = (q 0 (.), q 1 (.), q 2 (.), q 3 (.)) ∈ W 1,1 ([0, T ]; R 4 ) s.t. ( q1 , q2 , q3 )(t) = f (q(t), u(t)), for a.e. t ∈ [0, T ] q0 (t) = L(q(t), u(t)), for a.e. t ∈ [0, T ] u(t) ∈ U for a.e. t ∈ [0, T ] q 0 (0) = 0, q 1 (0) = 0, q 1 (T ) ≥ x T q 2 (0) = q 2 (T ), q 3 (0) = q 3 (T ) .

We consider the pseudo Hamiltonian (also referred to as 'unmaximized' Hamiltonian)

We are interested in establishing normality properties of the Copepod model when we minimize an efficiency cost like JE . 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:

(I) L(q, u) = u 2 1 + u 2 2 (this is the sub-Riemannian problem).