Robust stability of milling operations based on pseudospectral approach

The prediction of chatter-free machining parameters is not reliable for industrial applications without guaranteed robustness against modeling uncertainties. Measurement inaccuracies, ﬁtting uncertainties, simpliﬁcationsandmodelingerrorstypicallyleadtoamismatchbetweenthemathematicalmodel and the real physical system. This paper presents the robust stability analysis of milling operations basedonapseudospectralapproachtotakeintoaccounttheeﬀectofboundedparametricuncertaintiesbothinthecuttingcoeﬃcientsandinthemodalparameters.Inordertomakethepredictionsmoreaccurate,theoperationalmodalanalysisofthespindleduringrotationwasconducted.Thenaturalfre-quenciesanddampingratiosofthedominantvibrationmodeswereidentiﬁedfromtheimpacttestsoftherotatingtoolatdiﬀerentspindlespeeds.Theuncertaintiesoftheﬁttedparameterswereincludedinthecomputationofthepseudospectralradiusofthemonodromyoperatorofthetime-periodicsystem.Thesolveristestedinacasestudyandexperimentalchattertestsarealsoincludedfordemonstrationpurposes.


Introduction
The continuously increasing demand for competitive, fast, high-quality and economical production forces the manufacturing industry to push the operational conditions of machines to their limit. However, after a certain point often quality and reliability decreases, while deterioration and costs increase. In material removal processes, such as milling, turning, grinding, or drilling, one of the most destructive phenomenon that limits productivity the most is machine tool vibration [1].
The first mathematical equations describing the evolution of unstable vibrations appeared in the work of Tobias [2] and Tlusty [3] in the 1950s and 1960s, although the importance of the phenomenon was already recognized by Taylor at the beginning of the 20th century [4]. Their observations laid down the fundamentals of the so-called surface regeneration effect, which became the most widely accepted explanation for machine tool vibration. The geometry of the chip is modified by the relative motion between the tool and workpiece, which induces variation in the cutting force. The vibrations in the past are embedded in the surface of the workpiece, which affects every subsequent cut by the modified geometry of the chip that depends on present and past states, too. Dependencies on past events in dynamical models are described by delay-differential equations (DDEs). When the stationary solution of the DDE loses stability, the amplitude of vibration increases until the physical contact between the workpiece and the tool is lost. This large-amplitude selfexcited vibration is called chatter.
A milling operation is an intermittent cutting process, where the chip geometry varies periodically at every toothpassing interval. The surface regeneration couples with parametric excitation and results in a DDE with periodically varying coefficients, also called time-periodic delay-differential equation (TPDDE). The stationary solution in a stable operation is periodic (at most with the spindle revolution period), which can undergo secondary Hopf or period doubling bifurcations [5] as the machining parameters change.
The harmful unstable vibrations must be avoided during the operation in order to reach the required quality and maximize the production. The highest performance of machining centers, however, can only be reached if the machine is safely operated at high speeds close to stability boundaries. These conditions require a reliable and accurate prediction of dynamical behavior. The domain of chatter-free machining parameters are visualized by the so-called stability lobe diagrams. Nevertheless, deterministic models considered during calculations often fail to accurately predict dynamical behavior and therefore experimental cutting tests do not always match the expectations. One of the reason is the uncertainties in the dynamical model that affect the prediction of chatter-free machining parameters.
The need for reliable predictions considering variations in machining and dynamical parameters was already recognized by the engineering community, and many attempts were made to construct such stability lobe diagrams that take the effect of modeling errors into account. These diagrams are referred to as robust stability lobe diagrams, where the domain of robust stable machining parameters can retain the stability if bounded uncertainties in the model arise. In the simplest case, Monte Carlo simulations and statistical evaluations were performed in [6,7], which require significant computational effort and therefore their practical applica-Robust stability of milling operations based on pseudospectral approach tion without major improvements is limited. The advantage of these methods, however, is the statistical data provided after a large amount of computations. To obtain similar results without huge computational efforts, research papers also provide local sensitivity methods, e.g. [8,9,10], that are based on local derivatives (sensitivity) and extrapolations. A global sensitivity method is presented in [11], which approximates a multi-dimensional integral of the probability density function (PDF) by separated one-dimensional integrals. The Robust Chatter Prediction Method (RCPM) presented in [12] is also an alternative method for global analysis, which is based on the discretization of the PDF.
Robust methods, where the worst-case perturbations are sought, were also used in the literature. For instance, the Edge Theorem combined with the Zero Exclusion Principle presented in [13] provides robust stability boundaries for time-invariant models, but the computational effort is often large. A different robust method called Multi-frequency Solution with Structured Singular Value Analysis is based on the direct perturbation of the measured frequency response functions (FRFs), where parameter perturbations are replaced by FRF perturbations, see [14]. In this latter case the computational effort is significantly reduced, but the prediction can be very conservative.
This paper presents the application of a pseudospectral approach that considers uncertainties in the dynamical parameters of milling operations and provides robust stable machining domains to exploit the highest productivity. The method is based on an iterative solver, that computes the worst-case uncertainty and therefore no conservativeness is introduced. In order to make the computations effective, the calculation of the characteristic multipliers (Floquet multipliers) of the time-periodic DDE is carried out by discretizing the eigenvalue problem of the monodromy operator into a generalized eigenvalue problem (GEP), using a spectral method. The iterations converge to a local maximizer of the pseudospectral radius in the space of allowable perturbations. However, global optimality is achieved by incorporating a restarting strategy, where several dominant Floquet multipliers of the original system are used to initialize the iteration. The related problem of computing the globally rightmost point of pseudospectra of time-invariant systems is addressed in [15].
The paper is structured as follows. In Sec. 2 the dynamical model of milling operations is presented, then the stability of time-periodic DDEs is investigated in Sec. 3 by introducing the new formalism. The pseudospectral approach to access robust stability is detailed in Sec. 4. A case study to test the computation of robust boundaries along with experimental chatter tests are presented in Sec. 5. The calculated stability boundaries are compared to a structured singular value analysis in Sec. 6. The conclusions are presented in Sec. 7.

Modeling of milling operations
A regular cutting tool with uniform helix angle is presented in Fig. 1(a)  where ( , ⋅) ∈ ℝ 2 is the forcing vector at the tool-tip in the ( , )-plane, is the modeled number of degrees of freedom, [⋅] = diag(⋅) ∈ ℝ × are diagonal matrices, n, is the -th natural angular frequency ( = 1, … , ), is the corresponding relative damping ratio, ( ) = [ ( ), ( )] ⊤ denotes the displacement of the tool-tip in the ( , )-plane and ∈ ℝ 2× is the modal transformation matrix that connects the modal space and the tool-tip's motion as ( ) = ( ) [16].
The cutting force acting at the tool-tip is described by empirical force models, that highly depend on the tool geometry and material properties. Next, for sake of brevity, regular cutting tools with number of equally distributed cutting edges and with uniform helix angle are modeled. Thus, the angular position of the -th cutting edge ( = 1, … , ) along the axial direction reads where is the coordinate along the axial immersion and Ω is the spindle speed given in rpm. The cross-section of the cutter is presented in Fig. 1(b). The elementary cuttingforce components in tangential and radial directions acting on tooth at a disk element of width d are given as where ℎ ( , ) is the theoretical chip thickness cut by tooth at axial immersion , moreover ℎ (ℎ) and ℎ (ℎ) are the empirical cutting force characteristics assumed in the shifted linear form where and are cutting coefficients in the main direction (tangential), and and are cutting coefficients in radial direction according to the Shearing & Ploughing (S&P) cutting force model [17]. The indicator function ( , ) is a piecewise constant function, which gives whether the cutting edge is in contact with the material or not. ( ) can be written as where en is the angle where the elementary disk starts cutting and ex is the angular position where it leaves the workpiece, see Fig. 1(b). The actual chip thickness cut by tooth at axial immersion then can be calculated approximately as where = [ , 0] ⊤ is the feed per tooth in direction , and the time-delay (equal to the tooth-passing period) in case of equally distributed cutting edges is = 60∕( Ω). The resultant cutting force vector in the ( , )-plane concentrated theoretically at the tool-tip is calculated as where the transformation matrix is written as Inserting the regenerative forcing model into the governing equation, then assuming small perturbation ( ) about the periodic motion p ( ) of the stationary cutting, i.e., ( ) = p ( ) + ( ), the equation of the variational system omitting the periodic term is written as where the directional matrix ( ) = ( + ) is introduced as where the dependency on ( , ) for = 1, … , make ( ) discontinuous, but piecewise smooth, see [18]. The time-periodic DDE (10) can be rewritten into a first-order form where is the × identity matrix and is the × zero matrix.

Stability of time-periodic time-delay systems
We observe that the governing equation (12) is of the forṁ where ( ) ∈ ℝ is the vector of state variables at time , and , are ℝ × -valued periodic functions with main period = . In the present and in the following section we briefly illustrate a novel pseudospectral approach to assess the stability robustness of a system of time-periodic DDEs, which was introduced in a more general form in [19]. In this manuscript we adapt the description to systems of the form (13). We first explain the adopted approach for general matrix-valued functions and . At the end of the section, we explain how the method can be improved by exploiting the property that the coefficients in (12) are piecewise smooth functions. The solution of the infinite-dimensional time-delay system governed by (13) depends on the time history of ( ), therefore an initial function is required to determine the evolution of states in time. In order to generalize the notation, let  = ([− , 0], ℂ ) be the space of continuous functions on the interval = [− , 0] and let us introduce the notation ∶= ( + ), ∈ [− , 0] in order to denote the solution segment over the interval [ − , ]. The Cauchy problem associated with the system (13) is written as with 0 ∈  being the initial function. The solution segment of the system at any time associated with the initial function 0 is described by the evolution operator  ( , 0 ) ∶  → , defined through the relation The asymptotic stability of the trivial solution ( ) ≡ of (13) is determined by the spectral properties of the monodromy operator  ∶=  ( , 0) (setting 0 = 0), as it is proven by the Floquet theory [20]. The nonzero elements of the spectrum of  are called characteristic multipliers (or Floquet multipliers) and can also be defined by with  being the identity operator. Any ∈ ℂ, ≠ 0 is a Floquet multiplier if and only if there exists an eigenfunction ∈ , ≠ , such that In general, there are no explicit expressions for the evolution operator and the characteristic multipliers. However, several numerical methods exist for approximating the operator and the multipliers.
As examples, interested readers are referred to numerical methods such as the Semidiscretization technique [5], the Spectral Element Method [21], [22], Chebyshev polynomialbased collocation methods [23,24] or the Improved Chebyshev Collocation Method [17], just to mention a few. In order to provide a technique amendable for the adopted pseudospectrabased approach towards uncertainty, a different Chebyshev polynomial-based approach is presented: as we will see, the Floquet multipliers are computed as the eigenvalues of a generalized eigenvalue problem which preserves linearity w.r.t. the entries of the nominal matrices and .
For systems like (13), where = , the proposed approach approximates the solution on the interval [− , ] by a piecewise polynomial in the form + 1 is the number of Chebyshev-distributed points, from which we can define the mesh { ( ) } +1 =1 on any = [ , ] (and viceversa) by using the following transformations The solution is approximated on each interval based on a linear combination of Chebyshev polynomials ( ) of the first kind with degree as and ( ) ∈ ℂ are the coefficients to be determined. In order to formulate the solution in terms of the polynomials, the following conditions are imposed: • Collocation constraint for (17) for = 1, … , + 1 • Continuity in = 0 gives 2 (0) = 1 (0), i.e., • Collocation constraint for the TPDDE (13) giveṡ for = 1, … , + 1, which using (21) gives for all = 1, … , . Here, we have applied the fact that ′ ( (1) ) = 2∕ anḋ ( ) = −1 ( ), where ( ) is the Chebyshev polynomial of the second kind of degree .
The conditions above result a generalized eigenvalue problem in the form (1) ] ⊤ ∈ ℂ 2 ( +1) and , ∈ ℝ 2 ( +1)×2 ( +1) . The coefficient matrices in (26) can be expressed as follows Robust stability of milling operations based on pseudospectral approach When and are smooth functions, the solutions emanating from an eigenfunction are smooth on [− , 0] and [0, ]; hence, they can be well approximated by polynomials on each of these intervals, ensuring spectral convergence of the eigenvalue approximations as a function of the degree . However, if and are only piecewise smooth, as might be the case for (12), the discontinuities in and impact the smoothness of the solutions and negatively affect the convergence rate of polynomial approximations. However, this problem can be overcome by a choice of piecewise polynomial consistent with the location of the discontinuities, while at the intersections only continuity is imposed. We refer the reader interested in a more detailed description on interval splitting to [19,22] or [24].

Pseudospectral method for robust analysis
In this section we assume that the generic system of timeperiodic DDEs (13) is affected by a set of real-valued scalar uncertainties, that result the perturbed equation in the forṁ where Δ ∈ ℝ, = 1, … , , are the uncertainties, and ( ), ( ) ∈ ℝ ( = 0, 1, = 1, … , ) are timeperiodic scaling matrices defining the structure and the weight of each uncertainty Δ (this will be clarified in the numerical experiments). We indicate each vector of uncertainties with the compact notation = [Δ 1 , … , Δ ] ⊤ ∈ ℝ . In addition, we set an upper bound equal to 1 on the absolute value of each Δ , i.e., |Δ | ≤ 1 for = 1, … , : observe that this can be done without loosing of generality, since we can assign different weights to each uncertainty using the scaling matrix-valued functions , . For the more complex case of matrix-valued uncertainties we refer again the reader to [19].
Each uncertainty Δ perturbing the time-periodic matrix or necessarily affects the collocation constraints for the TPDDE defined in (25) for = 1, … , + 1: thus we redefine the generalized eigenvalue problem corresponding to (26) as where ( ) indicates the perturbation of induced by parametric uncertainty Δ . The construction of the perturbation matrix is given in the Appendix A.
In order to analyse the worst-case scenario for this class of perturbations we look for the Floquet pseudospectral radius , i.e., the eigenvalue with largest modulus that the perturbed generalized eigenvalue problem (32) can attain. If the Floquet pseudospectral radius has modulus less than one, then we are guaranteed that the zero solution of (31) is asymptotically stable despite the uncertainties. Hence, the Floquet pseudospectral radius allows to assess the stability robustness. The worst-case analysis is performed by applying a (projected) gradient-ascent method in the space of uncertainties, thereby optimizing the modulus of the dominant Floquet multiplier: this approach is conceptually similar to the state-of-the-art approach ( [15,25]) in the worstcase analysis of time-invariant functional-differential equations, where the rightmost eigenvalue generated by a class of bounded perturbations is maximized.
To compute the Floquet pseudospectral radius we solve the optimization problem where is the dominant eigenvalue of perturbed problem ( ; ) = , for some ∈ ℂ 2 ( +1) .
This problem can be solved by simply using a projected gradient method in the space of 1-bounded perturbations : in the following we first construct the gradient of the dominant eigenvalue of the perturbed generalized eigenvalue problem (32) w.r.t. to uncertainties Δ . Let , ∈ ℂ 2 ( +1) be the left and the right eigenvector of with unit norms and such that ∶= * > 0, where * is the complex conjugate of . The derivative of | | 2 w.r.t. Δ as follows for = 1, … , . Note, that ( ) depends affinely on the uncertainty Δ , and therefore ( ) ∕ Δ is a (sparse) matrix, which does not change during the iterations. Therefore the calculation of in (35) is numerically efficient. A fast way to compute is given in Appendix A. We can now use a (projected) gradient method to solve the optimization problem (33) using the following fundamental step at each iteration where is the stepsize at iteration and ( ) is the ascent direction defined as follows where ( ) is the derivative of the largest eigenvalue of the generalized eigenvalue problem (32) perturbed with uncertainties Δ ( ) . Observe that ( ) is the projection of the derivative ( ) on the feasible set, which is needed to avoid the violation of the norm constraint.

Implementation
We here provide some more technical details about the implementation of the described method: the algorithm is stopped at some iteration such that the norm of the (projected) derivative ( ) is below a prescribed tolerance for = 1, … , ; of course, being a gradient ascent method, this algorithm is not guaranteed to converge globally rather than locally. Sticking to the state of the art (see [25,26]), including the following strategies in order to compute the globally dominant Floquet multiplier has proved to be effective: • The adaptive stepsize guarantees the monotonicity of the iterations: at each step , the initial step −1 is reduced by a scaling factor 2 until monotonicity is guaranteed; vice versa, if −1 already guarantees monotonicity, an extra computation with a doubled stepsize is carried out: the stepsize that guarantees the largest eigenvalue is then adopted.
• The restarting strategy: the method is run multiple times, each time initializing Δ (0) , = 1, ..., , from the left and right eigenvectors (see [15] for the details) corresponding to several dominant eigenvalues of the unperturbed problem.

Extension to modal parameter uncertainties
The application of the pseudospectral method relies on a numerically efficient computation of the gradient in (35). In order to apply the strategies discussed in this paper, we need the governing equation (including the perturbations) to be represented in the form of (31). In order to make the method easily applicable to system (12), we introduce a transformation in this subsection.
The modal parametes n, and (and their uncertainty) show up nonlinearly in the system (12), however, the latter can be transformed to a linear form by introducing a slack variable ( ) = [ n, ] ( ), and increasing the dimension of the state vector. Then, the linearity is established and the system reads as follows Observe that the leading matrix of the system is not the identity matrix anymore, and that it is perturbed by the uncertainty affecting : however, the methodology illustrated in Sec. 4 can be trivially extended to account for uncertain leading matrices as well. Regarding the computation of Floquet multipliers as introduced in Sec. 3, it is easy to see that in the definition of the first matrix in 22 in (30) must be replaced by (see the left-hand-side of (24), which has to be modified by introducing ). Also note that this extension requires two pairs of scaling matrices associated with the perturbation of n, , but the implementation of the pseudospectral method is not affected.

Case study and experiment
Experimentally determined stability limits are often different from the ones predicted based on a dynamical model. It is known that natural frequencies and damping ratios measured on the rotating spindle may differ from the ones measured on the idle tool, and this effect can be responsible for variations in the stability lobe diagrams. However, this characteristics can be measured prior to the experiments or predicted based on preliminary chatter tests, see [28,29] and especially [30] for the so-called inverse stability solution, and the references therein.
An experiment is carried out in an NCT EmR-610Ms milling machine. The cutting forces were measured by a Kistler 9129AA multicomponent dynamometer and data were acquired by NI-9234 Input Modules and an NI cDAQ-9178 chassis. Other devices are shown in Fig. 3. First, cutting tests were performed at two levels of depth of cut ( p = 1, 2 mm) and five feed rates ( = 0.02, 0.04, 0.06, 0.08, 0.1 mm/tooth) with full-immersion. Then cutting parameters were estimated by minimizing the difference between the measurements and the theoretical cutting forces. Second, we first performed an experimental modal analysis (EMA) on the real cutting tool during idle conditions, by measuring the responses and the excitations in directions and at the tool tip. The measured frequency response functions (FRFs) are presented in Fig. 4. The fitted FRFs are given in the form Robust stability of milling operations based on pseudospectral approach  [27].
where Φ , = , , ∕̂ 2 n, are the static compliances ( , = , ) and hat-symbol on̂ n, refers to the reference value of the -th natural frequency measured in idle state. The identified dominant modes and modal parameters are given in Table 1.
The cutter was replaced by a cylindrical dummy tool and a thorough operational modal analysis (OMA) was performed from 0 rpm to 13000 rpm. During the OMA the tool was hit by a small bullet and vibrations were measured by a highprecision laser sensor, see Fig. 3 (similar test are performed in a hardware-in-the-loop experiment in [29]). Then natural frequencies and damping ratios corresponding to the vibration modes of the real tool were identified, and the difference was scaled proportionally to match the modal parameters of the real cutter. The fitted parameters in the range 5000 -13000 rpm can be seen in Fig. 5, where black dots denote the fitted data, and dashed curves are second-order fitted characteristics in the presented region. The parameters of the fitted  characteristics are given in Table 2. The uncertainty bound was determined by calculating the standard deviation ( ) of the fitted parameters around the mean, which gave the ±3 gray region in Fig. 5. The uncertainty region was assumed to be independent of the spindle speed Ω. This results in Δ ∈ ℝ, = 1, ..., 8, scalar parametric uncertainties, that can be included in the pseudospectral method. Since Δ is a normalized variable (|Δ | < 1), the weight is introduced to represent the actual scaling of the uncertain parameters. Let be one of the parameter under considerations (given in Table 2 for this specific case), then its perturbed value can be written as + Δ . The weights can be used to construct the scaling matrices and in (31). The list of the speed-dependent modal parameters and their uncertainties determined from experiments are given in Table 2. The cutting force characteristics were identified from preliminary cutting tests, where the tangential and radial cutting coefficients and were both assumed to have 10% relative uncertainty. The uncertainties of the modal parameters were igure 5: Change of modal parameters along the spindle speed. (a) Natural frequencies (b) Damping ratios. Black dots denote the fitted modal parameters, dashed curves are the fitted characteristics and gray region is the constant ±3 uncertainty region around the fitted mean.

Table 2
Change of the parameters along the spindle speed ( (Ω) = 0, + 1, Ω + 2, Ω 2 , Ω ∈ [5000, 13000] rpm) and their (additive) uncertainty ( is independent from Ω). chosen to be equal to the 3 deviation, see Fig. 5 and the last row of Table 2. Note, that variation in the modeshapes are difficult to measure and are often neglected (see [30]). In this study we also assume that the modeshapes do not change and no uncertainty is added. The stability lobe diagrams corresponding to machining parameters given in Table 3 can be seen in Fig. 6(b), where dashed curves denote the stability boundaries obtained by neglecting speed dependency (modal parameters in Table 1), and solid black curves are the boundaries obtained by the speed-dependent parameters in Table 2. Fig. 6(a) presents the strengths of the harmonics of the chatter vibrations that can be calculated from the dominant eigenvalue and the Fourier transform of the corresponding eigenvector, see the method presented in [18]. In panel (a) black denotes the dominant chatter frequency c,d , while gray shading indicates the other harmonics with less energy. This diagram is used to identify the dominant frequencies from measurements. The robust stability lobe diagram was determined using the pseudospectral method. During the calculations, ten parametric uncertainties were assumed with weights given in the last column of Table 2.
Experimental chatter tests were carried out with fine resolution from spindle speed 5000 rpm to 13000 rpm. Accelerometers were mounted onto the workpiece and spindle housing, moreover a microphone was also placed close to the machining area. Chatter was identified based on the spectra of the signals measured during cutting. Each measured point in panel (a) was marked by a filled circle if the operation was stable, by a cross if it was unstable, and by a triangle if chatter was not clearly identifiable. When the operation was unstable, the dominant chatter frequency was selected and it was marked in Fig. 6(a), also by a cross. For quantitative chatter indicators obtained from the cutting force and acceleration signals, we refer the reader to the method presented in [31].
Four sample spectra are given in Fig. 6(c-f), which were captured by an accelerometer mounted on the spindle housing in direction . Panels (c-d) are measured at spindle speed 7600 rpm, with depth of cut 0.75 and 1 mm. From the signal of the accelerometer the spectrum of the velocity was calculated and the dominant peaks were compared to the predicted chatter frequencies. The absolute value of the FFT of the signal is indicated by dark gray curves, while dashed and point-dotted vertical lines indicate the tooth-passing frequency f tp = 60∕( Ω), its half, and their integer multiplies. The dominant chatter peak and the shifted harmonics are colored by black. It can clearly be seen in panel (d) that the dominant chatter frequency occurs at 2204 Hz, which is close to the third vibration mode. A similar test is shown in panels (e-f), where tests were carried out at speed 12000 rpm, at depth of cuts were 1 and 1.5 mm. The dominant chatter frequency was identified at 724 Hz, close to the predicted first mode.
Although small uncertainties were assumed initially in the model, the robust stable maximum depth of cuts at the lowest points decrease by nearly 30%. On the other hand by measuring the vertical differences, a 0.15-1 mm wide uncertainty band can be observed, which is also significant.
Experiments show a qualitatively good agreement with the predictions, but not all of the points inside the robust stable domain were indeed stable during the chatter tests. This indicates that there are still modeling errors and uncertainties present, which can explain this discrepancy. Possible reasons can be the linear force characteristics used during calculations or the dynamical model with proportional damping which are simplifications of reality. On the other hand the run-out of the tool, nonlinearities and stochastic effects may also have a strong impact on the stability boundaries.

Comparison with the MFS-SSV method
There exist several statistical methods in the literature to approximate the probability of stability if the probability density function of the uncertain parameters is known. However, there are only a few methods that guarantee stability if only upper and lower bounds on parameter values are available. These techniques are all called robust methods, such as the approach inferred from the Edge Theorem [13], the structured singular value analysis [14] or the pseudospectral method. These approaches cannot determine the probability of stability, but the calculation is typically less time-consuming.
The Edge Theorem is applied to time-invariant models in [13], but the extension to time-periodic systems is not straightforward, because the characteristic equation cannot be determined in a similar way. It also requires an affine de-pendence of the (scalar) characteristic equation on the uncertain parameters. The Multi-frequency Solution with Structured Singular Value analysis (MFS-SSV, [14]) provides an alternative way to assess the robust stable region, but this method also has drawbacks. It is only applicable to modal parameter uncertainties, and the approximations of robust stable regions can be more conservative than others. In this section we assume that only the modal parameters are uncertain (Table 2.) and compare the results obtained by the pseudospectral method and the MFS-SSV approach.
The Multi-frequency Solution of [32] can directly be applied to the measured frequency response functions ( ) without parameter identification. This technique has been extended in [14] in order to consider additive uncertainties in the FRFs in the form ( ) + ( ), where ( ) is a complex-valued uncertainty. The structured singular value analysis assumes that the uncertainty can be restructured into a block-diagonal matrix̂ , such that the governing equation is written as wherê ,̂ ∈ ℂ 2(2 +1)×2(2 +1) are constructed according to [14], and is the highest number of harmonics considered during the expansion. The system is not robust stable if the determinant in (40) is zero for somê and c . The computational algorithm presented in [14] gives a lower bound on sucĥ , which can be rephrased in terms of sufficient conditions for robust stability of the system. Figure 7 presents two scenarios, which compare the above detailed methods. Red curves indicate the robust boundary obtained by the pseuduospectral method, where = 22 and [0, ] is split into 3 subintervals. The iterations have been restarted from the 3 largest Floquet multipliers at each point of the parameter space, the maximum iterations were set to 50, but iterations were stopped if the error reached |Δ ( +1) −Δ ( ) | ≤ 10 −3 (for all ). The pseudospectral radius was calculated on a 120 × 40 grid on the plane (Ω, p ). The computation time was approximately 6 hours. Blue curves bound the approximation of the robust stable domain obtained by the MFS-SSV approach, where the highest number of harmonics considered in the approximation was = 30, and the grid in the parameter space (Ω, p , c ) was 120×40× 200 (this method requires a sweep along a third dimension also). The calculation time was approximately 50 minutes, which includes the generation on the uncertainty envelopes of the frequency response functions for all spindle speeds. We applied the Multi-dimensional Bisection Method presented from [33] in order to reduce the computational time.
Panels (a) and (b) in Fig. 7 show the different robust boundaries at different uncertainty levels. There is a very important remark. The pseudospectral method is not conservative, since the iterations are kept inside the domain of the allowable perturbations (resulting in an inner approximation of the robust stable region), and at the same time the globally dominant Floquet multiplier is computed. Meantime, the MFS-SSV method is always conservative, because the uncertainty envelope is bounded from above, leading to an outer approximation of the robust stable region. In our case (Fig. 7) the gap between the two curves is very narrow, which proves the accuracy of both methods in the present study. Note also, that the MFS-SSV method can be significantly more conservative in case of large uncertainties (see the explanations in [14]), and the method cannot handle cutting force model perturbations either. As opposed to this, the pseudospectral method is more flexible and can provide a reliable robust boundary.

Conclusion
An accurate prediction of dynamical stability of machining operations is limited by model simplifications and parametric uncertainties. In order to reduce the gap between predictions and actual measured stability limits, variations and uncertainties in modal parameters of the tool-tip were determined by means of experimental and operational modal analyses performed on the idle and rotating spindle. To reach the highest stable depths of cuts and avoid chatter, parametric uncertainties were considered and a pseudospectral method was applied.
The solver was tested in a case study, and stability predictions were compared to experimental chatter tests. The traditional stability lobe diagrams without considering the uncertainties predicted the dominant chatter frequencies well,  Table 2.
but the maximum stable axial depth of cuts were inaccurate. Robust stability boundaries were determined to include the effect of inaccuracies in modal parameters and in the cutting coefficients. The obtained new boundaries gave more reliable predictions than conventional methods, which proves the effectiveness of robust calculation. Another property of the pseudospectral method is that the basic algorithm converges locally to the worst-case perturbation by keeping the uncertainties inside the prescribed bounded region, and therefore inner approximations of the robust stable regions are obtained. This can also be viewed as a drawback, because the iterations may not converge globally. However, by incorporating restarts from several largest Floquet multipliers of the original problem (three in the case-study) a robust algorithm for computing the globally dominant Floquet multipliers is obtained, allowing to compute the actual robust stable region. In order to test the accuracy, the solver is compared to a robust structured singular value based analysis, which showed good accuracy for both methods.
In order to make the computations numerically efficient, the pseudospectral approach is applied to perturbed generalized eigenvalue problems, where uncertainties linearly affect the coefficient matrices. In this paper we treat scalar uncertainties only, however, the method is also effective to matrix-valued uncertainties, see [19]. Let us recall the perturbed time-periodic DDE in the forṁ where Δ ∈ ℝ, = 1, … , , are the uncertainties, and ( ), ( ) ∈ ℝ ( = 0, 1, = 1, … , ) are timeperiodic scaling matrices. These scales propagate through the collocations constraints and affect the generalized eigenvalue problem in the form wherê ,̂ ∈ ℝ 2 ( +1) ( = 1, … , + 1) are new scaling matrices that can be constructed from and . The optimization problem is solved by using a projected gradient method. Let again , ∈ ℂ 2 ( +1) be the left and the right eigenvector of with unit norms and such that ∶= * > 0, and let us definê