# Dynamical Analysis and FPGA Implementation of a Novel Hyperchaotic System and Its Synchronization Using Adaptive Sliding Mode Control and Genetically Optimized PID Control

Authors: Karthikeyan Rajagopal; Laarem Guessas; Sundarapandian Vaidyanathan; Anitha Karthikeyan; Ashokkumar Srinivasan
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We announce a new 4D hyperchaotic system with four parameters. The dynamic properties of the proposed hyperchaotic system are studied in detail; the Lyapunov exponents, Kaplan-Yorke dimension, bifurcation, and bicoherence contours of the novel hyperchaotic system are derived. Furthermore, control algorithms are designed for the complete synchronization of the identical hyperchaotic systems with unknown parameters using sliding mode controllers and genetically optimized PID controllers. The stabilities of the controllers and parameter update laws are proved using Lyapunov stability theory. Use of the optimized PID controllers ensures less time of convergence and fast synchronization speed. Finally the proposed novel hyperchaotic system is realized in FPGA.


---

## Body

## 1. Introduction

The study of chaotic attractors arising in nonlinear dynamical systems has received much attention as this has potential applications in many branches of science and engineering. Chaotic systems are characterized as nonlinear dynamical systems which are irregular, aperiodic, and unpredictable and have sensitive dependence on the initial conditions [1, 2]. Lorenz discovered a chaotic system while he was studying 3D weather model [3]. Some important paradigms of 3D autonomous chaotic systems are Rössler system [4], Chen system [5], Lü system [6], and so forth. Some novel 3D autonomous chaotic systems are Zhu system [7], Vaidyanathan systems [8–10], Pehlivan system [11], Akgul system [12], Tacha system [13], and so forth.Ahyperchaotic system is mathematically defined as a chaotic system having more than one positive Lyapunov exponent implying that its dynamics are expended in many different directions simultaneously. Thus, a hyperchaotic system has more complex dynamical behaviours than a chaotic system [1, 2]. Rössler was the first to discover a hyperchaotic system [14]. There are many well-known hyperchaotic systems such as hyperchaotic Lorentz system [15], hyperchaotic Chen system [16], hyperchaotic Li system [17], hyperchaotic Vaidyanathan systems [18–20], and so forth.Sliding mode control is one of the best robust control algorithms as discussed in literature [21–23]. Most of the works were concentrated towards active sliding mode control assuming that the system parameters are known. Some authors [24, 25] have derived a methodology to synchronize systems assuming the parameters are unknown. Systems with uncertainty and external disturbance is proposed to realize chaos control [25] in which the bounds of the uncertainty and external disturbance are assumed to be unknown. FPGA implementation of the chaotic generators is more useful for real-time applications such as secure communications and cryptography [26]. A new methodology for implementing integer order integrators in FPGA using Xilinx System Generator toolbox was proposed [27].In this paper, we have derived a novel hyperchaotic system and described the qualitative properties of the system. The Lyapunov exponents of the proposed system are having two positive exponents showing that the system is hyperchaotic. The bifurcation and bicoherence contours of the hyperchaotic system are studied to investigate the parameter dependence and quadratic nonlinearities, respectively.The synchronization of hyperchaotic systems has many applications like secure communications [28, 29], cryptosystems [30, 31], and so forth. In this work, we also derive new results for the adaptive synchronization of the novel hyperchaotic systems with unknown parameters using sliding mode control and PID control [32] optimized by genetic algorithm [33, 34]. The main adaptive control results are established via Lyapunov stability theory [35]. Finally the novel hyperchaotic system is realized in FPGA using the Xilinx (Vivado) System Generator toolbox.


## 2. A Novel Hyperchaotic System

In this work, we propose a novel hyperchaotic system described by(1)x˙=ay-x+w,y˙=cy-xz,z˙=xy-bz+y,w˙=xz+rw,where x,y,z,w are the state variables and a,b,c,r are positive constant parameters.In this work, we show that system (1) ishyperchaoticwhen we take the parameters as(2)a=39,b=3,c=22,r=1.2.For numerical simulations, we take the initial values of the hyperchaotic system (1) as(3)x0=0.2,y0=0.2,z0=0.2,w0=0.2.Figures1–4 show the 3D projections of the novel hyperchaotic system (1) on (x,y,z),(x,y,w),(x,z,w), and (y,z,w) spaces, respectively.Figure 1
3D projection of the novel hyperchaotic system on(x,y,z) space.Figure 2
3D projection of the novel hyperchaotic system on(x,y,w) space.Figure 3
3D projection of the novel hyperchaotic system on(x,z,w) space.Figure 4
3D projection of the novel hyperchaotic system on(y,z,w) space.


## 3. Properties of the Novel Hyperchaotic System

In this section, we discuss the qualitative properties of the novel hyperchaotic system (1) introduced in Section 2. We suppose that the parameter values of system (1) are as in the hyperchaotic case for the system parameters (2).


### 3.1. Dissipativity

In vector notation, the 4D system (1) can be expressed as(4)X˙=fX=f1x,y,z,wf2x,y,z,wf3x,y,z,wf4x,y,z,w,where(5)f1x,y,z,w=ay-x+w,f2x,y,z,w=cy-xz,f3x,y,z,w=xy-bz+y,f4x,y,z,w=xz+rw.LetΩ be any region in R4 with a smooth boundary and also, Φ(t)=ΦtΩ, where Φt is the flow of the vector field f. Furthermore, let V(t) denote the hypervolume of Φ(t).By Liouville’s theorem, we have(6)V˙=∫Φt∇·fdxdydzdw.The divergence of the vector fieldf is easily calculated as(7)∇·f=∂f1∂x+∂f2∂y+∂f3∂z+∂f4∂w=-a+c-b+r=-μ,where μ=a-c+b-r=18.8>0.Substituting (7) into (6), we obtain the first-order differential equation(8)V˙t=-μVt.Integrating (8), we obtain the unique solution as(9)Vt=exp⁡-μtV0∀t≥0.Sinceμ>0, it follows that V(t)→0 exponentially as t→∞. This shows that the novel hyperchaotic system (1) is dissipative. Thus, the system limit sets are ultimately confined into a specific limit set of zero volume, and the asymptotic motion of the novel hyperchaotic system (1) settles onto a strange attractor of the system.


### 3.2. Equilibrium Points

The equilibrium points of the novel hyperchaotic system (1) are obtained by solving(10)f1x,y,z,w=ay-x+w=0,f2x,y,z,w=cy-xz=0,f3x,y,z,w=xy-bz+y=0,f4x,y,z,w=xz+rw=0.The parameter values are taken as in the hyperchaotic case (2).Then the equilibrium points of the novel hyperchaotic system (1) are found as(11)E0=0000,E1=-8.6394-16.303441.5161298.8957,E2=7.639414.416341.5161-264.2990.The Jacobian matrix of the novel hyperchaotic system (1) at any point X∈R4 is found as(12)JX=-aa01-zc-x0yx+1-b0z0xr=-393901-z22-x0yx+1-30z0x1.2.The eigenvalues ofJ0=J(E0) are easily calculated as(13)λ1=1.2,λ2=-3,λ3=22,λ4=-39.This shows thatE0 is a saddle point, which is unstable.The eigenvalues ofJ1=J(E1) are easily calculated as(14)λ1=0.6290,λ2=-13.9770,λ3,4=-2.7260±27.2250i.This shows thatE1 is a saddle-focus, which is unstable.The eigenvalues ofJ2=J(E2) are easily calculated as(15)λ1=0.6292,λ2=-12.5821,λ3,4=-3.4221±26.9621i.This shows thatE2 is a saddle-focus, which is unstable.Hence, all the three equilibrium points of the novel hyperchaotic system (1) are unstable.


### 3.3. Lyapunov Exponents and Kaplan-Yorke Dimension

We take the parameters of the novel 4D system (1) as in the hyperchaotic case (2).Then the Lyapunov exponents of the novel 4D system (1) are numerically found as(16)L1=1.0623,L2=0.2952,L3=0,L4=-20.1063.Since there are two positive Lyapunov exponents in (16), it is clear that the 4D system (1) is hyperchaotic.We note that the sum of the Lyapunov exponents of system (1) is negative. In fact,(17)L1+L2+L3+L4=-18.7488<0.This shows that the novel hyperchaotic system (1) is dissipative.Also, the Kaplan-Yorke dimension of the hyperchaotic system (1) is derived as(18)DKY=3+L1+L2+L3L4=3.0675,which is fractional.


### 3.4. Bicoherence

The bicoherence or the normalized bispectrum is a measure of the amount of phase coupling that occurs in a signal or between two signals. Both bicoherence and bispectrum are used to find the influence of a nonlinear system on the joint probability distribution of the system input. Phase coupling is the estimate of the proportion of energy in every possible pair of frequency componentsfj,fk. Bicoherence analysis is able to detect coherent signals in extremely noisy data, provided that the coherency remains constant for sufficiently long times, since the noise contribution falls off rapidly with increasing number of segments N.The autobispectrum of a chaotic system is given by Pezeshki et al. [36]. They derived the autobispectrum with the Fourier coefficients(19)Bω1,ω2=EAω1Aω2A∗ω1+ω2,where ωn is the radian frequency and A is the Fourier coefficients of the time series. The normalized magnitude spectrum of the bispectrum known as the squared bicoherence is given by(20)bω1,ω2=Bω1,ω22Pω1Pω2Pω1+ω2,where P(ω1) and P(ω2) are the power spectrums at f1 and f2.The motivation to study the bicoherence is twofold. First, the bicoherence can be used to extract information due to deviations from Gaussianity and suppress additive (colored) Gaussian noise. Second, the bicoherence can be used to detect and characterize asymmetric nonlinearity in signals via quadratic phase coupling or identify systems with quadratic nonlinearity. The bicoherence is the third-order spectrum. Whereas the power spectrum is a second-order statistics, formed fromX′f∗Xf, where Xf is the Fourier transform of xt, the bispectrum is a third-order statistics formed from Xfj∗Xfk∗X′fj+fk. The bispectrum is therefore a function of a pair of frequencies fj,fk. It is also a complex-valued function. The (normalized) square amplitude is called the bicoherence (by analogy with the coherence from the cross-spectrum). The bispectrum is calculated by dividing the time series intoM segments of lengthN_seg, calculating their Fourier transforms and biperiodogram, and then averaging over the ensemble. Although the bicoherence is a function of two frequencies the default output of this function is a one dimensional output, the bicoherence refined as a function of only the sum of the two frequencies.Figures5 and 6 depict the contours of the bicoherence of the states x and y for the novel hyperchaotic system (1). Figure 7 shows the bicoherence contours of all the states together. Shades in yellow represent the multifrequency components contributing to the power spectrum. From Figures 5, 6, and 7, the cross-bicoherence is significantly nonzero and nonconstant, indicating a nonlinear relationship between the states. In Figures 5 and 6 (autobicoherence) the yellow shades indicate that the nonlinear relationships of the states x and y are not narrow-band processes. The yellow shades and nonsharpness of the peaks as well as the presence of structure around the origin in Figure 7 (cross-bicoherence) indicate that the nonlinearity between the states x,y,z,w is not of the quadratic nonlinearity and hence may be because of nonlinearity of higher dimensions. The most two dominant frequencies (f1,f2) are taken for deriving the contour of bicoherence. The sampling frequency (fs) is taken as the reference frequency. Direct FFT is used to derive the power spectrum for individual frequencies and Hankel operator is used as the frequency mask. Hanning window is used as the FIR filter to separate the frequencies.Figure 5
Contour of bicoherence of statex for the novel hyperchaotic system.Figure 6
Contour of bicoherence of statey for the novel hyperchaotic system.Figure 7
Contour of bicoherence of all the states of the novel hyperchaotic system.For the bicoherence plots we used the pseudocode [37] given below.Analysis of the Input Data n = n e l e m e n t s x if (n<100) then Not enough data’ Is the segment length supplied? if(nelements(nseg)=0) then nseg=n/20 Is the sampling perioddt supplied? if(nelements(dt)=0) then dt=1.0Break the Time Series into M Segments of N Length Each m = n / n s e g if (m<1) then S i z e D A T A < n s e g n c u t = m ∗ n s e g d a t a = r e f o r m x 0 : n c u t - 1 , n s e g , m n p = + v e   f r e q u e n c i e s n f = n s e g / 2 Frequency resolution d f = 1.0 / d t ∗ n s e g   d f = 1.0 / ( d t ∗ n _ s e g ) Frequency array f = f i n d g e n n f + 1 ∗ d fCalculate the Periodogram for Each Segment Using FFT Command Calculate mean value (DC component) for each segment m e a n  x = t o t a l d a t a , 1 / f l o a t n s e g Subtract mean (DC component) d a t a = t e m p o r a r y d a t a - t r a n s p o s e r e b i n m e a n  x , m , n s e g , / s a m p l e Calculate the Fourier transform of each segment (row) this givesDFTfj,ti data = fft(data, 1, dimension = 1,/overwrite) Extract only positive frequenciesfj:j=0,1,2,…,nf D F T = d a t a 0 : n f , ∗ Calculate the index matrixj,k=j+k i n d x = i n d g e n n f + 1 j = r e b i n i n d x , n f + 1 , n f + 1 , / s a m p l e k = t r a n s p o s e j j k = j + kPrepare Arrays for Storing Components of Bicoherence b i c o h = | U | 2 / ( V ∗ W ) whereU=Xj∗Xk∗X′j+k is the complex value V = | X j ∗ X k | 2, W=|Xj+k|2 and Poij=|Xj|2 are the real values U = m a k e a r r a y n f + 1 , n f + 1 , / c o m p l e x , v a l u e = 0.0 V = m a k e a r r a y n f + 1 , n f + 1 , / f l o a t , v a l u e = 0.0 W = m a k e a r r a y n f + 1 , n f + 1 , / f l o a t , v a l u e = 0.0 P o i j = m a k e a r r a y n f + 1 , / f l o a t , v a l u e = 0.0 P o i n = 0.0Normalize So That | D F T | 2 Is in [ r m s / m e a n ] 2 H z - 1 Units n o r m = s q r t ( 2.0 ∗ d t / f l o a t ( n s e g ) / m e a n  x [ i ] 2 ) X j = t e m p o r a r y X j ∗ n o r m Calculate the matrixXj∗Xk X j X k = X j # X j Calculate the matrixXj+k X j k = X j j kCalculate the Triple Term X j ∗ X k ∗ X j + k X j X k X j k = X j X k ∗ c o n j X j kSum over the Ensembles of Each Segment U = t e m p o r a r y ( U ) + X j X k X j k V = t e m p o r a r y ( V ) + a b s ( X j X k ) 2 W = t e m p o r a r y ( W ) + a b s ( X j k ) 2Define Power Spectrum P o i j = t e m p o r a r y ( P o i j ) + a b s ( X j ) 2 and its Poisson noise level P o i n = P o i n + 2.0 / m e a n  x iNormalize the Sums to Get Means U = t e m p o r a r y U / f l o a t m V = t e m p o r a r y V / f l o a t m W = t e m p o r a r y W / f l o a t m P o i j = t e m p o r a r y P o i j / f l o a t m P o i n = P o i n / f l o a t mCalculate the Bicoherence from | U | 2 / ( V ∗ C ) U = a b s ( U ) 2 V W = V ∗ WMasking the Variables to Avoid/0 Errors mask = where (VW≠0, count) Prepare array for 2-dimensional bicoherence b i c o h = m a k e a r r a y n f + 1 , n f + 1 , / f l o a t , v a l u e = 0 normalizebispectrum2 to get bicoherence i f   c o u n t > 0   t h e n   b e g i n b i c o h m a s k = U m a s k / V W m a s k Plot the results.


### 3.5. Bifurcation

In order to understand the dynamical behaviour of the MNECS, the bifurcation plots are derived for three cases as follows.Case 1.
Fixb=3,c=22 and vary a between [2540].Case 2.
Fixa=39,c=22 and vary b between [03.5].Case 3.
Fixa=39,b=3 and vary c between [2025].Figures8(a), 8(b), and 8(c) show the bifurcation plots of the system for the parameters a,b,c, respectively. As it is observed from Figures 8(a), 8(b), and 8(c), for the values of 38≤a≤39.2,2.8≤b≤3.1,20.8≤c≤22.3, the bifurcation plots show denser points confirming the existence of two stable positive Lyapunov exponents (hyperchaos) showing the parameter dependence of system (1).Figure 8
Bifurcation plots versus (a)a; (b) b; (c) c.
(a) (b) (c)


## 3.1. Dissipativity

In vector notation, the 4D system (1) can be expressed as(4)X˙=fX=f1x,y,z,wf2x,y,z,wf3x,y,z,wf4x,y,z,w,where(5)f1x,y,z,w=ay-x+w,f2x,y,z,w=cy-xz,f3x,y,z,w=xy-bz+y,f4x,y,z,w=xz+rw.LetΩ be any region in R4 with a smooth boundary and also, Φ(t)=ΦtΩ, where Φt is the flow of the vector field f. Furthermore, let V(t) denote the hypervolume of Φ(t).By Liouville’s theorem, we have(6)V˙=∫Φt∇·fdxdydzdw.The divergence of the vector fieldf is easily calculated as(7)∇·f=∂f1∂x+∂f2∂y+∂f3∂z+∂f4∂w=-a+c-b+r=-μ,where μ=a-c+b-r=18.8>0.Substituting (7) into (6), we obtain the first-order differential equation(8)V˙t=-μVt.Integrating (8), we obtain the unique solution as(9)Vt=exp⁡-μtV0∀t≥0.Sinceμ>0, it follows that V(t)→0 exponentially as t→∞. This shows that the novel hyperchaotic system (1) is dissipative. Thus, the system limit sets are ultimately confined into a specific limit set of zero volume, and the asymptotic motion of the novel hyperchaotic system (1) settles onto a strange attractor of the system.


## 3.2. Equilibrium Points

The equilibrium points of the novel hyperchaotic system (1) are obtained by solving(10)f1x,y,z,w=ay-x+w=0,f2x,y,z,w=cy-xz=0,f3x,y,z,w=xy-bz+y=0,f4x,y,z,w=xz+rw=0.The parameter values are taken as in the hyperchaotic case (2).Then the equilibrium points of the novel hyperchaotic system (1) are found as(11)E0=0000,E1=-8.6394-16.303441.5161298.8957,E2=7.639414.416341.5161-264.2990.The Jacobian matrix of the novel hyperchaotic system (1) at any point X∈R4 is found as(12)JX=-aa01-zc-x0yx+1-b0z0xr=-393901-z22-x0yx+1-30z0x1.2.The eigenvalues ofJ0=J(E0) are easily calculated as(13)λ1=1.2,λ2=-3,λ3=22,λ4=-39.This shows thatE0 is a saddle point, which is unstable.The eigenvalues ofJ1=J(E1) are easily calculated as(14)λ1=0.6290,λ2=-13.9770,λ3,4=-2.7260±27.2250i.This shows thatE1 is a saddle-focus, which is unstable.The eigenvalues ofJ2=J(E2) are easily calculated as(15)λ1=0.6292,λ2=-12.5821,λ3,4=-3.4221±26.9621i.This shows thatE2 is a saddle-focus, which is unstable.Hence, all the three equilibrium points of the novel hyperchaotic system (1) are unstable.


## 3.3. Lyapunov Exponents and Kaplan-Yorke Dimension

We take the parameters of the novel 4D system (1) as in the hyperchaotic case (2).Then the Lyapunov exponents of the novel 4D system (1) are numerically found as(16)L1=1.0623,L2=0.2952,L3=0,L4=-20.1063.Since there are two positive Lyapunov exponents in (16), it is clear that the 4D system (1) is hyperchaotic.We note that the sum of the Lyapunov exponents of system (1) is negative. In fact,(17)L1+L2+L3+L4=-18.7488<0.This shows that the novel hyperchaotic system (1) is dissipative.Also, the Kaplan-Yorke dimension of the hyperchaotic system (1) is derived as(18)DKY=3+L1+L2+L3L4=3.0675,which is fractional.


## 3.4. Bicoherence

The bicoherence or the normalized bispectrum is a measure of the amount of phase coupling that occurs in a signal or between two signals. Both bicoherence and bispectrum are used to find the influence of a nonlinear system on the joint probability distribution of the system input. Phase coupling is the estimate of the proportion of energy in every possible pair of frequency componentsfj,fk. Bicoherence analysis is able to detect coherent signals in extremely noisy data, provided that the coherency remains constant for sufficiently long times, since the noise contribution falls off rapidly with increasing number of segments N.The autobispectrum of a chaotic system is given by Pezeshki et al. [36]. They derived the autobispectrum with the Fourier coefficients(19)Bω1,ω2=EAω1Aω2A∗ω1+ω2,where ωn is the radian frequency and A is the Fourier coefficients of the time series. The normalized magnitude spectrum of the bispectrum known as the squared bicoherence is given by(20)bω1,ω2=Bω1,ω22Pω1Pω2Pω1+ω2,where P(ω1) and P(ω2) are the power spectrums at f1 and f2.The motivation to study the bicoherence is twofold. First, the bicoherence can be used to extract information due to deviations from Gaussianity and suppress additive (colored) Gaussian noise. Second, the bicoherence can be used to detect and characterize asymmetric nonlinearity in signals via quadratic phase coupling or identify systems with quadratic nonlinearity. The bicoherence is the third-order spectrum. Whereas the power spectrum is a second-order statistics, formed fromX′f∗Xf, where Xf is the Fourier transform of xt, the bispectrum is a third-order statistics formed from Xfj∗Xfk∗X′fj+fk. The bispectrum is therefore a function of a pair of frequencies fj,fk. It is also a complex-valued function. The (normalized) square amplitude is called the bicoherence (by analogy with the coherence from the cross-spectrum). The bispectrum is calculated by dividing the time series intoM segments of lengthN_seg, calculating their Fourier transforms and biperiodogram, and then averaging over the ensemble. Although the bicoherence is a function of two frequencies the default output of this function is a one dimensional output, the bicoherence refined as a function of only the sum of the two frequencies.Figures5 and 6 depict the contours of the bicoherence of the states x and y for the novel hyperchaotic system (1). Figure 7 shows the bicoherence contours of all the states together. Shades in yellow represent the multifrequency components contributing to the power spectrum. From Figures 5, 6, and 7, the cross-bicoherence is significantly nonzero and nonconstant, indicating a nonlinear relationship between the states. In Figures 5 and 6 (autobicoherence) the yellow shades indicate that the nonlinear relationships of the states x and y are not narrow-band processes. The yellow shades and nonsharpness of the peaks as well as the presence of structure around the origin in Figure 7 (cross-bicoherence) indicate that the nonlinearity between the states x,y,z,w is not of the quadratic nonlinearity and hence may be because of nonlinearity of higher dimensions. The most two dominant frequencies (f1,f2) are taken for deriving the contour of bicoherence. The sampling frequency (fs) is taken as the reference frequency. Direct FFT is used to derive the power spectrum for individual frequencies and Hankel operator is used as the frequency mask. Hanning window is used as the FIR filter to separate the frequencies.Figure 5
Contour of bicoherence of statex for the novel hyperchaotic system.Figure 6
Contour of bicoherence of statey for the novel hyperchaotic system.Figure 7
Contour of bicoherence of all the states of the novel hyperchaotic system.For the bicoherence plots we used the pseudocode [37] given below.Analysis of the Input Data n = n e l e m e n t s x if (n<100) then Not enough data’ Is the segment length supplied? if(nelements(nseg)=0) then nseg=n/20 Is the sampling perioddt supplied? if(nelements(dt)=0) then dt=1.0Break the Time Series into M Segments of N Length Each m = n / n s e g if (m<1) then S i z e D A T A < n s e g n c u t = m ∗ n s e g d a t a = r e f o r m x 0 : n c u t - 1 , n s e g , m n p = + v e   f r e q u e n c i e s n f = n s e g / 2 Frequency resolution d f = 1.0 / d t ∗ n s e g   d f = 1.0 / ( d t ∗ n _ s e g ) Frequency array f = f i n d g e n n f + 1 ∗ d fCalculate the Periodogram for Each Segment Using FFT Command Calculate mean value (DC component) for each segment m e a n  x = t o t a l d a t a , 1 / f l o a t n s e g Subtract mean (DC component) d a t a = t e m p o r a r y d a t a - t r a n s p o s e r e b i n m e a n  x , m , n s e g , / s a m p l e Calculate the Fourier transform of each segment (row) this givesDFTfj,ti data = fft(data, 1, dimension = 1,/overwrite) Extract only positive frequenciesfj:j=0,1,2,…,nf D F T = d a t a 0 : n f , ∗ Calculate the index matrixj,k=j+k i n d x = i n d g e n n f + 1 j = r e b i n i n d x , n f + 1 , n f + 1 , / s a m p l e k = t r a n s p o s e j j k = j + kPrepare Arrays for Storing Components of Bicoherence b i c o h = | U | 2 / ( V ∗ W ) whereU=Xj∗Xk∗X′j+k is the complex value V = | X j ∗ X k | 2, W=|Xj+k|2 and Poij=|Xj|2 are the real values U = m a k e a r r a y n f + 1 , n f + 1 , / c o m p l e x , v a l u e = 0.0 V = m a k e a r r a y n f + 1 , n f + 1 , / f l o a t , v a l u e = 0.0 W = m a k e a r r a y n f + 1 , n f + 1 , / f l o a t , v a l u e = 0.0 P o i j = m a k e a r r a y n f + 1 , / f l o a t , v a l u e = 0.0 P o i n = 0.0Normalize So That | D F T | 2 Is in [ r m s / m e a n ] 2 H z - 1 Units n o r m = s q r t ( 2.0 ∗ d t / f l o a t ( n s e g ) / m e a n  x [ i ] 2 ) X j = t e m p o r a r y X j ∗ n o r m Calculate the matrixXj∗Xk X j X k = X j # X j Calculate the matrixXj+k X j k = X j j kCalculate the Triple Term X j ∗ X k ∗ X j + k X j X k X j k = X j X k ∗ c o n j X j kSum over the Ensembles of Each Segment U = t e m p o r a r y ( U ) + X j X k X j k V = t e m p o r a r y ( V ) + a b s ( X j X k ) 2 W = t e m p o r a r y ( W ) + a b s ( X j k ) 2Define Power Spectrum P o i j = t e m p o r a r y ( P o i j ) + a b s ( X j ) 2 and its Poisson noise level P o i n = P o i n + 2.0 / m e a n  x iNormalize the Sums to Get Means U = t e m p o r a r y U / f l o a t m V = t e m p o r a r y V / f l o a t m W = t e m p o r a r y W / f l o a t m P o i j = t e m p o r a r y P o i j / f l o a t m P o i n = P o i n / f l o a t mCalculate the Bicoherence from | U | 2 / ( V ∗ C ) U = a b s ( U ) 2 V W = V ∗ WMasking the Variables to Avoid/0 Errors mask = where (VW≠0, count) Prepare array for 2-dimensional bicoherence b i c o h = m a k e a r r a y n f + 1 , n f + 1 , / f l o a t , v a l u e = 0 normalizebispectrum2 to get bicoherence i f   c o u n t > 0   t h e n   b e g i n b i c o h m a s k = U m a s k / V W m a s k Plot the results.


## 3.5. Bifurcation

In order to understand the dynamical behaviour of the MNECS, the bifurcation plots are derived for three cases as follows.Case 1.
Fixb=3,c=22 and vary a between [2540].Case 2.
Fixa=39,c=22 and vary b between [03.5].Case 3.
Fixa=39,b=3 and vary c between [2025].Figures8(a), 8(b), and 8(c) show the bifurcation plots of the system for the parameters a,b,c, respectively. As it is observed from Figures 8(a), 8(b), and 8(c), for the values of 38≤a≤39.2,2.8≤b≤3.1,20.8≤c≤22.3, the bifurcation plots show denser points confirming the existence of two stable positive Lyapunov exponents (hyperchaos) showing the parameter dependence of system (1).Figure 8
Bifurcation plots versus (a)a; (b) b; (c) c.
(a) (b) (c)


## 4. Adaptive Synchronization of the Novel Hyperchaotic Systems Using Sliding Mode Control

In this section, we deploy the adaptive sliding mode control (ASMC) method to derive new results for the complete synchronization of the novel hyperchaotic systems using sliding mode control.The master system is given as(21)x˙m=aym-xm+wm,y˙m=cym-xmzm,z˙m=xmym-bzm+ym,w˙m=xmzm+rwm.The slave systems with the adaptive sliding mode controllers (ui) are defined as(22)x˙s=a^ys-xs+ws+ux,y˙s=c^ys-xszs+uy,z˙s=xsys-b^zs+ys+uz,w˙s=xszs+r^ws+uw.The parameters of the slave system are assumed to be unknown and hence replaced by their respective parameter estimates a^,b^,c^,r^.For the complete synchronization of the master and slave systems, the synchronization errorslimt→∞ei=0, where i=x,y,z,w and ei is defined as (23)ei=is-im.The error dynamics are derived by differentiating the synchronization errors(24)e˙i=i˙s-i˙m.Using the master system (21) and slave system (22) in the error dynamics (24),(25)e˙x=a^ys-xs+ws-aym-xm-wm+ux,e˙y=c^ys-xszs-cym+xmzm+uy,e˙z=xsys-b^zs+ys-xmym+bzm-ym+uz,e˙w=xszs+r^ws-xmzm-rwm+uw.


### 4.1. Control Design

The main control objective of this paper is to design adaptive sliding mode controllers (ux,uy,uz,uw) such that the synchronization errors limt→∞ei=0. Sliding mode control methodology is a simple approach to robust control and good at dealing with dynamic uncertainty [25, 38]. The control design procedure consists of two steps: first constructing a sliding surface which presents the desired dynamics [39] and second selecting a switching control law so as to verify sliding condition.The integral sliding surface [38] can be defined as(26)si=ei+ki∫eiτdτ,where i=x,y,z,w and ki>0.The time derivate of the sliding surface is given by(27)s˙i=e˙i+kiei.For any initial conditions, the problem of synchronization of master and slave systems is equivalent to that of the synchronization errors ei remaining on the surface si for all t>0. When the system operates in the sliding surface, it satisfies s=0 and s˙=0.The parameter estimation errors are given by(28)ea=a^-a,eb=b^-b,ec=c^-c,er=r^-r.Let the adaptive sliding mode controllers be defined as(29)ux=-a^ey-ws+wm-γxsgnsx-ρxsx-kxex,uy=-c^ey+xszs-xmzm-γysgnsy-ρysy-kyey,uz=-xsys+b^ez+xmym-ys+ym-γzsgnsz-ρzsz-kzez,uw=-xszs-r^ew+xmzm-γwsgnsw-ρwsw-kwew,where ρi>0,γi>0 are the sliding surface gains and ki is the controller gain for i=x,y,z,w.Let us define the parameter adaptive laws as(30)a^˙=-sxym-xm,b^˙=szzm,c^˙=-syym,r^˙=-swwm,where a^˙,b^˙,c^˙,r^˙ are the dynamics of the parameter estimates a^,b^,c^,r^.


### 4.2. Stability Analysis of the Controller

Theorem 1.
Considering that adaptive sliding mode control input law in (5) is used to control error system in (6) with update laws of parameters in (7), then error system in (6) is globally asymptotically stable.Proof.
To check the stability of the controlled system, let us consider the following Lyapunov candidate function:(31)V=12sx2+sy2+sz2+sw2+ea2+eb2+ec2+er2.The derivative of the Lyapunov function is given by (32)V˙=sxs˙x+sys˙y+szs˙z+sws˙w+eae˙a+ebe˙b+ece˙c+ere˙r.Using (27) and (25) in (32),(33)V˙=sxa^ys-a^xs-aym+axm+ws-wm+ux+kxex+syc^ys-xszs-cym+xmzm+uy+kyey+szxsys-b^zs+ys-xmym+bzm-ym+uz+kzez+swxszs+r^ws-xmzm-rwm+uw+kwew+eaa^˙+ebb^˙+ecc^˙+err^˙.After some mathematical manipulations we arrive at(34)V˙=sxa^ey+kxex+ws-wm+ux+syc^ey-xszs+xmzm+uy+kyey+szxsys-b^ez-xmym+ys-ym+uz+kzez+swxszs+r^ew-xmzm+uw+kwew+easxym-xm+a^˙+ebszzm+b^˙+ecsyym+c^˙+erswwm+r^˙,where si=-γisgn(si)-ρisi is the sliding reaching law [38], i=x,y,z,w, and ρi>0,γi>0 are the sliding surface gains.
Introducing the parameter update laws (30) and the adaptive sliding mode controllers (29) in (34), we simplify the Lyapunov first derivative as (35)V˙=-γxsx-γysy-γzsz-γwsw-ρxsx2-ρysy2-ρzsz2-ρwsw2.As can be seen from (35), Lyapunov first derivative is a negative semidefinite as ρi>0 and γi>0, hence, by Barbalat’s lemma in Lyapunov stability theory [36, 40], it follows that the synchronization errors ex(t),ey(t),ez(t),ew(t) converge exponentially to zero as t→∞.
For numerical simulations, we take the initial conditions of the master system (27) as(36)xm0=1.2,ym0=6.8,zm0=3.9,wm0=5.4.We take the initial conditions of the slave system (28) as(37)xs0=4.7,ys0=2.6,zs0=1.4,ws0=1.9.We take the initial conditions of the parameter estimates as(38)a^0=5.2,b^0=3.1,c^0=1.2,r^0=6.4.Figure 9 depicts the time-history of the synchronization errors and Figure 10 shows the synchronized states. Figure 11(a) shows the time-history of the states. Figure 11(b) shows the time-history of adaptive sliding mode (ASM) controllers.Figure 9
Time-history of the synchronization errors.Figure 10
Time-history of the parameter estimates.Figure 11
(a) Time-history of the synchronized states. (b) Time-history of the ASM controllers.
(a) (b)


## 4.1. Control Design

The main control objective of this paper is to design adaptive sliding mode controllers (ux,uy,uz,uw) such that the synchronization errors limt→∞ei=0. Sliding mode control methodology is a simple approach to robust control and good at dealing with dynamic uncertainty [25, 38]. The control design procedure consists of two steps: first constructing a sliding surface which presents the desired dynamics [39] and second selecting a switching control law so as to verify sliding condition.The integral sliding surface [38] can be defined as(26)si=ei+ki∫eiτdτ,where i=x,y,z,w and ki>0.The time derivate of the sliding surface is given by(27)s˙i=e˙i+kiei.For any initial conditions, the problem of synchronization of master and slave systems is equivalent to that of the synchronization errors ei remaining on the surface si for all t>0. When the system operates in the sliding surface, it satisfies s=0 and s˙=0.The parameter estimation errors are given by(28)ea=a^-a,eb=b^-b,ec=c^-c,er=r^-r.Let the adaptive sliding mode controllers be defined as(29)ux=-a^ey-ws+wm-γxsgnsx-ρxsx-kxex,uy=-c^ey+xszs-xmzm-γysgnsy-ρysy-kyey,uz=-xsys+b^ez+xmym-ys+ym-γzsgnsz-ρzsz-kzez,uw=-xszs-r^ew+xmzm-γwsgnsw-ρwsw-kwew,where ρi>0,γi>0 are the sliding surface gains and ki is the controller gain for i=x,y,z,w.Let us define the parameter adaptive laws as(30)a^˙=-sxym-xm,b^˙=szzm,c^˙=-syym,r^˙=-swwm,where a^˙,b^˙,c^˙,r^˙ are the dynamics of the parameter estimates a^,b^,c^,r^.


## 4.2. Stability Analysis of the Controller

Theorem 1.
Considering that adaptive sliding mode control input law in (5) is used to control error system in (6) with update laws of parameters in (7), then error system in (6) is globally asymptotically stable.Proof.
To check the stability of the controlled system, let us consider the following Lyapunov candidate function:(31)V=12sx2+sy2+sz2+sw2+ea2+eb2+ec2+er2.The derivative of the Lyapunov function is given by (32)V˙=sxs˙x+sys˙y+szs˙z+sws˙w+eae˙a+ebe˙b+ece˙c+ere˙r.Using (27) and (25) in (32),(33)V˙=sxa^ys-a^xs-aym+axm+ws-wm+ux+kxex+syc^ys-xszs-cym+xmzm+uy+kyey+szxsys-b^zs+ys-xmym+bzm-ym+uz+kzez+swxszs+r^ws-xmzm-rwm+uw+kwew+eaa^˙+ebb^˙+ecc^˙+err^˙.After some mathematical manipulations we arrive at(34)V˙=sxa^ey+kxex+ws-wm+ux+syc^ey-xszs+xmzm+uy+kyey+szxsys-b^ez-xmym+ys-ym+uz+kzez+swxszs+r^ew-xmzm+uw+kwew+easxym-xm+a^˙+ebszzm+b^˙+ecsyym+c^˙+erswwm+r^˙,where si=-γisgn(si)-ρisi is the sliding reaching law [38], i=x,y,z,w, and ρi>0,γi>0 are the sliding surface gains.
Introducing the parameter update laws (30) and the adaptive sliding mode controllers (29) in (34), we simplify the Lyapunov first derivative as (35)V˙=-γxsx-γysy-γzsz-γwsw-ρxsx2-ρysy2-ρzsz2-ρwsw2.As can be seen from (35), Lyapunov first derivative is a negative semidefinite as ρi>0 and γi>0, hence, by Barbalat’s lemma in Lyapunov stability theory [36, 40], it follows that the synchronization errors ex(t),ey(t),ez(t),ew(t) converge exponentially to zero as t→∞.
For numerical simulations, we take the initial conditions of the master system (27) as(36)xm0=1.2,ym0=6.8,zm0=3.9,wm0=5.4.We take the initial conditions of the slave system (28) as(37)xs0=4.7,ys0=2.6,zs0=1.4,ws0=1.9.We take the initial conditions of the parameter estimates as(38)a^0=5.2,b^0=3.1,c^0=1.2,r^0=6.4.Figure 9 depicts the time-history of the synchronization errors and Figure 10 shows the synchronized states. Figure 11(a) shows the time-history of the states. Figure 11(b) shows the time-history of adaptive sliding mode (ASM) controllers.Figure 9
Time-history of the synchronization errors.Figure 10
Time-history of the parameter estimates.Figure 11
(a) Time-history of the synchronized states. (b) Time-history of the ASM controllers.
(a) (b)


## 5. PID Controller Based on Genetic Algorithms (GAPID)

Many advanced and intelligent control algorithms have been developed, however, the PID-type controller [16] remains the most popular in industry for the simplicity of control law, few tuning parameters, less time convergence and fast synchronization speed. But finding appropriate parameters values for the PID controller is still a difficult task, so in practice control engineers still often use trial and error for the tuning process. So in practice control engineers still often use trial and error for the tuning process.PID control consists of three types of control, proportional, integral, and derivative control; they can be implemented in many forms but they are mostly used in feedback loops.In this paper we use four PID controllers (ux,uy,uz,uw) as given in (39) to synchronize the master system states (xm,ym,zm,wm) with slave system states (xs,ys,zs,ws).(39)ui=KPei+KI∫0teiτdτ+KDdeidt,where ui is the PID controller for i=x,y,z,w,ei is the error signal, and KP,KI,KD are the proportional, integral, and derivative gains of the PID controller for the synchronization of the master novel hyperchaos system (21) and the slave system (22). Tuning of PID controllers involves the selection of KP,KI,KD gains for better control performance which is defined with reference to the required performance index.In this work we will investigate the effectiveness of the PID controller gains when optimized using genetic algorithms (GA). In particular, many PID control based on GA schemes have been successfully applied to the control and the synchronization of chaotic systems [41], GA’s [33, 34] are a stochastic global search method that mimics the process of natural evolution. It starts with no knowledge of the correct solution and depends entirely on responses from its environment and evolution operators (i.e., reproduction, crossover and mutation) to arrive at the best solution. By starting at several independent points (population) and searching in parallel, the algorithm avoids local minima and converging to suboptimal solutions. A genetic algorithm is typically initialized with a random population consisting of between 20–100 individuals. This population is usually represented by a real-valued number or a binary string called a chromosome. How well an individual performs a task is measured and assessed by the objective function. The objective function assigns each individual a corresponding number called its fitness. The fitness of each chromosome is assessed and a survival of the fittest strategy is applied. The magnitude of the error is usually used to assess the fitness of each chromosome. Writing an objective function is the most difficult part of creating a genetic algorithm. An objective function could be created to find a PID controller that gives the smallest overshoot, fastest rise time or quickest settling time, however, in order to combine all of these objectives it was decided to design an objective function that will minimize the error of the controlled system. Each chromosome in the population is passed into the objective function one at a time. The chromosome is then evaluated and assigned a number to represent its fitness, the bigger its number the better its fitness. The genetic algorithm uses the chromosome’s fitness value to create a new population consisting of the fittest members. There are two important performance indices ISE (Integrated Squared Error) and IAE (Integrated Absolute Error) given as(40)ISE=∫0∞ei2τdτ,IAE=∫0∞eiτdτ.In this paper we use IAE as the objective function and the fitness functions is given as(41)fitness=1∑abseiwherei=x,y,z,w.Using the fitness function (41) and objective function (40) in the optimization toolbox of Matlab with the following options, we derive the optimal Kp,Ki,Kd values for each of the controllers (ux,uy,uz,uw).Variable bounds matrix =-10,10.Population size = 80; GA. Generally the bigger the population size, the better the final approximation.Number of generations =100.Selection function = stochastic uniform.Crossover fraction =0.8.Mutation function = Gaussian.Stopping criteria = error performance criterion.Length of the chromosome =12, decimal codage.With the same parameter values used in (2) and the same the initial conditions in (36) and (37), Table 1 shows the PID gain values after running the genetic algorithm (GA) solver in the optimization tool with the options cited above, we get the best solutions tracked over generations for the complete synchronization of the novel hyperchaotic system:Table 1
PID controller gain values optimized with GA.
PID controller K P K I K D u x 7.9839 −22.1195 17.4396 u y −75.1076 −0.5207 −2.2663 u z −15.8710 −65.0964 14.4902 u w −63.6645 −15.3596 −6.7482Figure12 depicts the time-history of the synchronization errors using genetically optimized PID controllers, Figure 13 shows the synchronized states, and Figure 14 shows the time-history of the PIDs controllers.Figure 12
Time-history of the synchronization errors.Figure 13
State synchronization(xm,xs) (control in action at t=0.07s).Figure 14
Time-history of the PIDs controller.As seen from Figures8 and 11, in the ASMC method of synchronization the errors converges at t=0.3s whereas the GAPID controllers based synchronization the errors converges at t=0.15s. The states synchronization Figures 9 and 12 and the controllers convergence plots shown in Figures 11 and 14 clearly indicates that GAPID controllers are efficient compared to the ASMC controllers. As can be seen from Figures 11 and 14, ASM controllers converging time (t=200mS) is slower than the GAPID controllers (t=160mS). Also the ASM controllers shoot to peak values of [-10001000] which when realized in electronic circuits need special attenuators to suppress the oscillations whereas in GAPID controllers the shoot values are reduced comparatively.


## 6. FPGA Implementation of the Novel Hyperchaotic System

In this section we discuss about the implementation of the novel hyperchaotic system in FPGA using the Xilinx (Vivado) System Generator toolbox in Simulink. Firstly we configure the available built in blocks of the System Generator toolbox. The Add/Subblocks are configured with zero latency and 32/16 bit fixed point settings. The output of the block is configured to rounded quantization in order to reduce the bit latency. For the multiplier block a latency of 3 is configured and the other settings are the same as in Add/Subblock. Next we will have to design the integer order integrator which is not a readily available block in the System Generator. Hence we implement the integrators using the mathematical relationdxi/dt=limh→0[xi(n+1)-xi(n)]/h and the value of h is taken as 0.001 and the initial conditions are fed in to the forward register. Figure 15 shows the Xilinx RTL schematics of the system using Virtex 7-XC7 chips and Figures 16, 17, 18, and 19 shows the 3D phase portraits of the novel hyperchaotic system using Xilinx System Generator. Comparing Figures 1, 2, 3, and 4 with Figures 16, 17, 18, and 19 one can clearly see that the FPGA implemented novel hyperchaotic system exhibits the same phase portraits for the initial conditions [0.20.20.20.2].Figure 15
Xilinx schematics of the system implemented in Virtex 7 (Device-xc7vx690tffg1761).Figure 16
3D projection of the novel hyperchaotic system on(x,y,z) space.Figure 17
3D projection of the novel hyperchaotic system on(x,y,w) space.Figure 18
3D projection of the novel hyperchaotic system on(x,z,w) space.Figure 19
3D projection of the novel hyperchaotic system on(y,z,w) space.


## 7. Conclusions

In this research work, we first described a novel hyperchaotic system and discussed its qualitative properties. Lyapunov exponents are derived to prove that the system is hyperchaotic. Bicoherence plots are drawn to prove the multifrequency terms dependency of the novel hyperchaotic systems. Next, identical hyperchaotic systems are synchronized using adaptive sliding mode control assuming that the system parameters are unknown and also genetically optimized PID controllers. The complete synchronization results derived in this paper have been proved using the Lyapunov stability theory. A comparison between the ASM controllers and GAPID controllers are discussed with various plots. Finally the 4D novel hyperchaotic system is implemented in FPGA implementation to show that the system is hardware realizable.



---

*Source: 7307452-2017-02-19.xml*