Issues on Applying One- and Multi-Step Numerical Methods to Chaotic Oscillators for FPGA Implementation

Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Department of Electronics, INAOE, Tonantintla, Puebla 72840, Mexico; ing.omargufe@gmail.com (O.G.-F.); xerk.kun@gmail.com (M.F.M.-L.) * Correspondence: etlelo@inaoep.mx; Tel.: +52-222-2470-517 † These authors contributed equally to this work.


Introduction
Chaos is a nonlinear and unpredictable behavior that can be modeled by ordinary differential equations (ODEs). In continuous-time, the minimum number of ODEs for autonomous chaotic oscillators is three, as for example in [1,2]. A dynamical system modeled by four or more ODEs can generate hyper-chaotic behavior, as for example in [3]. Although sensitivity to initial conditions does not necessarily yield chaos [4], the majority of authors agree that the main characteristic of a dynamical system that generates chaos is the high sensibility to initial conditions, which is associated with a high unpredictability in the evolution of the time series of the state variables. The chaotic time series can be used to estimate Lyapunov exponents, as already shown in the seminal work [5], and by using the software for TIme SEries ANalysis (TISEAN) introduced in [6]. Lyapunov exponents are quite useful to characterize the behavior of a dynamical system, and they quantify the exponentially fast divergence or convergence of nearby orbits that can be seen in phase space.
Nowadays, it is said that a system with one positive Lyapunov exponent (LE+) is defined to be chaotic, and a system with more than one LE+ is hyper-chaotic. Some engineering applications of chaotic oscillators can be found in [7], which provides guidelines on the implementation by using field-programmable gate arrays (FPGAs), and shows the design of random number generators (RNGs) and chaotic secure communication systems [8]. The applications based on (hyper-)chaotic oscillators can be enhanced by guaranteeing higher unpredictability of the chaotic time series. One way is finding the chaotic oscillator having the highest LE+ [2], and also one must take into account other dynamical characteristics such as entropy and Kaplan-Yorke dimension (D KY ). For a chaotic oscillator having three ODEs, one computes three Lyapunov exponents, where one must be positive, one zero and one negative. There are methods to evaluate the Lyapunov spectrum [9,10], the seminal one was introduced in [5], and herein we apply TISEAN [6].
Recent works show the usefulness of chaotic oscillators in different engineering problems [11][12][13], however, there is no information on the issues related to the implementation of the numerical methods in electronic systems. In this manner, this paper uses three representative chaotic and two hyper-chaotic oscillators as case studies, which are listed in Table 1, along with their associated name, ODEs and parameter values that are used herein to generate chaotic behavior. The five chaotic oscillators are case studies to evaluate LE+ and D KY from their chaotic time series that are generated by applying three one-step and three multi-step methods. Representative numerical methods are chosen to be implemented on a FPGA and their hardware resources are counted to show the challenges of minimizing hardware resources while guaranteeing the highest exactness and stability of the numerical simulations.
The three chaotic and two hyper-chaotic oscillators that are case studies in this paper are detailed in Section 2. Three one-step and three multi-step numerical methods are given in Section 3. The chaotic time series of each state variable of each (hyper-)chaotic oscillator are generated by applying all the numerical methods, and the LE+ and D KY of each state variable are evaluated in Section 4. The FPGA implementation of representative numerical methods is detailed in Section 5. Finally, the conclusions are summarized in Section 6.

Chaotic and Hyper-Chaotic Oscillators
The three chaotic (modeled by three ODEs) and two hyper-chaotic (modeled by four ODEs) oscillators that are case studies herein are given in Table 1. In this Table CO1 is the well-known Lorenz system, introduced in 1963 as a simplified mathematical model for atmospheric convection [14], and from which was accidentally discovered the property associated to the high sensitivity to initial conditions. This originated one of the main characteristics in chaos theory and this CO1 is widely used as a work-horse to verify simulation and hardware implementation issues. In phase space, the Lorenz attractor resembles a butterfly effect, which stems from the real-world implications, i.e., in any physical system, the prediction of the evolution of the chaotic trajectories of the state variables will always fail in the absence of perfect knowledge of the initial conditions. In this manner, although physical systems can be completely deterministic, their chaotic behavior makes them inherently unpredictable (https://en.wikipedia.org/wiki/Lorenz_ system#cite_note-lorenz-1).
The chaotic oscillator labeled as CO2 is another well-known system introduced by Otto Rössler in 1976, originally intended to behave similarly to the Lorenz attractor, but its dynamical behavior is simpler and has only one manifold. In the Rössler system, an orbit within the attractor follows an outward spiral around an unstable fixed point. From the mathematical model of CO2 given in Table 1, this spiral effect is seen in the x, y plane, and once the graph spirals out enough, the z-dimension shows the influence of a second fixed point causing rise and twist. After the introduction of the Rössler system, important news was that the original model was useful in modeling equilibrium in chemical reactions [15].
CO3 is based on a saturated nonlinear function series that can be approximated by a piecewise-linear (PWL) function. Considering that the PWL function has saturation levels k i , break-points B i and slope m, then Equation (1) can be used to generate two scrolls, and Equation (2) to generate three scrolls. In a general sense, the PWL function given in Equation (1) can be increased to generate an even number of scrolls, and Equation (2) to generate an odd number of scrolls, as shown in [16].
The hyper-chaotic oscillators labeled as HO4 and HO5, both have more than three ODEs in order to have more than one positive Lyapunov exponent, so that they present a more complex behavior than chaotic oscillators modeled by three ODEs. Table 1. Chaotic and hyper-chaotic oscillators.

One-Step and Multi-Step Methods
The mathematical models of the chaotic and hyper-chaotic oscillators given in Table 1 can be formulated as initial value problems of the typeẋ = f (x). The solution of the ODEs can be performed by applying one-step and multi-step methods. The former requires values evaluated in one step x i to evaluate the next step denoted by x i+1 , while the multi-step methods require two or more previous step values denoted as to evaluate x i+1 . Other classifications are predictor or explicit and corrector or implicit methods. The explicit methods require past steps to evaluate the current step at iteration i + 1, but the implicit methods require estimation of the value at the current step i + 1 and past values at steps x i , x i−1 , x i−2,... . In this manner, it is common to name predictorcorrector [18] to the implicit methods, and they require an explicit method to evaluate the functions at the current iteration.
The explicit methods are faster than the implicit ones, but they may present numerical instability and lower exactness than the implicit methods. There are some rules for choosing the explicit method that is used within an implicit one to evaluate the current step x i+1 [18]. The step-size can also be varied during the computation or it can be constant and can be estimated from the stability analysis of the method, but one must take care of choosing the correct step-size to avoid non-convergence [19]. The explicit or predictor is the weak part in an implicit method due to the inherent truncation error, so that it puts a condition on the exactness of the initial prediction and the step-size of the corrector [18]. To enhance FPGA implementations of the numerical methods, the challenge is the selection of a method that allows a large step-size. That way, the larger the step-size of a numerical method, the higher the operating frequency of the FPGA implementation, as shown in Section 5.
The solution of the five chaotic oscillators given in the previous section are solved herein by applying the three one-step methods given in Table 2, and the three multi-step methods given in Table 3. The one-step methods are labeled as Forward Euler (FE), Backward Euler (BE) and fourth-order Runge-Kutta (RK4). The multi-step methods are labeled as sixth-order Adams-Bashforth (AB6), fourth-order Adams-Moulton (AM4) and fourthorder Gear (G4). Table 2. One-step methods.

Method
Iterative Equation Table 3. Multi-step methods.

Chaotic Time Series, LE+ and D KY
In Table 1, CO1 is the well-known Lorenz system, therefore, we show the simulation results for CO2, CO3, HO4 and HO5. The step-size h for each numerical method is given in the upper corner of each figure. One can appreciate that in some cases h is decreased to generate the same behavior provided by the majority of methods. Although the time evolution of the chaotic series is different for each method, the LE+ and D KY are similar, and it can be improved by varying h, which is not a trivial task and requires the analysis of the eigenvalues associated to each Jacobian matrix of each equilibrium point of each chaotic oscillator. Figures 1-4 show some chaotic time series of the (hyper-)chaotic oscillators simulated by applying the six numerical methods and listing the step-size h. The six methods were programmed into MatLab, and afterwards described in hardware language for FPGA implementation. In this case, a large h is desired to increase the operation frequency of an FPGA implementation, as shown in the following section.  Table 1 with initial conditions x 0 = y 0 = z 0 = 0.01.  Table 1 with initial conditions x 0 = y 0 = 0.1, z 0 = 0.0.   Table 1 with initial conditions x 0 = y 0 = z 0 = w 0 = 0.2.   Table 1 with initial conditions The D KY is evaluated from the Lyapunov exponents [9], and for an n-dimensional system it is evaluated by Equation (3), where LE 1 , . . . , LE n are Lyapunov exponent values ordered from the highest to the lowest value.
The LE+ and D KY were evaluated by TISEAN [6], which is based on the method introduced in [20]. The parameters for TISEAN are different for each state variable and analysis is performed using 50,000 samples for each chaotic time series. The LE+ for each state variable of each oscillator is shown in Table 4, and ordered from the highest to the lowest value. The highest LE+ is from the state variable x of HO4 and simulated with the fourth-order Runge-Kutta method, so that it is labeled as x_HO4_RK4. The same labels were adopted for the evaluation of D KY , whose results are shown in Table 5.

FPGA Implementation Issues
The development of engineering applications like chaotic secure communication systems and lightweight cryptography have positioned chaotic oscillators as a hot topic for research in this century. Nowadays, one can find implementations of chaotic systems using either analog or digital electronics, as already shown in [21]. This paper shows the implementation of (hyper-)chaotic oscillators from the selection of a numerical method, and by using FPGAs, which can be programmed/configured in the field after manufacture, and allow fast prototyping at relatively low development cost while providing good performance, computational power and programming flexibility.
Lets us consider the Lorenz oscillator (CO1) given in Table 1. The ODEs can be discretized by applying the most simple method known as Forward Euler (FE), to give the equations given in Equation (4). It is easy to see that these equations can be implemented by using multipliers, adders and subtractors. In addition, each block can be implemented including a clock (clk) and a reset (rst) pin to control the iterative process. As the multiplier consumes more power, if the multiplication includes a constant, as h, σ, ρ, β, one can design single-constant-multipliers (SCMs), as shown in [21], which use shift registers and adders to reduce power consumption and hardware resources. In this manner, the block description of Equation (4) is shown in Figure 5. The registers have an enable (en) pin and the description is divided into the macro-blocks labeled as Function Evaluation and Integrator FE. A counter is added to control the number of clks required in the FPGA implementation to evaluate the current iteration at n + 1, which is saved in the registers to process the next iteration.
x f e n+1 = x n + h[σ(y n − x n )] y f e n+1 = y n + h[−x n z n + ρx n − y n ] z f e n+1 = z n + h[x n y n − βz n ] The discretization of CO1 by applying an implicit method like Backward Euler (BE) is given in Equation (5), where it can be appreciated that the predictor is the FE given in Equation (4) to evaluate x f e n+1 , x f e n+1 , x f e n+1 . The block description for FPGA implementation is more complex and it embeds the FE method as shown in Figure 6. One can infer that the hardware resources for the BE method almost double compared to FE.
x n+1 = x n + h[σ(y f e n+1 − x f e n+1 )] y n+1 = y n + h[−x f e n+1 z f e n+1 + ρx f e n+1 − y f e n+1 ] z n+1 = z n + h[x f e n+1 y f e n+1 − βz f e n+1 ] The application of other one-step and multi-step methods to discretize a (hyper-)chaotic oscillator is performed in a similar manner as for the FE and BE methods. For example, the application of the multi-step sixth-order Adams-Bashforth (AB6) method is more complex than FE or BE. It requires five past steps associated to f (n), f (n − 1), f (n − 2), f (n − 3), f (n − 4), f (n − 5) that can be evaluated by the 4th-order Runge-Kutta (RK4) method. In this manner, using the iterative equation associated to AB6 that is given in Table 3, the discrete equations of CO1 are given in Equation (6). Figure 7 shows the block description for the FPGA implementation of CO1. One can see the predictor RK4, function evaluation and integrator Adams-Bashforth blocks, which are designed as for the FE and BE methods. The evaluation of Equation (6) also requires a finite-state machine (FSM) to control the iterative process, a cumulative sum block to process the RK4 method and random access memories (RAMs) to save the past steps f (n), f (n − 1), f (n − 2), f (n − 3), f (n − 4), f (n − 5) that are required for the next iteration, and they are controlled by STP (StarT Prediction) and EOP (End Of Prediction). The predictor RK4 is disconnected after the first iteration, which is controlled by the FSM.
In all the previous cases, the FPGA synthesis can be performed by adopting computer arithmetic of fixed-point notation, where the number of bits depends on the amplitudes of the state variables, as detailed in [7], where one can also find guidelines on Very High Speed Integrated Circuit Hardware Description Language (VHDL) programming. In this paper, the fixed-point notation has the format 12.20. Table 6 summarizes the hardware resources for the implementation of CO1, CO2, CO3, HO4 and HO5 applying FE and using FPGA Cyclone IV EP4CGX150DF31C7 under the synthesizer of software "Quartus II 13.0". In the same table, the last two rows provide the number of clk cycles that are required to evaluate one iteration n, and the latency is given in nanoseconds when using a 50 MHz clk signal. Register rst EnR times the hardware resources than FE, it is more exact and allows a higher h [7]. AB6 requires the higher number of hardware resources compared to the other five methods. If one does not design an SCM, the VHDL description of AB6 will require more than the 720 available multipliers in the FPGA Cyclone IV, and therefore it may not be implemented on this FPGA, so that one must use an FPGA with more density resources. The hardware resources for the FPGA implementation of the remaining chaotic systems labeled as CO2, CO3, HO4 and HO5, have similar increases for each numerical method, the main difference being due to the number of ODEs and nonlinear functions. Figure 8 shows a representative case of the FPGA implementation of CO1 applying the one-step method BE, and Figure 9 shows the application of the multi-step method G4, considering h = 0.001 in both cases.

Conclusions
We have shown the issues with the FPGA implementation of chaotic and hyperchaotic oscillators from the selection of a one-step and multi-step numerical method. The challenge is the selection of the time-step h to increase the frequency of operation of the FPGA design. It was appreciated that each one-step or multi-step method requires different hardware resources, so that trade-offs arise among reducing hardware resources, improving exactness and maximum operation frequency. Another open problem is the selection of the best chaotic oscillator, which can be done by evaluating the LE+ and D KY . This last characteristic increases as the number of ODEs increases, so that according to the results provided by TISEAN, the hyper-chaotic oscillators have the higher LE+ and D KY values. The FPGA implementation of the Lorenz system CO1 showed good agreement on the time series generated by applying BE and G4 methods, and using 32 bits in fixed-point notation of 12.30. The exactness can also be accomplished through using more bits, so that one can enhance applications in chaotic secure communications and the Internet of Things (IoT) to guarantee security and privacy. In particular, the IoT application requires a connectivity protocol in which chaotic oscillators can be synchronized to mask the data being transmitted, like in the extremely lightweight publish/subscribe messaging transport known as MQTT (mqtt.org), which is ideal for connecting remote devices with a small code footprint and minimal network bandwidth.