# Mathematical modelling and optimal control of the seasonal coffee leaf rust propagation

Clotilde Djuikem, Agustín Gabriel Yabo, Frédéric Grognard, Suzanne Touzeau

► To cite this version:

Clotilde Djuikem, Agustín Gabriel Yabo, Frédéric Grognard, Suzanne Touzeau. Mathematical modelling and optimal control of the seasonal coffee leaf rust propagation. ADHS 2021 - 7th IFAC Conference on Analysis and Design of Hybrid Systems, IFAC, Jul 2021, Brussels, Belgium. pp.193-198, 10.1016/j.ifacol.2021.08.497. hal-03274865

HAL Id: hal-03274865

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

Submitted on 30 Jun 2021

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.
Mathematical modelling and optimal
control of the seasonal coffee leaf rust
propagation*

Clotilde Djuikem* Agustín Gabriel Yabo*
Frédéric Grognard * Suzanne Touzeau **,**

* Université Côte d'Azur, Inria, INRAE, CNRS, Sorbonne Université,
BIOCORE, Sophia Antipolis, France (e-mail: clotilde.djuikem@inria.fr,
agustin.yabo@inria.fr, frederic.grognard@inria.fr)

** Université Côte d'Azur, INRAE, CNRS, ISA, Sophia Antipolis,
France (e-mail: suzanne.touzeau@inria.fr)

**Abstract:** This paper deals with the problem of modelling and optimal control of the coffee leaf rust (CLR) in a coffee plantation. We first formulate and analyze a mathematical model for the propagation of CLR with impulsive effect, which represents the non-production period of coffee. Then, the CLR control is formulated and solved as an optimal control, showing how a biological control strategy can be implemented in a coffee plantation to reduce the damage caused by CLR, therefore increasing the coffee production. Numerical simulations are performed to illustrate and validate the theoretical results.

**Keywords:** Epidemiological modeling; Impulsive system; Floquet theory; Stability; Optimal control, *Hemileia vastatrix*;.

# 1. INTRODUCTION

Coffee is one of the most widely consumed beverages in the world, with their consumption in 2019/20 being estimated at 168 million 60-kilogram bags (International Coffee Organization, 2020). Its cultivation is widespread in almost all the countries of the tropic and is an important factor of social stability as it supports no less than twenty-five million small producers and their families worldwide. However, due to the climate of these areas, multiple coffee diseases have persisted for several years. Among these pests, the Coffee Leaf Rust (CLR) is considered the most important leaf disease known to this day for being one of the most devastating pathogens. It is caused by the fungus *Hemileia vastatrix*, first discovered in Sri Lanka in 1869, and nowadays present in the majority of the coffee-growing regions of the world. The pathogenic principle behind CLR is that it causes premature defoliation, which reduces the photosynthetic capacity of coffee leaves, inducing severe yield losses to the sector (up to 70%).

The dynamics of CLR is complex due to the multiple interactions between the plant host and the pathogen in the environment. A better understanding of such CLR dynamics can significantly improve the accuracy of the models, which can be then analyzed through theoretical studies and numerical simulations. Many efforts have been and are still being devoted to the modeling of the propagation of CLR (Vandermeer et al., 2018). However, none of them take into account the impulsive effect in the dynamics. Considering discrete events in the model—such as harvesting—could allow for a more comprehensive

representation of the plantation process, providing in turn a more precise framework for pest treatment. In this regard, several works have used ordinary impulsive differential equations for representing population dynamics in ecology and epidemiology in the recent years (Mailleret and Grognard, 2009; Tamen et al., 2017; Tankam-Chedjou et al., 2020; d'Onofrio, 2002).

The capacity to accurately model the time evolution of coffee plantations allows not only to effectively predict trends in the process, but also to act on them. In this context, optimal control theory has proven to be a powerful tool for investigating potential control strategies in pest treatment (Abbasi et al., 2020; Wong and Tong, 2012; Belbas and Schmidt, 2009). Such approaches are commonly based on the well-known PMP (Pontryagin's maximum principle) (Pontryagin, 2018): for a given cost function to maximize, this theory can provide necessary (and often sufficient) conditions for optimality of control strategies in systems of ordinary differential equations, partial differential equations and hybrid systems with given constraints (Dmitruk and Kaganovich, 2008).

The present contribution offers an impulsive perspective to the modeling of the propagation of CLR in a coffee plantation. We first formulate and analyze a non-controlled impulsive model for the propagation of CLR, for which we compute the periodic disease-free solution (PDFS) and a threshold parameter $R$. Then, we conduct an in-depth analysis of the global asymptotic stability of the PDFS when $R < 1$, and we perform numerical simulations to illustrate the theoretical results. Furthermore, we devise an extended model that considers a biological control for pest treatment, intended to maximize the coffee production

<footnote>* The financial support of this work are Inria and EPITAG, an Inria Associate team part of the LIRIMA (https://team.inria.fr/epitag/).</footnote>
while minimizing the damages undergone by the CLR. The latter objective is written as an Hybrid OCP (Optimal Control Problem), for which we propose an analogous problem in continuous form, which allows to obtain numerical solutions with Bocop (Team Commands, 2017). Later on, we provide a numerical solution of an optimal trajectory of the system for a relevant set of parameters and initial conditions. As a result, we show that a cost-effective inoculative biological control can successfully contain the damages on a coffee plantation.

## 2. THE MODEL AND ITS BASIC PROPERTIES

### 2.1 Mathematical Modeling of coffee leaf rust

In this section, we formulate a mathematical model for the propagation of the CLR in the coffee plantation. To do so, we consider the dynamics of fungus during the production season—which corresponds to the rainy season for some countries—through ordinary differential equations; and during the non-production period—which is the dry season for some countries—represented by the impulsive effect. More precisely, we consider that in a coffee plantation we can find: susceptible branches S, infected branches I and urediniospores U. During the production period with length $T > 0$, the recruitment of healthy branches occurs at rate $\Lambda$. Urediniospores are deposited on leaves of all branches at rate $\nu$, and a fraction $S/N$ lands on susceptible branches. When the conditions are favorable, urediniospores germinate with efficiency rate $\omega$ and the susceptible branches becomes infected branches. All branches undergo natural mortality with baseline rate $\mu$ and the infected branches have an additional mortality rate $d$ due to the disease. Urediniospores are produced by infected branches at rate $\gamma$ and lose their ability to infect coffee branches at constant rate $\mu_U$. Susceptible and infected branches produce the berries with constant rate $\delta_S$ and $\delta_I$, respectively. At the end of the production period, harvest occurs instantaneously and we consider that the dry season or non-production period corresponds to the time when the number of branches at the end of the production period switches, this switch due to the fact that harvesting reduces the number of susceptible branches, infected branches and urediniospores with rates $\varphi_S$, $\varphi_I$ and $\varphi_U$ respectively. Using the fact that urediniospores lose their ability to infect fastly, we assume that $\varphi_U$ is very close to 0. The model is given by the following hybrid system:

$$ 
\left\{
\begin{aligned}
\dot{S} &= \Lambda - \frac{\omega\nu U}{N}S - \mu S, \quad t \neq nT; \\
\dot{I} &= \frac{\omega\nu U}{N}S - (\mu + d)I, \quad t \neq nT; \\
\dot{U} &= \gamma I - (\nu + \mu_U)U, \quad t \neq nT; \\
\dot{B} &= \delta_S S + \delta_I I - \mu_B B, \quad t \neq nT; \\
S(nT^+) &= \varphi_S S(nT); \\
I(nT^+) &= \varphi_I I(nT); \\
U(nT^+) &= \varphi_U U(nT); \\
B(nT^+) &= 0.
\end{aligned}
\right.
\quad (1) $$

where $S(t)$, $I(t)$, $U(t)$, $N(t) = S(t) + I(t)$ and $B(t)$ represent the number of susceptible branches, infected branches, urediniospores, total branches, and berries, respectively at

time $t$. With $0 < \varphi_I \le \varphi_S < 1$ and $\varphi_U > 0$ but close to 0 due to the fact that, at the end of season (which is also the harvest period), the number of spores is reduced considerably.

Most of the parameter values are based in literature values calibrated from ground-truth data. In particular, the reduction rates $\varphi_S$, $\varphi_I$ and $\varphi_U$ (of S, I and U respectively) are set to 0.7/day, 0.4/day and 0.001/day, respectively. For more information on the model calibration and the validity of the parameters see Djuikem et al. (2021).

### 2.2 Basic properties of the model

Let $\mathbb{R}_+^4 = \{X \in \mathbb{R}^4, X \ge 0\}$. Denote by $f = (f_1, f_2, f_3, f_4)$ the map given by the right-hand side of (1), and consider the initial conditions

$$ (S(0), I(0), U(0), B(0)) \in \mathbb{R}_+^4, S(0) > 0. \quad (2) $$

The solutions of (1) are non-negatives. Indeed, suppose that one of the variables, that we will denote $x$, is equal to 0 at some instant, with all the others being non-negative. A quick analysis shows that $\dot{x} \ge 0$ cannot become negative. This implies that $x(t) \ge 0$. Then, the solutions are non-negatives. The solutions of (1) are piecewise continuous functions $X: \mathbb{R}_+ \to \mathbb{R}_+^4$, $X$ is continuous on $(nT, (n+1)T]$, $n \in \mathbb{N}$ and $X(nT^+) = \lim_{t \to nT^+} X(t)$. Clearly, the smooth properties of $f$ guarantee the existence and uniqueness of positive solutions of (1). By adding the first and second equations in (1), we see that the dynamics of total branches satisfies

$$ 
\left\{
\begin{aligned}
\dot{N} &= \Lambda - \mu N - dI, \quad t \neq nT \\
N(nT^+) &\le \varphi_S N(nT).
\end{aligned}
\right.
\quad (3) $$

From (3), and using the positivity of the variables of the model, we have $\dot{N} \le \Lambda - \mu N$ and $N(nT^+) \le \varphi_S N(nT)$. It follows that $0 \le \lim_{t \to \infty} N(t) \le \frac{\Lambda}{\mu}$. Then, solutions of system (1) are bounded.

## 3. MATHEMATICAL ANALYSIS

Herein, we present the mathematical analysis of system (1). Since the state variable of berries $B$ is not present in the other equations of system (1), we do not consider it in the mathematical analysis.

### 3.1 Periodic disease free solution and its stability

The periodic disease free solution (PDFS) occurs when $I = 0$ and $U = 0$. Replacing these values in the impulsive differential equation (IDE) (1), we obtain

$$ \dot{S} = \Lambda - \mu S, \quad t \neq nT; \quad S(nT^+) = \varphi_S S(nT). \quad (4) $$

The resolution of the first equation of (4) gives:

For $t \in (nT, (n+1)T]$

$$ S(t) = \frac{\Lambda}{\mu} + \left( S(nT^+) - \frac{\Lambda}{\mu} \right) e^{-\mu(t-nT)}. \quad (5) $$

This implies that

$$ S((n+1)T) = \frac{\Lambda}{\mu} + \left( S(nT^+) - \frac{\Lambda}{\mu} \right) e^{-\mu T}. $$
Using the impulsive condition $S((n+1)T^+) = \varphi_S S((n+1)T)$ yields

$$S((n+1)T^+) = \varphi_S \left( \frac{\Lambda}{\mu} + (S(nT^+) - \frac{\Lambda}{\mu})e^{-\mu T} \right). \quad (6)$$

The fixed point of (6) is given by

$$S(nT^+) = \frac{\Lambda \varphi_S (1 - e^{-\mu T})}{\mu (1 - \varphi_S e^{-\mu T})} > 0.$$

Substituting the value of $S(nT)$ into (5), for $t \in (nT, (n+1)T]$ one has

$$eS^P(t) = \frac{\Lambda}{\mu} \left[ 1 - \frac{(1 - \varphi_S)e^{\mu T}}{e^{\mu T} - \varphi_S} e^{-\mu(t-nT)} \right]. \quad (7)$$

Finally, the PDFS is $X^P(t)=(S^P(t), 0, 0)$, where $S^P(t)$ is defined as in (7).

## 3.2 Stability of the PDFS

We first study the local stability of $X^P(t)$ using small amplitude perturbation methods. Let us denote $\tilde{X}(t) = X(t) - X^P(t)$ where $X(t) = (S(t), I(t), U(t))^T$ and $\tilde{X}(t)$ is understood to be small amplitude perturbations. Substituting the expression of $\tilde{X}(t)$ in (1) gives

$$\left\{
\begin{aligned}
\dot{\tilde{S}} &= -\frac{\omega\nu\tilde{U}}{\tilde{N} + S^P(t)}(\tilde{S} + S^P(t)) - \mu\tilde{S}, \quad t \neq nT; \\
\dot{\tilde{I}} &= \frac{\omega\nu\tilde{U}}{\tilde{N} + S^P(t)}(\tilde{S} + S^P(t)) - (\mu + d)\tilde{I}, \quad t \neq nT; \\
\dot{\tilde{U}} &= \gamma\tilde{I} - (\nu + \mu_U)\tilde{U}, \quad t \neq nT; \\
\tilde{S}(nT^+) &= \varphi_S\tilde{S}(nT); \\
\tilde{I}(nT^+) &= \varphi_I\tilde{I}(nT); \\
\tilde{U}(nT^+) &= \varphi_U\tilde{U}(nT).
\end{aligned}
\right.
\quad (8)$$

where $\tilde{N}(t) = \tilde{S}(t) + \tilde{I}(t)$ and $N^P(t) = S^P(t)$. The linearization of (8) gives

$$\left\{
\begin{aligned}
\dot{\tilde{X}}(t) &= A\tilde{X}(t), \quad t \neq nT; \\
\tilde{X}(nT^+) &= \operatorname{diag}(\varphi_S, \varphi_I, \varphi_U)\tilde{X}(nT);
\end{aligned}
\right.
\quad (9)$$

where

$$A = \begin{pmatrix} -\mu & 0 & -\omega\nu \\ 0 & -(\mu+d) & \omega\nu \\ 0 & \gamma & -(\nu+\mu_U) \end{pmatrix}.$$

Solving (9) yields $\tilde{X}(t) = \Phi(t)\tilde{X}(0)$, where $\Phi(0) = I$ with $I$ being the identity matrix of dimension 3 and $\Phi(t)$ the fundamental matrix that satisfies $\frac{d\Phi(t)}{dt} = A\Phi(t)$. Then, one has

$$\Phi(t) = \begin{pmatrix} e^{-\mu t} & \psi_{12}(t) & \psi_{13}(t) \\ 0 & \psi_{22}(t) & \psi_{23}(t) \\ 0 & \psi_{32}(t) & \psi_{33}(t) \end{pmatrix},$$

where

$$\left\{
\begin{aligned}
\psi_{22}(t) &= \frac{1}{2\beta}[(\beta + k_1 - k_2)e^{\lambda_1 t} + (\beta - k_1 + k_2)e^{\lambda_2 t}], \\
\psi_{23}(t) &= \frac{\nu\omega}{\beta}[-e^{\lambda_1 t} + e^{\lambda_2 t}], \\
\psi_{33}(t) &= \frac{1}{2\beta}[(\beta - k_1 + k_2)e^{\lambda_1 t} + (\beta + k_1 - k_2)e^{\lambda_2 t}].
\end{aligned}
\right.
$$

with

$$\left\{
\begin{aligned}
k_1 &= \mu + d, & k_2 &= \nu + \mu_U, \\
\alpha &= k_1 + k_2, & \beta &= \sqrt{(k_1 - k_2)^2 + 4\gamma\omega\nu}, \\
\lambda_1 &= -\frac{\alpha}{2} - \frac{\beta}{2} & \text{and} \quad \lambda_2 &= -\frac{\alpha}{2} + \frac{\beta}{2}.
\end{aligned}
\right.$$

Using the impulsive conditions, one obtains the impulsive matrix $\tilde{X}(nT^+) = \operatorname{diag}(\varphi_S, \varphi_S, \varphi_U)\tilde{X}(nT)$. Finally, the solution of (9) becomes

$$\tilde{X}(t) = \operatorname{diag}(\varphi_S, \varphi_I, \varphi_U)\Phi(t)\tilde{X}(0).$$

From the above equation, the monodromy matrix is $M = \operatorname{diag}(\varphi_S, \varphi_I, \varphi_U)\Phi(T)$, i.e

$$M = \begin{pmatrix} \varphi_S e^{-\mu T} & \varphi_S \xi_1(T) & \varphi_S \xi_2(T) \\ 0 & \varphi_I \psi_{22}(T) & \varphi_I \psi_{23}(T) \\ 0 & \varphi_U \psi_{32}(T) & \varphi_U \psi_{33}(T) \end{pmatrix}.$$

The Floquet multipliers of the matrix $M$ is given by $\chi_1 = \varphi_S e^{-\mu T} < 1$ and the multipliers of the submatrix

$$M_1 = \begin{pmatrix} \varphi_I \psi_{22}(T) & \varphi_I \psi_{23}(T) \\ \varphi_U \psi_{32}(T) & \varphi_U \psi_{33}(T) \end{pmatrix}.$$

Then, the PDFS $X^P(t)$ is locally asymptotically stable if the eigenvalues of matrix $M_1$ stay in the unit cycle, which need to satisfy the following Jury conditions (Wang, 2012):

$$\left\{
\begin{array}{l}
-\operatorname{tr}(M_1) - \operatorname{det}(M_1) < 1, \\
\operatorname{det}(M_1) < 1, \\
\operatorname{tr}(M_1) - \operatorname{det}(M_1) < 1,
\end{array}
\right.
\quad (10)$$

where

$$\left\{
\begin{array}{l}
\operatorname{tr}(M_1) = \varphi_I \psi_{22}(T) + \varphi_U \psi_{33}(T), \\
\operatorname{det}(M_1) = \varphi_I \varphi_U (\psi_{22}(T)\psi_{33}(T) - \psi_{23}(T)\psi_{32}(T)),
\end{array}
\right.
$$

The conditions (10) hold if $\varphi_U \to 0$ and $\varphi_I \psi_{22}(T) < 1$. Finally, one has the following Lemma.

**Lemma 1.** The PDFS $X^P(t)$ is locally asymptotically stable provided that $\mathcal{R} < 1$, where

$$\mathcal{R} = \frac{\varphi_I}{2\beta} [(\beta + k_1 - k_2)e^{\lambda_1 T} + (\beta - k_1 + k_2)e^{\lambda_2 T}].$$

We have the following result about the global stability of the PDFS of (1).

**Theorem 2.** The PDFS $X^P(t)$ of (1) is globally asymptotically stable in region $G$ provided that $\mathcal{R} < 1$.

**Proof.** We will prove the global attractivity of the PDFS $X^P(t)$. Consider the subsystem following:

$$\left\{
\begin{aligned}
\dot{I} &= \frac{\omega\nu U}{N} S - (\mu + d)I, \quad t \neq nT; \\
\dot{U} &= \gamma I - (\nu + \mu_U)U, \quad t \neq nT; \\
I(nT^+) &= \varphi_I I(nT), \quad U(nT^+) = \varphi_U U(nT).
\end{aligned}
\right.
\quad (11)$$

The trivial solution of (11) is $(0, 0)$ and the solution $(S(t), I(t), U(t)) \in G$. The system (11) is cooperative because $(\frac{\partial I}{\partial U} = \frac{\omega\nu S}{N} > 0, \frac{\partial U}{\partial I} = \gamma > 0)$ and using the fact that $\frac{S(t)}{N(t)} < 1$, one has the following cooperative system:

$$\left\{
\begin{aligned}
\dot{I}_1 &= \omega\nu U_1 - (\mu + d)I_1, \quad t \neq nT; \\
\dot{U}_1 &= \gamma I_1 - (\nu + \mu_U)U_1, \quad t \neq nT; \\
I_1(nT^+) &= \varphi_I I_1(nT), \quad U_1(nT^+) = \varphi_U U_1(nT).
\end{aligned}
\right.
\quad (12)$$

Applying the Kamke's theorem (Coppel, 1965), one has that $I(t) \le I_1(t), U(t) \le U_1(t)$. Asymptotically, system
<figure><img src="image_1.png" /><figcaption>Fig. 1. Trajectories of (1) when $\gamma = 1.6$ (so that $\mathcal{R}_0 = 0.36$).</figcaption></figure>

(12) behaves like the system (11) (Hale, 1969). Then, using Lemma 1 when $\mathcal{R}_0 < 1$, one has

$$I_1(t), U_1(t) \to 0^+ \Rightarrow I(t), U(t) \to 0^+.$$

Since $U \to 0^+$, for any arbitrary positive $\varepsilon_U$ there exists a $t_0 > 0$ such that for $t \ge t_0$, $U(t) \le \varepsilon_U$. Using the positivity of the solutions and the fact that $S(t)/N(t) \le 1$ into the first equation of (1), one has:

$$\Lambda - \omega\nu\varepsilon_U - \mu S(t) \le \dot{S}(t) \le \Lambda - \mu S(t)$$

Applying the comparison theorem on the above differential inequalities, one has $S_1(t) \le S(t) \le S_2(t)$, where for $t \in [nT; (n+1)T]$,

$$\left\{ \begin{array}{l} S_1(t) = \frac{\Lambda - \omega\nu\varepsilon_U}{\mu} \left[ 1 - \frac{(1-\varphi_S)e^{\mu T}}{e^{\mu T} - \varphi_S} e^{-\mu(t-nT)} \right], \\ S_2(t) = \frac{\Lambda}{\mu} \left[ 1 - \frac{(1-\varphi_S)e^{\mu T}}{e^{\mu T} - \varphi_S} e^{-\mu(t-nT)} \right], \end{array} \right.$$

with $S_1$ and $S_2$ the solutions of the following impulsive differential equations:

$$\dot{S}_1 = \Lambda - \omega\nu\varepsilon_U - \mu S_1, \quad t \ne nT; \quad S_1(nT^+) = \varphi_S S_1(nT),$$

$$\dot{S}_2 = \Lambda - \mu S_2, \quad t \ne nT; \quad S_2(nT^+) = \varphi_S S_2(nT),$$

One can observe that $S_1(t) \to S^P(t)$ if $\varepsilon_U \to 0$ and $S_2(t) \to S^P(t)$, which implies that $S(t) \to S^P(t)$ asymptotically. In the other words, independently from the initial value $(S_0, I_0, U_0) \in G$, one has that

$$(S(t), I(t), U(t)) \to (S^P(t), 0, 0)$$

This concludes the proof.

Numerical results for the global stability of the PDFS and persistence are present in Figure 1 and 2 respectively. To do so, we consider the initial conditions $(S_0, U_0) \in [(500, 200), (1000, 400), (1500, 600), (2000, 800)]$, $I_0 = 0$ and $B_0 = 0$. Figure 1 shows that, for different initial conditions, the number of infected branches and urediniospores converge towards zero, which is the state free disease. Additionally, we can observe that the trajectories of susceptible branches and berries becomes stationary after three seasons, as proved in Theorem 2. Although the persistence has not been proved analytically, numerical

<figure><img src="image_2.png" /><figcaption>Fig. 2. Trajectories of (1) for different initial condition, when $\gamma = 2.1$ (so that $\mathcal{R}_0 = 5.7 > 1$).</figcaption></figure>

results in Figure 2 tend to show that the persistence of the CLR appears for (1) when $\mathcal{R}_0 = 5.7 > 1$. From Figures 1 and 2 one can deduce that, during the disease free period (when $\mathcal{R}_0 < 1$), the number of berries is 378910, during the persistence of disease state (when $\mathcal{R}_0 > 1$) the number of berries is 117698. Then, one can conclude that the yield lost due to the CLR is estimate at 68,93%. Thus, one can ask the following question, how to reduce the damage caused by the CLR in coffee plantation.

# 4. OPTIMAL CONTROL

## 4.1 Optimal problem

Many different methods are commonly used in CLR control, with varying levels of environmental impact. Herein, we choose the biological control, because it is a potentially powerful tool for managing coffee leaf rust that allows organic certification. This biological control uses natural enemies, e.g. mycoparasites such as *Lecanicillium lecanii* (Carrion and Rico-Gray, 2002; Vandermeer et al., 2009). The goal is to introduce the mycoparasite in the plantation so that it will consume the urediniospores. This method represents the inoculative biological control and is denoted by $v$. The new system with control $v(t)$ is the first system (1) where the dynamical equation of $U$ is replaced by

$$\dot{U} = \gamma I - (\nu + \mu_U)U - v, \quad t \ne nT; \quad (13)$$

where $v(t) \in [0, v_{max}]$, with $v_{max}$ the maximum value of the biological control. Given $k$ seasons of duration $T$, we are interested in maximizing the net revenue, which is directly proportional to the amount of grains of coffee at the end of each season $\sum_{j=1}^k B(jT)$ and inversely proportional to the usage of the control treatment throughout the whole period $\int_0^{kT} v dt$. Thus, we define the cost function to maximize representing the profit for a time period $[0, kT]$, given by

$$J(v) = \alpha_B \sum_{j=1}^{k} B(jT) - \alpha_v \int_{0}^{kT} v dt,$$
where $\alpha_B$ and $\alpha_v$ are the market prices per coffee grain and per spore, respectively. By using the relation

$$B(jT) = \int_{(j-1)T}^{jT} (\delta_S S + \delta_I I - \mu_B B) dt,$$

which takes into account that $B(jT^+) = 0$, we can rewrite the cost function as

$$J(v) = \int_{0}^{kT} (\alpha_B(\delta_S S + \delta_I I - \mu_B B) - \alpha_v v) dt.$$

Then, we define the Hybrid Optimal Control Problem

$$\left\{
\begin{array}{ll}
\text{maximize} & J(v) \\
\text{subject to} & \text{dynamics of system (1) and (13),} \\
& \text{initial conditions (2),} \\
& v(\cdot) \in \mathcal{V},
\end{array}
\right.
\qquad (\text{HOCP})$$

with $\mathcal{V}$ being the set of admissible controllers, which are Lebesgue measurable real-valued functions defined on the interval $[0, kT]$ and satisfying the constraint $v(t) \in [0, v_{max}]$. We notice that $v(t)$ appears linearly in the dynamics (13), and, since $S = 0$ is repulsive, $S + I > 0$ for any trajectory, which means that the dynamics cannot blow up in finite time (and so all trajectories remain in a compact set). An analytical computation of an hybrid optimal control can be derived through the Hybrid Maximum Principle (Dmitruk and Kaganovich, 2008), a generalization of the classical PMP (Pontryagin's Maximum Principle). For (HOCP), we define the state $x = (S, I, U, B)$ and adjoint state $p = (p_S, p_I, p_U, p_B)$, and so the Hamiltonian is written as

$$H(x, p, p_0, v) = \langle \lambda, f(x, v) \rangle - p_0 f^0(x, v)$$

where $f^0(x, v) = \alpha_B \dot{B} - \alpha_v v$, which corresponds to the integrand of the cost function $J(v)$, $f$ is the right-hand side of (13), and $p_0 \le 0$. According to the Maximum Principle, the optimal control $v(t)$ should maximize the Hamiltonian along the optimal trajectory. Since $H$ is linear in the control, the solution of (HOCP) becomes

$$v_{opt}(t) = \begin{cases} 0 & \text{if } \phi(x, p) < 0; \\ v_{max} & \text{if } \phi(x, p) > 0; \\ v_{sing}(t) & \text{if } \phi(x, p) = 0. \end{cases} \qquad (14)$$

where $\phi(x, p)$ is the switching function defined as $\phi(x, p) \doteq \partial H/\partial v$, and $v_{sing}(t)$ is called a singular control. According to (14), any solution of (HOCP) is a concatenation of bangs and singular arcs. A singular arc occurs when the switching function $\phi$ vanishes over some subinterval of time, and can be obtained by successively differentiating $\phi$ until the control $v$ can be computed explicitly in terms of the state and adjoint state. However, a more detailed description of the control law is often hard to obtain, due to the presence of the adjoint state $p$ in the singular arc $v_{sing}$ and in the switching function $\phi(x, p)$. One way to help understand the structure of the optimal control is to simulate numerical trajectories using direct methods, by approximating the problem by a finite dimensional optimization problem. Nevertheless, the use of optimal control solvers is usually restricted to continuous systems. Thus, in order to be able to apply such methods, we rewrite (HOCP) in a continuous form.

## 4.2 Problem in continuous form

(HOCP) can be rewritten in the continuous form by defining $k$ set of states $x_n = (S_n, I_n, U_n, B_n)$ and $k$ controls

$v_n(t)$, with dynamics defined for $t \in [0, T]$ as

$$\left\{
\begin{aligned}
\dot{S}_n &= \Lambda - \omega \nu U_n \frac{S_n}{S_n + I_n} - \mu S_n, \\
\dot{I}_n &= \omega \nu U_n \frac{S_n}{S_n + I_n} - (\mu + d) I_n, \\
\dot{U}_n &= \gamma I_n - (\nu + \mu_U) U_n - v_n, \\
\dot{B}_n &= \delta_S S_n + \delta_I I_n - \mu_B B_n,
\end{aligned}
\right.
\qquad (15)$$

for $n = 1, \dots, k$, with boundary conditions

$$\left\{
\begin{aligned}
S_1(0) &= S_0 > 0, & I_1(0) &= I_0 \ge 0, \\
U_1(0) &= U_0 \ge 0, & B_1(0) &= B_0 > 0 \\
S_{n+1}(0) &= \varphi_s S_n(T), & I_{n+1}(0) &= \varphi_I I_n(T), \\
U_{n+1}(0) &= \varphi_U U_n(T), & B_{n+1}(0) &= 0.
\end{aligned}
\right.
\qquad (16)$$

for $n = 1, \dots, k-1$. Then, the cost function becomes

$$J_c(v_1, \dots, v_k) = \alpha_B \sum_{n=0}^{k} B_n(T) - \alpha_v \sum_{n=0}^{k} \int_{0}^{T} v_n dt,$$

and the resulting continuous optimal control problem is

$$\left\{
\begin{array}{ll}
\text{maximize} & J_c(v_1, \dots, v_k) \\
\text{subject to} & \text{dynamics of system (15),} \\
& \text{boundary conditions (16),} \\
& v_n(\cdot) \in \mathcal{V}, \text{ for } n = 1, \dots, k.
\end{array}
\right.
\qquad (\text{COCP})$$

While this new OCP includes the boundary conditions relating the consecutive set of states $x_n$, there are no terminal conditions on (HOCP) which means that the set of admissible controllers is not empty. Thus, existence of a solution for (COCP)—and therefore for (HOCP)—holds by Filippov's theorem (Agrachev and Sachkov, 2013).

## 4.3 Numerical results

Figure 3 presents the results of a numerical simulation performed with Bocop of the optimal control and the optimal trajectory for problem (HOCP). In this Figure, we compare an open-loop strategy ($v \equiv 0$) with $R_0 > 1$, and the same initial conditions with the optimal control. For the latter case, the biological treatment successfully erradicates the pest. The optimal strategy is characterized by a bang control $v \equiv v_{max}$ during an initial period of time in the first season, followed by a singular arc in which the use of the pest treatment is gradually reduced until no further control is applied. Indeed, the optimal control exploits the fact that both $U$ and $I$ are greatly reduced at the end of each season due to harvesting, to economize on pest treatment. After that, no further control is required. The latter pattern has also been observed for different set of parameters and initial conditions defined in the previous section. In general, the treatment is applied rather for an initial interval of the season, followed by a decrease that takes the control to $v \equiv 0$. As a conclusion for the analyzed particular case, coffee producers should maximize the quantity of mycoparasitic (between given boundaries) roughly only during the beginning of the first year, while discontinuing the treatment in posterior seasons. Naturally, the instant at which the treatment should be turned off is strongly linked to the initial conditions, and so each case requires a specific analysis.

# 5. CONCLUSION

This paper presented a comprehensive, impulsive model for the transmission dynamics of the CLR within a coffee
<figure><img src="image_3.png" /><figcaption>Fig. 3. Trajectories of HCOP (13) when $\gamma = 2.1$ (so that $R_0 = 5.7$ without control). The maximum control $v_{max} = 6.5$. Initial conditions are set to $S(0) = 500$, $U(0) = 200$, $I(0) = B(0) = 0$. In the case of the optimal trajectory, $I = U = 0$ for all $t$.</figcaption></figure>

plantation. The model has been analyzed to gain insight into its qualitative dynamics. Using the Floquet theory, we have found that there exists a threshold parameter that determines the extinction and persistence of CLR in a coffee plantation. Later on, we presented a controlled model based on an inoculative biological control, for which we defined an OCP targeting the maximization of the net revenue of a coffee plantation. For such problem, we performed a numerical analysis for a relevant set of parameters and initial conditions, intended to provide guidance in how to control the damages caused by CLR. In a follow-up of this work, we will perform a more detailed pest treatment analysis by proposing different control strategies including the dynamics of the mycoparasite in the model.

## ACKNOWLEDGEMENTS

We thank Pr Samuel Bowong for comments that greatly improved the manuscript.

## REFERENCES

Abbasi, Z., Zamani, I., Mehra, A.H.A., Shafieirad, M., and Ibeas, A. (2020). Optimal control design of impulsive sqeiar epidemic models with application to covid-

19. *Chaos, Solitons and Fractals*, 139, 110054. doi: https://doi.org/10.1016/j.chaos.2020.110054.

Agrachev, A.A. and Sachkov, Y. (2013). *Control theory from the geometric viewpoint*, volume 87. Springer Science & Business Media.

Belbas, S. and Schmidt, W. (2009). Optimal control of impulsive volterra equations with variable impulse times. *Applied mathematics and computation*, 214(2), 353-369.

Carrion, G. and Rico-Gray, V. (2002). Mycoparasites on the coffee rust in Mexico. *Fungal Diversity*, 11(5), 49-60.

Coppel, W.A. (1965). *Stability and asymptotic behavior of differential equations*. Heath.

Djuikem, C., Grognard, F., Wafo, R.T., Touzeau, S., and Bowong, S. (2021). Modelling coffee leaf rust dynamics to control its spread. *Mathematical Modelling of Natural Phenomena*.

Dmitruk, A.V. and Kaganovich, A.M. (2008). The hybrid maximum principle is a consequence of pontryagin maximum principle. *Systems & Control Letters*, 57(11), 964-970.

d'Onofrio, A. (2002). Stability properties of pulse vaccination strategy in seir epidemic model. *Mathematical biosciences*, 179(1), 57-72.

Hale, J. (1969). *Ordinary differential equations*, 1969. Wiley-Interscience, New York.

International Coffee Organization (2020). Consumption of coffee. Trade statistics tables: https://cutt.ly/Lg8R5X9 (accessed 13/11/2020).

Mailleret, L. and Grognard, F. (2009). Global stability and optimisation of a general impulsive biological control model. *Mathematical Biosciences*, 221(2), 91-100.

Pontryagin, L.S. (2018). *Mathematical theory of optimal processes*. Routledge.

Tamen, A.T., Dumont, Y., Tewa, J.J., Bowong, S., and Couteron, P. (2017). A minimalistic model of tree-grass interactions using impulsive differential equations and non-linear feedback functions of grass biomass onto fire-induced tree mortality. *Mathematics and Computers in Simulation*, 133, 265-297.

Tankam-Chedjou, I., Touzeau, S., Mailleret, L., Tewa, J.J., and Grognard, F. (2020). Modelling and control of a banana soilborne pest in a multi-seasonal framework. *Mathematical Biosciences*, 322, 108324.

Team Commands, I.S. (2017). Bocop: an open source toolbox for optimal control. http://bocop.org.

Vandermeer, J., Hajian-Forooshani, Z., and Perfecto, I. (2018). The dynamics of the coffee rust disease: an epidemiological approach using network theory. *European Journal of Plant Pathology*, 150(4), 1001-1010.

Vandermeer, J., Perfecto, I., and Liere, H. (2009). Evidence for hyperparasitism of coffee rust (Hemileia vastatrix) by the entomogenous fungus, Lecanicillium lecanii, through a complex ecological web. *Plant Pathology*, 58(4), 636-641.

Wang, H. (2012). *Mathematical modeling I-preliminary*. Bookboon.

Wong, K.H. and Tong, W.M. (2012). Optimal control of switched impulsive systems with time delay. *The ANZIAM Journal*, 53(4), 292-307. doi: 10.1017/S1446181112000284.