A comparison of different approaches to detect the transitions from regular to chaotic motions in SMA oscillator

It is well known that dynamical systems that include devices based on shape memory alloys (SMA) can exhibit a wide spectrum of responses: periodic, quasi-periodic and chaotic motions. In view of the different types of possible applications, it is important to recognize the qualitative features of the system vibrations. To this end, various methods have been proposed in the literature and evaluated in different conditions. In this work, a comparison between some of the available methods is proposed, focusing attention on their ability to detect the regular–chaotic and chaotic–regular transitions. The specific system under consideration is a thermomechanical SMA oscillator with superelastic behavior subject to harmonic excitation. The diagnostic methods compared are 0–1 test, maximum Lyapunov exponent and the recurrence indicators. The obtained results show that each method used is suitable for distinguishing between the regular and non-regular response of the SMA oscillator, so one of them can be chosen, taking into account, for example, the length and a sampling of the collected data.


Introduction
Shape memory alloys (SMAs) are smart materials whose macroscopic behavior is strongly influenced by the occurrence, at the microscale, of various solid phase transitions, the more typical being forward and reverse Austenite-Martensite transformations [1]. Consequently, the dynamical response of oscillators where the restoring force is provided by SMA-based devices is characterized by thermomechanical hysteresis.
The nonlinear dynamics of Shape Memory Alloys oscillators has been studied by several authors, using different constitutive models [2,3] The occurrence of chaotic responses under harmonic excitation has been systematically observed within various ranges of excitation amplitude and frequency [4,5].
Different nonlinear methods have been used to analyze the regular features of the system responses, including, among the others, standard frequency analysis, the maximal Lyapunov exponent, the method of wandering trajectories, the 0-1 test , and recurrence statistics [4][5][6][7][8][9].
Unfortunately, the above mentioned analyses were limited to detect the possible chaoticity of single trajectories selected within suitable ranges of parameters [9]. The various methods proposed in the literature are more or less effective in discerning between regular and non-regular responses. However, there are cases in which there is a complex transition from the regions where regular and chaotic responses are observed and the performances of the various methods have not been tested in such conditions.
In this work, the bifurcations leading to transitions from regular motion to chaos are studied in detail by changing the chosen bifurcation parameter, several diagnostic methods are compared in this context. Among these methods, the largest Lyapunov exponent and recurrence analysis are based on the study of the distance between the trajectories of the phase space [10,11] In the case of numerical simulations, all the coordinates from the phase space are available. Otherwise, a single coordinate has to be used to reconstruct the remaining ones. By virtue of Takens's embedding theorem [12], the topological properties of the reconstructed attractor will be the same as the original one.
It should be noted that the method of the delay embedding coordinate has been successfully applied to experimental data if the level of noisiness is not too high. In this work we use a single variable time series for the 0-1 test [13], a two-dimensional projection of the phase space for the maximum Lyapunov exponent and an embedding D-dimensional reconstructed phase space for recurrence quantification analysis.

System description
A thermomechanical SMA oscillator subjected to harmonic excitation was considered in [4]. This is a model inspired by a phenomenological approach to the Martensite-Austenite structural phase transition using an internal variable indicating a content of the Martensite phase. We adopted that model for our analysis and performed the calculations using Matlab ode45 integration schema which enable determination of higher order statistics [14]. The system of dynamical equations that determine the motion and the temperature evolution can be expressed in dimensionless form as follows: where x, v are displacement and velocity, n the internal variable that represents the fraction of Martensitic phase and # the system temperature, while # e is the environment dimensionless temperature. Besides the equation of motion, the system (Eq. 1) includes an equation of evolution for the fraction coordinate n and a heat equation governing the evolution of the temperature driven by the thermomechanical coupling arising from the peculiarities of the solid phase transformations by the constitutive functions K and H, given by Eqs. 3-5 after [4]. The normalized excitation amplitude and frequency are denoted by A and x, whereas the remaining system parameters determine the features of the hysteresis loop as well as the evolution of temperature (Table 1). The thermomechanical coupling (in Eqs. 1, 2) is defined by the constitutive functions K F;R and H F;R [4] where the subscripts F and R indicate, respectively, the expressions relative to the Forward (Austenite ! Martensite), upper plateaus, and Reverse (Martensite ! Austenite), lower plateaus, transformations. The first function couple K F;R depends on the martensitic fraction n completed by the W R;F functions that are directly dependent on n where n 0 is the martensitic fraction n at the beginning of the last phase transformation that occurred before current time s.
Furthermore, H F;R depend on n and # variables 3 Numerical simulations To set up the stage for the comparison of different methods for the analysis of the response, a set of 17 simulations with non-dimensional excitation amplitude A ¼ 1 at different frequencies is considered. Signals are composed of 200000 points sampled at 200 points per period. Using the Runge-Kutta method of the fourth order (ode45 in Matlab), we adopted a small integration step, obtaining a relatively large number of points for the excitation period, while for our analysis the sampling was changed. Being a thermomechanical system, besides displacement there is also a temperature variable and the simulations are carried out in non-isothermal conditions. In the following, results relative to both the displacement and the temperature evolutions are presented. The comparison between displacement and temperature response is interesting since temperature oscillated at double frequency ( Fig. 1). Note that the basic system, uncoupled to internal variables, has a resonance at dimensionless frequency x 0 ¼ 1 (see Eq. 1). The trajectories to be analyzed have an excitation frequency belonging to the sequential values from the interval x 2 ð0:147; 0:307Þ, (close to x 0 =4) chosen in such a way to illustrate a typical transition from regular to chaotic and then a reverse transition back to regular motion. Analyzing this diagram one can observe periodic motion for To give qualitative characterization of the different system responses, four cases are presented in Figs. 2, 3, 4 and 5 with phase portraits, Poincaré sections and corresponding power spectra. Note that a periodic motion is identified, being represented by a single point on the Poincaré map. It is noticeable that one of the chosen cases ( Fig. 4) represents non-regular solutions by a large number of Poincaré points and additional continuous frequency background in the power spectrum renormalized by excitation frequency. In all the cases, apart of the excitation renormalized frequency, one can observe a few odd superharmonics. The nonregular solution of x ¼ 0:237 presented in Fig. 4 shows a fairly large amplitude in terms of displacement. Interestingly, the corresponding frequency is close to the subharmonic resonance 1/4. Such proximity of frequencies in nonlinear systems may favor bifurcations and an occurrence of a chaotic solution.

Data preparation
The methods of the largest Lyapunov exponent and the recurrence analysis are based on the study of the distance between the trajectories in the phase space. In the case of numerical simulations, one has all the coordinates from the phase space available. Otherwise, one can typically use a single coordinate and use it to reconstruct the remaining ones. By virtue of Takens's embedding theorem [12], the topological properties of the reconstructed attractor will be the same as the original one, provided the system is smooth enough. It should be noted that the method of delay-embedding coordinates has been successfully applied to experimental data if the level of noise is not too high. Let x(t) be a single-time series, then the reconstructed vector u(t) will have form: uðtÞ ¼ fxðtÞ; xðt þ sÞ; ::: where s denotes the time delay and D denotes the embedding dimension. The values of the embedding parameters should be properly selected. In the case of a time delay, one can use the autocorrelation function and choose its zero, but keep in mind that it only takes into account linear correlations. To take into account nonlinearities that occur in the system, the mutual information function (MI) [15] is used more often and the s value is chosen as its first minimum. This choice avoids too strong correlations between nearby points when s is too small, as well as too weak correlations between points when s is too large. The embedding dimension m can be determined empirically by visualizing the attractor in increasing dimensions, so that its structure is preserved. More often, however, the function of the false nearest neighbors (FNN) is used, which examines whether nearby points in a fixed space remain neighbors in the space of an increased dimension. Zero of the FNN function corresponds to the embedding dimension, because all false neighbors disappear and no further increase of the dimension is necessary [16]. In this work we use a one-dimensional time series for the 0-1 test, a two-dimensional projection of the phase space for the maximal Lyapunov exponent and a D-dimensional reconstructed phase space for the recurrence quantification analysis. More details on individual methods will be presented in the following chapters. Such a procedure will allow us to check different scenarios depending on the method of data availability, both from numerical simulations as well as from experience. The above methods have been checked many times in the literature and successfully used to identify the type of dynamic response of nonlinear systems [6].

0-1 test
The first approach to the evaluation of the dynamic character of the system's response considered in this comparison is the so-called 0-1 test. The method, developed by Gottwald and Melbourne [13,17,18], is based on the statistical and spectral properties of a single time series and can be used to analyze the nonlinear dynamics of both model dynamic and experimental systems [19]. As a first step, the considered one-dimensional time series is mapped to an auxiliary two-dimensional space (p, q) by means of transformation: xðjÞ cos ðjcÞ; q c ðnÞ ¼ X n j¼1 xðjÞ sin ðjcÞ; ð7Þ where x is a time series and c 2 0; p ð Þ is a constant corresponds to fixed frequency in Fourier decomposition of the time series x. In the present case, as the input data for the 0-1 test we used the x coordinate from the phase space sampled with the s ¼ 10 as the time delay. This allowed us to avoid the problem of oversampling due to the small numerical integration step [6]. In a similar way, Poincaré points or any other mapping can be used. Exemplary representation of the auxiliary trajectory in (p, q) plane are shown in Fig. 6.
Comparing the corresponding figures for individual x values, one can notice their periodic character for x 2 0:147; 0:177; 0:297 f g and similarity to random walk for x ¼ 0:237. In other words, the auxiliary trajectory in new coordinates can be bounded (periodic motion) or unbounded (non-periodic motion). This fact can be examined using the value of the mean square displacement function MSD : ½p c ðj þ nÞ À p c ðjÞ 2 þ ½q c ðj þ nÞ À q c ðjÞ 2 ; where n corresponds to the total number of points (in practice the above limit is assumed by taking n ¼ n max ) and n max \ \N (usually N ¼ n=10). The values of this function can increase asymptotically over time in case of chaotic motion or be limited in time for regular motions and it can be quantified by means of the parameter : which thus plays the role of measureing of the regularity of motion. Alternatively, the correlation method can be used to determine another estimate of the parameter K: where X ¼ 1; ::; n max f g , M c ¼ M c ð1Þ; ::; M c ðn max Þ f g . In the above, the covariance covðX; YÞ and variance varðXÞ, for vectors X and Y of n max elements, and the corresponding mean valuesX andŶ, respectively, are defined as covðX; YÞ ¼ 1 n max X n max n¼1 ðX(n) ÀXÞðY(n) ÀŶÞ; varðXÞ ¼ covðX; XÞ: As we have noticed before, the values of K c depend on the choice of the value of c. If the parameter c would take values that are multiple of the resonant frequency incorrect results are obtained regardless of the dynamics of the system. In the practical implementation of the 0-1 test, to overcome this issue the test is repeated for a large number of values of c, typically 100 chosen within the interval c 2 ð0; pÞ, and the final value of the parameter K is given by the median of the corresponding values obtained in this way. If K % 0 then the dynamics of the system is periodic, and if K % 1 then the dynamics of the system is regular.
The test has then been applied for all the trajectories in the chosen set and Fig. 7 shows the variation of the parameter K in terms of the bifurcation parameter x. The results show that the test is able to detect properly both the forward and reverse transition from regular to chaotic motions. The transition is detected in a rather sharp way although there are points where the test gives intermediate values of the parameter K corresponding to zones close to the bifurcation points. Figure 7 shows clear differences in the dynamic response characteristics of the SMA oscillator. We see that for x ¼ 0:147; . . .; 0:177; 0:257; . . .; 0:307 f g the dynamics of the system are periodic, and for x ¼ 0:187; ::; 0:247 f gare chaotic. The solution obtained for x ¼ 0:187 is not identified clearly, but it is located near the bifurcation point and corresponds to the socalled weak chaos. It should be noted that values closer to unity are obtained in the limit n ! 1, and in our case n ¼ 4000 points.

Maximal Lyapunov exponent
The second method considered in this comparison is based on the computation of the Maximal Lapunov exponent L MAX , which is the most commonly used measure of non-periodicity, based on the divergence of nearby trajectories in the phase space. Each dynamical system has Lyapunov's exponents corresponding to the number of phase coordinates, but the sign of the largest of them unambiguously distinguishes the regular dynamics from chaos. In this work, we analyze the dynamics of the non-linear system of SMA oscillators using the 2d projection of the space from higher dimension (due to the hysteretic effects). To calculate the maximal Lyapunov exponent, we  [20] and Rosenstein [21], which developed it at the same time.
The method is based on the analysis of the geometric properties of the trajectories. Let Y(t) be a reference point on the trajectory. For each point Y(t) it is necessary to determine first the nearest points that are located as close as the average period that has to be determined before, for example by means of a Fourier transform: Yðt 1 Þ; . . .; Yðt u Þ such as jYðtÞ À Yðt j Þj\; j ¼ 1. . .u: Then, it is possible to define the average distance between the reference point and its nearest neighbors in a certain trajectory segment s via: Finally, the distances defined as above are averaged over a suitable interval hence determining the value of a function S, whose linear increase corresponds to the value of LMAX: This type of analysis has been carried out for all the trajectories in the considered set and Fig. 9 shows all values of L MAX in chosen range of the bifurcation parameter.
From Fig. 9 it follows that the method is able to discern between periodic and non-periodic dynamics, but the differences in values are not very large, especially in consideration of the scale of the plot. This may be caused by either discontinuities in the system or the procedure for determining the final value of L MAX as the slope of the linear section of the SðsÞ plot, which was based on sections of equal length.

Recurrence plots
The third method considered in this comparison is based on recurrence plots (RP). The method was originally proposed by Eckmann [22] and it is also based on testing the distance between points on trajectories in the phase space. If the distance is small enough then two points of nearby trajectories are marked as recurrence points. This can be applied to either numerical solutions of a system or signals recorded in experimental stands. In the second case, it is necessary to make a proper reconstruction of the phase space coordinates [12]. By means of these methods, we determined the embedding dimension m ¼ 4 and the time delay s ¼ 10 for the x coordinate of the system. The recurrence plot is constructed from the distance matrix R with its element R ij given by [22]: where is the threshold value. Norm jj Á jj can be Euclidean, maximal or any other ( here we choose Euclidean distance norm). In a RP, the elements 0 and 1 of the matrix R are represented by white and black dots. It is assumed that periodic systems are characterized by points forming long lines parallel to the main diagonal, non-periodic systems i.e. points forming short diagonal lines with isolated single recurrence points. Sample plots of the recurrence matrix for different dynamic responses of the SMA system are shown in Fig. 10. As one would expect, qualitative analysis using RP matrix visualization differs significantly for periodic and non-periodic vibrations. For x 2 f0:147; 0:177; 0:297g long diagonal lines are present, whereas they are not present for the RP calculated for x ¼ 0:237.
Later, Webber and Zbilut [23] and Marwan et al. [11,24] extended the concept of recurrence plots and on the basis of it proposed certain measures whose values identify the dynamics of the system in a quantitative way, calling it Recurrence Quantification Analysis (RQA). This method was successfully used for failure diagnosis [25,26]. Its application to a dynamical system with a SMA element was also proposed by Iwaniec et al. [8]. Unfortunately, these preliminary studies did not make any general conclusions.
First of all the RQA analysis includes RR which describes the ability of the system to visit the neighborhood of previous states, ie. measuring the number of recurrence points according to the formula The quantifiers of the recurrence are calculated using the probability of occurrence of diagonal p(l) or vertical p(v) lines (in this work we will only focus on those based on statistics of diagonal lines): P ðlÞ denotes the histogram of l lengths for a fixed threshold . Based on the probability p(l) subsequently the determinism DET can be calculated: where l min denotes the minimal value which should be chosen for a specific dynamical system. The next measures based on diagonal line statistics are the mean diagonal line (L) and the longest diagonal line (LMAX): where N l denotes the number of diagonal lines in the RP. To exclude the main diagonal from the statistics, we have used l min ¼ 4. The determinism, DET, is a measure of the system's predictability and determines the ratio of the recurrence points forming diagonal lines to all recurrence points. It is worth noting that in the case of complex periodic systems all points will be included in diagonal lines. L refers to the predictability time of the dynamical system and LMAX to the longest diagonal line, respectively. It should be emphasized that the key parameter in RP analysis is , which has to be chosen carefully. On the one hand, the threshold value should be small enough to capture the dynamics of the system in a short time scale, and the reverse side can not be too big, because then all of the points on the trajectory would be identified as recurrences. One way to choose is to select it so that the number of recurrence points RR is a few percent of all the points in the RR matrix. As before, each time series has been sampled every five displacement points As the RP method can be used for short time series, we have divided the entire time series with a length of 200000 points into 50 series with a length of 4000 points each. We calculated the final values of the RQA discriminators as the median of all 50 values and summarized them in Figs. 11 and 12 respectively.
By comparing the values of the RQA measures based on the diagonal line statistics, you can see that the biggest differences in changes in the periodicity of the system are visible in the values L (only for displacement- Fig. 11) and LMAX (for both displacement and temperature-Figs. 11 and 12, respectively). However, the RR and DET indicators can be useful for identifying mainly the same periodic attractors. To sum up, the transition from regular to chaotic and regular vibrations can be identified with RQA discriminants for the displacement and the temperature time series.

Discussion and final conclusion
In this paper, the regularity features of the nonlinear oscillations of a pseudoelastic shape memory alloy oscillator have been studied by means of three methods. Firstly, numerical simulations corresponding to regular-chaotic and chaotic-regular transitions were carried out (Fig. 1). In the case of the 0-1 test, the statistical and spectral properties of a single time series were examined asymptotically. In this regard, several aspects should be noted: oversampling may be a problem, which, due to strong correlation over short intervals, may give ambiguously interpretable results. On the other hand, nonlinear dynamical systems are characterized by correlations in different time scales, hence as part of preprocessing, one can resample the time series so that the input data are not too weak or too correlated. Autocorrelation (linear relationships)  or mutual information (non-linear relationships) functions that help you select a sampling step can be helpful. Furthermore, the value of c 0 can be a multiple of the system forcing frequency, which leads to ambiguity, hence the repeated calculations for different values of c 0 , whose statistics give the final value of the indicator, were performed. In the case of the largest Lapunov exponent determined from time series there are several parameters to look out for. First of all, the average orbital period of the time series, which we consider as the distance of nearby points on the trajectory.
In the case of regular vibrations, it is quite easily determinable, but in the case of irregular vibrations, this characteristic period does not exist but instabilities of close trajectories can be still determined. Secondly, the time to study the distance between adjacent points also cannot be too long, because after a certain time of evolution we will observe an exponential increase in distance regardless of the type of dynamics (in this case we adopted half the average orbital period), which does not change the fact that this is not universal value. Thirdly, using linear regression to determine the value of the L max exponent, we need to select the appropriate fit interval, which we have chosen subjectively by analyzing the values of the function SðsÞ. Like the 0-1 test, this is a symptomatic method, so the input data must be long enough to reflect the dynamics in both short and long time scales.
In the case of recurrence analysis, the most important parameter (in addition to the embedding parameters: time delay and space dimension) is the smallest value that marks points as recurrent. These are reflected in diagonal line patterns including the average length L, and the maximal length LMAX indicators. It is worth noting that this method does not require steady states and long time series, unlike the other two methods: the 0-1 test and the maximal Lyapunov exponent.
In final conclusion, we report that all three methods gave consistent results. They were able to detect bifurcations and different types of attractors (Test 0-1 (Fig. 7), Maximal Lyapunov Exponent (Fig. 9) and Recurrence Quantification Analysis (Figs. 11 and 12). The quantitative analyzes were also confirmed qualitatively by means of phase portraits, Poincaré sections and power spectra. In this scope, one can choose an appropriate method based on, for example, the length of the time series: Recurrence Plots-short time series with not so many cycles or Maximal Lyapunov Exponent/0-1 Test for longer simulations. For the studied cases, the 0-1 test gives a more precise detection of the transition, whereas the maximal Lyapunov exponent is less clear. This could be caused by discontinuities in the mathematical description of the system. On the other hand, the dynamical response of the system is also affected by internal variables which model the hysteretic material properties. It is noteworthy that the results of the temperature time series analysis correspond to the result obtained from the displacement/velocity time series (Figs. 7,9,11,12). This may be relevant in the context of an