# A novel approach to feedback control with deep reinforcement learning

Kala Agbo Bidi, Jean-Michel Coron, Amaury Hayat, Nathan Lichtlé

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

## To cite this version:

Kala Agbo Bidi, Jean-Michel Coron, Amaury Hayat, Nathan Lichtlé. A novel approach to feedback control with deep reinforcement learning. Systems and Control Letters, 2025, 202, pp.106102. 10.1016/j.sysconle.2025.106102. hal-05339948

HAL Id: hal-05339948

https://hal.science/hal-05339948v1

Submitted on 31 Oct 2025

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

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

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

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

# A novel approach to feedback control with deep reinforcement learning

Kala Agbo dici<sup>a ID</sup>, Jean-Michel Coron<sup>a ID</sup>, Amaury Hayat<sup>b ID</sup>,*, Nathan Lichtlé<sup>b,c ID</sup>

## ARTICLE INFO

**Keywords:**

Asymptotic stability

Global stabilization

Reinforcement learning

Robust control

Sterile Insect Technique

## ABSTRACT

We present a novel approach to feedback control design by leveraging the power of deep reinforcement learning (RL). The goal is to blend the RL methodology with mathematical analysis to extract an explicit feedback control for systems where, because of constraints or limited measurements, the classical approaches of control theory (e.g. dynamical programming, optimal control-based feedback, backstepping, etc.) cannot be used.

We study a dynamical system of mosquito populations for biological pest control using the Sterile Insect Technique (SIT), a method traditionally applied in agriculture that involves releasing large numbers of sterile insects into the wild to reduce pest populations. Our goal is to derive a feedback control that globally stabilizes the system around the zero-mosquito equilibrium using only practical measurements, such as total male and female mosquito counts, rather than detailed counts of sterilized versus potent males or fecund versus unfecund females, which are often not accessible. This physical constraint presents challenges for classical methods control theory, as the full state cannot be measured. To address this, we apply deep reinforcement learning to suggest feedback laws for a discretized system that only rely on these accessible, real-world measurements, obtainable through methods like pheromone traps. Finally, we leverage the trained neural network to extract explicit feedback controls that stabilize the original continuous system over a wide range of initial conditions.

Many other dynamical systems arising from practical applications are subject to measurement constraints, which render the stabilization problem complex from a mathematical perspective. We believe that this approach could help in finding new solutions to these problems.

## 1. Introduction

Stabilization is a branch of mathematics that is very useful for many applications. Many controls have been developed in the last century to address the local or global stabilization of dynamic systems [1]. One can cite for instance optimal control-based feedback, dynamical programming [2,3], backstepping [4–6], and many others. However, in some cases, such methods cannot be used, for instance because the system suffers from physical constraints, such as limited measurements. Observers strategy to tackle this limitation have been developed but in many cases finding robust observers is still a challenge. In this paper, we introduce an new approach that employs deep reinforcement learning (RL) to suggest mathematical control strategies for dynamical systems. We illustrate it on a mathematical model of the Sterile Insect Technique (SIT) applied on mosquitoes population, namely (2.1)–(2.4) below. Controlling this system is particularly challenging because we only rely on infrequent, discrete measurements of the physically observable part of the state.

Our approach differs from classical methods and uses deep reinforcement learning (RL) to construct control feedback laws. Over the past few years, RL [7,8] has emerged as a powerful approach for control, with its ability to learn near-optimal decision-making strategies through interactions with an environment, and has demonstrated remarkable successes across a wide range of domains and applications with long-term horizons, high-dimensional partially-observable states. In particular, RL can perform well on certain tasks exhibiting states that are too large, too complex to model or too partially-observable to be efficiently solved with traditional methods such as optimal control, the cost being that RL controllers typically do not inherently have safety guarantees, explainability or optimality [8,9]. In robotics, RL has enabled machines to learn complex tasks such as locomotion, manipulation, and dexterous object handling [10–12]. In the realm of games, RL algorithms have achieved superhuman performance in challenging domains like Go, Chess, StarCraft [13,14], as well as in classical Atari games [15]. These remarkable achievements highlight the versatility and potential of RL as a general-purpose approach for solving complex

<footnote><sup>a</sup> Sorbonne Université, CNRS, Université Paris-Cité, Laboratoire Jacques-Louis Lions, LJLL, INRIA équipe CAGE, Paris, 75005, France</footnote>

<footnote><sup>b</sup> CERMICS, Ecole Nationale des Ponts et Chaussées, 6-8 avenue Blaise Pascal, Champs-sur-Marne, France</footnote>

<footnote><sup>c</sup> EECS, University of California, Berkeley, Berkeley, CA, USA</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: amaury.hayat@enpc.fr (A. Hayat).</footnote>
practical control problems in diverse domains. Akin to robotics, RL has been used to solve systems of ODEs or PDEs [16], which can sometimes be augmented by enforcing constraints on the training of the RL algorithm [17] or directly on the neural networks, for instance with physically-informed neural networks (PINNs) [18]. [19] does RL through symbolic regression using a basis of functions, which has the advantage of always being explicit and thus more interpretable, but in general is likely less expressive than state-of-the-art RL algorithms coupled with neural networks. Other approaches to symbolic regression using RL and neural-guided genetic programming can be found in [20], although it considers learning explicit forms of simple mathematical formulas, which does not consider the mathematical system solving part.

RL techniques, while powerful for decision-making, inherently provide control mechanisms that are discrete and numerical in nature. However, from a more rigorous mathematical point of view, these mechanisms often do not translate directly into analytical feedback control formulas. Recognizing this limitation, in our work, we blend RL methodologies with mathematical analysis to extract an explicit mathematical control. Importantly, in our approach we apply RL to solve the mathematical problem, rather than solely applying RL to the discretized system, a more conventional application of RL. This allows us to employ deep RL to architect control feedback laws which can be juxtaposed with existing controls that have been derived from more traditional methodologies.

Our example of application is the Sterile Insect Technique (SIT): it is a method of biological pest control that consists of releasing sterilized insects to reduce or eliminate a target population. Initially used in agriculture to control insect pests, it is now employed in the vector-born disease fight against mosquitoes that carry illnesses such as malaria and arboviruses [21,22] and there is a great interest both in research and in practice to understand which control to use for releasing the sterilized insects [23–25].

Several mathematical approaches have already been used in the literature to treat the SIT control problem applied to mosquitoes either for the complete system or for reduced models. Two reduced models have been considered: a two dimensional (2D) model obtained by assuming that the dynamics of males and eggs are fast so that these two populations can be assumed to be at equilibrium (see [25], $S_1$, page 231-232) or [26], (2)); and a three dimensional (3D) model obtained by overlooking the non-adult stages (see [27], (7a)-(7b)-(7c)). These mathematical approaches have led to the following stabilizing feedback controls:

* Stabilization using impulsive feedback controls for the 3D model: [27], Theorem 6] and [27, Theorem 7] for the case of sparse measurements. The case of vector migration is also considered in [28].

* Stabilization using optimal feedback control of the reduced model 2D model and the complete model: [25] and in particular [25, Remark 3.4].

* Stabilization using the backstepping method: [26] considers the 2D model while [29, Section 3.1] considers the complete model. See, for example, [30, Section I.2.2], [31, Pages 242–246] or [1, Section 12.5] for tutorial presentations of the backstepping method.

* Stabilization using properties of monotone systems [28,32].

* Stabilization using simple linear feedback laws: [29] proves stabilization for positively invariant subsets and conjectures it for the complete model.

These feedback controls above are constructed using classical tools in control theory such as control Lyapunov functions, the LaSalle invariance principle, the maximum principle, monotone dynamical systems, barrier functions, or the backstepping method. Although these feedback controls provide evidence in terms of robustness, their applications requires continuously measuring separately the different components

of the state of the model, which is difficult or impossible in practice. Thank to our approach, we are able to provide a control only based on the most accessible data to be measured: the adult population (adult females and adult males) using pheromone traps, which are already used in practice.

## 2. Mathematical framework of the problem

The SIT model in mosquitoes population is given by the following system of equations:

$$ \dot{E} = \beta_E F \left(1 - \frac{E}{K}\right) - (v_E + \delta_E)E, \quad (2.1) $$

$$ \dot{M} = (1 - v) v_E E - \delta_M M, \quad (2.2) $$

$$ \dot{F} = v v_E E \frac{M}{M + M_s} - \delta_F F, \quad (2.3) $$

$$ \dot{M}_s = u - \delta_s M_s, \quad (2.4) $$

where, at time $t$, $E(t) \ge 0$ is the mosquito density in aquatic phase, $M(t) \ge 0$ is the wild adult male density, $F(t) \ge 0$ is the density of adult females which have been fertilized, $M_s(t) \ge 0$ is the sterilized adult male density, and $u(t) \ge 0$, the control, is the density of sterilized males released at time $t$.

In system (2.1)–(2.4) we assume that all females mate as soon as they emerge from the pupal stage. The density of unfertilized females, i.e. the density of females that have mated with sterilized males, is denoted by $F_s(t)$. One has $F_s(t) = F(t)M_s(t)/M(t)$. Besides, we also assume that $\delta_s \ge \delta_M$, which is usually considered to be a biologically relevant assumption [25]. The interpretation of the parameters are given below [25]:

* $\beta_E > 0$ is the oviposition rate,

* $\delta_E, \delta_M, \delta_F > 0$ are the death rates for eggs, wild adult males and fertilized females respectively,

* $v_E > 0$ is the hatching rate for eggs,

* $v \in (0, 1)$ the probability that a pupa gives rise to a female (and $(1 - v)$ is, therefore, the probability to give rise to a male),

* $\delta_s > 0$ is the death rate of sterilized adults,

* $K > 0$ is the environmental capacity for eggs. It can be interpreted as the maximum density of eggs that females can lay in breeding sites. Since here the larval and pupal compartments are not present, it can be interpreted as $E$ representing all the aquatic compartments and this term $K$ representing a logistic law's carrying capacity for the aquatic phase (that also includes the effects of competition between larvae).

Typical values for these parameters as well as the values used in this work are given in Table A.5.

For the parameters given in Table A.5, when $u(t) = M_s(t) = 0$ for any $t \ge 0$, the system (2.1)–(2.3) has a unique globally asymptotically stable equilibrium ($E(t), M(t), F(t) \equiv (E^*, M^*, F^*)$ where $E^*, M^*$ and $F^*$ are large constant values. This corresponds to the situation where mosquitoes reproduce freely. The state ($E(t), M(t), F(t) \equiv (0, 0, 0)$ is also an equilibrium, albeit an unstable one. The mathematical problem is to find $u(t)$ of the form

$$ u(t) = f(M(t) + M_s(t), F(t) + F_s(t)), \quad (2.5) $$

where $f \in L^\infty(\mathbb{R}^2)$, such that the zero equilibrium $(0, 0, 0)$ is globally asymptotically stable and $M_s$ is asymptotically small, meaning there exists $c \in \mathbb{R}_+$ such that

$$ \lim_{t \to +\infty} \|u(t)\| = c < U^* \quad (2.6) $$

where

$$ U^* := \frac{K \beta_E v (1 - v) v_E^2 \delta_s}{4(\delta_E + v_E) \delta_F \delta_M} \left( 1 - \frac{\delta_F (v_E + \delta_E)}{\beta_E v v_E} \right)^2, \quad (2.7) $$
and the equilibrium (0, 0, 0, c/δs) of the system (2.1)–(2.4) is globally asymptotically stable (see Definition 2.1 below, where the notion of solutions of the closed-loop system is understood in the Filippov sense [29, Section 2.2]). Ideally, one would even like to be able to find, for any ε > 0, a control feedback law fε such that

$$ \lim_{t \to +\infty} u(t) = \varepsilon. \qquad (2.8) $$

**Definition 2.1.** The equilibrium (0, 0, 0, c/δs) of the system (2.1)–(2.4) with the feedback law (2.5) is globally asymptotically stable if, for any initial condition (E₀, M₀, F₀, Ms,₀), the (forward maximal) solutions (E, M, F, Ms) to the system (2.1)–(2.4) with the feedback law (2.5) are defined on [0, +∞) and for any ε > 0 there exists δ > 0 such that

$$ \begin{gathered} \|(E_0, M_0, F_0, M_{s,0} - c/\delta_s)\| \le \delta \\ \|(E(t), M(t), F(t), M_s(t) - c/\delta_s)\| \le \varepsilon, \quad \forall t \in [0, +\infty), \end{gathered} \qquad (2.9) $$

$$ \lim_{t \to +\infty} \|(E(t), M(t), F(t), M_s(t) - c/\delta_s)\| = 0, \qquad (2.10) $$

The form constraint (2.5) corresponds to a practical limitation due to the quantities that we are able to reliably measure in practice: M + Ms and F + Fs are the total number of males and females which are typically what can be measured in practice (see [21]), and in particular f cannot depend on E or M and Ms separately since these two quantities are difficult to measure. Nevertheless, we consider different variants in this work.

**Remark 2.1 (Constant Control).** The rationale behind the definition of (2.6) is that this is the critical value above which a constant control can stabilize the state (E*, M*, F*) = (0, 0, 0). Indeed, for a constant control u(t) ≡ ¯U, if ¯U > U* then the equilibrium (E*, M*, F*, Ms*) = (0, 0, 0, U*/δs) is globally asymptotically stable (see [25]).

**Remark 2.2 (Optimal Decay Rate).** Assume that E(0) ≤ K. Then, for any control u(t) ≥ 0 is, for every time t ≥ 0, we have E(t) ≤ K and

$$ E(t) \ge \bar{E}(t), \quad M(t) \ge \bar{M}(t), \quad F(t) \ge \bar{F}(t), \qquad (2.11) $$

where (Ẽ, M̃, F̃) is the solution to the Cauchy problem

$$ \dot{\bar{E}} = \beta_E \bar{F} \left(1 - \frac{\bar{E}}{K}\right) - (v_E + \delta_E) \bar{E}, \qquad (2.12) $$

$$ \dot{\bar{M}} = (1 - v) v_E \bar{E} - \delta_M \bar{M}, \qquad (2.13) $$

$$ \dot{\bar{F}} = -\delta_F \bar{F}, \qquad (2.14) $$

$$ (\bar{E}(0), \bar{M}(0), \bar{F}(0)) = (E(0), M(0), F(0)). \qquad (2.15) $$

It would be interesting to see if one can get with suitable output feedback laws (vanishing or small at the origin) a decay rate close to the one imposed by (2.11), i.e.

$$ E(t) \approx \bar{E}(t), \quad M(t) \approx \bar{M}(t), \quad F(t) \approx \bar{F}(t). \qquad (2.16) $$

(Note that it is possible to get (2.16) by taking u constant and large, depending on (E(0), M(0), F(0)).) This would be particularly useful in the case where the insect under study reproduces both sexually and asexually: indeed, this would give the best way to reduce the sexual reproduction part as much as possible by output feedback laws (vanishing or small at the origin).

In [29] a backstepping feedback control was built to stabilize this specific system at the origin. However, this control depends on the three variables (E, M, Ms) and not only on the feasibly observable quantities M + Ms and F + Fs. As of now, there is no known control depending only on M + Ms and F + Fs.

One possible approach could be to use the measurements of M + Ms and F + Fs to design an observer that reconstructs the rest of the state, enabling the application of backstepping control as an observer-based control. This has been successfully implemented when separate

measurements on M and Ms are available [33]. However, it remains a challenge when only M + Ms and F + Fs are measured. The main difficulty arises from the term M/(M + Ms), which is not globally Lipschitz and becomes singular at (0, 0). This can pose a problem when M + Ms is small and M is not directly measured. Consequently, certain common observer methods, such as extended Kalman filters for nonlinear systems and linearization techniques, are already inapplicable. To the best of our knowledge, no effective approach exists to address this challenge.

### 3. Approach and contributions of this paper

In this work we use a deep Reinforcement Learning (RL) approach to construct control feedback laws and compare them with the existing feedback controls that were deduced in previous work.

Our approach works in three steps:

1. We discretize the equations in a numerical scheme,

2. We train an RL model to obtain a numerical control feedback based on this numerical scheme,

3. We recover from the numerical control feedback an explicit mathematical control. We then perform several tests to ensure that the explicit control is efficient.

This is detailed in Section 4. We use this approach to construct three types of control feedback laws:

* A feedback control depending on M only,

* A feedback control depending on M + Ms (total number of males) only,

* A feedback control depending on M + Ms and F + Fs (total number of males and of females).

These different types of controls are studied in Sections 5, 6 and 7 respectively.

### 4. Reinforcement learning based control

Reinforcement Learning (RL) is a subfield of machine learning concerned with training agents to make decisions in an environment to maximize a long-term reward function. From a control perspective this can equivalently be seen as finding an optimal control for a cost function over the trajectories. We first define these concepts more formally and discuss the RL algorithm we use in Section 4.1, then explain how we apply this formalism to our specific problem in Section 4.2. After that, we discuss our approach in Section 4.3 and finally go through experiment details in Section 4.4.

#### 4.1. RL background

We model the environment using the common formalism of a Partially-Observable Markov Decision Process (POMDP) [34] $M = (S, A, T, R, \gamma, \mu, \Omega, \mathcal{O})$ where $S \subseteq \mathbb{R}^n$ is a set of states, $A \subseteq \mathbb{R}^m$ a set of actions, $T: S \times A \to \Delta S$ is the state transition function (ie. $T(s'|s, a)$ is the probability of transitioning to state $s'$ given state $s$ and action $a$), $R: S \times A \to \mathbb{R}$ is the reward function, $\gamma \in [0, 1)$ is the discount factor, $\mu \in \Delta(S)$ is the initial state distribution, $\Omega \subseteq \mathbb{R}^p$ is a set of observations of the hidden state, and $\mathcal{O}: S \to \Delta(\Omega)$ is the observation distribution (ie. $\mathcal{O}(o|s)$ is the probability of getting observation $o$ given current state $s$). Note that given a set $X$, $\Delta(X)$ denotes the set of probability distributions over $X$.

The goal for the agent is to learn a policy $\pi_\theta: \Omega \to \Delta(\mathcal{A})$ (stochastic in our case) mapping observations to actions, where $\theta$ are the parameters of the policy (typically the weights of a neural network
in the case of deep RL), which maximizes the expected discounted sum of rewards

$$J(\pi_{\theta}) = \mathbb{E}_{\tau \sim (\pi_{\theta}, M)} \left[ \sum_{t=0}^{\infty} \gamma^t r_t \right] \quad (4.1)$$

where the expectation is taken over all trajectories $\tau = (s_t, a_t, r_t)_{t \ge 0}$ generated by the current policy $\pi_{\theta}$ acting in the POMDP $M$. In other words, we aim to find optimal parameters $\theta^* = \arg\max_{\theta} J(\pi_{\theta})$. Note that maximizing this objective $J(\pi_{\theta})$ is analogous to the standard minimizing of the cost function $-J(\pi_{\theta})$ in control theory.

One common way to minimize this cost is to use policy gradient methods, which directly optimize the policy parameters by estimating the gradient of the expected cumulative reward. The basic policy gradient algorithm updates the policy parameters in the direction of the estimated gradient to increase the likelihood of actions that lead to higher rewards:

$$\theta_{\text{new}} \leftarrow \theta_{\text{old}} + \nabla_{\theta_{\text{old}}} \hat{J}(\pi_{\theta_{\text{old}}}), \quad (4.2)$$

where $\hat{J}$ is the sample mean of $J$, calculated by averaging the discounted sum of rewards over a number of trajectories. To improve stability and convergence, several techniques have been developed, such as Trust Region Policy Optimization (TRPO) [35] and Proximal Policy Optimization (PPO) [36] which we use in this work. TRPO limits the policy update step size by constraining the divergence between the new policy and the old policy. PPO introduces a simple surrogate objective function that includes a clipping mechanism to prevent large policy updates and improve sample efficiency, while still maximizing the likelihood of high-reward actions, and vice versa. In parallel, PPO also learns a value function, which is a second neural network that aims to estimate the value of each state, i.e. the total amount of rewards that can be expected when starting from a given state. A second term in the loss function drives the learning of this value function, by minimizing its mean squared error with respect to the empirical values of the state obtained from policy trajectories. Finally, a third term encourages exploration by rewarding higher-entropy policies.

## 4.2. Defining the POMDP

In this section, we define the specific POMDP problem formulation that we consider in this work: state, observations, actions and rewards, as well as the initial state distribution and state transition function that model the system. We present all the variants that we have considered, and the specific settings used for each controller will be specified in the respective feedback law parts in Sections 5–7.

**State space** The state for the POMDP are exactly the state of the SIT model, that is, $(E(t), M(t), F(t), M_s(t))^T$, all components being nonnegative and introduced in Section 2. Thus our state space can be formally written as $S = \mathbb{R}_{\ge 0}^4$.

**Observation space and distribution** To account for real-world partial observability constraints, the control does not have access to the full state, but only to some partial observation of it. In this work, we consider three different types of observations which we analyze and compare in their respective sections:

* Section 5 considers an observation consisting of only the number of wild males. Formally, the observation space is $\Omega = \mathbb{R}_{\ge 0}$ and the function mapping the state to observations (which in this work is deterministic) is $\mathcal{O}(E(t), M(t), F(t), M_s(t)) = M(t)$.

* Section 6 considers an observation consisting of the total number of males. Formally, $\Omega = \mathbb{R}_{\ge 0}$ and $\mathcal{O}(E(t), M(t), F(t), M_s(t)) = M(t) + M_s(t)$.

* Section 7 considers an observation consisting of the total number of males and the total number of females. Formally, $\Omega = \mathbb{R}_{\ge 0}^2$ and $\mathcal{O}(E(t), M(t), F(t), M_s(t)) = (M(t) + M_s(t), F(t) + \hat{F}_s(t))$, where $\hat{F}_s(t) = F(t)M_s(t)/M(t)$ was defined in Section 2.

All of the observations that we consider are quantities that we are able to measure in the real world. This is in contrast to the backstepping control of [29], which requires measurements of $E$, $M$ and $M_s$, which is currently not possible in the real world. Besides, these observations may grow very large and thus can span quite a large range. As such, to enable the neural network to observe large values while still being able to discriminate between smaller values, we input each observation into the control at different orders of magnitude. For instance, $M + M_s$ is inputted at several scales, namely $M + M_s$ becomes $\min(M + M_s, k)$ for $k$ ranging from 5 to 100K, and similarly for $F + F_s$. We normalize all observations so that they lie within [0, 1]. This has proved important to help with convergence during training: normalizing inputs is a common preprocessing technique, and inputting each observation at several scales further helps with training stability without technically adding more inputs to the control. We also consider adding memory of the past observations as an input to the control. The rationale behind is that we hope that with this additional information the neural network will be able to reconstruct in some sense the state and, in other words, internally build a kind of observer that would remain a black box. This approach however, makes it significantly more complex to convert into an explicit feedback.

**Action space** Our action space $\mathcal{A} = [-1, 1]$ corresponds to a single action $a(t) \in \mathcal{A}$ that is remapped to the range [0, 10K] and then directly inputted into the model Eqs. (2.1)–(2.4) through $u(t) = 5K(a(t) + 1)$. Having the neural network model output a normalized action is a common technique for more robust training, akin to normalization of the inputs (here we output a normalized action, then scale it up).

**Initial state distribution** We sample the initial state uniformly between 0 and 10K, which corresponds to physically-realistic values for the components of the state. Namely, $(E(0), M(0), F(0), M_s(0)) \sim \mu = \mathcal{U}([0, 10K]^4)$.

**State transition function** The state transition function $T$, that maps a current state and action to a next state, is deterministic in this work and implicitly defined by the ODE system (2.1)–(2.4): we discretize the state $x_t = (E(t), M(t), F(t), M_s(t))$, and at each time step compute the next state through a simple Euler update $x_{t+1} = x_t + \dot{x}_t dt$ (with an abuse of notation, since the update is technically done on each of the four individual components of the state), where the action $u_t = u(t)$ comes in the definition of $\dot{M}_s$. The value for $dt$ is indicated in Section 4.4. Thus, the transition function can be written as $T(x_t, u_t) = x_{t+1}$.

**Reward function** The optimization criterion is usually the most crucial part of the RL process. Our reward function, which we aim to maximize over time as per Eq. (4.1), takes the following form at time step $t$:

$$r_t = -c_1 (\|E(t)\|_2 + \|M(t)\|_2 + \|F(t)\|_2) - c_2(t) \|M_s(t)\|_2 \quad (4.3)$$

with

$$c_2(t) = \begin{cases} c_3 & \text{if } t < 0.9T, \\ c_3 + c_4 & \text{otherwise,} \end{cases} \quad (4.4)$$
where $T$ is the simulation horizon. In a continuous control framework this would correspond to maximizing the functional

$$J(u) = \int_{0}^{T} c_{1} \|E(t), M(t), F(t)\|_{2} + c_{2}(t) \|M_{s}(t)\|_{2} dt. \quad (4.5)$$

The specific values we use for $c_1$, $c_3$ and $c_4$ are indicated in Section 4.4, and are derived empirically and heuristically through hyperparameter tuning. The reason for the specific shape of this reward function Eq. (4.3) is as follows: ideally, we would simply penalize all the components of the state equally and have

$$r_t = -C (\|E(t)\|_2 + \|M(t)\|_2 + \|F(t)\|_2 + \|M_s(t)\|_2)$$

for some $C > 0$ in order to drive all the components of the state to zero. This is what we initially tried; however, for practical reasons, we found this objective to be harder to minimize due to the fact that $M_s$ can typically take much larger values than the 3 other component of the state. Indeed, one of the difficulties of this control problem is that $M_s$ should necessarily take high values to be able to bring $M$, $F$ and $E$ closer to 0, given (2.1)–(2.4), and this results in a delay between the control action and its effect on the mosquito population. A practical interpretation is that the females have no preferences between sterile males and fertile males, hence $M_s$ should be larger than $M$ to have an influence.

As a result, our choice of $c_2(t)$ corresponds to penalizing $M_s$ with less amplitude than we penalize the other 3 components of the state; nevertheless around the end of the simulation we increase the penalty on $M_s$ to encourage convergence to 0. This allows for a high action at the start to grow $M_s$ without excessive penalty, which in turn makes the other components of the state decrease, then a slowly decreasing action so that $M_s$ converges to 0. The main challenge is to not decrease $M_s$ too quickly, or the other components of the state would increase again. We empirically found that this reward design led to increased training stability: it is designed in a way to guide the controller's, which is initially random, to a reasonable behavior more quickly. We note that the rewards are also normalized by $K$ to lie within a reasonable range.

**Remark 4.1.** We note that the reward function during training has access to the full simulation state and may depend on components outside the control input to optimize performance. However, at inference, there is no more reward and the control relies solely on the provided inputs. Actually, one of the key challenges of this control problem is precisely that the reward to be maximized during training can depend on the full state, whereas the control itself cannot measure these states. This distinction highlights an important difficulty: the controller must learn to optimize a function it does not fully observe. In fact, the method would be less powerful if we only optimized cost functions that are entirely observable. Thus the RL control has no direct way of fully understanding the reward but still manages to optimize it successfully.

Additionally, in order to artificially reduce the horizon and make training more robust (RL usually suffers from overly long simulations, as it makes optimization a reward over time much more complex), we repeat each action several times, meaning that for each environment step we use the same action to run $n_{sims}$ simulation steps. Finally, the value for $\gamma$ is indicated in Section 4.4.

### 4.3. Method

In our work, we aim to train control policies using RL to regulate a dynamic system described by ordinary differential equations. To do that, we design an environment that simulates the behavior of the ODE system, allowing us to interact with it in a controlled manner. Our

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Diagram representing the (simplified) procedure by which we simulate our model in an environment that is used to train an RL agent, whose policy we then convert into an explicit control. The policy $\pi_\theta$ is modeled by a neural network with parameters $\theta$, which takes a state as an input and outputs an action (or a distribution over actions in the stochastic case). The neural network is then optimized to maximize the sum of rewards it obtains over simulations.</figcaption></figure>

approach is summarized in Fig. 1. We start by discretizing our system of Eqs. (2.1)–(2.4), and use those dynamics to create a simulation of our model. We use it to create an environment by implementing the observations, actions and rewards described in Section 4.2. Using this code, we train an RL agent that learns to maximize the objective function we assign it through many simulations, using the PPO algorithm mentioned in Section 4.1. Once the policy has converged, we can evaluate it on any simulation, in particular we can query a control $u$ for any current state $(E, M, F, M_s)^T$. Since we only trained policies with 1 or 2 observations, we can plot the action as a function of the input in 1D or 2D space. This allows us to perform a regression and empirically write a simple explicit control that has the same general shape as the neural network control (see Sections 5 and 7). Finally, we simulate this explicit control to ensure that it still stabilizes the system, and analyze its properties and its sensitivity with respect to the parameters of the system.

### 4.4. Experiment details

Here we present the experiment details in greater details. To implement the PPO algorithm we use to train our RL policies, we use Stable Baselines 3 [37] (version 1.6.2 in Python 3.8), a popular RL library that provides a collection of state-of-the-art algorithm implementations, as well as various tools for RL research. The models in Sections 5–7 are trained for 10 million environment timesteps (or 7 billion simulation timesteps) on 12 CPUs, which takes about 7 h. During each iteration, we collect 12288 (1024 per CPU) environment steps, then run 5 epochs of optimization with a batch size of 1024. The agent's policy is a fully-connected neural network with 2 hidden layers of 256 neurons each, with tanh non-linearities between each layer, outputting the mean and standard deviation of a normal distribution that is then used to sample the action. More formally, for a given observation vector $o_t$, the neural network policy outputs $(\mu_t, \sigma_t) = \pi_\theta(o_t)$ and the action is sampled as $a_t \sim N(\mu_t, \sigma_t)$. We train with a learning rate of $3 \times 10^{-4}$, discount coefficient of $\gamma = 0.99$, and all other hyperparameters are left to their default values.
We run each simulation for $T = 1000$ days, with a timestep $dt = 0.01$ days, and each action is repeated $n_{sims} = 500$ times, meaning that the environment horizon is 200 steps. For each simulation, the initial condition is uniformly sampled between 0 and $10K: E(0), M(0), F(0), M_s(0) \sim U(0, 10K)$. For our reward function, we use coefficients $c_1 = 0.1, c_3 = 0.001$ and $c_4 = 0.01$. For the model, we use the parameters given in Table A.5.

## 5. Feedback RL control using M

In the literature, one of the previous approaches consist in using a linear feedback control that only depends on $M(t)$ [29], that is $u(t) = f(M)$ where $f$ is linear instead of the control law (2.5). They observe that the following linear control feedback law

$$u(t) = \alpha M(t), \quad (5.1)$$

seems to stabilize the system as long as

$$\alpha > \frac{(\beta_E v v_E - (v_E + \delta_E) \delta_F) \delta_s}{(v_E + \delta_E) \delta_F}. \quad (5.2)$$

In this section, we use our RL procedure to deduce a potentially nonlinear feedback law that similarly only depends on $M(t)$.

We force the feedback control to go to 0 when the number of wild male goes to 0 by searching the control under the form

$$u(t) = \min \left( \min \left( f_{RL} \left( \frac{M(t)}{K} \right), \alpha_M \right) M(t), u_M \right), \quad (5.3)$$

where $\alpha_M$ is a chosen constant, $u_M$ is the maximal value of the control allowed which is dictated by physical constraints, and $f_{RL}$ is the function searched by the RL model. This $f_{RL}$ is searched using the procedure described in Section 4 and using the cost function

$$\begin{aligned} J(u) = & - \int_{0}^{T} r_{t} dt \\ r_{t} = & - \frac{\|E(t), M(t), F(t)\|_{2}}{K} \\ & + q_{1} \left( \frac{M_{s}(t)}{K} + \max \left( 0, \left( \frac{M_{s}(t)}{K} - 30 \right) \right)^{2} \right) \end{aligned} \quad (5.4)$$

where $q_1$ is a chosen constant, typically much smaller than 1, and $T$ is a chosen horizon. After training, the RL model converges and the optimal numerical $f_{RL}$ obtained has a relatively simple form, shown in Fig. 2, which happens to be exactly piecewise linear. This allows us to deduce the following nonlinear control feedback law for the system

$$u(t) = \min(\bar{u}(t), u_M) \quad (5.5)$$

with

$$\bar{u}(t) = \min \left[ \max \left( \alpha_1 - \alpha_2 \frac{M(t)}{K}, 0 \right), \alpha_M \right] M(t), \quad (5.6)$$

where $\alpha_1$ and $\alpha_2$ are positive constants. For the values in Table A.5 with $\alpha_M = 15, \alpha_1 \approx 16$ and $\alpha_2 \approx 12$.

**Remark 5.1.** The RL control (5.5) tends to be linear when $K \to +\infty$.

In Fig. 3 we present 100 numerical simulations of the closed-loop system when using this explicit control. Each simulation has a different initial condition taken uniformly at random with each component of the state having values in $[0, 10K]$. We can see that all the components of the state of the system converge quickly to 0 after 800 days, while the main components of the system (not sterile males $M$, females $F$ and eggs $E$) converge much faster to 0.

In Table 1 we show the average and variance of the components of the state $E, M, F$ and $M_s$ at different times over 10000 numerical simulations as well as the maximal absolute value of each component of the state.

<figure><img src="image_6.png" /><figcaption>Fig. 2. Value of the control $f_{RL}$ as a function of $M(t)/K$ when $\alpha_M = 15$ and parameters as in Table A.5.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 3. Representation of the state in a 100 simulations, with $\alpha_M = 15, q_1 = 0.004, T = 2000$ days, $u_M = 10K$, and the parameters of Table A.5. Each simulation corresponds to one color.</figcaption></figure>

<table><caption>**Table 1**
Average, variance and maximum of the different components over 100 simulations, with $\alpha_M > \alpha_1 = 13, q_1 = 0.004, u_M = 10K$ and the parameters of Table A.5.</caption><thead><tr><th></th><th>200 days</th><th>400 days</th><th>800 days</th><th>2000 days</th></tr></thead><tbody><tr><td>average |E| + |M| + |F|</td><td>8.7 10<sup>4</sup></td><td>3.3 10<sup>4</sup></td><td>1.3 10<sup>3</sup></td><td>1.5 10<sup>-2</sup></td></tr><tr><td>variance |E| + |M| + |F|</td><td>3.6 10<sup>7</sup></td><td>1.4 10<sup>7</sup></td><td>3.7 10<sup>4</sup></td><td>5.5 10<sup>-6</sup></td></tr><tr><td>average |M<sub>s</sub>|</td><td>1.6 10<sup>6</sup></td><td>1.5 10<sup>6</sup></td><td>8.4 10<sup>4</sup></td><td>1.0</td></tr><tr><td>variance |M<sub>s</sub>|</td><td>8.2 10<sup>9</sup></td><td>1.3 10<sup>10</sup></td><td>1.6 10<sup>8</sup></td><td>2.5 10<sup>-3</sup></td></tr><tr><td>maximum |E| + |M| + |F|</td><td>9.3 10<sup>4</sup></td><td>3.8 10<sup>4</sup></td><td>1.5 10<sup>3</sup></td><td>1.8 10<sup>-2</sup></td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Fig. 4. Evolution of the sterilized male population as a function of time when applying the feedback law $u(t) = \beta(M(t) + M_s(t))$ with $\beta = 0.118$ (top) or with $\beta = 0.116$ (bottom), over 3000 days.</figcaption></figure>

## 6. RL control using $M + M_s$

In practice, measuring $M$ is a challenge. In a wild population of mosquitoes, synthetic versions of female insect pheromones are released to attract and capture male insects. This allows for measuring $M + M_s$, however there is no easy way to distinguish between the wild male mosquitoes and the released sterilized male. Another possible measurement can be done by placing simple traps for adult mosquitoes in the wild, then differentiating them based on physical attributes such as their size to separately count total males $M + M_s$ and total females $F + F_s$.

This motivates the search, in practice, for a control that would only depend on $M + M_s$ and $F + F_s$. In this section, we consider using only total males $M + M_s$, since this quantity is more easily measured and we compare to the previous section where $M$ is available instead. Interestingly, the situation is much more complicated when $M + M_s$ is available instead of $M$, as we are presenting in this section.

In [29], it was conjectured that a linearly dependent feedback of $M + M_s$ stabilizes the dynamics at the origin. However this control lacks robustness (see Fig. 4) with respect to the parameters of the model or of the controller. Indeed the control has the linear form

$$u(t) = \beta(M(t) + M_s(t)), \qquad (6.1)$$

and can only work if $\beta$ satisfy

$$\left( \frac{\beta_E v v_E - (v_E + \delta_E) \delta_F}{\beta_E v v_E} \right) \delta_s \le \beta < \delta_s. \qquad (6.2)$$

With the parameters set in the Table A.5, this condition becomes

$$0.118 \le \beta < 0.12, \qquad (6.3)$$

which means that even a tiny imprecision on the model parameters renders the control inefficient.

We tested our RL control procedure in this framework, having only access to $M + M_s$ at the current time $t$. Even when allowed a nonlinear control, the model does not converge to an efficient control after 1000 iterations of training, with the same setup used to train the other controls. Fig. 5 illustrates the behavior of the trained control, which outputs an approximately constant control $u(t) \approx 200000$ yielding a stabilization of the populations of mosquitoes, except for the released sterilized male mosquitoes $M_s$ which the control does not manage to reduce with the limited information it has about the state of the system. This demonstrates the lack of robustness mentioned earlier.

Nevertheless, when allowing the control to depend not only on $M + M_s$ at the given time $t$ but also on previous times $s \le t$, that is to say enabling the control with memory of past values of $M + M_s$, the RL policy converges to what seems to be a robust numerical control. The rationale behind giving the model access to a memory of the measurements is

<figure><img src="image_9.png" /><figcaption>Fig. 5. Evolution of the population of mosquitoes, represented as $E(t) + M(t) + E(t)$ and $M_s(t)$, and of the control $u(t)$, as a function of time, over 10 simulations with random initial conditions, using an RL control that only has access to $M + M_s(t)$ at the current time $t$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Evolution of the population of mosquitoes, represented as $E(t) + M(t) + E(t)$ and $M_s(t)$, and of the control $u(t)$, as a function of time, over 100 simulations with random initial conditions, using an RL control that only has access to measurements of $M + M_s$ over the past 6 months, and every week obtains a new measurement and take a new action.</figcaption></figure>

to allow it to derive additional information on the state by internally deriving a kind of observer (that would remain a black box). We train the model by allowing it to measure the state of the system every 7 days and to keep the memory of the 26 most recent measurements (so 6 months of measurements). This measurement frequency corresponds to what is possible in practice [38]. Besides, the control also takes a single action every 7 days. Since this control has strictly more information than the control with only $M + M_s$ at time $t$, we expect it to perform better as long as the training procedure is stable. It is however not obvious that adding this memory of past observations would be enough to stabilize the system, but Fig. 6 shows numerically that it appears to be sufficient, across 100 numerical simulations of the closed-loop system using this numerical control with initial conditions chosen at random in $[0, 5K]^4$.
<figure><img src="image_11.png" /><figcaption>Fig. 7. Heatmap of the model's action $u(M + M_s, F + F_s)$ as a function of total males and total females, in linear scale (left) and logarithmic scale (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 8. Heatmap of the regression model's action $u_{reg}(M + M_s, F + F_s)$ as a function of total males and total females. A state-space trajectory is plotted in red, with the dot indicating initial state and the cross final state, comparing the no-noise case (left) with the case when a small noise $\mu \sim N(0,5)$ is added on top of the action (right).</figcaption></figure>

## 7. RL control using $M + M_s$ and $F + F_s$

Using past values of the state can be a challenge to find a mathematical formula from the numerical control. Indeed, when the control depends on the past state, the feedback that is searched is not anymore a function of a finite-dimensional vector but a functional on an infinite dimensional space containing portions of the trajectories (e.g. ($M + M_s(s)_{s \in [t-T, t]}$). This makes the symbolic regression a challenge. For this reason we try to find a control using $M + M_s$ and $F + F_s$ only at the current time, that is a control of the form (2.5). These two quantities can be measured in real life, and lead to a much simpler model than the model with memory that has 26 inputs.

A first control. In this framework, the RL model is trained as described in Section 4 and converges to a numerical control that we represent in Fig. 7 as a function of $M + M_s$ and $F + F_s$. We see that the plot of the control in linear scale is not really informative (see Fig. 7 left). However in log scale the expression of the control seems clearer (see Fig. 7 right) and clearly has two parts. In each of them the control seems to be close to a bang-bang control with a thin transition. With a simple regression we approximate this numerical control with

$$u_{\text{reg}}(M + M_s, F + F_s) = \begin{cases} u_{\text{reg}}^{\text{left}} & \text{if } M + M_s < 200, \\ u_{\text{reg}}^{\text{right}} & \text{otherwise,} \end{cases} \quad (7.1)$$

for $M + M_s, F + F_s > 0$, where

$$u_{\text{reg}}^{\text{left}} = \begin{cases} u_{\text{min}} & \text{if } \log \frac{200}{F+F_s} > 4, \\ u_{\text{max}} \left(4 - \log \frac{200}{F+F_s}\right) & \text{if } 4 \ge \log \frac{200}{F+F_s} > 3, \\ u_{\text{max}} & \text{otherwise, and} \end{cases} \quad (7.2)$$

$$u_{\text{reg}}^{\text{right}} = \begin{cases} u_{\text{min}} & \text{if } \log \frac{M+M_s}{F+F_s} > 4, \\ u_{\text{max}} \left(4 - \log \frac{M+M_s}{F+F_s}\right) & \text{if } 4 \ge \log \frac{M+M_s}{F+F_s} > 3, \\ u_{\text{max}} & \text{otherwise.} \end{cases} \quad (7.3)$$

Table 2 shows that this explicit control is still able to quickly stabilize the state over a wide range of initial conditions.

<figure><img src="image_13.png" /><figcaption>Fig. 9. Norm of the components of the state $\|E(t), M(t), F(t)\|_2$ (blue) and control $v_{\text{reg}}$ (red) as a function of time for different values of $u_{\text{min}}$ (0, 0.001, 1, and 5 respectively from top to bottom) and $u_{\text{max}} = 300000$, over 2000 days and with the same initial condition. When the minimum control $u_{\text{min}}$ is 0, the control occasionally (more and more rarely) outputs the maximum action (corresponding to the spikes) to prevent the norm of the state from going back up. However, when $u_{\text{min}} > 0$, after some time this minimal action is sufficient to stabilize the state around our target (i.e. $(E^*, F^*, M^*) = (0, 0, 0)$) and no more spikes are observed.</figcaption></figure>

During training and testing, the numerical control feedback law includes by default a small (gaussian) noise. This is made to promote some robustness of the control during the training and a good exploration. We tested the mathematical control we derived (given in (7.1)) with and without noise. To our surprise, the control with a small noise does seem to ensure the asymptotic stability, whereas the control without any noise does not seem to. Indeed, without noise, the control seems to have a cyclic behavior and never converges (see Fig. 8 (left)). When adding a small noise, however, the stability is restored (see Fig. 8 (right)). The explication to this apparent paradox is that having exactly $u_{\text{min}} = 0$ in one of the branches of the control given in (7.1) is seemingly too strong to allow the model to converge completely. Replacing this value with $u_{\text{min}} > 0$ for a small $u_{\text{min}}$ (typically $u_{\text{min}} = 5$) allows to stabilize the system without noise. In the system with noise, because the control $u$ has to be positive, the noise has the effect of increasing in average the effective value of $u_{\text{min}}$ of the control (7.1), which explains the apparent stabilization. Of course when setting $u_{\text{min}}$ to a small value, the equilibrium that is stabilized is not anymore $(0, 0, 0, 0)$ but $(0, 0, 0, u_{\text{min}}/\delta_s)$ which is very close to it, remains very acceptable in practice compared to the uncontrolled attraction point (especially as only the density of the sterile male mosquitoes does not converge to 0 and, moreover, this density converges to a small value) and answers the mathematical problem described in Section 2.

A simpler control. We decided to simplify the control found by the RL algorithm. We wanted to see if there is really a need for a different regime when there is only very few mosquitoes. The motivation behind this is that this region does not influence much the cost function that the RL algorithm tries to optimize and the control might be less precise on this part. This leads to the following simplified formula for the feedback control:

$$v_{\text{reg}}(M + M_s, F + F_s) = \begin{cases} u_{\text{min}} & \text{if } \log \frac{M+M_s}{F+F_s} > \alpha_2, \\ u_{\text{max}} & \text{otherwise,} \end{cases} \quad (7.4)$$

for $M + M_s, F + F_s > 0$, where $\alpha_2 = 4$ is a constant found by regression.

In this case again, the control with $u_{\text{min}} = 0$ and no noise does not seem to ever lead to the convergence of the state $(E(t), M(t), F(t), M_s(t))$, with a notable difference however: there is no purely cyclic behavior and $(E(t), M(t), F(t))$ converges rapidly to $(0, 0, 0)$. The obstacle to the convergence manifests in large peaks that appear in the control feedback and are increasingly spaced in time, as evidenced in Fig. 9. When choosing again a small $u_{\text{min}} > 0$ the convergence to the equilibrium is recovered and seems to work for arbitrarily small $u_{\text{min}} > 0$ (see Fig. 9).
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Components of the state and control $v_{\text{reg}}$ with $u_{\text{min}} = 5$ and $u_{\text{max}} = 300000$ over a duration of 1000 days for 1000 simulations with random initial conditions uniformly sampled in $[0, 20K]^4$ (note that the interval $[0, 5K]^4$ used during training was much smaller). Due to the larger values, we have reduced $dt$ to $10^{-3}$ for these simulations for numerical stability. Each color correspond to a different simulation.</figcaption></figure>

<caption>**Table 2**</caption>

Statistics over 100 simulations with random initial conditions in $[0, 10K]^4$ using control $u_{\text{reg}}$ (see (7.1)) with $u_{\text{min}} = 5$ and $u_{\text{max}} = 300000$ over a duration of 800 days.

<table><thead><tr><th></th><th>200 days</th><th>400 days</th><th>600 days</th><th>800 days</th></tr></thead><tbody><tr><td>average |E| + |M| + |F|</td><td>5.08 × 10<sup>4</sup></td><td>8.02 × 10<sup>3</sup></td><td>9.63 × 10<sup>4</sup></td><td>6.00 × 10<sup>-1</sup></td></tr><tr><td>variance |E| + |M| + |F|</td><td>4.54 × 10<sup>7</sup></td><td>4.01 × 10<sup>6</sup></td><td>1.39 × 10<sup>3</sup></td><td>2.00 × 10<sup>-2</sup></td></tr><tr><td>maximum |E| + |M| + |F|</td><td>5.90 × 10<sup>4</sup></td><td>1.05 × 10<sup>4</sup></td><td>1.49 × 10<sup>2</sup></td><td>8.00 × 10<sup>-1</sup></td></tr><tr><td>average |M<sub>s</sub>|</td><td>2.21 × 10<sup>6</sup></td><td>6.94 × 10<sup>5</sup></td><td>1.67 × 10<sup>4</sup></td><td>5.03 × 10<sup>1</sup></td></tr><tr><td>variance |M<sub>s</sub>|</td><td>5.31 × 10<sup>10</sup></td><td>1.31 × 10<sup>10</sup></td><td>4.14 × 10<sup>7</sup></td><td>8.12 × 10<sup>1</sup></td></tr><tr><td>maximum |M<sub>s</sub>|</td><td>2.47 × 10<sup>6</sup></td><td>8.22 × 10<sup>5</sup></td><td>2.58 × 10<sup>4</sup></td><td>7.06 × 10<sup>1</sup></td></tr></tbody></table>

We numerically demonstrate that this simplified control seems stabilizing over a wide range of initial conditions, much larger than at training time. Fig. 10 shows the evolution of the control and components of the state for randomly sampled initial conditions. The control appears effective, quickly stabilizing the state in every simulation. Table 3 shows corresponding statistics, illustrating that the state rapidly converges to 0 (with the exception of $M_s$ which takes longer to converge to ensure the other components of the state do not regrow) with very small variance. Besides, we can notice by comparing Tables 2 and 3 that the simplified control leads to quicker convergence of the state than the control defined in (7.1)–(7.3).

**Sensitivity of the control** We previously showed in Fig. 10 that the control appears effective in a wide range of initial conditions outside of the training distribution. We now investigate the sensitivity

<caption>**Table 3**</caption>

Statistics over 100 simulations with random initial conditions in $[0, 10K]^4$ using control $v_{\text{reg}}$ (see (7.4)) with $u_{\text{min}} = 5$ and $u_{\text{max}} = 300000$ over a duration of 800 days.

<table><caption>**Table 4**</caption><thead><tr><th></th><th>200 days</th><th>400 days</th><th>600 days</th><th>800 days</th></tr></thead><tbody><tr><td>average |E| + |M| + |F|</td><td>4.88 × 10<sup>4</sup></td><td>6.89 × 10<sup>2</sup></td><td>2.47 × 10<sup>0</sup></td><td>2.00 × 10<sup>-3</sup></td></tr><tr><td>variance |E| + |M| + |F|</td><td>7.58 × 10<sup>7</sup></td><td>7.85 × 10<sup>4</sup></td><td>1.31 × 10<sup>0</sup></td><td>1.61 × 10<sup>-6</sup></td></tr><tr><td>maximum |E| + |M| + |F|</td><td>5.91 × 10<sup>4</sup></td><td>1.13 × 10<sup>3</sup></td><td>4.37 × 10<sup>0</sup></td><td>6.00 × 10<sup>-3</sup></td></tr><tr><td>average |M<sub>s</sub>|</td><td>2.50 × 10<sup>6</sup></td><td>1.29 × 10<sup>5</sup></td><td>2.20 × 10<sup>3</sup></td><td>4.17 × 10<sup>1</sup></td></tr><tr><td>variance |M<sub>s</sub>|</td><td>3.07 × 10<sup>-11</sup></td><td>2.95 × 10<sup>9</sup></td><td>5.20 × 10<sup>6</sup></td><td>2.17 × 10<sup>-7</sup></td></tr><tr><td>maximum |M<sub>s</sub>|</td><td>2.50 × 10<sup>6</sup></td><td>2.48 × 10<sup>5</sup></td><td>1.08 × 10<sup>4</sup></td><td>4.17 × 10<sup>1</sup></td></tr></tbody></table>

Statistics over 10,000 simulations with random initial conditions in $[0, 10K]^4$ and random system parameters as described in (7.5), using control $v_{\text{reg}}$ (see (7.4)) with $u_{\text{min}} = 5$ and $u_{\text{max}} = 300000$ over a duration of 800 days.

<table><thead><tr><th></th><th>200 days</th><th>400 days</th><th>600 days</th><th>800 days</th></tr></thead><tbody><tr><td>average |E| + |M| + |F|</td><td>9.71 × 10<sup>2</sup></td><td>6.68 × 10<sup>0</sup></td><td>6.18 × 10<sup>-2</sup></td><td>1.02 × 10<sup>-4</sup></td></tr><tr><td>variance |E| + |M| + |F|</td><td>3.98 × 10<sup>3</sup></td><td>3.38 × 10<sup>2</sup></td><td>5.37 × 10<sup>0</sup></td><td>8.78 × 10<sup>-3</sup></td></tr><tr><td>maximum |E| + |M| + |F|</td><td>8.15 × 10<sup>4</sup></td><td>3.31 × 10<sup>4</sup></td><td>5.36 × 10<sup>2</sup></td><td>8.75 × 10<sup>-1</sup></td></tr></tbody></table>

of the control with respect to the system parameters (see Table A.5). Fig. 11 considers modifying one system parameter while leaving the others to their default values, and shows that the controlled system still converges for a wide range of parameter values, most of which are outside of the training distribution since the training happened with a single set of parameters.

Secondly, we consider randomly perturbing the system parameters simultaneously:

$$ 
\begin{aligned}
\hat{\rho}_E & \sim \mathcal{U}(0.65, 9.03) \\
\hat{\nu}_E & \sim \mathcal{U}(0.01, 0.35) \\
\hat{\delta}_E & \sim \mathcal{U}(0.03, 1.0) \\
\hat{\delta}_F & \sim \mathcal{U}(0.03, 0.1) \\
\hat{\delta}_M & \sim \mathcal{U}(0.09, 0.2) \\
\hat{\nu} & \sim \mathcal{U}(0.39, 1.0)
\end{aligned} 
\qquad (7.5) $$

where $\mathcal{U}(a, b)$ is the uniform distribution on interval $[a, b]$. Table A.5 shows more details about reasonable and realistic parameter ranges; the parameters tested above have mostly been chosen from the green converging points in Fig. 11. Table 4 shows that over 10,000 random simulations with parameters chosen randomly according to (7.5), all of them are converging. Overall, this shows that the control trained with given parameters appears to be efficient even for different parameter values within a reasonably large range, which are in line with the uncertainties on the parameters that could be experienced in real life.

Secondly, we consider a more refined model where $F_s$ is the solution of an ODE. Namely, we add the following equation to system (2.1)–(2.4) (which is the last equation of (1) in [39] that we consider with $\beta = +\infty$ and $\mu_F = 0$):

$$ \dot{F}_s = \nu \nu_E \frac{M_s}{M + M_s} - \delta_F F_s \qquad (7.6) $$

In this case, we investigate whether the previously-obtained control still works. We show numerically that control (7.4) appears to still be stabilizing for the extended system (2.1)–(2.4), (7.6), as evidenced by Fig. 12.

## 8. Conclusion and perspective

In this paper, we have explored various feedback control strategies for a population control problem involving mosquitoes. We have investigated the use of nonlinear feedback control strategies that depend only on measurements. We have first considered the case where only
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Sensitivity of the control to variations in the system parameters. Green points correspond to parameter values for which the system converges at inference time using control $v_{reg}$ (see (7.4)), while orange points correspond to values where the control fails to generalize. The vertical red line shows the unique parameter value at which the controller was trained on; all of the green points correspond to generalization with parameters that the control had never seen during training. When one parameter is changed, all of the others are left to their default values (see Table A.5). We note that for the $\delta_F$ and $\delta_M$ rates, all of the subsequent points up to 1 that are not shown for visibility reasons are also converging.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Components of the state and control $v_{reg}$ for system (2.1)–(2.4), (7.6) with $u_{min} = 5$ and $u_{max} = 300000$ over a duration of 1000 days for 100 simulations with random initial conditions uniformly sampled in $[0, 10K]^4$. Each color correspond to a different simulation. The control $v_{reg}$ takes one new action every 7 days. In all the simulations, $E, M, F$ and $F_i$ appear to converge to 0, $M_s$ converges to a low value (around 40) and the control converges to $u_{min}$.</figcaption></figure>

the wild male mosquito density is measured, then the case where the male mosquito density is measured, and finally the case where both the male and female mosquito densities are measured.

In particular, we have used reinforcement learning to develop a control strategy that depends only on the total number of male and sterile male mosquitoes, as well as the total number of female and sterile female mosquitoes, at a given time. This control strategy appears to be robust to parameter uncertainties and effective, and has the

<caption>**Table A.5**</caption>

Parameters for the system (2.1)-(2.4). This includes typical value ranges, which can be found in a population of Aedes polynesiensis in French Polynesia, as well as our chosen values in this work. Further details can be found in [39,42–48].

<table><thead><tr><th>Parameter name</th><th>Typical interval</th><th>Value in our work<sup>a</sup></th><th>Unit</th></tr></thead><tbody><tr><td>β<sub>E</sub></td><td>[7.46, 14.85]</td><td>8</td><td>Day<sup>-1</sup></td></tr><tr><td>ν<sub>E</sub></td><td>[0.005, 0.25]</td><td>0.25</td><td>Day<sup>-1</sup></td></tr><tr><td>δ<sub>E</sub></td><td>[0.023, 0.046]</td><td>0.03</td><td>Day<sup>-1</sup></td></tr><tr><td>δ<sub>F</sub></td><td>[0.033, 0.046]</td><td>0.04</td><td>Day<sup>-1</sup></td></tr><tr><td>δ<sub>M</sub></td><td>[0.077, 0.139]</td><td>0.1</td><td>Day<sup>-1</sup></td></tr><tr><td>δ<sub>s</sub></td><td>–</td><td>0.12</td><td>Day<sup>-1</sup></td></tr><tr><td>ν</td><td>–</td><td>0.49</td><td></td></tr><tr><td>K</td><td>–</td><td>50 000</td><td></td></tr></tbody></table>

<sup>a</sup> We consider a larger array of values when assessing the sensitivity of the control with respect to the parameters of the model (see the related paragraph at the end of Section 7).

potential to be applied in practice to control mosquito populations and prevent the spread of diseases.

Our results highlight the usefulness of machine learning and control theory in developing effective control strategies for complex dynamical systems, especially when the measurements are limited.

An interesting question would be to find a Lyapunov function for the closed-loop systems with these feedbacks. This is a challenging question, even when the control feedback is simple [29]. Note that other machine learning based approach are dedicated to this problem of finding a Lyapunov function [40,41]. Further research in this field could lead to even more powerful techniques for controlling populations of pests and disease vectors.

### CRediT authorship contribution statement

**Kala Agbo dici:** Writing – review & editing, Writing – original draft, Investigation. **Jean-Michel Coron:** Writing – review & editing, Writing – original draft, Investigation. **Amaury Hayat:** Writing – review & editing, Writing – original draft, Investigation. **Nathan Lichtlé:** Writing – review & editing, Writing – original draft, Investigation.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors would like to thank Luis Almeida for fruitful discussions about the SIT model. This work was granted access to the HPC resources of IDRIS under the allocation 2024-AD011014527R1 made by GENCI. The authors also wish to thank the ANR-Tremplin StarPDE ANR-24-ERCS-0010.

During the preparation of this work the author(s) used an LLM in order to correct spellings in the introduction. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the publication.

### Appendix. Typical system parameters

The parameters we use in this work for the (2.1)-(2.4) system can be found in Table 4, including their typical value ranges and additional literature regarding how those specific parameters are derived to match realistic values.
# Data availability

Data will be made available on request.

# References

[1] J.-M. Coron, Control and nonlinearity, in: Mathematical Surveys and Monographs, vol. 136, American Mathematical Society, Providence, RI, 2007, p. xiv+426, http://dx.doi.org/10.1090/surv/136.

[2] D.P. Bertsekas, Dynamic Programming and Optimal Control. Vol. I, fourth ed., Athena Scientific, Belmont, MA, 2017, p. xix+555.

[3] F.L. Lewis, D. Liu, Reinforcement Learning and Approximate Dynamic Programming for Feedback Control, John Wiley & Sons, 2013, http://dx.doi.org/10.1002/9781118453988.

[4] C.I. Byrnes, A. Isidori, New results and examples in nonlinear feedback stabilization, Systems Control Lett. 12 (5) (1989) 437-442, http://dx.doi.org/10.1016/0167-6911(89)90080-7.

[5] D.E. Koditschek, Adaptive techniques for mechanical systems, in: Proc. 5th. Yale University Conference, 1987, pp. 259-265.

[6] J. Tsinas, Sufficient Lyapunov-like conditions for stabilization, Math. Control Signals Systems 2 (4) (1989) 343-357, http://dx.doi.org/10.1007/BF02551276.

[7] R.S. Sutton, A.G. Barto, Reinforcement learning: an introduction, second ed., in: Adaptive Computation and Machine Learning, MIT Press, Cambridge, MA, 2018, p. xxii+526.

[8] D.P. Bertsekas, Reinforcement learning and optimal control, in: Athena Scientific Optimization and Computation Series, Athena Scientific, Belmont, MA, 2019, p. xiv+373, ©2019, Second printing with editorial revisions.

[9] B. Recht, A tour of reinforcement learning: The view from continuous control, 2018, URL https://arxiv.org/abs/1806.09460, arXiv:1806.09460.

[10] S. Levine, C. Finn, T. Darrell, P. Abbeel, End-to-end training of deep visuomotor policies, J. Mach. Learn. Res 17 (1) (2016) 1334-1373.

[11] S. Gu, E. Holly, T. Lillicrap, S. Levine, Deep reinforcement learning for robotic manipulation with asynchronous off-policy updates, in: 2017 IEEE International Conference on Robotics and Automation, ICRA, 2017, pp. 3389-3396, http://dx.doi.org/10.1109/ICRA.2017.7989385.

[12] A.R. Mahmood, D. Korenkevych, G. Vasan, W. Ma, J. Bergstra, Benchmarking reinforcement learning algorithms on real-world robots, 2018, URL https://arxiv.org/abs/1809.07731. arXiv:1809.07731.

[13] D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre, D. Kumaran, T. Graepel, et al., A general reinforcement learning algorithm that masters chess, shogi, and go through self-play, Science 362 (6419) (2018) 1140-1144, http://dx.doi.org/10.1126/science.aar6404.

[14] O. Vinyals, I. Babuschkin, W.M. Czarnecki, M. Mathieu, A. Dudzik, J. Chung, D.H. Choi, R. Powell, T. Ewalds, P. Georgiev, et al., Grandmaster level in StarCraft II using multi-agent reinforcement learning, Nature 575 (7782) (2019) 350-354, http://dx.doi.org/10.1038/s41586-019-1724-z.

[15] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, M. Riedmiller, Playing Atari with deep reinforcement learning, 2013, URL https://arxiv.org/abs/1312.5602. arXiv:1312.5602.

[16] S. Spielberg, R. Gopaluni, P. Loewen, Deep reinforcement learning approaches for process control, in: 2017 6th International Symposium on Advanced Control of Industrial Processes, ADCONIP, IEEE, 2017, pp. 201-206, http://dx.doi.org/10.1109/ADCONIP.2017.7983780.

[17] S. Chen, K. Saulnier, N. Atanasov, D.D. Lee, V. Kumar, G.J. Pappas, M. Morari, Approximating explicit model predictive control using constrained neural networks, in: 2018 Annual American Control Conference, ACC, IEEE, 2018, pp. 1520-1527, http://dx.doi.org/10.23919/ACC.2018.8431275.

[18] L. Lu, X. Meng, Z. Mao, G.E. Karniadakis, DeepXDE: A deep learning library for solving differential equations, SIAM Rev. 63 (1) (2021) 208-228, http://dx.doi.org/10.1137/19M1274067.

[19] J. Kubalík, E. Derner, J. Žegklitz, R. Babuška, Symbolic regression methods for reinforcement learning, IEEE Access 9 (2021) 139697-139711, http://dx.doi.org/10.1109/ACCESS.2021.3119000.

[20] T.N. Mundhenk, M. Landajuela, R. Glatt, C.P. Santiago, D.M. Faissol, B.K. Petersen, Symbolic regression via neural-guided genetic programming population seeding, 2021, URL https://arxiv.org/abs/2111.00053. arXiv:2111.00053.

[21] L. Almeida, M. Duprez, Y. Privat, N. Vauchelet, Mosquito population control strategies for fighting against arboviruses, Math. Biosci. Eng. 16 (6) (2019) 6274-6297, http://dx.doi.org/10.3934/mbe.20191313.

[22] N. Alphey, L. Alphey, M.B. Bonsall, A model framework to estimate impact and cost of genetics-based sterile insect methods for dengue vector control, PLoS One 6 (10) (2011) e25384, http://dx.doi.org/10.1371/journal.pone.0025384.

[23] P.-A. Bliman, M.S. Aronna, F.C. Coelho, M.A. da Silva, Ensuring successful introduction of wolbachia in natural populations of *Aedes aegypti* by means of feedback control, J. Math. Biol. 76 (2018) 1269-1300, http://dx.doi.org/10.1007/s00285-017-1174-x.

[24] P.-A. Bliman, Feedback control principles for biological control of dengue vectors, in: 2019 18th European Control Conference, ECC, IEEE, 2019, pp. 1659-1664, http://dx.doi.org/10.23919/ECC.2019.8795984.

[25] L. Almeida, M. Duprez, Y. Privat, N. Vauchelet, Optimal control strategies for the sterile mosquitoes technique, J. Differential Equations 311 (2022) 229-266, http://dx.doi.org/10.1016/j.jde.2021.12.002.

[26] A. Cristofaro, L. Rossi, Backstepping control for the sterile mosquitoes technique: stabilization of extinction equilibrium, 2023, Preprint.

[27] P.-A. Bliman, D. Cardona-Salgado, Y. Dumont, O. Vasilieva, Implementation of control strategies for sterile insect techniques, Math. Biosci. 314 (2019) 43-60, http://dx.doi.org/10.1016/j.mbs.2019.06.002.

[28] P.-A. Bliman, Y. Dumont, Robust control strategy by the sterile insect technique for reducing epidemiological risk in presence of vector migration, Math. Biosci. 350 (2022) http://dx.doi.org/10.1016/j.mbs.2022.108856, Paper No. 108856, 23.

[29] K. Agbo bidi, L. Almeida, J.-M. Coron, Global stabilization of a sterile insect technique model by feedback laws, J. Optim. Theory Appl. 204 (2) (2025) 30, http://dx.doi.org/10.1007/s10957-024-02566-4.

[30] M. Krstic, J. Kanellakopoulos, P. Kokotovic, Nonlinear and adaptive control design, in: Adaptive and Cognitive Dynamic Systems: Signal Processing, Learning, Communications and Control, John Wiley and sons, 1995, p. xiii+563.

[31] E.D. Sontag, Mathematical control theory, second ed., in: Texts in Applied Mathematics, vol. 6, Springer-Verlag, New York, 1998, p. xvi+531, http://dx.doi.org/10.1007/978-1-4612-0577-7, Deterministic finite-dimensional systems.

[32] P.-A. Bliman, Basic offspring number and robust feedback design for the biological control of vectors by sterile insect release technique, 2024, HAL, 04757871 URL https://inria.hal.science/hal-04757871.

[33] K. Agbo bidi, Feedback stabilization and observer design for sterile insect technique models, Math. Biosci. Eng. 21 (6) (2024) 6263-6288, http://dx.doi.org/10.3934/mbe.2024274, URL https://www.aimpress.com/article/doi/10.3934/mbe.2024274.

[34] M.T.J. Spaan, Partially observable Markov decision processes, in: M. Wiering, M. van Otterlo (Eds.), Reinforcement Learning: State-of-the-Art, Springer, Berlin, Heidelberg, 2012, pp. 387-414, http://dx.doi.org/10.1007/978-3-642-27645-3_12.

[35] J. Schulman, S. Levine, P. Moritz, M.I. Jordan, P. Abbeel, Trust region policy optimization, 2017, URL https://arxiv.org/abs/1502.05477. arXiv:1502.05477.

[36] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, O. Klimov, Proximal policy optimization algorithms, 2017, URL https://arxiv.org/abs/1707.06347. arXiv:1707.06347.

[37] A. Raffin, A. Hill, A. Gleave, A. Kanervisto, M. Ernestus, N. Dormann, Stable-BASELINES3: Reliable reinforcement learning implementations, J. Mach. Learn. Res. 22 (268) (2021) 1-8, URL http://jmlr.org/papers/v22/20-1364.html.

[38] R. Gato, Z. Menéndez, E. Prieto, R. Argilés, M. Rodríguez, W. Baldoquín, Y. Hernández, D. Pérez, J. Anaya, I. Fuentes, et al., Sterile insect technique: successful suppression of an *Aedes aegypti* field population in Cuba, Insects 12 (5) (2021) 469, http://dx.doi.org/10.3390/Insects12050469.

[39] M. Strugarek, H. Bossin, Y. Dumont, On the use of the sterile insect release technique to reduce or eliminate mosquito populations, Appl. Math. Model. 68 (2019) 443-470, http://dx.doi.org/10.1016/j.apm.2018.11.026.

[40] A. Alfarano, F. Charton, A. Hayat, Discovering Lyapunov functions with transformers, in: The 3rd Workshop on Mathematical Reasoning and AI at NeurIPS'23, 2023.

[41] A. Alfarano, F. Charton, A. Hayat, Lyapunov functions: a long-standing open problem in mathematics, with symbolic transformers, 2024, Preprint.

[42] H. Hughes, N.F. Britton, Modelling the use of Wolbachia to control dengue fever transmission, Bull. Math. Biol. 75 (2013) 796-818, http://dx.doi.org/10.1007/s11538-013-9835-4.

[43] F. Rivière, Ecologie De Aedes (Stegomyia) Polynesiensis, Marks, 1951 Et Transmission De La Filariose De Bancroft En Polynesie, Université Paris- Sud, Paris (France), 1988.

[44] T. Suzuki, F. Sone, Breeding habits of vector mosquitoes of filariasis and dengue fever in Western Samoa, Med. Entomol. Zoöl. 29 (4) (1978) 279-286, http://dx.doi.org/10.7601/mez.29.279.

[45] E.W. Chambers, L. Hapairai, B.A. Peel, H. Bossin, S.L. Dobson, Male mating competitiveness of a Wolbachia-introgressed *Aedes polynesiensis* strain under semi-field conditions, PLoS Neglected Trop. Dis. 5 (8) (2011) e1271, http://dx.doi.org/10.1016/j.actatropica.2013.11.024, URL https://www.sciencedirect.com/science/article/pii/S0001706X13003471.

[46] L.K. Hapairai, J. Marie, S.P. Sinkins, H.C. Bossin, Effect of temperature and larval density on *Aedes polynesiensis* (Diptera: Culicidae) laboratory rearing productivity and male characteristics, Acta Trop. 132 (2014) S108-S115, http://dx.doi.org/10.1016/j.actatropica.2013.11.024, URL https://www.sciencedirect.com/science/article/pii/S0001706X13003471.

[47] L.K. Hapairai, M.A.C. Sang, S.P. Sinkins, H.C. Bossin, Population studies of the filarial vector *Aedes polynesiensis* (Diptera: Culicidae) in two island settings of French Polynesia, J. Med. Entomol. 50 (5) (2013) 965-976, http://dx.doi.org/10.1603/me112246.

[48] L. Hapairai, Studies on *Aedes Polynesiensis* Introgression and Ecology to Facilitate Lymphatic Filariasis Control [Ph.D. thesis], Oxford University, UK, 2013.