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

,


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][2][3][4][5][6][7][8][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.
ere 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. e 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]. e 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. e 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. e 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. is design scheme showcases one of the many possible applications of the proposed chaotic system in practice. e rest of the paper is structured as follows. In Section 2, 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.

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][13][14][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. erefore, 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]. e aforementioned system is described as in which the parameter d has been added to the variable z inside the system's only nonlinear term, sinh(z + d), and c 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 (2) It is a structurally simple hyperjerk system with only one nonlinear term.
it is easy to verify that the system has an equilibrium at (− (sinh(d)/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 volume V 0 � V(t � 0) in state space occupied by a set of initial conditions remains unchanged. So, the volume contraction/expansion rate Λ � V − 1 (dV/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. e following expression is obtained: which is the same as in system (1) of Ref. [10]. erefore, system (2) is also a dissipative system. So, it converges in the index of the following form dV/ dt � V 0 e − at . is means that each volume, containing 2 Complexity the trajectories of system (2), will shrink to zero when t ⟶ ∞ with the exponential rate V 0 e − at . us, 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.
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, 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).
As an additional simulation, different phase portraits of the system are shown in Figure 2

Simulation Results
In this section, the analysis of the dynamical behavior of the proposed system (2) for initial conditions (x 0 , y 0 , z 0 , w 0 ) � (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)). e 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.
Also, studying the bifurcation diagram of Figure 3(a), another interesting phenomenon can be observed. is is antimonotonicity, which has been introduced by Dawson et al. [19]. is is a fundamental phenomenon in bifurcations for a large class of nonlinear dissipative systems [20][21][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 perioddoubling route (period-1 ⟶ period-2 ⟶ . . . ⟶ chaos) and exits from the chaos following the reverse perioddoubling route (chaos ⟶ . . . ⟶ period-2 ⟶ period-1). As a result, in the bifurcation diagram, the shape of a chaotic bubble is formed.
is phenomenon is very interesting because it describes the complex scenario of how a nonlinear system creates or destroys unstable periodic orbits by parameter alterations. e phenomenon of antimonotonicity is clearly presented in the bifurcation diagrams of x 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.
Next, the phenomenon of coexisting attractors is studied. For this reason, the continuation diagram is obtained by plotting the variable x 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. e 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. is 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 Complexity continuation bifurcation diagrams of Figure 5 are different, it results in coexisting attractors.
In more detail, from the continuation diagrams of Figure 5, 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.
Also, from the continuation diagrams of Figure 5, it is obvious that the system enters chaos (or exits chaos) through a crisis phenomenon as the value of c decreases (or increases). e occurrence of sudden qualitative changes of chaotic dynamics as a parameter varies is known from 1983.
at 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].

Synchronization
In this section, the problem of synchronization of the proposed chaotic system is considered. e 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. e 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][28][29][30][31][32]. Here, the method of active control [33][34][35] is used to synchronize the proposed system with a 4D Rössler type hyperchaotic system [36]. e method of active control consists of two parts. Initially, the error dynamics of the master and slave system are considered. en, the control input is constructed as u(t) � k(t) + a(t). e first control part k(t) 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. en, the part a(t) is constructed such that the error dynamics take the form _ e � A active e, where A active is a stable matrix, and so the error will converge to zero.
Here, the proposed hyperjerk system (2) is taken as the master system: where to avoid confusion, x, y, z, w have been replaced by x 1 , x 2 , x 3 , x 4 , and the hyperchaotic Rössler system from [36] is considered as the slave system: where u i are the control inputs to be designed to achieve synchronization. Taking the difference between the states of the two systems e i � y i − x i , the error dynamics are given as Now, the control inputs are computed as where A active i are the linear control terms to be computed in the following, and Now, substituting (9) into the error dynamics (8) In view of (11), we must choose A active i such that the error dynamics become asymptotically stable. One such choice is the following:

Complexity 5
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. e system parameter values for the two systems are a � 1.4, b � 3, c � 1.2, and d � 0.85 and b 1 � 0.25, b 2 � − 0.5, and b 3 � 0.05 and the initial conditions are (x 1 , x 2 , x 3 , x 4 ) 0 � (0.01, 0.01, 0.01, 0.01) and (y 1 , y 2 , y 3 , y 4 ) 0 � (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.

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. us, 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][40][41][42][43][44][45][46][47][48][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][49][50][51][52][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][55][56][57]. To avoid uneven coverage, a problem that was also addressed in [50,58], a deskewing is first performed to the chaotic bit sequence. is 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.

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 period t s is chosen and a series of N samples is obtained from the original system (2).
e sample value N is chosen according to the number of path commands that are desired. e 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 t s � 1 sec, then the system must be simulated for 5,000 seconds to obtain the required samples, while if t s � 0.1 sec, 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 sequence D i as follows: where [·] denotes the integer part of the argument and the subscript i denotes the sample of each state. e product  Complexity of the sampled states is used to generate a random chaotic sequence for the modulo tactic. is 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 s i , the following motion commands are generated: is 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. Figure 9 depicts the four states of the chaotic system, and Figure 10 depicts their product x 1 x 2 x 3 x 4 , for a simulation time of 500 sec, starting from initial conditions (x 1 (0), x 2 (0), x 3 (0), x 4 (0)) � (0.1, 0.1, 0.1, 0.1), for parameter values a � 1.4, b � 3, c � 1.4, and d � 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 LE 1 � 0.15266, LE 2 � 0, LE 3 � − 0.63290, and LE 4 � − 1.53962. e power spectrum of the product x 1 x 2 x 3 x 4 , is shown in Figure 11, from which it can be seen that the spectrum of the circuit is continuous. ere are peaks in the spectrum and the peaks are connected together, from which its chaotic behavior is further demonstrated. Figure 12 depicts the histogram for N � 1,000,000 values of s i , obtained with a simulation time of 100,000 sec and a sampling time of t s � 0.1 sec, using the same initial conditions and parameter values. It is seen that there is a uniform distribution among the different direction commands. is 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.
For the simulation of the robot movement, we consider a 100 × 100 grid, consisting of 100 2 distinct cells. e robot makes a move with each value of s i , 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,50] T , for 35,000 steps.
ese are generated from the proposed system for initial values chosen randomly from the interval [0, 1] and a sampling time of t s � 0.1 sec. is means that to generate the required samples, the system was simulated for 3,500 sec.
e 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. is gives a better visualization with respect to the number of times the robot revisits cells in the grid.

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. us, the motion sequence is adjusted as

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 values a � 1.4, b � 3, c � 1.4, and d � 0.85. For each simulation, an initial value for the chaotic system and an initial position are chosen randomly. en, a series of moves are generated, for three different sampling periods, and the coverage is then computed and finally averaged. e coverage (C) is computed by where M is the number of cells and I is the coverage situation of each cell.
I(i) � 1, cell i has been visited, 0, cell i has not been visited.
e simulation results are shown in Figure 18. Here, it is observed that the sampling time does not significantly affect the performance of the system.
is 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. us, 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.
e curve fitting performed for the coverage percent showed an exponential convergence rate of the form where the values of the parameters are shown in Table 1.   0  50  100  150  200  250  300  350  400 450 500    10 Complexity Also, Figure 19 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. e fitting performed in Matlab indeed resulted in a linear curve.
where the values of the parameters are given in Table 2.
Overall, it is seen from the simulations that movement in eight directions yields better results compared to movement in four. is 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. at 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.

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. e 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][60][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][65][66].

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. e synchronization of the system was then considered using a simple and effective technique of active control. e 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. us, 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.

Data Availability
No data were used to support this study