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

HAL
open science

# Feedback Classification and Optimal Control with Applications to the Controlled Lotka-Volterra Model

Bernard Bonnard, Jérémy Rouot

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

To cite this version:

Bernard Bonnard, Jérémy Rouot. Feedback Classification and Optimal Control with Applications to the Controlled Lotka-Volterra Model. 2024. hal-03917363v2

HAL Id: hal-03917363

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

Preprint submitted on 24 May 2024

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

Distributed under a Creative Commons Attribution 4.0 International License
# Feedback Classification and Optimal Control with Applications to the Controlled Lotka-Volterra Model

Bernard Bonnard, Jérémy Rouot

► To cite this version:

Bernard Bonnard, Jérémy Rouot. Feedback Classification and Optimal Control with Applications to the Controlled Lotka-Volterra Model. Optimization, 2024, 10.1080/02331934.2024.2392209. hal-03917363v4

HAL Id: hal-03917363

https://inria.hal.science/hal-03917363v4

Submitted on 4 May 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
ARTICLE TEMPLATE

Feedback Classification and Optimal Control with Applications to the
Controlled Lotka–Volterra Model

Bernard Bonnard<sup>a</sup> and Jérémy Rouot<sup>b</sup>

<sup>a</sup>Institut Mathématique de Bourgogne and Inria Sophia Antipolis, 9 rue Alain Savary, 21 000 Dijon

<sup>b</sup>Univ Brest, UMR CNRS 6205, Laboratoire de Mathématiques de Bretagne Atlantique, Brest

7

ARTICLE HISTORY

9

Compiled May 23, 2024

10

*Dedicated to Roger W. Brockett and Ivan Kupka.*

ABSTRACT

11
12
13
14
15
16
17

Let $M$ be a $\sigma$-compact $C^\infty$ manifold of dimension $n \ge 2$ and consider a single-input con-
trol system: $\dot{x}(t) = X(x(t)) + u(t) Y(x(t))$, where $X, Y$ are $C^\infty$ vector fields on $M$. We prove
that there exist an open set of pairs $(X, Y)$ for the $C^\infty$-Whitney topology such that they
admit singular abnormal rays so that the spectrum of the projective singular Hamiltonian
dynamics is feedback invariant. It is applied to controlled Lotka-Volterra dynamics where
such rays are related to shifted equilibria of the free dynamics.

KEYWORDS

18
19
20

Feedback classification; Nonlinear systems; Lotka–Volterra model; Optimal control; Direct numerical methods

1. Introduction

Consider a single-input affine control system on a $\sigma$-compact $C^\infty$ manifold of dimension $n \ge 2$ defined by

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

where X, Y are C∞ vector fields on M and the set of admissible controls U is the set of
bounded measurable mappings u : [0, T(u)] → ℝ, T(u) > 0.

Let (X, Y) and (X', Y') be two such pairs. They are called *feedback equivalent* if there
exist a $C^\infty$ diffeomorphism $\varphi$ on M and a feedback $u = \alpha(x) + \beta(x)u'$, where $\alpha, \beta$ are $C^\infty$
mappings and $\beta$ invertible so that

(i) $X' = \varphi * X + (\varphi * Y) \cdot \alpha,$

(ii) $Y' = (\varphi * Y) \cdot \beta,$
# Feedback Classification and Optimal Control with Applications to the Controlled Lotka–Volterra Model

Bernard Bonnard<sup>a</sup> and Jérémy Rouot<sup>b</sup>

## ARTICLE HISTORY

Compiled May 4, 2025

*Dedicated to Roger W. Brockett and Ivan Kupka.*

## ABSTRACT

Let $M$ be a $\sigma$-compact $C^\infty$ manifold of dimension $n \ge 2$ and consider a single-input control system: $\dot{x}(t) = X(x(t)) + u(t) Y(x(t))$, where $X, Y$ are $C^\infty$ vector fields on $M$. We prove that there exist an open set of pairs $(X, Y)$ for the $C^\infty$-Whitney topology such that they admit singular abnormal rays so that the spectrum of the projective singular Hamiltonian dynamics is feedback invariant. This result is applied to the controlled Lotka–Volterra dynamics, where such rays are related to shifted equilibria of the free dynamics.

## KEYWORDS

Feedback classification; Nonlinear systems; Lotka–Volterra model; Optimal control; Direct numerical methods

## 1. Introduction

Consider a single-input affine control system on a $\sigma$-compact $C^\infty$ manifold $M$ of dimension $n \ge 2$ defined by

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

where $X, Y$ are $C^\infty$ vector fields on $M$ and the set of admissible controls $\mathcal{U}$ is the set of bounded measurable mappings $u: [0, T(u)] \to \mathbb{R}, T(u) > 0$.

Let $(X, Y)$ and $(X', Y')$ be two such pairs. They are called *feedback equivalent* if there exist a $C^\infty$ diffeomorphism $\varphi$ on $M$ and a feedback $u = \alpha(x) + \beta(x)u'$, where $\alpha, \beta$ are $C^\infty$ mappings and $\beta$ invertible so that

$$ (i) \quad X' = \varphi * X + (\varphi * Y) \cdot \alpha, $$

$$ (ii) \quad Y' = (\varphi * Y) \cdot \beta, $$

<footnote><sup>a</sup>Institut Mathématique de Bourgogne and Inria Sophia Antipolis, 9 rue Alain Savary, 21 000 Dijon</footnote>

<footnote><sup>b</sup>Univ Brest, UMR CNRS 6205, Laboratoire de Mathématiques de Bretagne Atlantique, Brest</footnote>
where if $Z$ is a $C^\infty$ vector field $\varphi * Z$ is the image of $Z$ by $\varphi$ given in local coordinates by

$$ \varphi * Z = \left( \frac{\partial \varphi}{\partial x} \right)^{-1} (Z \circ \varphi). $$

This action defines a group structure $G_f$ on the set of triplets $\{g = (\varphi, \alpha, \beta)\}$ called the *feed-back group*.

It is well known from geometric linear system theory [23] that, restricting to linear autonomous controllable systems: $\dot{x}(t) = Ax(t) + u(t)b$, $A$ being a constant matrix and $b$ a constant vector, every such pairs $(A, b)$ and $(A', b')$ are feedback equivalent, restricting to linear diffeomorphisms and feedbacks: $u = \alpha \cdot x + \beta u'$, $\beta \neq 0$ constant.

The feedback equivalence for control affine systems was studied in the earliest reference [4] (see also [1]) in relation with the time-minimal geodesics corresponding to the so-called singular trajectories. More precisely consider the pair $(X, Y)$ and denote by $x(\cdot, x_0, u)$ the response to $u(\cdot) \in \mathcal{U}$ starting at time $t=0$ from $x_0$. Fixing $t_f$, the fixed extremity mapping at time $t_f$ is the map $E^{x_0, t_f} : \mathcal{U} \ni u \mapsto x(t_f, x_0, u) \in M$. If $t_f$ is free, the extremity mapping is the map $E^{x_0} : \mathcal{U} \ni u \mapsto x(\cdot, x_0, u)$. If we endow the set of controls with the $L^\infty$-norm topology both maps are Fréchet differentiable and they have singularities, that is the image of $L^\infty$ by the derivative is not of full rank $n$. Such a pair $(x(\cdot), u(\cdot))$ is called singular for the fixed time case and if the final time is not fixed such singular trajectories are called exceptional or abnormal.

Singular trajectories can be parametrized thanks to the Maximum Principle [16] as projection of singular extremal pairs $(z(\cdot), u(\cdot))$, $z = (x, p)$ where $p$ is non vanishing adjoint vector solutions of the Hamiltonian dynamics:

$$ \begin{aligned} \frac{dz}{dt}(t) &= H_X(z(t)) + u(t) H_Y(z(t)), \\ H_Y(z(t)) &= 0, \end{aligned} \qquad (2) $$

where $H_X = p \cdot X(x)$ and $H_Y = p \cdot Y(x)$ are respectively the Hamiltonian lifts of $X, Y, H_X + uH_Y$ is the pseudo Hamiltonian and additionally in the abnormal case, we have:

$$ H_X(z(t)) = 0. $$

Given a pair $(X, Y)$ the collinearity set is the set $\mathcal{C}$ of points $x \in M$ such that $X(x)$ and $Y(x)$ are linearly dependent. Clearly it is a feedback invariant. Taking a point $x_e$ in $\mathcal{C}$ there exist a constant $u_e$ so that $X(x_e) + u_e Y(x_e) = 0$ and the pair $(x_e, u_e)$ is called a forced equilibrium.

The contribution of this article is to construct if $n \ge 2$ a set of pairs $(X, Y)$ of codimension $n$ in the jets space such that there exist abnormal trajectories reduced to isolated point $x_0$, which can be lifted into lines $\ell$ in the projective space $PT_{x_0}^*M$ called abnormal rays and the whole spectrum of the projectivized linearized dynamics defined by (2) is feedback invariant.

This result can be illustrated by the following 2d-case, $x = (x_1, x_2)$,

$$ \begin{aligned} \dot{x}_1 &= \lambda x_1 + x_2^2 \\ \dot{x}_2 &= u \end{aligned} \qquad (3) $$

in which the singular line is the axis $x_2 = 0$ on which the singular dynamics is given by $\dot{x}_1 = \lambda x_1$ and the abnormal point is the equilibrium point 0. Denoting $p = (p_1, p_2)$ the
where if Z is a $C^\infty$ vector field $\varphi * Z$ is the image of Z by $\varphi$ given in local coordinates by

$$ \varphi * Z = \left( \frac{\partial \varphi}{\partial x} \right)^{-1} (Z \circ \varphi). $$

This action defines a group structure $G_f$ on the set of triplets $\{g = (\varphi, \alpha, \beta)\}$ called the *feed-back group*.

It is well known from geometric linear system theory [23] that, restricting to linear autonomous controllable systems: $\dot{x}(t) = Ax(t) + u(t)b$, A being a constant matrix and b a constant vector, every such pairs $(A, b)$ and $(A', b')$ are feedback equivalent, restricting to linear diffeomorphisms and feedbacks: $u = \alpha \cdot x + \beta u'$, $\beta \neq 0$ constant.

The feedback equivalence for control affine systems was studied in the earliest reference [4] (see also [1]) in relation with the time-minimal geodesics corresponding to the so-called singular trajectories. More precisely consider the pair $(X, Y)$ and denote by $x(\cdot, x_0, u)$ the response to $u(\cdot) \in \mathcal{U}$ (measurable and bounded controls) starting at time $t=0$ from $x_0$. Fixing $t_f$, the fixed extremity mapping at time $t_f$ is the map $E^{x_0, t_f} : \mathcal{U} \ni u \mapsto x(t_f, x_0, u) \in M$. If $t_f$ is free, the extremity mapping is the map $E^{x_0} : \mathcal{U} \ni u \mapsto x(\cdot, x_0, u)$. If we endow the set of controls with the $L^\infty$-norm topology both maps are Fréchet differentiable and they have singularities, that is the image of $L^\infty$ by the derivative is not of full rank $n$. Such a pair $(x(\cdot), u(\cdot))$ is called singular for the fixed time case and if the final time is not fixed such singular trajectories are called exceptional or abnormal.

Singular trajectories can be parametrized thanks to the Maximum Principle [16] as projection of singular extremal pairs $(z(\cdot), u(\cdot))$, $z = (x, p)$ where $p$ is non vanishing adjoint vector solutions of the Hamiltonian dynamics:

$$ \begin{aligned} \frac{dz}{dt}(t) &= H_X(z(t)) + u(t) H_Y(z(t)), \\ H_Y(z(t)) &= 0, \end{aligned} \qquad (2) $$

where $H_X = p \cdot X(x)$ and $H_Y = p \cdot Y(x)$ are respectively the Hamiltonian lifts of $X, Y, H_X + uH_Y$ is the pseudo Hamiltonian and, additionally, in the abnormal case, we have:

$$ H_X(z(t)) = 0. $$

Given a pair $(X, Y)$ the collinearity set is the set $\mathcal{C}$ of points $x \in M$ such that $X(x)$ and $Y(x)$ are linearly dependent. Clearly it is a feedback invariant. Taking a point $x_e$ in $\mathcal{C}$ there exist a constant $u_e$ so that $X(x_e) + u_e Y(x_e) = 0$ and the pair $(x_e, u_e)$ is called a forced equilibrium.

The contribution of this article is to construct, if $n \ge 2$, a set of pairs $(X, Y)$ of codimension $n$ in the jets space such that there exist abnormal trajectories reduced to isolated point $x_0$, which can be lifted into lines $\ell$ in the projective space $PT_{x_0}^*M$ called abnormal rays and the whole spectrum of the projectivized linearized dynamics defined by (2) is feedback invariant.

This result can be illustrated by the following 2d-case, $x = (x_1, x_2)$,

$$ \begin{aligned} \dot{x}_1 &= \lambda x_1 + x_2^2 \\ \dot{x}_2 &= u \end{aligned} \qquad (3) $$

in which the singular line is the axis $x_2 = 0$ on which the singular dynamics is given by $\dot{x}_1 = \lambda x_1$ and the abnormal point is the equilibrium point 0. Denoting $p = (p_1, p_2)$ the
adjoint vector one has $p_2 = 0$ since $H_Y = 0$ and the adjoint vector is given by $\dot{p}_1 = -\lambda p_1$.
Hence the ray at $T_0^*\mathbb{R}^2$ is $\ell: p_1(t) = e^{-\lambda t}p_1(0)$. The singular Hamiltonian dynamics reads

$$ \dot{x}_1 = \lambda x_1, \quad \dot{p}_1 = -\lambda p_1. $$

Clearly $\lambda$ is feedback invariant and the projectivized dynamics is given by the dynamics of $v = x_1/p_1$ with singular point at $v = 0$ with eigenvalue $2\lambda$ (the projectivized invariant).

Note that in the $n$-dimensional case, for each ray one gets $2n-3$ eigenvalues to classify the systems.

This result is applied to analyze vermin reduction for controlled Lotka–Volterra dynamics in the cone $M = \mathbb{R}_{\ge 0}^n$ where the equations read

$$ \frac{dx}{dt}(t) = (\operatorname{diag} x(t)) (Ax(t) + r + u(t)\mathcal{E}) \quad (4) $$

with $x := (x_1, \dots, x_n)^\top$ is the vector of interacting species, $\operatorname{diag} x$ denotes in short the diagonal matrix with diagonal coefficients $x_i$, $A = (a_{ij})$ is the matrix of interaction coefficients and $r = (r_1, \dots, r_n)^\top$ is the individual growth vector without interaction. The constant vector $\mathcal{E} = (\epsilon_1, \dots, \epsilon_n)^\top$ describes the effect of probiotic or antibiotic agents to reduce the population of the infecting agent, which can be taken as the $x_1$-population. The optimal control problem is a *Mayer problem*: $\min_{u(.)} x_1(t_f)$, which can be formulated in the dual form as: $\min_{u(.)} t_f$, $x_1(t_f) = d$, where $d$ is a desired amount of population $x_1$ at final time $t_f$, taken as a parameter. More generally we shall consider any Mayer problem: $\min_{u(.)} \phi(x(t_f))$ with $\phi$ is a smooth mapping from $M$ to $\mathbb{R}$.

The free motion was studied by the historical contributors Lotka and Volterra, generalizing the prey–predator model for populations dynamics [21]. The interest of the model is to get easily computable equilibria using linear analysis only. In the regular case where $A$ is invertible, the *interior equilibrium* is $x_e := -A^{-1}r$ and the dynamics is related in particular to interaction of $x_e$ with the trivial boundary equilibrium 0. A large amount of literature was devoted to analyze the dynamics in the 2d or 3d case starting from the seminal complete analysis in the prey–predator model. In particular in the case of competitive systems, the achievement in the 3d–case was the explicit analysis of the limit set of the origin to determine the boundary of the basin of repulsion of this point. Related references are [11, 24, 19] and in particular the May–Leonard model [15] will serve as case study.

An important very recent application of controlled Lotka–Volterra model comes from the control of complex microbiota, see the work by Jones et al [12] based on the Stein et al study [20] to model the C. difficile infection of the intestinal microbiote. It contains 11 interaction species and can admit up to $2^{11}$ equilibria. Hence the analysis is a huge task. More modestly we shall concentrate our computation to specific 2d and 3d models to illustrate our techniques. Besides the biological applications, the controlled Lotka–Volterra model allows algebraic computations of the forced equilibria, which are by construction shifted equilibria of the free equilibria and hence in fine one can compute in this case, the spectrum of corresponding ray abnormal solutions. We show that this information is relevant in an applicative context to analyze a 3d May and Leonard model corresponding to a certain microbiote with three species using antibiotic and probiotic. Recent contribution [7] presents a protocol based on a MPC tracking algorithm of a trajectory computed with a direct method. The first consists in reaching the invariant manifold – carrying simplex – in minimum time using antibiotic only. The second phase of the treatment focuses on the asymptotic dynamics in the 2d carrying simplex, where the aim is to reach an healthy region using probiotics only. This article presents an alternative MPC algorithm for the
adjoint vector one has $p_2 = 0$, since $H_Y = 0$, and the adjoint vector is given by $\dot{p}_1 = -\lambda p_1$. Hence the ray at $T_0^*\mathbb{R}^2$ is $\ell: p_1(t) = e^{-\lambda t}p_1(0)$. The singular Hamiltonian dynamics reads

$$ \dot{x}_1 = \lambda x_1, \quad \dot{p}_1 = -\lambda p_1. $$

Clearly $\lambda$ is feedback invariant and the projectivized dynamics is given by the dynamics of $v = x_1/p_1$ with singular point at $v = 0$ with eigenvalue $2\lambda$ (the projectivized invariant).

Note that in the $n$-dimensional case, for each ray one gets $2n-3$ eigenvalues to classify the systems.

This result is applied to analyze vermin reduction for controlled Lotka–Volterra dynamics in the cone $M = \mathbb{R}_{\ge 0}^n$ where the equations read

$$ \frac{dx}{dt}(t) = (\operatorname{diag} x(t)) (Ax(t) + r + u(t)\mathcal{E}) \quad (4) $$

with $x := (x_1, \dots, x_n)^\top$ is the vector of interacting species, $\operatorname{diag} x$ denotes in short the diagonal matrix with diagonal coefficients $x_i$, $A = (a_{ij})$ is the matrix of interaction coefficients and $r = (r_1, \dots, r_n)^\top$ is the individual growth vector without interaction. The constant vector $\mathcal{E} = (\epsilon_1, \dots, \epsilon_n)^\top$ describes the effect of probiotic or antibiotic agents to reduce the population of the infecting agent, which can be taken as the $x_1$-population. The optimal control problem is a *Mayer problem*: $\min_{u(\cdot)} x_1(t_f)$, which can be formulated in the dual form as: $\min_{u(\cdot)} t_f$, $x_1(t_f) = d$, where $d$ is a desired amount of population $x_1$ at final time $t_f$, taken as a parameter. More generally we shall consider any Mayer problem: $\min_{u(\cdot)} \phi(x(t_f))$ with $\phi$ is a smooth mapping from $M$ to $\mathbb{R}$.

The free motion was studied by the historical contributors Lotka and Volterra, generalizing the prey–predator model for populations dynamics [21]. The interest of the model is to get easily computable equilibria using linear analysis only. In the regular case where $A$ is invertible, the *interior equilibrium* is $x_e := -A^{-1}r$ and the dynamics is related in particular to interaction of $x_e$ with the trivial boundary equilibrium 0. A large amount of literature was devoted to analyze the dynamics in the 2d or 3d case starting from the seminal complete analysis in the prey–predator model. In particular in the case of competitive systems, the achievement in the 3d–case was the explicit analysis of the limit set of the origin to determine the boundary of the basin of repulsion of this point. Related references are [11, 24, 19] and in particular the May–Leonard model [15] will serve as case study.

An important very recent application of controlled Lotka–Volterra model comes from the control of complex microbiota, see the work by Jones et al [12] based on the Stein et al study [20] to model the C. difficile infection of the intestinal microbiote. It contains 11 interaction species and can admit up to $2^{11}$ equilibria. Hence the analysis is a huge task. More modestly we shall concentrate our computation to specific 2d and 3d models to illustrate our techniques. Besides the biological applications, the controlled Lotka–Volterra model allows algebraic computations of the forced equilibria, which are by construction shifted equilibria of the free equilibria and hence in fine one can compute in this case, the spectrum of corresponding ray abnormal solutions. We show that this information is relevant in an applicative context to analyze a 3d May and Leonard model corresponding to a certain microbiote with three species using antibiotic and probiotic. Recent contribution [7] presents a protocol based on a MPC tracking algorithm of a trajectory computed with a direct method. The first consists in reaching the invariant manifold – carrying simplex – in minimum time using probiotic only. The second phase of the treatment focuses on the asymptotic dynamics in the 2d carrying simplex, where the aim is to reach an healthy region using antibiotic only. This article presents an alternative MPC algorithm for the sec-
second phase in which we reach an healthy region in minimum time (see also [6, 8]).

This article is organized in four sections. In section 2, we make a recap of properties of singular trajectories and their parametrization, thanks to the Maximum Principle as an Hamiltonian dynamics in the projectivized cotangent bundle. We introduce the basic concepts of transversality theory in the jets space [14] to formulate in a neat geometric framework the main result of this article concerning feedback invariants. In section 3, we apply this theorem to classify controlled Lotka–Volterra dynamics in dimension 2 or 3. We point out some applications to controlled stability and time minimal syntheses. The final section 4 relates those calculations based on the Maximum Principle in the permanent case to the computations of optimal solution of vermin reduction in the sampled-data control frame.

## 2. A brief recap about singular trajectories and main result

This section is based on the results of [4] and [5], which are briefly presented to deduce the classification result. It combines the techniques from transversality theory in the jets space, see [14], and elementary geometric invariant theory [10].

### 2.1. Notations and general results

Let $M$ be a $\sigma$-compact $C^\infty$ (smooth) manifold of dimension $n \ge 2$. We introduce the following notations:

* $TM$ is the tangent space of $M$ and $T_x M$ is the tangent space at $x \in M$.

* $T^*M$ is the cotangent space of $M$ and $T_x^* M$ is the cotangent space at $x \in M$. The null section of $T^*M$ is denoted by $0$ and $(T^*M)_0 = T^*M \setminus \{0\}$. We note $PT^*M$ the projectivized cotangent space i.e. $PT^*M = T^*M/\mathbb{R}^*$ and $[z]$ is the equivalence class of $z$ in $PT^*M$.

* For any integer $N$, $J^N TM$ is the space of all $N$-jets of vector fields i.e. Taylor expansions up to order $N$ and $JTM = \cup_{N \ge 0} J^N TM$ is the jets space.

* $VF(M)$ is the vector space of all smooth vector fields on $M$ endowed with the Whitney topology.

* We denote by $z = (x, p)$ the canonical coordinates on $(T^*M, \omega)$, where $\omega$ is the Darboux form induced by the Liouville form.

* Take $(X, Y) \in VF(M)$, the Lie bracket is calculated in local coordinates with the convention:

$$ [X, Y](x) := \frac{\partial X}{\partial x}(x)Y(x) - \frac{\partial Y}{\partial x}(x)X(x). $$

* Given any smooth function $H$ defined on an open subset $\Omega$ of $T^*M$, $\tilde{H}$ denotes the Hamiltonian vector field defined by $H$ on $\Omega$, $\tilde{H} := (\partial_p H, -\partial_x H)$. Given $H_1, H_2$ on $M$, $\{H_1, H_2\}$ denotes the Poisson bracket:

$$ \{H_1, H_2\}(z) := dH_1(\tilde{H}_2(z)) = \omega(H_1, H_2). $$

* If $X \in VF(M)$, $H_X$ denotes the Hamiltonian lift: $H_X(x) := p \cdot X(x)$. Given $(X, Y) \in VF(M)$, one has $\{H_X, H_Y\} = H_{[X, Y]}$.
ond phase in which we reach an healthy region in minimum time (see also [6, 8]).

This article is organized in four sections. In section 2, we make a recap of properties of singular trajectories and their parametrization, thanks to the Pontryagin Maximum Principle as an Hamiltonian dynamics in the projectivized cotangent bundle. We introduce the basic concepts of transversality theory in the jets space [14] to formulate in a neat geometric framework the main result of this article concerning feedback invariants. In section 3, we apply this theorem to classify controlled Lotka–Volterra dynamics in dimension 2 or 3. We point out some applications to controlled stability and time minimal syntheses. The final section 4 relates those calculations based on the Maximum Principle in the permanent case to the computations of optimal solution of vermin reduction in the sampled-data control frame.

## 2. A brief recap about singular trajectories and main result

This section is based on the results of [4] and [5], which are briefly presented to deduce the classification result. It combines the techniques from transversality theory in the jets space, see [14], and elementary geometric invariant theory [10].

### 2.1. Notations and general results

Let $M$ be a $\sigma$-compact $C^\infty$ (smooth) manifold of dimension $n \ge 2$. We introduce the following notations:

* $TM$ is the tangent space of $M$ and $T_x M$ is the tangent space at $x \in M$.

* $T^*M$ is the cotangent space of $M$ and $T_x^* M$ is the cotangent space at $x \in M$. The null section of $T^*M$ is denoted by $0$ and $(T^*M)_0 = T^*M \setminus \{0\}$. We note $PT^*M$ the projectivized cotangent space i.e. $PT^*M = T^*M/\mathbb{R}^*$ and $[z]$ is the equivalence class of $z$ in $PT^*M$.

* For any integer $N$, $J^N TM$ is the space of all $N$-jets of vector fields i.e. Taylor expansions up to order $N$ and $JTM = \cup_{N \ge 0} J^N TM$ is the jets space.

* $VF(M)$ is the vector space of all smooth vector fields on $M$ endowed with the Whitney topology.

* We denote by $z = (x, p)$ the canonical coordinates on $(T^*M, \omega)$, where $\omega$ is the Darboux form induced by the Liouville form.

* Take $(X, Y) \in VF(M)$, the Lie bracket is calculated in local coordinates with the convention:

$$[X, Y](x) := \frac{\partial X}{\partial x}(x) Y(x) - \frac{\partial Y}{\partial x}(x) X(x).$$

* Given any smooth function $H$ defined on an open subset $\Omega$ of $T^*M$, $\tilde{H}$ denotes the Hamiltonian vector field defined by $H$ on $\Omega$, $\tilde{H} := (\partial_p H, -\partial_x H)$. Given $H_1, H_2$ on $M$, $\{H_1, H_2\}$ denotes the Poisson bracket:

$$\{H_1, H_2\}(z) := dH_1(\tilde{H}_2(z)) = \omega(H_1, H_2).$$

* If $X \in VF(M)$, $H_X$ denotes the Hamiltonian lift: $H_X(x) := p \cdot X(x)$. Given $(X, Y) \in VF(M)$, one has $\{H_X, H_Y\} = H_{[X, Y]}$.
• Finally, to each pair $(X, Y)$ of vector fields on $M$, we associate the single-input control affine system:

$$ \frac{dx}{dt}(t) = X(x(t)) + u(t) Y(x(t)), \quad u(t) \in \mathbb{R}, \quad x \in M. \qquad (5) $$

The study of time-minimal trajectory of (5) leads to introduce the extremal trajectories: $(z, u) : [0, T] \to T^*M \times \mathbb{R}, T > 0$, such that

(1) $z$ is absolutely continuous, $u$ is measurable and bounded,

(2) $z(t) \neq 0$ (0 being the null section) for all $t \in [0, T]$,

(3) $\frac{dz}{dt}(t) = \vec{H}_X(z(t)) + u(t) \vec{H}_Y(z(t))$ for a.e. $t \in [0, T]$,

(4) $\vec{H}_X(z(t)) + u(t) \vec{H}_Y(z(t)) = \max_{v \in \mathbb{R}} H_X(z(t)) + v \vec{H}_Y(z(t))$ for a.e. $t \in [0, T]$ or equivalently

$$ H_Y(z(t)) = 0, \quad \text{for all } t \in [0, T] \qquad (6) $$

since $z \mapsto H_Y(z)$ is continuous.

**Definition 2.1.** A curve $(z, u) : [0, T] \to T^*M \times \mathbb{R}$ satisfying the above conditions (1)–(4) is called a singular extremal and its projection $(\Pi_M(z, u), u) = (x, u)$ is a singular trajectory. Denoting $z = (x, p)$, $p$ is then the adjoint (non zero) vector.

**Proposition 2.2.** Let $(z, u) = (x, p, u)$ be a singular extremal on $[0, T]$. Then:

• The Fréchet derivative of the fixed time extremity mapping $E^{x(0),tf}$ along $(x(\cdot), u(\cdot))$ is given by the linear dynamics:

$$ \dot{\delta}_{1}x(t) = A(t)\delta_{1}x(t) + u(t)b(t) \qquad (7) $$

with $A = \frac{\partial X}{\partial x}(x(t)) + u(t) \frac{\partial Y}{\partial x}(x(t))$ and $b(t) = Y(x(t))$ with initial condition $\delta_{1}x(0) = 0$.

• The adjoint vector $p(.)$ is orthogonal to the image of the Fréchet derivative called the first order Pontryagin space vector.

One can easily calculate with the Maximum Principle many singular trajectories $(x, u)$. Indeed deriving twice with respect to $t$ the equation (6) we get for a.e. $t \in [0, T]$:

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

**Definition 2.3.** A singular extremal $(z, u)$ on $[0, T]$ is called of minimal order if $\mathcal{R} = \{t \in [0, T], \{H_X, H_Y\}, H_Y\}(z(t)) \neq 0$ is dense in $[0, T]$.

**Proposition 2.4.** Let $(z, u)$ be a singular extremal and $\mathcal{R}$ be a non empty set, then:

• $z$ restricted to $\mathcal{R}$ is smooth,

• the set $\Sigma' = \{z, H_Y(z) = \{H_Y, H_X\}(z) = 0\}$ is invariant for the singular dynamics, which is given on $\Sigma'$ by:

$$ \frac{dz}{dt}(t) = \vec{H}_X(z(t)) + u_s(z(t)) \vec{H}_Y(z(t)) \qquad (9) $$
• Finally, to each pair $(X, Y)$ of vector fields on $M$, we associate the single-input control affine system:

$$ \frac{dx}{dt}(t) = X(x(t)) + u(t) Y(x(t)), \quad u(t) \in \mathbb{R}, \quad x \in M. \qquad (5) $$

The study of time-minimal trajectory of (5) leads to introduce the extremal trajectories: $(z, u) : [0, T] \to T^*M \times \mathbb{R}, T > 0$, such that

(1) $z$ is absolutely continuous, $u$ is measurable and bounded,

(2) $z(t) \neq 0$ (0 being the null section) for all $t \in [0, T]$,

(3) $\frac{dz}{dt}(t) = \vec{H}_X(z(t)) + u(t) \vec{H}_Y(z(t))$ for a.e. $t \in [0, T]$,

(4) $\vec{H}_X(z(t)) + u(t) \vec{H}_Y(z(t)) = \max_{v \in \mathbb{R}} H_X(z(t)) + v \vec{H}_Y(z(t))$ for a.e. $t \in [0, T]$ or equivalently

$$ H_Y(z(t)) = 0, \quad \text{for all } t \in [0, T] \qquad (6) $$

since $z \mapsto H_Y(z)$ is continuous.

**Definition 2.1.** A curve $(z, u) : [0, T] \to T^*M \times \mathbb{R}$ satisfying the above conditions (1)–(4) is called a singular extremal and its projection $(\Pi_M(z, u), u) = (x, u)$ is a singular trajectory. Denoting $z = (x, p)$, $p$ is then the adjoint (non zero) vector.

**Proposition 2.2.** Let $(z, u) = (x, p, u)$ be a singular extremal on $[0, T]$. Then:

• The Fréchet derivative of the fixed time extremity mapping $E^{x(0),tf}$ along $(x(\cdot), u(\cdot))$ is given by the linear dynamics:

$$ \dot{\delta}_{1}x(t) = A(t)\delta_{1}x(t) + u(t)b(t) \qquad (7) $$

with $A = \frac{\partial X}{\partial x}(x(t)) + u(t) \frac{\partial Y}{\partial x}(x(t))$ and $b(t) = Y(x(t))$ with initial condition $\delta_{1}x(0) = 0$.

• The adjoint vector $p(\cdot)$ is orthogonal to the image of the Fréchet derivative called the first order Pontryagin space vector.

One can easily calculate with the Maximum Principle many singular trajectories $(x, u)$. Indeed deriving twice with respect to $t$ the equation (6) we get for a.e. $t \in [0, T]$:

$$ 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. \qquad (8) $$

**Definition 2.3.** A singular extremal $(z, u)$ on $[0, T]$ is called of minimal order if $\mathcal{R} = \{t \in [0, T], \{H_X, H_Y\}, H_Y\}(z(t)) \neq 0$ is dense in $[0, T]$.

**Proposition 2.4.** Let $(z, u)$ be a singular extremal and $\mathcal{R}$ be a non empty set, then:

• $z$ restricted to $\mathcal{R}$ is smooth,

• the set $\Sigma' = \{z, H_Y(z) = \{H_Y, H_X\}(z) = 0\}$ is invariant for the singular dynamics, which is given on $\Sigma'$ by:

$$ \frac{dz}{dt}(t) = \vec{H}_X(z(t)) + u_s(z(t)) \vec{H}_Y(z(t)) \qquad (9) $$
where $u_s(z)$ is the singular dynamics feedback:

$$u_s(z) = \frac{\{\{H_X, H_Y\}, H_X\}(z)}{\{\{H_Y, H_X\}, H_Y\}(z)}. \qquad (10)$$

**Proposition 2.5.** Let $(X, Y)$ be a pair such that the open subset $\Omega$ of all $z \in (T^*M)_0$ such that $\{z, \{\{H_X, H_Y\}, H_Y\}(z) \neq 0\}$ is not empty. Let $H_s: \Omega \to \mathbb{R}$ be the true Hamiltonian:

$$H_s := H_X + \frac{\{\{H_X, H_Y\}, H_X\}(z)}{\{\{H_Y, H_X\}, H_Y\}(z)} H_Y \qquad (11)$$

restricted to the set $\Sigma' : \{z, H_Y(z) = \{H_Y, H_X\}(z) = 0\}$. Then there exist an open set in $VF(M) \times VF(M)$ such that for any couple $(X, Y)$ in this set, $\Omega$ is open and dense and the set of all $z \in \Sigma'$ is a codimension 2 symplectic manifold of $\Omega$ for the induced symplectic form: $\omega_{|\Sigma'}$.

**Notations 2.6.** Restricting to such pairs $\mathcal{G}$ (called good pairs), denote:

* $\Sigma : \{z, H_Y(z) = 0\}$ the switching surface.

* $\tilde{H}_s$ the Hamiltonian dynamics restricted to $\Sigma' \subset \Sigma$ with $\Sigma' : H_Y(z) = \{H_Y, H_X\}(z) = 0$.

* Let $\lambda_s$ be the mapping $\mathcal{G} \ni (X, Y) \mapsto (\tilde{H}_s, \Sigma)$.

## 2.2. Action of the feedback group $G_f$ on the set of good pairs

We shall briefly recall the results of [4]. First of all, borrowed from elementary geometric invariant theory [10], we have.

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

The action of $(\varphi, \alpha, \beta) \in G_f$ can be lifted as the action of Matthieu symplectomorphisms $\tilde{\varphi}$ on $T^*M$ defined in canonical coordinates by $x = \varphi(y)$, $p = q^{\frac{\partial \varphi^{-1}}{\partial y}}$. The action of $(\varphi, \alpha, \beta)$ on $(\tilde{H}_s, \Sigma)$ being reduced to the action of $\varphi$ only. Also note that $\Sigma : H_X(z) = 0$ codes the distribution: $x \mapsto \text{span} Y(x)$.

**Proposition 2.8.** Restricting to good pairs $(X, Y)$, $\lambda_s$ is a covariant, that is the following diagram is commutative:

$$ \begin{array}{ccc} \mathcal{G} \ni (X, Y) & \xrightarrow{\lambda_s} & \lambda_s(X, Y) \\ \xleftarrow{G_f} & \circled{\overline{\quad}} & \xleftarrow{G_f} \\ \mathcal{G} \ni (X', Y') & \xrightarrow{\lambda_s} & \lambda_s(X', Y') \end{array} . $$

**Definition 2.9.** Let $(z, u)$ be a singular extremal on $[0, T]$ of minimal order. The singular extremal is called strict if the adjoint vector is unique up to a scalar (that is unique in $PT^*M$) on $[0, T]$. In the strict case, a singular trajectory $(x, u) = (\Pi_M(z), u)$ on $[0, T]$ is called

(1) Abnormal or exceptional if for every $t \in [0, T]$, $H_X(z(t)) = 0$,
(2) Hyperbolic if for every $t \in [0, T]$, $H_X(z(t)) \{\{H_Y, H_X\}, H_Y\}(z(t)) > 0$,

(3) Elliptic if for every $t \in [0, T]$, $H_X(z(t)) \{\{H_Y, H_X\}, H_Y\}(z(t)) < 0$.

According to the high order Maximum Principle [13], the hyperbolic trajectories are candidates to time-minimal control, while elliptic trajectories are candidates to time-maximal control, while abnormal (exceptional) can be both.

## 2.3. The 2d-case

The 2d-case is a specific situation but can be used to illustrate the general result avoiding technical difficulties. Without losing any generality, one can take $M = \mathbb{R}^2$ and we denote by $x = (x_1, x_2)$ the coordinates while $p = (p_1, p_2)$ denotes the non zero adjoint vector. We introduce the following determinantal sets called respectively the singular locus $\mathcal{S}$ and the collinearity locus $\mathcal{C}$:

$$ \mathcal{S} : \{x, \det(Y, [Y, X])(x) = 0\}, \quad \mathcal{C} : \{x, \det(Y, X)(x) = 0\}. $$

One will assume that $Y$ is not vanishing so that one can choose (local) coordinates with $Y = \frac{\partial}{\partial x_2}$. Moreover we assume that $\mathcal{S}$ and $\mathcal{C}$ are regular and intersect transversally at the point 0. If $D := \det(Y, [[Y, X], Y])$ is not vanishing when restricted to $\mathcal{S}$, the singular control is given by the feedback: $u_s(x) = -\frac{D'(x)}{D(x)}|_{\mathcal{S}}$, where $D' := \det(Y, [[Y, X], X])$. Moreover using $H_Y = 0$, the adjoint vector is such that $p_2 = 0$ identically.

We choose coordinates preserving $Y$ so that $\mathcal{S}$ coincides with the $x_1$-axis. In a neighbourhood of 0, using the action of the feedback group, the system reads:

$$ \begin{aligned} \dot{x}_1 &= \lambda x_1 - x_2^2 + o_{x_1}(x_2^2), \\ \dot{x}_2 &= u, \end{aligned} \qquad (12) $$

where $o_{x_1}(x_2^2)$ represents a term of order $\ge 3$ in the jets space of $(X, Y)$ along the singular line, identified to the $x_1$-axis (the singular dynamics being feedback equivalent to the linear dynamics).

Assuming $\lambda \neq 0$ and restricting the dynamics to singular line identified to $x_2 = 0$, straightforward computation gives us that the singular dynamics is $\dot{x}_1 = \lambda x_1$ and is foliated by the abnormal points $x_1 = 0$, the hyperbolic arc in $x_1 > 0$ and the elliptic arc in $x_1 < 0$.

The adjoint dynamics using $p_2 = 0$ is defined by the adjoint system $\dot{p}_1 = -\lambda p_1 + o_{x_1}(x_2)$. Hence linearized adjoint dynamics reads: $\dot{p}_1 = -\lambda p_1$.

In particular, the abnormal singular point $x = 0$ lifts into a ray in the projective bundle defined by $p_2 = 0$ and $p_1(t) = e^{-\lambda t} p_1(0)$. The linearized Hamiltonian dynamics takes the form

$$ \dot{x}_1 = \lambda x_1, \quad \dot{p}_1 = -\lambda p_1. $$

Clearly $\lambda$ is a feedback invariant and using the projective coordinates $v = x_1/p_1$, leads to the dynamics $\dot{v} = 2\lambda v$ so that $2\lambda$ is a *projectivized feedback invariant*.

From this analysis one deduces the following.

**Theorem 2.10.** In the 2d-case there exist a nonempty open set of pairs $(X, Y)$ for the Whitney topology such that:
(2) Hyperbolic if for every $t \in [0, T]$, $H_X(z(t)) \{\{H_Y, H_X\}, H_Y\}(z(t)) > 0$,

(3) Elliptic if for every $t \in [0, T]$, $H_X(z(t)) \{\{H_Y, H_X\}, H_Y\}(z(t)) < 0$.

According to the high order Maximum Principle [13], the hyperbolic trajectories are candidates to time-minimal control, while elliptic trajectories are candidates to time-maximal control, while abnormal (exceptional) can be both.

### 2.3. The 2d-case

The 2d-case is a specific situation but can be used to illustrate the general result avoiding technical difficulties. Without losing any generality, one can take $M = \mathbb{R}^2$ and we denote by $x = (x_1, x_2)$ the coordinates while $p = (p_1, p_2)$ denotes the non zero adjoint vector. We introduce the following determinantal sets called respectively the singular locus $\mathcal{S}$ and the collinearity locus $\mathcal{C}$:

$$ \mathcal{S} : \{x, \det(Y, [Y, X])(x) = 0\}, \quad \mathcal{C} : \{x, \det(Y, X)(x) = 0\}. $$

One will assume that $Y$ is not vanishing so that one can choose (local) coordinates with $Y = \frac{\partial}{\partial x_2}$. Moreover we assume that $\mathcal{S}$ and $\mathcal{C}$ are regular and intersect transversally at the point 0. If $D := \det(Y, [[Y, X], Y])$ is not vanishing when restricted to $\mathcal{S}$, the singular control is given by the feedback: $u_s(x) = -\frac{D'(x)}{D(x)}|_{\mathcal{S}}$, where $D' := \det(Y, [[Y, X], X])$. Moreover using $H_Y = 0$, the adjoint vector is such that $p_2 = 0$ identically.

We choose coordinates preserving $Y$ so that $\mathcal{S}$ coincides with the $x_1$-axis. In a neighbourhood of 0, using the action of the feedback group, the system reads:

$$ \begin{aligned} \dot{x}_1 &= \lambda x_1 - x_2^2 + o_{x_1}(x_2^2), \\ \dot{x}_2 &= u, \end{aligned} \qquad (12) $$

where $o_{x_1}(x_2^2)$ represents a term of order $\ge 3$ in the jets space of $(X, Y)$ along the singular line, identified to the $x_1$-axis (the singular dynamics being feedback equivalent to the linear dynamics).

Assuming $\lambda \ne 0$ and restricting the dynamics to singular line identified to $x_2 = 0$, straightforward computation gives us that the singular dynamics is $\dot{x}_1 = \lambda x_1$ and is foliated by the abnormal points $x_1 = 0$, the hyperbolic arc in $x_1 > 0$ and the elliptic arc in $x_1 < 0$.

The adjoint dynamics using $p_2 = 0$ is defined by the adjoint system $\dot{p}_1 = -\lambda p_1 + o_{x_1}(x_2)$. Hence linearized adjoint dynamics reads: $\dot{p}_1 = -\lambda p_1$.

In particular, the abnormal singular point $x = 0$ lifts into a ray in the projective bundle defined by $p_2 = 0$ and $p_1(t) = e^{-\lambda t} p_1(0)$. The linearized Hamiltonian dynamics takes the form

$$ \dot{x}_1 = \lambda x_1, \quad \dot{p}_1 = -\lambda p_1. $$

Clearly $\lambda$ is a feedback invariant and using the projective coordinates $v = x_1/p_1$, leads to the dynamics $\dot{v} = 2\lambda v$ so that $2\lambda$ is a *projectivized feedback invariant*.

From this analysis one deduces the following.

**Theorem 2.10.** In the 2d-case there exist a nonempty open set of pairs $(X, Y)$ for the Whitney topology such that:
(1) Equilibria of the singular dynamics are isolated and reduced to abnormal equilibria $x_0$ defining the singular lines in M with dynamics $\dot{x} = \lambda(x_0)x$.

(2) Each such point defines a ray $z(t) = e^{-\lambda(x_0)t} z_0$, $z_0 \in \ell$ ($\ell$ being a line in $T_{x_0}^*M$).

(3) The singular line is foliated into the abnormal equilibrium and hyperbolic, elliptic arcs.

(4) The eigenvalue $2\lambda(x_0)$ is a projective feedback invariant.

**Remark 1.** Introducing the clock one form: $\alpha = pdx$ outside the collinear set, the singular lines are the zero of $d\alpha = dp \wedge dx$.

## 2.4. The 3d-case

The previous planar case can be generalized to the three dimensional case, which is a very rich situation and can be treated similarly, paving the road to the general case.

Using (8), we introduce the determinantal mappings:

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

Assume that $D$ is non zero, one has the following proposition.

**Proposition 2.11.**

* Singular trajectories of minimal order are solutions of the dynamics:

$$ \frac{dx}{dt}(t) = X_s(x) := X(x) + u_s(x) Y(x), \qquad (14) $$

where the singular controls $u_s$ is given by the feedback:

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

* The sets $D'' = 0$, $DD'' > 0$ and $DD'' < 0$ are invariant for the singular dynamics and correspond respectively to abnormal, hyperbolic and elliptic trajectories.

* The adjoint vector in the projective space $PT^*\mathbb{R}^3$ is uniquely defined by the relations $H_Y = \{H_Y, H_X\} = 0$.

Introducing the clock form $\alpha = pdx$ defined by the relations: $H_X = 1$, $H_Y = \{H_Y, H_X\} = 0$, outside the abnormal locus $D'' = 0$, the singular trajectories are the characteristics of $d\alpha$.

Clearly we have:

**Lemma 2.12.** The feedback group $G_f$ acts on the singular dynamics by change of coordinates only.

From which we deduce:

**Proposition 2.13.** The singular points of (14) are abnormal equilibria and denoting by $J$ the Jacobian matrix of (14) at such points then the whole spectrum $\sigma(J)$ is feedback invariant. Each equilibrium point $x_0$ defines a ray solution in the projective cotangent space given by $z(t) = e^{-\lambda t}z(0)$, where $x_0$ is the canonical projection of $z(0)$, $\lambda \in \sigma(J)$.
(1) Equilibria of the singular dynamics are isolated and reduced to abnormal equilibria $x_0$ defining the singular lines in M with dynamics $\dot{x} = \lambda(x_0)x$.

(2) Each such point defines a ray $z(t) = e^{-\lambda(x_0)t} z_0$, $z_0 \in \ell$ ($\ell$ being a line in $T_{x_0}^*M$).

(3) The singular line is foliated into the abnormal equilibrium and hyperbolic, elliptic arcs.

(4) The eigenvalue $2\lambda(x_0)$ is a projective feedback invariant.

**Remark 1.** Introducing the clock one form: $\alpha = pdx$ outside the collinear set, the singular lines are the zero of $d\alpha = dp \wedge dx$.

## 2.4. The 3d-case

The previous planar case can be generalized to the three dimensional case, which is a very rich situation and can be treated similarly, paving the road to the general case.

Using (8), we introduce the determinantal mappings:

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

Assume that $D$ is non zero, one has the following proposition.

**Proposition 2.11.**

* Singular trajectories of minimal order are solutions of the dynamics:

$$ \frac{dx}{dt}(t) = X_s(x) := X(x) + u_s(x) Y(x), \quad (14) $$

where the singular controls $u_s$ is given by the feedback:

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

* The sets $D'' = 0$, $DD'' > 0$ and $DD'' < 0$ are invariant for the singular dynamics and correspond respectively to abnormal, hyperbolic and elliptic trajectories.

* The adjoint vector in the projective space $PT^*\mathbb{R}^3$ is uniquely defined by the relations $H_Y = \{H_Y, H_X\} = 0$.

Introducing the clock form $\alpha = pdx$ defined by the relations: $H_X = 1$, $H_Y = \{H_Y, H_X\} = 0$, outside the abnormal locus $D'' = 0$, the singular trajectories are the characteristics of $d\alpha$.

Clearly we have:

**Lemma 2.12.** The feedback group $G_f$ acts on the singular dynamics by change of coordinates only.

From which we deduce:

**Proposition 2.13.** The singular points of (14) are abnormal equilibria and denoting by $J$ the Jacobian matrix of (14) at such points then the whole spectrum $\sigma(J)$ is feedback invariant. Each equilibrium point $x_0$ defines a ray solution in the projective cotangent space given by $z(t) = e^{-\lambda t}z(0)$, where $x_0$ is the canonical projection of $z(0)$, $\lambda \in \sigma(J)$.
2.5. The general case $n \ge 3$

264 In this case, more technicality is necessary, but the result follows mainly from the proof of
265 Lemma 1 in the reference [5].

267 One needs the following.

268 **2.5.1. Ad–condition and the bad set of infinite codimension**

269 Given a pair (X, Y) of vector fields denote by adX the operator defined by: ad$^0$X(Y) = Y,
270 adX(Y) = [X, Y] and inductively: ad$^k$X(Y) = [ad$^{k-1}$X · Y, Y] for k ≥ 2. Denote by adH$_X$ the
271 induced operator on the Hamiltonians: adH$_X$ · H$_Y$ = {H$_X$, H$_Y$}.

272 **Definition 2.14.** For $N$ large enough we define the following subset of $J^N TM \times J^N TM$:

273 (1) $B'_{\ell}(N)$ is the subset of all couples $(j_X^N X, j_X^N Y)$ such that:
274 dim span $\{X(x), Y(x), [X, Y](x)\} \le 1$.

275 (2) $\hat{B}_{\ell}''$ is the subset of $J^N TM \times J^N TM \times \mathbb{R}$ of all triples $(j_X^N X, j_X^N Y, a)$ such that:

276 (i) $Y(x) \neq 0$

277 (ii) $X(x) = aY(x)$

278 (iii) $\dim \text{span}\{\text{ad}^k G_a(Y)(x), 0 \le k \le n-1, [[X, Y], Y](x)\} < n$, where $G_a = X - aY$.

279 (3) Denote by $B_{\ell}''(N)$ the canonical projection of $\hat{B}_{\ell}''(N)$ on $J^N TM \times J^N TM$. Let $B_{\ell}(N) = B_{\ell}' \cup B_{\ell}''(N)$.
280

281 From the proof of Lemma 1 in [5], one has:

282 **Proposition 2.15.** Let (X, Y) be a pair in VF(M) × VF(M) such that for all $x \in M$,
283 $(j_X^N X, j_X^N Y) \notin B_\ell(N)$. Then:

284 (1) Let $(z, u)$ be a singular extremal on $[0, T]$ such that $\dim \text{span}\{X(x(t)), Y(x(t))\} \le 1$, where $x$ is the canonical projection of $z$. Then $x(\cdot)$ is constant and is an abnormal singular arc reduced to a point $x_0$.

287 (2) The extremal $z(\cdot)$ is of minimal order and strict and there exist a line $\ell \in T_{x_0}^* M$ so that
288 $z(t) = e^{-\lambda t} z_0$, $z_0 \in \ell$ and the control $u$ is constant a.e. Hence $z(\cdot)$ is a ray solution.

289 (3) The point $x_0$ is a forced abnormal equilibrium contained in the collinear set and
290 let $(A, b)$ be the linearized dynamics (7) at $(x_0, u)$, then $(A, b)$ are constant and
291 $\text{span}\{b, Ab, \dots, A^{n-1}b\}$ is of codimension one. The line $\ell$ is orthogonal to this space.

292 **Theorem 2.16** (Main theorem for $n \ge 3$). There exist an open set of pairs $(X, Y)$ for the $C^\infty$–
293 Whitney topology such that:

294 (1) Every singular extremal pair $(z, u)$ on $[0, T]$ is with minimal order and strict.

295 (2) Every singular trajectory reduced to a point $x_0$ is abnormal and $x_0$ is a forced abnor-
296 mal equilibrium associated to a ray solution $z(t) = e^{-\lambda t}z(0)$, $z(0) \in \ell$, with $\ell$ a line in
297 $PT_{x_0}^* M$.

298 (3) Every abnormal equilibrium is isolated.

299 (4) The pair $(A, b)$ is not controllable and $\text{span}\{b, Ab, \dots, A^{n-1}b\} = n-1$. The unique un-
300 controllable mode satisfies the dynamics $\dot{x} = \lambda x$.

301 (5) The whole spectrum of the linearized dynamics in $PT_{x_0}^* M$ is feedback invariant.

302 **Proof.** The proof follows from [5] and standard linear geometric control theory. More pre-
303 cisely, from linear theory, if the pair $(A, b)$ is controllable, the pole placement theorem [23]
**2.5. The general case $n \ge 3$**

In this case, more technicality is necessary, but the result follows mainly from the proof of Lemma 1 in the reference [5].

One needs the following.

**2.5.1. Ad-condition and the bad set of finite codimension**

Given a pair $(X, Y)$ of vector fields denote by $\mathrm{ad}X$ the operator defined by: $\mathrm{ad}^0 X(Y) = Y$, $\mathrm{ad}X(Y) = [X, Y]$ and inductively: $\mathrm{ad}^k X(Y) = [\mathrm{ad}^{k-1} X \cdot Y, Y]$ for $k \ge 2$. Denote by $\mathrm{ad}H_X$ the induced operator on the Hamiltonians: $\mathrm{ad}H_X \cdot H_Y = \{H_X, H_Y\}$.

**Definition 2.14.** For $N$ large enough we define the following subset of $J^N TM \times J^N TM$:

(1) $B'_\ell(N)$ is the subset of all couples $(j_X^N X, j_X^N Y)$ such that: $\dim \mathrm{span}\{X(x), Y(x), [X, Y](x)\} \le 1$.

(2) $\hat{B}''_\ell$ is the subset of $J^N TM \times J^N TM \times \mathbb{R}$ of all triples $(j_X^N X, j_X^N Y, a)$ such that:
(i) $Y(x) \neq 0$

(ii) $X(x) = aY(x)$

(iii) $\dim \mathrm{span}\{\mathrm{ad}^k G_a(Y)(x), 0 \le k \le n-1, [[X, Y], Y](x)\} < n$, where $G_a = X - aY$.

(3) Denote by $B''_\ell(N)$ the canonical projection of $\hat{B}''_\ell(N)$ on $J^N TM \times J^N TM$. Let $B_\ell(N) = B'_\ell \cup B''_\ell(N)$.

From the proof of Lemma 1 in [5], one has:

**Proposition 2.15.** Let $(X, Y)$ be a pair in $VF(M) \times VF(M)$ such that for all $x \in M$, $(j_X^N X, j_X^N Y) \notin B_\ell(N)$. Then:

(1) Let $(z, u)$ be a singular extremal on $[0, T]$ such that $\dim \mathrm{span}\{X(x(t)), Y(x(t))\} \le 1$, where $x$ is the canonical projection of $z$. Then $x(\cdot)$ is constant and is an abnormal singular arc reduced to a point $x_0$.

(2) The extremal $z(\cdot)$ is of minimal order and strict and there exist a line $\ell \in T_{x_0}^* M$ so that $z(t) = e^{-\lambda t}z_0$, $z_0 \in \ell$ and the control $u$ is constant a.e. Hence $z(\cdot)$ is a ray solution.

(3) The point $x_0$ is a forced abnormal equilibrium contained in the collinear set and let $(A, b)$ be the linearized dynamics (7) at $(x_0, u)$, then $(A, b)$ are constant and span $\{b, Ab, \dots, A^{n-1}b\}$ is of codimension one. The line $\ell$ is orthogonal to this space.

**Theorem 2.16 (Main theorem for $n \ge 3$).** There exist an open set of pairs $(X, Y)$ for the $C^\infty$-Whitney topology such that:

(1) Every singular extremal pair $(z, u)$ on $[0, T]$ is with minimal order and strict.

(2) Every singular trajectory reduced to a point $x_0$ is abnormal and $x_0$ is a forced abnormal equilibrium associated to a ray solution $z(t) = e^{-\lambda t}z(0)$, $z(0) \in \ell$, with $\ell$ a line in $PT_{x_0}^* M$.

(3) Every abnormal equilibrium is isolated.

(4) The pair $(A, b)$ is not controllable and $\mathrm{span}\{b, Ab, \dots, A^{n-1}b\} = n-1$. The unique uncontrollable mode satisfies the dynamics $\dot{x} = \lambda x$.

(5) The whole spectrum of the linearized dynamics in $PT_{x_0}^* M$ is feedback invariant.

**Proof.** The proof follows from [5] and standard linear geometric control theory. More precisely, from linear theory, if the pair $(A, b)$ is controllable, the pole placement theorem [23]
asserts that one can assign every spectrum using a linear gain control: $u = kx$. On the opposite, the non controllable modes cannot be modified. In the strict case, there exist a single uncontrollable mode. □

### 2.5.2. Application to generic properties of pairs $(X, Y, \varphi)$.

Practically one aim of biological models is to describe equilibria and their stability property. From control point of view, one can choose $Y$ in a given class and the Mayer cost to be maximized defines a family of terminal manifolds $N(d)$ of codimension one given as the level sets: $\{x, \varphi(x) = d\}$. Take a point $x \in N$, we denote by $n(x)$ the normal vector to $N$ at $x$, and accessibility properties near the terminal point can be classified in a generic context for the $C^\infty$-Whitney topology on the triples: $(X, Y, N)$ using the so-called transversality condition. This amounts to maximize the scalar product: $n \cdot \dot{x} = n \cdot (X + uY)$ for all $u$ in a feasible interval, which can be taken as $[-1, +1]$.

This leads to stratify the final target $N(d)$ into:

* The switching locus $\Sigma: \{x \in N, n(x) \cdot Y(x) = 0\}$

* The singular locus $\mathcal{S}: \{x \in N, n(x) \cdot [Y, X](x) = 0\}$

* The exceptional locus $\mathcal{E}: \{x \in N, n(x) \cdot X(x) = 0\}$.

In particular this leads, for fixed pairs $(X, Y)$, to consider specific critical manifolds. In our study we consider manifolds $N$ for which: $\Sigma \cap \mathcal{S} \cap \mathcal{E} \neq \emptyset$, to analyze

* time-minimal control syntheses,

* controlled stability.

This will be studied in the next sections using algebraic computations of equilibria on the Lotka–Volterra models.

## 3. Controlled Lotka–Volterra model

### 3.1. A brief recap about controlled Lotka–Volterra related to microbiote control

**Definition 3.1.** A controlled Lotka–Volterra dynamics is a control system of the form

$$ \frac{dx}{dt}(t) = (\operatorname{diag} x(t)) (Ax(t) + r + u(t)\mathcal{E}) \quad (16) $$

with $x := (x_1, \dots, x_n)^\top \in \mathbb{R}_{\ge 0}^n$, $A = (a_{ij})$, $r = (r_1, \dots, r_n)^\top$ and $\mathcal{E} = (\epsilon_1, \dots, \epsilon_n)^\top$ are constant vectors, $u(t)$ represents the control intensity, which can be taken in $[0, 1]$. The free dynamics is called *regular* if the matrix $A$ is invertible. One can extend the dynamics to the whole $\mathbb{R}^n$ and the control intensity to the whole $\mathbb{R}$.

An *interior equilibrium* is a point $x$ in $\mathbb{R}_{>0}^n$ such that $Ax + r = 0$ i.e. $x = A^{-1}r$. To each Lotka–Volterra dynamics one can assign in the regular case up to $2^n$ equilibria by considering all the induced Lotka–Volterra models with extinction of at most one species $x_i$ (i.e. $x_i = 0$). For the controlled Lotka–Volterra model, the infecting agent population is denoted $x_1$, hence vermin reduction aims to minimize the $x_1$-population.

**Lemma 3.2.** Consider the dynamics (16) and let $\Omega = (K_1, \dots, K_n)^\top$ be an interior equili-
asserts that one can assign every spectrum using a linear gain control: $u = kx$. On the opposite, the non controllable modes cannot be modified. In the strict case, there exist a single uncontrollable mode. □

### 2.5.2. Application to generic properties of pairs $(X, Y, \varphi)$.

Practically one aim of biological models is to describe equilibria and their stability property. From control point of view, one can choose $Y$ in a given class and the Mayer cost to be maximized defines a family of terminal manifolds $N(d)$ of codimension one given as the level sets: $\{x, \varphi(x) = d\}$. Take a point $x \in N$, we denote by $n(x)$ the normal vector to $N$ at $x$, and accessibility properties near the terminal point can be classified in a generic context for the $C^\infty$-Whitney topology on the triples: $(X, Y, N)$ using the so-called transversality condition. This amounts to maximize the scalar product: $n \cdot \dot{x} = n \cdot (X + uY)$ for all $u$ in a feasible interval, which can be taken as $[-1, +1]$.

This leads to stratify the final target $N(d)$ into:

* The switching locus $\Sigma: \{x \in N, n(x) \cdot Y(x) = 0\}$

* The singular locus $\mathcal{S}: \{x \in N, n(x) \cdot [Y, X](x) = 0\}$

* The exceptional locus $\mathcal{E}: \{x \in N, n(x) \cdot X(x) = 0\}$.

In particular this leads, for fixed pairs $(X, Y)$, to consider specific critical manifolds. In our study we consider manifolds $N$ for which: $\Sigma \cap \mathcal{S} \cap \mathcal{E} \neq \emptyset$, to analyze

* time-minimal control syntheses,

* controlled stability.

This will be studied in the next sections using algebraic computations of equilibria on the Lotka–Volterra models.

## 3. Controlled Lotka–Volterra model

### 3.1. A brief recap about controlled Lotka–Volterra related to microbiote control

**Definition 3.1.** A controlled Lotka–Volterra dynamics is a control system of the form

$$ \frac{dx}{dt}(t) = (\operatorname{diag} x(t)) (Ax(t) + r + u(t)\mathcal{E}) \quad (16) $$

with $x := (x_1, \dots, x_n)^\top \in \mathbb{R}_{\ge 0}^n$, $A = (a_{ij})$, $r = (r_1, \dots, r_n)^\top$ and $\mathcal{E} = (\epsilon_1, \dots, \epsilon_n)^\top$ are constant vectors, $u(t)$ represents the control intensity, which can be taken in $[0, 1]$. The free dynamics is called *regular* if the matrix $A$ is invertible. One can extend the dynamics to the whole $\mathbb{R}^n$ and the control intensity to the whole $\mathbb{R}$.

An *interior equilibrium* is a point $x$ in $\mathbb{R}_{>0}^n$ such that $Ax + r = 0$ i.e. $x = A^{-1}r$. To each Lotka–Volterra dynamics one can assign in the regular case up to $2^n$ equilibria by considering all the induced Lotka–Volterra models with extinction of at most one species $x_i$ (i.e. $x_i = 0$). For the controlled Lotka–Volterra model, the infecting agent population is denoted $x_1$, hence vermin reduction aims to minimize the $x_1$-population.

**Lemma 3.2.** Consider the dynamics (16) and let $\Omega = (K_1, \dots, K_n)^\top$ be an interior equili-
339 *brium.* Then there exist coordinates such that the dynamics (16) can be written:

$$
\dot{v}_i = -(v_i + 1) \left( \sum_{j=1}^{n} a_{ij}^* v_j + u\epsilon_i \right), \quad i = 1, \dots, n. \qquad (17)
$$

340

341 **Proof.** Let $y_i$ be the dimensionless coordinates $y_i = x_i/K_i$, $i = 1, \dots, n$ so that the dynamics

342 (16) with $A \leftarrow -A$ written as $\dot{x}_i = x_i \left( r_i - \sum_{j=1}^{n} a_{ij} x_j + u\epsilon_i \right)$, $i = 1, \dots, n$, becomes

$$
\dot{y}_i = y_i \left( r_i - \sum_{j=1}^{n} (a_{ij} K_i) y_j + u \epsilon_i \right), \quad i = 1, \dots, n. \qquad (18)
$$

343 Denote $A^* = (a_{ij}^*) = (K_i a_{ij})$. By construction, the interior equilibrium is normalized to

344 $\Omega = (1, \dots, 1)$ so that: $r_i = \sum_{j=1}^{n} a_{ij}^*$. Hence (18) becomes

$$
\dot{y}_i = y_i \left( \sum_{j=1}^{n} a_{ij}^* (1 - y_j) + u \epsilon_i \right), \quad i = 1, \dots, n.
$$

345 Therefore if we set $v_i = y_i - 1$ the dynamics takes the form

$$
\dot{v}_i = -(v_i + 1) \left( \sum_{j=1}^{n} a_{ij}^* v_j - u\epsilon_i \right), \quad i = 1, \dots, n. \qquad (19)
$$

346 It can be written shortly as:

$$
\dot{v} = -(\mathrm{diag}(1 + v))(A^* v - u\mathcal{E}),
$$

347 where the equilibrium is normalized to 0.

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

348 Therefore this triggers to consider controlled Lotka-Volterra model of the form
349 $-(\mathrm{diag}(v+1))(Av-u\mathcal{E})$ for which we have the following Lemma.

350 **Lemma 3.3.** Consider the controlled Lotka–Volterra model $\dot{x} = -(\mathrm{diag}(x+1))(Ax-u\mathcal{E})$ with interior equilibrium $x=0$. Denote $X(x) = -(\mathrm{diag}(x+1))A$, $Y(x) = (\mathrm{diag}(x+1))\mathcal{E}$ so that $-A = \frac{\partial X}{\partial x}(0)$ and $Y(0) = \mathcal{E}$. Hence $\mathrm{span}\{\mathrm{ad}^k X \cdot Y(0), k=0, \dots, n-1\} = \mathrm{span}\{\mathcal{E}, A\mathcal{E}, \dots, A^{n-1}\mathcal{E}\}$. Therefore $(x, u) = (0,0)$ is an abnormal strict singular point if and only if the rank of the Kalman matrix $K = [\mathcal{E}, A\mathcal{E}, \dots, A^{n-1}\mathcal{E}]$ is strictly less than $n-1$.

355 **Corollary 3.4.** The pair $(A, b)$ is controllable if and only if $\mathcal{E}$ is a cyclic vector of $A$.

356 This leads to easily characterize every pair $(A, \mathcal{E})$ using a Jordan decomposition of $A$
357 such that $(0,0)$ is an (abnormal) singular trajectory.

358 **Example 3.5.** Take a matrix $A$ with Jordan blocks with equal eigenvalues e.g. the diagonal
359 matrix $A = \mathrm{diag}(\lambda_1, \lambda_1, \lambda_2)$. Then for every vector $\mathcal{E}$ the pair $(A, \mathcal{E})$ is not controllable.
brium. Then there exist coordinates such that the dynamics (16) can be written:

$$ \dot{v}_i = (v_i + 1) \left( \sum_{j=1}^{n} a_{ij}^* v_j + u \mathcal{E} \right), \quad i = 1, \dots, n. \qquad (17) $$

**Proof.** Let $y_i$ be the dimensionless coordinates $y_i = x_i/K_i$, $i = 1, \dots, n$ so that the dynamics (16) with $A \leftarrow -A$ written as $\dot{x}_i = x_i (r_i - \sum_{j=1}^n a_{ij} x_j + u \epsilon_i)$, $i = 1, \dots, n$, becomes

$$ \dot{y}_i = y_i \left( r_i - \sum_{j=1}^{n} (a_{ij} K_j) y_j + u \epsilon_i \right), \quad i = 1, \dots, n. \qquad (18) $$

Denote $A^* = (a_{ij}^*) = (K_i a_{ij})$. By construction, the interior equilibrium is normalized to $\Omega = (1, \dots, 1)$ so that: $r_i = \sum_{j=1}^n a_{ij}^*$. Hence (18) becomes

$$ \dot{y}_i = y_i \left( \sum_{j=1}^{n} a_{ij}^* (1 - y_j) + u \epsilon_i \right), \quad i = 1, \dots, n. $$

Therefore if we set $v_i = y_i - 1$ the dynamics takes the form

$$ \dot{v}_i = (v_i + 1) \left( \sum_{j=1}^{n} a_{ij}^* v_j + u \epsilon_i \right), \quad i = 1, \dots, n. \qquad (19) $$

It can be written shortly as:

$$ \dot{v} = (\mathrm{diag}(1 + v)) (A^* v + u \mathcal{E}), $$

where the equilibrium is normalized to 0. $\square$

Therefore this triggers to consider controlled Lotka-Volterra model of the form $-(\mathrm{diag}(v+1))(Av - u\mathcal{E})$ for which we have the following Lemma.

**Lemma 3.3.** Consider the controlled Lotka-Volterra model $\dot{x} = -(\mathrm{diag}(x+1))(Ax - u\mathcal{E})$ with interior equilibrium $x = 0$. Denote $X(x) = -(\mathrm{diag}(x+1))A$, $Y(x) = (\mathrm{diag}(x+1))\mathcal{E}$ so that $-A = \frac{\partial X}{\partial x}(0)$ and $Y(0) = \mathcal{E}$. Hence $\mathrm{span}\{\mathrm{ad}^k X \cdot Y(0), k = 0, \dots, n-1\} = \mathrm{span}\{\mathcal{E}, A\mathcal{E}, \dots, A^{n-1}\mathcal{E}\}$. Therefore $(x, u) = (0, 0)$ is an abnormal strict singular point if and only if the rank of the Kalman matrix $K = [\mathcal{E}, A\mathcal{E}, \dots, A^{n-1}\mathcal{E}]$ is strictly less than $n-1$.

**Corollary 3.4.** The pair $(A, b)$ is controllable if and only if $\mathcal{E}$ is a cyclic vector of $A$.

This leads to easily characterize every pair $(A, \mathcal{E})$ using a Jordan decomposition of $A$ such that $(0, 0)$ is an (abnormal) singular trajectory.

**Example 3.5.** Take a matrix $A$ with Jordan blocks with equal eigenvalues e.g. the diagonal matrix $A = \mathrm{diag}(\lambda_1, \lambda_1, \lambda_2)$. Then for every vector $\mathcal{E}$ the pair $(A, \mathcal{E})$ is not controllable.
### 3.2. Some examples of case studies

#### 3.2.1. The prey–predator model

This historical model [21], using the $(x, y)$ coordinates, takes the form

$$ \frac{dx}{dt} = x(\lambda_1 + \mu_1 y) + u\epsilon_1 y, \qquad \frac{dy}{dt} = y(\lambda_2 + \mu_2 x) + u\epsilon_2 y, \qquad (20) $$

which covers both the elliptic and hyperbolic case.

The point $\Omega = (-\lambda_2/\mu_2, -\lambda_1/\mu_1) := (K_1, K_2)$ is the interior equilibrium in the quadrant $x, y > 0$ provided that $\lambda_1\mu_1, \lambda_2\mu_2 < 0$. Introduce the dimensionless coordinates $x \leftarrow x/(-\lambda_2/\mu_2)$ and $y \leftarrow y(-\lambda_1/\mu_1)$ we obtain the system:

$$ \dot{x} = \lambda_1 x(1 - y) + u x \epsilon_1, \qquad \dot{y} = \lambda_2 y(1 - x) + u x \epsilon_2, $$

where the interior equilibrium is normalized to $\Omega = (1, 1)$.

The Jacobian matrix at $\Omega$ is

$$ J = \begin{pmatrix} 0 & -\lambda_1 \\ -\lambda_2 & 0 \end{pmatrix} $$

and is elliptic if $\lambda_1\lambda_2 < 0$ and hyperbolic if $\lambda_1\lambda_2 > 0$.

The collinearity locus is $\mathscr{C}: xy[\epsilon_2\lambda_1(1-y) - \epsilon_1\lambda_2(1-x)] = 0$ and contains $\Omega = (1, 1)$ by construction.

The singular locus is $\mathscr{S}: xy[\epsilon_1^2\lambda_2x - \epsilon_2^2\lambda_1y] = 0$ and defines the interior line for $xy \neq 0$ given by $y = \frac{\epsilon_1^2\lambda_2}{\epsilon_2^2\lambda_1}x$. Hence we have:

* *Elliptic case:* $\lambda_1\lambda_2 < 0$. The singular line is not contained in the cone $x, y > 0$.

* *Hyperbolic case:* $\lambda_1\lambda_2 > 0$. The singular line is contained in the cone of positive population $x, y > 0$.

The crucial point of the classification is to compute abnormal points, intersection of such two interior lines. At such point, the analysis of 3.2 applies.

A richer situation making interacting the elliptic and hyperbolic cases is to consider the model:

$$ \frac{dx}{dt} = (1 - u)X(x) + uY(x), \quad u \in [0, 1], $$

with

$$ X = \lambda_1 x(1 - y) \frac{\partial}{\partial x} + \lambda_2 y(1 - x) \frac{\partial}{\partial y}, \quad Y = x(\lambda'_1 + \mu'_1 y) \frac{\partial}{\partial x} + y(\lambda'_2 + \mu'_2 x) \frac{\partial}{\partial y}, $$

which leads to the system: $\dot{x} = X + u(Y - X)$ connecting the elliptic case to the hyperbolic case. It is an extension of models (16). Straightforward computations lead to an interior collinearity locus, which is defined by a quadratic mapping vs. a linear map in the prey–predator model, while the interior singular locus is defined by a cubic mapping vs. a linear map. More generally automatic computations lead to analyze all the 2d–cases.
## 3.2. Some examples of case studies

### 3.2.1. The prey–predator model

This historical model [21], using the $(x, y)$ coordinates, takes the form

$$ \frac{dx}{dt} = x(\lambda_1 + \mu_1 y) + u\epsilon_1 x, \qquad \frac{dy}{dt} = y(\lambda_2 + \mu_2 x) + u\epsilon_2 y, \qquad (20) $$

which covers both the elliptic and hyperbolic case.

The point $\Omega = (-\lambda_2/\mu_2, -\lambda_1/\mu_1) := (K_1, K_2)$ is the interior equilibrium in the quadrant $x, y > 0$ provided that $\lambda_1\mu_1, \lambda_2\mu_2 < 0$. Introduce the dimensionless coordinates $x \leftarrow x/(-\lambda_2/\mu_2)$ and $y \leftarrow y(-\lambda_1/\mu_1)$ we obtain the system:

$$ \dot{x} = \lambda_1 x(1 - y) + u x \epsilon_1, \qquad \dot{y} = \lambda_2 y(1 - x) + u y \epsilon_2, $$

where the interior equilibrium is normalized to $\Omega = (1, 1)$.

The Jacobian matrix at $\Omega$ is

$$ J = \begin{pmatrix} 0 & -\lambda_1 \\ -\lambda_2 & 0 \end{pmatrix} $$

and is elliptic if $\lambda_1\lambda_2 < 0$ and hyperbolic if $\lambda_1\lambda_2 > 0$.

The collinearity locus is $\mathscr{C} : xy[\epsilon_2\lambda_1(1 - y) - \epsilon_1\lambda_2(1 - x)] = 0$ and contains $\Omega = (1, 1)$ by construction.

The singular locus is $\mathscr{S} : xy[\epsilon_1^2\lambda_2x - \epsilon_2^2\lambda_1y] = 0$ and defines the interior line for $xy \neq 0$ given by $y = \frac{\epsilon_1^2\lambda_2}{\epsilon_2^2\lambda_1}x$. Hence we have:

* Elliptic case: $\lambda_1\lambda_2 < 0$. The singular line is not contained in the cone $x, y > 0$.

* Hyperbolic case: $\lambda_1\lambda_2 > 0$. The singular line is contained in the cone of positive population $x, y > 0$.

The crucial point of the classification is to compute abnormal points, intersection of such two interior lines. At such point, the analysis of 3.2 applies.

A richer situation making interacting the elliptic and hyperbolic cases is to consider the model:

$$ \frac{dx}{dt} = (1 - u)X(x) + uY(x), \quad u \in [0, 1], $$

with

$$ X = \lambda_1 x(1 - y) \frac{\partial}{\partial x} + \lambda_2 y(1 - x) \frac{\partial}{\partial y}, \quad Y = x(\lambda'_1 + \mu'_1 y) \frac{\partial}{\partial x} + y(\lambda'_2 + \mu'_2 x) \frac{\partial}{\partial y}, $$

which leads to the system: $\dot{x} = X + u(Y - X)$ connecting the elliptic case to the hyperbolic case. It is an extension of models (16). Straightforward computations lead to an interior collinearity locus, which is defined by a quadratic mapping vs. a linear map in the prey–predator model, while the interior singular locus is defined by a cubic mapping vs. a linear map. More generally automatic computations lead to analyze all the 2d–cases.
### 3.2.2. 3d-case studies

The 3d-case is a very complicated situation due to the complexity of the classification of the singular dynamics. This is not surprising in the context of dynamical systems due to the phenomenon of chaos. For Lotka–Volterra models the situation can be tamed if we consider the case of competitive 3d-models. Restricting to this case, we shall make a brief recap of this theory based on the presentation of [3] introducing the concept of carrying simplicia and the May–Leonard model [15] as a case study.

**Definition 3.6.** If the Lotka–Volterra dynamics takes the form: $\dot{x} = (\text{diag}x)(r - Ax)$, where $a_{ij}, r_i > 0$ for all $i, j = 1, \dots, n$ the system is called competitive.

**Property 3.7.** For a competitive Lotka–Volterra system, we have:

* The equilibrium $x_e = 0$ is an unstable node.

* Every positive trajectory in the $\mathbb{R}_{\ge 0}^n$-space is bounded.

The interesting case, in the non trivial case, is when the dynamics admits an unique equilibrium in the cone $\mathbb{R}_{>0}^n$. Without loss of generality, it can be taken as $\Omega = (1, \dots, 1)$ using Lemma 3.2.

Denoting in short by $X(x)$ the vector field $(\text{diag}x)(r - Ax)$ and by $\varphi_t$ be the (local) parameter group $\{\exp tX\}$, one can assume that $X$ is complete and let $\Lambda^+(x)$ denotes the $\omega$-limit set of a point $x$. Recall that $p \in \Lambda^+(x)$ if there exist a sequence $t_k \to +\infty$ such that $\exp t_k X(x) \to p$ when $k \to +\infty$. One can take [18] as a reference for the properties of such set in relation with stability analysis.

**Definition 3.8.** Since every positive trajectory is bounded in $\mathbb{R}_{\ge 0}^n$ the basin of repulsion of the unstable node $x_e = 0$ is bounded and its boundary is called the *carrying simplex* and is denoted by $\Pi$.

The following result holds [11].

**Theorem 3.9.** In the competitive case, every positive trajectory of the dynamics in $\mathbb{R}_{\ge 0}^n \setminus \{0\}$ is asymptotic to one trajectory in $\Pi$ and $\Pi$ is a Lipschitz submanifold transverse to all strictly positive direction and homeomorphic to the probability simplex: $\Delta_1: \sum_{i=1}^n x_i = 1$.

**Example 3.10.** In many cases, the carrying simplex can be computed and the dynamics on $\Pi$ describes the asymptotic behaviours of positive trajectories. The simplest example where it differs from $\Delta_1$ is given by the 2d-system:

$$ \dot{x} = x(1 - x - y/2), \quad \dot{y} = y(1 - 3x - y). $$

The only equilibria are the origin and the non interior equilibria: (1,0) and (0,1). The carrying simplex is the graph of the function $y(x) = (1-x)^2$ on $[0,1]$ and is a trajectory joining the saddle (0,1) to the attracting node (1,0).

Clearly the computation of this set leads to solve the global stability problem, see [24]. This opens the road to the concept of *controlled stability* of a chosen point in $\Pi$ using the controlled Lotka–Volterra model (see Section 4.2).

### 3.2.3. The May–Leonard model [15]

It is a basic model where the carrying simplex coincides for some parameters of the model with the probability simplex and where the dynamics can be investigated using *bifurcation*
### 3.2.2. 3d-case studies

The 3d-case is a very complicated situation due to the complexity of the classification of the singular dynamics. This is not surprising in the context of dynamical systems due to the phenomenon of chaos. For Lotka–Volterra models the situation can be tamed if we consider the case of competitive 3d-models. Restricting to this case, we shall make a brief recap of this theory based on the presentation of [3] introducing the concept of carrying simplicia and the May–Leonard model [15] as a case study.

**Definition 3.6.** If the Lotka–Volterra dynamics takes the form: $\dot{x} = (\text{diag}x)(r - Ax)$, where $a_{ij}, r_i > 0$ for all $i, j = 1, \dots, n$ the system is called competitive.

**Property 3.7.** For a competitive Lotka–Volterra system, we have:

* The equilibrium $x_e = 0$ is an unstable node.

* Every positive trajectory in the $\mathbb{R}_{\ge 0}^n$-space is bounded.

The interesting case, in the non trivial case, is when the dynamics admits an unique equilibrium in the cone $\mathbb{R}_{>0}^n$. Without loss of generality, it can be taken as $\Omega = (1, \dots, 1)$ using Lemma 3.2.

Denoting in short by $X(x)$ the vector field $(\text{diag}x)(r - Ax)$ and by $\varphi_t$ be the (local) parameter group $\{\exp tX\}$, one can assume that $X$ is complete and let $\Lambda^+(x)$ denotes the $\omega$-limit set of a point $x$. Recall that $p \in \Lambda^+(x)$ if there exist a sequence $t_k \to +\infty$ such that $\exp t_k X(x) \to p$ when $k \to +\infty$. One can take [18] as a reference for the properties of such set in relation with stability analysis.

**Definition 3.8.** Since every positive trajectory is bounded in $\mathbb{R}_{\ge 0}^n$ the basin of repulsion of the unstable node $x_e = 0$ is bounded and its boundary is called the *carrying simplex* and is denoted by $\Pi$.

The following result holds [11].

**Theorem 3.9.** In the competitive case, every positive trajectory of the dynamics in $\mathbb{R}_{\ge 0}^n \setminus \{0\}$ is asymptotic to one trajectory in $\Pi$ and $\Pi$ is a Lipschitz submanifold transverse to all strictly positive direction and homeomorphic to the probability simplex: $\Delta_1: \sum_{i=1}^n x_i = 1$.

**Example 3.10.** In many cases, the carrying simplex can be computed and the dynamics on $\Pi$ describes the asymptotic behaviours of positive trajectories. The simplest example where it differs from $\Delta_1$ is given by the 2d-system:

$$ \dot{x} = x(1 - x - y/2), \quad \dot{y} = y(1 - 3x - y). $$

The only equilibria are the origin and the non interior equilibria: $(1, 0)$ and $(0, 1)$. The carrying simplex is the graph of the function $y(x) = (1 - x)^2$ on $[0, 1]$ and is a trajectory joining the saddle $(0, 1)$ to the attracting node $(1, 0)$.

Clearly the computation of this set leads to solve the global stability problem, see [24]. This opens the road to the concept of *controlled stability* of a chosen point in $\Pi$ using the controlled Lotka–Volterra model (see Section 4.2).

### 3.2.3. The May–Leonard model [15]

It is a basic model where the carrying simplex coincides for some parameters of the model with the probability simplex and where the dynamics can be investigated using *bifurcation*
426 analysis. We follow the presentation of [3], see also [2].

427 **Definition 3.11.** The May–Leonard model has the dynamics

$$
\begin{align}
\dot{x} &= x(1 - x - \alpha y - \beta z) \nonumber \\
\dot{y} &= y(1 - \beta x - y - \alpha z) \tag{21} \\
\dot{z} &= z(1 - \alpha x - \beta y - z) \nonumber
\end{align}
$$

428 and we denote by $A = \begin{pmatrix} 1 & \alpha & \beta \\ \beta & 1 & \alpha \\ \alpha & \beta & 1 \end{pmatrix}$ the associated (circular) matrix.

429 **Property 3.12.** For α, β > 0 and α + β = 2, the carrying simplex Π coincides with the proba-
bility simplex x + y + z = 1 and eliminating in (21) the z-variable leads to the dynamics on
Π:
431

$$
\begin{align}
\dot{x} &= \frac{\alpha - \beta}{2} x(1 - x - 2y) \nonumber \\
\dot{y} &= -\frac{\alpha - \beta}{2} y(1 - 2x - y) \tag{22}
\end{align}
$$

432 with α + β = 2. The system is conservative (i.e. with zero divergence) and we have a canonical
433 Hamiltonian system associated to

$$
H(x, y) = \frac{\alpha - \beta}{2} (1 - x - y) xy
$$

434 so that in the open triangle $\mathcal{T} : \{(x, y) > 0, x + y < 1\}$ all the solutions are periodic.

435 More generally for arbitrary parameters α, β the eigenvalues of A can be determined,
436 see [15] for a discussion.

437 The equilibria points are the origin, three single population of the form (1, 0, 0) and three
438 two populations solutions of the form (1 − α, 1 − β, 0) / (1 − αβ). The interior equilibrium is
439 Ω = (1 , 1 , 1) / (1 + α + β), where the linearized dynamics is determined by the eigenvalues of
440 A.
441

441 Therefore the model is a good case study to apply our analysis in particular taking the
442 control direction Y tangent to the 2d-plane. This allows to get a complete classification of
443 all cases with respect to the parameters (ε₁, ε₂) and computations associated to this exam-
444 ple are given in Section 4.2

4. Applications and numerical results

446 **4.1. Discussion about the nonemptyness property of the good set in Theorem 2.16**

447 **4.1.1. Nonemptyness in the controlled Lotka–Volterra model**

448 Using example 3.5 from Section 3.1, consider the controlled Lotka–Volterra model:
449 −diag( x + 1) [ Ax − u E ] in R 3 with:

450 A = diag(λ₁, λ₁, λ₃),

• $\mathcal{E} = (\epsilon_1, \epsilon_2, \epsilon_3)^T$, $\epsilon_1 \neq \epsilon_2$.
analysis. We follow the presentation of [3], see also [2].

**Definition 3.11.** The May–Leonard model has the dynamics

$$ 
\begin{aligned}
\dot{x} &= x(1 - x - \alpha y - \beta z) \\
\dot{y} &= y(1 - \beta x - y - \alpha z) \\
\dot{z} &= z(1 - \alpha x - \beta y - z)
\end{aligned}
\qquad (21) $$

and we denote by $A = \begin{pmatrix} 1 & \alpha & \beta \\ \beta & 1 & \alpha \\ \alpha & \beta & 1 \end{pmatrix}$ the associated (circular) matrix.

**Property 3.12.** For $\alpha, \beta > 0$ and $\alpha + \beta = 2$, the carrying simplex $\Pi$ coincides with the probability simplex $x + y + z = 1$ and eliminating in (21) the z-variable leads to the dynamics on $\Pi$:

$$ 
\begin{aligned}
\dot{x} &= \frac{\alpha - \beta}{2} x(1 - x - 2y) \\
\dot{y} &= -\frac{\alpha - \beta}{2} y(1 - 2x - y)
\end{aligned}
\qquad (22) $$

with $\alpha + \beta = 2$. The system is conservative (i.e. with zero divergence) and we have a canonical Hamiltonian system associated to

$$ H(x, y) = \frac{\alpha - \beta}{2} (1 - x - y) xy $$

so that in the open triangle $\mathcal{T} : \{(x, y) > 0, x + y < 1\}$ all the solutions are periodic.

More generally for arbitrary parameters $\alpha, \beta$ the eigenvalues of $A$ can be determined, see [15] for a discussion.

The equilibria points are the origin, three single population of the form $(1, 0, 0)$ and three two populations solutions of the form $(1 - \alpha, 1 - \beta, 0)/(1 - \alpha\beta)$. The interior equilibrium is $\Omega = (1, 1, 1)/(1 + \alpha + \beta)$, where the linearized dynamics is determined by the eigenvalues of $A$.

Therefore the model is a good case study to apply our analysis in particular taking the control direction $Y$ tangent to the $2d$-plane. This allows to get a complete classification of all cases with respect to the parameters $(\epsilon_1, \epsilon_2)$ and computations associated to this example are given in Section 4.2

# 4. Applications and numerical results

## 4.1. Discussion about the nonemptiness property of the good set in Theorem 2.16

### 4.1.1. Nonemptiness in the controlled Lotka–Volterra model

Using example 3.5 from Section 3.1, consider the controlled Lotka–Volterra model: $-\mathrm{diag}(x+1) [Ax - u\mathcal{E}]$ in $\mathbb{R}^3$ with:

* $A = \mathrm{diag}(\lambda_1, \lambda_1, \lambda_3)$,

* $\mathcal{E} = (\epsilon_1, \epsilon_2, \epsilon_3)^T$, $\epsilon_1 \neq \epsilon_2$.
Hence since the matrix A contains two identical eigenvalues, one deduces that for every $\mathcal{E}$, the equilibrium point $x_e = 0$ is abnormal with $u_e = 0$ being the abnormal control.

The computation of the determinantal mappings $D, D'$ defined in Section 2.4 yields the singular control $u_s(x) = -D'(x)/D(x) = \lambda_3(x_1 - x_2)/(\varepsilon_1 - \varepsilon_2)$. By construction at $x = 0$, one has $u_s(0) = u_e = 0$. Moreover the spectrum of the linearized singular dynamics at 0 is $\{\lambda_3 - \lambda_1, -\lambda_3, -\lambda_1\}$, proving that the singular point 0 is isolated and that the set introduced in Theorem 2.16 is nonempty.

The calculation can be generalized to the $n$-dimensional Lotka–Volterra case.

### 4.1.2. Emptyness in the quadratic case

**Proposition 4.1.** Let $(Q, b)$ be a pair of vector fields on $\mathbb{R}^3$ with $Q$ being a quadratic homogeneous and $b$ being a constant vector. Then we have:

(1) $D = \det(b, [b, Q], [[b, Q], b]):$ linear form

$D' = \det(b, [b, Q], [[b, Q], Q]):$ cubic form

$D'' = \det(b, [b, Q], Q):$ cubic form.

Hence the singular dynamics

$$ \frac{dx}{dt}(t) = Q(x(t)) - \frac{D'(x(t))}{D(x(t))} b $$

is homogeneous and quadratic.

(2) Let be the time reparameterization defined by $Dd\tau = dt$ so that the singular dynamics reads

$$ \frac{dx}{d\tau}(\tau) = Q(x(\tau))D(x(\tau)) - D'(x(\tau))b $$

and is cubic. It can be projected on the projectivized space $P(\mathbb{R}^3)$. Moreover the map

$$ \hat{\lambda}_s : (Q, b) \mapsto QD - D'b $$

is a semi-covariant.

**Definition 4.2.** Let $\dot{x} = H(x)$ be a differential equation on $\mathbb{R}^n$ with $H$ being a cubic homogeneous vector field. A ray is a line $\ell$ so that the dynamics restricted to $\ell$ is given by $\dot{x}_1 = \lambda x_1^3$. In the generic case, $\lambda \neq 0$, it is an asymptotic direction but if $\lambda = 0$ it is a set of non isolated equilibria for the dynamics.

From which we deduce.

**Corollary 4.3.** For every pair $(Q, b)$ on $\mathbb{R}^3$, every abnormal equilibrium point is not isolated so that in this family the set described in Theorem 2.16 is empty.

Nevertheless, one can use for pair $(Q, b)$ the projectivized singular dynamics on $P(\mathbb{R}^3)$ to compute feedback invariants. Indeed every ray projects onto an equilibrium point, where we can compute the linearized dynamics. We refer to [4] for the application of this result to classify the controlled Euler dynamics in the attitude control problem, where $b$ describes the position of gas jet on the satellite.
Hence since the matrix A contains two identical eigenvalues, one deduces that for every $\mathcal{E}$, the equilibrium point $x_e = 0$ is abnormal with $u_e = 0$ being the abnormal control.

The computation of the determinantal mappings $D, D'$ defined in Section 2.4 yields the singular control $u_s(x) = -D'(x)/D(x) \sim_{x \to 0} \lambda_3(x_1 - x_2)/(\epsilon_1 - \epsilon_2)$. By construction at $x = 0$, one has $u_s(0) = u_e = 0$. Moreover the spectrum of the linearized singular dynamics at 0 is $\{\lambda_3 - \lambda_1, -\lambda_3, -\lambda_1\}$, proving that the singular point 0 is isolated and that the set introduced in Theorem 2.16 is nonempty.

The calculation can be generalized to the $n$-dimensional Lotka–Volterra case.

### 4.1.2. Emptyness in the quadratic case

**Proposition 4.1.** Let $(Q, b)$ be a pair of vector fields on $\mathbb{R}^3$ with $Q$ being a quadratic homogeneous and $b$ being a constant vector. Then we have:

(1) $D = \det(b, [b, Q], [[b, Q], b]):$ linear form

$D' = \det(b, [b, Q], [[b, Q], Q]):$ cubic form

$D'' = \det(b, [b, Q], Q):$ cubic form.

Hence the singular dynamics

$$ \frac{dx}{dt}(t) = Q(x(t)) - \frac{D'(x(t))}{D(x(t))} b $$

is homogeneous and quadratic.

(2) Let the time reparameterization defined by $Dd\tau = dt$ so that the singular dynamics reads

$$ \frac{dx}{d\tau}(\tau) = Q(x(\tau))D(x(\tau)) - D'(x(\tau))b $$

and is cubic. It can be projected on the projectivized space $P(\mathbb{R}^3)$. Moreover the map

$$ \hat{\lambda}_s : (Q, b) \mapsto QD - D'b $$

is a semi-covariant.

**Definition 4.2.** Let $\dot{x} = H(x)$ be a differential equation on $\mathbb{R}^n$ with $H$ being a cubic homogeneous vector field. A ray is a line $\ell$ so that the dynamics restricted to $\ell$ is given by $\dot{x}_1 = \lambda x_1^3$. In the generic case, $\lambda \neq 0$, it is an asymptotic direction but if $\lambda = 0$ it is a set of non isolated equilibria for the dynamics.

From which we deduce.

**Corollary 4.3.** For every pair $(Q, b)$ on $\mathbb{R}^3$, every abnormal equilibrium point is not isolated so that in this family the set described in Theorem 2.16 is empty.

Nevertheless, one can use for pair $(Q, b)$ the projectivized singular dynamics on $P(\mathbb{R}^3)$ to compute feedback invariants. Indeed every ray projects onto an equilibrium point, where we can compute the linearized dynamics. We refer to [4] for the application of this result to classify the controlled Euler dynamics in the attitude control problem, where $b$ describes the position of gas jet on the satellite.
## 4.2. The May-Leonard model

### 4.2.1. Geometric properties

Following Section 3.2.3, we study the asymptotic dynamics of the three dimensional May-Leonard model. This dynamics is characterized by the carrying simplex $\Pi: x + y + z = 1$. We then consider the restriction of (21) to the carrying simplex and the control directions are tangent to the carrying simplex. This procedure can be reproduced for general Lotka-Volterra systems, where we can select a set of species that satisfy a May and Leonard dynamics (see [8]). This leads to the system

$$ \dot{\mathbf{x}} = X(\mathbf{x}) + uY(\mathbf{x}), \quad \mathbf{x} = (x, y), \quad u \in [0, 1] $$

with

$$ X = x(1 - x - 2y)\frac{\partial}{\partial x} - y(1 - 2x - y)\frac{\partial}{\partial y}, \quad Y = \epsilon_1 x\frac{\partial}{\partial x} + \epsilon_2 y\frac{\partial}{\partial y}. \quad (23) $$

The singular set and the collinearity locus are given respectively by

$$ \mathcal{S}: \epsilon_2(\epsilon_1 + 2\epsilon_2)y + \epsilon_1(\epsilon_2 + 2\epsilon_1)x = 0 \quad \text{and} \quad \mathcal{C}: (\epsilon_1 + 2\epsilon_2)y + (\epsilon_2 + 2\epsilon_1)x = \epsilon_1 + \epsilon_2. $$

For $\epsilon_1 \neq \epsilon_2$, $\epsilon_1 \neq -2\epsilon_2$ and $\epsilon_2 \neq -2\epsilon_1$, their intersection $\mathcal{S} \cap \mathcal{C}$ is the point

$$ x_{se} = \left( \frac{\epsilon_2}{\epsilon_2 + 2\epsilon_1}, \frac{\epsilon_1 + \epsilon_2}{\epsilon_2 - \epsilon_1}, \frac{\epsilon_1}{\epsilon_1 + 2\epsilon_2}, \frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2} \right) = \left( \frac{1}{1 + 2\kappa}, \frac{1 + \kappa}{1 - \kappa}, \frac{\kappa}{\kappa + 2}, \frac{\kappa + 1}{\kappa - 1} \right), \quad \kappa = \frac{\epsilon_1}{\epsilon_2} $$

and the study of the components of $x_{se}$ as functions of $\kappa$ shows that $x_{se}$ is in the positive orthant if $\kappa \in (-\infty, -2[\cup]-1/2, 0]$ and is outside the triangle $\mathcal{T} = \{\mathbf{x} = (x, y) > 0, x+y<1\}$ for any value of $\kappa$ (recall that in the triangle $\mathcal{T}$ the free dynamics consists in periodic orbits).

For $x_e = (x_{1e}, x_{2e}) \in \mathcal{C}$, the value of the control $u_e(x_e)$ such that $X(x_e) + u_e Y(x_e) = 0$ is

$$ u_e(x_e) = \frac{1 - 3x_{1e}}{\epsilon_1 + 2\epsilon_2} $$

and the spectrum of the Jacobian matrix $J(x_e) := \frac{\partial}{\partial \mathbf{x}}(X(\mathbf{x}) + uY(\mathbf{x}))|_{\mathbf{x}=x_e}$ is

$$ \sigma(J(x_e)) = \left\{ 3u(\epsilon_1 + \epsilon_2) \pm \frac{1}{2\sqrt{3}} \sqrt{11u^2(\epsilon_1 + \epsilon_2)^2 + 4(\epsilon_1 u + 1)(\epsilon_2 u - 1)} \right\}. $$

In particular,

• at $u = 0$, $\sigma(J(x_e)) = \{\pm \frac{i}{\sqrt{3}}\}$,

• at $u = 1$, $\sigma(J(x_e)) = \{3(\epsilon_1 + \epsilon_2) \pm \frac{1}{2\sqrt{3}} \sqrt{11(\epsilon_1 + \epsilon_2)^2 + 4(\epsilon_1 + 1)(\epsilon_2 - 1)}\}$,

• at $u = u_e(x_{se})$, $\sigma(J(x_e)) = \left\{ \frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2}, -\frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2} \frac{3\epsilon_1\epsilon_2}{(2\epsilon_1 + \epsilon_2)(2\epsilon_2 + \epsilon_1)} \right\}$. Note that in this case, the first eigenvalue corresponds to the noncontrollable mode of the system, while the second eigenvalue is not meaningful since it can be replaced by any value via pole shifting.

The singular dynamics $\dot{z} = \tilde{H}_s(z)$ (restricted to $\Sigma$) has eigenvalues at $x = x_{se}$:

$$ \left\{ -\frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2}, \frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2} \right\} $$
## 4.2. The May-Leonard model

### 4.2.1. Geometric properties

Following Section 3.2.3, we study the asymptotic dynamics of the three dimensional May-Leonard model. This dynamics is characterized by the carrying simplex $\Pi: x + y + z = 1$. We then consider the restriction of (21) to the carrying simplex and the control directions are tangent to the carrying simplex. This procedure can be reproduced for general Lotka-Volterra systems, where we can select a set of species that satisfy a May and Leonard dynamics (see [8]). This leads to the system

$$ \dot{\mathbf{x}} = X(\mathbf{x}) + u Y(\mathbf{x}), \quad \mathbf{x} = (x, y), \quad u \in [0, 1] $$

with

$$ X = x(1 - x - 2y)\frac{\partial}{\partial x} - y(1 - 2x - y)\frac{\partial}{\partial y}, \quad Y = \epsilon_1 x\frac{\partial}{\partial x} + \epsilon_2 y\frac{\partial}{\partial y}. \quad (23) $$

The singular set and the collinearity locus are given respectively by

$$ \mathcal{S}: \epsilon_2(\epsilon_1 + 2\epsilon_2)y + \epsilon_1(\epsilon_2 + 2\epsilon_1)x = 0 \quad \text{and} \quad \mathcal{C}: (\epsilon_1 + 2\epsilon_2)y + (\epsilon_2 + 2\epsilon_1)x = \epsilon_1 + \epsilon_2. $$

For $\epsilon_1 \neq \epsilon_2$, $\epsilon_1 \neq -2\epsilon_2$ and $\epsilon_2 \neq -2\epsilon_1$, their intersection $\mathcal{S} \cap \mathcal{C}$ is the point

$$ x_{se} = \left( \frac{\epsilon_2}{\epsilon_2 + 2\epsilon_1}, \frac{\epsilon_1 + \epsilon_2}{\epsilon_2 - \epsilon_1}, \frac{\epsilon_1}{\epsilon_1 + 2\epsilon_2}, \frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2} \right) = \left( \frac{1}{1 + 2\kappa}, \frac{1 + \kappa}{1 - \kappa}, \frac{\kappa}{\kappa + 2}, \frac{\kappa + 1}{\kappa - 1} \right), \quad \kappa = \frac{\epsilon_1}{\epsilon_2} $$

and the study of the components of $x_{se}$ as functions of $\kappa$ shows that $x_{se}$ is in the positive orthant if $\kappa \in (-\infty, -2[\cup]-1/2, 0]$ and is outside the triangle $\mathcal{T} = \{\mathbf{x} = (x, y) > 0, x + y < 1\}$ for any value of $\kappa$ (recall that in the triangle $\mathcal{T}$ the free dynamics consists in periodic orbits).

For $x_e = (x_{1e}, x_{2e}) \in \mathcal{C}$, the value of the control $u_e(x_e)$ such that $X(x_e) + u_e Y(x_e) = 0$ is

$$ u_e(x_e) = \frac{1 - 3x_{1e}}{\epsilon_1 + 2\epsilon_2} $$

and the spectrum of the Jacobian matrix $J(x_e) := \frac{\partial}{\partial \mathbf{x}}(X(\mathbf{x}) + uY(\mathbf{x}))|_{\mathbf{x}=x_e}$ is

$$ \sigma(J(x_e)) = \left\{ 3u(\epsilon_1 + \epsilon_2) \pm \frac{1}{2\sqrt{3}} \sqrt{11u^2(\epsilon_1 + \epsilon_2)^2 + 4(\epsilon_1 u + 1)(\epsilon_2 u - 1)} \right\}. $$

In particular,

* at $u = 0$, $\sigma(J(x_e)) = \{\pm \frac{i}{\sqrt{3}}\}$,

* at $u = 1$, $\sigma(J(x_e)) = \{3(\epsilon_1 + \epsilon_2) \pm \frac{1}{2\sqrt{3}} \sqrt{11(\epsilon_1 + \epsilon_2)^2 + 4(\epsilon_1 + 1)(\epsilon_2 - 1)}\}$,

* at $u = u_e(x_{se})$, $\sigma(J(x_e)) = \left\{ \frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2}, \frac{-\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2} \frac{3\epsilon_1\epsilon_2}{(2\epsilon_1 + \epsilon_2)(2\epsilon_2 + \epsilon_1)} \right\}$. Note that in this case, the first eigenvalue corresponds to the noncontrollable mode of the system, while the second eigenvalue is not meaningful since it can be replaced by any value via pole shifting.

The singular dynamics $\dot{z} = \dot{H}_s(z)$ (restricted to $\Sigma$) has eigenvalues at $x = x_{se}$:

$$ \left\{ -\frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2}, \frac{\epsilon_1 + \epsilon_2}{\epsilon_1 - \epsilon_2} \right\} $$
and $\lambda = \frac{c_1+c_2}{c_1-c_2}$ is a feedback invariant of the control system.

### 4.2.2. Direct and semi-direct methods

We present two numerical schemes to solve the time-minimal control problem for the Lotka-Volterra class of models specifically a direct method and a semi-direct method.

The objective is to reach, from an initial position $x_0$ a terminal manifold of codimension one in minimum time, namely, this target is taken as a disk $N(\mathbf{x}) \le 0$ centered on the collinearity locus $\mathcal{C}$ in which local controllability is guaranteed provided that the forcing feedback is interior. This is an example of controlled stability.

There is no use comparing these methods in terms of computational time. The direct method computes open-loop controls, while the semi-direct method computes closed-loop controls, which can be used for real-time application since the trajectory is computed step by step. In this sense we use the direct method to have an upper bound of the value function and it is compared to the objective value associated to the semi-direct method presented below.

Both methods are implemented in the WOLFRAM LANGUAGE using the Brent's principal axis method of the Mathematica's routine FindMinimum.

**Direct method.** We consider the time-optimal control problem:

$$ (OCP) \qquad \begin{array}{ll} \min_{u(\cdot),T} & T \\ & \dot{\mathbf{x}}(t) = X(\mathbf{x}(t)) + u(t) Y(\mathbf{x}(t)), \quad u(t) \in [0,1], \text{ a.e. } t \in [0,T] \\ & \mathbf{x}(0) = x_0 \quad (\text{given}) \\ & N(\mathbf{x}(T)) \le 0 \end{array} $$

where $\mathbf{x} = (x, y)$ and $X, Y$ are given by (23).

We perform a discretization over the state and the control spaces for (OCP) to obtain a nonlinear finite dimensional optimization problem. The optimization variables are the values of the control at each time step and a primal-dual interior point method is used to solve numerically the optimization problem. The optimality conditions – written as a relaxation of the Karush-Kuhn-Tucker conditions – are determined using automatic differentiation.

It is usually a quite robust method with respect to the initialization compared to indirect methods based on the Pontryagin Maximum Principle. However it does not exploit the geometric structure of the optimal control. The BOCOP software<fnref n="2" /> provides an implementation of this method based on the IPOPT optimization solver<fnref n="3" />.

**Model predictive control method.** The semi-direct method is based on a sampled-data control formulation of the problem, which is adapted to medical protocols that can be used to cure the C. difficile infection. The optimal path is constructed iteratively, where at each iteration we solve an optimal control problem on a reduced time horizon, whereas the direct method discretizes the problem on the whole time interval.

In this sense, our method is closely related to model predictive control (MPC) widely used for control theory applications [17, 22]. The current state of the method $x_c$ is initial-

<footnote n="2">²www.bocop.org</footnote>

<footnote n="3">³www.coin-or.github.io/Ipopt</footnote>
ized to $x_0$ and is updated iteratively by solving the optimal control problems of the form

$$ (OCP') \quad \begin{array}{ll} \min_{u \in \mathbb{R}^h} & N(\mathbf{x}(t_h; u, x_c)) \\ & \dot{\mathbf{x}}(t) = X(\mathbf{x}(t)) + u_i Y(\mathbf{x}(t)), \quad u_i \in [0, 1], \text{ a.e. } t \in [t_i, t_{i+1}], \quad i = 0, \dots, h-1, \\ & \mathbf{x}(0) = x_c \end{array} $$

where the integer $h$ is the horizon, $0 = t_0 < \dots < t_h$ are given fixed times and $\mathbf{x}(.; u, x_c)$ is the state response associated to $(u_1, \dots, u_h)$ and starting at $x_c$ at $t=0$. The algorithm terminates when $|N(x_c)|$ is smaller than a given threshold.

To solve (OCP') numerically, we derive a finite dimensional optimization problem by constructing an approximation of the objective function $u = (u_0, \dots, u_{h-1}) \mapsto N(\mathbf{x}(t_h))$ via an approximation of $\mathbf{x}(t_h; u, x_c)$ by discretizing the differential constraint with a midpoint rule (the discretization of the state on $[0, t_h]$ is finer than the partition $0 < t_1 < \dots < t_h$). The approximation of $N(\mathbf{x}(t_h; u, x_c))$ – together with its derivatives with respect to $u_i, i = 1, \dots, h$ – can be computed offline using symbolic computations. Then we solve the optimization problem associated to (OCP') using a primal-dual interior point method. Once (OCP') is solved for the current value of $x_c$, we retrieve the values of $t_1$ and $u_1$ to update $x_c \leftarrow \mathbf{x}(t_1; u_1, x_c)$ and we iterate considering the resulting new instance of (OCP').

### 4.2.3. Numerical results

We apply the previous numerical methods to the specific May–Leonard model and our aim is to control the system in the carrying simplex $\Pi$.

From the computations in Section 4.2.1, we can choose $\epsilon_1, \epsilon_2$ to fix the positions of $x_{se}$ and $x_e$ (corresponding to $u_e = 1$) and so that the singular dynamics goes toward $\mathcal{C}$ and is hyperbolic. The target is taken as the circle $N$ centered on $\mathcal{C}$ and is reachable with an admissible control $u \in ]0, 1[$.

We fix the initial point to $x_0 = (0.1, 0.8)$. The direct method converges in about 200 iterations to a bang–bang–singular–bang control (see Fig. 1) and the corresponding trajectory reaches $N$ in less than 7.7 unit of time.

The semi-direct method is tuned with an horizon of $t_h = 2$ unit of time with $h = 3$ i.e. we compute three controls $u_1, u_2, u_3$ over this horizon. The resulting control seems to have the same structure as the control of the direct method. The trajectory reaches $N$ in about 12 unit of time (see Fig. 1). The depicted singular behavior depends on the size of the horizon $t_h$ and the number of controls on that horizon. A smaller horizon $t_h$ would typically keep the trajectory away from the singular arc.

## 5. Conclusion

In this article we describe feedback invariants to classify single–input control affine systems in relation with the time–minimal control problem. They can be explicitly calculated using algebraic computations in the jets spaces of the dynamics at a given point corresponding to an abnormal stationary geodesic. They correspond at such points to the spectrum of the projectivized Hamiltonian geodesic dynamics.

This result completes the computation of the feedback invariants related to the concept of conjugate points for both normal and abnormal geodesics. This gives a neat common geometric frame since those are related to the spectrum of the "projectivized" second–order intrinsic derivative, in relation with algebraic computations in the jets space of geodesics.
ized to $x_0$ and is updated iteratively by solving the optimal control problems of the form

$$ (OCP') \quad \begin{array}{ll} \min_{u \in \mathbb{R}^h} & N(\mathbf{x}(t_h; u, x_c)) \\ & \dot{\mathbf{x}}(t) = X(\mathbf{x}(t)) + u_i Y(\mathbf{x}(t)), \quad u_i \in [0, 1], \text{ a.e. } t \in [t_i, t_{i+1}], \quad i = 0, \dots, h-1, \\ & \mathbf{x}(0) = x_c \end{array} $$

where the integer $h$ is the horizon, $0 = t_0 < \dots < t_h$ are given fixed times and $\mathbf{x}(\cdot; u, x_c)$ is the state response associated to $(u_1, \dots, u_h)$ and starting at $x_c$ at $t=0$. The algorithm terminates when $|N(x_c)|$ is smaller than a given threshold.

To solve (OCP') numerically, we derive a finite dimensional optimization problem by constructing an approximation of the objective function $u = (u_0, \dots, u_{h-1}) \mapsto N(\mathbf{x}(t_h))$ via an approximation of $\mathbf{x}(t_h; u, x_c)$ by discretizing the differential constraint with a midpoint rule (the discretization of the state on $[0, t_h]$ is finer than the partition $0 < t_1 < \dots < t_h$). The approximation of $N(\mathbf{x}(t_h; u, x_c))$ – together with its derivatives with respect to $u_i, i = 1, \dots, h$ – can be computed offline using symbolic computations. Then we solve the optimization problem associated to (OCP') using a primal-dual interior point method. Once (OCP') is solved for the current value of $x_c$, we retrieve the values of $t_1$ and $u_1$ to update $x_c \leftarrow \mathbf{x}(t_1; u_1, x_c)$ and we iterate considering the resulting new instance of (OCP').

### 4.2.3. Numerical results

We implement the previous numerical methods to the specific May–Leonard model and our aim is to control the system in the carrying simplex $\Pi$.

More precisely, we choose, according to the computations in Section 4.2.1, $\epsilon_1, \epsilon_2$ such that

* $x_{se}$ and $x_e$ (corresponding to $u_e = 1$) belong to the nonnegative orthant,

* $x_{se}$ is a hyperbolic stable fixed point for the singular dynamics.

From the initial point $x_0 = (0.1, 0.8)$, we seek the minimum time to reach the circle $N$, centered at $x_c = (1.8, 0.058) \in \mathcal{C}$, with a control valued in $[0, 1]$.

The results of the direct method and model predictive control method, described in the previous section, are represented in Fig. 1. The direct method converges in about 200 iterations to a bang–bang–singular–bang control and the corresponding trajectory reaches $N$ in less than 7.7 unit of time.

The semi-direct method associated to the optimization problem (OCP') is tuned with an horizon of $t_h = 2$ unit of time with $h = 3$. The resulting control have the same structure as the control of the direct method that is bang–bang–singular–bang. The trajectory reaches $N$ in about 12 unit of time. The depicted singular behavior depends on the size of the horizon $t_h$ and the number of controls on that horizon. A smaller horizon $t_h$ would typically keep the trajectory away from the singular arc.

## 5. Conclusion

In this article we describe feedback invariants to classify single–input control affine systems in relation with the time–minimal control problem. They can be explicitly calculated using algebraic computations in the jets spaces of the dynamics at a given point corresponding to an abnormal stationary geodesic. They correspond at such points to the spectrum of the projectivized Hamiltonian geodesic dynamics.

This result completes the computation of the feedback invariants related to the concept
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Numerical simulation performed on the May–Leonard model (21) with $\alpha = \beta = 2$, $\varepsilon_1 = 1$, $\varepsilon_2 = -2.9$. The aim is to reach in minimum time the ball of radius $\sqrt{0.03}$ centered at $x_C = (1.8, 0.058)$ located on the collinearity locus. (left) The dashed trajectory starting from $x_0 = (0.1, 0.8)$ and obtained with a direct method is bang–bang–singular–bang. It reaches the terminal circle $N$ in less than 7.7 unit of time. The continuous MPC trajectory, obtained with an horizon of $t_h = 2$ and with three controls ($h = 3$), seems to reproduce the singular behavior. It reaches $N$ in about 12 unit of time. (right) Time evolution of the control for the direct and MPC methods.</figcaption></figure>

Our result is briefly applied to the controlled Lotka–Volterra, where abnormal stationary geodesics are shifted equilibria of the free dynamics, which can be determined using linear calculations only. This gives an algebraic frame in relation with population dynamics control. Another applications are for quadratic systems in relation for instance with the attitude control problem of a rigid spacecraft.

A final numerical application shows the relation of the study to the time–minimal control problem, combining Pontryagin Maximum Principle with direct and semi–direct numerical frame. Additionally it provides an example of controlled stability method popular in biological models. More generally our analysis leads to structurally stable results based on Lie algebraic computations in the jets space of the geodesics dynamics. It opens the road to fine applications translating for instance results for the free dynamics of the May–Leonard model [2] to controlled stability [9] or time optimal results, in the controlled case.

Preliminary results in this direction are obtained in [6, 8] to control the infection of a complex microbiote with three species using probiotic and antibiotic treatment. The analysis of section 4.2 leads to analyze the effect of the therapy in the final phase of the treatment (see [7]). Additional applications concern the controlled Euler equation where the classification of the singular dynamics initialized in [4] is crucial in the attitude control problem.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Numerical simulation performed on the May–Leonard model (21) with $\alpha = \beta = 2$, $\varepsilon_1 = 1$, $\varepsilon_2 = -2.9$. The aim is to reach in minimum time the ball of radius $\sqrt{0.03}$ centered at $x_C = (1.8, 0.058)$ located on the collinearity locus. (left) The dashed trajectory starting from $x_0 = (0.1, 0.8)$ and obtained with a direct method is bang–bang–singular–bang. It reaches the terminal circle $N$ in less than 7.7 unit of time. The continuous MPC trajectory, obtained with an horizon of $t_h = 2$ and with three controls ($h = 3$), seems to reproduce the singular behavior. It reaches $N$ in about 12 unit of time. (right) Time evolution of the control for the direct and MPC methods.</figcaption></figure>

of conjugate points for both normal and abnormal geodesics. This gives a neat common geometric frame since those are related to the spectrum of the "projectivized" second-order intrinsic derivative, in relation with algebraic computations in the jets space of geodesics.

Our result is briefly applied to the controlled Lotka–Volterra, where abnormal stationary geodesics are shifted equilibria of the free dynamics, which can be determined using linear calculations only. This gives an algebraic frame in relation with population dynamics control. Another applications are for quadratic systems in relation for instance with the attitude control problem of a rigid spacecraft.

A final numerical application shows the relation of the study to the time–minimal control problem, combining Pontryagin Maximum Principle with direct and semi–direct numerical frame. Additionally it provides an example of controlled stability method popular in biological models. More generally our analysis leads to structurally stable results based on Lie algebraic computations in the jets space of the geodesics dynamics. It opens the road to fine applications translating for instance results for the free dynamics of the May–Leonard model [2] to controlled stability [9] or time optimal results, in the controlled case.

Preliminary results in this direction are obtained in [6, 8] to control the infection of a complex microbiote with three species using probiotic and antibiotic treatment. The analysis of section 4.2 leads to analyze the effect of the therapy in the final phase of the treat-
600 **Funding**

601 This work benefited from the support of the program PEPS "Jeunes chercheurs et jeunes
602 chercheuses" of Insmi 2022.

603 **References**

604 [1] A. AGRACHEV, I. ZELENKO, *On feedback classification of control-affine systems with one and two-dimensional inputs*. SIAM J. Control Optim., **46**, no.4 (2007), pp.1431–1460.

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

609 [3] S. BAIGENT, *Lotka–Volterra Dynamical Systems*. In: S. Bullett, T. Fearn, F. Smith, eds. Dynamical and Complex Systems, LTCC Advanced Mathematics Series: Chapter **5**. Singapore: World Scientific, 2017, 227 pages.

612 [4] B. BONNARD, *Feedback equivalence for nonlinear systems and the time optimal control problem*. SIAM J. on Control and Optim., **29** (1991), pp. 1300–1321.

614 [5] B. BONNARD, I. KUPKA, *Generic properties of singular trajectories*. Ann. Inst. H. Poincaré Anal. Non Linéaire, **14**, no.2 (1997), pp. 167–186.

616 [6] B. BONNARD, J. ROUOT, *Optimal Control of the Controlled Lotka-Volterra Equations with Applications - The Permanent Case*. SIAM J. Appl. Dyn., **22**, no. 2 (2023), pp. 2761–2791.

619 [7] B. BONNARD, J. ROUOT, *Optimal Control of the Lotka–Volterra Equations with Applications*. Ivan Kupka’s legacy. A tour through controlled dynamics, in AIMS Applied math. **12** (2024), pp. 15–33.

622 [8] B. BONNARD, J. ROUOT, C. J. SILVA, *Geometric optimal control of the generalized Lotka–Volterra model of the intestinal microbiome*. Optimal Control Applications and Methods, **45**, no. 2 (2024), pp. 544–574.

625 [9] F.H. CLARKE, Yu. S. LEDYAEV, L. RIFFORD, R.J. STERN, *Feedback stabilization and Lyapunov functions*. SIAM J. Control Optim., **39**, no.1 (2000), pp.25–48.

627 [10] J.A. DIEUDONNÉ, J.B. CARRELL, *Invariant Theory, Old and New*. Academic Press, New York, 1971, 85 pages.

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

631 [12] E.W. JONES, P.S. CLARKE, 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.

634 [13] 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.

636 [14] J. MARTINET, *Singularities of smooth functions and maps*. London Mathematical Society Lecture Note Series, **58**. Cambridge University Press, Cambridge-New York, 1982, 256 pages.

639 [15] R.M. MAY, W.J. LEONARD, *Nonlinear aspects of competition between three species*. SIAM J. Appl. Math., **29** (1975), pp. 243–253.

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

643 [17] J. B. RAWLINGS, D. Q. MAYNE, M. M. DIEHL, *Model Predictive Control: Theory, Computation, and Design*. Nob Hill Publishing, LLC, 2020, 770 pages.

644
ment (see [7]). Additional applications concern the controlled Euler equation where the classification of the singular dynamics initialized in [4] is crucial in the attitude control problem.

## Funding

This work benefited from the support of the program PEPS "Jeunes chercheurs et jeunes chercheuses" of Insmi 2022.

## References

[1] A. AGRACHEV, I. ZELENKO, *On feedback classification of control-affine systems with one and two-dimensional inputs*. SIAM J. Control Optim., **46**, no.4 (2007), pp.1431–1460.

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

[3] S. BAIGENT, *Lotka–Volterra Dynamical Systems*. In: S. Bullett, T. Fearn, F. Smith, eds. Dynamical and Complex Systems, LTCC Advanced Mathematics Series: Chapter **5**. Singapore: World Scientific, 2017, 227 pages.

[4] B. BONNARD, *Feedback equivalence for nonlinear systems and the time optimal control problem*. SIAM J. on Control and Optim., **29** (1991), pp. 1300–1321.

[5] B. BONNARD, I. KUPKA, *Generic properties of singular trajectories*. Ann. Inst. H. Poincaré Anal. Non Linéaire, **14**, no.2 (1997), pp. 167–186.

[6] B. BONNARD, J. ROUOT, *Optimal Control of the Controlled Lotka–Volterra Equations with Applications – The Permanent Case*. SIAM J. Appl. Dyn., **22**, no. 2 (2023), pp. 2761–2791.

[7] B. BONNARD, J. ROUOT, *Optimal Control of the Lotka–Volterra Equations with Applications*. Ivan Kupka’s legacy. A tour through controlled dynamics, in AIMS Applied math. **12** (2024), pp. 15–33.

[8] B. BONNARD, J. ROUOT, C. J. SILVA, *Geometric optimal control of the generalized Lotka–Volterra model of the intestinal microbiome*. Optimal Control Applications and Methods, **45**, no. 2 (2024), pp. 544–574.

[9] F.H. CLARKE, Yu. S. LEDYAEV, L. RIFFORD, R.J. STERN, *Feedback stabilization and Lyapunov functions*. SIAM J. Control Optim., **39**, no.1 (2000), pp.25–48.

[10] J.A. DIEUDONNÉ, J.B. CARRELL, *Invariant Theory, Old and New*. Academic Press, New York, 1971, 85 pages.

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

[12] E.W. JONES, P.S. CLARKE, 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.

[13] 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.

[14] J. MARTINET, *Singularities of smooth functions and maps*. London Mathematical Society Lecture Note Series, **58**. Cambridge University Press, Cambridge-New York, 1982, 256 pages.

[15] R.M. MAY, W.J. LEONARD, *Nonlinear aspects of competition between three species*. SIAM J. Appl. Math., **29** (1975), pp. 243–253.
645

646

[18] N. ROUCHE, J. MAWHIN, *Equations Différentielles Ordinaires*. Masson, Paris, **2** (1973), 266 pages.

647

648

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

649

650

651

[20] R.R. STEIN, V. BUCCI, N.C. TOUSSAINT, C.G. BUFFIE, G. RÄTSCH, E.G. PAMER, et al., *Ecological modelling from time-series inference: insight into dynamics and stability of intestinal microbiota*. PLoS Comp. Biology, **9** no. 12 (2013).

652

653

[21] V. VOLTERRA, *Leçons sur la théorie mathématique de la lutte pour la vie*. Les Grands Classiques Gauthier-Villars. Éditions Jacques Gabay, Sceaux, 1990, 215 pages.

654

655

[22] Y. WANG, S. BOYD, *Fast model predictive control using online optimization*, IEEE Transactions on control systems technology, **18** no. 2, (2010), 267–278.

656

657

[23] W.M. WONHAM, *Linear Multivariable Control*. Springer-Verlag, New York, 1985, 348 pages.

658

659

[24] E.C. ZEEMAN, M.L. ZEEMAN, *From local to global behavior in competitive Lotka-Volterra systems*. Trans. Amer. Math. Soc., 355 (2003), pp. 713–734.
[16] L.S. PONTRYAGIN, V.G. BOLTYANSKII, R.V. GAMKRELIDZE, E.F. MISHCHENKO, *The math-ematical theory of optimal processes*. Oxford, Pergamon Press, 1964, 362 pages.

[17] J. B. RAWLINGS, D. Q. MAYNE, M. M. DIEHL, *Model Predictive Control: Theory, Com-putation, and Design*. Nob Hill Publishing, LLC, 2020, 770 pages.

[18] N. ROUCHE, J. MAWHIN, *Equations Différentielles Ordinaires*. Masson, Paris, **2** (1973), 266 pages.

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

[20] R.R. STEIN, V. BUCCI, N.C. TOUSSAINT, C.G. BUFFIE, G. RÄTSCH, E.G. PAMER, et al., *Ecological modelling from time-series inference: insight into dynamics and stability of intestinal microbiota*. PLoS Comp. Biology, **9** no. 12 (2013).

[21] V. VOLTERRA, *Leçons sur la théorie mathématique de la lutte pour la vie*. Les Grands Classiques Gauthier-Villars. Éditions Jacques Gabay, Sceaux, 1990, 215 pages.

[22] Y. WANG, S. BOYD, *Fast model predictive control using online optimization*, IEEE Transactions on control systems technology, **18** no. 2, (2010), 267–278.

[23] W.M. WONHAM, *Linear Multivariable Control*. Springer-Verlag, New York, 1985, 348 pages.

[24] E.C. ZEEMAN, M.L. ZEEMAN, *From local to global behavior in competitive Lotka-Volterra systems*. Trans. Amer. Math. Soc., 355 (2003), pp. 713–734.