Strange Attractors Generated by Multiple-Valued Static Memory Cell with Polynomial Approximation of Resonant Tunneling Diodes

This paper brings analysis of the multiple-valued memory system (MVMS) composed by a pair of the resonant tunneling diodes (RTD). Ampere-voltage characteristic (AVC) of both diodes is approximated in operational voltage range as common in practice: by polynomial scalar function. Mathematical model of MVMS represents autonomous deterministic dynamical system with three degrees of freedom and smooth vector field. Based on the very recent results achieved for piecewise-linear MVMS numerical values of the parameters are calculated such that funnel and double spiral chaotic attractor is observed. Existence of such types of strange attractors is proved both numerically by using concept of the largest Lyapunov exponents (LLE) and experimentally by computer-aided simulation of designed lumped circuit using only commercially available active elements.


Introduction
A general property of chaos is long-time unpredictability; i.e., random-like evolution of dynamical system even if the describing mathematical model does not contain stochastic functions or parameters. Because of its nature, chaotic behavior was often misinterpreted as noise. The first mention of this kind of the complex solution was in [1] where Lorenz noticed the extreme sensitivity of autonomous deterministic dynamics to tiny changes of the initial conditions. After this very milestone, chaos started to be reported in many distinct scientific fields as well as daily life situations. Chaotic motion has been observed in chemical reactions [2], classical mechanics [3], hydrodynamics [4], brain activity [5], models of biological populations [6], economy [7] and, of course, in many lumped circuits.
Two basic vector field mechanisms are required for evolution of chaos: stretching and folding. The first mechanism is responsible for exponential divergence of two neighboring state trajectories and second one bounds strange attractor within a finite state space volume. Pioneering work showing the presence of robust chaotic oscillation within dynamics of simple electronic circuit is [8]. So far, the so-called Chua´s oscillator was subject of laboratory demonstrations, deep numerical investigations and many research studies [9][10][11]. Several interesting strange attractors associated with different vector field local geometries have been localized within the dynamics of three-segment piecewise-linear Chua systems [12][13][14]. However, the inventors of Chua´s oscillator built it intentionally to construct a vector field capable of generating chaotic waveforms. Progress in computational power together with development of the parallel processing allows chaos localization in standard functional blocks of radiofrequency subsystems such as in harmonic oscillators [15,16], frequency filters [17,18], phase-locked loops [19], power [20] and dc-dc [21] converters, etc. From a practical point of view, chaos represents an unwanted operational regime that needs to be avoided. It can be recognized among where the state vector is x = (v 1 , v 2 , i) T , C 1 and C 2 is parasitic capacitance of first and second RTD, L and R is summarized (RTDs are connected in series) lead inductance and resistance, respectively. Details about modeling of high frequency RTD including typical values of parasitic elements can be found in [25]. We can express both nonlinear functions (for k = 1, 2) as: Entropy 2018, 20, x 2 of 23 It can be recognized among regular behavior because of the specific features in the frequency domain: continuous and broadband frequency spectrum. However, an approach that is more sophisticated is to derive a set of describing differential equations and utilize the concept of LLE to find regions of chaotic solutions [22]. Searching for chaos in a mathematical model that describes simplified real electronic memory block is also topic of this paper. Three programs were utilized for the numerical analysis of MVMS: Matlab 2015 for the search-through-optimization algorithm including CUDA-based parallelization, Mathcad 15 for graphical visualization of results and Orcad Pspice 16 for circuit verification. The content of this paper is divided into four sections with the logical sequence: model description, numerical analysis, circuit realization and verification; both through simulation and measurement.

Dynamical Model of Fundamental MVMS
Basic mathematical model of MVMS [23,24] is given in Figure 1 and can be described by three first-order ordinary differential equations in the following form: where the state vector is x = (v1, v2, i) T , C1 and C2 is parasitic capacitance of first and second RTD, L and R is summarized (RTDs are connected in series) lead inductance and resistance, respectively. Details about modeling of high frequency RTD including typical values of parasitic elements can be found in [25]. We can express both nonlinear functions (for k = 1, 2) as: Thus, AVC of each RTD is a cubic polynomial that should form an N-type curve with a negative segment. Fixed points xe are all solutions of the problem dx/dt = 0. For further simplicity, let's assume that R = 0 Ω. We can determine global conditions and position for its existence within state space as each solution of system of the nonlinear algebraic equations, namely: 1 ( − 1 ) 3 + 1 ( − 1 ) + 1 = 2 ( − − 2 ) 3 + 2 ( − − 2 ) + 2 = − = 1 ( − 1 ) 3 + 1 ( − 1 ) + 1 Vector field geometry depends on the eigenvalues; i.e., roots of a characteristic polynomial. It can be calculated as det(s·E-J) = 0 where E is the unity matrix and J is the Jacobi matrix: Characteristic polynomial in symbolical form becomes: Thus, AVC of each RTD is a cubic polynomial that should form an N-type curve with a negative segment. Fixed points x e are all solutions of the problem dx/dt = 0. For further simplicity, let's assume that R = 0 Ω. We can determine global conditions and position for its existence within state space as each solution of system of the nonlinear algebraic equations, namely: Vector field geometry depends on the eigenvalues; i.e., roots of a characteristic polynomial. It can be calculated as det(s·E-J) = 0 where E is the unity matrix and J is the Jacobi matrix:  2 2 ) − 6 1 2 − 6 1 2 2 ] + 1 2 + 3 1 1 2 2 + 3 1 2 2 2 + 3 1 2 2 + 3 1 2 2 − 6 1 1 2 + 6 1 2 2 + 1 2 (9 4 − 18 1 3 − 18 2 3 + 9 1 2 2 2 + 9 1 2 2 + 9 2 2 2 − 36 1 2 2 − 18 1 2 2 + 18 1 2 2 ) + 2} + 1 + 2 + 3 2 2 + 3 1 1 2 + 3 2 2 2 + 3 2 ( 1 + 2 ) − 6 2 ( 2 + 2 ) − 6 ( 1 1 − 2 2 ) = 0, Obviously, symbolical expressions for the individual eigenvalues are very complicated and cannot further contribute to the better understanding of a vector field configuration and chaos evolution; check well-known Cardan rules. In situation, where xe and ye coordinate of equilibrium point is close to offset voltages represented by d1 and d2, characteristic polynomial simplifies into the relation: and the eigenvalues depend only on linear part of polynomial approximation of AVCs of RTDs. Until very recently, analysis of MVMS was focused only on a high-frequency modeling of RTDs, influence of a pulse driving force on overall stability [26] and global dynamics [27] and specification of the boundary planes [28]. However, existence of chaos has not been uncovered and examined.

Numerical Results and Discussion
Numerical values of MVMS parameters can be obtained by the optimization technique described in [29]. In this case, the mathematical model of MVMS was considered piecewise-linear. Such kind of a vector field allows better understanding of chaos evolution, allows partial analytic solution and makes linear analysis generally more powerful. However, situation when AVCs of both RTDs are approximated by the polynomial functions is closer to reality. Thus, our problem stands as follows: couple of three-segment piecewise-linear functions needs to be transformed into the cubic (or higherorder if necessary) polynomial functions without losing robust chaotic solution having numerically close metric fractal dimension (Kaplan-Yorke is preferred over capacity because of rapid and precise calculation). For more details, readers should consult [30] where the inverse problem has been successfully solved. Finally, chaotic attractor like the so-called Rossler attractor [31] was localized.
Searching within smooth vector field (1) and by considering default normalized values C1 = 11 F, C2 = 37 F, L = 100 mH and R = 0 Ω leads to the following optimal values of the cubic polynomials (2): 1 = 648.5, 1 = −23.1, 1 = −0.13, 1 = 0.3, 2 = 58.1, 2 = −21.9, 2 = −0.65, 2 = 0.5. (7) Adopting these values and fourth-order Runge-Kutta numerical integration process we get the reference chaotic orbits provided in Figure 2 Obviously, symbolical expressions for the individual eigenvalues are very complicated and cannot further contribute to the better understanding of a vector field configuration and chaos evolution; check well-known Cardan rules. In situation, where x e and y e coordinate of equilibrium point is close to offset voltages represented by d 1 and d 2 , characteristic polynomial simplifies into the relation: and the eigenvalues depend only on linear part of polynomial approximation of AVCs of RTDs. Until very recently, analysis of MVMS was focused only on a high-frequency modeling of RTDs, influence of a pulse driving force on overall stability [26] and global dynamics [27] and specification of the boundary planes [28]. However, existence of chaos has not been uncovered and examined.

Numerical Results and Discussion
Numerical values of MVMS parameters can be obtained by the optimization technique described in [29]. In this case, the mathematical model of MVMS was considered piecewise-linear. Such kind of a vector field allows better understanding of chaos evolution, allows partial analytic solution and makes linear analysis generally more powerful. However, situation when AVCs of both RTDs are approximated by the polynomial functions is closer to reality. Thus, our problem stands as follows: couple of three-segment piecewise-linear functions needs to be transformed into the cubic (or higher-order if necessary) polynomial functions without losing robust chaotic solution having numerically close metric fractal dimension (Kaplan-Yorke is preferred over capacity because of rapid and precise calculation). For more details, readers should consult [30] where the inverse problem has been successfully solved. Finally, chaotic attractor like the so-called Rossler attractor [31] was localized.
Searching within smooth vector field (1) and by considering default normalized values C 1 = 11 F, C 2 = 37 F, L = 100 mH and R = 0 Ω leads to the following optimal values of the cubic polynomials (2): Adopting these values and fourth-order Runge-Kutta numerical integration process we get the reference chaotic orbits provided in Figure 2 146. This is an interesting geometric configuration of the vector field: a funnel-type strange attractor generated by saddle-node equilibria with the stability index 1, 0 and 0; saddle-focuses are completely missing. Developed optimization algorithm can be utilized to maximize unpredictability of a dynamical flow and increase system entropy. Starting with values (7) the following set of numerical values was obtained:  Figure 4. Note that local geometries of a vector field are not affected since, in a closed neighborhood of the fixed points, dynamical movement is given by three eigenvectors as in the case of a funnel chaotic attractor. Also note that one direction of the flow (along the eigenvector associated with λ 1 ) is strongly attracting; this nature is evident from Figure 5.
Interesting fragments of the bifurcation diagrams are depicted in Figure 6. are not affected since, in a closed neighborhood of the fixed points, dynamical movement is given by three eigenvectors as in the case of a funnel chaotic attractor. Also note that one direction of the flow (along the eigenvector associated with λ1) is strongly attracting; this nature is evident from Figure 5.
Interesting fragments of the bifurcation diagrams are depicted in Figure 6.   are not affected since, in a closed neighborhood of the fixed points, dynamical movement is given by three eigenvectors as in the case of a funnel chaotic attractor. Also note that one direction of the flow (along the eigenvector associated with λ1) is strongly attracting; this nature is evident from Figure 5.
Interesting fragments of the bifurcation diagrams are depicted in Figure 6.        Here, transient motion has been removed and plane x = d1 has been used for individual slices. These are plotted against parameters ak, bk and ck of polynomial approximations of AVC of k-th RTD. Parameters a1 and a2 can burn or bury chaotic attractor since these affect eigenvalues in "outer" parts of the active vector field (given by size of the state attractor) while geometry around "middle" fixed point remains almost unchanged. Figure 7 shows numerical calculation of a gained energy with small time step with respect to the state space location; red regions mark large increment while dark blue stands for a small evolution.
As stated before, chaotic solution is sensitive to the changes of parameters ak, bk, ck and dk, where k = 1, 2. If we consider these values variable we create eight-dimensional hyperspace in which chaos    Here, transient motion has been removed and plane x = d1 has been used for individual slices. These are plotted against parameters ak, bk and ck of polynomial approximations of AVC of k-th RTD. Parameters a1 and a2 can burn or bury chaotic attractor since these affect eigenvalues in "outer" parts of the active vector field (given by size of the state attractor) while geometry around "middle" fixed point remains almost unchanged. Figure 7 shows numerical calculation of a gained energy with small time step with respect to the state space location; red regions mark large increment while dark blue stands for a small evolution.
As stated before, chaotic solution is sensitive to the changes of parameters ak, bk, ck and dk, where k = 1, 2. If we consider these values variable we create eight-dimensional hyperspace in which chaos Here, transient motion has been removed and plane x = d 1 has been used for individual slices. These are plotted against parameters a k , b k and c k of polynomial approximations of AVC of k-th RTD. Parameters a 1 and a 2 can burn or bury chaotic attractor since these affect eigenvalues in "outer" parts of the active vector field (given by size of the state attractor) while geometry around "middle" fixed point remains almost unchanged. Figure 7 shows numerical calculation of a gained energy with small time step with respect to the state space location; red regions mark large increment while dark blue stands for a small evolution.
Maximal merit of LLE is 0.089 for a value set (7) and 0.103 for (8) and associated Kaplan-Yorke dimension [33] is 2.016 and 2.021 respectively. Specification of sufficiently large "chaotic" area is important also from practical viewpoint; as will be clarified later.
As stated before, chaotic solution is sensitive to the changes of parameters a k , b k , c k and d k , where k = 1, 2. If we consider these values variable we create eight-dimensional hyperspace in which chaos alternates with periodic orbits or fixed-point solution. Two-dimensional subspaces hewed-out from such hyperspace are provided in Figure 8. Each graph is composed of 101 × 101 = 10,201 points, calculation routine deals with time interval t ∈ (100, 10 4 ), random initial conditions inside basin of attraction and Gram-Smith orthogonalization [32]. Here, dark blue represents a trivial solution, light blue a limit cycle, green color stands for weak chaos and yellow marks strong chaotic behavior. Note that LLE for set of values (7) can be found in each visualized plot.
Maximal merit of LLE is 0.089 for a value set (7) and 0.103 for (8) and associated Kaplan-Yorke dimension [33] is 2.016 and 2.021 respectively. Specification of sufficiently large "chaotic" area is important also from practical viewpoint; as will be clarified later.

Circuitry Realization of MVMS-Based Chaotic Oscillators
Design of analog equivalent circuit is common way how to prove existence of structurally stable strange attractors within the dynamics of a prescribed set of ordinary differential equations. Realization of such the so-called chaotic oscillators is a simple and straightforward task that we can solve by using several approaches [34][35][36][37][38], both using discrete components and in integrated form. A favorite method that allows us to utilize commercially available active elements follows the concept of analog computers. Thus, only three building blocks are required for circuit construction: inverting summing integrator, summing amplifier and, in the case of a polynomial nonlinearity, four-segment analog multiplier. Fully analog circuit implementation is shown in Figure 9.  a1-b1, a1-c1, a1-a2, a1-b2, a1-c2, b1-c1, b1-a2, b1-b2, b1-c2, c1-a2, c1-b2, c1-c2, a2-b2, a2-c2, b2-c2.

Circuitry Realization of MVMS-Based Chaotic Oscillators
Design of analog equivalent circuit is common way how to prove existence of structurally stable strange attractors within the dynamics of a prescribed set of ordinary differential equations. Realization of such the so-called chaotic oscillators is a simple and straightforward task that we can solve by using several approaches [34][35][36][37][38], both using discrete components and in integrated form. A favorite method that allows us to utilize commercially available active elements follows the concept of analog computers. Thus, only three building blocks are required for circuit construction: inverting summing integrator, summing amplifier and, in the case of a polynomial nonlinearity, four-segment analog multiplier. Fully analog circuit implementation is shown in Figure 9. Note that this circuit synthesis requires many active devices: two TL084, two AD844 and a single four channel four quadrant analog multiplier MLT04. Supply voltage is symmetrical ±15 V but, for MLT04, this voltage is lowered to ±5 V. Majority of the analog realizations of the chaotic systems with a polynomial nonlinearity utilize AD633, i.e., single channel multiplier. It is possible also in our case but with the cost of eight active devices. Dynamical range for correct operation of MLT04 is only ±2 V. However, prescribed strange attractor is smaller in v1-v2 dimension. Advantage of this circuit is that individual MVMS parameters can be adjusted independently using potentiometers. Note that this circuit synthesis requires many active devices: two TL084, two AD844 and a single four channel four quadrant analog multiplier MLT04. Supply voltage is symmetrical ±15 V but, for MLT04, this voltage is lowered to ±5 V. Majority of the analog realizations of the chaotic systems with a polynomial nonlinearity utilize AD633, i.e., single channel multiplier. It is possible also in our case but with the cost of eight active devices. Dynamical range for correct operation of MLT04 is only ±2 V. However, prescribed strange attractor is smaller in v 1 -v 2 dimension. Advantage of this circuit is that individual MVMS parameters can be adjusted independently using potentiometers. Theoretically, using different decomposition of the polynomial functions, total number of the active Entropy 2018, 20, 697 9 of 23 elements can be lowered to four. Proposed chaotic oscillator is uniquely described by the following set of the differential equations: where V C1 , V C2 and V bias are independent dc voltage sources and K = 0.4 is internally-trimmed scaling constant of the analog multiplier cells MLT04. Fundamental time constant of this chaotic oscillator is chosen to be τ = R·Cv = 10 3 ·10 −7 = 100 µs.
Of course, it is also possible to build both analog networks provided in Figure 1 directly. Instead of nonlinear two-ports we must construct a couple of resistors with polynomial AVC; systematic design towards these network elements can be found in [39,40]. Circuitry realization of original MVMS with state vector x = (v 1 , v 2 , i) T is demonstrated by means of Figure 10; i.e., the state variables are voltages across grounded capacitors and current flowing through the inductor. Note that both polynomial function (2) need to be rewritten into the form i = f (v) = a·v 3 + b·v 2 + c·v + d. Thus, a new set of the dimensionless ordinary differential equations to be implemented as lumped analog electronic circuit as follows: Entropy 2018, 20, x 9 of 23 Theoretically, using different decomposition of the polynomial functions, total number of the active elements can be lowered to four. Proposed chaotic oscillator is uniquely described by the following set of the differential equations: ) + ] where VC1, VC2 and Vbias are independent dc voltage sources and K = 0.4 is internally-trimmed scaling constant of the analog multiplier cells MLT04. Fundamental time constant of this chaotic oscillator is chosen to be τ = R·Cv = 10 3 ·10 −7 = 100 μs.
Of course, it is also possible to build both analog networks provided in Figure 1 directly. Instead of nonlinear two-ports we must construct a couple of resistors with polynomial AVC; systematic design towards these network elements can be found in [39,40]. Circuitry realization of original MVMS with state vector x = (v1, v2, i) T is demonstrated by means of Figure 10; i.e., the state variables are voltages across grounded capacitors and current flowing through the inductor. Note that both polynomial function (2) need to be rewritten into the form i = f(v) = a·v 3 + b·v 2 + c·v + d. Thus, a new set of the dimensionless ordinary differential equations to be implemented as lumped analog electronic circuit as follows:  Now assume that impedance and frequency norm are 10 5 and 10 4 , respectively. Such values lead to the nominal inductance 1H. This simplified concept of chaotic oscillators is given in Figure 10 and described by following set of the ordinary differential equations: where a small value R S can still model lead to resistances of both RTDs that are parts of MVMS in Figure 1.
Note that this kind of realization utilizes second generation current conveyors implemented by using ideal voltage-controlled voltage-source E and current-controlled current-source F. A positive variant of this active three-port element is commercially available as the AD844 while a negative variant is EL2082 (only one negative device is required). In practice, the inductor should be substituted by the synthetic equivalent; i.e., active floating gyrator (Antoniou's sub-circuit) with a capacitive load, check Figure 11. In this case, the number of the active elements raises to eight: a single TL084, six AD844s and a single MLT04. Dynamical behavior is uniquely determined by the following mathematical model: where R in represents input resistance of current input terminal of AD844. Its typical value is 50 Ω and, due to the high values of the polynomial coefficients and in the case of small impedance norm chosen, it cannot be generally neglected. On the other hand, we must avoid output current saturation of each AD844 and consider frequency limitations of each active device. Thus, choice of both normalization factors is always a compromise. Thanks to the symmetry inside floating Antoniou's structure R g5 = R g3 , R g7 = R g1 , R g6 = R g2 and C g2 = C g1 . Behavior of this chaotic oscillator is extremely sensitive to the working resistors connected in the nonlinear two-terminal devices. Thus, calculated values were specified by Orcad Pspice optimizer where fitness functions (several should be defined to create tolerance channel) are absolute difference between polynomials in (11) and actual input resistance of designed circuit. Corresponding dc sweep analysis were estimated for input voltages from 0 V to 2 V with step 10 mV. Last circuitry implementation is provided in Figure 12; namely dual network to the original MVMS. Impedance norm is chosen to be 10 3 and frequency normalization factor is 10 5 . To get the reasonable values of the resistors further impedance rescaling is possible. Set of describing differential equations can be expressed as follows: Last circuitry implementation is provided in Figure 12; namely dual network to the original MVMS. Impedance norm is chosen to be 10 3 and frequency normalization factor is 10 5 . To get the Entropy 2018, 20, 697 12 of 23 reasonable values of the resistors further impedance rescaling is possible. Set of describing differential equations can be expressed as follows: where r Tk is trans-resistance of k-h ideal current-controlled voltage-source and I 1 is independent dc current source.
Entropy 2018, 20, x 12 of 23 where rTk is trans-resistance of k-h ideal current-controlled voltage-source and I1 is independent dc current source. Polynomial nonlinearity is implemented by ideal multiplication using block MULT. Current flowing through both inductors can be sensed via small resistor RS1 and RS2 respectively. However, these resistors represent error terms that are inserted into describing differential equations; similarly as a parasitic shunt resistor Rp. Such parasitic properties change global dynamics, can boost system dimensionality (not in this situation) and desired chaotic behavior can eventually disappear.
Unfortunately, active elements where output voltage is controlled by input current are not offthe-shelf components. However, trans-resistance amplifier can be constructed using single standard operational amplifier and feedback resistor. To do this, input is fed directly to − terminal, node + is connected to ground, resistor between − and OUT. Node OUT also represents output of a designed transimpedance amplifier. Equivalently, the AD844 can also do the trick. Input can be connected to − terminal, ground to + terminal, resistor between C and ground and output to OUT. Considering the latter case, number of the active devices becomes seven: one MLT04 and six AD844.
Colored plots provided in Figure 8 demonstrate that region of chaos around discovered values (7) is wide enough to provide a structurally stable strange attractor; both funnel and double-scroll type. The same analysis was performed for set (8); this strange attractor should be also observable.

Circuit Simulation, Experimental Verification and Comparison
As mentioned before, AVC of both nonlinear resistors in Figure 11 should be as precise as possible to reach desired state space attractor. To fulfil this requirement, build-in Orcad Pspice Polynomial nonlinearity is implemented by ideal multiplication using block MULT. Current flowing through both inductors can be sensed via small resistor R S1 and R S2 respectively. However, these resistors represent error terms that are inserted into describing differential equations; similarly as a parasitic shunt resistor R p . Such parasitic properties change global dynamics, can boost system dimensionality (not in this situation) and desired chaotic behavior can eventually disappear.
Unfortunately, active elements where output voltage is controlled by input current are not off-the-shelf components. However, trans-resistance amplifier can be constructed using single standard operational amplifier and feedback resistor. To do this, input is fed directly to − terminal, node + is connected to ground, resistor between − and OUT. Node OUT also represents output of a designed transimpedance amplifier. Equivalently, the AD844 can also do the trick. Input can be connected to − terminal, ground to + terminal, resistor between C and ground and output to OUT. Considering the latter case, number of the active devices becomes seven: one MLT04 and six AD844. Figure 8 demonstrate that region of chaos around discovered values (7) is wide enough to provide a structurally stable strange attractor; both funnel and double-scroll type. The same analysis was performed for set (8); this strange attractor should be also observable.

Circuit Simulation, Experimental Verification and Comparison
As mentioned before, AVC of both nonlinear resistors in Figure 11 should be as precise as possible to reach desired state space attractor. To fulfil this requirement, build-in Orcad Pspice optimizer can be adopted; see Figure 13 where optimal AVC of one nonlinear two-terminal device is reached. Of course, the operational regime of this nonlinear resistor needs to be limited; in our case to input voltages into range starting with −0.5 V and ending with 2 V. Note that up to thirteen fitness functions were supposed to cover predefined operational range. Also note that optimization was stopped while some objective functions still have nonzero error (expressed in terms of the differential percentages). Important is the matter of similarity between simulated and numerically integrated strange attractor.  (middle table).
Implementation of the chaotic oscillator based on the integrator block schematic is more robust; i.e., the desired strange attractor is less vulnerable to the passive component matching: fabrication series (E12) and component tolerances (1%) were considered for design. We are also experiencing superior parasitic properties of the voltage-mode active elements: both integrated circuits TL084 and MLT04 have very high input and very low output resistances. Values of resistors that form external network connected to AD844 were chosen such that parasitic input and output impedances can be neglected. This is the main reason why this kind of MVMS realization was picked and forwarded into practical experiments and undergoes laboratory measurement.
Generally, parasitic properties of the active elements play important role in a design process of analog chaotic oscillators and, of course, should be minimized. Besides input and output impedances also roll-off effects of transfer constants need to be analyzed. Several publications have been devoted to reveal and study problems associated with the parasitic properties of the specific active devices and how these affect global dynamics; for example [41].
In simulation profile of transient response, final time was set to 1 s and maximal time step 10 μs. This setup is kept for each simulation mentioned in this section. Circuit simulation associated with Figure 9 is provided in Figure 14. To transform the funnel into a double-scroll chaotic attractor we should change value of the resistors Rz3, Rz1 and R10. Simulation results associated with direct realization of original MVMS provided in Figure 10 is demonstrated by means of Figure 15.
Computer-aided analysis of chaotic oscillator with idealized controlled sources is showed in Figure  16.
Existence of observable strange attractors and numerically expected routing-to-chaos scenarios within dynamics of the fundamental MVMS has been proved experimentally; by its construction on the bread-board (see Figure 17) and consequent laboratory measurement using the analog Of course, the operational regime of this nonlinear resistor needs to be limited; in our case to input voltages into range starting with −0.5 V and ending with 2 V. Note that up to thirteen fitness functions were supposed to cover predefined operational range. Also note that optimization was stopped while some objective functions still have nonzero error (expressed in terms of the differential percentages). Important is the matter of similarity between simulated and numerically integrated strange attractor.
Implementation of the chaotic oscillator based on the integrator block schematic is more robust; i.e., the desired strange attractor is less vulnerable to the passive component matching: fabrication series (E12) and component tolerances (1%) were considered for design. We are also experiencing superior parasitic properties of the voltage-mode active elements: both integrated circuits TL084 and MLT04 have very high input and very low output resistances. Values of resistors that form external network connected to AD844 were chosen such that parasitic input and output impedances can be neglected. This is the main reason why this kind of MVMS realization was picked and forwarded into practical experiments and undergoes laboratory measurement.
Generally, parasitic properties of the active elements play important role in a design process of analog chaotic oscillators and, of course, should be minimized. Besides input and output impedances also roll-off effects of transfer constants need to be analyzed. Several publications have been devoted to reveal and study problems associated with the parasitic properties of the specific active devices and how these affect global dynamics; for example [41].
In simulation profile of transient response, final time was set to 1 s and maximal time step 10 µs. This setup is kept for each simulation mentioned in this section. Circuit simulation associated with Figure 9 is provided in Figure 14. To transform the funnel into a double-scroll chaotic attractor we should change value of the resistors R z3 , R z1 and R 10 . Simulation results associated with direct realization of original MVMS provided in Figure 10 is demonstrated by means of Figure 15. Computer-aided analysis of chaotic oscillator with idealized controlled sources is showed in Figure 16. oscilloscope. Due to its simpler realization and increased robustness, only circuitry illustrated by means of Figure 9 was decided for a real measurement. Selected chaotic waveforms in a time domain are provided by means of Figure 18. Independent voltage source Vbias = 750 mV in series with Rbias = 1 kΩ can be replaced by the positive supply voltage +5 V and the fixed series resistance Rbias = 6.7 kΩ. Analogically, combinations Vc1 = −130 mV, Rx = 1 kΩ and Vc2 = −650 mV, Ry = 1 kΩ can be replaced by the negative supply voltage −5 V, Rx = 38.5 kΩ and the same voltage −5 V, Ry = 7.7 kΩ respectively.

Figure 14.
Orcad Pspice circuit simulation associated with chaotic circuit given in Figure 9: selected plane projection v3-v2 (blue, upper graph), v2-v1 (red) and v3-v1 (blue) of the chaotic attractor, generated chaotic signal v1 (red) and v2 (blue) in the time domain, chaotic waveform v1 (red) and v2 (blue) in the frequency domain. Note that significant frequency components are concentrated in audio range.

Figure 14.
Orcad Pspice circuit simulation associated with chaotic circuit given in Figure 9: selected plane projection v 3 -v 2 (blue, upper graph), v 2 -v 1 (red) and v 3 -v 1 (blue) of the chaotic attractor, generated chaotic signal v 1 (red) and v 2 (blue) in the time domain, chaotic waveform v 1 (red) and v 2 (blue) in the frequency domain. Note that significant frequency components are concentrated in audio range.
Entropy 2018, 20, x 15 of 23 Figure 15. Orcad Pspice circuit simulation associated with network given in Figure 10: selected plane projection v1-v2 (blue, upper plot), v1-iL (red) and v2-iL (blue) of chaotic attractor, generated signal v1 (red) and v2 (blue) in time domain, chaotic waveform v1 (red) and v2 (blue) in frequency domain. Figure 15. Orcad Pspice circuit simulation associated with network given in Figure 10: selected plane projection v 1 -v 2 (blue, upper plot), v 1 -i L (red) and v 2 -i L (blue) of chaotic attractor, generated signal v 1 (red) and v 2 (blue) in time domain, chaotic waveform v 1 (red) and v 2 (blue) in frequency domain. Figure 16. Orcad Pspice circuit simulation associated with network given in Figure 10: plane projection vZ-iL2 (blue, upper plot), iL2-iL1 (red) and vZ-iL1 (blue) of the observed strange attractor, generated chaotic signal iL2 (red) and vZ (blue) in the time domain, chaotic waveform iL2 (red) and vZ (blue) in the frequency domain. Note that generated chaos is not affected by saturation levels of the active devices. Figure 16. Orcad Pspice circuit simulation associated with network given in Figure 10: plane projection v Z -i L2 (blue, upper plot), i L2 -i L1 (red) and v Z -i L1 (blue) of the observed strange attractor, generated chaotic signal i L2 (red) and v Z (blue) in the time domain, chaotic waveform i L2 (red) and v Z (blue) in the frequency domain. Note that generated chaos is not affected by saturation levels of the active devices.
Existence of observable strange attractors and numerically expected routing-to-chaos scenarios within dynamics of the fundamental MVMS has been proved experimentally; by its construction on the bread-board (see Figure 17) and consequent laboratory measurement using the analog oscilloscope. Due to its simpler realization and increased robustness, only circuitry illustrated by means of Figure 9 was decided for a real measurement. Selected chaotic waveforms in a time domain are provided by means of Figure 18. Independent voltage source V bias = 750 mV in series with R bias = 1 kΩ can be replaced by the positive supply voltage +5 V and the fixed series resistance R bias = 6.7 kΩ. Analogically, combinations V c1 = −130 mV, R x = 1 kΩ and V c2 = −650 mV, R y = 1 kΩ can be replaced by the negative supply voltage −5 V, R x = 38.5 kΩ and the same voltage −5 V, R y = 7.7 kΩ respectively.  Thus, there is no need to introduce three additional independent dc voltage sources into oscillator. However, to trace route-to-chaos scenarios, MVMS parameters c1 and c2 represented by voltages Vc1 and Vc2 have been considered as the variables. In practice, two voltage dividers based on  Thus, there is no need to introduce three additional independent dc voltage sources into oscillator. However, to trace route-to-chaos scenarios, MVMS parameters c1 and c2 represented by voltages Vc1 and Vc2 have been considered as the variables. In practice, two voltage dividers based on Thus, there is no need to introduce three additional independent dc voltage sources into oscillator. However, to trace route-to-chaos scenarios, MVMS parameters c 1 and c 2 represented by voltages V c1 and V c2 have been considered as the variables. In practice, two voltage dividers based on potentiometers followed by two voltage buffers (remaining part of fifth integrated circuit TL084) did the trick. Since extreme values of V c1 , V c2 voltages can be managed observed Monge projections given in Figure 19 (plane v 3 vs. v 2 ) originate both inside and outside prescribed bifurcation scheme pictured in Figure 6 (third and sixth plot). Of course, very small steps of the hand-swept parameters cannot be captured by the oscilloscope. Different plane projections, namely v 3 vs. v 1 , are shown in Figure 20. The provided screenshots are centered on captured state attractors not the origin of the state coordinates. Of course, the mentioned plane projections are not in full one-to-one correspondence with the theoretical results given in Figure 14, simply because transfer functions of the nonlinear two-ports cannot be defined precisely using imperfect discrete resistors; it does not get better shape even if the working resistors are replaced by the standard hand potentiometers. Finally, the sequence of periodic and chaotic windows with continuous change of the individual coefficients of the polynomials has been confirmed; roughly proving the correctness of bifurcation diagrams in Figure 6. During laboratory experiments, unexpected and a very interestingly-shaped strange attractors were identified (see Figure 21). Unfortunately, numerical values of the mathematical model parameters were not found, and reference state trajectories cannot be created. Initial conditions need not to be imposed on the outputs of the inverting integrators (capacitors need not to be pre-charged) since neighborhood of zero belongs into the basin of attraction for both strange attractors.

Conclusions
Existence of robust chaotic attractors in the smooth vector field of MVMS have been demonstrated in this paper. This work represents a significant extension of the discoveries discussed in [29]. Both Entropy 2018, 20, 697 21 of 23 referred chaotic attractors are self-excited attractors; hidden attractors were not sought after and, in the case of analyzed MVMS dynamics, remain a mystery. The proposed mathematical model: (1) together with the nonlinear functions (2) and parameters (7) or (8) can be also considered as a new chaotic dynamical system. This is still an up-to-date problem for many engineers and the topic of many recent scientific papers. However, algebraically much simpler dynamical systems that exhibit chaotic attractors exist; reading [42][43][44][45][46][47] is recommended.
Three different circuit realizations are presented. Each was verified by circuit simulation and the most robust implementation undergoes experimental measurement. Two designed oscillators utilize off-the-shelf active elements and can serve for various demonstrations.
The fundamental motivation of this work is to show that structurally stable strange attractors can be observed in smooth dynamical system naturally considered as nonlinear but with non-chaotic limit sets. However, designed autonomous chaotic oscillators can serve as the core circuits in many practical applications such as cryptography [48], spread-spectrum modulation techniques [49], useful signal masking [50,51], random number generators [52,53], etc.
This work also leaves several places for future research. For example, following interesting questions needs to be answered: Are there some hidden attractors? Are multi-scroll chaotic attractors observable if many RTD will be connected appropriately? Can MVMS generate chaos if values of the internal parameters are close to microelectronic memory cell fabricated in common technology?
Funding: Research described in this paper was financed by the National Sustainability Program under Grant LO1401. For the research, infrastructure of the SIX Center was used.

Conflicts of Interest:
The author declares no conflict of interest.