# Optimal Control Approach for Implementation of Sterile Insect Techniques

Pierre-Alexandre Bliman, Daiver Cardona-Salgado, Yves Dumont, Olga Vasilieva

► To cite this version:

Pierre-Alexandre Bliman, Daiver Cardona-Salgado, Yves Dumont, Olga Vasilieva. Optimal Control Approach for Implementation of Sterile Insect Techniques. Journal of Mathematical Sciences, 2024, 279 (5), pp.607-622. hal-02418508

HAL Id: hal-02418508

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

Submitted on 18 Dec 2019

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.
# Optimal Control Approach for Implementation of Sterile Insect Techniques

Pierre-Alexandre Bliman¹, Daiver Cardona-Salgado², Yves Dumont³,⁴,⁵, Olga Vasilieva⁶*

¹ Sorbonne Université, Université Paris-Diderot SPC, Inria, CNRS

Laboratoire Jacques-Louis Lions, équipe Mamba, Paris, France

² Universidad Autónoma de Occidente, Cali, Colombia

³ CIRAD, Umr AMAP, Pretoria, South Africa

⁴ AMAP, Univ Montpellier, CIRAD, CNRS, INRA, IRD, Montpellier, France,

⁵ University of Pretoria, Department of Mathematics and Applied Mathematics, South Africa

⁶ Universidad del Valle, Cali, Colombia

## Abstract

Vector or pest control is essential to reduce the risk of vector-borne diseases or crop losses. Among the available biological control tools, the Sterile Insect Technique (SIT) is one of the most promising. However, SIT-control campaigns must be carefully planned in advance in order to render desirable outcomes. In this paper, we design SIT-control intervention programs that can avoid the real-time monitoring of the wild population and require to mass-rear a minimal overall number of sterile insects, in order to induce a local elimination of the wild population in the shortest time. Continuous-time release programs are obtained by applying an optimal control approach, and then laying the groundwork of more practical SIT-control programs consisting of periodic impulsive releases.

Keywords: pest control, vector control, sterile insect technique, optimal control, periodic impulsive control, open-loop control.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>A sex-structured entomological model involving sterile males</td><td>3</td></tr><tr><td>3</td><td>Optimal control approach</td><td>5</td></tr><tr><td>4</td><td>Numerical simulations and discussion</td><td>9</td></tr><tr><td>4.1</td><td>Preliminaries</td><td>9</td></tr><tr><td>4.2</td><td>Continuous-time optimal release programs</td><td>9</td></tr><tr><td>4.3</td><td>Practical applications: suboptimal impulsive release programs</td><td>11</td></tr><tr><td>4.4</td><td>Discussion of results</td><td>12</td></tr><tr><td>5</td><td>Conclusions</td><td>15</td></tr></table>

<footnote>*Corresponding author: olga.vasilieva@correounivalle.edu.co</footnote>
# 1 Introduction

The presence and abundance of mosquito species is strongly correlated with progressive dissemination and persistence of various vector-borne diseases in human populations, and thus constitutes a major issue to health-care authorities in many countries. In the absence of effective vaccine against malaria (transmitted by *Anopheles* mosquito species), encephalitis and West Nile virus (transmitted by *Culex* species), as well as dengue fever, Chikungunya and Zika viruses (spread by *Aedes* species), the control of these infectious diseases is primarily centered on suppression of the local mosquito populations through the use of chemical substances (larvicides, insecticides).

However, due to the negative environmental impact of chemical substancies on non-target species and the detected resistance developed by mosquitoes to insecticides and larvicides, other methods have been proposed and evaluated. One of them is the so-called Sterile Insect Technique (SIT) that relies on massive releases of male insects sterilized by ionizing radiations or using incompatible strains of *Wolbachia* symbiont. After mating with sterilized or *Wolbachia*-carrying males, female mosquitoes duly lay eggs that never hatch into larvae, and the latter induces a progressive suppression of the target wild population (see [5] regarding the detailed description of SIT and its applications).

Similarly, pest, like fruit flies, are responsible of considerable losses in crops, and orchards, around the world. For instance, *Bactrocera dorsalis*, a highly invasive species, was found in South Africa in 2010 [13]. It was first monitored in Réunion island in 2017¹, and was also detected for the first time in Europe, southern Italy, in 2018. Unfortunately, to control such a pest, very few tools, apart from insecticides, exist.

SIT-control interventions have been modeled and studied theoretically and numerically in a large number of papers using temporal and spatiotemporal models with discrete, continuous or hybrid time scaling while applying open-loop, closed-loop or mixed modeling approaches to derive continuous-time and/or periodic impulsive release programs (see, e.g., [1, 2, 4, 9, 11, 18] and references therein). However, the proposed programs for SIT-control required either to go on with excessively abundant releases for a considerable time or to monitor the current sizes of wild populations in order to reach an eventual elimination. Several scholars have also applied the dynamic optimization approach seeking to minimize the overall number of sterile insects to be released during the SIT-control campaign while anticipating a certain fixed duration of the control intervention [6, 14].

The principal goal of this paper is to design the SIT-control intervention programs that exhibit a better tradeoff between different objectives set by the decision-makers. The primary objective of SIT-control campaigns is to ensure the elimination or drastic reduction in the number of wild insects within some target locality in a finite time. Other (secondary) objectives are related to reducing the underlying costs of SIT-control interventions, which are usually expressed by key indicators such as the overall duration of the SIT-control campaign, the number of releases to be performed, and the number of sterilized males to be mass-reared for releases. Some SIT-control programs require for periodic assessments of the wild population sizes in order to determine the sizes of subsequent releases. Such measurements constitute extra costs of the SIT-control programs and must be accounted for. Therefore, it is essential to identify the better tradeoff between different secondary objectives and to design a SIT-control program that requires lower overall costs. In other words, we are seeking to design the control intervention programs that require not only a moderate control effort (expressed through the cumulative number of sterile insects to be released), but also induce local elimination of wild insects in a reasonable time and without compelling for real-time assessments of the wild population sizes.

For that purpose, we make use of an entomological model initially proposed and thoroughly analyzed in [2] (its key features are given in Section 2) and then formulate an optimal control problem with free terminal time and a multi-criteria objective functional (Section 3). This problem is reduced, by applying the traditional technique [10], to an optimality system with a time-optimality condition that is solved numerically to obtain the optimal release programs in continuous time while considering different priorities in the decision-making. Furthermore, we adopt the continuous-time optimal release programs for practical implementation and propose “almost optimal” (further referred to as *suboptimal*) SIT-control programs consisting of periodic impulsive releases of sterile insects. For the numerical simulations we consider the case of *Aedes spp* mosquito population, like in [2]. The results of all numerical simulations together with an underlying discussion that highlights the performance of the purely open-loop suboptimal release programs as well as their comparison

<footnote n="1">¹See the report given in http://daaf.reunion.agriculture.gouv.fr/Detection-d-une-nouvelle-mouche</footnote>
with mixed open/closed-loop impulsive release programs [2] are presented in Section 4, and Section 5 provides some final remarks and conclusions.

## 2 A sex-structured entomological model involving sterile males

Let $M(t)$ and $F(t)$ be the current population sizes of males and females present in some target locality at the time $t \ge 0$. Without loss of generality, we assume in the sequel that time is measured in days. Suppose that sterile males are being released in this locality at the per-day rate $u(t) \ge 0$ and denote their population size $M_S(t)$ at each $t \ge 0$. The release rate can be modeled by a piecewise continuous function $u(\cdot) \in PC$ defined for all $t \ge 0$. This function represents the decision or control variable and expresses the number of sterile males $M_S$ introduced at the day $t$ in the target locality. In other words, $u(t), t \ge 0$ defines an external control action over the natural population dynamics of the wild subpopulations, $M$ and $F$.

Using the modeling framework proposed by Bliman et al [2], the population dynamics of wild insects ($M$ and $F$) and their interaction with sterile males ($M_S$) can be described by the following three-dimensional ODE system:

$$ \dot{M} = r \rho \frac{FM}{M + \gamma M_S} e^{-\beta(M + F)} - \mu_M M, \quad M(0) = M^0 > 0 \quad (1a) $$

$$ \dot{F} = (1 - r) \rho \frac{FM}{M + \gamma M_S} e^{-\beta(M + F)} - \mu_F F, \quad F(0) = F^0 > 0 \quad (1b) $$

$$ \dot{M}_S = u(t) - \mu_S M_S, \quad M_S(0) = M_S^0 \ge 0. \quad (1c) $$

All the parameters of the model (1) are positive, and their detailed descriptions are provided in Table 1. The model assumes that all females $F$ are equally able to mate, with either wild males $M$ or sterile males $M_S$. However, viable offspring is produced by wild females only after their successful mating with wild males. Therefore, the recruitment rate of wild insects (cf. positive terms in Eqs. (1a)-(1b)) depends on the probability $M/(M + \gamma M_S)$ of effective contacts (or successful matings) between wild males and females. It is worth noting that sterilization of mass-reared male mosquitoes may affect their mating competitiveness [15], and the latter is expressed in the model by the coefficient $0 < \gamma \le 1$ denoting the relative mating efficiency or reproductive fitness of sterilized males (compared to the wild ones).

Eqs. (1a)-(1b) state that new individuals are recruited as either males or females in accordance with a primary sex ratio $r \div (1 - r)$, $r \in (0, 1)$ in viable offsprings that survive to adulthood. Here $\rho$ stands for the mean number of eggs deposited by a single female in average per day, while direct and/or indirect competition effect at different stages (larvae, pupae, adults) are included through the parameter $\beta$. Thus, lower values of $\beta$ imply that a larger fraction of eggs may survive to adulthood, whereas its higher values express stronger competition and fewer breeding sites.

<table><caption>Table 1: *Aedes spp* mosquito parameters of the model (1) with values borrowed from [2]</caption><thead><tr><th>Parameter</th><th>Description</th><th>Unit</th><th>Value (<i>Aedes</i> spp)</th></tr></thead><tbody><tr><td><i>r</i></td><td>Primary sex ratio</td><td>–</td><td>0.5</td></tr><tr><td><i>ρ</i></td><td>Average number of eggs deposited per female per day</td><td>day<sup>-1</sup></td><td>4.55</td></tr><tr><td><i>β</i></td><td>Characteristic of the competition effect per individual</td><td>–</td><td>3.57 × 10<sup>-4</sup></td></tr><tr><td><i>γ</i></td><td>Coefficient of mating efficiency of sterile males</td><td>–</td><td>1</td></tr><tr><td>μ<sub>M</sub></td><td>Natural mortality rate for male mosquitoes</td><td>day<sup>-1</sup></td><td>0.04</td></tr><tr><td>μ<sub>F</sub></td><td>Natural mortality rate for female mosquitoes</td><td>day<sup>-1</sup></td><td>0.03</td></tr><tr><td>μ<sub>S</sub></td><td>Natural mortality rate for sterile males</td><td>day<sup>-1</sup></td><td>0.04</td></tr></tbody></table>

Parameters $\mu_M$, $\mu_F$, and $\mu_S$ in (1) represent, respectively, the natural sex-specific mortality rates for wild males, wild females, and sterile males. Scientific evidence attests that mating females have higher average longevity than mating males [12], whereas sterilization of males may reduce their average lifespan [15]. Therefore, we suppose throughout the paper that:

$$ \mu_S \ge \mu_M \ge \mu_F. \quad (2) $$
Finally, Eq. (1c) states that sterile males $M_S$ are “exogenously” released at a time-dependent instantaneous rate $u(t)$. Here we assume that $u: \mathbb{R}_+ \to [0, u_{\max}]$ is a piecewise continuous function, and $u_{\max} > 0$ is the maximum capacity of releases. For a known $u(t)$, Eq. (1c) can be explicitly solved:

$$M_S(t; u(\cdot)) = e^{-\mu_S t} \left[ M_S^0 + \int_0^t e^{\mu_S \xi} u(\xi) d\xi \right] \quad (3)$$

According to (3), $M_S(t; u(\cdot))$ remains nonnegative and bounded from above whenever $u(t)$ is nonnegative and bounded from above. Additionally, the right-hand sides of Eqn. (1a)-(1b) are decreasing in $M_S(t)$. On the other hand, Bliman et al [2] have shown that, in absence of sterile males (that is, with $M_S(t) = 0$), the sex-structured entomological model<fnref n="2" /> has bounded trajectories that belong to an absorbing set

$$\mathcal{D} = \{ (M(t), F(t)) : 0 \le M(t) \le C, 0 \le F(t) \le C, \forall t \ge 0 \} \quad (4)$$

for some $C > 0$.

Since the right-hand sides of the dynamical system in (1) are Lipschitz-continuous, the initial-value problem (1) has a unique solution for all piecewise continuous and bounded $u(t) \in [0, u_{\max}], t \ge 0$.

In the context of this model, the upper bound of $u(t)$ (expressed by the maximum release capacity $u_{\max}$) plays an essential role, and its value must be defined in accordance with the model's parameters (presented in Table 1). A recent study performed by Bliman et al [2] can shed some light on an adequate choice of $u_{\max}$.

In particular, the authors of [2] had thoroughly analyzed the case when $u(t) = \Lambda$, $\Lambda = \text{const} > 0$ for all $t \ge 0$. Under such settings, Eq. (1c) becomes completely decoupled from the system (1) and its steady-state value $M_S^\# = \frac{\Lambda}{\mu_S}$ can be replaced in Eqs. (1a)-(1b) thus leading to the reduced ODE system

$$\dot{M} = r\rho \frac{FM}{M + \gamma M_S^\#} e^{-\beta(M + F)} - \mu_M M \quad (5a)$$

$$\dot{F} = (1-r)\rho \frac{FM}{M + \gamma M_S^\#} e^{-\beta(M + F)} - \mu_F F \quad (5b)$$

for naturally persistent wild populations. Let us recall that pest or mosquito populations are referred to as “naturally persistent” when one individual (male or female) produces, on average, more than one individual. It fact, Bliman et al [2] introduced two positive constants, $N_F$ and $N_F$, that represent the so-called basic offspring numbers related to the wild female and male populations, respectively. In the sequel, we will be dealing with naturally persistent insect populations. Therefore, we suppose that for the parameters of the model (defined in Table 1) it holds that

$$N_F := \frac{(1-r)\rho}{\mu_F} > 1, \quad N_M := \frac{r\rho}{\mu_M} > 1. \quad (6)$$

For the reduced system (5), the following fundamental results has been established in [2].

**Theorem 1** (Existence of positive equilibria for the reduced SIT system (5), [2]). Assume that $N_F > 1$. Then there exists $\Lambda^{crit} > 0$ such that system (5) admits

* two positive distinct equilibria if $0 < \Lambda < \Lambda^{crit}$,

* one positive equilibrium if $\Lambda = \Lambda^{crit}$,

* no positive equilibrium if $\Lambda > \Lambda^{crit}$.

<footnote n="2">²This model is obtained from Eqs. (1a)-(1b) when $M_S(t) = 0$ for all $t \ge 0$.</footnote>
The value of $\Lambda^{crit}$ is uniquely determined by the formula

$$ \Lambda^{crit} := 2 \frac{\mu_S \phi^{crit}(\mathcal{N}_F)}{\beta \gamma (1 + \frac{\mathcal{N}_F}{\mathcal{N}_M})}, \qquad (7) $$

where $\phi^{crit} := \phi^{crit}(\mathcal{N}_F)$ is the unique positive solution to the transcendental equation

$$ 1 + \phi \left( 1 + \sqrt{1 + \frac{2}{\phi}} \right) = \mathcal{N}_F \exp \left[ - \frac{2}{1 + \sqrt{1 + \frac{2}{\phi}}} \right]. \qquad (8) $$

**Theorem 2** (Stability of the insect-free equilibrium, [2]). If the reduced system (5) admits no positive equilibrium (that is, if $\Lambda > \Lambda^{crit}$), then the trivial equilibrium $(0,0) \in \mathbb{R}^2$ of (5) is globally exponentially stable.

Formal proofs of Theorems 1 and 2 are given in [2]. These theorems jointly provide a threshold value $\Lambda_{crit} > 0$ of the constant release rate such that for all constant release rates $u(t) = \Lambda > \Lambda_{crit}, t \ge 0$ the trajectories of the SIT system (1) engendered by any nonnegative set of initial conditions $(M^0, F^0, M_S^0)$ converge globally exponentially toward its unique equilibrium $E_0^\sharp = (0, 0, M_S^\sharp)$ free of wild insects, and the latter ultimately ensures an eventual elimination of the wild population.

Therefore, it is reasonable to take $u_{max}$ larger than $\Lambda_{crit}$.

In the following section, we apply the dynamic optimization approach in order to design a time-dependent release program under which local elimination of the wild population can be reached in minimum time. Such program is characterized by the shape of piecewise continuous function $u(t)$ that expresses for each $t$ an instantaneous release rate of sterile males.

## 3 Optimal control approach

Our primary goal is to define the time-dependent release rate $u^*(t)$ that ensures local elimination of the wild population in minimum time $T^* \in (0, \infty)$ while also minimizing the total (cumulative) number of sterile males to be released in the target locality during the whole campaign, that is, during the period $[0, T^*]$.

Let us suppose that the terminal time of control action $0 < T < \infty$ is set free, and then define several objectives of the decision-making:

1. Reaching elimination of the wild population in finite time.

2. Minimizing the overall time of control action.

3. Minimizing the underlying costs of control action.

The first objective can be formalized by imposing an endpoint condition

$$ F(T) = \varepsilon, \qquad (9) $$

where $\varepsilon \to 0^+$ is specified by the decision-maker. Strictly speaking, the endpoint condition (9) should be of the form $F(T^*) = 0$. However, given the exponential nature of state equations (1a)-(1b), the trajectory of $F(t)$ may only approach zero asymptotically when $t \to \infty$ but it cannot reach this value in finite time. Therefore, reaching elimination of wild females at some finite $0 < T < \infty$ can be modeled by (9) with some small enough $\varepsilon > 0$. Elimination of wild females inevitably entails the elimination of wild male insects, since wild males are progenies of wild females. Therefore, condition (1) ensures the elimination of the total wild population (males and females).

Time minimization (the second objective) can be expressed in the standard way by the identity $T = \int_0^T dt$, while the costs associated with a certain release program $u(t), t \in [0, T]$ (the third objective) are conveyed
through the cumulative number of sterile males (CNSM) released during the whole campaign, which is exactly assessed by the following formula:

$$ \mathrm{CNSM}(u) = \int_{0}^{T} u(t)dt. \qquad (10) $$

Thus, the above-mentioned three objectives can be altogether expressed by the following performance index (or objective functional):

$$ \mathcal{J}(u, T) = A_1(F(T) - \varepsilon)^2 + \int_0^T \left[ A_2F(t) + A_3 + \frac{1}{2}A_4u^2(t) \right] dt \qquad (11) $$

The nonnegative coefficients $A_i, i = 1, 2, 3, 4$ in (11) should be adequately chosen in order to reflect the desirable priorities of decision-making. Namely, $A_1$ defines the top priority of reaching elimination in finite time$^3$, $A_2$ characterizes the importance of “elimination intensity” of wild females $F(t)$ during the release campaign, $A_3$ expresses the time appreciation ($T = \int_0^T dt$), and $A_4$ refers to the costs of control effort (i.e., mass-rearing of sterile males).

It is worthwhile to point out that we assume no linear relationship between the coverage of control interventions and their respective costs. However, we do assume that the marginal cost of control action, i. e. $A_4u(t)$, is proportional to the control effort at each $t \in [0, T]$ which, in its turn, is expressed by the number sterile males to be released at each $t$. For that reason, the integrand in (11) is assumed quadratic with respect to control variable $u$. Under such assumption, the “direct” costs of a certain release program $u(t), t \in [0, T]$ can be expressed by $A_4 \times \mathrm{CNSM}(u)$, where $\mathrm{CNSM}(u)$ is given by (10), while other “indirect” costs can be assessed by the length $T$ of the SIT-control intervention. We are aware that this assumption is not very appealing but, at this stage, it is the most convenient one.

Our goal is to find an optimal release program $u^*(t) \in [0, u_{\max}], t \in [0, T^*]$ and the minimum time $T^* \in (0, \infty)$ that jointly minimize the performance index (11) subject to the wild population dynamics given by (1). In other words, we seek to solve the following optimal control problem:

$$ \min_{\substack{0 \le u(t) \le u_{\max} \\ 0 < T < \infty}} \mathcal{J}(u, T) = A_1(F(T) - \varepsilon)^2 + \int_0^T \left[ A_2F(t) + A_3 + \frac{1}{2}A_4u^2(t) \right] dt \qquad (12a) $$

subject to

$$ \dot{M} = r\rho \frac{FM}{M + \gamma M_S} e^{-\beta(M + F)} - \mu_M M, \qquad M(0) = M^0 > 0 \qquad (12b) $$

$$ \dot{F} = (1 - r)\rho \frac{FM}{M + \gamma M_S} e^{-\beta(M + F)} - \mu_F F, \qquad F(0) = F^0 > 0 \qquad (12c) $$

$$ \dot{M}_S = u(t) - \mu_S M_S, \qquad M_S(0) = M_S^0 \ge 0. \qquad (12d) $$

with

$$ u(\cdot) \in PC, \quad 0 \le u(t) \le u_{\max}, \quad \forall t \in [0, T]. \qquad (12e) $$

Existence of solution to the optimal control problem (12) is justified by the following result.

**Theorem 3.** Assume that $N_F > 1$ and $u_{\max} > \Lambda_{crit}$. Then optimal control problem (12) has a solution $(u^*, T^*)$ such that

$$ \min_{\substack{0 \le u(t) \le u_{\max} \\ 0 < T < \infty}} \mathcal{J}(u, T) = \mathcal{J}(u^*, T^*). $$

<footnote>$^3$Here $A_1$ stands for a penalty for violation of the endpoint condition (9), and therefore we suppose that $A_1 \gg A_i, i = 2, 3, 4$.</footnote>
Proof. To prove this theorem, we employ the classical approach based on the Filippov-Cesari existence theorem (thoroughly described in [7, 17]). Roughly speaking, we have to show that our optimal control problems (12) fulfils the following sufficient conditions for existence of an optimal solution:

(i) Solution $\mathbf{X}(t) := (M(t), F(t), M_S(t))'$ of the dynamical system (12b)-(12d) is well-defined and unique for each admissible $u(t) \in [0, u_{max}]$, $t \ge 0$, while the set of all $\mathbf{X}(t)$ is non-empty and bounded $\forall u(t) \in [0, u_{max}]$ and $\forall t \ge 0$.

(ii) The sets of all initial and terminal states $\mathbf{X}(0) = (M(0), F(0), M_S(0))$, $\mathbf{X}(T) = (M(T^*), F(T^*), M_S(T^*))$ are closed and bounded in $\mathbb{R}_+^3$.

(iii) The control set $[0, u_{max}]$ is closed, bounded and convex in $\mathbb{R}$.

(iv) The right-hand sides of the dynamical system (12b)-(12d) are linear in the control variable $u$.

(v) The terminal-state function $\phi(\mathbf{X}) := A_1(F - \varepsilon)^2$ is continuous in its arguments.

(vi) The integrand of (12a) is convex in $u$ and satisfies the coercivity condition $A_2F(t) + A_3 + \frac{1}{2}A_4u_1^2(t) \ge \mathcal{K}_1\|u\|^\alpha - \mathcal{K}_2$ with some constants $\mathcal{K}_1 > 0$, $\alpha > 1$, and $\mathcal{K}_2$.

Items (i)-(ii) have been already corroborated in Section 2 where it was established that dynamical system (12b)-(12d) has a unique solution for any $u(t) \in [0, u_{max}]$ that exists for all $t \ge 0$. Moreover, the trajectories $M(t), F(t), M_S(t)$ remain bounded for all $t \ge 0$ since they belong to the compact absorbing set

$$ \mathcal{Z} := \left\{ (M(t), F(t)) \in \mathcal{D}, 0 \le M_S(t) \le \frac{u_{max}}{\mu_S}, \forall t \ge 0 \right\} $$

if engendered by $(M(0), F(0), M_S(0)) \in \mathcal{Z}$. In the above definition of $\mathcal{Z}$, $\mathcal{D}$ refers to (4).

The credibility of items (iii)-(v) is beyond doubt, while item (vi) is fulfilled with $\mathcal{K}_1 = \frac{1}{2}A_4$, $\alpha = 2$ and $\mathcal{K}_2 = 0$ for all $t \ge 0$. This completes the proof. $\square$

Formal solution of the optimal control problem (12) can be achieved by applying the Pontryagin maximum principle that constitutes the necessary condition for optimality of $(u^*, T^*)$, which exists in virtue of Theorem 3. In particular, we are interested in the variant of maximum principle applicable to optimal control problems with free terminal time, concisely described by Lenhart and Workman [10]. The Hamiltonian function $H(M, F, M_S, u, \lambda_1, \lambda_2, \lambda_3) : \mathcal{Z} \times \mathbb{R}^3 \mapsto \mathbb{R}$ associated with the optimal control problem (12) is defined as

$$ H(\mathbf{X}, u, \lambda) = -A_2F - A_3 - \frac{1}{2}A_4u^2 + \lambda_1 \left[ r\rho \frac{FM}{M + \gamma M_S} e^{-\beta(M + F)} - \mu_M M \right] \\ + \lambda_2 \left[ (1 - r)\rho \frac{FM}{M + \gamma M_S} e^{-\beta(M + F)} - \mu_M M \right] + \lambda_3 [u(t) - \mu_S M_S] \quad (13) $$

where $\lambda := (\lambda_1, \lambda_2, \lambda_3)'$ can be viewed as a vector of Lagrange multipliers linked to differential constraints (12b)-(12d). On the other hand, $\lambda_i = \lambda_i(t)$, $i = 1, 2, 3$ are time-varying real adjoint functions that express the marginal variations in the value of objective functional $\mathcal{J}(u, T)$ induced by changes in the current values of state variables $(M(t), F(t), M_S(t))$ (in the “component-by-component” sense). Thus, the current values of $\lambda_i = \lambda_i(t)$, $i = 1, 2, 3$ reflect additional benefits or costs associated with changes in $(M(t), F(t), M_S(t))$ and they are necessary elements of the Pontryagin maximum principle [10], which is formulated as follows.

Let $(u^*, T^*)$ be an optimal pair for (12) in the sense that $u^*(t)$ is a piecewise continuous real function with domain $[0, T^*]$ and range $[0, u_{max}]$ and

$$ \mathcal{J}(u^*, T^*) \le \mathcal{J}(u, T) $$
for all other controls u and times T. Let $\mathbf{X}^*(t) = (M^*(t), F^*(t), M_S^*(t))'$ be the corresponding state defined for all $t \in [0, T^*]$. Then there exists a piecewise differentiable adjoint function $\lambda: [0, T^*] \mapsto \mathbb{R}^3$ satisfying the adjoint differential equation

$$ \dot{\lambda} = -H_{\mathbf{X}}(\mathbf{X}^*, u^*, \lambda) \quad (14) $$

with transversality conditions

$$ \lambda_1(T^*) = 0, \quad \lambda_2(T^*) = -\phi_{\mathbf{X}}|_{\mathbf{X}=\mathbf{X}(T^*)} = -2A_1(F^*(T^*) - \varepsilon), \quad \lambda_3(T^*) = 0 \quad (15) $$

while $0 < T^* < \infty$ fulfills the time-optimality condition

$$ H(\mathbf{X}^*(T^*), u^*(T^*), \lambda(T^*)) = 0. \quad (16) $$

Furthermore, the Hamiltonian (13) has a critical point (maximum$^4$) at $u = u^*$, and therefore

$$ H(\mathbf{X}^*(t), u^*(t), \lambda(t)) \ge H(\mathbf{X}^*(t), u(t), \lambda(t)) \quad (17) $$

for any $u(t) : [0, T^*] \mapsto [0, u_{\max}]$ and almost for all $t \in [0, T^*]$.

According to [10], condition (17) can be written more precisely as

$$ \left. \begin{array}{ll} u^*(t) = 0 & \text{if } H_u < 0 \\ 0 < u^*(t) < u_{\max} & \text{if } H_u = 0 \\ u^*(t) = u_{\max} & \text{if } H_u > 0 \end{array} \right\} \quad (18) $$

that leads to a more compact form of $u^*(t)$ also known as the *characterization* of the optimal control:

$$ u^*(t) = \max \left\{ 0, \min \left\{ \frac{1}{A_4} \cdot \lambda_3(t), u_{\max} \right\} \right\}. \quad (19) $$

It is worth pointing out that the optimality conditions (18) are rather meaningful and provide interesting insights regarding the benefits and costs of control interventions. Namely, the condition

$$ H_u = -A_4 u + \lambda_3 = 0 $$

implies that, under optimal release program $u^*$ and at each $t$, the marginal cost of control intervention (given by the term $A_4 u$) must be equal to its marginal benefit (expressed by the term $\lambda_3$). If the marginal cost of $u^*$ is higher than its marginal benefit (that is, $H_u < 0$ in (18)) then it is optimal not to implement this program at all, and therefore we set $u^*(t) = 0$. On the other hand, if the marginal cost of $u^*$ is lower than its marginal benefit (that is, $H_u > 0$ in (18)), then it is optimal to use all available resources, and therefore we set $u^*(t) = u_{\max}$.

Using the characterization of $u^*$ given by (19), the original optimal control problem (12) can be reduced to a two-point boundary value problem. This boundary value problem (also known as *optimality system*) is composed by six differential equations with six endpoint conditions, namely:

* three original equations with underlying initial conditions (12b)-(12d) where $u(t)$ is replaced by its characterization (19);

* three adjoint equations (14) with transversality conditions (15) where $u(t)$ is also replaced by its characterization (19).

The optimal time $0 < T^* < \infty$ is defined by the time-optimality condition (16).

Due to non-linearity, complexity, and high dimension of the optimality system described above, it can only be solved numerically, and the forthcoming section is fully devoted to this issue, with a particular focus on *Aedes spp* parameters, given in Table 1.

<footnote n="4">⁴It is straightforward to verify that $H_{uu} = -A_4 < 0$ for all u.</footnote>
# 4 Numerical simulations and discussion

## 4.1 Preliminaries

Our numerical simulations are focused on the "worst scenario", that is, supposing that wild mosquitoes have equilibrium densities at the commencement of control intervention ($t = 0$), while sterile males are not present in the target locality ($M_S(0) = 0$). Bliman et al [2] had deduced that the positive steady state $E^\ddagger = (M^\ddagger, F^\ddagger)$ of the sex-structured system (1) (or (5)) with $M_S(t) \equiv 0$ has coordinates

$$M^\ddagger := \frac{N_M}{N_F + N_M} \frac{1}{\beta} \ln N_F \quad \text{and} \quad F^\ddagger := \frac{N_F}{N_F + N_M} \frac{1}{\beta} \ln N_F, \qquad (20)$$

where the basic offspring numbers $N_F$, $N_M$ satisfy the relationships (6). Thus, initial conditions for the dynamical system (12b)-(12d) are defined as

$$M(0) = M^\ddagger, \quad F(0) = F^\ddagger, \quad M_S(0) = 0.$$

For the parameter values given in the last column of Table 1, we have that $N_F \approx 75.83$, $N_M \approx 56.87$, while the mosquito densities at equilibrium $E^\ddagger = (M^\ddagger, F^\ddagger)$ are $M^\ddagger \approx 5.19 \times 10^3$ and $M^\ddagger \approx 6.93 \times 10^3$ individuals per hectare (indv/ha).

We assume that elimination of wild population is reached when the endpoint condition (9) is fulfilled with $\varepsilon = 10^{-1}$ meaning that no female mosquitoes are left in the target locality ($F^*(T^*) < 1$).

The value of $u_{\max}$ must be set greater than $\Lambda_{\text{crit}} \approx 1.29 \times 10^3$ (see Theorems 1 and 2 in Section 2). On the one hand, $u_{\max}$ should not be too high in order to avoid extra-large releases and to maintain feasible the overall cost of the control intervention. On the other hand, $u_{\max}$ must be large enough in order to potentiate a faster elimination of wild populations. Therefore, we set

$$u_{\max} = 2.5 \times 10^3$$

as the maximum release capacity in all ensuing simulations. This choice is also linked to the production capacity of the sterile insect factory, and the proper value of $u_{\max}$ must be defined by practitioners.

The coefficients $A_i, i = 1, 2, 3, 4$ in (12a) must account for scaling of all quantities involved in the objective functional. In this study, we adopt the following scaling for $A_i, i = 1, 2, 3, 4$:

$$A_1 = \frac{P_1}{F^\ddagger}, \quad A_2 := \frac{P_2}{F^\ddagger}, \quad A_3 = \frac{P_3}{365}, \quad A_4 = \frac{P_4}{u_{\max}}, \qquad (21)$$

where $F^\ddagger$ refers to the steady-state value of wild females (20), time is scaled to years (= 365 days), and the control effort is normalized by the maximal release capacity $u_{\max}$. In the gauged values (21), the coefficients $P_i, i = 1, 2, 3, 4$ define “pure” priorities of the decision-making, which are dimensionless.

Our top priority is to reach elimination of wild mosquitoes in minimum time; therefore, the largest value is assigned to $P_1$ thus seeking to fulfill the endpoint condition (9) and to obtain that $A_1 \phi(\mathbf{X}^*(T^*)) = 0$. On the other hand, the lowest priority can be assigned to $P_4$ that stands for the rearing (unit) cost of one sterile males.

The intensity of wild mosquitoes elimination (expressed by $P_2$ in the objective functional (12a)) is also essential, and a higher value of $P_2$ must help to continue elimination even when the current value of $F(t)$ becomes relatively small. The role of time minimization ($P_3$) should be further explored by considering three options: (1) $P_3 = 0$, (2) $P_4 < P_3 < P_2$, and (3) $P_3 > P_2$.

Summarizing the above-mentioned considerations, we consider the following values of “pure” priority coefficients:

$$P_1 = 10^{10}, \quad P_2 = 10^4, \quad P_3 \in \{0, 10^3, 10^5\}, \quad P_4 = 1. \qquad (22)$$

## 4.2 Continuous-time optimal release programs

To solve numerically the optimality systems described in the closing part of Section 3 together with time optimality condition (16) for free terminal time, we employ the GPOPS-II solver<fnref n="5" /> designed for the MATLAB

<footnote n="5">⁵For more information regarding the GPOPS-II solver please visit http://gpopsc2.com/. A concise description of this solver and its capacities is also available in [3, Appendix B].</footnote>
<figure><img src="image_1.png" /><figcaption>Figure 1: *Upper row:* Optimal release programs $u^*(t)$, $t \in [0, T^*]$ for $P_3 \in \{0, 10^3, 10^5\}$. *Lower row:* Time evolution of mosquito populations under respective $u^*(t)$: wild males $\log_{10} M(t)$ (blue-colored curves), wild females $\log_{10} F(t)$ (red-colored curves), and sterile males $\log_{10} M_S(t)$ (yellow-colored curves).</figcaption></figure>

platform. This solver implements an adaptive numerical technique based on the *direct orthogonal collocation*, which is also known as *Radau pseudospectral method* [8, 16].

Additionally, GPOPS-II solver automatically scales all input intervals $[0, T]$ to the range $[-1, 1]$, thus it suffices to hold its standard numerical tolerance of $10^{-5}$ with regards to internal scaling.

Figure 1 presents the results of numerical solutions of the optimal control problem (12) for three sets of “pure” priority values assigned in (22). The upper row of Figure 1 displays very similar structure of the optimal release programs $u^*(t)$ for all three cases, what naturally shows off their robustness. However, as $P_3$ increases, the overall number of sterile insects needed for the program’s implementation grows larger, while the overall time $T^*$ of the control intervention slowly shortens down. This provides interesting insights for decision-makers by illustrating the existence of a certain tradeoff between the overall duration of the elimination campaign and its underlying costs.

Additionally, it is worthwhile to note that optimal release programs are suspended exactly at $t = T^*$, that is, when the endpoint condition (9) is finally met. The latter is induced by the transversality conditions (15), according to which no additional benefit is expected after $T^*$.

Even though all three optimal release programs $u^*(t)$, $t \in [0, T^*]$ displayed in Figure 1 have a very clear and straightforward structure, their tangible implementation may become unfeasible since the continuous-time releases are hardly practicable.

In the following subsection, we propose a more realistic type of release programs which are inspired by the optimal solutions $u^*(t)$ and, due to that reason, they are further referred to as “suboptimal” strategies based on impulsive releases of sterile males.
## 4.3 Practical applications: suboptimal impulsive release programs

From the practical standpoint, it is more realistic to consider periodic releases of sterile male mosquitoes every $n\tau$ days, where $n \in \mathbb{N}$ until meeting the endpoint condition (9). In fact, Bliman et al [2] studied an extension of SIT-control model (1), namely

$$ \dot{M} = r\rho \frac{FM}{M + \gamma M_S} e^{-\beta(M + F)} - \mu_M M, \quad M(0) = M^0 > 0, \quad (23a) $$

$$ \dot{F} = (1 - r)\rho \frac{FM}{M + \gamma M_S} e^{-\beta(M + F)} - \mu_F F, \quad F(0) = F^0 > 0, \quad (23b) $$

$$ \dot{M}_S = -\mu_S M_S, \quad \forall t \in \bigcup_{n \in \mathbb{N}} (n\tau, (n+1)\tau), \quad (23c) $$

$$ M_S(n\tau^+) = \tau\Lambda_n + M_S(n\tau^-), \quad n \in \mathbb{N}, \quad (23d) $$

where the piecewise continuous function $M_S(t)$ has jumps at $t = n\tau$, $M_S(n\tau^+)$ denote its right and left limits at $t = n\tau$, and $\Lambda_n$ stands for a “per day” number of sterile male mosquitoes to be released in the target locality. The relationship (23d) models that at $t = n\tau$, a lump quantity $\tau\Lambda_n$ of sterile insects is emitted to the wild population at once, whereas on a union of open intervals $(n\tau, (n+1)\tau)$, $n \in \mathbb{N}$ the mosquito populations $M(t), F(t)$, and $M_S(t)$ evolve according to Eqs. (23a)-(23c).

It is worth pointing out that Bliman et al [2] have proposed and studied different strategies for implementation of impulsive releases of sterile males that are engendered by the choice of $\Lambda_n$ in (23d). Namely, the open-loop (or feedforward) SIT-control strategies for impulsive releases have been designed by assuming $\Lambda_n = \Lambda$ constant for all $n \in \mathbb{N}$. Alternatively, choosing each $\Lambda_n$ in accordance with periodic measurements (either synchronized with releases or more sparse) of the wild population size has enabled the design of the closed-loop (or feedback) periodic impulsive strategies for SIT-control. Moreover, the combination of two above-mentioned approaches has resulted in the design of mixed (open/closed-loop) strategies for periodic impulsive SIT-control, and this control mode have actually rendered the best outcome expressed through a shorter overall time needed to reach elimination, smaller cumulative number of sterile insects to be released, and fewer number of releases to be effectively carried out during a whole SIT-control campaign.

On the other hand, mixed (open/closed-loop) strategies for periodic impulsive SIT-control still require for real-time assessments of wild population sizes $M(t)$ and $F(t)$, and this protocol is time-consuming and costly, even with sparse measurements carried out every $p\tau$ days where $p = 2, 3, \dots$. When $p = 1$, it is understood that the measurements are synchronized with the releases (the reader may find further details regarding the synchronized and sparse measurements in [2]).

In this paper, we propose another type of periodic impulsive SIT-control strategies that are expressly feedforward (or open-loop) and thus require no real-time assessments of the wild population sizes. Their design is performed on the grounds of continuous-time numerical solutions $u^*(t), t \in [0, T^*)$ to the optimal control problem (12) (obtained in Subsection 4.2), and for this reason we call them “suboptimal”. Contrary to the continuous case (Section 3) and the one considered by Bliman et al [2], we cannot rigorously prove the existence of an optimal solution, neither to show the convergence of $F(t)$ to $\epsilon$ in a finite time. Nevertheless, we explore this case numerically.

First, we define the extension of the optimal control program $u^*(t)$ as

$$ \hat{u}^*(t) = \begin{cases} u^*(t), & \text{if } t \in [0, T^*) \\ 0, & \text{if } t > T^* \end{cases} \quad (24) $$

For a chosen period $\tau$, we then construct a sequence of impulses $\{U_n\}$ using the following rule

$$ U_n = \max_{t \in [n\tau, (n+1)\tau]} \hat{u}^*(t), \quad n \in \mathbb{N}. \quad (25) $$

This sequence $\{U_n\}$ converges to zero and reaches this value in some finite time $\hat{T}^* < \infty$. In view of the relationship (25), it is expected that the overall number of released sterile insects during the whole SIT-control campaign be a bit higher under suboptimal impulsive releases than under continuous-time optimal
<figure><img src="image_2.png" /><figcaption>Figure 2: *Upper row:* suboptimal release programs {$U_n$} for $P_3 \in \{0, 10^3, 10^5\}$ with $\tau = 7$ days. *Lower row:* Time evolution of mosquito populations under respective {$U_n$}: wild males $\log_{10} M(t)$ (blue-colored curves), wild females $\log_{10} F(t)$ (red-colored curves), and sterile males $\log_{10} M_S(t)$ (yellow-colored curves).</figcaption></figure>

release programs $u^*(t)$. On the other hand, the endpoint condition (9) may be reached a bit sooner under suboptimal impulsive releases than under continuous-time optimal release programs $u^*(t)$, and therefore $\hat{T}^* \le T^*$. This is exactly what numerical experiments illustrate — see results displayed in Figures 2 and 3 for $\tau = 7$ days and $\tau = 14$ days, respectively.

As shown in the upper rows of Figures 2 and 3, the release pick-values ($\tau u_{\max}$ per hectare) are known in advance, and the decision-makers may choose the frequency of releases ($\tau$ days) in accordance with the mass-rearing capacities of sterile insects available at situ.

In the following subsection, we compare our results displayed by Figures 2 and 3 for feedforward suboptimal impulsive SIT-control programs with mixed open/closed-loop impulsive SIT-control programs designed in [2].

## 4.4 Discussion of results

It is worth noting that our simulations have been performed using exactly the same numerical values of the model's parameters (provided in Table 1) as in the work authored by Bliman et al [2]. Therefore, our results are comparable with those obtained in [2].

It should be pointed out that suboptimal impulsive SIT-control programs designed for $P_3 = 10^3$ (see middle columns in Figures 2 and 3) do not make much sense. Namely, they require the same time $\hat{T}^*$ to reach elimination of the wild mosquitoes as those designed for $P_3 = 0$ (see left columns in Figures 2 and 3) while
<figure><img src="image_3.png" /><figcaption>Figure 3: *Upper row:* suboptimal release programs {$U_n$} for $P_3 \in \{0, 10^3, 10^5\}$ with $\tau = 14$ days. *Lower row:* Time evolution of mosquito populations under respective {$U_n$}: wild males $\log_{10} M(t)$ (blue-colored curves), wild females $\log_{10} F(t)$ (red-colored curves), and sterile males $\log_{10} M_S(t)$ (yellow-colored curves).</figcaption></figure>

demanding to release more sterile insects. Therefore, these two suboptimal release programs are excluded from further consideration, and we focus on two key options with regards to priorities for decision-making, namely: (1) disregarding the time appreciation $P_3 = 0$; (2) encouraging the time appreciation $P_3 = 10^5$.

Table 4.4 summarizes the key features of simulation data for suboptimal (open-loop) and mixed (open/closed-loop, [2]) impulsive SIT-control programs, such as the cumulative number of sterile males needed for successful implementation of the SIT-based campaign, the number of weeks to reach elimination of wild populations, and the number of effective (nonzero) releases to be performed during the whole SIT campaign.

It should be noted that the cumulative number of sterile males needed for successful implementation of an impulsive suboptimal SIT-control program is calculated as $\tau \sum_{n=0}^{\infty} U_n$. In the first column of Table 4.4, the control gain is defined as $\left(\frac{1}{k} - 1\right)$ and the parameter $k$ is chosen to satisfy $0 < k < \frac{1}{N_F}$. Thus, smaller values of $k$ implies a larger control gain, while its larger values imply a smaller control gain$^6$. Furthermore, $p=1$ stands for the measurements that are synchronized with the releases, while $p=4$ expresses the sparser measurements to be carried out every $4\tau$ days.

According to the second column of Table 4.4, the suboptimal impulsive release programs require to mass-rear a lesser number of sterile mosquitoes than mixed open/closed-loop strategies with an exception of the

<footnote>$^6$The detailed explanation regarding the choice of $k$ is provided in [2], where two values of $k$ are considered: $k = 0.2/N_F$ and $k = 0.2/N_F$ that correspond to the control gain of the order ≈ 378 and ≈ 76, respectively.</footnote>
<table><caption>Table 2: Summary of simulation data for suboptimal (open-loop) and mixed (open/closed-loop) impulsive SIT-control programs</caption><thead><tr><th>Impulsive release program</th><th>Cumulative number of released sterile males</th><th>Number of weeks needed to reach elimination</th><th>Number of nonzero releases</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center;">Frequency of releases = 7 days</td></tr><tr><td>Suboptimal without time appreciation (P<sub>3</sub> = 0)</td><td>434.820</td><td>72</td><td>72</td></tr><tr><td>Suboptimal with time appreciation (P<sub>3</sub> = 10<sup>5</sup>)</td><td>457.500</td><td>70</td><td>70</td></tr><tr><td>Mixed with synchronized measurements and larger control gain (p = 1, k = 0.2/N<sub>F</sub>, [2])</td><td>450.668</td><td>72</td><td>72</td></tr><tr><td>Mixed with synchronized measurements and smaller control gain (p = 1, k = 0.99/N<sub>F</sub>, [2])</td><td>457.489</td><td>246</td><td>246</td></tr><tr><td>Mixed with sparse measurements and larger control gain (p = 4, k = 0.2/N<sub>F</sub>, [2])</td><td>534.849</td><td>65</td><td>53</td></tr><tr><td>Mixed with sparse measurements and smaller control gain (p = 4, k = 0.99/N<sub>F</sub>, [2])</td><td>450.077</td><td>69</td><td>53</td></tr><tr><td colspan="4" style="text-align: center;">Frequency of releases = 14 days</td></tr><tr><td>Suboptimal without time appreciation (P<sub>3</sub> = 0)</td><td>442.480</td><td>74</td><td>37</td></tr><tr><td>Suboptimal with time appreciation (P<sub>3</sub> = 10<sup>5</sup>)</td><td>464.050</td><td>72</td><td>36</td></tr><tr><td>Mixed with synchronized measurements and larger control gain (p = 1, k = 0.2/N<sub>F</sub>, [2])</td><td>465.187</td><td>72</td><td>36</td></tr><tr><td>Mixed with synchronized measurements and smaller control gain (p = 1, k = 0.99/N<sub>F</sub>, [2])</td><td>427.701</td><td>136</td><td>68</td></tr><tr><td>Mixed with sparse measurements and larger control gain (p = 4, k = 0.2/N<sub>F</sub>, [2])</td><td>499.497</td><td>66</td><td>25</td></tr><tr><td>Mixed with sparse measurements and smaller control gain (p = 4, k = 0.99/N<sub>F</sub>, [2])</td><td>449.099</td><td>74</td><td>28</td></tr></tbody></table>

one bearing $p = 1, k = 0.99/N_F$ (that is, a mixed strategy with synchronized measurements performed every $\tau = 14$ days and a smaller control gain for the feedback mode). However, this release program needs the longest time to reach the elimination of wild mosquitoes besides numerous measurements of wild population sizes it requires to perform.

Analyzing the data from the third and fourth columns of Table 4.4, it becomes clear that suboptimal release programs render about the same benefits as the mixed open/closed-loop programs based on the synchronized measurements while not requiring for assessments of wild population sizes and needing a lesser number of sterile insects. On the other hand, mixed open/closed-loop programs based on sparse measurements display better results than suboptimal open-loop programs in terms of the overall time needed to reach elimination and the number of effective (nonzero) releases. However, their advantages are counterpoised by the greater number of sterile insects needed for successful implementation of the SIT-control campaign and the extra costs for performing real-time assessments of wild population sizes.

In the end, the ultimate choice of the release program for SIT-control campaign must be made by practitioners after evaluating the realistic costs related to implementation of the SIT-control campaign, such as the costs for mass-rearing of larger/smaller cohorts of sterile insects, logistics costs for accomplishment of a single release, as well as the extra costs for taking a single real-time measurement of the wild population sizes.
# 5 Conclusions

In this work, we have applied the dynamic optimization approach for design of feedforward (open-loop) continuous-time programs for SIT-control, which do not require to assess the sizes of wild populations in real time. We have also proposed their more realistic suboptimal variants for practical implementation in the field.

All designed programs have a very clear structure and exhibit monotonicity with respect to the quantity of sterile insects to be released at each day $t$ (in the case of continuous-time optimal programs) or at the commencement of each period of $\tau$ days (in the case of impulsive suboptimal programs). The fact that release sizes gradually decrease during the SIT-control campaign may actually help in the adequate planning of the infrastructure and underlying logistics of the SIT-control interventions.

Another advantage of open-loop impulsive release program is the anticipated knowledge of the release pick-values ($\tau u_{\max}$) that enables the proper choice of the release frequency $\tau$ in accordance with the mass-rearing capacities of sterile insects available at situ.

We have also explored the impact of the exogenous parameter $P_3$ (expressing the time appreciation) on the structure of the designed release programs and their respective outcomes. From the latter, we have detected the existence of a certain tradeoff between the cumulative number of the released sterile insects needed to reach elimination and the overall duration of the SIT-control campaign. Namely, by increasing (significantly) the value of $P_3$ (ceteris paribus) the overall time $T^*$ (and also $\hat{T}^*$) can be (slightly) reduced on the cost of increasing the overall quantity of sterile males to be released during the SIT-based campaign (cf. left and right columns in Figures 1-3). However, the general structure of both optimal and suboptimal release programs exhibit robustness and resilience to variations in $P_3$.

In this work, we have tried to keep the value of another exogenous parameter – the daily mass-rearing capacity of sterile insects $u_{\max}$ – at a moderate realistic level of 2,500 individuals per hectare. However, a temperate enhancement (reduction) of $u_{\max}$ does not affect much the general structure of optimal and suboptimal release programs. The only difference we detected through numerical experiment (that are left beyond the scope of this paper) consists in less (more) initial releases to be performed at the maximum release capacity $u_{\max}$ per day (or $\tau u_{\max}$ every $\tau$ days).

Finally, the methodology presented in this paper can be easily adjusted to different biological characteristics of other mosquito or pest species (such as $\rho, \beta, \mu_M, \mu_F$) and of sterile insects ($\gamma, \mu_S$), while accounting for the total area of the target locality$^7$ and the underlying mass-rearing capacities to produce a sufficient quantity of sterile males every $\tau$ days. Therefore, the outcomes of this study could virtually help in the design and implementation of successful SIT-control campaigns.

# Acknowledgements

Support from the Colciencias and ECOS-Nord Program (Colombia: Project CI-71089; France: Project C17M01) is kindly acknowledged. DC and OV were supported by the inter-institutional cooperation program MathAmsud (18-MATH-05, MOVECO project). YD is partially endorsed by the GEMDOTIS project, funded by the call ECOPHYTO 2018 (Action 27). YD is also (partially) supported by the DST/NRF SARCHI Chair in Mathematical Models and Methods in Biosciences and Bioengineering at the University of Pretoria (grant 82770).

# References

[1] R. Anguelov, Y. Dumont, and J. Lubuma. Mathematical modeling of sterile insect technology for control of *Anopheles* mosquito. *Computers and Mathematics with Applications*, 64:374–389, 2012.

[2] P.-A. Bliman, D. Cardona-Salgado, Y. Dumont, and O. Vasilieva. Implementation of control strategies for sterile insect techniques. *Mathematical Biosciences*, 314:43–60, 2019.

<footnote>$^7$It must be stressed again that all our simulations were carried out for a “standardized” area comprising 1 hectare.</footnote>
[3] D. Campo-Duarte, O. Vasilieva, D. Cardona-Salgado, and M. Svinin. Optimal control approach for establishing *wMelPop Wolbachia* infection among wild *Aedes aegypti* populations. *Journal of Mathematical Biology*, 76(7):1907–1950, 2018.

[4] Y. Dumont and J.M. Tchuenche. Mathematical studies on the sterile insect technique for the Chikungunya disease and *Aedes albopictus*. *Journal of mathematical Biology*, 65(5):809–854, 2012.

[5] V.A. Dyck, J. Hendrichs, and A.S. Robinson, editors. *Sterile Insect Technique: Principles and Practice in Area-Wide Integrated Pest Management*. Springer, Dordrecht, The Netherlands, 2005.

[6] K. Fister, M. McCarthy, S. Oppenheimer, and C. Collins. Optimal control of insects through sterile insect release and habitat modification. *Mathematical biosciences*, 244(2):201–212, 2013.

[7] W. Fleming and R. Rishel. *Deterministic and stochastic optimal control*, volume 1 of *Applications of Mathematics*. Springer Verlag, New York, 1975.

[8] D. Garg, M. Patterson, C. Francolin, C. Darby, G. Huntington, W. Hager, and A. Rao. Direct trajectory optimization and costate estimation of finite-horizon and infinite-horizon optimal control problems using a Radau pseudospectral method. *Computational Optimization and Applications*, 49(2):335–358, 2011.

[9] M. Huang, X. Song, and J. Li. Modelling and analysis of impulsive releases of sterile mosquitoes. *Journal of biological dynamics*, 11(1):147–171, 2017.

[10] S. Lenhart and J. Workman. *Optimal Control Applied to Biological Models*. CRC Press, Taylor & Francis Group, Boca Raton, FL, 2007.

[11] J. Li and Z. Yuan. Modelling releases of sterile mosquitoes with different strategies. *Journal of biological dynamics*, 9(1):1–14, 2015.

[12] J. Liles. Effects of mating or association of the sexes on longevity in *Aedes aegypti* (L.). *Mosquito News*, 25(4):434–439, 1965.

[13] A. Manrakhan, J.H. Venter, and V. Hattingh. The progressive invasion of *Bactrocera dorsalis* (Diptera: Tephritidae) in South Africa. *Biological Invasions*, 17(10):2803–2809, 2015.

[14] L. Multerer, T. Smith, and N. Chitnis. Modeling the impact of sterile males on an *Aedes aegypti* population with optimal control. *Mathematical biosciences*, 311:91–102, 2019.

[15] C. Oliva, M. Jacquet, J. Gilles, G. Lemperiere, P.-O. Maquart, S. Quilici, F. Schooneman, M. Vreysen, and S. Boyer. The sterile insect technique for controlling populations of *Aedes albopictus* (Diptera: Culicidae) on Reunion Island: mating vigour of sterilized males. *PloS One*, 7(11):e49414, 2012.

[16] M. Patterson and A. Rao. GPOPS-II: A MATLAB software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming. *ACM Transactions on Mathematical Software (TOMS)*, 41(1):1, 2014.

[17] A. Seierstad and K. Sydsaeter. *Optimal control theory with economic applications*, volume 24 of *Advanced Textbooks in Economics*. Elsevier North-Holland, Inc., Amsterdam, 1987.

[18] M. Strugarek, H. Bossin, and Y. Dumont. On the use of the sterile insect release technique to reduce or eliminate mosquito populations. *Applied Mathematical Modelling*, 68:443–470, 2019.