Optimizing the Maximum Lyapunov Exponent of Fractional Order Chaotic Spherical System by Evolutionary Algorithms

: The main goal of this work is to optimize the chaotic behavior of a three-dimensional chaotic-spherical-attractor-generating fractional-order system and compare the results with its novel hyperchaotic counterpart. The fractional-order chaotic system is a smooth system perturbed with a hyperbolic tangent function. There are two major contributions in this investigation. First, the maximum Lyapunov exponent of the chaotic system was optimized by applying evolutionary algorithms, which are meta-heuristics search algorithms, namely, the differential evolution, particle swarm optimization, and invasive weed optimization. Each of the algorithms was populated with one hundred individuals, the maximum generation was ﬁve hundred, and the total number of design variables was eleven. The results show a massive increase of over 5000% in the value of the maximum Lyapunov exponent, thereby leading to an increase in the chaotic behavior of the system. Second, a hyperchaotic system of four dimensions was constructed from the inital chaotic system. The dynamics of the optimized chaotic and the new hyperchaotic systems were analyzed using phase portraits, time series, bifurcation diagrams, and Lyapunov exponent spectra. Finally, comparison between the optimized chaotic systems and the hyperchaotic states shows an evidence of more complexity, ergodicity, internal randomness, and unpredictability in the optimized systems than its hyperchaotic counterpart according to the analysis of their information entropies and prediction times.


Introduction
The subject of chaos and its application continues to attract research interests since it was made popular by Edward Norton Lorenz in 1963. His work on deterministic chaos [1] provided the foundation for chaos theory, which set the stage for enormous research focusing on the behavior of nonlinear dynamical systems, which display high sensitivity to initial values [2][3][4][5]. Basically, a nonlinear dynamical system has a state space, the coordinates of which show the state at any given time. The future states of any nonlinear dynamical system are crucial and they can be determined through a systematic approach. Beginning from an initial value, the entire future states of a nonlinear dynamical system, called trajectory or orbit, can be ascertained by iterating the system for a considerable number of times with a small time step advancement in each iteration.
Nowadays, chaos is increasingly becoming a subject of multidisciplinary research. Hence, it is being widely applied in security, communication, electronics, medicine, robotics, and so on. The dynamics of chaotic systems is crucial to the unpredictability property [6], and measuring the chaotic behavior is often done by computing the Lyapunov exponents (LEs) [7][8][9][10][11][12][13]. The LE measures the specific average rate of convergence or divergence of nearby orbits in the phase space. As a matter of fact, in a system that shows exponential optimization of the spherical-attractor-generating fractional-order system by EA, which is compared with a novel four-dimensional (4-D) hyperchaotic system modelled from the 3-D system. Therefore, the contributions arising from this investigation are highlighted as follows: (i) Optimizing the chaotic behavior of the 3-D spherical-attractor-generating system by maximizing the MLE with the application of DE, PSO and IWO metaheuristics. The results showed a significantly larger MLEs than the non-optimized system, which is evident in the phase space portraits, time series and the entropy of the optimized systems. (ii) Construction of a hyperchaotic system of 4-D. The hyperchaotic system was created with a state feedback controller added to the second equation of the original 3-D system. The analyses of the hyperchaotic system revealed that it possesses rich dynamics, exhibiting three different states, namely, hyperchaotic, chaotic, and periodic.
This article is organized into sections. Section 2 presents the 3-D chaotic dynamical system considered in this work. Section 3 contains the modelling of a new 4-D hyperchaotic system. In Section 4, we introduce the EA optimization algorithms applied to optimized the 3-D system. The procedure for evaluating the LEs is presented in Section 5. The results of this investigation, including the MLE optimization and comparison of the optimized systems with the novel hyperchaotic system, are shown in Section 6. In Section 7, we present the discussion of the results. Lastly, Section 8 contains the conclusion drawn from this investigation.

Chaotic Dynamical System Considered
This section describes the model and dynamical behavior of the 3-D fractional-order system that will be optimized in this investigation. The authors in [28] proposed a smooth, quadratic, and autonomous chaotic system with a hyperbolic tangent function. The chaotic system was constructed based on the Shilnikov criterion [29,30] and it can generate spherical attractors. The constructed chaotic system is presented in the next system of Equation (1), with its fractional-order counterpart given in (2).
x = a 1 x − a 2 y + a 3 z + 2 1 − e −200 sin y 1 + e −200 sin y ẏ = −dxz + b + eẋ z = c 1 xy + c 2 yz + c 3 z + c The original values of the parameters are a 1 = −4.1, a 2 = 1.2, a 3 = 13.45, b = 0.161, c = 3.5031, c 1 = 2.76, c 2 = 0.6, c 3 = 13.13, d = 1.6, and e = 0. The numerical evidences of the system have (x 0 , y 0 , z 0 ) = (−0.04, −15.8, −1.4) as the initial condition. The chaotic system has one equilibrium point, which is (x e , y e , z e ) = (0.7217, −2.5698, 0.1394), while the eigenvalues at this point are (λ 1 , λ 2 , λ 3 ) = (−0.3150, 3.9016 + 6.1032i, 3.9016 − 6.1032i). For the fractionalorder system (2), which is the focus of this investigation, the numerical integration was performed using optimized predictor-corrector Adams-Bashforth-Moulton (ABM) scheme with an integration step-size of h = 0.001. All numerical simulations were done in Matlab version 2016b. The phase diagrams of System (2) from the data of the last 100 seconds, to remove transient states, are presented in Figure 1. Applying a long simulation time of t = 2000 s to system (2), chaos is kept in the spherical attractor as seen in the phase portraits in Figure 1, in which we can see a permanent chaotic behavior in the system. In the next section, System (2) is considered for generating hyperchaos while Section 6 shows the optimization of its chaotic behavior via the MLE. In Figure 2 are found the bifurcation diagrams with the peak of state x against varying parameter e, together with its corresponding LE spectrum for System (2).

Novel Hyperchaotic System
Hyperchaotic systems usually possess more complex dynamics than the traditional chaotic systems. We intend to compare the 3-D fractional-order system, whose chaotic behavior will be optimized through the MLE in this work, with a fractional-order hyperchaotic system, to see how the optimized systems measure up against the hyperchaotic system. To achieve this, we present in this section a new 4-D non-linear system modelled from System (1). Presently, one popular method for constructing hyperchaotic systems from an existing 3-D chaotic system is to add state feedback controllers, linear or nonlinear, to the system. In general, the following two conditions must be satisfied to generate hyperchaos from a system: (1) The autonomous system should have a phase space of at least four dimensions, meaning that the number of the coupled first-order autonomous ordinary differential equations should be a minimum of four. (2) The coupled equations should have at least two terms responsible for instability, and a minimum of one should have a nonlinear function. In this work, we created the hyperchaotic system by including a periodic driving signal in the second equation of System (1), as shown in the next equation: where the values of a i , c i (1 ≤ i ≤ 3), d, b, c, and e remain the same as in the chaotic version, and r and g are parameters determining the dynamics of the hyperchaotic system. Hence, the controller w transformed the System (1) into a four-order autonomous system having four LEs. The fractional-order counterpart is modelled using the Caputo's derivative as follows:Ḋ We observed that the hyperchaotic system is invariant with respect to the following transformation: (x, y, z, w) → (x, y, z, −w) The flow means that the hyperchaotic system is symmetric with respect to x, y, and z planes. Therefore, we note that if (x, y, z, w) is a solution in the context that g is definite, there is a twin solution (x, y, z, −w) when g is negative. Furthermore, we define the vector field on the right-hand side of System (3) as follows: Consequently, the divergence of the vector field (6) is and the hyperchaotic system is dissipative if (a 1 + yc 2 + c 3 ) < 0. Therefore, The volume of the hyperchaotic system reduces to zero as t tends to infinity at an exponential rate of a 1 + yc 2 + c 3 .
To obtain the equilibrium points of the novel hyperchaotic system, we equate each of the system's equations to zero and evalute them to find the values of x, y, z, and w. Hence, at a 1 = −4.1, a 2 = 1.2, a 3 = 13.45, b = 0.161, c = 3.5031, c 1 = 2.76, c 2 = 0.6, c 3 = 13.13, d = 1.6, e = 0, and r = 15, the proposed hyperchaotic system has one equilibrium point: (x e , y e , z e , w e ) = (−0.3874, 0, −0.2668, 1.5705) (9) The eigenvalues at the equilibrium point P(x e , y e , z e , w e ) were obtained by evaluating the characteristic equation where J P is the Jacobian matrix for the equilibrium point P(x e , y e , z e , w e ), I is an identity matrix of the same size as J P , and λ is the eigenvalue. After some mathematical steps, the resulting characteristic equation is: Then, the eigenvalues are λ 1 = 0, λ 2 = 0, λ 3 = −4.1, and λ 4 = 13.13. Hence, the equilibrium point P(x e , y e , z e , w e ) is unstable, which is an indication of chaos in the proposed hyperchaotic system.

Optimization Algorithms
Enhancing the chaotic behavior of the fractional-order system (2) by optimization involves applying metaheuristics to explore the huge search space for the new system parameters and fractional-order that will increase the value of the MLE. In this work, three evolutionary algorithms, namely DE, PSO, and IWO, were applied to maximize the MLE of the chaotic system. Generally, the choice of DE, PSO, and IWO relies in the fact that they need only the value of the objective function to work, and it is not necessary to get any information about the derivative of the function.

Differential Evolution
The DE [31,32] is one of the most popular evolutionary algorithms that has been widely applied in systems optimization, especially chaotic systems. Some of the advantages of this algorithm are that it uses a few control parameters and has the ability to find the true global optimization of a multimodal search space, regardless of the initial parameter values. The DE uses a set of solution vectors N called individuals as a population in each generation G. Basically, the population of randomly created solution vectors is successfully improved by applying genetic operators in the following order: mutation, crossover and selection. There exists several variants of DE but the one used in this investigation is rand/1/bin, where the base vector at mutation are randomly selected, only one vector difference is used to create the mutant vector, and binomial recombination is used during crossover. The following equations in (12)- (14) represent the three genetic operators, where V i,G is the mutant vector, random indexes r1, r2, r3 {1, 2, ·, N}, and r1 = r2 = r3 = i, while F ∈ [0, 1] is the scaling factor controlling the amplification of the differential variation (X r2,G − X r3,G ).
whereby X i,G+1 represents the admitted vector into the next generation, while U i,G and X i,G are the trial vector and parent vector, respectively.

Particle Swarm Optimization
Just like the DE, the PSO is another popularly applied algorithm for systems optimization. Some of the benefits of the PSO are that it is a swarm intelligent algorithm that is based on a simple concept, and it is robust to control parameters and easy to implement. The PSO was inspired by the social behavior of birds and fish to solve optimization problems in a cooperative and intelligent framework using a simple mathematical model [33]. It uses a set of vectors N, which are potential solutions called particles, as a population for each generation G. The particles move through the problem space while following the latest optimal particles. In the problem space, each particle preserves its coordinates, which are linked to the best solution (fitness) it has ever achieved. The stored fitness value is often called pbest. In the event that a particle has the entire population as its topological neighbours, then the best value is a global best, which is often called gbest. Moreover, the velocity of each particle is changed at each time step toward its pbest and the best location, and it is weighted by a random term. However, separate random numbers are generated for acceleration toward pbest and the best location. The core of the PSO is the velocity and position update, represented by the next two models, respectively.
where k denotes the current generation, w is the inertia weight or constriction coefficient, c 1 and c 2 represent the cognitive and social coefficients, r 1 and r 2 are uniformly distributed random numbers in [0, 1], P bj is the personal best of a particle (coordinates of the best solution achieved so far by a particular particle), and P g represents the global best, which is the overall best solution obtained in the swarm at generation k.

Invasive Weed Optimization
The IWO has similar advantages as the PSO in that it is also a swarm intelligent algorithm with a simple concept, easy to implement, and robust to control parameters as well. It is a stochastic global optimization method inspired by the spreading and colonizing strategy of weeds [34]. The core of the IWO algorithm is reproduction and spatial dispersal of seeds. At the reproduction stage, seeds grow to become flowering plants and produce seeds depending on their fitness obtained from the objective function. Next is the spatial dispersal stage, during which the seeds produced are randomly dispersed over the search space. These operations are represented in the following models in (17) and (18).
where f i is the fitness of the i-th weed, f wst and f bst are the worst and the best fitness in the weed population, respectively. S min and S max represent the minimum and maximum number of seeds, respectively.
where itr max denotes the maximum iteration (generation), itr represents the current iteration, n is the nonlinear modulation index, and σ in and σ f n are the previously defined standard deviations.

Evaluation of the Lyapunov Exponents
The numerical evaluation of the LEs was by Benettin-Wolf algorithm, which employs the Continuous Gram-Schmidt orthogonalization (CGSO) procedure as seen in [35]. The procedure for evaluating the n Lyapunov exponents of an n-dimensional dynamical system requires finding the evolution of small perturbation to an orbit over a long time. The evaluation of the LEs is summarized in the following steps: while the initial condition of the variational system is set to I, an n × n identity matrix. (iii) The integration of the original and variational systems are done until the orthonormalization period K is reached. (iv) The variational system is then orthonormalized using the Continuous Gram-Schmidt orthogonalization.
(v) Next, the algorithm obtains and gathers in time the logarithm of the norm of each Lyapunov vector in the variational system. (vi) The next integration begins with the new orthonormalized vectors as the initial conditions. (vii) Steps (iii) to (vi) are repeated until the integration period T is reached. (viii) The n Lyapunov exponents are obtained by evaluating: where m = 1, 2, . . . , n, T is the integration period, K is the orthonormalization period, and ω i n are the orthonormal vectors. The largest of L 1 , L 2 , . . . , L n , is called the maximum Lyapunov exponent.
To maximize the MLE in this work, the optimization algorithms have to compute the Lyapunov exponents each time a likely optimal solution parameters are found, which might have a major impact on the running time.

Results
We present in this section the results of optimizing the chaotic behavior of the fractional-order system (2) by the evolutionary algorithms introduced in Section 4. In addition, analyses of the dynamic behavior of the hyperchaotic system (4) are presented and both the optimized systems and hyperchaotic system are compared.

MLE Optimization
The DE, PSO, and IWO algorithms were given the same conditions: randomly populated with 100 individuals, maximum generation of 500, and five individual optimization runs were performed. We considered two factors in choosing the search spaces for the design variables. First, we decided that the values of the parameters should not be discrete but continuous, real numbers of four decimal places in our work, that will vary continuously with no gaps within a specified range. Second, defining the lower and upper bounds of the search space for each design variable was based on the parametric studies through initial bifurcation analysis with each of the parameters. In doing this, we ensured that the search spaces are not overly restricted. The following are the search spaces of the parameters: −6.0000 ≤ a 1 ≤ 0.0000, −10.0000 ≤ a 2 ≤ 30.0000, 0.0000 ≤ a 3 ≤ 50.0000, 0.0000 ≤ b ≤ 2.0000, 0.0000 ≤ c ≤ 15.0000, 0.0000 ≤ c 1 ≤ 20.0000, 0.0000 ≤ c 2 ≤ 10.0000, 0.0000 ≤ c 3 ≤ 20.0000, 0.0000 ≤ d ≤ 5.0000, −12.0000 ≤ e ≤ 2.0000, and 0.0001 < q ≤ 1.0000. The specifications of the computer used for the optimization and the specific parameters of each algorithm are given next: In this optimization, the best results, in terms of the MLEs from the optimized parameters and fractional order, were obtained in 500 generations for all the three algorithms. Therefore, any mention of the DE, PSO, and IWO-optimized fractional-order systems in the rest of this paper refers to the best results obtained in the respective algorithms. The optimized MLEs in the five runs for each algorithm were higher than the non-optimized systems, with increase in the MLE having a corresponding increase in the chaotic behavior. In Table 1, we present the global optimized parameters for the highest MLE in the five runs from each algorithm. As seen in the Table, the MLEs in all the cases are higher than the non-optimized system. The effectiveness of the applied EAs was noticed in the variance among the MLEs in the five runs, as the maximum variation was approximately ±0.05. Furthermore, Table 1 shows the equilibrium points, eigenvalues, information entropies, and the instability analysis results. The instability analysis, examined from the eigenvalues, is determined according to the instability theorem in [19]. To get the equilibrium point, we equate each of the three equations of the chaotic system to zero and by simplification and substitutions, we arrived at the following equations to compute z, y, and x, respectively: and where (22) gives one real root for x. Hence, the chaotic system has one equilibrium point, denoted by P(x e , y e , z e ). To compute the eigenvalues, the characteristic polynomial of the Jacobian matrix at the equilibrium point P(x e , y e , z e ) is in the form: where p 1 = 1.3804, p 2 = −719.8723 and p 3 = −4549.8224 for DE-optimized system, p 1 = 1.3593, p 2 = −719.5002 and p 3 = −4611.0113 for PSO-optimized system, and p 1 = 1.3586, p 2 = −752.0069 and p 3 = −4907.0811 for IWO-optimized system. The information entropy H, which describes the degree of uncertainty in an information source, is defined as: where s is the information source and P(s i ) = Pr(s = s i ) is the probability of the ith value of s. We note that the LEs of a chaotic system can change if there is a slight change in the value of a parameter, the fractional order, or initial condition. For example, the MLE of the non-optimized fractional-order system computed in this work when q = 0.9999 is almost the same as that of the integer order when q = 1 using the same parameters. The slight difference of 0.0024 is caused by the change in the q value. In Figure 3, the evolution of the phase diagrams in 3-D planes obtained from the MLEs of the DE-optimized, PSOoptimized, and IWO-optimized systems are presented, simulated with a time-step h = 0.001 and integration period T = 2000 s. The increased chaotic behavior in the optimized systems can be observed in the phase diagrams. Moreover, the improved chaoticity in the optimized systems can also be observed in the time series presented in Figure 4. It is noted that the usual spherical shape seen in the non-optimized system has been greatly altered. This is due to the increased chaos in the optimized systems. Furthermore, the dynamics of the non-optimized and optimized systems was examined by plotting their bifurcation diagrams with the peak of state x against varying parameters a 1 , a 2 , c 1 , d, and fractional-order q, as shown in Figures 5-8. The LE spectra have a good agreement with the bifurcation diagrams.
The maximization of the MLE in this investigation further shows the effectiveness of the EAs as global optimization algorithm [6,7,15,16,19,36,37]. Table 1 reveals that the three algorithms obtained system parameter and fractional order values that produced better MLEs than the non-optimized fractional-order system. The higher MLEs obtained is an indication that the optimized systems have more unpredictability than the non-optimized system. (c) Figure 3. Phase diagrams of the optimized fractional-order chaotic system (2) as presented in Table 1:    Figure 5. Bifurcation diagrams and Lyapunov exponents spectra of the optimized fractional-order system (2), using parameters a 1 , a 2 , c 1 , d and fractional-order q, for the non-optimized system.  Figure 6. Bifurcation diagrams and Lyapunov exponents spectra of the optimized fractional-order system (2), using parameters a 1 , a 2 , c 1 , d and fractional-order q, for the DE-optimized system.  Figure 7. Bifurcation diagrams and Lyapunov exponents spectra of the optimized fractional-order system (2), using parameters a 1 , a 2 , c 1 , d and fractional-order q, for the PSO-optimized system.  Figure 8. Bifurcation diagrams and Lyapunov exponents spectra of the optimized fractional-order system (2), using parameters a 1 , a 2 , c 1 , d and fractional-order q, for the IWO-optimized system.

Optimized Systems against Hyperchaotic System
In this segment, we compare the optimized systems with the new hyperchaotic fractional-order system (4). We begin by examining the dynamics of the hyperchaotic system, considering the bifurcation diagram, LE spectrum, and the phase portraits. The initial value for simulating the hyperchaotic system is (x 0 , y 0 , z 0 , w 0 ) = (−0.04, −15.8, −1.4, −0.12), integration step time h = 0.001, and integration time interval I = [0, 1000]. Unlike the parameters of the original 3-D system, parameters r and g are responsible for the dynamic behavior of the new hyperchaotic system. Hence, they were selected as the bifurcation parameters to observe their effect on the new system.
For the purpose of identification, L 1 , L 2 , L 3 , and L 4 denote the four LEs of System (4). The peak of state x and the LEs were separately plotted against the varied value of parameter −10 ≤ r ≤ 20 and −10 ≤ g ≤ 20 to generate bifurcation diagrams and LE spectra. In Figure 9a is found the bifurcation diagram and the corresponding LE spectrum when g = 5 and r varies. In addition, Figure 9b contains the bifurcation diagram and LE spectrum, respectively, when r = 15 and g varies, while Table 2 shows selected LEs. It is observed in the table that the system exhibits three different behaviors with changing values of parameter r and g, namely hyperchaotic, chaotic, and periodic. The LE spectra revealed that the hyperchaotic state has two positive, one negative, and one zero LEs with L 1 , L 2 > 0, L 3 = 0 and L 4 < 0, chaotic state has one positive, two negative, and one zero LEs with L 1 > 0, L 2 , L 4 < 0 and L 3 = 0, and the periodic state has three negative and one zero LEs with L 1 , L 2 , L 4 < 0 and L 3 = 0. (b) Figure 9. Dynamics of fractional-order hyperchaotic system (4) showing the bifurcation diagram (left) and Lyapunov exponent spectrum (right), when (a) g = 5 and r varies; (b) r = 15 and g varies. Table 2. Selected Lyapunov exponents for the hyperchaotic fractional-order system of (4) when (i) g = 5 was fixed and r varies, and (ii) when r = 15 was fixed and g varies. The Kolmogorov-Sinai entropy h ks , which is the sum of the positive LEs, was computed. The highest h ks of the hyperchaotic state when r varies, obtained at r = 15.25, is shown in (25), while for when g varies, obtained at g = 5, the highest h ks is in (26).
To further validate the dynamics of the new hyperchaotic system, the geometric illustration of the trajectories in the phase plane were plotted and projected on both 3-D and 2-D planes using the data of the final 100 s to remove transient states. The various projections of the hyperchaotic state when r = 15.25 and g = 5 are presented in Figure 10. We noted in our investigation that the integer-order hyperchaotic system (3) and its fractional-order counterpart produced similar dynamics shown in the bifurcation and LE spectra plots when both parameters r and g were varied. In this work, we applied information entropy and prediction time as the basis for comparing the optimized chaotic systems and the hyperchaotic system. Particularly in a chaotic dynamical system, the information entropy, defined in (24), quantifies the ergodicity and randomness inherent in the system, and the larger the information entropy, the better the chaos in the system. On the other hand, the prediction time can be used to measure the security strength of dynamical systems, depending on the positive Lyapunov exponents. Generally, a dynamical system with shorter prediction time provides higher security for encryption algorithms [38]. In this work, the prediction times of the non-optimized and optimized chaotic systems and the three hyperchaotic states were computed. The prediction time is defined as: where K 1 is the Kolmogorov-Sinai entropy. Table 3 shows the information entropies for state x(t) of the non-optimized and optimized systems in Table 1, juxtaposed with that of the hyperchaotic states from the new hyperchaotic system (4) that has the best three Kolmogorov-Sinai entropies. They are named hyperchaotic 1, 2, and 3. In addition, the table presents the prediction times of the systems. In the table, it can be seen that the information entropies of the optimized fractional-order systems are higher than those of the hyperchaotic states. This is an indication of more complexity, ergodicity, and internal randomness in the optimized systems than the hyperchaotic states. On the other hand, the optimized chaotic systems have shorter prediction times than the non-optimized chaotic and the hyperchaotic systems. The prediction times of the hyperchaotic system are only shorter than the non-optimized chaotic system. Specifically, DE-optimized system has the largest entropy and the shortest prediction time. Table 3. Information entropies and prediction times of the non-optimized and optimized chaotic systems, and the best three hyperchaotic states of the new hyperchaotic system.

System
Entropy H(s)

Discussion
In the investigation on optimizing the chaotic behavior of the fractional-order spherical system, the best MLEs were obtained with optimized fractional orders 0.9400 < q < 0.9420 (see Table 1). Although the DE, PSO and IWO algorithms produced almost the same optimized MLEs, the DE has overall best MLE of 1.0808, which is a whopping 5, 806% increase over the non-optimized system. The structure of the strange attractors obtained from the optimized systems (see Figure 3) is a clear indication of the increased irregularity and non-periodicity in the systems [39,40], as a result of the optimization of the MLE. It is observed that the values of parameters a 3 , c 1 and e are the same in the three optimized results. The maximization of the MLE further affirms the EAs as an effective optimization algorithm. The DE, PSO, and IWO algorithms obtained system parameter and fractional order values that produced higher MLE than the non-optimized chaotic system.
At this juncture, it is relevant to note that the global optimization by IWO is the fastest in terms of the average running time. The reason that can be adduced to this is that IWO performed the least computational activities than DE and PSO. Overall, DE required 163 h, PSO 167 h, and IWO 144 h, to complete an optimization run of 500 generations. The asymptotic analysis results presented in Table 1 shows that the equilibrium points of the DE, PSO and IWO-optimized systems are saddle points, just like the non-optimized system. Specifically, the equilibrium points of the DE, PSO and IWO-optimized systems are saddle points of index 2. From the instability analysis, we can observe that the DE-optimized system can be asymptotically stable for fractional order q < 0.9150, PSO-optimized system for q < 0.9145, and IWO-optimized system for q < 0.9150 as well.
The computed information entropy to measure the complexity of the systems shows consistency with the optimized MLEs, having higher entropy than the three best hyperchaotic states, with the DE-optimized system producing the best value of 6.9334 while the non-optimized has the lowest with 0.9772 (see Table 3). This means that the DE-optimized system exhibits the most randomness in it. Consequently, the optimized systems are more complex and have more internal randomness and greater unpredictability than the non-optimized 3-D system and the hyperchaotic system, which is beneficial to encryption systems and random number generators. Moreover, the shorter prediction times in the optimized systems, 0.6413, 0.6433 and 0.6501, respectively, for the DE, PSO and IWO-optimized systems, resulted from the increase in the MLE (see Table 3). A system that has a low prediction time is regarded as safer for designing security systems [38]. Hence, a better secure communication system is guaranteed in the optimized fractional-order systems.
In the information provided in Table 4, we can see the comparison of our investigation on the chaotic behavior optimization with some other works in [6,7,15,16,36,37], which were performed on global optimization of dynamical behavior and parameter estimation in chaotic systems. It can be seen that we used more algorithms than others to determine their performance for the optimization problem. In our investigation, the fractional order was part of the design variables, making a total of eleven. Furthermore, this investigation used a wider search space for each of the design variables. Therefore, the considerably high MLEs obtained for the chaotic spherical attractor can be mostly accredited to the higher population size, higher number of generations, and wider search space. In addition, unlike the compared works, we studied the bifurcation and LE spectra of the optimized chaotic systems, varying the parameters and the fractional order, and both were in agreement with each other for each parameter and the fractional order. It should be noted that the sum of the LEs of the hyperchaotic systems as the parameters were varied is negative, meaning that the hyperchaotic system is dissipative [41]. As a result, the phase space volume contracts along a trajectory with a strange attractor formed, such as the hyperchaotic attractors projected in Figure 10. It can be seen in the hyperchaotic attractors that the spherical shape observed in chaotic systems (1) and (2) has been replaced with spiral-like shape of multiple scrolls. Hence, there is an evidence of hyperchaos in the new systems.
Finally, the optimization results provided in Table 1 are very dependent not only on the time averaging, but also on the numerical method that was used to solve the dynamical systems. However, the length of the time simulation is enough to observe the attractors in the phase space and for this reason, we can conclude on the appropriateness of applying evolutionary algorithms for the optimization of chaotic systems through the MLE. We note that Matlab is time-consuming, but this can be mitigated by programming the algorithms in C or Python language. In addition, one can save more computing time if the numerical method is also improved by paying attention to the selection of the highest time-step, which can be estimated according to the numerical method and eigenvalues. Another important point is that the optimized results can be used to design integrated chaotic systems, as shown in [2].

Conclusions
This work considers a 3-D fractional-order nonlinear, smooth, autonomous system, having a hyperbolic tangent function and capable of generating spherical chaotic attractors. First and foremost, the investigation shows the MLE maximization of the fractional-order chaotic system using three evolutionary algorithms, namely, differential evolution, particle swarm optimization, and invasive weed optimization. To achieve this, a total of eleven design variables, including ten parameters and one fractional order, were optimized. The result of the optimization shows improved chaotic behavior resulting from over 5000 percent increase in the MLE value, compared to the non-optimized system. In addition, the optimization shows that the chaotic behavior of the fractional-order system considered in this work is multifaceted in relation to the parameters and fractional-order. The main advantage of using the DE, PSO, and IWO for this work is that they only need the value of the objective function to work, and it is not necessary to get any information about the derivative of the function.
Furthermore, we modelled a new 4-D hyperchaotic system from the initial 3-D system and examined its dynamical behavior. Based on the four Lyapunov exponents, the hyperchaotic system has multistates, namely, periodic, chaotic, and hyperchaotic. The DE, PSO and IWO-optimized systems were compared with the non-optimized chaotic system and the three best hyperchaotic states using information entropy and prediction time. The results show that the DE-optimized system has the best MLE, information entropy, and prediction time. Therefore, we can conclude that the DE-optimized system has the best unpredictability and internal randomness. This investigation offers two major contributions. First and to the best of our knowledge, optimization of the chaotic behavior of the considered 3-D dynamical system was performed for the first time. The MLE was maximized, which increases the chaotic behavior of the system, evidenced in the phase space portraits, time series, and the information entropy analysis. Second, a novel hyperchaotic system was created using a state feedback controller. The new hyperchaotic system possesses rich dynamics and it is symmetric and dissipative. Lastly, based on the contributions highlighted above, the significance of this investigation is that the optimized systems and the hyperchaotic system are promising tools for developing cryptosystems, secure communication systems, and random number generators.

Funding:
The authors wish to thank the Instituto Politecnico Nacional, Secretaría de Investigación y Posgrado, for its support provided through the project SIP 20220014.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Particle Swarm Optimization