# Optimal Control of the Controlled Lotka-Volterra Equations with Applications - The Permanent Case

Bernard Bonnard, Jérémy Rouot

► To cite this version:

Bernard Bonnard, Jérémy Rouot. Optimal Control of the Controlled Lotka-Volterra Equations with Applications - The Permanent Case. SIAM Journal on Applied Dynamical Systems, 2023, 22 (4), pp.2761-2791. 10.1137/22M151978X. hal-03757060v2

HAL Id: hal-03757060

https://inria.hal.science/hal-03757060v2

Submitted on 4 Jul 2023

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 CC BY 4.0 - Attribution - International License
# Optimal Control of the Controlled Lotka-Volterra Equations with Applications. The Permanent Case*

Bernard Bonnard† and Jérémy Rouot‡

**Abstract.** In this article motivated by the control of complex microbiota in view to reduce the infection by a pathogenic agent, we introduce the theoretical frame from optimal control to analyze the problem. Two complementary approaches can be applied in the analysis: one is the so-called permanent case, where no digital constraints are concerning the control (taken as a measurable mapping) versus the sampled-data control case taking into account the logistic constraints, e.g. frequency of the medical interventions. The model is the n-dimensional Lotka-Volterra equation controlled using either probiotics or antibiotic agents or transplantation and bactericides. In this article, we concentrate to the permanent case associated to probiotic or antibiotic agent. The Maximum Principle is used to parameterize the geodesics and the optimal synthesis boils down to analyze mainly the singular trajectories and their concatenation with bang arcs.

**Key words.** Optimal control in the permanent case, biomathematics and population dynamics, geometric control theory.

**MSC codes.** 49K15, 92B05, 93C10, 93C15, 92D25

**1. Introduction.** The book by Vito Volterra [28] "Leçons sur la théorie mathématique pour la lutte pour la vie" leads to the Lotka-Volterra model to predict in a general frame the evolution on interacting biological species. The problem was studied independently by Lotka which makes the connection with chemical networks. The relation with control systems was already present in the original predator-prey model with two species set by Umberto D'Ancona to explain the evolution of the species in relation with reduction of the fishing activity during the first World War. The original memoir [28] starts with an interesting discussion about the evolution of the species in relation with integrability properties of the conservative model. It validates D'Ancona observations and opens the road to analyze different problems of populations dynamics. This leads to extend the model to the non conservative 2d-case in any dimension.

Recently our attention was attracted by the works of Jones et al. [17] based on the model by Stein et al. [27] in an attempt to model and cure a gut mouse microbiota infected by the C. difficile bacteria which leads to a 11-dimensional controlled Lotka-Volterra model, using either fecal transplantation or one antibiotic agent. The problem is analyzed using 2d-reduced (projection) non conservative Lotka-Volterra model. This model validates the effect of an antibiotic agent prior to infection and followed by a single fecal injection to cure the infection by constructing a separatrix which allows to decide about success or failure of the

<footnote>*Submitted to the editors DATE.</footnote>

<footnote>**Funding:** This work benefited from the support of the FMJH Program PGMO and from the support of EDF; Thales, Orange and by the program PEPS "Jeunes chercheurs et jeunes chercheuses" of Insmi.</footnote>

<footnote>†Institut Mathématique de Bourgogne, Université de Dijon, 9 rue Alain Savary, 21078 Dijon and Inria Sophia Antipolis Equipe McTao (bernard.bonnard@u-bourgogne.fr).</footnote>

<footnote>‡Laboratoire de Mathématiques de Bretagne Atlantique, UMR 6205, 6 avenue Le Gorgeu, 29238 Brest (jeremy.rouot@univ-brest.fr).</footnote>
procedure due to delay in the transplantation or insufficient dosing.

To understand this approach let us introduce the Lotka–Volterra dynamics and analyze its limits, already discussed in [28]. The system is given by the equations

$$ (1.1) \qquad \frac{dx}{dt}(t) = (\operatorname{diag}(x(t)))(Ax(t) + r), $$

where $x = (x_1, x_2, \dots, x_n)^\top$ belongs to the positive quadrant $x_i \ge 0$ and is the vector of populations of the interacting species, $\operatorname{diag}x$ denotes in short the diagonal matrix with coefficients $x_i$, $A = (a_{ij})$ is the matrix containing the interacting coefficients between the species and $r = (r_1, r_2, \dots, r_n)^\top$ is the vector describing the individual growth rate without interaction. Assuming the matrix $A$ invertible, the system possesses a unique interior equilibrium given by $x_0 = -A^{-1}r$. The model is only valid provided that each species is contained in an interval $[\varepsilon, M]$, $\varepsilon > 0$ and hence in the model only the so-called persistent trajectories defined for positive times and contained in $[\varepsilon, M]^n$ have physical signification. Therefore the differential equation represents an analytic continuation of the dynamics to the whole euclidean space. But the model in fine describes the interaction between the interior equilibrium and equilibria of a hierarchy of reduced Lotka–Volterra dynamics associated to extinction of the different species, e.g. extinction of $x_1$ leads to analyze the reduced model of $n-1$ interacting species where we substitute in the dynamics $x \to x = (x_2, \dots, x_n)^\top$ and again we can compute the interior equilibria of the reduced dynamics. Therefore this leads to interpret the model as a system of interaction between a given number $k$ of isolated equilibria with $k \le 2^n$.

From this point of view the analysis of [17] is precisely to describe a policy in the experimental setting where a mouse is treated prior to infection by antibiotic and a final fecal injection in order to cure the C. difficile infection. The 2d-reduced model is with a saddle interior point with stability and unstability domain delimited by separatrixes and the success of the therapy is to reach a stability domain of an healthy equilibrium.

At this time the connection with the network of chemical species pointed by Lotka has to be made in relation with intense research activities in the seventies. They were motivated by the analysis of chemical batch reactors and realized mainly by Feinberg-Horn-Jackson, see for instance the recent book [15]. The study is related to the graph of the reactions using the concept of deficiency.

Similarly concerning the Lotka–Volterra equations there was an intense research activities due in particular to Zeeman, Smale, Hirsch to analyze the dynamics in the frame of the theory of dynamical systems [29, 16, 26], the limit of the studies being related to complex chaotic behaviors, see [3] for a seminal presentation of the computational complexity of the problem.

The analysis from geometric optimal control viewpoint of the problem of reducing the infection of a complex microbiote using the Lotka–Volterra model leads to introduce the polysystem $D = \{X, Y\}$, where $X$ is the vector field $(\operatorname{diag}x)(Ax + r)$ and $Y$ is an additional vector field associated to a specific treatment and to study the action of the pseudo-semigroup $S(D)$ generated by concatenation of the positive orbits of the vectors fields denoted respectively $\varphi_t = \exp tX, \psi_t = \exp tY$. The accessibility set at time $t$ denoted $A(x_0, t)$ is the orbit $S(D)(x_0)$ when the total time is $t$ and its boundary contains the extremities of the time minimal geodesics. It can be evaluated for small time using Lie brackets computations of $X$ and $Y$. This boundary can have a complicated structure but an intense research activities
at the end of the eighties was devoted to stratify this set under some suitable assumptions in relation with the problem of computing a regular synthesis, see in particular [10], [12]. In this research program the Maximum Principle [24] is a very powerful tool to analyze the optimal solutions since it gives a parametrization of the boundary of the accessibility set, where the geodesics dynamics is the projection of an n-Hamiltonian dynamics that is a differential system in dimension 2n and this system is not smooth.

Let us formulate this principle to analyze for instance the effect of an antibiotic or probiotic treatment only, whence the control system takes the form

$$ \frac{dx}{dt}(t) = X(x(t)) + u(t)Y(x(t)) $$

where the control $u(\cdot)$ in the permanent case belongs to the set $\mathcal{U}$ of bounded measurable mappings valued in the convex set $U = [-1, +1]$, where $u = -1$ corresponds to zero dosing that is $X - Y = (\mathrm{diag}x)(Ax + r)$ while to maximal dosing corresponds the dynamics $(\mathrm{diag}x)(Ax + r) + 2Y(x)$, where $Y(x)$ is related to the action of the antibiotic or probiotic agent.

Consider the problem of reaching in minimum time a terminal manifold $N$ related to reduce the infected agent $x_1$, e.g. reach a small population $x_1$ in minimal time (a dual formulation of the problem is min $x_1$ for fixed final time).

**Theorem 1.1.** The Maximum Principle tells us that if $(x(\cdot), u(\cdot))$ is a time minimal control-trajectory pair on $[0, t_f]$, then there exists a non vanishing adjoint n-dimensional vector $p(\cdot)$ such that the triple $(z(\cdot), u(\cdot))$, $z = (x, p)$ satisfies the equations that we introduce next.

Denoting $H(z, u) := H_X(z) + uH_Y(z)$ the pseudo-Hamiltonian, where if $Z$ is a vector field, $H_Z(z) = p \cdot Z(x)$ denotes the Hamiltonian lift and $M$ is the maximized or true Hamiltonian $\max_{u \in [-1, +1]} H(z, u)$. Then for almost every $t$ one has

$$ (1.2) \qquad \frac{dx}{dt}(t) = \frac{\partial H}{\partial x}(x(t), p(t), u(t)), \quad \frac{dp}{dt}(t) = -\frac{\partial H}{\partial p}(x(t), p(t), u(t)). $$

Moreover the optimal control satisfies the maximization condition a.e.

$$ (1.3) \qquad H(z(t), u(t)) = M(z(t)), $$

and $M$ is a nonnegative constant.

At the terminal time $t_f$ the transversality condition holds

$$ (1.4) \qquad p(t_f) \perp T^*_{x(t_f)} N. $$

From the maximization condition one deduces that an optimal control is the concatenation of:

* Regular subarcs where $u(t) = \mathrm{sign}H_Y(z(t))$ a.e..

* Singular subarcs defined by the implicit relation $H_Y(z(t)) = 0$.

Such singular arcs define the geodesics solutions when relaxing the control bound to the whole $u \in \mathbb{R}$. They form an Hamiltonian flow constrained to the switching set $\Sigma: H_Y(z) = 0$ in which they filled in general a subset of codimension one. Hence they play an important role in our study, see for instance the reference [4].
The main objective of this article is to analyze geometrically this dynamics in the case of the controlled Lotka–Volterra equations and from the control optimal point of view it is the analogue of the program of classifying the geometric dynamics for the Lotka–Volterra model. The limit being essentially the same that is to handle the curse of dimension and the computational complexity. Making the connection of the problem of minimizing the infection with the problem of maximizing the production of one species for chemical network, our objective is to use in this frame a series of articles started in the eighties, see for instance [6, 21] in the frame of geometric optimal control [1, 4], aiming to optimize the production of batch chemical reactors and recently pursued in [7].

Moreover in practice our objective is to compute an approximation of the time minimal synthesis that is to reach the terminal manifold for every initial condition determining the closed loop optimal control: $x \to u^*(x)$. Regularity conditions have to be satisfied in order to define the solutions and they are related to the regularity properties of the time minimal value function.

The article is organized as follows. In section 2 the controlled Lotka–Volterra model is introduced and the optimal control problem presented in the frame of permanent controls. The Maximum Principle leads to the classification of the geodesics in the context of geometric optimal control and singularity theory using the seminal earliest contributions of [19, 20, 14]. The section 3 is based on the series of articles [6, 21, 7], dealing with a terminal manifold of codimension one. They are rather technical and our contribution being to introduce two main concepts. The first one is the notion of Whitney chart to determine in an appropriated coordinates system the time minimal synthesis in a neighborhood of the terminal manifold using the construction of semi-normal forms to estimate the switching and cut loci up to all cases of codimension two for the $C^\infty$-Whitney topology. Secondly the concept of unfolding from singularity theory is introduced in our control frame to reduce the classification to 2d and 3d cases and using their description (dictionary) of the aforementioned references. The final section is devoted to the analysis and the classification of singular arcs in the 2d and 3d cases to deduce the time minimal syntheses. In this program the computation can be automatized in the 2d case using two classical invariants: the collinear set and the singular locus. In dimension 3 the problem is intricate due to the complexity of the singular dynamics. But a program of computations in the general case can be outlined based on gluing Whitney charts as an alternative of programs as for example in [2] to derive patchy feedbacks to approximate the time minimal function, our construction differing by the dominant use of the feedback singular control.

## 2. Model of Controlled Lotka–Volterra Equation and Optimal Control.

### 2.1. A quick tour in the 2d–Lotka–Volterra predator-prey model.
The original Lotka–Volterra model analyzed in [28] describes, in the frame of (conservative) integrable dynamics, the interaction between two species. More precisely it was constructed to explain the evolution of the averaged populations of two fishing species in relation with diminution of the fishing activity observed by D’Ancona during the first World War and succeeds to explain the observation. From dynamical point of view it concerns the case of a center, the mechanical analogue being the oscillating non linear pendulum and it was extended to the case of an integrable saddle by Volterra himself. The limit of the model are clearly indicated in the
memoir and a parallel can be made with the 2d-reduced model of the gut microbiome that we
are analyzing in this paper.

We conserve the same notation than in the memoir. In the oscillatory situation the prey-predator population $N = (N_1, N_2)$ satisfies the dynamics:

$$
(2.1) \qquad \frac{dN_1}{dt} = N_1(\varepsilon_1 - \gamma_1 N_2), \quad \frac{dN_2}{dt} = -N_2(\varepsilon_2 - \gamma_2 N_1)
$$

where $(\varepsilon_1, \varepsilon_2, \gamma_1, \gamma_2) > 0$, which can be written in a general form as

$$
(2.2) \qquad \frac{dN_1}{dt} = N_1(\lambda_1 + \mu_1 N_2) \quad \frac{dN_2}{dt} = N_2(\lambda_2 + \mu_2 N_1)
$$

for real parameters $\lambda_i, \mu_i, i = 1, 2$.

This leads to the equation

$$
(2.3) \qquad \mu_2 \frac{dN_1}{dt} + \lambda_2 \frac{1}{N_1} \frac{dN_1}{dt} - \mu_1 \frac{dN_2}{dt} - \lambda_1 \frac{1}{N_2} \frac{dN_2}{dt} = 0
$$

and integrating one gets

$$
\mu_2 N_1 + \lambda_2 \ln N_1 - (\mu_1 N_2 + \lambda_1 \ln N_2) = \text{constant.}
$$

Hence

$$
N_1^{\lambda_2} e^{\mu_2 N_1} = C N_2^{\lambda_1} e^{\mu_1 N_2},
$$

where $C$ is a constant depending upon the initial conditions $(N_1(0), N_2(0))$.

Volterra describes the solution using the auxiliary curves

$$
(2.4) \qquad Y = N_1^{-\varepsilon_2} e^{\gamma_2 N_1}, \quad X = N_1^{\varepsilon_1} e^{-\gamma_1 N_2}
$$

so that the solution can be locally either represented as a graph $Y = CX$ or $X = CY$.

One denotes

$$
\Omega := (K_1, K_2) = \left( \frac{\varepsilon_2}{\gamma_2}, \frac{\varepsilon_1}{\gamma_1} \right)
$$

the interior equilibrium.

Exact formulae give the time evolution of the population $t \mapsto N_i(t)$. But from physical point of view, what matters is the averaged population $\langle N_i \rangle$ and a simple computation gives the following.

**Theorem 2.1.** The averaged populations are given by:

$$
\langle N_i \rangle := \frac{1}{T} \int_0^T N_i(t) \mathrm{d}t = K_i, \quad i = 1, 2,
$$

which are not depending upon the initial conditions but only from the equilibrium $\Omega$.
<figure><img src="image_2.png" /><figcaption>Figure 1. Volterra memoir. Left: center-case. Right: saddle case.</figcaption></figure>

The Lotka–Volterra model can be normalized in the same category using the dimensionless variables:

$$n_1 := \frac{N_1}{K_1}, \quad n_2 := \frac{N_2}{K_2}$$

and we get

$$(2.5) \qquad \frac{dn_1}{dt} = \varepsilon_1 n_1 (1 - n_2), \quad \frac{dn_2}{dt} = -\varepsilon_2 n_2 (1 - n_1)$$

and the linearized system at $\Omega$ takes the form

$$(2.6) \qquad \frac{\delta n_1}{dt} = -\varepsilon_1 \delta n_2, \quad \frac{\delta n_2}{dt} = \varepsilon_2 \delta n_1.$$

The geometric construction of the dynamics in the memoir is based on the auxiliary curves (2.4) to compute $Y = CX$.

This opens the path to treat in the same frame the case where $\varepsilon_1\varepsilon_2 < 0$ so that $\Omega$ is a saddle point to get the form $XY = C$. Volterra gives a geometric representation of the two cases in a single figure, see Fig. 1.

Let us analyze in the oscillatory case the role of the fishing activity, introducing control in the model. The system takes the form:

$$(2.7) \qquad \frac{dN_1}{dt} = (\varepsilon_1 - \alpha\lambda - \gamma_1 N_2)N_1, \quad \frac{dN_2}{dt} = -(\varepsilon_2 + \beta\lambda - \gamma_2 N_1)N_2,$$

where $\alpha, \beta \ge 0$ are the modes of destruction and $\lambda \ge 0$ is the intensity.

Assuming $\varepsilon_1 - \alpha\lambda > 0$ so that the population is still oscillating, the averaged values of

$N_1, N_2$ become

$$\frac{\varepsilon_2 + \beta\lambda}{\gamma_2}, \frac{\varepsilon_1 - \alpha\lambda}{\gamma_1} \quad \text{versus} \quad \frac{\varepsilon_2}{\gamma_2}, \frac{\varepsilon_1}{\gamma_1}$$
<figure><img src="image_3.png" /><figcaption>Figure 2. Jones et al. model.</figcaption></figure>

without fishing activity (see Theorem 2.1). This confirms the observations by D’Ancona.

This simple model can be generalized to analyze, in the same frame, controlled populations dynamics introducing dissipation in the model associated to non zero coupling coefficients, to get a 2d-control model of the form

$$ (2.8) \qquad \frac{dx}{dt} = (\operatorname{diag} x) [(r + Ax) + u \epsilon], $$

where A is the interaction matrix, $\epsilon$ is the sensitivity vector and the control u is the intensity.

Hence, in this frame the role of the constant control $u = 1$ is to shift the interior equilibrium and the spectrum of the linearized dynamics.

The previous discussion clarifies the construction by Jones et al. of a 2d-Lotka-Volterra model to describe a complex microbiome with an interior saddle point and the interaction between the C. difficile population $x_1$ and the healthy microbiome population aggregated into a single population $x_2$. Noting again $\Omega = (K_1, K_2)$ the interior saddle equilibrium and normalizing: $x_1 \to \frac{x_1}{K_1}$, $x_2 \to \frac{x_2}{K_2}$, this leads to a system with four equilibria: $\Omega$, the origin and two non interior equilibria $O_1, O_2$ associated respectively to the reduced dynamics with $x_1$ or $x_2$ equal zero. This can be interpreted as gluing together two prey-predators models (see Fig. 2).

The model simulates the following medical protocol:

$$ \begin{array}{l} \text{day 0} \quad \rightarrow \quad \text{day } k \quad \rightarrow \quad \text{day } k' > k, \\ \text{(healthy mouse} \quad \text{(C. difficile infection)} \quad \text{(fecal transplantation)} \\ \text{& antibiotic)}, \end{array} $$

where the control actions are the choice of the day of the transplantation and the composition of the fecal injection.

The final state posterior to the fecal injection is either the healthy state $O_1$ or the infected state $O_2$.

A separatrix $\Sigma$ is constructed in the 2d-model indicating the success or the failure of the transplantation therapy. This yields decision algorithms based on the computations of the parameters of the reduced models using the observations. The infection has to be compensated by the fecal transplantation. The mathematical limit of the model being that each species
population shall satisfies $M \ge x_i \ge \varepsilon$, in particular if the $x_2$ population becomes too small the model is not valid. Also from medical point of view a probiotic injection versus antibiotic injection increases the healthy population aiming to struggle against infection.

**2.2. Controlled Lotka–Volterra model and optimal control.** Next we introduce the definitions and concepts in a general frame.

**Definition 2.2.** Let $x = (x_1, \dots, x_n)^\top \in \mathbb{R}_+^n$ be the state of interacting species, $x_1$ being the infected agent, $x' = (x_2, \dots, x_n)^\top$ being the state of healthy agents. The dynamics is described by

$$ (2.9) \qquad \frac{dx}{dt} = (\mathrm{diag}x)(Ax + r) $$

the matrix $A = (a_{ij})$ being the matrix of coefficients of interaction and $r = (r_1, \dots, r_n)^\top$ is the vector of individual growth rate. We denote by $M^+ = \mathbb{R}_+^n$ the invariant domain $x_i > 0$ and $M^\sim$ the union of $M^+$ with its boundary. The dynamics is called regular if $A$ is invertible. The interior equilibrium is the point $\Omega = (K_1, \dots, K_n)$ given by $x = -A^{-1}r$.

We note $n = (n_1, \dots, n_n)^\top$ the dimensionless coordinates so that $\Omega$ is identified to $(1, \dots, 1)^\top$.

One can associate to (2.9) an hierarchy of dynamics replacing $x \to x' = (x_2, \dots, x_n)^\top$ and this leads to up to $2^n$ equilibria for the dynamics in the physical space, which can be easily computed by recurrence.

The dynamics can be compactified using Poincaré compactification, identifying $\mathbb{R}^n$ to the hyperplane $(x, z = 1)$ in $\mathbb{R}^{n+1}$ to define the system:

$$ (2.10) \qquad \frac{dx}{dt} = (\mathrm{diag}x)(Ax + r), \quad \frac{dz}{dt} = 0, $$

where the right-member can be homogenized to define an homogeneous vector field which can be projected on the $n$-sphere $S^n$.

Each equilibrium of the hierarchy of dynamics can be classified according to the $L$- (linear) stability status associated to the linearized dynamics at equilibrium.

Our study is related to the interaction of $k$ non interior equilibria interacting with the interior equilibrium and one can introduce a model reduction consisting in a polynomic dynamics of the form: $\frac{dx}{dt} = P(x)$.

Introducing the $\ln$-coordinates $x = e^y$ so that the system takes the form:

$$ (2.11) \qquad \frac{dy}{dt} = (Ae^y + r), $$

where $y \in \ln M^+$ and the non interior equilibria are at the infinity.

We denote by $X(x)$ the vector field defined by (1.1) where $x$ can be taken as the normalized coordinates.

**2.2.1. Antibiotic or probiotic agent.** For a single antibiotic or probiotic agent the control system takes the form

$$ (2.11) \qquad \frac{dx}{dt} = X(x) + uY(x), $$
with $Y(x) = (\mathrm{diag}x)(\varepsilon_1, \dots, \varepsilon_n)^\top$ is the sensitivity vector and the control $u(t)$ describes the dosing regimen. Note that $u$ can be restricted to domain $[0, 1]$ using normalizations.

This notation can be applied to the hierarchy of models so that this leads to an hierarchy of control problems which can be analyzed independently. For instance, prior to infection one can analyze the effect of probiotics agents where $\varepsilon_i > 0$, while posterior to infection one can consider the effect of antibiotic treatment.

The set of admissible controls fits in the frame of permanent control, where $u(\cdot)$ is a measurable mapping on $[0, t_f]$ valued in $[0, 1]$. Each measurable bounded mapping can be approximated by a sequence of piecewise constant mappings in the $L^\infty$-topology and accessibility can be studied restricting to this class. This is the point of view of geometric control which leads to introduce the polysystem: $D = \{X + uY; u \text{ constant in the interval } [0, 1]\}$.

**2.2.2. Fecal transplantation and bactericide.** In the case of transplantation or bactericide the control system takes the form:

$$ \frac{dx}{dt} = X(x) + uY(x), $$

where $Y = (v_1, \dots, v_n)^\top$ is a constant vector field and $u$ takes its values in the whole $\mathbb{R}^+$, where its action being to get Dirac pulse $\delta(t - t')$, defined as the limit of controls sequence: $u_n = \lambda n$ on $[t', t' + 1/n]$.

**2.2.3. Optimal control problem in the permanent case vs the sampled-data frame.** Resuming the previous discussion, one can write the general control system in the form:

$$ \frac{dx}{dt}(t) = X(x(t)) + \sum_{ant.,pro.} u_i(t)Y_i(x(t)) + \sum_{transp.,bac.} u'_i(t)Y'(x(t)), $$

where the two sums are related respectively to probiotic, antibiotic agents and transplantations and bactericides. Moreover discontinuity in the state $x = (x_1, \dots, x_n)^\top \to x' = (x'_1, x_2, \dots, x_n)^\top$ can be understood as the jump action $x \to x'$. Hence this leads in the general case to a mixture of permanent and sampled-data control systems. The first action is related to permanent control, but due to logistic medical constraints it can fit in the sampled-data frame, e.g. a finite number of medical interventions at some predefined times: $0 < t_1 < \dots < t_k < t_f$ to modify the treatment. The second action fits in the sampled-data frame, since in particular it corresponds to invasive therapies.

In this article we concentrate to the permanent case associated to a single antibiotic or probiotic agent.

**OCP in the permanent case.** The problem is either to reduce the $x_1$-infection or to increase the production or ratio of healthy agents, prior to infection. This leads to consider in a dual formulation problems of the form:

* Reach in minimum time $t_f$ a given terminal manifold $N$ of codimension one for the control system written as: $\frac{dx}{dt} = X(x) + uY(x)$, $Y$ being associated to a specific treatment.

In this formulation candidates as minimizers are selected using the Maximum Principle stated in the introduction.
**2.3. Singular trajectories and time optimal control.** In this section, we make a brief recap of the properties of singular trajectories, crucial in our analysis, for full details see [4].

**Definition 2.3.** Let $(X, Y)$ the pairs of $C^\omega$ (real analytic) vector fields on $M$ associated to control system: $\frac{dx}{dt} = X(x) + uY(x)$. The feedback group $G_f$ is the set of triples $(\varphi, \alpha, \beta)$, where $\varphi$ is a local diffeomorphism and $u = \alpha(x) + \beta(x)v$, $\beta \neq 0$ is a feedback, the group structure being induced by the actions:

* local diffeomorphism $\varphi: (X, Y) \to (\varphi * X, \varphi * Y)$, where $\varphi * Z = d\varphi(Z \circ \varphi^{-1})$ denotes the image of a vector field $Z$.

* feedback: $u = \alpha(x) + \beta(x)v: (X, Y) \to (X + Y\alpha, Y\beta)$.

The control system can be lifted on the cotangent bundle $T^*M$ with symplectic structure defined by $d\omega$, where $\omega = pdx$ is the Liouville form. The Hamiltonians $H_X(z) = p \cdot X(x)$, $H_Y(z) = p \cdot Y(x)$, where $z = (x, p)$ are the symplectic coordinates, are the Hamiltonian lifts of $X, Y$. The system lift takes the form: $\frac{dz}{dt} = H_X(z) + uH_Y(z)$, and $H_X(z) + uH_Y(z)$ is the pseudo or non maximized Hamiltonian. One can lift every local diffeomorphism $\varphi$ into a Mathieu symplectomorphism $\varphi$ defined by: $x = \varphi(y)$, $p = q(\frac{\partial \varphi}{\partial y})^{-1}$, where $p, q$ are row vectors. This induces an action of $G_f$ on the pairs $(H_X, H_Y)$.

**2.3.1. Computations of singular extremals.** Relaxing the control bound to $u \in \mathbb{R}$, from the Maximum Principle the candidates as time minimizers are the so-called singular extremals control-trajectory pairs written shortly $(z, u)$ solutions of the constrained Hamiltonian dynamics: $\frac{dz}{dt}(t) = H(z(t))$, $H$ being the Hamiltonian vector field, the constraints coming from the maximization condition: $\frac{\partial H}{\partial u} = H_Y(z) = 0$.

Hence they are solutions contained in the switching set $\Sigma: H_Y(z(t)) = 0$. They can be computed, deriving this equation with respect to $t$. Introducing the Poisson bracket of $H_{Z_1}(z) = p \cdot Z_1(x)$, $H_{Z_2}(z) = p \cdot Z_2(x)$, by $\{H_{Z_1}, H_{Z_2}\}(z) = p \cdot [Z_1, Z_2](x)$.

Hence we deduce:

$$ (2.12) \qquad \begin{aligned} H_Y(z(t)) &= \{H_Y, H_X\}(z(t)) = 0, \\ \{\{H_Y, H_X\}, H_X\}(z(t)) &+ u(t)\{\{H_Y, H_X\}, H_Y\}(z(t)) = 0. \end{aligned} $$

We introduce the following.

**Definition 2.4.** The Generalized Legendre-Clebsch condition (GLC) along a singular extremal $(z(\cdot), u(\cdot))$ on $(0, t_f]$ is given by:

$$ \frac{\partial}{\partial u} \frac{d^2}{dt^2} \frac{\partial H}{\partial u}(z(t)) = \{\{H_Y, H_X\}, H_Y\}(z(t)) \neq 0, $$

for every $t$ in $[0, t_f]$. The switching surface is $\Sigma: H_Y(z) = 0$ and we denote $\Sigma'$ the subset: $H_Y(z) = \{H_Y, H_X\}(z) = 0$. Then outside the collinear set of $Y(x)$ and $[Y, X](x)$, if $\{\{H_Y, H_X\}, H_Y\} \neq 0$, the restriction of the symplectic form $d\omega$ to $\Sigma'$ defines a symplectic manifold $(M', \omega|_{M'})$.

This gives the following.

**Proposition 2.5.** Assume that the GLC-condition holds along $(z(t), u(t))$ then the extremal is called of minimal order. We have:
1. The singular control $u_s(\cdot)$ is the dynamic feedback: $u_s(z) = -\frac{\{H_Y, H_X\}, H_X\}(z) \{H_Y, H_X\}, H_Y\}(z)$.

2. Introduce the true Hamiltonian $H_s(z) := H_X(z) + u_s(z) H_Y(z)$, the singular extremals of minimal order are the solutions of $H_s(z)$ contained in the set $\Sigma'$.

They are the solutions of $H_s(z)$ restricted to the symplectic manifold $M'$.

Higher-order singular extremals can be determined at any order using the following algorithm. If

$$ \{\{H_Y, H_X\}, H_Y\}(z(t)) = 0 \text{ and } \{\{H_Y, H_X\}, H_X\}(z(t)) = 0 $$

then, deriving both relations one gets:

$$ \{\{\{H_Y, H_X\}, H_Y\}, H_X\}(z(t)) + u(t) \{\{\{H_Y, H_X\}, H_Y\}, H_Y\}(z(t)) = 0, $$

$$ \{\{\{H_Y, H_X\}, H_X\}, H_X\}(z(t)) + u(t) \{\{\{H_Y, H_X\}, H_X\}, H_Y\}(z(t)) = 0. $$

If the control cannot be derived from the previous equations, we repeat the derivation procedure.

### 2.3.2. Singular extremals as feedback invariants.

**Definition 2.6.** Let $E, F$ be two vector spaces and $G$ a group acting linearly on $E, F$. An homomorphism $\chi: G \to \mathbb{R} \setminus \{0\}$ is called a character. Let $\chi$ be a character, a semi-invariant of weight $\chi$ is a map $\lambda: E \to \mathbb{R}$ such that for all $x \in E, g \in G, \lambda(g \cdot x) = \chi(g)\lambda(x)$. It is called an invariant if $\chi = 1$. A map $\lambda: E \to F$ is a semi-covariant of weight $\chi$ if for all $x \in E, g \in G, \lambda(g \cdot x) = \chi(g)g \cdot \lambda(x)$. It is called a covariant if $\chi = 1$.

The following theorem is unveiled in [4, Theorem 13, p.103]

**Theorem 2.7.** Denote by $\lambda_s$ the map which associates to the pair $(X, Y)$ the Hamiltonian vector field $H_s$ restricted to $\Sigma'$. Then it is a covariant for the respective actions of the feedback group. In particular singular extremals are feedback invariants.

### 2.3.3. High order Maximum Principle [18].

**Proposition 2.8.** Assume $p$ is oriented using the convention of the Maximum Principle along the singular extremal $z(\cdot): H_X(z(t)) \ge 0$. Then a necessary time optimality condition on $[0, t_f]$ is given by

$$ \frac{\partial}{\partial u} \frac{d^2}{dt^2} \frac{\partial H}{\partial u}(z(t)) = \{\{H_Y, H_X\}, H_Y\}(z(t)) \ge 0. $$

**Definition 2.9.** The singular extremal $z(t) = (x(t), p(t))$, $t \in [0, t_f]$ is called strict if $p$ is unique up to a scalar.

**Corollary 2.10.** Assume the strict case. Then the singular trajectories projections of singular extremals of minimal order are stratified according to the following:

* Hyperbolic case: $H_X \cdot \{\{H_Y, H_X\}, H_Y\}(z) > 0$,

* Elliptic case: $H_X \cdot \{\{H_Y, H_X\}, H_Y\}(z) < 0$,

* Abnormal or exceptional case: $H_X(z) = 0$.

### 2.3.4. Applications.
**2d-case.** Singular extremals satisfy $H_Y = \{H_Y, H_X\} = 0$ so that singular trajectories are located on the set $\mathcal{S}$: $\det(Y(x), [Y, X](x)) = 0$. Outside the collinear set $\mathcal{C}$: $\det(Y(x), X(x)) = 0$ one can take $Y, X$ as a frame and writing $[[Y, X], Y](x) = \alpha(x)X(x) + \beta(x)Y(x)$. The singular control is given by: $[[Y, X], X](x) + u_s(x) [[Y, X], Y](x) = 0$. Hyperbolic case corresponds to $\alpha(x) > 0$ and elliptic case to $\alpha(x) < 0$.

**3d-case.** The 3d-case is already a very rich situation to analyze the singular extremals and the program goes as follows.

Introduce the following determinants:

$$ (2.13) \qquad \begin{aligned} D(x) &= \det(Y(x), [Y, X](x), [[Y, X], Y](x)), \\ D'(x) &= \det(Y(x), [Y, X](x), [[Y, X], X](x)), \\ D''(x) &= \det(Y(x), [Y, X](x), X(x)), \end{aligned} $$

and using the relations

$$ (2.14) \qquad \begin{aligned} H_Y(z) &= \{H_Y, H_X\}(z) = 0, \\ \{\{H_Y, H_X\}, H_X\}(z) + u_s\{\{H_Y, H_X\}, H_Y\}(z) &= 0, \end{aligned} $$

we can eliminate $p$ and the singular control is given by the feedback:

$$ u_s(x) = -\frac{D'(x)}{D(x)}. $$

**Lemma 2.11.** The action of the feedback group implies the following changes on $D, D'$ and $D''$:

$$ D^{\phi*X,\phi*Y}(x) = \det\left(\frac{\partial\phi}{\partial x}\right) D^{X,Y}(\phi^{-1}(x)), \quad D^{X+\alpha Y,\beta Y}(x) = \beta^4 D^{X,Y}(x), $$

$$ D^{\phi*X,\phi*Y}(x) = \det\left(\frac{\partial\phi}{\partial x}\right) D^{X,Y}(\phi^{-1}(x)), \quad D^{X+\alpha Y,\beta Y}(x) = \beta^3 (D^{X,Y}(x) + \alpha D^{X,Y}(x)), $$

$$ D^{\phi*X,\phi*Y}(x) = \det\left(\frac{\partial\phi}{\partial x}\right) D^{X,Y}(\phi^{-1}(x)), \quad D^{X+\alpha Y,\beta Y}(x) = \beta^2 D^{X,Y}(x). $$

In particular the surfaces $D = 0$ and $D'' = 0$ are feedback invariant.

Defining the vector field:

$$ (2.15) \qquad X_s(x) := X(x) + u_s(x)Y(x), $$

we have:

**Proposition 2.12.** Singular trajectories of minimal order stratified the dynamics into:

* Hyperbolic arcs in $DD'' > 0$,

* Elliptic arcs in $DD'' < 0$,

* Exceptional or abnormal arcs in $D'' = 0$.
From Lemma 2.11, we obtain the following proposition.

**Proposition 2.13.** For the action of the feedback group $G_f = \{(\varphi, \alpha, \beta)\}$ reducing to $\varphi$-changes of coordinates on vector fields $X_s$, the map $\lambda_s : (X, Y) \to X_s$ is a covariant. Hence this allows to generate feedback invariants using the dynamics (2.15).

Equilibria of this dynamics split into two types:

* If $D \neq 0$, they are given by the solutions of $X_s(x) = 0$.

* If $D = 0$, one can reparameterize the dynamics to get the vector field $D(x)X(x) - D'(x)Y(x)$ and additional (non isolated) singular points are located on $D(x) = D'(x) = 0$.

Exceptional trajectories are contained in the invariant set $D''(x) = 0$ for the dynamics.

### 3. Extremals classification and local time minimal syntheses near a terminal manifold of codimension one.

**3.1. Introduction and definitions.** In this section we consider the time minimal control problem for the system: $\frac{dx}{dt} = X(x) + uY(x)$, $|u| \le 1$, with terminal manifold $N$ of codimension one. We denote by $H(z, u) = H_X(z) + uH_Y(z)$ the pseudo-Hamiltonian and $M(z) = \max_{|u| \le 1} H(z, u)$ the true or maximized Hamiltonian. In this setting the control domain is taken as $U = [-1, +1]$, where extreme control values lead to introduce the vector fields $X(x) - Y(x)$ (no medical treatment) or $X(x) + Y(x)$, that is the maximal dosing regimen.

**Definition 3.1.** The extremals are concatenation of regular extremals for which almost everywhere

$$u(t) = \operatorname{sign}H_Y(z(t))$$

and singular extremals if

$$H_Y(z(t)) = 0$$

holds identically.

An extremal is called exceptional if the maximized Hamiltonian is such that $M(z) = 0$. A BC-extremal is an extremal satisfying the transversality condition. A switching time is an instant such that the extremal control is discontinuous. A bang-bang extremal is a regular extremal with a finite number of switches.

Since the control domain is $U = [-1, +1]$, a singular extremal is called strictly feasible (admissible) if $|u_s| < 1$ and saturating at time $t$ if $|u_s(t)| = 1$. A regular or singular extremal is called strict if $p$ is unique up to a scalar. In the strict case singular extremals can be classified into hyperbolic, elliptic and exceptional extremals. We denote by $\sigma_+$, $\sigma_-$ bang arcs associated respectively to $u = +1$ or $u = -1$, $\sigma_s$ is a singular arc associated to $u_s$. We denote by $\sigma_1\sigma_2$ an arc $\sigma_1$ followed by $\sigma_2$.

**Definition 3.2.** Taking an open set $V$ of $M$, the problem is called geodesically complete on $V$ if for each $x_0, x_1 \in V$ there exists a time minimizing geodesic in $V$ joining $x_0$ to $x_1$. Fixing the target to $N (=N \cap V)$, a time minimal synthesis is a discontinuous feedback $x \to u^*(x)$ so that the solutions of $\frac{dx}{dt} = X(x) + u^*(x)Y(x)$ are well defined and $u^*(x)$ is the optimal feedback solution to steer $x$ to the target in minimum time.
**3.2. Small time classification of regular extremals.** In this section we recall some basic properties of regular extremals, see [19] but also [20] for the analysis of the Fuller phenomenon as a recommended reading.

The surface $\Sigma: H_Y(z) = 0$ is called the switching surface and we denote by $\Sigma'$ the set $H_Y(z) = \{H_Y, H_X\}(z) = 0$. Let $z(\cdot) = (x(\cdot), p(\cdot))$ be a reference extremal on $[0, t_f]$. We note $\Phi(t) := H_Y(z(t))$ the switching function coding the switching times.

Deriving twice $\Phi$ with respect to time, one gets:

$$ \frac{d\Phi}{dt}(t) = \{H_Y, H_X\}(z(t)), $$

$$ (3.1) \qquad \frac{d^2\Phi}{dt^2}(t) = \{\{H_Y, H_X\}, H_X\}(z(t)) + u(t)\{\{H_Y, H_X\}, H_Y\}(z(t)). $$

**Definition 3.3.** The time $t$ is called an ordinary switching time if $\Phi(t) = 0$ and $\frac{d\Phi}{dt}(t) \neq 0$.

**Lemma 3.4.** Assume $t$ be an ordinary switching time, then near $z(t)$ every extremal projects onto:

* $\sigma_+\sigma_-$ if $\frac{d\Phi}{dt}(t) > 0$,

* $\sigma_-\sigma_+$ if $\frac{d\Phi}{dt}(t) < 0$.

**Definition 3.5.** Let $z(\cdot)$ be a bang extremal on $[0, t_f]$ with $u = \varepsilon \in \{-1, +1\}$. We note by $\frac{d^2\Phi}{dt^2}$ the expression (3.1) in which $u \equiv \varepsilon$. The point $z(t)$ is called a fold point if $\Phi(t) = \frac{d\Phi}{dt}(t) = 0$ and $\frac{d^2\Phi}{dt^2}(t) \neq 0$. Assume that $\Sigma'$ is a regular surface of codimension two. We have three cases:

* Parabolic case: $\frac{d^2\Phi_+}{dt^2}(t) \frac{d^2\Phi_-}{dt^2}(t) > 0$,

* Hyperbolic case: $\frac{d^2\Phi_+}{dt^2}(t) > 0$ and $\frac{d^2\Phi_-}{dt^2}(t) < 0$,

* Elliptic case: $\frac{d^2\Phi_+}{dt^2}(t) < 0$ and $\frac{d^2\Phi_-}{dt^2}(t) > 0$.

This leads to:

**Proposition 3.6.** In the neighborhood of a fold point every extremal projects onto:

* In the parabolic case: $\sigma_+\sigma_-\sigma_+$ or $\sigma_-\sigma_+\sigma_-$.

* In the hyperbolic case $\sigma_\pm\sigma_s\sigma_\pm$.

* In the elliptic case, every extremal is bang-bang but the number of switches is not uniformly bounded.

This is illustrated in Fig. 3.

Note that in the elliptic case, there is a foliation by a cylinder of the regular extremal dynamics and the number of switches is related to the distance to $\Sigma'$ versus the Fuller phenomenon, where the sequence of regular arcs is not bang-bang and converges to $\Sigma'$.

**Application to the 3d-case.** Consider the case where $Y$, $[Y, X]$ and $X$ form a frame i.e.

$$ D'' = \det(Y, [Y, X], X) $$

is not vanishing. The problem is strict for singular extremals since $Y$ and $[Y, X]$ are independent. Hyperbolic trajectories are small time minimizing trajectories, while elliptic trajectories
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Fold case.</figcaption></figure>

are small time maximizing, this up to the first conjugate time $t_{1c}$ computed in [5] if they are strictly admissible (even in the limit case with no constraints on the control). In the parabolic case, they can be absent or not feasible. Consider the separating case where the singular arc is exceptional (abnormal) and assume that it is strictly admissible (for instance relaxing the control bound to the whole $\mathbb{R}$). In this case using again [5] an exceptional arc $\sigma_s$ is time minimizing and time maximizing, up to the first conjugate time $t_{1c}$. Such a point is absent in the 3d-case. Such an arc can be lifted into two extremals ($\pm p, \sigma_s$) and it corresponds either to an hyperbolic or elliptic situation in $\Sigma'$. One contribution of [21] is to analyze the time minimal syntheses near the terminal manifold in this situation.

**3.3. General concepts of regular synthesis with a terminal manifold of codimension one.** Take a triple $(X, Y, N)$ and let $x_0 \in N$ which can be identified to 0 while $N$ is the plane $x_1 = 0$. Let $U$ be a neighborhood of 0, which divides the space into neighborhoods $V$ and $W$ contained respectively in $x_1 < 0$ and $x_1 > 0$ so that $U = V \cup W$. The problem is to compute the time minimal regular synthesis to steer each point of $U$ to the terminal manifold.

$N$ can be taken locally as $f^{-1}(0)$ where $f$ is a submersion from $U$ into a neighborhood of 0 in $\mathbb{R}$. The set of triples $(X, Y, f)$ is endowed with the $C^\infty$-Whitney topology and we denote by $j^k X$ (resp. $j^k Y, j^k f$) the $k$-jet of $X$ (resp. $Y, f$) obtained by taking the Taylor expansion at $x_0 = 0$ up to order $k$. We say that the triple $(X, Y, f)$ has at 0 a singularity of codimension $i$ if $(j^k X, j^k Y, j^k f)$ belongs to a semi-algebraic set of codimension $i$ in the jet space.

Our aim is to make a short presentation of the results of [6, 21, 7] to classify local syntheses up to an homeomorphim preserving the target $N$ for all cases of codimension $\le 2$, by considering two cases occurring in the application that we introduce next.

**Definition 3.7.** When $Y$ is everywhere tangent to the target $N$, the control $u$ is indirect and this case is called the flat case. In the non flat case, the action of the control is direct and the
set of points where $Y$ is tangent to $N$ is of codimension $\ge 1$.

**3.3.1. Stratified synthesis.** Our aim is to describe the local time minimal synthesis in a neighborhood of $N$ by estimating at any order the different strata. Actually, the optimal control feedback $u^*(x)$ is not always defined on the whole subset $V$ of $U$ in the domain $x_1 < 0$, since for some $x \in V$ the target $N$ is not accessible.

In fact we can reduce our study to two cases:

• The case when the convex cone $C$ generated by $\{X \pm Y\}$ is strict and the set of admissible directions points towards the space $x_1 > 0$.

• The exceptional case where the set of admissible directions are tangent to the terminal manifold.

This leads to introduce the exceptional locus in the construction of the stratification of $N$.

**Definition 3.8.** Let $n$ be the normal to $N$ oriented toward $x_1 > 0$. The exceptional locus (restricted to $N$) $\mathcal{E}$ is the set of points of $N$ such that: $n(x) \cdot Y(x) = n(x) \cdot X(x) = 0$.

The second part of the stratification amounts to introduce the singular locus.

**Definition 3.9.** The singular locus (restricted to $N$) $\mathcal{S}$ is the set of points of $N$ such that: $n(x) \cdot Y(x) = n(x) \cdot [Y, X](x) = 0$.

**Definition 3.10.** A stratified synthesis amounts to find a partition of $V$ (or a partition of $U$ in the exceptional case) into $V^+$ (resp. $U^+$) where $u^*(x) = +1$ and $V^-$ (resp. $U^-$) where $u^*(x) = -1$ and a stratified surface separating $V^+$ and $V^-$ (resp. $U^+, U^-$) with three kind of strata:

• Switching locus. It is the closure of the set or ordinary switching points and forming the set $W_\#$, $\# \in \{-1, +1\}$, where $W_+$ is associated to $\sigma_+\sigma_-$ and $W_-$ to $\sigma_-\sigma_+$.

• Cut locus. Let $\sigma : [t_f, 0] \to M$ be a minimizing curve, integrating backwards from $N$ so that $\sigma(0) \in N, t_f < 0$. The cut-locus is the closure of the set of points where optimality is lost. It is denoted $C$ and contains the splitting locus $L$ where the optimal feedback is not unique.

• Switching singular locus. It is foliated by optimal singular arcs and is denoted $\Gamma_s$. Recall that if $u_s \in [-1, +1[$ the singular arc is strictly feasible but it can be saturated if $u_s^*(x) \in \{-1, +1\}$.

To estimate the different strata we use semi-normal forms restricting the action of the feedback group to local diffeomorphisms $\varphi$ preserving $0$ and feedbacks $u \to -u$ so that $\sigma_+$ and $\sigma_-$ can be inverted in the classification.

We can choose local coordinates to normalize a reference trajectory to $t \mapsto (t, 0, \dots, 0)$, the terminal manifold $N$ and the controlled vector field $Y$ to compute the optimal synthesis in a neighborhood of $N$. More precisely, we introduce the following important concept.

**Definition 3.11.** Let $x_0 \in N$ which can be identified to $0$. A Whitney chart is a pair $(U, \varphi)$ where $U$ is a neighborhood of $x_0$ and $\varphi$ a system of coordinates $(x, y, y_1, \dots, y_{n-3}, z)$ such that:

1. $Y = \frac{\partial}{\partial z}$,

2. $N$ is the surface of codimension one parameterized by $(ks^2 + s^2O(|w_1, w_2, \dots, w_{n-3}, w, s|), w, w_1, \dots, w_{n-3}, s)$,

3. The time minimal synthesis in $U$ is $C^0$ described with foliations by 2d or 3d syntheses
542

543

544

545

with triples $(X^b(x), Y^b(x), N^b(x))$, $x \in \mathbb{R}^2$ or $\mathbb{R}^3$ where the 3d-cases occur only in the exceptional case if $n \ge 3$.

Moreover in $U$, one can construct the stratified optimal synthesis, where each strata can be estimated at any order.

546

547

**Definition 3.12.** In the previous construction the restriction from $n$ to $n-1$ decomposing $x = (x', \lambda)$, $N = \cup_{\lambda} N'(\lambda)$, where $\lambda \in \mathbb{R}$, is a parameter is called a one dimensional unfolding.

548

The previous definitions will be clarified in the examples we present next.

549

550

551

552

553

**3.3.2. Main points of the geometric and analytic construction.** One can restrict our presentation to the 2d-case. Let $n(x)$ be the normal to $N$, oriented towards $x_1 > 0$. We denote by $N^{\perp}$ the symplectic lift of $N$: $N^{\perp} = \{(x, p); x \in N, p = n(x)\}$. The stratification of $N$ by the conditions $n(x) \cdot Y(x) = 0$ and $n(x) \cdot [Y, X](x) = 0$ selects ordinary switching points or fold points classified in Proposition 3.6.

554

555

556

557

558

* Ordinary switching locus $K$: One can restrict for simplicity the analysis to the 2d non flat case. Let $x_0 \in N$ such that $n \cdot Y(x_0) = 0$ and both $n \cdot X(x_0)$ and $n \cdot [Y, X](x_0)$ non zero. This leads to compute a switching locus $K$ terminating at $x_0$. Such a switching locus is part of $W$ if the corresponding bang-bang extremal crosses the switching locus versus reflects on the switching locus. This leads to estimate the slope of $K$.

559

560

561

562

563

564

565

566

567

* Singular locus $\Gamma_s$: Consider again the non flat case. Let $x_0$ such that $(n(x_0), x_0))$ is a fold point hence $n \cdot Y(x_0) = n \cdot [Y, X](x_0) = 0$. Moreover assume that $n \cdot X(x_0) \ne 0$ (non exceptional case) and that the singular arc $\sigma_s$ terminating at $x_0$ is strictly admissible, such an arc being small time minimizing for the problem with fixed extremities (assuming the condition $n \cdot [[Y, X], Y](x_0) \ne 0$). One can choose coordinates $(x, y)$ such that $Y = \frac{\partial}{\partial y}$ and $N$ is given by $(\frac{1}{2}ks^2, s)$, while $X$ can be normalized to $(1 - y^2X_1(x, y)) \frac{\partial}{\partial x} + (u - u_s(x) + yX_2(x, y)) \frac{\partial}{\partial y}$, where $\sigma_s$ is identified to $t \to (t, 0)$ and $u_s$ is the singular control. Note that such normalizations were introduced in [5] in a more general context.

568

569

570

571

In order to decide about optimality one can compare by direct computations in the normalized coordinates the curvature of the boundary of the accessibility set along the singular arc with the curvature of $N$ given by $k$. This leads to the cases described in Fig. 4. In the first case, there exists a cut locus $C$ with slope $-u_s(0)$.

572

573

574

The 3d-case with an hyperbolic point terminating at $x_0 = 0$ can be obtained by constructing a semi-normal form in coordinates $(x, y, z)$ by unfolding the case $k > 0$ and this leads to Fig. 5.

575

576

577

578

579

The central picture corresponds to Fig. 4 with $k > 0$. Left we have represented a Whitney neighborhood containing the optimal switching locus $W_+$ intersecting the hyperbolic singular locus. They can be easily computed using the stratification of $N^{\perp}$ into sectors with ordinary switching points and fold points while the singular arc terminating at 0 corresponds to an hyperbolic fold point.

580

581

582

This gives the procedure to compute the optimal syntheses using [6] in the non exceptional case as unfolding of 2d-cases. The procedure fails in the exceptional case for $n \ge 3$, where the following example from [21] shows why it cannot be reduced to 2d-foliations.
<figure><img src="image_5.png" /><figcaption>Figure 4. Optimal synthesis in the 2d-hyperbolic case.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Unfolding in the 2d-hyperbolic case with parameter $w_0$.</figcaption></figure>

**3.3.3. An exceptional 3d-case not 2d-reducible.** One take a flat case so that $N$ can be identified to $(0, w, s)$ in $(x, y, z)$ coordinates and $Y = \frac{\partial}{\partial z}$. In the construction the main point is to take a bang arc $\sigma_-$, which is optimal in the domain $x \ge 0$ with a contact of order 3 at 0 with the surface $N$. Hence this gives birth in the domain $x > 0$ to arcs $\sigma_-$ intersecting three times the target $N$, thanks to contact analysis.

Such a situation occurs for instance for the model which in fine is a $C^0$-normal form describing the situation:

$$ \frac{dx}{dt} = z, \quad \frac{dy}{dt} = b, \quad \frac{dz}{dt} = 1 + u + y, $$

the target $N$ given by $(w, s) \to (0, w, s)$.

One considers the situation with $b > 0$, where each point of the neighborhood $U$ can be steered to the target.

In the domain $x < 0$, every time optimal trajectory is of the form $\sigma_+$ and the contact of $\sigma_+$ at 0 with $N$ is of order 2.

Optimal arcs $\sigma_-(t)$ included in $x \ge 0$ are satisfying:

$$ (3.2) \qquad x_-(t) = t(s + wt/2 + bt^2/6 + \dots), $$
<figure><img src="image_7.png" /><figcaption>Figure 6. Synthesis exceptional case when $b < 0$. $\sigma_-$ has a contact of order 3 with $N$.</figcaption></figure>

where the weight of $s$ is one and the weight of $w$ is two, neglected having weights greater than 3.

The curve $\Xi$ is the set of points $(0, w, s)$ of $N$ such that: $\frac{x_-(t)}{t} = \frac{d}{dt} \left( \frac{x_-(t)}{t} \right)$ have a common zero and is given using (3.2) by: $w^2 \sim \frac{8bs}{3}$.

The optimal synthesis is represented on Fig. 6. Optimal arcs $\sigma_-$ in the domain $x \ge 0$ are cutting twice $N$ in the subsets of the target denoted 2 and 2' and have two subarcs which are optimal in $x > 0$, but the subarc in $x < 0$ is not optimal.

**4. Computations and preliminary results on the Controlled Lotka–Volterra model.** The aim of this section is to present the geometric study of the controlled Lotka–Volterra model:

$$ \frac{dx}{dt}(t) = X(x(t)) + u(t)Y(x(t)), $$

with $X(x) = (\text{diag}x)(Ax+r)$ and $Y(x) = (\text{diag}x)\epsilon$, where $x = (x_1, \dots, x_n)^\top$ is the population species, $x_1$ represents the infected species and $\epsilon = (\epsilon_1, \dots, \epsilon_n)^\top$.

The system can be written in $\ln$-coordinates: $y = \ln x$ and it takes the form:

$$ (4.1) \qquad \frac{dy}{dt} = (Ae^y + r) + u \epsilon. $$

**4.1. Equilibria and the collinear set.** The collinear set $\mathcal{C}$ is one of the main feedback invariant, related to computations of the free equilibria of the system for the hierarchy of models with no treatment $u = 0$ or maximal dosing treatment $u = 1$, but also for all intermediate dosing.

This set is a one dimensional algebraic variety and is the projection of the set:

$$ (4.2) \qquad \{(x_e, u_e) \in \mathbb{R}^n \times \mathbb{R}; \exists u_e, X(x_e) + u_e Y(x_e) = 0\}. $$

At such a point introduce the Jacobian matrix:

$$ J(x_e, u_e) = \frac{\partial}{\partial x|_{(x_e, u_e)}} (X(x) + u_e Y(x)). $$

The following dimensionless coordinates are useful in the computations.
The system can be written:

$$621 \qquad (4.3) \qquad \frac{dx_i}{dt} = x_i r_i - x_i \sum_{j=1}^{n} a_{ij} x_j + u x_i \varepsilon_i$$

and denotes by $x^*$ the free equilibrium given by:

$$624 \qquad r_i - \sum_{j=1}^{n} a_{ij} x_j^* = 0, i = 1, \dots, n.$$

One sets $v_i = \frac{x_i}{x_i^*}$ so that the dynamics takes the form:

$$626 \qquad \frac{dv_i}{dt} = v_i r_i - v_i \sum_{j=1}^{n} a_{ij} x_j^* v_j + u v_i \varepsilon_i,$$

and denoting $a_{ij}^* = a_{ij} x_j^*$, one has $r_i = \sum_{j=1}^n a_{ij}^*$ since the interior equilibrium is normalized to $\Omega = (1, \dots, 1)$.

If we set: $x_i = v_i - 1$, the dynamics is given by:

$$630 \qquad \frac{dx_i}{dt} = (x_i + 1) \sum_{j=1}^{n} a_{ij}^* - (x_i + 1) \sum_{j=1}^{n} a_{ij}^* (x_j + 1) + u (x_i + 1) \varepsilon_i.$$

Hence we have:

**Proposition 4.1.** In the dimensionless coordinates the controlled Lotka–Volterra model is given:

$$632 \qquad \frac{dx_i}{dt} = -(x_i + 1) \sum_{j=1}^{n} a_{ij} x_j + u (x_i + 1) \varepsilon_i,$$

so that $\Omega$ is identified to 0 and the Jacobian matrix at 0 is $-A$.

**4.1.1. Computations in the 2d–case in the dimensionless coordinates.** We consider the regular 2-dimensional dynamics given in Proposition 4.1.

**Collinear set and classification of equilibria.** The collinear set is one of the main feedback invariant related to the computations of free equilibria with no treatment $u = 0$ and forced equilibria with maximal dosing $u = 1$.

The collinear set is given by the determinantal variety: $\det(X(x), Y(x)) = 0$:

$$642 \qquad (x_1 + 1)(x_2 + 1)(x_1 \kappa_2 - x_2 \kappa_1) = 0,$$

where $\kappa_1 = \varepsilon_1 a_{22} - \varepsilon_2 a_{12}$ and $\kappa_2 = \varepsilon_2 a_{11} - \varepsilon_1 a_{21}$.

Alternatively, it can be viewed as the one dimensional algebraic variety projection of the set

$$646 \qquad \{(x, u_e) \in \mathbb{R}^2 \times \mathbb{R}, \exists u_e \text{ such that } X(x) + u_e Y(x) = 0\}.$$
The condition $u_e \in [0, 1]$ selects a segment of persistent equilibria located on the line:

$$C := \left\{ x_2 = u_e(x_1) \frac{\kappa_2}{\det A}, u_e(x_1) = x_1 \frac{\det A}{\kappa_1} \in [0, 1] \right\}.$$

In particular, introduce $x_e = (x_{1e}, x_{2e}) \in C$, the origin $x_e = 0$ is associated to the control $u_e = 0$ while $x_e = (\kappa_1/\det A, \kappa_2/\det A)$ is associated to the control $u_e = 1$.

For $u_e \in [0, 1]$, define the Jacobian matrix:

$$J(x_e, u_e) = \frac{\partial}{\partial x} (X(x) + u_e Y(x))|_{x=x_e},$$

we have:

Lemma 4.2. Let $x_e \in C$ associated to the control $u_e$, the spectrum of $J$ is:

$$\mathrm{spec}(J(x_e, u_e)) = \left\{ \left( k(x_e) \pm \sqrt{k(x_e)^2 + k'(x_e)} \right) / 2 \right\},$$

where $k(x) = -x_1 a_{22} \kappa_2 / \kappa_1 - a_{22} - a_{11}(x_1 + 1)$ and $k'(x) = -4 \det A(x_1 + 1)(1 + x_1 \kappa_2 / \kappa_1)$.

## 4.2. Computation of the collinearity locus and properties. Construction of a normal form in ln-coordinates.

**Computations about $C$ in the $n$-dimensional case.** We have the following algorithm, taking the system represented in the $x$-coordinates.

* Step 1. The collinear set is the projection of the algebraic curve defined by: There exists $u_e$ constant such that $X(x_e) = -u_e Y(x_e)$. This gives $n$-equations depending upon $(n+1)$ variables $(x_e, u_e)$.

* Step 2. Take such a pair $(x_e, u_e)$ so that $x_e$ is a forced equilibrium for $u = u_e$ and they form a set with extreme points associated to $u_e = 0$ and $u_e = 1$, when restricting to $u_e \in [0, 1]$.

The linear dynamics at the point $x_e$ is characterized by the Jacobian matrix:

$$J(x_e, u_e) = \frac{\partial}{\partial x} \bigg|_{(x_e, u_e)} (X(x) + u_e Y(x))$$

the spectrum being $\sigma(J) = (\lambda_1, \dots, \lambda_n)$ with associated generalized eigenspaces $E_{\lambda_i}, i = 1, \dots, n$.

The linear stability is determined by this spectrum, thanks to Lyapunov linear stability theory.

* Step 3. From control point of view we have three cases:

1. $u_e \notin [0, 1]$: $u_e$ is not feasible,

2. $u_e \in ]0, 1[$: $u_e$ is strictly feasible,

3. $u_e = \pm 1$: $u_e$ being feasible but saturating.

One can discuss the linear controllability of the pair $(J(x_e, u_e), b)$ where $b = Y(x_e)$.

* Kalman condition: rank $[b, Jb, \dots, J^{n-1}] = n$ and the singular point $x_e$ is regular. If the control $u_e$ is strictly feasible local controllability holds [22].
• If rank $[b, Jb, \dots, J^{n-1}] = n-k < n$, then the singular point is a singular trajectory (reduced to a point) associated to $u_e$ and $k$ is the codimension of the singularity.

From linear controllability theory, one can construct a normal form, at a given equilibrium pair $(x_e, u_e)$.

We take $\ln$-coordinates so that $X(x)$ takes the form $X^{\flat}(y) = (Ae^y + r)$ and the controlled vector field $Y(x)$ becomes the constant vector field $Y^{\flat} = \epsilon$.

Let $(y_e, u_e)$ be the selected forced equilibrium in the $\ln$-coordinates and let $z = y - y_e$, $v = (u - u_e)$ so that the system takes the form:

$$ \frac{dz}{dt} = J(z) + R(z) + v \epsilon, $$

where $J$ is the Jacobian matrix at $(y_e, u_e)$.

One can find coordinates such that the linear dynamics decomposes into

$$ \frac{dz_1}{dt} = J_{11}z_1 + J_{12}z_2 + v\epsilon $$

$$ \frac{dz_2}{dt} = J_{21}z_2 $$

where the restriction to the controllable space $z_2 = 0$ is given by the dynamics:

$$ \frac{dz_1}{dt} = J_{11}z_1 + v \epsilon. $$

The pair $(J_{11}, \epsilon)$ can be set in Brunovsky canonical form:

$$ J_{11} = \begin{bmatrix} 0 & \mathrm{Id}_{n-k-1} \\ -a_1 & -a_2 \ \dots \ -a_{n-k} \end{bmatrix}, \quad \epsilon = \begin{bmatrix} 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix}, $$

where the coefficients $(a_1, \dots, a_{n-k})$ are the coefficients of the characteristic polynomial of $J_{11}$.

Finally this leads to the construction of a normal form:

$$ \frac{dz}{dt} = J(z) + R(z) + v \epsilon, $$

where the pair $(J(x_e), u_e), \epsilon)$ is in linear canonical form. Note in particular that the sensitivity vector $\epsilon$ is normalized to $(0, \dots, 0, 1)^T$ and $R(z)$ is the non linear part of the dynamics, related to singular trajectories computations.

**4.3. Singular extremals.** Lie brackets can be computed in the original coordinates but the computations are simpler in $\ln$-coordinates, since the control vector field $Y^{\flat}$ is constant. This comes from the following.

**Proposition 4.3.**
1. If $X = (\operatorname{diag}x)X_1(x)$, $Y = (\operatorname{diag}x)Y_1(x)$, where $X_1, Y_1$ are polynomial then the iterated Lie brackets are in the same category of polynomial vector fields of the form $(\operatorname{diag}x)P(x)$.

2. Let $\varphi$ be the diffeomorphism: $x = e^y$. Denote by $X^b, Y^b$ the images of $X, Y$ by this diffeomorphism. Then by invariance of the Lie bracket: $[X^b, Y^b](y) = d\varphi^{-1}[X, Y](e^y)$.

**4.3.1. Computations in the 2d-case.** We start by computing the $\ln$-coordinates:

$$ \bullet X^b = (r_1 + \sum_{j=1,2} a_{1j} e^{x_j}, r_2 + \sum_{j=1,2} a_{2j} e^{x_j}), $$

$$ \bullet Y^b = (\varepsilon_1, \varepsilon_2), $$

$$ \bullet [X^b, Y^b] = (\sum_{j=1,2} \varepsilon_j a_{1j} e^{x_j}, \sum_{j=1,2} \varepsilon_j a_{2j} e^{x_j}). $$

Hence the singular locus: $S: \det([X, Y](x), Y(x)) = 0$ is given by:

$$ (4.4) \qquad x_1x_2(\varepsilon_1x_1\kappa_2 - \varepsilon_2x_2\kappa_1) = 0 $$

which is stratified into $x_1 = x_2 = 0$ and a permanent straight-line $L: \varepsilon_1x_1\kappa_2 = \varepsilon_2x_2\kappa_1$.

Moreover:

$$ \bullet [[X^b, Y^b], Y^b] = (\sum_{j=1,2} \varepsilon_j^2 a_{1j} e^{x_j}, \sum_{j=1,2} \varepsilon_j^2 a_{2j} e^{x_j}). $$

Outside the collinear set $C$, $X, Y$ form a frame and writing:

$$ [[Y, X], X](x) = \alpha(x)X(x) + \beta(x)Y(x), $$

so that:

$$ \bullet \text{ Hyperbolic (feasible) subarcs are such that } \alpha(x) > 0, $$

$$ \bullet \text{ Elliptic subarcs are such that } \alpha(x) < 0. $$

$$ \bullet \text{ At the persistent point intersection of } S \text{ and } C, \text{ one gets an exceptional point.} $$

**4.3.2. Computations in the 3d-case.** The Lie brackets computations are as before, except that the index $j$ goes from 1 to 3. So that

$$ D = \det(Y^b, [Y^b, X^b], [[Y^b, X^b], Y^b]) $$

is homogeneous and quadratic with respect to the variable $e^{x_i}$.

Moreover the exceptional locus is given by the relation:

$$ D'' = \det(Y^b, [Y, X]^b, X^b) = 0, $$

this set being a quadratic non homogeneous variety with respect to the variables $e^{x_i}$.

Note that the computation of the Lie bracket: $[[Y^b, X^b], X^b]$ is more complex and formal computations are necessary.

This fixes the limit of the computational complexity in the $n$-dimensional case.

## 4.4. Applications.

**4.4.1. A 2d-working case.** In this section, the general techniques from [6, 21] will be applied to analyze a 2d-case study which occurs in the Lotka–Volterra model, see also [11] as a complementary study for the fixed end point problem.
Lie brackets and feedback invariants. The first step is to compute the collinear set $\mathcal{C}$ defined by $\det(X(x), Y(x)) = 0$, which takes the form in the original coordinates:

$$ (4.5) \qquad x_1x_2(x_1\kappa_2 - x_2\kappa_1 + \varepsilon_2r_1 - \varepsilon_1r_2) = 0, $$

while the singular locus is $\mathcal{S}: x_1x_2(\varepsilon_1x_1\kappa_2 - \varepsilon_2x_2\kappa_1) = 0$ from (4.4) and Lie brackets of length 3 are:

$$ [[Y, X], Y](x) = -x_1(\varepsilon_1^2 x_1 a_{11} - \varepsilon_2^2 x_2 a_{12}) \frac{\partial}{\partial x_1} - x_2(\varepsilon_1^2 x_1 a_{21} + \varepsilon_2^2 x_2 a_{22}) \frac{\partial}{\partial x_2}, $$

$$ \begin{aligned} [[Y, X], X](x) = & -x_1(\varepsilon_1 x_1 (r_1 a_{11} + x_2 a_{12} (a_{11} - a_{21})) + \varepsilon_2 x_2 a_{12} (x_1 (a_{21} - a_{11}) + r_2)) \frac{\partial}{\partial x_1} \\ & -x_2(\varepsilon_2 x_2 (r_2 a_{22} + x_1 a_{21} (a_{22} - a_{12}))) + \varepsilon_1 x_1 a_{21} (x_2 (a_{12} - a_{22}) + r_1)) \frac{\partial}{\partial x_2}. \end{aligned} $$

The geometric situation that we analyze in our working example is the exceptional case where we consider the intersection of the collinear locus with the singular locus, which corresponds from (4.5) and (4.4) to intersection of two straight-lines. It corresponds to a generic interaction between an hyperbolic arc and an elliptic arc.

Constructing a semi-normal form. The second step is to construct a semi-normal form for the system. The construction is detailed in [6] but computation are simple in the 2d-case. It consists to choose coordinates such that the intersection is taken as the origin (0,0), Y is identified to the constant vector field $Y = \frac{\partial}{\partial x_2}$ (this amounts mainly to choose ln-coordinates) while the reference singular direction is identified to the straight line $(0x_1)$.

Expanding X in the jet space at (0,0), this leads to analyze the control system:

$$ (4.6) \qquad \dot{x}_1 = -\lambda x_1 + \alpha x_2^2, \qquad \dot{x}_2 = u - u_e $$

with $u_e \in ]-1, 1[, |u| \le 1, \lambda > 0$ and $\alpha > 0$.

Properties of the system. Computing Lie brackets in the new coordinates show relevant simplification with respect to the previous formulae:

$$ X(x) = (-\lambda x_1 + \alpha x_2^2) \frac{\partial}{\partial x_1} - u_e \frac{\partial}{\partial x_2}, \quad Y(x) = \frac{\partial}{\partial x_2}, $$

$$ \begin{aligned} [Y, X](x) = -2\alpha x_2 \frac{\partial}{\partial x_1}, & \quad [[Y, X], Y](x) = -2\alpha \frac{\partial}{\partial x_1}, \end{aligned} $$

and the singular line is given by

$$ x_2 = 0. $$

Restricting to $x_2 = 0$, one has:

$$ X|_{x_2=0} = -\lambda x_1 \frac{\partial}{\partial x_1}, \quad [[Y, X], Y]|_{x_2=0}(x) = -2\alpha \frac{\partial}{\partial x_1}. $$

Hence

$$ [[Y, X], Y]|_{x_2=0}(x) = \frac{2\alpha}{\lambda x_1} X|_{x_2=0}(x) $$

for $x_1 \ne 0$.

Therefore, we obtain the following lemmas.
771 Lemma 4.4. 1. The origin (0,0) is an abnormal singular arc reduced to a point and
the arc $x_1 > 0$ is hyperbolic and the arc $x_1 < 0$ is elliptic.

772
773
774
2. The singular control along the line $x_2 = 0$ is given by $u = u_e$ and is constant and strictly admissible since $u_e \in ]-1, +1[$.

775 Lemma 4.5. The collinear set C given by $\det(X, Y) = 0$ is the parabola $x_1 = \alpha x_2^2/\lambda$.

Clock form. To analyze a 2d-time minimal problem with fixed extremities the standard
technique is to introduce the clock form $\omega = pdx$ defined outside the collinear set by:

$$
p \cdot X(x) = 1, \quad p \cdot Y(x) = 0.
$$

Computing one has $\omega = \frac{1}{-\lambda x_1 + \alpha x_2^2} dx_1$ so that

$$
d\omega = \frac{2\alpha x_2}{(-\lambda x_1 + \alpha x_2^2)^2} dx_1 \wedge dx_2.
$$

One can decompose $\mathbb{R}^2 \setminus (C \cup S)$ in four domains:

• domain A: $\dot{x}_1 < 0 \cap x_2 > 0$,

• domain B: $\dot{x}_1 > 0 \cap x_2 > 0$,

• domain C: $\dot{x}_1 < 0 \cap x_2 < 0$,

• domain D: $\dot{x}_1 > 0 \cap x_2 < 0$.

On each domain one can compare the time along arcs $\gamma_1, \gamma_2$ joining respectively $x_0$ to $x_1$
where $\gamma_1 = \sigma_+\sigma_-$, $\gamma_2 = \sigma_+\sigma_+$ using Stokes theorem. One has:

Lemma 4.6. In domain A and D for such arcs the time minimal policy is $\sigma_{-}\sigma_{+}$ while in domain B and C the optimal policy us $\sigma_{+}\sigma_{-}$.

*Proof.* Take the case of domain A, one has:

$$
\int_{\gamma_1} \omega - \int_{\gamma_2} \omega = \int_{\gamma_1 \vee (-\gamma_2)} d\omega > 0,
$$

hence the time along $\gamma_1$ is longer than the time along $\gamma_2$. The discussion is similar for the
other cases. ■

*Integrating the extremal curves.* The adjoint system takes the form with $p = (p_1, p_2)$:

$$
\dot{p}_1 = \lambda p_1, \quad \dot{p}_2 = -2\alpha x_2 p_1.
$$

Denoting for $u = \pm 1$, $\beta = u - u_e$, one gets:

Lemma 4.7. The extremal system is characterized by:

• $x_2(t) = x_2(0) + \beta t$,

$$
\bullet x_1(t) = e^{-\lambda t} \left( x_1(0) + \int_0^t e^{\lambda s} (x_2(0) + \beta s)^2 ds \right)
$$

• $p_1(t) = e^{\lambda t} p_1(0)$,

• $p_2(t) = -2\alpha p_1(0) \int_0^t e^{\lambda s} (x_2(0) + \beta s) ds$

and they belong to the polynomial exponential category.
The integrals in the expressions of $x_1(t)$ and $p_2(t)$ can be evaluated using:

$$ \int_{0}^{t} se^{\lambda s} = \frac{te^{\lambda t}}{\lambda} - \frac{1}{\lambda^2} (e^{\lambda t} - 1), \qquad \int_{0}^{t} s^2 e^{\lambda s} = \frac{t^2 e^{\lambda t}}{\lambda} - \frac{2}{\lambda} \int_{0}^{t} se^{\lambda s} ds. $$

**Lemma 4.8.** The switching function is $\Phi(t) = p_2(t)$ so that for $u = \pm 1$ and one has: $\ddot{p}_2(t) = -2\alpha p_1(\lambda x_2 + u - u_e)$ where $p_1(t)$ is of constant sign given by the sign of $p_1(0)$.

**Geometric discussion of the synthesis with a terminal manifold N of codimension one.** Next we present the discussion of the time minimal synthesis with a terminal manifold N of codimension one by gluing cases discussed in [6].

One takes N as a circle with radius d centered at 0 where the synthesis amounts to glue the hyperbolic and elliptic situation. The circle intersects the hyperbolic arc $\sigma_s$ at $x_1 = d$ and the elliptic arc $\sigma_s$ at $x_1 = -d$.

The BC-extremals curves can be parameterized by Lemma 4.7 with $(p_1(0), p_2(0)) = \pm n(0)$ where $n(0)$ is the normal to the circle, $n(0) = (x(s), y(s))$. Geodesics curves are integrated backwards, geodesics in the interior to the circle are associated to $n(0)$ and geodesics exterior to the circle are associated to $-n(0)$.

Since $Y = \frac{\partial}{\partial x_2}$, from [6] we can at once deduced the synthesis outside the circle, near the hyperbolic point $(d, 0)$ and the elliptic point $(-d, 0)$, using the curvature of N only at such point which are the images of the curve $(-1/2ks^2, s)$, $k > 0$ at the hyperbolic point and $(1/2ks^2, s)$ at the elliptic point. They are represented on Fig. 7.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** 2d-syntheses.</figcaption></figure>

The main properties are

* **hyperbolic case:** the singular arc is optimal and the optimal policy is $-1$ for $x_2 > 0$ and $+1$ for $x_2 < 0$.

* **elliptic case:** The time minimal synthesis is defined by the stratification $W = W_{-} \cup W_{+}$ of the switching locus and there exists a cut locus C terminating at $(-d, 0)$.

To complete the analysis one must glue the two syntheses along the target N and consider also geodesics interior to the circle.

Using the analysis of [6] to glue different Whitney charts, we obtain the time minimal synthesis in a neighborhood V of the origin, which is schematically represented in Fig. 8.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Schematic representation of the time minimal synthesis gluing hyperbolic and elliptic case with N being the unit circle. The cut locus was numerically computed with parameters $\lambda = \alpha = 1$ and $u_e = 1/2$.</figcaption></figure>

**4.4.2. Complexity of the singular flow in the 3d-case.** From dynamical system point of view, the complexity of the Lotka–Volterra model is related to the existence of persistent equilibrium point, which leads to complicated dynamics related to the hierarchy of dynamics associated to the hierarchy of at most $2^n$–equilibria.

Hence, in this section one analyses the same question regarding to existence of permanent equilibria for the singular dynamics associated to the time minimal control problem.

First one must extent the previous 2d–result concerning the existence of (exceptional) singular arc reduced to a point, in the context of controlled Lotka–Volterra model.

**Proposition 4.9.** *Consider the controlled Lotka–Volterra model in the n–dimensional case. Then there exists (isolated) exceptional arcs reduced to a point.*

*Proof.* Consider the pair $(X, Y)$, and choose ln-coordinates so that $Y(x) = \frac{\partial}{\partial x_n}$ with $X(x) = \sum_{i=1}^{n} X_i(x) \frac{\partial}{\partial x_i}$, the collinear set is defined by the $(n-1)$–equations $X_i(x) = 0$, $i = 1, \dots, n-1$ and let $\lambda = -u_e$ so that $X_n(x) = -\lambda$. Let $\tilde{X}(x) = \sum_{i=1}^{n-1} X_i(x) \frac{\partial}{\partial x_i} + (X_n(x) + u_e) \frac{\partial}{\partial x_n}$ and by construction there exists $x_e$ on the collinear set so that $\tilde{X}(x_e) = 0$. Denote by $\tilde{J}$ the Jacobian matrix of $\tilde{X}$ at $x = x_e$ and let $\sigma(\tilde{J})$ be its spectrum.

Denoting by $\mathrm{ad}\tilde{X} \cdot Y = [\tilde{X}, Y]$. Then at $x_e$ the matrix with columns of iterated Lie brackets

$$K = (Y, \mathrm{ad}\tilde{X} \cdot Y, \dots, \mathrm{ad}^{n-1}\tilde{X} \cdot Y)$$

coincides with the Kalman matrix

$$(b, \tilde{J}b, \dots, \tilde{J}^{n-1}b)$$

where $b$ is the constant vector $Y(x)$. The singular point is exceptional if and only if $\mathrm{rank}K \le n-1$.
At isolated point, the condition rank = n − 1 can be realized for the controlled Lotka–Volterra model (see Example 4.10 for n = 3).

Example 4.10. In the dimensionless coordinates (see Proposition 4.1), take: A = diag( λ 1 , λ 1 , λ 3 ) , X ( x ) = − diag( x + 1) A x, Y ( x ) = diag( x + 1) ( ε 1 , ε 2 , ε 3 ) T so that the persistent equi- librium is located at x$_{e}$ = 0. The columns of the Kalman matrix K defined in the proof of Proposition 4.9 are

$$
\begin{align*}
Y(0) &= (\varepsilon_1, \varepsilon_2, \varepsilon_3)^\top, \\
[X, Y](0) &= (-\varepsilon_1\lambda_1, -\varepsilon_2\lambda_1, -\varepsilon_3\lambda_1)^\top, \\
[X, [X, Y]](0) &= (\varepsilon_1\lambda_1^2, +\varepsilon_2\lambda_1^2, \varepsilon_3\lambda_1^2)^\top,
\end{align*}
$$

hence rank K < 3 and 0 is a singular exceptional point.

Note that the Jacobian matrix of the singular flow, defined for ε 1 ̸ = ε 2 , evaluated at 0 is

$$
J = \begin{pmatrix}
\frac{\varepsilon_1 \lambda_3}{\varepsilon_1 - \varepsilon_2} - \lambda_1 & \frac{\varepsilon_1 \lambda_3}{\varepsilon_2 - \varepsilon_1} & 0 \\
\frac{\varepsilon_2 \lambda_3}{\varepsilon_1 - \varepsilon_2} & \frac{\varepsilon_1 \lambda_3}{\varepsilon_2 - \varepsilon_1} - \lambda_1 + \lambda_3 & 0 \\
\frac{\varepsilon_3 \lambda_3}{\varepsilon_1 - \varepsilon_2} & -\frac{\varepsilon_3 \lambda_3}{\varepsilon_1 - \varepsilon_2} & -\lambda_3
\end{pmatrix}
$$

and its spectrum { λ 3 − λ 1 , − λ 3 , − λ 1 } is resonant.

3d–case. In the 3d–case, the singular trajectories are solutions of the vector field: ˙ x = X$_{s}$ ( x ) = X ( x ) − u$_{s}$Y ( x ) , where the singular control feedback is u$_{s}$ = − D $^{′}$( x ) /D ( x ) with

$$
D(x) = \det(Y(x), [Y, X](x), [[Y, X], Y](x)), \quad D'(x) = \det(Y(x), [Y, X](x), [[Y, X], X](x)).
$$

Moreover exceptional trajectories are located on the exceptional locus D ′′ ( x ) = 0 with

$$
D''(x) = \det(Y(x), [Y, X](x), X(x)).
$$

Computing in the original coordinates leads to complicated expressions:

$$
\begin{align*}
D(x)/x_1 x_2 x_3 &= (\varepsilon_1^2 x_1 a_{21} + \varepsilon_1 (\varepsilon_2 (x_2 a_{22} - x_1 a_{11}) + \varepsilon_3 x_3 a_{23}) - \varepsilon_2 (\varepsilon_2 x_2 a_{12} + \varepsilon_3 x_3 a_{13})) \\
&\quad + (\varepsilon_1^2 x_1 a_{31} + \varepsilon_2^2 x_2 a_{32} + \varepsilon_3^2 x_3 a_{33}) + (\varepsilon_1^2 x_1 a_{11} + \varepsilon_2^2 x_2 a_{12} + \varepsilon_3^2 x_3 a_{13}) (\varepsilon_2^2 x_2 a_{32} + \varepsilon_3 \varepsilon_2 (x_3 a_{33} - x_2 a_{22}) \\
&\quad - \varepsilon_3^2 x_3 a_{23} + \varepsilon_1 x_1 (\varepsilon_2 a_{31} - \varepsilon_3 a_{21})) - (\varepsilon_1^2 x_1 a_{21} + \varepsilon_2^2 x_2 a_{22} + \varepsilon_3^2 x_3 a_{23}) \\
&\quad + (\varepsilon_1^2 x_1 a_{31} + \varepsilon_1 (\varepsilon_2 x_2 a_{32} + \varepsilon_3 (x_3 a_{33} - x_1 a_{11})) - \varepsilon_3 (\varepsilon_2 x_2 a_{12} + \varepsilon_3 x_3 a_{13})),
\end{align*}
$$

$$
\begin{align*}
D'(x)/x_1 x_2 x_3 &= (-\varepsilon_1^2 x_1 a_{21} + \varepsilon_1 (\varepsilon_2 (x_1 a_{11} - x_2 a_{22}) - \varepsilon_3 x_3 a_{23}) + \varepsilon_2 (\varepsilon_2 x_2 a_{12} + \varepsilon_3 x_3 a_{13})) \\
&\quad (\varepsilon_2 x_2 (x_1 a_{12} a_{31} - a_{32} (x_1 a_{21} + x_3 (a_{23} - a_{33}) + r_2)) - \varepsilon_1 x_1 (r_1 a_{31} + x_3 (a_{13} - a_{33}) a_{31} \\
&\quad + x_2 (a_{12} a_{31} - a_{21} a_{32})) + \varepsilon_3 x_3 (-r_3 a_{33} + x_1 a_{31} (a_{13} - a_{33}) + x_2 a_{32} (a_{23} - a_{33}))) \\
&\quad + (\varepsilon_2^2 (-x_2) a_{32} + \varepsilon_3 \varepsilon_2 (x_2 a_{22} - x_3 a_{33}) + \varepsilon_3^2 x_3 a_{23} + \varepsilon_1 x_1 (\varepsilon_3 a_{21} - \varepsilon_2 a_{31})) \\
&\quad (-\varepsilon_1 x_1 (r_1 a_{11} + x_2 a_{12} (a_{11} - a_{21}) + x_3 a_{13} (a_{11} - a_{31})) + \varepsilon_2 x_2 (x_3 a_{13} a_{32} - a_{12} (x_1 (a_{21} - a_{11} \\
&\quad + x_3 a_{23} + r_2)) - \varepsilon_3 x_3 (a_{13} (x_1 (a_{31} - a_{11}) + x_2 a_{32} + r_3) - x_2 a_{12} a_{23}))
\end{align*}
$$
$$ 
\begin{aligned}
& - (\varepsilon_1^2 (-x_1) a_{31} + \varepsilon_1 (\varepsilon_3 (x_1 a_{11} - x_3 a_{33}) - \varepsilon_2 x_2 a_{32}) + \varepsilon_3 (\varepsilon_2 x_2 a_{12} + \varepsilon_3 x_3 a_{13})) \\
& (\varepsilon_1 x_1 (x_3 a_{23} a_{31} - a_{21} (x_3 a_{13} + x_2 (a_{12} - a_{22}) + r_1)) + \varepsilon_2 x_2 (-r_2 a_{22} + x_1 a_{21} (a_{12} - a_{22}) \\
& + x_3 a_{23} (a_{32} - a_{22}))) + \varepsilon_3 x_3 (x_1 a_{13} a_{21} - a_{23} (x_1 a_{31} + x_2 (a_{32} - a_{22}) + r_3))),
\end{aligned}
$$

$$ 
\begin{aligned}
D''(x)/x_1 x_2 x_3 = & (-\varepsilon_1^2 x_1 a_{21} + \varepsilon_1 (\varepsilon_2 (x_1 a_{11} - x_2 a_{22}) - \varepsilon_3 x_3 a_{23}) + \varepsilon_2 (\varepsilon_2 x_2 a_{12} + \varepsilon_3 x_3 a_{13})) \\
& (x_1 a_{31} + x_2 a_{32} + x_3 a_{33} + r_3) + (-\varepsilon_2^2 x_2 a_{32} + \varepsilon_3 \varepsilon_2 (x_2 a_{22} - x_3 a_{33}) + \varepsilon_3^2 x_3 a_{23} + \varepsilon_1 x_1 (\varepsilon_3 a_{21} \\
& - \varepsilon_2 a_{31})) (x_1 a_{11} + x_2 a_{12} + x_3 a_{13} + r_1) + (\varepsilon_1^2 x_1 a_{31} + \varepsilon_1 (\varepsilon_2 x_2 a_{32} + \varepsilon_3 (x_3 a_{33} - x_1 a_{11})) \\
& - \varepsilon_3 (\varepsilon_2 x_2 a_{12} + \varepsilon_3 x_3 a_{13})) (x_1 a_{21} + x_2 a_{22} + x_3 a_{23} + r_2).
\end{aligned}
$$

**Proposition 4.11.** If $D(x) \neq 0$, the equilibria of the singular dynamics $\dot{x} = X_s(x)$ are exceptional trajectories reduced to a point. At such a point $x_e$, the spectrum of $J = \frac{\partial X_s}{\partial x}(x_e)$ is a feedback invariant. Moreover the dynamics is foliated by the invariant set $D''(x) = 0$ and $D(x)D''(x) > 0$ or $< 0$ associated respectively to hyperbolic and elliptic arcs. The singular feedback $u_s$ acts as a geometric pole placement of the dynamics on the collinear set.

*Proof.* The proof is clear following the construction detailed in the proof of Proposition 4.9.

**4.4.3. The 4d-case.** This gives the road to the n-dimensional case.

The singular exceptional control can be expressed as a feedback using the relation

$$ 
\begin{aligned}
H_X(z) = H_Y(z) &= \{H_X, H_Y\}(z) = 0 \\
\{\{H_X, H_Y\}, H_X\}(z) + u_{se} \{\{H_X, H_Y\}, H_Y\}(z) &= 0
\end{aligned}
$$

and this leads to

$$ u_{se}(x) = -\frac{D'(x)}{D(x)}, $$

where

$$ 
\begin{aligned}
D(x) &= \det(X(x), Y(x), [Y, X](x), [[Y, X], Y](x)), \\
D'(x) &= \det(X(x), Y(x), [Y, X](x), [[Y, X], X](x)).
\end{aligned}
$$

Similarly to the 3d-case, the singular exceptional dynamics: $\dot{x} = X(x) + u_{se}Y(x)$ can be used to generate feedback invariants.

Remaining singular dynamics are parameterized by a dynamic feedback $u(x, \lambda)$ depending upon a one dimensional coefficient.

We refer to [9] for a general approach to compute feedback invariants for the dynamics.

**5. Conclusion.** In this article we have presented the general techniques from geometric control to analyze in the permanent case, the optimal control problem related to vermin reduction in a complex microbiote modelled by the Lotka–Volterra equations.

Our analysis is based on a series of articles classifying the time minimal syntheses for a single-input affine system with terminal manifold of codimension one developed for chemical networks [6, 7, 21]. Using the concepts of Whitney chart and unfolding the explicit computations in a neighbourhood of the terminal manifold can be reduced to problems in dimension 2
<table><caption>**Table 1**</caption><tbody><tr><td>Barnesiella (Bar.)</td><td></td><td>0.3680</td><td>Akkermansia (Akk.)</td><td></td><td>0.2297</td></tr><tr><td>undefined genus of Lachnospiraceae (Und. Lac.)</td><td></td><td>0.3102</td><td>Coprobacillus (Cop.)</td><td></td><td>0.8300</td></tr><tr><td>undefined genus of unclassified Mollicutes (Und. Mol.)</td><td></td><td>0.4706</td><td>Clostridium difficile (C. diff.)</td><td></td><td>0.3918</td></tr><tr><td>unclassified Lachnospiraceae (Uncl. La.)</td><td></td><td>0.3561</td><td>Enterococcus (Ent.)</td><td></td><td>0.2907</td></tr><tr><td>Blautia (Bla.)</td><td></td><td>0.7089</td><td>undefined genus of Enterobacteriaceae (Und. En.)</td><td></td><td>0.3236</td></tr><tr><td>Other</td><td></td><td>0.5400</td><td></td><td></td><td></td></tr><tr><td></td><td>Bar.</td><td>Und. Lac.</td><td>Uncl. Lac.</td><td>Other</td><td>Bla.</td><td>Und. Mol.</td><td>Akk.</td><td>Cop.</td><td>Und. En.</td><td>Ent.</td><td>C. diff.</td></tr><tr><td>Bar.</td><td>-0.205</td><td>0.098</td><td>0.167</td><td>-0.164</td><td>-0.143</td><td>0.019</td><td>-0.515</td><td>-0.391</td><td>-0.268</td><td>0.008</td><td>0.346</td></tr><tr><td>Und. Lac.</td><td>0.062</td><td>-0.104</td><td>-0.043</td><td>-0.154</td><td>-0.187</td><td>0.027</td><td>-0.459</td><td>-0.413</td><td>-0.196</td><td>0.022</td><td>0.301</td></tr><tr><td>Uncl. Lac.</td><td>0.143</td><td>-0.192</td><td>-0.101</td><td>-0.139</td><td>-0.165</td><td>0.013</td><td>-0.504</td><td>-0.772</td><td>-0.206</td><td>-0.006</td><td>0.292</td></tr><tr><td>Other</td><td>0.224</td><td>0.138</td><td>0.000</td><td>-0.831</td><td>-0.223</td><td>0.220</td><td>-0.205</td><td>-1.009</td><td>-0.400</td><td>-0.039</td><td>0.666</td></tr><tr><td>Bla.</td><td>-0.180</td><td>-0.051</td><td>0.000</td><td>-0.054</td><td>-0.708</td><td>0.016</td><td>-0.507</td><td>0.553</td><td>0.106</td><td>0.224</td><td>0.157</td></tr><tr><td>Und. Mol.</td><td>-0.111</td><td>-0.037</td><td>-0.042</td><td>0.041</td><td>0.261</td><td>-0.422</td><td>-0.185</td><td>-0.432</td><td>-0.264</td><td>-0.061</td><td>0.164</td></tr><tr><td>Akk.</td><td>-0.126</td><td>-0.185</td><td>-0.122</td><td>0.380</td><td>0.400</td><td>-0.160</td><td>-1.212</td><td>1.389</td><td>-0.096</td><td>0.191</td><td>-0.379</td></tr><tr><td>Cop.</td><td>-0.071</td><td>0.000</td><td>0.080</td><td>-0.454</td><td>-0.503</td><td>0.169</td><td>-0.562</td><td>-4.350</td><td>-0.207</td><td>-0.223</td><td>0.443</td></tr><tr><td>Und. Ent.</td><td>-0.374</td><td>0.278</td><td>0.248</td><td>-0.168</td><td>0.084</td><td>0.033</td><td>-0.232</td><td>-0.395</td><td>-0.384</td><td>-0.038</td><td>0.314</td></tr><tr><td>Ent.</td><td>-0.042</td><td>-0.013</td><td>0.024</td><td>-0.117</td><td>-0.328</td><td>0.020</td><td>0.054</td><td>-2.096</td><td>0.023</td><td>-0.192</td><td>0.111</td></tr><tr><td>C. diff.</td><td>-0.037</td><td>-0.033</td><td>-0.049</td><td>-0.090</td><td>-0.102</td><td>0.032</td><td>-0.181</td><td>-0.303</td><td>-0.007</td><td>0.014</td><td>-0.055</td></tr></tbody></table>

(top) Growth rates $a_{ij}$ of each microbial population $i$ of the CDI model. (bottom) Interactions between pairwise microbial populations of the CDI model. Both tables are excerpted from [27].

or 3 and a dictionary of the time minimal syntheses is described in [6, 21], up to codimension 2 cases. Global syntheses can be described by gluing different Whitney charts.

It can be applied to the controlled Lotka-Volterra model to analyze either the problem of reducing the infection using an antibiotic agent or to reinforce the body prior to infection using a probiotic agent. A case study is given to construct a global optimal synthesis by gluing distinct Whitney charts.

Our article shows the parallel between the analysis of the free dynamics in the frame of dynamical systems and the dynamics of the Hamiltonian dynamics deduced from the Maximum Principle which parameterizes the extremals candidates as minimizers.

The complexity of the Hamiltonian dynamics for nonlinear control systems is related to the existence and their complexity of the singular extremals dynamics. Our contribution is to present preliminary analysis of this complexity, in the frame of the controlled Lotka-Volterra model. It is shown to be related to the collinear set on which are located the equilibrium of the free dynamics, where no treatment is applied, and the forced equilibrium associated to maximal treatment.

An additional step in our analysis will be to analyze the problem in the sampled-data control frame in relation with the permanent case, see [8] for such an analysis.

REFERENCES

[1] A.A. AGRACHEV, Y.L. SACHKOV, *Control theory from the geometric viewpoint*. Encyclopaedia of Mathematical Sciences, 87. Control Theory and Optimization, II. Springer-Verlag, Berlin, 2004, 412 pages.

[2] F. ANCONA, A. BRESSAN, *Nearly time optimal stabilizing patchy feedbacks*. Ann. Inst. H. Poincaré Anal. Non Linéaire **24** no. 2 (2007), pp. 279–310.
940

941

[3] S. BAIGENT, *Geometry of carrying simplices of 3-species competitive Lotka-Volterra systems*. Nonlinearity, **26** no. 4 (2013), pp. 1001-1029.

942

943

[4] B. BONNARD, M. CHYBA, *The role of singular trajectories in control theory*. Springer Verlag, New York, 2003, 357 pages.

944

945

[5] B. BONNARD, I. KUPKA, *Théorie des singularités de l'application entrée/sortie et optimalité des trajets toires singulières dans le problème du temps minimal*. Forum Math. **5** no. 2 (1993), pp. 111-159.

946

947

[6] B. BONNARD, G. LAUNAY, M. PELLETIER, *Classification générique de synthèses temps minimales avec cible de codimension un et applications*. Annales de l'I.H.P. Analyse non linéaire **14** no.1 (1997), pp. 55-102.

948

949

950

951

[7] B. BONNARD, J. ROUOT, *Towards Geometric Time Minimal Control without Legendre Condition and with Multiple Singular Extremals for Chemical Networks*. Advances in Nonlinear Biological Systems, Modeling and Optimal Control, AIMS on applied Maths **11** (2021), pp. 1-34.

952

953

[8] B. BONNARD, J. ROUOT, C. SILVA, *Geometric Optimal Control of the Generalized Lotka-Volterra Model of the Intestinal Microbiome*. Preprint 2023: hal-03861565.

954

955

[9] B. BONNARD, J. ROUOT, *Feedback Classification and Optimal Control with Applications to the Controlled Lotka-Volterra Model*. Preprint 2023: hal-03917363.

956

957

[10] V.G. BOLTYANSKII, *Sufficient conditions for optimality and the justification of the dynamic programming method*. SIAM J. Control **4** (1966), pp. 326-361.

958

959

[11] U. BOSCAIN, B. PICCOLI, *Optimal syntheses for control systems on 2-D manifolds*. Springer Science & Business Media **43** 2003, 275 pages.

960

961

[12] P. BRUNOVSKÝ, *Existence of regular synthesis for general control problems*. J. Differential Equations **38** no. 3 (1980), pp. 317-343.

962

963

[13] W.L. CHOW, *Über Systeme von linearen partiellen Differentialgleichungen erster Ordnung*. Mathematische Annalen, **117** (1939) pp. 98-105.

964

965

[14] I. EKELAND, *Discontinuités de champs hamiltoniens et existence de solutions optimales en calcul des variations*. Inst. Hautes Études Sci. Publ. Math. **47** (1977), pp. 5-32.

966

967

[15] M. FEINBERG, *Foundations of chemical reaction network theory*. Applied Mathematical Sciences, 202. Springer, Cham, 2019, 454 pages.

968

969

[16] M.W. HIRSCH, *Systems of differential equations which are competitive or cooperative: III. Competing species*. Nonlinearity, **1**, (1988), pp.51-71.

970

971

972

[17] E.W. JONES, P. S. CLARCKE, J. M. CARSLON, *Navigation of outcome in a generalized Lotka-Volterra model of the microbiome*. Advances in Nonlinear Biological Systems, Modeling and Optimal Control, AIMS on applied Maths **11** (2021), pp. 97-117.

973

974

[18] A.J. KRENER, *The high order maximal principle and its application to singular extremals*. SIAM J. Control Optim. **15** no.2 (1977), pp. 256-293.

975

976

[19] I. KUPKA, *Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case*. Trans. Amer. Math. Soc. **299** no.1 (1987), pp. 225-243.

977

978

[20] I. KUPKA, *The ubiquity of Fuller's phenomenon*. Nonlinear controllability and optimal control, 313-350, Monogr. Textbooks Pure Appl. Math., 133, Dekker, New York, 1990.

979

980

981

[21] G. LAUNAY, M. PELLETIER, *The generic local structure of time-optimal synthesis with a target of codimension one in dimension greater than two*. Journal of Dynamical and Control Systems **3**, no. 165 (1997).

982

983

[22] E.B. LEE, L. MARKUS, *Foundations of optimal control theory*. John Wiley & Sons, Inc., New York-London-Sydney 1967, 576 pages.

984

985

[23] A.J. LOTKA, *Elements of mathematical biology*. Dover Publications, Inc., New York, N.Y., 1958, 465 pages.

986

987

[24] L.S. PONTRYAGIN, V.G. BOLTYANSKII, R.V. GAMKRELIDZE, E.F. MISHCHENKO, *The mathematical theory of optimal processes*. Oxford, Pergamon Press, 1964, 362 pages.

988

989

990

[25] H. SCHÄTTLER, U. LEDZEWICZ, *Optimal control for mathematical models of cancer therapies. An application of geometric methods*. Interdisciplinary Applied Mathematics, 42. Springer, New York, 2015, 496 pages.

991

992

[26] S. SMALE, *On the differential equations of species in competition*. Journal of Mathematical Biology, **3** (1976), pp. 5-7.

993

[27] R.R. STEIN, V. BUCCI, N.C. TOUSSAINT, C.G. BUFFIE, G. RÄTSCH, E.G. PAMER, ET AL., *Ecological*
994

modelling from time-series inference: insight into dynamics and stability of intestinal microbiota.

995

PLos Comp. Biology, **9** no. 12 (2013).

996

[28] V. VOLTERRA, *Leçons sur la théorie mathématique de la lutte pour la vie.* Les Grands Classiques

997

Gauthier-Villars. Éditions Jacques Gabay, Sceaux, 1990, 215 pages.

998

[29] E.C. ZEEMAN, M.L. ZEEMAN, *From local to global behavior in competitive Lotka Volterra systems*. Trans.

999

Amer. Math. Soc., 355 (2003), pp. 713–734.