Poincaré maps for detecting chaos in fractional-order systems with hidden attractors for its Kaplan-Yorke dimension optimization

1.   Introduction

Chaos theory is the branch of mathematics that deals with complex behavior emerging in dynamical systems ^[1]. Chaotic systems possess special features such as extreme sensitivity to initial conditions, random-like unpredictable behavior and ergodicity. Improving these characteristics enhance applications in almost every field of science and engineering, such as: Biology ^[2], astrophysics ^[3], mechanics ^[4], economics ^[5], secure communications ^[6], cryptography ^[7], robotics ^[8], control ^[9] and so on. With the study of fractional calculus, it has been found that Fractional-Order (FO) derivatives serve as an extra-degree freedom for enhancing chaos-based applications, this thanks to the FO derivatives that have inheritance or memory properties that preserve complex phenomena with more accuracy ^[10,11]. For instance, in cryptography and secure communications, FO chaotic systems can generate larger keyspace than their integer counterpart ^[12,13].

According to the literature, attractors of dynamical systems can be divided into two categories depending on the localization of their equilibrium points in the corresponding basin of attraction ^[14]. On the one hand, chaotic systems with self-excited attractors are those whose basin of attraction intersects with an unstable equilibrium. On the other hand, chaotic systems with hidden attractors have a basin of attraction that does not intersect with any open neighborhoods of equilibria. Since the first hidden attractor was found in the classical Chua's circuit ^[15], much more attention has been shifted towards this new type of attractors. Firstly, the dynamics of hidden attractors are highly complex to analyze since they can lead to unexpected behavior. For instance, complex phenomena such as "Multistability" ^[16], which refers to when two or more different attractors coexist under the same parameters are likely to appear in dynamical systems with hidden attractors. The multistability phenomena depend on the initial conditions to switch between one attractor to a totally different one, and this property can be used as an additional source of randomness for potential applications. For instance, but not limited to, they might enhance the performance on cryptography ^[17], and secure communications ^[18].

There exist few studies related to FO chaotic systems with hidden attractors, but FO derivatives can generate several families of hidden chaotic attractors in their commensurate and incommensurate models where multistability phenomena can also be observed ^[19,20]. Therefore, the research effort oriented to FO chaotic systems with hidden attractors is an area of opportunity to enhance chaos-based applications. One of the critical factors that guarantees the success of those applications, is the complexity degree of the chaotic attractor that can be quantified by its Kaplan-Yorke dimension ($D_{KY}$). In this manner, the application of metaheuristics, such as: Differential Evolution (DE) and Accelerated Particle Swarm Optimization (APSO) algorithms, has gained interest in recent years for solving optimization problems of FO chaotic systems ^[21]. In ^[22] the optimization through metaheuristics of the $D_{KY}$ is achieved on four FO chaotic systems with self-excited attractors, a similar approach is presented in ^[23]. To the best of the authors' knowledge, there are no works that explore the optimization of FO chaotic systems with hidden attractors.

Motivated by the above discussion, this paper aims to introduce the optimization of FO chaotic systems with hidden attractors. One of the challenges in the optimization process, is that hidden attractors have small basin of attraction ^[24,25,26], and thus, tiny variations in the parameters, and fractional-order derivatives or initial conditions can easily shift the chaotic behavior to a periodic dynamics. Hence, careful considerations needs to be taken to ensure that the optimization process is performing correctly. The literature reports several methods for detecting Chaos in dynamical systems, such as phase portraits, bifurcation diagrams, equilibrium points analysis through Shilnikov theorem, Lyapunov exponents (LEs), Poincaré maps, among others. It is well-known that LEs provide quantitative information over the stretching and folding of the chaotic attractors. However, accurate computation of Lyapunov exponents requires the analysis of large time series by specialized library packages such as TISEAN ^[27], which can be high computational expensive, but LEs are necessary for calculating the $D_{KY}$ of a chaotic attractor. For such problem, one can introduce constraints that ensures the chaotic motion for its LEs computation through TISEAN. In this manner, Poincaré maps provide a qualitative understating of the dynamical characteristics of a $n$-dimensional system into a $(n-1)$ representation ^[28]. This is done by the intersection analysis of a hyper-plane transverse to the flow of the dynamical system, then the Poincaré map corresponding to a chaotic system should display enough distribution over a 2-D plane. In this regard, the contribution of this work is to detect chaotic motion in FO chaotic systems with hidden attractors through the distribution quantification of Poincaré maps. Moreover, to save execution time, the optimization algorithms are implemented by using the Numba just-in-time (jit) compiler in the Python programming language.

The rest of the paper is organized as follows: Section 2 summarizes the background of this work. DE and APSO algorithms are described in Section 3. Section 4 shows the numerical analysis of the FO chaotic systems. Section 5 details the optimization process whereas the discussion of the optimization results is presented in Section 6. Finally, the authors' conclusions are given in Section 7.

2.   Mathematical background

2.1.   Grünwald-Letnikov method

Although one can cite various definitions of FO derivatives such as as the Grünwald-Letnikov ^[29], Riemman Liouville ^[30] and Caputo definitions ^[31], these mathematical interpretations do not allow a direct understanding of the derivative solution. Therefore, the solution of FO differential equations are often obtained through numerical approximations accordingly the used fractional derivative definition ^{[32,33,34,35]}. For instance, the Adams Bashforth Moulton (ABM) method is a predictor-corrector scheme highly recognized as a valuable fractional calculus tool due to the good approximation that it produces ^[36]. Drawbacks of this scheme are related to the complexity for representing the memory properties of FO derivatives, that requires from long time computation. The Grünwald Letnikov (GL) method is another numerical scheme very suitable for solving FO systems. It can be considered as a generalization of the Euler method for solving ordinary differential equations ^[37]. The GL method is of low computational cost and has approximately the same order of accuracy and good match of numerical solutions as ABM ^[38]. Therefore, this work uses the Grünwald-Letnikov method for solving the FO differential equations that comprises the FO chaotic systems that are the case study. The mathematical description is given below.

According to the explicit numerical approximation of the Grünwald-Letnikov derivative ^[39], the FO derivative $0 < q < 1$ for a discrete function $f(t_k)$, can be described by (2.1), $h$ is the integration step-size, and $C^q_j$ are binomial coefficients which are recursively calculated by (2.2).

The solution of a FO system in the form $D^q_{t_k}x(t_k) = f(x(t_k))$, can be obtained by (2.3), which ideally requires infinite memory length for numerical simulation since the summation term depends on the discretized time $t_k$.

2.2.   Poincaré maps

One of the features of chaos behavior is the existence of a dense set of periodic orbits, which implies that any periodic trajectory of the orbit visits an arbitrarily small neighborhood of a non-periodic one ^[40]. In this manner, and since Poincaré maps preserve many of the orbits of dynamical systems, it is a viable tool for detecting chaos in both integer and FO systems ^[41,42]. Below is described the elaboration of a Poincaré map.

Lets consider a FO system of the form given in (2.4), where $0 < q < 1$, $x(t) = [x_1(t), x_2(t), ..., x_n(t)]^T\in\mathbb{R}^n$ and $f\big(x(t)\big) = [f_1\big(x(t)\big), f_2\big(x(t)\big), ..., f_n\big(x(t)\big)]^T$.

The Poincaré map of the system given in (2.4) is obtained as follows:

$(1)$ Without loss of generality, let $\Sigma\in\mathbb{R}^{n-1}$ be the Poincaré section of system (2.4) such that $\Sigma\cap x(t) = [x_1(t), x_2(t), ..., x_{n-1}(t)]^T$.

$(2)$ Next, the Poincaré section is settled on a given $x_n = \sigma$.

$(3)$ Finally, all states $[x_1(t), x_2(t), ..., x_{n-1}(t)]^T$ are captured each time the dynamical flow intersects the Poincaré section.

The captured points by the above procedure portrait the Poincaré map of system (2.4), then, the behavior of the system is obtained by analyzing the distribution of the distinct points of the Poincaré map. If the Poincaré map consists neither of a small number of points or filling a continuous line but instead a large and irregular distribution, then it is a strong indicator of chaos ^[43].

2.3.   Kaplan Yorke dimension

The Kaplan Yorke dimension $D_{KY}$ describes the fractal dimension of an attractor by means of its Lyapunov exponents as follows ^[44]:

herein, the Lyapunov exponents are sorted from highest to lowest $\lambda_1\geq\lambda_2\geq\cdots\geq\lambda_{j+1}$ and $j\in\mathbb{N}$ is the largest integer for which $\lambda_1 + \lambda_2 +\cdots + \lambda_j\geq0$. The $D_{KY}$ gives a meaningful measure to quantify the complexity of a chaotic attractor ^[45,46]. Therefore, it is the interest of this optimization work to determine a set of parameters and fractional orders that leads to a higher $D_{KY}$.

3.   DE and APSO algorithms

3.1.   Differential Evolution algorithm

The Differential Evolution (DE) algorithm belongs to the family of algorithms based on biological evolution, which uses genetic cross and mutation operators to generate new candidate solutions ^[47]. The structure of DE is illustrated in the chart flow of Figure 1, whose description is as follows:

$(1)$ In the first step, the initial population is generated within the range of interest.

$(2)$ Next, the objective function is evaluated for each individual in the population

$(3)$ In the next step, the mutation operation is carried out:

● (3.1) First, three randomly individuals $X_1$, $X_2$ and $X_3$ are selected.

● (3.2) Then, the following operation is performed:

where $\mu\in\mathbb{R}$ is a scaling parameter.

$(4)$ In the next step, individuals are crossed according to the condition

Specifically, if the random number (rand) is less than or equal to the cross parameter $Cr$, then the term $F_i (j)$ takes the value $w_i (j)$, otherwise it will be $X_i (j)$. Next, the first generation of individuals generates new individuals through mutation and cross operations. Subsequently, the objective function is evaluated in the new individual.

$(5)$ Steps 2–5 are repeated for a number of $k$ generations.

$(6)$ In the last step, the lowest cost function is printed.

3.2.   Accelerated Particle Swarm Optimization algorithm

In the Accelerated Particle Swarm Optimization (APSO) algorithm, the tentative solutions are called particles. This algorithm consists of disturbing the behavior of the particles in such a way that they can move within a search space and finding the solution to the optimization problem ^[48]. The algorithm is inspired by the flight of birds to find food in an unknown field with the lowest possible energy expenditure, which can be translated into finding the minimum value of a function. The structure of the APSO is illustrate in the chart flow of Figure 2, and its description is given as follows:

$(1)$ In the first step, the initial population is generated, meaning the particles position within a range of interest.

$(2)$ Next, the objective function is evaluated for each particle

$(3)$ In the next step, the position of the particles is disturbed:

● (3.1) First, the best particles are selected

● (3.2) Then, the position is updated

Specifically, the new population is generated from the change in position of the particles. The value for $\epsilon_n$ is randomly initialized within the range $[0, 1]$ whereas $\alpha$ and $\beta$ are estimated experimentally.

$(4)$ Steps 2–5 are repeated for a number of $k$ generations.

$(5)$ In the last step, the lowest cost function is printed.

4.   Fractional order chaotic systems

This section shows the numerical simulations and Poincaré maps of the FO chaotic systems that are case study. For instance, the authors in ^[19] introduced the 3D FO chaotic system given in (4.1), where $0 < q_i < 1$, $(i = 1, 2, 3)$ and $a\in\mathbb{R}$. According to expression (4.2), it can be seen that by setting the parameter $a > 1/4$, the 3D FO chaotic system does not posses an equilibrium point ($E_p = (\bar{x}, \bar{y}, \bar{z})$), and thus, any existing basin of attraction is considered as hidden.

By taking into account (2.3), the solution of the 3D FO chaotic system by means of the Grünwald-Letnikov method is given in (4.3), whereas Figure 3 shows the simulation results. Figure 4 shows the LEs spectra by varying the fractional derivative.

To demonstrate the usefulness of Poincaré maps in detecting chaotic behavior, a Poincaré section $\Sigma\in\mathbb{R}^2$ on $y = 2$ is defined as depicted in Figure 5 (a), then, by following the steps of Section 2.2, the Poincaré map of Figure 5 (b) is obtained. It can be seen that the displayed points resembles the $xz$ phase portrait of Figure 3 (b), the markers illustrates if the flow intersects the Poincaré section on $y > 2$ (blue, in-direction) and if it intersects on $y < 2$ (red, out-direction).

It is worthy to mention that the number of points portrayed in the Poincaré map is in relation to the amount of used data. However, the number of points itself says nothing about the dynamical behavior. For instance, lets us consider the dynamical behavior of the system (4.1) illustrated in Figure 6, then, a Poincaré section $\Sigma\in\mathbb{R}^2$ on $y = 2$ is defined as depicted in Figure 7 (a). The obtained Poincaré map in Figure 7 (b) displays only two points, indicating the flow intersects twice in different in-out directions of the Poincaré section, and thus, resembling a periodic behavior. Nevertheless, this is only apparent, due to the fact of the non-existence of exactly periodic solutions in FO systems ^[49,50], thereby, the importance of the distribution quantization of the Poincaré map. The next step is the distribution quantization of the Poincaré map. Herein, a simple procedure is proposed as follows:

$(1)$ Verify that the solution of a system (4.1) remains bounded but does not converge to a fixed point.

$(2)$ Next, the first 10 % samples of the time series are deleted to avoid the initial transient.

$(3)$ Then, the Poincaré map is portrayed with the obtained points when the flow intersects on $y < 2$ (out-direction).

$(4)$ In the next step, only the points that satisfies the condition (4.4) are considered, where $(u_j, v_j)$ and $(u_k, v_k)$ denotes the coordinates of the $P_j$ and $P_k$ points respectively, $p$ is the total number of points intersecting the Poincaré section and $\epsilon = 0.1$, this ensures to preserve only the points distributed over the plane and remove those located closer than $(\epsilon, \epsilon)$,

$(5)$ Finally, the number of points remaining in the Poincaré map is co-related with the dynamical behavior.

For instance, Figures 8 and 9 illustrates the hidden attractors and its corresponding Poincaré maps when system (4.1) behaves quasi-periodic. After applying 1–5 from above procedure, the obtained results in Table 1 shows that the remaining points are 10 and 3 respectively. On the other hand, when system (4.1) behaves chaotic such as in Figure 5, the remaining points in the Poincaré map must portrait $\approx 100$ points.

The second case study is a 4D FO chaotic system introduced in ^[51] as a 4D integer chaotic system with coexisting hidden attractors. Later on, in ^[52] the FO version given in (4.5) was proposed, where $0 < q_i < 1$, $(i = 1, 2, 3)$ and $a\in\mathbb{R}$. According to expression (4.6), it can be seen that for $a\neq 0$, the 4D FO chaotic system does not present an equilibrium point ($E_p = (\bar{x}, \bar{y}, \bar{z}, \bar{w})$), and thus, any existing basin of attraction is considered as hidden.

By taking into account (2.3), the solution of the 4D FO chaotic system by means of the Grünwald-Letnikov method is given in (4.7), whereas Figure 10 shows the simulation results. Figure 11 shows the LEs spectra by varying the fractional derivative.

In the same manner, a Poincaré map is elaborated for analyzing the behavior of the 4D FO chaotic system. Herein, the Poincaré section $\Sigma\in\mathbb{R}^2$ is settled on $z = 0$ as depicted in Figure 12 (a), then, by following the steps of Section 2.2, the Poincaré map of Figure 12 (b) is obtained.

The distribution quantization of the Poincaré map of the 4D FO chaotic system is quite similar to the quantization of the Poincaré map of the 3D FO chaotic system. Herein, only the obtained points when the flow intersects on $z > 0$ (in-direction) are considered with a separation of $\epsilon = 0.1$. This procedure is repeated for different parameters and fractional order for the system (4.5). The obtained results with a number of data of $1.20\times 10^{5}$ are shown in Table 2. It can be seen that for chaotic behavior to exist in the system (4.1), its corresponding Poincaré map must portrait $\approx 40-60$ points.

5.   Description of the optimization process

This section details the maximization process of the $D_{KY}$ for the 3D and 4D FO chaotic systems. Given the nature of the DE and APSO, the optimization problem is given as follows:

where $X\in\mathbb{R}^n$ are the parameters and fractional order derivative values of the FO chaotic system of interest and $f(X)$ returns the $D_{KY}$.

The optimization process is performed applying DE and APSO who are programmed in Python, because it is a multi-paradigm programming language highly versatile for creating all kind of applications and, therefore, it has become one of the most used in recent years. However, Python based applications pays a penalty in the form of speed execution, thereby complex projects can benefit from the use of the Numba JIT compiler ^[53], which allows to translate Python functions to optimized machine code and thus obtaining better speed performance ^[54]. Moreover, the implementation of the Numba JIT compiler only requires from slight modifications in the original Python code and, therefore, maintaining the flexibility of the language. Table 3 shows a comparison of the speed performance between a normal Python code and one Python code using the Numba JIT compiler. In both executions, the argument is to solve the proposed FO chaotic system and elaborate its corresponding Poincaré maps for the chaotic behavior validation. It can be seen, that the optimized Python code performs $\approx 375$ times faster than the non-optimized. Thereby the interest of using this compiler in the optimization process.

Algorithms 1 and 2 show the pseudocodes for the implementation of DE and APSO, respectively. It can be seen that the implementation is quite similar in both algorithms. The Numba JIT compiler is used in most of the Python code with the exception of the Lyapunov exponents (LEs) computation $D_{KY}$ computation, in this step, the computation of the Lyapunov exponents $\lambda_i$, $(i = 1, 2, ..., n)$ is made with the TISEAN package by removing the first $10\%$ samples of the time series.

6.   Optimization results

The DE and APSO were executed for a 100 population and for 50–100 generations. Table 4 shows the ranges of interest and the configuration of the metaheuristics for the optimization of the 3D and 4D FO chaotic systems, respectively. The values of the parameters and fractional-order derivatives in the search space, were encoded within three decimal places. For convenience, the initial population is set with parameters and fractional-order derivatives that lead to chaotic behavior and, therefore, the elaboration of Poincaré maps is not necessary for this step. In Tables 5 and 6 are shown the optimization results of the 3D and 4D FO chaotic systems by DE and APSO, respectively. It can be seen that in all cases, the obtained results lead to higher $D_{KY}$ than those shown in Tables 1 and 2, thus validating the optimization process. The best results for the 3D FO chaotic system are those obtained APSO by using 50 generations whereas the best results for the 4D FO chaotic system are obtained with DE by using 100 generations.

7.   Conclusions

In this paper, the optimization of the $D_{KY}$ of a 3D and a 4D FO chaotic systems was performed by applying two metaheuristics, DE and APSO. The main contribution of this work is the use of Poincaré maps for validating the chaotic behavior in FO chaotic systems with hidden attractors. Herein a simple numerical procedure was proposed for the analysis of points distribution in Poincaré maps. It has been shown that with this simple procedure, it is possible to differentiate a chaotic attractor from a periodic or quasi-periodic attractor. Moreover, the use of the Numba JIT compiler in the Python programming language has been an essential tool to save execution time, while preserving the versatility in the implementation of the metaheuristics. This last part is essential and it has been demonstrated that it facilitates the optimization process of FO chaotic systems with hidden attractors. Finally, it is worthy to mention that Poincaré maps can be used in both FO and integer chaotic systems regardless of the attractor nature (hidden or self-excited). The proposed approach can be applied in FO chaotic systems with self-excited attractor as done in integer chaotic systems with hidden and self-excited attractors.

Acknowledgments

Daniel Clemente-López acknowledges CONACYT for PhD scholarship with CVU no. 859036.

Conflict of interest

The authors declare no conflicts of interests.