Puu System of Fractional Order and Its Chaos Suppression

Abstract

In this paper, the fractional-order variant of Puu’s system is introduced, and, comparatively with its integer-order counterpart, some of its characteristics are presented. Next, an impulsive chaos control algorithm is applied to suppress the chaos. Because fractional-order continuous-time or discrete-time systems have not had non-constant periodic solutions, chaos suppression is considered under some numerical assumptions.

1. Introduction

The impulsive control concept has a long history that is based on the mathematical foundation of impulsive differential equations. Most examples of impulsive phenomena can be found in mechanical systems with impacts where a sudden change in their states appears, but also in population dynamics, biotechnology processes, chemistry, engineering, medicine, spacecraft optimal control, and so on. Impulsive systems can be studied via the mathematical tool based on impulsive differential or discrete equations. In the last few decades, differential equations with impulses can be found in e.g., nanoelectronic devices, chaotic spread-spectrum communication systems, or electrical engineering applications, and so on (see, e.g., [1,2]). On the other side, there are classes of systems like biological systems, economical systems where discrete time models seem to be more realistic than the continuous ones.

There exist several kinds of impulses [3]—for example, systems where the impulses are applied at fixed time-moments (see, e.g., one of the first references [4]) and also systems where the impulses are applied at variable times [5].

Due to the memory and hereditary properties of the fractional derivatives (see [6,7]), the discrete-time or continuous-time systems of fractional order (FO) are more suitable than integer-order systems.

While the definition of fractional derivative for continuous-time real functions has been formulated in the late 19th Century by Liouville, Grunwald, Letnikov, and Riemann, the first definition of a fractional difference operator was made by Diaz and Olser in 1974 [8]. Nowadays, differential or difference equations of FO represent useful models in viscoelasticity, mechatronics, seismology, aerodynamics, electrical circuits, biophysics, biology, blood flow phenomena, chemistry, control theory, etc (see, e.g., [9]).

In this paper, the fractional order (FO) variant of the Puu’s system is numerically analyzed and compared with its integer order (IO) counterpart. Using the results on non-existence of exact periodic solutions of FO discrete systems, one consider the periodicity as computationally approached. Moreover, using an impulsive algorithm, the chaotic motion can be suppressed.

2. Puu’s Fractional Order System

In 1939, Paul Samuelson [10] introduced the standard principle of acceleration to consumption used in one of the first formal mathematical models for business cycles, while Sir John Hicks (1950) later improved the Samuelson model [11], introducing the “floor” and “roof” to depreciation levels.

Starting from the standard Samuelson-Hicks model, in 1989, Puu and Sushko [12] developed the discrete dynamical income system of integer order (IO) with cubic nonlinearity of integer order, modeled by the following cubic initial value problem (IVP)

$$x(n+1)=ax\left(n\right)-(a+1)x^3\left(n\right),x\left(0\right)=x_0,$$

(1)

with $a>0$ a real parameter.

It is usual to have multiple attractors for a dynamical system. Each of them can be considered as the attractor for a given initial condition within its own attraction basin. In the bifurcation diagram vs. $a\in [0,3]$, presented in Figure 1a, the Puu system of IO (1) reveals two symmetric regions representing regular and chaotic attractors generated by two symmetric (negative and positive) initial conditions within two attraction basins (Due to its symmetry, the coexistence of symmetric attractors is a feature of the Puu system. Like for the cubic logistic map $x\to kx(1-x^2)$, the oddness of the map (1) induces an equivariance that is a ${\mathbb{Z}}_2$-symmetry [13,14].). Therefore, attractors corresponding to $a\in [0,3]$ in the upper (bottom) region present non-symmetric new generated points. The successive bifurcations of these points lead to chaos via the standard period-doubling cascade. In addition, note the sudden crisis—at $a=2.6$, the previous two non-symmetric chaotic attractors (red and blue, respectively) collide and give birth to a symmetric chaotic attractor that covers both positive and negative values of x. Details on similar dynamics but related to the cubic logistic map can be found in [15].

Let us consider the FO Puu discrete system. Following the approach in [16] (see also [17], where the synchronization of the FO Puu is considered), the Caputo-like Puu system of FO can be modeled by the following general IVP of FO:

$${\Delta }_{*}^{q}x\left(k\right)=f\left(x(k+q-1)\right),k\in {\mathbb{N}}_{1-q},x\left(0\right)=x_0,$$

(2)

where ${\Delta }_{*}^q$ is the Caputo-like delta difference of order $0<q\le 1$, and $N_{1-q}=\{1-q,2-q,\dots ,\}$ represents the isolated time scale. With the cubic right-hand side $f\left(x\right)=x-(a+1)x^3$, the IVP (2) reads

$${\Delta }_{*}^{q}x\left(k\right)=ax(k+q-1)-(a+1)x^3(k+q-1),k\in {\mathbb{N}}_{1-q},x\left(0\right)=x_0,$$

(3)

Theorem. 

[16,18] The IVP (3) admits the following discrete integral

$$x\left(n\right)=x\left(0\right)+\frac{1}{\Gamma \left(q\right)}\sum _{j=1-q}^{n-q}\frac{\Gamma (n-j)}{\Gamma (n-j-q)}[ax(j+q-1)-(a+1)x^3(j+q-1)].$$

(4)

In (4), $\Gamma$ represents the Gamma function.

Numerically implemented, via the substitution $i=j+q$, the integral (4) can be expressed as follows [16]:

$$x\left(n\right)=x\left(0\right)+\frac{1}{\Gamma \left(q\right)}\sum _{i=1}^n\frac{\Gamma (n-i+q)}{\Gamma (n-i+1)}[ax(i-1)-(a+1)x^3(i-1)],n=1,2,\dots$$

(5)

The local Lyapunov exponent $\lambda$, can be approximated numerically at the first n iterations with the following relation [19]:

$$\lambda =\frac{1}{n}ln\left(\right|\tau (n-1)\left|\right),$$

(6)

where the tangent map $\tau \left(n\right)$, obtained by the linearization of (5) along the orbit $x\left(n\right)$, is

$$\tau \left(n\right)=\tau \left(0\right)+\frac{1}{\Gamma \left(q\right)}\sum _{j=1}^n\frac{\Gamma (n-i+q)}{\Gamma (n-i+1)}[a\tau (j-1)-3(a+1)x^2(j-1)\tau (j-1)],\tau \left(0\right)=1.$$

Like the solution x, the tangent map $\tau$ at the moment n also presents the so-called discrete memory effect (or time history) i.e., the numerical determined value at the moment n depends on all previous values. Because in relations (5) and (6), the term $R:={\sum }_{j=1}^n\frac{\Gamma (n-i+q)}{\Gamma (n-i+1)}$ presents divergency problems, to deal numerically with R for large values of n, one can use the relation $\frac{\Gamma (n-i+q)}{\Gamma (n-i+1)}=e^{ln(\Gamma (n-i+q\left)\right)-ln(\Gamma (n-i+1\left)\right)}$. In this way, the numerical experiments can be extended for n up to a couple thousand.

The bifurcation diagrams of the Puu system of FO vs. q, generated within $x_0=\pm 0.085$ and $x_0=\pm 0.8$, for $a=1.27$ and $a=1.14$ (Figure 1b,c, respectively), reveal symmetrical transitions to chaos and, also, contrary to its IO part, attractors’ coexistence.

Hereafter, only the two merged “positive” regions (blue and cyan) located most entirely in the positive axis are considered. The other two “negative” symmetric regions (red and light blue) can be generated with $x_0=-0.085$ and $x_0=-0.8$, respectively. As can be seen, both “positive” and “negative” attractors corresponding to $q\in (0,1]$ maintain the odd symmetry existing in the integer order system and the bifurcations in the positive (negative) axis are non-symmetric.

Remark 1.

As for the case of fractional differential equations, contrary to integer order difference equations, fractional difference equations do not admit exact periodic but only asymptotically periodic solutions (see, e.g., [20,21,22], respectively). Therefore, notions like “stable cycle”, “chaos control”, and even “chaos” (as consisting in infinitely many unstable periodic motions), in continuous-time or discrete-time systems of FO seem to not be adequate. Up to a considered enough tiny computational error, the apparent periodic orbits are called here “computationally periodic orbits” (CPOs). Therefore, in this paper, the chaos suppression means to obtain stable CPOs.

To study some orbit $x\left(n\right)$, the relations (5) and (6) together with time series, histograms, and the ‘0–1’ test [23] will be used. Because we are interested in chaos suppression, the value of the parameter a is chosen so that the system evolves chaotic for all initial conditions. Thus, for $a=1.3$, the system behaves chaotically for all initial conditions.

One of the tools generated by the ‘0–1’ test [23], which can be easily numerically implemented starting from a time series related to the considered dynamical continuous or discrete system, is the asymptotic growth rate, K, which gives important information to distinguish chaotic behavior from regular behavior. Note that the ‘0-1’ test does not need the system equations, but only a time series of the system. If $K\approx 1$, then the underlying dynamics are chaotic, while, if $K\approx 0$, the system behaves along some stable CPO. It is commonly assumed that about 2000 elements in the considered time series, after transients are discarded, are enough. In this paper, after the first 1000 iterations neglected, 1500 elements give enough accurate results.

While the integer order system develops, as usual, stronger chaos once the bifurcation parameter a increases (Figure 1a), the FO system presents chaos extinction with the increase of a (Figure 1b,c).

Hereafter, the case $q=0.5$ will be considered.

Another interesting property of the FO Puu system, contrary to its integer order counterpart, is the coexistence of stable CPOs with chaotic attractors. For example, in Figure 2, two zoomed areas are presented (Figure 2b,c) from the bifurcation diagram for $a\in [0,2]$ in Figure 2a. Both zoomed areas reveal, for some parameter values of a, the coexistence of a chaotic attractor with a stable CPO (dotted lines at $a=1.27$ and $a=1.14$).

Consider $a=1.27$ within a stable CP window of period-5 for $a\in [1.22,1.32]$ (red plot) coexisting with its chaotic counterpart (blue plot) in Figure 2b. The coexistence of the underlying stable CPO and the chaotic attractor is studied in Figure 3. Figure 3a,b present superimposed, the time series entirely and the last 100 iterations, respectively; in Figure 3c,d, superimposed, the Lyapunov exponent $\lambda$ and the asymptotic growth rate, K, are represented, while in Figure 3e,f plot the histograms of the coexisting attractors. The chaotic attractor generated from $x_0=0.8$ (blue plot in Figure 3a,b is characterized by the evolution of the time series, the K value, which is close to 1, and the positiveness of the superimposed Lyapunov exponent $\lambda$ (Figure 3d) and the histogram in Figure 3f. The stable CPO revealed after about 900 neglected iterations from the considered 2500 iterations, generated from $x_0=0.085$ (red plot in Figure 3a,b), is characterized by the evolution of the time series, the K value, and the superimposed $\lambda$, in which both are close to 0 (Figure 3c), and the discrete five peaks bars in the histogram (Figure 3e) indicating the branches of the cycle, numbered 1–5 in the zoomed time series in Figure 3b. Note that the zero value of $\lambda$ seems to indicate a kind of Neimark–Sacker bifurcation, generating the stable CPO of period-5 which coexists with the chaotic attractor. In addition, as can be seen in Figure 3c,d, this kind of bifurcation seems to be not singular, but this subject does not represent the purpose of this paper.

3. Chaos Suppression in the FO Puu System

Being a model proposed in economy, where chaos is usually present, chaos suppression could represent an important task. A simple and efficient way used in this paper is to apply constant periodic impulses $\gamma$ every $\Delta$ steps. The numerical form of the impulsive algorithm is

$$x(n+1)=\left\{\begin{array}{c}f\left(x\right(n\left)\right),n\in \mathbb{N},\hfill \\ (1+\gamma )x(n+1),\mathrm{if}mod(n,\Delta )=0,\hfill \end{array}\right.$$

(7)

where $f\in C(\mathbb{R},\mathbb{R})$ (here the cubic map (1)), $\Delta \in {\mathbb{N}}^{*}$ and the impulse $\gamma$ some relative small real number.

Numerically, the control algorithm (7) reads as follows: after the new value of x is calculated at the moment $n+1$, $x(n+1)$, if n is multiple of $\Delta$ steps, $x(n+1)$ is perturbed with $1+\gamma$.

The algorithm is suitable to systems where the state variable is accessible to perturbations and has been successfully utilized for discontinuous systems of FO [24], continuous fractional-order systems of IO [25], discrete systems of FO [26], and discrete systems of IO [27]. The stability of impulsive fractional difference equations is studied in [28]. Some analytical study, such as boundness and periodicity, that applied to a discrete economical supply and demand system can be found in [29].

Consider the case $q=0.5$ and $a=1.3$ when the both coexisting chaotic attractors of the Puu’s system of FO are chaotic (see Figure 2b). This means that for whatever initial conditions $x_0$ the system evolves along one of the two chaotic attractors (red or blue).

To not modify the system structure, $\gamma$ values are chosen to be relatively small, here of $1e-2$ order. The effectiveness of the algorithm (7) can be deduced form the bifurcation diagram of the controlled system vs. $\gamma$ for $\gamma \in [-0.1,0.1]$ and fixed $\Delta$. The stable windows reveal the $\lambda$ values generating stable CPOs. While the chaotic orbits are identified by the positiveness of the Lyapunov exponent, $\lambda$, the chaotic evolution of the underlying time series, dense set of bars in histograms, and $K\approx 1$, the COPs have negative (or zero) Lyapunov exponent, computationally periodic behavior in the time series, discrete bars in histograms and $K\approx 0$.

The case $\Delta =1$, when the system is impulsed at every step, is presented in Figure 4a,d for $x_0$ chosen in the two mentioned attraction basins, $x_0=0.8$ and $x_0=0.085$, respectively. As can be seen in the bifurcation diagrams, the impulsive algorithm applied within the two attraction basins preserves the existence of the two previously existing chaotic attractors, but also generates stable CP windows in the $\gamma$. Obviously, it is desirable to find those $\gamma$ values for which both coexisting chaotic attractors are stable. However, for $\Delta =1$, the periodic windows cannot suppress simultaneously the chaos in both chaotic attractors (see dotted green and black lines). For example, for $\gamma =-0.0375$ (dotted black lines in Figure 4a,d), the chaotic attractor corresponding to $x_0=0.8$ of the impulsed system is stabilized and evolves along a stable CPO of period-5 (see the five intersections of the green dotted line in Figure 4a, the time series in Figure 4b and the history in Figure 4c), but the coexisting attractor corresponding to $x_0=0.085$ of the impulsed system is chaotic for the same value of $\gamma$ (see Figure 4d–f). The stability within the periodic windows is revealed by the Lyapunov exponent (blue plot) and the asymptotic growth rate K (red plot) which are approximately 0.

Concluding, for $\Delta =1$, it is possible that, depending on the initial conditions, the system evolves regularly but also chaotically. This impediment can be avoided if $\Delta =2$, when the system is impulsed only every two steps. Now, both chaotic attractors can be controlled for a large interval of $\gamma$ values (Figure 5). Thus, for $\gamma =-0.09$, the impulsed system evolves, for whatever initial condition, along a stable CPO: a stable CPO of period-8, for $x_0=0.8$ (Figure 5a–c), and also a stable CPO of period-4, for $x_0=0.085$ (Figure 5d–f).

Chaos suppression can also be realized with higher values of $\Delta$. For example, for $\Delta =3$ (Figure 6), if one sends an impulse to the system at every three steps, for $\gamma =-0.086$ (see dotted green line), and one obtains two stable CPOs of multiple periods.

Next, let the averaged energy of the impulsed system, $\overline{E}$, which can be determined with the following relation:

$$\overline{E}=\frac{{\sum }_{n=1}^N{x}^2\left(n\right)}{N},$$

for N sufficiently large (see e.g., [30] for the averaged energy applied to a generalized logistic map).

The variation of the averaged energy vs. $\gamma$, for the considered cases $\Delta =1,2,3$, for $N=2000$, is presented in Figure 7. Note that, if the system is impulsed rarer, i.e., $\Delta =3$, the energy of the impulsed system (red) is smaller than that of the not impulsed system, ${\overline{E}}_0$ (horizontal dotted line) for most of $\gamma$ values, while, if the system is impulsed more often, i.e., $\Delta =1$ (black) or $\Delta =2$ (blue), the necessary energy for the controlled system (vertical dotted lines) is bigger than ${\overline{E}}_0$. In addition, for positive values of $\gamma$ (not considered in this work), $\overline{E}<{\overline{E}}_0$, for all considered values for $\Delta$. As expected, when $\gamma$ is negative, the system increases the energy.

4. Conclusions

In this paper, the FO variant of the Puu system has been considered and some of its properties are characterized comparatively with its IO counterpart. Contrary to the IO variant, the FO variant presents significant symmetries. In addition, contrary to the FO continuous-time systems, the chaos in the Puu system of FO reduces once the fractional order q increases.

Because, as usual for FO systems, the system does not admit periodic orbits, instead “periodic orbit”, the notion of “computationally periodic orbit” has been introduced. In this way, the door of everything that involves “periodicity”, such as stable/unstable cycles, chaos, chaos control, remains open for study. Note that another approach of this kind of orbits would be done via “almost periodicity” notion (see, e.g., [31]).

Using the impulsive algorithm (7), which perturbs periodically the state variable with the value $1+\gamma$, chaos can be suppressed. The values of $\gamma$ can be determined from the bifurcation diagram of the impulsed system vs. $\gamma$. In this paper, the perturbations $\gamma$ are negative, but, depending on the considered system, positive values can also be used (see the mentioned references). Note that, in nature, these kinds of perturbations can be shocks, and can appear in natural disasters, ecology, or can be used in ecosystems management, harvesting, etc.

As for other examples of nonlinear discrete FO systems modeled by Caputo’s derivative, chaos in Puu’s system of FO vanishes when the bifurcation parameter approaches 1. Therefore, the problem of chaos suppression makes sense only for relatively smaller values of the fractional order q, where the considered system behaves chaotically.