A new approach of the partial control method in chaotic systems
aa r X i v : . [ n li n . C D ] F e b A new approach of the partial control method in chaotic systems
Rub´en Cape´ans, Juan Sabuco, and Miguel A. F. Sanju´an
1, 31
Nonlinear Dynamics, Chaos and Complex Systems Group,Departamento de F´ısica, Universidad Rey Juan Carlos,M´ostoles, Madrid, Tulip´an s/n, 28933, Spain Institute for New Economic Thinking at the Oxford Martin School,Mathematical Institute, University of Oxford,Walton Well Road, Eagle House OX2 6ED, Oxford, UK. Department of Applied Informatics, Kaunas University of Technology,Studentu 50-415, Kaunas LT-51368, Lithuania (Dated: February 19, 2019)
Abstract
We present here a new approach of the partial control method, which is a useful control techniqueapplied to transient chaotic dynamics affected by a bounded noise. Usually we want to avoid theescape of these chaotic transients outside a certain region Q of the phase space. For that purpose,there exists a control bound such that for controls smaller than this bound trajectories are kept in aspecial subset of Q called the safe set. The aim of this new approach is to go further, and to computefor every point of Q the minimal control bound that would keep it in Q . This defines a specialfunction that we call the safety function, which can provide the necessary information to computethe safe set once we choose a particular value of the control bound. This offers a generalizedmethod where previous known cases are included, and its use encompasses more diverse scenarios. Keywords: chaos control, transient chaos, time series. . INTRODUCTION Transient chaos is a behaviour found in nonlinear systems where trajectories behavechaotically in a certain region Q of the phase space, before eventually escaping to an externalattractor. In some occasions, this escape involves a highly undesirable state and thereforethe application of some control scheme is required to prevent it.Different control methods have been proposed in the literature [1–4] to achieve this goal.However, these methods sometimes fail dramatically in presence of noise due to the expo-nential growth of small perturbations in chaotic dynamics. To deal with real systems wherethe presence of noise can be unavoidable, it has been recently proposed the partial controlmethod [5, 6]. This method is applied on chaotic maps and it is based on the followingscheme: q n +1 = f ( q n ) + ξ n + u n | ξ n | ≤ ξ and | u n | ≤ u , with ξ > u > . (1)Here, the term f ( q n ) represents the action of the map, while the terms ξ n and u n represent thedisturbance and control acting on the nth iteration of the map. Both, the disturbance andcontrol are bounded so that | ξ n | ≤ ξ and | u n | ≤ u , with ξ > u >
0. These constraints area consequence of the limitations of both the disturbance ands control in most applications.One of the remarkable findings of the partial control method is that, controlled trajectoriesexist for values u < ξ . This means that the changes in the dynamics of the system inducedby the disturbances can be counteracted with the application of a smaller amount of control.Such a counterintuitive result was proven in several paradigmatic systems like the H´enonmap [6], the Duffing oscillator [6] or the Lorenz system [7] as well as other models in thecontext of ecology or cancer dynamics [8, 9].To implement this method, it is necessary to know the map f ( q ) and the disturbancebound ξ . Then specify the region Q where we want to keep the trajectories, and setthe control bound u that we want to apply. By using an algorithm called the SculptingAlgorithm [6] it is possible to found the subset of points q ∈ Q that can be controlled underthe scheme (1). This subset is called the safe set and its shape depends on the choice ofthe bound u .When we compute the safe set, there is a minimum u < ξ for which the safe set exists.That is, the safe set can be computed only for values of u larger than this minimun u .2 igure 1. Partial control method.
In this figure the slope-3 tent map is represented. The mapis affected by a uniform disturbance bounded by ξ = 0 .
04. The small dots help to visualize themagnitude and distribution of the disturbance. (a) An uncontrolled trajectory that escapes fromthe interval [0,1] after a few iterations is shown. (b) The partial control method was applied withthe control bound u = 0 . ,
1] when a control | u n | ≤ .
025 is applied at every iteration tothe x variable. An example of this technique is shown in Fig. 1, where an uncontrolled trajectory and acontrolled one are compared in the case of the slope-3 tent map. This map is given by: x n +1 = x n + ξ n + u n for x n ≤ − x n ) + ξ n + u n for x n > (2)where as an example, we consider a disturbance bound ξ = 0 . Q = [0 ,
1] after a few iterations. InFig. 1(b) the partial control method is applied. For a u = 0 .
025 we found the safe setdisplayed at the bottom of the figure. By forcing the trajectory to pass through this set,the orbit is kept in the interval [0 ,
1] by using at every iteration a control | u n | ≤ . II. A NEW APPROACH
With the aim to extend the applications of the partial control method, a new approachhas been developed. Now, the maps considered are more general and have the following3 [ i = ] f( q [ i = ]) U k u [ i = , j = ] U k [j=9] q [ j = ] Pair ( u[i=3, j=9], U k [j=9] ) = ( ) i n d e x region Q=[0,1]
Figure 2.
Definition of the function U k . In this figure a region Q (where we want to keep thedynamics) and a possible function U k (in blue) are represented. We assumed that the dynamics inthis region has escapes and the control is applied to avoid these escapes. To implement the controltechnique, a grid covering the region Q (in this case N=11 points) was taken. The controlleddynamics is given by q n +1 = f ( q n ) + u n . We identify the starting and arrival point as q [ i ] = q n and q [ j ] = q n +1 respectively. The control corresponding to a point q [ i ] to go to the point q [ j ], isdenoted as u [ i, j ], while the value U k [ j ] represents the control bound for the point q [ j ] to remainin Q the next k iterations. Therefore the pair of values (cid:0) u [ i, j ] , U k [ j ] (cid:1) can be read as (presentcontrol, future control) . Each possible pair, represents a choice of control. This approach evaluatesall choices of control and takes the one with the minimum control bound. As an example, we haveillustrated the starting point q [ i = 3] and the choice of control (cid:0) u [ i = 3 , j = 9] , U k [ j = 9] (cid:1) whichreaches the arriving point q [ j = 9]. form: q n +1 = f ( q n , ξ n ) + u n , (3)where ξ n is a disturbance term (random perturbation), belonging to a bounded distribution.However here, the bound of the disturbance distribution is allowed to be space-dependent,and can act over the variables or parameters of the map. The u n term is the control appliedto the variables of the map with the aim of keeping the trajectory in the desirable region Q .To explain the goal of this approach, suppose that we start with the initial condition4 ∈ Q , and in order to sustain the trajectory in Q during the next k iterations, we applya sequence of control magnitudes ( | u | , | u | , .., | u k | ). However this choice of controls is notunique, and certain strategy should be followed in order to keep these controls low. Whatwe pursue here is to find a control strategy that minimizes the upper bound (the maximum)of these controls.To do that we define in the region Q a special function that we name U k . The value U k ( q ) of this function represents the minimum control bound needed to sustain a trajectory(starting in q ) in the region Q during k iterations. This means that the sequence of k controlsapplied to this trajectory satisfy the condition max( | u | , | u | , .., | u k | ) ≤ U k ( q ). This boundis minimal, so no other controlled trajectory exists with a smaller control bound.The U k , U k − , ..U functions implicitly define the control strategy, so we will focus hereon finding these functions. This finding is not trivial due to the chaotic dynamics presentin the region Q . However, it is possible to obtain them, following an iterative procedure sothat with the initial knowledge of U we can obtain U , U .. etc. To explain the procedure,we consider first the particular case where no disturbances are present in the controlleddynamics, and then we extend the reasoning to the case where the disturbances appear. A. Computing the functions U k in absence of disturbances. When no disturbances affect the system, the controlled map has the form q n +1 = f ( q n ) + u n . We use a grid on Q of N points, and the index i = 1 : N , to identify the startingpoint q [ i ] ≡ q n . Alternatively, we use the index j = 1 : N to denote the arrival point q [ j ] ≡ q n +1 . The controlled map in this grid takes de form q [ j ] = f (cid:0) q [ i ] (cid:1) + u [ i, j ]. Thisis illustrated in Fig. 2, where we have considered the interval [0 ,
1] as the region Q , andwe have selected a grid of N = 11 points. We show an iteration of the map, where thepoint q [ i = 3] maps (control included) to the point q [ j = 9]. The particular control used isrepresented as u [ i = 3 , j = 9]. In the same figure, we also display a hypothetical function U k and its value in the arrival point U k [ j = 9]. The value u [ i, j ] represents the current controlcorresponding to the point i to reach the point j, while the value U k [ j ] represents the controlbound corresponding to the point j to remain in Q for the next k iterations.To illustrate the computation of the U k , the slope-3 tent map shown in Fig. 3 will beused as an example. The region Q selected is the interval [0 , igure 3. Computing the function U → U in the tent map. (a) In the figure a tent-like mapwith no disturbance is represented. We have selected the interval [0 ,
1] in the region Q . The initialfunction U [ i ] = 0, ∀ i , is indicated in blue. Every new value U [ i ] can be computed individually,and the procedure to compute U [ i = 3] is shown in the figure. For that, we need to know theimage f ( q [ i = 3] ) and later, we compute all possible controls u [ i = 3 , j ], which are representedwith horizontal arrows. Afterwards, we build the corresponding pairs ( u [ i = 3 , j ] , U [ j ] ) indicatednear the arrows respectively. The upper bound (or maximum) of each pair is indicated in bold.As we want to minimize the new control bound U [ i = 3], the pair with a minimum upper boundhas to be selected, that is, U [ i = 3] = min ≤ j ≤ N (cid:0) max ( u [ i = 3 , j ] , U [ j ] ) (cid:1) = 0 .
02. (b) The resultingfunction U [ i ], ∀ i is drawn (the values indicated are approximate). escape after one iteration. The idea is to compute recursively the functions U → U → U → ... → U k . Taking into account that U [ i ] represents control bound needed by q [ i ] tokeep its trajectory in Q during 0 iterations, it follows that U [ i ] = 0, ∀ i . This function isdisplayed in blue in Fig. 3(a). For visual convenience, both the tent map and the U functionare represented using the same axes. In the following, we will use this joint representation6 igure 4. Computing the function U → U in the tent map. (a)The previous function U [ i ] is indicated in blue. Every new value U [ i ] can be computed individually, and we show theprocedure to compute the value U [ i = 3]. For that, it is necessary to know the image f ( q [ i ] )and then, compute all possible controls u [ i = 3 , j ], which are represented with horizontal arrows.Then, we build all the pairs ( u [ i = 3 , j ] , U [ j ] ) indicated near the arrows respectively. Theupper bound (or maximum) of each pair is indicated in bold. As we want to minimize the newcontrol bound U [ i = 3], the pair with a minimum upper bound has to be selected, that is, U [ i = 3] = min ≤ j ≤ N (cid:0) max ( u [ i = 3 , j ] , U [ j ] ) (cid:1) . (b) The resulting function U [ i ], ∀ i is drawn (thevalues indicated are approximate). In a similar way, functions U , U ... etc. can be obtained. when the scale axis overlap.To explain how to compute U [ i ], we take for instance, the point q [ i = 3] shown inFig. 3(a). This point maps into f ( q [ i = 3]) and then, all possible controls u [ i = 3 , j ] arecomputed, which are shown in the figure with the horizontal arrows at the bottom. For eachcontrol, the corresponding pair ( u [ i = 3 , j ] , U [ j ] ) is also indicated. This pair can be read as (present control, future control) , so that the pair that minimizes the overall control will be the7 igure 5. Functions U → U → U → U → U in the slope-3 tent map. (a) The slope-3tent map where the region Q selected is the interval [ − . , . U (in blue) was computed. (b) The successive functions U k (starting with U ) thathave been computed to obtain U . pair with the minimum bound. In this case, the pair ( u [ i = 3 , j = 6] , U [ j = 6]) = (0 . , . U [ i = 3] = 0 .
02. This value represents the minimumupper control bound for just one iteration. In general, the values of the function U can befound as U [ i ] = min ≤ j ≤ (cid:0) max ( u [ i, j ] , U [ j ] ) (cid:1) .The resulting function U is displayed in Fig. 3(b). It can be seen that the central pointsof Q maps outside Q , and therefore they need a big control to return to Q in just oneiteration. Therefore a central peak appears in the function U .Once we have U , the function U can be computed following the same process (see Fig. 4).Taking again the initial point q [ i = 3], the action of the tent map f ( q [ i = 3] ) is shown inthe figure. Then a control u [ i = 3 , j ] is applied. All possible pairs ( u [ i = 3 , j ] , U [ j ]) areindicated. In this case, the pair ( u [ i = 3 , j = 5] , U [ j = 5]) = (0 . , .
15) marked in redhas the minimum bound (0.15). This value represents the minimum upper control boundfor 2 iterations. Therefore U [ i = 3] = 0 .
15. In general the values of the function U can befound as U [ i ] = min ≤ j ≤ (cid:0) max ( u [ i, j ] , U [ j ] ) (cid:1) . In Fig. 4(b) the function U is shown.8quivalently, we compute U , U ... etc. In general, in absence of any disturbance, wehave the following recursive formula to compute the functions U k : U k +1 [ i ] = min ≤ j ≤ N (cid:0) max ( u [ i, j ] , U k [ j ] ) (cid:1) i ≡ index of the starting point q [ i ] , i = 1 : N. where N = total number of grid points .j ≡ index of the arrival point q [ j ] , j = 1 : N (4)starting with U [ i ] = 0 ∀ i . Note that the values of u [ i, j ] remain unchanged for everyiteration of the algorithm, so they only need to be calculated once. In Fig. 5, we displaythe process for the slope-3 tent map. The region Q has been selected to be the interval[ − . , . U → U → U → U → U are shown. B. Computing the functions U k in presence of disturbances The extension of the recursive algorithm in the case of systems affected by disturbancesis rather straightforward. Now the dynamics is given by q n +1 = f ( q n , ξ n ) + u n , where ξ n isthe disturbance term belonging to a bounded distribution.In Fig. 6, we illustrate the case of a map affected by a bounded disturbance distribution.The main complication here is that, due to the disturbance, the same point has multipledisturbed images. This number can be infinite and therefore, a discretization must be takento perform the computations (see the red dots in Fig. 6). Given a point q [ i ], we denote thegrid of possible images as f ( q [ i ] , ξ [ s ]), where s = 1 : M i is the index of every individualdisturbance. The number of disturbed images M i can take different values depending onthe particular point q [ i ]. The control corresponding to the point q [ i ] and affected by thedisturbance ξ [ s ], to reach the point q [ j ], is denoted as u [ i, s, j ].Now, to compute the functions U k in presence of disturbances, we follow a similar rea-soning as in the case where there are no disturbances. However, in this case we also have to9 [ i ] f ( q [ i ] , [ s ] ) y = x q n+1 q n U k [ s ] escape u [ i , s , j ] U k [j] q [ j ] Figure 6.
Scheme of a map affected by a bounded disturbance distribution.
The extensionof the algorithm in the case of maps affected by a bounded disturbance distribution, is ratherstraightforward. In this case, given a point q [ i ], to compute the upper control bound U k +1 [ i ], wehave to consider all disturbed images f ( q [ i ] , ξ [ s ]). Then compute all the corresponding controlbounds as in the case of no disturbances, and finally extract the maximum control among themall. evaluate all the disturbed images for a given point q [ i ], and take the maximum among themall to obtain an overall upper control bound. Therefore, the recursive formula in presenceof disturbances is given by: 10 igure 7. Safety function U ∞ for different maps affected by a disturbance. This figureshows how the safety function (in blue) changes depending on the map and the disturbance affectingit. In all cases, the convergence of the safety functions was achieved with 15 iterations or less ofthe algorithm. The maps represented are the following: (a,e) Tent map. (b,f) Logistic map.(c,g) Asymmetric tent map. (d,h) Map with two symmetric hills. The horizontal grey line at x n +1 = 1 . q that map above this line, escape directly from the region Q = [ − . , . ξ = 0 .
05. In contrast, for the maps at the bottom (e,f,g,h), the disturbance bound is ξ = 0 .
2. Note that the safety functions for the bottom maps take larger values, due to the largerdisturbances affecting them. k +1 [ i ] = max ≤ s ≤ M i (cid:16) min ≤ j ≤ N (cid:0) max ( u [ i, s, j ] , U k [ j ] ) (cid:1) (cid:17) i ≡ index of the starting point q [ i ] , i = 1 : N. where N = total number of grid points .s ≡ index of the disturbance ξ [ s ] , s = 1 : M i . where M i = number of disturbed images corresponding with q [ i ] .j ≡ index of the arrival point q [ j ] , j = 1 : N (5)starting with U [ i ] = 0 , ∀ i . Note that in this iterative formula, the u [ i, s, j ] values remainunchanged every iteration of the algorithm. Thus, they only need to be calculated once. III. THE SAFETY FUNCTION U ∞ AND THE SAFE SETS.
In this section, we study the important case where the goal of the controller is to keepthe trajectory in the region Q forever with the smallest control bound. To do that, it isrequired to find U ∞ (that we call the safety function ) and therefore iterate infinite timesthe algorithm. However, if the algorithm converges for a given iteration k so that U k +1 = U k ,then it follows that U ∞ = U k and the iterative process is finished. We do not intend hereto explore the necessary mathematical conditions to achieve the convergence. Our findingis that for the analyzed transient chaotic maps, the algorithm converges in a few iterations.In the next sections some examples supporting this point will be provided.To show different examples of safety functions and how the disturbances affect them, werepresent in Fig. 7 the safety function U ∞ (in blue) for different maps. The maps at the top(a,b,c,d) are affected by the same disturbance bound ξ = 0 .
05. The maps at the bottom(e,f,g,h) are the same respectively, but affected by a bigger disturbance bound ( ξ = 0 . U ∞ has larger values in this case, since a larger control boundis needed to sustain trajectories affected by larger disturbances.To control a trajectory by means of the safety function, we need to specify first the uppercontrol bound u that we want to apply. This value must be chosen so that u ≥ min( U ∞ ).12 igure 8. Extracting the safe sets from the safety function.
Once the safety function U ∞ iscomputed, the safe set is the set of points q ∈ Q that satisfy U ∞ ( q ) ≤ u , where u is the controlbound that we want to apply. In panel (a) we draw the safety function (in blue) correspondingto the map of Fig. 7(a), the tent map affected by an uniform disturbance with bound ξ = 0 . u = 0 .
03 is shown at the bottom (green bars). In panel (d) atrajectory is controlled by applying every iteration a control | u n | ≤ u that forces the trajectoryto pass through the safe set. The controls | u n | applied are shown in panel (g), and for conveniencewe only display the first 100 iterations. Panels (b,e,h) are equivalent but taking instead a controlbound u = 0 .
07. In panels (c,f,i) the control bound is u = 0 .
15. Note that the larger the u value, the larger the safe set, and therefore the trajectory is allowed to explore more points of the Q = [ − . , .
1] region. u is allowed. The set of points q for which U ∞ ( q ) ≤ u constitutes what we call the safe set . Only this set of points can be controlled forever byapplying controls | u n | ≤ u , where in each iteration of the map, u n is chosen to force thetrajectory to pass through the safe set. Very often, the choice of the control u n is not uniqueand therefore multiple controlled trajectories are possible. This makes the method veryflexibleIn Figs. 8(a-b-c) the safety functions corresponding to the maps of the Fig. 7(a) are shown.Different control bounds u were taken and at the bottom the respective safe sets have beendrawn. For each value u , a particular controlled trajectory is shown in Figs. 8(d-e-f). Forclarity, we display only 100 iterations of the trajectory. The control applied in every iterationof these trajectories is represented at the bottom of Figs. 8(g-h-i) respectively. A. Application to the tent map affected by asymmetric disturbances
In the previous examples, we have considered maps where the disturbance ξ n affectingthe trajectories were uniformly bounded so that | ξ n | ≤ ξ ∀ x ∈ Q . However there is noimpediment to apply the algorithm in case of non-uniform disturbance bounds. To showan example, we consider the slope-3 tent map affected by a non-uniform disturbance. Thesystem is given by: x n +1 = x n + ξ n (4 x n −
3) + u n for x n ≤ − x n ) + ξ n (4 x n −
3) + u n for x n > , (6)where the term ξ n (4 x n −
3) models the asymmetric disturbance distribution (see Fig. 9).This particular choice of disturbance was made on purpose to show the particular shape ofthe function U ∞ . For this map, the fixed point x ∗ = 0 .
75 is affected by a zero disturbance,and therefore it needs zero control since f ( x ∗ ) = x ∗ . For this reason, we expect that thesafety function evaluated in the fixed point takes the value U ∞ ( x ∗ ) = 0.We have chosen a uniform grid of 1000 points in the region Q = [ − . , . U ∞ , which is shown in Fig. 9. We can observe that U ∞ has a minimum in the fixed point x = 0 .
75. This minimum control is virtually zero, aswe expected. In the right panel of Fig. 9 different controlled trajectories are displayed forincreasing control values u . Note that with the control bounds u = 0 .
12 and u = 0 .
06 thetrajectory behaves chaotically (affected by the disturbances), while in the case of u = 0 . igure 9. Asymmetric disturbance affecting the map.
On the left panel we show the slope-3tent map affected by an asymmetric disturbance. In particular the disturbance was set to have azero value in the fixed point x ∗ = 0 .
75 highlighted with a circle. This choice was made on purpose toshow that the safety function U ∞ (in blue) has a minimum in this point, since no control is neededto keep the trajectory in the fixed point due to the zero disturbance affecting it. Three differentcontrol bounds u (red horizontal bounds) have been tested. On the right panel we represent thecorresponding controlled trajectories. Depending on the control bound the qualitative behaviorchanges drastically. A small control keeps the trajectory around the fixed point, while larger u values let the trajectory explore other points of the region [-0.1,1.1]. the trajectory remains in the fixed point. This interesting result could be used by thecontroller to change the qualitative behavior of the trajectory, just varying the control value u . B. Application to the H´enon map
In order to compute a two-dimensional safety function, we use here the H´enon map,defined as: x n +1 = a − by n − x n y n +1 = x n . (7)15 igure 10. Uncontrolled trajectory in the H´enon map.
The H´enon map for the parametervalues a = 2 .
16 and b = 0 . ξ = 0 .
1. The bluedot is the initial condition and thr reds dots describe a chaotic transient path.an eventually escapes.The dot marked with a cross is the last iteration of the trajectory in the square Q = [ − , × [ − , Q . This map shows transient chaos for a wide range of parameters a and b . Here we havechosen the parameter values a = 2 .
16 and b = 0 .
3. For these values, the trajectories withinitial conditions in the square [ − , × [ − ,
4] have a short chaotic transient, before finallyescaping this region towards infinity (see Fig. 10).In this example, we consider a situation where the variables ( x, y ) are affected by auniform and bounded disturbance ( ξ xn , ξ yn ) so that k ξ xn , ξ yn k = p ( ξ xn ) + ( ξ yn ) ≤ ξ . To keepthe orbits in Q = [ − , × [ − , u xn , u yn ) also bounded k u xn , u yn k≤ u .The controlled dynamics of the system is then given by: x n +1 = a − by n − x n + ξ xn + u xn y n +1 = x n + ξ yn + u yn . (8)We have applied the extended partial control algorithm with a disturbance bound ξ =0 .
10, obtaining the safety function U ∞ shown in Fig. 11. The logarithm of U ∞ has beenplotted for a better visualization. The minimum of U ∞ is found at the value 0.07. In thefigure it has been represented a controlled trajectory (red dots) obtained by setting a controlbound u = 0 .
08. The controlled trajectory remains in the square [ − ,
4] forever.16 R log(U ꝏ ) Figure 11.
The 2D safety function for the H´enon map.
Taking a uniform disturbancebounded by ξ = 0 . Q = [ − , × [ − , U ∞ has been computed. The algorithm takes 13 iterations to converges in agrid of 2000 × .
08. The logarithm of U ∞ is shown here to enhance the visualization. We represent a controlled trajectory (in red) with acontrol bound u = 0 .
08. This trajectory never abandons the square Q = [ − , × [ − , C. Application to a time series from an ecological system.
In this example, we have worked with an ecological model that describes the interactionbetween 3 species: resources, consumers and predators. The interest of this model lies in thefact that, for some choices of parameters, transient chaos appears involving the extinction ofone of the species. Without no control, the system evolves from a situation where the threespecies coexist towards a state where just two species survive, while predators get extinct.The model that we have used is an extension of the McCann-Yodzis model [10] proposedby Duarte et al. [10, 11], which describes the dynamics of the population density of aresource species R , a consumer C and a predator P . The resulting model is given by the17 igure 12. Dynamics of the extended McCann-Yodzis (Eqs. 9) . Depending on the values ofthe parameters different dynamics are possible. (a) Before the boundary crisis ( K = 0 . σ = 0),there are two possible attractors depending on the initial conditions: one chaotic attractor wherethe three species coexist, and one limit cycle where only the resources and consumers coexist. (b)The case treated here, for values ( K = 0 . σ = 0 . b ). The predators eventually get extinct. following set of nonlinear differential equations: dRdt = R (cid:18) − RK (cid:19) − x c y c CRR + R dCdt = x c C (cid:18) y c RR + R − (cid:19) − ψ ( P ) y p CC + C (9) dPdt = ψ ( P ) y p CC + C − x p P. Depending on the parameters values, different dynamical behaviours can be found (seeFig. 12). Following [11] we have fixed the model parameters : x c = 0 . y c = 2 . x p = 0 . y p = 2 . R = 0 . C = 0 . K = 0 .
99 and σ = 0 .
07. For these valuestransient chaos appears, and the predators eventually get extinct as shown in Figs. 12(b)and 12(c).With the aim of avoiding the extinction, we have computed the safety function. To dothat, first we have discretized the dynamics to obtain a map. It is straightforward to builda map taking a Poincar´e section that intersects the flow. In this case, we have chosen the18 igure 13.
Building the map from several trajectories.
It is possible to discretize thedynamics of the ecological model by taking a Poincar´e section. In this case, we have chosen thesection with C = 0 .
24 as shown on the left panel. With the set of points ( R n , C n , P n ) intersectingthe plane, it is possible to build a return map of the form P n +1 = f ( P n ) as represented on the rightpanel. As the trayectories escape towards values P → R n and C n in the Poincar´e section remain practically constant sothat only the values P n will be controlled. plane C = 0 .
24 as shown in Fig. 13(a). For this Poincar´e section the intersection of the planeand the flow, gives us a set of points ( R n , C n , P n ) that is approximately one-dimensional.Note that C n has a constant value equal to 0.24, and the variable R n is practically constant.Therefore it is possible to construct a return map of the form ( P n , P n +1 ) and control thesystem just perturbing the variable P n . Due to the finite escape time of the transient chaotictrajectories, several trajectories were simulated (displayed with different colors in Fig. 13)to obtain a representative return map.We consider here two different cases. First, a situation where the trajectories are affectedby continuous noise in the variables. Second, the case where a continuous noise is affectingthe parameter K of the system. We want to point out here the difference between themeanings of disturbance and noise. In our convention, the disturbance term only appearsin the map and represents the amount of uncertainty measured in this map. In this sense,the disturbance is the product of the accumulated noise along the trajectory during one19teration of the map. The controlled scheme is given by: P n +1 = f ( P n , ξ n ) + u n , (10)where ξ n is a particular disturbance whose bound ξ may be space-dependent.In the first scenario, the trajectories were obtained by using a RK4 integrator with aGaussian noise affecting the variables ( R, C, P ). In Fig. 14(a) the return map obtained via3000 intersections of the trajectories with the Poincar´e section is shown. With these pointsit is possible to reconstruct the map including the disturbance. Note that in this sense,noise removal techniques are useless here since we want to include the disturbances (theaccumulated noise measure in the map). To do that, different statistical techniques canbe used. One very powerful is the bootstrapping technique that allows the estimation ofthe sampling distribution of almost any statistic using random sampling methods. Howeverfor simplicity, we use here a quantile regression technique to estimate the upper and lowerbounds of the map. Taking the quantile values 0.01 (lower bound) and 0.99 (upper bound)we obtain the two red curves shown in Fig. 14(a). The gap between the two curves containsthe disturbed points corresponding to each P n value. We can see that the disturbance gapis rather uniform in this case.In order to avoid the extinction of predators, the region Q selected to keep the trajectoryis the interval [0 . , . U ∞ shown in Fig. 14(b). The minimum of thisfunction corresponds to the value 0.010. Taking a control bound u = 0 .
011 a trajectorywas controlled using the corresponding safe set. Only the variable P n needs to be controlledsince R n and C n remain practically constant. In Fig. 14(c), 500 iterations of the controlledtrajectory are displayed. Every time the Poincar´e section is crossed, a suitable control | u n | ≤ .
011 is applied. As a result, the extinction of the predators is avoided and the 3species coexist in a stable chaotic regime.In the second situation, we consider a small Gaussian noise affecting the parameter K of the system. This noise affects continuously K and it has been included in the integra-tor. Proceeding in a similar way to the previous case, we obtain the return map shownin Fig. 15(a). It can be appreciated that, in comparison with the first scenario, the dis-turbance interval (gap between red lines) is smaller and less uniform. Therefore the U ∞ a)(b) (c) u u =0.011 Figure 14.
Continuous noise affecting the variables. (a) Return map obtained by means of3000 intersections of the trajectories with the Poincar´e section. A continuous noise is affecting thevariables (
R, C, P ) and it arises in the return map as a stripe. Red lines represent the quantileregression calculated for quantiles 0.01 and 0.99. The gap between the red lines represent thedisturbance bound. In this case the gap is rather uniform in all the map. (b) Taking the region Q as the interval [0 . , . U ∞ (in blue) has been computed obtaining aminimum value of 0 . u = 0 . | u n | ≤ .
011 is applied to putthe orbit again in the nearest point P n with U ∞ ( P n ) ≤ . function, which is shown in Fig. 15(b) is quite different. In the Fig. 15(c) a controlled tra-jectory is displayed, for which we have used a control bound u = 0 . u n exceeds the control bound u . However, due to the Gaussian noise (not bounded)affecting the parameter K , it may happen that at certain iteration we need an extra control.For example, if we work with a map affected by a normal disturbance distribution, and webound it with a three-sigma interval, the safety function U ∞ obtained and the upper bound u selected, will be valid the 99 .
7% of the times. The rest of iterations (0 , a)(b) (c) u =0.006 Figure 15.
Continuous noise affecting the parameter K . (a) Return map obtained by meansof 3000 intersections of the trajectories with the Poincar´e section. A continuous noise is affectingthe parameter K and it arises as a stripe in the return map. Red lines represent the quantileregression calculated for quantiles 0.01 and 0.99. The gap between the red lines represent thebound of the disturbance. The gap in this case is not uniform, since some points P n are affectedby bigger disturbances than others. (b) Taking the region Q as the interval [0 . , . U ∞ (in blue) has been computed obtaining a minimum value of 0 . u = 0 . | u n | ≤ .
006 is applied to put the orbit again in the nearest point P n , with U ∞ ( P n ) ≤ . we know how safe is every point q ∈ Q , this suitable control can be chosen efficiently. IV. CONCLUSIONS
We have presented here a new algorithm in the context of the partial control method.This method is applied to maps of the form q n +1 = f ( q n , ξ n )+ u n , where ξ n is the disturbanceand u n the control. Given a region Q where the dynamics presents an escape, the methodcalculates directly the minimum control bound needed to sustain a trajectory in the region Q forever. To do that, we have introduced the safety function U ∞ that can be computed22hrough a recursive algorithm. This function characterizes every state q ∈ Q and tell us howmuch effort is required to control it. Once the safety function is computed, we only needto pick a bound u ≥ min( U ∞ ). Controlled trajectories are possible by applying a suitablecontrol | u n | ≤ u every iteration.The new partial control algorithm has been proven in the one-dimensional tent map andthe two-dimensional H´enon map, under a non-uniform and a uniform disturbance boundrespectively. We have also applied the control method to a continuous ecological systemwhere one of the species eventually gets extinct via a boundary crisis. Two different scenarioswere studied, a continuous noise affecting the variables, and a continuous noise affectingone parameter of the system. In both cases the safety function U ∞ was obtained and thetrajectories controlled, avoiding the extinction.We show that the use of the safety functions U ∞ makes this partial control approach veryrobust and specially useful in the case of experimental time series. Although the methodwas presented here to avoid undesirable escapes in chaotic transient dynamics, we believethat this method can be extended, under minor modifications, to other interesting scenarios. ACKNOWLEDGMENTS
This work was supported by the Spanish State Research Agency (AEI) and the EuropeanRegional Development Fund (FEDER) under Project No. FIS2016-76883-P. [1] Schwartz IB, Triandaf I. 1996 Sustainning chaos by using basin boundary saddles.
Phys. Rev.Lett.
77, 4740-4743.[2] Dhamala M, Lai YC. 1999 Controlling transient chaos in deterministic flows with applicationsto electrical power systems and ecology.
Phys. Rev. E
59, 1646-1655.[3] Bertsekas DP. Infinite-time reachability of state-space regions by using feedback control. 1972
IEEE Trans. Autom. Control
17, 604-613.[4] Bertsekas DP and Rhodes IB. 1971 On the minimax reachability of target set and targettubes.
Automatica
7, 233-247.[5] Sabuco J, Zambrano S, Sanju´an MAF, Yorke JA. 2012 Dynamics of partial control.
Chaos
6] Sabuco J, Zambrano S, Sanju´an MAF, Yorke JA. 2012 Finding safety in partially controllablechaotic systems.
Commun. Nonlinear Sci. Numer. Simul.
17, 4274-4280.[7] Lorenz E. 1963 Deterministic nonperiodic flow.
J. Atmos. Sci.
20, 130-141.[8] Lop´ez AG, Sabuco J, Seoane JM, Duarte J, Janu´ario C, Sanju´an MAF. 2014 Avoiding healthycells extinction in a cancer model.
J. Theor. Biol.
Ecol. Complex . 19, 1-8.[10] McCann K, Yodzis P. 1995 Bifurcation structure of a three-species food chain model.
Theor.Popul. Biol.
48, 93-125.[11] Duarte, J., Janu´ario, C., Martins, N., Sardany´es., J., 2009. Chaos and crises in a model forcooperative hunting:a symbolic dynamics approach.
Chaos
58, 863-883.58, 863-883.