Mathematical Model and FPGA Realization of a Multi-Stable Chaotic Dynamical System with a Closed Butterfly-Like Curve of Equilibrium Points

Abstract

This paper starts with a review of three-dimensional chaotic dynamical systems equipped with special curves of balance points. We also propose the mathematical model of a new three-dimensional chaotic system equipped with a closed butterfly-like curve of balance points. By performing a bifurcation study of the new system, we analyze intrinsic properties such as chaoticity, multi-stability, and transient chaos. Finally, we carry out a realization of the new multi-stable chaotic model using Field-Programmable Gate Array (FPGA).

1. Introduction

Chaotic dynamical systems are applicable in a wide range of domains such as encryption [1,2,3], secure communication [4,5,6], robotics [7,8,9], lasers [10], chemical reactors [11], memristors [12,13,14], etc. Chen et al. [1] reported an optical image encryption algorithm in gyrator transform domain. Zhang [2] devised a new unified image encryption system with identical encryption and decryption algorithms. Hua et al. [3] developed a new color image encryption algorithm using a 2-D logistic tent modular map. Samimi et al. [4] designed a chaotic secure communication system with the help of a brain emotional learning-based intelligent controller. Abbasi et al. [5] proposed a chaotic evolutionary amino acid codons Model for image encryption. Alemami et al. [6] discussed chaos-based multiple encryption algorithms for speech recordings. Chaotic systems find use in secure communication devices and engineering applications featuring the implementations of analog or digital electronics [15].

Many types of chaotic dynamical systems have been studied in the literature such as jerk systems [16], butterfly attractors [17], multi-scroll systems [18,19,20], no-equilibrium chaotic systems [21,22,23], etc.

In the recent years, many papers have discussed the modeling of chaotic dynamical systems equipped with special curves having an infinite number of balance points [24].

Table 1. Chaotic systems with special curves of balance points.

Chaotic System	The Curve of Balance Points
Vaidyanathan et al. [25]	Hyperbola
Sambas et al. [26]	Circle
Sambas et al. [27]	Apple
Vaidyanathan et al. [28]	Axe
Mamat et al. [29]	Conch
Sambas et al. [30]	Peanut

lists various dynamical systems exhibiting chaotic behavior with special curves of balance points such as hyperbola [25], circle [26], apple [27], axe [28], conch [29] and peanut [30].

A main contribution of this work is the finding of a new chaotic dynamical system with a butterfly-like curve of balance points. We perform a bifurcation analysis of the new chaotic system and analyze its intrinsic properties such as multi-stability and transient chaos. Multi-stability refers to the coexistence of attractors of the new chaotic system, generated with different initial states of the system under fixed values of the parameters [31,32,33]. Transient chaos phenomenon occurs in a dynamical system when trajectories starting from randomly chosen initial conditions seem chaotic up to a certain time interval and then collapses into a periodic orbit or balance point with the evolution of time [34,35,36].

Furthermore, we carry out field-programmable gate array (FPGA) realization of the new chaotic system with a butterfly-like curve of balance points. Chaotic systems can be implemented in FPGA using numerical analysis and computer algorithms [37,38,39,40]. Realizations of chaotic systems using FPGA are useful for engineering applications such as generators of random numbers [41], secure devices [42], etc.

2. Mathematical Model of a Chaotic System with a Butterfly-Like Curve of Balance Points

The chaotic dynamical systems with special curves of balance points mentioned in Table 1 are particular cases of a general dynamical model stated as follows:

$$\left\{\begin{array}{l}\dot{u}=w\\ \dot{v}=w\phi (u,v,w)\\ \dot{w}=\psi (u,v)\end{array}\right.$$

(1)

It is easy to observe that the balance points of the dynamical system (1) lie on the curve $S$ on the $(u,v)-$plane, which can be defined as follows:

$$S=\left\{(u,v,w)\in R^3|\psi (u,v)=0,w=0\right\}$$

(2)

In this research work, we define the mappings $\phi$ and $\psi$ as given below:

$$\left\{\begin{array}{l}\phi (u,v,w)=-au-b\left|v\right|-cuw\\ \psi (u,v)=u^4+v^4-\left|u\right|\left|v\right|-u^2\end{array}\right.$$

(3)

where $a,b,c$ are real, constant, parameters. We assume that $a,b,c$ are all positive.

Substituting (3) into (1), we obtain a new 3-D mathematical model as given below:

$$\left\{\begin{array}{l}\dot{u}=w\\ \dot{v}=w(-au-b|v|-cuw)\\ \dot{w}=u^4+v^4-\left|u\right|\left|v\right|-u^2\end{array}\right.$$

(4)

The balance points of the system (4) belong to the following curve $S$ on the $(u,v)-$plane.

$$S=\left\{(u,v,w)\in R^3|u^4+v^4-|u\left|\right|v|-u^2=0,w=0\right\}$$

(5)

The closed curve $S$ is symmetrical about both $u$ and $v$ coordinate cases in $R^3.$ In the $(u,v)-$plane, $S$ looks like a butterfly as illustrated in Figure 1.

The new dynamical model (4) exhibits a chaotic attractor if we take $(a,b,c)=(14,0.5,1.4)$ and the initial phase vector $(u(0),v(0),w(0))=(0.05,0.05,0.05)$ as the Lyapunov exponents spectrum can be calculated using Wolf’s algorithm [43] as ${\rho }_1=0.0144,$ ${\rho }_2=0$ and ${\rho }_3=-0.0155$.

The presence of ${\rho }_1>0$ in the Lyapunov exponents spectrum (${\rho }_1,{\rho }_2,{\rho }_3$) establishes the existence of a chaotic attractor of the system (4) [44]. Since ${\rho }_1+{\rho }_2+{\rho }_3=-0.0011<0,$ we also conclude that the system (4) is a dissipative chaotic system [44].

The Kaplan-Yorke dimension [44] of a general n-D dissipative system with its Lyapunov exponents arranged in non-increasing order (${\rho }_1\ge {\rho }_2\ge \dots \ge {\rho }_n$) is defined as follows:

$$K_D=j+\frac{1}{\left|{\rho }_{j+1}\right|}{\displaystyle \sum _{i=1}^j{\rho }_i}$$

(6)

where $j$ is the largest index such that ${\displaystyle \sum _{i=1}^j{\rho }_i}\ge 0$.

Using the Formula (6), we determine the Kaplan-Yorke dimension of the system (4) as follows:

$$K_D=2+\frac{{\rho }_1+{\rho }_2}{\left|{\rho }_3\right|}=2+\frac{0.0144+0}{0.0155}=2.9290$$

(7)

The large value of $K_D$ near three signifies the high complexity of the model (4).

The 2-D signal plots of the chaotic model (4) with butterfly-shaped curve of balance points for $(a,b,c)=(14,0.5,1.4)$ and the initial state $(u(0),v(0),w(0))=(0.05,0.05,0.05)$ are represented in Figure 2, Figure 3 and Figure 4.

3. Dynamic Study of the Chaotic System with Butterfly-Molded Curve of Balance Points

In this section, we shall detail the bifurcation analysis, multi-stability, coexisting attractors and transient chaos of the new chaotic system (4) with butterfly-shaped equilibrium curve.

3.1. Bifurcation Study

In Section 2, it was shown that the system (4) is chaotic for $(a,b,c)=(14,0.5,1.4).$ For the bifurcation analysis, we take $b=0.5$ and $c=1.4$. We consider the initial state as $u(0)=0.05,$ $v(0)=0.05,$ and $w(0)=0.05$. We vary $a$ in the range of $[10,20].$ The bifurcation range of the signal $u$ and the corresponding Lyapunov Exponents (LE) values is shown in Figure 5.

Figure 5 signifies robust chaos in the system (4) in the whole bifurcation region.

Next, we take $a=14$ and $c=1.4$. We consider the initial state as $u(0)=0.05,$ $v(0)=0.05,$ and $w(0)=0.05$. We vary $b$ in the range of $[0,2].$ The bifurcation range of the state variable $y$ and the corresponding Lyapunov exponents is shown in Figure 6.

Figure 6 signifies robust chaos in the system (4) in the whole bifurcation region.

Next, we take $a=14$ and $b=0.5$. We consider the initial state as $u(0)=0.05,$ $v(0)=0.05,$ and $w(0)=0.05$. We vary $c$ in the range of $[1,2]$. The bifurcation range of the state variable $w$ and the corresponding Lyapunov exponents is shown in Figure 7.

Figure 7 also signifies robust chaos in the system (4) in the whole bifurcation region.

Usually when analyzing the dynamics of a nonlinear system with parameter changes, it will encounter stable state, periodic chaos, and period-doubling bifurcation. However, for our system (4), only a chaotic state appears. The system (4) showing robust chaos is very advantageous in chaotic-based engineering applications, which can provide stable chaotic signals without external interference.

3.2. Multi-Stability

Multi-stability for a dynamical system undergoing chaotic behavior means the coexistence of two or more attractors with different initial states but fixed values of the system parameters [31,32,33].

We fix $a=14,$ $b=0.5$ and vary $c$ in the region of $[1,2]$. The initial states are selected as $U_0=(0.05,0.05,0.05)$ and $V_0=(-0.05,-0.05,-0.05)$. The coexisting bifurcation plot for the chaotic system (4) is represented in Figure 8.

The blue color trajectory of the system (4) corresponds to the initial state $U_0,$ while the red color trajectory of the system (4) corresponds to the initial state $V_0.$

When $c=1$ and $c=1.5,$ the chaotic dynamical system (4) exhibits two types of coexisting chaotic attractors in Figure 9a,b, respectively.

3.3. Transient Chaos

Transient chaos means that a dynamical system is in a chaotic state at the beginning, and after time evolution, it eventually tends to a point or a periodic orbit [34,35,36]. When fixing $a=14,$ $b=0.5,$ $c=10,$ and $U(0)=[0.05,0.05,0.05],$ the system (4) exhibits transient chaos as represented in Figure 10. By calculating the Lyapunov exponents of transient chaos of the system (4) in interval [0, 225], we obtain the values of Lyapunov exponents as ${\rho }_1=0.0318,$ ${\rho }_2=0,$ and ${\rho }_3=-0.0363.$

4. FPGA Realization of the Chaotic Dynamical System with a Butterfly-Like Curve of Balance Points

Here, we realize the new chaotic dynamical system (4) using FPGA. Chaotic systems are implemented in FPGA using finite-difference numerical schemes and computer algorithms [37,38,39,40].

The new dynamical system (4) with a closed butterfly-like curve of balance points can be represented in $(x,y,z)-$coordinate system as follows:

$$\left\{\begin{array}{l}\dot{x}=z\\ \dot{y}=z(-ay-b\left|y\right|-cxz)\\ \dot{z}=x^4+y^4-\left|x\right|\left|y\right|-x^2\end{array}\right.$$

(8)

For FPGA implementation of the chaotic system (8), we need to discretize it first.

With the help of a forward Euler integration method [45], we obtain a discrete-time representation of the chaotic system as given below.

$$\left\{\begin{array}{l}x[n+1]=x[n]+hz[n]\\ y[n+1]=y[n]+\left\{z[n](-ay[n]-b\mathrm{abs}(y[n])-cx[n]z[n])\right\}\\ z[n+1]=z[n]+h\left\{x^4[n]+y^4[n]-\mathrm{abs}(x[n])\mathrm{abs}(y[n])-x^2[n]\right\}\end{array}\right.$$

(9)

In the forward Euler scheme (9), $n$ and $n+1$ are the current and next counters respectively, while $h$ is the discretization step size.

Also, $\mathrm{abs}(y)=\left|y\right|$ is the absolute value function, which is defined as follows:

$$\mathrm{abs}(y)=\left\{\begin{array}{ccc}y& \mathrm{if}& y\ge 0\\ -y& \mathrm{if}& y<0\end{array}\right.$$

(10)

The chaotic behavior of the system (8) can be obtained by selecting the coefficient parameters as $(a,b,c)=(14,0.5,1.4)$ and the initial phase as $(x_0,y_0,z_0)=(0.05,0.05,0.05).$

A Python code for implementing the discretized chaotic system (9) is shown in Algorithm 1.

Algorithm 1 Python code for the new discrete-time chaotic system (9)

Input: Initial states x0, y0, and z0. System parameters a, b, and c. Output: xn, yn, and zn    j = 1    step = 500,000    x = x0    y = y0z = z0 while j < step do    xn.insert(j, x + h ∗ z)    yn.insert(j, y + h ∗ (z ∗ (−a ∗ y − b ∗ abs(y) − c ∗ x ∗ z))    zn.insert(j, z + h ∗ (pow(x,4) + pow(y,4) −abs(x) ∗ abs(y) − pow(x,2)))    x = xn[j]    y = yn[j]    z = zn[j]    j = j + 1 end

Figure 11, Figure 12 and Figure 13 represent the simulated orbits for our chaotic oscillator system (7) using the Euler method (Algorithm 1) in the $(x,y),$ $(y,z)$ and $(x,z)$ planes, respectively with $h=0.001$.

Figure 14 shows the top-level hardware design of the discrete-time chaotic system (8).

The chaotic system (8) has three inputs for the initial phase signals $x0,$ $y0$ and $z0$. In addition to the clock and reset signals, the chaotic system (9) has the system parameters $a,$ $b$ and $c$ as the other inputs. As shown in Figure 14, the top-level block diagram contains three sub-systems, which are implemented according to Equations (9) and (10).

Each state variable in any of the chaotic sub-system has 32-bits. It is described using fixed-point representation which consists of 1-bit for the sign, 7-bits for the integer part, and 24-bits for the fractional part (see Figure 15).

The Hardware Description Language (VHDL) code for the top level of the chaotic system (9) is devised as in [37].

The three sub-circuits are implemented using basic mathematical operations (Adder, Subtractor, Multiplier, Multiplexer, and Single Constant Multiplier). The basic blocks of these operations are shown in Figure 16. The VHDL codes for the adder/subtractor, and multiplier are taken as represented in [37].

The Adder, Subtractor, and Multiplier blocks are sequential blocks with both clock (clk) and reset (rst) signals. The Single Constant Multiplier (SCM) performs a multiplication of the input by a constant with the use of only adders, subtractors, and shift operations. The Multiplexer (MUX) block chooses one of the inputs depending on the control signal (ctrl) and passes it to the output. In Figure 16, all inputs and outputs of the hardware blocks have a bus width of 32-bits.

In our hardware design, MUX is used to implement the absolute value function required for the state variables y and z. Hence, Equation (9) should be modified to use the multiplexer (MUX).

In the discrete equation of state variable $y,$ the absolute value function $\left|y\right|$ is used, which is defined in Equation (11).

Equation (10) shows the state variable y after including MUX.

$$\left\{\begin{array}{l}y[n+1]=y[n]-hz[n]\left((a+b)y[n]+cx[n]\cdot z[n]\right)\mathrm{if}y\ge 0\\ y[n+1]=y[n]-hz[n]\left((a-b)y[n]+cx[n]\cdot z[n]\right)\mathrm{if}y<0\end{array}\right.$$

(11)

In the discrete equation of state variable $z,$ two absolute functions are used, viz. $\left|x\right|$ and $\left|y\right|.$ If both state variables have the same sign, then their multiplication is positive. Otherwise, it is negative. To check if the two state variables have the same sign, we need to XOR the sign bit of the state variable $x$$(S_x)$ with the sign bit of the state variable $y$$(S_y),$ i.e., $S_x\mathrm{XOR}{S}_y.$

Equation (12) shows the state variable z after including MUX.

$$\left\{\begin{array}{l}z[n+1]=z[n]+h\left(x^4[n]+y^4[n]-x[n]\cdot y[n]-x^2[n]\right)\mathrm{if}{S}_x\mathrm{XOR}{S}_y=0\\ z[n+1]=z[n]+h\left(x^4[n]+y^4[n]+x[n]\cdot y[n]-x^2[n]\right)\mathrm{if}{S}_x\mathrm{XOR}{S}_y=1\end{array}\right.$$

(12)

To implement the sub-system $X,$ the discrete equation of the signal $x$ in (9) is realized as represented in Figure 17. It uses SCM to multiply the signal $z$ by the constant $h$ and then adds the result to the signal $x$ to get the output $x=x+hz.$

To implement the sub-system $Y,$ the discrete equations of the signal $y$ in (11) are realized as represented in Figure 18. It uses four SCM modules, two multipliers, two adders, one subtractor, and one multiplexer. The SCM modules multiplies the input signals with a constant. The output signals are forwarded to the next modules. The multiplexer receives two signals produced by the adders $(a+b)y+cxz,$$(a-b)y+cxz$ and passes one of them according to the signal $y.$

To implement the sub-system $Z,$ the discrete equations of the signal $z$ in (12) are realized as represented in Figure 19. This sub-system is the most complex one. It includes seven multipliers, three adders, two subtractors, one SCM, and one multiplexer.

The signal $y4$ or $x4$ is produced using three multipliers in two stages. The multiplexer receives two signals produced by an adder $x4+y4-x2+xy$ and a subtractor $x4+y4-x2-xy.$ It passes one of them according to the result of $S_x\mathrm{XOR}{S}_y.$ The output of the multiplexer is multiplied by $h$ and added to the signal $z$ to produce the final output.

Our FPGA realization was carried out with the Zed-Board FPGA prototyping board from Xilinx. The system specifications of the board can be listed as Zynq-7000 all programmable SoC FPGA, 512 MB DDR3Memory, and some other peripherals. We synthesized the generated VHDL codes using Xilinx Vivado design suite and executed on the FPGA prototyping board.

Table 2. FPGA resources use of sub-systems and complete chaotic oscillator system.

Resources	Available	Sub-System X	Sub-System Y	Sub-System Z	Complete System
		Use (%)	Use (%)	Use (%)	Use (%)
FF	106,400	107 (0.10%)	172 (0.16%)	209 (0.2%)	390 (0.37%)
LUT	53,200	137 (0.26%)	877 (1.65%)	434 (0.82%)	1295 (2.43%)
I/O	200	99 (49.5%)	99 (49.5%)	99 (49.5%)	99 (49.5%)
DSP48	220	0 (0%)	8 (3.64%)	20 (9.09%)	28 12.73%)

lists the FPGA resources use for each sub-system and for the complete chaotic oscillator system. As expected, the sub-system X uses less resources as it is less complex than the others. The sub-system Y uses more resources in addition to the Digital Signal Processors (DSPs) as it has more calculations. The sub-system Z is the most complex one. It requires 20 DSPs to perform seven multiplications and other operations.

The complete implementation of discrete-time chaotic system (8) requires 390 Flip-Flops, 1295 Look-up Tables, 99 Input/outputs, and 28 DSPs. The use of the complete system does not exceed 0.37% of the total Flip-Flops, 2.43% of the total Look-up Tables, 12.73% of the DSPs available on the FPGA. The operation frequency of our chaotic oscillator system on the FPGA is 85.39 MHz.

5. Conclusions

Chaotic systems with special curves of balance points are studied with keen interest in the recent years. In this paper, we gave a concise review of three-dimensional chaotic systems equipped with open and closed curves of balance points. We presented a new three-dimensional chaotic system equipped with a closed butterfly-like curve of balance points in the (u,v)-plane. Using bifurcation analysis and Lyapunov exponent diagrams, we found many intrinsic properties such as chaoticity, multi-stability, and transient chaos. Finally, we carried out a realization of the new multi-stable chaotic model using FPGA.