Robust Stability Limit of Delayed Dynamical Systems

Computation of the stability limit of systems with time delay is essential in many research and industrial applications. Most of the computational methods consider the exact model of the system, and do not take into account the uncertainties. However, the stability charts are highly sensitive to the change of input parameters, such as eigenfrequency and time-delay. Furthermore, the computation of the dense stability lobe structure is numerically intensive; however, in some cases an envelope of these stability boundaries would be sufficient. A method has been developed to determine the robust stability limits of delayed dynamical systems, which is insensitive to fluctuation of system parameters. It is shown, that the resultant robust stability limits form the lower envelopes of the stability lobe structure.


Introduction
The determination of the stability of dynamical systems with time delay is of high importance for many industrial and research applications.The investigation of the dynamics of turning processes is such a typical industrial problem.In order to reach high process efficiency the material removal rate (MRR) has to be maximized.The major limitation for increasing the MRR is the so-called chatter vibration.The origin of these self-excited vibrations is the surface regeneration effect, which can be described by linear delayed differential equations (DDEs) [1,2].Another example is human balancing, which is often modelled by a simple inverted pendulum in the presence of reflex delay [3,4,5], which is a relevant issue to human motion control [6].A similar problem is the remote control of periodic robotic motions, when the delay in the information transmission system is not negligible [7,8].
The most important qualitative property of the corresponding dynamical system is the stability of the equilibrium or periodic motions.This is usually presented in the form of the socalled stability chart, that identifies those ranges of parameters where the linear system is stable.
In many cases, especially for machining operations at low spindle speed, the computation of the stability boundary requires very high computational effort and unnecessarily high resolution, due to the dense and sharp line segments of the stability boundary (see Fig. 1).In these cases, the computation of the lower envelope of the dense stability lobe structure would be adequate.
In case of time-domain computations [10,13,15,16,17], really high degree of discretization is required for proper results.In frequency domain computations [14,18,19,20], the determination of the very dense lobe structure in the unstable domain is unnecessary and also high resolution is required to determine the chatter frequency parameter, which is also a required computational step.
Traditional computational methods consider exact models of the system, and do not take into account the uncertainty of the input parameters; nevertheless, the results of these computations are highly sensitive to the change of input parameters such as eigenfrequency and time delay.In engineering practice, these dynamic parameters of a system can only be determined with a given uncertainty [22,23], moreover, parameters can change during operation [24].Some methods can be found in literature to determine the robust stability in which stability boundaries are computed for a large set of parameters, and the intersection of the stability ranges is determined.These algorithms lead to numerically intensive computations with many limitations, hence the "robust formulation cannot accommodate more than two varying parameters due to increasing model complexity" [25].
In the present study, a perturbation method is introduced to determine the robust stability limit.In Section 2, a well-known mathematical model, the delayed oscillator is used to introduce the method.In Section 3, the computation is applied to turning processes with process damping effect.The robust stability limits are determined in analytical form for both cases.

Stability of the delayed oscillator
Let us consider the second order delayed oscillator in the following traditional form [26]: where the dimensionless time delay τ = 2π .The stability chart is published in [26,27,28].
In the followings, the main steps of the stability computation based on the D-subdivision method [9,29] are summarised briefly.
The characteristic equation D(λ) of the dimensionless delayed differential equation Eq. ( 1) can be found by substituting the trial function x(t) = e λt as follows: The stability boundaries can be determined if the critical value of the root λ = iω c is substituted: In this case, a co-dimension 2 problem is defined by the real and imaginary part of the characteristic equation: in the three dimensional parameter space (δ, b, ω c ).The so-called Multi-Dimensional Bisection Method (MDBM) [30,31] is designed for these types of root-finding problems.This robust technique is able to find the submanifolds of the roots of a system of non-linear equations in any high dimension and co-dimension.The MDBM can be used for the determination of multiple boundary curves and it can even find the closed curves of stable and unstable islands in the stability chart automatically.The roots of Eq. ( 4)-( 5) are determined by the MDBM and are presented in Fig. 2 for different values of the damping ratio κ.In the top view, the resultant lines determine the stable area.
For the undamped system κ = 0, the stability boundaries are straight lines with slope +1 and -1.For κ > 0, the stability boundaries are not straight lines only.The δ = κ line is associated to saddle-node instability at ω c = 0, all other boundary curves represent Hopf instabilities (see [32]).

Parameter uncertainty
Equation ( 1) can be considered as a dimensionless form of the governing equation of a mass-spring system with delayed control [7,9,26,32].The previous computation method considers the exact model of the mechanical system, and does not take into account the uncertainty of the input parameters, such as eigenfrequency and feed-back delay.These parameters have large influence on the dimensionless time delay, which is the main source of the stability problem.Hence, small differences in these parameters could lead to significant change in the stability chart.To represent the effect of uncertainty in the delay parameter τ, numerous stability charts were calculated for a set of time delays in the range of [0,4π] and these were plotted together in Fig. 3.The intersection of the stability areas can be used as an approximation of the robust stability region [25].See the resultant graph in Fig. 3.

Robust stability limit
The robust stability limit should be computed for a continuous variation of the time delay, however, the computation for obtaining the results presented in Fig. 3 was time consuming even for 30 different τ values only.
It was found, that this set of boundary curves can be connected to form one surface by means of introducing an additional parameter.Let us define the regenerative phase shift ϕ := τω c and consider this as an independent extra (time-delay perturbation) parameter in the exponential term in Eq. (3): The usage of the MDBM is essential to solve the resultant co-dimension 2 problem in the extended 4 dimensional parameter space (δ, b, ω c , ϕ).The resultant surface is plotted in Fig. 4. 2) The envelope of the surface, which forms the robust stability limit is also plotted with a thick black line.
The robust stability limit is defined by the envelope of this surface, where the surface segments are parallel to the ω c axis in the 3D representation in Fig 4 .Based on the previous fact, it can be deduced, that in the vicinity of these parameter points, the real part of the roots λ of the characteristic equation (2) does not change as a function of the perturbation parameter.This condition can be described as follows: Re 0 The left hand side can be determined by the implicit derivation of Eq. ( 6) [9] d d Now, the extra condition Eq. ( 7) for robust stability can be written as follows: ( ) During the numerical implementation, it is better to use a rearranged form to eliminate the division in which the critical value λ = iω c is already considered.
Figure 5 shows the traditional stability chart and the robust stability limit which is determined by the MDBM as a codimension 3 problem formed by Eq. ( 4), ( 5) and ( 11) in the extended 4 dimensional parameter space (δ, b, ω c , ϕ).In case of 65 grid points along each parameter dimension, the time of computing the traditional stability limits was 1.3 s, while the robust stability limit was computed in 4.4 s (Matlab 2013b; Intel Core i7-2620M CPU 2.70GHz, 4GB Memory).Compared to the much higher complexity of the robust stability computation, this slight growth in the computational time is negligible.The high computation efficiency of the MDBM is described in details in [30] and it is characterised by the efficiency number.

Robust stability of the turning process
In the following section, the above method is applied to orthogonal turning with process damping effect (see Fig. 6) [2,21].The governing equation of the system is given by: where the parameters are: mass m, damping c, process damping coefficient C, stiffness k, cutting coefficient k 1 , chip width w, feed per revolution f 0 , time delay τ = 2π ⁄ Ω and spindle speed Ω.The number of parameters can be reduced by using dimensionless parameters (natural frequency n k m and dimensionless time n t tω =  .Rearranging Eq. ( 12), and separating the static tool deflection, the timedelayed governing equation which describe the stability properties is written as follows [9]: , the computational results obtained by the MDBM are plottedin Fig. 7.In the next step, it is important to note, that the c τω  component is substituted by the perturbation parameter ϕ only in the exponential term, which represents the time-delay.) is plotted in Fig. 8.The robust stability limit is shown in Fig. 8, too, for which the condition in Eq. ( 11) is also considered.The final result of the stability computation should be given in the form shown in Fig. 9, where the traditional stability chart and the robust stability limit are represented in the same graph.This combined chart could be well applied in industrial practice, while the safe range of parameters is presented by the robust stability limit together with the stability pockets of the traditional instability-lobe structure.In these pockets large MRR [33] can be achieved, but one has to be aware of the risk of parameter uncertainty which can affect the position of the lobes.

Analytic results
Analytical result can be evaluated for the turning model with process damping.The corresponding system of equations is given by: This system can be solved analytically, and the result which is valid for positive chip width values is: ( In case of zero process damping coefficient 0 C =  , Eq. ( 16) provides a straight line along the spindle speed Ω .This line connects the minimum points of the lobes, which is at 2 (1 ) w ξ ξ = +  as shown in [9].

Conclusion
In the present study, I have developed a perturbation method for calculating the robust stability limits of delayed differential equations.The method eliminates the essential parameter sensitivity of the stability charts.The determined curve forms the lower envelope of the lobe structure, thus it can be used in case of inaccurate input parameters, such as: natural frequency or time-delay.In the proposed computation method, an additional perturbation parameter must be introduced, which increases the dimension of the parameter space by one, nevertheless, the applied Multi-Dimensional Bisection Method can overcome this problem without significant increase of computational time.
In addition, in case of the turning process, the presented method also speeds-up the computation at lower spindle speed range, moreover, smaller resolution in the parameter space is sufficient.
The robust stability limit can be determined analytically in closed form equations in special mechanical models, as it is shown in this paper for the delayed oscillator and for the single degree of freedom turning model with traditional process damping effect.
In the future, we plan to extend this method to dynamical systems with multiple and/or distributed delays.

Fig. 1
Fig. 1 Stability diagrams for low spindle speeds range a) tool geometry optimization [17]b) measured and computed stability limits with process damping[21].

Fig. 2
Fig. 2 Stability chart of the delayed oscillator system and their 3D view with the critical angular frequency.Black, blue and red lines denote κ = 0,0.1,0.2 and the corresponding stable areas are shaded with gray, bright blue and pink, respectively.

Fig. 3
Fig. 3 Stability chart of the delayed oscillator for a set of different time delays τ = [0,4π] (black lines).The original chart (τ = 2π) is denoted by red thick lines.In the δ − b plane the robust stability area is shaded with green.(κ=0.2)

Fig. 4
Fig. 4 Surface of the connected stability boundaries of the delayed oscillator system and their 3D view with the extra independent regenerative phase shift parameter ϕ. (κ = 0.2) The envelope of the surface, which forms the robust

Fig. 5
Fig. 5 Stability boundaries of the exact model of the delayed oscillator (black line) and its robust stability limit (colored line) and their 3D view.In the δ − b plane the area of robust stability and stable area of the exact model are shaded by green and blue, respectively.(κ = 0.2)

Fig. 7
Fig.7 Stability boundaries and their 3D view with the dimensionless chatter frequency.The stable area is shaded with blue.(ξ = 0.05, 0 001 C = . ) The resultant surface in the extended 4 dimensional parameter space ( c w ω φ Ω, , ,  