Stabilization of structured populations via vector target oriented control

TL;DR

This paper introduces a new control method for stabilizing structured population models, effectively managing chaos and stabilizing periodic solutions in complex discrete systems like the LPA and delayed Ricker models.

Contribution

The paper develops a novel target oriented control approach for structured population models, enabling stabilization of states and periodic solutions in higher-order and delayed systems.

Findings

Successfully stabilizes chaotic dynamics in structured populations.

Extends control techniques to higher-order and delayed difference equations.

Demonstrates effectiveness on LPA and delayed Ricker models.

Abstract

In contrast with unstructured models, structured discrete population models have been able to fit and predict chaotic experimental data. However, most of the chaos control techniques in the literature have been designed and analyzed in a one-dimensional setting. Here, by introducing target oriented control for discrete dynamical systems, we prove the possibility to stabilize a chosen state for a wide range of structured population models. The results are illustrated with introducing a control in the celebrated LPA model describing a flour beetle dynamics. Moreover, we show that the new control allows to stabilize periodic solutions for higher order difference equations, such as the delayed Ricker model, for which previous target oriented methods were not designed.

Equations

$${\mathbf{x}}_{n+1}=\mathbf{f}({\mathbf{x}}_n),n=0,1,2,\dots ,$$

$$x_{n+1}=f(cT+(1-c)x_n),T\ge 0,c\in [0,1).$$

$$\varphi (x)=cT+(1-c)x$$

$$x_{n+1}=cT+(1-c)f(x_n),T\ge 0,c\in [0,1),$$

$${\mathbf{x}}_{n+1}=\mathbf{f}(c\mathbf{T}+(1-c){\mathbf{x}}_n),\mathbf{T}\in D,c\in [0,1),$$

$${\mathbf{x}}_{n+1}=c\mathbf{T}+(1-c)\mathbf{f}({\mathbf{x}}_n),\mathbf{T}\in D,c\in [0,1).$$

$$B=\underset{\mathbf{x}\in S,j=1,\dots ,d}{\max }\sum _{i=1}^d\left|\frac{\partial f_j}{\partial x_i}(\mathbf{x})\right|\text{mbox}andC=\underset{\mathbf{x}\in S,j=1,\dots ,d}{\max }\sum _{i=1}^d\left|\frac{\partial f_i}{\partial x_j}(\mathbf{x})\right|$$

$$c^{\ast }=\max \left\{0,1-\frac{1}{A}\right\}.$$

$$\mathbf{K}=\mathbf{f}(\mathbf{K}).$$

$$\Vert \mathbf{f}(\mathbf{x})-\mathbf{K}\Vert \le L\Vert \mathbf{x}-\mathbf{K}\Vert ,\mathbf{x}\in D.$$

$${\mathbf{x}}_{n+1}=\mathbf{f}((1-c){\mathbf{x}}_n),c\in [0,1),$$

$${\mathbf{x}}_{n+1}=(1-c)\mathbf{f}({\mathbf{x}}_n),c\in [0,1),$$

$$\left\{\begin{array}[]{l}L_{n+1}=bA_{n}\exp(-c_{el}L_{n}-c_{ea}A_{n}),\\ P_{n+1}=(1-\mu_{L})L_{n},\\ A_{n+1}=P_{n}\exp(-c_{pa}A_{n})+A_{n}(1-\mu_{A}),\end{array}\right.$$

$$\mathbf{f}(\mathbf{x})=(b{x}_3\exp (-c_{el}{x}_1-c_{ea}{x}_3),(1-{\mu }_L)x_1,x_2\exp (-c_{pa}{x}_3)+x_3(1-{\mu }_A)).$$

$$\left\{\mathbf{x}\in {\mathbb{R}}_{+}^d:0<\Vert \mathbf{x}\Vert <\max \{\text{mbox}{e}^{-1}b/c_{ea},\Vert \mathbf{T}\Vert \}\right\}.$$

$$J\mathbf{f}(\mathbf{x})=\left(\begin{array}{ccc}-b{c}_{el}{x}_3{\text{mbox}e}^{-c_{ea}{x}_3-c_{el}{x}_1}& 0& b(1-c_{ea}{x}_3){\text{mbox}e}^{-c_{ea}{x}_3-c_{el}{x}_1}\\ 1-{\mu }_L& 0& 0\\ 0& {\text{mbox}e}^{-c_{pa}{x}_3}& 1-{\mu }_A-c_{pa}{x}_2{\text{mbox}e}^{-c_{pa}{x}_3}\end{array}\right),$$

$$u_{n+1}=u_n\exp (r-u_{n-d+1}),$$

$${\mathbf{x}}_{n+1}=\mathbf{f}({\mathbf{x}}_n)$$

$$P$$

$$P=c\Vert \mathbf{f}({\mathbf{x}}^{\ast })-\mathbf{T}\Vert .$$

$${\mathbf{x}}_{n+1}=C\mathbf{T}+(I_d-C)\mathbf{f}({\mathbf{x}}_n),\mathbf{T}\in {\mathbb{R}}^d,C=diag\{c_j\},c_j\in [0,1),$$

$$\mathbf{g}(\mathbf{x})=c\mathbf{T}+(1-c)\mathbf{f}(\mathbf{x}).$$

$$\Vert \mathbf{g}(0)\Vert =\Vert c\mathbf{T}+(1-c)\mathbf{f}(0)\Vert >\Vert 0\Vert .$$

$$\Vert \mathbf{g}(\mathbf{x})\Vert \ge \Vert \mathbf{x}\Vert ,\mathbf{x}\in {\mathbb{R}}_{+}^d,\Vert \mathbf{x}\Vert =m.$$

$$\Vert \mathbf{g}(\mathbf{x})\Vert$$

$$=c\Vert \mathbf{T}\Vert +(1-c)\Vert \mathbf{f}(\mathbf{x})\Vert \le c\Vert \mathbf{x}\Vert +(1-c)\Vert \mathbf{x}\Vert =\Vert \mathbf{x}\Vert .$$

$$\left\{\mathbf{x}\in {\mathbb{R}}_{+}^d:0<\Vert \mathbf{x}\Vert <\max \{M,\Vert \mathbf{T}\Vert \}\right\}$$

$${\mathbf{x}}^{\ast }=c\mathbf{T}+(1-c)\mathbf{f}({\mathbf{x}}^{\ast }).$$

$$R=\{\mathbf{f}(\mathbf{K})+\alpha (\mathbf{K}-\mathbf{f}(\mathbf{K})):\alpha >1\},$$

$$\mathbf{g}(\mathbf{K})=c_{\mathbf{K}}{\mathbf{T}}_{\mathbf{K}}+(1-c_{\mathbf{K}})\mathbf{f}(\mathbf{K})$$

$${\varphi }_2\left({\varphi }_1(\mathbf{f}(\mathbf{x}))\right)$$

Full text

††journal: Bulletin of Mathematical Biology

1 Introduction

In 1970ies, R. May showed in may1974biological ; May1976 that certain one-dimensional discrete population models commonly used in theoretical ecology, can exhibit chaos for some parameter values. This fact immediately started two extremely interesting lines of research. The first line, and probably the hardest, is related to finding chaotic models fitting experimental data and, which is more important, predicting experiment outcomes. May himself, along with Hassel and Lawton, observed that there are species for which, in order to fit experimental data to a simple one-dimensional model, we have to consider parameters that generate chaotic dynamics hassell1976patterns . The task of predicting experimental data using these simple one-dimensional models appears to be hard morris1990problems and, as far as we know, is an open problem. However, more elaborated models, e.g. including age structure or the interaction between different species, were able to fit and predict experimentally the detected chaos Costantino389 ; kot1993population . The second line of research deals with whether and how this complexity can be controlled sole1999controlling . Several strategies have been proposed for controlling chaos in population dynamics, e.g. braverman2012stabilization ; capeans2014less ; Dattani ; desharnais2001chaos ; hilker2005control ; liz2010control ; parthasarathy1995controlling ; sah2013stabilizing ; segura2016adaptive ; tung2014comparison . Most of the control techniques proposed for discrete dynamical systems have been introduced and studied for scalar maps. However, the outcome of the study of experimental data and model-data matching process points out that stabilizing multidimensional maps is more of an interest and practical application.

In the present paper, our purpose is to stabilize nonlinear chaotic dynamical systems given by the first order vector difference equation

[TABLE]

where $\mathbf{f}\colon D\to D$ is a continuous function, $D$ is a convex subset of $\mathbb{R}^{d}$, and $\mathbf{x}_{0}\in D$ is the initial condition. This type of systems is suitable to describe multi-species models, as well as single species populations with a structure (e.g., an age structure given by different age stages, such as juvenile and adults, or a spatial structure given by a population living in different patches connected by dispersal), as well as physical models. System (1) can exhibit chaotic behavior Marotto , and chaos control of multidimensional systems is usually a harder problem than for scalar maps. See, for example, the recent paper LizPotzsche on multidimensional prediction based control, where the challenges to stabilize the two-dimensional Hénon map Henon1976 were outlined.

Here, we consider a natural extension of a method called target oriented control (TOC). This method was introduced in Dattani for first order difference equations, i.e. when $d=1$. TOC establishes a target state and implements an increase or a decrease of the state variable each time step, depending on whether its value exceeds or is below the target state

[TABLE]

In a population, to apply TOC we fix a target population size $T$ and cull/restock a fraction $c$ of the difference between the target and the actual population sizes. This fraction $c$ measures the control intensity. Therefore, TOC is a two-parameter control method in which the controller chooses $T$ and $c$. It is very interesting that recent experiments with fruit flies support that TOC, as other two-parameter control methods, has a better performance than one-parameter control techniques in enhancing simultaneously different ecological stability properties tung2016simultaneous . Indeed, in such experiments TOC concurrently reduced population fluctuations, decreased extinction probability and increased effective population size.

Some of these stabilization properties were explained mathematically in Chaos2014 ; TPC . For instance, it was proved that if the control intensity $c$ is close enough to one, TOC can provide global stabilization of a positive equilibrium for a wide class of smooth maps. In TPC , it was also noticed that if we describe the linear transformation of the variable

[TABLE]

and consider the modified target oriented control (MTOC)

[TABLE]

then the global (local) asymptotic stability of the equilibrium $K_{c}$ of (2) is equivalent to the global (local) asymptotic stability of the equilibrium $P_{c}=\phi(K_{c})$ of (4). In other words, the stability of the equilibrium is not altered by switching the time of control application: before or after the reproduction period.

Here, we consider the natural extensions of TOC and MTOC to multidimensional systems: Vector Target Oriented Control (VTOC)

[TABLE]

as well as Vector Modified Target Oriented Control (VMTOC)

[TABLE]

Our main results give sufficient conditions for the local and global stabilization of either an equilibrium or some other target vector state and estimate the minimum control intensity $c^{\ast}$ to attain such stabilization using these new methods. Regular (5) and modified (6) vector target oriented controls are topologically conjugate (see Lemma 5 in the Appendix), which implies that from a stability perspective both systems have the same properties. Since the results are focused on stability properties, from now on we restrict ourselves to VMTOC without loss of generality.

The paper cushing2002chaos is perhaps the best example of the intersection of the two lines of research described in the first paragraph of this introduction. There, Cushing et al. put forward a control method to stabilize fluctuations of an insect population with three age stages: larvae, pupae and adults. The method was based on the analysis of the chaotic attractor of the theoretical model known to describe the population dynamics. In order to illustrate our results, we consider the same model, called LPA model, showing that, at least theoretically, VMTOC globally stabilizes an equilibrium if the control intensity is high enough. The main advantages of VMTOC, compared to the method presented in cushing2002chaos , are its simplicity, since VMTOC does not need any information about the chaotic attractor, and flexibility, since any age-stage configuration can be stabilized.

We also use the LPA model to illustrate that the selection of the target can have important consequences, not previously reported for target oriented control. We show that depending on the selection of the target vector $\mathbf{T}$, an increase of the control intensity $c$ may not always be stabilizing due to the presence of bubbles.

VMTOC is not the first extension of target oriented control. Indeed, TOC has been recently extended to higher order difference equations braverman2015stabilization . Equations of this type naturally arise when considering populations with non-overlapping generations but with multi-seasonal interactions levin1976note . Since higher order difference equations can be rewritten as a vector map, our results here are also applicable to this particular case. If the target $\mathbf{T}$ is a vector with equal coordinates, the stabilization scheme of braverman2015stabilization can be obtained as a particular case of the results of the present paper. However, stabilizing a certain vector state in $\mathbb{R}^{d}$ with non-equal coordinates corresponds to a $d$-periodic orbit stabilization for the original higher order difference equation. This opens the possibility to stabilize periodic orbits in delayed population models, which seemed not possible with the approach presented in braverman2015stabilization .

Let us note that the method developed in the present paper allows to effectively deal not only with chaotic but also with multi-stable systems, as well as control oscillation amplitudes.

The paper is organized as follows. In Section 2 some auxiliary results are collected: the fact that we can combine two or more VMTOCs to obtain a control of the same type, and that, with an appropriate combination of the target vector $\mathbf{T}$ and the control intensity $c$, we can get any vector in the interior of the domain of $\mathbf{f}$ as an equilibrium (in fact, there is an infinite number of such $c,\mathbf{T}$ pairs). Moreover, we present sufficient conditions for VMTOC to have at least one non-trivial equilibrium for all control intensities $c\in(0,1)$. In Section 3 we justify the possibility of local and global stabilization of an equilibrium of the uncontrolled system or of any prescribed vector in $D$. Section 4 illustrates the results with two different applications: LPA model and the delayed Ricker model. The discussion section summarizes the results obtained and further possible developments. The auxiliary result on the equivalence of VTOC and VMTOC is postponed to the appendix.

2 Calculus of VMTOCs

Since we are interested in the capacity of VMTOC to stabilize an equilibrium, we begin by showing that under quite general conditions on $\mathbf{f}$ and $D$, such an equilibrium exists. Indeed, if $\mathbf{f}:D\to D$ is continuous and the set $D$ is convex and compact, then Brouwer Fixed-point Theorem implies the existence of at least one equilibrium point of VMTOC in $D$ for every $c\in(0,1)$ and any $\mathbf{T}\in D$. However, $D$ is not always compact. For example, in population models, where each component of $\mathbf{x}$ corresponds to a certain population size, it is natural to assume that $D={\mathbb{R}}^{d}_{+}$, ${\mathbb{R}}_{+}=[0,\infty)$ but also that, due to the competition for resources, the inequality $\|\mathbf{f(x)}\|\leq\|\mathbf{x}\|$ holds when $\|\mathbf{x}\|$ is large, where $\|\cdot\|$ is an arbitrary fixed norm in $\mathbb{R}^{d}$. The next result shows that in this situation VMTOC has a nontrivial equilibrium for every $c\in(0,1)$ as well.

Lemma 1

Assume that $\mathbf{f}\colon\mathbb{R}^{d}_{+}\to\mathbb{R}^{d}_{+}$ is continuous and there exists a positive constant $M$ such that $\mathbf{\|f(x)\|\leq\|x\|}$ holds for any $\mathbf{x}\in\mathbb{R}^{d}_{+}$ with $\|\mathbf{x}\|\geq M$. Then for any target $\mathbf{T}\in\mathbb{R}^{d}_{+}\setminus\{\mathbf{0}\}$ and $c\in(0,1)$, there exists at least one equilibrium $\mathbf{x}^{\ast}\in\mathbb{R}^{d}_{+}$ of VMTOC satisfying $0<\|{\mathbf{x}}^{\ast}\|<\max\{M,\|\mathbf{T}\|\}.$ Moreover, if all the components of $\mathbf{T}$ are non-zero, then for all $c\in(0,1)$ the equilibria of VMTOC have all their components positive.

Remark 1

It is interesting to note the practical consequences of Lemma 1. A controller can select the target $\mathbf{T}$ depending on the aimed result. If the system (1) models the interactions among different species, and we aim to eradicate one of them while keeping the others, then it is natural to select $\mathbf{T}$ with a zero component, and such that the fixed point belongs to the boundary of $\mathbb{R}^{d}_{+}$. Whereas if the system (1) models the interactions between different age stages of the same species, and we aim to reduce the fluctuations, then $\mathbf{T}$ will belong to the interior of $\mathbb{R}^{d}_{+}$ and, by Lemma 1, the equilibria of VMTOC as well.

Remark 2

Note that the original map $\mathbf{f}$ in Lemma 1 does not need to have a fixed point in $\mathbb{R}^{d}_{+}\setminus\{\mathbf{0}\}$, for example, the unique fixed point of $\mathbf{f}(\mathbf{x})=\alpha\mathbf{x}$, $\alpha\in(0,1)$, is $\mathbf{0}$. However, with any nontrivial target $\mathbf{T}$, the equilibrium of the controlled map becomes non-trivial.

Next, we consider the possibility to have an arbitrary vector state in the interior of $D$ as an equilibrium in controlled model (6).

Lemma 2

Let $\mathbf{f}\colon D\to D$ be a continuous function, where $D$ is convex. Then for any $\mathbf{K}\in D$, $\mathbf{K}\not\in\partial D$ there exist $c_{\mathbf{K}}\in(0,1)$ and $\mathbf{T_{K}}\in D$ such that $\mathbf{K}$ is an equilibrium of $\mathbf{g}(\mathbf{x})=c_{\mathbf{K}}\mathbf{T_{K}}+(1-c_{\mathbf{K}})\mathbf{f}(\mathbf{x})$.

As easily follows from the proof of Lemma 2, there are infinitely many pairs $(c_{\mathbf{K}},\mathbf{T_{K}})$. For any $c_{\mathbf{K}}\in(0,1)$, $\mathbf{T_{K}}$ is unique, but depending on the geometry of $D$, $\mathbf{T_{K}}$ can be an arbitrary point on the ray starting at $\mathbf{K}$ and $\mathbf{K}-\mathbf{f(K)}$-directed. The closer is $c_{\mathbf{K}}$ to 1, the smaller the distance between $\mathbf{T_{K}}$ and $\mathbf{K}$ is. In other words, we could choose either a larger target and a weaker control, or a smaller target but a tighter control. Of course, this election can affect the stability of the equilibrium $\mathbf{K}$ of $\mathbf{g}$. In the following section we will consider such stability. But before that, we present our last auxiliary result on the effect of combining two different VMTOCs.

Lemma 3

Let $D$ be convex. Then a combination of two VMTOCs is a VMTOC.

3 Stabilization of an equilibrium

3.1 Local stabilization

When a sufficiently strong control is implemented, VMTOC locally asymptotically stabilizes any finite equilibrium.

Theorem 3.1

Assume that $\mathbf{f}=(f_{1},\dots,f_{d})\colon D\to D$ is continuously differentiable and that, for a fixed compact set $S\subset D$, VMTOC with a target $\mathbf{T}\in D$ has at least one equilibrium in $S$ for every $c\in(0,1)$. Define $A=\min\{B,C\}$, where

[TABLE]

Then, all equilibria of VMTOC in $S$ are locally asymptotically stable for $c\in(c^{*},1)$ with

[TABLE]

Remark 3

If $\mathbf{T}$ is chosen as an equilibrium $\mathbf{K}$ of $\mathbf{f}$, then $\mathbf{K}$ is also an equilibrium of VMTOC for every $c\in(0,1)$. In such a case, in Theorem 3.1 we can take $S=\{\mathbf{K}\}$, and in the definition of $c^{*}$ we can replace $A$ by $\rho(J\mathbf{f}(\mathbf{K}))$, where $\rho$ denotes the spectral radius of a matrix.

3.2 Global stabilization

First, we consider the case when an equilibrium of the uncontrolled system is to be stabilized, that is, it satisfies

[TABLE]

In this section, we assume that there exists $L>0$ such that

[TABLE]

Note that (9) yields that $\mathbf{K}$ is a fixed point of $\mathbf{f}$.

Our first result gives not only a sufficient global stabilization condition but also an estimate of the control intensity necessary to achieve it, which depends on the constant $L$ in condition (9).

Theorem 3.2

Assume that for a continuous function $\mathbf{f}\colon D\to D$ equality (8) and inequality (9) hold, where $D\subseteq\mathbb{R}^{d}$ is convex, $\mathbf{K}\in D$ and $\mathbf{T=K}$.

If $L<1$, then all solutions of (6) with an initial condition $\mathbf{x}_{0}\in D$ converge to $\mathbf{K}$ for any $c\in[0,1)$.

If $L\geq 1$, then there exists $c^{*}\in(0,1)$ such that for $c\in(c^{*},1)$ all solutions of (6) with an initial condition $\mathbf{x}_{0}\in D$ converge to $\mathbf{K}$.

If we select the zero target $\mathbf{T}=\mathbf{0}$ in (5), the resulting control will be the proportional feedback (PF)

[TABLE]

assuming the reduction of the state variable, which is proportional to the size of this variable gm . Proportional reduction models constant effort harvesting or culling processes. Switching the variable transformation $\psi(\mathbf{x})=(1-c)\mathbf{x}$ with the map $\mathbf{f}$, we get a modified proportional feedback method (MPF)

[TABLE]

in which the control occurs after the process modeled by $\mathbf{f}$ takes place (e.g. reproduction). Applying Theorem 3.2 to (11), we obtain a result on the stabilization of the origin using MPF control.

Corollary 1

Assume that $D\subseteq\mathbb{R}^{d}$ is convex, $\mathbf{0}\in D$ and that for a continuous function $\mathbf{f}\colon D\to D$, the inequality $\|\mathbf{f(x)}\|\leq L\|\mathbf{x}\|$ holds for $\mathbf{x}\in D$. Then for $c\in(c^{*},1)$ with $c^{*}=\min\{0,1-\frac{1}{L}\}$, the origin is an attractor for any sequence starting with $\mathbf{x}_{0}\in D$ and satisfying the controlled equation (11).

In Theorem 3.2 and Corollary 1, the constant $L$ in (9) is used to estimate the control intensity necessary to stabilize globally an equilibrium: the smaller $L$ is, the sooner is the global stability attained. Therefore, it is interesting to have easy ways to calculate $L$ for a given map. If $\mathbf{f}$ is globally Lipschitz continuous with $\mathbf{K}\in D$ being an equilibrium of $\mathbf{f}$, we can take $L$ as the global Lipschitz constant of $\mathbf{f}$, though this is not necessarily a minimal $L$. It is also possible to estimate $L$ if $\mathbf{f}$ is a locally Lipschitz continuous bounded function.

Lemma 4

Let $\mathbf{f}\colon D\to D$ be a locally Lipschitz continuous bounded function and $\mathbf{K}\in D$ be a fixed point of $\mathbf{f}$. Then there exists $L\geq 1$ such that condition (9) holds for any $\mathbf{x}\in D$.

So far we have justified the possibility to stabilize only a fixed point of the original map with VMTOC. Let us aim to stabilize an arbitrary interior point in $D$ with VMTOC.

Theorem 3.3

Suppose $\mathbf{f}\colon D\to D$ is a continuous function which is either globally Lipschitz, or locally Lipschitz and bounded, and $D$ is convex. Then for any $\mathbf{K}\in D$, $\mathbf{K}\not\in\partial D$ there exists a VMTOC for which $\mathbf{K}$ is a globally asymptotically stable equilibrium.

4 Applications

4.1 LPA model

Consider the following age structured model designated to describe the changes in the densities of the flour beetle Tribolium castaneum life stages

[TABLE]

where $L_{n}$ is the number of feeding larvae, $P_{n}$ is the number of nonfeeding larvae, pupae and callow adults, and $A_{n}$ is the number of sexually mature adults at stage $n$, whereas $b>0$ is the number of larval recruits per adult per unit of time in the absence of cannibalism, $c_{el},c_{ea},c_{pa}\in\mathbb{R}_{+}$ account for the cannibalism of eggs by both larvae and adults and the cannibalism of pupae by adults, and $\mu_{L},\mu_{A}\in(0,1)$ are the larval and adult rates of mortality. We refer to desharnais2001chaos for a more detailed explanation of the model.

System (12) is known as the LPA (larvae-pupae-adults) model. Depending on the values of the parameters, the dynamics predicted by the LPA model can be different (e.g. stable equilibria, periodic cycles, chaotic oscillations) and, very interestingly, these different types of dynamics have been demonstrated with flour beetle populations in the laboratory.

Here, we set $b=10.45$, $c_{el}=0.01731$, $c_{ea}=0.01310$, $c_{pa}=0.35$, $\mu_{L}=0.200$, and $\mu_{A}=0.96$. For these parameter values the LPA model exhibits chaotic oscillations and there exists a strange attractor desharnais2001chaos . Figure 1 shows this chaotic attractor and the oscillations for all age stages.

In order to stabilize this chaotic LPA model, Desharnais et al proposed the in-box control method desharnais2001chaos . This method has two steps: first, detecting the region of the attractor more sensible to small perturbations by computing the local Lyapunov exponents bailey1996local , and, second, modifying the population (by adding a fixed number of individuals) when the population is inside that region. The in-box method was able to stabilize chaotic populations of the flour beetle in the laboratory, and the experimental data obtained was predicted correctly by numerical simulations. The in-box method is an empirically tested control of chaos strategy in age-structured population dynamics, see also fryxell2005evaluation ; sah2013stabilizing ; tung2016simultaneous . However, the in-box method is not easy to apply. A controller needs to detect the region of the attractor more sensible to small perturbations and then to establish an action (remove or add certain type of individuals) that sends the population out of that region. Here, we show that, at least theoretically, VMTOC can be used to stabilize age-structured populations in a simpler way.

LPA model (12) can be written as $\mathbf{x}_{n+1}=\mathbf{f}(\mathbf{x}_{n})$ with $\mathbf{x}=(x_{1},x_{2},x_{3})$

[TABLE]

Next, note that $f_{1}(\mathbf{x})\leq\mbox{e}^{-1}b/c_{ea}\approx 295,44$, $f_{2}(\mathbf{x})\leq 0.8x_{1}$, and $f_{3}(\mathbf{x})\leq\max\{x_{2},x_{3}\}$ for each $\mathbf{x}\in\mathbb{R}^{3}_{+}$. Fixing the $\max$-norm $\|\mathbf{x}\|=\max_{1\leq j\leq n}|x_{j}|$, we have that $\mathbf{\|f(x)\|\leq\|x\|}$ for any $\mathbf{x}$ with $\|\mathbf{x}\|\geq\mbox{e}^{-1}b/c_{ea}$. Lemma 1 guarantees that applying VMTOC to the LPA model, independently of the target $\mathbf{T}\in\mathbb{R}^{3}_{+}\setminus\{\mathbf{0}\}$ and the control intensity $c\in(0,1)$, there exists at least one equilibrium in the open domain

[TABLE]

Numerically, one finds that $\mathbf{K}\approx(28.0120,22.4096,4.6251)$ is a fixed point of the LPA model. Applying VMTOC with target $\mathbf{T}=(28.0120,22.4096,4.6251)$ to the LPA model 12 has the effect illustrated in Figure 2.

This agrees with Theorem 3.1 and Remark 3. Indeed, the Jacobian of $\mathbf{f}$ is

[TABLE]

and calculating numerically the value of the spectral radius of $J\mathbf{f}(\mathbf{T})$, we obtain $\rho(J\mathbf{f}(\mathbf{T}))\approx 1.3803$. Therefore, using Theorem 3.1 and Remark 3, we can guarantee that $\mathbf{K}$ is asymptotically stable for MVTOC if $c$ is greater than $c^{*}=1-\frac{1}{1.3803}\approx 0.2756$.

For the scalar case, in TPC the effect of changing the size of the target on TOC and MTOC was studied. There, the focus was mainly on the size of the equilibrium, showing, for example, that this size increases as the size of the target increases. Here, we present an effect of changing the target which has not been previously reported. To this end, we select three different targets $\mathbf{T}_{1}=(30,30,200)$, $\mathbf{T}_{2}=(30,200,30)$, and $\mathbf{T}_{3}=(200,30,30)$ and use the same Euclidean norm. Note that choosing one of them means that the controller designs one of the age stages to be prevalent in comparison to the other two.

Our results show that independently of the target, VMTOC will stabilize an equilibrium if the control intensity is large enough. However, the responses are different, not only from the point of view of the sizes of the different age stages at the equilibrium (which one should expect), but also the effect of increasing $c$ on the stability of the equilibrium. For targets $\mathbf{T}_{2}$ and $\mathbf{T}_{3}$ (see Figure 3), increasing $c$ has no negative effect on the stability: if for $c=\hat{c}$ the controlled system has a stable equilibrium, then it has a stable equilibrium for any $c\in(\hat{c},1)$. For the other target, $\mathbf{T}_{1}$, this is not true: for a certain $c=\hat{c}$, the controlled system has a stable equilibrium for $c\in(\hat{c},c_{1})$, while for $c\in(c_{1},c_{2})$, the adult population is still stable but both pupae and larvae have a stable two-cycle, which has the form of bubbling, see Figure 4. Increasing further the control intensity further removes these oscillations. Figure 4 illustrates it, for $c=0.1$ the three age stages tend to an equilibrium value, but increasing the control intensity to any $c$ in approximately $(0.17,0.36)$ causes oscillations in larvae and pupae populations; if $c>0.36$, the three age stages tend again to their equilibrium values.

4.2 Higher order equations revisited

In braverman2015stabilization , we showed that using a fixed target is possible to stabilize an equilibrium of a chaotic higher order difference equation. Higher order difference equations arise, for instance, when studying multi-seasonal interactions in a population with non-overlapping generations levin1976note .

Although increasing the control intensity in the method discussed in braverman2015stabilization stabilizes an equilibrium, it does not follow the characteristic route from chaos of folding-period bifurcations that many stabilizing strategies show. Indeed, TOC itself presents this route for one-dimensional models—see the bifurcation diagram of Figure 1 in Dattani and compare with the bifurcation diagrams in Figures 1, 2 and 3 in braverman2015stabilization for higher-dimension models. Those bifurcation diagrams indicate that the stabilization of an equilibrium in higher order equations happens through a Neimark-Sacker bifurcation while in the one-dimensional case happens through a period-folding bifurcation. Thus, if the aim is stabilization of periodic orbits of higher order equations, direct application of TOC presented in braverman2015stabilization is not useful.

Let us illustrate with a numerical example that applying a periodic target, i.e. using VMTOC with a target $\mathbf{T}$ with non-equal components, can stabilize a periodic orbit for higher-order equations. We consider the delay Ricker equation in the form of peran2015global

[TABLE]

where $d\geq 2$ is a fixed natural number determining the time lag in the intraspecific regulatory mechanisms of the population. Equation (13) can be rewritten as the system

[TABLE]

with $\mathbf{x}=(x_{1},x_{2},\dots,x_{k})$ and $\mathbf{f}(\mathbf{x})=(x_{1}\exp(r-x_{k}),x_{1},\dots,x_{k-1})$.

Let us fix $r=2$, $d=2$ and choose the target $\textbf{T}=(1,3)$. Figure 5B shows the effect of applying VMTOC to system (14). There, the population evolves without any control during the first 20 generations, then VMTOC with control intensity $c=0.4$ is applied. We can observe how a period-two orbit is stabilized. Figure 5A corresponds to the population dynamics without control, showing that the uncontrolled population has phases of low population density followed by phases with higher-density.

5 Discussion

Compared to other (especially one-parameter) methods, target oriented control (TOC) has the advantage that it allows independent choice of both the stabilized state and the control intensity. For instance, TOC can stabilize a point which is not an equilibrium of the original map. The vector modifications VTOC and VMTOC inherit these two properties of scalar TOC: a variety of originally unstable states can be stabilized. Certainly the minimal control sufficient for stabilization would depend on the choice of the target $\mathbf{T}$ and the state to be stabilized. However, as we can observe, the ability to stabilize an arbitrary state can be quite cost-involved. The cost efficiency of TOC is discussed in Dattani and analyzed numerically. If we describe the eventual cost-per-step $P$ as the “average” control perturbation

[TABLE]

then, once a state $\mathbf{x}^{\ast}$ is stabilized, the eventual cost-per-step becomes

[TABLE]

If the state to be stabilized is a fixed point $\mathbf{K}=\mathbf{f}(\mathbf{K})$ of the original vector map, the cost-per-step may be quite high in the transient period but eventually tends to zero. From the above, the closer we can choose the target to the image of the stabilized state, more efficiency can be achieved. However, the minimal stabilizing $c$ also depends on the choice of $\mathbf{x}^{\ast}$ and $\mathbf{T}$. The cost optimization of VMTOC is a separate question, and its solution is not in the framework of the current paper. Also, simplified estimate (15) does not take into account the following factors:

• •

the costs of culling and restocking can be different;

• •

the costs may strongly depend on either age stage or patch location.

Moreover, culling and restocking of certain age-stages or locations can be problematic. This brings us to the discussion of some further VMTOC generalization. Originally, TOC included two parameters: the target $T$ and the control intensity $c$. We considered VMTOC in $\mathbb{R}^{d}$ with $d+1$ parameters involved. This allowed us to obtain results on the possibility of stabilization for a variety of states, and to estimate the minimal control intensity $c$ in each case. However, it is quite natural also to consider a $2d$-parameter method with a diagonal control matrix $C$, where the controls $c_{j}=c_{jj}$ are applied to each $j$-th stage. With this modification, VMTOC has the form

[TABLE]

where $I_{d}$ is the $d\times d$ identity matrix and $j=1,2,\dots,d$. If some stages cannot be controlled, the corresponding $c_{j}$ will be zero.

A higher order difference equation is a particular case of the one-dimensional vector equation, and VMTOC allows to stabilize periodic orbits of higher order equations, see the example of stabilization for a two-orbit of the delayed Ricker model above. Let us note that the same technique applies to stabilizing a periodic $k$-orbit of a vector map. Indeed, denote by $\mathbf{f}^{k}$ the $k$-th iterate of $\mathbf{f}$. Then stabilizing a periodic $k$-orbit of $\mathbf{f}$ is equivalent to stabilizing a certain state of $\mathbf{f}^{k}$. All the results of the present paper apply to this case. We can also consider a control type similar to (16), with either $c_{j}=c$ or zeros on the main diagonal of $C$. In the case of $c_{1}=\dots=c_{d}=c$, all the other $c_{j}$ for stabilization of $\mathbf{f}^{k}$ being zeros, the result is the pulse control Chaos2014 . If all $c_{j}=0$, except $c_{k},c_{k+d},\dots$, this corresponds to a control of only one age stage or patch. While it may definitely be problematic to achieve stabilization goals with this limited type of control, it is still an interesting question whether controlling one stage only (for example, juveniles), we can reduce the risk of extinction and population fluctuations.

Appendix A Appendix

The next result shows that VTOC and MVTOC are topologically conjugate, thus they have the same dynamics. We recall that two maps $\phi$ and $\psi$ are topologically conjugate if there is a homeomorphism $h$ such that $\phi\circ h=h\circ\psi$.

Lemma 5

Assume $\mathbf{f}\colon D\to D$, with $D$ convex and $\mathbf{T}\in D$. Then the difference equations (5) and (6) are topologically conjugate.

Proof

We are going to show that the maps defining the difference equations (5) and (6) are topologically conjugate. We begin by defining such maps.

Consider the map $\varphi(\mathbf{x})=c\mathbf{T}+(1-c)\mathbf{x}$ from $D$ to $\varphi(D)$. Moreover, since $D$ is convex $\varphi(D)\subset D$ and the map $\mathbf{f}\circ\varphi\colon D\to D$ is well defined. Clearly, map $\mathbf{f}\circ\varphi$ defines the recurrent relation in equation (5). On the other hand, note that after the first iterate the solutions of (6) belong to $\varphi(D)$. Therefore, we have that after the first iterate the map that defines the recurrence given by (6) is $\varphi\circ\mathbf{f}\colon\varphi(D)\to\varphi(D)$.

It is easy to check that $\varphi$ is a homeomorphism from $D$ to $\varphi(D)$ and obviously $(\varphi\circ\mathbf{f})\circ\varphi=\varphi\circ(\mathbf{f}\circ\varphi)$. Hence, $\mathbf{f}\circ\varphi$ and $\varphi\circ\mathbf{f}$ are topologically conjugate.

Proof

of Lemma 1 An equilibrium $\mathbf{x}^{\ast}$ of VMTOC is a fixed point of the map

[TABLE]

Since $\mathbf{T}\in\mathbb{R}^{d}_{+}\setminus\{\mathbf{0}\}$ and $\mathbf{f}\colon\mathbb{R}^{d}_{+}\to\mathbb{R}^{d}_{+}$, we have for every $c\in(0,1)$ that $c\mathbf{T}+(1-c)\mathbf{f(0)}\in\mathbb{R}^{d}_{+}\setminus\{\mathbf{0}\}$ and

[TABLE]

Hence, by the continuity of $\mathbf{g}$ and the norm, for each fixed $c\in(0,1)$ is possible to find $0<m<M$ such that

[TABLE]

On the other hand, we have for any $\mathbf{x}\in\mathbb{R}^{d}_{+}$ with $\|\mathbf{x}\|\geq\max\{M,\|\mathbf{T}\|\}$,

[TABLE]

Thus, by Krassnosel’skiĭ Fixed-point Theorem for cone-compressing operators (see e.g. (granas2013fixed, , Theorem 7.12)), the map $\mathbf{g}$ has at least one fixed point ${\mathbf{x}}^{\ast}$ in the set

[TABLE]

for every $c\in(0,1)$.

Finally, the last statement of the lemma follows from noticing that an equilibrium $\mathbf{x}^{*}$ of VMTOC satisfies

[TABLE]

∎

Proof

of Lemma 2 If $\mathbf{K}=\mathbf{f}(\mathbf{K})$, we take $\mathbf{T}_{\mathbf{K}}=\mathbf{K}$ and any $c_{\mathbf{K}}\in(0,1)$. Let $\mathbf{K}\neq\mathbf{f}(\mathbf{K})$. Consider the set of points

[TABLE]

which is the ray in the direction $\mathbf{K}-\mathbf{f}(\mathbf{K})$ starting at $\mathbf{K}$. Using that $\mathbf{K}\not\in\partial D$, we have that $R\cap D\neq\emptyset$. Therefore, there exists $\alpha>1$ such that $\mathbf{T_{K}}:=\alpha\mathbf{K}+(1-\alpha)\mathbf{f(K)}\in D$. Fixing $c_{\mathbf{K}}=\alpha^{-1}\in(0,1)$ and noticing that

[TABLE]

concludes the proof.∎

Proof

of Lemma 3 We note that it is sufficient to consider only controls $c\in(0,1)$, since for $c=0$ we have the original map. Let $\phi_{1}(\mathbf{x})=c_{1}\mathbf{T}_{1}+(1-c_{1})\mathbf{x}$ be applied first, and next another argument transformation $\phi_{2}(\mathbf{x})=c_{2}\mathbf{T}_{2}+(1-c_{2})\mathbf{x}$. Then

[TABLE]

where

[TABLE]

Since $c_{1},c_{2}\in(0,1)$, also $c=c_{1}(1-c_{2})+c_{2}\in(0,1)$. On the other hand, the vector $\mathbf{T}=\alpha\mathbf{T}_{1}+(1-\alpha)\mathbf{T}_{2}$, where $\alpha=c_{1}(1-c_{2})/c\in(0,1)$, therefore by the convexity of $D$ the target $\mathbf{T}\in D$.∎

Proof

of Theorem 1 Recall that, by linearization, any equilibrium $\mathbf{p}_{c}$ of VMTOC is asymptotically stable if the spectral radius of the Jacobian matrix $J\mathbf{g}$ of $\mathbf{g(x)}=c\mathbf{T}+(1-c)\mathbf{f(x)}$ at $\mathbf{p}_{c}$ is smaller than $1$.

By a well-known bound for the eigenvalues of a matrix (see e.g. (horn2012matrix, , Corollary 6.1.5)), we have that any eigenvalue $\lambda$ of $J\mathbf{g}(\mathbf{p}_{c})$ satisfies

[TABLE]

Therefore, for $c\in(c^{*},1)$ the eigenvalues of $J\mathbf{g}(\mathbf{p}_{c})$ have modulus smaller than $1$. ∎

Proof

of Theorem 2 If $L$ in (9) satisfies $L\in(0,1)$ and $\mathbf{x}_{n}$ is a solution of (6) with an initial condition $\mathbf{x}_{0}\in D$ and $c\in[0,1)$, then for $n\in\mathbb{N}$

[TABLE]

For any $\varepsilon\in(0,\|\mathbf{x}_{0}-\mathbf{K}\|)$, we have

[TABLE]

so $\mathbf{x}_{n}\to\mathbf{K}$ as $n\to\infty$.

Next, let $L\geq 1$. Denote

[TABLE]

and assume that $c\in(c^{*},1)\subseteq(0,1)$. Then $(1-c)L<1$, denote $\theta=(1-c)L\in(0,1)$. We have $1-c=\theta/L$, $c\mathbf{T}=c\mathbf{K}$ and

[TABLE]

By induction,

[TABLE]

Therefore, $\|\mathbf{x}_{n+1}-\mathbf{K}\|<\varepsilon<\|\mathbf{x}_{0}-\mathbf{K}\|$ for any $\displaystyle n>\left.\ln\left(\frac{\varepsilon}{\|\mathbf{x}_{0}-\mathbf{K}\|}\right)\right/\ln\theta$. Thus, $\lim\limits_{n\to\infty}\mathbf{x}_{n}=\mathbf{K}$. Moreover, $\|\mathbf{x}_{n}-\mathbf{K}\|$ decays at least geometrically, which concludes the proof. ∎

Proof

of Lemma 4 Let $M>0$ be an upper bound of $\mathbf{f}$:

[TABLE]

Since $\mathbf{f}$ is locally Lipschitz continuous, for $\mathbf{x}$ in the intersection of the ball $\|\mathbf{x-K}\|\leq\|\mathbf{K}\|$ with $D$, there is a constant $\tilde{L}>0$ such that $\|\mathbf{f(x)-K}\|\leq\tilde{L}\|\mathbf{x-K}\|$. Next, let $\mathbf{x}\in D$, $\|\mathbf{x-K}\|>\|\mathbf{K}\|$. Then, by (19),

[TABLE]

Finally, choosing

[TABLE]

we obtain that inequality (9) is satisfied for any $\mathbf{x}\in D$. ∎

Proof

of Theorem 3 By Lemma 2 there exist $c_{\mathbf{K}}\in(0,1)$ and $\mathbf{T_{K}}\in D$ such that $\mathbf{K}$ is an equilibrium of $\mathbf{g}(\mathbf{x})=c_{\mathbf{K}}\mathbf{T_{K}}+(1-c_{\mathbf{K}})\mathbf{g}(\mathbf{x})$. By Lemma 4, $\mathbf{f}$ satisfies (9) with some constant $L_{1}$ instead of $L$. Then

[TABLE]

thus condition (9) holds for $\mathbf{g}$ with $L=(1-c_{\mathbf{K}})L_{1}$.

Now, applying Theorem 3.2, we obtain that there exists $c^{*}\in[0,1)$ such that for $\hat{c}\in(c^{*},1)$ and $\mathbf{T=K}$, all solutions of (6) with $\mathbf{g}$ instead of $\mathbf{f}$, $\hat{c}$ instead of $c$ and $\mathbf{x}_{0}\in D$ converge to $\mathbf{K}$. By Lemma 3, a combination of two VMTOCs is a VMTOC. Thus, if we choose $c=c_{\mathbf{K}}(1-\hat{c})+\hat{c}$, where $\hat{c}\in(c^{*},1)$, $c_{\mathbf{K}}$ is defined above, and

[TABLE]

we get that $\mathbf{K}$ is a global attractor of the combination of VMTOCs.∎

Acknowledgements.

The authors are grateful to the anonymous referees whose valuable comments contributed to the presentation of the results of the paper.