## Abstract

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*. 2003*a*; 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 is(2.1)where *τ* represents the time delay in the feedback of the position variable. Section 2*a* 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.

### (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,(2.2)where and are state vectors (Guckenheimer & Holmes 1983). Expanding the nonlinear terms on the right side about a fixed point solution () and retaining only linear terms results in(2.3)where is a 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é mappings(2.4)where *m* is the number of periods a sample is collected. The matrices and are of dimension and the elements of these matrices are given by(2.5)The matrix 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 can be found in a least-squares sense(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 can be compared to the eigenvalues of in equation (3.6). Section 2*a*(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 . 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 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 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 using(3.1)Here, the trial functions or polynomials (), which are given by equations (A 1*a*)–(A 1*d*) in the appendix, are written as a function of the local time (). The local time is defined to vary from zero to the maximum time for each element () and the superscript for the polynomial coefficients denotes the *n*th 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 by(3.2)where the integration time for each element is and 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 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 period(3.3)where the above matrices are written for two elements and the terms inside the above matrices are given by , , , and(3.4)(3.5)Equation (3.3) can be written in more compact form or in the form of a discrete system dynamic map(3.6)The stability of the system can then be determined directly from the monodromy operator (), 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.

## 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*. 2003*b*; Mann *et al*. 2005*a*). 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*. 2005*b*; 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 and , the steady-state region of the time series is used to extract the fixed points , 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.

### (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.

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 () 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 .

### (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 model(4.1)where *ζ* is the damping ratio, *ω* is the circular natural frequency and *b* is the axial depth of cut. The terms of and are compact notation for(4.2)(4.3)where *h* represents the feed per tooth and is a Heaviside step function (i.e. for and is otherwise zero). The terms and 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 motion(4.4)substituting equation (4.4) into equation (4.1) and subtracting off the equation for only periodic motion gives(4.5)With the introduction of a common approximation, (see Tlusty 2000) the expression for final undefined term can be written as(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)(4.7)(4.8)(4.9)(4.10)(4.11)(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.

A nonlinear least-square optimization routine was used for the parameter identification process with the following values taken as initial conditions for the routine , *ω*=930 rad s^{−1} and *Γ*_{1}=5.0×10^{8} 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 , *ω*=921.1 rad s^{−1} and *Γ*_{1}=2.16×10^{8} 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.

### (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).

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.

## Acknowledgments

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

## Footnotes

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

- © 2006 The Royal Society