A New Chaotic System with a Self-Excited Attractor: Entropy Measurement, Signal Encryption, and Parameter Estimation

In this paper, we introduce a new chaotic system that is used for an engineering application of the signal encryption. It has some interesting features, and its successful implementation and manufacturing were performed via a real circuit as a random number generator. In addition, we provide a parameter estimation method to extract chaotic model parameters from the real data of the chaotic circuit. The parameter estimation method is based on the attractor distribution modeling in the state space, which is compatible with the chaotic system characteristics. Here, a Gaussian mixture model (GMM) is used as a main part of cost function computations in the parameter estimation method. To optimize the cost function, we also apply two recent efficient optimization methods: WOA (Whale Optimization Algorithm), and MVO (Multi-Verse Optimizer) algorithms. The results show the success of the parameter estimation procedure.

As we know, modeling of the real world chaotic systems has received great attention in recent decades [50][51][52][53][54][55]. Choosing proper values for model parameters is essential in chaotic systems, since they are very sensitive, both to model parameters and initial conditions. A slight change in the parameters of the chaotic system may cause important bifurcation in its behavior, because of the butterfly effect of the chaotic system [56]. Therefore, the parameter estimation problem of chaotic system models is a complex problem [57][58][59].
There are some widely used methods for the parameter estimation of the chaotic systems which are based on optimization methods [60][61][62]. In these methods, the problem of the parameter estimation is generally formulated as a cost function based on an error function between a time series obtained from a real system and a time series obtained from a known model with unknown parameters of that system. The goal of the parameter estimation method will then be to find the best values of the unknown parameters of the model which minimize the cost function. In addition, the optimization approaches have been used algorithms for this problem to find the best values of the unknown parameters as quickly as possible. They are algorithms such as genetic [63], particle swarm optimization [64], and evolutionary programming [65]. However, approaches that utilize cost function based on the error function seem to bear major limitations because of the butterfly effect of the chaotic systems [57][58][59].
It was remarked that the state space would be a proper domain to analyze the chaotic systems rather than the time-domain. The time series generated by the chaotic systems have random-like behavior in the time-domain, but they are ordered in the state space. They can show specific topologies in the state space named strange attractors. In this paper, we use a non-conventional metric as a useful cost function for the parameter estimation method. Accordingly, we model the attractor distribution of a real chaotic system by a parametric model named the Gaussian mixture model (GMM). It can provide flexible and probabilistic modeling for data distributions. GMM is also a commonly used parametric model in the pattern recognition and machine learning domain [66]. For example, in the speech recognition field, a set of GMMs was introduced to model phone attractors in a reconstructed phase space (RPS) in which the RPS is a time-independent domain similar to the state space [67][68][69]. The phone classification results showed that the GMM could be a useful model to capture the structure and topology of the speech attractors in the RPS. In addition, models of Gaussian mixture were recently used as the parameter identification method for some chaotic systems [70][71][72].
Here, to optimize the cost function, two recent efficient optimization methods are applied, including the WOA (Whale Optimization Algorithm), and MVO (Multi-Verse Optimizer) algorithm. Also, for testing the parameter estimation method in the chaotic systems, a real circuit is utilized based on a new chaotic system in this paper. All the data (time series) are obtained from the circuit that is designed based on the new chaotic system.
The contributions of this paper are described as: • A new 3D chaotic system with saddle equilibriums is proposed by a set of ordinary differential equations. • Dynamical properties of the 3D chaotic system are then reported that exhibit its dynamics.

•
The electronic circuit implementation of the 3D chaotic system is studied and used to present a random number generator (RNG), and its signal encryption is then introduced as an engineering application. • 1D and 2D parameter estimation of the electronic circuit is done by a GMM based cost function. • The cost function is optimized using two new efficient optimization methods called the WOA and the MVO algorithms. • By comparing the experimental data with numerically generated time series, the best-fitting parameters are found because the circuit had (almost) the same dynamics as the 3D chaotic system.
The structure of the paper is organized as follows: in the next section we introduce and analyze the new chaotic system with saddle equilibriums. In Section 2, we investigate it carefully through bifurcation analysis, spectrum of Lyapunov exponents, and its entropy. Section 3 deals with the circuit implementation of this new system and a real circuit application based on mobile RNG design. In the next section, the cost function based on the GMM is introduced. Two meta-heuristic optimization algorithms (WOA and MVO) are presented in Section 5. Results of the cost function and the parameter estimation of the new chaotic system using the WOA and MVO methods are available in Section 6. Finally, Section 7 is the conclusion of the paper.

A New Chaotic System and Its Analysis
In this section we introduce a new 3D system which can show chaotic behavior. Consider a system described with the following ordinary differential equations: This system is in the chaotic state when a = 4.0, b = 1.0, c = 1.0, d = 1.0, e = 1.0, f = 4.0 and g = 1.0. Different projections of the phase portrait for this system are plotted in Figure 1, which shows its strange attractor in 2D state spaces. System (1) is a new offset-boostable one [1][2][3][4] in which the variable z can be boosted with a direct constant in the first dimension. The cost function is optimized using two new efficient optimization methods called the WOA and the MVO algorithms.  By comparing the experimental data with numerically generated time series, the best-fitting parameters are found because the circuit had (almost) the same dynamics as the 3D chaotic system.
The structure of the paper is organized as follows: in the next section we introduce and analyze the new chaotic system with saddle equilibriums. In Section 3, we investigate it carefully through bifurcation analysis, spectrum of Lyapunov exponents, and its entropy. Section 4 deals with the circuit implementation of this new system and a real circuit application based on mobile RNG design. In the next section, the cost function based on the GMM is introduced. Two meta-heuristic optimization algorithms (WOA and MVO) are presented in Section 6. Results of the cost function and the parameter estimation of the new chaotic system using the WOA and MVO methods are available in Section 7. Finally, Section 8 is the conclusion of the paper.

A New Chaotic System and Its Analysis
In this section we introduce a new 3D system which can show chaotic behavior. Consider a system described with the following ordinary differential equations: This system is in the chaotic state when = 4.0, = 1.0, = 1.0, = 1.0, = 1.0, = 4.0, and = 1.0. Different projections of the phase portrait for this system are plotted in Figure 1, which shows its strange attractor in 2D state spaces. System (1) is a new offset-boostable one [1][2][3][4] in which the variable z can be boosted with a direct constant in the first dimension. This system has fixed-points in every ( * , * , * ) which satisfy the following equation, According to Equation (2), the system (1) has two equilibria in = (0.7321, 3.4641, 0) and = (0.7321, −3.4641, 0). The Jacobian matrix of the system (1) is and the corresponding eigenvalues for A and B are y y x x z z This system has fixed-points in every (x * , y * , z * ) which satisfy the following equation, According to Equation (2), the system (1) has two equilibria in A = (0.7321, 3.4641, 0) and B = (0.7321, −3.4641, 0). The Jacobian matrix of the system (1) is Therefore, both equilibria are saddle-foci. Thus, the attractor is self-excited.

Bifurcation Analysis
In this part, we investigate the behaviors of the system (1) with respect to changing parameter g. In part (A) of Figure 2 the bifurcation diagram of the system is shown and in part (B) of this figure Lyapunov exponents can be observed. It is important to be careful about numerical calculation of Lyapunov exponents, since improper use of usual methods may cause some issues [14,15,[73][74][75][76][77]. We have used the algorithm of [78] for computation of Lyapunov exponents.
Therefore, both equilibria are saddle-foci. Thus, the attractor is self-excited.

Bifurcation Analysis
In this part, we investigate the behaviors of the system (1) with respect to changing parameter . In part (A) of Figure 2 the bifurcation diagram of the system is shown and in part (B) of this figure Lyapunov exponents can be observed. It is important to be careful about numerical calculation of Lyapunov exponents, since improper use of usual methods may cause some issues [14,15,[73][74][75][76][77]. We have used the algorithm of [78] for computation of Lyapunov exponents. As can be seen in Figure 2A, changing parameter causes a familiar period doubling route to chaos. In addition, positive values of the Lyapunov exponents in Figure 2B show that the underlying system is the chaotic system. As can be seen in Figure 2A, changing parameter g causes a familiar period doubling route to chaos. In addition, positive values of the Lyapunov exponents in Figure 2B show that the underlying system is the chaotic system.

Entropy Analysis
There are many techniques to evaluate the system complexity from data. One of the most famous method which had been used since 1991 is Approximate Entropy (ApEn) [79]. ApEn can be applied to short and noisy data with outliers [80]. Therefore, many systems can be categorized by means of complexity [81]. Consider the N data sample u(1), u(2), . . . , u(N) with the vector sequence x(1), x(2), . . . , x(N − m + 1) ∈ R m which can be defined as: where m is an integer and determines the dimension of x(i) as the length of compared run of data. Then, for each i in the 1 ≤ i ≤ N − m + 1, the following equation is defined: where J is the number of correct vectors in Equation (7), the number of vectors that the distance (infinity norm or maximum norm) between them and x(i) is lower than r, and r is also a tolerance threshold that is defined by the product of a constant C to the standard deviation of data.
Then, the ApEn can be written as: The estimation of Equation (11) for N data sample is as follows, It can be derived that the ApEn values determine the similarity between chosen window and the sliding window of the data. Therefore, m determines the length of the window to be compared, and r is the tolerance threshold for accepting similar pattern between two windows. Figure 3 represents the ApEn diagram of the system (1) with respect to parameter g.

Entropy Analysis
There are many techniques to evaluate the system complexity from data. One of the most famous method which had been used since 1991 is Approximate Entropy (ApEn) [79]. ApEn can be applied to short and noisy data with outliers [80]. Therefore, many systems can be categorized by means of complexity [81]. Consider the N data sample (1), (2), … , ( ) with the vector sequence (1), (2), … , ( − + 1) ∈ which can be defined as: where m is an integer and determines the dimension of ( ) as the length of compared run of data. Then, for each in the 1 ≤ ≤ − + 1, the following equation is defined: where J is the number of correct vectors in Equation (7), the number of vectors that the distance (infinity norm or maximum norm) between them and ( ) is lower than r, and r is also a tolerance threshold that is defined by the product of a constant C to the standard deviation of data.
Then, the ApEn can be written as: The estimation of Equation (11) for N data sample is as follows, It can be derived that the ApEn values determine the similarity between chosen window and the sliding window of the data. Therefore, m determines the length of the window to be compared, and r is the tolerance threshold for accepting similar pattern between two windows. Figure 3 represents the ApEn diagram of the system (1) with respect to parameter .

Real Circuit Design of the New Chaotic System as a Mobile RNG and Its Application for Signal Encryption
Random number generator (RNG) algorithms produce a sequence of numbers with properties of randomness and they are a research subject since a few decades. Chaotic systems are commonly used in the random numbers generation algorithms because they are complex and very sensitive. In this section, a mobile RNG design is implemented based on the introduced chaotic system (1) and then signal encryption application is realized with the RNG.
The micro-computer based mobile RNG can be used in many fields especially in encryption studies with low cost and high performance. It is aimed at encryption of multimedia data (audio, image, video, text etc.) with the realized mobile RNG to be flexible and user friendly.
As far as we know, random number generators require high cost hardware like computers and FPGA in order to successfully pass the universal tests [82][83][84][85][86]. In this paper, the design of a microcomputer-based mobile RNG and a signal encryption application with the designed RNG is realized without needing hardware such as FPGA, computers, etc. Therefore, "Raspberry Pi 3" is used here as hardware which supports 64-bit processing capability. Since the "Raspberry Pi 3" card has 64-bit processing capability, it can generate very sensitive decimal numbers; thus, randomness of these generated numbers is very high. BCM2837 SoC (system-on-chip) 64-bit ARMv8 quad core Cortex A53 processor running @1.2GHz produced by Broadcom is available on the card. The general view of "Raspberry Pi 3" is as given as in Figure 4.

Real Circuit Design of the New Chaotic System as a Mobile RNG and Its Application for Signal Encryption
Random number generator (RNG) algorithms produce a sequence of numbers with properties of randomness and they are a research subject since a few decades. Chaotic systems are commonly used in the random numbers generation algorithms because they are complex and very sensitive. In this section, a mobile RNG design is implemented based on the introduced chaotic system (1) and then signal encryption application is realized with the RNG.
The micro-computer based mobile RNG can be used in many fields especially in encryption studies with low cost and high performance. It is aimed at encryption of multimedia data (audio, image, video, text etc.) with the realized mobile RNG to be flexible and user friendly.
As far as we know, random number generators require high cost hardware like computers and FPGA in order to successfully pass the universal tests [82][83][84][85][86]. In this paper, the design of a microcomputer-based mobile RNG and a signal encryption application with the designed RNG is realized without needing hardware such as FPGA, computers, etc. Therefore, "Raspberry Pi 3" is used here as hardware which supports 64-bit processing capability. Since the "Raspberry Pi 3" card has 64-bit processing capability, it can generate very sensitive decimal numbers; thus, randomness of these generated numbers is very high. BCM2837 SoC (system-on-chip) 64-bit ARMv8 quad core Cortex A53 processor running @1.2GHz produced by Broadcom is available on the card. The general view of "Raspberry Pi 3" is as given as in Figure 4. Our proposed circuit is used as an entropy source for RNG. Then, the NIST-800-22 tests are performed on random numbers to evaluate the performance of the designed RNG. In the next step, a signal encryption application is realized as an example application in "Raspberry Pi 3". Also, an electronic circuit implementation of the chaotic circuit is done in OrCAD-PSpice and on the oscilloscope.

Micro-Computer-Based Mobile RNG Design
As before mentioned, the "Raspberry Pi 3" board is used as a micro-computer for RNG design and encryption application. The chaotic system of (1) is also utilized in the RNG design. The RNG design steps are given in Algorithm 1 as a pseudo code. Our proposed circuit is used as an entropy source for RNG. Then, the NIST-800-22 tests are performed on random numbers to evaluate the performance of the designed RNG. In the next step, a signal encryption application is realized as an example application in "Raspberry Pi 3". Also, an electronic circuit implementation of the chaotic circuit is done in OrCAD-PSpice and on the oscilloscope.

Micro-Computer-Based Mobile RNG Design
As before mentioned, the "Raspberry Pi 3" board is used as a micro-computer for RNG design and encryption application. The chaotic system of (1) is also utilized in the RNG design. The RNG design steps are given in Algorithm 1 as a pseudo code. Entering parameters and initial condition of the chaotic system 3: Determination of the value of ∆h 4: Sampling with determination ∆h value 5: while (least 1 M. Bit data) do 6: Solving the chaotic system with RK4 7: Convert float to binary number (32 bit) 8: Select the bits (LSB-16 bit) from 32 bit binary number 9: end while 10: The implementation of NIST After entering parameters and initial condition of the system (1), the outputs are discretized with the RK4 differential equation solving method. Then, float numbers are obtained and converted into 32 bits binary numbers. Later, the RNG design is executed with obtained binary numbers. The last 16 bits of the outputs (x, y and z variables) are used in the design. The NIST-800-22 statistical tests are also used to prove the success of the RNG design [87]. The NIST-800-22 tests consist of 16 different tests such as monobit, serial and discrete Fourier transform tests. The p-values of the test should be greater than 0.001 in order to be counted as successful in NISTS-800-22 tests.
Our experiments show that the random numbers generated from x, y and z outputs successfully passed all the tests with the last 16 bits. The NIST-800-22 tests results are given in Table 1. The ready tested random numbers that pass all of the NIST-800-22 tests can be used in applications that require high security such as cryptology, data hiding, watermarking, etc. To obtain the random numbers, the pins x, y, and z GPIO (General purpose input/output) are utilized as shown in Figure 5. They are the 37th pin for x output, the 35th pin for y output, and the 38th pin for z output from "Raspberry Pi 3". The generated x, y and z outputs (first 50 bits) are shown in Figure 6 as real-time oscilloscope outputs. The 35th, 37th, and 38th GPIO pins give x, y and z outputs in Figure 5, respectively. They are used for real-time oscilloscope outputs.  The generated x, y and z outputs (first 50 bits) are shown in Figure 6 as real-time oscilloscope outputs. The 35th, 37th, and 38th GPIO pins give x, y and z outputs in Figure 5, respectively. They are used for real-time oscilloscope outputs. The generated x, y and z outputs (first 50 bits) are shown in Figure 6 as real-time oscilloscope outputs. The 35th, 37th, and 38th GPIO pins give x, y and z outputs in Figure 5, respectively. They are used for real-time oscilloscope outputs.

Signal Encryption Application Using "Raspberry Pi 3"
In this section, a signal encryption application with RNG that was generated from the proposed chaotic system is realized in "Raspberry Pi 3". The steps of the encryption and decryption process are given in Algorithm 2. In the encryption application, a signal that consists of 512 bits is used and shown (first 50 bits) in Figure 7 as the real-time oscilloscope outputs. Algorithm 2. Chaos based encryption and decryption algorithm pseudo code. 1: Getting ready to test random numbers for keys 3: Getting signal data to be encrypted 4: for i = 1 for all original data 5 random number bit xor original data bit 6: end 7: Encrypted data 8: for i = 1 for all encrypted data 9: random number bit xor encrypted data bit 10: end 11: Decrypted data 12: End For the encryption process, the 'XOR' operator is used. Figure 8 shows the first 50 bits of the encrypted signal as real-time oscilloscope outputs. Since the encryption process is performed for each bit, the size of the encrypted data is also 512. The same keys generated from the chaotic system are needed for decryption. With these keys, the original data can be obtained, again. The first 50 bits of the decrypted signal are shown in Figure  9 as the real-time oscilloscope outputs. As can be seen, comparing Figures 7 and 9, for the first 50 bits, there is no deformation. Getting ready to test random numbers for keys 3: Getting signal data to be encrypted 4: for i = 1 for all original data 5: random number bit xor original data bit 6: end 7: Encrypted data 8: for i = 1 for all encrypted data 9: random number bit xor encrypted data bit 10: end 11: Decrypted data 12: End For the encryption process, the 'XOR' operator is used. Figure 8 shows the first 50 bits of the encrypted signal as real-time oscilloscope outputs. Since the encryption process is performed for each bit, the size of the encrypted data is also 512.

Signal Encryption Application Using "Raspberry Pi 3"
In this section, a signal encryption application with RNG that was generated from the proposed chaotic system is realized in "Raspberry Pi 3". The steps of the encryption and decryption process are given in Algorithm 2. In the encryption application, a signal that consists of 512 bits is used and shown (first 50 bits) in Figure 7 as the real-time oscilloscope outputs. Algorithm 2. Chaos based encryption and decryption algorithm pseudo code. 1: Getting ready to test random numbers for keys 3: Getting signal data to be encrypted 4: for i = 1 for all original data 5 random number bit xor original data bit 6: end 7: Encrypted data 8: for i = 1 for all encrypted data 9: random number bit xor encrypted data bit 10: end 11: Decrypted data 12: End For the encryption process, the 'XOR' operator is used. Figure 8 shows the first 50 bits of the encrypted signal as real-time oscilloscope outputs. Since the encryption process is performed for each bit, the size of the encrypted data is also 512. The same keys generated from the chaotic system are needed for decryption. With these keys, the original data can be obtained, again. The first 50 bits of the decrypted signal are shown in Figure  9 as the real-time oscilloscope outputs. As can be seen, comparing Figures 7 and 9, for the first 50 bits, there is no deformation. The same keys generated from the chaotic system are needed for decryption. With these keys, the original data can be obtained, again. The first 50 bits of the decrypted signal are shown in Figure 9 as the real-time oscilloscope outputs. As can be seen, comparing Figures 7 and 9, for the first 50 bits, there is no deformation.
In the implemented method, a cryptoanalyser who wants to crack the encrypted data must know exactly all of the parameters and initial values of the chaotic system used in the encryption. Also, encrypted data will be not decrypted without "Raspberry Pi 3". In the implemented method, a cryptoanalyser who wants to crack the encrypted data must know exactly all of the parameters and initial values of the chaotic system used in the encryption. Also, encrypted data will be not decrypted without "Raspberry Pi 3".

Electronic Circuit Implementation of the Chaotic System in OrCAD-PSpice and on the Oscilloscope
In this part, we design an electronic circuit based on system (1) in OrCAD-PSpice ( Figure 10) and on the board ( Figure 11). The circuit includes simple electronic elements such as resistors, multipliers, capacitor, and opamps. Note that PSPICE simulation of chaotic circuits is quite trivial. In the literature, such systems are implemented with integrated circuit technology [88].   Figure 9. Decrypted Signal Data (first 50 bits).

Electronic Circuit Implementation of the Chaotic System in OrCAD-PSpice and on the Oscilloscope
In this part, we design an electronic circuit based on system (1) in OrCAD-PSpice ( Figure 10) and on the board (Figure 11). The circuit includes simple electronic elements such as resistors, multipliers, capacitor, and opamps. Note that PSPICE simulation of chaotic circuits is quite trivial. In the literature, such systems are implemented with integrated circuit technology [88]. In the implemented method, a cryptoanalyser who wants to crack the encrypted data must know exactly all of the parameters and initial values of the chaotic system used in the encryption. Also, encrypted data will be not decrypted without "Raspberry Pi 3".

Electronic Circuit Implementation of the Chaotic System in OrCAD-PSpice and on the Oscilloscope
In this part, we design an electronic circuit based on system (1) in OrCAD-PSpice ( Figure 10) and on the board (Figure 11). The circuit includes simple electronic elements such as resistors, multipliers, capacitor, and opamps. Note that PSPICE simulation of chaotic circuits is quite trivial. In the literature, such systems are implemented with integrated circuit technology [88].   The OrCAD-PSpice simulation outputs, which are two-dimensional phase portraits of the chaotic system, are seen in Figures 12 and 13, respectively. As can be seen from the ORCAD-PSpice outputs in Figure 12 and oscilloscope outputs in Figure 13, the results are similar.

Parameter Estimation of the Chaotic System
In this section, we introduce the parameter estimation method used for the chaotic circuit. This method utilizes a cost function which was adopted for the chaotic systems. The cost function of the parameter estimation method is based on a similarity metric using a parametric model of strange attractors in the state space. It was shown that this cost function could yield better results than the conventional error-based cost function over the time-domain [71]. The time-independent property of the state space is a sufficient reason to use this cost function because the state space can show complex behaviors of the strange attractor of chaotic systems [89].
As before mentioned, the utilized cost function is based on the attractor modeling; therefore, we need a model to represent the distribution of the attractor points in the state space. As a smooth  The OrCAD-PSpice simulation outputs, which are two-dimensional phase portraits of the chaotic system, are seen in Figures 12 and 13, respectively. As can be seen from the ORCAD-PSpice outputs in Figure 12 and oscilloscope outputs in Figure 13, the results are similar. The OrCAD-PSpice simulation outputs, which are two-dimensional phase portraits of the chaotic system, are seen in Figures 12 and 13, respectively. As can be seen from the ORCAD-PSpice outputs in Figure 12 and oscilloscope outputs in Figure 13, the results are similar.

Parameter Estimation of the Chaotic System
In this section, we introduce the parameter estimation method used for the chaotic circuit. This method utilizes a cost function which was adopted for the chaotic systems. The cost function of the parameter estimation method is based on a similarity metric using a parametric model of strange attractors in the state space. It was shown that this cost function could yield better results than the conventional error-based cost function over the time-domain [71]. The time-independent property of the state space is a sufficient reason to use this cost function because the state space can show complex behaviors of the strange attractor of chaotic systems [89].
As before mentioned, the utilized cost function is based on the attractor modeling; therefore, we need a model to represent the distribution of the attractor points in the state space. As a smooth  The OrCAD-PSpice simulation outputs, which are two-dimensional phase portraits of the chaotic system, are seen in Figures 12 and 13, respectively. As can be seen from the ORCAD-PSpice outputs in Figure 12 and oscilloscope outputs in Figure 13, the results are similar.

Parameter Estimation of the Chaotic System
In this section, we introduce the parameter estimation method used for the chaotic circuit. This method utilizes a cost function which was adopted for the chaotic systems. The cost function of the parameter estimation method is based on a similarity metric using a parametric model of strange attractors in the state space. It was shown that this cost function could yield better results than the conventional error-based cost function over the time-domain [71]. The time-independent property of the state space is a sufficient reason to use this cost function because the state space can show complex behaviors of the strange attractor of chaotic systems [89].
As before mentioned, the utilized cost function is based on the attractor modeling; therefore, we need a model to represent the distribution of the attractor points in the state space. As a smooth

Parameter Estimation of the Chaotic System
In this section, we introduce the parameter estimation method used for the chaotic circuit. This method utilizes a cost function which was adopted for the chaotic systems. The cost function of the parameter estimation method is based on a similarity metric using a parametric model of strange attractors in the state space. It was shown that this cost function could yield better results than the conventional error-based cost function over the time-domain [71]. The time-independent property of the state space is a sufficient reason to use this cost function because the state space can show complex behaviors of the strange attractor of chaotic systems [89].
As before mentioned, the utilized cost function is based on the attractor modeling; therefore, we need a model to represent the distribution of the attractor points in the state space. As a smooth parametric model, a Gaussian mixture model (GMM) can model the chaotic attractor geometry in the state space [67]. The GMM is a parametric probability density function represented by a weighted sum of Gaussian component densities [90]. It can model the distribution of the attractor points in the state space based on its powerful characteristics [91]. So far, metrics such as Kullback-Leibler divergence (also called relative entropy) were defined to measure distance between GMMs [90]. In addition, similarity-based metrics such as likelihood functions have been used to measure distance between a time series and a GMM. This idea was recently used as phone classification methods by parametric models of the distribution points of the speech signal in a high-dimension domain named RPS [67][68][69].
The GMMs have also been used for parameter estimation of some chaotic systems [70,71]. They were utilized similar to the task of the phone classification method. Suppose we have a chaotic system with a known model and its trajectory was recorded. We can then generate a GMM for the strange attractor of the chaotic system in the state space. Utilizing a distance-like metric over a likelihood function, we can compute dissimilarity between the learned GMM model of the real system attractor (with unknown parameters) and a distribution of a new attractor obtained by a system's model (with known parameters) in the state space to complete the parameter estimation method. Therefore, the score of the distance-based metric will be equal to the cost function of the parameter estimation method.

The GMM Computation as a Cost Function
A GMM with M mixtures is a weighted sum of M individual Gaussian densities. Each Gaussian density as a component of the GMM is represented by three main factors, mixture weight, mean vector, and covariance matrix. Therefore, they can be shown by a set of parameters, λ, as follows, and p(v|λ ) is the conditional probability of a D-dimensional single observation vector v given the GMM of λ. The p(v|λ ) can show a likelihood score. It expresses how probable the observed vector v is for the GMM of λ. In Equation (13) |.| is the determinant operator, exp(.) denotes the exponential function, and M is the number of mixtures (Gaussian components). In addition, for m-th mixture, w m ∈ [0, 1] is an scalar and named m-th mixing coefficient or mixture weight, µ m is the m-th D-dimensional mean vector, and Σ m is the m-th D × D covariance matrix. The mean vector and covariance matrix of a Gaussian component can show the center and the shape of points distribution around those of the component. It should be noted that the mixing coefficients w m are constrained to sum to 1, i.e., ∑ M m=1 w m = 1. As a problem which depends on the complexity of the data distribution, there is no analytical solution to determine the optimum number of GMM mixtures, M, needed for modeling of the attractors. Therefore, it is common to use a trial-and-error method to choose an adequate value of M. In our attractor modeling problem, to obtain a proper GMM model of the attractor in the state space, we evaluate some values of M. Generally, we need a higher value of M for attractor modeling if it has a very complex dynamic in the state space. One should note that while a higher number of mixtures can increase the performance of the cost function, it also increases the computational cost.
To find the similarity score between the attractor of a real system and the state space points of a specific model obtained from a chaotic system with known parameters (for example chaotic system (1)), the likelihood score can be calculated. Therefore, the parameter estimation of a known chaotic system with unknown parameters can be performed using the following two phases; a learning phase, here named "phase A", which includes fitting the GMM to the attractor of the real system, and an evaluation phase, named "phase B", to select the best values of parameters for the known chaotic model which causes the maximum similarity score or equally minimum distance score (cost function) over the learned GMM. Following are those phases in details:

Phase A
The first phase of the parameter estimation approach is the learning phase to find the GMM parameters, λ in Equation (13). The GMM learns the attractor's distribution of a real system, e.g., a chaotic circuit. Suppose S = {s 1 , s 2 , . . . , s N } is an N × D matrix consisting of N-samples of the time series of the real data in the D-dimensional state space. Therefore, each sample is a D-dimensional observation vector. To find the GMM parameters, an iterative expectation-maximization (EM) algorithm is utilized as follows:

Initialization Step
Initialize the mean vector µ m , covariance matrix Σ m and mixing coefficients w m in Equation (13) and evaluate the initial value of the logarithm of the likelihood score obtained from the input time series as follows,

Expectation
Step Evaluate values of r(s i , m), named responsibility of i-th sample of S given the m-th Gaussian component, using the current values of the GMM parameters:

Maximization Step
Re-estimate the parameters of the GMM utilizing the estimated values of the responsibilities as follows:

Likelihood Computation
Step Evaluate the logarithm of the likelihood score in Equation (14) and check for convergence criterion. If the convergence criterion is not satisfied, return to Section 5.2.2.

Phase B
The second phase is finding the best parameters of the known model of the chaotic system (with unknown parameters) using the learned GMM in the phase A. Here, the search space will be formed from a set of acceptable values of the model parameters. Now we suppose that the values of the parameters a&b of the system (1) are unknown. Then, for each pair of parameters (a, b), the chaotic system of (1) will be simulated, and a trajectory T(a, b) = (t 1 , t 2 , t 3 , . . . , t K |a, b ) with K samples will be obtained where each t K is D-dimensional measured data in the state space. Finally, using an average point-by-point log-likelihood score obtained from the learned GMM, λ, a similarity-based score is computed as follows, where T (a,b) is a matrix whose rows are composed from the state space vectors of the system trajectory with the model's parameters (a, b), and K is the number of the state space point. The parameter estimation method of the model is accomplished by computing Equation (20) and selecting the parameters of the model that can obtain the best similarity-based score, which here means the maximum score. If we use the negative of the similarity-based score, then the parameter estimation becomes a cost function minimization. Therefore, the best parameter selection, (a, b) * , would be conducted by the following criteria, J(.), based on the negative of mean log-likelihood score, Equation (21) shows the utilized cost function and (a, b) is the set of the system parameters (1). Here, λ is the learned GMM of the real system attractor obtained from the phase A. The objective of the parameter estimation method is to determine the parameters of the system, (a, b) when the cost function is minimized to result in the minimum value of J ((a, b)). The minimum value of the cost function guarantees the best solution with the proper parameters.

The GMM of Chaotic Circuit
Based on the observation vector v of the chaotic circuit, in this work, D = 3 is selected for the dimension of the state space according to system (1). Using the prepared real training data for the attractor by the chaotic circuit as a real system, the GMM will be specialized in order to model the geometry of that attractor. Figure 14 shows the attractor of the chaotic system in a three-dimensional state space with its GMM modeling using 256 Gaussian components, M = 256, where every three-dimensional ellipsoid corresponds to one of the Gaussian components. average point-by-point log-likelihood score obtained from the learned GMM, , a similarity-based score is computed as follows, where ( , ) is a matrix whose rows are composed from the state space vectors of the system trajectory with the model's parameters ( , ), and K is the number of the state space point. The parameter estimation method of the model is accomplished by computing Equation (20) and selecting the parameters of the model that can obtain the best similarity-based score, which here means the maximum score. If we use the negative of the similarity-based score, then the parameter estimation becomes a cost function minimization. Therefore, the best parameter selection,( , ) * , would be conducted by the following criteria, (. ), based on the negative of mean log-likelihood score, Equation (21) shows the utilized cost function and ( , )is the set of the system parameters (1). Here, is the learned GMM of the real system attractor obtained from the phase A. The objective of the parameter estimation method is to determine the parameters of the system, ( , ) when the cost function is minimized to result in the minimum value of (( , )). The minimum value of the cost function guarantees the best solution with the proper parameters.

The GMM of Chaotic Circuit
Based on the observation vector v of the chaotic circuit, in this work, D = 3 is selected for the dimension of the state space according to system (1). Using the prepared real training data for the attractor by the chaotic circuit as a real system, the GMM will be specialized in order to model the geometry of that attractor. Figure 14 shows the attractor of the chaotic system in a three-dimensional state space with its GMM modeling using 256 Gaussian components, = 256, where every threedimensional ellipsoid corresponds to one of the Gaussian components.

Optimization Algorithm
There are four major categories to which different kinds of optimization methods belong: Enumerative methods, Calculus-Based methods, Heuristic methods, and Meta-heuristic methods [62]. Meta-heuristic optimization algorithms, which cover a wide range of problems, are becoming more and more popular in engineering applications [92,93]. Nature can be named as one of the most important sources of inspiration for new meta-heuristic algorithms. On this subject, black and white holes in cosmology and Humpback whales in the sea aid in constructing the MVO (Multi-Verse Optimizer) and WOA (Whale Optimization Algorithm) meta-heuristic algorithm [94,95]. We now introduce these methods:

The Whale Optimization Algorithm
The WOA algorithm is based on the hunting behavior of Humpback whales which can encircle the recognized location of prey. The WOA algorithm assumes that the current best candidate solution is the target prey or is close to the optimum. The next step is about the attacking strategy which is the bubble-net strategy. Putting it all together, the proposed WOA method includes three major steps in the simulation: the search for prey, encircling prey, and the bubble-net foraging behavior of humpback whales. For complete details see [94].

Multi-Verse Optimizer: A Nature-Inspired Algorithm for GlobalOptimization
Another novel nature-inspired algorithm is Multi-Verse Optimizer (MVO). Cosmology (white hole, black hole, and wormhole) is the main inspiration of this algorithm. As mentioned before, every search process in the optimization algorithm consists of two phase: exploration and exploitation. The MVO supports this by white and black holes in order to respond to the exploration phase and wormholes for the exploitation phase. Further details are described in [95].

Experimental Results
In this section, some simulations are done to investigate the acceptability of the parameter estimation method of the chaotic circuit. We have used a fourth-order Runge-Kutta method with a step size of 10 ms and a total of 30, 000 samples corresponding to a time of 300 s. Here, we assume that the original chaotic system of (1) should be estimated by minimization of the GMM-based cost function.

Optimization Algorithm
There are four major categories to which different kinds of optimization methods belong: Enumerative methods, Calculus-Based methods, Heuristic methods, and Meta-heuristic methods [62]. Meta-heuristic optimization algorithms, which cover a wide range of problems, are becoming more and more popular in engineering applications [92,93]. Nature can be named as one of the most important sources of inspiration for new meta-heuristic algorithms. On this subject, black and white holes in cosmology and Humpback whales in the sea aid in constructing the MVO (Multi-Verse Optimizer) and WOA (Whale Optimization Algorithm) meta-heuristic algorithm [94,95]. We now introduce these methods:

The Whale Optimization Algorithm
The WOA algorithm is based on the hunting behavior of Humpback whales which can encircle the recognized location of prey. The WOA algorithm assumes that the current best candidate solution is the target prey or is close to the optimum. The next step is about the attacking strategy which is the bubble-net strategy. Putting it all together, the proposed WOA method includes three major steps in the simulation: the search for prey, encircling prey, and the bubble-net foraging behavior of humpback whales. For complete details see [94].

Multi-Verse Optimizer: A Nature-Inspired Algorithm for GlobalOptimization
Another novel nature-inspired algorithm is Multi-Verse Optimizer (MVO). Cosmology (white hole, black hole, and wormhole) is the main inspiration of this algorithm. As mentioned before, every search process in the optimization algorithm consists of two phase: exploration and exploitation. The MVO supports this by white and black holes in order to respond to the exploration phase and wormholes for the exploitation phase. Further details are described in [95].

Experimental Results
In this section, some simulations are done to investigate the acceptability of the parameter estimation method of the chaotic circuit. We have used a fourth-order Runge-Kutta method with a step size of 10 ms and a total of 30,000 samples corresponding to a time of 300 s. Here, we assume that the original chaotic system of (1) should be estimated by minimization of the GMM-based cost function.
First, using some 1D parameter estimation methods, different number of the GMM's components, = (64,96,128,192,256), are used to show the sufficiency of the cost function. The experimental results of the cost function versus the values of the parameters a&b are depicted in Figures 15 and 16, respectively.   As can be seen, all of the cost functions show convex functions around the desired point. Therefore, they are acceptable for the parameter estimation methods. Specifically, they show the effect of changing the parameter of the model as a monotonically trend along with a global minimum at the exact expected value of the desired parameters ( = 4.00, = 1.00). Therefore, the GMMbased cost function has the desired ideal properties for the parameter estimation problem. Moreover, Figures 15 and 16 show the effect of increasing the number of GMM components, M, used in the GMM modeling. In this case, = 256 represents better performance to identify the parameters a&b. In Figures 17 and 18, a contour plot of the cost function and its "cost surface" are respectively shown for the chaotic system (1) with = 256 along with variation in the parameters, a&b. They show dissimilarity between the real system attractor and each model attractor for a 2D parameter estimation problem. The minimum value of the point on those plots gives the parameters for the best model. As can be seen, all of the cost functions show convex functions around the desired point. Therefore, they are acceptable for the parameter estimation methods. Specifically, they show the effect of changing the parameter of the model as a monotonically trend along with a global minimum at the exact expected value of the desired parameters (a = 4.00, b = 1.00). Therefore, the GMM-based cost function has the desired ideal properties for the parameter estimation problem. Moreover, Figures 15 and 16 show the effect of increasing the number of GMM components, M, used in the GMM modeling. In this case, M = 256 represents better performance to identify the parameters a&b.
In Figures 17 and 18, a contour plot of the cost function and its "cost surface" are respectively shown for the chaotic system (1) with M = 256 along with variation in the parameters, a&b. They show dissimilarity between the real system attractor and each model attractor for a 2D parameter estimation problem. The minimum value of the point on those plots gives the parameters for the best model. As can be seen, all of the cost functions show convex functions around the desired point. Therefore, they are acceptable for the parameter estimation methods. Specifically, they show the effect of changing the parameter of the model as a monotonically trend along with a global minimum at the exact expected value of the desired parameters ( = 4.00, = 1.00). Therefore, the GMMbased cost function has the desired ideal properties for the parameter estimation problem. Moreover, Figures 15 and 16 show the effect of increasing the number of GMM components, M, used in the GMM modeling. In this case, = 256 represents better performance to identify the parameters a&b. In Figures 17 and 18, a contour plot of the cost function and its "cost surface" are respectively shown for the chaotic system (1) with = 256 along with variation in the parameters, a&b. They show dissimilarity between the real system attractor and each model attractor for a 2D parameter estimation problem. The minimum value of the point on those plots gives the parameters for the best model. As can be seen in Figures 17 and 18, the global minimum of the cost function is in the right place ( = 4.00 and = 1.00). Furthermore, the surface of the cost function is almost convex near the best parameters, which makes it an easy case for any optimization approach that moves downhill.
In order to examine the efficiency of the cost function in the parameter estimation, two mentioned meta-heuristic optimization methods are applied. All the basic parameters, such as maximum number of iterations (50) and number of search agent (25), are the same in both algorithms. For further details about the algorithms and their particular parameters, see [96]. Comparison between the performances of MVO and WOA optimization algorithm is shown in Figure 19. Based on the results of Figure 19, the MVO optimization method showed a superior performance in comparison with the WOA algorithms. In addition, Figure 20 represents the process of finding the best parameters using the WOA algorithm performed once for every 10 iterations. As can be seen, the individuals converge to the optimum area ( = 4.00 and = 1.00). As can be seen in Figures 17 and 18, the global minimum of the cost function is in the right place (a = 4.00 and b = 1.00). Furthermore, the surface of the cost function is almost convex near the best parameters, which makes it an easy case for any optimization approach that moves downhill.
In order to examine the efficiency of the cost function in the parameter estimation, two mentioned meta-heuristic optimization methods are applied. All the basic parameters, such as maximum number of iterations (50) and number of search agent (25), are the same in both algorithms. For further details about the algorithms and their particular parameters, see [96]. Comparison between the performances of MVO and WOA optimization algorithm is shown in Figure 19. As can be seen in Figures 17 and 18, the global minimum of the cost function is in the right place ( = 4.00 and = 1.00). Furthermore, the surface of the cost function is almost convex near the best parameters, which makes it an easy case for any optimization approach that moves downhill.
In order to examine the efficiency of the cost function in the parameter estimation, two mentioned meta-heuristic optimization methods are applied. All the basic parameters, such as maximum number of iterations (50) and number of search agent (25), are the same in both algorithms. For further details about the algorithms and their particular parameters, see [96]. Comparison between the performances of MVO and WOA optimization algorithm is shown in Figure 19. Based on the results of Figure 19, the MVO optimization method showed a superior performance in comparison with the WOA algorithms. In addition, Figure 20 represents the process of finding the best parameters using the WOA algorithm performed once for every 10 iterations. As can be seen, the individuals converge to the optimum area ( = 4.00 and = 1.00). Based on the results of Figure 19, the MVO optimization method showed a superior performance in comparison with the WOA algorithms. In addition, Figure 20 represents the process of finding the Entropy 2018, 20, 86 18 of 23 best parameters using the WOA algorithm performed once for every 10 iterations. As can be seen, the individuals converge to the optimum area (a = 4.00 and b = 1.00).

Conclusions
In this paper, a new chaotic system has been investigated carefully through bifurcation, the largest Lyapunov exponent, ApEn, and stability analysis. Then, an engineering application of that system was proposed using a random number generator and its signal encryption application. After that, a GMM-based cost function was utilized in the parameter estimation of the chaotic circuit designed from the chaotic system. The cost function was based on the minimization of dissimilarity between the phase portrait obtained from the real system and that obtained from the model of the chaotic system. In order to minimize the cost function and to obtain the correct parameters, we used two new efficient optimization methods, the Whale Optimization Algorithm (WOA), and Multi-Verse Optimizer (MVO) algorithm. The MVO optimization method showed superior performance in comparison with the WOA algorithm.

Parameter b
Parameter a a b c d Figure 20. Process of finding the best parameters using the WOA algorithm. (a-d) represent the first, 10th, 20th, and 30th iteration, respectively.

Conclusions
In this paper, a new chaotic system has been investigated carefully through bifurcation, the largest Lyapunov exponent, ApEn, and stability analysis. Then, an engineering application of that system was proposed using a random number generator and its signal encryption application. After that, a GMM-based cost function was utilized in the parameter estimation of the chaotic circuit designed from the chaotic system. The cost function was based on the minimization of dissimilarity between the phase portrait obtained from the real system and that obtained from the model of the chaotic system. In order to minimize the cost function and to obtain the correct parameters, we used two new efficient optimization methods, the Whale Optimization Algorithm (WOA), and Multi-Verse Optimizer (MVO) algorithm. The MVO optimization method showed superior performance in comparison with the WOA algorithm.