## Abstract

Chaotic dynamical systems are essentially nonlinear and are highly sensitive to variations in initial conditions and process parameters. Chaos may appear both in natural (e.g., heartbeat rhythms and weather fluctuations) and human-engineered (e.g., thermo-fluid, urban traffic, and stock market) systems. For prediction and control of such systems, it is often necessary to be able to distinguish between non-chaotic and chaotic behavior; several methods exist to detect the presence (or absence) of chaos, specially in noisy signals. A dynamical system may exhibit multiple chaotic regimes, and apparently, there exist no methods, reported in open literature, to classify these regimes individually. This paper demonstrates an application of standard hidden Markov modeling (HMM), which is a commonly used supervised method, as a technique to classify multiple regimes from a time series of dynamical systems, where classified regimes could be chaotic or non-chaotic. The proposed HMM-based method of regime classification has been tested using numerical data obtained from several well-known chaotic dynamical systems (e.g., Hénon, forced Duffing, Rössler, and Lorenz attractor). It is apparently well-suited to serve as a bench mark for the development of alternative data-driven methods to enhance the performance (e.g., accuracy and computational speed) of regime classification in chaotic dynamical systems.

## 1 Introduction

Edward Lorenz introduced the concept of chaos in the form of the Butterfly Effect [1], where he suggested that the atmospheric phenomena are highly sensitive to exogenous perturbations. Chaotic dynamical systems are essentially nonlinear and are usually deterministic with an appearance of randomness, but having underlying patterns, feedback loops, and a structure of *repetitive order*, which are very sensitive to initial conditions and perturbations in the form of parametric changes and/or impulse inputs.

Chaos has been studied in great details and initially observed in several natural phenomena like weather fluctuations and climate changes [2]. The presence of chaos has also been seen in several other natural processes like heartbeat rhythms, which has helped physicians and scientists to make predictions on patient health [3], and anthropological studies to mathematically understand human evolution [4]. Chaotic nature also shows up in human-engineered processes such as urban traffic [5] and stock market [6]. Chaos is also seen in physical systems [7] starting from the simple double-pendulum going all the way up to thermo-fluid processes, electronic circuits, and even the fairground Tilt-a-Whirl.

Chaotic systems do follow a deterministic map (often with additive noise), where the initial conditions, process parameters, excitation inputs, and exogenous disturbances may significantly affect the time evolution of system responses. The *map* of a chaotic dynamical system may show significantly different nature for different values of process parameters. The state to which the system converges for a given set of parameters is called an *attractor* for that regime, which remains locally unchanged irrespective of the initial conditions or exogenous perturbations. There are essentially three types of attractors: (i) fixed-point attractors, where the system function converges to a single value, (ii) limit cycles, where the system shows a periodic nature, and (iii) strange attractors, where the system shows complex dynamical behavior that may be locally unstable and yet globally stable, implying that the system output may diverge but never significantly deviate from the attractor. This third kind of attractor is of much interest in the field of chaos.

It is often important to be able to distinguish noise from chaos and also to discriminate between the behaviors of a strange attractor and a limit cycle, which is often a challenge because the differences between these states may not be very apparent in the signal form, because of the presence of *hidden* dynamics and dependencies. Much work has been reported in literature in an attempt to make these distinctions based on time series data from various sources.

The 0–1 test method [8] has been introduced in the last decade to classify a given deterministic dynamical system as chaotic or non-chaotic. This method has been used extensively to study chaos in dynamical systems, including experimental plants (e.g., Ref. [9]). Djurovíc and Rubez̆ić [10] proposed a different approach by using short-time Fourier transform (STFT) to distinguish between chaotic signals and periodic oscillations. Recently, the concept of a simple decision tree structure is introduced by Toker et al. [11] to distinguish between non-chaotic and chaotic signals, which modifies the standard 0–1 test to show its applicability to several standard chaotic systems as well as to biological data.

Nearly all of the above techniques involve binary detection tests to study time series data from a system and decide whether it is chaotic or not. Few other researchers have attempted to classify the type of chaos but not with much rigor, because chaos occurs in complex dynamical processes, and it is not easy to generate a single test method to classify the various types of chaos occurring in the signal.

This study proposes regime classification in chaotic dynamical systems in the setting of a standard supervised data-driven method, called hidden Markov modeling (HMM) [12]. Several other well-known data-driven methods are potentially applicable for regime classification; examples are symbolic time series analysis (STSA) [13] and deep neural networks (DNN) using long short-term memory (LSTM) [14]. However, HMMs provide a good trade off with relatively smaller training times compared to that in LSTM-based neural networks with similar accuracy and may have higher accuracy than traditional STSA methods. Thus, in this study, the authors have chosen to demonstrate the efficacy of an HMM classifier in the context of performing the difficult task of discriminating multiple regimes in chaotic systems. To the best of authors’ knowledge, there is no reported literature that reports classification of different regimes of a chaotic dynamical system.

The primary objective here is not only to detect the presence of chaos but also to classify the operational regimes of the system, after the classifier is trained with the knowledge of the regimes (i.e., in a supervised fashion). As a classifier, HMMs are capable of identifying the current regime of the system, regardless of whether the regime is chaotic or non-chaotic. In the event that the dynamical system is needed to be controlled, this information could help in choosing a better control action to steer the system to a preferred state. Efficacy of the HMM-based regime classifier is demonstrated on several well-known chaotic dynamical systems that are sensitive to initial conditions and parametric perturbations and that have very different structures of attractors depending on their model parameters.

**Contributions**: Major contributions of the work reported in this study are as follows:

*Development of a HMM-based regime classifier for chaotic dynamical systems*: The underlying algorithm takes time-series signals as inputs for regime classification and does not require the information on whether the classified regimes are chaotic or not. The output is the identified regime based on previously trained regime class data.*Efficacy demonstration of the above HMM-based regime classifier*: The classification algorithm has been validated on time series data, generated from four chaotic dynamical systems—Hénon map, forced Duffing, Rössler, and Lorenz attractor.

**Organization**: The study is organized in five sections. Section 2 introduces the essential concepts of chaos along with a mathematical description of four different chaotic system models that are investigated here. Section 3 provides a brief background of HMM that is the backbone of the regime classification algorithm in this study. Section 4 presents the results for various cases of each of the four chaotic systems. Section 5 summarizes and concludes the reported work along with a few recommendations for future research.

## 2 Chaotic Dynamical Systems

Chaotic dynamical systems are essentially deterministic; they are sensitive to the initial conditions, only for chaotic regimes, and to system parameters for (possibly) all regimes. This study focuses on chaotic dynamical system models with parameters that can be varied to obtain multiple regimes of operation, including both chaotic and non-chaotic. This section presents the four chaotic dynamical systems used to generate time-series data from various operating regimes, obtained by varying the parameters (see Table 1 for details).

Case | Chaotic system | Fixed parameter values | Varying parameter values | Initial conditions for training data | Initial conditions for testing data | HMM error |
---|---|---|---|---|---|---|

I | Henon | b = 0.3 | a = 1 to 1.4 in increments of 0.05 | (x_{0}, y_{0}) = (1, 1) | (x_{0}, y_{0}) = (1, 1) | 0.97% |

I-A | Henon | b = 0.3 | a = 1 to 1.4 in increments of 0.05 | (x_{0}, y_{0}) = (1, 1) | (x_{0}, y_{0}) = (0, 0) | 12.24% |

II | Duffing | γ = 0.2, δ = 0.5,α = 1, β = −1 | ω = 0 to 2 in increments of 0.5 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0,0)$ | 0% |

II-A | Duffing | γ = 0.2, δ = 0.5,α = 1, β = −1 | ω = 0 to 2 in increments of 0.5 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0.1,0)$ | 0% |

III | Duffing | α = −1, β = 1,δ = 0.3, ω = 1.2 | γ = 0.2 to 0.60 in increments of 0.1 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0,0)$ | 4.48% |

IV | Duffing | α = −1, β = 1,δ = 0.3, ω = 1.2 | γ = 0.2, 0.28, 0.29, 0.37, 0.50 & 0.65 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0,0)$ | 1.11% |

V | Duffing | α = 1, β = 1,γ = 22, ω = 5 | δ = 0.1 to 0.35 in increments of 0.05 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0,0)$ | 30.77% |

VI | Rossler | b = 2, c = 4 | a = 0.25 to 0.55 in increments of 0.05 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | 1.69% |

VI-A | Rossler | b = 2, c = 4 | a = 0.25 to 0.55 in increments of 0.05 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (0, 1, 1) | 1.69% |

VII | Rossler | b = 2, c = 4 | a = 0.3547, 0.398, 0.410,0.458, 0.503, and 0.550 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | 1.17% |

VIII | Rossler | a = 0.2, c = 5.7 | b = 0.2 to 1.8 in increments of 0.4 | (x_{0}, y_{0}, z_{0}) = (0, 0, 0) | (x_{0}, y_{0}, z_{0}) = (0, 0, 0) | 0.85% |

IX | Rossler | a = 0.1, b = 0.1 | c = 5 to 45 in increments of 5 | (x_{0}, y_{0}, z_{0}) =(0.2, 0.3, 0.2) | (x_{0}, y_{0}, z_{0}) =(0.2, 0.3, 0.2) | 0.10% |

X | Lorenz | σ = 10, β = 8/3 | ρ = 50 to 250 in increments of 20 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | 0% |

X-A | Lorenz | σ = 10, β = 8/3 | ρ = 50 to 250 in increments of 20 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (0, 1, 1) | 0% |

XI | Lorenz | σ = 10, β = 8/3 | ρ = 13, 14, 15, 28, 99.96, 100.3,155, 193, 228, 229, 230.5, and 240 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | 0% |

Case | Chaotic system | Fixed parameter values | Varying parameter values | Initial conditions for training data | Initial conditions for testing data | HMM error |
---|---|---|---|---|---|---|

I | Henon | b = 0.3 | a = 1 to 1.4 in increments of 0.05 | (x_{0}, y_{0}) = (1, 1) | (x_{0}, y_{0}) = (1, 1) | 0.97% |

I-A | Henon | b = 0.3 | a = 1 to 1.4 in increments of 0.05 | (x_{0}, y_{0}) = (1, 1) | (x_{0}, y_{0}) = (0, 0) | 12.24% |

II | Duffing | γ = 0.2, δ = 0.5,α = 1, β = −1 | ω = 0 to 2 in increments of 0.5 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0,0)$ | 0% |

II-A | Duffing | γ = 0.2, δ = 0.5,α = 1, β = −1 | ω = 0 to 2 in increments of 0.5 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0.1,0)$ | 0% |

III | Duffing | α = −1, β = 1,δ = 0.3, ω = 1.2 | γ = 0.2 to 0.60 in increments of 0.1 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0,0)$ | 4.48% |

IV | Duffing | α = −1, β = 1,δ = 0.3, ω = 1.2 | γ = 0.2, 0.28, 0.29, 0.37, 0.50 & 0.65 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0,0)$ | 1.11% |

V | Duffing | α = 1, β = 1,γ = 22, ω = 5 | δ = 0.1 to 0.35 in increments of 0.05 | $(x0,x0\u02d9)=(0,0)$ | $(x0,x0\u02d9)=(0,0)$ | 30.77% |

VI | Rossler | b = 2, c = 4 | a = 0.25 to 0.55 in increments of 0.05 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | 1.69% |

VI-A | Rossler | b = 2, c = 4 | a = 0.25 to 0.55 in increments of 0.05 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (0, 1, 1) | 1.69% |

VII | Rossler | b = 2, c = 4 | a = 0.3547, 0.398, 0.410,0.458, 0.503, and 0.550 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | 1.17% |

VIII | Rossler | a = 0.2, c = 5.7 | b = 0.2 to 1.8 in increments of 0.4 | (x_{0}, y_{0}, z_{0}) = (0, 0, 0) | (x_{0}, y_{0}, z_{0}) = (0, 0, 0) | 0.85% |

IX | Rossler | a = 0.1, b = 0.1 | c = 5 to 45 in increments of 5 | (x_{0}, y_{0}, z_{0}) =(0.2, 0.3, 0.2) | (x_{0}, y_{0}, z_{0}) =(0.2, 0.3, 0.2) | 0.10% |

X | Lorenz | σ = 10, β = 8/3 | ρ = 50 to 250 in increments of 20 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | 0% |

X-A | Lorenz | σ = 10, β = 8/3 | ρ = 50 to 250 in increments of 20 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (0, 1, 1) | 0% |

XI | Lorenz | σ = 10, β = 8/3 | ρ = 13, 14, 15, 28, 99.96, 100.3,155, 193, 228, 229, 230.5, and 240 | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | (x_{0}, y_{0}, z_{0}) = (1, 1, 1) | 0% |

Figure 1 shows the bifurcation diagram for the Rössler system in Eq. (3) by varying the parameter *c*, while the other system parameters are kept fixed at *a* = 0.1 and *b* = 0.1. It is seen that as *c* is varied, the system behavior changes drastically, with different attractor basins for different ranges of *c*. Non-chaotic attractors cause the system to come to respective equilibria, irrespective of initial conditions; for chaotic attractors, the initial conditions may play a major role in the system evolution while the chaotic attractor remains the same irrespective of the initial conditions. Figure 2 shows the phase plots of signals for nine selected values of the parameter *c* along with a mention of the type of system regime that it belongs to. It is seen that changes of the parameter *c* from 8.5 to values in the range of 8.7–9 alters the system behavior from periodic with period 4, to periodic with period 8 to chaotic. This demonstrates possible sensitivity of a chaotic system to certain parameters. The plots in Figs. 1 and 2 show that it is not always easy to distinguish between the regimes by traditional analysis of time-series data, and perhaps that is why it is not reported in the current literature.

## 3 Hidden Markov Modeling

While the details of HMM are extensively reported in technical literature (e.g., Refs. [12,19]), this section briefly introduces the underlying concept for completeness of this study. HMMs have found applications primarily in speech recognition [20], time series classification [21] and even in image classification [22]. In all of these problems, the HMM method has shown high classification accuracy. Since HMM is a well-known and well-reported technique, the mathematical details are not presented here in detail for brevity; interested readers are referred to the cited references.

In the training phase of a HMM [12], time series data from each regime, i.e., classes $k=1,\u2026,K$, are used to learn the HMM. Subsequently, in the testing phase, the learned HMMs are compared to the HMM trained from the unknown regime. A succinct description of HMM is outlined below.

During training, the *Baum-Welch algorithm* [12,19] has been used to learn the HMM models of $K$ classes. For each class (regime), an ensemble of time-series windows is obtained from the available time-series data. This is done by first downsampling the long data with a downsampling rate *DS*, and then taking data windows of length *L* shifted by a distance *S*; the parameters, *DS*, *L*, and *S* are user-specified [23]. Then, the *Baum-Welch algorithm* [12,19] is applied to train the HMM *λ*_{k} which is a triplet *λ*_{k} = {*A*_{k}, *B*_{k}, *π*_{k}} [12], where

*A*≜ [*a*_{ij}] is the*N*×*N*state-transition probability matrix where*N*is the assumed finite number of hidden states belonging to a set of symbols $\u2208Q$:where $\u2211jaij=1\u2200i$.$aij=p(zt+1=qj|zt=qi):qi,qj\u2208Q$*B*≜ [*b*_{j}(*y*_{t})] is the probability density of the observation given the state:$bj(yt)=p(yt|zt=qj)$*π*≜ [*π*_{i}] is the probability distribution of the initial state*z*_{1}:*π*_{i}=*p*(*z*_{1}=*q*_{i}), where*π*is a 1 ×*N*vector with $\u2211i\pi i=1$.

This procedure is repeated for each of the $K$ classes.

*DS*,

*L*, and

*S*, as in the training phase, are provided as inputs to the algorithm. Then, the

*Forward Procedure*[12,19] of HMM is used to compute the log likelihood (

*L*

_{k}) of a given window of unknown time series data belonging to each of the $K$ classes. The final decision, as to which class the unknown data belongs, is made by selecting the class with the largest log likelihood as follows:

In this study, a continuous HMM formulation has been used, where the emission is assumed to follow a Gaussian mixture model with *M* Gaussian components and *N* hidden states. The algorithms and theory for all the above HMM procedures are available in the literature (e.g., Refs. [12,19]).

Although dynamical systems may have multiple regimes of operation, a real-life system would probably show only a few chaotic and periodic regimes. The HMM-based classification should be applicable to those cases as a subset of the dynamical systems demonstrated in this study.

## 4 Results and Discussions

This section presents the results of analysis for each of the four chaotic dynamical systems (see Eqs. (1)–(4)). Multiple cases are considered, each having a different set of fixed and varying parameters. For each case, time-series data are obtained numerically for the system parameters using a set of initial conditions of the states. In some cases, the initial conditions are varied to examine robustness of the HMM-based classifier to capture the attractor behavior when the initial conditions are perturbed in the chaotic system.

In the Hénon map (see Eq. (1)), the time evolution is accomplished by forward integration. For Duffing, Rössler and Lorenz systems (see Eqs. (2)–(4)), the signals are generated using the fourth-order Runge–Kutta formulation at a fixed time-step size of 1 × 10^{−3} s. Data of 5000 s have been generated corresponding to each parameter set, for both training and testing phases, yielding 5,000,000 data points per condition. Furthermore, the signals are artificially contaminated with 10% noise (amplitude ratio formulation), which corresponds to a signal-to-noise ratio (SNR) of 20 dB.

For both HMM training and testing (see Sec. 3), the data are first downsampled by *DS* = 100. This is necessary because the time-step of integration is intentionally made very small, which leads to oversampling. The processed data are then windowed with a moving window formulation (see Sec. 3) with a window size of *L* = 1000 and a shift value of *S* = 100 data points. The parameters *DS*, *L*, and *S* are held constant across all training and test cases.

The HMM parameters, namely, number of Gaussian components (*M*) and number of hidden states (*N*), also need to be assigned. The optimal values for each system of equations are slightly different. To allow for easy comparison, a common set of values are taken and then kept unchanged for all tested cases; these parameters are *M* = 4 and *N* = 30. Higher values of *N* degrade the performance due to over-fitting and needing more data for proper training. Also computational time increases with increasing *N*.

Table 1 lists the results for each case, which show both fixed and varying parameters for individual systems. The initial conditions for generation of training and testing data are also given. These parameters and their ranges are taken from cited research publications (e.g., Refs. [1,13,15–17,24,25]) and are reported here for reproducibility of the results.

For the same map, multiple cases have been presented with different fixed and varying parameters. The last (right hand) column in Table 1 provides the overall classification accuracy for various classes, where the varying data are used for both training and testing. The total error is defined as the ratio of the number of data windows erroneously classified, to the total number of data windows analyzed.

Majority of the cases are trained and tested using data from a range of values of the varying parameter, as obtained from literature or from the bifurcation diagrams of the dynamical systems. For cases I, II, VI, and X in Table 1, a second analysis is done, where the initial conditions for the testing data are made different from those of the training data. These cases are appended with “-A” in the first (left hand) column of Table 1. Certain cases, namely, IV, VII, and XI, have been analyzed using specific values of the varying parameter. These values are taken from literature as regions of interest and have been analyzed for illustrating the strengths and weaknesses of HMM-based classification. In general, it is seen that HMM-based classification yields low error rates (e.g., less than 2%). However, it is important to know why HMM-based classification performs somewhat poorly for cases I-A, III, and V, as explained below.

The poor performance in Case I-A is explained by the fact that, as these systems are chaotic, changing the initial conditions causes the system behavior to change significantly. This possibly leads to incorrect classification, especially for the regions that have *strange attractors*.

For case III, the chosen range of values of *γ* is such that several values correspond to signals coming from the same regime with possibly a minor difference in the signal form (e.g., *γ* = 0.3 and *γ* = 0.4 have very similar signals); this leads to poor classification. However, in the next case IV, the same fixed parameters are used, with the values of the varying parameter chosen to match regions where the signal nature varies. Here, a better classification accuracy is achieved.

In case V, the signals are from a region where the parameter *δ* ranges from 0.1 to 0.35. For *δ* = 0.1 to 0.3, the time signals are very similar, and the signal bifurcates to a different regime when *δ* ≈ 0.35, as seen in Figs. 3 and 4. Here too, HMM-based classification is unable to capture the difference between the signals when *δ* = 0.1 to 0.3. Thus, in this region of interest [13], the HMM method is apparently insufficient.

## 5 Summary, Conclusions, and Future Work

This study has demonstrated the applicability of the HMM-based time-series analysis in an attempt to classify regimes of operation in chaotic dynamical systems, which is not apparently explored much in literature. The proposed data-driven regime classification does not need the knowledge of whether the classified regime is chaotic or non-chaotic; this property has been demonstrated by considering various operational cases from four standard chaotic maps, namely, Hénon, Duffing, Rössler, and Lorenz attractor systems. The study has been conducted on global ranges as well as local ranges of interest near bifurcation zones.

The proposed HMM-based method of regime classification yields good accuracy, in general, and it performs sub-par only for a few specific cases, the possible causes for which are discussed in Sec. 4. Thus, this method, coupled with the standard tests for chaotic nature could help not only to identify chaos but also to classify chaotic “sub”-regimes.

The above results evince that the proposed HMM-based method is expected to be well-suited to serve as a bench mark for development of alternative data-driven methods of regime classification in chaotic dynamical systems.

While there are many areas of research to enhance the work reported in this study, the following topics are suggested to be conducted in the near future.

Development of other data-driven methods (such as Symbolic Time Series Analysis [13,26]) that can improve upon the accuracy of HMM-based classification, especially in

*closeby*regions of interest, typically seen just prior to bifurcation.Enhancement of robustness of data-driven regime classification to initial conditions of chaotic systems.

Application of the proposed classification method to real-life systems that display multiple chaotic regimes.

Development of a data-driven method to identify chaos type with no or partial knowledge of system dynamics (e.g., unsupervised or semi-supervised classification).

## Acknowledgment

The reported work has been supported in part by the U.S. Air Force Office of Scientific Research (AFOSR) under Grant No. FA9550-15-1-0400 in the area of dynamic data-driven application systems (DDDAS). Any opinions, findings and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the sponsoring agencies.

## Conflict of Interest

There are no conflicts of interest.