Royal Society Publishing

An empirical approach for delayed oscillator stability and parametric identification

Brian P Mann , Keith A Young


This paper investigates a semi-empirical approach for determining the stability of systems that can be modelled by ordinary differential equations with a time delay. This type of model is relevant to biological oscillators, machining processes, feedback control systems and models for wave propagation and reflection, where the motion of the waves themselves is considered to be outside the system model. A primary aim is to investigate the extension of empirical Floquet theory to experimental or numerical data obtained from time-delayed oscillators. More specifically, the reconstructed time series from a numerical example and an experimental milling system are examined to obtain a finite number of characteristic multipliers from the reduced order dynamics. A secondary goal of this work is to demonstrate a benefit of empirical characteristic multiplier estimation by performing system identification on a delayed oscillator. The principal results from this study are the accurate estimation of delayed oscillator characteristic multipliers and the utilization the empirical results for parametric identification of model parameters. Combining these results with previous research on an experimental milling system provides a particularly relevant result—the first approach for identifying all model parameters for stability prediction directly from the cutting process vibration history.


1. Introduction

Systems governed by delay differential equations are relevant to many fields of science and engineering. For instance, recent literature has examined applications in economics, robotics, biology and manufacturing processes (see Baker et al. 1998; Mann et al. 2004; Wahi & Chatterjee 2004). The qualitative investigation of these dynamical systems often includes a stability analysis—where compact representations of the stable and unstable parameter domains are often presented in the form of stability charts.

A well-known outcome of introducing a time delay into a dynamical system is the growth of the phase space from a finite dimension to an infinite dimension (see Stépán 1989; Hale & Verduyn 1993). However, recent approaches for examining the stability behaviour of time-delayed systems have shown that it is possible to ascertain the stability characteristics of the resulting infinite dimensional system from a finite dimensional phase space (see Stépán 2001; Davies et al. 2002). A common goal in these analyses is to transform the original set of equations into the form of a dynamic map that maps the system behaviour over a single time delay (Insperger et al. 2003a; Mann et al. 2003). However, it is potentially useful to determine the stability behaviour solely from experimental or numerical data.

The primary goal of this paper is to investigate an empirical method for characteristic multiplier estimation applicable to systems with time delays. The presented methodology differs from previous works on systems without delays (see Bayly & Virgin 1993) since the dimension of the phase space or truncated dimension in this case, is not known a priori. Since additional considerations are required in the examination of experimental data, such as the presence of measurement noise, a numerical study of the delayed damped Mathieu equation is undertaken prior to examining an experimental system.

Empirical characteristic multiplier estimates are obtained via two different methods: (i) a perturbation is applied to the steady-state motion and the stability of the underlying behaviour is inferred from the resulting transient and (ii) the coupled tool–workpiece transient interactions at cutter entry are analysed for an experimental milling system. In each instance, empirical estimates are compared with the characteristic multiplier estimates of a semi-analytical approach based upon temporal finite element analysis (Bayly et al. 2003; Mann et al. 2004).

A secondary goal of this paper is to demonstrate the potential benefit of characteristic multiplier estimation in the common problem of system identification. This idea is demonstrated on a non-smooth delayed oscillator, an experimental milling system, to provide a particularly interesting result—the first method capable of identifying the necessary model parameters for stability prediction directly from the cutting process vibration history.

The work presented in this paper is organized as follows. Section 2 introduces a system for numerical study and describes an empirical approach for characteristic multiplier estimation. The accuracy of the largest characteristic multiplier is then examined by comparisons with a semi-analytic approach. The final stage of this work focuses on an experimental study of a delayed oscillator system and describes the steps taken for parametric identification of model parameters.

2. Empirical characteristic multiplier estimation

The stability of systems governed by time-periodic differential equations is important to various fields of science and engineering. For instance, recent literature has described applications in high-speed milling (Altintas & Budak 1995; Tlusty 2000; Balachandran 2001), quantum mechanics and rotating helicopter blades (Sinha et al. 1996). Since it is commonly difficult to provide real-time control, it is imperative to consider the influence of a time delay on the system stability. The delayed damped Mathieu equation, which has been the focus of several recent works (Insperger & Stépán 2002), is studied here since it is a representative system with both delayed feedback and parametric excitation. Furthermore, this system is considered prior to studying an experimental system to provide noise-free numerical data.

The expression for the delayed damped Mathieu equation isEmbedded Image(2.1)where τ represents the time delay in the feedback of the position variable. Section 2a describes an empirical approach for investigating the stability behaviour of the above system. Empirical characteristic multiplier estimates are then compared to the characteristic multiplier estimates of a semi-analytical technique in figure 1. The semi-analytical technique transforms equation (2.1) into a discrete dynamical system; this results in a mapping solution over a single time delay and the system stability can be determined directly from the dynamic map characteristic multipliers.

Figure 1

Numerically generated results for a single time series: (a) the time series with markers indicating the location of state selection and (b) the reconstructed Poincaré map from numerical time series. (c) The largest empirical characteristic multiplier (open circle) and semi-analytical characteristic multiplier (dashed line) that have been estimated as a function of δ. All results are presented for Embedded Image, ϵ=1 and Embedded Image, and (a) and (b) additionally use a fixed value of Embedded Image.

(a) The empirical technique

This section describes an empirical approach that can be applied to experimental or numerical data to infer stability. Characteristic multiplier estimates are obtained by applying a single perturbation during numerical simulation. The empirical approach begins by assuming a general nonlinear system that can be written as a discrete time system or dynamic map,Embedded Image(2.2)where Embedded Image and Embedded Image are Embedded Image state vectors (Guckenheimer & Holmes 1983). Expanding the nonlinear terms on the right side about a fixed point solution (Embedded Image) and retaining only linear terms results inEmbedded Image(2.3)where Embedded Image is a Embedded Image Jacobian matrix. Collecting periodic samples of the perturbed motion, at the frequency of the time delay, results in a Poincaré section. Successive periodic samples provide a Poincaré map which can be written as m Poincaré mappingsEmbedded Image(2.4)where m is the number of periods a sample is collected. The matrices Embedded Image and Embedded Image are of dimension Embedded Image and the elements of these matrices are given byEmbedded Image(2.5)The matrix Embedded Image can be found directly from equation (2.4) once a sufficient number of Poincaré points have been recorded (i.e. m=d satisfies the uniqueness criteria). However, better estimates are obtained for the over constrained case (m>d) and the expression for Embedded Image can be found in a least-squares senseEmbedded Image(2.6)The explanation for this improved result is that the least-squares approach provides a reduction in either experimental noise or numerical imprecision (see Virgin 2000).

At this point, it is important to note that the correct choice for the truncated system dimension has not been discussed. Once the choice for the state vector dimension has been obtained, the characteristic multipliers from empirical approach can be compared to those obtained via the semi-analytical method. More specifically, the eigenvalues of Embedded Image can be compared to the eigenvalues of Embedded Image in equation (3.6). Section 2a(i) discusses the methods employed to obtain the correct reduced order model.

(i) Reduced order dynamics

Two questions about the empirical approach have remained unanswered at this point. The first question is the proper choice of the state space vector dimension. The second question is the correct choice for the additional states of the system since only a single position state is likely to be available in an experiment and one additional velocity state would be available from numerical data. To help and facilitate the discussion that follows, an example of time series from the delayed damped Mathieu equation is shown in figure 1.

The initial quest to answer the above state selection and reduced order dimension questions was investigated by the application of delayed embedding techniques and standard methods for dimension estimation from nonlinear time-series analysis (Nayfeh & Balachandran 1995). For instance, a pseudo phase space—which is considered topologically equivalent to the original phase space—can be reconstructed from a single measurement. This can be accomplished by finding a suitable time lag in the data, such as the first minimum of a mutual information diagram and reconstructing the pseudo phase space via delayed reconstruction techniques (Abarbanel 1996). The time-series dimension estimation of the data shown in figure 1 was performed with a false nearest neighbour routine (Abarbanel 1996); this resulted in an attractor dimension of one for a time series of length Embedded Image. Since this approach determines the attractor dimension, as opposed to the state space dimension for a delayed oscillator (Kantz & Schreiber 2004), the correct choice for the state space dimension was not clear. When only the short duration transient motions are examined with these methods, difficulties occur in estimating the proper time lag and a suitable choice for the system dimension due to the short-time history.

A second argument for choosing the additional state vectors can be made from the semi-analytical technique presented in §3. The use of the interpolated Hermite polynomials makes the coefficients of the assumed solution represent the velocity and displacement at evenly spaced points along the interval of the time delay. The final question of obtaining a proper lower dimensional representation of the higher dimensional data was investigated via signal analysis techniques. The two methods applied to examine the experimental data include: (i) the rank of the data matrices and (ii) singular value decomposition which has also been referred to as proper orthogonal decomposition, principal component analysis or the Karhunen–Loeve decomposition (see Chatterjee 2000).

The first approach for choosing the system dimension included embedding the d-dimensional state vector with equally spaced displacement data over the time interval τ. Poincaré points, shown by the three different markers for the time series in figure 1, were collected over m periods, where m was taken to be much larger than d and the rank of the Embedded Image matrix was computed while increasing the dimension of the state vector. The proper choice for the system dimension was set to the dimension at which the rank of the matrix stopped growing with an increase in the assumed system dimension. For the time series of interest, the rank of the data matrix was found to be three. Figure 1 shows the reconstructed Poincaré map using this approach and a comparison of the largest characteristic multiplier estimates from the empirical approach and the semi-analytical approach as a function of δ.

Since the empirical data in the afore mentioned example was void of experimental noise, it is imperative to note that erroneous estimates of the correct dimension were found for the introduction of experimental noise. Therefore, a second method was used to obtain the system dimension; this approach embedded the d-dimensional state vector with equally spaced displacement data over the interval of the time delay. Poincaré points were collected over m periods, where m was taken to be much larger than d and singular value decomposition of the Embedded Image matrix was computed while increasing the dimension of the state vector. The proper choice for the truncated system dimension was taken to be at the point, where a levelling off in the magnitude of the singular values occurred. Although these criteria produced identical results for the noise-free empirical data examined here, it is noted that very different results were obtained for the experimental data shown in the latter half of this paper. Therefore, the singular value approach was found to be necessary when choosing the system dimension in the presence of experimental noise.

3. Semi-analytic characteristic multiplier estimation

This section presents a semi-analytical approach for investigating the stability of the delayed damped Mathieu equation. The system is examined over a single period, where the period for the system corresponds to the time delay (τ), by dividing the time period into a finite number of temporal elements. For instance, the motion of delay oscillator is approximated as a linear combination of polynomials within each time interval or time element usingEmbedded Image(3.1)Here, the trial functions or polynomials (Embedded Image), which are given by equations (A 1a)–(A 1d) in the appendix, are written as a function of the local time (Embedded Image). The local time is defined to vary from zero to the maximum time for each element (Embedded Image) and the superscript for the polynomial coefficients denotes the nth period. Substitution of the assumed solution (equation (3.1)) into the governing equation for the system leads to a non-zero error. The error from the assumed solution is weighted by multiplying a set of test functions and setting the result of the integrated weighted error to zero. This provides two equations per element (Peters & Idzapanah 1988), which are given byEmbedded Image(3.2)where the integration time for each element is Embedded Image and Embedded Image is a different constant for each element (i.e. the elapsed time from the beginning of first element to the start of the current element for the current delay period). The introduction of Embedded Image is required to maintain continuity of parametric excitation term at the beginning and end of each element. It is noted that the equations for each element are linear in the coefficients of the assumed solution. Combining the resulting equations for each element, a global matrix equation can be obtained that relates the states of the system in the current period to the states of the system in the previous periodEmbedded Image(3.3)where the above matrices are written for two elements and the terms inside the above matrices are given by Embedded Image, Embedded Image, Embedded Image, Embedded Image andEmbedded Image(3.4)Embedded Image(3.5)Equation (3.3) can be written in more compact form or in the form of a discrete system dynamic mapEmbedded Image(3.6)The stability of the system can then be determined directly from the monodromy operator (Embedded Image), which provides a mapping over a single time delay. If any monodromy operator eigenvalue, which are also referred to as the characteristic multipliers for the dynamic map, has magnitude greater than one, the system is considered unstable. Figure 2 is an example of the stability chart generated using the semi-analytical method along with characteristic multiplier trajectories in the complex plane. The characteristic multiplier trajectories are examples of a flip bifurcation, when negative real characteristic multipliers penetrate the unit circle and a saddle node bifurcation that penetrates the unit circle with values that are positive and real.

Figure 2

Stability chart for the delayed damped Mathieu equation, shown in graph (a), for δ versus ϵ parameter space; graphs (a)–(c) were generated for Embedded Image and Embedded Image. The two graphs at the bottom provide characteristic multiplier trajectories in the complex plane for a flip bifurcation (b) and a saddle node bifurcation (c).

4. Parametric identification from milling data

Machining operations are a common example, where a time delay can result in unstable behaviour. For instance, a primary limiting factor in material removal applications is the relative oscillations between a cutting tool and workpiece system (Davies et al. 2002; Insperger et al. 2003b; Mann et al. 2005a). More specifically, cutting force models become non-ideal energy sources that must capture force modulations created from relative motions between the cutting tool and workpiece system (Mann et al. 2005b; Olgac & Sipahi 2005). The most common approach is to prescribe a cutting force model that is a function of the uncut chip area. This type of model couples the relative system of oscillations to the cutting forces and chip load variations; these motions can cause dynamic cutting forces and self-excitation of the machine-tool structural modes that lead to instability. Unless avoided, these unstable regions may cause large dynamic loads on the machine spindle and table structure, damage to the cutting tool and a poor surface finish. Therefore, a key strategy is to apply analytical and numerical methods to predict parameter domains where instabilities exist.

This section investigates empirical characteristic multiplier estimation, parameter identification and stability prediction for an experimental system that can presumably be modelled as a time delay system. The previously presented ideas are implemented for an experiment that is designed to be a single degree of freedom. Figure 3 shows an example for cutting test, where two regions of transient motion are found. Unlike the cutter entry transient, which consists of coupled tool–workpiece motion, the exit transient consist of partial tool–workpiece engagement and post-cutting motions of the system. This makes it difficult to utilize the post-cutting oscillations in parameter estimation. In the study that follows, the coupled cutter-workpiece transients at cutter entry are used to obtain Embedded Image and Embedded Image, the steady-state region of the time series is used to extract the fixed points Embedded Image, and the post-cutting transients are used to obtain a first approximation for the system damping ratio (ζ) and natural frequency (ω); the estimated post-cutting parameters are then utilized as initial conditions for a nonlinear optimization routine.

Figure 3

Experimental measurement of flexure oscillations during a cutting test performed at the following parameters b=2.9 mm, h=0.1016 mm rev−1 and Ω=3810 r.p.m. (see §4a for parameter definitions).

(a) Milling experiment description

Experiments were performed on a flexure system designed to be an order of magnitude more compliant than the cutting tool (see figure 4). The monolithic, uni-directional flexure was machined from aluminium and instrumented with a single non-contact, eddy current displacement transducer. In comparison to the compliant direction of the flexure, the values of stiffness in the perpendicular directions were more than 20 times greater, as was the stiffness of the tool. The displacement transducer output was anti-alias filtered and sampled (16-bit precision, 12 800 samples s−1) with data acquisition hardware connected to a laptop computer. A periodic 1/tooth pulse was obtained with the use of a laser tachometer to sense a black–white transition on the rotating tool holder.

Figure 4

Schematic of the single degree of freedom milling system that provided an experimental system with time delays due to self-excited oscillations.

The compliant direction of flexure was aligned with the tool feed direction during the cutting tests. Aluminium (7075-T6) test specimens, of width 6.35 mm and length 100 mm, were mounted on the flexure and milled with a single flute 19.050 mm diameter end mill. The spindle speed (Embedded Image) and depth of cut (b) were changed for each cutting tests to determine the onset of unstable vibrations. Specimens were up-milled, meaning the tool was positioned behind the workpiece and progressed from right to left in figure 4, at a constant feed of 0.1016 mm rev−1 and a fraction of the spindle period in the cut of Embedded Image.

(b) Milling process model and semi-analytic dynamic map formulation

Parametric identification of the system model parameters requires a method to estimate the characteristic multipliers from the governing equations of motion. Although the analysis presented here is reminiscent to the approach previously described in Mann et al. (2004), the details are included here for completeness. The dynamics of the system under consideration can be represented as a single degree of freedom modelEmbedded Image(4.1)where ζ is the damping ratio, ω is the circular natural frequency and b is the axial depth of cut. The terms of Embedded Image and Embedded Image are compact notation forEmbedded Image(4.2)Embedded Image(4.3)where h represents the feed per tooth and Embedded Image is a Heaviside step function (i.e. Embedded Image for Embedded Image and is otherwise zero). The terms Embedded Image and Embedded Image represent cutting coefficients that have been divided by the flexure first modal mass (i.e. the final units are pressure per unit mass). A variational system can be formed by writing the solution as a perturbation about the τ-periodic motionEmbedded Image(4.4)substituting equation (4.4) into equation (4.1) and subtracting off the equation for only periodic motion givesEmbedded Image(4.5)With the introduction of a common approximation, Embedded Image (see Tlusty 2000) the expression for final undefined term can be written asEmbedded Image(4.6)Following the same discretization method as shown for the delayed damped Mathieu equation of §3, a dynamic map can be obtained from equation (4.5) in the same form as equations (3.6) and (3.3). However, the terms inside the matrices of equation (3.3) are given by equations (4.7)–(4.12)Embedded Image(4.7)Embedded Image(4.8)Embedded Image(4.9)Embedded Image(4.10)Embedded Image(4.11)Embedded Image(4.12)

(c) Experimental system parameter identification

This section describes a novel approach taken to identify the unknown milling model parameters directly from the piecewise continuous delayed oscillator time series. This approach utilizes the empirically estimated characteristic multipliers in the identification process. Parameters identified from the empirical approach are then compared with parameters identified from separate modal and dynamometer tests. Next, a stability chart is constructed for both parameter sets to illustrate the validity of the results.

An example of data series, embedded in a four dimension phase space, is shown in figure 5 along with the normalized singular values. The rapid decay in the trend of the singular values by the fourth value was used as an indicator for the dimension selection. Parametric identification of the unknown milling model parameters was performed by minimizing the error between the largest semi-analytical characteristic multiplier and the largest empirical characteristic multiplier for three different cutting tests. The testing conditions and empirical characteristic multiplier values from each test are listed in table 1.

Figure 5

Experimental cutting test data performed at 2850 r.p.m. and b=2.9 mm. The three views are: (a) stroboscopic samples of the chosen four states; (b) singular values of the experimental data matrix Embedded Image and (c) a mapping of data to a mod τ torus and sectioning to provide a Poincaré map of a stable fixed point.

View this table:
Table 1

Empirical characteristic multipliers and testing conditions.

A nonlinear least-square optimization routine was used for the parameter identification process with the following values taken as initial conditions for the routine Embedded Image, ω=930 rad s−1 and Γ1=5.0×108 N m−2 kg−1. Since the optimization routine was a gradient-based search algorithm, perturbations to these initial conditions were also investigated. The convergent result of the optimization routine provided estimated parameters of Embedded Image, ω=921.1 rad s−1 and Γ1=2.16×108 N m−2 kg−1. Upon comparing these values to the parameters obtained from separate modal and dynamometer cutting tests, see table 2, the difference in the estimated damping ratio was the only concern. However, a comparison stability charts for each parameter combination, shown in figure 6, reveals only subtle differences between the two stability borders.

View this table:
Table 2

Comparison of estimated parameters from the empirical method to parameters obtained from separate dynamometer cutting tests and experimental modal tests.

Figure 6

Predicted stability boundaries from empirical approach (solid line) and separate modal/dynamometer tests (dotted line) versus experimental results. Symbols denote: (i) open circle is a clearly stable case; (ii) open inverted triangle is an unstable cutting test and (iii) plus is a borderline unstable case (i.e. not clearly stable or unstable).

(d) Milling experiment results

Experimental stability results have been superimposed onto stability predictions made using the estimated parameters from a single cutting tests (see figure 6). The results from the experimental cutting tests show strong agreement with the predicted results. Additionally, the stability predictions are nearly identical to those obtained when separate modal test and dynamometer cutting force measurements are used for parameter estimation.

Raw displacement measurements were periodically sampled at the tooth passage frequency to create 1/tooth passage displacement samples and Poincaré sections shown in displacement versus delayed displacement coordinates; these plots are also shown with the power spectral density (PSD) of the continuously sampled displacement in figure 7. The time-series reconstructions shown in this section were performed using standard embedding procedures outlined in Abarbanel (1996). Tests were declared stable if the 1/tooth passage sampled displacement or τ-periodic samples, approached a fixed point value (see figure 7, case D).

Figure 7

Experimental up-milling measurement data for cases (a, b, c, d) shown in figure 6. Each row contains a x-axis 1/tooth displacement plot, a Poincaré section shown in displacement (Embedded Image) versus delayed displacement (Embedded Image) coordinates and the tooth passing frequency is marked by open circle on the power spectrum (PSD). Unstable cutting processes are shown in cases: (a) Ω=3050 r.p.m., b=2.7 mm, (b) Ω=3565 r.p.m., b=2.1 mm, (c) Ω=3635 r.p.m., b=2.0 mm, and (d) Ω=3810 r.p.m., b=2.9 mm is a stable process. Cases (a) and (b) result from unstable Neimark–Sacker or secondary Hopf bifurcations and case (c) shows a flip bifurcation.

Unstable behaviour predicted by complex characteristic multipliers with a magnitude greater than one, shown by cases A and B, corresponds to Neimark–Sacker or secondary Hopf bifurcation. These post-bifurcation test results show the 1/tooth-passage displacement samples are incommensurate with the tooth passage frequency and quasiperiodic motions can be observed in the Poincaré sections.

Unstable period-doubling behaviour is predicted when the dominant characteristic multiplier of the dynamic map model is negative and real with a magnitude greater than one. Experimental evidence of this type of post-bifurcation behaviour is shown by case C of figure 7.

5. Conclusions

This paper investigates an empirical method for determining the truncated system characteristic multipliers for two systems that contain a single time delay. It is shown that the characteristic multipliers can be accurately obtained for both numerical and experimental data. The technique differs from previous research focused on systems without delays since the dimension of these systems is already known to be finite. In this paper, the authors assume the delay oscillator infinite dimensional phase space can be closely approximated by a finite dimensional system and methods for obtaining the truncated dimension are explored.

Empirical estimates for the system characteristic multipliers are obtained using either perturbations in the steady-state motion or the initial transients of the system, such as in the case of the experimental milling system. Characteristic multiplier estimates are compared with the results of a semi-analytical approach to perform parameter identification. The presented results are thought to be applicable to a wide variety of systems with time delays. In particular, the presented methods provide the first known methodology for parameter estimation and stability predictions directly from the milling process vibration history.


Support from US National Science Foundation CAREER Award (CMS-0348288) is gratefully acknowledged.


    • Received September 12, 2005.
    • Accepted January 13, 2006.


View Abstract