{'fraction_non_alphanumeric': 0.028766274023558587, 'fraction_numerical': 0.007398222773300269, 'mean_word_length': 4.919388379204893, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Autoencoder Based Iterative Modeling and Multivariate Time-Series Subsequence Clustering Algorithm

This paper introduces an algorithm for the detection of change-points and the identification of the corresponding subsequences in transient multivariate time-series data (MTSD). The analysis of such data has become more and more important due to the increase of availability in many industrial fields. Labeling, sorting or filtering highly transient measurement data for training condition based maintenance (CbM) models is cumbersome and error-prone. For some applications it can be sufficient to filter measurements by simple thresholds or finding change-points based on changes in mean value and variation. But a robust diagnosis of a component within a component group for example, which has a complex non-linear correlation between multiple sensor values, a simple approach would not be feasible. No meaningful and coherent measurement data which could be used for training a CbM model would emerge. Therefore, we introduce an algorithm which uses a recurrent neural network (RNN) based Autoencoder (AE) which is iteratively trained on incoming data. The scoring function uses the reconstruction error and latent space information. A model of the identified subsequence is saved and used for recognition of repeating subsequences as well as fast offline clustering. For evaluation, we propose a new similarity measure based on the curvature for a more intuitive time-series subsequence clustering metric. A comparison with seven other state-of-the-art algorithms and eight datasets shows the capability and the increased performance of our algorithm to cluster MTSD online and offline in conjunction with mechatronic systems.


Motivation
In the applications of machine diagnosis of mechatronic systems and the subfield condition-based maintenance (CbM), all supervised methods rely on high-quality labeled data [50,30]. But also, the unsupervised methods are more robust and provide better results when a "good" distribution of the data preexists. Unsupervised machine learning methods generally require independent and identically distributed (i.i.d) data to produce acceptable results. Providing data with this distribution in an industrial measurement setup is impractical, especially when continuously changing conditions cause time dependency. An option for a mechatronic system with different operating points is to measure the main operating points separately and create a diagnosis method for each of those individually. This requires a well-structured design and execution of experiments with a measurement labeling process. In a real world development environment for mechatronic system, where measurements are taken either automatically and/or manually and by many individuals and in different hardware and software development stages, providing consistently labeled and categorized data is a challenge.
Automatically labeling and categorizing multivariate time-series (MTS) data is therefore not only an alleviation but might be crucial for a successful CbM approach. As described above, labeled and categorized data is essential for training a model to represent a mechatronic system in a data driven approach. In the automotive sector where a lot of measurements occur at different operating points this is especially important. Some of these measurements are being recorded on a test bench in standard environment conditions with predefined operating points and for a given time (e.g., Worldwide harmonized Light Duty Test Cycle (WLTC)). Others can be in idle mode during waiting or preparation time for a longer trip or other measurements. Also, very transient episodes are existent (e.g., Real Driving Emissions (RDE)). All of these measurements do not necessarily have the same calibration of the underlying mechatronic system. To train a robust model of the mechatronic system, component group or a single component, a big effort has to be put in the design of the experiments alone, not to mention the experiments themselves. Therefore, a method of automatically labeling existing measurements is of advantage. Afterwards an automatic sorting of the labeled time sequences by statistical methods is possible, to enable a data driven mechatronic diagnosis approach.
Using advanced unsupervised approaches for CbM allows the data to be unlabeled (otherwise supervised methods could be used). In this case the labeling refers to the label of the condition (mechanical degradation) of the monitored system. When trying to diagnose mechatronic systems that have many operating points and are free to transfer in between those or are capable of totally transient operation modes, then a robust diagnosis of the actual condition of the mechatronic system is extremely challenging. An early and reliable (robust) diagnosis of a mechatronic system prevents accidents, enables optimal maintenance and increases uptime of machinery. Without the knowledge of the current condition of the system, fault prevention can only be done by predetermined maintenance intervals. Motivation is therefore to monitor the health condition of the mechatronic system as close as possible, resulting in the task of separating discrete sensory data into uniquely identifiable and recognizable segments or subsequences. This is beneficial to the performance of anomaly detection, because if all normally occurring subsequences are identified, the detection of abnormal or faulty subsequences is straightforward.

Introduction
When monitoring the health condition of a mechatronic system it is state of the art, to manually calibrate specific release conditions, during which the condition monitoring is enabled. This is done to exclude operating points which are very rare, too transient or are just not feasible for drawing conclusions about the condition of the mechatronic system. But even restricting the conditions on where to diagnose the machine (which already is reducing the probability of diagnosing the machine at all due to an operating state which is by chance outside the release conditions) cannot always help to improve the fault detection, identification and quantification of its magnitude. For example, in a mechatronic system with a complex nonlinear dependency of its subcomponents and its time dependency, "going in" or "going out" of the release conditions can result in very different system behavior. Comparing these two states does not lead to reliable conclusions for the mechatronic systems health condition.
Therefore, it is crucial for any monitoring strategy what kind of data sequences are used for training/calibration and validation. Manually screening, labeling and sorting data into comparable sequences is time-consuming, error-prone and cumbersome. Additionally, this is a decision process which requires expert and domain knowledge.
The new algorithm which we introduce in this work is capable of generating subsequence models from online streaming data which is processed sequentially. Any coherent subsequence that is identified can be recognized (clustered) if occurring again. Depending on multiple sensitivity calibration parameters, time-varying data points are associated and identified as a subsequence. The parameters determine the volatility or the strength of the affiliation required to be recognized as one time varying subsequence. These subsequence models can also be applied efficiently offline onto large existing datasets. During this prediction phase, the algorithm provides a vector of subsequence labels which were recognized as one from the training data. Depending on the calibration, it can also provide a label for unknown data which represents a phase where no pattern could be recognized. Otherwise, it finds the best fitting subsequence and labels it as that. The approach published in this work is currently only based on MTS input but could be adapted for a univariate input. It is a multivariate time-series sub-sequence discovery and identification method.
Our contribution is a new algorithm for online subsequence clustering of multivariate time-series data (MTSD) called "Autoencoder Based Iterative Modeling and Subsequence Clustering Algorithm (ABIMCA)" and a new metric to evaluate cluster algorithms focused on this task "Multivariate Time-Series Sub-Sequence Clustering Metric (MT3SCM)". We compare our algorithm with • seven other state-of-the-art algorithms • eight datasets from which six are publicly available and two are provided with our codebase • three widely used unsupervised clustering metrics • our own metric (MT3SCM) and its four components while varying the use of default algorithm parameters with optimized parameters on each algorithm and dataset via random grid search.

Related work
In this section we define the terminology used and its semantics to categorize our work within the large bibliography existing in this field and provide a selected list of related works and their ascendancy to this paper.

Terminology and Semantics
Numerous possibilities have been described for achieving our main goal of segmenting discrete time-series sensory data. Most approaches can be sorted into the following partially overlapping categories: time-series analysis [12], pattern recognition [10], temporal knowledge discovery [53], motif discovery [40], change-point detection [62], data clustering [31] or anomaly detection [19]. All those terms refer to methods or algorithms which could be used directly or indirectly to achieve our goal. Explicit description on each term or category can be found by the interested reader in the stated references.
To limit the scope of this work, we focus on data clustering which can be separated into six subcategories by the following groups of two: univariate and multivariate data, online and offline algorithms, variant and invariant data. The term "clustering" implies an unsupervised method. The equivalent supervised method would be called "classifi-

Multivariate Univariate
Online Offline Variant Invariant cation". In the relevant literature more and other distinctions are made, depending on the specific field and context. First, we describe time-series data (TSD) since this is the data format we focus on in this paper. Afterwards we explain the differences between the subcategories.
"A time-series is a sequence of observations taken sequentially in time." [12] We denote a time-series data point as an observation with exactly one connected timestamp. The timestamp is not a variable or feature.
univariate -multivariate If a single value or a scalar is the only variable of the data, then the data is univariate. It is the most basic format data can have. Considering TSD, a single temperature sensor with a timestamp would be univariate. Univariate TSD could also be interpreted as multivariate data of two dimensions, when taking the timestamp as another variable or feature.
online -offline The differentiation between online and offline algorithms or analysis of data is crucial. Offline refers to data analysis that is applied to all the data at once. Measurement data, for example, is available in one or multiple files or can be accessed via a previously filled database. Offline algorithms can therefore iterate and optimize their result based on a criterion applied to known data. Online analysis on the other hand, is applied sequentially. The algorithm needs to be able to function with a criterion that generalizes well on unknown data. One chunk or piece of data can be applied on the online algorithm without knowing the rest of the data. This approach cannot be as robust and accurate as an offline analysis, which is why most methods found in the literature are offline algorithms. Ideally the online algorithm learns every time new data is provided. For the sake of completeness, however, it should be noted that some online algorithms have to be pretrained offline and some algorithms referred to as offline, can be used sequentially. This depends on the underlying methods used. An offline algorithms' purpose is not to be used online.
Online algorithms, nevertheless, can be used offline.

dependent -independent TSD (time variant -time invariant)
Depending on the field of study the specific terminology of time dependency can differ. We want to emphasize on the common accepted assumption "An intrinsic feature of a time-series is that, typically, adjacent observations are dependent" ([12, S. 1]). Time dependency characterizes TSD, where a consecutive observation has some connection with its predecessor. In some fields a connection is not necessarily given for two data points in a database that are the closest to each other regarding their timestamp. The dependency of adjacent observations is self-evident, when collecting sensor values of a mechatronic system from an experimental rig, for example. We therefore use the term "time dependency" in the context of a dynamical system and extend it to a time-variant system in the terminology of control systems engineering. Independent TSD would be where the variance and the average are invariant along the time (stationary) and ergodic.
We position our work in the subcategory of online clustering of dependent multivariate time-series data. As of now we refer to a time-series as a sequence of dependent observations with a constant sample-rate. A "discrete series of consecutive data points" as a subset of this time-series is synonymously referred to as pattern, motif, sequence, operating state, state change, between change points, subsequence, episode or segment, among others. In this paper we will use the term subsequence (using terminology of [4]).
The algorithms we use for comparison are listed in Table 1. All are online clustering algorithms, that can be used for time-series clustering, implementations are publicly available in the Python programming language (see library column in Table 1) and are well established and tested. The Balanced Iterative Reducing and Clustering using Hierarchies (BIRCH) [67] algorithm is based on a clustering features (CF) tree with the CF as a triple of the number of data points, linear sum and the squared sum. This CF tree is built dynamically. It was also one of the earliest algorithms capable of online clustering. The Bayesian Online Change-point Detection (BOCPD) [1] algorithm as the name suggests, uses Bayesian methods to detect change-points (CPs) online. Since this algorithm only detects CPs, we manipulated the result to be able to interpret every CP as the beginning of a new cluster. This algorithm starts in our comparison with the limitation of not being able to recognize a previously seen cluster. The Stream Clustering Framework (CluStream) [2] algorithm is based on extended CF from BIRCH, following a k-means algorithm. The Density-based Stream Clustering (DBSTREAM) [28] algorithm is based on the Self Organizing density-based clustering over data Stream (SOStream) [29] and uses a shared density graph to capture the density between micro-clusters. The Density-Based Clustering over an Evolving Data Stream with Noise (DenStream) [15] algorithm is an extension of Density-Based Spatial Clustering of Applications with Noise (DBSCAN) [25] which uses a damped window model of CF to create core-micro-clusters and outlier-micro-clusters. The Mini-Batch K-Means (MiniBatchKMeans) [56] algorithm proposes "the use of mini-batch optimization for kmeans clustering" ( [56]) to improve the k-means optimization problem. The STREAMKMeans [45] algorithm uses an adaption of the original STREAM algorithm from [46]. Replacing the k-median subroutine LSEARCH by an incremental k-means algorithm. More information and comparison of most of the used algorithms can be found in [17,26] As described in section 2, the focus of this publication is on the online multivariate time dependent subsequence clustering using recurrent neural network (RNN) based Autoencoder (AE). The number of algorithms within this scope is limited compared to the number of clustering algorithms in general. The following approaches use at least some of those prerequisites. [32] emphasize the term segmentation for an offline sliding window and bottom-up algorithm.
Others are converting the time-series into a Markov chain (MC) and then using a Bayesian method to cluster the MCs [51]. Referring them as episodes. Here the data needs to be discretized into bins of equal length. Reference [64] use manually selected characteristics (e.g., kurtosis, skewness and frequency) for clustering univariate TSD.
Others are using the Augmented Dickey-Fuller test to evaluate time-series stationarity and perform a segmentation based on this [58]. In [43] dynamic latent variables from a vector autoregression (VAR) model in combination with a principal component analysis (PCA) is used for segmenting industrial TSD. Reference [18] show the advantage of an embedding approach as well, by introducing a PCA and a Vanilla-AE CP detection method with the restriction of focusing on multivariate power grid data.
Focused on transfer learning, [38] introduce an adversarial approach for domain adaption using a stacked AE. Offline convolutional sparse AE used for supervised sequence classification was done by [8] and adapted by [9] for unsupervised motif mining. Other AE based papers are, for example, by [20] who use a mixture of AEs for image and text clustering. Stacked AE and k-means for offline clustering is done by [11] without considering time dependency. Showing the combination of GRU-based AE and MTS for anomaly detection is done by [27]. Reference [21] apply the sliding windows approach on CNN-based AE for Anomaly Detection of Industrial Robots. A similar approach using AE for MTS segmentation is published in [37]. The focus there is on change point detection using latent space variables only and no clustering or identification of the subsequences is done.
Clustering is strongly depending on the data and task provided, as Jain et al. expressed in 1999: ". . . ; each new clustering algorithm performs slightly better than the existing ones on a specific distribution of patterns." ([31, S. 268]) Therefore, we try to apply the algorithms on multiple different MTS datasets and compute different metrics for comparison. Large efforts are made for making datasets available to the scientific community and the public to improve comparability and reproducibility by universities or governmental institutions [24,22]. For this paper we focus on data with multivariate quantitative features with continuous values. For a list of the datasets see Table 3 and a brief description is given in section 5 on page 7. Evaluating the performance of a clustering algorithm can be done with two different approaches. If external knowledge about the ground truth of each data point and its cluster is known, then so-called external measures can be applied. If no ground truth is available, internal measures need to suffice. Many external measures exist, like the well-known F1-score (based on the effectiveness measure by [63]). With the large number of data available and working in the context of transient machine behavior with the focus on finding internal states of the system, acquiring or providing the ground truth is time-consuming, error-prone and cumbersome (as described in section 2). "The definition of clusters depends on the user, the domain, and it is subjective." ([3, S. 30]). We therefore use internal measures for comparing our approach. Those internal measures commonly rely on a similarity measure of the actual data which being clustered. Using the provided label from the cluster algorithm many possible combinations of, for example, the density of all points in a cluster and the distance to the center of the nearest other cluster are available. Thorough work on metric comparison and similarity measures has been done [35,57,3]. Most of those measures are based on simple distances and densities computed for each data point but do not take time dependency into consideration. Because of this, we found, that for the use case described in this paper, the commonly used clustering evaluation measures are not well suited for "time-series clustering evaluation measures".
In section 6 we introduce an approach for similarity measures which considers time dependency in combination with well-established clustering metrics (see Table 2). Table 2: Metrics used for time-series clustering comparison. Implementations used from [49] valuation value metric silhouette [54] Ratio of distance to its own cluster and distance to the nearest cluster center {s ∈ R : −1 ≤ s ≤ 1} the higher, the better calinski-harabasz [13] Ratio of between-cluster variance and the within-cluster variance {s ∈ R : 0 ≤ s ≤ ∞} the higher, the better davies-bouldin [23] Ratio of cluster size and between-cluster distance {s ∈ R : 0 ≤ s ≤ ∞} the lower, the better

Definitions and Restrictions
In the following section we define in more detail our data, together with the restrictions of our environment. Considering online clustering, we can refer to our TSD as continuously incoming data or streaming data. This data is considered multivariate when the dimension (number of sensors or features) of the data stream d > 1. When we denote one value of one feature as x, we have at time step t the following feature vector: Whereas the natural numbers include zero {0, 1, 2, . . .} = N. A complete measurement sequence with n number of time steps using x from Equation (1) as so, X ∈ R d×n . With time dependency consideration, it is reasonable to denote a sliding window of the streaming data, considering Equation (1) and n the number of samples already collected as: Let's also assume, that within the measurement data X there exist subsequences S j which satisfy our requirements of non-overlapping and variable length. For the indexes of our subsequences, we denote where u is the number of identified subsequences in X. The subsequence then is a continuous sampling from X for a small period of time steps with consecutive data points and the length of m. The subsequence length is usually much smaller than the length of the full measurement data m n.
For each subsequence with the index j we have a first time step index q j = q j,start and a last q j + m = q j,end for which we do not allow overlapping This results in our uniquely identified non-overlapping set of subsequences Clustering these uniquely identified subsequences results in recognizing reoccurring subsequences and combining them into a subset of all subsequences which results in the following cluster set C with the cluster index or unique cluster label For the output of a clustering algorithm at time t we denote the scalar value y t as our label or designated subsequence identification. For evaluation purposes a clustering for a time-series produces a label array y for all time steps: Furthermore, it is a requirement, that the streaming data provided can be applied to a numerical differentiation algorithm. Therefore, a constant sample rate is necessary and in case of strong noise, filtering or smoothing of the data should be applied by a preprocessing step. Also, there mustn't be missing values and extreme outliers need to be removed. In our use case we assume that some knowledge about the incoming data exists, so that an estimate of the variance and the mean of the variable can be performed for standardization.

Datasets
All datasets used for comparison in this work are described briefly in this section and listed in Table 3. They all contain quantitative features with continuous values. For further use of the datasets, no missing values exist, the data is continuous and was standardized for the algorithms but not for the metric computations. No other preprocessing like smoothening or filtering was performed.
The bee-waggle dataset [47] contains movement of bees in a hive captured with a vision-based tracker. The first two features are the x and y coordinates of the bee added with the sine and the cosine function applied to the heading angle.
The cmapss dataset is a "dataset of run-to-failure trajectories for a small fleet of aircraft engines under realistic flight conditions" [7] with 18 features.
The eigen-worms dataset [66] contains measurements of worm motion. Preprocessing extracted six features, which represent the amplitudes along six previously identified base shapes of the worms The hydraulic dataset [55] is obtained from a hydraulic test rig with measuring 17 process values such as pressures, volume flows and temperatures.
Lorenz-attractor refers to a synthetic dataset which is calculated using a system of the three coupled ordinary differential equations which represent a hydrodynamic system:Ẋ = s(Y − X);Ẏ = rX − Y − XZ;Ż = XY − bZ with parameters used s = 10, r = 28 and b = 2.667 (see Figure 2). "In these equations X is proportional to the intensity of the convective motion, while Y is proportional to the temperature difference between the ascending and descending currents, similar signs of X and Y denoting that warm fluid is rising, and cold fluid is descending." [41] The mocap or The Motion Capture Database (MOCAP) dataset [16] contains 93 features from human motion captured with markers. The occupancy dataset [14] is a measurement of sensory data in an office with the following sensors: temperature, humidity, the derived humidity ratio, light and CO2.

Multivariate Time-Series Sub-Sequence Clustering Metric (MT3SCM)
As emphasized in section 2 and section 3, to our knowledge none of the existing clustering metric exists, that takes the time space variations like curvature, acceleration or torsion in a multidimensional space into consideration. We believe using these curve parameters, is an intuitive method to measure similarities between mechatronic system state changes or subsequences in MTSD in general (in regard to the restrictions in section 4).
Our MT3SCM score consists of three main components.
The weighted curvature consistency (cc w ), the silhouette location based (s L ) and the silhouette curve-parameter based (s P ). When making the attempt of clustering TSD, it is subjective and domain specific. Nevertheless, we try to take the intuitive approach of treating MTSD as space curves and use the parameterization as a similarity measure. This is done in two different ways. First, we create new features by computing the curve parameters sample by sample (e.g., curvature, torsion, acceleration) and determine their standard deviation for each cluster. Our hypothesis is, that with a low standard deviation of the curve parameters inside a cluster, the actions of a mechatronic system in this cluster are similar. We call this the curvature consistency (cc) (see Equation (24) used in line 14 in algorithm 1). The second procedure is to apply these newly computed features, which are computed to scalar values per subsequence, onto a well-established internal clustering metric, the silhouette score [54] (see Table 2).
The computation of the cc comprises the calculation of the curvature κ and the torsion τ at every time step t with x t .
Whereas e 1 is the unit tangent vector (or first Frenet vector), e 2 is the unit normal vector (or second Frenet vector) and e 3 is the unit binormal vector (or third Frenet vector) which are defined as: From which we can also derive the speed v = ẋ t and the acceleration a = ẍ t . Figure 3 shows exemplarily the curvature κ, torsion τ , speed v and acceleration a for the first part of the thomas-attractor dataset. Afterwards the cc is calculated per cluster i ∈ I, by taking the empirical standard deviation for each curve parameter (exemplarily for κ in Equation (20) with the set of subsequence indexes J i within our cluster i). The arithmetic mean (Equation (21)) of the standard deviations for the curvature κ, torsion τ and the acceleration a results in the final cc per cluster (see Equation (22)).
The cc w , is directly derived from the cc per cluster, by weighting it with the number of data points per cluster i ∈ I Equation (24).
The calculation of the scores s P and s L is different to the standard estimation of the silhouette score, which is shown in Equation (25) and originally based on every data point of the time-series X and the assigned cluster label array y: Our s P is the silhouette score derived from our previously computed curve parameters per subsequence per cluster as well as the standard deviation of those and the number of data points per subsequence. with The s L uses the silhouette score based on the median valuex dij of a subsequences original feature space per feature d. with The main idea of this approach is to combine three main parts inside one metric. First incentive is to reward a low standard deviation of the curve parameters in between a cluster (accomplished by cc). Second, to benchmark the clusters spatial separation based on the new feature space (curve parameters, accomplished by s P ). And third, to benchmark the clusters spatial separation based on the median of the subsequence in the original feature space (accomplished by s L ). The proposed algorithm for this new metrics computation is described in algorithm 1.

Evaluation
For computational tests, we manually created a "perfect" synthetic dataset with respect to our metric (see Figure 4). Figure 4 (a) shows the original synthetic dataset, where the subsequences in cluster 1 are a helix along the increasing x axis. For cluster 2 the subsequences are a straight movement, with quadratic decreasing distances along the y axis. Cluster 3 is representing a helix along the decreasing x axis but with a different resolution than cluster 1. Cluster 4 is, along with cluster 2, a straight movement with quadratic increasing distances along the y axis. This cycle is repeated six times. Figure 4  K ← GETCURVEPARAMETERSFORALLDATA(X) 5: y unique ← FINDUNIQUECLUSTERIDS(y) 6: for i in y unique do 7: for j in s do 10: cc w ← WEIGHTEDAVERAGE(C, n pc ) Compute weighted average curvature consistency (cc w ) from cc i with number of points per cluster 18: s L ← SILHOUETTECOMPUTATION(L, y unique ) Compute the silhouette coefficient using the center positions of each identified subsequence 19: s P ← SILHOUETTECOMPUTATION(P, y unique ) Compute the silhouette coefficient with the curve parameters 20: score ← (cc w + s L + s P )/3 21: return score The final score 22: end procedure To show the benefit of the new feature space, we applied the agglomerative clustering 10 not on the original lorenzattractor dataset but on the newly computed feature space based on curvature, torsion and acceleration (see Figure 5) The metric values for Figure 5 (b) show a high cc w and a decent s P value for the low number of 10 clusters specified.
To further evaluate our metric, we used the lorenz-attractor and the thomas-attractor dataset (see Table 3) and applied an agglomerative clustering, a time-series k-means clustering as well as a random subsequence clustering. Varying the number of clusters and some algorithm specific parameters. Afterwards the metrics calinski-harabasz, daviesbouldin and silhouette scores were computed and compared to our new metric MT3SCM. From these results we derived a correlation matrix (see Figure 6). The cc and the cc w are clearly related due to their direct combination. The positive correlation between the internal components to the overall MT3SCM score is obvious. We see a clear positive correlation to the silhouette score which is evident due to the internal use of this metric. Interestingly, the correlation between the cc w and the s P is negative. This is due to the types of datasets and algorithms we used. Because with higher number of clusters we theoretically expect a better cc because of the lower standard deviation by chance. On the other hand, the more clusters exist, the more likely a similar curve parameter between the clusters exists and therefore creates a new feature space with overlapping clusters, which results in a low s P score. This can be retraced within the subfigures of Figure 7. The low correlation between the calinski-harabasz and the davies-bouldin scores supports our point that the available clustering metrics are not well suited to be used for time-series clustering evaluation measures. Figure 7 shows examples where the agglomerative clustering was applied on the lorenz-attractor dataset (part of the data used for the correlation matrix Figure 6). It can be seen that the agglomerative clustering on the original dataset is not an optimal cluster algorithm, when comparing the metrics to Figure 5 (b). Comparing Figure 5 (b) and Figure 7 (d) we can see a similar MT3SCM score but very different standard metrics scores. The similar MT3SCM score is based on the much higher number of clusters and equally distributed subsequence length in Figure 7 (d), which results in a high cc w value as well as a good spatial separation (s L ), which is compensating the low s P value due to the similar curve parameters of the clusters. Figure 5 (b) however, also has a very high cc w value with a good s P value reaching a similar MT3SCM score but with a fifth of the number of clusters. How our metric handles random clustering with critical scenarios, is shown in Figure 8. The Python code and a more detailed evaluation are publicly available at [33] Table 4 for metric comparison of the following subplots. (a) New curve parameters feature space computed from the Lorenzattractor dataset with labels from agglomerative clustering (b) Lorenz-attractor dataset with labels from agglomerative clustering on the new curve parameters feature space (c) New feature space from the curve parameters extracted from each subsequence. Standard metrics scores computed with new feature space (d) New feature space from the centers (median value) of each subsequence. Standard metrics scores computed with new feature space

Conclusion
We have described a more suitable similarity measure for dependent TSD. After showing how to compute our metric and evaluated on different datasets its use case and effectiveness. Further we will use this metric in addition to the standard metrics to evaluate our proposed online time-series clustering algorithm which is described in section 7 on the next page Figure 6: Own metric (MT3SCM) correlation analysis. Own metric and its four subcomponents (curvature consistency (cc), weighted curvature consistency (ccw), silhouette location based (sL), silhouette curve-parameter based (sP )) correlation to calinskiharabasz, davies-bouldin and silhouette score for random, agglomerative and k-means clustering on lorenz and thomas-attractor dataset.  In this section we describe the concept of our time-series clustering approach in detail. Afterwards, we apply our algorithm onto the datasets described in section 5 on page 7 and present the results.

Method
As described in [39] a key component in a time-series clustering algorithm is the similarity function to quantify the clustering criteria. Common similarity functions used are distance measures like euclidean distance or some kind of correlation coefficients like Pearson's correlation coefficient. Those are also used for static data clustering algorithms. More suitable for time-series clustering are similarity functions like Dynamic Time Warping (DTW) distance, short time-series (STS) distance [44] or considering space curves like we introduced in section 6.
In this work we analyzed an approach which is data driven, based on unsupervised machine learning algorithms and has online capabilities (see Figure 11). Our approach uses a RNN based AE to generate scores which are used as similarity measures. Specifically, the experiments in this work were performed using a pytorch [48] implementation of a bidirectional one-layer gated recurrent unit (GRU) RNN with a hidden size of the input dimensions minus one h = d − 1. Other prerequisites regarding the dataset and preprocessing are described in section 5 on page 7 and section 4 on page 6.  Table 4 for metric comparison of the following subplots. Subfigures (a) (b) and (c) all have a similarly low MT3SCM score compared to Figure 5 (b) but considerably good standard metric scores. Subfigure (d) can achieve a relatively high MT3SCM score due to the high number of clusters and the resulting good ccw and sL value which compensates the low sP value.  Table 4 for metric comparison of the following subplots (a) Own metric and all of its subcomponents are around zero, as desired. Calinskiharabasz value is low and davies-bouldin is high, which also indicate a "bad" clustering (b) Longer random subsequences also generate a MT3SCM result around zero. Calinski-harabasz and davies-bouldin scores are stronger influenced by the subsequence length (c) Own metric and all of its subcomponents are around zero, as desired. Calinski-harabasz value is low and davies-bouldin is high, which also indicate a "bad" clustering (d) As seen for the lorenz-attractor data in (b), longer subsequences have a high impact on calinski-harabasz and davies-bouldin scores The main procedure of the approach is as follows: The incoming data is taken as a sliding window W t at the current time t with length ζ of past time steps and number of features d. This matrix W t ∈ R d×ζ is used for the input of, what we call, the Base Autoencoder (BAE). The key element of our algorithm is, that this BAE's parameters are not constant but being adapted iteratively with a stochastic gradient descent (SGD) optimization method for each new incoming sliding window. For this training of the BAE, we use a slight adaption of the sparse AE loss function L from [52] with a basic regularization term or sparsity penalty Ω where h = f (x) is the encoders output or latent space. The sparsity penalty we denote as: with the penalty factor λ = 1e−10 and the latent center constant c lc = 0.5. The mean squared error (MSE) is which results in the final loss computation L(x, h, r) where the first part is the MSE between the input matrix W t and the reconstruction W t and the second part is the penalty of the latent space deviation.
To determine if a subsequence is recognized at the current time step, we denote the scoring function SF as follows whereas the weighting factor in current implementation is c f w = 1 and latent center constant is c lc = 0.5. It utilizes the reconstruction error as well as the deviation of the latent space with a doubled emphasis on the latent space deviation due to its dependence in the loss function as well as the scoring function which includes the loss again (see Equation (34)). In combination with a threshold, the score is used to determine when a recognizable subsequence is present.
The bottom row of Figure 9 shows, that a subsequence is present, when the BAE's score (blue line) is below the horizontal black line (s b <= η). If a subsequence is present, a copy of the BAE is made and its parameters are frozen and associated with this specific pattern of a subsequence. These copies of the BAE, which we call Subsequence Autoencoder (SAE), are used to recognize previously seen subsequences using the same scoring function. A concept drawing of the approach is shown in Figure 11. The algorithm is described in pseudocode in algorithm 2.
The functionality can be retraced considering Figures 9 and 10. This example shows the algorithm applied on a threedimensional synthetic data set. The input data consists of four different operation points with small white noise. The sequence of the four subsequences is repeated once. The other rows are described in the caption of Figure 9. It can be seen that the algorithm needs a few time steps to adapt to the current subsequence until it is recognized as such.
Recognizing a previously identified subsequence, however, is almost instantaneous. It is apparent, how the calibration of the thresholds η (horizontal black line) and ζ (horizontal gray line) are crucial. The necessary time steps to adapt to a current subsequence can be altered by the calibration of the learning rate α and the number of BAE's training cycles per time step ω. A faster recognition of a subsequence has the drawback of the algorithm being very sensitive and therefore identifying even small changes of the input as a new subsequence. A strategy could be to calibrate the algorithm first to be rather insensitive and cluster the time-series in major subsequences. These can then be further clustered with a more sensitive calibration. This procedure can be repeated until the required degree of granularity is achieved. For a streaming application multiple runs of this procedure could also be applied in parallel and combined into a cluster tree.
Algorithm 2 ABIMCA. Using the following parameters: α: learning rate, ω: number of base model training cycle per step, η: subsequence detection score threshold, ρ: subsequence recognition score threshold, window length: ζ. Also, we denote θ S as a list of subsequence model parameters, s as an array of scores for all existing subsequence models, W t the reconstruction of the sliding window input, 1: procedure ABIMCA(W t ) Sliding window data W t ∈ R d×ζ 2: c S ← 0 Initialize subsequence counter 4: θ b ← SPARSE(sparsity = 0.1) Initialize base model parameters 5: while W t do New input available 6: for j ← 1, ω do Iterative base model train loop 7: 10: Update base model parameters 11: end for 12: s ←GETSUBSEQUENCESCORES(W t , θ S ) 15: if min(s) < ρ then Subsequence recognized 16: S ID ← arg min(s) Set ID to recognized subsequence index 17: else No subsequence recognized 18: if s b <= η then Score below new subsequence threshold 19: θ S .append(θ b ) Append current base model parameters to list of subsequence models parameters 20: c s ← c s + 1 increase subsequence counter 21: S ID ← c s Set ID to new subsequence index 22: else In transition 23: S ID ← 0 Set ID to unknown 24: end if 25: end if 26: end while 27: end procedure

Evaluation
For the evaluation study of our algorithm, we chose eight different MTS datasets (see Table 3) from which six are publicly available and two are provided with our codebase [34], seven other state-of-the-art algorithms (see Table 1) and three widely used unsupervised clustering metrics (see Table 2). Each algorithm has been applied to each dataset with default parameters. Additionally, we performed a hyperparameter search for each algorithm based on a random grid search of 300 samples. The parameter boundaries for this hyperparameter search are listed in Table 15. Overall, 19 264 experiments were run.  For a better overview of the results, we chose to compare every algorithm to the "MiniBatchKMeans" algorithm and counted the number of times they performed better. Table 5 shows the results for the hyperparameter search and the number of outperformances of each algorithm compared to the "MiniBatchKMeans" algorithm. Table 6 shows the same results with default parameter settings for each algorithm. We can see, that in sum and in two of the metrics our algorithm beats state-of-the-art algorithms. The full list of results is attached in Tables 7 to 14 on pages 21-20 As cited before, every algorithm performs differently on the specific distribution of patterns and the hyperparameter search was a simple random grid search of "only" 300 samples, so these results unlikely represent the optimal solution for each algorithm on each dataset. Nevertheless, we demonstrate, that the algorithm we present in this work, is highly effective of detecting subsequences online in a MTS.

Limitations and Discussion
For the evaluation of the segmentation of MTSD, we introduced a new metric which is based on space-curve parameters in the feature space. Due to the wide variety of fields, use cases and applications, this falls into place for some applications and uses cases but not for all. The calculation of these space-curve parameters is sensitive to outliers and smoothness and questionable for steady-state conditions or non-moving point clouds in the feature space. We have implemented specific numerical boundary limits for computing the derivatives of the data in these states, but it needs to be considered and evaluated if this is compatible with the application. Because of the outlier sensitivity we use the mean value of these parameters as well as their standard deviation. A low-pass filter for very noisy data should be considered before applying it to the metric. Attention is called for, when the data is scaled or standardized. This effects the actual space curve parameters, since a constant curvature is likely not constant anymore after scaling. The metric also tends to reward short subsequences who only occur once. Due to the mean value of the curve parameters the subsequence separation appears to be good, but the variance of one large subsequence is high. This needs to be compensated or prevented more and will be part of future analysis and improvement of the metric. It might also be reasonable to introduce weighting factors for the three parts of our metric in Equation (12) to consider domain specific emphasis on a rather spatial or curve parameter separation requirement.
Regarding our clustering method, calibration of the main thresholds (η, ρ) needs special attention. In combination with the learning rate, they mainly influence the "sensitivity" of the segmentation process. Using only the reconstruction error or only the representation deviation can be beneficial for different use cases and complexities (e.g., changing the weighting factor in Equation (34)). Advantages are that no additional information about the data for offline clustering needs to be saved. All necessary information for clustering data after training is the parameters of the subsequence specific AEs. It is a completely unsupervised method which can cluster online data. In the context of CbM the once identified subsequence AEs can be used for deviation quantification of the underlying system. This can be used for deterioration analysis and maintenance strategies. Further investigations for improving the ABIMCA method would be to explore different kind of AEs like feedforward neural network (FNN), CNN or a combination of such. Also, a VAE could be reasonable depending on the underlying process. Future work should analyze the effect of reducing the latent space dimension by multiple factors of the input dimension (when input dimension is very high). This could reduce computation costs and improve representation learning without performance loss. A detailed analysis of the optimal default parameters or a generic automatic calibration depending on some statistics of the expected input could increase performance and decrease calibration efforts.

Conclusion
In this paper we have introduced the Autoencoder Based Iterative Modeling and Subsequence Clustering Algorithm (ABIMCA) which is a deep learning method to separate multivariate time-series data (MTSD) into subsequences. It is beneficial in a variety of fields, to cluster MTSD into smaller segments or subsequences in an unsupervised manner. The ability to filter measurement data based on specific subsequences can improve downstream development products such as anomaly detection or machine diagnosis in condition-based maintenance (CbM) strategies. Our algorithm is specifically useful for MTSD generated by a mechatronic system in a transient environment. It can be used offline as well as online for streaming data. It utilizes recurrent neural network (RNN) based Autoencoders (AE) by iteratively training a Base Autoencoder (BAE), generating a segmentation score and saving the intermediate parameters of the BAE to recognize previously identified subsequences. By comparing our algorithm with seven other algorithms on eight different publicly available datasets using four different unsupervised metrics (from which we introduced one ourselves), we have shown that our algorithm outperforms state-of-the-art algorithms. Our unsupervised metric introduced (Multivariate Time-Series Sub-Sequence Clustering Metric (MT3SCM)), is an attempt to use a more intuitive similarity measure based on the curvature and other space-curve parameters of the spanned feature space. Additionally, all our code is open source and publicly available for benchmarking.