# Optimal Darwinian Selection of Microorganisms with Internal Storage

Walid Djema, Térence Bayen, Olivier Bernard

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

## To cite this version:

Walid Djema, Térence Bayen, Olivier Bernard. Optimal Darwinian Selection of Microorganisms with Internal Storage. Processes, 2022, 10 (3), pp.461. 10.3390/pr10030461. hal-03699709

HAL Id: hal-03699709

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

Submitted on 20 Jun 2022

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Optimal Darwinian Selection of Microorganisms with Internal Storage

Walid Djema<fnref n="1" />*, Térence Bayen<fnref n="2" /> and Olivier Bernard <fnref n="1" />,<fnref n="3" />

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

check for
updates

**Citation:** Djema, W.; Bayen, T.; Bernard, O. Optimal Darwinian Selection of Microorganisms with Internal Storage. *Processes* **2022**, *10*, 461. https://doi.org/10.3390/pr10030461

Academic Editors: Philippe Bogaerts and Alain Vande Wouwer

Received: 19 January 2022

Accepted: 14 February 2022

Published: 24 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** In this paper, we investigate the problem of species separation in minimal time. Droop model is considered to describe the evolution of two distinct populations of microorganisms that are in competition for the same resource in a photobioreactor. We focus on an optimal control problem (OCP) subject to a five-dimensional controlled system in which the control represents the dilution rate of the chemostat. The objective is to select the desired species in minimal-time and to synthesize an optimal feedback control. This is a very challenging issue, since we are dealing with a ten-dimensional optimality system. We provide properties of optimal controls allowing the strain of interest to dominate the population. Our analysis is based on the Pontryagin Maximum Principle (PMP), along with a thorough study of singular arcs that is crucial in the synthesis of optimal controls. These theoretical results are also extensively illustrated and validated using a direct method in optimal control (via the Bocop software for numerically solving optimal control problems). The approach is illustrated with numerical examples with microalgae, reflecting the complexity of the optimal control structure and the richness of the dynamical behavior.

**Keywords:** optimal control; modelling; microalgae; chemostat; nonlinear control; Pontryagin's principle; singular control; Droop model; photobioreactor

## 1. Introduction

The interaction between species coexisting in an ecosystem is complex and affected by external factors. Depending on their environment, some species will dominate, while others, less adapted, will progressively decline. This Darwinian pressure, when it can be manipulated [1], provides the opportunity of guiding the evolution of species of interest. This concept can be applied to artificial ecosystems to select individuals with a desired trait. Here, we focus on microalgae, unicellular photosynthetic microorganisms with promising potential for industrial applications [2,3]. The great biodiversity of microalgae opens the door for a large range of applications [4]. They are grown for their pigments, antioxidants or essential fatty acids [5], and, over the longer term, their efficient way of producing proteins, bricks for green chemistry, biofuel and CO₂ mitigation [2,6,7]. To date, microalgae do not have the place they deserve in biotechnology (see, e.g., [8–10]) and many optimization steps must be carried out to improve the economic and environmental performances of these processes at a large scale [6,11]. Currently, only wild organisms sampled in nature are used on an industrial scale. One of the key challenges is to improve the productivity of these strains.

Species in agriculture have been improved after centuries of selection and hybridization. The objective of this work is to develop an alternative approach adapted to microorganisms to select, on a shorter time scale, more productive microalgae strains, by Darwinian

<footnote n="1">¹ Université Côte d'Azur, Inria, INRAE, CNRS, Sorbonne Université, Biocore Team, 06902 Sophia Antipolis—Valbonne, France; olivier.bernard@inria.fr</footnote>

<footnote n="2">² Laboratoire de Mathématiques d'Avignon (EA 2151), Avignon Université, 84018 Avignon, France;
terence.bayen@univ-avignon.fr</footnote>

<footnote n="3">³ Laboratoire d'Océanographie de Villefranche (LOV), CNRS, Sorbonne University, IMEV, 06230 Villefranche-sur-Mer, France</footnote>

<footnote>* Correspondence: walid.djema@inria.fr</footnote>
pressure. The idea is based on the competitive exclusion principle in a continuous reactor [12] stating that the species which more efficiently uses the available resources will win the competition. The conditions for a bacterial or microalgal species to win the competition for a limiting substrate have been well established, and the outcome of the competition is known to depend on the minimum requested substrate to support a growth rate equal to the dilution rate [12,13].

Experimental works for maintaining a long-term selection process to favour individuals of interest have already been carried out [14,15]. Since the experiments can last several months or several years, these approaches are costly in time. There is a margin of improvement by applying optimal control theory [16] to enhance the selection process for N strains competing for the same resource (the control parameter being the dilution rate).

One main issue is to decrease the operating time when the species of interest starts to dominate. Several works addressed the question of improving the selection process in minimal time in the case of the chemostat system with Monod’s laws [17,18]. Microalgae are more complicated microorganisms better represented by the Droop model, taking into account the internal accumulation of the limiting nutrient [19,20]. Such a model for two strains in competition leads to a five-dimensional problem. The minimal time selection problem with this model is the main focus of the paper. It has been tackled in [21,22] after a simplification allowing for reducing the model dimension. This necessitated to oversimplify the initial dynamics which can play a role in the minimal time selection.

Optimal control [23] strategies ensuring the domination of the strain of interest are derived using the Pontryagin maximum principle [24]. Since the system is affine w.r.t. the control, we obtain various possible structures for an optimal control, namely, the concatenation of several bang arcs or of a bang arc with a singular arc of first order satisfying Legendre–Clebsch’s condition. The paper is structured as follows: in Section 2, we introduce the model and present the optimal control problem. We also prove the reachability of the target set. In Section 3, we make explicit the necessary conditions provided by the Pontryagin Maximum Principle and we introduce properties of the switching function. A thorough study of singular arcs is provided in Section 4 thanks to geometrical control theory. The paper is concluded with numerical simulations of optimal strategies using a direct method in Section 5.

## 2. The Optimal Control Problem (OCP)

### 2.1. Droop Model and Main Assumptions

We consider the Droop model [19]. This emblematic variable yield model represents the growth rate of microorganisms which can intracellularly store nutrients. When two strains are in competition, it results in a five-dimensional system. The growth of each strain depends on the intracellular quota-storage (q-variable) of the limiting nutrient (s-variable). More precisely, when two species/strains, of biomass concentrations $x_1$ and $x_2$ are competing for one limiting nutrient s in the bioreactor, the Droop model reads as follows:

$$ 
\left\{ 
\begin{aligned}
\dot{s} &= (s_{in} - s)D(t) - \rho_1(s)x_1 - \rho_2(s)x_2, \\
\dot{q}_1 &= \rho_1(s) - \mu_1(q_1)q_1, \\
\dot{x}_1 &= [\mu_1(q_1) - D(t)]x_1, \\
\dot{q}_2 &= \rho_2(s) - \mu_2(q_2)q_2, \\
\dot{x}_2 &= [\mu_2(q_2) - D(t)]x_2,
\end{aligned}
\right.
\quad (1) $$

where $q_i$ is the quota storage of the i-th species and $s_{in}$ is the input substrate concentration. The dilution rate $D(\cdot)$ is a bounded non-negative control function such that $D(t) \in [0, D_{max}]$, where $D_{max} > 0$ is the maximal admissible value of the dilution rate, above the maximum actual growth rates [22] (this will be made more precise in Section 2.3) of the two species as shown in Figure 1. In addition, for $i = 1, 2$, $\rho_i$ is a non-negative function representing the rate of substrate absorption, i.e., the uptake rate of the free nutrient
s, and $\mu_i$ is also a non-negative function representing the growth rate of the i-th species (see [25]).

Following for instance [25], we suppose that the uptake rates $\rho_i(s)$ are expressed as,

$$ \rho_i(s) = \frac{\rho_{m_i} s}{K_i + s'} \qquad (2) $$

which corresponds to Michaelis-Menten's kinetics. Here, the parameters $K_i$ and $\rho_{m_i}$ are positive, $i = 1,2$. In addition, we assume that, for $i = 1,2$, there is $k_i > 0$ such that cell division does not occur if $k < k_i$. Concerning the Droop model, the kinetics $\mu_i$ are defined by

$$ \left| \begin{array}{ll} \mu_i(q_i) = 0, & 0 \le q_i \le k_i, \\ \mu_i(q_i) = \mu_{i\infty} \left(1 - \frac{k_i}{q_i}\right), & k_i \le q_i. \end{array} \right. \qquad (3) $$

In addition, for $i = 1,2$, let $M_i$ stand for,

$$ M_i := \sup_{s \in [0, s_{in}]} \rho_i(s) = \rho_i(s_{in}) < \rho_{m_i}, $$

and let $\bar{q}_1, \bar{q}_2$ be such that $\mu_i(\bar{q}_i)\bar{q}_i = M_i, i = 1,2$ (observe that $\bar{q}_1, \bar{q}_2$ are uniquely defined). Thus, one has,

$$ \rho_i(s_{in}) = \mu_i(\bar{q}_i)\bar{q}_i. $$

System (1) satisfies the following invariance property.

**Proposition 1.** For every $q_{m_1} \ge \bar{q}_1$, and for every $q_{m_2} \ge \bar{q}_2$, the set

$$ \Omega := (0, s_{in}) \times [k_1, q_{m_1}] \times \mathbb{R}_+^* \times [k_2, q_{m_2}] \times \mathbb{R}_+^*, \qquad (4) $$

is forward invariant by (1).

**Proof.** First, observe that, for $i = 1,2$, $x_i$ never vanishes whenever $x_i^0 = x_i(0) > 0$. Now, $(0, s_{in})$ is clearly invariant by the dynamics of $s(\cdot)$ since $\dot{s} \ge 0$ (resp. $\dot{s} \le 0$) whenever $s = 0$ (resp. $s = s_{in}$). Similarly, for $i = 1,2$:

$$ \begin{array}{ll} q_i = k_i & \Rightarrow \dot{q}_i = \rho_i(s) > 0, \\ q_i = q_{m_i} & \Rightarrow \dot{q}_i \le M_i - \mu_i(q_{m_i})q_{m_i} \le 0, \end{array} $$

where the last inequality follows from the choice of $M_i$ and the fact that $q_{m_i} \ge \bar{q}_i, i = 1,2$. This ends the proof. $\square$

The parameter $q_{m_i}$ represents the maximum internal storage quota. Since $\Omega$ is invariant by (1) (Proposition 1), we notice that $\bar{q}_i$ stands for the effective maximum internal storage quota for $s \in (0, s_{in})$. Thus, in the sequel, we consider without loss of generality that $q_{m_i} = \bar{q}_i$ for $i = 1,2$. In the sequel, we also assume that $\rho_1, \rho_2$ fulfill the following hypothesis:

**Assumption 1.** The affinity of species 1 for the substrate is higher than the one of species 2, i.e., $K_{s_2} > K_{s_1}$, or equivalently:

$$ \frac{\rho_2''}{\rho_2'} > \frac{\rho_1''}{\rho_1'}. \qquad (5) $$

From (2), we deduce that, for $s \ge 0$,

$$ \rho_2''(s)\rho_1'(s) - \rho_1''(s)\rho_2'(s) = \frac{2K_1K_2\rho_{m_1}\rho_{m_2}(K_2 - K_1)}{(K_1 + s)^3(K_2 + s)^3}. $$

It means that species 1 with the lower $K_i$ absorbs nutrients slightly faster.
We are now in a position to formulate the OCP of interest.

## 2.2. Statement of the Optimal Control Problem (OCP)

In this work, we suppose that the first species (with biomass concentration $x_1$) is the one of interest. Our aim is to compute the best feeding strategy, that is, the optimal (dilution rate) control function $D(\cdot)$, such that $x_1$ becomes predominant in the photobioreactor in minimal-time. This can be formulated and quantified in terms of the ratio between the two competing species. Intuitively, we wish to find an adequate control strategy (if possible optimal) $D(\cdot)$ for which, at the end of the process, we have $\frac{x_1}{x_2} \gg 1$.

Firstly, the set of admissible controls is defined as,

$$ \mathcal{D} := \{D : [0, +\infty) \to [0, D_{\max}] ; D(\cdot) \in \mathcal{L}_{\text{loc}}^{\infty}(\mathbb{R}_{+})\}, $$

where $\mathcal{L}_{\text{loc}}^{\infty}(\mathbb{R}_{+})$ is the space of locally integrable functions on every compact on $\mathbb{R}_{+}$ and $D_{\max} > 0$ is the maximum pump feeding capacity. In practice, $D_{\max}$ is designed above the maximum growth rates of the coexisting species (see Section 2.3).

To handle the selection process between the two species, let us define a subset $\mathcal{T}$ of $\Omega$ as,

$$ \mathcal{T} := \{ X := (s, q_1, x_1, q_2, x_2) \in \Omega ; x_2 \le \varepsilon x_1 \}. $$

We choose the parameter $\varepsilon > 0$ such that $\varepsilon \ll 1$ in such a way to quantify the contamination rate of the interesting strain $x_1$. Whenever a trajectory reaches the target set $\mathcal{T}$, this means that the biomass of the first species is significantly greater than the other one when reaching the target $\mathcal{T}$ at the terminal time (if possible).

**Objective 1.** The optimal control problem (OCP) can then be stated: determine a dilution-based control strategy $D(\cdot)$ in such a way that trajectories of (1) starting from an initial condition within the set $\Omega$ reach the target set $\mathcal{T}$ in minimal-time, i.e.,

$$ \inf_{D \in \mathcal{D}} t_f^D \quad \text{s.t.} \quad X(t_f^D) \in \mathcal{T} \quad \text{and} \quad X^0 \in \Omega, \qquad (6) $$

where $X(\cdot)$ is the unique solution of (1) associated with $D(\cdot) \in \mathcal{D}$ such that $X(0) = X^0 \in \Omega$ and $t_f^D \in (0, +\infty]$ is the first entry time of $X(\cdot)$ into the target set. In the sequel, we will use the simpler notation $t_f$ instead of $t_f^D$.

In other words, for every positive initial conditions $X^0 = (s^0, q_1^0, x_1^0, q_2^0, x_2^0)$ such that $q_i^0 > k_i$, we are seeking an admissible control strategy $D = D_{X_0} \in \mathcal{D}$, steering the trajectory $X(t)$ of the system (1) from $X^0$ to the target set $\mathcal{T}$ in minimal-time, for a fixed $D_{\max}$ (Figure 1) and a given contamination rate $\varepsilon \ll 1$. Note that, if one is able to synthesize such an optimal control for every $X^0 \in \Omega$, then one is able to construct an optimal feedback control over $\Omega$ as $X_0 \mapsto D_{X_0}(0)$. Such an optimal control problem falls into the class of minimal-time control problems governed by a mono-input affine controlled system, for which the synthesis of an optimal feedback control, thanks to geometric control theory, is a crucial (but also delicate) issue. In particular, handling the high dimension of the Droop model in competition and its resulting optimality system is challenging. Note also that the linearity of the problem w.r.t. $D$ (in contrast, for instance, with strictly convex cost functionals) leads to technicalities because singular arcs usually occur in this setting, see Section 4.

## 2.3. Basic Properties

We now introduce the so-called actual growth rates. These key functions will have an important role in the optimal separation strategy. For that, let us firstly start with the following observations:

* The mapping $\rho_i : [0, s_{in}] \to [0, M_i]$ is one-to-one with $M_i < \rho_{m_i}$;
* For every $q_i \in [k_i, \bar{q}_i]$, one has $q_i\mu_i(q_i) \in [0, M_i]$ and $q_i \mapsto q_i\mu_i(q_i)$ is one-to-one from $[k_i, \bar{q}_i]$ into $[0, M_i]$;

* It follows that the composition $\delta_i: [k_i, \bar{q}_i] \to [0, s_{in}]$ such that

$$ \delta_i(q_i) := \rho_i^{-1}(q_i\mu_i(q_i)) = \frac{K_i q_i \mu_i(q_i)}{\rho_{m_i} - q_i \mu_i(q_i)}, \quad q_i \in [k_i, \bar{q}_i], $$

is well-defined and is also one-to-one;

* Hence, the mapping $\delta_i^{-1}: [0, s_{in}] \to [k_i, \bar{q}_i]$ such that

$$ \delta_i^{-1}(s) := \frac{k_i K_i \mu_{i\infty} + (\rho_{m_i} + \mu_{i\infty} k_i)s}{\mu_{i\infty}(K_i + s)}, \quad s \in [0, s_{in}], $$

is well-defined over $[0, s_{in}]$ with values in $[k_i, \bar{q}_i]$ and is one-to-one.

From these observations, one can immediately check that the mappings $\delta_i, i = 1, 2$ and $\delta_i^{-1}$ are increasing.

Indeed, for $i = 1$, one can write $\delta_1^{-1}(s) = c_1 + \frac{c_2}{s+c_3}$ with $c_1, c_3 > 0$ and

$$ c_2 := \delta_i^{-1}(s) = C + K_{s_1} \left[ \frac{k_{k_1} \mu_{1\infty}}{k_1 \mu_{1\infty} + \rho_{m_1}} - 1 \right] < 0, $$

which implies that $\delta_i^{-1}$ is increasing.

The actual growth rate of species $i$ is then defined as the mapping $\mu_i \circ \delta_i^{-1}$,

$$ \mu_i(\delta_i^{-1}(s)) = \frac{\rho_{m_i} \mu_{i\infty}s}{k_i K_i \mu_{i\infty} + (\rho_{m_i} + \mu_{i\infty} k_i)s}, \quad s \in [0, s_{in}], \quad i = 1, 2. $$

The resulting generic functions are illustrated in Figure 1: Let us now define,

$$ \Delta(s) := \mu_1(\delta_1^{-1}(s)) - \mu_2(\delta_2^{-1}(s)), \quad s \in [0, s_{in}]. \qquad (7) $$

Throughout the paper, we suppose that $\Delta$ satisfies the following assumption.

**Assumption 2.** There is a unique $\hat{s} \in (0, s_{in})$ such that $\Delta(s) > 0$ for every $s \in (0, \hat{s})$ and $\Delta(s) < 0$ for every $s \in (\hat{s}, s_{in})$. In addition, $\Delta$ has a unique maximum $s_c \in [0, \hat{s}]$.

Taking into account that $\Delta(\hat{s}) = 0$, the inequalities satisfied by $\Delta$ according to Assumption 2 can also be written:

$$ \Delta(s) (s - \hat{s}) < 0, \quad s \in (0, s_{in}) \setminus \{\hat{s}\}. $$

If we assume that $s$ is regulated to $s^*(t) = s_c$ using an appropriate control $D$, we notice that the $q$-variables are regulated to some unique $q_{ic} \in [k_i, q_{mi}]$, for $i = 1, 2$. The unique point $(s_c, q_{1c}, q_{2c})$, where $s_c \in [0, s_{in}]$, plays a crucial role in the optimal control strategy of (OCP) as discussed in Section 5. Finally, $D_{max}$ (the maximum dilution rate) is assumed to be large enough in order to drive competition between the two species. More precisely, we assume that $D_{max}$ satisfies the hypothesis:

**Assumption 3.** The maximal value of the dilution rate $D_{max}$ satisfies

$$ \forall s \in [0, s_{in}], \quad D_{max} > \max(\mu_1(\delta_1^{-1}(s)), \mu_2(\delta_2^{-1}(s))). $$

These assumptions will ensure reachability (as detailed in the next section) of the target $\mathcal{T}$, and establishes a generic framework where both species may win the competition for sufficiently large time (considering for instance various constant control parameters $D$ that favor species 1 or 2). Thus, under these considerations, we ensure the well-posedness
of the optimal control problem of interest. In this general framework, the objective is then to determine the optimal D that steers the trajectories in minimal-time to the desired target.

<figure><img src="image_5.png" /><figcaption>Figure 1. The illustrated absorption functions $\rho_i$ and growth rates $\mu_i$ lead to a generic form of the actual growth rates $s \mapsto \mu_i(\delta_i(s))$ where both species may win the competition. The maximum dilution rate $D_{\max}$ is a fixed constant value above the maximum of the functions $s \mapsto \mu_i(\delta_i(s))$ for $i = 1,2$, as stated in Assumption 3.</figcaption></figure>

## 2.4. Reachability of the Target

Our next aim is to show that the target is reachable from every initial condition. First, let us recall that the set

$$ \mathcal{M} := \{X \in \Omega; x_1q_1 + x_2q_2 + s = s_{in}\} \qquad (8) $$

is an invariant and attractive manifold for (1) for a given persistently exciting control (i.e., an admissible control function $D(\cdot)$ such that $\int_0^{+\infty} D(t) dt = +\infty$).

**Proposition 2.** For every initial condition $X^0 \in \Omega$, there exists an admissible control $D(\cdot)$ and a time $t_e \ge 0$ such that $X(t_e) \in \mathcal{T}$, where $X(\cdot)$ is the unique solution of (1), starting from $X^0$, associated with $D(\cdot)$.

**Proof.** Let $s^\dagger \in (0, \hat{s})$ and $X^0 \in \Omega$. Without any loss of generality, we may assume that $s(0) = s^\dagger$. Indeed, observe that $s = s^\dagger$ is not a steady-state of $\hat{s}$ whenever $D = 0$ over $\mathbb{R}_+$ or $D = D_{\max}$ over $\mathbb{R}_+$. Thus, if we apply $D = 0$ (in that case $\hat{s} < 0$) or $D = D_{\max}$ (in that case $\hat{s} > 0$), then $s(t) = s^\dagger$ is reached in a finite horizon. Consider the feedback control function,

$$ D^{\dagger}(x_1, x_2) := \frac{\rho_1(s^{\dagger})x_1 + \rho_2(s^{\dagger})x_2}{s_{in} - s^{\dagger}}, $$

in such a way that the unique solution of (1) associated with this control satisfies $s(t) = s^\dagger$ for every time $t \ge 0$ (Cauchy-Lipschitz’s Theorem). We claim that there exists $t_1 \ge 0$ large
enough such that $D^{\dagger}(x_1(t), x_2(t)) \in [0, D_{max}]$ for every time $t \ge t_1$. Indeed, when $t \to +\infty$, taking into account (8), it follows:

$$D^{\dagger}(x_1(t), x_2(t)) \sim \bar{D}(t) := \frac{\rho_1(s^{\dagger})x_1(t) + \rho_2(s^{\dagger})x_2(t)}{x_1(t)q_1(t) + x_2(t)q_2(t)}.$$

Now, when $t \to +\infty$, it is easy to see that $q_1(t)$ converges to $\frac{\rho_1(s^{\dagger})}{\mu_1^{\infty}} + k_1$, since $q_1$ satisfies the linear ODE $\dot{q}_1 + \mu_1^{\infty}q_1 = \rho_1(s^{\dagger}) + \mu_1^{\infty}k_1$. At steady-state, we thus have

$$\rho_1(s^{\dagger}) = q_1\mu_1(q_1).$$

In conclusion, when $t \to +\infty$,

$$\bar{D}(t) \sim \frac{q_1\mu_1(q_1((t))x_1(t) + q_2(t)\mu_2(q_2(t))x_2(t))}{x_1(t)q_1(t) + x_2(t)q_2(t)},$$

or, equivalently, using that, at steady state, $s^{\dagger} = \rho_1^{-1}(q_1\mu_1(q_1))$ that is, $q_1 = \delta_i^{-1}(s^{\dagger})$, we end up with

$$\bar{D}(t) \sim \frac{\mu_1(\delta_1^{-1}(s^{\dagger}))q_1(t)x_1(t) + \mu_2(\delta_2^{-1}(s^{\dagger}))q_2(t)x_2(t)}{x_1(t)q_1(t) + x_2(t)q_2(t)}.$$

Thanks to Assumption 3, this last expression is upper bounded by $D_{max}$ for $t$ large enough, which proves our claim. Finally, posit $y_i := \ln(x_i)$ and observe that

$$\dot{y}_1 - \dot{y}_2 = \Delta(s^{\dagger}) > 0.$$

It follows that $y_1(t) - y_2(t) \to +\infty$ when $t \to +\infty$, which implies that $\lim_{t \to +\infty} \frac{x_2(t)}{x_1(t)} = 0$. This ends the proof. $\square$

## 2.5. Motivation of Studying the OCP

Thanks to Proposition 2, the target set is reachable from any initial condition, thus the existence of an optimal control of (6) is standard (namely because the dynamics is affine w.r.t. the control): it is an application of the Fillipov Theorem, see, e.g., [26–28]. Considering such a control as in the proof of Proposition 2 then indeed allows the system to let the species of interest dominate the reactor, but this process can be long (see, for instance, Example 1). Another possible strategy is to use a constant control $D$. Following [12], depending on the value of $D$, species 1 may win the competition, i.e.,

$$\lim_{t \to +\infty} x_1(t) > 0 \quad ; \quad \lim_{t \to +\infty} x_2(t) = 0.$$

In that case, this (simple) strategy indeed allows for reaching the target. However, this convergence is asymptotic and depends on the value of $D$ as in the competitive exclusion principle. Roughly speaking, if $D > \mu_1(\tilde{q}_1)$ (where $\tilde{q}_1$ and $s^{\dagger}$ are such that $\tilde{q}_1 := \delta_1^{-1}(s^{\dagger})$ and $(\mu_1 \circ \delta_1^{-1})(s^{\dagger}) = (\mu_2 \circ \delta_2^{-1})(s^{\dagger})$), then species 2 wins the competition, whereas, if $D < \mu_1(\tilde{q}_1)$, species 1 wins the competition. We refer to [12] for more details about the asymptotic behavior of (1) for a constant control $D$. Thus, the target set may not always be reachable with a constant control $D$. Our objective in this paper is precisely to propose a methodology to compute a control strategy to reach the target set $\mathcal{T}$ faster, playing on the control $D(\cdot)$ as illustrated in Figure 2.

**Example 1.** Let us consider the following parameters: $\rho_{1m} = 0.8$, $\rho_{2m} = 0.95$, $K_{s1} = 1$, $K_{s2} = 1.4$, $k_1 = 1.1$, $k_2 = 1.4$, $\mu_{1\infty} = 1.8$, $\mu_{2\infty} = 1.7$, $s_{in} = 10$. The contamination rate is fixed to $\varepsilon = 0.05$. The initial conditions are given by: $s^0 = 2$, $q_i^0 = 2.5$ and $x_i^0 = 1$. As illustrated in Figure 2, the target $\mathcal{T}$ is reached after $t_f = 46.774$ days using the control $D(t) = D_{opt}$, while $x_1$ dominates the culture after $t_f = 62.68$ days using the constant control $D = 0.48$. Let us also point out
that, using an arbitrary constant control $D \in [0, D_{\max}]$, we are not even sure that $x_1$ wins the competition (trajectories do not reach the target in that case).

<figure><img src="image_6.png" /><figcaption>Figure 2. It is possible to find a constant feeding strategy $D \in [0, D_{\max}]$ that could favor one species or the other. However, this is a risky time-consuming process, while the target is reached much faster using an optimized control $D_{opt}$, thus saving several days of costly cultures.</figcaption></figure>

### 3. Necessary Conditions on Optimal Controls

We start this section by generalizing previously obtained results [22,29] characterizing optimal solutions of (6). For that, we apply the PMP which allows for obtaining necessary conditions satisfied by optimal controls of (6). We denote by $X = (s, q_1, x_1, q_2, x_2)$ and $\lambda = (\lambda_s, \lambda_{q_1}, \lambda_{x_1}, \lambda_{q_2}, \lambda_{x_2})$, respectively, the state and adjoint variables (also called co-state or covector). The Hamiltonian associated with the optimal control problem

$$H = H(s, q_1, x_1, q_2, x_2, \lambda_s, \lambda_{q_1}, \lambda_{x_1}, \lambda_{q_2}, \lambda_{x_2}, \lambda^0, D),$$

is given by

$$H = \lambda^0 - (\rho_1(s)x_1 + \rho_2(s)x_2)\lambda_s + (\rho_1(s) - \mu_1(q_1)q_1)\lambda_{q_1} + \mu_1(q_1)x_1\lambda_{x_1} + (\rho_2(s) - \mu_2(q_2)q_2)\lambda_{q_2} + \mu_2(q_2)x_2\lambda_{x_2} + D[(s_{in} - s)\lambda_s - x_1\lambda_{x_1} - x_2\lambda_{x_2}].$$

Let $X^0 \in \Omega \setminus \mathcal{T}$ and let $(X(\cdot), u(\cdot))$ be an optimal pair such that $X(\cdot)$ reaches the set $\mathcal{T}$ in a time $t_f \ge 0$. Thanks to the PMP, there exist an absolutely-continuous map $\lambda: [0, t_f] \to \mathbb{R}^5$ and $\lambda^0 \le 0$ such that:

* The pair $(\lambda(\cdot), \lambda^0)$ is non-trivial, i.e., $(\lambda(\cdot), \lambda^0) \ne (0, 0)$.

* The covector satisfies

$$\dot{\lambda}(t) = -\nabla_X H(X(t), \lambda(t), \lambda^0, D(t)) \quad \text{a.e. } t \in [0, t_f]. \qquad (9)$$

* The Hamiltonian maximization condition writes

$$D(t) \in \arg\max_{\xi \in [0, D_{\max}]} H(X(t), \lambda(t), \lambda^0, \xi) \quad \text{a.e. } t \in [0, t_f]. \qquad (10)$$

* At the terminal time, the transversality condition writes:

$$\lambda(t_f) \in -N_{\mathcal{T}}(X(t_f)). \qquad (11)$$

Here, $N_{\mathcal{T}}(X) = \{p \in \mathbb{R}^5; \forall Y \in \mathcal{T}, p \cdot (Y - X) \le 0\}$ stands for the normal cone to $\mathcal{T}$ at some point $X \in \mathcal{T}$, see [28]. The adjoint Equation (9) is equivalent to:

$$\left\{ \begin{array}{l} \dot{\lambda}_s = (\rho'_1(s)x_1 + \rho'_2(s)x_2)\lambda_s - \rho'_1(s)\lambda_{q_1} - \rho'_2(s)\lambda_{q_2} + D\lambda_s, \\ \dot{\lambda}_{q_1} = \mu_{1\infty}\lambda_{q_1} - \mu'_1(q_1)x_1\lambda_{x_1}, \\ \dot{\lambda}_{x_1} = \rho_1(s)\lambda_s - \mu_1(q_1)\lambda_{x_1} + D\lambda_{x_1}, \\ \dot{\lambda}_{q_2} = \mu_{2\infty}\lambda_{q_2} - \mu'_2(q_2)x_2\lambda_{x_2}, \\ \dot{\lambda}_{x_2} = \rho_2(s)\lambda_s - \mu_2(q_2)\lambda_{x_2} + D\lambda_{x_2}. \end{array} \right. \qquad (12)$$
We define an extremal as a quadruplet $(X(\cdot), \lambda(\cdot), \lambda^0, D(\cdot))$ such that $(\lambda(\cdot), \lambda^0)$ is non zero and such that (1) and (9)-(11) is verified. Whenever $\lambda^0 = 0$, we say that the extremal is abnormal ; if $\lambda^0 \neq 0$, then we say that the extremal is normal. Because (6) is autonomous (i.e., the system does not depend explicitly on time), the Hamiltonian computed along an extremal is constant. In addition, since the terminal time is not fixed, we classically obtain, following optimal control theory, that $H = 0$.

The transversality condition is crucial for obtaining properties on optimal controls by reasoning backward in time from the terminal time $t = t_f$. We shall next extend earlier results [22] by taking into account explicitly the fact that $\mathcal{T}$ is a half-space of $\mathbb{R}^5$ and exploiting that $X(t_f)$ belongs to the set $E := \{X \in \mathbb{R}^5; x_2 - \varepsilon x_1 = 0\}$ (the boundary of the target set). Then, condition (11) can be transformed more explicitly as follows. At $X(t_f)$, the normal cone to $\mathcal{T}$ writes

$$N_{\mathcal{T}}(X(t_f)) = \mathbb{R}_+(0,0, -\varepsilon, 0, 1).$$

Therefore, inclusion (11) is then equivalent to

$$\lambda_s(t_f) = \lambda_{q_1}(t_f) = \lambda_{q_2}(t_f) = 0, \qquad (13)$$

together with

$$\lambda_{x_1}(t_f) + \varepsilon\lambda_{x_2}(t_f) = 0, \qquad (14)$$

and the inequalities $\lambda_{x_1}(t_f) \geq 0$ and $\lambda_{x_2}(t_f) \leq 0$. Actually, one has $\lambda_{x_1}(t_f) > 0$ and $\lambda_{x_2}(t_f) < 0$. Suppose indeed that $\lambda_{x_1}(t_f) = 0$. Then, (14) would imply $\lambda_{x_2}(t_f) = 0$, and, since (9) is linear w.r.t. $\lambda$, we would obtain $\lambda \equiv 0$ over $[0, t_f]$. Using the constancy of $H$, we would also obtain $\lambda^0 = 0$ contradicting the PMP. We can then conclude that

$$\lambda_{x_1}(t_f) > 0 \text{ and } \lambda_{x_2}(t_f) < 0. \qquad (15)$$

Throughout the paper, we suppose that only normal extremals occur, i.e., $\lambda^0 < 0$, and, without any loss of generality, we may assume that $\lambda^0 = -1$ (up to a renormalization of the necessary conditions that are linear w.r.t. $\lambda$).

**Remark 1.** Abnormal extremals are not generic. They correspond to the optimal path reaching the target set in some particular subset of the target set and are such that $\lambda^0 = 0$ (and thus $\lambda(\cdot) \neq 0$). From the conservation of $H$, this implies that $\mu_1(q_1(t_f)) = \mu_2(q_2(t_f))$ in such a way that $\lambda(t_f)$ is not uniquely defined in contrast with the normal case (see below). At such a singular point, the value function (the minimal time as a function of $X^0$) is also non-differentiable.

Going back to the normal case, i.e., $\lambda^0 = -1$, the covector $\lambda$ at $t = t_f$ can be completely determined (thanks to the conservation of $H$) as follows:

$$\begin{align}
\lambda_s(t_f) &= \lambda_{q_1}(t_f) = \lambda_{q_2}(t_f) = 0, \nonumber \\
\lambda_{x_1}(t_f) &= \frac{1}{x_1(t_f)(\mu_1(q_1(t_f)) - \mu_2(q_2(t_f)))} > 0, \tag{16} \\
\lambda_{x_2}(t_f) &= -\frac{\lambda_{x_1}(t_f)}{\varepsilon} < 0. \nonumber
\end{align}$$

Notice that the quantity $\mu_1(q_1(t_f)) - \mu_2(q_2(t_f))$ is non-zero along a normal extremal. As a consequence, the transversality condition (11) coupled with the conservation of $H$ are equivalent to (16). The computation of $\lambda$ at $t = t_f$ is useful to integrate the state adjoint system backward in time from the target set.

We now wish to exploit the Hamiltonian condition (10). It is of particular interest to introduce the *switching function*, which allows us to determine the optimal control $D$
according to the sign of the switching function. For that, let us denote by $\tilde{\Phi}$ the switching function:

$$ \tilde{\Phi} := -(x_1\lambda_{x_1} + x_2\lambda_{x_2}) + (s_{in} - s)\lambda_s, $$

associated with the control function $D(\cdot)$.

$$ \begin{aligned} \tilde{\Phi}(t) > 0 & \Rightarrow D(t) = D_{max}, \\ \tilde{\Phi}(t) < 0 & \Rightarrow D(t) = 0, \end{aligned} \qquad (17) $$

for a.e. $t \in [0, t_f]$.

In the case where $\tilde{\Phi} > 0$, resp. $\tilde{\Phi} < 0$ over a time interval $[t_1, t_2]$, we say that the optimal control $u$ is of bang type (denoted by $B_+$, resp. $B_-$). When the control $D(\cdot)$ is non-constant in every neighborhood of a time $t_c \in (0, t_c)$, we say that $t_c$ is a switching time, and one must have $\tilde{\Phi}(t_c) = 0$. Next, when the switching function $\tilde{\Phi}$ vanishes over a time-interval $[t_1, t_2]$, we state that a singular arc occurs. In this case, the corresponding trajectory is singular over $[t_1, t_2]$, and such an arc will be denoted by $S$. Singular arcs are essential to optimize the time to steer an initial condition to the target set. Now, we are ready to state some main features of the switching function $\tilde{\Phi}$, and then investigate properties of the singular paths.

**Lemma 1.**

(i) The function $\tilde{\Phi}$ is continuously differentiable over $[0, t_f]$ and, moreover,

$$ \dot{\tilde{\Phi}} = (s_{in} - s)(\rho_1'(s)[x_1\lambda_s - \lambda_{q_1}] + \rho_2'(s)[x_2\lambda_s - \lambda_{q_2}]). \qquad (18) $$

(ii) At the terminal time $t = t_f$, it holds:

$$ \tilde{\Phi}(t_f) = \dot{\tilde{\Phi}}(t_f) = 0. \qquad (19) $$

**Proof.** By differentiating $\tilde{\Phi}$ w.r.t. $t$, we find that

$$ \dot{\tilde{\Phi}} = -\dot{s}\lambda_s + (s_{in} - s)\dot{\lambda}_s - [\dot{x}_1\lambda_{x_1} + \dot{x}_2\lambda_{x_2}]. $$

Using (1) and (9), we obtain (18), which proves (i). For proving (ii), note that $x(t_f) \in E$ and $\lambda(t_f) \in E^\perp$, thus $\tilde{\Phi}(t_f) = 0$. Using (13), we also obtain $\dot{\tilde{\Phi}}(t_f) = 0$, which ends the proof. $\square$

**Remark 2.** Following the formalism of geometric control theory, $\tilde{\Phi}$ never involves $D$ explicitly, but $D$ is present in the expression of $\tilde{\Phi}^{(2k)}$, $k \ge 1$, see [27]. If $k \ge 1$ is the first integer for which the control is present in the expression of $\tilde{\Phi}^{(2k)}$, we usually say that the singular arc is of order $k$.

At this step, an optimal control is a concatenation of bang and singular arcs:

$$ B_{\pm}, B_{\pm}B_{\mp}, B_{\pm}S, B_{\pm}B_{\mp}S, B_{\pm}B_{\mp}SB_{\pm}, \dots $$

with possibly infinitely many crossing times (in particular, if there is a singular arc of order 2 [27]). The occurrence and properties of singular arcs as well as the various (possible) structure for an optimal control of (6) will be precisely the matter of the next section. The goal is to reduce (if possible) the number of possible structures for an optimal control.

**4. Singular Arcs and Insights into Optimal Solutions**

**4.1. Legendre–Clsebsch’s Necessary Condition and Computation of the Singular Control**

The analysis of singular arcs requires to compute $\tilde{\Phi}$. Indeed, the Hamiltonian condition does not give any information about an optimal control during a singular phase. Thanks to
the next computations, we shall also be able to deduce an expression of the singular control during a singular arc.

Doing so, let us define $\theta_1, \theta_2 : [0, T] \to \mathbb{R}$ as

$$ \theta_1 := \rho'_1(s)[x_1\lambda_s - \lambda_{q_1}] + \rho'_2(s)[x_2\lambda_s - \lambda_{q_2}], $$

$$ \theta_2 := \rho''_1(s)[x_1\lambda_s - \lambda_{q_1}] + \rho''_2(s)[x_2\lambda_s - \lambda_{q_2}]. $$

Hereafter, to simplify the layout, we do not write the dependency of $\theta_i, s, x_i, \lambda_s, \lambda_{q_i}$ w.r.t. the time. To shorten the notation, we also do not write explicitly the dependency of certain functions w.r.t. some variables. Using (12)–(18), one can write

$$ \dot{\tilde{\Phi}} = (s_{in} - s)\theta_1 \quad \text{and} \quad \dot{\lambda}_s = \theta_1 + D\lambda_s. \qquad (20) $$

**Lemma 2.**

(i) The derivative of $\theta_1$ can be expressed as:

$$ \dot{\theta}_1 = \theta_2 \dot{s} + [\rho'_1 x_1 + \rho'_2 x_2]\theta_1 + \lambda_s[x_1 \mu_1 \rho'_1 + x_2 \mu_2 \rho'_2] \\ - \rho'_1[\mu_{1\infty} \lambda_{q_1} - \mu'_1(q_1) x_1 \lambda_{x_1}] - \rho'_2[\mu_{2\infty} \lambda_{q_2} - \mu'_2(q_2) x_2 \lambda_{x_2}]. \qquad (21) $$

(ii) The second derivative of $\tilde{\Phi}$ fulfills the equality:

$$ \ddot{\tilde{\Phi}} = [-\theta_1 + (s_{in} - s)\theta_2]\dot{s} + (s_{in} - s)[\rho'_1 x_1 + \rho'_2 x_2]\theta_1 + (s_{in} - s)\lambda_s[x_1 \mu_1 \rho'_1 + x_2 \mu_2 \rho'_2] \\ - (s_{in} - s)\rho'_1[\mu_{1\infty} \lambda_{q_1} - \mu'_1(q_1) x_1 \lambda_{x_1}] - (s_{in} - s)\rho'_2[\mu_{2\infty} \lambda_{q_2} - \mu'_2(q_2) x_2 \lambda_{x_2}]. \qquad (22) $$

**Proof.** By differentiating $\theta_1$ w.r.t. $t$, we have

$$ \dot{\theta}_1 = \theta_2 \dot{s} + \rho'_1[\dot{x}_1 \lambda_s + x_1 \dot{\lambda}_s - \dot{\lambda}_{q_1}] + \rho'_2[\dot{x}_2 \lambda_s + x_2 \dot{\lambda}_s - \dot{\lambda}_{q_2}]. $$

Using (20) and (12), we obtain (21). Using that $\ddot{\tilde{\Phi}} = -\dot{s}\theta_1 + (s_{in} - s)\dot{\theta}_1$, we obtain (22). □

Note that these computations have been verified thanks to a symbolic computation software. The next step is to establish whether Legendre–Clebsch’s condition is verified or not along a singular arc. Recall that this condition is necessary for optimality and that it can be stated as follows (see, e.g., [27,30,31]).

**Theorem 1** (Legendre–Clebsch’s condition [27]). Let $I = [t_1, t_2]$ be such that the trajectory is singular over $[t_1, t_2]$. Then, one has:

$$ \ddot{\tilde{\Phi}}|_{D} \geq 0, \qquad (23) $$

which is fulfilled over $I = [t_1, t_2]$.

Using the expression of the derivative $\ddot{\tilde{\Phi}}$ given in (22), we provide in the next lemma the expression of the second derivative, $\ddot{\tilde{\Phi}}|_D$.

**Lemma 3.** Let $I = [t_1, t_2]$ be such that the trajectory is singular over $[t_1, t_2]$. Then, one has:

$$ \ddot{\tilde{\Phi}}|_D = (s_{in} - s)^2 \theta_2 = (x_1 \lambda_s - \lambda_{q_1}) \frac{(\rho''_1 \rho'_2 - \rho''_2 \rho'_1)}{\rho'_2}. \qquad (24) $$

**Proof.** In (22), the only term involving the control $D$ is related to $\dot{s}$. We obtain (24) using that $\theta_1 \equiv 0$ over $[t_1, t_2]$. □

**Proposition 3.** Along a singular arc that occurs over a time interval $[t_c, t_f]$, it holds that:

$$ \lambda_s = 0, $$
over $[t_c, t_f]$ and Legendre-Clebsch's condition (with a strict inequality) is fulfilled over $[t_c, t_f)$.

**Proof.** Because the trajectory is singular over $[t_c, t_f]$, one has $\theta_1 \equiv 0$ over this interval; thus, $\lambda_s(\cdot)$ satisfies:

$$ \dot{\lambda}_s = D\lambda_s, \quad \lambda_s(t_f) = 0. $$

It follows that $\lambda_s \equiv 0$ over $[t_c, t_f]$. In a left neighborhood of $t = t_f$, one has from (12) $\dot{\lambda}_{q_1} < 0$; thus, since $\lambda_{q_1}$ vanishes at $t = t_f$, we necessarily have $\lambda_{q_1} > 0$ in a left neighborhood of $t_f$. Because $\lambda_s$ is zero over $[t_c, t_f]$, we deduce that $x_1\lambda_s - \lambda_{q_1} = -\lambda_{q_1} < 0$ over $[t_c, t_f)$ (at $t = t_f$, $\lambda_{q_1}$ vanishes at $t_f$, as well as $\ddot{\Phi}_{|D}$). Over $[t_c, t_f]$, we note that

$$ \dot{\lambda}_{x_1} = \lambda_{x_1}(D - \mu_1(q_1)), $$

hence $\lambda_{x_1}$ does not vanish over $[t_c, t_f]$. Suppose that $\lambda_{q_1}$ vanishes over $[t_c, t_f]$ at a time $t' \in [t_c, t_f]$. Then, one must have $\dot{\lambda}_{q_1}(t') \ge 0$ since $\lambda_{q_1} > 0$ over $(t', t_f)$. However, at $t = t'$, the adjoint equation implies that

$$ \dot{\lambda}_{q_1}(t') = -\mu_1'(q_1(t'))x_1(t')\lambda_{x_1}(t') < 0 $$

because $q_1(t') > k_1$, $\mu_1' > 0$ over $(k_1, +\infty)$, $x_1 > 0$, and $\lambda_{x_1} > 0$ over $[t_c, t_f]$. This is a contradiction and thus $\lambda_{q_1}$ does not vanish over $t_c, t_f$. Assumption 1 implies that $\rho_1''\rho_2' - \rho_2''\rho_1' < 0$, thus $\ddot{\Phi}_{|D} > 0$ over $[t_c, t_f)$ as desired. We can then conclude that Legendre-Clebsch's condition (with a strict inequality) is fulfilled over the whole interval $[t_c, t_f)$. $\square$

A consequence of the previous proposition is that, when a singular arc occurs over some time interval $[t_c, t_f]$, then it is of order 1. Based on this proposition, we shall only consider singular arcs of first order in the remaining of the paper. If Legendre-Clebsch's condition holds true, the singular arc is said to be of turnpike type [26]. The expression defining the singular control can then be derived using (22). Next, let $\zeta(X, \lambda)$ be defined by:

$$ \zeta(X, \lambda) := \theta_2[\rho_1 x_1 + \rho_2 x_2] - \lambda_s(x_1\mu_1\rho_1' + x_2\mu_2\rho_2') - \rho_1'[\mu_{1\infty}\lambda_{q_1} - \mu_1'(q_1)x_1\lambda_{x_1}] \\ - \rho_2'[\mu_{2\infty}\lambda_{q_2} - \mu_2'(q_2)x_2\lambda_{x_2}]. \quad (25) $$

We now give an expression of the singular control as a feedback of the state and covector.

**Proposition 4.** Suppose that an extremal is singular over $[t_1, t_2]$ and that (23) is verified over $[t_1, t_2]$ with a strict inequality. Then, the singular control $D_s$ is given by

$$ D_s(X, \lambda) := \frac{\zeta(X, \lambda)}{(s_{in} - s)\theta_2}, \quad (26) $$

where we recall that $\theta_2 = \rho_1''(s)[x_1\lambda_s - \lambda_{q_1}] + \rho_2''(s)[x_2\lambda_s - \lambda_{q_2}]$ and $\zeta$ is given by (25).

**Proof.** This expression follows from (22) in which $s$ is replaced by $(s_{in} - s)D - \rho_1 x_1 - \rho_2 x_2$ and $\theta_1 \equiv 0$ (since $\ddot{\Phi} = \ddot{\Phi} = 0$). $\square$

**Corollary 1.** If the singular arc occurs over some time interval $[t_c, t_f)$, expression (26) simplifies into

$$ D_s(X, \lambda) := \frac{\zeta(X, \lambda)}{(s_{in} - s)\theta_2}, \quad (27) $$

where $\zeta$ is given by

$$ \zeta(X, \lambda) := \theta_2[\rho_1 x_1 + \rho_2 x_2] - \rho_1'[\mu_{1\infty}\lambda_{q_1} - \mu_1'(q_1)x_1\lambda_{x_1}] - \rho_2'[\mu_{2\infty}\lambda_{q_2} - \mu_2'(q_2)x_2\lambda_{x_2}], \\ \text{and, in this case, } \theta_2 \text{ simplifies also into } \theta_2 = -\rho_1''\lambda_{q_1} - \rho_2''\lambda_{q_2} \text{ because } \lambda_s \equiv 0. $$
**Remark 3.**

(i) Note that Legendre–Clebsch’s condition (with a strict inequality) is equivalent to $\theta_2 > 0$ (and that this condition is always verified over $[t_c, t_f]$) with $t_c$ close to $t_f$.

(ii) In view of the general expression giving the singular control, see (26), there is no guarantee a priori that the singular control $D_s$ is always with values in $[0, D_{max}]$, i.e., that the singular arc is always admissible (even if Legendre–Clebsch’s condition is verified). This can bring additional difficulties; however, we may discard this point by choosing $D_{max}$ large enough.

(iii) Notice that (27) is at least active at $t = t_f$ in the case where a singular arc $D_s$ steers the model trajectories towards the target $\mathcal{T}$, since the transversality conditions ensure that $\lambda_s(t_f) = 0$.

## 4.2. About the Occurrence of a Terminal Singular Arc at the Terminal Time

The aim of this section is to discuss the possibility of having a singular arc over some time interval $[t_f - \tau, t_f]$ (with $\tau > 0$) and the structure of optimal controls. Our main questioning is as follows:

Does any optimal trajectory contain a singular arc over some time interval $[t_f - \tau, t_f]$?

To analyze this point, let us summarize properties of the switching function at $t = t_f$ (that are consequences of transversality conditions associated with the codimension 1 target):

* The switching function and its derivative vanish at $t = t_f$:

$$ \dot{\Phi}(t_f) = \Phi(t_f) = 0. \qquad (28) $$

* The second derivative of the switching function satisfies:

$$ \ddot{\Phi}_{|D}(t_f) = 0. \qquad (29) $$

* In addition, Legendre–Clebsch’s condition (23) is always satisfied along a singular arc defined in a left neighborhood of the terminal time $t = t_f$.

The necessary conditions (28) and (29) are a very good indication for the occurrence of a singular arc and are thus strong arguments to answer positively to the above question. Thus, we could now wonder whether or not conditions (28) and (29) are sufficient to ensure the occurrence of a singular arc in some time interval $[t_f - \tau, t_f]$. It appears that this question is complex and falls into the setting of geometric optimal control theory. As far as we know, such conditions are not equivalent to the occurrence of a singular arc over some time interval $[t_f - \tau, t_f]$ (this may depend, in particular, on the initial condition). It is, however, worth mentioning that these conditions (in particular (28)) are commonly used numerically to implement a singular arc in shooting methods [32]. In our context of Droop model, it is very interesting to notice that singular arcs are the cornerstone of the optimal control, in particular for a large set of initial conditions that are biologically meaningful (typically for heterogeneous cultures where $x_1^0/x_2^0 \approx 1$). However, the answer to the above question is not always true and depends on the initial condition (as it has been confirmed using direct optimization methods, see Section 5). Indeed, as illustrated in Example 2—Section 5, when the initial conditions $x_i^0$ are taken very close to the target $(x_2(t_f)/x_1(t_f) = \varepsilon)$, and $\mu_1(q_1^0) - \mu_2(q_2^0) > 0$, the singular arc does not appear or appears marginally at $t = t_f$ to satisfy the transversality conditions.

Recall that $\dot{\Phi}(t_f) = \Phi(t_f) = 0$. Hence, the sign of $\dot{\Phi}$ depends on $\ddot{\Phi}(t_f)$ that is computed in the next lemma.

**Lemma 4.** *At the terminal time, one has*

$$ \mu_1(q_1(t_f)) - \mu_2(q_2(t_f)) > 0. \qquad (30) $$
In addition, the second derivative of the switching function exists at $t = t_f$ and:

$$ \ddot{\Phi}(t_f) = \frac{[s_{in} - s(t_f)][\rho'_2(s(t_f))\mu'_2(q_2(t_f)) - \rho'_1(s(t_f))\mu'_1(q_1(t_f))]}{\mu_2(q_2(t_f)) - \mu_1(q_1(t_f))}. \quad (31) $$

**Proof.** Inequality (30) follows from the expression of $\lambda(t_f)$ and the transversality condition (16). The expression of $\ddot{\Phi}(t_f)$ in (31) follows from (22). $\square$

We can now define the function:

$$ t \mapsto \zeta(t) := \rho'_2(s(t))\mu'_2(q_2(t)) - \rho'_1(s(t))\mu'_1(q_1(t)). $$

From the previous lemma, we deduce the behavior of an optimal path near the terminal time:

* First, the target set can only be reached at some point $X(t_f)$ such that (30) is fulfilled.

* In addition, if $\zeta(t_f) \neq 0$, then, in a left neighborhood of $t = t_f$, the optimal control $D(\cdot)$ is of bang type and satisfies

$$ D(t) = \mathrm{sign}(\zeta(t)). $$

* If a singular arc occurs in a left neighborhood of $t = t_f$, then one must have $\zeta(t_f) = 0$, i.e., a singular arc reaches the target in the subset of $\mathcal{T}$ defined as:

$$ \mathcal{T}' := \{ X = (s, q_1, x_1, q_2, x_2) \in \mathcal{T} ; \mu_1(q_1) - \mu_2(q_2) > 0 \text{ and} \\ \rho'_2(s)\mu'_2(q_2) - \rho'_1(s)\mu'_1(q_1) = 0 \}. $$

### 4.3. Toward an Optimal Synthesis Characterizing the Optimal Solutions

Reducing the number of switching times is in general non-tractable for nonlinear optimal control problems governed by a system in dimension greater than three. Nevertheless, thanks to the properties of singular arcs, we obtained previously and of the switching function at $t = t_f$, we can expect a limited number of possible structures for an optimal control as we formulate in the next conjecture.

**Conjecture 1.** Every initial condition in $\Omega$ is steered optimally to the target set via a control $D$ that has a finite number of switching times. In addition, for almost every initial condition, an optimal control presents the following structure:

$$ B_{\pm}B_{\mp}S. $$

For a large set of initial conditions in some subset $\Omega^\dagger \subset \Omega$ (far from the target), there is a single bang arc and a terminal singular arc, whereas, for some initial conditions close to the target set, no singular arc occurs (i.e., $S$ is of zero duration).

This conjecture has been verified numerically for a large number of initial conditions (see Section 5). Our argumentation to confirm this conjecture is as follows.

* From the PMP, we have seen that, for every initial condition in $\Omega$, an optimal control is a concatenation of bang arcs $B_{\pm}$ and singular arcs $S$.

Moreover, since the switching function $\Phi$ satisfies the strong requirements $\ddot{\Phi}(t_f) = \Phi(t_f) = 0, \ddot{\Phi}|_D(t_f) = 0$

(from the transversality conditions) as well as Legendre-Clebsch's condition, we conjecture that, for almost all initial conditions, an optimal control is singular in a left neighborhood of the terminal time. This implies in particular that the number of switchings is finite since we proved that any terminal singular arc is of first order. In addition, the number of switchings is minimal in general (apart when chattering
occurs, see [27], which is not the case here). Thus, an optimal control should be of type $B_{\pm}B_{\mp}S$ with one or two bang arcs before the terminal singular arc.

As we have seen, we also must have $\xi(t_f) = 0$, which only involves state variables at the terminal time. This surprising condition mixing the first derivative of the basic functions in the Droop model, with the $s$ variable on one side and the $q_i$ variables on the other side, is, however, hard to interpret biologically.

* For some marginal—but admissible—initial conditions outside of $\Omega^\dagger$, the structure of an optimal control of (6) may be of bang type for almost all $t \in [0, t_f)$, or $[0, t_f - \tau]$, with very small $\tau > 0$ (see, e.g., Example 2 in the next section). This is the case when typically $x_1^0 \gg x_2^0$, i.e., the initial condition is very close to the target set $\mathcal{T}$, with in addition $\mu_1(q_1^0) - \mu_2(q_2^0) > 0$. Thus, the requirement $\mu_1(q_1(t_f)) - \mu_2(q_2(t_f)) > 0$ is easily satisfied. Thus, in this particular situation, it comes as no surprise that the fastest path to reach the target $\mathcal{T}$ is the one exploiting the fact that $\dot{x}_1(0) \gg \dot{x}_2(0)$ (since $x_1^0 \gg x_2^0$) along with $D(t) = 0$, since it maximizes $\dot{x}_1(t)$ (we recall that $\dot{x}_i = (\mu_i(q_i) - D(t))x_i$).

It is worth noticing that, when no singular arc occurs, the strategy mainly consists of “pushing” $x_1$ and $x_2$ as quickly as possible towards the target $\mathcal{T}$ using $D = 0$, when the initial conditions $x_1^0$ and $x_2^0$ are very close to $\mathcal{T}$. Nevertheless, this strategy may not be the optimal one whenever $q_1^0$ and $q_2^0$ are “far” from satisfying $\mu_1(q_1(t_f)) - \mu_2(q_2(t_f)) > 0$. This is typically the case illustrated in Example 3 in the next section.

For an optimal control of type $B_{\pm}S$, the occurrence of a singular arc is related to the so-called turnpike phenomenon that we now explain in this framework. For a large subset of initial conditions $S^\dagger \subset \Omega$ that are biologically the most relevant, the structure of the optimal control is bang-singular $B_{\pm}S$. The singular arc is the control $D_s$ given in (27) that reaches the target $\mathcal{T}$. Moreover, this singular phase coincides with optimal trajectories $(s(t), q_1(t), q_2(t))$ that stay most of the time close to the critical point $(s_c, q_{1c}, q_{2c})$ defined in Section 2.3 (related to the actual growth rates and the function $\Delta(s) = \mu_1(\delta_1^{-1}(s)) - \mu_2(\delta_2^{-1}(s)))$. Observe, for instance, the trajectories $s$, $q_1$ and $q_2$ in Figure 3a. We also believe that the concatenation of bang arcs before the major singular phase exclusively aims at moving $(s^0, q_1^0, q_2^0)$ towards $(s_c, q_{1c}, q_{2c})$. Then, the singular arc $D_s$ takes over at a switching-time instant $t = t_s$ and ensures that the associated singular trajectory, denoted $(s_s(t), q_{1s}(t), q_{2s}(t))$, satisfies the so-called turnpike inequality (see, e.g., [33]),

$$ \|s_s(t) - s_c\| + \|q_{1s}(t) - q_{1c}\| + \|q_{2s}(t) - q_{2c}\| \le a_1(e^{-a_2t} + e^{a_2(t-t_f)}), \quad (32) \\ a_1, a_2 > 0, \text{ for all } t \in [t_s, t_f]. $$

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

<figure><img src="image_8.png" /><figcaption>Figure 3. (a) Optimal state in Droop model, (b) the resulting optimal co-state trajectories, which satisfy in particular all the transversality conditions of the PMP.</figcaption></figure>

This is, for instance, the case for optimal controls illustrated in Figure 4 (of type $B_-S$) and in Figure 5a (of type $B_-B_+S$). The inequality (32) usually holds when the time interval
$[0, t_f]$ is not excessively short [33–35], which is the generic case in the Droop model (1) associated with (6). Indeed, in practice, the most significant biological experiments aim to separate species and select $x_1$ starting from an homogeneous culture (a well-balanced initial culture with $x_1^0/x_2^0 \approx 1$) or even from $x_2^0 \gg x_1^0$ with the challenging issue of selecting the minority species ($x_1$), which is not naturally promoted. In these cases, Droop's kinetics ensure that the minimum selection time $t_f$ cannot be excessively short and therefore singular arcs as well as the *turnpike-type* behavior appear systematically in the optimal strategy of (6).

<figure><img src="image_9.png" /><figcaption>**Figure 4.** The optimal control in Example 2 ($s^0 = 4, q_i^0 = 1.9, x_i^0 = 0.3$) is of *bang(0)-singular type*.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5.** (a) The optimal state trajectories $s, q_1$ and $q_2$ (associated with the optimal control $D$ in Figure 6) get closer over time to the critical point ($s_c, q_{1c}, q_{2c}$). The target $\mathcal{T}$, with $\varepsilon = 0.08$, is reached quickly, without resorting to the singular arc. The transversality conditions are satisfied, and in particular $\lambda_s(t_f) = \lambda_{q_1}(t_f) = \lambda_{q_2}(t_f) = 0$, as illustrated in (b).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Evolution of the quantity $\rho_2'(s(t))\mu_2'(q_2(t)) - \rho_1'(s(t))\mu_1'(q_1(t))$ along the optimal trajectories given in Figure 3a.</figcaption></figure>
## 5. Direct Optimization and Numerical Results

In this section, a direct-optimization approach is performed in order to solve (6) and illustrate the different cases discussed in Section 4.3. The numerical direct methods that we use throughout this paper are implemented in Bocop [36] (an optimal control toolbox), and they have the characteristic to transform the studied (6) into a nonlinear programming problem (NLP) in finite-dimension [37], through the discretization step of the control and the state variables [38]. Numerical results are organized as follows:

* An optimal control of type $B_{-}S$ is developed throughout Example 1.

* An optimal control of type $B_{-}$ is developed throughout Example 2.

* An optimal control of type $B_{-}B_{+}S$ is developed throughout Example 3.

In all the numerical examples, we consider the model parameters given in Table 1, with the settings in Table 2.

<table><caption>**Table 1.** The model parameters used in Section 5.</caption><thead><tr><th></th><th>ρ<sub>m<sub>i</sub></sub></th><th>K<sub>i</sub></th><th>μ<sub>i∞</sub></th><th>k<sub>i</sub></th></tr></thead><tbody><tr><td>Species/strain 1</td><td>7</td><td>0.3</td><td>1.7</td><td>1.75</td></tr><tr><td>Species/strain 2</td><td>8</td><td>0.6</td><td>1.8</td><td>1.80</td></tr></tbody></table>

<table><caption>**Table 2.** Model and OCP settings. The contamination rate $\epsilon$ characterizes the target set $\mathcal{T}$.</caption><thead><tr><th>s<sub>in</sub></th><th>Control D</th><th>Contamination Rate ε</th></tr></thead><tbody><tr><td>6</td><td>[0, D<sub>max</sub> = 1.5]</td><td>0.08</td></tr></tbody></table>

Assumption 3 is verified when $D_{max} = 1.5$, namely because $D_{max}$ is precisely chosen above the maximum actual growth rates of the species. The contamination rate in all the examples is fixed to a significantly small value, $\epsilon = 0.08$.

In Bocop, the state variables (and even the time, in minimal-time OCPs) of the Droop model (1) are discretized with a Lobatto scheme based on Runge-Kutta methods of type Lobatto-IIIC of order 6, which uses an implicit trapezoidal rule. The main settings used in Bocop are given in Table 3.

<table><caption>**Table 3.** Bocop settings used in Section 5.</caption><thead><tr><th>Discretization Method</th><th>Lobatto IIIC (Implicit, 4-Stage, Order 6)</th></tr></thead><tbody><tr><td>Time steps</td><td>130</td></tr><tr><td>NLP tolerance</td><td>&lt;10<sup>-14</sup></td></tr></tbody></table>

**Example 2.** In the first example, we consider the Droop model resulting from the parameters in Table 1 and Figure 7, associated with the settings in Tables 2 and 3, and the initial conditions given in Table 4.

<table><caption>**Table 4.** The initial conditions used in Example 2.</caption><thead><tr><th>s<sup>0</sup></th><th>q<sup>0</sup><sub>1</sub></th><th>x<sup>0</sup><sub>1</sub></th><th>q<sup>0</sup><sub>2</sub></th><th>x<sup>0</sup><sub>2</sub></th></tr></thead><tbody><tr><td>4</td><td>1.9</td><td>0.3</td><td>1.9</td><td>0.3</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>Figure 7. Growth and absorption (uptake) functions, respectively $\mu_i$ and $\rho_i$, $i = 1,2$, produced using the model parameters given in Table 1, and used throughout the Examples 1–3. We note that the value of $s$ that maximizes the function $\Delta(s)$ is given $s_c = 0.092$. This point maximizes the difference between the actual growth rates (as discussed in Section 2.3) and defines the unique static solution $(s_c, q_{1c}, q_{2c})$. All the assumptions that ensure the well-posedness of the generic (6) are satisfied in this case. In particular, we highlight that $\Delta(s)$ can be positive and negative, and both species may win the competition using an appropriate control $D$ (see Section 2.4).</figcaption></figure>

In this example, we have a well-balanced initial culture since $x_1^0/x_2^0 = 1$ (see Section 4.3).

The direct optimization method allows us to determine the optimal control $D$, given in Figure 4, that steers the model trajectories towards $\mathcal{T}$ (with $\varepsilon = 0.08$) in minimal-time $t_f = 18.3526$ days.

We check and analyze the evolution of the switching function $\tilde{\Phi}$, its derivatives, and the co-state of the substrate $s$ (Figure 8) in order to characterize the switching instant $t_s \in (0, t_f)$. This time-instant coincides with $\lambda_s = 0$ (since the singular arc is the one reaching the target $\mathcal{T}$), $\tilde{\Phi} = \dot{\tilde{\Phi}} = 0$ (thus activating $D_s$, according to the PMP). We also notice that the condition $\ddot{\tilde{\Phi}}(t_f) = 0$ is also satisfied. The optimal state and co-state trajectories are depicted in Figure 3, where we notice that $s(t), q_1(t)$ and $q_2(t)$ evolve around the static critical point $(s_c, q_{1c}, q_{2c})$ for almost all $t \in [t_s, t_f]$, see the turnpike-like property discussed in Section 4.3.

The optimal control strategy aims to maximize the difference between the actual growth rates (the function $\Delta(s)$) as illustrated in Figure 9. The initial arc bang(0) drives $s^0$ towards $s_c$ (we recall that $\dot{s} = (s_{in} - s)D - \rho_1(s)x_1 - \rho_2(x_2)$, $s^0 = 4$, $s_{in} = 6$, and, $s_c = 0.092$). The quantity $\mu_1(q_1(t)) - \mu_2(q_2(t))$ is maximized, with a delayed-dynamics, as a consequence of maximizing $\Delta(s)$. At the final time of $t_f = 18.3526$ days, we have $\mu_1(q_1(t_f)) - \mu_2(q_2(t_f)) > 0$. Finally, we check in Figure 6 that, at $t = t_f$, we have

$$ \frac{\rho_2'(s)}{\rho_1'(s)}(t_f) = \frac{\mu_1'(q_1)}{\mu_2'(q_2)}(t_f). $$
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Characterization of the switching instant and of the singular arc in Example 2. (a) The switching function and its derivatives. (b) The co-state of the substrate s.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** (a) Evolution of $\Delta(s(t))$ and $\mu_1(q_1(t)) - \mu_2(q_2(t))$ along the optimal trajectories. The optimal control aims to maximize the function $\Delta$. We also check in (b), using the optimal model trajectories given in Figure 3, that $s(t) + q_1(t)x_1(t) + q_2(t)x_2(t)$ converges towards $s_{in}$ (see Section 2.4).</figcaption></figure>

The behavior of the optimal control and optimal trajectories described in Example 2 is definitively the most compelling one (with *bang(0)-singular* or *bang(Dmax)-singular* arcs) from a biological standpoint, since it is the one that systematically appears when the final time is not extremely short. Indeed, in practice, initial conditions start more commonly sufficiently “far” from the target $T$, leading to a final time that allows the singular arc and the *turnpike*-like behaviors to hold.

**Example 3.** Now, let us consider the initial conditions in Table 5.

<table><caption>**Table 5.** The initial conditions used in Example 3.</caption><thead><tr><th>s<sup>0</sup></th><th>q<sub>1</sub><sup>0</sup></th><th>x<sub>1</sub><sup>0</sup></th><th>q<sub>2</sub><sup>0</sup></th><th>x<sub>2</sub><sup>0</sup></th></tr></thead><tbody><tr><td>1</td><td>2.5</td><td>1.2</td><td>2</td><td>0.1</td></tr></tbody></table>

It is worth noticing that the initial conditions in Example 3 are intuitively favourable for reaching the target $T$ in a very short time, since $\mu_1(q_1^0) - \mu_2(q_2^0) > 0$ and $x_2^0/x_1^0 = 0.083$, while $\varepsilon = 0.08$ (very close to the target). The optimal control in this case is mainly a *bang(0)* over time, as illustrated in Figure 10 (see Section 4.3 for more details). The optimal state trajectories and co-state trajectories (that satisfy the transversality conditions) are illustrated in Figure 5.
<figure><img src="image_15.png" /><figcaption>**Figure 10.** The optimal control given in (a) is *bang(0)* almost everywhere over $[0, t_f)$. It achieves species separation (b) and selects the species $x_1$ at $t_f = 1.1810$ days, since $x_1(t_f) = 1.3425$ and $x_2(t_f)/x_1(t_f) = \varepsilon = 0.08$.</figcaption></figure>

**Example 4.** In the last example, let us consider the initial conditions in Table 6.

<table><caption>**Table 6.** The initial conditions used in Example 4.</caption><thead><tr><th>s<sup>0</sup></th><th>q<sub>1</sub><sup>0</sup></th><th>x<sub>1</sub><sup>0</sup></th><th>q<sub>2</sub><sup>0</sup></th><th>x<sub>2</sub><sup>0</sup></th></tr></thead><tbody><tr><td>5</td><td>1.8</td><td>2.2</td><td>5</td><td>0.18</td></tr></tbody></table>

In this situation, we notice that initial conditions are still favourable for reaching the target $\mathcal{T}$ in a very short time because $x_2^0/x_1^0 = 0.0818$ (with $\varepsilon = 0.08$, so $x_i^0$ are very close to the target). However, we also note that $\mu_1(q_1(0)) - \mu_2(q_2(0)) < 0$, while $\mu_1(q_1(t_f)) - \mu_2(q_2(t_f))$ should be positive at $t = t_f$ (see Section 3). Thus, the issue of minimal-time separation is slightly more complex than Example 3, and we obtain a structure for the optimal control of *bang(0)-bang(D<sub>max</sub>)-singular* type, as illustrated in Figure 11. The corresponding model trajectories and optimal co-states are provided in Figure 12.

<figure><img src="image_16.png" /><figcaption>**Figure 11.** The optimal control given in (a) is *bang(0)-bang(D<sub>max</sub>)-singular* over $[0, t_f)$, where $t_f = 4.5541$ days. A first switching instant from *bang(0)* to *bang(D<sub>max</sub>)* occurs around 0.3 days, then $t_s$ occurs when $\lambda_s = 0$, as illustrated in (b), starting the singular arc that steers the model trajectories towards $\mathcal{T}$, with $\varepsilon = 0.08$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 12.** The optimal model trajectories (a), and the optimal co-state trajectories (b) that satisfy the transversality conditions.</figcaption></figure>

**6. Conclusions**

The minimal-time OCP for the competition between two microbial species accumulating nutrients is a key issue. Progressing along this problem will definitely help for experiments that nowadays last more than 6 months [14,15]. However, the competition described by the Droop model turns out to be significantly more complicated than for the Monod model in dimension 2 [17]. We have improved the preliminary results about this problem to be found in [22,29] in order to provide an optimal synthesis depending on the initial condition. In particular, we applied the PMP, we discussed the structure of the optimal control, and we identified the singular arc steering optimal paths to the desired target set in a minimal amount of time. This study also highlights the *turnpike* behavior [33] although an exact verification in our case is not possible in our setting since the problem is affine w.r.t. entries in contrast with [35] that, in general, requires coercive hypotheses on the Hamiltonian w.r.t. entries. As usual in optimal control problems that are affine in the control, the study of singular arcs via geometric methods is a crucial issue.

This study also raised a mathematical (open) question outside the scope of this paper on the existence or not of a terminal singular arc whenever the switching function, its derivative, and second derivative w.r.t. the input vanish on a terminal manifold of codimension 1 (see, e.g., [39]).

Future work will focus on the determination of closed loop (sub-optimal) controllers to be applied for bioreactor control subject to uncertainties that are inherently present in biological systems. Finally, the proposed strategy must now be tested experimentally to assess the gain in experimental time it can offer.

**Author Contributions:** Conceptualization, W.D., T.B. and O.B.; methodology, W.D., T.B. and O.B.; software, W.D.; validation, W.D., T.B. and O.B.; formal analysis, W.D., T.B. and O.B.; investigation, W.D. and T.B.; resources, W.D., T.B. and O.B.; data curation, W.D., T.B. and O.B.; writing—original draft preparation, W.D., T.B. and O.B.; writing—review and editing, W.D., T.B. and O.B.; visualization, W.D., T.B. and O.B.; supervision, W.D., T.B. and O.B.; project administration, W.D., T.B. and O.B.; funding acquisition, W.D., T.B. and O.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors would like to thank the *Bio-MSA* (Graine-ADEME) project, funded by the *Agence Nationale de la Recherche (ANR)* in France, for partially supporting this work.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.
**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors are grateful to J.-B. Caillau and J.-B. Pomet for fruitful exchanges about singular arcs.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Masci, P.; Bernard, O.; Grognard, F. Continuous selection of the fastest growing species in the chemostat. *IFAC Proc. Vol.* **2008**, *41*, 9707–9712. [CrossRef]

2. Angermayr, S.A.; Hellingwerf, K.J.; Lindblad, P.; de Mattos, M.J.T. Energy biotechnology with cyanobacteria. *Curr. Opin. Biotechnol.* **2009**, *20*, 257–263. [CrossRef] [PubMed]

3. Ducat, D.C.; Way, J.C.; Silver, P.A. Engineering cyanobacteria to generate high-value products. *Trends Biotechnol.* **2011**, *29*, 95–103. [CrossRef] [PubMed]

4. Wijffels, R.H.; Barbosa, M.J.; Eppink, M.H. Microalgae for the production of bulk chemicals and biofuels. *Biofuels Bioprod. Biorefin. Innov. Sustain. Econ.* **2010**, *4*, 287–295. [CrossRef]

5. Mata, T.M.; Martins, A.A.; Caetano, N.S. Microalgae for biodiesel production and other applications: A review. *Renew. Sustain. Energy Rev.* **2010**, *14*, 217–232. [CrossRef]

6. Thangavel, M.; Pugazhendhi, A. Utilization of algae for biofuel, bio-products and bio-remediation. *Biocatal. Agric. Biotechnol.* **2019**, *17*, 326–330.

7. Morales, M.; Collet, P.; Lardon, L.; Hélias, A.; Steyer, J.P.; Bernard, O. Life-cycle assessment of microalgal-based biofuel. In *Biofuels from Algae*; Elsevier: Amsterdam, The Netherlands, 2019; Chapter 20, pp. 507–550. ISBN 9780444641922.

8. Cervantes-Gaxiola, M.E.; Hernández-Calderón, O.M.; Rubio-Castro, E.; Ortiz-del-Castillo, J.R.; González-Llanes, M.D.; Rios-Iribe, E.Y. In silico study of the microalgae-bacteria symbiotic system in a stagnant pond. *Comput. Chem. Eng.* **2020**, *135*, 106740. [CrossRef]

9. De-Luca, R.; Trabuio, M.; Barolo, M.; Bezzo, F. Microalgae growth optimization in open ponds with uncertain weather data. *Comput. Chem. Eng.* **2018**, *117*, 410–419. [CrossRef]

10. Malek, A.; Alamoodi, N.; Almansoori, A.S.; Daoutidis, P. Optimal dynamic operation of microalgae cultivation coupled with recovery of flue gas CO₂ and waste heat. *Comput. Chem. Eng.* **2017**, *105*, 317–327. [CrossRef]

11. Nhat, P.V.H.; Ngo, H.H.; Guo, W.S.; Chang, S.W.; Nguyen, D.D.; Nguyen, P.D.; Bui, X.T.; Zhang, X.B.; Guo, J.B. Can algae-based technologies be an affordable green process for biofuel production and wastewater remediation? *Bioresour. Technol.* **2018**, *256*, 491–501. [CrossRef]

12. Smith, H.L.; Waltman, P. *The Theory of the Chemostat, Dynamics of Microbial Competition*; Press of Cambridge University: Cambridge, UK, 1995.

13. Grognard, F.; Masci, P.; Benoît, E.; Bernard, O. Competition between phytoplankton and bacteria: Exclusion and coexistence. *J. Math. Biol.* **2015**, *70*, 959–1006. [CrossRef] [PubMed]

14. Bonnefond, H.; Grimaud, G.; Rumin, J.; Bougaran, G.; Talec, A.; Gachelin, M.; Boutoute, M.; Pruvost, E.; Bernard, O.; Sciandra, A. Continuous selection pressure to improve temperature acclimation of Tisochrysis lutea. *PLoS ONE* **2017**, *12*, e0183547.

15. Gachelin, M.; Boutoute, M.; Carrier, G.; Talec, A.; Pruvost, E.; Guihéneuf, F.; Bernard, O.; Sciandra, A. Enhancing PUFA-rich polar lipids in Tisochrysis lutea using adaptive laboratory evolution (ALE) with oscillating thermal stress. *Appl. Microbiol. Biotechnol.* **2021**, *105*, 301–312. [CrossRef]

16. Nolasco, E.; Vassiliadis, V.S.; Kähm, W.; Adloor, S.D.; Al Ismaili, R.; Conejeros, R.; Espaas, T.; Gangadharan, N.; Mappas, V.; Scott, F.; et al. Optimal control in chemical engineering: Past, present and future. *Comput. Chem. Eng.* **2021**, *155*, 107528. [CrossRef]

17. Bayen, T.; Mairet, F. Optimization of the separation of two species in a chemostat. *Automatica* **2014**, *50*, 1243–1248. [CrossRef]

18. Bayen, T.; Mairet, F. Optimization of strain selection in evolutionary continuous culture. *Int. J. Control* **2017**, *90*, 2748–2759. [CrossRef]

19. Droop, M.R. Vitamin B₁₂ and marine ecology. (IV). The kinetics of uptake growth and inhibition in Monochrysis lutheri. *J. Mar. Biol. Assoc.* **1968**, *48*, 689–733. [CrossRef]

20. Droop, M.R. Some thoughts on nutrient limitation in algae. *J. Phycol.* **1973**, *9*, 264–272. [CrossRef]

21. Djema, W.; Giraldi, L.; Bernard, O. An Optimal Control Strategy Separating Two Species of Microalgae in Photobioreactors. In *Proceedings of the Conference on Dynamics and Control of Process Systems, including Biosystems*—12th DYCOPS, Florianopolis, Brazil, 23–26 April 2019; Volume 52, pp. 910–915.

22. Djema, W.; Bernard, O.; Giraldi, L. Separating Two Species of Microalgae in Photobioreactors in Minimal Time. *J. Process Control (JPC)* **2020**, *87*, 120–129. [CrossRef]

23. Andrés-Martínez, O.; Biegler, L.T.; Flores-Tlacuahuac, A. An indirect approach for singular optimal control problems. *Comput. Chem. Eng.* **2020**, *139*, 106923. [CrossRef]

24. Pontryagin, L.S. *Mathematical Theory of Optimal Processes*; Springer: Berlin/Heidelberg, Germany, 1964.

25. Caperon, J.; Meyer, J. Nitrogen-limited growth of marine phytoplankton—II. Uptake kinetics and their role in nutrient limited growth of phytoplankton. *Deep Sea Res. Oceanogr. Abstr.* **1972**, *19*, 619–632. [CrossRef]
26. Boscain, U.; Piccoli, B. *Optimal Syntheses for Control Systems on 2D Manifolds*; Springer: Berlin/Heidelberg, Germany, 2004; Volume 43.

27. Schättler, H.; Ledzewicz, U. *Geometric Optimal Control*; Springer: New York, NY, USA, 2012.

28. Vinter, R. *Optimal Control, Systems and Control: Foundations and Applications*; Birkhauser: Basel, Switzerland, 2000.

29. Djema, W.; Bernard, O.; Bayen, T. Optimal control separating two microalgae species competing in a chemostat. In Proceedings of the 59th IEEE Conference on Decision and Control (CDC), IEEE Proceedings 2020, Jeju, Korea, 14–18 December 2020; pp. 2368–2373.

30. Bayen, T.; Cots, O. Tangency Property and Prior-Saturation Points in Minimal Time Problems in the Plane. *Acta Appl. Math.* **2020**, *170*, 515–537. [CrossRef]

31. Bonnard, B.; Chyba, M. *Singular Trajectories and Their Role in Control Theory*; Springer: Berlin/Heidelberg, Germany, 2003; Volume 40.

32. Aronna, M.S.; Bonnans, J.F.; Martinon, P. A Shooting Algorithm for Optimal Control Problems with Singular Arcs. *J. Optim. Theory Appl.* **2013**, *158*, 419–459. [CrossRef]

33. Djema, W.; Giraldi, L.; Maslovskaya, S.; Bernard, O. Turnpike features in optimal selection of species represented by quota models. *Automatica* **2021**, *132*, 109804. [CrossRef]

34. Porretta, A.; Zuazua, E. Long time versus steady state optimal control. *SIAM J. Cont. Opti.* **2013**, *51*, 4242–4273. [CrossRef]

35. Trélat, E.; Zuazua, E. The turnpike property in finite-dimensional nonlinear optimal control. *J. Differ. Equ.* **2015**, *258*, 81–114. [CrossRef]

36. Bonnans, F.J.; Giorgi, D.; Grelard, V.; Heymann, B.; Maindrault, S.; Martinon, P.; Tissot, O.; Liu, J. *BOCOP: An Open Source Toolbox for Optimal Control—A Collection of Examples*; Technical Reports, Project-Team Commands; Inria, Saclay: Palaiseau, France, 2017.

37. Biegler, L.T. *Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes*; Series MPS-SIAM on Optimization (Book Number 10); SIAM: Philadelphia, PA, USA, 2010; p. 415.

38. Betts, J.T. *Practical Methods for Optimal Control and Estimation Using Nonlinear Programming*, 2nd ed.; Advances in Design & Control; SIAM: Philadelphia, PA, USA, 2010; Volume 19, p. 427.

39. Bonnard, B.; Rouot, J. *Towards Geometric Time Minimal Control without Legendre Condition and with Multiple Singular Extremals for Chemical Networks*; Applied Mathematics, Advances in Nonlinear Biological Systems: Modeling and Optimal Control; American Institute of Mathematical Sciences: San Jose, CA, USA, 2020; pp. 1–34.