Phase space warping: nonlinear time-series analysis for slowly drifting systems

D Chelidze, J.P Cusumano


A new general dynamical systems approach to data analysis is presented that allows one to track slowly evolving variables responsible for non-stationarity in a fast subsystem. The method is based on the idea of phase space warping, which refers to the small distortions in the fast subsystem's phase space that results from the slow drift, and uses short-time reference model prediction error as its primary measurement of this phenomenon. The basic theory is presented and the issues associated with its implementation in a practical algorithm are discussed. A vector-tracking version of the procedure, based on smooth orthogonal decomposition analysis, is applied to the study of a nonlinear vibrating beam experiment in which a crack propagates to complete fracture. Our method shows that the damage evolution is governed by a scalar process, and we are able to give real-time estimates of the current damage state and identify the governing damage evolution model. Using a final recursive estimation step based on this model, the time to failure is continuously and accurately predicted well in advance of actual failure.

1. Introduction

An important class of contemporary problems in experimental nonlinear dynamics involves the study of slowly drifting systems in which the slow dynamics cannot be treated as merely some undesirable ‘non-stationarity’ because it is, on the contrary, the very important phenomenon that one desires to analyse. One such problem, namely the study of damage evolution, performance degradation and failure prediction in engineered systems, has been the primary motivation for the work presented here.

The field of machinery condition monitoring has traditionally emphasized the development of robust threshold-based damage detectors capable of classifying and ‘predicting’ impending failures. However, we are interested in developing methods capable of actually tracking the evolution of damage states (or equivalently, ‘health states’) of systems in real time. Such methods are necessary if one wishes to track damage, or other unwanted system changes, as it evolves from its earliest stages. Furthermore, a reliable method for estimating the damage states is a prerequisite for true prognostics, which we interpret in the strictest sense as providing continuously updated estimates of the remaining life. More fundamentally, however, such experimental methods will provide never-before-seen time-series datasets that will shed light on the still poorly understood multi-scale dynamics of failure processes, which currently can only be studied empirically by statistical approaches based on post-failure specimen evaluation.

In this paper, we present a new approach to track the evolution of a slowly evolving ‘hidden’ state using only the measurements from a ‘fast’ subsystem for which it is the source of non-stationarity on long time-scales. More specifically, we consider hierarchical dynamical systems of the form,Embedded Image(1.1a)Embedded Image(1.1b)Embedded Image(1.1c)where Embedded Image is the fast dynamic variable (the directly observable state), Embedded Image is the slow dynamic variable which is assumed not to be directly accessible, the parameter vector Embedded Image is a function of Embedded Image, t is the time, the rate constant Embedded Image is defined as the time-scale separation between the fast dynamics and the slow ‘drift’, and y is a scalar quantity derived only from the fast variable x via the smooth function h.

Models of the form of equations (1.1a)–(1.1c) have been proposed for systems with evolving damage (Cusumano & Chatterjee 2000), in which Embedded Image represents the slowly evolving damage variable, and in this paper, we take this as the primary interpretation. Indeed, in the remainder of this paper, we will treat ‘slow variable’ and ‘damage variable’ as synonymous. However, we emphasize that this situation arises in a variety of other contexts as well, and we expect the approach described here to be applicable quite generally for systems with adequate time-scale separation.

One might be tempted to eliminate the function Embedded Image in the fast subsystem of equation (1.1a) and treat Embedded Image itself as a slowly drifting parameter vector. However, by using this notation, we wish to emphasize that the effect of damage on the fast subsystem will be mediated by its effect on the constitutive parameters required to describe it. Obviously, if the ‘damage’ does not affect any of the model parameters at the scale of dynamical observation, it cannot be said to have an effect on the system. Thus, it can be concluded that either the model was inadequate for the purpose of studying damage evolution or the proposed damage process was, in reality, harmless. Beyond this practical issue, the appearance of the damage variable Embedded Image only through the function Embedded Image is a consequence of the hierarchical nature of equations (1.1a)–(1.1c), which are not merely hierarchical in time (due to the time-scale separation induced by ϵ), but also often in space. In particular, for material damage processes, the slow variable Embedded Image can be viewed as a generalized state variable characterizing the material of the fast subsystem at the microstructural level, and Embedded Image is a mapping from the microstructural scale to the macrostructural scale of the fast subsystem.

All diagnostic methods rely on some multivariate statistic, or feature vector, in which the changes are monitored and used as input to a prognostic layer when it exists. The heuristic basis of the work describe here is that, given dynamical data from any functioning system, one expects intuitively that the most information one might hope to extract from the data would be a model capable of generating the observed time-series itself. Hence, this work uses nonlinear time-series analysis to generate a dynamical feature vector based on the short-time reference model prediction (STRMP) error, and discusses conditions under which this particular type of feature vector can function as a tracking function for an underlying drift process.

Dynamical systems based time-series analysis has matured to the point that its basic methods, which fall under the general framework of phase space reconstruction, are now well-described in widely available texts (e.g. Tong 1990; Weigend & Gershenfeld 1994; Abarbanel 1996; Kantz & Schreiber 2004), though the field continues to develop rapidly. However, there are two significant limitations to standard time-series analysis approaches when one intends to study systems as in equations (1.1a)–(1.1c). First, time-series analysis usually focuses on the temporal ‘mesoscale’, corresponding to the problem-specific time-scale at which the data should be sampled, what we are here calling the ‘fast’ subsystem of equations (1.1a)–(1.1c). Indeed, a significant part of the craft of time-series modelling resides on the various methods for identifying this very time-scale.

Components of the data that are significantly faster than this scale are treated as ‘noise’, and are typically removed by low-pass filtering or averaging, whereas components that are significantly slower are treated as non-stationarity and are eliminated, if possible, by high-pass filtering or differentiation. In regard to high-frequency noise, our approach will not be different. However, in this case, the non-stationarity is driven by the slow process of equations (1.1a)–(1.1c), and it is precisely the drifting variable that we wish to examine using only measurements at the mesoscale of equations (1.1a)–(1.1c). Thus, it is necessary to modify standard time-series analysis methods to make this possible.

A second limitation of standard nonlinear time-series approaches, closely related to the first, is that they have been attractor based, a fact that reflects the origins of developments of the last few decades in the study of chaotic systems. Hence, standard approaches have tended to emphasize the estimation of invariant quantities, such as Lyapunov exponents or generalized dimensions. However, while these quantities (which are just other kinds of feature vectors) may be adequate for detecting sudden changes in a system, they do not provide continuous tracking of the evolving slow variable. Hence, our approach focuses not on the long-time asymptotic system behaviour, but rather on the short-time dynamics of the fast subsystem.

To study systems of the form of equations (1.1a)–(1.1c) experimentally, the next section introduces the concept of phase space warping (PSW), an approach to phase space reconstruction aimed at studying the small distortions that occur in the fast subsystem's vector field as a result of an underlying slowly evolving process. After summarizing the basic theory, the fundamental issues that must be addressed for algorithmic implementation are discussed, and the algorithm is described. Next, we show the results by applying the method to a system with evolving material damage. In Chelidze et al. (2002) and Cusumano et al. (2002), a version of the method capable of tracking scalar slow variables has been applied to the study of a system in which the potential energy is perturbed by a battery-powered electromagnet: ‘failure’ of the system in that case corresponded to complete discharge of the battery. It was shown that the tracking metric output by the algorithm was related in a 1–1, approximately linear fashion to the scalar generalized damage variable, which in that case was the open circuit battery voltage.

Here, we apply the method to a vibrating beam nonlinear oscillator in which a crack propagates to complete fracture. The analysis presented uses a new version of the PSW algorithm that can track vector slow variables and that makes use of the recently developed method of the smooth orthogonal decomposition (SOD; Chelidze & Zhou 2006). Using this approach, we are able to show that the damage evolution is dominated by a scalar process, which is the first time such a result has been obtained. Again, the tracking metric is shown to provide a 1–1 relationship with an independent measurement of the damage.

Using empirical damage evolution models and recursive filtering, the tracking metric can be used to predict the remaining useful life. This approach has been applied to the battery experiment in Chelidze & Cusumano (2004), and here we apply it to the fracture experiment. We demonstrate that accurate real-time estimates of current damage state and time to failure (TTF) can be made well in advance of actual failures.

2. Phase space warping

In equations (1.1a)–(1.1c), we assume that U and V are the compact subsets and that the phase space is the Cartesian product Embedded Image, where T is the manifold of which t is an element. We also assume that T itself is a compact manifold such as, e.g. a p-torus corresponding to the vector field f in equation (1.1a) being p-quasi-periodic in time. The key issue is that the dynamics of equation (1.1a) must take place in a region of the extended fast phase space Embedded Image that is diffeomorphic (via delay coordinate embedding) with a compact subset Embedded Image where d is the embedding dimension.

Our goal is to use only experimental measurements of the fast variable Embedded Image to track, and ultimately predict, the slow variable ϕ, which is not directly measurable. In general, equations (1.1a)–(1.1c) have solution Embedded Image and Embedded Image, where Embedded Image is the initial condition, and Embedded Image is the prediction time.

We wish to compare the fast dynamics starting at some current time Embedded Image to what the dynamics would have been had it been in the original or ‘reference’ condition. Letting the constant Embedded Image to represent the reference value of the damage variable, we refer to the solution Embedded Image as the reference model. We then define the time Embedded Image ahead STRMP error starting at time Embedded Image asEmbedded Image(2.1)The arguments in the last line are arranged to emphasize that Embedded Image is a map that takes the current value of the slow variable, Embedded Image, into Embedded Image, and that we will treat Embedded Image, Embedded Image and Embedded Image as experimentally manipulable parameters (in addition to the ‘fixed’ parameters Embedded Image and ϵ). We refer to the vector Embedded Image as a tracking function, if it is an injective mapping (preferably linear) from V into Embedded Image, i.e. if every point in its range corresponds to a unique point Embedded Image. The conditions under which this is true are presented below.

It is of course true that equation (2.1) can be related to changes in the fast vector field f, since taking the derivative with respect Embedded Image givesEmbedded Image(2.2)An alternative, therefore, would be to estimate Embedded Image, which is a measure of the distortion, or warping, in the vector field f caused by changes in the slow variable, which motivated giving the name PSW for our approach. However, it is difficult to measure vector fields directly in experiments, and thus in what follows, we use the flow form of the tracking function in equation (2.1), which gives the total amount of PSW at any given point Embedded Image.

We are only considering a fixed, finite prediction time and so we expect the solutions to be smooth with respect to both initial conditions and parameters, even in the presence of bifurcations. We require that the prediction time be ‘short’, i.e. Embedded Image. For any given initial time Embedded Image, the short prediction time allows us to treat the fast subsystem equation (1.1a) as ‘quasi-static’, since the damage variable Embedded Image will be approximately constant during this time interval. More precisely, we can use regular perturbations to expand the terms in equation (2.1) in a power series about ϵ=0, as Embedded Image where Embedded Image. Substitution into equation (2.1) givesEmbedded Image(2.3)

We further expand the leading term in equation (2.3) in a Taylor series about Embedded Image, which then gives the STRMP error asEmbedded Image(2.4)in which we have suppressed the fact that the derivative matrices are evaluated at Embedded Image. We expect smaller values of ϵ to allow relatively larger prediction times while maintaining quasi-stationarity. In practice, however, the prediction time is indeed very small, so that the Embedded Image remainder term in equation (2.4) is expected to be truly negligible.

Equation (2.4) shows that Embedded Image is to leading order an affine transformation of the slowly drifting variable,Embedded Image(2.5)where Embedded Image is Embedded Image, and Embedded Image is Embedded Image. We have dropped the zero subscript from Embedded Image since the above discussion is true for any future value of the slow variable over sufficiently short-time intervals. However, the zero subscript is still required on t and x since the matrices of equation (2.5) will depend on the selected initial point in the extended fast-time phase space, Embedded Image, used to compute the STRMP error.

A necessary condition for the transformation of equation (2.5) to be a tracking function is that the matrix Embedded Image has maximal rank. This requirement on C can be interpreted as a linear observability condition for the slow variable. (However, note that even when the linear condition fails, higher order observability may be possible.) Thus, in the neighbourhood of a point in the extended phase space Embedded Image for which C has maximal rank, for ϵ and Embedded Image both sufficiently small, we can conclude that the STRMP error allows us to unambiguously track the changes in the slow variable Embedded Image using only measurements of the fast variable x. Furthermore, we see that under ideal circumstances, this tracking function can be expected to be approximately linear.

Finally, we note that given the fixed values of Embedded Image, Embedded Image, Embedded Image and Embedded Image, we can take the time derivative of equation (2.5) so thatEmbedded Image(2.6)From equation (2.2), we see that this slow evolution can also be interpreted as giving the rate at which PSW, i.e. distortion in the fast vector field f, is occurring, and that this rate is controlled to leading order by the slow vector field of equation (1.1b). Given the Embedded Image nature of the vector field, we expect the solution to equation (2.6) to be well approximated by a slow flow differential equation Embedded Image on the time-scale of Embedded Image, whereEmbedded Image(2.7)which does not depend on x. We will not attempt to use this fact analytically: not only would such a calculation typically be very cumbersome, we rarely expect to have sufficient a priori knowledge about the system to do so. However, it is useful when interpreting results, and when determining the necessary model structure of the drift process a posteriori, which needed to use the tracking output for prognostics.

3. Implementation of PSW algorithms

In practice, the existence of a tracking function is used as a hypothesis. Analytical derivation is not only typically unfeasible, but data-based estimation, via regression, is also impossible because independent measurements of the true damage variable are rarely available—indeed, that is the entire motivation for the method.

In general, under the assumptions presented earlier, the output of the tracking function cannot, strictly speaking, be used as an observer for the damage state variable, but instead gives the damage state to within an unknown, approximately affine transformation. This situation is analogous to that found with delay coordinate embedding. The embedding process creates a new observable or ‘coordinate’ that can be shown to be smoothly related to the true state variable given the satisfaction of certain hypotheses (see Sauer et al. 1991, for a complete discussion), but this ‘coordinate transformation’ is not, in general, explicitly available. However, the situation with PSW is somewhat better; in cases where independent measurements of the slow variable are available, the tracking function can, in effect, be ‘calibrated’, so that the exact transformation can be determined.

(a) Tracking metrics

In experiments, the tracking function provided by the STRMP error, equation (2.5), is difficult to apply directly because C and c are parameterized by the position in the extended fast phase space, (Embedded Image). For most applications, it is difficult or impossible to repeatedly start the fast subsystem equation (1.1a) from these same initial conditions. Furthermore, it is not necessarily desirable to do so even when it is possible, since one might accidentally use a point for which the singular values of C are small, yielding a tracker with very poor sensitivity.

Thus, one should use many values of Embedded Image from some ensemble Embedded Image to deal with the problem of initial state repeatability, as well as to increase the robustness of the method. Then for every fixed Embedded Image, the STRMP error Embedded Image can be thought of as a random variable determined by random maps of the form of equation (2.5), with each map corresponding to a single element of Embedded Image. Given such an ensemble, one can examine the multivariate statistics Embedded Image in order to obtain information about the structure of the slow phase space. We refer to statistics obtained in this way as tracking metrics.

For example, one might use the ensemble of tracking functions to replace the map of equation (2.5) by Embedded Image, where the angled brackets indicate the (possibly weighted) ensemble average over Embedded Image. A more general approach, that makes better use of all the information contained in the multivariate statistics of Embedded Image, is presented in §3d. Furthermore, if it is reasonable to assume, given the physics of the problem, that the slow variable is both monotonic in time and scalar, then an even simpler approach is to consider only the magnitude of Embedded Image, which gives to leading orderEmbedded Image(3.1)where Embedded Image and Embedded Image. Different versions of this scalar-tracking metric were used in the original versions of the PSW method (Chelidze et al. 2002; Cusumano et al. 2002; Chelidze & Cusumano 2004), as described in more detail in §3c.

(b) Application to experimental time-series

Experimental time-series of the fast state variable x are obtained indirectly via the function equation (1.1c), which represents the measurement system. In this paper, we assume that the data are sampled with a uniform time interval Embedded Image, and that N point data records are collected over the time period Embedded Image. Under these assumptions, during the collection of any given data record, Embedded Image is constant to Embedded Image. Over the course of an experiment, a total of Embedded Image data records are collected.

The measured scalar time-series Embedded Image are used to reconstruct a fast-time trajectory for each record of data using delay coordinate embedding (Sauer et al. 1991). In particular, the scalar time-series is converted to a series of vectors Embedded Image (for Embedded Image), where Embedded Image, τ is a suitable delay and d is the appropriate embedding dimension. Embedding parameters, τ and d, are usually determined using the first minimum of the average mutual information (Fraser & Swinney 1986) and method of false nearest neighbours (Kennel et al. 1992), respectively.

Of course, we do not have access to the analytical solution X used to define the STRMP error via equation (2.1). Instead, in the reconstructed phase space, the state vectors are assumed to be governed by a map of the formEmbedded Image(3.2)where Embedded Image is generally nonlinear and needs to be determined from the data. The experimental tracking function is then defined for any point in the reconstructed phase space, following equation (2.1), asEmbedded Image(3.3)where as before, Embedded Image is the reference state of the damage (slow) variable and Embedded Image is its current value. However, the first term on the right-hand side of equation (3.3) is already available from the current data record, and so we have, after equation (3.2),Embedded Image(3.4)for every point on the reconstructed phase space trajectory. Thus, one only needs to construct the experimental reference model Embedded Image for the ‘undamaged’ system, that is using data from the reference condition (usually the first data record).

A variety of nonlinear models can be used to estimate the map P in equation (3.4). In work to date, we have successfully used locally linear maps of the formEmbedded Image(3.5)where Embedded Image is a Embedded Image matrix and Embedded Image is a Embedded Image vector. The parameters of the local linear models are determined for each point y in the current data record using regression on the Embedded Image nearest neighbours of that point and their future states for data taken in the reference condition. Then the PSW tracking function can be written asEmbedded Image(3.6)in which Embedded Image is the modelling error and Embedded Image is the estimated tracking function, which can be determined experimentally and is given byEmbedded Image(3.7)It is clear from equation (3.6) that, in general, the quality of the estimated tracking function depends strongly on the modelling error being small. Furthermore, we see that an important consequence of having to estimate the dynamics of the reference model is that one obtains another source of fluctuations in tracking functions associated with each point in the reconstructed phase space. In addition to the variations caused by the actual dependency of the tracking function on the chosen point in the fast phase space, there is also a variation caused by local changes in the quality of the reference model fit. The reference model fit error, in turn, has two sources: one comes from the type of model used and the accuracy with which it can represent the actual dynamics and the other is experimental noise. All of these fluctuations must be addressed for successful implementation of the PSW approach.

(c) Scalar tracking

In the previous work, the authors have used some form of the scalar-tracking metric shown in equation (3.1). Several ensemble-averaging strategies were advanced. For example, in Chelidze et al. (2002) a scalar-tracking metricEmbedded Image(3.8)was proposed, with weighting function Embedded Image given byEmbedded Image(3.9)where Embedded Image is the radius of the ball centred on y containing all the nearest neighbours used in estimating the parameters of equation (3.5) and Embedded Image is the point-wise dimension about y in the reference data set. This type of averaging compensates for fluctuations in the accuracy of the reference model built using a fixed number of nearest neighbours to y; tracking functions estimated using neighbourhoods with smaller radii are emphasized over those with larger radii. By scaling the radii using the point-wise dimension, this is effectively the same as favouring those regions of the reconstructed phase space with a higher density of points.

In subsequent work (Chelidze & Cusumano 2004), the ensemble-averaged scalar-tracking function was estimated usingEmbedded Image(3.10)where F was a simple Kalman filter based on a single step constant predictor. In the Kalman filter algorithm, changes in accuracy of the reference model, even in the absence of noise, were treated as ‘process noise’, whereas experimental noise was treated as ‘measurement noise’. The Kalman filter was found to provide more effective reduction of the model error Embedded Image.

(d) Vector tracking: smooth orthogonal decomposition

While scalar tracking has been shown to be very effective in several cases to date, there are several reasons not to be satisfied with it, in general. The most obvious problem is that a scalar-tracking metric can only track a scalar slow process. Indeed, if one does not have a priori knowledge of the slow process dimensionality, but instead wishes to determine it experimentally, a scalar-tracking metric will be of little use. While one might imagine performing phase space reconstruction on a scalar-tracking time-series to study dimensionality, there would be several potential problems with such an approach, not the least of which is that the tracking data evolves on the slowest time-scale, and thus one will typically not have sufficient data to estimate dimensionality in this way.

The maximal rank requirement for tracking discussed in the derivation of equation (2.5) suggests that one can track a slow variable with any dimension less than or equal to that of the fast variable. Therefore, an examination of the tracking function data that makes use of its full dimensionality could be used to provide an unprecedented tool for studying the dimensionality of slow processes such as damage evolution. In fact, the method described below may make it possible to track damage processes with dimension exceeding that of the fast variable.

The basic idea behind the multivariate tracking procedure is to partition the reconstructed reference phase space into a fixed number of local neighbourhoods and evaluate the expected value of the tracking function in each of them. A similar approach has been used by Chelidze (2004) to track a two-dimensional battery discharge process. An advantage of this partitioning method is that, in addition to addressing the dimensionality issue, it directly addresses the phase space variability of the tracking function (discussed in detail in §3a) since the behaviour of the tracking functions in each neighbourhood is monitored throughout each experiment. Furthermore, the approach leads to an inherently spatial analysis that allows one, in principle, to identify how regions of the fast phase space are coupled to the drift process.

We begin by defining Embedded Image disjoint regions Embedded Image, where the subscript of Embedded Image emphasizes that the regions are defined with respect to points in the reference phase space. Then the average of the estimated tracking function Embedded Image is found in the ith region byEmbedded Image(3.11)where Embedded Image is the number of points in Embedded Image. All of the local average vectors in data record j are combined into one multi-dimensional feature vector:Embedded Image(3.12)where the semicolon in the above expression indicates column-wise concatenation of each column vector Embedded Image, giving an Embedded Image vector Embedded Image for each data record, where Embedded Image and Embedded Image.

For a small change in slow-time variables, based on the above discussion and equation (2.5), it is conjectured that there is an affine projection Embedded Image that maps the proposed feature vector onto the damage state in the jth data record:Embedded Image(3.13)where V is an Embedded Image matrix, and v is an Embedded Image vector. Affine projection parameters, V and v, can be determined in a least-squares sense if independent measurements of the damage state are available (Chelidze 2004). However, as discussed previously, in practical situations, direct measurements of slow-time variables are not available and so the parameters of equation (3.13) cannot be obtained. Thus, the feature vector Embedded Image has to be used directly to infer any observable changes in the slow-time variables.

For this purpose, the estimated feature vectors Embedded Image for all Embedded Image data records are stacked as row vectors into a tracking matrix Y, i.e.Embedded Image(3.14)Thus, each of the Embedded Image columns of Y represents the time-series on the slow time-scale of each component of the average tracking metric Embedded Image, in each region Embedded Image. The time between each row of Y, Embedded Image, is of the same order as Embedded Image, where we recall that N is the number of data points in each record and Embedded Image is the sampling time of the fast data. Thus, Embedded Image can be though of as the effective slow sampling time used for measurements of the slow process.

Using the tracking matrix, the drift identification problem reduces to the statistical analysis of the slow time-series in its columns. Classical multivariate data analysis tools such as proper orthogonal decomposition (POD; Mardia et al. 1979) are an obvious choice for the study of Y. However, POD may not always be suitable for extracting slow-time trajectories from the tracking matrix. The main reason for this is that the POD tries to maximize the variance (or ‘signal energy’) in data projections. However, considerable signal variance can be attributed to possible bifurcations in the fast-time behaviour, which will not smoothly track the slow variables.1 Therefore, the identified proper orthogonal coordinates (POCs) can be expected to have jumps due to the bifurcations.

In order to correct this problem, a new multivariate analysis procedure has been developed that is specifically tailored for extracting information about smooth, deterministic processes: the SOD. The SOD analysis can be viewed as a constrained version of the POD (Chelidze & Zhou 2006).

In the POD, proper orthogonal modes (POMs) identify principal directions of the variance ellipsoid in the multi-dimensional data space. When the data is projected onto the POMs, we obtain POCs that maximize the signal's ‘energy content’ (i.e. the variance) in each mode. In contrast, with the SOD, the projected coordinates are required to maximize their total variance, but also to minimize the local variation (i.e. to maximize local smoothness). The first mention of this idea can be attributed to the optimal tracking method that was proposed for scalar damage process identification by (Chatterjee et al. 2002). This method relies on the existence of underlying deterministic behaviour of the damage accumulation process, but does not require its model. It assumes only that the damage accumulation process follows some smooth trend.

Given the matrix Y, we are looking for its linear projection Embedded Image onto a slow time-series Embedded Image that evolves smoothly and has maximum total variance. In other words, we seek the solution to the following constrained maximization problem:Embedded Image(3.15)where D is a slow-time differential operator that can be approximated by a Embedded Image matrix:Embedded Image(3.16)Using Euler–Lagrange equations, one can easily show that equation (3.15) is equivalent to the following generalized eigenvalue problem:Embedded Image(3.17)The eigenvector, or smooth orthogonal mode (SOM), q corresponding to the largest eigenvalue λ of equation (3.17) yields the projection of the matrix Y that maximizes the smoothness and overall variance of the smooth orthogonal coordinate (SOC) time-series Embedded Image.

For practical purposes, the solution to equation (3.17) is obtained using the generalized singular value decomposition (GSVD) of the matrices Y and Embedded Image (Golub & Van Loan 1996). In particular, these matrices are decomposed as:Embedded Image(3.18)where I is the identity matrix, the matrices U and V are unitary, X is a square matrix, and C and S are non-negative diagonal matrices. Then the SOMs (Embedded Image for Embedded Image) are the column vectors of Embedded Image, the smooth orthogonal values (SOVs) are given by Embedded Image, and the SOC time-series Embedded Image are given by the columns of UC. Some important properties of the SOD are described in Chelidze & Zhou (2006). In particular, it is shown that the procedure is invariant under linear transformations of data, while the POD is not.

In the context of tracking applications, the end result of the SOD analysis can be interpreted as follows. The transformed tracking matrix Embedded Image, with columns Embedded Image, can be sorted column-wise according to Embedded Image. Thus, the jth row of Embedded Image represents the original multi-dimensional feature vector for record j, Embedded Image, expressed in terms of the basis Embedded Image, i.e.Embedded Image(3.19)We can estimate the dimension of the underlying drift process by identifying those Embedded Image corresponding to the dominant Embedded Image, say Embedded Image. Then, using only the first k columns of Embedded Image, equation (3.19) gives a k-dimensional tracking metricEmbedded Image(3.20)that is expected to be within a 1–1, nearly affine transformation of the true slow variable.

(e) Failure prognosis

Application of the ideas outlined in the previous sections results in a damage tracking time-series over the slow time-scale. Following the discussion at the end of §2, one can consider candidate damage models of the formEmbedded Image(3.21)where Embedded Image is related to the long-time average of Embedded Image, the slow vector field, via equation (2.7). We emphasize that we are not assuming that averaging can be applied analytically; rather, we use the concept of averaging only to justify the autonomous form of equation (3.21). A suitable model structure must be found from first principles, or empirically from prior application of the tracking algorithm.

Given a damage model equation (3.21), the TTF can, in principle, be estimated using the tracking metric time-series, although in practice this may be difficult. In the easiest case, when a scalar-tracking metric is adequate, i.e. when there is only a single dominant SOV, the tracking metric of equation (3.20) becomes Embedded Image, and so equation (3.21) takes the formEmbedded Image(3.22)

This scalar case will be adequate for the results presented in the next section.

Then, given a failure limit Embedded Image, the equation (3.22) can be integrated to give the failure time asEmbedded Image(3.23)Nonlinear recursive estimation, such as the extended Kalman filter or unscented filtering (Julier & Uhlmann 1997), can be used to get robust estimates of Embedded Image (Chelidze & Cusumano 2004). The different equations used to define the estimator treat the tracking metric time-series as a series of observations from which one wishes to estimate the actual damage state using the ‘sensor model’ of equation (2.5). Thus, a side benefit of using the recursive filter with a specific damage model is that one is able to obtain an estimate of the actual damage state consistent with the model.

4. Experimental application

For the experiments described here, we used a two-well magneto-elastic oscillator, modified as described in (Chelidze et al. 2002). In the system (refer to figure 1), a clamped-free beam is restricted to a single degree of freedom by steel bar stiffeners (each 192.1×5.2×12.8 mm). The beam itself is made of 209.6×0.7×12.7 mm high carbon steel bar. Top 25.5 mm of the beam is left clear and acts as hinge. Two rare-earth magnets near its free end provide a two-well potential. They are placed 30 mm apart with a clearance of 3 mm between the tip of the beam and the magnets. The natural frequencies of small oscillation in the front and back energy wells are 7.85±0.01 and 8.58±0.01 Hz, respectively. The beam displacement is measured by a strain gauge mounted close to the clamped end. A shallow notch is machined on one side of the beam across its width, below the strain gauge and just above the stiffeners. The system is mounted on a shaker and is forced at 6 Hz. Care was taken in assembly to ensure that only bending vibrations were excited. The damage accumulates slowly and the experiment is run until complete fracture of the beam. Strain gauge output is low-pass filtered with 50 Hz cut-off frequency, sampled at 100 Hz sampling frequency, digitized (16 bit A/D), and stored on a computer.

Figure 1

Configuration of the two-well oscillator experimental system.

Several experiments were conducted. As expected, a range of failure times is observed in the experiments: in some instances, the beam broke rapidly, after as little as half an hour and in others, experiments were run for several days without failure. However, a detailed examination of the statistical properties of the observed failure times is beyond the scope of the present paper, and is reserved for future work. Here, we aim only to demonstrate the capability of the approach to determine fundamental properties of damage evolution dynamics in specific cases. Thus, after some trial and error, two experiments were run in which the beams failed in approximately 2.5 h. These provided sufficient and manageable amounts of data to apply our technique, for trials with approximately similar outcomes, as required for the purposes of this study. The data from these two trials were used for further analysis. In both experiments, the estimated embedding dimension for the reference dataset was found to be 5, and the delay time was estimated to be 5ts. The first 214 data points were used for the reference data set, and Embedded Image nearest neighbours were used for the local linear model parameter estimation. We have considered data records of 104 point size. The results obtained for both experiments were very similar, confirming the reproducibility of the calculations. Therefore, in the interest of brevity, results of only one experiment will be described in detail in what follows.

After going through the embedding and modelling process, we split the data into Embedded Image overlapping records of Embedded Image size. The overlap between the consecutive records is 90% of data. The estimated tracking function Embedded Image was calculated using the single step (Embedded Image) STRMP error of the reference model for each point y in the data record. We have also recorded the corresponding distances Embedded Image, which are defined for equation (3.9).

The reference phase space was partitioned into 27 equi-probable hypercuboids. In particular, first the reference phase space was cut along the second dimension into three parts so that each part contained approximately the same number of points. Then, this process was repeated to obtain hypercuboids along the third and fourth dimensions. We used the central coordinates of the reconstructed phase space for partitioning to avoid the edge effects.

The average of the estimated tracking function was calculated in each of the partitions for each data record and combined into the multi-dimensional feature vector, equation (3.12). These vectors were assembled into the 924×135 tracking matrix Y, equation (3.14). Each column of Y was mean subtracted in preparation for multivariate analysis.

The SOD analysis showed that the data were dominated by one scalar coordinate, as shown by the results in figure 2. In the semi-log plot of figure 2a, it is evident that the first generalized eigenvalue is about 1.5 orders of magnitude larger then the rest. Plotting the corresponding SOC, Embedded Image, as in figure 2b shows smooth power law type behaviour, even though the actual load history at the notch is quite complex, consisting of many chaotic/periodic transitions. By comparison, the SOCs corresponding to the next four significant SOVs (Embedded Image through Embedded Image) are shown in figure 3. All of them show an increase in amplitude at the end of the experiment, after about 1.5 h where the process begins to run away as indicated by the plot of figure 2b. However, in general the signals for the higher SOCs have significantly larger fluctuations and do not show the monotonic trend observed in the first SOC. Hence, we have a strong indication that the damage evolution process is one-dimensional, and we take Embedded Image to be our tracking metric.

Figure 2

SOD tracking analysis results: (a) the first 20 largest generalized eigenvalues (SOVs); (b) the SOC time-series (Graphic) corresponding to the largest SOV.

Figure 3

Plots of the second through the fifth most dominant SOCs (i.e. Graphic through Graphic) versus time.

We remark that for these particular experiments, a POD analysis yields approximately the same results as the SOD analysis. We speculate that the similarity in this case could be attributed to the dominance of only one slow coordinate in the data. Another contributing factor for this similarity could be that the total change in the fast-time system parameters is relatively small, at least prior to total failure. This is due to the fact that the major part of the restoring force of the beam is provided by the magnetic potential, which is constant, whereas fatigue damage and crack growth at the notch effects the elastic potential, which constitutes a relatively small part of the total restoring force. The similarity between the SOD and POD results in these experiments is contrary to what was observed in Chelidze (2004), in which the magnetic potential itself was altered and significant and repeated bifurcations were observed in the fast-time response. In this electromagnetic ‘damage’ case, the SOD was shown to provide a more consistent identification of drifting variables since it is robust with respect to the bifurcations. Nevertheless, even though a POD analysis could be used for the analysis of this paper, it was observed that the SOD was still superior in that the separation between the first two POVs were not as pronounced as for the SOVs, and thus the SOD gives a sharper indication of the slow process dimensionality.

In these experiments, we had no independent measure of accumulating fatigue damage and crack growth at the notch. However, in another set of experiments (Chelidze & Liu 2005), digital images of the beam cross-section near the notch were used to obtain an independent scalar measure of damage, which was shown to be in approximately linear 1–1 relationship with the identified tracking metric. This provides further support for taking the first SOC as our tracking metric.

Examination of figure 2b suggests that the tracking metric satisfies a Paris-type power law familiar from the study of macroscopic crack growth,Embedded Image(4.1)which after integrating as in equation (4.1) givesEmbedded Image(4.2)where β is some constant of integration. In fact, for carbon steel, which is used in our experiments, we expect α=2. Therefore, the reciprocal of the tracking metric should be a linear function of time, which is indeed the case as shown by the blue line in the left plot of figure 4, starting a little over half way through the experiment. The linear fit to this reciprocal over the time interval of 1.4–2.5 h (shown with a red line) yields the appropriate value for ϵ=60.5 h−1.

Figure 4

Damage model estimation and time to failure (TTF) prediction. (a) The reciprocal of the dominant SOC versus time given by the black jagged line, indicating that the variable evolves according to a power-2 Paris-type power law. The linear fit to this data over the time interval from 1.4 to 2.5 h is shown by the straight grey line. (b) The estimated TTF versus time, given by the black line, with the actual TTF (known a posteriori) overlayed in grey.

So far, the PSW tracking results have been used to determine the dimensionality of the slow damage process, identify the appropriate damage model, and determine its parameters. If we assume that this model and its parameters are known a priori, we can use it for failure prognosis. In our case, we can use equation (4.2) with α=2 to obtain a measurement of the TTF for some predefined failure value Embedded Image:Embedded Image(4.3)

Then the Kalman process and the measurement equations will simply be:Embedded Image(4.4)where t is the remaining TTF, Embedded Image is the current measurement of the TTF from equation (4.3), Embedded Image is the time interval between the measurements of Embedded Image, w is the process noise and v is the measurement noise. The noise terms w and v are assumed to be white, independent, mean zero Gaussian random variables with standard deviations Q and R, respectively. The measurement noise standard deviation R is estimated from equation (4.3) asEmbedded Image(4.5)where Embedded Image is taken to be the corresponding average of the standard deviation of Embedded Image, which is defined for equation (3.9).

For the TTF estimation, we applied a recursive Kalman filter to equations (4.3) and (4.4), with the final (failure) value of Embedded Image taken to be Embedded Image. Six hours was given as an initial guess for the TTF and the process noise standard deviation was set at Q=1. The results of this calculation are shown in the right plot of figure 4. Here, the estimated TTF is shown by the blue line and the actual TTF (known a posteriori) shown as a red line.

Figure 4 demonstrates that the TTF estimate converges to the true value more than 1 h before the total failure of the beam. The time required for convergence may be related to the difficulty of carrying out the required estimation when the tracking metric time-series is almost flat (refer to figure 2). However, we also hypothesize that it indicates the difference between the crack nucleation phase, for which the simple power law model does not work, and the actual crack propagation phase, for which it does. The fact that the α=2 Paris' law model works at all is remarkable given that the crack loading is decidedly not periodic.

5. Summary and conclusions

In this paper, we have presented a new type of nonlinear time-series analysis for slowly drifting dynamical systems based on the concept of PSW. The small distortions, or ‘warping’, in a fast subsystem's phase space caused by the slowly drifting variable are tracked using the STRMP error of a nonlinear reference model. The reference model is estimated using data in the reconstructed phase space of the fast-time subsystem while it is in its reference, or undamaged state. We have discussed the conditions under which the STRMP error will be a tracking function, equal to within an affine map to the hidden slow variable. Using this approach, one can generate a tracking metric that allows one to continuously track changes in a slowly evolving hidden variable. We have shown how a vector-tracking metric can be obtained from the raw STRMP error vector data using the SOD analysis, which provides a measure of how the STRMP error varies across different regions of the fast subsystem's phase space.

The approach was applied to the study of an experimental system consisting of a vibrating beam nonlinear oscillator with a crack that grows to complete failure. The vector PSW tracking analysis provided evidence that the damage evolution is governed by a scalar process, and by looking at the slow time-series of the resulting scalar-tracking metric, we demonstrated that the governing damage evolution equation has the form of a power-2 Paris-type damage growth law. Finally, the identified damage evolution model was used to predict the TTF using an additional recursive estimation step. It was shown that the TTF estimates converged to the true values well in advance of the actual complete fracture of the beam.

The PSW method has been demonstrated on a system involving material damage evolution; however, the formulation of the approach (which is really a family of techniques more than a single method) is quite general, and we expect it to be applicable to the study of a wide variety of both natural and human-made slowly drifting systems.


D. C. would like to acknowledge the support of this work by the NSF grant no. 0237792. The work of J. C. is supported in part by the Air Force Research Laboratory, via Miltec contract no. 04-C-0015.


  • One contribution of 15 to a Theme Issue ‘Exploiting chaotic properties of dynamical systems for their control’.

  • Even though the STRMP error at a point is insensitive to bifurcations, the ensemble of points available for computing the tracking metric can vary significantly in systems that undergo stability transitions. This is the source of the possible additional signal variance referred to here.


    View Abstract