# Optimal Separation of Two Microbial Species Competing for Two Substitutable Resources. Species Selection in Minimum-Time

Walid Djema, Térence Bayen, Jean-Luc Gouzé

► To cite this version:

Walid Djema, Térence Bayen, Jean-Luc Gouzé. Optimal Separation of Two Microbial Species Competing for Two Substitutable Resources. Species Selection in Minimum-Time. Journal of Optimization Theory and Applications, 2025, 205 (3), pp.61. 10.1007/s10957-025-02672-x. hal-04984791

HAL Id: hal-04984791

https://inria.hal.science/hal-04984791v1

Submitted on 10 Mar 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal Separation of Two Microbial Species Competing for Two Substitutable Resources
## Species Selection in Minimum-Time

Walid Djema<fnref n="1" />,\*

Térence Bayen<fnref n="2" />

Jean-Luc Gouzé<fnref n="3" />

March 9, 2025

## Abstract

We address a selection problem between two microbial species by controlling the dilution rate in a continuous stirred tank reactor. Unlike earlier studies, we focus on the competition between species for two substitutable substrates instead of a single limiting resource, with each species capable of dominating the culture. The selection of the most suitable strain is formulated as a minimum-time optimal control problem. After examining the reachability of the target set, we conduct a thorough analysis of optimal solutions and singular arcs, using Pontryagin's Maximum Principle and geometric control theory. Notably, unlike studies involving a single substrate, the target is not always reachable, and we construct a counterexample to prove it. In the case where the target can be reached, we highlight that optimal solutions exhibit turnpike-like phenomena by combining the approach with a direct optimization method, and these properties differ from those in the single-substrate case. Although there is no natural static optimization problem associated with a minimum-time control problem, our analysis enables us to characterize the turnpike point as the solution to a static optimization problem associated with an auxiliary Lagrange-type optimal control problem. These results are extensively illustrated throughout this work using direct optimization.

**Keywords:** Pontryagin's Maximum Principle; Singular arcs; Turnpike phenomenon; Chemostat system; Strain selection; Microbial competition.

## 1 Introduction

Microorganisms are essential for life; for instance, cyanobacteria in oceans are responsible for fixing 50% of anthropogenic CO₂ and form the basis of the aquatic food chain. They also play a crucial role in numerous biotechnological applications, such as wastewater treatment, production of high-value pharmaceutical products, and biofuel extraction. Modeling and analyzing these organisms are essential for understanding their behavior, achieving effective control, and optimizing yields [21, 31]. However, the complexity of manipulating these living organisms, which are difficult to manage, limits their use in biotechnological applications. Additionally, the microbial species used are often wild strains not specifically adapted for industrial use, generally resulting in unsatisfactory yields, which can be improved by

<footnote n="1">¹Inria, Université Côte d'Azur (UniCA), GreenOwl team, Centre Inria d'Université Côte d'Azur, Valbonne. France.
walid.djema@inria.fr</footnote>

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

<footnote n="3">³Inria, Université Côte d'Azur (UniCA), Centre Inria d'Université Côte d'Azur, Macbes team, Valbonne. France.
jean-luc.gouze@inria.fr</footnote>

<footnote>\*Corresponding author: walid.djema@inria.fr</footnote>
using more efficient microorganisms. For instance, in algae biofuel production, selecting more transparent microalgae enhances light penetration, photosynthesis, and growth. Additionally, choosing species with higher lipid content can significantly boost biofuel extraction [26]. Consequently, the concept of selecting advantageous species in continuous cultures is gaining traction (see, e.g., [10]). This paradigm extends beyond microalgae and encompasses a wide range of microbial populations used in biotechnology. Mathematical optimization and optimal control provide powerful tools to tackle the challenge of species selection, having already proven effective in various biological and biotechnological contexts. The selection of the most suitable microbial species in practice relies on implementing control laws to drive evolution and optimize specific criteria. Optimal control theory provides a rigorous framework for understanding species dynamics and optimizing microbial selection in biotechnological applications.

This work aims to address the issue of selecting the species with the best traits between two competitors to enhance yields in a given biotechnological application. For that, we consider, in a continuous culture system (chemostat), two microbial species initially coexist, each with the potential to dominate the culture over time. By adjusting the bioreactor's dilution rate (the control variable), we aim to quickly select the most advantageous species. Thus, we aim to select the desirable microbial species competing with another in minimum time, which defines our objective function. This means that, by the end of the culture, the ratio between the concentrations of the two species must reach a prescribed value, with the goal of achieving this state as quickly as possible using an admissible dilution-control strategy. Similar issues have been investigated in previous works, such as [4, 5, 19, 25], where species compete over a single limiting substrate. In [16, 18, 19], a single limiting substrate is also considered, but additional dynamics are included to model the storage of that substrate within each species, known as the Droop model.

The main novelty of the model studied in this work is that it takes into account two limiting substrates instead of a single resource. Specifically, the competing species grow on substitutable nutrient resources that are supplied to the culture. This change in competition rules can enhance growth opportunities for both species and make it more challenging to eliminate the undesirable species due to their access to multiple subsistence resources. More precisely, we consider the evolution of two microbial species in a continuous stirred tank reactor, growing on two different but substitutable nutrients (for more details about substitutable and complementary resources, see [24]). Following the Monod model, originally designed for bacterial cultures (see, e.g., [24, 31]), our central objective is to investigate the problem of species separation or selection—an issue addressed through geometric optimal control. In broader biotechnological contexts, several other relevant models and applications are worth mentioning. For instance, in [20], the competition between two bacterial strains in co-culture is regulated to stabilize their relative abundance. Additionally, the study in [1] explores the asymptotic behavior of a distinct model involving two competing species growing on a limited substrate, which exists in two forms: insoluble and soluble.

Our framework leads us to study a minimum-time control problem governed by a mono-input affine control system with respect to the input. More precisely, the system involves four states (two species and two nutrient sources) and one positive and bounded control (the dilution rate), with the target set to be optimally reached defined as a half-space in Euclidean space. The target set is primarily related to the ratio between the two species. Throughout the paper, we aim to fully characterize the best admissible control strategy that steers the model trajectories from initial conditions (in the four-dimensional space) to a final condition defined as a targeted ratio between the two microbial species. The time necessary to reach the target, when it is possible to do so, is what we seek to minimize. The study of the minimum-time problem is divided into several parts and we now outline the structure of the paper, highlighting the key points and the results obtained.

# 2 Key Findings and Paper Structure

## 1. Reachability of the target set.
In Section 3, we state the optimal control problem (OCP) and analyze the reachability of the target set. We ensure the well-posedness of the studied OCP by considering generic model parameters such that each species can dominate the culture. Then, we prove that under additional conditions on the model parameters, the target set can be reached in finite time from any initial
condition, ensuring the existence of optimal controls. When these assumptions on the model parameters are not satisfied, the reachability of the target cannot be guaranteed. In particular, we highlight a case where the target is not reachable by providing a counterexample, achieved through a careful analysis of the underlying dynamical system using slow-fast dynamical features. These properties indicate that the system is not generally controllable, which will impact optimal trajectories (see point 4 below). In summary, the study of reachability in the model involving two substrates is much more complex, whereas for a single nutrient source, we have proven the attainability of the target from all initial conditions.

**2. Necessary optimality conditions.** To study the behavior of optimal trajectories when the target is reachable, we apply the Pontryagin's Maximum Principle (PMP) [29], in Section 4 to the minimum-time problem. We discuss the occurrence of normal and abnormal solutions ([11]), analyze the transversality conditions, and then, based on the Hamiltonian maximization condition, deduce that the optimal control may consist of a concatenation of bang and singular arcs (see, e.g., [9]). In our case, only first-order singular arcs occur, thus excluding chattering phenomena ([30, 33]). While chattering can arise in certain biotechnological applications (see, e.g., [12, 22]), it is not observed in the species selection problem. Finally, we derive the analytical expression of the first-order singular control using the second derivative of the switching function with respect to time. Additionally, we express the general form of the second-order (Legendre-Clebsch) conditions. Note that the minimum-time synthesis is generally very complex to obtain when dealing with a non-linear control system in dimensions three or higher. Therefore, we combine our results with direct approaches based on direct optimization.

**3. Direct optimization method.** In Section 5, we complement our theoretical study by implementing a direct optimization method using bocop (see [8]). This software relies on the solver IPOPT (Interior Point OPTimizer) to numerically solve a nonlinear programming problem (NLP ([6, 7])) resulting from the discretization of the minimum-time problem. The direct optimization indicates that optimal controls are a concatenation of one or two bang arcs together with a singular arc for which the state and costates are nearly constant. Using the expression of the singular control obtained in Section 4.3, we corroborate the numerical solution with the optimality conditions derived in Section 4. We numerically verify that the numerical singular control coincides with the first-order singular control computed theoretically, and we also verify that the Legendre-Clebsch necessary condition is fulfilled along the singular arc, confirming the optimality of the singular arc. The last part of the trajectory thus corresponds to a turnpike-type behavior, affecting the dynamics of the nutrients, although it is less evident in the biomass of the species. It is important to note that the observed behavior differs from the classic definition of the turnpike feature (in the spirit of [32]), primarily because the turnpike feature is originally defined in problems of Lagrange type (i.e., for which there is an explicit objective function to optimize). The basic assumption is that the turnpike appears when the final time is sufficiently large, which contradicts the minimum-time problem that seeks to minimize the final time. This phenomenon deserves further study, which we conduct in the next section.

**4. Analysis of the turnpike-like phenomenon.** We conclude the paper with a deeper investigation into the turnpike-like phenomenon (Section 6.1). Our objective is to characterize the values around which the optimal solution remains close for most of the time until reaching the target. Specifically, the singular arc approaches particular values in the nutrient variable plane and stays nearly constant. Strong evidence leads us to associate this steady-state solution with some static optimization problem. However, unlike an optimal control problem (OCP) of Lagrange type, there is no obvious way to associate a static optimal control problem with a minimum-time problem (since the cost function is constantly equal to one). Our initial observations show that, for a given fixed set of parameters where the target is reachable, the turnpike point does not depend on initial conditions but is strongly influenced by the yield values, which correspond to the conversion rates between nutrients and microbial biomass. More importantly, we propose an approach that establishes a link between the minimum-time OCP and an auxiliary OCP in the Lagrange form. This makes it possible to define a static-OCP within the auxiliary-OCP framework, for which we determine the solution. The results from the direct optimization method show that the singular arc of the minimum-time problem is significantly close to the solution of the static-OCP associated with the auxiliary OCP. These results are demonstrated in a specific case where we determine the static-OCP for a particular equilibrium point of the species, in which the second species vanishes. Additionally,
it is important to compare the solutions and turnpike-like behavior in problems involving one or two nutrients. In the second case, which is the focus of this work, the turnpike point is not simply the expected maximization of the differences between the growth rates, as it is with one substrate (see, [19, 4, 18]), due to controllability issues. Ultimately, the study of the two-substrate case proved to be more challenging but also much more intriguing than the single-substrate case.

# 3 Model Presentation and Reachability of the Target Set

## 3.1 Presentation of the Control System

Throughout this paper, we consider the evolution of two microbial species in a continuous stirred tank reactor that grow on two different but substitutable nutrients. Following the Monod model, initially designed for bacterial cultures (see, e.g., [24, 31]), this leads to the four-dimensional control system,

$$ egin{cases} 
\dot{m} = (m_{in} - m) D(t) - \frac{1}{Y_{1m}} \mu_{1m} (m) x_1 - \frac{1}{Y_{2m}} \mu_{2m} (m) x_2, \\ 
\dot{n} = (n_{in} - n) D(t) - \frac{1}{Y_{1n}} \mu_{1n} (n) x_1 - \frac{1}{Y_{2n}} \mu_{2n} (n) x_2, \\ 
\dot{x}_1 = [\mu_{1m}(m) + \mu_{1n}(n) - D(t)] x_1, \\ 
\dot{x}_2 = [\mu_{2m}(m) + \mu_{2n}(n) - D(t)] x_2. 
\end{cases} 
\quad (1) $$

In this dynamical system, $m$, $n$ and $x_1$, $x_2$ represent respectively the substrate concentrations and the species concentrations. Additionally, $m_{in}$ and $n_{in}$ are positive parameters and denote the input substrate concentrations while $Y_{1m}$, $Y_{2m}$, $Y_{1n}$, and $Y_{2n}$ are positive and represent the yield coefficients. The dilution rate is $D(\cdot)$ and it is a control variable in $L^\infty(\mathbb{R}_+, [0, D_{max}])$ where $D_{max}$ is the maximum feeding rate. The functions $\mu_{jm}$ and $\mu_{jn}$ for $j = 1, 2$ are the kinetics and they are supposed to be of Monod type (see [28]):

$$ \mu_{jm}(m) := \frac{\mu_{jm}^{\max} m}{K_{jm} + m} \quad ; \quad \mu_{jn}(n) := \frac{\mu_{jn}^{\max} n}{K_{jn} + n}, \quad (2) $$

where $\mu_{jm}^{\max}$, $\mu_{jn}^{\max}$, $K_{jm}$ and $K_{jn}$ are positive parameters. Next, let $R := [0, m_{in}] \times [0, n_{in}]$ and let $\Delta : R \to \mathbb{R}$ be defined as

$$ \Delta(m, n) := \mu_{1m}(m) - \mu_{2m}(m) + \mu_{1n}(n) - \mu_{2n}(n). $$

The function $\Delta$ will play an important role in this paper. We also define the combined growth rates of the species as

$$ \mu_1(m, n) := \mu_{1m}(m) + \mu_{1n}(n); \quad \mu_2(m, n) := \mu_{2m}(m) + \mu_{2n}(n). $$

In the sequel, we suppose that the following hypotheses are verified. The function $\Delta$ takes both positive and negative values over $R$ and it has a unique maximum $(m^*, n^*)$ in the interior of $R$. The maximum dilution rate $D_{max}$ satisfies

$$ \forall (m, n) \in R, \quad D_{max} > \max(\mu_1(m, n), \mu_2(m, n)). $$

These hypotheses will be used to study the reachability of the target set (see Remark 3.2) and to ensure that the optimal control problem is properly posed.

The functions $\mu_1$, $\mu_2$, and $\Delta$, together with the point $(m^*, n^*)$, are depicted in Fig. 1. Parameter values are given in Table 1.

We can now proceed to the formulation of the optimal control problem of interest.

## 3.2 The Minimum-Time Control Problem

Our aim is to drive the system to a state where one of the two species (the one of interest) is dominant over the other. Without any loss of generality, we choose to select the first microbial species $x_1$. So, at
<figure><img src="image_2.png" /><figcaption>Figure 1: Plot of the growth functions $\mu_j$ for $j = 1, 2$ (model parameters are specified in Table 1) which are also used in Examples 1 and 2 (Section 5). Observe that $\Delta$ has a unique maximum $(m^*, n^*) \approx (0.06, 0.23)$ and that it is with positive and negative values (Assumption 3.1).</figcaption></figure>

the end of the process, we wish for the ratio $\frac{x_1}{x_2}$ to be significantly large, that is,

$$ \frac{x_1}{x_2} \gg 1, $$

and our aim is to fulfill this inequality in a minimum amount of time. Since $R$ is forward invariant under the dynamics satisfied by $(m, n)$, initial conditions $X^0 = (m_0, n_0, x_1^0, x_2^0) \in \mathbb{R}^4$ will be chosen in the set

$$ \Omega := R \times (0, +\infty) \times (0, +\infty) \subset \mathbb{R}^4. $$

Similar to the single-substrate case ([4, 17, 18]), we define the target set with respect to the two species as the half-space,

$$ \mathcal{T}_{\varepsilon} := \{ X := (m, n, x_1, x_2) \in \Omega ; x_2 \le \varepsilon x_1 \}. $$

The minimum-time control problem then reads as follows:

$$ \inf_{D(\cdot) \in [0, D_{\max}]} t_f \quad \text{s.t.} \ X(t_f) \in \mathcal{T}_{\varepsilon}, \qquad (\text{OCP}) $$

where $X(\cdot)$ is the unique solution of (1) associated with the control $D(\cdot)$ and such that $X(0) = X^0 \in \Omega \setminus \mathcal{T}_{\varepsilon}$. The existence of an optimal control of (OCP) relies on the reachability of the target that will be studied in Section 3.3. We now wish to introduce an auxiliary optimal control problem that will be useful to understand the optimal synthesis together with the turnpike phenomenon (see Section 6.1). By integrating (1), we have for every $T > 0$,

$$ \ln \left( \frac{x_1(T)}{x_2(T)} \right) - \ln \left( \frac{x_1(0)}{x_2(0)} \right) = \int_0^T \Delta(m(t), n(t)) \, dt. \qquad (3) $$
Given an arbitrary $T > 0$, we introduce the following Lagrange optimal control problem (in the spirit of [17]):

$$ \sup_{D(\cdot) \in [0, D_{\max}]} J_T(D) := \int_0^T \Delta(m(t), n(t)) \, dt, \qquad (\text{OCP}_T) $$

where $X(\cdot) = (m(\cdot), n(\cdot), x_1(\cdot), x_2(\cdot))$ is the unique solution of (1) over $[0, T]$ associated with the control $D(\cdot)$ and such that $X(0) = X^0$. Suppose that (OCP) has a solution and let $(X(\cdot), D(\cdot))$ be an optimal pair such that $D(\cdot)$ steers (1) from $X^0$ to the target set in a time $T$. Then, $D(\cdot)$ is a solution to (OCP$_T$). By contradiction, suppose that there is an admissible control $\tilde{D}$ associated with $\tilde{X}(\cdot) := (\tilde{m}(\cdot), \tilde{n}(\cdot), \tilde{x}_1(\cdot), \tilde{x}_2(\cdot))$ such that $J_T(\tilde{D}) < J_T(D)$. Using (3), we get

$$ \frac{\tilde{x}_1(T)}{\tilde{x}_2(T)} > \frac{x_1(T)}{x_2(T)} = \frac{1}{\varepsilon}. $$

The preceding inequality implies that $\tilde{x}_2(T) < \varepsilon\tilde{x}_1(T)$, so, $\tilde{X}(T)$ belongs to the interior of the target set $\mathcal{T}_\varepsilon$ at time $T$. Since $X_0 \notin \mathcal{T}_\varepsilon$, we deduce that $\tilde{X}$ enters the target set at some time strictly less than $T$. Following [3], let $\tilde{T}$ be the first entry time of $\tilde{X}$ into the target set so that $\tilde{X}(\tilde{T}) \in \mathcal{T}_\varepsilon$. Clearly, $\tilde{X}$ is admissible for (OCP), moreover $\tilde{T} < T$ which contradicts the optimality of the pair $(X(\cdot), D(\cdot))$. This ends the proof. (i) The preceding lemma suggests that a solution to (OCP) may take advantage to reaching the point $(m^*, n^*)$ for which $\Delta$ is maximal over $R$. However, we shall see in Section 6.1 that the situation is more intricate and that optimal solutions of (OCP) will involve a turnpike-like phenomenon for the sub-system satisfied by $(m, n)$ near another point.

(ii) The fact that $\Delta$ changes of sign (see Assumption 3.1) does not necessarily imply that the target set is reachable from $X^0$, but, it indicates that $t \mapsto \ln(\frac{x_1(t)}{x_2(t)})$ is non necessarily strictly monotone (as it would be the case if $\Delta$ is of constant sign) so that it is not possible to favor one species over the other for every admissible control.

(iii) Assumption 3.1 implies that if $D = 0$, then $\dot{x}_i > 0$ for $i = 1, 2$ whereas if $D = D_{\max}$, then $\dot{x}_i < 0$ for $i = 1, 2$. This property ensures that each species can be controlled, however, it is not sufficient to ensure reachability of the target set as we will see in the next section.

### 3.3 Sufficient Conditions for the Reachability of the Target

In this section, we study the reachability of the target set. We shall give sufficient conditions on the data ensuring that it can be reached from any initial condition which will guarantee the existence of an optimal control for (OCP). For convenience, we set:

$$ \alpha(m) := \mu_{1m}(m) - \mu_{2m}(m); \quad \beta(n) := \mu_{1n}(n) - \mu_{2n}(n), $$

so that $\Delta(m, n) = \alpha(m) + \beta(n)$. Now, we state an additional assumption. There are $\bar{m} > 0$ and $\kappa_1 > 0$ such that

$$ \forall n \in [0, n_{in}], \quad \alpha(\bar{m}) + \beta(n) \geq \kappa_1, \qquad (4) $$

or there are $\bar{n} > 0$ and $\kappa_2 > 0$ such that:

$$ \forall m \in [0, m_{in}], \quad \alpha(m) + \beta(\bar{n}) \geq \kappa_2. \qquad (5) $$

This assumption is illustrated on Fig. 2. Roughly speaking, it means that the set $\{\Delta > 0\}$ contains either a vertical or a horizontal line segment.

We now prove that under the previous hypothesis, every initial condition $X_0 \in \Omega$ can be steered to the target set $\mathcal{T}_\varepsilon$ in finite time. Doing so, given $\bar{m} \in (0, m_{in})$, we define a function $\psi_{\bar{m}} : \mathbb{R}_+ \times \mathbb{R}_+ \times (0, m_{in}) \to \mathbb{R}$ as,

$$ \psi_{\bar{m}}(x_1, x_2, m) := \frac{\frac{1}{Y_1}\mu_{1m}(\bar{m})x_1 + \frac{1}{Y_2}\mu_{2m}(\bar{m})x_2}{m_{in} - m}. $$

We have the following property related to the stabilization of (1) using the preceding feedback control law. Consider the closed-loop system (1) with $D = \psi_{\bar{m}}(x_1, x_2, m)$. Then, for every initial condition
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of Assumption 3.3 under parameters given in Table 1 (Fig. (a)) and given by $\mu_{1n}^{\max} = 2$, $\mu_{1m}^{\max} = 3$, $K_{1n} = 0.2$, $K_{1m} = 0.35$, $\mu_{2n}^{\max} = 2$, $\mu_{2m}^{\max} = 3$, $K_{2n} = 0.35$, $K_{2m} = 0.25$ (Fig. (b)).</figcaption></figure>

$(m_0, n_0, x_1^0, x_2^0) \in \Omega, \lim_{t \to +\infty} m(t) = \bar{m}.$ If we set $w := m - \bar{m}$, we find that $w$ satisfies $\dot{w}(t) = -f(t)w(t)$ for all $t \ge 0$ where,

$$f(t) := \frac{\mu_{1m}^{\max} K_{1m} x_1(t)}{Y_{1m}(K_{1m} + \bar{m})(K_{1m} + m(t))} + \frac{\mu_{2m}^{\max} K_{2m} x_2(t)}{Y_{2m}(K_{2m} + \bar{m})(K_{2m} + m(t))}.$$

The result then follows using Barbalat's Lemma, classical reasonings related to dynamical systems, and considering separately the cases $w(0) > 0$ and $w(0) < 0$. The same property also holds with the substrate $n$ in place of $m$ and it is as follows: given $\bar{n} \in (0, n_{in})$, we define in the same way a function $\phi_{\bar{n}}: \mathbb{R}_+ \times \mathbb{R}_+ \times (0, n_{in}) \to \mathbb{R}$ as,

$$\phi_{\bar{n}}(x_1, x_2, n) := \frac{\frac{1}{Y_{1n}}\mu_{1n}(\bar{n})x_1 + \frac{1}{Y_{2n}}\mu_{2n}(\bar{n})x_2}{n_{in} - n}.$$

Consider the closed-loop system (1) with $D = \phi_{\bar{n}}(x_1, x_2, n)$. Then, for every initial condition $(m_0, n_0, x_1^0, x_2^0) \in \Omega, \lim_{t \to +\infty} n(t) = \bar{n}$. The proof being similar, it is omitted.

It is worth noticing that given some initial condition $X^0 \in \Omega$, then $(x_1, x_2)$ is not necessarily bounded over $\mathbb{R}_+$, which means that both feedback controls $\psi_{\bar{m}}$ and $\phi_{\bar{n}}$ may be unbounded. However, following Assumption 3.1, given a fixed time horizon $T > 0$, then, it is always possible to choose $D_{\max}$ large enough in order to ensure that $\psi_{\bar{m}}$ and $\phi_{\bar{n}}$ are always with values within $[0, D_{\max}]$ in such a way that such feedback controls are always admissible over a fixed time period $[0, T]$. Our main result in this section is to prove that the target set can be reached from any initial condition in $\Omega$ provided that Assumption 3.3 is fulfilled. Under Assumption 3.3, the target set can be reached in finite time from any initial condition in $\Omega$. Suppose first that (4) is fulfilled. Then, Lemma 3.3 implies that the unique solution to (1) with $D = \psi_{\bar{m}}(x_1, x_2, m)$ such that $X(0) = X_0$ satisfies $m(t) \to \bar{m}$ as $t \to +\infty$. This implies that there exists $t_0 \ge 0$ such that for every $t \ge t_0$, one has

$$\alpha(m(t)) \ge \alpha(\bar{m}) - \frac{\kappa_1}{2}.$$

If we set $z := \ln(x_1/x_2)$, we have for $t \ge t_0$:

$$\dot{z}(t) = \alpha(m(t)) + \beta(n(t)) \ge \alpha(\bar{m}) - \frac{\kappa_1}{2} + \beta(n(t)).$$
Assumption 3.3 thus implies that for $t \ge t_0$, one has $\dot{z}(t) \ge \frac{\kappa_1}{2} > 0$. We deduce that $z(t) \to +\infty$ as $t \to +\infty$ or, equivalently, that $\frac{x_1(t)}{x_2(t)} \to +\infty$ as $t \to +\infty$. Hence, the target set is necessarily reached in a finite horizon. The same proof can be done with $\phi_{\bar{n}}$ in place of $\psi_{\bar{m}}$ if (5) is fulfilled in place of (4).

Note that such a strategy is suboptimal for reaching the target set, that is why, we investigate in Section 4 the synthesis of an optimal control reaching the target set in minimum time when it is reachable.

## 3.4 Designing a Counterexample Demonstrating Target Inaccessibility

In this part, we point out that it is not always possible for the control system (1) to reach the target set (which is a set of co-dimension 1 in $\mathbb{R}^4$) from certain initial conditions. This is due to the lack of controllability of the system (with a single input and a four-dimensional state). Next, we wish to explain how the parameters defining (1) are chosen in order to highlight the non-reachability of the target set.

### 3.4.1 Slow-fast Dynamics Hypothesis for Model Reduction

The first step for constructing the counterexample is to choose parameters such that the system has time-scale separation properties, and can be reduced to a two-dimensional system, applying the Tikhonov Theorem for slow-fast systems (see, [23, Thm. 11.1]). For this purpose, we choose yield coefficients $Y_{ij}$ small enough, of the form $Y_{ij} = \eta Y'_{ij}$, where $\eta$ is small with respect to the other parameters and $Y'_{ij} > 0$. Thanks to this change of variables, $m$ and $n$ become fast variables, and we can apply the Tikhonov Theorem taking $\eta$ small enough (for conciseness, we have omitted the details regarding the verification of the hypotheses of Tikhonov's Theorem in our setting). This result is also valid for singular perturbations over an infinite interval of time (see, [23, Thm. 11.2]). It follows that the states $m$ and $n$ rapidly converge to their quasi-steady states values $m^\ddagger$ and $n^\ddagger$, which are significantly smaller than $m_{\text{in}}$, $n_{\text{in}}$, and other parameters. The algebraic equations for $m^\ddagger$ and $n^\ddagger$ are given by:

$$m^{\ddagger} = \frac{m_{\text{in}} D}{\left( \frac{1}{Y_{1m}} \frac{\mu_{1m}^{\text{max}}}{K_{1m}} x_1 + \frac{1}{Y_{2m}} \frac{\mu_{2m}^{\text{max}}}{K_{2m}} x_2 \right)}, \quad n^{\ddagger} = \frac{n_{\text{in}} D}{\left( \frac{1}{Y_{1n}} \frac{\mu_{1n}^{\text{max}}}{K_{1n}} x_1 + \frac{1}{Y_{2n}} \frac{\mu_{2n}^{\text{max}}}{K_{2n}} x_2 \right)}$$

We are now interested in the reduced approximate system obtained using the quasi-steady states $m \approx m^\ddagger$ and $n \approx n^\ddagger$ and substituting their values in the slow variables $(x_1, x_2)$. We denote the variables of the reduced system by $(z_1, z_2)$ and obtain a two-dimensional driftless mono-input affine control system $\dot{z}(t) = D(t)F(z(t))$ given by:

$$\begin{aligned} \dot{z}_1 &= D \left[ \frac{\mu_{1m}^{\text{max}}}{K_{1m}} \frac{m_{\text{in}}}{\left( \frac{1}{Y_{1m}} \frac{\mu_{1m}^{\text{max}}}{K_{1m}} z_1 + \frac{1}{Y_{2m}} \frac{\mu_{2m}^{\text{max}}}{K_{2m}} z_2 \right)} + \frac{\mu_{1n}^{\text{max}}}{K_{1n}} \frac{n_{\text{in}}}{\left( \frac{1}{Y_{1n}} \frac{\mu_{1n}^{\text{max}}}{K_{1n}} z_1 + \frac{1}{Y_{2n}} \frac{\mu_{2n}^{\text{max}}}{K_{2n}} z_2 \right)} - 1 \right] z_1, \\ \dot{z}_2 &= D \left[ \frac{\mu_{2m}^{\text{max}}}{K_{2m}} \frac{m_{\text{in}}}{\left( \frac{1}{Y_{1m}} \frac{\mu_{1m}^{\text{max}}}{K_{1m}} z_1 + \frac{1}{Y_{2m}} \frac{\mu_{2m}^{\text{max}}}{K_{2m}} z_2 \right)} + \frac{\mu_{2n}^{\text{max}}}{K_{2n}} \frac{n_{\text{in}}}{\left( \frac{1}{Y_{1n}} \frac{\mu_{1n}^{\text{max}}}{K_{1n}} z_1 + \frac{1}{Y_{2n}} \frac{\mu_{2n}^{\text{max}}}{K_{2n}} z_2 \right)} - 1 \right] z_2, \end{aligned} \quad (6)$$

where $D(\cdot) \in [0, D_{\text{max}}]$ is the control function. Given $z^0 \in (0, +\infty)^2$, we denote by $\Gamma_+(z^0)$ the positive semi-orbit of (6) under a constant control in $(0, D_{\text{max}})$. Actually, this set should depend on the value of the control, but, it is easily seen that $\Gamma_+(z^0)$ is independent of $D(\cdot)$ (indeed, one can proceed to a time scaling $t' = \int_0^t D(\tau) d\tau$ in (6) given any admissible control $D(\cdot)$ such that $D(t) > 0$ for a.e. $t \ge 0$). This justifies the notation $\Gamma_+(z^0)$.

The system (6) is not controllable in the following sense: for every pair of points $z^0, z^1 \in (0, +\infty)^2$ such that $z^1 \notin \Gamma_+(z^0)$, then it is not possible to steer (6) from $z^0$ to $z^1$ in finite time.

By the preceding remark on the time scaling, the reachable set from $z^0$ is $\Gamma^+(z^0)$, which concludes the proof.

Since (6) is a driftless system with a single input, the control $D(\cdot)$ cannot change the geometric phase portrait of the system. It can only introduce new equilibria (when $D(t) = 0$ for every $t \ge 0$), but semi-orbits in the $z$-plane remain the same.
### 3.4.2 Dynamic Behavior of the Reduced Model: Phase Portrait

Note that for a sufficiently small contamination rate $0 < \varepsilon \ll 1$, reaching the target means that $z_2$ must be very small with respect to $z_1$. Additionally, the reduced system admits $\bar{z}_2 = 0$ as a steady-state, for which either the steady-state $\bar{z}_1$ is also zero or $\bar{z}_1$ is given by $\bar{z}_1 = Y_{1m}m_{in} + Y_{1n}n_{in}$. From the z-dynamics, we deduce that the trivial steady-state $(\bar{z}_1, \bar{z}_2) = (0, 0)$ is clearly non-attractive. We can also make the equilibrium point

$$ (\bar{z}_1, \bar{z}_2) = (Y_{1m}m_{in} + Y_{1n}n_{in}, 0) $$

unstable, under certain conditions on the model parameters, as follows. The approximated dynamics $(\tilde{z}_1, \tilde{z}_2)$ of the reduced system $(z_1, z_2)$ around the steady state $(\bar{z}_1, \bar{z}_2) = (Y_{1m}m_{in} + Y_{1n}n_{in}, 0)$ are given by:

$$ \dot{\tilde{z}}_1 = D(t) \left[ \frac{Y_{1m}m_{in}}{Y_{1m}m_{in} + Y_{1n}n_{in}} + \frac{Y_{1n}n_{in}}{Y_{1m}m_{in} + Y_{1n}n_{in}} - 1 \right] \tilde{z}_1, $$

$$ \dot{\tilde{z}}_2 = D(t) \left[ \frac{\mu_{2m}^{\max} K_{1m}}{\mu_{1m}^{\max} K_{2m}} \frac{Y_{1m}m_{in}}{Y_{1m}m_{in} + Y_{1n}n_{in}} + \frac{\mu_{2n}^{\max} K_{1n}}{\mu_{1n}^{\max} K_{2n}} \frac{Y_{1n}n_{in}}{Y_{1m}m_{in} + Y_{1n}n_{in}} - 1 \right] \tilde{z}_2. $$

Therefore, it follows that if the condition,

$$ \frac{\mu_{2m}^{\max} K_{1m}}{\mu_{1m}^{\max} K_{2m}} \frac{Y_{1m}m_{in}}{Y_{1m}m_{in} + Y_{1n}n_{in}} + \frac{\mu_{2n}^{\max} K_{1n}}{\mu_{1n}^{\max} K_{2n}} \frac{Y_{1n}n_{in}}{Y_{1m}m_{in} + Y_{1n}n_{in}} > 1, \quad (7) $$

is satisfied, then $(\bar{z}_1, \bar{z}_2) = (Y_{1m}m_{in} + Y_{1n}n_{in}, 0)$ is unstable. In particular, we expect that if the condition (7), derived for the local approximation $\tilde{z}$, is satisfied, the equilibrium $\bar{z}_2 = 0$ of the z-dynamics is thus unstable, which means that $z_2$ (and by extension, $x_2$) will not approach the domain making the target set $\mathcal{T}_\epsilon$ attainable, for a sufficiently small $\epsilon$. It should be noted that this condition (7) is not incompatible with the assumptions made on the system parameters, namely condition (3.1). Indeed, locally (under the previous approximation of the growth function), assumption (3.1) is satisfied if

$$ \left( \frac{K_{2m} \mu_{1m}^{\max}}{K_{1m} \mu_{2m}^{\max}} - 1 \right) \left( \frac{K_{2n} \mu_{1n}^{\max}}{K_{1n} \mu_{2n}^{\max}} - 1 \right) < 0. \quad (8) $$

In Appendix A, we provide a numerical example that satisfies both conditions: local instability of the point $(\bar{z}_1, 0)$ and a sign change in $\Delta$. This example particularly highlights the validity of the approximation of the $x_i$-dynamics by $z_1$ and $z_2$ under the slow-fast hypothesis. Indeed, in this numerical example, the trajectories of $x_i$ are well approximated by those of $z_i$, as illustrated in Figs. 18-19 in Appendix A. The phase portrait associated with the trajectories of the reduced model in this numerical example, with the target $\mathcal{T}_\epsilon$ defined for $\epsilon = 0.1$, starting from various initial conditions, is illustrated in Figure 3, where we observe that none of the trajectories reach or converge towards the target $\mathcal{T}_\epsilon$ and do not pass through it. This example shows that the target cannot be reachable for all parameter sets.

# 4 Necessary Optimality Conditions

## 4.1 Application of the PMP

In what follows, we assume the existence of an optimal control for OCP. Note that since (OCP) is a minimum-time control problem governed by a mono-input affine system w.r.t. the control, the existence follows straightforwardly from the reachability of the target set (see Assumption 3.3 and Section 3.3). To study OCP, we use the PMP which will allow us to obtain necessary conditions on optimal controls. In particular, it allows us to compute singular arcs which are intermediate strategies in comparison with bang-bang controls (i.e., controls taking extremal values of the admissible control set). The application of the PMP to (OCP) is as follows. Recall that the state variable is denoted by $X := (m, n, x_1, x_2)$ and
<figure><img src="image_4.png" /><figcaption>Figure 3: Phase portrait of the reduced model trajectories with the target set, starting from various initial conditions (gray points). The figure shows that none of the trajectories converge towards or reach $\mathcal{T}_\epsilon$ (all trajectories converge towards the red point). The model parameters in this example are given in Appendix A.</figcaption></figure>

let us write the adjoint variable (also called covector) as $\Lambda := (\lambda_m, \lambda_n, \lambda_{x_1}, \lambda_{x_2}) \in \mathbb{R}^4$. The Hamiltonian $H = H(X, \Lambda, \lambda^0, D)$ associated with (OCP) is the map $H: \mathbb{R}^4 \times \mathbb{R}^4 \times \mathbb{R} \times \mathbb{R} \to \mathbb{R}$ defined as

$$ H = \lambda^0 + \varphi D - \lambda_m \left( \frac{1}{Y_{1m}} \mu_{1m}(m) x_1 + \frac{1}{Y_{2m}} \mu_{2m}(m) x_2 \right) + \lambda_{x_1} \mu_1(m, n) x_1 + \lambda_{x_2} \mu_2(m, n) x_2 - \lambda_n \left( \frac{1}{Y_{1n}} \mu_{1n}(n) x_1 + \frac{1}{Y_{2n}} \mu_{2n}(n) x_2 \right), \quad (9) $$

where $\varphi$ is defined as:

$$ \varphi := \lambda_m(m_{in} - m) + \lambda_n(n_{in} - n) - \lambda_{x_1}x_1 - \lambda_{x_2}x_2. $$

According to the PMP, if $(X(\cdot), D(\cdot))$ is an optimal pair such that $X(\cdot)$ reaches the target set $\mathcal{T}_\epsilon$ at some terminal time $t_f \ge 0$, then, there exist an absolutely continuous function $\lambda: [0, t_f] \to \mathbb{R}^4$ (the covector) and $\lambda^0 \le 0$ such that the following conditions are satisfied:

* Non-triviality of the multipliers:

$$ (\lambda(\cdot), \lambda^0) \neq 0; \quad (10) $$

* Adjoint equation:

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

where $\nabla_X H$ denotes the gradient of the Hamiltonian with respect to $X$;

* Hamiltonian maximization condition:

$$ D(t) \in \operatorname*{argmax}_{\xi \in [0, D_{\max}]} H(X(t), \lambda(t), \lambda^0, \xi) \quad \text{a.e. } t \in [0, t_f]; \quad (12) $$

* Transversality condition (at the terminal time $t_f$):

$$ \lambda(t_f) \in -N_{\mathcal{T}_\epsilon}(X(t_f)), \quad (13) $$

where $N_{\mathcal{T}_\epsilon}(X) := \{p \in \mathbb{R}^4; p \cdot (Y - X) \le 0, \forall Y \in \mathcal{T}_\epsilon\}$ stands for the normal cone (in the sense of convex analysis) to $\mathcal{T}_\epsilon$ at some point $X \in \mathcal{T}_\epsilon$.
An extremal is defined as a quadruple $(X(\cdot), \lambda(\cdot), \lambda^0, D(\cdot))$ satisfying (1)-(10)-(11)-(12). When $\lambda^0 = 0$, the extremal is said abnormal whereas if $\lambda^0 \neq 0$, then it is said normal. Since (OCP) is autonomous, the Hamiltonian remains constant almost everywhere over $[0, t_f]$. Moreover, the terminal time being free, the PMP implies that the Hamiltonian is zero over $[0, t_f]$. Using (1), the adjoint equation can be equivalently written:

$$ 
\left\{ 
\begin{aligned}
\dot{\lambda}_m &= \left( \frac{1}{Y_{1m}} \mu'_{1m}(m)x_1 + \frac{1}{Y_{2m}} \mu'_{2m}(m)x_2 \right) \lambda_m - \mu'_{1m}(m)x_1 \lambda_{x_1} \\
&\quad - \mu'_{2m}(m)x_2 \lambda_{x_2} + D\lambda_m, \\
\dot{\lambda}_n &= \left( \frac{1}{Y_{1n}} \mu'_{1n}(n)x_1 + \frac{1}{Y_{2n}} \mu'_{2n}(n)x_2 \right) \lambda_n - \mu'_{1n}(n)x_1 \lambda_{x_1} \\
&\quad - \mu'_{2n}(n)x_2 \lambda_{x_2} + D\lambda_n, \\
\dot{\lambda}_{x_1} &= \frac{1}{Y_{1m}} \mu_{1m}(m)\lambda_m + \frac{1}{Y_{1n}} \mu_{1n}(n)\lambda_n - \mu_1(m, n)\lambda_{x_1} + D\lambda_{x_1}, \\
\dot{\lambda}_{x_2} &= \frac{1}{Y_{2m}} \mu_{2m}(m)\lambda_m + \frac{1}{Y_{2n}} \mu_{2n}(n)\lambda_n - \mu_2(m, n)\lambda_{x_2} + D\lambda_{x_2}.
\end{aligned}
\right. 
\quad (14) 
$$

Thanks to the transversality conditions and the expression of $\mathcal{T}_\varepsilon$, we obtain that at $t = t_f$:

$$ \lambda_m(t_f) = \lambda_n(t_f) = 0 \text{ and } \lambda_{x_1}(t_f) = \rho\varepsilon; \lambda_{x_2}(t_f) = -\rho, $$

for some $\rho \ge 0$. Actually, note that $\rho > 0$ (otherwise, from the linearity of the adjoint equation w.r.t. $\Lambda$, we would have $\Lambda \equiv 0$ implying that $\lambda^0 = 0$ from $H = 0$ and a contradiction). Thus, we have:

$$ \lambda_{x_1}(t_f) > 0; \lambda_{x_2}(t_f) < 0; \lambda_{x_1}(t_f)x_1(t_f) + \lambda_{x_2}(t_f)x_2(t_f) = 0. $$

The next lemma allows to exclude abnormal extremals depending on the terminal value of the state. If $\Delta(m(t_f), n(t_f)) \neq 0$, then the extremal is normal. Suppose that $\lambda^0 = 0$. Using that $H = 0$, $\lambda_{x_1}(t_f)x_1(t_f) + \lambda_{x_2}(t_f)x_2(t_f) = 0$ and $\lambda_m(t_f) = \lambda_n(t_f) = 0$, we obtain

$$ 0 = H = \lambda_{x_1}(t_f)x_1(t_f)\mu_1(m(t_f), n(t_f)) + \lambda_{x_2}(t_f)x_2(t_f)\mu_2(m(t_f), n(t_f)), $$

implying that $\Delta(m(t_f), n(t_f)) = 0$ and a contradiction. This ends the proof. Thanks to this lemma, if an abnormal extremal occurs, then we necessarily have $\Delta(m(t_f), n(t_f)) = 0$. Since the measure of $\{\Delta = 0\}$ is zero in $R$, we believe that generically, only normal extremals occur. In the rest of the paper, we suppose that only normal extremals occur for (OCP). From Lemma 4.1, we can thus conclude that at $t = t_f$, the adjoint vector satisfies the following (since $H = 0$):

$$ \lambda_{x_1}(t_f) = \frac{1}{x_1(t_f)[\mu_1(m(t_f), n(t_f)) - \mu_2(m(t_f), n(t_f))]}, \quad (15) $$

and $\lambda_{x_2}(t_f) = -\frac{\lambda_{x_1}(t_f)}{\varepsilon}$.

## 4.2 Switching Function

Following the expression of the Hamiltonian given by (9), we call *switching function* the map $t \mapsto \varphi(t)$ defined over $[0, t_f]$ as :

$$ \varphi(t) := \lambda_m(t)(m_{in} - m(t)) + \lambda_n(t)(n_{in} - n(t)) - \lambda_{x_1}(t)x_1(t) - \lambda_{x_2}(t)x_2(t). $$

The behavior of this function is crucial in order to derive the expression of an optimal control. From (12), the control law can be expressed as follows,

$$ D(t) \in \left\{ 
\begin{array}{ll}
\{D_{\max}\} & \text{if } \varphi(t) > 0, \\
\{0\} & \text{if } \varphi(t) < 0, \quad \text{a.e. } t \in [0, t_f]. \\
(0, D_{\max}) & \text{if } \varphi(t) = 0.
\end{array}
\right. $$

We call a bang arc $B_\pm$ any time interval where $\varphi > 0$, resp. $\varphi < 0$ (hence, $D = D_{\max}$, resp. $D = 0$). A singular arc is a time interval $[t_1, t_2]$ such that $\varphi(t) = 0$ for every $t \in [t_1, t_2]$. Note that the singular control
(i.e., the control that is associated with the singular arc) cannot be derived using (12), but it can be found by differentiating $\varphi$ an even number of times (see [30]). A switching time is any instant $t_s \in (0, t_f)$ at which the control switches from a bang arc $B_+$ to a bang arc $B_-$ or from a singular arc to a bang arc (or vice versa). Synthesizing an optimal control is in general difficult, and it relies on the properties of $\varphi$ and on the number of switching times. Whenever chattering arises (see [33]), then, this number is infinite over $[0, t_f]$ and this means that the singular arc is (at least) of second order [30]. We shall see that this situation does not occur in the present application model. However, knowing in advance the number of switchings is a difficult task from a theoretical viewpoint, that is why, the application of the PMP will also be combined with a direct method (see Section 5) which provides an insight into the structure of an optimal control together with the number of switching times. The conditions obtained via the PMP will then serve to check the optimality of the numerical solution. We now turn to properties of $\varphi$ that will be used in order to characterize singular arcs.

The function $\varphi$ is continuously differentiable over $[0, t_f]$ and satisfies $\varphi(t_f) = 0$. Moreover, one has

$$ \dot{\varphi} = (m_{\text{in}} - m) \left( \mu'_{1m}(m)x_1 \left[ \frac{1}{Y_{1m}}\lambda_m - \lambda_{x_1} \right] + \mu'_{2m}(m)x_2 \left[ \frac{1}{Y_{2m}}\lambda_m - \lambda_{x_2} \right] \right) + (n_{\text{in}} - n) \left( \mu'_{1n}(n)x_1 \left[ \frac{1}{Y_{1n}}\lambda_n - \lambda_{x_1} \right] + \mu'_{2n}(n)x_2 \left[ \frac{1}{Y_{2n}}\lambda_n - \lambda_{x_2} \right] \right), \quad (16) $$

for every $t \in [0, t_f]$. By differentiating $\varphi$ w.r.t. $t$, we find that

$$ \dot{\varphi} = (m_{\text{in}} - m) \dot{\lambda}_m - \dot{m}\lambda_m + (n_{\text{in}} - n) \dot{\lambda}_n - \dot{n}\lambda_n - [\dot{x}_1\lambda_{x_1} + x_1\dot{\lambda}_{x_1} + \dot{x}_2\lambda_{x_2} + x_2\dot{\lambda}_{x_2}]. $$

Since $\varphi$ does not involve the control $D(t)$, $\varphi$ is of class $C^1$. Moreover, using (1) and (14), we obtain (16). Now, combining the expression of the target set with the transversality conditions gives us $\varphi(t_f) = 0$.

Following the preceding computations, we can check that

$$ \dot{\varphi}(t_f) = \lambda_{x_2}(t_f)x_2(t_f)\alpha(m(t_f), n(t_f)), $$

with $\alpha(m, n) := (m_{\text{in}} - m)(\mu'_{1,m}(m) - \mu'_{2,m}(m)) + (n_{\text{in}} - n)(\mu'_{1,n}(n) - \mu'_{2,n}(n))$. If $\alpha(m(t_f), n(t_f)) = 0$, then, $\dot{\varphi}(t_f) = 0$, and we can expect that a normal extremal ends up with a singular arc in the spirit of [4] and [18] since we would then $\varphi(t_f) = \dot{\varphi}(t_f) = 0$. However, at this step, we can only prove that $\varphi(t_f) = 0$ which does not necessarily imply that the extremal ends up with a singular arc (although this is what we find numerically).

## 4.3 Computation of the Singular Control and Legendre-Clebsch Condition

The analysis of singular arcs and the computation of the singular control are based on the computation of $\ddot{\varphi}$. Hereafter, for conciseness, we shall sometimes omit the explicit dependence on time $t$ of $m(\cdot)$, $n(\cdot)$, $x_j(\cdot)$, $\lambda_m(\cdot)$, $\lambda_n(\cdot)$, $\lambda_{x_j}(\cdot)$ for $j = 1, 2$. Additionally, the dependence of certain functions on specific variables will not be explicitly indicated to simplify the notation. Let $\theta_m, \theta_n, \bar{\theta}_m, \bar{\theta}_n : [0, T] \to \mathbb{R}$ be defined as:

$$ \theta_m(t) := \theta_{1m}(t) + \theta_{2m}(t); \quad \theta_n(t) := \theta_{1n}(t) + \theta_{2n}(t), $$

$$ \bar{\theta}_m(t) := \bar{\theta}_{1m}(t) + \bar{\theta}_{2m}(t); \quad \bar{\theta}_n(t) := \bar{\theta}_{1n}(t) + \bar{\theta}_{2n}(t), $$

where for $j = 1, 2$, the functions $\theta_{jm}$ and $\theta_{jn}$ are defined as:

$$ \theta_{jm} := \mu'_{jm}(m)x_j \left[ \frac{1}{Y_{jm}}\lambda_m - \lambda_{x_j} \right]; \quad \theta_{jn} := \mu'_{jn}(n)x_j \left[ \frac{1}{Y_{jn}}\lambda_n - \lambda_{x_j} \right], $$

and where for $j = 1, 2$, the functions $\bar{\theta}_{jm}$ and $\bar{\theta}_{jn}$ are defined as:

$$ \bar{\theta}_{jm} := \mu''_{jm}(m)x_j \left[ \frac{1}{Y_{jm}}\lambda_m - \lambda_{x_j} \right]; \quad \bar{\theta}_{jn} := \mu''_{jn}(n)x_j \left[ \frac{1}{Y_{jn}}\lambda_n - \lambda_{x_j} \right]. $$
Using (14) and (16), we get $\dot{\lambda}_m = \theta_m + D\lambda_m$, $\dot{\lambda}_n = \theta_n + D\lambda_n$ and

$$ \dot{\varphi} = (m_{\text{in}} - m) \theta_m + (n_{\text{in}} - n) \theta_n. \quad (17) $$

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

$$ \begin{aligned} \dot{\theta}_m ={}& \bar{m}\bar{\theta}_m + \mu'_{1m}x_1 \left[ \mu_{1n} \left( \frac{1}{Y_{1m}}\lambda_m - \frac{1}{Y_{1n}}\lambda_n \right) + \frac{1}{Y_{1m}}\theta_m \right] \\ & + \mu'_{2m}x_2 \left[ \mu_{2n} \left( \frac{1}{Y_{2m}}\lambda_m - \frac{1}{Y_{2n}}\lambda_n \right) + \frac{1}{Y_{2m}}\theta_m \right]. \end{aligned} \quad (18) $$

(ii) The second order derivative of $\varphi$ can be expressed as:

$$ \begin{aligned} \ddot{\varphi} ={}& [-\theta_m + (m_{\text{in}} - m)\bar{\theta}_m]\dot{m} + [-\theta_n + (n_{\text{in}} - n)\bar{\theta}_n]\dot{n} \\ & + (m_{\text{in}} - m)\mu'_{1m}x_1 \left[ \mu_{1n} \left( \frac{1}{Y_{1m}}\lambda_m - \frac{1}{Y_{1n}}\lambda_n \right) + \frac{1}{Y_{1m}}\theta_m \right] \\ & + (m_{\text{in}} - m)\mu'_{2m}x_2 \left[ \mu_{2n} \left( \frac{1}{Y_{2m}}\lambda_m - \frac{1}{Y_{2n}}\lambda_n \right) + \frac{1}{Y_{2m}}\theta_m \right] \\ & + (n_{\text{in}} - n)\mu'_{1n}x_1 \left[ \mu_{1m} \left( \frac{1}{Y_{1n}}\lambda_n - \frac{1}{Y_{1m}}\lambda_m \right) + \frac{1}{Y_{1n}}\theta_n \right] \\ & + (n_{\text{in}} - n)\mu'_{2n}x_2 \left[ \mu_{2m} \left( \frac{1}{Y_{2n}}\lambda_n - \frac{1}{Y_{2m}}\lambda_m \right) + \frac{1}{Y_{2n}}\theta_n \right]. \end{aligned} \quad (19) $$

By differentiating $\theta_m$ w.r.t. $t$, we find that

$$ \begin{aligned} \dot{\theta}_m ={}& \bar{\theta}_m \dot{m} + \mu'_{1m} \dot{x}_1 \left[ \frac{1}{Y_{1m}} \lambda_m - \lambda_{x_1} \right] + \mu'_{1m} x_1 \left[ \frac{1}{Y_{1m}} \dot{\lambda}_m - \dot{\lambda}_{x_1} \right] \\ & + \mu'_{2m} \dot{x}_2 \left[ \frac{1}{Y_{2m}} \lambda_m - \lambda_{x_2} \right] + \mu'_{2m} x_2 \left[ \frac{1}{Y_{2m}} \dot{\lambda}_m - \dot{\lambda}_{x_2} \right]. \end{aligned} $$

Using (14) and (17), we obtain (18). Finally, (19) follows from the equality

$$ \ddot{\varphi} = (m_{\text{in}} - m)\dot{\theta}_m - \dot{m}\theta_m + (n_{\text{in}} - n)\dot{\theta}_n - \dot{n}\theta_n. $$

Following geometric optimal control theory (see, e.g., [30]), the expression of $\ddot{\varphi}$ is affine w.r.t. the control, so that the singular control can be obtained solving $\ddot{\varphi} = 0$ as we now do. Given a normal extremal, suppose that a singular arc occurs over a sub-interval $[t_1, t_2] \subset [0, t_f]$. Additionally, suppose that $t \mapsto (m_{\text{in}} - m(t))^2\bar{\theta}_m(t) + (n_{\text{in}} - n(t))^2\bar{\theta}_n(t)$ is non-zero over $[t_1, t_2]$. Then, the singular control $D_s$ can be expressed as,

$$ D_s(X, \lambda) := - \frac{\zeta_m + \zeta_n}{(m_{\text{in}} - m)^2 \bar{\theta}_m + (n_{\text{in}} - n)^2 \bar{\theta}_n}, \quad (20) $$

where

$$ \begin{aligned} \zeta_m(X, \lambda) :={}& \mu'_{1m}x_1 \left[ \mu_{1n} \left( \frac{1}{Y_{1m}}\lambda_m - \frac{1}{Y_{1n}}\lambda_n \right) + \frac{1}{Y_{1m}}\theta_m \right] (m_{\text{in}} - m) \\ & + \mu'_{2m}x_2 \left[ \mu_{2n} \left( \frac{1}{Y_{2m}}\lambda_m - \frac{1}{Y_{2n}}\lambda_n \right) + \frac{1}{Y_{2m}}\theta_m \right] (m_{\text{in}} - m) \\ & + \left( \frac{1}{Y_{1m}}\mu_{1m}x_1 + \frac{1}{Y_{2m}}\mu_{2m}x_2 \right) (\theta_m - (m_{\text{in}} - m)\bar{\theta}_m), \end{aligned} \quad (21) $$

and where $\zeta_n(X, \lambda)$ is defined similarly. This expression follows from (19) in which $\dot{m}$ and $\dot{n}$ are replaced according to (1) and $\dot{\varphi} = (m_{\text{in}} - m)\theta_m + (n_{\text{in}} - n)\theta_n \equiv 0$ (since $\dot{\varphi} = \dot{\varphi} = 0$). We now recall the *Legendre-Clebsch necessary optimality condition* which is a fundamental issue in geometric optimal
control in order to exclude a singular arc or to relate optimality of a singular arc together with the sign of the denominator defining the singular control $D_s$.

Following [30], recall that if a normal extremal is optimal and if a singular arc occurs over some time interval $[t_1, t_2] \subset [0, T]$, then, one must have

$$ \ddot{\varphi}_{|D} \geq 0, \qquad (22) $$

over $[t_1, t_2]$ where $\ddot{\varphi}_{|D}$ represents the derivative of $\ddot{\varphi}$ w.r.t. $D$. In our setting, following the expression of $\ddot{\varphi}$ given in the preceding lemma, this condition becomes:

$$ \ddot{\varphi}_{|D} = (m_{\text{in}} - m)^2 \bar{\theta}_m + (n_{\text{in}} - n)^2 \bar{\theta}_n \geq 0, \qquad (23) $$

over $[t_1, t_2]$. Whenever the inequality in (23) is strict, we say that the strict Legendre-Clebsch condition is verified. This condition will be used as follows in Section 5. We shall compute numerically the extremal trajectory (in particular the adjoint vector) so that we will be able to check if (23) is fulfilled over a singular arc. This will then serve as a certificate for optimality. Note that we cannot guarantee in advance that $D_s$ is always admissible, i.e., that $D_s(t) \in [0, D_{\text{max}}]$ along any singular arc and for almost all $t \in [t_1, t_2]$. Nevertheless, it is always possible to assume that $D_{\text{max}}$ is large enough to satisfy the upper bound. Additionally, we will see that $D_s$ is a particular arc that coincides with a turnpike feature and remains within $[0, D_{\text{max}}]$ whenever it is activated.

# 5 Direct Optimization and Link with the PMP

## 5.1 Description of the Method and Numerical Results

We now solve (OCP) numerically for two different sets of initial conditions using a direct method (see, e.g., [14]). This approach employs the Bocop toolbox ([8]), which transforms the optimal control problem into a finite-dimensional nonlinear programming problem (NLP) through the discretization of the dynamical control system. The resulting NLP is then solved using the IpOPT routine. For an overview of nonlinear programming techniques, we refer to [7], and for practical aspects of direct numerical methods, we refer to [6].

The direct optimization approach is a natural choice for solving the species separation OCP in this study, as it allows us to compute optimal trajectories without requiring prior knowledge of the control function, unlike indirect methods. However, to assess the validity of the obtained solutions, we subsequently compare the results obtained via direct methods with the first- and second-order PMP conditions derived in the previous section. Additionally, while alternative approaches such as those based on the Hamilton-Jacobi-Bellman (HJB) equation could provide globally optimal solutions—unlike our local approach—the high dimensionality of our optimality system, combined with the lack of controllability previously discussed through the counterexample, makes solving the HJB equation computationally challenging and even intractable. Given these challenges, the direct method remains a practical and effective choice for obtaining high-quality numerical solutions, despite their local nature. Its reliability is further supported by our numerical experiments, which confirm that the structure of the computed solutions remains consistent across different initial conditions. This numerical consistency, along with the theoretical validation through the first- and second-order PMP conditions from the previous section, strengthens our confidence in the robustness and accuracy of this approach. Furthermore, several heuristic and alternative optimization algorithms could be explored for solving the species separation OCP, but are beyond the scope of the paper. However, a detailed discussion on these numerical approaches, along with relevant references, can be found in [15].

In our direct optimization approach, the model parameters are specified in Tables 1 and 2 (in particular, $D_{\text{max}}$ is chosen to satisfy Assumption 3.1). The configuration settings for Bocop, including the discretization scheme and solver parameters, are provided in Table 3.

**Example 1:** for the initial condition $(m_0, n_0, x_1^0, x_2^0) = (1, 1, 1, 1)$ the solution found by bocop is of type $B_S$. The optimal control together with the associated trajectory are shown in Fig. 4 and 5.
<table><caption>Table 1: Model parameters related to the dynamics of the two species.</caption><thead><tr><th></th><th>Y<sub>jm</sub></th><th>Y<sub>jn</sub></th><th>K<sub>jm</sub></th><th>K<sub>jn</sub></th><th>μ<sup>max</sup><sub>jm</sub></th><th>μ<sup>max</sup><sub>jn</sub></th></tr></thead><tbody><tr><td>Species 1</td><td>1.1</td><td>1.3</td><td>0.2</td><td>0.15</td><td>1.98</td><td>2</td></tr><tr><td>Species 2</td><td>1.2</td><td>1.4</td><td>0.35</td><td>0.25</td><td>2.8</td><td>1.9</td></tr></tbody></table>

<table><caption>Table 2: Input substrate concentrations, maximal control, contamination rate $\varepsilon$.</caption><thead><tr><th>m<sub>in</sub></th><th>n<sub>in</sub></th><th>D<sub>max</sub></th><th>Contamination Rate ε</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>5</td><td>0.1</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Example 1: plot of the optimal control (of type $B_S$).</figcaption></figure>

**Example 2:** for the initial condition $(m_0, n_0, x_1^0, x_2^0) = (1, 1, 4, 3)$ the solution found by bocop is of type $B_B S$. The optimal control together with the associated trajectory are depicted on Fig. 6 and 7.

The numerical simulations reveal the turnpike-like behavior in the optimal controls. In our setting, this means that over a first time period, bang arcs occur and that the rest of the trajectory is composed of a singular arc playing the role of the turnpike. Note that we also found no terminal bang arcs (in contrast with the description of the turnpike phenomenon in [32]), probably because (OCP) is a minimum-time problem that does not involve an additional Mayer cost function. Finally, we would like to insist that we cannot apply the result of [32] since the problem is linear w.r.t. the control (however, the turnpike-like phenomenon is still present for solutions to (OCP)).

## 5.2 Verification of the Legendre-Clebsch Necessary Condition

We now turn to the verification of the optimality of the singular control obtained numerically in Examples 1-2. Doing so, we will check if the singular control (computed numerically) coincides with the singular control value $D_s$ computed theoretically by the PMP and if Legendre-Clebsch necessary condition is verified along the singular arc. The key point consists in recovering the adjoint vector $\lambda$ over $[0, T]$ from the Lagrange multipliers associated with the equality constraints in the NLP problem arising from the discretization of the control system. Note that these multipliers can be easily obtained thanks to bocop. Following this approach, we proceed in three steps to match the necessary conditions obtained from the PMP together with the numerical solution:

* switching function and its derivative : on Fig. 8 (related to Example 1) and on Fig. 11 (related to Example 2), we compute the switching function and its derivative thanks to the Lagrange multipliers and we verify that $\varphi = \dot{\varphi} = 0$ over the singular arc ;
<table><caption>Table 3: bocop settings.</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>1000</td></tr><tr><td>NLP Tolerance</td><td>&lt;1.0e-14</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Example 1: optimal trajectories $(m, n)$ and $(x_1, x_2)$ (above) and their respective adjoint vectors $(\lambda_m, \lambda_n)$ and $(\lambda_{x_1}, \lambda_{x_2})$ (below).</figcaption></figure>

• singular control : on Fig. 9 (related to Example 1) and on Fig. 12 (related to Example 2), we compute the singular control $D_s$ thanks to the Lagrange multipliers and we check that it coincides with the optimal control computed numerically;

• Legendre-Clebsch necessary condition : on Fig. 10 (related to Example 1) and on Fig. 13 (related to Example 2), we check the Legendre-Clebsch condition thanks to the computation of $\dot{\varphi}|_D$ (the adjoint vector is replaced by the Lagrange multipliers).

Thanks to these numerical comparisons valid for the two examples, we can deduce the following points. First, based on Figs. 8-11, we verify that $\varphi = \dot{\varphi} = \ddot{\varphi} = 0$ along the singular arc. Second, Fig. 9 and 12 indicate that the numerical optimal control coincides with the singular control (on the singular phase only). Finally, the Legendre-Clebsch condition is verified numerically over the singular arcs in both examples (see Fig. 10 and Fig. 13) which indicates that the singular control is indeed optimal. From these numerical comparisons, we deduce that the numerical solution indeed satisfies the necessary optimality
<figure><img src="image_7.png" /><figcaption>Figure 6: Example 2 : plot of the optimal control (of type $B_B_S$).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Example 2 : optimal trajectories $(m, n)$ and $(x_1, x_2)$ (above) and adjoint vectors $(\lambda_m, \lambda_n)$ and $(\lambda_{x_1}, \lambda_{x_2})$ (below).</figcaption></figure>

conditions obtained from the PMP which thereby confirms its optimality.

# 6 Characterization of the Turnpike

## 6.1 Study of the Static Optimization Problem

Following [32], we can associate to the Lagrange optimal control problem (OCP$_T$) a static optimization problem. This problem amounts to maximize the integrand in (OCP$_T$) over the steady-states of (1) parametrized by the constant control $D$ and it reads as follows:

$$ \max_{(m,n,x_1,x_2,D) \in \mathbb{R}^5} \Delta(m,n), \qquad (24) $$
<figure><img src="image_9.png" /><figcaption>Figure 8: Example 1: plot of $\varphi$ and $\dot{\varphi}$ (Fig. left) and plot of $\varphi$, $\dot{\varphi}$, $\ddot{\varphi}$ (Fig. right) using the Lagrange multipliers.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Example 1: plot of the numerical optimal control (in blue) together with the singular control $D_s$ (in orange) computed thanks to the PMP and to the Lagrange multipliers.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Example 1: verification of the Legendre-Clebsch Condition $\ddot{\varphi}_D \ge 0$ (along the singular arc).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Example 2: plot of $\varphi$ and $\dot{\varphi}$ (Fig. left) and plot of $\varphi$, $\dot{\varphi}$, $\ddot{\varphi}$ (Fig. right) using the Lagrange multipliers.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Example 2: plot of the numerical optimal control (in blue) together with the singular control $D_s$ (in orange) computed thanks to the PMP and to the Lagrange multipliers.</figcaption></figure>

under the equality and inequality constraints :

$$ \left\{ \begin{array}{l} (m_{\text{in}} - m) D - \frac{1}{Y_{1m}} \mu_{1m}(m) x_1 - \frac{1}{Y_{2m}} \mu_{2m}(m) x_2 = 0, \\ (n_{\text{in}} - n) D - \frac{1}{Y_{1n}} \mu_{1n}(n) x_1 - \frac{1}{Y_{2n}} \mu_{2n}(n) x_2 = 0, \\ [\mu_{1m}(m) + \mu_{1n}(n) - D] x_1 = 0, \\ [\mu_{2m}(m) + \mu_{2n}(n) - D] x_2 = 0, \\ m \in [0, m_{\text{in}}], n \in [0, n_{\text{in}}], x_1 \ge 0, x_2 \ge 0, D \in [0, D_{\text{max}}]. \end{array} \right. \qquad (25) $$

Note also that $(\text{OCP}_T)$ always has a solution since from (25), the pair $(m, n)$ takes values within a compact subset of $\mathbb{R}$. Here, we restrict our attention to solving this problem accurately for the data of Example 1 and 2. We find that the optimal solution to (24) is the point:

$$ (\hat{m}, \hat{n}, \hat{x}_1, \hat{x}_2, \hat{D}) = (0.13, 0.13, 2.1, 0, 1.7). $$

In particular, the optimal steady-state is such that $\hat{x}_2 = 0$ and $\hat{x}_1 > 0$ (which is meaningful in view of the target set to reach). Moreover, we find that $\Delta(\hat{m}, \hat{n}) = 0.3 > 0$. We used the routine optimization
<figure><img src="image_14.png" /><figcaption>Figure 13: Example 2 : verification of the Legendre-Clebsch Condition $\phi_D \ge 0$ (along the singular arc).</figcaption></figure>

(scipy.optimize) together with the method trust-constr which handles constraints. The numerical optimization has been done via a trust-region method. This numerical computation can be confirmed by inspecting more into details (24). In view of the target to reach, we impose that $x_1 \neq 0$. We thus have two cases:

• if $\hat{x}_1 > 0$ and $\hat{x}_2 = 0$, we get (thanks to the two first equations in (25)) that $n = f(m)$ for some function $f: \mathbb{R} \to \mathbb{R}$. The value of $\hat{m}$ is recovered by maximizing $m \mapsto \Delta(m, f(m))$ over $[0, m_{in}]$. Numerically, we also find that the Jacobian matrix of the system at this point is a Hurwitz matrix.

• if $\hat{x}_1 \neq 0$ and $\hat{x}_2 \neq 0$, then we obtain that $\Delta$ must be zero at the steady-state which is non optimal (since the first case provides a positive value of $\Delta$);

A general analysis of (24) should be (in theory) be carried out for all possible values of the pair $(x_1, x_2)$, in particular when $x_1$ is zero. However, in that case, one should also take into account stability properties of the (1) at the steady-state (in such a way that this steady-state can be reached). A theoretical analysis of the above (non-convex) optimization problem with arbitrary values for the parameters is involved and out of scope of the paper.

## 6.2 Link Between the Singular Arc and the Static Optimization Problem

Direct numerical simulations in Section 5 have revealed a turnpike-like behavior of optimal trajectories and that the triple $(m, n, D)$ remains close to a constant value over the singular phase that we call $(m^\dagger, n^\dagger, D^\dagger)$. In view of Lemma 3.2, a solution to (OCP) that reaches the target set in a time $T$ is necessarily a solution to $(OCP_T)$ (in which $T$ is fixed to this value). Now, if $T$ is large (typically, if $\varepsilon$ becomes very small), we can expect a solution to this problem to remain close to the point $(\hat{m}, \hat{n}, \hat{x}_1, \hat{x}_2, \hat{D})$ according to the turnpike phenomenon (although we cannot apply any result of [32] since $(OCP_T)$ is linear w.r.t. the control and the control $D(\cdot)$ is bounded). Coming back to the direct optimization method, this is exactly what we find. By extrapolating the constant values observed in Examples 1 and 2 along the turnpike, we observe that,

$$(m^\dagger, n^\dagger, D^\dagger) \approx (\hat{m}, \hat{n}, \hat{D}).$$

This observation shows the following interesting feature related to optimal solutions to (OCP). Even though there is no natural way to associate a static optimization problem to a minimum-time control problem like (OCP), optimal solutions to this problem possess a turnpike-like’s phenomenon around a constant value. Now, the corresponding point can be characterized as the optimal solution to a static optimization problem (namely Problem (24)) associated to an auxiliary Lagrange optimal control problem (namely Problem $(OCP_T)$). In particular, the point $(m^\dagger, n^\dagger)$ differs from the point $(m^*, n^*)$ which corresponds to the maximum of $\Delta$ over $R$ (see Fig. 14 below). But by its definition, $(m^*, n^*)$ does not
take into account the control system, that is why, it is not surprising that optimal trajectories $(m(\cdot), n(\cdot))$ do not necessarily converge to this point. Besides, reaching this point by solutions to (1) would also require controllability of the system.

## 6.3 Computation of $(m^\dagger, n^\dagger)$ for Other Values of the Parameters

In this section, we pursue the analysis of the turnpike phenomenon arising in (OCP) (see, Fig. 14 below) by computing the optimal solution for various values of the yield coefficients, first, when $Y_{1m} = Y_{1n}, Y_{2m} = Y_{2n}$ (see Fig. 15 and Fig. 16) and then in the case of arbitrary values for the yield coefficients (see Fig. 17). In the case where $Y_{1m} = Y_{1n}, Y_{2m} = Y_{2n}$, we have depicted on Fig. 16 several trajectories corresponding

<figure><img src="image_15.png" /><figcaption>Figure 14: Definition of the point $(m^\dagger, n^\dagger)$. In the gray region, the turnpike like solution is close to point $(m^\dagger, n^\dagger)$. The theoretical static optimal point is $(m^*, n^*)$.</figcaption></figure>

to the initial conditions $(m_0, n_0, x_1^0, x_2^0) \in [0.1, 1] \times [0, 1] \times \{1\} \times \{1\}$. We observe that all trajectories (even those such that $\Delta < 0$ at time $t = 0$) are such that $(m(\cdot), n(\cdot))$ approaches $(m^\dagger, n^\dagger)$ over a significant sub-interval of $[0, T]$.

For arbitrary yield coefficients $Y_{1m}, Y_{2m}, Y_{1n}, Y_{2n}$, we have plotted the optimal solution of (OCP) starting from (0.85,0.85,1,1) in the $(m, n)$ plane on Fig. 17 together with the new point $(m^\dagger, n^\dagger)$ for each value of $(Y_{1m}, Y_{2m}, Y_{1n}, Y_{2n})$. Note that the yield coefficients have been chosen in the interval $[0.1, 2]$. The variation of the point $(m^\dagger, n^\dagger)$ is of utmost importance from a practical viewpoint since they are correlated with the considered species of interest. We see that the point $(m^\dagger, n^\dagger)$ has a high variation as yield coefficients vary (see Fig. 17). This highlights the importance of having found a characterization of this point together with the ability of finding an optimal control of (OCP) for each possible value of the parameters defining the system.

Interestingly, in Fig. 17, many optimal trajectories exhibit multiple crossings at the interface between the two subsets $\{\Delta > 0\}$ and $\{\Delta < 0\}$. Such a behavior is non-intuitive and was also not observed in separation problems involving only one substrate (such as in [17, 18]). Indeed, the analysis carried out in [17, 18] revealed that optimal trajectories involve only one such crossing.

# 7 Conclusions

We fully characterized the optimal solution of a minimum-time optimal control problem aimed at separating two species competing for two substitutable nutrients in a chemostat system. The addition of a
<figure><img src="image_16.png" /><figcaption>Figure 15: Plot of $(m(\cdot), n(\cdot), D(\cdot))$ for those solutions to (1) starting at the initial conditions $(m_0, n_0, x_1^0, x_2^0) \in [0.1, 1] \times [0, 1] \times \{1\} \times \{1\}$. Each trajectory approaches $(m^*, n^*)$ over a significant sub-interval of $[0, T]$ as predicted by (24). Numerical simulations have been realized thanks to the direct optimization method.</figcaption></figure>

second nutrient source and the introduction of yields $Y_{ij}$ have considerable practical interest, making the model more intriguing to analyze while rendering the solution significantly different from those involving a single substrate. The solution we determined is not a mere generalization of the single-substrate case. Among the major differences, we highlight the following.

1. The attainability of the target is not systematic. We proved that for certain model parameter values, the target is unreachable due to a loss of controllability when using dilution to control a system with two different substrates.

2. The analysis of the PMP is significantly more complex and relies heavily on inputs from direct methods to obtain a plausible global description of the optimal solution.

3. The turnpike-like phenomenon observed in this problem is distinct and does not generalize the findings from the single-substrate case, whether in Monod or Droop models (such as in [17, 18]). In these previous cases, the turnpike solution always coincided with the point (noted * in the paper) that maximizes the difference between growth rates in static conditions. This point is unattainable in the two-substrate case, and another turnpike point (noted † in the paper) successfully steers the model trajectories towards the target when reachable. By formulating an auxiliary OCP and solving its associated static OCP, we characterized the dagger point.

In summary, this study reveals that, unlike a Lagrange optimal control problem, the minimum-time control problem we studied does not permit associating a turnpike computed as the solution to a
<figure><img src="image_17.png" /><figcaption>Figure 16: **Fig. left:** plot of the trajectories computed in Fig. 15 in the $(m, n)$ plane (under the same settings as in Fig. 15). **Fig. right:** plot of the corresponding functions $t \mapsto \Delta(m(t), n(t))$ for each trajectory. The figures highlight the fact that trajectories approach $(m^\dagger, n^\dagger)$ which differs from $(m^*, n^*)$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Plot of the points $(m^\dagger, n^\dagger)$ for various values of the yield coefficients $Y_{1m}, Y_{2m}, Y_{1n}, Y_{2n} \in [0.1, 2]$ together with the optimal solution starting from $(m^0, n^0) = (0.85, 0.85)$.</figcaption></figure>

static optimization problem. Nevertheless, considering the target set to reach, we could introduce such a problem by defining an auxiliary Lagrange optimal control problem. By combining numerical and theoretical approaches, we observed that optimal solutions exhibit a turnpike-like phenomenon near the optimal solution to the static optimization problem. This result highlights the importance of the turnpike phenomenon, which emerges in many optimal control problems, even in cases where classical theory does not apply (due to the lack of ellipticity in the Hamiltonian function). In particular, the characterization
of the turnpike point is more involved in the context of minimum-time control problems, as highlighted by this study.

This work opens up several interesting theoretical and practical research directions that we aim to explore in future studies. One immediate extension would be to consider more general growth rate models incorporating mortality and maintenance terms. Such modifications would fundamentally alter key properties of our framework, including the structure of the function $\Delta$, which plays a central role in our approach. Indeed, this shift could impact key system properties, such as invariance principles in the species competition model. Additionally, introducing growth functions with inhibition effects, such as Haldane-type kinetics, could further complicate the analysis by leading to a different formulation of the static OCP, which was a crucial step in characterizing the optimal solution. Another extension would be to consider the selection problem with two control inputs, such as the inflows of both substrates, a common practice in bioprocess engineering. This enhances controllability, potentially improving the turnpike point, but also increases complexity with interwoven bang and singular phases for both controls (see, e.g., [2], and also [13, 27]). Beyond these aspects, microbial species selection problems offer rich numerical, theoretical, and biological challenges, providing valuable insights into mathematical modeling and bioprocess optimization.

**Acknowledgements.** This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence Université Côte d'Azur under reference number ANR-15-IDEX-01. It also benefited from the collaborative facilities between Inria France and Inria Chile. The authors would like to warmly thank Dominique Y. Lagos (Federico Santa María Technical University, Valparaíso, Chile) for her help in plotting the data and producing the illustrative figures during her Master-1 research visit, and both Dominique Y. Lagos and Claudia Alvarez-Latuz for helpful discussions on the optimization problem.

# A Model Parameters in the Counterexample - Instability Condition and Numerical Validation of the Slow-fast Hypothesis

In the proposed counterexample, we consider the model parameters below, which satisfy both the instability conditions of the equilibrium point $(\bar{z}_1, 0)$ and the assumptions on the model (Assumption 1).

<table><thead><tr><th></th><th>Y<sub>jm</sub></th><th>Y<sub>jn</sub></th><th>K<sub>jm</sub></th><th>K<sub>jn</sub></th><th>μ<sup>max</sup><sub>jm</sub></th><th>μ<sup>max</sup><sub>jn</sub></th></tr></thead><tbody><tr><td>Species/strain 1</td><td>0.02</td><td>0.02</td><td>0.4</td><td>0.4</td><td>2.8</td><td>3.25</td></tr><tr><td>Species/strain 2</td><td>0.02</td><td>0.02</td><td>0.5</td><td>0.5</td><td>4.5</td><td>3.0</td></tr></tbody></table>

We also consider the initial conditions given by,

<table><thead><tr><th>Initial Conditions</th><th>x<sub>1</sub><sup>0</sup> = z<sub>1</sub><sup>0</sup></th><th>x<sub>2</sub><sup>0</sup> = z<sub>2</sub><sup>0</sup></th><th>m<sup>0</sup></th><th>n<sup>0</sup></th></tr></thead><tbody><tr><td></td><td>0.5</td><td>0.5</td><td>0.025</td><td>0.025</td></tr></tbody></table>

To ensure that the model satisfies the required instability and feasibility conditions, the parameters are chosen as follows.

1. The condition that $\Delta$ changes its sign is satisfied, since

$$ \left( \frac{K_{2m}\mu_{1m}^{\max}}{K_{1m}\mu_{2m}^{\max}} - 1 \right) \left( \frac{K_{2n}\mu_{1n}^{\max}}{K_{1n}\mu_{2n}^{\max}} - 1 \right) = -0.0787 < 0 \qquad (26) $$

2. The condition for $(\bar{z}_1, 0)$ to be unstable is also satisfied, since

$$ \frac{\mu_{2m}^{\max} K_{1m} Y_{1m} m_{in}}{\mu_{1m}^{\max} K_{2m} Y_{1m} m_{in} + Y_{1n} n_{in}} + \frac{\mu_{2n}^{\max} K_{1n} Y_{1n} n_{in}}{\mu_{1n}^{\max} K_{2n} Y_{1m} m_{in} + Y_{1n} n_{in}} = 1.012 > 1 \qquad (27) $$

For these parameters, we verify that the z-dynamics provides a satisfactory approximation of the full model dynamics, as illustrated in Figures 18-19.
<figure><img src="image_19.png" /><figcaption>Figure 18: Comparison between the full and reduced dynamics. Fig. left: dynamics of the original complete system satisfied by $(m, n, x_1, x_2)$. Fig. right: dynamics of the reduced system satisfied by $(z_1, z_2)$. The dynamics of $z_1$ and $z_2$, which are approximations of $x_1$ and $x_2$, are obtained under the assumption that $m$ and $n$ are fast variables, rapidly reaching a quasi-steady state.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Comparison between the full and reduced dynamics. We observe that the dynamics satisfied by $(z_1, z_2)$, see (6), provide a good approximation of the slow dynamics of $(x_1, x_2)$, thereby validating the hypothesis of the Tikhonov slow-fast reduction in the two-species chemostat model for the chosen parameter set in the counterexample.</figcaption></figure>

## References

[1] Almuashi, H.H., Almuallem, N.A. and El Hajji, M. The Effect of Leachate Recycling on the Dynamics of Two Competing Bacteria with an Obligate One-Way Beneficial Relationship in a Chemostat. Mathematics, 12(23), p.3819 (2024)

[2] Asswad, R., Djema, W., Bernard, O., Gouzé, J.L. Cinquemani, E. Optimization of microalgae biosynthesis via controlled algal-bacterial symbiosis. Proceedings of the 63rd IEEE Conference on Decision and Control-CDC, (arXiv preprint arXiv:2410.08575) (2024)

[3] Bardi, M., Dolcetta, I. C. Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations, (Vol. 12). Boston: Birkhäuser (1997)

[4] Bayen, T., Mairet, F. Optimization of the separation of two species in a chemostat. Automatica J. IFAC, vol. 50, 4, pp. 1243-1248 (2014)
[5] Bayen, T., Mairet, F. Optimisation of strain selection in evolutionary continuous culture. Internat. J. Control, 90(12), pp. 2748-2759 (2017)

[6] Betts, J. T. Practical methods for optimal control and estimation using nonlinear programming. Society for Industrial and Applied Mathematics (2010)

[7] Biegler, L. T. Nonlinear programming: concepts, algorithms, and applications to chemical processes. Society for Industrial and Applied Mathematics (2010)

[8] Bonnans, F.J., Giorgi, D., Grelard V., Heymann, B., Maindrault S., Martinon, P., Tissot, O., Liu, J. Bocop – A collection of examples. Technical Report, Inria (2017)

[9] Bonnard, B. and Chyba, M. Singular trajectories and their role in control theory. Vol. 40. Springer Science & Business Media (2003)

[10] Bonnefond, H., Lie, Y., Lacour, T., Saint-Jean, B., Carrier, G., Pruvost, E., Sciandra, A. Dynamical Darwinian selection of a more productive strain of Tisochrysis lutea. Algal Research, 65, 102743 (2022)

[11] Boscain, U., Piccoli, B. Optimal syntheses for control systems on 2-D manifolds. Volume 43. Springer Science & Business Media (2003)

[12] Caillau, J.B., Djema, W., Gouzé, J.L., Maslovskaya, S. and Pomet, J.B. Turnpike property in optimal microbial metabolite production. Journal of Optimization Theory and Applications, 194(2), pp. 375-407 (2022)

[13] Cherkaoui-Dekkaki, O., Djema, W., Raissi, N., Gouzé, J.L. and El Khattabi, N. Optimal control of waste recovery process. International Journal of Dynamics and Control, 12(12), pp.4244-4262 (2024)

[14] Caillau, J.B., Ferretti, R., Trélat, E. and Zidani, H. An algorithmic guide for finite-dimensional optimal control problems. In Handbook of Numerical Analysis (Vol. 24, pp. 559-626). Elsevier (2023)

[15] Della Marca, R. and d'Onofrio, A. Volatile opinions and optimal control of vaccine awareness campaigns: Chaotic behaviour of the forward-backward sweep algorithm vs. heuristic direct optimization. Communications in Nonlinear Science and Numerical Simulation, 98, p.105768 (2021)

[16] Djema, W., Bernard, O. and Bayen, T. Optimal control separating two microalgae species competing in a chemostat. Proceedings of the 59th IEEE Conference on Decision and Control (CDC), pp. 2368-2373 (2020)

[17] Djema, W., Bernard, O., Giraldi, L. Separating two species of microalgae in photobioreactors in minimal time. Journal of Process Control, 87, pp. 120-129 (2020)

[18] Djema, W., Bayen, T. and Bernard, O. Optimal Darwinian Selection of Microorganisms with Internal Storage. Processes, 10(3), 461 (2022)

[19] Djema, W., Giraldi, L., Maslovskaya, S. and Bernard, O. Turnpike features in optimal selection of species represented by quota models. Automatica, 132, p.109804 (2021)

[20] Gutiérrez Mena, J., Kumar, S. and Khammash, M. Dynamic cybergenetic control of bacterial co-culture composition via optogenetic feedback. Nature Communications, 13(1), p.4808 (2022)

[21] Harmand, J., Lobry, C., Rapaport, A., Sari, T. The chemostat: Mathematical theory of microorganism cultures. Vol. 1. John Wiley & Sons (2017)

[22] Innerarity Imizcoz, J., Djema, W., Mairet, F., Gouzé, J.L. Optimal resource allocation in microorganisms under periodic nutrient fluctuations. Journal of Theoretical Biology, 595, 111953 (2024)
[23] Khalil, H. Nonlinear systems, third edition. Prentice Hall (2002)

[24] León, J.A., Tumpson, D.B. Competition between two species for two complementary or substitutable resources. Journal of Theoretical Biology, 50(1), pp. 185-201 (1975)

[25] Masci, P., Bernard, O., Grognard, F. Continuous selection of the fastest growing species in the chemostat. IFAC Proceedings Volumes, 41(2), pp. 9707-9712 (2008)

[26] Mata, T. M., Martins, A. A., Caetano, N. S. Microalgae for biodiesel production and other applications: a review. Renewable and Sustainable Energy Reviews, 14(1), pp. 217-232 (2010)

[27] Mazel, P., Djema, W., Grognard, F. Optimal control for a combination of cancer therapies in a model of cell competition. Proceedings of the 63rd IEEE Conference on Decision and Control-CDC (2024)

[28] Monod, J. Recherches sur la Croissance des Cultures Bactériennes. Hermann, Paris (1942)

[29] Pontryagin, L.S., Boltyanskiy V.G., Gamkrelidze, R.V., Mishchenko, E.F. Mathematical theory of optimal processes. The Macmillan Company (1964)

[30] Schättler, H., Ledzewicz, U. Geometric optimal control: theory, methods and examples, (Vol. 38). New York: Springer (2012)

[31] Smith, H. L., Waltman, P. The theory of the chemostat: dynamics of microbial competition, (Vol. 13). Cambridge University Press (1995)

[32] Trélat, E., Zuazua, E. The turnpike property in finite-dimensional nonlinear optimal control. J. Differential Equations, 258(1), pp. 81-114 (2015)

[33] Zelikin, M. I., Borisov, V. F. Theory of chattering control: with applications to astronautics, robotics, economics, and engineering. Springer Science & Business Media (2012)