# Optimal pest control for banana production under varying infestation density

Frank Kemayou, Walid Djema, Samuel Bowong, Frédéric Grognard, Suzanne

Touzeau

► To cite this version:

Frank Kemayou, Walid Djema, Samuel Bowong, Frédéric Grognard, Suzanne Touzeau. Optimal pest control for banana production under varying infestation density. 2025. hal-05247381

HAL Id: hal-05247381

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

Preprint submitted on 9 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Optimal pest control for banana production under
varying infestation density

Frank Kemayou<sup>1,2,3,4*</sup>, Walid Djema<sup>5</sup>, Samuel Bowong<sup>1,4</sup>,
Frédéric Grognard<sup>3</sup>, Suzanne Touzeau<sup>2,3</sup>

¹ University of Douala, Department of Mathematics and Computer science, Cameroon

² Université Côte d'Azur, INRAE, ISA, France

³ Université Côte d'Azur, Inria, INRAE, CNRS, MACBES, France

⁴ IRD, UMMISCO, Sorbonne University, Bondy, France

⁵ Université Côte d'Azur, Inria, INRAE, CNRS-Sorbonne Université, GREENOWL, France

* Corresponding author: kemayoufranck00@gmail.com

September 2025

Abstract

Banana and plantain production faces significant caused by many pests, the most harmful being *Radopholus similis*. This nematode infests and spends most of its life in the plant roots, which makes it difficult to control. Conventional methods face efficiency, cost, and environmental limitations, making it crucial to design optimized strategies. Our aim is to tackle this issue by formulating an optimal control problem on an original model of banana-pest interactions.

We propose a four-compartment model describing the host and pest dynamics. Its originality is that both infected roots and infesting nematodes are considered, which allows variable nematode densities in the roots. By employing time-scale separation, we reduce the model dimension and reveal a scenario where both stable pest-free and endemic equilibria coexist (backward bifurcation). We then include a control variable which reduces the root infestation rate. We formulate an optimal control problem aiming at maximizing the profit during a cropping season. The existence of an optimal control and necessary optimality conditions are established. We solve the optimal control problem by a direct method implemented in Bocop. For numerical settings of interest, the optimal solutions exhibit a bang-bang quasi-periodic structure. We extensively test whether candidate suboptimal strategies, based on periodic or constant control, can approximate the optimal performance. Numerical comparisons reveal that constant control performs poorly, whereas periodic control closely approaches the optimal solution. This confirms the overyielding effect and supports the realistic implementation of periodic control.

**Keywords:** Epidemiological model, Quasi-periodic control, Bang-bang control, Direct optimization

AMS keywords: 49K15 92D45 34C23

# 1 Introduction

Banana is the fourth most important crop in the world after rice, wheat and maize, as well as the second most traded fruit after citrus [14]. It is a major staple food in developing countries and it contributes to food security by producing fruit all year round, with an estimate of 119.83 million tonnes, including 40 million tonnes of plantain in 2020 [14]. In West and Central Africa, about 70 million people depend on these fruits for a large proportion of their daily carbohydrate intake [2]. Moreover, banana production is an important source of income and employment [14, 3].

Despite its significant production, banana is far from covering the increasing food needs of constantly growing populations [14]. Indeed, the banana sector is experiencing drops in yields due to nematodes, the most harmful being *Radopholus similis* [19], which causes up to 60% crop losses in many countries around the world [11]. It feeds and spends most of its life in the banana roots [7]. This particular aspect of nematode biology makes it difficult to study their behavior and to control them. Also, nematicides, the currently prevailing control, have a high environmental and financial cost (many are gradually being banned in Europe). Alternative control methods include the use of biostimulants, improved cultural practices, use of nematode tolerant or resistant banana cultivars, and biological control using fungi and bacterial hyperparasites [15].

Our objectives in this work are twofold: first, to formulate a mathematical model describing banana- nematode interactions based on the life cycle of *Radopholus Similis* presented in [23]; second, to propose and implement an optimal control strategy based on the use of biostimulants, which reduce nematode infestation [26]. The goal is to maximize the banana yield, while minimizing the control costs (i.e. to maximize profit), as
well as minimizing the nematode infestation at the end of a cropping season to avoid significant losses for the
next season.

Optimal control has previously been used to design efficient control strategies for various plants pests.
For instance, optimal fungicide applications using quadratic control functionals not easily linked to profit were
determined for the banana black spot disease [1] and Bayoud disease, caused by the fungus *Fusarium oxysporum*,
which spreads through contact with palm roots [4]. Similarly, optimal chemical insecticide applications have
been established to effectively control cocoa mirids, notably Sahlbergella singularis [24] and to control insects
responsible for the spread of viral plant diseases through an optimal combination of these chemical insecticides
and natural predators [8]. Finally, optimal strategies based on the release of an entomopathogenic fungus using
a profit related linear control functional were also established for effective control of the coffee berry borer [13].

To our knowledge, three models have been proposed for R. similis. The first is a cohort-based model to simulate the population dynamics of two major plant-parasitic nematode species of banana: *Radopholus similis* and *Pratylenchus coffeae* [25]. This model is fairly complex and requires extensive data to be calibrated. The second is a model designed to analyze the impact of climatic factors on the growth of *Radopholus similis*, without considering the direct interaction with banana plants [12]. It latter also incorporates biological and chemical controls, used to regulate this nematode population according to climatic conditions. The third is a multi-seasonal host-parasite model, coupling banana roots and *Radopholus similis*, in which two control strategies, against this pest are considered: pesticides and fallows [23]. However, our work differs from these three previous works not only by the optimal control application, but also by the simultaneous integration of both infected roots and infesting nematodes, with a dynamic variation of nematode densities in the roots.

This paper is organized as follows. Section 2 is devoted to the plant-pest model: we formulate the model in 2.1, derive its basic properties in 2.2, reduce it via Tikhonov's slow-fast system theory in 2.3, and analyse its equilibria in 2.4-2.5. Section 3 deals with control: we present the optimal control problem in 3.1, provide existence and necessary optimality conditions in 3.2, solve it numerically using the Bocop software in 3.3, and compare it with suboptimal controls in 3.4. Finally, we conclude the paper in Section 4.

2 Mathematical modeling and analysis of banana–pest interactions

**2.1 Formulation of the epidemiological model**

Herein, we formulate a mathematical model describing the infestation dynamics of banana roots by nematodes
during a cropping season. To build the model, we make the modelling assumptions described below.

1. The nematode population is divided into two compartments: free nematodes in the soil ($N_F$) and infesting nematodes in the roots ($N_I$).

2. The root biomass of banana plants is subdivided into two compartments: healthy roots ($S$) and infected roots ($I$). We define by $\mathcal{X} = \frac{N_I}{I}$ the infestation density.

3. Healthy roots have a constant recruitment ($\Lambda$).

4. Healthy roots ($S$) become infected ($I$) either when free nematodes enter the roots (rate $\beta_F$), or when infesting nematodes migrate in the roots (rate $\beta_I$). Free nematodes ($N_F$) then become infesting nematodes ($N_I$, with conversion parameter $c_\beta$); infesting nematodes remain infesting nematodes. Furthermore, reinfection may occur when some free nematodes enter infected roots (fraction $\eta$). As a consequence, the infestation density $\mathcal{X}$ may vary.

5. Some infesting nematodes leave the roots and migrate to the soil to become free nematodes ($N_F$), at a rate which increases with the infestation density (baseline rate $\gamma$).

6. Infesting nematodes ($N_I$) feed on the plant roots with a Holling type II functional response (rate $\rho$, half-saturation $a$).

7. Infesting nematodes reproduce inside the roots according to a logistic-like function. The growth rate is the sum of two terms: a constant intrinsic growth rate ($r$) and a feeding-dependent rate (with conversion parameter $c_\rho$). The carrying capacity ($K$) limits the infestation density ($\mathcal{X}$), rather than the number of infesting nematodes.

8. All roots undergo natural mortality (rate $\mu$).

9. Nematodes undergo natural mortality, which is higher for free nematodes than for infesting nematodes (rates $\nu_F > \nu_I$). Furthermore, we assume that intrinsic growth is higher than natural mortality for infesting nematodes ($r > \nu_I$), ensuring that nematodes strive in the roots (no Allee effect).
10. We include pest control, based on the use of biostimulants, which reduce nematode infestation rates $(\beta_F, \beta_I)$. The corresponding control variable is $u \in \mathcal{U} = \{u \in L^\infty(\mathbb{R}_+) : 0 \le u(t) \le u_{\max} \le 1, \forall t \in \mathbb{R}_+\}$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Scheme of the epidemiological model describing banana root infestation by nematodes. State variables: healthy roots (S), infected roots (I), free nematodes ($N_F$), and infesting nematodes ($N_I$). Solid arrows represent flows, dashed arrows interactions. The full model (black and grey arrows and boxes), given by equations (1), and the reduced model (black arrows and boxes only), given by equations (3), are represented.</figcaption></figure>

Under the hypotheses listed above, we obtain the compartmental model depicted in Figure 1 and represented by the following system of nonlinear ordinary differential equations

$$ \dot{S} = \Lambda - ((1-u)\beta_F N_F + (1-u)\beta_I N_I)S - \mu S \quad (1a) $$

$$ \dot{I} = ((1-u)\beta_F N_F + (1-u)\beta_I N_I)S - \mu I - \rho N_I \frac{I}{\alpha + I} \quad (1b) $$

$$ \dot{N}_F = -(1-u)c_\beta\beta_F(S+\eta I)N_F + \left(\gamma + \gamma \frac{N_I}{K I}\right) N_I - \nu_F N_F \quad (1c) $$

$$ \begin{aligned} \dot{N}_I = (1-u)c_\beta\beta_F(S+\eta I)N_F & - \left(\gamma + \gamma \frac{N_I}{K I}\right) N_I \\ & + \left(r + c_\rho \rho \frac{I}{\alpha + I}\right) N_I \left(1 - \frac{N_I}{K I}\right) - \nu_I N_I \end{aligned} \quad (1d) $$

As they represent biological quantities, more precisely root biomasses and nematode population sizes, the state variables should be non-negative. Moreover, mathematically, system (1) is not defined when the infected root biomass $I = 0$. Biologically, there are no infected roots without infesting nematodes and vice versa, so both $I$ and $N_I$ should be positive. System (1) is then defined in

$$ \Delta = \mathbb{R}_{\ge 0} \times \mathbb{R}_{> 0} \times \mathbb{R}_{\ge 0} \times \mathbb{R}_{> 0}. $$

The parameters used in model (1) are positive and described in Table 1. The parameter values used for the numerical simulations were partly based on the literature [7, 24], partly set to realistic values.

## 2.2 Basic properties of the epidemiological model

For system (1) to be mathematically and biologically meaningful, all state variables should remain in $\Delta$ and bounded for initial conditions in $\Delta$. This is stated more precisely in the following proposition.

**Proposition 1.** System (1) admits a unique solution that stays in $\Delta$ at all time $t \ge 0$. Moreover, the bounded region

$$ \Omega = \{(S, I, N_F, N_I) \in \Delta, S + I \le A, N_F + N_I \le B\} $$

with $A = \max\{S(0) + I(0), \frac{\Delta}{\mu}\}$ and $B = \max\{N_F(0) + N_I(0), \frac{KA(r+c_\rho\rho)}{4\min\{\nu_F, \nu_I\}}\}$, is positively invariant for system (1).

*Proof.* To prove the positivity of system (1), we move to the $(S, I, N_F, \mathcal{X})$ system, which is equivalent on $\Delta$. We substitute $N_I$ by $\mathcal{X}I$ in the first three equations of system (1) and we replace the last equation with the dynamics of $\mathcal{X}$. Hence, we obtain the following system well defined in $\mathcal{X} = 0$
<table><caption>Table 1: Parameters of the full model (black and grey), given by equations (1), and the reduced model (black only), given by equations (3).</caption><thead><tr><th>Symbol</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>Λ</td><td>Recruitment of healthy roots</td><td>300 g.day<sup>-1</sup></td></tr><tr><td>β<sub>F</sub></td><td>Infection rate by free nematodes</td><td>10<sup>-7</sup> nematode<sup>-1</sup>.day<sup>-1</sup></td></tr><tr><td>β<sub>I</sub></td><td>Infection rate by infesting nematodes</td><td>10<sup>-7</sup> nematode<sup>-1</sup>.day<sup>-1</sup></td></tr><tr><td>η</td><td>Reinfection fraction</td><td>0.0002</td></tr><tr><td>c<sub>β</sub></td><td>Infection conversion factor</td><td>100 nematodes.g<sup>-1</sup></td></tr><tr><td>γ</td><td>Nematode migration rate</td><td>0.02 day<sup>-1</sup></td></tr><tr><td>ρ</td><td>Root consumption rate</td><td>2 × 10<sup>-4</sup> g.nematode<sup>-1</sup>.day<sup>-1</sup></td></tr><tr><td>α</td><td>Root consumption half-saturation</td><td>60 g</td></tr><tr><td>c<sub>ρ</sub></td><td>Root consumption conversion factor</td><td>400 nematodes.g<sup>-1</sup></td></tr><tr><td>r</td><td>Nematode intrinsic growth rate</td><td>0.15 day<sup>-1</sup></td></tr><tr><td>K</td><td>Nematode carrying capacity</td><td>1000 nematodes.g<sup>-1</sup></td></tr><tr><td>μ</td><td>Mortality rate of roots</td><td>0.05 day<sup>-1</sup></td></tr><tr><td>ν<sub>F</sub></td><td>Mortality rate of free nematodes</td><td>0.0495 day<sup>-1</sup></td></tr><tr><td>ν<sub>I</sub></td><td>Mortality rate of infesting nematodes</td><td>0.045 day<sup>-1</sup></td></tr></tbody></table>

$$ \dot{S} = \Lambda - ((1-u)\beta_F N_F + (1-u)\beta_I X I)S - \mu S \quad (2a) $$

$$ \dot{I} = ((1-u)\beta_F N_F + (1-u)\beta_I X I)S - \mu I - \rho X \frac{I^2}{\alpha + I} \quad (2b) $$

$$ \dot{N}_F = -(1-u)c_\beta\beta_F(S+\eta I)N_F + \left(\gamma + \gamma \frac{X}{K}\right) X I - \nu_F N_F \quad (2c) $$

$$ \dot{X} = \left( r + c_{\rho}\rho \frac{I}{\alpha + I} \right) X \left( 1 - \frac{X}{K} \right) - \left( \gamma + \gamma \frac{X}{K} \right) X - (\nu_I - \mu)X + (c_{\beta}\beta(S + \eta I) - \beta S X) \frac{N_F}{I} - \beta_I S X^2 + \rho X^2 \frac{I}{a + I}, \quad (2d) $$

Suppose that variable $Y_i = 0$ with $Y_i \in \{S, N_F, X\}$ at some time, while the other variables are non-negative. A quick analysis shows that in this case $\dot{Y}_i \ge 0$, so that $Y_i$ cannot become negative. In addition, suppose that $I \to 0$ at some time while the other variables are non-negative, then $\dot{I} \ge 0$ since $X \le \frac{B}{I}$. Hence, solutions of system (2) are non-negative and consequently, those of system (1).

To ensure that any solution remains in $\Delta$, we still need to prove that $I$ and $N_I$ cannot reach zero in finite time. For that, we first need to show the boundedness of the variables. Adding equations (1a) and (1b), and using the comparison theorem, we obtain

$$ \dot{S} + \dot{I} \le \Lambda - \mu(S + I) \quad \Rightarrow \quad S + I \le \mathcal{A} = \max \left\{ S(0) + I(0), \frac{\Lambda}{\mu} \right\}. $$

Adding equations (1c) and (1d), we have

$$ \dot{N}_F + \dot{N}_I = \left( r + c_{\rho}\rho \frac{I}{\alpha + I} \right) N_I \left( 1 - \frac{N_I}{K I} \right) - \nu_I N_I - \nu_F N_F. $$

Using the fact that $N_I(1 - \frac{N_I}{K I}) \le \frac{K I}{4}$ provided that $I \ge 0$, we have

$$ \dot{N}_F + \dot{N}_I \le \frac{K A}{4} (r + c_{\rho}\rho) - \min\{\nu_F, \nu_I\} (N_F + N_I). $$

Using the comparison theorem, we obtain

$$ N_F + N_I \le \mathcal{B} = \max \left\{ N_F(0) + N_I(0), \frac{K A (r + c_{\rho}\rho)}{4 \min\{\nu_F, \nu_I\}} \right\}. $$

Having shown that the variables are non negative and upper-bounded, provided that $I$ remains positive, we now need to show that $I$ and $N_I$ cannot reach zero in finite time. Using equation (1b), we obtain

$$ \dot{I}(t) \ge -\mu I - \frac{\rho}{\alpha} \mathcal{B} I, $$
which after integration gives

$$I(t) \geq I(0)e^{-(\mu+\frac{\epsilon}{\alpha}\mathcal{B})t} > 0.$$

So $I$ cannot reach zero in finite time and remains positive. Using equation (1d), we have

$$\begin{aligned} \dot{N}_I(t) & \geq \left(r - \nu_I - \gamma + c_\rho \rho \frac{I}{\alpha + I}\right) N_I - \frac{1}{KI} \left(r + \gamma + c_\rho \rho \frac{I}{\alpha + I}\right) N_I^2, \\ & \geq (r - \nu_I - \gamma) N_I - \frac{1}{KI} (r + \gamma + c_\rho \rho) N_I^2. \end{aligned}$$

The comparison theorem applied to this last inequality gives us

$$N_I(t) \geq \frac{N_I(0)}{\omega_1 N_I(0) + \omega_2} > 0,$$

where

$$\omega_1 = \int_0^t \frac{1}{KI(s)} (r + \gamma + c_\rho \rho) e^{(r-(\nu_I+\gamma))(t-s)} ds \quad \text{and} \quad \omega_2 = e^{(r-(\nu_I+\gamma))t}.$$

So $N_I(t)$ cannot reach zero in finite time either. Therefore, the right hand-side of system (1) is Lipschitz continuous along all solutions, which implies that this system admits a unique solution that stays in $\Delta$ at all time. $\square$

## 2.3 Model reduction: a slow-fast analysis

Considering that root infestation is fast compared to nematode development, we use the change of variables $N = N_F + N_I$ and Tikhonov's theorem on slow-fast systems to reduce the model dimension. Hence, we obtain the following nonlinear reduced system (see Appendix A for more details)

$$\dot{S} = \Lambda - (1-u)\beta NS - \mu S, \quad (3a)$$

$$\dot{I} = (1 - u)\beta NS - \mu I - \rho N \frac{I}{\alpha + I}, \quad (3b)$$

$$\dot{N} = \left( r + c_{\rho} \rho \frac{I}{\alpha + I} \right) N \left( 1 - \frac{N}{KI} \right) - \nu N, \quad (3c)$$

with simplified parameter notations $\beta_I = \beta$ and $\nu_I = \nu$. It is defined in

$$\Delta' = \mathbb{R}_{\ge 0} \times \mathbb{R}_{> 0} \times \mathbb{R}_{> 0}.$$

As for system (1), we assume that nematode intrinsic growth is higher than natural mortality ($r > \nu$).

**Proposition 2.** System (3) admits a unique solution that stays in $\Delta'$ at all time $t \geq 0$. Also the bounded region

$$\Omega' = \{(S, I, N) \in \Delta', S + I \leq \mathcal{A}, N \leq \mathcal{C}\},$$

where $\mathcal{A} = \max\{S(0) + I(0), \frac{\Delta}{\mu}\}$ and $\mathcal{C} = \max\{N(0), K(S(0) + I(0)), \frac{K\Delta}{\mu}\}$, is positively invariant for system (3).

*Proof.* The principle is the same as in Proposition 1. $\square$

**Remark 1.** From proposition 2, system (3) is mathematically and biologically well-posed. Then, it is sufficient to consider its dynamics in $\Omega'$.

## 2.4 Asymptotic analysis of the reduced model

As the reduced model in system (3) is a good approximation of the full epidemiological model in system (1), for simplicity reasons, we carry out the equilibrium study on this reduced model. The results can be extended to the full epidemiological model. In the following, we assume that $u = 0$, which corresponds to a case with no pest control, or a constant control with a lower infestation rate.
### 2.4.1 Existence of endemic equilibria

According to proposition 2, we can look for equilibria in $\Omega'$, where $I$ and $N$ are positive. Therefore, they would necessarily be endemic equilibria, in which the pest ($N$) and the plant ($S+I$) coexist. We obtain the following result.

**Proposition 3.**

1. There is no equilibrium for system (3) when $c_\rho = 0$ and $\mathcal{R} \le 1$, where $\mathcal{R} = \frac{\Lambda\beta K}{\mu^2}(1 - \frac{\nu}{r})$.

2. When $c_\rho = 0$ and $\mathcal{R} > 1$, system (3) admits a unique endemic equilibrium.

3. When $c_\rho \ne 0$, system (3) admits 0 or 2 or 4 endemic equilibria if $\mathcal{R} < 1$, and 1 or 3 endemic equilibria if $\mathcal{R} > 1$.

*Proof.* See Appendix B. $\square$

### 2.4.2 Asymptotic behavior towards the pest-free situation

System (3) is not defined if there is no pest ($N = I = 0$), so it does not formally admit a pest-free equilibrium (PFE). However, we are still interested in its behavior when $N$ and $I$ tend to zero. Thanks to the same state variable change as for system (2), replacing the nematode population size $N$ by the infestation density $\mathcal{X} = \frac{N}{I}$, we obtain the following system

$$ \dot{S} = \Lambda - \beta \mathcal{X}IS - \mu S, \quad (4a) $$

$$ \dot{I} = \beta \mathcal{X}IS - \mu I - \rho \mathcal{X} \frac{I^2}{\alpha + I}, \quad (4b) $$

$$ \dot{\mathcal{X}} = \left( r + c_{\rho} \rho \frac{I}{\alpha + I} \right) \mathcal{X} \left( 1 - \frac{\mathcal{X}}{K} \right) - (\nu - \mu) \mathcal{X} - \beta S \mathcal{X}^2 + \rho \frac{I}{\alpha + I} \mathcal{X}^2, \quad (4c) $$

which is defined when the infected root biomass $I = 0$.

Now, looking for PFE for system (4), i.e. equilibria that verify $I = 0$, we obtain

$$ \zeta^1 = \left( \frac{\Lambda}{\mu}, 0, 0 \right) \text{ and } \zeta^2 = \left( \frac{\Lambda}{\mu}, 0, \mathcal{X}^0 \right), \text{ with } \mathcal{X}^0 = K\mu \frac{r - \nu + \mu}{r\mu + K\beta\Lambda}. $$

To study the local stability of these equilibria, we check the sign of the eigenvalues of the Jacobian matrix at each equilibrium. The Jacobian matrix at the first equilibrium $\zeta_1$ is

$$ \mathcal{J}_1 = \begin{pmatrix} -\mu & 0 & 0 \\ 0 & -\mu & 0 \\ 0 & 0 & r - \nu + \mu \end{pmatrix}. $$

As we assume that $r > \nu$ (item 9 in Section 2), we have $r - \nu + \mu > 0$, so $\zeta_1$ is unstable.

The Jacobian matrix at the second equilibrium $\zeta_2$ gives

$$ \begin{aligned} \mathcal{J}_2 &= \begin{pmatrix} -\mu & -\beta \frac{\Lambda}{\mu} \mathcal{X}^0 & 0 \\ 0 & \beta \frac{\Lambda}{\mu} \mathcal{X}^0 - \mu & 0 \\ -\beta(\mathcal{X}^0)^2 & \frac{c_\rho \rho}{\alpha} \mathcal{X}^0 \left(1 - \frac{\mathcal{X}^0}{K}\right) + \frac{\rho}{\alpha}(\mathcal{X}^0)^2 & r \left(1 - 2\frac{\mathcal{X}^0}{K}\right) - (\nu - \mu) - 2\beta \frac{\Lambda}{\mu} \mathcal{X}^0 \end{pmatrix} \\ &= \begin{pmatrix} -\mu & -\beta \frac{\Lambda}{\mu} \mathcal{X}^0 & 0 \\ 0 & \frac{r\mu^2}{r\mu + K\beta\Lambda}(\mathcal{R} - 1) & 0 \\ -\beta(\mathcal{X}^0)^2 & \frac{c_\rho \rho}{\alpha} \mathcal{X}^0 \left(1 - \frac{\mathcal{X}^0}{K}\right) + \frac{\rho}{\alpha}(\mathcal{X}^0)^2 & -(r - \nu + \mu) \end{pmatrix} \end{aligned} $$

Hence, $\zeta_2$ is locally asymptotically stable if $\mathcal{R} < 1$ and unstable if $\mathcal{R} > 1$.

## 2.5 Bifurcation analysis of the reduced model

### 2.5.1 Forward bifurcation

In this subsection we study the behavior of system (3) when there is a forward bifurcation. For this, we use $\beta$ as bifurcation parameter, as its variations cover all $\mathcal{R}$ values. We choose the parameters given in Table 1. Initial conditions are

$$ x(0) = \left\{ S(0) = \frac{\Lambda}{\mu}, I(0) > 0, N(0) = 900 \right\}. \quad (5) $$
Figure 2 presents the bifurcation diagram and dynamics of infected roots for the forward bifurcation case. It also illustrates the existence of the endemic equilibria and their stability. Subplot (a) presents a forward transcritical bifurcation and Subplot (b) the dynamics of infected roots for parameters in Table 1.

For $\mathcal{R} = 8.4$ (dotted line), subplot (a) shows that system (3) has one endemic equilibrium which is stable, as well as an unstable pest-free situation, and subplot (b) shows the convergence of infected roots to their endemic value.

<figure><img src="image_3.png" /><figcaption>Figure 2: Bifurcation diagram and dynamics of infected roots for the forward bifurcation case. Subplot (a) presents the $(\mathcal{R}, I)$-bifurcation diagram of system (3). Subplot (b) shows the dynamics of infected roots $I^*$ for $\mathcal{R} = 8.4$. The underlying bifurcation parameter is $\beta$. All remaining parameters are chosen in Table 1. Initial conditions are chosen in (5) with $I_0 = 1500$ (solid curves) and $I_0 = 50$ (dashed curves).</figcaption></figure>

### 2.5.2 Backward bifurcation

In this subsection we study the behavior of system (3) when there is a backward bifurcation. We also use $\beta$ as bifurcation parameter. We choose the parameters given in Table 1 except $\Lambda = 10$ and $r = 0.05$. We obtain Figure 3, which illustrates the existence of the endemic equilibria and their stability. Subplot (a) presents a backward transcritical bifurcation and subplot (b), the dynamics of infected roots for $\beta = 1.25 \times 10^{-6}$, which yields $\mathcal{R} = 0.5$ (dotted line).

<figure><img src="image_4.png" /><figcaption>Figure 3: Bifurcation diagram and dynamics of infected roots for the backward and forward bifurcation case. Subplot (a) presents the $(\mathcal{R}, I)$-bifurcation diagram of system (3). Subplot (b) shows the dynamics of infected roots $I^*$ for $\beta = 1.25 \times 10^{-6}$, which yields $\mathcal{R} = 0.5$. The underlying bifurcation parameter is $\beta$. All remaining parameters are given in Table 1 except $\Lambda = 10$ and $r = 0.05$. Initial conditions are given in (5) with $I_0 = 12$ (dashed curves) and $I_0 = 40$ (solid curves).</figcaption></figure>

For $\mathcal{R} = 0.5$, Subplot (a) shows that system (3) has two endemic equilibria (EE), one of which is stable
and the other unstable. It also shows that the system has a pest-free situation (PFE) which is stable and coexists with the stable EE. Subplot (b) shows the convergence of infected roots to their stable endemic value (for $I_0 = 40$) or towards the pest-free situation (for $I_0 = 12$).

<figure><img src="image_5.png" /><figcaption>Figure 4: Extinction and persistence of nematodes in the banana plantation. Simulation of model (3) for $\mathcal{R} = 0.5$, using initial conditions (5), with $I_0 = 2$ (blue curves), $I_0 = 12$ (green curves), $I_0 = 30$ (black curves) and $I_0 = 40$ (red curves). All remaining parameters are given in Table 1 except $\Lambda = 10$ and $r = 0.05$.</figcaption></figure>

Figure 4 shows the convergence of system (3) towards the PFE (blue and green curves) or towards the stable endemic equilibrium (black and red curves), according to different initial values of $I_0$. At the EE, subplot (a) shows that healthy roots are notably reduced, as expected with the strong increase of the nematode population. It is interesting to note in subplot (d) that the infestation density tends to a positive value, which is higher for the EE than for the PFE.

# 3 Optimization and optimal control strategies

## 3.1 Statement of the optimal control problem (OCP)

Applying biostimulants can reduce nematode infestation [26] and hence its negative impact on banana yield. The effect of this control method is introduced in system (3) through control variable $u$. A proxy of the banana yield is given by the integral of the healthy root biomass during a cropping season [23]. The objective is hence to maximize this yield minus the control costs, so as to maximize profit. Moreover, to ensure a reasonable yield for the next season, we penalize the final infected root biomass. Indeed, banana plants reproduce through suckers, which are offshoots of the parent plant. Denoting by $t_f$ the duration of the cropping season, we obtain the following objective function

$$ \begin{aligned} \mathcal{J}(u) &= \int_{0}^{t_f} (B_u u(t) - B_S S(t)) dt + B_I I(t_f), \\ &= \mathcal{J}_u(u) + \mathcal{J}_S(u) + \mathcal{J}_I(u), \end{aligned} \quad (6) $$

where,

$$ \mathcal{J}_u(u) = \int_{0}^{t_f} B_u u(t) dt, \quad \mathcal{J}_S(u) = -B_S \int_{0}^{t_f} S(t) dt \quad \text{and} \quad \mathcal{J}_I(u) = B_I I(t_f). $$
This type of problem is called a Bolza optimal control [20]. $\mathcal{J}_u$ is the total cost of the biostimulant, with $B_u$ its unit price. $\mathcal{J}_S$ represents the opposite of the cropping yield, with $B_S$ its unit sale price. Finally, $\mathcal{J}_I$ is the cost associated with infected roots remaining after harvest, weighted by $B_I$, which reflects their economic impact. Hence, $-\mathcal{J}$ represents the profit over the entire cultivation period, penalized by the infestation at the end of the cropping season.

The set of admissible controls is defined as follows

$$ \mathcal{U} = \{u \in L^{\infty}([0, t_f]) : 0 \le u(t) \le u_{\max}, \forall t \in [0, t_f]\}. \quad (7) $$

Hence, the main optimal control problem (OCP) is stated below.

**Problem 1.** Find an admissible control $u^* \in \mathcal{U}$ defined in (7) and the corresponding state variable $x^* = (S^*, I^*, N^*)$ that solve the following OCP:

$$ \left\{ \begin{array}{l} \min_{u \in \mathcal{U}} \mathcal{J}(u), \\ \text{with } \mathcal{J}(u) \text{ given by (6).} \end{array} \right. $$

**Remark 2.** The existence of an optimal control for the Problem 1 is guaranteed by Fleming and Rishel [10] since $\mathcal{U}$ is a closed convex set and $S$ is bounded.

## 3.2 Necessary optimal conditions: application of Pontryagin's Maximum Principle

In this subsection, we establish the first-order necessary optimality conditions using the Pontryagin's Maximum Principle (PMP) [20]. These optimality conditions are used to compute the optimal control $u^*$ of Problem 1. This principle transforms the optimization problem into a problem consisting in determining the minimum with respect to $u^*$ of the corresponding Hamiltonian. Let $\lambda = (\lambda_1, \lambda_2, \lambda_3)$ be the adjoint variables, $\lambda^0 \in \mathbb{R}$. The Hamiltonian of Problem 1 is defined as follows

$$ \begin{aligned} \mathcal{H}(x, u, \lambda, \lambda^0) = & \lambda^0 [B_u u - B_S S] + \lambda_1 [\Lambda - (1-u)\beta N S - \mu S] \\ & + \lambda_2 \left[ (1-u)\beta N S - \mu I - \rho N \frac{I}{\alpha + I} \right] \\ & + \lambda_3 \left[ \left( r + c_\rho \rho \frac{I}{\alpha + I} \right) N \left( 1 - \frac{N}{K I} \right) - \nu N \right]. \end{aligned} \quad (8) $$

Let $u^*(\cdot)$ be an optimal control and $x^* = (S^*, I^*, N^*)$ the associated trajectory on $[0, t_f]$ starting from $x_0 = (S_0, I_0, N_0)$. Then, according to the PMP, we have the following conditions

* There exists $\lambda^0 \ge 0$ and an absolutely continuous non-zero vector mapping $\lambda: [0, t_f] \to \mathbb{R}^3$ called the adjoint vector, which satisfies the following adjoint equation $\dot{\lambda} = -\frac{\partial \mathcal{H}}{\partial x}(x(t), u(t), \lambda(t), \lambda^0)$ $\forall t \in [0, t_f]$, i.e

$$ \left\{ \begin{array}{l} \dot{\lambda}_1 = \lambda_0 B_S + (1-u)\beta N(\lambda_1 - \lambda_2) + \mu\lambda_1, \\ \dot{\lambda}_2 = \left(\mu + \rho\alpha\frac{N}{(\alpha+I)^2}\right)\lambda_2 - \left(c_\rho\rho\alpha\frac{N}{(\alpha+I)^2}\left(1-\frac{N}{KI}\right) + \frac{1}{K}\left(\frac{N}{I}\right)^2\right)(r+c_\rho\rho\frac{I}{\alpha+I})\lambda_3, \\ \dot{\lambda}_3 = (1-u)\beta S(\lambda_1 - \lambda_2) + \rho\frac{I}{\alpha+I}\lambda_2 - \left((r+c_\rho\rho\frac{I}{\alpha+I})\left(1-2\frac{N}{KI}\right) - \nu\right)\lambda_3. \end{array} \right. $$

* The transversality conditions satisfy $\lambda(t_f) = \nabla\psi(x(T))\lambda^0$, that is

$$ \lambda_1(t_f) = \lambda_3(t_f) = 0 \quad \text{and} \quad \lambda_2(t_f) = B_I \lambda^0. \quad (9) $$

* Control $u(\cdot)$ satisfies the minimization condition

$$ u^*(t) \in \arg \min_{\omega \in \mathcal{U}} \mathcal{H}(x(t), \omega(t), \lambda(t), \lambda^0), \quad \forall t \in [0, t_f]. \quad (10) $$

If $\lambda^0 = 0$, we have $\lambda(t_f) = 0$. This leads to $\dot{\lambda}_i(t_f) = 0$ for $i = 1, 2, 3$, which implies that $\lambda_i(t_f) = 0$ for all $t \in [0, t_f]$ and therefore contradicts the first condition of PMP. Hence $\lambda^0 \ne 0$, and we choose $\lambda^0 = 1$. Now, as $\mathcal{H}$ is linear in $u$, its minimization depends on the sign of

$$ \partial_u \mathcal{H}(x(t), u(t), \lambda(t)) = B_u - \beta N S(\lambda_2 - \lambda_1). $$

By introducing the switching function

$$ \phi(t) = B_u - \beta N S(\lambda_2 - \lambda_1), \quad (11) $$

we characterise the optimal control $u^*$ as follows

$$ u^*(t) = \begin{cases} 0 & \text{if } \phi(t) > 0, \\ \text{undefined} & \text{if } \phi(t) = 0, \\ u_{\max} & \text{if } \phi(t) < 0. \end{cases} $$
**Remark 3.** At the end of the season, $\phi(t_f) = B_u - B_I \beta N(t_f) S(t_f)$. Hence, if $B_I = 0$, then $\phi(t_f) > 0$, which implies $u^*(t_f) = 0$.

The switching function $\phi(t)$ may vanish at isolated time points $t \in (0, t_f)$, indicating a switch between bang-bang regimes. If instead $\phi(t)$ vanishes identically on a nonempty interval $[t_1, t_2] \subset [0, t_f]$ with $t_1 < t_2$, a singular arc arises, associated with a singular control $u_{\text{sing}}^*$. To characterize a possible first-order singular arc, one follows standard geometric control theory. If $\phi(t) = 0$ on $[t_1, t_2]$, then necessarily $\dot{\phi}(t) = 0$ on this interval, and the control variable does not appear explicitly. Differentiating once more yields

$$\ddot{\phi}(t) = \mathcal{D}(x, \lambda) + (1-u)\mathcal{F}(x, \lambda) = 0,$$

where the analytical expressions of $\mathcal{D}$ and $\mathcal{F}$ are provided in Appendix C. This leads to the formal expression of a first-order singular control,

$$u_{\text{sing}}^*(t) = 1 + \frac{\mathcal{D}(x, \lambda)}{\mathcal{F}(x, \lambda)}. \qquad (12)$$

It should be emphasized that, due to the complexity of $\mathcal{D}$ and $\mathcal{F}$ (Appendix C), there is no a priori guarantee that $u_{\text{sing}}^*(t)$ remains within the admissible bounds $[0, u_{\text{max}}]$. Although higher-order singular arcs cannot be excluded in theory, their analysis would be technically demanding: for instance, second-order singularities would require differentiating the switching function up to $\phi^{(4)}(t)$. As confirmed later by the numerical results, however, singular arcs do not play a significant role in the optimal strategy. We therefore restrict attention here in the theoretical analysis to the first-order case only. Finally, a candidate singular arc can only be optimal if it satisfies the second-order optimality conditions, namely the Legendre-Clebsch condition [21]:

$$(-1)^h \partial_u \left[ \frac{d}{dt^2} \partial_u \mathcal{H}(x(\cdot), u(\cdot), \lambda(\cdot)) \right] \geq 0. \qquad (13)$$

In our problem this reduces to

$$(-1)^1 \partial_u [\ddot{\phi}(t)] = \mathcal{F}(x, \lambda). \qquad (14)$$

Since the sign of $\mathcal{F}(x, \lambda)$ depends on the adjoint variables, it cannot be determined analytically, and the Legendre-Clebsch condition can only be verified numerically if a singular arc is detected.

It is worth noting that singular arcs are not uncommon in biological and ecological optimal control problems, where both first- and second-order cases have been shown to play a significant role (see, e.g., [16] for a second-order singular arc and [9] for a first-order one). However, such regimes are in general difficult to implement in practice, as they are open-loop controls which typically depend on the co-state variables. By contrast, when the switching function vanishes only at isolated time points and then changes sign, the resulting strategy is of the bang-bang type (see, e.g., [18]). Bang-bang structures are easier to interpret and implement in practice, since they rely on clearly defined maximal and minimal control actions. In the next section, we will show that for the parameter set considered, the optimal solution indeed exhibits a bang-bang quasi-periodic structure, although different parameter values could in principle lead to alternative control patterns. In cases where singular arcs arise under other conditions or datasets, the analytical expression and optimality test derived above provide a framework for assessing their relevance.

## 3.3 Direct optimization

To investigate the structure of the optimal solutions, we solve our OCP using a direct optimization method. Direct methods are numerical approaches that solve the optimal control problem numerically without a priori applying the PMP. These methods involve selecting a discretization for the state and control variables, as well as a numerical scheme to discretize the control system and the cost function. This process transforms the discretized optimal control problem into a family of nonlinear optimization problems in finite dimensions [6]. Bocop is such a software dedicated to solving optimization and optimal control problems based on the IPOPT optimization library [5].

To obtain the numerical solution of system 3, we implement and solve this optimal control problem on Bocop, choosing the Lobatto discretisation method (implicit, 4 stages, order 6), 1000 time steps, and a $10^{-20}$ tolerance. Model parameter values are given in Table 1 and initial conditions in

$$x(0) = \{S(0) = 10^2, I(0) = 10^2, N(0) = 10^3\} \qquad (15)$$

We solve the OCP on finite time horizon $t_f = 300$ days, which corresponds to a cropping season, with $u_{\text{max}} = 0.8$ (unless otherwise specified) and the following weight parameters

$$B_u = 7800, B_S = 2, B_I = 180 \text{ or } 0. \qquad (16)$$
<figure><img src="image_6.png" /><figcaption>Figure 5: Optimal control with penalty. Simulation of system (3) without pest (dashed gray line), without control (black curves) and with control (red curves) using parameters given in Table 1 and equations (15)-(16), with $B_I = 180$. We obtain $\mathcal{J}(u^*) = -2.023 \times 10^6$, $\mathcal{J}_u(u^*) = 1.014 \times 10^6$, $\mathcal{J}_S(u^*) = -3.066 \times 10^6$ and $\mathcal{J}_I(u^*) = 0.029 \times 10^6$.</figcaption></figure>

The numerical results thus found are presented in Figure 5.

Figure 5a) shows that in the absence of control (black curves), healthy root production declines rapidly after initial growth, dropping by more than half compared with the pest-free scenario (dashed gray line). For the first 50 days or so, the red curve, which represents the controlled scenario, follows a similar trajectory to that observed in the absence of control, but then begins to oscillate throughout the season, remaining at a high level to maximise the yield. Figures 5b) and 5c) illustrate the dynamics of infected roots and the total nematode population. In the absence of control (black curves), the number of infected roots and the total nematode population initially increase rapidly and then stabilise at a high level. In contrast, with control (red curve), infected roots and the total nematode population oscillate at much lower levels throughout the season, demonstrating that optimal control limits the spread of infection and reduces nematode proliferation. Finally, Figure 5d) shows the optimal control scheme applied. The control oscillates regularly between two levels: 0 and $u_{max} = 0.8$. This type of “bang-bang” control means that at certain times no control effort is applied ($u = 0$), while at other times maximum control ($u = 0.8$) is applied. This regular alternation makes it possible to adjust infestation management effectively throughout the season. The penalized profit $-\mathcal{J}(u^*) = 2.120 \times 10^6$ is nearly double the profit obtained in the uncontrolled case $-\mathcal{J}(0) = 1.681 \times 10^6$.

In summary, this bang-bang control strategy, where control is applied quasi-periodically (approximately every 10 days), maintains high levels of healthy roots, while limiting infected roots and the nematode population. This oscillation mechanism appears to be a particular feature of our optimal control problem, leading to an effective balance between crop production, pest control and cost efficiency.

Figure 6 shows that when the final penalty on infected roots is lifted ($B_I = 0$), the infected roots and the total nematode population increase sharply at the end, while healthy roots decrease. This indicates that, without an additional constraint at the end of the cropping season, the optimal solution consists in stopping the control towards the end of the season, as it reduces the control costs with a limited decrease of the yield (corresponding to the integral of the healthy root biomass), compromising root health for the next cropping
<figure><img src="image_7.png" /></figure>

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

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

<figure><img src="image_10.png" /><figcaption>Figure 6: Optimal control without penalty. Simulation of system (3) without pest (dashed gray line), without control (black curves) and with control (red curves) using parameters given in Table 1 and equations (15)-(16), with $B_I = 0$. We obtain $\mathcal{J}(u^*) = -2.120 \times 10^6$, $\mathcal{J}_u(u^*) = 0.964 \times 10^6$, $\mathcal{J}_S(u^*) = -3.084 \times 10^6$ and $\mathcal{J}_I(u^*) = 0$.</figcaption></figure>

season. The introduction of this penalty promotes a balance that keeps infestation at controllable levels and preserves yield in the long term.

When $u_{\max}$ is reduced to 0.5, the control behavior becomes simpler and less reactive. The options are now limited to three scenarios, whatever the cost of control: either it remains constant at “bang(0)” when the control cost $B_u$ is high, or constant at 0.5 “bang(0.5)” when the control cost $B_u$ is low, or it oscillates according to the pattern “bang(0)–bang(0.5)–bang(0)” when the control cost $B_u$ is intermediate, as illustrated in Figure 7. These simpler, more predictable strategies are due to the reduced and relatively poor efficiency of the control method in this case. As depicted in Figure 7, the healthy roots are notably lower than in the pest-free scenario, the infestation is conversely higher.

**Remark 4.** *From our numerical resolution of the optimal control problem using direct methods (Bocop), we draw the following observations.*

1. The optimal control systematically exhibits a bang-bang structure.

2. The resulting optimal trajectories, both for the states and the control, are quasi-periodic.

3. Despite a numerical exploration, we did not observe the presence of singular arcs in this problem.

## 3.4 Suboptimal controls

In this section, we conduct a detailed numerical comparison between the optimal bang-bang quasi-periodic solution and candidate suboptimal strategies, including an ideal bang-bang control for its effectiveness and regularity.
<figure><img src="image_11.png" /><figcaption>Figure 7: Optimal control with reduced on control upper bound. Simulation of system (3) without pest (dashed gray line), without control (black curves) and with control (red curves) using parameters given in Table 1 and equations (15)-(16), with $B_I = 0$ and $u_{max} = 0.5$. We obtain $\mathcal{J}(u^*) = -1.813 \times 10^6$, $\mathcal{J}_u(u^*) = 0.871 \times 10^6$, $\mathcal{J}_S(u^*) = -2.684 \times 10^6$ and $\mathcal{J}_I(u^*) = 0$.</figcaption></figure>

### 3.4.1 Suboptimal constant controls

In order to further understand the problem, we will examine if efficient constant controls can be identified despite the quasi-periodic structure of the OCP solution. We consider three approaches detailed below.

1. Firstly, we solve the static version of the problem at steady state, the so-called static-OCP. Indeed, when the time horizon is sufficiently large, the optimal control often tends towards this solution [27]. Hence, we solve

$$ \min_{u \in [0, u_{\max}]} \bar{\mathcal{J}}(u) = B_u u - B_S S, \qquad (17) $$

under the constraint that the system is at equilibrium

$$ \left\{ \begin{array}{l} \Lambda - (1-u)\beta NS - \mu S = 0, \\ (1-u)\beta NS - \mu I - \rho N \frac{I}{a+I} = 0, \\ (r + c_{\rho} \rho \frac{I}{a+I}) N (1 - \frac{N}{KI}) - \nu N = 0. \end{array} \right. \qquad (18) $$
The static-OCP then reduces to (see Appendix D for more details)

$$
\begin{equation}
\left\{
\begin{array}{l}
\displaystyle \min_{u \in [0, u_{\max}]} \bar{\mathcal{J}}(u) = B_u u - B_S S, \\[1em]
\text{under the constraints:} \\[1em]
\displaystyle S = \frac{\Lambda}{\beta(1-u)N + \mu}, \\[1em]
\displaystyle N = KI \left( \frac{(r+c_{\rho}\rho - \nu)I + (r-\nu)\alpha}{(r+c_{\rho}\rho)I + r\alpha} \right), \\[1em]
\displaystyle a_4 I^4 + a_3 I^3 + a_2 I^2 + a_1 I + a_0 = 0.
\end{array}
\right.
\tag{19}
\end{equation}
$$

Under the assumption that Problem 19 has a unique solution $(\bar{x}, u_s) = (\bar{S}, \bar{I}, \bar{N}, u_s)$, the definition of the static-OCP is extended to include the co-state variables using the Hamiltonian equation (8). Indeed, according to the Lagrange multipliers rule, there exists $\bar{\lambda} = (\bar{\lambda}_1, \bar{\lambda}_2, \bar{\lambda}_3)$ such that

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial \mathcal{H}}{\partial \lambda_i}(\bar{x}, u_s, \bar{\lambda}, 1) = 0, \\[2ex]
\displaystyle -\frac{\partial \mathcal{H}}{\partial x_i}(\bar{x}, u_s, \bar{\lambda}, 1) = 0, \\[2ex]
\displaystyle \frac{\partial \mathcal{H}}{\partial u}(\bar{x}, u_s, \bar{\lambda}, 1) = 0.
\end{array}
\right.
\tag{20}
$$

The solution of the static-OCP, $u_s$, is given by Bocop.

2. Secondly, we explore whether the average of the bang-bang periodic optimal control, denoted by $u_a$, provides an appropriate suboptimal solution. Indeed, a singular control often appears as an average of a bang-bang control, an intermediate means to avoid abrupt phases with an instantaneous change, and more practical to implement.

3. Finally, we solve the constant control optimization problem, which corresponds to the control value in $[0, u_{\max}]$ denoted by $u_c$ for which the minimum of the objective function is reached along the dynamic system.

Applying the value of $u_s$ found in Bocop, as well as those of $u_a$ and $u_c$ for $B_u = 7800$, $B_S = 2$, $B_I = 0$ and $u_{\max} = 0.8$, we compare the profits obtained when solving the OCP and with the three constant control strategies. By definition, the OCP solution performs globally better than the other strategies. However, more detailed comparisons can be drawn from the dynamics in Figure 8 and the corresponding objective function values in Table 2: the yield is proportional to the integral of healthy roots in Figure 8(a) and is given by $-\mathcal{J}_S$ in Table 2; the control cost is proportional to the integral of $u$ in Figure 8(b) and is given by $\mathcal{J}_u$ in Table 2.

<table><caption>Table 2: Objective function values for the OCP and the suboptimal strategies. Components of function $\mathcal{J}$ defined in (6) are detailed: $\mathcal{J}_u$ corresponds to the control cost, $-\mathcal{J}_S$ to the yield. Constant control and periodic strategies are illustrated in Figure 8. $B_u = 7800$, $B_S = 2$, $B_I = 0$ and $u_{\max} = 0.8$, parameters are given in Table 1 and initial conditions in (15).</caption>
  <thead>
    <tr>
      <th rowspan="2">Strategy</th>
      <th rowspan="2">Control u</th>
      <th colspan="3">10<sup>-6</sup> criterion values</th>
    </tr>
    <tr>
      <th>J<sub>u</sub></th>
      <th>J<sub>S</sub></th>
      <th>J</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>OCP</td>
      <td>u<sup>*</sup> varying</td>
      <td>0.964</td>
      <td>-3.084</td>
      <td>-2.120</td>
    </tr>
    <tr>
      <td>periodic-OCP</td>
      <td>u<sub>p</sub> varying</td>
      <td>0.930</td>
      <td>-3.047</td>
      <td>-2.117</td>
    </tr>
    <tr>
      <td>constant-OCP</td>
      <td>u<sub>c</sub> 0.620</td>
      <td>1.451</td>
      <td>-3.155</td>
      <td>-1.704</td>
    </tr>
    <tr>
      <td>static-OCP</td>
      <td>u<sub>s</sub> 0.662</td>
      <td>1.549</td>
      <td>-3.242</td>
      <td>-1.685</td>
    </tr>
    <tr>
      <td>average-OCP</td>
      <td>u<sub>a</sub> 0.412</td>
      <td>0.963</td>
      <td>-2.521</td>
      <td>-1.558</td>
    </tr>
  </tbody>
</table>

The static-OCP and constant-OCP solutions are comparable, their control values being close. They both maintain very good yields, even better than the OCP, but at very high control costs. Conversely, the control cost of the average-OCP and the OCP are the same (by definition of the average-OCP), but the yield is strongly degraded for the average-OCP. As a consequence, the three constant control strategies are markedly less efficient than the OCP solution, the average-OCP being the worse.

Hence, our numerical results suggest the presence of an overyielding-like feature induced by the OCP. In the ecological literature, overyielding usually refers to situations where a periodic environment achieves higher
<figure><img src="image_12.png" /><figcaption>Figure 8: Comparison between optimal and suboptimal strategies. Simulation of system (3): without pest (dashed gray line); without control (black curves); with control, solution of OCP (red curves), static-OCP (green curves), average-OCP (yellow curves), constant-OCP (magenta curves), periodic-OCP (dashed blue curves). Parameters are given in Table 1 and equations (15)-(16), with $B_I = 0$. The different objective function values are given in Table 2.</figcaption></figure>

performance than a constant one (for instance, day-night light cycles compared to constant illumination). In our case, the periodicity is not externally imposed but instead generated endogenously by the optimal bang-bang control input, which consistently outperforms constant or averaged strategies.

In view of the results of this subsection, it seems appropriate to consider whether a periodic control might yield better performance than constant controls.

### 3.4.2 Suboptimal periodic control

In this subsection, we consider a periodic control that is similar to the optimal solution of Figure 6, that is starting and ending with $u = 0$, and presenting the same number of upper-bangs at $u_{max}$ as the optimal control, through four parameters, $\tau_0$, $\tau_1$, $\tau_2$, $\tau_3$ as described below.

* If $t < \tau_0$, $u = 0$.

* Once $t \ge \tau_0$, the control alternates periodically between:

- $u = u_{max}$ for a duration $\tau_1$ (upper-bang);

- $u = 0$ for a duration $\tau_2$ (lower-bang);

- the cycle ends after $n$ upper-bangs and $n - 1$ lower-bangs.

* After this cycle, $u = 0$ for a duration $\tau_3$.

* Parameters are constrained by: $\tau_0 + n\tau_1 + (n - 1)\tau_2 + \tau_3 = t_f$.

Hence, the problem is mathematically formulated as follows

$$ \begin{cases} \min_{\tau_0, \tau_1, \tau_2, \tau_3} \mathcal{J}(u), \\ \text{with } \mathcal{J}(u) \text{ defined in equation (6).} \end{cases} \qquad (21) $$

To solve this problem, we use the global Differential Evolution method [22], available in the python SciPy library (scipy.optimize.differential_evolution).

For $n = 5$ upper-bangs and $t_f = 300$ days, as in Figure 6 (and the same parameter values), using [0, 300] for the parameter bounds, we obtain the following optimal values [$\tau_0, \tau_1, \tau_2, \tau_3$] = [51.7, 29.9, 16.3, 33.7] (in days). The suboptimal solution is illustrated in Figure 8 and the objective function value is given in Table 2.

As could be expected, the solutions of the OCP and periodic-OCP are very close, both in terms of control and of objective function value. Moreover, the periodic solution is easier for farmers to implement. This approach therefore represents an economically efficient and accessible solution.
# 4 Conclusions

In this paper, we developed a mathematical model based on the life cycle of the nematode *Radopholus similis* in order to study its infestation dynamics in the banana roots during a cropping season. The originality of the model is that it takes into account infected roots and infesting nematodes, which allows variable nematodes densities in the roots. After consideration that the root infestation is very fast compared to nematode development, we reduced the dimension of our original system using the slow-fast system theory. In the theoretical analysis of the reduced system, we derived a threshold $\mathcal{R}$ which classically determines the stability of the pest-free situation and identified two cases of bifurcations. The first case is a forward bifurcation which corresponds to a stable pest-free situation when $\mathcal{R} < 1$, and a stable endemic equilibrium when $\mathcal{R} > 1$. The second case is a backward bifurcation in which, on top of the stable pest-free situation, there is a stable endemic equilibrium when $\mathcal{R} < 1$. This means that $\mathcal{R} < 1$ does not ensure pest control. Burrowing nematodes being striving pests, realistic parameter values correspond to $\mathcal{R} > 1$, a situation in which pest control is necessary.

To reduce nematode damages in the plantation, we designed an optimal control problem based on the use of biostimulants, which reduce nematode infestation. The objective of this control problem was to maximise the profit during a cropping season, while minimising the infected roots at the end of the season, to ensure reasonable yields for the next season. We showed that an optimal control exists, which we characterised using Pontryagin's Maximum Principle. In our numerical simulations, we obtained bang-bang quasi-periodic solutions that are consistent with our theoretical results. We also showed that the biostimulant application in the plantation not only effectively reduced the dynamics of infected roots and the total nematode population by more than half, but also doubled the yield and the penalized profit compared to the uncontrolled case. The application of suboptimal constant controls was notably less efficient in terms of yield and profit, but the periodic control was very similar to the optimal control, with a simpler structure. These results highlighted an overyielding effect.

# Data Availability

All parameter values used for simulations are available within the paper. The python codes used to generate all the figures presented in the paper, as well as the Bocop codes obtained from solving the optimal control problems, are available in the public repository: (https://github.com/fkemayou/Myproject).

# Acknowledgements

This work is supported by EPITAG, an Inria Associate team part of the LIRIMA, supported by INRAE (https://team.inria.fr/epitag/). The authors also thank Jean Baptiste Caillau and Pierre Martinon for helpful discussions about the Bocop software.

# Appendix A. Model reduction

Banana roots being the main resource of the nematodes (the ratio between the population density in the soil and in the roots expressed as number of nematodes per gram is generally less than 1/100 [7]), we can then assume that root infestation by free nematodes occurs at a scale faster than nematode demography. In other words, the corresponding epidemiological parameter is superior to the demographic parameters. One can then apply singular perturbation theory to approximate this system by a lower dimensional system. The theorem supporting our following analysis is given in Khalil ([17]; Theorem 11.1; page 434). We hence introduce a scale parameter such that

$$c_{\beta} = \frac{c'_{\beta}}{\varepsilon}, \quad \text{with} \quad 0 < \varepsilon \ll 1. \qquad (22)$$

With the following change in variables: $N = N_F + N_I$, system (1) and (22) can be rewritten as

$$\dot{S} = \Lambda - ((1-u)\beta_F N_F + (1-u)\beta_I(N - N_F))S - \mu S \qquad (23a)$$

$$\dot{I} = ((1-u)\beta_F N_F + (1-u)\beta_I(N - N_F))S - \mu I - \rho(N - N_F) \frac{I}{\alpha + I} \qquad (23b)$$

$$\dot{N} = \left( r + c_{\rho} \rho \frac{I}{\alpha + I} \right) (N - N_F) \left( 1 - \frac{(N - N_F)}{KI} \right) - \nu_I (N - N_F) - \nu_F N_F \qquad (23c)$$

$$\varepsilon \dot{N}_F = -(1-u)c'_{\beta} \beta_F (S + \eta I) N_F + \varepsilon \left( \gamma + \gamma \frac{(N - N_F)}{KI} \right) (N - N_F) - \varepsilon \nu_F N_F \qquad (23d)$$
Setting $x = (S, I, N)$ and $z = N_F$, this system takes the form of an autonomous singular perturbation model:

$$ \dot{x} = f(x, z, \varepsilon) \qquad (24a) $$

$$ \varepsilon \dot{z} = g(x, z, \varepsilon) \qquad (24b) $$

where, $x$ are the slow variable and $z$ the fast variables. The fast system has a unique non-negative quasi steady state when $\varepsilon = 0$

$$ z^* = N_F^* = h(x) = 0. \qquad (25) $$

By substituting equation (25) in the slow system, we obtain the following nonlinear reduced system

$$ \dot{S} = \Lambda - (1 - u)\beta_I N S - \mu S \qquad (26a) $$

$$ \dot{I} = (1 - u)\beta_I N S - \mu I - \rho N \frac{I}{\alpha + I} \qquad (26b) $$

$$ \dot{N} = \left( r + c_{\rho} \rho \frac{I}{\alpha + I} \right) N \left( 1 - \frac{N}{K I} \right) - \nu_I N \qquad (26c) $$

The boundary-layer system is given by

$$ \frac{dy}{d\tau} = g(t, x, y + h(t, x), 0) = -(1-u)c'_{\beta}\beta_{F}(S + \eta I)y, \quad \text{with} \quad dt = \varepsilon d\tau. \qquad (27) $$

As (27) is linear, the origin is trivially exponentially stable uniformly in $(t, x)$.

Moreover $f, g, h$, the Jacobian of the reduced model (3) and $[\partial g(t, x, z, 0)/\partial z]$ are smooth functions. Hence, we conclude that when $\varepsilon$ is “small enough”, the reduced model (3) is a good approximation of the original system (1), with $z = z^* = h(x)$ given by equation (25).

To illustrate the effect of time scale separation, in Figure(9), we simulate and compare the solutions of the reduced system (3) and the complete system (23), that we have shown to be equivalent to system (1)-(22), using different epsilon values, namely $\varepsilon = \{1/10, 1/100, 1/1000\}$ and initial conditions in

$$ x_1(0) = \{S(0) = 6000, I(0) = 500, N(0) = 10^5, N_F(0) = 5 \times 10^4\}. \qquad (28) $$

When $\varepsilon = 1/10$, the dynamics of both systems differ, especially for the fast variable $N_F$: its quasi steady state value $N_F^* = 0$ is not a good approximation of its dynamics in the complet model. When epsilon is smaller, the dynamics of the slow variables $S, I$ and $N$ are perfectly approximated by the reduced system. When $\varepsilon = 1/100$, the dynamics of the fast variable $N_F$ is fairly well approximated and very well approximated when $\varepsilon = 1/1000$.

## Appendix B. Existence of endemic equilibria

We present here the proof of Proposition 3. Let us first consider the case when $c_\rho = 0$. Setting the right hand side of system (3) to 0, equations (3c) and (3a) give

$$ \bar{N} = K \left(1 - \frac{\nu}{r}\right) \bar{I} \quad \text{and} \quad \bar{S} = \frac{\Lambda}{\beta K \left(1 - \frac{\nu}{r}\right) \bar{I} + \mu}. \qquad (29) $$

By substituting equation (29) in equation (3b), we obtain

$$ m_2(\bar{I})^2 + m_1 \bar{I} + m_0 = 0, \qquad (30) $$

where

$$ \left\{ \begin{array}{l} m_2 = \mu\beta K \left(1 - \frac{\nu}{r}\right) + \rho\beta K^2 \left(1 - \frac{\nu}{r}\right)^2, \\ m_1 = \mu^2(1 - \mathcal{R}) + \mu\alpha\beta K \left(1 - \frac{\nu}{r}\right) + \rho\mu K \left(1 - \frac{\nu}{r}\right), \\ m_0 = \alpha\mu^2(1 - \mathcal{R}). \end{array} \right. $$

**Remark 5.** We assume in the model that $r > \nu$. So we always have $m_2 > 0$. Moreover, the sign of $m_0$ depends on $\mathcal{R}$ and $m_1 > 0$ if $\mathcal{R} < 1$.

Using the Descartes' rule for polynomial (30), we obtain the number of positive solutions, as shown in in Table 3. Moreover, there is no positive solution when $\mathcal{R} = 1$, as $m_0 = 0$ and $m_1 > 0$.

Now, let us compute this equilibrium in the case $c_\rho \neq 0$. Setting the right hand side of system (3) to 0, equations (3a) and (3c) give

$$ \bar{S} = \frac{\Lambda}{\beta \bar{N} + \mu} \quad \text{and} \quad \bar{N} = K \bar{I} \left( \frac{(r + c_{\rho} \rho - \nu) \bar{I} + (r - \nu) \alpha}{(r + c_{\rho} \rho) \bar{I} + r \alpha} \right). \qquad (31) $$
<figure><img src="image_13.png" /><figcaption>Figure 9: Simulations of system (3) (red curves) for $\mathcal{R} = 8.4$ and system (23), using $\epsilon = 1/10$ (black curves), $\epsilon = 1/100$ (green curves) and $\epsilon = 1/1000$ (blue curves). All parameter values are given in Table 1 and initial conditions in (28).</figcaption></figure>

<table><caption>Table 3: Descartes’ rule of signs for I-polynomial (30)</caption><thead><tr><th>m<sub>2</sub></th><th>m<sub>1</sub></th><th>m<sub>0</sub></th><th>R</th><th>Number of positive solutions</th></tr></thead><tbody><tr><td>+</td><td>+</td><td>+</td><td>R &lt; 1</td><td>No solution</td></tr><tr><td>+</td><td>+</td><td>-</td><td>R &gt; 1</td><td>1 solution</td></tr><tr><td>+</td><td>-</td><td>-</td><td>R &gt; 1</td><td>1 solution</td></tr></tbody></table>

Multiplying equation (3b) by $(\beta\bar{N} + \mu)(\alpha + \bar{I})[(r + c_{\rho}\rho)\bar{I} + r\alpha]^2$ and using equation (31), we obtain

$$f(\bar{I}) + g(\bar{I}) + h(\bar{I}) = 0, \qquad (32)$$

where

$$ \left\{ \begin{aligned} f(\bar{I}) &= \Lambda\beta K[(r + c_{\rho}\rho - \nu)\bar{I} + (r - \nu)\alpha](a + \bar{I})[(r + c_{\rho}\rho)\bar{I} + r\alpha], \\ g(\bar{I}) &= -\mu(\alpha + \bar{I})\left[\beta K \bar{I} \left(\frac{(r+c_{\rho}\rho-\nu)\bar{I}+(r-\nu)\alpha}{(r+c_{\rho}\rho)\bar{I}+r\alpha}\right) + \mu\right] [(r + c_{\rho}\rho)\bar{I} + r\alpha]^2, \\ h(\bar{I}) &= -\rho\beta K^2(\bar{I})^2 \left(\frac{(r+c_{\rho}\rho-\nu)\bar{I}+(r-\nu)\alpha}{(r+c_{\rho}\rho)\bar{I}+r\alpha}\right)^2 [(r + c_{\rho}\rho)\bar{I} + r\alpha]^2 \\ &\quad -\mu\rho K \bar{I} \left(\frac{(r+c_{\rho}\rho-\nu)\bar{I}+(r-\nu)\alpha}{(r+c_{\rho}\rho)\bar{I}+r\alpha}\right) [(r + c_{\rho}\rho)\bar{I} + r\alpha]^2. \end{aligned} \right. $$
After development, we obtain

$$
\left\{
\begin{aligned}
f(\bar{I}) &= \Lambda\beta K(r + c_{\rho}\rho)(r + c_{\rho}\rho - \nu)(\bar{I})^3 + \Lambda\beta K\alpha(r + c_{\rho}\rho)(r + c_{\rho}\rho - \nu)(\bar{I})^2 \\
&\quad + [\Lambda\beta Kr\alpha(r + c_{\rho}\rho - \nu) + \Lambda\beta K\alpha(r + c_{\rho}\rho)(r - \nu)](\bar{I})^2 \\
&\quad + [\Lambda\beta K\alpha^2 r(r + c_{\rho}\rho - \nu) + \Lambda\beta K\alpha^2(r + c_{\rho}\rho)(r - \nu) + \Lambda\beta Kr\alpha^2(r - \nu)]\bar{I} \\
&\quad + \Lambda\beta Kr\alpha^3(r - \nu), \\
g(I^*) &= -\mu\beta K(r + c_{\rho}\rho)(r + c_{\rho}\rho - \nu)(\bar{I})^4 - [\mu\alpha\beta K(r + c_{\rho}\rho)(r + c_{\rho}\rho - \nu)](\bar{I})^3 \\
&\quad - [\mu^2(r + c_{\rho}\rho)^2 + \mu\alpha\beta Kr(r + c_{\rho}\rho - \nu) + \mu\alpha\beta K(r + c_{\rho}\rho)(r - \nu)](\bar{I})^3 \\
&\quad - [\mu\alpha^2\beta K(r + c_{\rho}\rho)(r - \nu) + \mu^2\alpha(r + c_{\rho}\rho)^2 + \mu\alpha^2\beta Kr(r - \nu)](\bar{I})^2 \\
&\quad - [2\mu^2(r + c_{\rho}\rho)r\alpha + \mu\alpha^2\beta Kr(r + c_{\rho}\rho - \nu)](\bar{I})^2 - \mu\alpha^3\beta Kr(r - \nu)\bar{I} \\
&\quad - [2\mu^2(r + c_{\rho}\rho)r\alpha^2 + \mu^2r^2\alpha^2]\bar{I} - \mu^2r^2\alpha^3, \\
h(I^*) &= -\rho\beta K^2(r + c_{\rho}\rho - \nu)^2(\bar{I})^4 - 2\rho\beta K^2(r + c_{\rho}\rho - \nu)(r - \nu)\alpha(\bar{I})^3 \\
&\quad - \rho K(r + c_{\rho}\rho - \nu)\mu(r + c_{\rho}\rho)(\bar{I})^3 - \rho K(r + c_{\rho}\rho - \nu)\mu r\alpha(\bar{I})^2 \\
&\quad - [\rho\beta K^2(r - \nu)^2\alpha^2 + \rho K(r - \nu)\alpha\mu(r + c_{\rho}\rho)](\bar{I})^2 - \rho K(r - \nu)\alpha^2\mu r\bar{I}.
\end{aligned}
\right.
$$

Therefore, equation (32) is equivalent to

$$a_4(\bar{I})^4 + a_3(\bar{I})^3 + a_2(\bar{I})^2 + a_1\bar{I} + a_0 = 0, \quad (33)$$

where

$$
\left\{
\begin{aligned}
a_4 &= -\mu K \beta (r + c_{\rho} \rho - \nu) (r + c_{\rho} \rho) - \rho \beta K^2 (r + c_{\rho} \rho - \nu)^2, \\
a_3 &= \mu^2 r (r + c_{\rho} \rho) (\mathcal{R} - 1) + c_{\rho} \rho (r + c_{\rho} \rho) \frac{r \mu^2}{r - \nu} \left(\mathcal{R} - \frac{r - \nu}{r}\right) \\
&\quad - \mu \alpha \beta K (r + c_{\rho} \rho - \nu) (r + c_{\rho} \rho) - \mu K \beta r \alpha (r + c_{\rho} \rho - \nu) - \mu K \beta \alpha (r - \nu) (r + c_{\rho} \rho) \\
&\quad - 2 \rho \beta K^2 \alpha (r + c_{\rho} \rho - \nu) (r - \nu) - \rho K \mu (r + c_{\rho} \rho - \nu) (r + c_{\rho} \rho), \\
a_2 &= 2 \mu^2 r \alpha (r + c_{\rho} \rho) (\mathcal{R} - 1) + c_{\rho} \rho \alpha (r + c_{\rho} \rho) \frac{r \mu^2}{r - \nu} \left(\mathcal{R} - \frac{r - \nu}{r}\right) + \mu^2 r^2 \alpha (\mathcal{R} - 1) \\
&\quad + r \alpha c_{\rho} \rho \frac{r \mu^2}{r - \nu} \left(\mathcal{R} - \frac{r - \nu}{r}\right) - \mu \alpha^2 \beta K r (r + c_{\rho} \rho - \nu) - \mu \alpha^2 \beta K (r - \nu) (r + c_{\rho} \rho) \\
&\quad - \mu \beta K r \alpha^2 (r - \nu) - \rho K \mu r \alpha (r + c_{\rho} \rho - \nu) - \rho \beta K^2 (r - \nu)^2 \alpha^2 \\
&\quad - \rho K (r - \nu) \alpha \mu (r + c_{\rho} \rho), \\
a_1 &= 2 \mu^2 r^2 \alpha^2 (\mathcal{R} - 1) + r \alpha^2 c_{\rho} \rho \frac{r \mu^2}{r - \nu} \left(\mathcal{R} - \frac{r - \nu}{r}\right) + \mu^2 r \alpha^2 (r + c_{\rho} \rho) (\mathcal{R} - 1) \\
&\quad - \mu \alpha^3 \beta K r (r - \nu) - \rho K \mu r \alpha^2 (r - \nu), \\
a_0 &= \mu^2 r^2 \alpha^3 (\mathcal{R} - 1).
\end{aligned}
\right.
$$

**Remark 6.** If $\mathcal{R} < \frac{r-\nu}{r} < 1$, then $a_0, a_1, a_2, a_3, a_4 < 0$. Hence, equation (33) has no solution.

Now, using the Descartes sign rules for I-polynomial (33), the number of positive solutions is given in Table 4. Hence, system (3) admits 0 or 2 or 4 endemic equilibria when $\mathcal{R} < 1$, and 1 or 3 endemic equilibria when $\mathcal{R} > 1$.

## Appendix C. Singular control

Differentiating the switching function 11 with respect to time $t$ along the singular control and using the fact that $\lambda_2 - \lambda_1 = \frac{B_u}{\beta NS}$, yields

$$
\begin{aligned}
\dot{\phi}(t) = & -B_u \left(r + c_{\rho} \rho \frac{I}{\alpha + \bar{I}}\right) \left(1 - \frac{N}{K\bar{I}}\right) + B_u \nu - B_u \frac{\Lambda}{S} - \beta NS \left(\rho a \frac{N}{(\alpha + \bar{I})^2} \lambda_2 - B_S\right) \\
& + \beta NS \left(c_{\rho} \rho \alpha \frac{N}{(\alpha + \bar{I})^2} \left(1 - \frac{N}{K\bar{I}}\right) + \frac{1}{K} \left(\frac{N}{\bar{I}}\right)^2 \left(r + c_{\rho} \rho \frac{I}{\alpha + \bar{I}}\right)\right) \lambda_3
\end{aligned}
$$

Now, differentiating $\phi$ a second time since the expression of its first differential does not explicitly depend on the control, we obtain
<table><caption>Table 4: Descartes’ rule of signs for I-polynomial 33</caption><thead><tr><th>a<sub>4</sub></th><th>a<sub>3</sub></th><th>a<sub>2</sub></th><th>a<sub>1</sub></th><th>a<sub>0</sub></th><th>R</th><th>Number of positive solutions</th></tr></thead><tbody><tr><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>R &lt; 1</td><td>No solution</td></tr><tr><td>-</td><td>-</td><td>-</td><td>-</td><td>+</td><td>R &gt; 1</td><td>1 solution</td></tr><tr><td>-</td><td>-</td><td>-</td><td>+</td><td>-</td><td>R &lt; 1</td><td>0 or 2 solutions</td></tr><tr><td>-</td><td>-</td><td>-</td><td>+</td><td>+</td><td>R &gt; 1</td><td>1 solution</td></tr><tr><td>-</td><td>-</td><td>+</td><td>-</td><td>-</td><td>R &lt; 1</td><td>0 or 2 solutions</td></tr><tr><td>-</td><td>-</td><td>+</td><td>-</td><td>+</td><td>R &gt; 1</td><td>1 or 3 solutions</td></tr><tr><td>-</td><td>-</td><td>+</td><td>+</td><td>-</td><td>R &lt; 1</td><td>0 or 2 solutions</td></tr><tr><td>-</td><td>-</td><td>+</td><td>+</td><td>+</td><td>R &gt; 1</td><td>1 solution</td></tr><tr><td>-</td><td>+</td><td>-</td><td>-</td><td>-</td><td>R &lt; 1</td><td>0 or 2 solutions</td></tr><tr><td>-</td><td>+</td><td>-</td><td>-</td><td>+</td><td>R &gt; 1</td><td>1 or 3 solutions</td></tr><tr><td>-</td><td>+</td><td>-</td><td>+</td><td>-</td><td>R &lt; 1</td><td>0 or 2 ou 4 solutions</td></tr><tr><td>-</td><td>+</td><td>-</td><td>+</td><td>+</td><td>R &gt; 1</td><td>1 or 3 solutions</td></tr><tr><td>-</td><td>+</td><td>+</td><td>-</td><td>-</td><td>R &lt; 1</td><td>0 or 2 solutions</td></tr><tr><td>-</td><td>+</td><td>+</td><td>-</td><td>+</td><td>R &gt; 1</td><td>1 or 3 solutions</td></tr><tr><td>-</td><td>+</td><td>+</td><td>+</td><td>-</td><td>R &lt; 1</td><td>0 or 2 solutions</td></tr><tr><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>R &gt; 1</td><td>1 solution</td></tr></tbody></table>

$$
\begin{aligned}
\ddot{\phi}(t) = & -B_u \frac{c_\rho \rho \alpha}{(\alpha + I)^2} \left(1 - \frac{N}{KI}\right) \left(\beta(1-u)NS - \mu I - \rho N \frac{I}{\alpha + I}\right) \\
& + B_u \frac{1}{KI^2} \left(r + c_\rho \rho \frac{I}{\alpha + I}\right) \left(I \frac{dN}{dt} - \beta(1-u)N^2 S + \mu I N - \rho N^2 \frac{I}{\alpha + I}\right) \\
& + \frac{\Lambda}{S^2} (\Lambda - (1-u)\beta NS - \mu S) - \beta \left[S \frac{dN}{dt} + (\Lambda - (1-u)\beta NS - \mu S)N\right] \\
& \times \left(\frac{\rho \alpha N}{(\alpha + I)^2} \lambda_2 - \left(\frac{c_\rho \rho \alpha N}{(\alpha + I)^2} \left(1 - \frac{N}{KI}\right) + \frac{1}{K} \left(\frac{N}{I}\right)^2 \left(r + c_\rho \rho \frac{I}{\alpha + I}\right)\right) \lambda_3 - B_S\right) \\
& - \frac{\rho \alpha \beta NS}{(\alpha + I)^3} \left[(\alpha + I) \frac{dN}{dt} - 2 \left((1-u)\beta NS - \mu I - \rho N \frac{I}{\alpha + I}\right) N\right] \lambda_2 - \frac{\rho \alpha N}{(\alpha + I)^2} \beta NS \\
& + \left(1 - \frac{N}{KI}\right) \frac{c_\rho \rho \alpha}{(\alpha + I)^3} \beta NS \left((\alpha + I) \frac{dN}{dt} - 2 \left((1-u)\beta N^2 S - \mu I N - \frac{\rho N^2 I}{\alpha + I}\right)\right) \lambda_3 \\
& - c_\rho \rho \alpha \frac{N}{(\alpha + I)^2} \frac{1}{KI^2} \beta NS \left(I \frac{dN}{dt} - (1-u)\beta N^2 S + \mu I N + \rho N^2 \frac{I}{\alpha + I}\right) \lambda_3 \\
& + \left[c_\rho \rho \alpha \frac{N}{(\alpha + I)^2} \left(1 - \frac{N}{KI}\right) + \frac{1}{K} \left(\frac{N}{I}\right)^2 \left(r + c_\rho \rho \frac{I}{\alpha + I}\right)\right] \\
& \times \left[(1-u)\beta S(\lambda_1 - \lambda_2) + \rho \frac{I}{\alpha + I} \lambda_2 - \left(\left(r + c_\rho \rho \frac{I}{\alpha + I}\right) \left(1 - \frac{2N}{KI}\right) - \nu\right) \lambda_3\right] \beta NS \\
& + \frac{2}{KI^2} \left(I \frac{dN}{dt} - (1-u)\beta N^2 S + \mu I N + \rho N^2 \frac{I}{\alpha + I}\right) \left(\frac{N}{I}\right) \left(r + c_\rho \rho \frac{I}{\alpha + I}\right) \lambda_3 \beta NS \\
& + \frac{1}{K} \left(\frac{N}{I}\right)^2 \frac{c_\rho \rho \alpha}{(\alpha + I)^2} \left((1-u)\beta NS - \mu I - \rho N \frac{I}{\alpha + I}\right) \lambda_3 \beta NS
\end{aligned}
$$

As $\dot{\phi}(t) = 0$ for all $t \in [t_1, t_2]$, we have

$$
\begin{aligned}
& \beta NS \left( \frac{\rho \alpha N}{(\alpha + I)^2} \lambda_2 - \left( \frac{c_\rho \rho \alpha N}{(\alpha + I)^2} \left(1 - \frac{N}{KI}\right) + \frac{1}{K} \left(\frac{N}{I}\right)^2 \left(r + c_\rho \rho \frac{I}{\alpha + I}\right) \right) \lambda_3 - B_S \right) \\
& = -B_u \left( \frac{1}{N} \frac{dN}{dt} + \frac{\Lambda}{S} \right)
\end{aligned}
\qquad (34)
$$

Then, substituting equation (34) and the expression of $\lambda_2 - \lambda_1$ into the above expression of $\ddot{\phi}(t)$, we obtain

$$
\ddot{\phi}(t) = \mathcal{D}(x, \lambda) + (1 - u)\mathcal{F}(x, \lambda)
$$
where

$$
\begin{aligned}
\mathcal{D}(x, \lambda) &= \frac{B_u c_\rho \rho \alpha}{(\alpha + I)^2} \left[ \mu I + \rho N \frac{I}{\alpha + I} \right] \left( 1 - \frac{N}{KI} \right) \\
&\quad + \frac{B_u}{KI^2} \left[ I \frac{dN}{dt} + \mu N I + \rho N^2 \frac{I}{\alpha + I} \right] \left( r + c_\rho \rho \frac{I}{\alpha + I} \right) + \frac{\Lambda}{S^2} (\Lambda - \mu S) \\
&\quad - \left( \frac{\rho \alpha N}{(\alpha + I)^2} \lambda_2 - B_S \right) \left[ \beta S \frac{dN}{dt} + \beta (\Lambda - \mu S) N \right] \\
&\quad + \left( \frac{c_\rho \rho \alpha N}{(\alpha + I)^2} \left( 1 - \frac{N}{KI} \right) + \frac{1}{K} \left( \frac{N}{I} \right)^2 \left( r + c_\rho \rho \frac{I}{\alpha + I} \right) \right) \lambda_3 \\
&\quad \times \left[ \beta S \frac{dN}{dt} + \beta (\Lambda - \mu S) N \right] - \frac{\rho \alpha \beta N S}{(\alpha + I)^3} \left[ (\alpha + I) \frac{dN}{dt} + 2 \left( \mu I + \rho N \frac{I}{\alpha + I} \right) N \right] \lambda_2 \\
&\quad - \frac{\rho \alpha}{(\alpha + I)^2} \left[ -B_u \left( \frac{dN}{dt} + \frac{\Lambda N}{S} \right) + \beta N^2 S (\mu \lambda_2 + B_1) \right] \\
&\quad + \left( 1 - \frac{N}{KI} \right) \frac{c_\rho \rho \alpha}{(\alpha + I)^3} \beta N S \left( (\alpha + I) \frac{dN}{dt} + 2 \left( \mu I N + \rho N^2 \frac{I}{\alpha + I} \right) \right) \lambda_3 \\
&\quad - \frac{c_\rho \rho \alpha}{(\alpha + I)^2} \frac{1}{KI^2} \beta N^2 S \left( I \frac{dN}{dt} + \mu I N + \rho N^2 \frac{I}{\alpha + I} \right) \lambda_3 \\
&\quad + \frac{c_\rho \rho \alpha \beta N^2 S}{(\alpha + I)^2} \left( 1 - \frac{N}{KI} \right) \left[ \rho \frac{I}{\alpha + I} \lambda_2 - \left( \left( r + c_\rho \rho \frac{I}{\alpha + I} \right) \left( 1 - \frac{2N}{KI} \right) - \nu \right) \lambda_3 \right] \\
&\quad + \frac{2 \beta N^2 S}{KI^3} \left( I \frac{dN}{dt} + \mu I N + \rho N^2 \frac{I}{\alpha + I} \right) \left( r + c_\rho \rho \frac{I}{\alpha + I} \right) \lambda_3 \\
&\quad - \frac{1}{KI^2} \frac{c_\rho \rho \alpha}{(\alpha + I)^2} \beta N^3 S \left( \mu I + \rho N \frac{I}{\alpha + I} \right) \lambda_3 \\
&\quad + \frac{\beta N^3 S}{KI^2} \left( r + c_\rho \rho \frac{I}{\alpha + I} \right) \left[ \rho \frac{I}{\alpha + I} \lambda_2 - \left( \left( r + c_\rho \rho \frac{I}{\alpha + I} \right) \left( 1 - \frac{2N}{KI} \right) - \nu \right) \lambda_3 \right],
\end{aligned}
$$
$$
\mathcal{F}(x, \lambda) = -2 \frac{c_\rho \rho \alpha}{(\alpha + I)^2} \beta N S \left( 1 - \frac{N}{KI} \right) \left[ B_u + \frac{\beta N^2 S \lambda_3}{\alpha + I} \right] - 2 B_u \frac{\Lambda \beta N}{S} \\
- B_u \beta \frac{dN}{dt} - \frac{2}{KI^2} \beta N^2 S \left( r + c_\rho \rho \frac{I}{\alpha + I} \right) \left[ B_u + \frac{\beta N^2 S \lambda_3}{I} \right] \\
+ 2 \frac{\rho \alpha}{(\alpha + I)^2} \beta^2 N^3 S^2 \left[ \frac{\lambda_2}{\alpha + I} + \frac{c_\rho N \lambda_3}{KI^2} \right].
$$

Hence, as $\ddot{\phi}(t) = 0$, we obtain the expression of the singular control as follows

$$ u_{\text{sing}}^{*}(t) = 1 + \frac{\mathcal{D}(x, \lambda)}{\mathcal{F}(x, \lambda)}. $$

## Appendix D. static-OCP

Trivially, first and third equations of system (18) give

$$ S = \frac{\Lambda}{\beta(1-u)N + \mu} \quad \text{and} \quad N = KI \left( \frac{(r + c_{\rho}\rho - \nu)I + (r - \nu)\alpha}{(r + c_{\rho}\rho)I + r\alpha} \right). \qquad (35) $$

Now, combining the second equation of system (18) and equation (35), and after development, we obtain

$$ a_4 I^4 + a_3 I^3 + a_2 I^2 + a_1 I + a_0 = 0, \qquad (36) $$
where

$$
\left\{
\begin{aligned}
a_4 &= -\mu K \beta(1-u)(r+c_\rho\rho-\nu)(r+c_\rho\rho) - \rho\beta(1-u)K^2(r+c_\rho\rho-\nu)^2, \\
a_3 &= \mu^2 r(r+c_\rho\rho)(\mathcal{R}-1) + c_\rho\rho(r+c_\rho\rho)\frac{r\mu^2}{r-\nu}\left(\mathcal{R}-\frac{r-\nu}{r}\right) \\
    &\quad - \mu\alpha\beta(1-u)K(r+c_\rho\rho-\nu)(r+c_\rho\rho) - \mu K \beta(1-u)r\alpha(r+c_\rho\rho-\nu_I) \\
    &\quad - \rho K \mu(r+c_\rho\rho-\nu)(r+c_\rho\rho) - \mu K \beta(1-u)\alpha(r-\nu)(r+c_\rho\rho) \\
    &\quad - 2\rho\beta(1-u)K^2\alpha(r+c_\rho\rho-\nu)(r-\nu_I), \\
a_2 &= 2\mu^2 r\alpha(r+c_\rho\rho)(\mathcal{R}-1) + c_\rho\rho\alpha(r+c_\rho\rho)\frac{r\mu^2}{r-\nu}\left(\mathcal{R}-\frac{r-\nu}{r}\right) + \mu^2 r^2 \alpha(\mathcal{R}-1) \\
    &\quad + r a c_\rho \rho \frac{r\mu^2}{r-\nu}\left(\mathcal{R}-\frac{r-\nu}{r}\right) - \mu\alpha^2\beta(1-u)Kr(r+c_\rho\rho-\nu_I) \\
    &\quad - \mu\alpha^2\beta(1-u)K(r-\nu)(r+c_\rho\rho) - \mu\beta(1-u)Kr\alpha^2(r-\nu) \\
    &\quad - \rho K(r-\nu)\alpha\mu(r+c_\rho\rho) - \rho K \mu r \alpha(r+c_\rho\rho-\nu) - \rho\beta(1-u)K^2(r-\nu)^2\alpha^2, \\
a_1 &= 2\mu^2 r^2 \alpha^2(\mathcal{R}-1) + r\alpha^2 c_\rho\rho \frac{r\mu^2}{r-\nu}\left(\mathcal{R}-\frac{r-\nu}{r}\right) + \mu^2 r\alpha^2(r+c_\rho\rho)(\mathcal{R}-1) \\
    &\quad - \mu\alpha^3\beta(1-u)Kr(r-\nu) - \rho K \mu r \alpha^2(r-\nu), \\
a_0 &= \mu^2 r^2 \alpha^3(\mathcal{R}-1) \quad \text{with} \quad \mathcal{R} = \frac{\Lambda\beta(1-u)K}{\mu^2}\left(1-\frac{\nu}{r}\right).
\end{aligned}
\right.
$$

Hence, we obtain the following formulation of the static-OCP

$$
\begin{array}{ll}
\displaystyle \min_{u \in [0, u_{\max}]} & \bar{J}(u) = B_u u - B_S S, \\
\text{under the constraints:} \\
\displaystyle S = & \frac{\Lambda}{\beta(1-u)N + \mu}, \\
\displaystyle N = & KI \left( \frac{(r + c_{\rho}\rho - \nu)I + (r - \nu)\alpha}{(r + c_{\rho}\rho)I + r\alpha} \right), \\
\displaystyle a_4 I^4 + a_3 I^3 + a_2 I^2 + a_1 I + a_0 = & 0.
\end{array}
$$

References

[1] Agouanet, F.P., Tewa, J.J., Etoua, M.R.: Mathematical modelling and optimal control of banana black sigatoka disease. In: CARI 2020-Colloque Africain sur la Recherche en Informatique et Mathématiques Appliquées (2020). URL https://hal.science/hal-02926131

[2] Akinro, A., Olufayo, A., Oguntunde, P.: Crop water productivity of plantain (Musa sp) in a humid tropical environment. Journal of Engineering Science & Technology Review 5(1), 19–25 (2012)

[3] Asha, M.K., Mwengi, A.I., M'Vubu, A.R.N.: Analyse de la performance économique de la production de banane plantain et son impact socio-économique sur les ménages agricoles dans le village Mbubu au Kongo Central en RD Congo. ESI Preprints 32, 243–243 (2024). doi:10.19044/esj.2024.v20n25p81

[4] Baala, Y., Rachik, M.: Discrete mathematical modeling and optimal control for Bayoud disease of date palm. Communications in Mathematical Biology and Neuroscience 2021, 68 (2021). doi:10.28919/cmbn/5859

[5] Bonnans, F., Martinon, P., Giorgi, D., Grélard, V., Maindrault, S., Tissot, O., Liu, J.: Bocop 2.0.5 - user guide (2017)

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

[7] Chabrier, C.: Survie et dissémination du nématode Radopholus similis (Cobb) Thorne dans les sols bruns-rouilles à halloysites (nitisols): effets de l'état hydrique et des flux hydriques. Ph.D. thesis, Université des Antilles-Guyane (2008). URL https://theses.hal.science/tel-00697213

[8] Chen-Charpentier, B., Jackson, M.: Optimal control of plant virus propagation. Mathematical Methods in the Applied Sciences 43(14), 8147–8157 (2020). doi:10.1002/mma.6244
[9] Djema, W., Bayen, T., Bernard, O.: Optimal darwinian selection of microorganisms with internal storage. Processes **10**(3), 461 (2022). doi:10.3390/pr10030461

[10] Fleming, W.H., Rishel, R.W.: Deterministic and stochastic optimal control, *Applications of Mathematics*, vol. 1. Springer (2012)

[11] Fogain, R.: Effect of Radopholus similis on plant growth and yield of plantains (Musa, aab). Nematology **2**(2), 129-133 (2000). doi:10.1163/156854100509015

[12] Fotso, S., Kolaye, G., Ntahomvukiye, J., Bowong, S., Taffouo, V.: Modelling the influence of climatic factors on the population dynamics of Radopholus similis: Banana-plantain pest. Acta Biotheoretica **70**(3), 21 (2022). doi:10.1007/s10441-022-09444-4

[13] Fotso Fotso, Y., Touzeau, S., Tsanou, B., Bowong, S., Grognard, F.: Modelling and optimal strategy to control coffee berry borer. Mathematical Methods in the Applied Sciences **44**(18), 14,569-14,592 (2021). doi:10.1002/mma.7726

[14] Fulgence, T.D.E., Zadjéhi, K.E.B., Emmanuel, I.A., Stanislas, K.T., Philippe, G.: Évaluation de la variabilité génétique d'une collection de bananiers plantain (Musa sp.) cultivés en Côte d'Ivoire à partir d'une caractérisation agro-morphologique. European Scientific Journal **18**(40), 366 (2022). doi:10.19044/esj.2022.v18n40p366

[15] Gowen, S.R., Quénéhervé, P., Fogain, R.: Nematode parasites of bananas and plantains. In: Plant parasitic nematodes in subtropical and tropical agriculture, pp. 611-643. CABI Publishing (2005)

[16] Imizcoz, J.I., Djema, W., Mairet, F., Gouzé, J.L.: Optimal cell growth under periodic environment. In: 26th International Symposium on Mathematical Theory of Networks and Systems (2024). URL https://inria.hal.science/hal-04724783

[17] Khalil, H.K., Grizzle, J.W.: Nonlinear systems, *Information and System Sciences Series*, vol. 3. Prentice Hall (2002)

[18] Ledzewicz, U., Maurer, H., Schättler, H.: Bang-bang optimal controls for a mathematical model of chemo-and immunotherapy in cancer. Discrete and Continuous Dynamical Systems-B **29**(3), 1481-1500 (2024). doi:10.3934/dcdsb.2023141

[19] Mouliom Péfoura, A., Fogain, R., Cohan, P., Messiaen, S., Tomekpé, K.: Pests and diseases of bananas and plantains in the Western and Central Africa: Impact and integrated control strategies. In: International Conference on IPM in Sub-Saharan Africa, Kampala, Uganda, pp. 26–37 (2002)

[20] Pontryagin, L., Boltyanskiy, V., Gamkrelidze, R.V., Mishchenko, Y.: Mathematical theory of optimal processes. Wiley (1962)

[21] Schättler, H., Ledzewicz, U.: Geometric optimal control: theory, methods and examples, *Interdisciplinary Applied Mathematics*, vol. 38. Springer (2012)

[22] Storn, R., Price, K.: Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization **11**, 341–359 (1997). doi:10.1023/A:1008202821328

[23] Tankam-Chedjou, I., Touzeau, S., Mailleret, L., Tewa, J.J., Grognard, F.: Modelling and control of a banana soilborne pest in a multi-seasonal framework. Mathematical Biosciences **322**, 108,324 (2020). doi:10.1016/j.mbs.2020.108324

[24] Tapi, M.D., Yakam, A.N., Wafo, R.T., Bowong, S.: Mathematical modelling and optimal control of production losses caused by Miridae. Mathematical Modelling of Natural Phenomena **18**, 28 (2023). doi:10.1051/mmnp/2023030

[25] Tixier, P., Risède, J.M., Dorel, M., Malézieux, E.: Modelling population dynamics of banana plant-parasitic nematodes: A contribution to the design of sustainable cropping systems. Ecological Modelling **198**(3-4), 321-331 (2006). doi:10.1016/j.ecolmodel.2006.05.003

[26] Torres-Asuaje, P.E., Cotes-Prado, A.M., Echeverría-Beirute, F., Blanco-Rojas, F.A., Sandoval-Fernández, J.A., Segura-Mena, R.A., Palomares-Rius, J.E.: Ensilaged biostimulants promoting root health and control of Radopholus similis in banana (Musa aaa) cv. Grande Naine. European Journal of Plant Pathology **165**(3), 465-474 (2023). doi:10.1007/s10658-022-02617-4

[27] Trélat, E., Zuazua, E.: The turnpike property in finite-dimensional nonlinear optimal control. Journal of Differential Equations **258**(1), 81–114 (2015). doi:10.1016/j.jde.2014.09.005