High speed CNC system design. Part II: modeling and identification of feed drives

International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

High speed CNC system design. Part II: modeling and identification of feed drives Kaan Erkorkmaz, Yusuf Altintas

*

The University of British Columbia, Department of Mechanical Engineering, 2324 Main Mall, Vancouver, BC, Canada V6T 1Z4 Received 30 October 2000; accepted 18 December 2000

Abstract Accurate modeling and identification of the feed drives’ dynamics is an important step in designing a high performance CNC. This paper presents a method for identifying the dynamic parameters, as well as the friction characteristics of machine tool drives. The inertia and viscous friction are estimated through an unbiased least squares scheme. The friction model is developed by observing the disturbance torque through a Kalman filter, while jogging the axes under closed loop control at various speeds. The overall axis model is used in designing a high speed feed drive control system, which has been presented in Part III of this paper. As verification of the identified friction model, contouring test results without and with friction compensation are also presented.  2001 Elsevier Science Ltd. All rights reserved.

1. Introduction The design of a high performance feed drive control system requires accurate knowledge of the axis dynamics. The feedback controllers need to be designed to impose the same closed loop response on all axes, in order to avoid contouring errors in linear motion [1]. On the other hand, the feedforward control techniques employed to widen the axis tracking bandwidth typically rely on cancelling as much of the open or closed loop dynamics as feasible in a stable manner, and make use of the future reference points known in advance in order to compensate for the phase lag introduced by the remaining axis dynamics [2–4]. Pritschow and Philipp [2] have shown that the performance of such feedforward controllers can significantly degrade in the presence of parameter mismatch, therefore the axis parameters have to be accurately known in order for the feedforward control scheme to be successful. As a solution to this problem, Tung and Tomizuka * Corresponding author. Tel.: +1-604-822-2182; fax: +1-604-822-2403. E-mail addresses: [email protected] (K. Erkorkmaz), [email protected] (Y. Altintas).

0890-6955/01/$ - see front matter  2001 Elsevier Science Ltd. All rights reserved. PII: S 0 8 9 0 - 6 9 5 5 ( 0 1 ) 0 0 0 0 3 - 7

1488 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

[5] have proposed designing the feedforward controller based on an experimentally identified model of closed loop axis dynamics. Their identification was done using a frequency weighted least squares algorithm with a unity d.c. gain constraint. Although a robust feedback control law is designed to include some sort of disturbance rejection mechanism, such as simple integral action or a more sophisticated disturbance observer, it takes a certain amount of time before a discontinuous change in friction force can be estimated and compensated for. This phenomenon is significant in machine tool applications where the contour accuracy degrades at corners and arc quadrants, as one of the axes changes its direction of motion. A counter measure for this problem is to use the readily available trajectory information in the interpolator for feedforward compensation of expected friction forces [6–8], before they have a chance to degrade the tracking accuracy. However, this technique requires accurate knowledge of the friction characteristics, whereas an incorrect compensation may also worsen the contour following properties. This paper presents a method for modeling and identifying the dynamics of feed drives. Together with the classical servo model in [1,9] the dominant nonlinear effects of friction [6] are also considered, and the least squares parameter identification [10] is formulated to include the Coulomb friction as an extra parameter, enabling the parameters of the linear model, to be estimated without bias. Then, the friction model is refined by jogging the machine tool axes under closed loop control at various velocities, and monitoring the actual friction using a Kalman filter [11] based disturbance observer. The identified axis model has been used in designing the high speed feed drive control system, which has been presented in Part III of this article. The remaining of this paper is organized as follows: Modeling of the linear dynamics, as well as the nonlinear friction effects are presented in Section 2. This is continued by the parameter identification in Section 3. High speed contouring test results, which verify the experimentally identified friction model are presented in Section 4, and the contributions are summarized in Section 5. 2. Modeling of feed drives In this section, a simple linear model for feed drives is derived. Noise and quantization properties, as well as an adequate disturbance model for observing the friction is also included into the formulation. Then the expected nonlinear friction characteristic is summarized from [6], and this friction model is combined with the linear dynamics previously derived. 2.1. Linear dynamics The simplified linear dynamics of a typical feed drive, can be represented as shown in Fig. 1. Here, u⬘ (V) is the control signal applied to the input of the current amplifier, denoted with the gain Ka (A/V). The amplifier produces the current i (A) in the motor armature, resulting in the motor torque Tm (Nm), which is linearly proportional to i with the motor torque constant Kt (Nm/A). In addition to the torque produced by the current running through the armature, the motor shaft is also subject to a disturbance torque Td (Nm) which contains the effect of friction in the guideways, ball screw, and bearings, as well as the cutting forces [9]. The difference between Tm and Td is used to actuate the mechanical system consisting of the equivalent inertia

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1489

Fig. 1. Linear dynamics of a feed drive.

J (kg m2) and viscous damping B (kg m2/s) of the corresponding axis, with respect to the motor shaft. This in return leads to an angular velocity w (rad/s) in the motor shaft, and through the lead screw gain rg (mm/rad), to axis velocity x˙ (mm/s), and its integral which is axis position x (mm). If there are any gear ratios in the transmission train from the motor shaft to the axis position, their effect would also be included into rg. For linear motor driven feed drives, rg=1, J and B are replaced by the translational axis mass and viscous damping values. It should be noted that a velocity control loop is not yet closed around the system, hence the amplifier is configured to operate in current mode, meaning that it tries to regulate the armature current i in a proportional manner to the control signal u⬘ (i.e. i=Kau⬘). Although there is physically a current feedback loop between the amplifier and the motor armature, its typical bandwidth is around a few kiloHertz, allowing the modeling of the current amplifier as a linear gain, as the discrete-time sampling frequency in most CNCs is around 1–2 kHz. Furthermore, although back e.m.f. is physically present in the motor, having a current regulation loop in the amplifier compensates for the torque loss due to the back e.m.f., allowing its effect to be neglected in the model. This would not have been valid if the amplifier were voltage controlled, and not current. Considering Fig. 1, the axis position x(s) can be written in terms of the control signal u⬘(s) and disturbance torque Td(s) in the Laplace domain as, rg 1 (K K u⬘(s)⫺Td(s)) x(s)⫽ · s Js+B t a

(1)

The identification techniques presented in the paper use the control signal input u⬘(s) to excite the axis dynamics, and the motor mounted encoder and tachometer to measure the motor shaft angular position and velocity. For modeling convenience, the tachometer measurements have been calibrated to yield the motor shaft velocity in rad/s, and for axis control convenience, the position measurements have been calibrated to correspond directly to the table position in mm. The disturbance torque Td(s) can be written as a combination of two components as, Td(s)⫽Tf(s)⫹Tc(s)

(2)

where Tf(s) and Tc(s) correspond to the effects of overall axis friction and cutting forces respectively, reflected on the motor shaft. To facilitate easier disturbance estimation, Td can be assumed to act at the control signal stage, as d(s)=Td(s)/(KtKa). Furthermore, defining the gain Kw=KtKa/J and pole pw=⫺B/J, the shaft velocity w can be expressed in terms of control input u⬘ and equivalent input disturbance d as,

1490 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

w(s)⫽

Kw (u⬘(s)⫺d(s)) s−pw

(3)

and the axis position x can be written as, rg Kw rg (u⬘(s)⫺d(s)) x(s)⫽ ·w(s)⫽ · s s s−pw

(4)

The above equations can be expressed in state space notation as,

冋 册 冋 册

冋 册

冋 册 冋 册

x˙ (t) x(t) u⬘(t) 0 rg 0 ⫽Ac ⫹[Bc⫺Bc]· , where Ac⫽ , Bc ˙ (t) 0 pw w w(t) d(t) Kw

(5)

The axis control signal u⬘ is generated by a D/A converter, hence the velocity expression in Eq. (3) can be rewritten in the discrete-time domain with a zero order hold at the input stage as [12], Kw Kwd (u⬘(k)⫺d(k)), where Kwd⫽ (1⫺ePwTs), pwd⫽ePwTs w(k)⫽ z−pwd −pw

(6)

where k is the sample counter, Ts (s) is the discrete time sampling period, Kwd is the discrete time transfer function gain, pwd is the discrete time pole and z is the forward shift operator. The discrete time version of the state space expression in Eq. (5) can be obtained as [12].

冋

x(k+1) w(k+1)

册 冋 册 ⫽Ad

x(k)

w(k)

⫹[Bd⫺Bd]

冋 册 u⬘(k) d(k)

冕 T

, where Ad⫽eAcTs, Bd⫽ eAcl dl·Bc

(7)

0

On top of the linear dynamics covered in the above equations, there will be various nonlinearities such as quantization errors, saturations, as well as sources of noise. As most axis controllers are designed to operate in the linear characteristic range of the drives, it can be assumed that the control signals will be below the saturation limits, which can be accomplished by setting limits to the acceleration and jerk profiles of the reference trajectory (see Part I of this article). On the other hand, the quantization and measurement noise errors will always be effective, and will dominate the achievable tracking accuracy. Hence, they need to be considered in both the identification, as well as the controller design stages. The main sources of error are shown in Fig. 2, where the desired control signal u is applied to the amplifier as u⬘ after being quantized by the DAC which has a resolution of du. The true axis position x is measured as xm using an encoder with a resolution of dx, and the motor shaft velocity w is measured as wm, through a tachometer with a measurement noise variance of Rw. The smallest difference in control voltage that can be generated by the DAC is du, hence the actual control signal u⬘ will lie in a range of ⫺du/2 to du/2 of the commanded control signal u. The quantization noise u˜ which can be defined as,

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1491

Fig. 2. Simplified axis dynamics including quantization and measurement noise.

u⬘(k)⫽u(k)⫹u˜ (k)⇒u˜ (k)⫽u⬘(k)⫺u(k)

(8)

has a uniform probability distribution function (p(u)) between ⫺du/2 and du/2 (Fig. 3), which can be expressed as, p(u˜ )⫽

再

1/du, −du/2ⱕu˜ ⱕdu/2 0,

(9)

otherwise

Considering that u has zero mean (E[u]=0) the variance (Ru˜ ) of this distribution can be calculated as,

冕 ⬁

Ru˜ ⫽E[(u˜ ⫺E[u˜ ])2]⫽E[u˜ 2]⫽

冕

du/2

p(u˜ )·u˜ 2·du˜ ⫽

⫺⬁

(du)2 1 2 ·u˜ ·du˜ ⫽ du 12

(10)

⫺d/2

Similarly, the position measurement noise x˜ defined as, xm(k)⫽x(k)⫹x˜ (k)⇒x˜ (k)⫽xm(k)⫺x(k)

(11)

will also have a uniform distribution between ⫺dx/2 and dx/2 with a variance of, (dx)2 Rx˜ ⫽ 12

(12)

Fig. 3.

Probability density function of control signal error.

1492 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

˜ is defined as, The tachometer measurement noise w ˜ (k)⇒w ˜ (k)⫽wm(k)⫺w(k) wm(k)⫽w(k)⫹w

(13)

and is assumed to have a normal distribution, with a variance Rw˜ estimated by computing the variance of the tachometer measurements (in rad 2/s2) while the axis is powered but not moving. In order to design an unbiased state estimator, the disturbance characteristic should also be considered in the dynamic model. In axis control applications the dominant source of disturbance is friction, which does not change very abruptly during steady state feed motion. Hence, it is reasonable to assume that the disturbance d(k) consists of piecewise constant signals with zero mean white noise perturbation wd(k) as below, d(k⫹1)⫽d(k)⫹wd(k)

(14)

where the variance of the disturbance perturbation Rwd will be a tuning parameter to determine the “activeness” of the Kalman Filter. An augmented discrete time state space model that includes the linear dynamics, input and measurement noise, as well as the disturbance model can be obtained by combining Eqs. (7, 8, 11, 13) and (14) as,

冤 冥冤 冥 x(k+1)

冋 册

x(k)

w(k+1) =A w(k) +B[u(k)]+W d(k+1)

冋 册 xm(k)

d(k)

冤 冥冋

册

x(k)

x˜ (k) =C w(k) +V ˜ (k) wm(k) w d(k)

u˜ (k)

wd(k)

冧

(15)

where,

冋

A⫽

Ad 0

−Bd 0 1

册 冋册 冋 册 , B⫽

Bd

0

, C⫽

1 0 0 0 1 0

冤 冥 Bd 0

, W⫽

冋 册

0 , V⫽

0

1

1 0

0 1

(16)

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1493

In the above model, the disturbance d(k) appears as a state, as it is to be estimated by the Kalman filter, both in friction identification, as well as feedback control. A and B are the augmented system and input matrices respectively. W governs how the process noise vector [u˜ wd]T affects the state transition. C is the output matrix, and V relates the measurement noise ˜ ]T to the measured position and velocity [xm wm]T, which is also referred to as the output [x˜ w vector. The above state space model is used in the Kalman filter design, both for friction identification, as well as feedback control purposes. 2.2. Friction Armstrong et al. [6] have presented an excellent survey on the physics behind the friction phenomenon, as well as compensation techniques of dealing with it. Some of the significant characteristics of friction are summarized below, and incorporated into the linear axis dynamics. The typical friction characteristics for lubricated metallic surfaces in contact can be described by the Stribeck curve, as shown in Fig. 4. If a tangential force is applied to the surfaces, it will first work to elastically deform the asperity junctions (local spikes which are in contact between the two surfaces). This phenomenon is referred to as presliding displacement or microslip. If the tangential force exceeds a certain threshold, the junctions will break, causing sliding to start. This threshold is referred to as the static friction value. It is noted in tribology literature that static friction level will be a function of dwell time, which is the duration that the surfaces are at rest before sliding occurs. Once the breakaway occurs, a film of lubricant will not be able to build up between the contact surfaces at very low velocities. In this case, sliding will occur between solid boundary layers of lubricant that are stuck to the metal surfaces. This regime of the Stribeck curve, is referred to as boundary lubrication. As the sliding velocity between the two surfaces increases, more lubricant

Fig. 4.

The relation between sliding velocity and friction force for two lubricated metallic surfaces in contact.

1494 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

is drawn into the contact zone, which allows a lubricant film to be formed. At this stage, the film is not thick enough to completely separate the two surfaces, and the contacts at some asperities still affect the friction force. This regime is named as partial fluid lubrication. As partial fluid lubrication increases, solid to solid contact between the boundary layers decreases, which results in the reduction of friction force with increasing velocity. Partial fluid lubrication is inherently an unstable regime; with increasing velocity, the lubricant film gets thicker, hence reducing the friction force, and causing the velocity to increase further. This regime is difficult to model, as it involves the interaction of elasto-hydrodynamic phenomena with surface roughness properties. Furthermore, the changes in friction force do not immediately follow the changes in velocity or load conditions, but there is a phase lag between the two, referred to as frictional memory which may be in the order of milliseconds to seconds. After sliding velocity reaches a certain level, a continuous fluid film is formed which completely separates the two surfaces. In this regime, referred to as full fluid lubrication, the viscosity of the lubricant is dominant on the friction force. If the two surfaces are conformal (i.e. they have same radius of curvature) hydrodynamic lubrication occurs. If the surfaces are nonconformal, elastohydrodynamic lubrication occurs, which usually results in larger friction forces and more wear compared to hydrodynamic lubrication. Neglecting the time dependent properties of friction, such as the increase of static friction with dwell time, and the phase lag between varying load and velocity conditions and its outcome in friction force, the friction characteristics similar to those shown in Fig. 4 can be incorporated into the axis dynamics in Fig. 1, which results in the block diagram shown in Fig. 5. In this case, as the equations of motion are written according to the motor shaft, the friction is considered to be a part of the disturbance torque Td, as in Eq. (2). Note that the effect of viscous friction has been included into the “friction dynamics” block. The friction torque Tf is considered to be dependent on axis velocity w, and actuation torque T which is what remains of the motor torque Tm after a part of it has been used to overcome the effect of cutting forces Tc. Due to the self-locking characteristic of the ball screw, Tc will only be felt when the axis is in motion, whereas this is not the case for direct drives. Velocity dependence arises from the nature of the Stribeck curve (Tst(w)). Actuation torque dependence has a decision making role in determining whether the static friction torque will be overcome to initiate motion, when the axis is at rest. Hence, the expression for friction torque Tf may be written as,

冦

0

if w=0 & T=0

T

if w=0 & T⬎0 & T⬍T +stat

T

if w=0 & T⬍0 & T −stat⬍T

Tf(w,T)⫽ T +stat

if w=0 & T⬎0 & T⬎T +stat

T −stat

if w=0 & T⬍0 & T −stat⬎T

(17)

T +st(w) if w⬎0

T −st(w) if w⬍0

The first case of Eq. (17) attributes to no actuation and when the axis at rest. The second and

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1495

Fig. 5. Interaction of friction with axis dynamics.

third cases deal with the actuation torque T not being large enough to overcome the static friction Tstat. The fourth and fifth cases handle the beginning of motion when the axis is at rest and the actuation torque is large enough to overcome the friction torque. The last two cases account for friction properties when the axis is in motion. In Eq. (17), T +st(w) refers to the part of the Stribeck curve that corresponds to positive velocity and T −st(w) to negative velocity (similar to f +s and f −s in Fig. 4). The static friction values T +stat and T −stat are obtained from these two curves as, T +stat⫽ lim T +st(w), T −stat⫽ lim T −st(w)

(18)

w→0−

w→0+

In practice T +st(w) and T −st(w) can be characterized analytically with sufficient closeness as, +

+

−

−

冎

T +st(w)=T +state−w/⍀1+T +coul(1−e−w/⍀2)+T +viscw T −st(w)=T −state−w/⍀1 +T −coul(1−e−w/⍀2)+T −viscw

(19)

where Tstat (Nm) is the static friction torque, Tcoul (Nm) is the Coulomb friction torque, and Tvisc [(Nm)/(rad/s)] is the viscous damping coefficient for the corresponding direction of motion (positive or negative). The velocity constants ⍀1 and ⍀2, determine the velocity spacing between the different friction regions. The above model has been used for experimental friction identification in this work.

3. Identification Using the models derived in Section 2, the x and y axis dynamics of a vertical CNC machining center were identified, and the experimental results obtained are presented in the following. 3.1. Unbiased estimation of inertia and viscous friction The discrete time axis dynamics given in Eq. (6) has only two parameters Kwd and pwd which are simple to identify using the Least Squares (LS) technique, if the disturbance d(k) is neglected. In a straight forward manner, assuming that d(k)⬵0, the difference equation relating to Eq. (6) may be written as, w(k)⫽Pwd·w(k⫺1)⫹Kwd(u⬘(k⫺1)⫺d(k⫺1))

(20)

1496 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

and assuming that w(k)⬵wm(k), Eq. (20) can be expressed for all N samples collected with a sufficiently exciting input u⬘(k) as, Y⫽⌽q⫹E

(21)

where,

冤 冥 冤 w(2)

Y⫽

w(3) …

, ⌽⫽

w(N)

w(1)

u⬘(1)

w(2)

u⬘(2)

…

…

冥

, q⫽

冋 册 pwd

Kwd

(22)

w(N−1) u⬘(N−1)

In the equation above, Y is the output vector, ⌽ is the regressor matrix, q is the parameter vector to be estimated, and E is the prediction error vector. The objective in the parameter estimation process is to find the parameter estimate vector qˆ that minimizes the cost function, 1 V(qˆ )⫽ (Y⫺⌽qˆ )T(Y⫺⌽qˆ ) 2

(23)

which penalizes the difference between the predicted output ⌽qˆ and the measured output Y. The optimal parameter estimate vector to minimize V(qˆ ) is given as [10], qˆ ⫽(⌽T⌽)−1⌽TY (24) The estimates for pwd and Kwd can be extracted from the first and second elements of qˆ , and using the definitions in Eq. (6) (Kw=(KtKa)/J, pw=⫺B/J), the physical inertia J and damping values B can be estimated as, (pˆ wd−1)KtKaTs ˆ (1−pˆ wd)KtKa , B⫽ Jˆ ⫽ ˆ Kwd ln(pˆ wd) Kˆ wd

(25)

Although the above approach is easy to implement, it is not realistic to totally neglect the nonlinear effect of friction d(k) in the parameter estimation. Furthermore, the measurement noise, even if it is white, may result in a colored noise sequence as it passes through the process dynamics, in the estimation process. Both of these factors will cause bias to the parameter estimates. Recalling from Eq. (13) that the measurement of axis velocity wm(k) is assumed to be corrupted ˜ (k), and from Eq. (8) that the effective input signal u⬘(k) is with uncorrelated Gaussian noise w also corrupted with uniformly distributed noise u˜ (k) as, ˜ (k) wm(k)=w(k)+w u⬘(k)=u(k)+u˜ (k)

冎

(26)

Eq. (20) may be rewritten in terms of the measured velocity values as, ˜ (k)⫺pwdw ˜ (k⫺1)⫹u˜ (k⫺1) wm(k)⫽pwdwm(k⫺1)⫹Kwdu(k⫺1)⫺Kwdd(k⫺1)⫹w

(27)

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1497

The Kwdd(k⫺1) term represents the friction, which varies with velocity w. Unless its effect is considered in the parameter estimation, it will result in a significant bias to pωd, which translates ˜ truly is a white noise process, to the estimated viscous friction B. On the other hand, even if w as can be seen from Eq. (27) that it enters the estimation process through the filter (1⫺pwdz−1), due to the fact that the true values of w are being substituted by the available noisy measurements wm·u˜ is still a random uncorrelated sequence, which would not result in a bias in the estimates. Of these bias sources, the disturbance term Kwdd(k⫺1) is the most severe one that needs to be considered in the identification process. In the following, the simple least squares estimation technique will be reformulated to also include the effect of Coulomb friction. Although the colored ˜ (k) would best be handled using a Maximum Likelihood Estimator noise caused by (1⫺pwd·z−1)w (MLE), in this article LS estimation has been used favoring its simplicity and its ability to cope with the Kwdd(k⫺1) term, which is the major source of bias. It can be assumed that Coulomb friction is the dominant source of nonlinear friction during motion, and that it enters the plant through the disturbance channel d as,

冦

if w(k)=0

0

df(w(k))⫽ d

+ f

if w(k)⬎0

(28)

d −f if w(k)⬍0

Noting that there is no machining being done during the identification tests (i.e. Tc=0 ), d(k)⫽df(k)

(29)

In Eq. (28), d +f and d −f correspond to the positive and negative velocity Coulomb friction values as, d +f⫽T +coul/(KtKa), d −f ⫽T −coul/(KtKa).

(30)

As the model in Eq. (28) is to be included into a parameter estimation process where the noisy measurements of velocity wm, are being used instead of the true values, a dead band at ⍀d, just above the maximum noise level can be used to distinguish between the true axis motion and noise. Defining a signum function with deadband ⍀d as,

冦

0

s(wm, ⍀d)⫽ 1

if 兩wm兩ⱕ⍀d if wm⬎⍀d

−1 if wm⬍−⍀d

(31)

1498 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

the friction model in Eq. (28) can be rewritten as, df(wm(k))⫽PV(wm(k))·d +f⫹NV(wm(k))·d −f

(32)

where, 1 PV: Positive Velocity= ·s(wm(k))·(1+s(wm(k))) 2

冧

(33)

1 NV: Negative Velocity=− ·s(wm(k))·(1−s(wm(k))) 2

˜ (k⫺1)+u˜ (k⫺1)) term, and including ˜ (k)⫺pwdw Eq. (27) can now be rewritten neglecting the (w the friction model as,

冤 冥 pwd

wm(k)⫽[wm(k−1) u(k−1) −PV(wm(k−1)) −NV(wm(k−1))]

Kwd

Kwdd +f

(34)

Kwdd −f

In a similar approach used for the simple LS problem before, Eq. (34) can be rewritten for N samples of data collected as, Y0⫽⌽0·q0⫹E0

(35)

where Y0 is the output vector, q0 is the parameter vector, ⌽0 is the matrix of regressors, and E0 is the prediction error vector. Y0, ⌽0, and q0 are given as, Y0⫽[wm(1) wm(2) … wm(N)]T, q0⫽[pwd Kwd Kwdd +f Kwdd −f ]T

(36)

and

⌽0⫽

冤

wm(1)

u(1)

−PV(wm(1))

−NV(wm(1))

wm(2)

u(2)

−PV(wm(2))

−NV(wm(2))

…

…

…

…

wm(N−1) u(N−1) −PV(wm(N−1)) −NV(wm(N−1))

The parameter vector estimate can be obtained as, qˆ 0⫽((⌽0)T⌽0)−1(⌽0)TY0

冥

(37)

(38)

The first entry for qˆ 0 will be pˆ wd and the second Kˆ wd. dˆ +f and dˆ f can be obtained by dividing the third and fourth entries with Kˆ wd respectively,

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1499

qˆ (3) qˆ (4) dˆ +f⫽ ˆ , dˆ −f ⫽ ˆ Kwd Kwd

(39)

Noting that Kw=KtKa/J and pw=⫺B/J, and using the expressions for Kwd and pwd given in Eq. (6), the estimates for the physical inertia and viscous friction can be obtained using Eq. (25). Also, the estimated values for positive and negative Coulomb friction torque T +coul and T −coul can be obtained as, (40) Tˆ +coul⫽KtKadˆ +f, Tˆ −coul⫽KtKadˆ −f It should be noted that dˆ +f and dˆ −f (hence Tˆ +coul and T −coul) are not the final estimates for friction, but they are merely to take the bias off of Kwd and pwd (hence J and B). Prior to the identification tests, the amplifier gains Ka for the x and y axes were measured by monitoring the motor current signals as, Kax=6.4898 A/V and Kay=7.5768 A/V. Both axes are driven by Glentek GM6060-50 d.c. brush motors, for which the catalogue torque gain was given as Ktx=Kty=0.4769 Nm/A. The ball screw pitch length of both axes was hp=10 mm, resulting the in the gear ratio to be rg=10/(2π)=1.5915 mm/rad. The tachometer signals were calibrated by removing their offsets and scaling them to yield consistent velocity measurements in rad/s with the shaft encoder readings, while the encoders were calibrated to yield axis position values in mm. Identification tests were conducted on the x and y axes of the machining center which involved inputting a series of step inputs, as shown in Fig. 6. For each different test, the amplitude of the input signal was scaled with the following factors, Ku=1, 2.5, 5, 7.5, 10, 12.5, 15 and the parameters J, B, T +coul and T −coul were estimated for the cases where sufficient excitation was achieved to overcome the static friction (i.e. Kuⱖ5). Both the simple LS parameter estimation, as well as the unbiased LS technique explained above, were used to allow comparison between the two approaches. The velocity dead band was chosen as 0.5 rad/s to exceed the detected noise level when the axes were at rest. The parameter estimates obtained for different scalings of the input signal are shown in Fig. 7. Considering the figure, it is seen that the unbiased LS technique results in smaller values of inertia and damping (J and B) compared to the simple LS technique. This is because the simple LS approach does not account for the presence of Coulomb friction in its formulation, hence the effect of the real friction causes the inertia and damping estimates to appear larger than their true values. This is especially significant for the estimates of viscous damping, where the values

Fig. 6.

Input signal used in identification tests.

1500 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

Fig. 7. X and Y axis parameter estimates obtained using simple least squares and unbiased least squares identification techniques.

obtained with simple LS are around three times larger than those obtained with the unbiased algorithm. As feedforward control and Kalman filter based disturbance and friction compensation techniques rely on the drive parameters being known with a reasonable accuracy, such an error is not desirable as it will hinder further successful identification, as well as tracking control performance. It can also be observed from Fig. 7 that the unbiased LS yields much less variability between estimates obtained at different input signal scalings, especially after values of Kuⱖ7.5 V/V. Also the Coulomb friction estimates obtained with the unbiased LS algorithm are quite consistent for different values of Ku. The significantly biased values obtained at Ku = 5 V/V are due to static friction still being dominant, where neither of the parameter estimation algorithms were designed to consider its effect. However, it can be said that for all of the other input signal scalings tried, the unbiased LS algorithm yields reliable estimates of the feed drive parameters. As less significance of static friction is expected at higher control signal scalings, the parameter estimates corresponding to the case where Ku=15 V/V have been chosen. Hence, the simple LS technique yields,

冎

Jˆ LSx=7.97×10−3 kg.m2, Jˆ LSy=10.14×10−3 kg.m2 Bˆ LSx=0.0945 kg.m2/s, Bˆ LSy=0.1225 kg.m2/s

(41)

On the other hand, the estimates of the unbiased LS technique are, Jˆ x=7.65×10−3 kg.m2, Jˆ y=9.79×10−3 kg.m2 Bˆ x=0.0296 kg.m2/s, Bˆ y=0.0385 kg.m2/s Tˆ +couly=3.170 Nm Tˆ +coulx=2.877 Nm, Tˆ −coulx=−2.363 Nm,

Tˆ −couly=−3.093 Nm

冧

(42)

The performance of the two estimation algorithms have been compared by reconstructing the velocity signal using the two different models, and comparing them to the measured velocity obtained for the same input. This has been has been done for the Ku=15 V/V case, at which the

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1501

parameters were estimated. Fig. 8 shows the predicted velocity response, as well as prediction error of the two models. As can be seen, the axis model obtained with unbiased LS yields much less prediction error compared to the model obtained with simple LS. For the rest of the identification work in this article, the values of inertia and damping in Eq. (42) will be used. The Coulomb friction will be reidentified during the overall friction identification process which considers the whole shape of the Stribeck curve. 3.2. Identification of friction characteristics using a kalman filter Following the identification of the linear axis dynamics, detailed friction identification tests were conducted by jogging the axes back and forth using a pole placement controller at various speeds. The effect of the overall disturbance d in Eq. (4), was estimated using a Kalman Filter designed for the state space model in Eq. (15). Then, considering the variation of disturbance with respect to axis velocity w, the parameters of the friction model in Eq. (19) were identified. In the following, the Kalman filter and control law equations will briefly be presented, which can be found in more detail in Part III of this article. This will be followed by the friction identification results. The Kalman filter is a state observer designed to yield the most accurate state estimate possible for a dynamic system with known parameters, where the system inputs, as well as the measurements are corrupted with noise of known statistical properties. The state space model in Eq. (15) can be used to estimate the axis position x, velocity w and disturbance d simultaneously as follows,

冤 冥

xˆ (k) ˆ (k) ⫽(I⫺KobsC)A w ˆd(k)

Fig. 8.

冤 冥

冋 册

xˆ (k−1) x (k) ˆ (k−1) ⫹(I⫺KobsC)B[u(k⫺1)]⫹Kobs m w wm(k) dˆ (k−1)

(43)

Velocity response prediction of models identified through simple and unbiased least squares techniques.

1502 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

ˆ (k), and dˆ (k) are the estimates of position, velocity, and disturbance respectively. where xˆ (k), w ˆ xˆ (k⫺1), w(k⫺1), and dˆ (k⫺1) are the estimates from the previous control interval. xm(k) and wm(k) are the quantized position and noisy velocity measurements obtained through the encoder and tachometer respectively, and u(k⫺1) is the control signal applied at the previous sampling interval. The matrices A, B, and C are given in Eq. (16). Kobs is the Kalman gain matrix, designed to yield the minimum estimation error for the known A, B, C matrices, as well as the ˜ (k), respectively in noise properties, i.e. variances R−u ,Rwd,R−x , and R−w of u˜ (k), wd(k), x˜ (k), and w Eq. (15). Using the previously obtained inertia (Jx, Jy) and viscous friction (Bx, By) estimates in Eq. (42), the corresponding continuous time velocity transfer function gain Kw and pole pw are obtained as, Kw=(KtKa)/J and pw=⫺B/J. These are then used to obtain the continuous time state space representation in Eq. (5), which is discretized using Eq. (7) with the sampling period Ts=1 ms. Finally, the form in Eq. (15) is obtained by using the definitions in Eq. (16). This results in the following numerical values for the A, B, C, W, and V matrices,

冤

1 0.001589 −0.000322

冤

1 0.001590 −0.000294

冥 冤 冥 冤 0.000322

Ax⫽ 0 0.997407 −0.403843 , Bx⫽ 0.403843 , Wx⫽ 0.403843 0 0 0

1

0

0

冥 冤 冥 冤 0.000294

1

Cx⫽Cy⫽

1

冋 册 1 0 0

0 1 0

, Vx⫽Vy⫽

0

冥

0.000294 0

Ay⫽ 0 0.997563 −0.368776 , By⫽ 0.368776 , Wy⫽ 0.368776 0 0 0

冥

0.000322 0

0

1

冋 册 1 0

0 1

The DAC used has 16 bits, its output is in the range ⫺10%+10 V hence its resolution is du=20/216=0.305176×10−3 V). The variance for the difference u˜ between the commanded and actual control signal can be calculated from Eq. (10) as, Ru˜ ⫽

(du)2 ⫽7.761⫻10−9 V2 12

(44)

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1503

Similarly, for the 1024 line quadrature encoders mounted on the ball screw with 10 mm of pitch length, the resolution is dx=10/(4×1024)=2.4414×10−3 mm. Hence the variance of position measurement noise x˜ can be calculated from Eq. (12) as, (dx)2 ⫽4.9671⫻10−7 mm2 Rx˜ ⫽ 12

(45)

˜ was estimated by monitoring On the other hand, the variance of velocity measurement noise w the noise in the tachometer signals when the motors were powered, but the axes were at rest. The measured velocity noise variances are, Rw˜ x⫽7.8691⫻10−3 rad 2/s2, Rw˜ y⫽43.1189⫻10−3 rad 2/s2

(46)

It was seen that the velocity noise exhibited significant variability between different measurements, hence changing the effective values of variance. To account for this uncertainty, the x and y axis velocity noise variance values were taken as, Rw˜ x⫽Rw˜ y⫽50⫻10−3 rad 2/s2

(47)

which is higher than the observed values, but reasonably close to the actual values in designing the Kalman filter. Finally, the disturbance perturbation variance Rwd was tuned in tracking tests where the estimated disturbance dˆ would be used to cancel out the true friction. Higher values of Rwd resulted in faster estimates, better disturbance rejection, but more sensitivity to noise. The following value of Rwd resulted in both good state estimation as well as tracking performance, Rwdx⫽Rwdy⫽0.01 V2

(48)

The optimal filter gains for the above values were calculated as,

冤

0.346114

Kobsx⫽ 65.070465

0.000646 0.349254

冥 冤

−25.794230 −0.308201

0.345170

, Kobsy⫽ 64.385433

0.000640 0.331971

−27.886652 −0.310235

冥

(49)

As mentioned earlier, the algorithm used to compute the Kalman gains is presented in Part III of this article. The resulting observer poles (eig[(I⫺KobsC)A]) in the discrete-time domain (Ts=1 ms) and their corresponding natural frequencies and damping ratios are, Pobsx⫽0.6779, 0.6944⫾j0.2863, fn⫽61.88, 77.12 Hz @ z⫽1,0.59 pobsy⫽0.6833, 0.7065⫾j0.2817, fn⫽60.61, 74.46 Hz @ z⫽1,0.58

1504 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

In order to be able to jog the axes in a controlled manner, a pole placement state feedback controller was designed as,

冉冋 册 冋 册冊

u(k)⫽Kc

xref(k)

wref(k)

⫺

xˆ (k) ˆ (k) w

(50)

where xref(k) and wref(k) are the position and velocity reference values respectively, xˆ (k) and ˆ (k) are the position and velocity estimates obtained by the Kalman filter. Kc has been chosen w so that the closed loop poles (eig(A⫺BKc)) have a natural frequency of 30 [Hz] and a damping ratio of z=0.7. Hence the resulting controller gains are, Kcx⫽[48.4441 0.6065], Kcy⫽[53.0507 0.6646]

(51)

Details of the pole placement controller are also presented in Part III of this article. The reference position and velocity can both be obtained from the trajectory generator, or if only reference position values are available, the reference velocity is obtained by a simple backward differencing operation as, 1 xref(k)−xref(k−1) wref(k)⫽ · rg Ts

(52)

With the state estimator and controller ready, the friction in the drives was observed using the Kalman filter, while jogging the axes back and forth at various speeds ranging from 0.1 to 100 mm/s. More tests were conducted in the lower velocity region (i.e. ±0.1, ±0.3, ±0.5, ±1.0, ±3.0, ±5.0, ±10.0 mm/s) to better identify the boundary and partial fluid lubrication characteristics, whereas four high speed points (±25, ±50, ±75, ±100 mm/s) were sufficient in characterizing the viscous friction range. Typical test results for the x axis at 5 and 75 mm/s axis speeds are shown in Fig. 9. As can be seen, position dependency is not significant in the disturbance estimate graphs. For each test, by computing the average disturbance estimate for the positive and negative directions of motion and scaling them in terms of friction torque (Tf=KtKadf), one can obtain the

Fig. 9.

Estimated friction for along the X axis for different jogging velocities.

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1505

plots of observed friction versus axis velocity for the x and y axes, as shown with the “×” marks in Fig. 10. If the viscous friction estimate B used in the observer design had exactly matched the true viscous friction B in the system, then no viscous effect would have been existent in the observed friction values. However, it can be seen that the observed friction values at ±15.7, ±31.4, ±47.1, and ±62.8 rad/s are lined up in a slope. A positive slope, as can be seen for the x axis case, indicates that the damping estimate was smaller than the true viscous friction B, and the opposite holds for the negative slope, as is in the y axis case. By using simple linear regression, it is possible to estimate the slopes of these “high speed” points as, ⌬B+x =0.0025 kg.m2/s,

⌬B−x =0.0024 kg.m2/s

⌬B+y =−0.0086 kg.m2/s, ⌬B−y =−0.0111 kg.m2/s

冎

(53)

where the superscripts “+” and “⫺” denote positive and negative velocities respectively. It is not normally expected for the lubricant to exhibit different viscosity values for different directions of motion. Hence, the positive and negative direction slopes can be combined into a single value by averaging them,

冎

⌬B¯ x=(⌬B+x +⌬B−x )/2=0.0025 kg.m2/s ⌬B¯ y=(⌬B+y +⌬B−y )/2=−0.0098 kg.m2/s

(54)

The observed extra viscous friction ⌬B can be lumped back into the previous damping estimates B as, Bˆ ⬘x=Bˆ x+⌬B¯ x=0.0321 kg.m2/s Bˆ ⬘y=Bˆ y+⌬B¯ y=0.0287 kg.m2/s

冎

(55)

This adjustment to the viscous friction estimates is shown in Fig. 10 as “initially identified B” corresponding to B, and “corrected B” to B⬘. By removing the extra viscous effect ⌬B·w from the observed friction points, which is now contained in the adjusted damping value B⬘, it is possible to obtain the corrected friction values,

Fig. 10.

Observed and modeled friction characteristics for the X and Y axes.

1506 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

shown with the dots in Fig. 10. The Coulomb friction values for the positive and negative directions of motion are found by averaging the values of corrected friction, that correspond to the speeds of ±15.7, ±31.4, ±47.1, and ±62.8 rad/s with the correct signs. Which result in, Tˆ +coulx=2.2054 Nm, Tˆ −coulx=−1.6562 Nm Tˆ +couly=3.0578 Nm, Tˆ −couly=−2.7386 Nm

冎

(56)

The static friction values are taken as the friction values recorded at the smallest jogging speed (i.e. 0.063 rad/s),

冎

Tˆ +statx=3.647 Nm, Tˆ −statx=−2.7021 Nm Tˆ −staty=3.6523 Nm, Tˆ −staty=−3.2183 Nm

(57)

The analytical model shown in Eq. (19) was used to describe the friction characteristics of the x and y axes as, +

+

−

−

冎

T +s(w)=T +state−w/⍀1+T +coul(1−e−w/⍀2 )+T +viscw

(58)

T −s (w)=T −state−w/⍀1+T −coul(1−e−w/⍀2)+T −viscw

The Coulomb and static friction values have already been determined in Eqs. (56) and (57). The viscous friction coefficients Tvisc+ and Tvisc− will be equal to the corrected viscous friction values, Tˆ +viscx⫽Tˆ −viscx⫽Bˆ ⬘x, Tˆ +viscy⫽Tˆ −viscy⫽Bˆ ⬘y

(59)

The velocity constants ⍀+1, ⍀−1 , ⍀+2, and ⍀−2 are estimated by minimizing the cost function,

冘 N2

1 [dT(k)]2 V⫽ 2k⫽N

(60)

1

that penalizes the error dT(k) between predicted and measured friction values, +/−

+/−

−w(k)/⍀1 −w(k)/⍀2 ⫹T +/− )] dT(k)⫽dT[w(k)]⫽T⬘f[w(k)]⫺[T +/− stat·e coul·(1⫺e

(61)

where, w(k) is the array of the jogging velocities used in the tests (i.e. w(k)=±0.063, ±0.189, ±0.314,…, ±47.124, ±62.832 rad/s for k=1, 2, …, 11). T⬘f[w(k)] is the corrected friction torque observed at axis the velocity w(k) (i.e. estimated average minus the viscous term). ⍀1, and ⍀2, are the velocity constants to be estimated, and the “+/⫺” superscript denotes the relevant direction of motion. Noting that Eq. (61) cannot be expressed as a linear combination ⍀1 and ⍀2, these parameters were obtained by testing a range of values from 0 to 10 rad/s (⫺10 rad/s for the negative velocity case) with increments of 0.1 rad/s, and choosing the ones that resulted in the lowest value of the cost function in Eq. (60). The points to be included into the

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1507

cost function evaluation were chosen as N1=4 and N2=11 (i.e. the range corresponding to velocities between 0.628 and 62.832 rad/s). The parameters obtained at the end of this optimization were, ˆ −1x=−2.0 rad/s ˆ +1x=2.0 rad/s, ⍀ ⍀ ˆ +2x=2.0 rad/s, ⍀ ˆ −2x=−2.0 rad/s ⍀ ˆ −1y=−3.6 rad/s ˆ +1y=4.1 rad/s, ⍀ ⍀ ˆ −2y=−4.2 rad/s ˆ +2y=5.4 rad/s, ⍀ ⍀

冧

(62)

The resulting friction model, from which the viscous component has been removed, is shown in Fig. 10 with the label “Model”. As can be seen, the analytical expression is sufficiently close in representing the experimentally observed friction values. The identified feed drive parameters, as well as the design values of the Kalman filter and the controller have been summarized in Table 1.

Table 1 Summary of identified feed drive and controller design parameters Identified/designed parameters Feed drive Amplifier Ka (A/V) Motor Kt (Nm/A) Inertia Jˆ (kg.m2) Viscous friction Bˆ (kg.m2/s) Lead screw rg (mm/rad) Noise and quantization DAC resolution du (V) Encoder resolution dx (µm) Tach. noise var. Rw (rad 2/s2) Friction Static friction T +/− stat (Nm) Coluomb friction T +/− coul (Nm) (rad/s) Velocity constant 1 ⍀+/− 1 (rad/s) Velocity constant 2 ⍀+/− 2 Closed loop observer (Kalman filter) poles Natural frequency f 1,2 0 (Hz) Damping ratio z1,2 0 Pole placement controller closed loop poles (Hz) Natural frequency f 1,2 c Damping ratio z1,2 c Discrete time sampling period Ts (s)

X axis 6.4898 0.4769 7.65×10⫺3 0.0321 1.5915

Y axis 7.5768 0.4769 9.79×10⫺3 0.0287 1.5915

3.0518×10⫺4 2.441 50×10⫺3

3.0518×10⫺4 2.441 50×10⫺3

3.647, ⫺2.702 2.205, ⫺1.656 2.0, ⫺2.0 2.0, ⫺2.0

3.652, ⫺3.218 3.058, ⫺2.739 4.1, ⫺3.6 5.4, ⫺4.2

61.88, 77.12 1, 0.59

60.61, 74.46 1, 0.58

30 0.7

30 0.7 0.001

1508 K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509

4. Experimental verification of the friction model in contouring tests As an experimental verification of the identified friction model, two circular contouring tests were conducted using the pole placement control scheme with the Kalman filter for state estimation. In both tests, a zero phase error tracking controller (ZPETC [3]) was used in feedforward to widen the overall axis tracking bandwidth. This resulted in the tracking errors to be mainly dependent on friction and parameter mismatch, and not the limited feedback control bandwidth. In the first test, only the pole placement control law with feedforward action was used. In the second test, the friction torque was compensated in feedforward using the identified model. In neither case, was on-line disturbance cancellation used, leaving only feedforward friction compensation to cope with the adverse effect of friction on the tracking accuracy. A circular reference trajectory with a radius of 35 mm was traveled at a feedrate of 170 mm/s with a tangential acceleration of 2000 mm/s2 and a jerk limit of 50 000 mm/s3. The reference trajectory generation and axis controller design details can be found in Parts I and III of this article. The tracking and contour error profiles for the two tests are shown in Fig. 11. Considering Fig. 11(a), where no feedforward friction compensation has been used, the axis tracking errors reach a value of 20 µm, and there is steady state error before and after following the circle (i.e. prior to 0.4 s and after 1.8 s) due to the feedback controller not having integral action. The contour error is also significant, and undergoes sudden changes as the quadrants of the circle are crossed. On the other hand in Fig. 11(b), it is seen that by injecting the predicted friction values into the control signal, both the tracking and contour errors are reduced significantly, to a maximum value of about 12 µm. Furthermore, there is now no steady state error, in spite of the absence of integral action. In fact, the results in Fig. 11(b) are quite comparable with those obtained in Part III of this article, where on-line disturbance cancellation is used on top of feedforward friction compensation to improve the robustness against parameter variations and unmodeled disturbances. These results demonstrate that the identified friction model can be used to successfully represent the actual friction, and counteract its degrading effect on the tracking and contouring accuracy.

Fig. 11. Recorded tracking and contouring errors for control (a) without and (b) with feedforward friction compensation.

K. Erkorkmaz, Y. Altintas / International Journal of Machine Tools & Manufacture 41 (2001) 1487–1509 1509

5. Conclusions This paper has presented the details of a necessary initial step in high speed CNC design, which is the modeling and identification of the feed drives’ dynamics. A linear servo model has been combined with a practical friction model, and the joint formulation has allowed unbiased estimation of the axis parameters like inertia and viscous friction. Accurate knowledge of these parameters enables successful state estimation, desired closed loop response realization, and feedforward control that cancels out the correct dynamics. The parameters of the friction model have been identified using a Kalman filter, and verified in high speed contouring trials where the modeled friction has been used to immediately counteract the actual friction, in order to prevent it from degrading the tracking and contouring performance. Hence, the axis dynamics and friction characteristics required to design the high performance CNC controller in Part III of this article have been accurately modeled and identified. References [1] Y. Koren, Computer Control of Manufacturing Systems, McGraw-Hill, New York, 1983. [2] G. Pritschow, W. Philipp, Research on the efficiency of feedforward controllers in direct drives, Annals of CIRP 41 (1) (1992) 411–415. [3] M. Tomizuka, Zero phase error tracking algorithm for digital control, Journal of Dynamic Systems, Measurement, and Control 109 (1987) 65–68. [4] M. Weck, G. Ye, Sharp corner tracking using the IKF control strategy, Annals of CIRP 39 (1) (1990) 437–441. [5] E.D. Tung, M. Tomizuka, Feedforward tracking controller design based on the identification of low frequency dynamics, Journal of Dynamic Systems, Measurement and Control 115 (3) (1993) 348–356. [6] H.B. Armstrong, P. Dupont, A. Canudas De Wit, Survey of models, analysis tools and compensation methods for the control of machines with friction, Automatica 30 (7) (1994) 1083–1138. [7] H.S. Lee, M. Tomizuka, Robust motion controller design for high-accuracy positioning systems, IEEE Transactions on Industrial Electronics 43 (1) (1996) 48–55. [8] K. Erkorkmaz, High speed contouring control for machine tool drives, M.A.Sc. thesis, The University of British Columbia, Department of Mechanical Engineering, Vancouver, 1999. [9] Y. Altintas, Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design, Cambridge University Press, Cambridge, 2000. [10] L. Ljung, System Identification: Theory for the User, 2nd ed., Prentice-Hall of Canada, Canada, 1998. [11] R.E. Kalman, A new approach to linear filtering and prediction problems, Journal of Basic Engineering 82 (1960) 35–44. [12] K.J. Astrom, B. Wittenmark, Computer-Controlled Systems: Theory and Design, 3rd ed., Prentice-Hall, Englewood Cliffs, NJ, 1997.