Maximizing the Chaotic Behavior of Fractional Order Chen System by Evolutionary Algorithms

This paper presents the application of three optimization algorithms to increase the chaotic behavior of the fractional order chaotic Chen system. This is achieved by optimizing the maximum Lyapunov exponent (MLE). The applied optimization techniques are evolutionary algorithms (EAs), namely: differential evolution (DE), particle swarm optimization (PSO), and invasive weed optimization (IWO). In each algorithm, the optimization process is performed using 100 individuals and generations from 50 to 500, with a step of 50, which makes a total of ten independent runs. The results show that the optimized fractional order chaotic Chen systems have higher maximum Lyapunov exponents than the non-optimized system, with the DE giving the highest MLE. Additionally, the results indicate that the chaotic behavior of the fractional order Chen system is multifaceted with respect to the parameter and fractional order values. The dynamical behavior and complexity of the optimized systems are verified using properties, such as bifurcation, LE spectrum, equilibrium point, eigenvalue, and sample entropy. Moreover, the optimized systems are compared with a hyper-chaotic Chen system on the basis of their prediction times. The results show that the optimized systems have a shorter prediction time than the hyper-chaotic system. The optimized results are suitable for developing a secure communication system and a random number generator. Finally, the Halstead parameters measure the complexity of the three optimization algorithms that were implemented in MATLAB. The results reveal that the invasive weed optimization has the simplest implementation.


Introduction
The nonlinear dynamical system is a branch in science, especially physics and mathematics, which is fundamental in the modelling of physical systems. Because of groudbreaking researches by pinoneering scientists in nonlinear dynamical system, a concept called chaos theory has emerged and, in the last few decades, the subject of chaos has attracted great interest from researchers. Initial studies on chaos have shown that the phenomenon has some unique properties in terms of sensitivity to initial conditions, complexity, ergodicity, transitivity, and determinism, which have been exploited in crucial areas, such as telecommunication [1][2][3], robotics [4][5][6], medicine [7][8][9], and so on. Therefore, investigations have led to the development of many chaotic systems, and there are several activities in areas, such as electronic realization, stability, synchronization, and chaos-based communication [10][11][12]. discontinuous fractional order Shimizu-Morioka system, while the dynamical behavior of a new hyperchaotic fractional order Rabinovich system was studied in [33].
Moreover, some more recent works on fractional order chaotic systems, their applications, and implementations are contained in [34][35][36][37][38][39][40]. The authors in [34] presented an analog approach towards implementing the fractional order derivates of a chaotic system, utilizing amplifiers and other passive circuit elements. In addition, the analog method includes the use of field-programmable analog arrays (FPAAs) to minimize the mismatch that may arise when using discrete devices. The work presented in [35] shows the application of a fractional order chaotic oscillator for secure communication systems based on the synchronization of two fractional order chaotic systems in a master-slave configuration. The secure communication system was applied to encrypt RGB and grayscale images. The work further shows fractional order chaotic systems as reliable for implementing random number generators for cryptographic applications. The authors presented the experimental field-programmable gate arrays (FPGAs) realization of four different fractional order chaotic systems in [36]. In the work, the Grünwald-Letnikov (GL) method was applied to solve the fractional order derivatives, while, in [37], the authors presented a work on numerical analysis of the Chen oscillator, whereby the computation of the LEs was a part of the results. The authors introduced an encryption, compression, and transmission scheme in [38], based on a fractional order chaotic system in conjunction with discrete wavelet transform (DWT) and quadrature phase shift keying (QPSK). The encrypted image file was wirelessly transmitted through a very small distance (7 cm). A novel multistable fractional order chaotic system, having an unstable equilibrium point of self-excited attractor and hidden attractors with no equilibrium point, was proposed in [39] with its digital signal processor (DSP) implementation, while the synchronization of two variable fractional orders chaotic systems using fuzzy modeling were considered in [40] with application in secure communications. In [41], the authors included a work on the optimization of the positive LE and Kaplan-Yorke Dimension of fractional order chaotic systems by differential evolution and particle swarm optimization. Here, the fractional order was fixed, while the system parameter values were varied within the search space, looking for the combination of values that increase MLE. This is probably the first attempt in this regard, because a rigorous search of the literature did not yield other works.
In the present investigation, the fractional order Chen system is considered with the primary objective of maximizing its MLE using EAs [42][43][44]. This effort will add to the pool of literatures that are available in the research area of optimizing the MLEs of fractional order chaotic systems. Specifically, the global optimization procedure presented in this work employs meta-heuristics to search for the values of the system parameters and fractional order that maximize the MLE. The dynamic properties and complexity of the optimized fractional order Chen chaotic system are theoretically and numerically verified. Bifurcation, LE spectrum, equilibrium point, eigenvalue, and sample entropy are the properties considered. The behavior of a system can change dramatically with a change in the values of its parameters. This is called bifurcation. The parameter values at which these changes occur are called bifurcation points. The spectrum of the LEs is a set of real numbers that are sorted in non-increasing order and equal to the dimension of the phase space. Knowledge of the complete LE spectrum gives more detailed analysis of a dynamical system. Geometrically, the equilibria are points in a system's phase space that are calculated by solving the function f (x) = 0. The stability of a dynamical system can be verified through the evaluation of the eigenvalues of the system. For instance, a saddle point is an equilibrium point in which at least one eigenvalue is located in the stable region and one eigenvalue is in the unstable region. An index 1 saddle point is the one in which only one eigenvalue is unstable and the others are stable. A saddle point of index 2 has one stable eigenvalue and two unstable eigenvalues.
The need to develop a more unpredictable and reliable fractional order chaotic Chen system for secure communication systems is the motivation behind this work [45][46][47][48][49]. The global optimization in this work is carried out using three single objective EAs, namely: differential evolution (DE), particle swarm optimization (PSO), and invasive weed optimization (IWO), and the results are then compared. In the light of this, the following contributions to the state of the art are hereby highlighted: (i) the optimization of the MLE of the fractional order chaotic Chen system by DE, PSO, and IWO, whereby the system's fractional order q is not fixed, but it is considered to be a design variable and optimized alongside the conventional design variables, which are the system parameters. This was done because a slight change in the value of the fractional order q can also lead to a big change in the LEs; and, (ii) optimizing the system parameters and fractional order q gives up to an 80% increase in the value of the MLE over the non-optimized system. The highest optimized MLE is obtained from the DE optimization. Consequently, the optimized fractional order chaotic Chen systems are more complex and unpredictable, which makes them suitable for developing random number generators and secure communication systems for cryptographic applications.
This article is organized, as follows: Section 2 presents the theoretical framework, describing the fractional order chaotic Chen oscillator, the computation of the LEs, the applied EAs, and the complexity and instability analyses. Section 3 shows the results that were obtained in this work. The discussion of the results can be found in Section 4. Finally, the conclusion that is derived from the investigation is given in Section 5.

Theoretical Framework
This section describes some concepts and techniques that are relevant to this investigation. The following concepts are analyzed: fractional order chaotic Chen oscillator, LEs, the EAs under consideration, and time series complexity and instability of nonlinear systems.
In modeling fractional order systems, there are several definitions for the fractional derivative of order q > 0, but the Caputo's definition is the most commonly used. Caputo's definition includes initial conditions for the function f as well as its integer order derivatives. Therefore, using the Caputo's derivative, the initial value problem (IVP) of a fractional order system is defined, as follows: where D q * is the Caputo's differential operator of order 0 < q ≤ 1 and q = [q 1 , q 2 , q 3 , . . . , q n ] T , f : R n → R n is a Lipschitz continuous nonlinear function, x 0 ∈ R n is the initial condition, and T > 0. If q 1 = q 2 = q 3 = . . . = q n , then system (1) is called to have commensurate order; otherwise, non-commensurate order. D q * is defined by: where m has the smallest integer that is larger than q, x (m) (t) is the m-order derivative in the usual sense, and Γ is the Euler's Gamma function. Numerical approximation methods, such as the Grünwald-Letnikov, Riemmann-Liouville, predictor-corrector Adams-Bashforth-Moulton (ABM), and optimized predictor-corrector ABM can be applied to evaluate the fractional differential equation in (1).

Fractional Order Chaotic Chen Oscillator
The Chen oscillator is a continuous nonlinear and autonomous dynamical system, and it is an initial value problem, either in integer or fractional order. Guanrong Chen and Tetsushi Ueta proposed the oscillator [50,51], and its fractional order version is modeled by Caputo's fractional derivative, as follows: where D q * is the Caputo's differential operator of order 0 < q ≤ 1; x, y, and z are the system dependent variables; and, a, b, and c are the system parameters with traditional values of 35, 3, and 28, respectively. We considered commensurate order in our work, hence, q = q 1 = q 2 = q 3 . The fractional order Chen system (3) has the same equilibrium points as its equivalent integer order. Hence, the three equilibria are [37]: The Jacobian matrix of the fractional order Chen system (3) is:

Computation of Lyapunov Exponents
The LEs of continuous dynamical systems provide the most quantitative description of the presence of a deterministic non-periodic flow. The computation begins with determining their variational equation.

Definition 1.
Given a three-dimensional continuous fractional order system of the form: where D q * is the Caputo's differential operator. If e 0 = (x 0 , y 0 , z 0 ) is the initial condition, a small perturbation to e 0 in x, y, and z directions will lead to the evolution of the initial perturbed condition vector in the direction of a different point, say, x , y , and z , respectively. The slopes of system (8) in each direction, which are described by the Jacobian matrix J, provide an indication of the evolution of the perturbation after a finite time t. The variation or perturbation in the three directions is defined by: Generally, the variational equation in (9) can be written in the form: where β represents the matrix solution of system (8), β ∈ R n×n , D v is the n × n Jacobian of f , and I is the n × n identity matrix. Because β(0) is non-singular, β(t) is non-singular for t ≥ 0. The MLE is computed following the algorithm that is described in [45]. The methodology is based on the work shown in [52,53], in which a Gram-Schmidt orthogonalization procedure for computing LEs was proposed.

Definition 2.
Given j linearly independent set of vectors {α 1 , α 2 , . . . , α j } in R n , the Gram-Schmidt orthogonalization procedure generates an orthonormal set {ζ 1 , ζ 2 , . . . , ζ j } of vectors spanning the same subspace as {α 1 , α 2 , . . . , α j }. The volume of the parallelepiped that is covered by {α 1 , α 2 , . . . , α j } is: The variational Equation (10) is integrated from {e 0 , Z 0 }, until a period T is reached, to obtain: The corresponding matrix of orthonormal vectors Z 1 is calculated, the variational equation is integrated, and e 2 and A 2 are computed. The integration and orthonormalization continue until K is reached. A, Z, and e are the n × n matrix, orthonormal vectors matrix, and initial condition, respectively. At k-th step of the process, the set of orthonormal vectors is {ω k 1 , ω k 2 , . . . , ω k j }. For a suitable value of T, the LEs are obtained by evaluating: until a maximum iteration is attained or when convergence occurs.

Description of the Evolutionary Algorithms
In optimization, the function that is to be optimized is generally of the form: g(X) : R D → R. The global optimization problem of maximizing (or minimizing) the objective function g(X) is by optimizing the values of parameters X = {x 1 , x 2 , . . . , x D }, X ∈ R D , where X is a vector comprising D objective function parameters. The parameters are bounded by lower and upper boundary constraints, i.e., . Most of the EAs use mechanisms that are inspired by biological evolution, such as reproduction, mutation, recombination, and selection. The use of meta-heuristics search is a norm in global optimization problems because of the huge search space that is usually involved. The optimization algorithms that are applied in this investigation, DE, PSO, and IWO, are evolutionary and based on a meta-heuristic. An initialization of population is the first of the main activities of every EA. The others are random variation of individuals, evaluation of cost function, and selection. The last three activities are normally conducted in an iterative manner. Equation (15) shows the initialization of population of N individuals and V decision variables, where l and u are the lower and upper bounds of the decision variables, respectively, i = 1, 2, . . . , N and j = 1, 2, . . . , V.

Differential Evolution
The DE was introduced and described in detail by the authors in [54,55]. It uses a set of vectors N called individuals as a population for each generation G. Each vector represents the potential solution for the global optimization problem. Currently, there exist several variants of DE, but the one used in this investigation is rand/1/bin. DE uses the difference of solution vectors to create new candidate solutions, whereby the distribution of the population and its orientation are hidden in the differences in population members. As the population size increases, the distribution behind the DE sampling technique tends to a multi-variable Gaussian (normal) distribution, which is very common in the context of EAs. Specifically, the population of randomly created N solution vectors is successfully improved by applying genetic operators in the following order: mutation, crossover, and, finally, selection. Mutation is the main novelty of DE. This operation enables DE to explore the search space and maintain diversity. The mutation is explained in the following definition.
Definition 3. For a given parameter vector X i,G , i.e., an individual i and a generation G, three vectors (X r1,G X r2,G X r3,G ) are randomly selected, and a mutant vector V i,G is then created, such that: where random indexes r1, r2, r3 {1, 2, . . . , N}, r1 = r2 = r3 = i, and F ∈ [0, 1] is the scaling factor that controls the amplification of the differential variation (X r2,G − X r3,G ).
The crossover operator is applied to the parent vector X i,G and mutant vector V i,G , in order to increase the diversity of the perturbed parameter vectors. Two crossover techniques are available, namely, binomial and exponential crossover. In binomial crossover, which is used in this investigation, the trial vector U i,j,G is formed, as follows: where ] is the crossover probability, and D is the number of decision variables. jrand is randomly chosen in order to ensure that U i,G obtains a minimum of one parameter from V i,G . The selection operator shown in Equation (18), for maximization problem, works by comparing the parent vector with the trial vector. The vector with the better fitness value is admitted to the next generation.

Particle Swarm Optimization
The PSO was developed by James Kennedy and Russel Eberhart, and described in [56,57]. PSO was inspired by the social behavior of birds and fish to solve optimization problems in a cooperative and intelligent framework. Unlike other EAs, PSO has no evolution operators, such as crossover and mutation. PSO uses a set of vectors N (that are called particles) as a population for each generation G. Potential solutions (particles) fly through the problem space by following the current optimum particles. Each particle keeps track of its coordinates in the problem space that are associated with the best solution (fitness) that it has achieved so far. The stored fitness value is called pbest. When a particle takes all the population as its topological neighbours, the best value is a global best and it is called gbest. The velocity of each particle, weighted by a random term, changes at each time step toward its pbest and best location. The core of the PSO is the particle's velocity and position update. Definition 4. Given a swarm of N dimensional population and D dimensional particles, the position and velocity of the j-th particle are represented by x j = x j1 , x j2 , . . . , x jD and v j = v j1 , v j2 , . . . , v jD , where j = 1, 2, . . . , N, the following equations give the velocity update and position update of the particle, respectively: v j,k+1 = w k v j,k + c 1 r 1 (P bj,k − x j,k ) + c 2 r 2 (P g,k − x j,k ) (19) x where k is the current generation, w is the inertia weight, c 1 and c 2 are 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 (the coordinates of the best solution obtained so far by a specific particle), and P g represents the global best, which is the overall best solution that was obtained in the swarm at generation k.

Invasive Weed Optimization
IWO is inspired by the spreading and colonizing strategy of weeds and it is described in [58]. Weeds represent the feasible solutions for the problem and population is the set of all weeds. A finite number of seeds (initial solutions) are randomly spread over the search area. This is the initialization stage. The reproduction stage follows, whereby the seeds grow to become flowering plants and produce seeds depending on their fitness obtained from the objective function [50]. Equation (21) represents the number of seeds N s to be produced by the plants, where f i is the fitness of the i-th weed, f wst and f bst are the worst and best fitness in the weed population, and S min and S max represent the minimum and maximum number of seeds. The spatial dispersal stage is next, during which the seeds produced are randomly dispersed over the search space by normally distributed random numbers with a mean equal to zero, but with a varying standard deviation σ itr . This is achieved by the following Equation (22): where itr max is the maximum number of iterations (generations), itr is the current iteration, σ in and σ f n are the previously defined standard deviations, and n is the nonlinear modulation index. When the maximum number of plants is reached, the competitive exclusion technique is applied, whereby only the plants with high fitness can survive and produce seeds, and others are eliminated. The process continues until the maximum number of iterations is reached. The plant with the best fitness value is the optimized solution.

Complexity Analysis and Instability of Equilibria
The field of information theory studies the complexity in data series by analyzing its entropy. The entropy is a quantity that measures the degree of randomness in a system or data series to ascertain not only whether it is random or not, but also how random it is. As the level of randomness in a system increases, so does its complexity and unpredictability. There are several algorithms for measuring the randomness that exist in data series, chaotic time series inclusive. Examples of such algorithms incude approximate entropy, sample entropy, and permutation entropy. In this work, sample entropy [59,60] is applied to measure the time series complexity of the fractional order chaotic Chen system. Sample entropy measures the probability of generating new patterns in signals. The sample entropy is computed, as follows: where m is the embedded dimension, r is the tolerance, N is the time series data, and A and B are template vector pairs.
The instability of an autonomous system, as in (1), is examined from the eigenvalues λ 1 , λ 2 , . . . , λ n of the Jacobian matrix J = ∂ f ∂x evaluated at each equilibrium point according to the following theorem [61]: Theorem 1. The equilibrium points of the autonomous fractional order system (1) are asymptotically stable if and only if σ = q − 2| arg(λ) min |/π (24) is strictly negative.

Remark 1.
If σ ≤ 0 and the critical eigenvalues satisfying τ = 0 have a geometric multiplicity of one, then the equilibrium point is stable.

Remark 2.
If σ > 0, then the equilibrium point is unstable and the system may exhibit chaotic behavior.
As a matter of fact, the system parameters that are used in the computation of the LE are the primary component for computing the eigenvalues upon which Theorem 1 is based. It is also safe to deduce that, if σ > 0, then there may be a positive LE and chaotic attractor.

Results
The DE and PSO algorithms were randomly populated with 100 individuals, while IWO has an initial population size of 10 and maximum of 100. The number of generations was arbitrarily set from 50 to 500 with a step of 50. The values of the system parameters a, b, and c, as well as fractional order q in the search space, were kept within four decimal places. The initial condition is chosen in the basin of attraction of the system. For the step-size value, some tests were carried out on the simulation program, starting with a trial value and continuing to reduce the value until convergence occurs. The underlying numerical method adopted is the fast optimized ABM predictor-corrector [62]. The numerical evaluation of the LEs is by the Continuous Gram-Schmidt orthogonalization (CGSO) procedure, as justified in [62]. The configuration of the computer that was employed in the simulations and the specific parameters of each algorithm are given next:  Table 1 presents the optimized parameters and fractional order by the DE, PSO, and IWO, as well as the optimized LEs, against the non-optimized system. In addition, the table shows the equilibrium points, eigenvalues, sample entropy, and the asymptotic results applying Theorem 1. We note that the MLE of the non-optimized system that is computed in this work when q = 0.9800 is almost the same as the 2.0192 reported in [63] for integer order when q = 1 using the same parameters. It can be seen that the slight difference of 0.0101 is caused by the q value. Generally, the LEs can be altered if there is a change in the value of a parameter or the fractional order.  Figure 1 presents the evolution of the phase diagram in 3D planes of the best MLE obtained from DE, PSO, and IWO, respectively, which were simulated with a time-step h = 0.001 and integration period T = 75 s. Additionally, the chaos in the optimized system by each EA is compared in the time series presented in Figure 2.  Furthermore, the dynamics of the non-optimized and optimized systems are examined with bifurcation diagrams of state x as well as the LE spectra against the varied values of parameters a, b, and c, and fractional order q. Figures 3-6 display the results. The corresponding LE spectra in the figures show good agreement with the bifurcation diagrams.

Comparison with Hyper-Chaotic Fractional Order System
The optimized parameters of Chen system that were obtained in this investigation are compared with a hyper-chaotic Chen system on the basis of their prediction times [64]. The hyper-chaotic system is a class of dynamical system with more than one positive Lyapunov exponent. They are usually more complicated than the traditional chaotic systems. The integer order of the compared hyper-chaotic Chen oscillator was described in [65] and the fractional order equivalent is represented, as follows: where x, y, z, and w are the system dependent variables, a, b, c, and r are the system parameters with conventional values of 35, 3, 28, and 35 respectively, D q * is the Caputo's differential operator of order 0 < q ≤ 1, and φ is a variable parameter that determines the chaotic attractor. Figure 7 presents the LE spectrum of model (25), while the LEs of some selected states of the hyper-chaotic system are listed in Table 2. Generally, because hyper-chaotic system has two positive LEs, its prediction time is expected to be shorter than the original chaotic system [66]. The prediction time is given by: where K 1 is the Kolmogorov-Sinai entropy (KS-Entropy). The sum of the positive LEs is the upper bound of the KS-Entropy [67]. The KS-Entropy is calculated, as follows [67]: where L n are the LEs. The prediction time of the hyper-chaotic system is only shorter than the non-optimized system. Overall, DE-optimized system have the shortest prediction time.

Complexity of Optimization Codes
The complexity of the implemented algorithms in MATLAB were measured using the Halstead metrics, and the results are compared in Table 3. The procedure for determining the complexity begins with the counting of the distinct operators and operands that make up the software, where n 1 is the number of distinct operators, n 2 is the number of distinct operands, N 1 is the total number of operators, and N 2 is the total number of operands. Other parameters are derived from these four primary parameters. It should be noted that the metrics only measure the optimization algorithm, not including the numerical method to solve the fractional order system and numerical computation of the LEs. IWO implementation is the least complex of the three algorithms, while PSO is the most complex, as shown, most importantly, by the program length, which is N 1 + N 2 , and volume, which is a measure of the amount of codes written, based on the vocabulary (n 1 + n 2 ) and length N. Other Halstead parameters in the table that show IWO to be the simplest are the values of the difficulty involved in writing and maintaining the codes, mental effort required to implement the algorithm in MATLAB, and the estimated time that is required to fully implement the algorithm, which is less than 4 h, as compared to DE of over 5 h. Additionally, the IWO's estimated amount of delivered bugs is the least.

Discussion
The highest MLEs are obtained with optimized fractional order values that are above 0.7900 (Table 1). Although the three algorithms produce optimized MLEs within the same range of values (3.6000 ≤ MLE ≤ 3.6500), DE gives the overall best MLE. While the optimized parameters a, b, and q are slightly different in the three optimized results, the values of parameter c are the same. In the asymptotic analysis results shown in Table 1, it is seen that the equilibrium points EP 0 , EP 1 , and EP 2 of the DE, PSO, and IWO optimized systems are saddle points, just like the non-optimized system. For each of them, EP 0 is a saddle point of index 1, while EP 1 and EP 2 are the saddle points of index 2. According to the instability measure shown in Equation (24), at the saddle point of index 1 the systems are unstable for all the values of fractional order q, while the saddle points of index 2 become asymptotically stable below the following q values: 0.7964 for DE-optimized, 0.7981 for PSO-optimized, 0.7966 for IWO-optimized, and 0.8244 for non-optimized. Additionally, the measured complexity of the chaotic time series of the optimized systems and the non-optimized one shows consistency with the MLEs, with the DE-optimized system having the highest sample entropy and the non-optimized the lowest. This means that the DE-optimized system exhibits the most randomness in it. The improved randomness in the optimized systems makes them more suitable for developing a random number generator. Table 4 compares this investigation with some other works in [13,18,19,68,69], which have been carried out for the global optimization of dynamical behavior and parameter estimation in chaotic systems. It could be seen that our investigation is the only one that measures the complexity of the implemented optimization algorithms. It uses more algorithms than the others to compare their performance in the optimization problem. Additionally, unlike others, the work is based on a chaotic system of fractional order.  In our work, the fractional order is a part of the design variables, making a total of four. Furthermore, we use wider search spaces for the design variables. Hence, the higher MLEs in the present investigation can be mostly accredited to the higher population size, higher number of generations, and wider search spaces. The bifurcation and LE spectra are studied while separately varying each of the parameters, including the fractional order, unlike the compared works presented in Table 4. The maximization of the MLE in this investigation further shows the effectiveness of EAs as a global optimization algorithm [13,18,19,68,69]. Our work shows that the three algorithms obtain system parameter and fractional order values that produced better MLE than the non-optimized fractional order Chen system. The higher MLEs obtained is an indication that the optimized fractional order Chen oscillators have greater unpredictability than the non-optimized one. In addition, the smaller prediction times obtained for the optimized systems are credited to the increased MLE value. A system with shorter prediction time is regarded as safer for designing security systems [66,70]. Therefore, an improved secure communication system is guaranteed in the optimized fractional order Chen system.

Conclusions
This work demonstrates the application of DE, PSO, and IWO evolutionary algorithms in optimizing the values of parameters a, b, and c, and fractional order q of the chaotic Chen system. The meta-heuristics produce optimized MLEs that fall with a range of values (3.6000 ≤ MLE ≤ 3.6500), but with completely different parameter and fractional order values. This underlines the fact that chaotic behavior in a fractional order Chen system is multifaceted with respect to the parameter and fractional order values. The dynamics and complexity of the optimized systems were verified with bifurcation, LE spectrum, asymptotic, and sample entropy analyses of the results. The optimized MLEs in this investigation are higher than the non-optimized system, which is an indication of greater unpredictability in the optimized systems, with the DE giving the best optimal MLE. Therefore, the most reliable secure communication system is guaranteed with the DEoptimized system. The optimized systems are compared with a Chen hyper-chaotic system using the prediction time. The optimized systems were found to have a shorter prediction time than the hyper-chaotic system, which is also promising for improved secure communication system. Additionally, the complexity of the three implemented algorithms is determined and the Halstead metrics reveal that IWO has the simplest implementation.