# Analysis, Synchronization, and Robotic Application of a Modified Hyperjerk Chaotic System

Authors: Lazaros Moysis; Eleftherios Petavratzis; Muhammad Marwan; Christos Volos; Hector Nistazakis; Salman Ahmad
Journal: Complexity
 (2020)

---

## Abstract

In this work, a novel hyperjerk system, with hyperbolic sine function as the only nonlinear term, is proposed, as a modification of a hyperjerk system proposed by Leutcho et al. First, a dynamical analysis on the system is performed and interesting phenomena concerning chaos theory, such as route to chaos, antimonotonicity, crisis, and coexisting attractors, are studied. For this reason, the system’s bifurcation diagrams with respect to different parameter values are plotted and its Lyapunov exponents are computed. Afterwards, the synchronization of the system is considered, using active control. The proposed system is then applied, as a chaotic generator, to the problem of chaotic path planning, using a combination of sampling and a modulo tactic technique.


---

## Body

## 1. Introduction

Over the last 60 years, chaotic systems have integrated almost every scientific discipline and have found applications in numerous fields, including, but not limited to, biology, engineering, finance, robotics, circuits, cryptography, secure communications, and many more (see, for example, [1–9] and the references therein). Due to their deterministic nature and high sensitivity to initial conditions, chaotic systems provide an efficient method of adding complexity and increased security to a design. In addition, the observance of a plethora of chaos-related behaviors in complex systems, natural or artificial, indicates that a wide spectrum of systems can be tagged as having chaotic behavior or having the potential of becoming chaotic when certain conditions regarding their parameters are met.For the above reasons, there is an ongoing research in creating novel chaotic systems and studying them in order to unmask their intricate dynamical behavior as well as their potential implementability in chaos-related applications. There are three main approaches in designing new chaotic systems from existing ones: firstly, by modifying an existing term in the differential equations that describe the system to a higher order nonlinear term, for example, by changing a quadratic power to cubic or a product term to an exponential term; secondly, by making adjustments to the existing terms without changing the type of nonlinearity; and thirdly, by adding more nonlinear terms to the system of differential equations.The approach considered here is the second one, where the system’s nonlinear term is slightly adjusted, giving a new chaotic system which is then applied to the problems of synchronization and to chaotic path planning. First, a novel hyperjerk system is proposed, as a modification of a system proposed by Leutcho et al. [10]. The proposed system has only one nonlinear term, yet it yields a plethora of interesting chaotic phenomena, which are analyzed and discussed, by considering its bifurcation diagrams, spectrum of Lyapunov exponents, and continuation diagrams. It is seen that the system has chaotic behavior for a wide spectrum of parameter values and also presents phenomena like antimonotonicity and coexisting attractors. Afterwards, the synchronization of the system with a 4D hyperchaotic Rössler system is considered, using the method of active control. The simulation results show an efficient synchronization time. At last, the proposed chaotic system is used as a chaotic motion generator and is applied to the problem of chaotic path planning. The system is first sampled and is then combined with a modulo tactic to generate a motion sequence for a robot moving in eight or four distinct directions in a given grid. Extensive simulations are performed to compare the coverage percent and the mean number of visits per cell to the number of movements performed, and the results are discussed. This design scheme showcases one of the many possible applications of the proposed chaotic system in practice.The rest of the paper is structured as follows. In Section2, the novel hyperjerk system is proposed and its structural characteristics are described. In Section 3, extensive simulations are performed to analyze the dynamical behavior of the system, including bifurcation diagrams, Lyapunov exponents, and continuation diagrams. In Section 4, the synchronization of the system is considered using active control. In Section 5, the application of the proposed chaotic system to the problem of chaotic path planning is proposed for a robot moving in eight or four directions. Section 6 concludes the paper.


## 2. The Proposed Dynamical System

In order to enrich the type of chaotic systems and meet the needs of the applications, it is necessary to improve the existing chaotic systems to get systems with higher nonlinearity and complexity [11]. As mentioned in the introduction, the optimization of chaotic systems is usually made by the following three easy modifications to the nonlinear term of the existing system. In the first approach, the nonlinear term is modified to a higher order nonlinear term, such as changing the product term to an exponential or logarithmic function [12–15], while in the second approach, simple adjustments to the nonlinear terms without affecting its order [16, 17] are made. As a third approach, more nonlinear terms can be added to the system.Therefore, in this work, the second approach has been adopted in order to modify a known hyperjerk system, which has been proposed by Leutcho et al. [10]. The aforementioned system is described as(1)x˙=y,y˙=z,z˙=dw,w˙=−aw−by−cx−γsinhz,while the modified system is considered as(2)x˙=y,y˙=z,z˙=w,w˙=−aw−by−cx−sinhz+d,in which the parameter d has been added to the variable z inside the system’s only nonlinear term, sinh(z + d), and γ has been set to 1.With this simple modification, the new system presents the following features:(1)
It is also a hyperjerk system described by equation (6) of Ref. [18], which presents the general form of hyperjerk systems as(3)d4xdt4+a0d3xdt3+a1d2xdt2+a2dxdt=gx,dxdt,d2xdt2,wherea0 = a, a1 = 0, a2 = b, and gx,dx/dt,d2x/dt2=−cx−sinhd2x/dt2+d.(2)
It is a structurally simple hyperjerk system with only one nonlinear term.(3)
By solving(4)y=0,z=0,w=0,−aw−by−cx−sinhz+d=0,it is easy to verify that the system has an equilibrium at−sinhd/c,0,0,0, in contrast to the system of Ref. [10] which has an equilibrium located at the origin of system coordinates.(4)
Chaotic systems can be categorized as either dissipative or conservative. A dissipative system has a property that trajectories are attracted to a bounded region of state space by a strange attractor and its divergence is negative. In contrast, a conservative chaotic system does not have a state space attractor and its divergence is zero implying that a volumeV0 = V(t = 0) in state space occupied by a set of initial conditions remains unchanged.So, the volume contraction/expansion rateΛ=V−1dV/dt of the hyperjerk system modeled by (2) at any given point (x, y, z, w) of the state space is computed in order to gain information about the existence of attractive sets. The following expression is obtained:(5)Λ=dx˙/dx+dy˙/dy+dz˙/dz+dw˙/dw=−a<0,which is the same as in system (1) of Ref. [10]. Therefore, system (2) is also a dissipative system. So, it converges in the index of the following form dV/dt = V0e−at. This means that each volume, containing the trajectories of system (2), will shrink to zero when t ⟶ ∞ with the exponential rate V0e−at. Thus, all the orbits of system (2) are finally confined to a specific subset of zero volume, and the asymptotic motion of system (5) will be fixed on an attractor of the system.(5)
System (2) is invariant under the transformation (x, y, z, w; d) ⟶ (−x, −y, −z, −w; −d) owing to the presence of the specific hyperbolic sine nonlinearity. In more detail, if the initial conditions (x, y, z, w) are a solution of system (2) for a given value of parameter d, then (−x, −y, −z, −w) is also a solution for the value of parameter −d and for the same parameter set. In Figure 1, this phenomenon, under the aforementioned transformation, is shown in two different cases of attractors (periodic and chaotic attractors).Figure 1
(a) Periodic attractors for (x, y, z, w)0 = (0.1, 0.2, 0.3, 0.4), d = 0.3 and (x, y, z, w)0 = (−0.1, −0.2, −0.3, −0.4), d = −0.3, while a = 1.4, b = 0.3, and c = 1.4. (b) Chaotic attractors for (x, y, z, w)0 = (0.1, 0.2, 0.3, 0.4), d = 0.8 and (x, y, z, w)0 = (−0.1, −0.2, −0.3, −0.4), d = −0.8, while a = 1.4, b = 0.3, and c = 1.4.
(a)
(b)As an additional simulation, different phase portraits of the system are shown in Figure2, for initial conditions x0,y0,z0,w0=0.1,0.2,0.3,0.4, for parameter values a=1.4,b=0.3,c=1.4,andd=0.85, and a simulation time of 800 sec.Figure 2
Phase portraits for initial conditionsx0,y0,z0,w0=0.1,0.2,0.3,0.4, parameter values a=1.4,b=0.3,c=1.4,andd=0.85, and simulation time of 800 sec.


## 3. Simulation Results

In this section, the analysis of the dynamical behavior of the proposed system (2) for initial conditions (x0, y0, z0, w0) = (0.1, 0.2, 0.3, 0.4) is performed, by using well-known tools of nonlinear theory, such as bifurcation diagram and Lyapunov exponents. For the investigation of the proposed system’s dynamical behavior, system (2) is integrated numerically using the classical fourth-order Runge–Kutta integration algorithm. For each set of the chosen parameters used in this work, the time step is always Δt = 0.001 and the calculations are performed using variables and parameters in extended precision mode. For each parameter setting, the system is integrated for sufficiently long time and the transient is discarded.Ιn more detail, for revealing the dynamics of system (2), regarding the value of parameter d, the system’s bifurcation diagram of x versus d has been plotted (Figure 3(a)). The bifurcation diagram is obtained by plotting the variable x when the trajectory cuts the plane y = 0 with dy/dt < 0, as the control parameter d increases in very small steps in the range 0.3 ≤ d ≤ 1.4. As illustrated in Figure 3(a), system (2) inserts to a chaotic region, which is interrupted by small periodic windows, through a period-doubling route, as the parameter d increases. It should be pointed out that if we choose the inverse direction (dy/dt > 0), a different set of points in the diagram will be produced. However, the system’s behavior will be the same. Furthermore, the spectrum of Lyapunov exponents (LEs) versus the parameter d is depicted in Figure 3(b). From this diagram, the system’s chaotic behavior as the control parameter increases can be confirmed, as the bifurcation diagram coincides with the spectrum of the Lyapunov exponents. Consequently, it is known from the nonlinear theory in periodic regions that the system only has the maximum Lyapunov exponent equal to zero, while in chaotic regions, the system has a positive maximum Lyapunov exponent.Figure 3
(a) Bifurcation diagram ofx versus d and (b) the spectrum of system’s Lyapunov exponents, for a = 1.4, b = 0.3, and c = 1.4 and initial conditions (x0, y0, z0, w0) = (0.1, 0.2, 0.3, 0.4).
(a)
(b)Also, studying the bifurcation diagram of Figure3(a), another interesting phenomenon can be observed. This is antimonotonicity, which has been introduced by Dawson et al. [19]. This is a fundamental phenomenon in bifurcations for a large class of nonlinear dissipative systems [20–22], where periodic orbits are not only created but also destroyed, when one increases the control parameter monotonically (smoothly) in any neighborhood of a homoclinic tangency value, leading to reversals of period-doubling cascades. So, antimonotonicity is due to tangential intersections between the stable and unstable manifolds of a system. According to this phenomenon, the system enters chaos with a period-doubling route (period-1 ⟶ period-2 ⟶ … ⟶ chaos) and exits from the chaos following the reverse period-doubling route (chaos ⟶ … ⟶ period-2 ⟶ period-1). As a result, in the bifurcation diagram, the shape of a chaotic bubble is formed. This phenomenon is very interesting because it describes the complex scenario of how a nonlinear system creates or destroys unstable periodic orbits by parameter alterations.The phenomenon of antimonotonicity is clearly presented in the bifurcation diagrams ofx versus d (Figure 4), for various values of parameter c, by keeping a = 1.4 and b = 0.3. In more detail, for c = 1.1, displays the so-called “primary bubble,” due to the fact that system (2) has only one period-doubling. However, as the value of c increases, bifurcation diagrams of system (2) become more complex, as smaller bubbles belonging to higher period orbits appear on the primary bubble (Figures 4(b)–4(d)). Finally, for c > 1.18 small chaotic bubbles appear in the bifurcation diagrams (Figures 4(e) and 4(f)), which gradually increase, as parameter c also increases. So, more complex structures are formed in the bifurcation diagram, such as the extended chaotic region which appears in the bifurcation diagram of Figure 3(a), for c = 1.4.Figure 4
Bifurcation diagrams ofx versus d for (a) c = 1.1, (b) c = 1.15, (c) c = 1.17, (d) c = 1.178, (e) c = 1.18, and (f) c = 1.2 with a = 1.4 and b = 0.3.
(a)
(b)
(c)
(d)
(e)
(f)Next, the phenomenon of coexisting attractors is studied. For this reason, the continuation diagram is obtained by plotting the variablex when the trajectory cuts the plane y = 0 with dy/dt < 0, as the control parameter, which in this case the parameter c has been chosen, is increased (or decreased) in tiny steps in the range of 1.87 ≤ c ≤ 1.92. The difference between the bifurcation diagram and the continuation bifurcation diagram is that in the bifurcation diagram, each simulation performed for a value of the control parameter starts from the same initial conditions, while in a continuation diagram, the initial conditions in each simulation are chosen as the final value that the system had in its previous simulation. So, in the continuation bifurcation diagram, the final states at each parameter serve as the initial conditions for the next value of the control parameter. This is considered for the cases where the simulations are performed also for increasing and decreasing values of parameter c. In the graph of Figure 5, two sets of data, one for increasing values (black) of c and one for decreasing values (red) of c, are superimposed. From this comparison, the range of values of control parameter c, in which the system develops coexisting attractors, are localized. When the continuation bifurcation diagrams of Figure 5 are different, it results in coexisting attractors.Figure 5
Continuation bifurcation diagrams of system (2) for increasing (black) and decreasing (red) values of c in the range of 1.87 ≤ c ≤ 1.92, for a = 1.4, b = 0.3, and d = 0.9.In more detail, from the continuation diagrams of Figure5, it is easy to verify the coexistence of attractors of system (2) in the range of 1.87 ≤ c ≤ 1.92. In the case of decreasing the control parameter, the system passes from periodic to chaotic behavior for lower value of c in contrary to the case of increasing the value of the control parameter. So, chaotic and period-3 attractors coexist with period-1, as presented in Figure 6.Figure 6
Coexisting attractors for (a)c = 1.897 and (b) c = 1.9, with a = 1.4, b = 0.3, and d = 0.9 and different initial conditions of system (2).
(a)
(b)Also, from the continuation diagrams of Figure5, it is obvious that the system enters chaos (or exits chaos) through a crisis phenomenon as the value of c decreases (or increases). The occurrence of sudden qualitative changes of chaotic dynamics as a parameter varies is known from 1983. That year, Grebogi et al. [23] reported that such changes may result from the collision of an unstable periodic orbit and a coexisting chaotic attractor. Since then, crisis phenomena have been observed in many dynamical systems [7, 24, 25].


## 4. Synchronization

In this section, the problem of synchronization of the proposed chaotic system is considered. The control of chaotic systems is overall considered a challenging task, given their complex dynamics and their sensitivity to initial conditions and parameter changes. A problem closely related to chaos control is synchronization, which appears when two chaotic systems are coupled, with the secondary system called slave system being driven by a master system. The aim here is to have the trajectory of the two systems converge, as time approaches infinity. Synchronization was first introduced by Pecora and Carroll [26], and over the years, it has found applications in cryptography, secure communications, and many more [4, 5, 27–32]. Here, the method of active control [33–35] is used to synchronize the proposed system with a 4D Rössler type hyperchaotic system [36].The method of active control consists of two parts. Initially, the error dynamics of the master and slave system are considered. Then, the control input is constructed asut=kt+at. The first control part kt is constructed in order to discard all terms that cannot be written in terms of the error so that the error dynamics is transformed to a linear system. Then, the part at is constructed such that the error dynamics take the form e˙=Aactivee, where Aactive is a stable matrix, and so the error will converge to zero.Here, the proposed hyperjerk system (2) is taken as the master system:(6)x˙1=x2,x˙2=x3,x˙3=x4,x˙4=−ax4−bx2−cx1−sinhx3+d,where to avoid confusion, x,y,z,w have been replaced by x1,x2,x3,x4, and the hyperchaotic Rössler system from [36] is considered as the slave system:(7)y˙1=b3y1+b2y2+u1,y˙2=2+y2y3+u2,y˙3=−y2+y4+u3,y˙4=y1+y3+b1y4+u4,where ui are the control inputs to be designed to achieve synchronization. Taking the difference between the states of the two systems ei=yi−xi, the error dynamics are given as(8)e˙1=b2e2+b3e1+b2−1x2+b3x1+u1,e˙2=3+e2e3+e2x3+x2e3−x3+u2,e˙3=−e2−e4−x2−2x4+u3,e˙4=e1+e3+1+cx1+x3+b1e4+y4+bx2+ax4+sinhx3+d+u4.Now, the control inputs are computed as(9)ui=Aactivei+ki,where Aactivei are the linear control terms to be computed in the following, and(10)k1=−b2−1x2−b3x1,k2=−3−e2e3−e2x3−x2e3+x3,k3=+x2+2x4,k4=−1+cx1−x3−a+b1x4−bx2−sinhx3+d.Now, substituting (9) into the error dynamics (8) results in(11)e˙1=b3e1+b2e2+Aactive1,e˙2=Aactive2,e˙3=−e2−e4+Aactive3,e˙4=e1+e3+b1e4+Aactive4.In view of (11), we must choose Aactivei such that the error dynamics become asymptotically stable. One such choice is the following:(12)Aactive=Aactive1Aactive2Aactive3Aactive4=−b3+3−b2000−3−1001−31−10−1−b1+3e1e2e3e4.Overall, using active control (9), the error dynamics converge to zero and the trajectories of the slave system synchronize with the master system. Figure 7 shows the simulation results for synchronization, and Figure 8 shows the synchronization error, which converges to zero. The system parameter values for the two systems are a=1.4,b=3,c=1.2,andd=0.85 and b1=0.25,b2=−0.5,andb3=0.05 and the initial conditions are x1,x2,x3,x40=0.01,0.01,0.01,0.01 and y1,y2,y3,y40=0.4,−0.4,0.9,0.9. Using a similar approach with small modifications, the synchronization between two hyperjerk systems (2) with different parameters can also be considered.Figure 7
Synchronization between hyperjerk system (2) and hyperchaotic Rössler systems, for parameter values a=1.4,b=3,c=1.2,andd=0.85 and b1=0.25,b2=−0.5,andb3=0.05 and initial conditions x1,x2,x3,x40=0.01,0.01,0.01,0.01 and y1,y2,y3,y40=0.4,−0.4,0.9,0.9.
(a)
(b)
(c)
(d)Figure 8
Synchronization errors.
(a)
(b)
(c)
(d)


## 5. Application to Chaotic Path Planning

An emerging application of chaotic systems is the chaotic path planning [8, 37]. In this problem, the aim is to utilize a chaotic system in order to design a motion trajectory for a robot moving in a given space, with the aim of covering the complete grid [38]. Using a chaotic system to generate the trajectory of the robot leads to a movement that is unpredictable and looks arbitrary to an outside observer. Thus, chaotic path planning is an effective approach for grid coverage that is attractive to applications that consider obstacles, surveillance, search and rescue missions, and domestic applications [39–49].For chaotic path planning, there are several approaches used in the literature, depending on the robot motion capabilities that are considered. For differential motion robots, the design consists of replacing the linear and directional velocities of the robot by a combination of the states of a chosen chaotic system [48–53]. For robots moving in discrete directions, that is, four or eight, the approach that is mostly considered consists of first creating a chaotic random bit generator using a chaotic system and then applying the resulting bit sequence in generating the chaotic motion of the system [54–57]. To avoid uneven coverage, a problem that was also addressed in [50, 58], a deskewing is first performed to the chaotic bit sequence. This approach yields good results but has the disadvantage that it requires a redundant number of iterations of the chaotic system in order to produce the desired number of movements. Here, in order to avoid this, an alternative approach is considered. By using a modulo operator on the chaotic system, a chaotic series of motions can be generated that yields an even grid coverage with satisfying percentage.


### 5.1. Motion in Eight Directions

Since the motion of the robot consists of discrete direction commands, we first need to consider the sampling of the chaotic system. So initially, a sampling periodts is chosen and a series of N samples is obtained from the original system (2).(13)Di=Xi⋅ts,i=0,…,N−1.The sample valueN is chosen according to the number of path commands that are desired. The choice of sampling period can vary, and depending on the sampling period, the original system may need to be simulated for a different duration in order to obtain the required samples. For example, if N = 5,000 and ts=1sec, then the system must be simulated for 5,000 seconds to obtain the required samples, while if ts=0.1sec, then a simulation for 500 seconds is enough. We will see in the following that the sampling period does not have a significant effect on the path planning generator.After obtaining the required samples from the continuous time system, a modulo tactic is applied to the sequenceDi as follows:(14)si=106⋅d1id2id3id4imod8,where ⋅ denotes the integer part of the argument and the subscript i denotes the sample of each state. The product d1id2id3id4i of the sampled states is used to generate a random chaotic sequence for the modulo tactic. This was chosen after some appropriate testing, but it is possible to use other combinations of the sampled states with equally good results. So, (14) generates a sequence of integers ranging from 0 to 7. Based on si, the following motion commands are generated:(15)mi=up,si=0,up−right,si=1,right,si=2,down−right,si=3,down,si=4,down−left,si=5,left,si=6,up−left,si=7.This command sequence can then be applied to the autonomous robot exploring a given area. Here, it is assumed that the robot can move in eight different directions. Figure9 depicts the four states of the chaotic system, and Figure 10 depicts their product x1x2x3x4, for a simulation time of 500 sec, starting from initial conditions x10,x20,x30,x40=0.1,0.1,0.1,0.1, for parameter values a=1.4,b=3,c=1.4,andd=0.85. For these sets of parameters and initial conditions, system (2) is chaotic, as it is verified by the calculation of system’s Lyapunov exponents, which are LE1 = 0.15266, LE2 = 0, LE3 = −0.63290, and LE4 = −1.53962. The power spectrum of the product x1x2x3x4, is shown in Figure 11, from which it can be seen that the spectrum of the circuit is continuous. There are peaks in the spectrum and the peaks are connected together, from which its chaotic behavior is further demonstrated.Figure 9
States of the chaotic system for initial valuesx10,x20,x30,x40=0.1, 0.1, 0.1, 0.1 and parameter values a = 1.4, b = 3, c = 1.4, and d = 0.85.Figure 10
Product of the states for initial valuesx10,x20,x30,x40=0.1, 0.1, 0.1, 0.1 and parameter values a = 1.4, b = 3, c = 1.4, and d = 0.85.Figure 11
Power spectrum (FFT) of the product signalx1x2x3x4 with a = 1.4, b = 3, c = 1.4, and d = 0.85.Figure12 depicts the histogram for N = 1,000,000 values of si, obtained with a simulation time of 100,000 sec and a sampling time of ts=0.1sec, using the same initial conditions and parameter values. It is seen that there is a uniform distribution among the different direction commands. This is a first indication that the coverage performance will give satisfactory results, since an uneven distribution of the motion commands would clearly lead to the robot moving in a skewed way on the grid. A similar approach using the modulo operator was considered in [46], where the discrete time logistics map was used. In that work, the chaotic system was discrete, so the use of sampling and its effect on the coverage performance were not explored.Figure 12
Histogram for 1,000,000 values of the proposed modulo tactics.For the simulation of the robot movement, we consider a100×100 grid, consisting of 1002 distinct cells. The robot makes a move with each value of si, according to the above rule (15). If the generated direction is unacceptable, like moving outside the defined limits or meeting an obstacle, then the robot remains in its place and awaits for the next motion command. In Figure 13, the motion of the robot is shown, starting from initial position 50,50T, for 35,000 steps. These are generated from the proposed system for initial values chosen randomly from the interval [0, 1] and a sampling time of ts=0.1sec. This means that to generate the required samples, the system was simulated for 3,500 sec. The system was solved in Matlab using ode45, with a tolerance of 10−5. It is seen that a coverage of 75% is achieved, which is very satisfying. Figure 14 depicts a color-coded graph showing the number of visits in each cell. This gives a better visualization with respect to the number of times the robot revisits cells in the grid.Figure 13
Grid coverage for 35,000 steps.Figure 14
Color-coded grid coverage for 35,000 moves.


### 5.2. Motion in Four Directions

In the case that the movement of the robot is limited to only four directions, the same technique can be performed using a modulo 4 operator instead. Thus, the motion sequence is adjusted as(16)si=106⋅di,1di,2di,3di,4mod4,mi=up,si=0,right,si=1,down,si=2,left,si=3.Figure15 depicts the histogram for N = 1,000,000 values of si, obtained with a simulation time of 100,000 sec and a sampling time of ts=0.1sec and parameter values the same as before. Also, Figures 16 and 17 show the path coverage, starting from initial position 50,50T, for 35,000 steps. The initial values are chosen randomly from the interval [0, 1], and the sampling time is ts=0.1sec.Figure 15
Histogram for 1,000,000 values of the modulo 4 tactics.Figure 16
Grid coverage for 35,000 steps for motion in 4 directions.Figure 17
Color-coded grid coverage for 35,000 moves for motion in 4 directions.


### 5.3. Coverage

To study the connection between the grid coverage and the number of performed moves, we considered the average of 10 simulations, performed for a varying number of steps, for parameter valuesa=1.4,b=3,c=1.4,andd=0.85. For each simulation, an initial value for the chaotic system and an initial position are chosen randomly. Then, a series of moves are generated, for three different sampling periods, and the coverage is then computed and finally averaged. The coverage (C) is computed by(17)C=1M∑i=1MIi,where M is the number of cells and I is the coverage situation of each cell.(18)Ii=1,cellihas been visited,0,cellihas not been visited.The simulation results are shown in Figure18. Here, it is observed that the sampling time does not significantly affect the performance of the system. This is due to the high sensitivity of the modulo command to changes in the values of the first six decimal digits of the computed value. Thus, it is seen that this approach is very robust and can be potentially used in combination with many different techniques and other systems as well. The curve fitting performed for the coverage percent showed an exponential convergence rate of the form(19)yx=y0+Aex/R,where the values of the parameters are shown in Table 1.Figure 18
Coverage versus number of iterations.Table 1
Curve fitting for coverage.
Moves
Sampling (sec)
y
0
A
R
R-square
8
1
100.1143
−94.03029
−34054.97
0.99632
0.1
99.4387
−96.42309
−31546.65
0.99663
0.01
99.19168
−97.76927
−30040.11
0.99571
4
1
97.29593
−95.88355
−39574.59
0.99556
0.1
100.35984
−96.48292
−44331.05
0.99493
0.01
97.5818
−96.348
−39907.35
0.99662Also, Figure19 shows the average amount of visits per visited cells for varying number of steps. Here, it is seen that the number of visits changes linearly as the number of steps rises. The fitting performed in Matlab indeed resulted in a linear curve.(20)yx=p1x+p2,where the values of the parameters are given in Table 2.Figure 19
Mean number of visits per visited cells.Table 2
Linear curve fitting for number of visits per visited cell.
Moves
Sampling (sec)
p
1
p
2
R-square
8
1
7.967·10−5
2.515
0.9974
0.1
7.948·10−5
2.483
0.9954
0.01
7.922·10−5
2.446
0.9939
4
1
7.591·10−5
3.519
0.9932
0.1
7.479·10−5
3.62
0.9918
0.01
7.516·10−5
3.583
0.9920Overall, it is seen from the simulations that movement in eight directions yields better results compared to movement in four. This was expected, since movement in four directions is much more limited and has a higher possibility (1 in 4) of revisiting a cell in each iteration, compared to movement in eight directions, where the possibility of revisiting the previous position is 1 in 8. That is why the coverage rate is higher for eight directions, while the number of visits per visited cell is higher for the case of four directions.


### 5.4. Discussion

For the application to path planning, it is seen that the modulo technique provides a robust approach to generating a chaotic path that gives even coverage of the terrain. The path planning generator algorithm is very easy to implement in any programming environment, using a 16 digit accuracy. It can also be easily adjusted to different scenarios, like grids with obstacles, or labyrinth-like environments. Also, this chaotic path planning technique can be potentially applied to robots with different moving capabilities [59–61] and can also be combined with other advanced techniques, like the inverse pheromone approach [46, 56], multiagent and swarm robotics [62, 63], and more [64–66].


## 5.1. Motion in Eight Directions

Since the motion of the robot consists of discrete direction commands, we first need to consider the sampling of the chaotic system. So initially, a sampling periodts is chosen and a series of N samples is obtained from the original system (2).(13)Di=Xi⋅ts,i=0,…,N−1.The sample valueN is chosen according to the number of path commands that are desired. The choice of sampling period can vary, and depending on the sampling period, the original system may need to be simulated for a different duration in order to obtain the required samples. For example, if N = 5,000 and ts=1sec, then the system must be simulated for 5,000 seconds to obtain the required samples, while if ts=0.1sec, then a simulation for 500 seconds is enough. We will see in the following that the sampling period does not have a significant effect on the path planning generator.After obtaining the required samples from the continuous time system, a modulo tactic is applied to the sequenceDi as follows:(14)si=106⋅d1id2id3id4imod8,where ⋅ denotes the integer part of the argument and the subscript i denotes the sample of each state. The product d1id2id3id4i of the sampled states is used to generate a random chaotic sequence for the modulo tactic. This was chosen after some appropriate testing, but it is possible to use other combinations of the sampled states with equally good results. So, (14) generates a sequence of integers ranging from 0 to 7. Based on si, the following motion commands are generated:(15)mi=up,si=0,up−right,si=1,right,si=2,down−right,si=3,down,si=4,down−left,si=5,left,si=6,up−left,si=7.This command sequence can then be applied to the autonomous robot exploring a given area. Here, it is assumed that the robot can move in eight different directions. Figure9 depicts the four states of the chaotic system, and Figure 10 depicts their product x1x2x3x4, for a simulation time of 500 sec, starting from initial conditions x10,x20,x30,x40=0.1,0.1,0.1,0.1, for parameter values a=1.4,b=3,c=1.4,andd=0.85. For these sets of parameters and initial conditions, system (2) is chaotic, as it is verified by the calculation of system’s Lyapunov exponents, which are LE1 = 0.15266, LE2 = 0, LE3 = −0.63290, and LE4 = −1.53962. The power spectrum of the product x1x2x3x4, is shown in Figure 11, from which it can be seen that the spectrum of the circuit is continuous. There are peaks in the spectrum and the peaks are connected together, from which its chaotic behavior is further demonstrated.Figure 9
States of the chaotic system for initial valuesx10,x20,x30,x40=0.1, 0.1, 0.1, 0.1 and parameter values a = 1.4, b = 3, c = 1.4, and d = 0.85.Figure 10
Product of the states for initial valuesx10,x20,x30,x40=0.1, 0.1, 0.1, 0.1 and parameter values a = 1.4, b = 3, c = 1.4, and d = 0.85.Figure 11
Power spectrum (FFT) of the product signalx1x2x3x4 with a = 1.4, b = 3, c = 1.4, and d = 0.85.Figure12 depicts the histogram for N = 1,000,000 values of si, obtained with a simulation time of 100,000 sec and a sampling time of ts=0.1sec, using the same initial conditions and parameter values. It is seen that there is a uniform distribution among the different direction commands. This is a first indication that the coverage performance will give satisfactory results, since an uneven distribution of the motion commands would clearly lead to the robot moving in a skewed way on the grid. A similar approach using the modulo operator was considered in [46], where the discrete time logistics map was used. In that work, the chaotic system was discrete, so the use of sampling and its effect on the coverage performance were not explored.Figure 12
Histogram for 1,000,000 values of the proposed modulo tactics.For the simulation of the robot movement, we consider a100×100 grid, consisting of 1002 distinct cells. The robot makes a move with each value of si, according to the above rule (15). If the generated direction is unacceptable, like moving outside the defined limits or meeting an obstacle, then the robot remains in its place and awaits for the next motion command. In Figure 13, the motion of the robot is shown, starting from initial position 50,50T, for 35,000 steps. These are generated from the proposed system for initial values chosen randomly from the interval [0, 1] and a sampling time of ts=0.1sec. This means that to generate the required samples, the system was simulated for 3,500 sec. The system was solved in Matlab using ode45, with a tolerance of 10−5. It is seen that a coverage of 75% is achieved, which is very satisfying. Figure 14 depicts a color-coded graph showing the number of visits in each cell. This gives a better visualization with respect to the number of times the robot revisits cells in the grid.Figure 13
Grid coverage for 35,000 steps.Figure 14
Color-coded grid coverage for 35,000 moves.


## 5.2. Motion in Four Directions

In the case that the movement of the robot is limited to only four directions, the same technique can be performed using a modulo 4 operator instead. Thus, the motion sequence is adjusted as(16)si=106⋅di,1di,2di,3di,4mod4,mi=up,si=0,right,si=1,down,si=2,left,si=3.Figure15 depicts the histogram for N = 1,000,000 values of si, obtained with a simulation time of 100,000 sec and a sampling time of ts=0.1sec and parameter values the same as before. Also, Figures 16 and 17 show the path coverage, starting from initial position 50,50T, for 35,000 steps. The initial values are chosen randomly from the interval [0, 1], and the sampling time is ts=0.1sec.Figure 15
Histogram for 1,000,000 values of the modulo 4 tactics.Figure 16
Grid coverage for 35,000 steps for motion in 4 directions.Figure 17
Color-coded grid coverage for 35,000 moves for motion in 4 directions.


## 5.3. Coverage

To study the connection between the grid coverage and the number of performed moves, we considered the average of 10 simulations, performed for a varying number of steps, for parameter valuesa=1.4,b=3,c=1.4,andd=0.85. For each simulation, an initial value for the chaotic system and an initial position are chosen randomly. Then, a series of moves are generated, for three different sampling periods, and the coverage is then computed and finally averaged. The coverage (C) is computed by(17)C=1M∑i=1MIi,where M is the number of cells and I is the coverage situation of each cell.(18)Ii=1,cellihas been visited,0,cellihas not been visited.The simulation results are shown in Figure18. Here, it is observed that the sampling time does not significantly affect the performance of the system. This is due to the high sensitivity of the modulo command to changes in the values of the first six decimal digits of the computed value. Thus, it is seen that this approach is very robust and can be potentially used in combination with many different techniques and other systems as well. The curve fitting performed for the coverage percent showed an exponential convergence rate of the form(19)yx=y0+Aex/R,where the values of the parameters are shown in Table 1.Figure 18
Coverage versus number of iterations.Table 1
Curve fitting for coverage.
Moves
Sampling (sec)
y
0
A
R
R-square
8
1
100.1143
−94.03029
−34054.97
0.99632
0.1
99.4387
−96.42309
−31546.65
0.99663
0.01
99.19168
−97.76927
−30040.11
0.99571
4
1
97.29593
−95.88355
−39574.59
0.99556
0.1
100.35984
−96.48292
−44331.05
0.99493
0.01
97.5818
−96.348
−39907.35
0.99662Also, Figure19 shows the average amount of visits per visited cells for varying number of steps. Here, it is seen that the number of visits changes linearly as the number of steps rises. The fitting performed in Matlab indeed resulted in a linear curve.(20)yx=p1x+p2,where the values of the parameters are given in Table 2.Figure 19
Mean number of visits per visited cells.Table 2
Linear curve fitting for number of visits per visited cell.
Moves
Sampling (sec)
p
1
p
2
R-square
8
1
7.967·10−5
2.515
0.9974
0.1
7.948·10−5
2.483
0.9954
0.01
7.922·10−5
2.446
0.9939
4
1
7.591·10−5
3.519
0.9932
0.1
7.479·10−5
3.62
0.9918
0.01
7.516·10−5
3.583
0.9920Overall, it is seen from the simulations that movement in eight directions yields better results compared to movement in four. This was expected, since movement in four directions is much more limited and has a higher possibility (1 in 4) of revisiting a cell in each iteration, compared to movement in eight directions, where the possibility of revisiting the previous position is 1 in 8. That is why the coverage rate is higher for eight directions, while the number of visits per visited cell is higher for the case of four directions.


## 5.4. Discussion

For the application to path planning, it is seen that the modulo technique provides a robust approach to generating a chaotic path that gives even coverage of the terrain. The path planning generator algorithm is very easy to implement in any programming environment, using a 16 digit accuracy. It can also be easily adjusted to different scenarios, like grids with obstacles, or labyrinth-like environments. Also, this chaotic path planning technique can be potentially applied to robots with different moving capabilities [59–61] and can also be combined with other advanced techniques, like the inverse pheromone approach [46, 56], multiagent and swarm robotics [62, 63], and more [64–66].


## 6. Conclusion

In this work, a novel hyperjerk system was proposed as a modification of the system proposed by Leutcho et al. [10]. A thorough dynamical analysis showed a rich chaotic behavior for the system. The synchronization of the system was then considered using a simple and effective technique of active control. The application of the proposed system to chaotic path planning was then proposed and studied. After appropriate sampling of the continuous chaotic system and a combination with a modulo tactic, the results showed a satisfying coverage rate. Thus, it is seen that the proposed chaotic system is suitable for various chaos and synchronization related applications. Future research directions can address the circuit realization of the system, a continuation of the parameter analysis to uncover more chaotic phenomena, the application of synchronization to the problem of secure communication, the implementation of the path planning to an experimental robotic setup, the use of different chaotic systems in path planning [67], a comparative study between the different chaotic path planning methods that exist, and the application of the proposed system to chaotic random bit generation and image encryption.



---

*Source: 2826850-2020-04-30.xml*