Controlling Complex Dynamics in a Protected Area Discrete-Time Model

Abstract

This paper investigates how the introduction of user fees and defensive expenditures changes the complex dynamics of a discrete-time model, which represents the interaction between visitors and environmental quality in an open-access protected-area (OAPA). To investigate this issue more deeply, we begin by studying in great detail the OAPA model, and then we introduce the user fee (𝛽) and the defensive expenditures (𝜌𝛽) specifically directed towards at the protection of the environmental resource. We observed that some values of 𝛽 can generate a chaotic regime from a stable dynamic of the OAPA model. Finally, to eliminate the chaotic regime, we design a controller by OGY method, assuming the user fee as a controller parameter.

1. Introduction

Empirical analysis has shown that tourists are willing to pay more for environmental management, if they believe that the money they pay will be allocated for biodiversity conservation and protected area management (see [1, 2]). Consequently, the funds for maintaining public goods can be increased by fees payed by visitors of the protected areas (PAs).

Several works in economic literature analyze the effects of ecological dynamics generated by economic activity and environmental defensive expenditures. In particular, [3, 4] analyze the stabilizing effect of ecological equilibria in an optimal control context in which ecological dynamics are represented by predator-prey equations.

More recently, economists, social and political scientists have started to develop and adapt chaos theory as a way of understanding human systems. Specifically, [5–10] have considered chaos theory as a way of understanding the complexity of phenomena associated with tourism.

In [11] a three-dimensional environmental defensive expenditures model with delay is considered. The model is based on the interactions among visitors 𝑉, quality of ecosystem goods 𝐸, and capital 𝐾, intended as accommodation and entertainment facilities, in PA. The visitors’ fees are used partly as a defensive expenditure and partly to increase the capital stock.

Based on the continuous environmental model of [11], in this paper we analyze a discrete-time model with no capital stock and with no time delay. We aim at analyzing how the dynamics change when switching from an open-access protected area (OAPA) regime, where there are no services or facilities, to a PA regime with visitor fees used for environmental protection.

This paper is organized as follows. In Section 2, we present the discrete-time model that embodies the user fees and defensive expenditures. In Section 3, the dynamics of an open-access protected area, that is, without the user fee and defensive expenditures, is studied, including stable fixed points, periodic motions, bifurcations (flip-flop and Neimark-Sacker bifurcations), and chaos. Section 4 deals with the control of chaotic motion and the process of control is achieved an appropriate determination of user fees and defensive expenditures.

2. The Mathematical Model with User Fee and Defensive Expenditures

The model refers to a generic protected area and describes the interplay between two state variables: the size 𝑉(𝑡) of the population of visitors of the protected area at time 𝑡 and an index 𝐸(𝑡) measuring the quality of environmental resources of the protected area. The dynamic of 𝑉(𝑡) is assumed to be described by the differential equation:𝑑𝑉(𝑡)(𝑡)𝑑𝑡=−𝑏−𝑐𝑉+𝑑𝐸. (2.1) According to such equation, the time evolution of 𝑉(𝑡) depends on three factors: (i) −𝑏 represents the negative effect of the fee (𝑏) that visitors have to pay to enter the protected area; (ii) −𝑐𝑉 is the negative effect due to congestion; (iii) 𝑑𝐸 (𝑑 is the parameter that presents attractiveness associated with high environmental quality) is the positive effect of environmental quality on visitors’ dynamics. 𝑏, 𝑐, and 𝑑 are strictly positive parameters.

The dynamic of the environmental quality index 𝐸(𝑡) is assumed to be given by:𝑑𝐸𝑑𝑡=𝑟₀(1−𝐸(𝑡))𝐸(𝑡)−𝑎𝑉²(𝑡)+𝑞𝑏𝑉(𝑡), (2.2) where the time evolution of environmental quality is described by a logistic equation (see [12]). According to (2.2), visitors generate a negative impact on environmental quality (this effect is represented by −𝑎𝑉²); however, visitors also generate a positive effect, in that, a share 𝑞 of the revenues 𝑏𝑉 deriving from the fees is used for environmental protection (this effect is represented by 𝑞𝑏𝑉). 𝑟₀ and 𝑎 are strictly positive parameters, while 𝑞 is a parameter satisfying 0≤𝑞≤1.

Euler’s difference scheme for the continuous system (2.1)-(2.2) takes the form (see [13]):𝑉(𝑡+Δ𝑡)−𝑉(𝑡)(𝑡)(𝑡),𝐸()(𝑡)Δ𝑡=−𝑏−𝑐𝑉+𝑑𝐸𝑡+Δ𝑡−𝐸Δ𝑡=𝑟₀(1−𝐸(𝑡))𝐸(𝑡)−𝑎𝑉²(𝑡)+𝑏𝑞𝑉(𝑡), (2.3) where Δ𝑡 denotes the time step. As Δ𝑡→0, the discrete system converges to the continuous system. Roughly speaking, a discrete system can give rise to the same dynamics as a continuous system if the Δ𝑡 is small enough. However, it may generate qualitatively different dynamics if Δ𝑡 is large. In this sense, the discrete system with Δ𝑡>0 generalizes the corresponding continuous system. In the following, we first simplify the discretized system (2.3). Notice that a variable 𝑤(𝑡) in continuous time can be written as 𝑤(𝑡_𝑛) in discrete time. Set 𝑡_𝑛=Δ𝑡⋅𝑛(𝑛=1,2,…); then, given Δ𝑡>0, the variable can be expressed as follows;𝑤(𝑡_𝑛)=𝑤(Δ𝑡⋅𝑛)=𝑤_𝑛and𝑤(𝑡_𝑛+Δ𝑡)=𝑤(Δ𝑡⋅(𝑛+1))=𝑤_𝑛+1. Thus the discretized dynamic system (2.3) can be written as: 𝑥_𝑛+1=𝑥_𝑛−𝑏Δ𝑡−𝑐Δ𝑡𝑥_𝑛+𝑑Δ𝑡𝑦_𝑛,𝑦_𝑛+1=𝑦_𝑛+𝑟₀Δ𝑡1−𝑦_𝑛𝑦_𝑛−𝑎Δ𝑡𝑥²_𝑛+𝑏Δ𝑡𝑞𝑥_𝑛. (2.4)

The length of each period is equal to Δ𝑡. For notational convenience, replacing 𝑛 with 𝑡 and posing 𝑟=𝑟₀Δ𝑡, 𝛼=𝑎Δ𝑡, 𝛽=𝑏Δ𝑡, 𝛾=𝑐Δ𝑡, 𝜌=𝑞Δ𝑡, and 𝜎=𝑑Δ𝑡,we obtain the following discrete-time system:𝑥_𝑡+1=𝑥_𝑡−𝛽−𝛾𝑥_𝑡+𝜎𝑦_𝑡,𝑦_𝑡+1=𝑦_𝑡+𝑟1−𝑦_𝑡𝑦_𝑡−𝛼𝑥²_𝑡+𝛽𝜌𝑥_𝑡, (2.5) where 𝑥_𝑡 and 𝑦_𝑡 represent, respectively, the size 𝑉 of the population of visitors and the value of the quality index 𝐸 at time 𝑡; the parameters 𝛼, 𝛽, 𝛾, 𝜎, 𝜌, and 𝑟 have the same meaning of the corresponding parameters 𝑎, 𝑏, 𝑐, 𝑑, 𝑞, and 𝑟₀, in the system (2.3).

3. The Dynamic Behavior of an Open-Access PA Model

In this section, we analyze the dynamics of our model under the assumption of free-access in the protected area; in this context, visitors do not have to pay a fee to visit the area, and system (2.5) becomes𝑥_𝑡+1=𝑥_𝑡−𝛾𝑥_𝑡+𝜎𝑦_𝑡,𝑦_𝑡+1=𝑦_𝑡+𝑟1−𝑦_𝑡𝑦_𝑡−𝛼𝑥²_𝑡. (3.1)

To compute the fixed points of (3.1), we have to solve the nonlinear system of equations:𝑥=𝑥−𝛾𝑥+𝜎𝑦,𝑦=𝑦+𝑟1−𝑦𝑦−𝛼𝑥². (3.2)

Proposition 3.1. The system (3.1) always present two fixed points:(a)𝑃₁=(𝑥^∗₁,𝑦^∗₁)=(0,0); (b)𝑃₂=(𝑥^∗₂,𝑦^∗₂)=(𝑟𝛾/𝜎/(𝛼+𝑟(𝛾/𝜎)²),(𝛾/𝜎)𝑥^∗).

Now we study the stability of these fixed points. The local stability of a fixed point (𝑥^∗,𝑦^∗) (it denotes (𝑥^∗₁,𝑦^∗₁) or (𝑥^∗₂,𝑦^∗₂)) is determined by the modules of the eigenvalues of the characteristic equation evaluated at the fixed point.

The Jacobian matrix of the system (3.1) evaluated at (𝑥^∗,𝑦^∗) is given by𝐽=−𝛾+1𝜎−𝜃₁1+𝜃₂, (3.3) where 𝜃₁(𝜎)=2𝛼𝑥^∗ and 𝜃₂(𝜎)=𝑟(1−2𝑦^∗). The characteristic equation of the Jacobian matrix 𝐽 can be written as𝜆²+𝑝(𝜎)𝜆+𝑞(𝜎)=0, (3.4) where 𝑝(𝜎)=𝛾−𝜃₂(𝜎)−2 and 𝑞(𝜎)=(1+𝜃₂(𝜎))(1−𝛾)+𝜎𝜃₁(𝜎). In order to study the moduli of the eigenvalues of the characteristic equation (3.4), we first give the following lemma, which can be easily proved.

Lemma 3.2. Let 𝐹(𝜆)=𝜆²+𝑝𝜆+𝑞. Suppose that 𝐹(1)>0,𝜆₁ and𝜆₂are two roots of𝐹(𝜆)=0. Then: (i)|𝜆₁|<1 and |𝜆₂|<1(sink) if and only if 𝐹(−1)>0 and 𝑞<1;(ii)|𝜆₁|<1 and |𝜆₂|>1(or |𝜆₁|>1 and (|𝜆₂|<1)(saddle) if and only if 𝐹(−1)<0;(iii)|𝜆₁|>1 and |𝜆₂|>1(source) if and only if 𝐹(−1)>0 and 𝑞>1;(iv)𝜆₁=−1 and |𝜆₂|≠1(flip-flop bifurcation) if and only if 𝐹(−1)=0 and 𝑝≠0,2;(v)𝜆₁ and 𝜆₂ are complex and |𝜆₁|=|𝜆₂|=1(Neimark-Sacker bifurcation) if and only if 𝑝²−4𝑞<0 and 𝑞=1.

From Lemma 3.2, it follows the following.

Proposition 3.3. The fixed point 𝑃₁=(0,0) is always unstable, while the fixed point 𝑃₂, varying 𝜎, can be a sink, a source, or a saddle (see Figure 1).

Figure 1 

The parameter values are 𝛼=0.12, 𝛽=0, 𝛾=0.375, 𝜌=0, and 𝑟=2.8.

Figure 1 shows the values of 𝐹(−1), 𝑞−1, 𝑝²−𝑞, defined in Lemma 3.2, as functions of the parameter 𝜎.

We fix 𝛼=0.12, 𝛾=0.375, and 𝑟=02.8, and assume that 𝜎 can vary. Smaller values of 𝜎 (see Figure 1) give rise to real eigenvalues, while higher values of it give rise to complex eigenvalues.

According to Lemma 3.2, when the parameter 𝜎 belongs to the interval (0,𝜎_ff) (dash-dot line) we are in the situation described in point (ii) of Lemma 3.2; when 𝜎=𝜎_ff=0.656407, a flip-flop bifurcation occurs; when 𝜎∈(𝜎_ff,𝜎_NS), we are in the context described in (i) (solid line); at the value 𝜎_NS=1.416516, a Neimark-Saker bifurcation takes place; finally, for 𝜎>𝜎_NS, the fixed point becomes unstable.

Such results are illustrated in Figures 2 and 3, which show that some remarkable phenomena occur.

(a)

(b)

Figure 2 

Bifurcation diagrams for the state variable 𝑥 (a) and for the state variable 𝑦 (b), varying 𝜎. The parameter values are 𝛼=0.12, 𝛽=0, 𝛾=0.375, 𝜌=0, and 𝑟=2.8.

(a) 𝜎=0.165

(b) 𝜎=1.415

(c) 𝜎=1.4165

(d) 𝜎=1.42

(e) 𝜎=1.477

Figure 3 

Phase plot with the parameters of Figure 1. (a) chaotic trajectory; (b) the stable fixed point before the Neimark-Sacker bifurcation occurs; (c) the Neimark-Sacker bifurcation; (d) the stable invariant closed curve around the fixed point created after bifurcation; (e) chaotic trajectory.

Figure 3(a) shows a strange attractor appearance posing 𝜎=0.165. If the value of 𝜎 increases, we obtain the attractive fixed point showed in Figure 3(b): both variables of the dynamic system approach a unique fixed point independently from the initial state. The fixed point is characterized by the coordinates 𝑥^∗₂=2.223 and 𝑦^∗₂=0.6949. The eigenvalues of the Jacobian matrix evaluated at such point are 𝜆=0.226651±𝑖0.7155 with |𝜆|=0.7635.

Continuing to increase the value of 𝜎, a Neimark-Sacker bifurcation takes place. For the parameter value 𝜎=1.4165, the fixed point has coordinates 𝑥^∗₂=2.3432 and 𝑦^∗₂=0.6205, and the associated pair of complex conjugate eigenvalues are 𝜆=0.47498±𝑖0.8799 with |𝜆|=1.000; this shows that the eigenvalues belong to the unit circle, and the stability properties of the equilibrium change through a Neimark-Sacker bifurcation. Figure 3(c) illustrates the phase plot corresponding to the bifurcation value of 𝜎.

Continuing to increase the value of 𝜎, we see what happens for 𝜎=1.42. The coordinates of the fixed point are 𝑥^∗₂=2.3456 and 𝑦^∗₂=0.61937, and the associated eigenvalues are 𝜆=0.4782±𝑖0.8819. The modulus of the complex conjugate eigenvalues is |𝜆|=1.0032, and so we can conclude that the fixed point becomes unstable, and an invariant closed curve arises around such point, which is shown in Figure 3(d).

As 𝜎 is further increased, a strange attractor is generated by successive stretching and folding. The fixed point has coordinates 𝑥^∗₂=2.3657 and 𝑦^∗₂=0.6006, and the corresponding eigenvalues are 𝜆=0.53066±𝑖0.88656, with |𝜆|=1.0332. The strange attractor is generated by the breaking of the invariant circles and the appearance of twelve chaotic (not shown in this figure) regions changes as they are linked into a single-chaotic attractor.

4. Controlling through β by OGY Method

In the preceding section we showed that, according to other works in the literature (see [14, 15]) environmental defensive expenditures may generate chaotic behavior which, in turn, may jeopardize environmental sustainability of economic activity. In this section, we show how chaos can be ruled out from the dynamics of our model by an appropriate choice of the visitors fee 𝛽 and of the defensive expenditure 𝜌𝛽. We are interested in modifying the dynamic behavior of the OAPA model by introducing the visitors fee 𝛽 and the defensive expenditure 𝜌𝛽. As it was shown in Figure 2, at the value 𝜎=1.2, the OAPA model presents a stable fixed point. Figure 4 shows the bifurcation diagram of the system (2.5), where parameter 𝛽 varies in the interval [0,0.8] and parameter 𝜌 is posed equal to 0.2. We can obtain both stable dynamics and chaotic dynamics. In fact, starting from a stable fixed point of the OAPA system, for values of 𝛽∈[0,0.42), the system (2.5) admits a stable fixed point, while for 𝛽>0.42 chaotic dynamic occurs.

(a)

(b)

Figure 4 

Bifurcation diagrams for the state variable 𝑥 (a) and for the state variable 𝑦 (b), varying 𝛽. The parameter values are 𝛼=0.12, 𝜎=1.2, 𝛾=0.375, 𝜌=0.2, and 𝑟=2.8.

In this section, we describe a method that allows to stabilize this chaotic dynamic. In order to achieve this goal, the so-called OGY method (see [16]) is used.

The OGY method was successfully used in several studies, both in economics and physics (see e.g., [17, 18]). As it is summarized in [18–20], the OGY method is based on the following assumptions:(a₁)a chaotic solution of a nonlinear dynamic system may have even an infinite number of unstable periodic orbits;(a₂)in a neighborhood of a periodic solution; the evolution of the system can be approximated by an appropriate local linearization of the equation of motion;(a₃)small perturbations of the parameter 𝑝 of the system can shift the chaotic orbit toward the so-called stable manifold of the chosen periodic orbit;(a₄)The points belonging to the stable manifold approach the periodic solution in the course of time;our goal is to find a “good” way to approach the periodic unstable orbit by proper changes of the parameter if the starting point is in a neighbourhood of the periodic unstable orbit.Let us assume that the model can be described as 𝑧_𝑛+1=𝑓𝑧_𝑛,𝑝, (4.1) where 𝑛=1,2,…, 𝑝 is real parameter, 𝑧_𝑛=(𝑥_𝑛,𝑦_𝑛)∈ℜ²,𝑓=(𝑓₁,𝑓₂);(a₅)Suppose that we have a fixed point 𝑧₀=(𝑥₀,𝑦₀) corresponding to a fixed parameter value 𝑝₀ such that 𝑧₀=𝑓𝑧₀,𝑝₀, (4.2) and such fixed point is unstable;(a₆)assume that the Jacobian matrix has two eigenvalues 𝜆₁,𝜆₁ satisfying|𝜆₁|<1<|𝜆₂|.

Then it follows from (𝑎₂) that, starting sufficiently close to 𝑧₀ and 𝑝₀, we can approximate the right-hand side of (2.5) by the first-degree terms of its Taylor expansion around 𝑧₀ and 𝑝₀ Then, by (𝑎₃), modifying 𝑝 we try to shift the chaotic orbit toward a stable manifold.

Thanks to the OGY method, the goal of approaching a stable manifold may be achieved as follows. Let 𝑧_𝑛 and 𝑝_𝑛 be close enough to 𝑧₀ and 𝑝₀ as required in (𝑎₂). Then, the next point of the orbit is determined by (4.1):𝑧_𝑛+1=𝑓𝑧_𝑛,𝑝_𝑛. (4.3) Our aim is to determine 𝑝_𝑛, that is, how to control the system that orbit approaches the unstable fixed point.

From the above results we get the following theorem.

Theorem 4.1. Under the assumptions (a₁)–(a₆), there is a value for 𝑝_𝑛 such that trajectory of the recurrence map (4.1) is shifted towards the stable manifold.

We fix the parameters 𝛼=0.12, 𝛾=0.375, 𝜎=1.2, 𝜌=0.2,𝑟=2.8, and 𝛽=0.745; in such context the system exhibits a chaotic attractor. We take 𝛽 as the control parameter which allows for external adjustment but is restricted to lie in a small interval |𝛽−𝛽₀|<𝛿, 𝛿>0, around the nominal value 𝛽₀=0.745. The system becomes𝑥_𝑡+1=𝑥_𝑡−𝛽−0.375𝑥_𝑡+1.2𝑦_𝑡,𝑦_𝑡+1=𝑦_𝑡+2.81−𝑦_𝑡𝑦_𝑡−0.12𝑥²_𝑡+0.2𝛽𝑥_𝑡. (4.4)

Following [21], we consider the stabilization of the unstable period one orbit 𝑃₂=(𝑥^∗,𝑦^∗)=(1.21738,1.00126). The map (4.4) can be approximated in the neighborhood of the fixed point by the following linear map:𝑥_𝑡+1−𝑥^∗𝑦_𝑡+1−𝑦^∗𝑥≅𝐴^𝑡−𝑥^∗𝑦_𝑡−𝑦^∗+𝐵𝛽−𝛽₀, (4.5) where𝑥𝐴=𝜕𝑓^∗,𝑦^∗𝜕𝑥_𝑡𝜕𝑓𝑥^∗,𝑦^∗𝜕𝑦_𝑡𝜕𝑔𝑥^∗,𝑦^∗𝜕𝑥_𝑡𝜕𝑔𝑥^∗,𝑦^∗𝜕𝑦_𝑡,𝜕𝑓𝑥^∗,𝑦^∗𝑥𝜕𝛽𝜕𝑔^∗,𝑦^∗𝜕𝛽. (4.6) are the Jacobian matrixes with respect to the control state variable (𝑥_𝑡,𝑦_𝑡) and to the control parameter 𝛽. The partial derivatives are evaluated at the nominal value 𝛽₀ and at (𝑥^∗,𝑦^∗). In our case we get𝑥_𝑡+1𝑦−1.21738_𝑡+1≅𝑥−1.001260.6251.2−0.14317−1.80708𝑦^𝑡−1.21738_𝑡+−1.00126−10.24347𝛽−0.75. (4.7) Next, we check whether the system is controllable. A controllable system is one for which a matrix 𝐻 can be found such that 𝐽−𝐵𝐻 has any desired eigenvalues. This is possible if rank(𝐶)=𝑛, where 𝑛 is the dimension of the state space, and𝐶=𝐵∶𝐽𝐵∶𝐽²𝐵∶⋯∶𝐽^𝑛−1𝐵. (4.8) In our case it follows that𝐶=(𝐵∶𝐽𝐵)=−1−0.33280.24347−0.29681, (4.9) which obviously has rank 2, and so we are dealing with a controllable system. If we assume a linear feedback rule (control) for the parameter of the form:𝛽−𝛽₀𝑥(𝑡)=−𝐻−𝑥^∗𝑦(𝑡)−𝑦^∗, (4.10) where 𝐻∶=[ℎ₁ℎ₂], then the linearized map becomes𝑥_𝑡+1−𝑥^∗𝑦_𝑡+1−𝑦^∗≅()𝑥𝐽−𝐵𝐻^𝑡−𝑥^∗𝑦_𝑡−𝑦^∗, (4.11) that is,𝑥_𝑡+1𝑦−1.21738_𝑡+1≅−1.001260.625−ℎ₁1.2−ℎ₂−0.14317+0.2437ℎ₁−1.80708+0.2347ℎ₂𝑥_𝑡𝑦−1.21738_𝑡−1.00126, (4.12) which shows that the fixed point will be stable provided that 𝐴−𝐵𝐻 is that all its eigenvalues have modulus smaller than one. The eigenvalues 𝜇₁ and 𝜇₂ of the matrix 𝐴−𝐵𝐻 are called the “regulator poles,” and the problem of placing these poles at the desired location by choosing 𝐻 with 𝐴,𝐵 given is called the “pole-placement problem”. If the controllability matrix’ from (4.8) is of rank 𝑛, 𝑛=2 in our case, then the pole-placement problem has a unique solution. This solution is given by𝐻=𝛼₂−𝑎₂𝛼₁−𝑎₁𝑇⁻¹, (4.13) where 𝑇=𝐶𝑊, and𝑎𝑊=¹1=101.1820110. (4.14) Here, 𝑎₁and𝑎₂ are the coefficients of the characteristic polynomial of 𝐽, that is|𝐽−𝜆𝐼|=𝜆²+𝑎₁𝜆+𝑎₂=𝜆²+1.1820𝜆−0.9576; (4.15)𝛼₁,and𝛼₂ are the coefficients of the desired characteristic polynomial of 𝐽−𝐵𝐻, that is(𝐽−𝐵𝐻)−𝜇𝐼=𝜇²−𝛼₁𝜇+𝛼₂⟹𝛼₁=−𝜇₁+𝜇₂⟹𝛼₂=𝜇₁𝜇₂. (4.16) From (4.13); we get that𝐻=𝜇₁𝜇₂+0.9576𝜇₁+𝜇₂=−1.1820−0.64437−2.64657−0.023824.00933−0.6444𝜇₁𝜇₂−0.5889+0.02382𝜇₁+0.02382𝜇₂−2.647𝜇₁𝜇₂−7.274−4.009𝜇₁−4.009𝜇₂. (4.17) Since the 2−𝐷 map is nonlinear, the application of linear control theory will succeed only in a sufficiently small neighborhood 𝑈 around (𝑥^∗,𝑦^∗). Taking into account the maximum allowed deviation from the nominal control parameter 𝛽₀ and (4.10), we obtain that we are restricted to the following domain:𝐷_𝐻=𝑥(𝑡),𝑦(𝑡)∈ℜ²∶𝐻𝑥(𝑡)−𝑥^∗𝑦(𝑡)−𝑦^∗<𝛿. (4.18) This defines a slab of width 2𝛿/|𝐻| and thus we activate the control (4.10) only for values of (𝑥_𝑡,𝑦_𝑡) inside this slab, and choose to leave the control parameter at its nominal value when (𝑥_𝑡,𝑦_𝑡) is outside the slab.

Any choice of regular poles inside the unit circle serves our purpose. There are many possible choices of the matrix 𝐻. In particular, it is very reasonable to choose all the desired eigenvalues to be equal to zero, and in this way the target would be reached at least after 𝑛 periods, and therefore, a stable orbit is obtained out of the chaotic evolution of the dynamics.

In Figures 5(a) and 5(b), we show the time series of the chaotic trajectory starting from the point (𝑥₀,𝑦₀)=(0.9,0.8) which we have chosen to control. In contrast, Figures 5(c) and 5(d) present the controlled orbit converging to the stabilized fixed point when the feedback matrix 𝐻 is chosen such that the eigenvalues of (𝐽−𝐵𝐻) are 𝜇₁=𝜇₂=0. This implies that 𝜇₁+𝜇₂=0, 𝜇₁𝜇₂=0 and so 𝐻=(−0.5889,−7.274). For this control strategy, we have also chosen 𝛿=0.1.

(a) Original chaotic orbit of the variable 𝑥(𝑡)

(b) Original chaotic orbit of the variable 𝑦(𝑡)

(c) Controlled chaotic orbit of the variable 𝑥(𝑡)

(d) Controlled chaotic orbit of the variable 𝑦(𝑡)

Figure 5 

Comparison between the original and the controlled orbit.

5. Conclusion

In this paper we studied a discrete-time model that describes the interaction between visitors and the environment resource, in an open-access protected area (OAPA). It was shown that by varying the parameter that indicates the preferences of visitors with reference to the environmental quality, complex dynamics may occur (flip-flop bifurcation, Neimark-Sacker bifurcation, and chaotic dynamics). Furthermore, we analyzed the impact that user fees and environmental defensive choices can have on the OAPA dynamics when it presents an attractive fixed point. Finally, we have applied the OGY control technique (with user fee 𝛽 as control parameter) and we have shown that the aperiodic and complicated motion arising from the dynamics of the model can be easily controlled by small perturbations in their parameters and turned into a stable steady state.