Optimal Spraying in Biological Control of Pests

TL;DR

This paper applies optimal control theory to determine the most effective spraying strategies for pest management in agroecosystems, aiming to minimize or eradicate pest populations.

Contribution

It introduces two novel optimal control formulations for pest spraying, providing insights into effective biological control strategies.

Findings

Optimal control solutions suggest effective spraying policies.

Strategies can potentially eradicate pests.

Results inform sustainable pest management practices.

Abstract

We use optimal control theory with the purpose of finding the best spraying policy with the aim of at least to minimize and possibly to eradicate the number of parasites, i.e., the prey for the spiders living in an agroecosystems. Two different optimal control problems are posed and solved, and their implications discussed.

Tables

Table 1. Table 1. Model parameters: the reference values are taken from Ezio:spiders:JNAIAM:2008 .

Symbol	Description	Value
$r$	insects’s natural birth rate	1
$e$	parasites’s natural birth rate	2.5
$a$	spider’s natural mortality	3.1
$b$	hunting rate of spiders on vineyard insects	1.2
$c$	hunting rate of spiders on wood insects	0.2
$h$	intensity of spraying	0.7
$q$	fraction of insecticide that lands on vineyards	0.9
$k$	“efficiency” constant at which prey biomass is turned into new spiders	1
$W$	carrying capacity for insects in the woods	5
$V$	carrying capacity for insects in the vineyards	1000
$K$	unwanted killing effect of spraying on spiders	0.01

Equations

$$\{\begin{array}{l}\dot{f}(t)=rf(t)\left(1-\frac{f(t)}{W}\right)-cs(t)f(t)-h(1-q)u(t),\\ \dot{s}(t)=s(t)\left(-a+kbv(t)+kcf(t)\right)-hKqu(t),\\ \dot{v}(t)=ev(t)\left(1-\frac{v(t)}{V}\right)-bs(t)v(t)-hqu(t).\end{array}$$

$$f(0)=f_0,s(0)=s_0,v(0)=v_0,\text{with}{f}_0,s_0,v_0\ge 0.$$

$$\Omega :=\left\{u(\cdot )\in L^{\infty }(0,T)|0\le u(t)\le 1,\forall t\in [0,T]\right\}.$$

$$\mathcal{J}(u(\cdot ))={\int }_0^T\left[v(t)+\frac{\xi }{2}{u}^2(t)\right]dt.$$

$$\mathcal{J}(u^{\ast }(\cdot ))=\underset{\Omega }{\min }\mathcal{J}(u(\cdot )).$$

$$\dot{f}=\frac{\partial H}{\partial {\lambda }_1},\dot{s}=\frac{\partial H}{\partial {\lambda }_2},\dot{v}=\frac{\partial H}{\partial {\lambda }_3}$$

$${\dot{\lambda }}_1=-\frac{\partial H}{\partial f},{\dot{\lambda }}_2=-\frac{\partial H}{\partial s},{\dot{\lambda }}_3=-\frac{\partial H}{\partial v},$$

$$\begin{gathered} H=H(f,s,v,u,{\lambda }_1,{\lambda }_2,{\lambda }_3)=v+\frac{\xi }{2}{u}^2+{\lambda }_1\left(rf\left(1-\frac{f}{W}\right)-csf-h(1-q)u\right) \\ +{\lambda }_2\left(s\left(-a+kbv+kcf\right)-hKqu\right)+{\lambda }_3\left(ev\left(1-\frac{v}{V}\right)-bsv-hqu\right) \end{gathered}$$

$$H(f^{\ast }(t),s^{\ast }(t),v^{\ast }(t),u^{\ast }(t),\lambda (t))=\underset{0\le u\le 1}{\min }H(f^{\ast }(t),s^{\ast }(t),v^{\ast }(t),u,\lambda (t))$$

$$\{\begin{array}{l}{\dot{\lambda }}_1(t)=-{\lambda }_1(t)\left(r\left(1-\frac{f^{\ast }(t)}{W}\right)-\frac{r{f}^{\ast }(t)}{W}-c{s}^{\ast }(t)\right)-{\lambda }_2(t)s^{\ast }(t)kc,\\ {\dot{\lambda }}_2(t)={\lambda }_1(t)c{f}^{\ast }(t)-{\lambda }_2(t)\left(-a+kb{v}^{\ast }(t)+kc{f}^{\ast }(t)\right)+{\lambda }_3(t)b{v}^{\ast }(t),\\ {\dot{\lambda }}_3(t)=-1-{\lambda }_2(t)s^{\ast }(t)kb-{\lambda }_3(t)\left(e\left(1-\frac{v^{\ast }(t)}{V}\right)-\frac{e{v}^{\ast }(t)}{V}-b{s}^{\ast }(t)\right),\end{array}$$

$$u^{\ast }(t)=\min \left\{\max \left\{0,\frac{h\left({\lambda }_1(t)-{\lambda }_1(t)q+{\lambda }_2(t)Kq+{\lambda }_3(t)q\right)}{\xi }\right\},1\right\}.$$

$$f(0)=3.1,s(0)=3.7,v(0)=2.2,$$

$${\int }_0^{T}v(t)dt.$$

$$\mathcal{J}={\int }_0^{T}1dt.$$

$$u^{\ast }(t)=\min \left\{\max \left\{0,u_{sing}^{\ast }(t)\right\},1\right\},$$

$$\begin{array}{r}{u}_{sing}^{\ast }(t)=-\frac{A(t)}{2B(t)}\end{array}$$

$$\begin{array}{rl}A& =q{b}^2{s}^2{W}^2{V}^2({\lambda }_3-{\lambda }_{2}k)+W^2{V}^{2}q(e-bs)+2q{\lambda }_{3}ev{W}^{2}bsV+2{\lambda }_1{r}^2{f}^2{V}^2+W^2{V}^2(r^2{\lambda }_1\\ & +2csq{\lambda }_{1}r+c^2{s}^{2}q{\lambda }_{2}k)+2csW{V}^2{\lambda }_{1}rf+W^2{V}^2(-2qe{\lambda }_{3}bs+r{\lambda }_{2}skc+Kq{a}^2{\lambda }_2-{\lambda }_{1}csa\\ & -c^2{s}^{2}q{\lambda }_1)+2qVW({\lambda }_1{r}^{2}fV-{\lambda }_3{e}^{2}vW)+W^2{V}^2(-2cs{\lambda }_{1}r-c^2{s}^2{\lambda }_{2}k-Kqbv)\\ & -2({\lambda }_1{r}^2{f}^{2}q{V}^2+q{\lambda }_3{e}^2{v}^2{W}^2)+W^2{V}^2(-r^{2}q{\lambda }_1+c^2{s}^2{\lambda }_1+q{e}^2{\lambda }_3)\\ & -2qev{W}^{2}V-2{\lambda }_1{r}^{2}f{V}^{2}W-W^2{V}^{2}skc{\lambda }_{3}bv+W^2{V}^{2}Kqa({\lambda }_{1}cf-2{\lambda }_{2}kbv-2{\lambda }_{2}kcf+{\lambda }_{3}bv)\\ & +{\lambda }_1{W}^2{V}^{2}cskbv-2csqW{V}^2{\lambda }_{1}rf+W^2{V}^2{\lambda }_{2}sk(qeb-rqc)+2{\lambda }_{2}sk(rfq{V}^{2}cW-evq{W}^{2}bV\\ & -rf{V}^{2}cW)+W^2{V}^{2}Kq(-kbv{\lambda }_{1}cf+k^2{b}^2{v}^2{\lambda }_2+2k^{2}bv{\lambda }_{2}cf-k{b}^2{v}^2{\lambda }_3-k{c}^2{f}^2{\lambda }_1+k^2{c}^2{f}^2{\lambda }_2\\ & -kcf{\lambda }_{3}bv+sc{\lambda }_{3}bv)+W^2{V}^2{\lambda }_1(-qskbcf+csqa-csqkbv)+q{\lambda }_{3}bs{W}^2{V}^2(-a+kcf)\\ & +Kq(c{f}^{2}W{V}^2{\lambda }_{1}r+b{v}^2{W}^{2}V{\lambda }_{3}e+{\lambda }_{2}k(-bv{W}^2{V}^{2}e+b{v}^2{W}^{2}Ve-cf{W}^2{V}^{2}r+c{f}^{2}W{V}^{2}r)),\end{array}$$

$$\begin{array}{rl}B& =hWV[{\lambda }_{2}KWVkcq(1-q)+{\lambda }_{1}V(2rq-r{q}^2-cKqW+cK{q}^{2}W-r)\\ & +KWV{q}^{2}b({\lambda }_{2}k-{\lambda }_3)-q^2{\lambda }_{3}eW],\end{array}$$

$$\{\begin{array}{l}{\dot{\lambda }}_1(t)=-{\lambda }_1(t)\left(r\left(1-\frac{f^{\ast }(t)}{W}\right)-\frac{r{f}^{\ast }(t)}{W}-c{s}^{\ast }(t)\right)-{\lambda }_2(t)s^{\ast }(t)kc,\\ {\dot{\lambda }}_2(t)={\lambda }_1(t)c{f}^{\ast }(t)-{\lambda }_2(t)\left(-a+kb{v}^{\ast }(t)+kc{f}^{\ast }(t)\right)+{\lambda }_3(t)b{v}^{\ast }(t),\\ {\dot{\lambda }}_3(t)=-{\lambda }_2(t)s^{\ast }(t)kb-{\lambda }_3(t)\left(e\left(1-\frac{v^{\ast }(t)}{V}\right)-\frac{e{v}^{\ast }(t)}{V}-b{s}^{\ast }(t)\right).\end{array}$$

$$H(f^{\ast }(t),s^{\ast }(t),v^{\ast }(t),u^{\ast }(t),\lambda (t),{\lambda }^0)=\underset{0\le u\le 1}{\min }H(f^{\ast }(t),s^{\ast }(t),v^{\ast }(t),u,\lambda (t),{\lambda }^0)$$

$$\begin{array}{rl}H& =H(f,s,v,u,\lambda ,{\lambda }^0)\\ & ={\lambda }^0+{\lambda }_1\left(rf\left(1-\frac{f}{W}\right)-csf-h(1-q)u\right)\\ & +{\lambda }_2\left(s\left(-a+kbv+kcf\right)-hKqu\right)\\ & +{\lambda }_3\left(ev\left(1-\frac{v}{V}\right)-bsv-hqu\right).\end{array}$$

$$u^{\ast }(t)=\min \left\{\max \left\{0,u_{sing}^{\ast }\right\},1\right\},$$

$$h(1-q)-{\lambda }_{2}hKq-{\lambda }_{3}hq=0,$$

$$a<ckW\text{ and }{c}^{2}kWe<br(ckW-a)$$

Full text

1. Introduction

Fairly recent researches have highlighted the role that generalist predators, i.e., predators that feed on several species, have on pests in agriculture, so that their presence in agroecosystems should be fostered RL . For instance, it is well known that toward the end of the nineteenth century, the European vineyards Vitis vinifera suffered from the accidental introduction of the plant louse phylloxera Pinney:HW . Later vineyard agroecosystems in every continent became affected by this spread of the pest. This fact caused relevant economic damages worldwide Winemaking . The crisis was overcome by the use of native American grapevines, grafting the European vines onto Eastern American roots, such as Vitis labrusca, which are resistant to phylloxera Ox:Comp:Wine . The role of spiders to combat these and other parasites affecting vines, such as black measles, little-leafs, nematodes and red ticks, and infesting fruit orchards, e.g., grapevine beetles, grape-berry moths, climbing cutworms, black rot and mildew, has been recently elucidated in the literature CD98 , IBBB , NM , W . The findings of these researches indicate that vineyards contain a very diversified and abundant spider community, which depends heavily on the availability of the grass lying in between the vine rows and possibly nearby bushes or small woods, and as a whole helps in keeping down the parasites population RB . However, spiders have different hunting strategies and locations, for which different insects are most commonly hunted by the various spider species MC . Spiders can be distinguished in two very broad sets, the web-builder species, which are in general stantial, except for a single “flight” due to the action of the wind when they are young. Then they build the web and wait for the prey. The second set is composed by wanderer species, like wolf spiders, that crawl on the ground in search for food. For instance, the nocturnal wandering spiders Clubiona brevipes, C. corticalis and C. leucaspis (Clubionidae) hunt non-flying Aphids and larvae of Lepidoptera, while the diurnal wandering species Ballus depressus (Salticidae) are effective against Cicadellidae; adults and larvae of Hymenoptera and Lepidoptera are prey of the ambush species Philodromus aureolus (Philodromidae) and Diaea dorsata (Thomisidae) IBB06 , Ma . The ultimate result is that we can consider spiders as specialists, i.e., feeding essentially only on one species, thus implying that parasites can be controlled only by ensuring the persistence of a wide variety of diverse spider species within the agroecosystem. Mathematical models for the understanding of these complex systems have been proposed and analysed in the past years: in VIBITB a space-free description of the dispersal phenomenon of web-builder spiders is presented; wanderer spiders are instead considered in VIBCB and its extensions CIVa , CIVb .

Another very common strategy to fight parasites is represented by the use of pesticides, although it is now widely recognized that they have high economical costs and negative impact in terms of both environment damages as well as human health. Besides, pests acquire resistance to the chemicals, which may lead to unwanted large parasites outbreaks RD-B .

The action of parasites has been also considered within the modeling framework of previous researches, but always considering very simple control strategies, as the focus was mainly on extracting the relationships between the various populations living in the environment, while the human role within the ecosystem model was considered marginal. In fact, either an instantaneous effect of the spray, VIBITB , or an exponentially decay of the poison, CIVb , VIBCB , were considered, but in both cases the control was taken to be a constant function.

Optimal control theory has a long history of being applied to problems in biomedicine (see, e.g., Eisen1979 , Ledzewicz2002 , Lenhart2007 and references cited therein) and recently to epidemic models for infectious diseases (see, e.g., Behncke2000 , Gaff2009 , Sofia2010 , Sofia2013 , PaulaSilvaTorresT2014 , SilvaTorresNACO , SilvaTorresMBS2013 , Silva:Torres:DCDS-A:2015 for population compartmental models, HattafHIV2012 , HattafHIVdelay2012 for HIV cellular-level models, and IbrahimIJDC2016 for numerical methods applied to optimal control problems). However, to our knowledge, little attention has been given to agroecosystem models.

In this paper, we apply optimal control theory to the activities related to an agroecosystem, and study specifically the possible spraying policies of the fruit orchards. Two optimal control problems are identified and solved in this context, providing optimal strategies for the realistic implementation of pest eradication: the complete elimination of the pests, the minimum time to achieve the previous goal and the minimization of the parasites in a given time frame.

The text is organized as follows. In Section 2, we recall the model and extend it as a general control system. In Section 3, two optimal control problems are proposed and solved: subsection 3.1 addresses the pest eradication problem, subsection 3.2 contains the time optimal control problem. We end with Section 4 of discussion and Section 5 of conclusions and future work.

2. Mathematical model with control function

We consider the mathematical model from Ezio:spiders:JNAIAM:2008 , which considers interactions between three populations: insects living in open fields and woods at time $t$ ($f(t))$; parasites living in the vineyards ($v(t)$); and wanderer spiders living in the whole environment ($s(t)$). Both insects $f(t)$ and parasites $v(t)$ can constitute prey for the spiders $s(t)$. The time unit is the day and populations are counted in numbers.

Let $W$ and $V$ denote the carrying capacity for insects in the woods and vineyards, respectively. Given the situation of highly exploited lands, it is assumed that the vineyard carrying capacity much exceeds the one of nearby green areas, so that $V\gg W$.

We assume that insect and parasite populations ($f$ and $v$) are the only food source for spiders $w$, so that when the former are absent the latter will die exponentially fast. Predation is accounted for by mass action law terms with suitable signs. Captured prey contribute to spider reproduction, however not the whole prey is consumed but only a fraction of it; this is expressed via the efficiency constant $0<k<1$.

Insects and parasites natural birth rates are denoted by $r$ and $e$, respectively. The spiders’ natural mortality is denoted by the parameter $a$. The hunting rates of spiders on vineyard and woods are denoted by $b$ and $c$, respectively. The parameter $h$ represents the intensity of spraying, i.e., the amount of poison released, $K$ models its unwanted killing effect on spiders, $q$ represents the fraction of insecticide that lands on target, i.e., within the vineyards, and $1-q$ the corresponding fraction that due to the action of the wind, water and possibly other causes, is dispersed in the woods. All the parameters are nonnegative constants unless otherwise specified.

We rewrite the most general form of the model in Ezio:spiders:JNAIAM:2008 including human spraying of the fields, by replacing its particular controls by a general function. Specifically, we add to the dynamic system a control variable, $u$, which is a function of time and not a very specific control as assumed in the previous literature Ezio:spiders:JNAIAM:2008 , CIVb , VIBCB , VIBITB . It represents the human intervention via spraying with insecticides, as follows:

[TABLE]

3. Problem formulation

In this section, we formulate two optimal control problems where the control system is given by (2.1). These problems are subject to initial conditions

[TABLE]

Note that for $u\equiv 0$ the control system (2.1) reduces to the model for the ecosystem without human intervention. The set $\Omega$ of admissible control functions is defined as

[TABLE]

The optimal control problems are the following:

1. (1)

minimize the number of parasites ($v(t)$) with human intervention (use of insecticides), taking or not into account the cost of insecticides (see subsection 3.1); 2. (2)

minimize the time $T$ for which the parasites of the vineyard will be eradicated, that is, $v(T)=0$ with $f(T)$ and $s(T)$ free and $s(T)>0$ (see subsection 3.2).

For the two problems, we compare the results with human intervention (with control, that is, $u\neq 0$) to the case where there is no human intervention (no control/no insecticides, that is, $u\equiv 0$).

3.1. Minimizing the number of parasites

Consider the following objective functional defined by a sum of two terms:

[TABLE]

The first term in our cost functional tells us that we want to minimize the number of pests. The constant $\xi$ in the second term is a measure of the relative cost of the intervention associated to the control $u$. One applies control measures that are associated with some implementation costs that we also intend to minimize. By considering the cost with the control in a quadratic form, we are being consistent with previous works in the literature (see, e.g., PaulaSilvaTorresT2014 , SilvaTorresMBS2013 ). Moreover, a quadratic structure in the control has mathematical advantages. Roughly speaking, it implies that the Hamiltonian attains its minimum over the control set at a unique point (given by (3.7)). The optimal control problem consists of determining the triple $\left(f^{*}(\cdot),s^{*}(\cdot),v^{*}(\cdot)\right)$, associated to an admissible control $u^{*}(\cdot)\in\Omega$ on the time interval $[0,T]$, satisfying (2.1), the initial conditions (3.1), and minimizing the cost functional (3.2), that is,

[TABLE]

The existence of an optimal control $u^{*}$ and associated $(f^{*},s^{*},v^{*})$ comes from the convexity of the integrand of the cost functional (3.2) with respect to the control $u$ and the Lipschitz property of the state system with respect to state variables $(f,s,v)$ (see, e.g., Cesari_1983 , Fleming_Rishel_1975 for existence results of optimal solutions).

According to the Pontryagin Maximum Principle Pontryagin_et_all_1962 , if $u^{*}(\cdot)\in\Omega$ is optimal for the problem (2.1), (3.3) with the initial conditions given by (3.1) and fixed final time $T$, then there exists a nontrivial absolutely continuous mapping $\lambda:[0,T]\to\mathbb{R}^{3}$, $\lambda(t)=\left(\lambda_{1}(t),\lambda_{2}(t),\lambda_{3}(t)\right)$, called adjoint vector, such that

[TABLE]

and

[TABLE]

where the function $H$ defined by

[TABLE]

is called the Hamiltonian, and the minimality condition

[TABLE]

holds almost everywhere on $[0,T]$. Moreover, one has the transversality conditions $\lambda_{i}(T)=0$, $i=1,\ldots,3$.

{thrm}

The problem (2.1), (3.3) with fixed initial conditions (3.1) and fixed small final time $T$, admits a unique optimal solution $(f^{*}(\cdot),s^{*}(\cdot),v^{*}(\cdot))$ associated with an optimal control $u^{*}(\cdot)$ on $[0,T]$. Moreover, there exist adjoint functions, $\lambda_{1}(\cdot)$, $\lambda_{2}(\cdot)$ and $\lambda_{3}(\cdot)$, such that

[TABLE]

with transversality conditions $\lambda_{i}(T)=0$, $i=1,\ldots,3$. Furthermore,

[TABLE]

Proof.

Existence of an optimal solution $(f^{*},s^{*},v^{*})$ associated to an optimal control $u^{*}$ comes from the convexity of the integrand of the cost functional $J$ with respect to the control $u$ and the Lipschitz property of the state system with respect to the state variables $\left(f,s,v\right)$ (see, e.g., Cesari_1983 , Fleming_Rishel_1975 ). System (3.6) is derived from the Pontryagin maximum principle (see (3.4), Pontryagin_et_all_1962 ) and the optimal controls (3.7) come from the minimization condition (3.5). For small final time $T$, the optimal control pair given by (3.7) is unique due to the boundedness of the state and adjoint functions and the Lipschitz property of systems (2.1) and (3.6) (see SLenhart_2002 and references cited therein). ∎

In this article the numerical results were obtained using the PROPT Matlab Optimal Control Software PROPT . We considered the initial conditions

[TABLE]

and the parameter values from Table 1.

If we consider the model without controls, i.e., (2.1) with $u\equiv 0$, then the behavior for $t\in[0,50]$ is shown in Figure 1.

Next, we consider human intervention, that is, the situation when spraying with insecticides is considered ($u\neq 0$).

3.1.1. Insecticides cost not accounted for

We begin by considering $\xi=0$, that is, the goal is simply to minimize

[TABLE]

The solution to this optimal control problem is shown in Figures 2 and 3 for $T=50$; and in Figures 4 and 5 for $T=150$.

When we do not consider the costs associated to the control (the spray), the number of pests in the vineyard $v(t)$ can attain the value six for $T$ near 50 days (see Figure 2). However, between days 20 and 30 the number $v(t)$ is lower than the one observed in Figure 1, when no control is applied. Since there is no cost associated with the spray, the optimal control attains the upper bound often (see Figure 3).

In Figures 4 and 5 we consider $\xi=0$ and $T=150$. We observe that the number of preys in the vineyard attains a value around 2 during almost 100 days, which is lower than the one observed when no spray is applied (compare Figures 1 and 4).

However, this is associated with a control that takes constantly the maximum value for more than 50 days (see Figure 5).

3.1.2. Inclusion of the cost of insecticides

Proceeding differently from Section 3.1.1, now we also account for the cost of insecticides ($\xi\neq 0$) in the objective functional to be minimized. More precisely, we set $\xi=50$ in (3.2). We observe that the coexistence equilibrium is attained at $t\simeq 20$ (see Figure 6). The optimal way of spraying $u^{*}(t)$, $0\leq t\leq 50$, is illustrated in Figure 7.

Comparing the range of oscillations of $(f(\cdot),s(\cdot),v(\cdot))$ for $u\equiv 0$ (no spraying with insecticide) with those of $(f^{*}(\cdot),s^{*}(\cdot),v^{*}(\cdot))$ corresponding to the optimal control $u^{*}$, we observe that in case of human intervention (the case with control $u^{*}$) the range of oscillations at $t\simeq 26$ is similar to the case without control at $t\simeq 150$ (Figures 6 and 8).

In contrast with Section 3.1.1, now we take into account the costs associated with the spray. For this reason the optimal control never attains the maximum allowed value, being always less than $0.4$. Moreover, for $t>20$ days, the optimal control takes values close to zero, which means that very small quantities of insecticides are applied (see Figure 7).

We observe that if no insecticide is applied, one needs almost 150 days to attain the values for $v(t)$ that we have at the end of 25 days with insecticides (compare Figures 6 and 8).

3.2. A time-optimal control problem

We now consider the objective functional in the form

[TABLE]

Our aim is to minimize the time $T$ for which the number of parasites living in a vineyard will vanish: $v(T)=0$. More precisely, the optimal control problem consists in finding the control $u^{*}$ that attains the final point ($f(T)=f_{T}$, $s(T)=s_{T}$, $v(T)=v_{T}$) in minimal time, where $f_{T}$, $s_{T}>0$ are free and $v_{T}=0$. In this section, we make the following controllability assumption:

$(H)$ The final point $\left(f_{T},s_{T},0\right)$, with $f_{T}$, $s_{T}>0$, can be steered from the initial point (3.1).

{thrm}

Under assumption $(H)$, the optimal control problem (2.1), (3.9) admits a solution $(f^{*}(\cdot),s^{*}(\cdot),v^{*}(\cdot))$ on $[0,T]$, associated to a control $u^{*}(\cdot)$, where $T>0$ is the minimal time. Moreover,

[TABLE]

where $u^{*}_{sing}$ denotes the singular control given by

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

Proof.

The existence of a solution $(f^{*}(\cdot),s^{*}(\cdot),v^{*}(\cdot))$, associated to a control $u^{*}(\cdot)$ on $[0,T]$, where $T>0$ is the minimal time, comes from the assumption $(H)$ (see, e.g., [Cesari_1983, , Chapter 9] for optimal control existence theorems). According to the Pontryagin Maximum Principle Pontryagin_et_all_1962 , there exist a real number $\lambda^{0}\geq 0$ and an adjoint function $\lambda(\cdot)=(\lambda_{1}(\cdot),\lambda_{2}(\cdot),\lambda_{3}(\cdot))\,:\,[0,T]\rightarrow\mathbb{R}^{3}$, with $\left(\lambda_{1}(\cdot),\lambda_{2}(\cdot),\lambda_{3}(\cdot),\lambda^{0}\right)\neq 0$, such that (3.10) holds together with the minimality condition

[TABLE]

holding almost everywhere on $[0,T]$. Here

[TABLE]

Moreover, $\min_{0\leq u\leq 1}H(f^{*}(t),s^{*}(t),v^{*}(t),u,\lambda(t),\lambda^{0})=0$ for every $t\in[0,T]$. From (3.11) we have

[TABLE]

where $u^{*}_{sing}$ denotes a singular control that can be obtained by differentiating twice $\frac{\partial H}{\partial u}=0$, i.e.,

[TABLE]

and replacing each derivative $\dot{\lambda}_{i}$, $i=1,2,3$, by (3.10), and $(\dot{f}^{*},\dot{s}^{*},\dot{v}^{*})$ by (2.1). ∎

For numerical simulations, we considered $e=0.3$ and the rest of the parameter values from Table 1, which satisfy the conditions for the feasibility and stability of the equilibrium $(f_{4},s_{4},0)$, that is,

[TABLE]

(see [Ezio:spiders:JNAIAM:2008, , pag. 197]). As for the initial conditions, we considered the same (3.8) as in the previous section. Figure 9 shows the optimal control, and Figure 10 the optimal state variables $(f^{*},s^{*},v^{*})$. In Figure 11 we can observe the solutions of model (2.1) without control, i.e., $u\equiv 0$ with initial conditions (3.8) and parameter values from Table 1 with the exception of $e=0.3$.

4. Discussion

Comparing the no-control policy with the control for minimization of the total number of pests, but not accounting for the costs of the pesticides, we find that in the latter case both spiders and pests experience outbursts during the last part of the simulation, higher than the peaks of the oscillations arising naturally in the absence of the control, compare Figures 1 and 2.

On the same problem, we replicated the simulation over a larger timeframe, $T=150$. In this case we note that the resulting control is different. Further, it allows again population peaks toward the final time, but their sizes now are comparable to the ones of the system without control, compare Figures 3 and 1. For earlier times instead, the populations remain at lower values. This apparently indicates that the longer the timespan of the control application, the lower the possible populations outbursts are.

Taking into account the costs of spraying, it is immediately seen that all the populations after a relatively initial short transient, up to $t=20$, are now confined to small oscillations about the average value that they attain when the control is not administered, compare Figures 6 and 1. Note also that the corresponding control is mainly used at the early times, and is essentially “negligible” after time $t=20$, see Figure 7. Further, the optimal spraying policy uses all the time only at most $40\%$ of the available insecticide. While in the previous cases, Figure 3 and, above all, Figure 5, the control is maximal and used constantly for a good portion of the time interval, here it is administered only at a few selected instants, compare the peaks of Figure 7 with the flat portion of Figure 5. The same result is achieved in Figure 8, where no control is used, but it takes much more time, in fact populations fall within the same range only at time $t=150$, rather than $t=20$ as found above.

Finally, comparing Figures 10 and 11, we see that the minimal time under spraying is about a hundredfold smaller than without the insecticide usage. The optimal control, Figure 9, needs to be used only for the first $0.6$ units of time, time at which the pests are eradicated, after which the system can keep to evolve naturally.

5. Conclusion and future work

In this paper, we have considered two optimal control problems for a model of pest management in agroecosystems. More precisely, we studied problems of how to spray pests with the aim to eradicate the parasites. We proved necessary optimality conditions and obtained numerical solutions for each of the problems using the optimal control solver PROPT PROPT . The problems considered in Section 3 are challenging, both from analytical and numerical points of view. Indeed, in Section 3.1.1, we consider a problem formulation that falls into the category of $L^{1}$-optimization. Numerical algorithms are especially poor in locating the optimal solutions for $L^{1}$-type problems, which is the main reason for the abundance of papers written on optimization with $L^{2}$-type objectives MyID:353 . In this work we have used the state of the art on optimal control solvers, provided by PROPT (see our code in Appendix A). It remains an open problem to prove that PROPT candidates for local minimum are indeed (local) minimizers. We are not aware of any kind of solvers that guarantee minimality, either local or global, for our specific problems. The difficulties are related with the occurrence of singular arcs, which may not be located by numerical procedures.

For future work, it would be interesting to validate the obtained numerical results with the necessary optimality results given by Theorems 3.1 and 3.2. Moreover, to prove sufficient optimality conditions and/or global minimality remain also nontrivial open questions.

{acknowledgement}

This work was partially supported by project TOCCATA, reference PTDC/EEI-AUT/2933/2014, funded by Project 3599 – Promover a Produção Científica e Desenvolvimento Tecnológico e a Constituição de Redes Temáticas (3599-PPCDT) and FEDER funds through COMPETE 2020, Programa Operacional Competitividade e Internacionalização (POCI), and by national funds through Fundação para a Ciência e a Tecnologia (FCT) and CIDMA, within project UID/MAT/04106/2013 (Silva and Torres); the post-doc fellowship SFRH/BPD/72061/2010 (Silva); and by the project “Metodi numerici in teoria delle popolazioni” of the Dipartimento di Matematica “Giuseppe Peano” (Venturino). The authors are grateful to two referees for useful comments and suggestions.

Appendix A PROPT Matlab codes

We present here our PROPT Matlab codes so that any interested reader can replicate the numerical results reported in this work. We begin with the code for the optimal control problem investigated in Section 3.1, where the parameter values are given in Table 1.

% Problem setup

toms t tf = 150; p = tomPhase(’p’, t, 0, tf, 101); setPhase(p);

tomStates x1 x2 x3

tomControls u1

x = [x1; x2; x3];

u = [u1];

% Initial conditions

x0i = [3.1; 3.7; 2.2];

x0 = icollocate({x1==x0i(1),x2==x0i(2),x3==x0i(3)});

% Box constraints and boundary

uL = [0]; uU = [1];

cbox = {0 <= icollocate(x1) <= 1000 0 <= icollocate(x2) <= 1000 0 <= icollocate(x3) <= 1000 0 <= icollocate(u1) <= 1 };

cbb = {collocate(uL <= u1 <= uU)

initial(x == x0i)};

% Parameters of the model

a = 3.1; b = 1.2; c= 0.2; e = 2.5; r = 1; W = 5; V = 1000; k = 1; h = 0.7; K = 0.01; q = 0.9;

% ----- control system -------------------

% f = x1; s = x2; v = x3;

ceq = collocate({ dot(x1) == r.x1.(1 - x1/W) - cx2.x1 - h(1-q)u1; dot(x2) == x2.(-a + kbx3 + kcx1) - hKqu1; dot(x3) == ex3.(1 - x3/V) - bx2.x3 - hqu1});

% ----- cost functional ------------------

objectiveJ1 = integrate(x3 + 50*u1.^2);

% objectiveJ1 = integrate(x3);

% ------ Solve the problem --------

options = struct;

options.name = ’system predator prey’;

solution = ezsolve(objectiveJ1, {cbb, cbox, ceq}, x0, options);

For the time-optimal control problem of Section 3.2, the following code was developed:

toms t toms t_f p = tomPhase(’p’, t, 0, t_f, 50); setPhase(p);

tomStates x1 x2 x3

tomControls u

% Box constraints

cbox = {0 <= t_f <= 50 0 <= collocate(u) <= 1};

% Boundary constraints

cbnd = {initial({x1 == 3.1; x2 == 3.7; x3 == 2.2}) final({x3 == 0})};

% Parameters of the model

a = 3.1; b = 1.2; c= 2; e = 0.3; r = 1; W = 5; V = 1000; k = 1; h = 0.7; K = 0.01; q = 0.9;

% ODEs and path constraints

ceq = collocate({dot(x1) == r.x1.(1 - x1/W) - cx2.x1 - h(1-q)u; dot(x2) == x2.(-a + kbx3 + kcx1) - hKqu; dot(x3) == ex3.(1 - x3/V) - bx2.x3 - hqu});

% Objective

objective = t_f;

% Solve the problem

options = struct; options.name = ’minimum time system predator prey’; options.prilev = 1;

solution = ezsolve(objective, {cbox, cbnd, ceq}, options);