Fundamental limits for nondestructive measurement of a single spin by Faraday rotation


Riccardo Messina

Weiliang Jin

Alejandro W Rodriguez

come

Exact formulas for radiative heat transfer between planar bodies under arbitrary temperature profiles:

Modified asymptotics and sign-flip transitions


Riccardo Messina

Weiliang Jin

Alejandro W Rodriguez

come

Strongly coupled near-field radiative and conductive heat transfer between planar bodies


Ghassen Dridi

Kaipeng Liu

Stéphane Guérin

Optimal robust quantum control by inverse geometric optimization

We develop an inverse geometric optimization technique that allows the derivation of optimal and robust exact solutions of low-dimension quantum control problems driven by external fields: we determine in the dynamical variable space optimal trajectories constrained to robust solutions by Euler-Lagrange optimization; the control fields are then derived from the obtained robust geodesics and the inverted dynamical equations. We apply this method, referred to as robust inverse optimization (RIO), to design optimal control fields producing a complete or half population transfer and a NOT quantum gate robust with respect to the pulse inhomogeneities. The method is versatile and can be applied to numerous quantum control problems, e.g. other gates, other types of imperfections, Raman processes, or dynamical decoupling of undesirable e↵ects.

Introduction.-Designing solutions of a system driven by external controls that are optimal with respect to practical costs such as area of the controls, energy or duration is a well known problem [START_REF] Glaser | Training Schrödinger's cat: quantum optimal control[END_REF]. The necessary conditions of optimality were established by Pontryagin via a maximum principle [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]. Based on this approach, various optimal quantum problems have been solved from low- [START_REF] Khaneja | [END_REF][4][START_REF]Introduction to Quantum Control and Dynamics[END_REF][START_REF] Bonnard | Optimal Control with Applications in Space and Quantum Dynamics[END_REF] to large-dimensional [START_REF] Rice | Optimal Control of Molecular Dynamics[END_REF][START_REF] Brif | [END_REF] systems.

Solutions that additionally feature robustness has become a major issue in quantum physics, especially in quantum information processing, where ultra highfidelity solutions are required (typically with relative errors below 10 4 ) [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]. Small imperfections in the design can cause fatal deviations of the performance. Robustness can be specifically taken in to account using adiabatic [START_REF] Vitanov | [END_REF]11], composite [12][13][14], combined [15] or shortcut to adiabaticity [16][17][18][19] techniques. However, these methods are not optimal and usually cost non-necessary energy and time.

Combining robustness constraints with the optimization methods has thus become a major challenge [20]. Numerical optimal control techniques based on time discretization with thousands of parameters to be optimized [21] can lead to very di↵erent results depending on the algorithm and the initial condition used. Alternative techniques involving from a few tens [22] to a few [23] parameters to be optimized have been developed, but they do not provide global optimal solutions in principle since they are based on restricted parametrizations. A recent proposal using Pontryagin's maximum principle (PMP) in an extended Hilbert space [24] allows an elegant integration of the robustness constraints, but leads to complicated systems to solve, only tractable for very simple targets, typically complete population transfers. A geometrical approach has been shown to provide optimal single-qubit phase gates [25]. All these methods use a direct optimization procedure, i.e. with the dynamical ⇤ sguerin@u-bourgogne.fr equations as constraints, which makes complicated the simultaneous integration of the robustness constraints.

In this Letter, we propose an alternative method of optimization based on inverse engineering: We apply in a first stage the Euler-Lagrange optimization [START_REF] Bertsekas | Constrained Optimization and Lagrange Multiplier Methods[END_REF] constrained by robustness integrals and by the boundaries ensuring exact fidelity in the dynamical variable space without invoking the dynamical equations. The constrained Euler-Lagrange optimization leads thus to a robust and exact geodesic. In a second stage, we derive the control field parameters from the geodesic and the dynamical equation, i.e. the time-dependent Schrödinger equation for the present quantum control problem, formulated in an inverted way, in which we express the controls from the dynamical variables.

We describe the technique by applying it to determine complete and partial transfers, as well as single-qubit gates that features all the desired properties with respect to the control pulse area (chosen as the cost): exactfidelity, optimality and robustness. The general applicability of the method is finally discussed with an emphasis on specific important quantum control problems.

Inverse engineering method in the dynamical variable space.-We can consider, without loss of generality, traceless Hamiltonians H = H 0 + V (in units such that ~= 1), where

H 0 = 1 2 ✓ ⌦ ⌦ ◆ , (1)
models the qubit {|0i, |1i} with the control parameters: the pulsed Rabi frequency ⌦ ⌘ ⌦(t) (considered positive without loss of generality) and the detuning ⌘ (t), and gathers unknown (time-independent) parameters representing the error in the description of the model. An error with respect to (i) the pulse amplitude (or area for a fixed pulse duration), i.e. the pulse inhomogeneities, is modeled with V = ⌦ x /2 and = ↵ corresponds to the relative deviation of the pulse amplitude; (ii) the detuning is modeled with V = z /2 and = corresponds to the deviation of the detun-ing. The solution | 0 (t)i = U 0 (t, t i )| 0 (t i )i of the time dependent Schrödinger equation (TDSE) associated to the qubit without errors, i~@ @t U 0 (t, t i ) = H 0 U 0 (t, t i ) with U 0 (t, t i ) the propagator from the initial time t i to time t, is conveniently parameterized with three angles, the mixing angle ✓ ⌘ ✓(t) 2 [0, ⇡], a relative ' ⌘ '(t) 2 [ ⇡, ⇡] and a global phase ⌘ (t) 2 [0, 2⇡], as

| 0 (t)i = ✓ cos ✓ 2 e i'/2 sin ✓ 2 e i'/2 ◆ e i /2 . ( 2
)
The TDSE can be equivalently rewritten as

✓ = ⌦ sin ', (3a)
' = + ⌦ cos ' cot ✓, (3b)
˙ = ⌦ cos ' sin ✓ = ✓ cot ' sin ✓ . ( 3c
)
The inverse-engineering method consists in determining the Hamiltonian elements (the controls) from the dynamical variables by inverting the TDSE:

H 0 = i~( @ @t U 0 (t, t i ))U † 0 (t, t i ), i.e. from inversion of Eqs. (3): = ' ˙ cos ✓, (4a)
⌦ = q ✓2 + ˙ 2 sin 2 ✓ = | ˙ | q ė ✓ 2 + sin 2 e ✓ (4b) with ė ✓ ⌘ d e ✓ d .
Since the optimization cost (the pulse area) and the integrals of robustness [see Eqs. ( 6)] depend on ✓(t) and (t), we can consider them as the dynamical variables providing a geometric representation of the problem, and the third variable '(t) is given by (3c), cot ' = ˙ sin ✓/ ✓, from which we obtain ' = ( ✓ ˙ sin ✓ ¨ ✓ sin ✓ ˙ ✓2 cos ✓)/( ✓2 + ˙ 2 sin 2 ✓). In the right part of Eq. (4b), we have assumed that one can write ✓(t) as a function of (t): e ✓( ) ⌘ ✓(t). In general we are led to consider multiple functions e ✓ j ( ), each depending on the time interval, in order to get the optimal solution (see for instance Fig. 3). We note that the pulse area from the initial t i to the final t f times (denoting i ⌘ (t i ), f ⌘ (t f ) and assuming a monotonic (t), such that ˙ (t) > 0) can be written as

A = Z t f ti dt ⌦(t) = Z f i d q ė ✓ 2 + sin 2 e ✓, (5)
which does not depend on the time-dependance of (t), but only on e ✓( ) (and its derivative). In this representation, it is thus relevant to consider trajectories e ✓( ) in the parameter space formed by the angles , ✓ and e '( ) ⌘ '(t) is given by (3c) for a given trajectory. We could have considered alternatively a representation with multiple functions e j (✓). However we will see that the chosen representation e ✓( ) is more convenient for the technical determination of the geodesics.

As a simplification of the presentation, we consider in this Letter robustness with respect to the pulse area (or identically to the pulse amplitude or inhonomegenities for a given time of interaction). It is shown below to allow one to consider the problem in the parameter space without invoking specific time parametrization: a robust optimal solution corresponds to a special trajectory e ✓( ) (the same when we consider optimization with respect to the pulse area, pulse energy or pulse duration for a given bound of its amplitude). The construction of the actual controls, i.e. the Rabi frequency ⌦ and the detuning , from (4) necessitates to use a specific time parametrization, (t), which can be chosen at will for the optimization with respect to the pulse area. On the other hand, the optimization with respect to the pulse duration corresponds to a specific time parametrization [see Eq. ( 8)].

We denote | (t)i the state of the complete dynamics including the error, solution of the TDSE i~@ @t | (t)i = H | (t)i. The single-shot shaped pulse method [17] allows one to define trajectories, in the dynamical variables space, resistant to errors. It can be formulated by a perturbative expansion of (t f ) with respect to ,

h T | (t f )i = 1 + O 1 + O 2 + O 3 + • • • ,
O 1 = i Z t f ti h 0 (t)|V (t)| 0 (t)idt ⌘ i Z t f ti e(t)dt, (6a) O 2 = ( i) 2 Z t f ti dt Z t ti dt 0 [e(t)e(t 0 ) + f (t) f (t 0 )], = 1 2 h Z t f ti dt e(t) i 2 Z t f ti dt f (t) Z t ti dt 0 f (t 0 ), (6b)
with e = 1 2 ( cos ✓ ↵ ˙ sin 2 ✓), and f = 1

2 h sin ✓ + ↵ ⇣ 1 2 ˙ sin 2✓ i ✓⌘i e i .
The other terms can be determined from a symbolic diagram [17]. We denote at a certain order n:

h T | ↵, (t f )i n = 1 + O 1 + O 2 + • • • + O n .
One notes the remarkable property that, assuming a function e ✓( ) and a monotonic (t) (such that ˙ (t) > 0), when one considers the robustness with respect to solely ↵ (i.e. = 0), then the integrals O n do not depend on the particular time-parametrization of (t) since

Z t ti e(t)dt = 1 2 ↵ Z t ti dt ˙ sin 2 ✓ = 1 2 ↵ Z i d sin 2 e ✓, (7a)
Z t ti f (t)dt = 1 2 ↵ Z i d ⇣ 1 2 sin 2 e ✓ i d e ✓ d ⌘ e i . ( 7b
)
The opposite situation of a decreasing (t) ( ˙ (t) < 0) would add a minus sign in the right hand sides of Eqs. [START_REF] Rice | Optimal Control of Molecular Dynamics[END_REF].

One can thus design an optimal robust trajectory e ✓( ) when = 0. The robust optimal time T min is determined from the particular time-parametrization (t) of the same trajectory e ✓( ) leading to a flat pulse (of chosen amplitude ⌦ 0 ):

T min = 1 ⌦ 0 Z f i d q ė ✓ 2 + sin 2 e ✓. (8)
We notice that robustness with respect to could be considered, from the trajectory e ✓( ) previously derived by exploiting the time-dependence of (t) (anticipated di↵erent from the one derived above for the time optimization).

The method, referred in short to as robust inverse optimization (RIO), is applied below for deriving first the optimal robust solution of two typical examples of population transfer: complete transfer and half coherent superposition (referred to as half transfer). We next consider a more complex target: the optimal robust quantum gate.

RIO for optimal robust population transfer. For the case of a population transfer to a target state | T i, the final global phase is not a priori fixed and is not robust since it is irrelevant. The figure of merit up to the third order reads

|h T | ↵, (t f )i| 2 = 1 + Õ2 , e O 2 = Z t f ti f (t)dt 2 . ( 9
)
The main result of the Letter is the following: The problem of optimal nullification up to the third order can be formulated as an optimization problem: finding the trajectory e ✓( ) that minimizes the pulse area ( 5)

A( e ✓) = Z f i d q ė ✓ 2 + sin 2 e ✓ ⌘ Z f i d L 0 ( , e ✓) (10)
under the constraints e O 2 = 0 rewritten for convenience as

1 ( e ✓) = 1 4 Z f i d (sin 2 e ✓ 2 e ✓) sin ⌘ Z f i d ' 1 ( , e ✓) = 1 2 (✓ f cos f ✓ i cos i ), (11a)
2 ( e ✓) = 1 4

Z f i d (sin 2 e ✓ 2 e ✓) cos ⌘ Z f i d ' 2 ( , e ✓) = 1 2 (✓ i sin i ✓ f sin f ) (11b)
with the initial state characterized by the angles (✓ i ⌘ ✓(t i ), ' i ⌘ '(t i ), i ) and the final target state (✓ f ⌘ ✓(t f ), ' f ⌘ '(t f ), f ). This can be solved by the Lagrange multiplier method extended to the function space as follows: The trajectory e ✓( ) is solution of

grad A( e ✓) + 1 grad 1 ( e ✓) + 2 grad 2 ( e ✓) = 0, ( 12
)
with j , j = 1, 2, the Lagrangian multipliers associated to the two constraints, where the gradient is defined according to the Euler-Lagrange equation (which is zero without constraint):

grad A( e ✓) = @L 0 @ e ✓ d d ✓ @L 0 @ ė ✓ ◆ , (13)
and similarly for grad j ( e ✓), j = 1, 2. We obtain the di↵erential equation

ë ✓ = 2 ė ✓ 2 cotan e ✓ + sin e ✓ cos e ✓ + ( 1 sin 2 cos ) ė ✓ 2 + sin 2 e ✓ 3/2 . ( 14
)
FIG. 1. Optimal robust geodesic e ✓( ) in the dynamical variable space ( , ✓) determined from numerical solution of ( 14 The optimal robust trajectory e ✓ opt ( ), solution of ( 14), is obtained for the set of values of 1 and 2 , which satisfies (11) (and we select the trajectory of smallest pulse area in case of more than one solution). We remark that the value of f results from this solution.

As a first example, we consider the complete population transfer from the ground state:

✓ i = 0, ✓ f = ⇡, i = ' i . Equation (3c), cot ' = ė sin ✓ = sin e ✓/ ė ✓, implies ' i = ' f = ⇡/2.
We can show that ė ✓ and ë ✓ are initially and finally infinite, and that the trajectory is symmetric about ✓ = ⇡/2 (which implies 1 = 1 sin f 2 cos f ). The obtained trajectory and the corresponding driving parameters (for a flat time-optimum pulse) are shown in Fig. 1. We notice that the optimal robust trajectory e ✓( ) is a more convenient representation since it forms a function contrary to the inverse one e (✓). We remark that we recover the robust optimal solution that has been derived in [24] by the PMP method in an extended Hilbert space.

We next focus on a (more complex) typical example of partial population transfer: the half transfer of internal (relative) phase ' 0 (up to an irrelevant global phase 0 /2), targeted from the ground state |0i:

| (t f )i = | T i ⌘ 1 p 2 ✓ e i'0/2 e i'0/2 ◆ e i 0 /2 . ( 15
)
This imposes the boundaries e ✓ opt ( ) and the corresponding driving parameters (for a flat pulse) are shown in Fig. 2. From this trajectory, we determine ė f = 0, which, from Eq. (3c), gives the optimal relative phase ' 0 = ⇡/2.

✓ i = 0, ✓ f = ⇡/2, i = ' i = ⇡/2, ' f = ' 0 ,
RIO for robust quantum gate -Achieving a quantum gate requires the additional robust control of the global phase . A traceless hamiltonian can only generate the SU(2)-type gate, which is taken as the targeted propa-

gator: U 0 = ✓ a b ⇤ b a ⇤ ◆ with |a| 2 + |b| 2 = 1.
The figure of merit usually adopted to determine the fidelity of a quantum gate is defined as F = 1 2 Tr(U † 0 U ) , where U is the actual propagator. Up to the third order, it involves then the real part of the integrals O 1 , O 3 (which are both zero), and O 2 :

F = 1 + R(O 2 ) = 1 1 2 h Z t f ti dt e(t) i 2 1 2 Z t f ti f (t)dt 2 , (16)
which is one when both R

t f
ti dt e(t) = 0 and

R t f
ti f (t)dt = 0. Robustness at third order can be thus expressed in this case (assuming ˙ (t) > 0) by ( 11) and the integral R f i d sin 2 e ✓ = 0. Since the argument of the latter is positive, the integral cannot be 0 and we conclude that (t) cannot be monotonic. We have thus to consider two (continuous) functions

⇢ e ✓ + ( ) for ˙ 0, = [ i , m ] e ✓ ( ) for ˙ < 0, = [ f , m [ (17)
with ✓ m = e ✓ + ( m ) = e ✓ ( m ) and the integrals to be

= 1 2 Z m i d e i (sin 2 e ✓ + 2 e ✓ + ) i(✓ f e i f ✓ i e i i ) 1 2 Z m f d e i (sin 2 e ✓ 2 e ✓ ). (18a) 0
We consider the typical NOT-type gate:

U NOT = ✓ 0 e i e i 0 ◆
, which means robust control to state |1i from the ground state |0i, i.e. ideally:

| (t i )i = |0i, | (t f )i = | T i ⌘ e i |1i
, with the phase  = (' f + f )/2, which should be additionally controlled in a robust way. This control implies the boundaries

✓ i = 0, ✓ f = ⇡, i = ' i = ' f = ⇡/2. The symmetric form f = i (giving  = ⇡/2), e ✓ ( ) = 2✓ m e ✓ + ( ), ✓ f = 2✓ m ,
for which the integral (18a) is automatically satisfied, is optimal. It leads to Eqs. (11) for the half transfer where f is replaced by m . This means that the optimal robust NOT gate is achieved when two consecutive optimal robust half transfers are achieved. The resulting trajectory is shown in Fig. 3. It corresponds to a minimum time T min ⇡ 8.1/⌦ 0 (for a flat pulse of Rabi frequency ⌦ 0 ) and to a pulse area of 2.58⇡.

Discussions and conclusions. -The technique of inverse optimization we have developed allows the design of optimal and robust solutions of quantum control problems of the general form H = H 0 + V , where robustness is meant with respect to . Its applicability necessitates the knowledge of the parametrization of the dynamics generated by H 0 and robustness is considered by perturbation of V . This includes low-dimensional dynamical symmetries for H 0 , typically SU(2), SU(3) [START_REF] Torrontegui | [END_REF] and SU(4) [28], for which dynamical invariants can be derived. But this does not limit the applicability to two-, three-or four-level systems; higher dimensions with specific symmetries can be considered [29]. The RIO method can treat robust optimization of the following multi-level problems:

(i) Stimulated Raman exact passage (STIREP) [18,19] featuring a SU(2) symmetry in the resonant case, or more generally a SU(3) symmetry;

(ii) Two-qubit gate, represented as a four-level problem in its simplest form: Following Ref. [30], one can compensate the error in the phase of a two-qubit controlled-PHASE gate (e.g. implemented in an ion trap) using interactions of the form T 1 ⌘ 1 2 x ⌦ x , T 2 ⌘ 1 2 x ⌦ y and T 3 ⌘ 1 2 1 1 ⌦ z , which feature SU(2) symmetry; (iii) Qudit gate (with an arbitrary dimension d), at the heart of quantum Fourier transform (a key ingredient of many quantum algorithms), in a multi-pod configuration with some overlapping controls [31] or with circulant symmetries [32].

The perturbation V is not limited to imperfections of the driving pulse, but can also concern the leakeage to undesirable states [33] or to a lossy environnement (for instance leading to dephasing noise [34]), where the latter problem takes the general form L = L 0 + V via Lindbladians. The application of the method can be then interpreted as a dynamical decoupling inverse optimization. Robustness and dynamical decoupling can also be treated simultaneously.

  FIG. 1. Optimal robust geodesic e ✓( ) in the dynamical variable space ( , ✓) determined from numerical solution of (14) corresponding to 1 ⇡ 1.11505, 2 ⇡ 0.30473 (leading to f = 5⇡/3). Inset: Resulting detuning and dynamics of the populations Pj, j = 1, 2, for robust timeoptimal control [obtained for a flat pulse of Rabi frequency ⌦0 according to (8)] showing the complete population transfer with the optimal time Tmin ⇡ 5.84/⌦0. The detuning has the form of a complete period of the elliptic cosine = 0 cn(4K(m)t/Tmin + K(m), m), t 2 [0, Tmin], with K(m) the complete elliptic integral of the first kind, m = 0.235, and 0 = 8K(m) p m/Tmin ⇡ 1.114⌦0.

  FIG. 2. Optimal robust geodesic e ✓( ) corresponding to 1 ⇡ 1.69741, 2 ⇡ 0.64653 (leading to f ⇡ 1.48⇡). Inset: Detuning and dynamics of the populations Pj, j = 1, 2, for robust time-optimal control (for a flat pulse of Rabi frequency ⌦0) showing the half superposition. We obtain the optimal time Tmin ⇡ 4.05/⌦0. The detuning has the form of three-quarters of the elliptic cosine period = 0 cn(3K(m)t/Tmin + K(m), m) with m = 0.4, and 0 = 6K(m) p m/Tmin ⇡ 1.66⌦0.

FIG. 3 .

 3 FIG. 3. Optimal robust geodesic e ✓( ) corresponding to the NOT gate. It corresponds to two consecutive symmetric trajectories of Fig. 2.

  where O n denotes the error term of total order n: O n ⌘ O( n ) and | T i the target state. The first two terms read

ACKNOWLEDGMENTS

We acknowledge support from the Investissements d'Avenir programs, project ISITE-BFC /IQUINS (ANR-15-IDEX-03), the EUR-EIPHI Graduate School (17-EURE-0002) and from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 765075 (LIM-QUET). KL acknowledges additional support from the CSC (China Scholarship Council). SG acknowledges discussions with D. Sugny.


Walid Djema
email: walid.djema@inria.fr

Keywords: Optimal control, turnpike properties, nonlinear systems, photobioreactor, strains selection, microalgae

The paper focuses on a generic optimal control problem (OCP) deriving from the competition between two microbial populations in continuous cultures. The competition for nutrients is reduced to a two-dimensional dynamical nonlinear-system that can be derived from classical quota models. We investigate an OCP that achieves species separation over a fixed time-window, suitable for a large class of empirical growth functions commonly used in quota models. Using Pontryagin's Maximum Principle (PMP), the optimal control strategy steering the model trajectories is fully characterized. Then, we provide sufficient conditions for the existence of a turnpike property associated with the optimal control and state-trajectories, as well as their respective co-state trajectories. Indeed, we prove that for a sufficiently large time, the optimal strategy achieving strain separation remains most of the time exponentially close to an optimal steady-state defined from an associated simpler static-OCP. This turnpike feature is based on the hyperbolicity of the linearized Hamiltonian-system around the solution of the static-OCP. The obtained theoretical results are then illustrated on microalgae, described by the Droop model in dimension 5. The optimal strategy is numerically computed in Bocop (open source toolbox for optimal control ) with direct optimization methods.

To cite this work: W. Djema, L. Giraldi, S. Maslovskaya, O. Bernard, 'Turnpike features in optimal selection of species represented by quota models'., Automatica, 2021.

Introduction

Natural selection drives the fate of the species that are permanently competing for limited survival sources [START_REF] Norberg | Biodiversity and ecosystem functioning: a complex adaptive systems approach[END_REF]. The individuals with a lower fitness in local conditions will eventually disappear from the ecosystem due to the competition for survival [START_REF] Williams | Adaptation and natural selection: A critique of some current evolutionary thought[END_REF]. Combined with mutations that may appear naturally over time, the Darwinian adaptation characterizing living organisms ensures that individuals which are better tailored to their local environment are continuously emerging and then progressively dominating their surroundings. In some instances, these evolution phenomena may be guided and redirected towards different fates. It is for instance possible to modify the outcome of species competition, leading to the emergence of enhanced organisms, by controlling and modulating the natural selection pressure [START_REF] Swenson | Artificial ecosystem selection[END_REF]. Such a selection process was then tamed for generating naturally more efficient microogranisms [START_REF] Almeida | Selection and application of microorganisms to improve oil recovery[END_REF]. Within this scope, some microalgae strains have been recently improved by setting-up a Darwinian selection pressure within artificial cultures such as photobioreactors (see, e.g., [START_REF] Gachelin | Enhancing PUFA-rich polar lipids in Tisochrysis lutea using adaptive laboratory evolution (ALE) with oscillating thermal stress[END_REF][START_REF] Huertas | Estimating the capability of different phytoplankton groups to adapt to contamination: herbicides will affect phytoplankton species differently[END_REF][START_REF] Bonnefond | Continuous selection pressure to improve temperature acclimation of Tisochrysis lutea[END_REF]). Microalgae attract considerable attention since they have a great potential to produce less impacting biofuels ( [START_REF] Chisti | Biodiesel from microalgae[END_REF][START_REF] Wijffels | An outlook on microalgal biofuels[END_REF]). In addition to reducing atmospheric CO 2 , they are also considered as promising elemental bricks for green chemistry and profitable new protein sources. Moreover, they are even used to improve wastewater treatment (see also, [START_REF] Walsh | Algal food and fuel coproduction can mitigate greenhouse gas emissions while improving land and water-use efficiency[END_REF][START_REF] Odjadjare | Potential biotechnological application of microalgae: a critical review[END_REF]). However, despite their large interest, microalgae are still in their infancy and hitherto mainly grown for cosmetic or pharmaceutical industries. Their exploitation on a larger scale, particularly for biofuel production, is conditioned by the overcoming of various technical challenges, including the necessity of enhancing the strain productivity since mainly wild species are currently cultivated. This explains the worldwide race during the last decade for selecting more suitable microalgae strains (see, e.g., [START_REF] Aravantinou | Selection of microalgae for wastewater treatment and potential lipids production[END_REF][START_REF] Rodolfi | Microalgae for oil: Strain selection, induction of lipid synthesis and outdoor mass cultivation in a low-cost photobioreactor[END_REF][START_REF] Yoo | Selection of microalgae for lipid production under high levels carbon dioxide[END_REF]).

In this paper, the objective is to develop an effective control strategy ensuring that the most suitable strains will dominate the continuous artificial culture in finite time. For instance, in biofuel production, desirable microalgae strains can be those which grow faster, with higher lipid content, and which are also well adapted to local environment and seasons.

We recall that when n distinct strains coexist in a chemostat with limited substrate, a wellknown principle -the competitive exclusion principle CEP [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF]-predicts that the species requiring the smallest amount of nutrients to grow at a rate equal to the dilution rate wins the competition, while the other species will be asymptotically washed out. However, depending on the culture conditions, the winner of the competition is not in general necessarily the one with the highest lipid content or unrivalled productivity. Therefore, our goal is to determine an optimal control strategy so that the species of interest is the one that dominates in finite time. To achieve our objectives, we use an approach that is based on optimal control theory (Pontryagin's principle [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]) and turnpike features (see, e.g., [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]). There are other works for optimal selection of microorganisms in minimal-time ( [START_REF] Bayen | Optimization of the separation of two species in a chemostat[END_REF][START_REF] Bayen | Optimization of strain selection in evolutionary continuous culture[END_REF]), based on the simpler Monod-type dynamics describing the growth of bacteria. The microalgae behavior is however more complex since the growth of these microorganisms does not depend on the substrate concentration in the photobioreactor (like in the Monod model), but rather on the quantity of nutrients stored inside the cells [START_REF] Bernard | Modelling of microalgae culture systems with applications to control and optimization[END_REF][START_REF] Bernard | Hurdles and challenges for modelling and control of microalgae for CO 2 mitigation and biofuel production[END_REF]. The Droop model ( [START_REF] Droop | Vitamin B 12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri[END_REF]) is widely used to represent algal growth, since it takes into account the additional dynamics governing the internal cell-quota storage [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF][START_REF] Droop | 25 years of algal growth kinetics, a personal view[END_REF]. Here we use a general form of quota models, extending the classical Droop's formulation to include different forms of growth rate. In particular, we suppose in this work that an ad'hoc controller is able to track the desired nutrient concentration in the medium. As a consequence, the classical 5 dimensional quota model (studied for instance in [START_REF] Djema | Optimal control separating two microalgae species competing in a chemostat[END_REF][START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF]) reduces to a system of dimension 2. In addition, this work generalizes the approach in [START_REF] Djema | An Optimal Control Strategy Separating Two Species of Microalgae in Photobioreactors[END_REF][START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF], not only because it considers a generic class of systems relevant for microalgae (using different growth rates), but also because there are no Inria assumption on the initial conditions, which are often assumed to be already confined to an attractive positively invariant set ([4, 18, 20]). The main novelty of our work is to prove the existence of a turnpike property characterizing the optimal solutions. The turnpike has been recently revisited in the optimal control literature and it is attracting more and more attention (see, e.g., [START_REF] Trélat | Integral and measure-turnpike properties for infinite-dimensional optimal control systems[END_REF][START_REF] Porretta | Long time versus steady state optimal control[END_REF][START_REF] Zaslavski | Turnpike theory of continuous-time linear optimal control problems[END_REF][START_REF] Trélat | Integral and measure-turnpike properties for infinite-dimensional optimal control systems[END_REF][START_REF] Lohéac | Minimal controllability time for finite-dimensional control systems under state constraints[END_REF][START_REF] Gugat | On the turnpike phenomenon for optimal boundary control problems with hyperbolic systems[END_REF]). The turnpike property states that the optimal solutions consist of three pieces: two transient short-time arcs at the boundaries of the time interval, and a long-time behavior remaining exponentially close to a static solution (see, e.g., [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]). From a practical standpoint, turnpike features are also advantageous for deriving sub-optimal controls that are easy to implement on real industrial systems.

The paper is organized as follows. Sect. 2 is devoted to the statement of the general optimal control problem (OCP), formulated for a generic two-dimensional dynamical system, which includes competition between algal species. In Sect. 3 we state the main theorem dealing with the general OCP. Its proof is provided in Sect. 4, which relies on Pontryagin's maximum principle (PMP, [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]) and also on the existence of a turnpike property along the optimal solutions. In Sect. 5, we focus on an application involving a variable yield model that includes the classical Droop's model. Thus, based on the results in Sect. 3-4, we solve the problem of microalgae selection under several growth-rates. Finally, in Sect. 5.5, the optimal control is determined in different cases and the turnpike feature is illustrated using a direct optimization method implemented in Bocop (optimal control toolbox, [START_REF] Bonnans | BOCOP: an open source toolbox for optimal control -A collection of examples[END_REF]).

Optimal control problem (OCP)

We consider a generic nonlinear two-dimensional dynamical system, suitable for competition between two populations, in which the dynamics are coupled through the feeding control-flow. Typically, the competition may hold between two species (e.g., phytoplankton, bacteria, yeasts,...) for a limiting substrate (w), whose concentration is controlled in a continuous bioreactor, can be represented as follows:

ẏ1 = -f 1 (y 1 )y 1 + g 1 (w), ẏ2 = -f 2 (y 2 )y 2 + g 2 (w), (1)
where the f i (•) are regular nonnegative real-valued and increasing functions representing the growth rate of the i-th species with respect to the intracellular limiting nutrient y i . The functions g i (•) are also regular nonnegative real-valued and increasing functions, standing for the flux of limiting nutrient into the cell.

We highlight that the system (1) typically appears as a subsystem of quota models, where y i is the internal quota of the i-th species, i.e. the intracellular amount of the element for which the involved species are competing. The Droop's model is extensively studied later in Sect. 5 as the most famous example of this class of systems. However, the model (1) can also be found in more complicated metabolic dynamical systems ( [START_REF] Wiechert | Modeling and simulation: tools for metabolic engineering[END_REF]) where the intracellular element y i is diluted by growth.

The resulting system (1) is positive, i.e., for positive initial conditions the trajectories remain positive.

In general, the dynamics of the system (1) is at higher dimension (at least 5 in the case of Droop's model for instance). However, the reduced system is derived by assuming that the limiting substrate has been controlled by an ad'hoc controller tracking the concentration w(t). Indeed, the way this subsystem can be extracted from the full model is developed in Sect. 5 in the typical case of Droop's model describing microalgae.

Without loss of generality, we consider throughout this paper that species 1 is the one of interest and must become predominant in the population in finite time. Our objective is to investigate the control strategies for this competition between two initially co-existing populations by controlling the nonnegative and bounded input w over a fixed time-interval [0, T f ].

Statement of the general OCP

The set of admissible controls associated with system (1) is defined by w ∈ W, s.t., W is a subset of L ∞ loc (R + ), the space of locally integrable functions on every compact set on R + . More precisely, we set,

w ∈ W = w : [0, T f ] → [0, w m ] | w(•) ∈ L ∞ loc (R + )
, where w m is a sufficiently large positive constant (the maximum influent substrate concentration in the bioreactor).

Then, we define now the main OCP that reads,

Maximize w∈W C [0,T f ] = T f 0 f 1 (y 1 ( )) -f 2 (y 2 ( ))d
, where, y = (y 1 , y 2 ), y(•) is the solution of (1) corresponding to the control w(•).

(

) 2
In what follows, the terminal conditions of the states y i , i = 1, 2, are assumed to be free and their initial conditions are positive. Besides, let us define (y * , w * ), the steady state of (1), where, y * > 0 and w * ∈ [0, w m ]. It follows that, f i (y * i )y * i = g i (w * ), for i = 1, 2, and since the map, y → f i (y i )y i , is increasing, we get, w * = δ i (y * ), where,

δ i (y) = g -1 i (f i (y i )y i ). Thus, y * i (w * ) = δ -1 i (w * ), i = 1, 2, for all fixed w * ∈ [0, w m ].
Then now, we define for all admissible y * i > 0 and w * ∈ [0, w m ] the function,

∆(w * ) = f 1 (δ -1 1 (w * )) -f 2 (δ -1 2 (w * )). (3)
Throughout this work, we assume that the maximizer of ∆ is the unique point w ∈ (0, w m ), as illustrated in Fig. 1. The cases where w is zero or w m are pointless and can be trivially addressed outside this framework. Let us remark that under these considered assumptions, OCP (2) has a maximizer using Filippov Theorem (see, e.g., [START_REF] Cesari | Optimization -Theory and applications[END_REF]). In fact, the generic form of the functions f i (δ -1 i (•)), involving one intersection point over (0, w m ) as in Fig. 1, ensures the existence of control functions maximizing and minimizing the criterion C [0,T f ] , as a direct consequence of the competitive exclusion principle (CEP), (see, e.g., [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF]).

Pontryagin's maximum principle (PMP)

The Hamiltonian associated with OCP (2) is defined by,

H =Φ(λ y1 , λ y2 , w) -λ y1 f 1 (y 1 )y 1 -λ y2 f 2 (y 2 )y 2 + λ 0 [f 1 (y 1 ) -f 2 (y 2 )] , (4)
where,

Φ(λ y1 , λ y2 , w) = λ y1 g 1 (w) + λ y2 g 2 (w). (5)
It is classical to consider that λ 0 = 1 since OCP (2) is stated as a maximization problem. The PMP implies the existence of absolutely-continuous λ yi , which are the co-state of y i , satisfying, λyi = -∂H/∂y i , for i = 1, 2. In addition, since the final states y i (T f ) are free, it follows that the transversality conditions are given by,

λ yi (T f ) = 0, for i = 1, 2. (6)
Finally, the PMP implies that the control function w ∈ W satisfies, for almost all t ∈ [0, T f ], max

w(t)∈[0,wm] H(y i (t), λ yi (t), w(t)). (7)
Inria ) provide the generic form of the functions f i and g i for i = 1, 2 (with arbitrary units), studied in this paper. (c) illustrates the resulting f i (δ -1 i (w)), in their typical general form studied throughout this work, where both situations

f 1 (δ -1 1 (w)) > f 2 (δ -1 2 (w)) and f 1 (δ -1 1 (w)) < f 2 (δ -1 2 (w)) hold over [0, w m ].
The resulting function ∆ in (d) satisfies the hypothesis w ∈ (0, w m ).

Introducing the static-OCP

The static-OCP (see, e.g., [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]) associated with the OCP (2), for the system at steady state is:

Maximize

w∈[0,wm] ∆(w) = f 1 (y 1 (w)) -f 2 (y 2 (w)), under the constraints f i (y i )y i = g i (w), i = 1, 2. (8)
In other words, the static-OCP seeks to maximize the function ∆ defined in (3) when system (1) is at steady state. Under the assumption that w ∈ (0, w m ) is the unique point maximizing ∆, as in Fig. 1, the definition of the static-OCP is extended to include the co-state variables using the Hamiltonian equation. Indeed, according to the Lagrange multipliers rule ( [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]), there exists, (λ y1 , λ y2 ) ∈ R 2 , s.t., the solution (y 1 , y 2 , w) of the static-OCP (8) implies that w, y 1 , y 2 , λ y1 , λ y2 also satisfies:

         ∂H ∂λy i (y 1 , y 2 , λ y1 , λ y2 , -1, w) = 0, -∂H ∂yi (y 1 , y 2 , λ y1 , λ y2 , -1, w) = 0, ∂H ∂w (y 1 , y 2 , λ y1 , λ y2 , -1, w) = 0. (9)
More precisely, using H given in (4), we deduce that the steady-state satisfying the static-OCP

(8) is solution of,              -f i (y i )y i + g i (w) = 0, i = 1, 2, λ y1 (f 1 (y 1 ) + f 1 (y 1 )y 1 ) -f 1 (y 1 ) = 0, λ y2 (f 2 (y 2 ) + f 2 (y 2 )y 2 ) + f 2 (y 2 ) = 0, g 1 (w)λ y1 + g 2 (w)λ y2 = 0, (10)
In the general case illustrated in Fig. 1, we define the unique static optimal steady-state that satisfies [START_REF] Bonnefond | Continuous selection pressure to improve temperature acclimation of Tisochrysis lutea[END_REF], denoted by,

E static = w, y 1 , y 2 , λ y1 , λ y2 . (11)
In Sect. 5, an explicit determination of E static is given for typical f i and g i functions in quota models (in the Droop's case).

3 Turnpike main result : Theorem 1

Let us consider the perturbed trajectories around the static optimal steady-state E static defined in [START_REF] Bonnans | BOCOP: an open source toolbox for optimal control -A collection of examples[END_REF]. More precisely, we perform the classical change of coordinates, δy i (t) = y i (t) -y i , δλ yi (t) = λ yi (t) -λ yi , for i = 1, 2, and we consider the shifted control, δw(t) = w(t) -w. Thus, as in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF], we can check that, δw(t) = -H -1 ww (H wy δy(t) + H wλ δλ y (t)), where, δy = (δy 1 , δy 2 ) T and δλ y = (δλ y1 , δλ y2 ) T . In what follows, we set z = (δy, δλ y ) T .

Using the Hamiltonian (4), we notice that,

H yw = ∂ 2 H ∂y∂w E static = (0, 0),
and thus around the unique point E static we end up with,

δ ẏ(t) δ λy (t) = Ψ δy(t) δλ y (t) + Λ 1 (δy(t), δλ y (t)) Λ 2 (δy(t), δλ y (t)) , (12)
for all t ∈ [0, T f ], where Λ i satisfies,

Λ i (δy, δλ y ) (δy, δλ y ) ---------→ (δy,δλy) → 0 0. (13)
and,

Ψ = H λy -H λw H -1 ww H wλ -H yy -H yλ , (14)
where,

H ww = ∂ 2 H ∂w 2 E static = λ y1 g 1 (w) + λ y2 g 2 (w), Inria H yy = diag{ -λ y1 (2f 1 (y 1 ) + y 1 f 1 (y 1 )) + f 1 (y 1 ), -λ y2 (2f 2 (y 2 ) + y 2 f 2 (y 2 )) -f 2 (y 2 )}, H λy = diag {-f 1 (y 1 ) -y 1 f 1 (y 1 ), -f 2 (y 2 ) -y 2 f 2 (y 2 )} , H wλ = (g 1 (w), g 2 (w)),
and thus we end up with,

H λw H -1 ww H wλ =    g 2 1 (w) λy 1 g 1 (w)+λy 2 g 2 (w) g 1 (w)g 2 (w) λy 1 g 1 (w)+λy 2 g 2 (w) g 1 (w)g 2 (w) λy 1 g 1 (w)+λy 2 g 2 (w) g 2 2 (w) λy 1 g 1 (w)+λy 2 g 2 (w)    .
For notational convenience, we set the following matrices and constants,

-H λw H -1 ww H wλ H yy = b 1 b 2 b 3 b 4 , b 23 = b 2 b 3 , a i = - ∂(f (y i )y i ) ∂y i E static , and
H λy = diag{a 1 , a 2 }.
Consequently, we can show that the eigenvalues κ j , j = 1, . . . , 4, of the matrix Ψ defined in ( 14) are given by,

κ 1 = 1 √ 2 a 2 1 + a 2 2 -b 1 -b 4 + ς 1/2 , κ 2 = -κ 1 , κ 3 = 1 √ 2 a 2 1 + a 2 2 -b 1 -b 4 -ς 1/2 , κ 4 = -κ 3 , (15)
where,

ς = (a 2 1 -a 2 2 ) -(b 1 -b 2 ) 2 -4b 23 .
Now, we are ready to state the main result of the paper.

Theorem 1. Let T f > 0 be a large fixed final-time and E static be the unique global solution of the static-OCP (8), s.t. w ∈ (0, w m ) maximizes ∆ in (3).

x The solution of OCP (2) with fixed initial conditions y i (0

) ∈ R * + , i = 1, 2, corresponds to the singular control function w s ∈ W, s.t., w s (t) = {0, w m } , for all t ∈ [0, T f ], (16)
i.e., the optimal control maximizing H in (4) is neither a bang 0 nor a bang w m .

y Moreover, if a 1 = a 2 , H ww < 0, and, Re(κ j ) = 0, for j = 1, . . . , 4,

where κ j are the eigenvalues of the matrix Ψ given in [START_REF] Chisti | Biodiesel from microalgae[END_REF], then there exist positive constants , c 1 , c 2 and a time-instant T * > 0, s.t., for all T f > T * , if,

2 i=1 y i (0) -y i + λ yi < , (18)
then, the optimal solutions of OCP (2) satisfy, for all t ∈ [0, T f ], the turnpike inequality,

2 i=1 y i (t) -y i + λ yi (t) -λ yi + w s (t) -w ≤ c 1 e -c2t + e -c2(T f -t) , (19)
where c 1 does not depend on T f and the optimal control is the singular w s in x.

Remark 1. Assuming that T f is sufficiently large is classical in turnpike framework ( [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]). The proof of Theorem 1 provides an estimate on the lower-bound of T f . This bound is satisfied in biological applications, e.g., metabolite production ( [START_REF] Caillau | The turnpike property in maximization of microbial metabolite production[END_REF][START_REF] Yegorov | Optimal control of bacterial growth for the maximization of metabolite production[END_REF]), algal growth ( [START_REF] Djema | Optimal control separating two microalgae species competing in a chemostat[END_REF][START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF]) due to the selection time which is much larger than the time-constants characterizing the system dynamics.

Remark 2. The conditions [START_REF] Djema | Optimization of Microalgae Selection: Highlighting turnpike features[END_REF] ensure that the real part of the eigenvalues of the matrix Ψ are non-zero (i.e. they are not on the imaginary axis), which leads to a hyperbolic behavior around the static steady-state. In practice, the conditions (17) are satisfied for a large class of systems describing microalgae growth, as illustrated in Sect. 5.5.

Remark 3. In the proof of the second part of Theorem 1, we use the approach given in [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF][START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF].

Notice that the main theorem from [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF] cannot be applied to our class of optimal control problems. Firstly, the control is bounded in the considered class of optimal control problem whereas it is unbounded in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]. Secondly, the matrix -H yy + H yw H -1 ww H wy is not positive definite in our case. For instance, choosing f i (y i ) = aiyi-bi yi , it follows that 2f 1 (y 1 ) + y 1 f 1 (y 1 ) = 0, and then we deduce that, -H yy + H yw H -1 ww H wy = -H yy = diag{f 1 (y 1 ), -f 2 (y 2 )}, which does not satisfy the positive definite condition in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]. This point is illustrated numerically in Sect. 5.5 using Droop's functions. Therefore, the hypotheses of the second part of Theorem 1 are slightly different from the hypotheses of the main theorem in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF].

The proof of Theorem 1 is divided in two parts. The first step concerns the statement x and derives from the PMP. The second step is based on the hyperbolicity of the linearized system ż = Ψz + Λ(z), where z = (δy, δλy) T . We end up the proof of point y using hyperbolic system properties, similarly to the arguments used in [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF][START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF], and which are adapted to the case of Theorem 1 as detailed in Appendix D.

Proof of Theorem 1

x Using (4)-( 5), we note that the PMP implies that the admissible control function satisfies, for almost all t ∈ [0, T f ], the maximization condition, max

w(t)∈[0,wm] H(y i , λ yi , w) = max w(t)∈[0,wm] Φ(λ yi , w). ( 20
)
At first sight, we deduce from ( 20) that any extremal solution should be a concatenation of bang and singular arcs over [0, T f ]. Indeed, using [START_REF] Bayen | Optimization of strain selection in evolutionary continuous culture[END_REF] we notice that in the case where λ y1 (t) ≥ 0 and λ y2 (t) ≥ 0, the control satisfying [START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF] is the bang w(t) = w m . Similarly, if λ y1 (t) ≤ 0 and λ y2 (t) ≤ 0, then it follows that the control satisfying [START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF] is the bang w(t) = 0. And finally, when λ y1 (t)λ y2 (t) < 0 at any t ∈ [0, T f ], we deduce from the fact that g i is a positive increasing function that the control satisfying [START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF] is the singular control denoted by w s (t), which maximizes Φ and which is neither 0 nor w m . However, using the co-states dynamics given by,

λy1 = λ y1 (f 1 (y 1 ) + f 1 (y 1 )y 1 ) -f 1 (y 1 ), λy2 = λ y2 (f 2 (y 2 ) + f 2 (y 2 )y 2 ) + f 2 (y 2 ), (21)
Inria where f i are continuously increasing, we deduce that,

λ y1 = f 1 (y 1 ) f 1 (y 1 ) + f 1 (y 1 )y 1 > 0,
and

λ y2 = -f 2 (y 2 ) f 2 (y 2 ) + f 2 (y 2 )y 2 < 0.
Then, since λ yi is an equilibrium point of ( 21) and using the transversality conditions ( 6), we deduce that,

0 < λ y1 (t) < λ y1 , λ y2 < λ y2 (t) < 0, (22)
for all t ∈ [0, T f ), where the latter inequalities derive from the fact that ∂(f i (y i )y i )/∂y i > 0 and f i (y i ) > 0 in [START_REF] Driver | Potential of bioenergy production from microalgae[END_REF]. In particular, [START_REF] Dochain | On-line estimation and adaptive control of bioreactors[END_REF] imply that for all t ∈ [0, T f ),

λ y1 (t)λ y2 (t) < 0. (23)
Consequently, we conclude from ( 23) that the control w s (t) is the one that satisfies the PMP and thus steers over the time window [0, T f ] the trajectories from y 1 (0), y 2 (0) to y 1 (T f ), y 2 (T f ). This concludes the proof of x. y Now, let us focus on the eigenvalues of the matrix Ψ defined in [START_REF] Cesari | Optimization -Theory and applications[END_REF]. Firstly, let us observe that,

det(Ψ -κI 4 ) = det H λy -κI 2 -H λw H -1 ww H wλ -H yy -H yλ -κI 2 ,
where I p is the p × p identity matrix. We notice that H yy and -H yλ -κI 2 are diagonal and thus they are commutative. Therefore, it follows that,

det(Ψ -κI 4 ) = det (H λy -κI 2 )(-H λy -κI 2 ) -H λw H -1 ww H wλ H yy .
Consequently, det(Ψ -κI 4 ) = 0, leads to,

(a 2 1 -κ 2 -b 1 )(a 2 2 -κ 2 -b 4 ) -b 23 = 0, (24)
where the constants a's and b's are defined before Theorem 1. From the equation [START_REF] Droop | Some thoughts on nutrient limitation in algae[END_REF] we readily get the eigenvalues κ j , j = 1, . . . , 4, given in [START_REF] Chisti | Biodiesel from microalgae[END_REF]. Under the suitable assumptions expressed in [START_REF] Djema | Optimization of Microalgae Selection: Highlighting turnpike features[END_REF], we ensure that the studied optimality system is hyperbolic, i.e., all eigenvalues of Ψ have nonzero real parts. Next, by Theorem 6 of [START_REF] Molinari | The time-invariant linear-quadratic optimal control problem[END_REF], it follows that there exist symmetric 2 × 2 matrices P -, P + maximal and minimal solutions of the Riccati equation,

H yλ P + P H λy -H yy -P H λw H -1 ww H wλ P = 0. (25)
Moreover, the corresponding matrices

Ω + = H λy -H λw H -1 ww H wλ P + and Ω -= H λy - H λw H -1
ww H wλ P -are such that eigenvalues of Ω + have strictly positive real part and eigenvalues of Ω -have strictly negative real part, i.e., Ω + is unstable and Ω -is stable. We denote by,

V := I I P -P + ( 26
)
where I is the 2 × 2 identity matrix. Straightforward calculations show that,

Ψ = V -1 ΨV = diag{Ω -, Ω + }.

RR n°9399

We denote by z the new variable s.t. V z = z, and then it follows from (12) that,

ż1 (t) = Ω -z1 (t) + θ 1 (z 1 (t), z2 (t)) , ż2 (t) = Ω + z2 (t) + θ 2 (z 1 (t), z2 (t)) , (27)
where the θ i , i = 1, 2, satisfy the o-little property as in [START_REF] Caperon | Nitrogen-limited growth of marine phytoplankton -II. Uptake kinetics and their role in nutrient limited growth of phytoplankton[END_REF]. At this juncture, we point out that under different hypothesis, we have rewritten the optimality conditions in the form [START_REF] Gugat | On the turnpike phenomenon for optimal boundary control problems with hyperbolic systems[END_REF] similarly to [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF] and [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF]. It can be shown (see Lemmas 2-3 in Appendix D), that there exists r > 0 s.t. for all (z 0 1 , zf 2 ) satisfying, z0

1 = z1 (0), zf 1 = z2 (T f ), (z 0 1 , zf 2 ) ∈ B(0, r), the corresponding solution (z 1 (t, (z 0 1 , zf 2 ), z2 (t, (z 0 1 , zf 2 )) of system (27) admits for any t ∈ [0, T f ] the bound,      z1 (t) ≤ C κ z0 1 e -κt + e -κ(T f -t) zf 2 Θ 1 ( z1 c 0 ) , z2 (t) ≤ C κ zf 2 e -κ(T f -t) + e -κt z0 1 Θ 2 ( z2 c 0 ) . ( 28
)
Next step is to verify that there exist T * > 0 and > 0 in (18) s.t., for any T f > T * and y(0) satisfying ( 18), (z 0 1 , zf 2 ) satisfy,

z0 1 + z2 (0) = y 0 -y, P -z1 (T f ) + P + zf 2 = -λ y . (29)
For that, we use the fixed point argument as in [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF], by defining the map

F : B(0, r) → R 4 ,s.t., F z0
1 zf 2 = y 0 -y -z2 (0) -P -1 + (P -z1 (T f ) + λ y )
, and showing that F admits the unique fixed point (z 0 1 , zf 2 ). This step is detailed Appendix D.4. It is worth mentioning that at this stage of the proof, we also get the lower-bound T * of T f and from [START_REF] Djema | An Optimal Control Strategy Separating Two Species of Microalgae in Photobioreactors[END_REF]. Finally, using (28) for T f > T * , we obtain for all t ∈ [0, T f ],

δy(t) = z1 (t, (z 0 1 , zf 2 )) + z2 (t, (z 0 1 , zf 2 )) ≤ C κ z0 1 e -κt + e -κ(T f -t) zf 2 Θ 1 (r) + C κ zf 2 e -κ(T f -t) + e -κt z0 1 Θ 2 (r) .
Up to defining a bigger constant C κ , we deduce that,

δy(t) ≤ C κ e -κt + e -κ(T f -t) , ∀t ∈ [0, T f ]. (30)
Note that C κ is independent of T f . Similarly, we get,

δλ y (t) ≤ C κ e -κt + e -κ(T f -t) , ∀t ∈ [0, T f ]. (31)
Note that there exists R = R(r) s.t., δy + δλ y ≤ R.

Then, using additional calculations for δw(t), we get,

δw(t) = w(y, λ y ) -w ≤ sup B((ȳ, λy),R) ∇ y w δy + sup B((ȳ, λy),R) ∇ λy w δλ y , (32)
for all t ∈ [0, T f ]. Combining (30), ( 31) and (32) leads us to the exponential turnpike estimate given in [START_REF] Djema | Optimal control separating two microalgae species competing in a chemostat[END_REF] and thus concludes the proof of Theorem 1.

Inria 5 Application to a higher dimensional dynamical system

Algal competition using Droop's model

In this section, we focus on a model of competition between two microalgae populations in dimension 5, with the aim of achieving species separation over a fixed time window. We tweak this issue using a system reduction in order to rewrite the optimization objective in term of the generic OCP [START_REF] Aravantinou | Selection of microalgae for wastewater treatment and potential lipids production[END_REF]. We stick to conventional biological notations following the classical Droop's model in order to remain consistent with earlier works (e.g., [START_REF] Djema | Optimal control separating two microalgae species competing in a chemostat[END_REF][START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF]). Let x 1 and x 2 be the biomass concentrations of two different strains of microalgae, initially coexisting in a bioreactor with one limiting substrate s, s.t., s(t) ∈ [0, s in ], for all t ≥ 0, where s in is the constant substrate influent concentration. Each microalgae species uptakes the free nutrient s in the bioreactor. The intracellular amount of this limiting element is the cell quota q i , with i = 1, 2.

A general class of quota models, including the classical Droop's model ( [START_REF] Droop | Vitamin B 12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri[END_REF]), writes,

     ṡ = (s in -s)D - 2 i=1 ρ i (s)x i , qi = ρ i (s) -µ i (q i )q i , ẋi = [µ i (q i ) -D] x i . ( 33
)
The function µ i in [START_REF] Norberg | Biodiversity and ecosystem functioning: a complex adaptive systems approach[END_REF] stands for the growth rate of the i-th species, while ρ i is its nutrient uptake rate of free nutrient. The dilution rate D and the substrate input s in can generally be manipulated in the bioreactor. The functions ρ i and µ i are nonnegative and increasing bounded functions, s.t., 0 ≤ ρ i (s) ≤ ρ mi , and, 0 ≤ µ i (q i ) ≤ µ mi , where ρ mi and µ mi are positive constants. Typically in Droop's model, the uptake rate ρ i (s) is expressed in terms of Michaelis-Menten kinetics,

ρ i (s) = ρ mi s K si + s , (34)
where K si is a positive constant. The most classical form -known as Droop's growth-for µ i is given for i = 1, 2 by,

µ i (q i ) = µ i∞ 1 - Q i min q i , q i ≥ Q i min . (35)
The minimum intracellular quota Q i min > 0 defines the quota under which cell division cannot occur (thus, µ i (q i ) = 0 for all q i ≤ Q i min ). Other forms have been used to describe the growth rate such as the Caperon-Meyer equation ( [START_REF] Caperon | Nitrogen-limited growth of marine phytoplankton -II. Uptake kinetics and their role in nutrient limited growth of phytoplankton[END_REF]),

µ i (q i ) = µ i∞ (q i -Q i min ) q i -Q i min + k qi , q i ≥ Q i min , k iq > 0, (36)
another possible form that we also use in this work is,

µ i (q i ) = µ i∞ q ni i q ni i + k ni qi , n i > 1. (37)
Notice that the choices (36)-(37) lead to nonlinear terms µ i (q i )q i in [START_REF] Norberg | Biodiversity and ecosystem functioning: a complex adaptive systems approach[END_REF]. In all these cases, the resulting quota is bounded for all t ≥ 0, i.e., Q i min ≤ q i (t) ≤ q mi , where q mi is the maximum internal storage rate ( [START_REF] Bernard | Hurdles and challenges for modelling and control of microalgae for CO 2 mitigation and biofuel production[END_REF]). In the sequel, the quota variables q i will stand for y i in model ( 1) and we will prove that:

(i) OCP (2) achieves microalgae separation, through maximization of relative concentration

x 1 (T f )/x 2 (T f ),
(ii) Theorem 1 applies for different functions µ i (given by ( 35)-( 37)), i.e., turnpike-type behaviors hold and characterize microalgae selection, under various types of growth rates.

Controlling the substrate concentration

Now, let us assume that the free substrate concentration s is controlled by an ad'hoc controller.For instance, by playing with the concentration of the influent substrate s in and with the dilution rate D, the concentration s can follow a desired (smooth enough) reference trajectory s(t) = s * (t).

Different control designs can be used such as linearizing control ( [START_REF] Dochain | On-line estimation and adaptive control of bioreactors[END_REF]), asymptotic control ( [START_REF] Mailleret | Nonlinear adaptive control for bioreactors with unknown kinetics[END_REF]) or other approaches that demonstrated their efficiency on microalgae (see, e.g., [START_REF] Bernard | Modelling of microalgae culture systems with applications to control and optimization[END_REF][START_REF] Tebbani | CO2 Biofixation by Microalgae: Modeling, estimation and control[END_REF]). By assuming a perfect tracking of the substrate towards w(t), the system (33) then reduces to,

qi = ρ i (w) -µ i (q i )q i , ẋi = [µ i (q i ) -D] x i , ( 38
)
This system is triangular, and the dynamics of x i results from the dynamics of q i . Consequently, we can focus exclusively on the quota-dynamics q i , for i = 1, 2.

Maximizing the relative biomass density

Species separation is achieved when maximizing the relative density of the species of interest x 1 at T f , i.e., finding the optimal control maximizing the quantity,

C(T f ) = x 1 (T f ) x 2 (T f ) . ( 39
)
Let us state the following key result to end up with a functional that does not depend on the states x i , but still maximizes the relative concentration at the final-time.

Proposition 1. Any solution of the class of systems [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF] maximising

C [0,T f ] defined in OCP (2)
is equivalently maximizing CT f in [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF], where the functions f i in OCP (2) stand for the growth rates µ i in system [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF].

Proof: Indeed, the cost function in OCP (2) writes,

C [0,T f ] = T f 0 [µ 1 (q 1 (t)) -D] -[µ 2 (q 2 (t)) -D] dt. ( 40
)
Using the state dynamics, µ i (q i ) -D = ẋi xi , we get,

C [0,T f ] = [ln(x 1 (t)) -ln(x 2 (t))] T f 0 = ln x 1 (T f )x 2 (0) x 1 (0)x 2 (T f ) . ( 41
)
Thus, maximizing C [0,T f ] indeed consists in maximizing CT f . The advantage of considering the criterion C [0,T f ] instead of the natural criterion CT f is that it is independent of the dynamics of the states x i , for i = 1, 2.

Thus, the optimization problem can now be reduced to an OCP as in (2), involving the cost C [0,T f ] and associated with the reduced-system (in the form (1)),

q1 = -µ 1 (q 1 )q 1 + ρ 1 (w), q2 = -µ 2 (q 2 )q 2 + ρ 2 (w), (42)
Inria where the growth rate µ i is given by one of the functions ( 35)-( 37), the uptake rate ρ i is given by [START_REF] Rodolfi | Microalgae for oil: Strain selection, induction of lipid synthesis and outdoor mass cultivation in a low-cost photobioreactor[END_REF], and the control w is equivalent to the substrate s in [START_REF] Norberg | Biodiversity and ecosystem functioning: a complex adaptive systems approach[END_REF] once regulated by the ad'hoc controller. This is especially true since x in Theorem 1 states that the solution of the OCP (2) corresponds to the singular control function with no bangs over the time-window [0, T f ], and consequently the optimal control w # (t) has thus the required smoothness to be tracked in the original system [START_REF] Norberg | Biodiversity and ecosystem functioning: a complex adaptive systems approach[END_REF] through the ad'hoc controller. To sum up, we can see that finding a control for system [START_REF] Norberg | Biodiversity and ecosystem functioning: a complex adaptive systems approach[END_REF] approaching as close as possible the solution w # (t) of the OCP studied in Sect. 2-3 will provide an excellent suboptimal control for the OCP of the ful-system (33) associated with the criterion (39).

Optimal synthesis of species separation

Now, we are ready to apply the results obtained in Sect. 2-3 to the problem of species selection. Firstly, the static-OCP associated with [START_REF] Trélat | Integral and measure-turnpike properties for infinite-dimensional optimal control systems[END_REF] and

C [0,T f ] states, Maximize w∈[0,wm] ∆(w) = µ 1 (q 1 (w)) -µ 2 (q 2 (w)),
where

-µ i (q i )q i + ρ i (w) = 0 for i = 1, 2 , ( 43
)
under the steady-state constraints -µ i (q i )q i + ρ i (w) = 0 for i = 1, 2. Similarly to Section 2, the solution of ( 43) -including the co-states-is denoted by, E static = q 1 , q 2 , λ q1 , λ q2 , w . This optimal steady-state satisfies ( 9)-( 10), using this time the quota-coordinates in [START_REF] Trélat | Integral and measure-turnpike properties for infinite-dimensional optimal control systems[END_REF], under the standard assumption of Section 2, i.e., w ∈ (0, w m ) as in Fig. 1, which is generally satisfied for ρ i and µ i in ( 34)-( 37), as illustrated later in numerical simulations. A case in point is the classical Droop's model associated with ( 34)-( 35), in which case the function ∆ rewrites,

∆(w) = µ 1 (δ -1 1 (w)) -µ 2 (δ -1 2 (w)), (44)
where, δ i ( ) = ρ -1 i (µ i ( ) ), and since, ρ -1 i ( ) = Ki ρim-, for all ∈ [0, ρ im ), we obtain,

δ i ( ) = K i -Qi min σi-
, where, σ i = ρim µi∞ + Q i min . This leads to,

δ -1 i ( ) = σ i + K i Q i min + K i ,
and consequently we end up in the Droop's case with,

µ i (δ -1 i ( )) = ρ im σ i + K i Q i min , i = 1, 2. (45)
Using [START_REF] Wiechert | Modeling and simulation: tools for metabolic engineering[END_REF] we can compute successive derivatives of ∆(w) in the Droop's case, i.e.,

∆ (w) = ρ 1m K 1 Q 1 min (K 1 Q 1 min + σ 1 w) 2 - ρ 2m K 2 Q 2 min (K 2 Q 2 min + σ 2 w) 2 ,
and, ∆ (w) = κ 1 (w)κ 2 (w), where,

κ i (w) = - 2K i µ 2 i∞ Q i min ρ im (µ i∞ Q i min + ρ im ) (K i µ i∞ Q i min + (µ i∞ Q i min + ρ im )w) 3 .
Thus, by looking for the points that satisfy ∆ (w) = 0 and ∆ (w) < 0 in the case where µ i is given by [START_REF] Odjadjare | Potential biotechnological application of microalgae: a critical review[END_REF], we deduce that the steady-state E static = w, q 1 , q 2 , λ q1 , λ q2 , satisfies,

q i = 1 µ i∞ (ρ i (w) + µ i∞ Q i min ) , for i = 1, 2, RR n°9399 λ q1 = - Q 1 min q 2 1 ,
and,

λ q2 = Q 2 min q 2 2 .
Numerically, we check that only one numerical value w belongs to (0, w m ) for a sufficiently large w m , leading then to the configuration in Fig. 1 (Fig. 4 in Droop's case). We refer to Appendices A-B for more details on the practical determination of E static in the Droop's case.

In a similar way, we can determine the analytic expression of E static for different growth-rates (( 36)-( 37)).

Therefore, under suitable conditions a unique E static exists and satisfies w ∈ (0, w m ), and thus in light of Theorem 1 we derive the following statement as a direct consequence.

Claim 1 (Application of Theorem 1). Let fix any biological parameters for the functions ρ i in (34) and µ i in (35)-( 37), i = 1, 2, in system (33), s.t., there exists a unique E static with w ∈ (0, w m ), solution of the static-OCP (43) then, x the optimal control w # (t) for all t ∈ [0, T f ] is singular, i.e., it does not involve bang arcs 0 and w m over [0, T f ]. y Moreover, if the corresponding conditions (17) and (18) are fulfilled, the singular optimal control (substrate concentration) w # (t) and its related quota q i (t), as well as their respective costates λ qi (t), i = 1, 2, exhibit over the fixed time-horizon [0, T f ] an exponential turnpike behavior characterized in [START_REF] Djema | Optimal control separating two microalgae species competing in a chemostat[END_REF].

Note that the overlapping of favorable features (namely the facts that w solution of the static-OCP is remote from the bounds 0 and w m , the optimal control is inherently singular, i.e., w s (t) ∈ (0, w m ) for all t ∈ [0, T f ], and also the exponential turnpike keeping w s (t) almost constant most of the time) results in strong smoothness of the optimal control and the optimal trajectories, as illustrated in the next section. Nonetheless, this smoothness does not imply that the resulting ad'hoc controller of system (33) should be smooth as well, as shown in [START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF] where the obtained dilution-based control is of type bang-singular.

To sum up the methodology provided in this paper: the first step is to determine the optimal point E static s.t. w ∈ (0, w m ), solution of the static-OCP corresponding to a given dataset of biological parameters and functions in the system [START_REF] Norberg | Biodiversity and ecosystem functioning: a complex adaptive systems approach[END_REF]. Next, we compute the matrix Ψ defined around E static as in [START_REF] Cesari | Optimization -Theory and applications[END_REF]. Then, the conditions stated in y-Theorem 1 ensure the existence of a turnpike-type behavior of the optimal solution of the reduced system (1). The numerical computation of the eigenvalues κ j , j = 1, . . . , 4, which are defined in [START_REF] Chisti | Biodiesel from microalgae[END_REF], is illustrated in different cases of the functions µ i in Sect. 5.5.

Finally, it is worth noting that the conditions [START_REF] Djema | Optimization of Microalgae Selection: Highlighting turnpike features[END_REF] in Theorem 1 appear to be systematically satisfied in numerical simulations, leading accordingly to a turnpike behavior in all the evaluated cases.

Numerical synthesis through direct optimization methods

The direct-optimization approach performed in this section illustrates the results of Theorem 1 applied to the quota-model in Sect. 5 (System (42)) for the different growth functions (( 35)-( 37)) illustrated in Fig. 2.

In this perspective, the state variables and the control function of system (42) are discretized, by setting a subdivision of the time-horizon [0, T f ] with fixed T f , and the substrate-control w is discretized over W, s.t. w is piecewise constant on the time-interval subdivision, with values Inria belonging to [0, w m ]. The state dynamics q 1 and q 2 (42) is discretized with Lobatto scheme ( [START_REF] Betts | Practical methods for optimal control and estimation using nonlinear programming[END_REF]), i.e. RK method of type Lobatto-IIIC, based on an implicit trapezoidal rule, of order 6. suitable for stiff-nonlinear problems, through a scheme of ordinary differential equations (ODE) decomposition (e.g., implicit/explicit Euler schemes, Runge-Kutta methods, etc.). Thus, numerical direct methods transform OCP (2) (stated in Sect. 2, associated with System (42)) through the discretization step, into a nonlinear programming problem (NLP) in finite-dimension of the form, Maximize {w1,...,wp}∈[0,wm] C [0,T f ] (w 1 , . . . , w p ),

q 1,2 k+1 = f i (t k , q 1,2 k , w k ), k = 0, . . . , N -1. (46)
When applied to an ODE of the form, ż = f(t, z), Lobatto Runge-Kutta methods are expressed as,

Z ki = z k + h k s j=1 a ij f(t k + c j h k , Z kj ), with, i = 1, . . . , s, z k+1 = z k + h k s j=1 b j f(t k + c j h k , Z kj ), (47)
where the stage value s defines the order of the Lobatto's method (the order is given by 2s + 2 in Lobatto-III family), and the coefficients a ij , b j and c j are given according to s, verifying a set of conditions (see, e.g., [START_REF] Hairer | Solving ordinary differential equations[END_REF], Chap. VI.5). In this case, the stage s is 4 and the coefficients a ij , b j and c j in [START_REF] Williams | Adaptation and natural selection: A critique of some current evolutionary thought[END_REF] are given in Tab. 1.

Using the numerical values of the biological parameters in Fig. 2, we performed several direct optimizations in Bocop, based on the settings given in Tab. 2. with appropriate settings (tolerance for NLP solver: 10 -34 , 2000 time-steps, see Tab. 2). No issue is encountered regarding the initialization and starting points for solving the studied problem. In all the following examples, the conditions of Theorem 1 are satisfied, and the optimal solutions satisfy the turnpike property [START_REF] Djema | Optimal control separating two microalgae species competing in a chemostat[END_REF] as stated in Claim 1. [START_REF] Caperon | Nitrogen-limited growth of marine phytoplankton -II. Uptake kinetics and their role in nutrient limited growth of phytoplankton[END_REF][START_REF] Droop | Vitamin B 12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri[END_REF]). From left to right: absorption functions ρ i in [START_REF] Rodolfi | Microalgae for oil: Strain selection, induction of lipid synthesis and outdoor mass cultivation in a low-cost photobioreactor[END_REF] with, ρ 1m = 9.1 (µmol.L -1 .d -1 ), K 1 = 0.08 (µmol.L -1 ) ρ 2m = 9.6(µmol.L -1 .d -1 ), K 2 = 0.25(µmol.L -1 ). Growth rate µ i in [START_REF] Odjadjare | Potential biotechnological application of microalgae: a critical review[END_REF] with, µ 1∞ = 1.6day -1 ,

Q 1 min = 0.95 µmol/L, µ 2∞ = 2.2 day -1 , Q 2 min = 1.5 µmol/L. Growth rate µ i in (36) with, µ 1∞ = 0.75day -1 , Q 1 min = 1.6µmol/L, k q1 = 0.9µmol/L, µ 2∞ = 0.9day -1 , Q 2 min = 2µmol/L, k q2 = 0.6 µmol/L. Growth rate µ i in (37) with, µ 1∞ = 1.5 day -1 , n 1 = 4, k q1 = 3 µmol/L, µ 2∞ = 1.15 day -1 , n 2 = 3, k q2 = 2 µmol/L.
Example 1. We set T f = 15 days and we consider that the growth functions µ i are in the Droop's form [START_REF] Odjadjare | Potential biotechnological application of microalgae: a critical review[END_REF]. In this case, the sufficient conditions given in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF] are not satisfied (see Remark 3). The computations performed in Appendices A-B give the analytical expression of the singular control in this case.

The function

∆(w) = µ 1 (δ -1 1 (w)) -µ 2 (δ -1 2 (w)) is given in Fig. ( 3
). In addition, the static steady-state E static and the eigenvalues of the matrix Ψ are given in Table . Notice also that Figure The functions µ i (δ -1 i (w)) and ∆(w) using the data given in Fig. 2 with µ i as in [START_REF] Odjadjare | Potential biotechnological application of microalgae: a critical review[END_REF].

Example 2. We set T f = 30 days and we consider that the growth functions µ i are in the form [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF], with the biological parameters given in Fig. 2. Results are in Table 3. The turnpike holds according to Theorem 1, as described in Claim 1 and illustrated in Fig. 5. Notice that the sufficient conditions for the turnpike existence given in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF] In all these examples (Fig. 456), the optimal solutions of the OCP (2), settled in relatively large time [0, T f ], consist of three pieces: the first and the last ones are relatively short arcs, and the middle piece being a relatively long-time arc staying close to the optimal steady-state solution of the associated static-OCP (8).

Conclusion

The optimal control has been fully-characterized through a complete study combining Pontryagin's principle and turnpike frameworks relying upon classical hyperbolic system properties. The turnpike theory ( [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF][START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF]) considerably simplifies the understanding of the control strategy and provides a clear insight into the optimal control structure. By investigating the OCP and its associated static version, the main result in this paper proves that the optimal control is always singular and provides sufficient conditions ensuring the existence of turnpike behaviors characterizing the optimal trajectories. We had to upgrade the existing turnpike results so that they can apply to the specificity of our models. The considered class of systems naturally appears in the framework of competition between two microbial or algal species, but a similar structure is likely to play a role for other classes of problems. In general, solving an OCP in dimension 5 is a challenging issue (see [START_REF] Djema | Optimal control separating two microalgae species competing in a chemostat[END_REF][START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF]). Thus, the microalgae optimal selection problem was reduced in this paper considering the core dynamics in dimension 2. This reduction approach using a back-stepping philosophy might also prove to be efficient in other contexts (see, e.g., [START_REF] Yegorov | Optimal control of bacterial growth for the maximization of metabolite production[END_REF][START_REF] Caillau | The turnpike property in maximization of microbial metabolite production[END_REF]), and this will be the main focus of our future work.

B The analytic expression of the singular static-control, solution of the static-OCP (43)

The singular controls (candidate) obtained in the case of the growth function µ i in Droop's form [START_REF] Odjadjare | Potential biotechnological application of microalgae: a critical review[END_REF], are given by, To give more details on the proof of the bound (28) for solutions of [START_REF] Gugat | On the turnpike phenomenon for optimal boundary control problems with hyperbolic systems[END_REF] we can follow ( [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF][START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]) using classical theory on hyperbolic system. The following technical Lemmas are given for a general system [START_REF] Gugat | On the turnpike phenomenon for optimal boundary control problems with hyperbolic systems[END_REF] of any dimension n ≥ 2.

D.1 Exponential estimate in the stable case

Figure 1 :

 1 Figure1: (a-b) provide the generic form of the functions f i and g i for i = 1, 2 (with arbitrary units), studied in this paper. (c) illustrates the resulting f i (δ -1 i (w)), in their typical general form studied throughout this work, where both situationsf 1 (δ -1 1 (w)) > f 2 (δ -1 2 (w)) and f 1 (δ -1 1 (w)) < f 2 (δ -1 2 (w)) hold over [0, w m ].The resulting function ∆ in (d) satisfies the hypothesis w ∈ (0, w m ).

Figure 2 :

 2 Figure 2: Various functions used in Droop and Caperon-Meyer models ([START_REF] Caperon | Nitrogen-limited growth of marine phytoplankton -II. Uptake kinetics and their role in nutrient limited growth of phytoplankton[END_REF][START_REF] Droop | Vitamin B 12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri[END_REF]). From left to right: absorption functions ρ i in[START_REF] Rodolfi | Microalgae for oil: Strain selection, induction of lipid synthesis and outdoor mass cultivation in a low-cost photobioreactor[END_REF] with, ρ 1m = 9.1 (µmol.L -1 .d -1 ), K 1 = 0.08 (µmol.L -1 ) ρ 2m = 9.6(µmol.L -1 .d -1 ), K 2 = 0.25(µmol.L -1 ). Growth rate µ i in[START_REF] Odjadjare | Potential biotechnological application of microalgae: a critical review[END_REF] with, µ 1∞ = 1.6day -1 , Q 1 min = 0.95 µmol/L, µ 2∞ = 2.2 day -1 , Q 2 min = 1.5 µmol/L. Growth rate µ i in (36) with, µ 1∞ = 0.75day -1 , Q 1 min = 1.6µmol/L, k q1 = 0.9µmol/L, µ 2∞ = 0.9day -1 , Q 2 min = 2µmol/L, k q2 = 0.6 µmol/L. Growth rate µ i in (37) with, µ 1∞ = 1.5 day -1 , n 1 = 4, k q1 = 3 µmol/L, µ 2∞ = 1.15 day -1 , n 2 = 3, k q2 = 2 µmol/L.

Figure 4 :Example 3 .Figure 5 :

 435 Figure 4: The optimal control and optimal trajectories given by Bocop in Example 1.

Figure 6 :

 6 Figure 6: The optimal control and optimal trajectories given by Bocop in Example 3.

w * 1 , 2 = 2 .

 122 ± µ1∞µ2∞ (K2µ2∞Q2 min (µ1∞Q1 min + ρ1m) -K1µ1∞Q1 min (µ2∞Q2 min + ρ2m) K1K2Q1 minQ2 minρ1mρ2m+ K1K2µ1∞µ2∞Q1 minQ2 min (-µ2∞ρ1mρ2m + µ1∞ (µ2∞Q2 minρ1m -µ2∞Q1 minρ2m + ρ1mρ2m)) K2µ 2 2∞ Q2 min (µ1∞Q1 min + ρ1m) 2 ρ2m -K1µ 2 1∞ Q1 minρ1m (µ2∞Q2 min + ρ2m) Let us consider, F 1 : [a, b] → R + ∈ C 0 ([a, b], R + ), F 2 : [a, b] → (0, +∞) ∈ C 0 ([a, b], R + ), F 3 : [a, b] → R + ∈ C 0 ([a, b], R + ). Let, ∀y : [a, b] → R + ∈ C 0 ([a, b], R + ), s.t., y(t) ≤ F 1 (t) +

t a F 2 F 2

 22 (t)F 3 ( )y( )d , ∀t ∈ [a, b]. Then, we get for all t ∈ [a, b],y(t) ≤ F 2 (t)e t a F2( )F3( )d t a e -m a F2( )F3( )d F 1 (m)F 3 (m)dm + F 1 (t).(52)The backward version is given for all t ∈ [a, b] by,∀y : [a, b] → R + ∈ C 0 ([a, b], R + ), s.t., y(t) ≤ F 1 (t) + b t (t)F 3 ( )y( )d .Then, we get for all t ∈ [a, b],y(t) ≤ F 2 (t)e b t F2( )F3( )d b t e -b t F2( )F3( )d F 1 (m)F 3 (m)dm + F 1 (t).(53)D Technical detail for the proof of Theorem 1

Lemma 1 . 1 2 . 2 2 . 3 7 3 8 4 10 5 6 20 A The singular control expression 24 B 25 C Gronwall's Lemma 25 D 25 D. 1

 1122237810620242525251 ([36]). Let Ψ ∈ M n (R) be s.t. its spectrum satisfies Re (spec( Ψ)) ⊂ (-∞, 0), andlet θ ∈ C 1 (R n , R n ) satisfies θ( )/ ----→ → 00. Then, for any κ ∈ 0, -sup Re (spec( Ψ)) , there exist C κ > 0 and r κ > 0, s.t., for anyT f > 0, if z(t) satisfies, ż = Ψz + θ(z), z ≤ rκ ,(54)for all t ∈ [0, T f ], then it follows that,z(t) ≤ C κ e -κt z(0) , ∀t ∈ [0, T f ] .Statement of the general OCP . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Pontryagin's maximum principle (PMP) . . . . . . . . . . . . . . . . . . . 6 Introducing the static-OCP . . . . . . . . . . . . . . . . . . . . . . . . . . . . Turnpike main result : Theorem 1 Proof of Theorem 1 Application to a higher dimensional dynamical system 13 5.1 Algal competition using Droop's model . . . . . . . . . . . . . . . . . . . . 13 5.2 Controlling the substrate concentration . . . . . . . . . . . . . . . . . . . . 14 5.3 Maximizing the relative biomass density . . . . . . . . . . . . . . . . . . . 14 5.4 Optimal synthesis of species separation . . . . . . . . . . . . . . . . . . . . 15 5.5 Numerical synthesis through direct optimization methods . . . . . . . . 16 Conclusion The analytic expression of the singular static-control, solution of the static-OCP (43) Technical detail for the proof of Theorem 1 Exponential estimate in the stable case . . . . . . . . . . . . . . . . . . . . 25 D.2 Existence, uniqueness and continuity : classical results for hyperbolic systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 D.3 Exponential estimate of the hyperbolic system (proving (28)) . . . . . . 27 D.4 Projection of the terminal conditions (proving (29)) . . . . . . . . . . . . 28

Table 1 :

 1 The values of the coefficients corresponding to s = 4 in Lobatto IIIC.

	Final time (days)	T f = 15 (Example 1), T f = 20 (Example 2), T f = 50 (Example 3)
	Discretization	Lobatto IIIC (implicit, 4-stage, 6th order)
	Time steps	2000
	Max. iterations	2000
	NLP solver tolerance	10 -34

Table 2 :

 2 The numerical settings used in Bocop, associated with the biological parameters given in Fig.2.in this case H yy = diag{0.2997, -1.0443}. According to Theorem 1 (see Claim 1), the turnpike property (19) holds as shown in Fig.4.

	1.5
											0.4
	1
											0.2
	0.5
											0
	0
	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5

Table 3 :

 3 are not satisfied, since H yy = diag{0.116993, -0.939631}. Values of the static steady-state E static and eigenvalues of the matrix Ψ defined at E static for the 3 examples.

		Example 1	Example 2	Example 3
	w	0.021	0.015	0.074
	q 1	2.164	0.015	3.91
	q 2	1.84	0.015	2.69
	λ q1	0.44	0.322	0.173
	λ q2	-0.20	-0.1377	2.690
	κ 1,2 2.332 ±0.81i	1.74 ±1.09i	±3.02
	κ 3,4 -2.332 ±0.81i -1.74 ±1.09i	±0.82
				Inria

RR n°9399

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements: This work was supported by the ADEME BioMSA and Maximic ANR projects in France, and also benefited from the support of the ANR PhotoBiofilm Explorer (ANR-20-CE43-0008). The authors acknowledge J.-B. Caillau, J.-B. Pomet and J.-L. Gouzé for the stimulating discussions on this topic.

A The singular control expression

We choose g i (w) = ρ im w/(K i + w), with K 1 = K 2 . In the case where λ y1 λ y2 < 0, it is useful to perform the change of variable,

Thus, maximizing [START_REF] Djema | Separating Two Species of Microalgae in Photobioreactors in Minimal Time[END_REF] is equivalent to maximizing the function,

We can readily check that ∂ Φ/∂v = 0 is equivalent to,

The corresponding discriminant is then given by,

and therefore, we get the solutions,

the two solutions possibly belonging to (0, v m ) are,

Therefore we can readily determine the singular controls w s , w s 1 and w s 2 , bounded by 0 and w m , using the expressions ( 49), ( 50) and [START_REF] Yegorov | Optimal control of bacterial growth for the maximization of metabolite production[END_REF]. Consequently, in the case of classical absorption function in Droop's model, we deduce that the optimal controls satisfying (20) can be a concatenation of bangs (when λ y1 λ y2 > 0) and singular controls w s (when λ y1 λ y2 < 0), s.t.,

On the other hand, x of Theorem 1 shows that λ y1 λ y2 < 0, and therefore only singular controls are admissible.

Inria

Proof. The proof relies essentially on the fact that for any κ ∈ (0, -sup Re (spec(Ψ))) there exists a Cκ > 0 s.

D.2 Existence, uniqueness and continuity : classical results for hyperbolic systems

Now, let us fix

Let us consider the system,

Following ( [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF]), we define for later use the constant K > 0, associated with Ω -and Ω + , by,

Now we are ready to state and prove the following lemma.

Lemma 2. ( [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF]) For any T f > 0 and any r > 0, there exists ∈ (0, r) s.t. the following statements hold.

• For any

Proof. The first point follows from standard properties of affine systems. Its proof relies on Picard theorem (see, e.g., [START_REF] Zaidenberg | Picard's theorem and hyperbolicity[END_REF]), which implies that for any

Inria Thus, we are left to show that the map Ξ defined by Ξ(h 1 , h 2 ) = (z 1 , z2 ) admits a fixed point in the space (B C 0 (0, ), • ∞ ), where B C 0 (0, ) is the ball in the space C 0 ([0, T f ], R 2n ), • ∞ of radius . This can be shown with chosen in such way that,

The second statement of Lemma 2 can be shown using Ascoli-Arzela's theorem (see, e.g., [START_REF] Dieudonné | Foundations of modern analysis[END_REF]).

D.3 Exponential estimate of the hyperbolic system (proving (28))

Based on the previous lemmas, we give more details on the proof of the estimate [START_REF] Hairer | Solving ordinary differential equations[END_REF]. The following results hold for a general system [START_REF] Gugat | On the turnpike phenomenon for optimal boundary control problems with hyperbolic systems[END_REF] of any dimension, associated with a more general form of matrix V , that we denote here by, V = ν 1 ν 2 ν 3 ν 4 , where ν 1 and ν 2 are invertible. Now, let us prove the following statement.

Lemma 3. ([36]

) For any κ ∈ (0, min{-max{Re (spec (Ω -))}, min{Re (spec (Ω + ))}) there exists rκ ∈ (0, ∞) independent from T f ∈ (0, ∞) and functions

then, we get the following estimates for all t ∈ [0, T f ],

Proof. First, let us highlight that the statement in Lemma 3 holds for linear systems, i.e., in the particular case when θ 1 = θ 2 ≡ 0 in system (56). Indeed, in the latter case the resulting system consists of two independent systems and Lemma 1 can be applied in forward time for the system associated with Ω -and in reverse time for the system associated with Ω + . Next, when θ 1 and θ 2 are not identically 0, they admit the property Jac θ i (0) = 0, for i = 1, 2, where Jac θ i stands for the Jacobian of θ i , which is in our case continuous near 0. Therefore, for any > 0 there exists δ > 0 such that Jac θ i (m) < when m < δ. Then, using the Taylor expansions of θ i , i = 1, 2, around 0, for each θ i we can defined

0. More precisely, functions Γ i 1 , Γ i 2 are defined in the following way for i = 1, 2,

and,

Then, we can express θ i in terms of Γ i 1 and Γ i 2 as follows, r). These inequalities can be applied together with several successive applications of the Gromwall's lemma leading to the desired estimates (57) where Θ 1 = Γ 1 1 and Θ 2 = Γ 2 2 .

D.4 Projection of the terminal conditions (proving (29))

Here we provide calculations used to show [START_REF] Huertas | Estimating the capability of different phytoplankton groups to adapt to contamination: herbicides will affect phytoplankton species differently[END_REF] in proof of Theorem 1 (Section 4). The main idea is again inspired by [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF], by applying Lemmas 2-3 to ( 27), we have that there exists r κ > 0 s.t., for all (z 0 1 , zf 2 ) ∈ B(0, rκ K ), the corresponding (z 1 , z2 ) is bounded z1 (t) + z2 (t) ≤ r κ for all t ∈ [0, T f ] and thus, it satisfies [START_REF] Hairer | Solving ordinary differential equations[END_REF]. Let us define F : B(0

, where ν j appear in the matrix of V (generalized form of ( 26)). Now, let us show that there exists r satisfying rκ K > r > 0, T * > 0, and (r), s.t., if T f > T * and the initial and final conditions satisfy [START_REF] Djema | An Optimal Control Strategy Separating Two Species of Microalgae in Photobioreactors[END_REF], then F is contracting on B(0, r), that is, F : B(0, r) → B(0, r). By definition we have,

By applying (28), we get

We take into account (z 0 1 , zf 2 ) ∈ B(0, r) and ( 18) which leads to the further simplifications

Let us define T * according to, C κ ( ν 2 + ν 3 ) e -κT * = 1 4( ν1 -1 + ν4 -1 ) , and r > 0 for which there holds (from Lemma 3), Θ 2 (Kr) ≤ 1 4 ( ν 1 -1 + ν 4 -1 ) C κ ν 2 , and, Θ 1 (Kr) ≤ 1 4 ( ν 1 -1 + ν 4 -1 ) C κ ν 3 .

Finally, we take T f ≥ T * and satisfying ≤ r 4 ( ν 1 -1 + ν 4 -1 ) .

This implies F (z 0 1 , z f 2 ) ≤ r. and proves that F is contracting. By Brouwer's fixed-point theorem, F admits the unique fixed point (z 0 1 , zf 2 ). Now, it sufficient to notice that F (z 0 1 , zf 2 ) = (z 0 1 , zf 2 ) is equivalent to [START_REF] Huertas | Estimating the capability of different phytoplankton groups to adapt to contamination: herbicides will affect phytoplankton species differently[END_REF] in the case of the typical form of V .

Inria