Stability of turning processes subjected to digital PD control

Stabilization of turning processes with a digital proportionalderivative feedback controller is analyzed. A one-degree-offreedom model of the turning process is considered. The control force is assumed to be acting directly on the tool. The sampling effect and the delay of the digital controller are involved in the model. The governing equation is a periodic delaydifferential equation, which includes a continuous point delay due to the regenerative effect of the material removal process and a discrete delay (i.e., a term with piecewise constant argument) due to the sampling effect of the controller. The principal period of the system is the sampling period. The stability analysis is performed using different implementations of the semidiscretization method. A series of stability diagrams are presented for different proportional and derivative control gains.


Introduction
In the last two decades machining technology went through a rapid development.The tendency shows an increase in the speed of manufacturing, hence a decrease in manufacturing time.High-speed manufacturing developed the need for the dynamical investigation of the cutting process.
Experience has shown, that at certain manufacturing parameters, self-excited vibration of the tool arises.The most commonly accepted explanation for this phenomena is the so-called regenerative effect which was first discussed by Tobias [11], Tlusty [10] and Kudinov [4].The point of this regenerative effect is that the cutting force is a function of the chip width which depends on not only the actual position of the tool but also on its position one turn before.Therefore, the mathematical description of the cutting processes results in Delay-Differential Equations (DDEs).Stability properties of the machining process are depicted by so-called stability lobe diagrams which plot the maximum stable axial depths of cut versus the spindle speed.
During cutting, self-excited vibrations of the tool are copied on the surface of the workpiece, which makes it wavy and worsens its quality.The time of the machining operation is influenced by the depth of cut, the spindle speed and the feed rate.The maximum available depths of cut for certain spindle speeds can be obtained from the stability lobe diagrams.One way to reduce the manufacturing time is the increase of the depth of cut.Several methods exist to increase the stable domains in the stability diagrams and to suppress regenerative chatter.Such methods are the vibration absorber [8], impedance modulation [6] and spindle speed variation [7].In each type of the cutting processes the increment of the damping ratio results in a growth of the stability domain.However, this increment is limited by the physical properties of the machine tool.
In control technology proportional-derivative (PD) controllers are widely used to stabilize linear unstable systems.Here, the proportional gain P works as an artificial stiffness and the differential gain D is a kind of artificial damping.PD control requires the measurement of the tool position and velocity.If the measurement and the feedback is continuous and delay-free, then the system can be described by Ordinary Differential Equa-tions (ODEs), and the stability analysis can be performed by standard techniques, like the Routh-Hurwitz criteria.However, in real applications, mostly digital controllers are used, where both the measurement and the feedback are non-continuous and delayed.
is usually modelled with a simple, linear into account the regenerative effect hence it leads to a DDE.For a non-linear linearised system is well known from the literature, this is This study deals with the theoretical investigation of the stability of the turning process under the effect of a digital PD controller.The application point of the control force is at the end of the cutting tool.The model is a linear 1 DOF model.The stability of this system is analysed numerically using the semidiscretization method [2].

Mechanical Model
The 1 DOF mechanical model of turning is shown by Fig. 1.The tool is assumed to be ideally elastic relative to the tool holder, with dominant mode parallel to the feed direction (x).The modal mass of the workpiece is m, the stiffness and damping are k and c respectively.The governing equation is The cutting force is given in the form where w is the depth of cut, h(t) is the chip thickness at time instant t and q is the cutting-force exponent.Parameter K depends on further technological parameters, and it is considered to be constant in the present analysis.Note that other formulas for the cutting force are also used in the literature [1], [3].If the spindle speed is constant, then the time delay within one turn is τ = 60/Ω, where spindle speed Ω is given in [rpm].The chip where v f is the feed velocity.If the control force is delay-free and continuous, than it can be given as where P is the proportional and D is the differential gain.Thus the governing equation can be written in the form This non-linear equation can be linearised around its trivial equilibrium point x 0 .The linearised equation reads where ξ = x − x 0 describes the displacement around the equilibrium point.After the introduction of the specific cutting force coefficient H = Kw(v f τ) q−1 /m, and specific control gains k p = P/m and k d = D/m, the governing equation can be rewritten in the form It can easily be seen that k p increases the stiffness and k d increases the damping of the system.Without control the stability diagram of ( 7) is a well known lobe diagram, which can be seen on Fig. 2. The boundaries of the stable region can be determined by the D-subdivision method [9].In parametric form these boundaries are

Digital control
The application of digital control with Zero-Order-Holder (ZOH) model [5] gives the equation where t ∈ [t j , t j+1 ), with ∆t = t j+1 − t j is the sampling period of the controller.This equation is a DDE with parametric excitation in the time delay, since the term ξ(t j−1 ) can be written as ξ(t − ρ(t)), where ρ(t) = ∆t − t j + t, t ∈ [t j , t j+1 ) is a sawtooth-like time-periodic time delay.The period of the parametric excitation is the sampling period ∆t of the controller.The first-order representation reads where The system is analysed by the semi-discretization method, which gives a finite dimensional approximation for the infinite dimensional problem.The sketch of the discretization for the case of zeroth-order one-point method is shown in Fig. 3.

One-point methods
One-point methods approximate the delayed value of the function with values taken from one discrete past time instant.The approximation of (11) for the time interval t ∈ [t i , t i+1 ) can be given generally as where D(t) is a weighting matrix which depends on the method and the order of the approximation.Short hand notation is used for x(t i−r ) = x i−r and respectively for the similar terms.
Here h = t i+1 − t i is the step size of the discretization, which is determined as τ = t i − t i−r = rh, where integer r is the resolution of the time delay.Integer p is the resolution of the sampling period determined as ph = ∆t.Integer κ = r/p is the ratio of the time delay and the sampling period.
The sketch of the zeroth-and the first-order one-point methods are shown in Fig. 4. The initial condition for (13) equation is x(t i ) = x i , which provides the continuity of the displacement and velocity functions at time instant t = t i .Using the method variation of constants the solution for (13) is Hence the relation between the two end points of the discretization step is where Equation ( 15) implies the following difference equation . .
In simpler form this reads as where X i is an augmented state vector, and G 1 is the coefficient matrix for this first step.If the control force were zero, the investigation of the eigenvalues of G 1 would be enough to determine the approximate stability chart of (10).However from Fig. 3 it is visible that the upper difference equation describes only the special case when the control force changes its value.To determine the approximate stability chart of (10) for a general case, the investigation of this periodic DDE is needed for a whole time period of the parametric excitation.The difference equation between the two ends of one sampling period is where Φ = G p G p−1 ...G k ...G 2 G 1 is the transition matrix (monodromy matrix) over the principal period ∆t = ph.Hence the condition of stability is that the magnitude of each eigenvalue of Φ must be less then 1, formally Until the end of the sampling period the control force keeps its value, thus the approximate differential equation for the second discretization step is where t ∈ [t i+1 , t i+2 ).After solving this differential equation similarly to (13), the difference equation between the endpoints of the second discretization step has the form . .
In this coefficient matrix (denoted by G 2 ) only the upper line is changed compared to G 1 : sub-matrix R C "jumped" one column right.With the induction of this phenomena the behaviour of G matrices for different discretization steps is shown by Fig. 5.

Zeroth-order
This method uses only the delayed function values as past values.The weighting matrix has the form (23)

First-order
It is visible from the structure of G matrices that the first time derivatives of ξ are introduced to the matrix only because the velocity is present in the control force.Without these derivatives G matrices can be reduced to half sizes, hence the presence of the derivatives drastically increase the computational time.With higher order methods computational time can be saved on an indirect way, by increasing the convergence of the numerical approximation.In this case, it is reasonable to apply higher order methods because of the "unused" first derivatives.This first order method uses not only the delayed function values but also their derivatives at past time instants.The weighting matrix is (24)

Two-point methods
Two-point methods take the past values from two past time instants for the approximation of the delayed function.The sketch of the applied two-point methods is shown in Fig. 6.The approximation of (10) for the time interval t ∈ [t i , t i+1 ) generally can be given as where the value of D 1 (t) and D 2 (t) weighting matrices depend on the order of the approximation and on the weighting between the past values.Similarly to (13) the solution of (25) can be determined by the method of variation of constants.The relation between the two endpoints of the discretization step is where The coefficient matrices G for the two-point methods have a similar form as in the case of one point methods.The only difference is that on the right end of the top row one more sub-matrix appears.The behaviour of matrices G is shown in Fig. 7.

Zeroth-order
This method simply takes the average of the function values in two past time instants.The weighting matrices are (28)

First-order
This method joins the function values in two past time instants with a strait line.Thus it uses linear interpolation between the two past values of the function.The weighting matrices are (29)

Second-order
This method approximates the function values between two past time instants by using not only the past values of the function but also the first derivatives of it.The approximation is done by the linear weighting of the tangentials calculated from the past values and derivatives.The weighting matrices are (30)

Third-order
This method approximates the function values between two past time instants with a third order polynomial.The coefficients of the polynomial are calculated by using the past values and derivatives.The weighting matrices are (31)

Results
Fig. 8 presents the stability boundaries for a turning process without any control.The negative depth of cut values are also presented.The exact stability boundaries are shown by black color.It can be seen that the higher the order of the approximation, the more accurate the stability boundaries.
Figs 9-10 present a series of stability diagrams for different (both negative and positive) proportional and derivative control gains for κ = 2 and 20.In these diagrams the exact stability boundaries of the turning process without any control are presented by gray color.Note that κ = r/p = τ/∆t, that is, it describes the ratio of the regenerative delay τ and the sampling period ∆t.
The stability diagrams were obtained numerically by analyzing the eigenvalues of the transition matrix Φ for a series of fixed parameters.

Conclusions
It can be seen in Fig. 8 that by increasing the number of steps of the semi-discretization method the calculated map approximates the real stability chart with higher accuracy.Higher order methods give better approximation at a given resolution of the time delay.It can be observed that the stability regions change with the ratio of the time delay and sampling period (κ).With certain control parameters for given κ values the stability region of the turning process can be increased by digital control.
Overall, it can be stated that the relation between the stability properties and the values of the proportional and the derivative gains is not trivial.For certain parameter combination, even negative control gains may stabilize the system, while for other parameters, positive control gains may destabilize the system.

Fig. 2 .
Fig. 2. Stability domain of the turning process without control force for ζ=0.05

Fig. 4 .
Fig. 4. Sketch of the zeroth-and the first-order one-point methods

Fig. 5 .
Fig. 5. Top row of G matrices for one point methods

Fig. 6 .
Fig. 6.Sketch of the applied two-point methods

Fig. 7 .
Fig. 7. Top row of G matrices for tw-point methods

Fig. 8 .
Fig. 8. Stability chart for the turning process without control force obtained by different approximation methods for different r values, ζ=0.05