Pre-compensation of contour errors in five-axis CNC machine tools

International Journal of Machine Tools & Manufacture 74 (2013) 1–11

Contents lists available at ScienceDirect

International Journal of Machine Tools & Manufacture journal homepage: www.elsevier.com/locate/ijmactool

Pre-compensation of contour errors in five-axis CNC machine tools Ke Zhang a, Alexander Yuen b, Yusuf Altintas b,n a b

Key Laboratory of Mathematics Mechanization, Chinese Academy of Sciences, Beijing 100190, China Manufacturing Automation Laboratory, The University of British Columbia, Vancouver, BC, Canada V6T1Z4

art ic l e i nf o

a b s t r a c t

Article history: Received 18 February 2013 Received in revised form 8 July 2013 Accepted 11 July 2013 Available online 20 July 2013

This paper presents an analytical prediction and compensation of contouring errors in five-axis machining of splined tool paths. The position commands are first fitted to piecewise quintic splines while respecting velocity, acceleration and jerk continuity at the spline joints. The transfer function of each servo drive is kept linear by compensating the disturbance effect of friction with a feed-forward block. Using the analytically represented five-axis, splined tool path, splined tracking errors and kinematic model of the five-axis machine tool, contouring errors are predicted ahead of axis control loops. The contouring errors are decoupled into three linear and two rotary drives, and the position commands are modified before they are sent to servo drives for execution. The proposed method has been experimentally demonstrated to show significant improvement in the accuracy of contouring fiveaxis tool paths. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Five-axis Tracking error Contour error Compensation

1. Introduction Five-axis CNC machining centers are widely used in industry in machining dies, molds and aerospace parts with sculptured surfaces. While low feeds and speeds are used in roughing operations, high feeds and spindle speeds are needed in finish machining of these parts with ball end mills. Although the cutting loads are quite low in finish machining, the feeds, i.e. the material removal rates, are constrained mainly by the contouring accuracy of the feed drives and controller. The drives have a limited frequency bandwidth, and they lag behind the position commands causing tracking errors which increase proportional to the axis feed speeds. When three translational and two rotary drives of a five-axis machine try to follow a curved path, the tracking errors of the drives lead to significant deviations from the commanded path which is described as contouring errors of the machine. The programmed feeds must be kept at conservative levels in order to avoid having contour errors that are larger than the tolerance of the part. The minimization of contouring errors has been studied in several directions in the past [1]. Machine tool designers install stronger motors with high torque constants to increase the bandwidth of the drives, but at the expense of higher cost and larger mass on the machine. Significant research has been directed towards improving the bandwidth of the drives by developing advanced axis control laws. Tomizuka [2] developed Zero Phase Error Tracking Control (ZPETC)

n

Corresponding author. Tel.: +1 604 822 5622; fax: +1 604 822 2403. E-mail addresses: [email protected], [email protected] (Y. Altintas).

0890-6955/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijmachtools.2013.07.003

system which cancels the stable poles of the drive's closed loop transfer function. ZPETC minimizes the phase lag between the reference position input and actual output of the drive, which yields reduced tracking errors. However, the algorithm is sensitive to the changes in the transfer function of the drive that occurs during positioning of the machine within the work envelope. Weck and Ye [3] removed the high frequency content of the commands generated by ZPETC, but at the expense of phase distortion. Renton and Elbestawi [4] developed a servo control law that uses the axis performance envelope as well as instantaneous position, velocity, and acceleration information of the target path in the presence of disturbances. Altintas et al. [5] presented an adaptive sliding mode control technique for high speed feed drives. The proposed control system is robust against uncertainties in the drive's parameters, maximizes the bandwidth within physical limitations, and compensates for external disturbances such as friction and cutting force. Erkorkmaz and Altintas [6] provided a systematic approach for designing a control law which provides a high tracking bandwidth as well as adequate disturbance rejection and parameter variation robustness. Unless it is actively damped, the bandwidth of the drive is still limited by the first natural frequency of the feed drive structure regardless of the control law used in rigid body mode. Kamalzadeh et al. [7] and Okwudire et al. [8] expanded the sliding mode controller to reject the external disturbances while actively damping the torsional, lateral and axial vibrations of the ball screw drives. The active damping expands the bandwidth and reduces the tracking errors of drives. Alternative methods are based on prediction and compensation of contouring errors. Koren [9] proposed a cross-coupled controller (CCC)

2

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

S(t) splined axis position e(t) predicted tracking error using splined axis position e ¼ ½ex ; ey ; ez ; ea ; ec T five-axis tracking error vector P ref , Oref reference tool tip position and orientation P act , Oact actual tool tip position and orientation εp , εo tool tip and orientation contour errors P comp , Ocomp pre-compensated tool tip position and orientation qcomp pre-compensated five-axis position command

Nomenclature P ¼ ½P x ; P y ; P z T tool tip position O ¼ ½Oi ; Oj ; Ok T tool orientation q ¼ ½X; Y; Z; A; CT five-axis position command Ts sampling period Ns number of interpolation steps tf total path travel time qref ¼ ½X ref ; …; C ref T reference five-axis position command

to minimize the contouring errors in biaxial machines. Later, Koren and Lo [10] introduced a variable gain CCC to adjust the controller as the path curvature changes. However, CCC acts in real time at the axis level, and the controller cannot act fast enough to avoid large contouring errors at sharp curvatures. Altintas and Sencer [11] proposed a sliding mode controller that predicts and compensates the contouring errors along five-axis tool path. Instead of controlling individual axis, their multi input–multi output controller controls all five drives to minimize the contouring error. Ernesto and Farouki [12] compensated the inertia and viscous damping of the drives by a priori modifications of the commanded path geometry. Altintas and Khoshdarregi [13] presented a vibration avoidance and contouring error compensation algorithm for feed drives by shaping the trajectory ahead of axis control commands. The successful pre-compensation of tracking and contouring errors require accurate knowledge of servo drive dynamics and their effects on the contouring errors along the tool path ahead of machine motion. This paper presents an analytical prediction of tracking errors for three translational and two rotary drives of fiveaxis machines when they are used in contour machining of sculptured surfaces. The tool tip position and tool orientation errors are predicted by combining a five-axis splined tool path and the analytically expressed tracking errors of drives. The trajectory commands are modified in the feed-forward blocks to compensate the predicted contouring errors before sending them to the servo drives. The proposed method is experimentally validated on a fiveaxis CNC research machine. The rest of this paper is organized as follows. Section 2 gives the tracking error prediction model. Section 3 gives the contour error pre-compensation method. Section 4 gives the simulation and experimental results, and the paper is concluded in Section 5.

workpiece coordinate (P-system) to the machine tool coordinates (M-system) shown in Fig. 1 is carried out using the following inverse kinematics transformation: 8 A ¼ arccosðOk Þ > > > > > C ¼ arctanðOi =Oj Þ > < X ¼  cos ðCÞP x  sin ðCÞP y > > > > Y ¼ cos ðAÞ sin ðCÞP x  cos ðAÞ cos ðCÞP y  sin ðAÞP z  sin ðAÞLac;z > > : Z ¼ sin ðAÞ sin ðCÞP  sin ðAÞ cos ðCÞP þ cos ðAÞP þ cos ðAÞL x

y

z

ð4Þ where Lac;z and LTya;z are offsets of the rotary table. Conversely, the transformation from the M-system to the P-system is carried out using the direct kinematics transformation: 8 Oi ¼ sin ðAÞ sin ðCÞ > > > > > O > j ¼ sin ðAÞ cos ðCÞ > > > < Ok ¼ cos ðAÞ P x ¼  cos ðCÞX þ sin ðCÞ cos ðAÞY þ sin ðAÞ sin ðCÞZ sin ðCÞ sin ðAÞLTya;z > > > > > P y ¼  sin ðCÞX cos ðCÞ cos ðAÞY sin ðAÞ cos ðCÞZ þ cos ðCÞ sin ðAÞLTya;z > > > > : P z ¼  sin ðAÞY þ cos ðAÞZ cos ðAÞLTya;z Lac;z

ð5Þ The discrete axis position commands at each control interval (Ts) are expressed, i.e. for x-axis, as X ref ¼ ½Xð0Þ; Xð1Þ; Xð2Þ; …; XðN s Þ;

2.1. Fitting a quintic spline to axis position commands

ð1Þ

and the unit tool orientation vector is given by ∥O∥ ¼ 1:

ð2Þ

Using the inverse kinematics transformation shown in Eq. (4), the reference tool path variables ðP x ; P y ; P z ; Oi ; Oj ; Ok Þ are transformed into position commands of the physical feed-drives of the machine tool as follows: qðtÞ ¼ ½XðtÞ; YðtÞ; ZðtÞ; AðtÞ; CðtÞT ;

ð6Þ

where Ns is the number of discrete steps and t f ¼ N s T s is the total path travel time. The objective of this section is to predict the tracking error of each axis analytically as described in Fig. 2.

The tool tip position at time (t) is expressed in the workpiece coordinate system (P-system) by the vector (Fig. 1) as follows:

OðtÞ ¼ ½Oi ðtÞ; Oj ðtÞ; Ok ðtÞT ;

þ LTya;z

Since discontinuities in the higher derivatives of the trajectory can result in an undesired vibration response from the machine

2. Tracking error prediction

PðtÞ ¼ ½P x ðtÞ; P y ðtÞ; P z ðtÞT ;

ac;z

ð3Þ

where ðX; Y; Z; A; CÞ are the axis coordinates in the machine tool coordinate system (M-system). The transformation from the

Fig. 1. Five-axis machine tool configuration.

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

Input: axis position commands

Spline fitting data:

Fit

quintic spline to the data with minimum jerk

Within the tolerance

?

Add the point with maximum error to the spline fitting data

No

Yes

3

the transition points of spline segments (see Appendix A for the details). The fitting errors (Fig. 3) from the reference positions X ref to spline S(t) are calculated in each iteration for each axis. If the maximum error is less than the tolerance ε, the fitting process is completed. Otherwise, as shown in Fig. 2, the reference point which has the maximum error is inserted to the spline fitting points, and the process is repeated until the errors are within tolerance ε. Finally, a C3 quintic spline for each of the five axes is obtained with n sub-spline segments (Eq. (7)).

quintic spline

Solve the differential equation for tracking error

Output: analytical tracking error

2.2. Analytical solution of tracking error

Fig. 2. Tracking error prediction model.

Fig. 3. Axis position fitting.

tool, the jerk continuity of the pre-compensated trajectory is preserved with a quintic spline [14,15]. This will result in a jerk continuous solution for tracking error, and subsequently a jerk continuous solution for contour error which will lead to smoother trajectory. The axis position commands within a tool path segment are first fitted to a C3 quintic spline iteratively within a tolerance ε, i.e., within an encoder resolution. The time stamped position commands (Eq. (6)) are subdivided into n coarse segments but with equal time length as shown in Fig. 3, where each segment may initially contains, i.e. 50 reference position commands (n ¼ ⌊N s =50⌋). The objective is to describe the position commands by a small number of piecewise, analytical splines, without having to consider all points in the reference trajectory (Eq. (6)). If n þ 1 equally time spaced x-axis positions (x0 ; x1 ; …; xn ) are initially selected from Xref (Eq. (6)), the quintic spline for segment (i) can be expressed piecewise as Si ðtÞ ¼ b5;i t 5 þ b4;i t 4 þ b3;i t 3 þ b2;i t 2 þ b1;i t þ b0;i ; S∈fSX ; SY ; SZ ; SA ; SC g

A five-axis machine tool with three linear and two rotary drives controlled by PID (Proportional Derivative Integral) controllers are considered as shown in Fig. 4. For the simplicity of the tracking error prediction model, only the rigid body machine dynamics is considered, and the cutting force disturbances are neglected since ball screw drives reduce their effect significantly. The Coulomb friction in the drives is identified and compensated in the feedforward block during experiments. In Fig. 4, the tracking error e ¼ Xx is the difference between the commanded and actual axis positions. The current amplifier has a gain of ka and the torque constant of the motor is kt : The equivalent inertia and viscous damping reflected at the motor shaft are J and B, respectively. The ball screw transmission gain from angular to linear motion is r g : The transfer function between the actual position of the drive (x) and the reference input (X) can be expressed in the Laplace domain as follows (Fig. 4): Kðkd s2 þ kp s þ ki Þ xðsÞ ¼ 3 : XðsÞ Js þ ðB þ Kkd Þs2 þ Kkp s þ Kki

ð9Þ

where K ¼ ka kt r g : The three poles of the transfer function (Eq. (9)) are denoted as λ1 , λ2 , and λ3 respectively, which are assumed to be all distinct. The transfer function between the tracking error (eðsÞ ¼ XðsÞxðsÞ) can be evaluated from Eq. (9) as follows: eðsÞ Js3 þ Bs2 ¼ 3 : XðsÞ Js þ ðB þ Kkd Þs2 þ Kkp s þ Kki

ð10Þ

which can be expressed as an equivalent differential equation in time domain as follows: 



_ þ Kki eðtÞ ¼ JX ðtÞ þ BX€ ðtÞ: € þ Kkp eðtÞ J e ðtÞ þ ðB þ Kkd ÞeðtÞ

t∈½0; hi ;

ð11Þ

ð7Þ

where time interval in spline segment i is hi ¼ t i t i1 , i ¼ 1; …; n and 0 ¼ t 0 o t 1 o … o t n ¼ t f . The coefficients in Eq. (7) are to be determined for each segment hi for n number of sub-spline segments. Once the spline coefficients are determined using the method explained below, the error between each reference position point and its corresponding value predicted by the spline is checked. The point with the highest error in the entire trajectory (Eq. (6)) is added to the spline fitting data, and the process is iteratively continued until all points are predicted by the spline within the tolerance ε. Depending on the curvature, a trajectory with 2000 points may be represented by a spline fitted to only 80∼100 points in the set. The spline parameters are identified by minimizing the integral of square of axis jerk as proposed in [14] Z tf  J¼ S ðtÞ2 dt ð8Þ

The linear differential equation (Eq. (11)) with constant coefficients has homogenous and particular solutions [16]. Since the reference positions X(t) are expressed by piecewise quintic polynomials as shown in Eq. (7), the right hand side of Eq. (11) is a cubic polynomial which yields to a particular solution in the following form: e^ i ðtÞ ¼ d3;i t 3 þ d2;i t 2 þ d1;i t þ d0;i

at each time segment ½0; hi , i ¼ 1; …; n. Substituting the undetermined solution (Eq. (12)) into Eq. (11), and equating the coefficients of terms 1, t, t 2 , t 3 on each side of the equation, the

0

while respecting the continuity of jerk, acceleration and velocity at

ð12Þ

Fig. 4. Block diagram for the PID controller and x-axis drive.

4

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

following needs to hold: 3 2 32 3 2 d0;i 2b2;i B þ 6b3;i J 6J Kki Kkp 2ðB þ Kkd Þ 7 6 76 7 6 2Kkp 6ðB þ Kkd Þ 76 d1;i 7 6 6b3;i B þ 24b4;i J 7 Kki 6 0 7; 6 76 7¼6 6 7 6 0 7 6 7 3Kkp 0 Kki 4 54 d2;i 5 4 12b4;i B þ 60b5;i J 5 d3;i 20b5;i B 0 0 0 Kki ð13Þ which yields to the solution of coefficients ðd0;i ; d1;i ; d2;i ; d3;i Þ. By combining homogenous (μ1;i eλ1 t þ μ2;i eλ2 t þ μ3;i eλ3 t ) and particular solutions, the general solution of Eq. (11) can be expressed at each time interval ½0; hi  as follows: ei ðtÞ ¼ e^ i ðtÞ þ μ1;i eλ1 t þ μ2;i eλ2 t þ μ3;i eλ3 t ;

t∈½0; hi ;

ð14Þ

where μ1;i , μ2;i and μ3;i are undetermined integration constants on ith time interval. For preserving continuity up to the second derivative of e(t), each two consecutive segments ei(t) and eiþ1 ðtÞ need to have the same value at the junctions, as well as their first and second derivatives. Consequently, the following equations should hold for i ¼ 1; …; n1: 8 > < ei ðhi Þ ¼ eiþ1 ð0Þ; e_ i ðhi Þ ¼ e_ iþ1 ð0Þ; ð15Þ > : e€ ðh Þ ¼ e€ ð0Þ: i i iþ1 Substituting Eq. (14) into Eq. 3 2 2 3 2 μ1;i eλ1 hi μ1;iþ1 1 6μ 7 6 6λ λ2 hi 7 6 7 μ e þ ¼ 4 2;iþ1 5 4 2;i 5 4 1 λ3 hi μ3;iþ1 λ21 μ3;i e

(15) 1 λ2 λ22

3 31 2 ^ e i ðhi Þe^ iþ1 ð0Þ 6 7 λ3 5 4 e^_ ðh Þe^_ ð0Þ 7 5; 1

i

λ23

i

iþ1

ð16Þ

e^€ i ðhi Þe^€ iþ1 ð0Þ

which holds for the segments i ¼ 1; …; n1. With the initial _ € conditions eð0Þ ¼ eð0Þ ¼ eð0Þ ¼ 0 on the first time interval ½0; h1 , Eq. (14) yields the initial integration constants as follows: 2 3 2 3 2 3 μ1;1 1 1 1 1 e^ 1 ð0Þ 6μ 7 6 7 6 7 ð17Þ 4 2;1 5 ¼ 4 λ1 λ2 λ3 5 4 e^_ 1 ð0Þ 5: 2 2 2 μ3;1 λ1 λ2 λ3 e^€ ð0Þ 1

The integration constants μ1;i , μ2;i and μ3;i at the remaining segments can be obtained recursively from Eqs. (16) and (17), which give the analytical expression of tracking errors e(t) for each axis (Eq. (14)). Note that the right hand side of Eq. (11) is continuous since X(t) is C3, which implies that the left hand side of Eq. (11) is also continuous. Since the continuity up to the second derivative of e(t) has been preserved, the predicted tracking error function e(t) is also C3. The tracking error prediction model proposed here can be generalized to deal with any other linear systems, such as using P– PI controllers.

3. Pre-compensation of contouring errors Due to tracking errors of individual servo drives, two types of contour errors occur in five-axis motion control systems [11]. The first one is the normal deviation of the tool tip from the desired tool path, called the tool tip contour error. The second one is defined as deviation from the desired tool orientation in spherical coordinates and controlled only by the rotary axes, which is called the orientation contour error. 3.1. Prediction of contour errors The contouring error has two components: the tool tip contour error εp and tool orientation error εo : The geometric reference tool tip positions and tool orientations are represented by C3 splines PðuÞ; 0≤u≤1 and OðwÞ; 0≤w≤1 respectively in this paper, which are generated by using the method proposed in [17]. Since the axis tracking errors are modeled analytically in the previous section (Eq. (14)), it is possible to predict the contouring errors in five-axis machining. Let the reference and actual tool tip positions are denoted as P ref and P act ; respectively, as shown in Fig. 5(a). P n is the point on the reference but splined tool path PðuÞ which is closest to the actual tool tip position (P act Þ: The corresponding tool tip contour error vector is defined as P n P act , which should be perpendicular to the tangent vector P′ðuÞ ¼ ðd=duÞPðuÞ of the reference tool path at P n . The point P n can be found by solving the following equation numerically: P′ðuÞ  ðPðuÞP act Þ ¼ 0;

ð18Þ

u∈½0; 1;

which yields the corresponding spline parameter un , i.e. P n ¼ Pðun Þ. The tool tip contour error is then evaluated as follows: εp ¼ ∥Pðun ÞP act ∥:

ð19Þ

Similarly, the reference and actual tool orientations are denoted as Oref and Oact , respectively, as shown in Fig. 5(b). Let On is the closest point from Oact to the reference orientation contour OðwÞ in spherical _ coordinates. The orientation contour error is the arc length of Oact On , i. e. the angular change from Oact to On , which is arccos ðOact  On Þ. The plane determined by Oact , On , and the origin O should be perpendicular to the tangent vector O′ðwÞ ¼ ðd=dwÞOðwÞ of the orientation contour at On . The equivalent condition is that, Oact is perpendicular to the tangent vector O′ðwÞ at On , hence the point On can be found by solving O′ðwÞ  Oact ¼ 0;

ð20Þ

w∈½0; 1; n

which yields the corresponding orientation spline parameter w , i.e. On ¼ Oðwn Þ. The tool orientation contour error becomes εo ¼ arccosðOact  Oðwn ÞÞ:

Fig. 5. Tool tip position and orientation contour errors.

ð21Þ

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

3.2. Compensation of contour errors It can be observed from the inverse and direct kinematics transformations shown in Eqs. (4) and (5) that the translations are coupled with rotations in five-axis machining. However, the tool orientation vector is not affected from the changes in translational motions (Eq. (5)) [18], hence the orientation contour error is precompensated ahead of compensating tool tip errors. The proposed strategy in compensating tool tip and tool orientation errors is outlined in Fig. 6, and detailed as follows. Tracking error prediction

Direct kinematics

Tracking error prediction

Inverse kinematics

Orientation contour error compensation

Inverse kinematics Tool tip contour error compensation

Direct kinematics

5

The tracking errors of all five drives have been analytically modeled in Section 2 from the closed loop transfer functions as follows: eðtÞ ¼ ½ex ðtÞ; ey ðtÞ; ez ðtÞ; ea ðtÞ; ec ðtÞT :

ð22Þ

The actual positions (q1 ) of the drives can be predicted as follows: q1 ¼ qref e←qref ¼ ½X ref ; Y ref ; Z ref ; Aref ; C ref T

ð23Þ

where qref is the reference five-axis position command vector. Using the direct kinematics transformation (Eq. (5)) for q1 , the predicted tool orientation is obtained and denoted as O1 . In Fig. 5 (b), Oref is the reference tool orientation; On is calculated from Eq. (20); Oact ¼ O1 ; and Ocomp is the pre-compensated tool orientation that needs to be determined using a “spherical parallelogram”. As shown in Fig. 5(b), Oref , Oact , On , and Ocomp are four vertexes of the spherical parallelogram, whose center is denoted _ as Oc . The_center is defined as the midpoint of the circle arc Oref On (and Ocomp Oact ) and obtained from the following equation: Oc ¼

Oref þ On Ocomp þ Oact ¼ : ∥Oref þ On ∥ ∥Ocomp þ Oact ∥

ð24Þ

As shown in Fig. 5(c), the length of the projection of Oact on Oc is ∥Ocomp þ Oact ∥=2, which can also be obtained from inner product Oact  Oc since Oact and Oc are both unit vectors. The relationships yield to

Fig. 6. Contouring errors pre-compensation strategy.

Ocomp ¼ ∥Ocomp þ Oact ∥Oc Oact ¼ 2ðOact  Oc ÞOc Oact :

10

0

20 Px [mm]

40

−40

−60

60

0

−20

Py [mm]

Fig. 7. Five-axis tool path used in the experiments.

1

Ok

Pz (mm)

0.5 0 0

0 −0.5

40

−20

P

y

20 −40

(m

m)

−1 −1

)

mm

0

Px (

−20

−1

0

0

Oj

1

Oi

1

300 200 100 0 −100

Axis position (rad)

2 X−axis Y−axis Z−axis

400

A−axis C−axis

1.5 1 0.5 0 −0.5

0

0.5

1

1.5

Time (s)

2

ð26Þ

Since the five-axis positions will change due to coupled kinematics after modifying the tool orientations, the tracking error prediction needs to be updated after finding tool tip contour error pre-compensation. Using the inverse kinematics transformation Eq. (4) with the reference tool tip position P ref and the precompensated tool orientation Ocomp , the new five-axis reference positions q2 are obtained. Again, taking q2 as the input of the

0 −20

ð25Þ

From Eqs. (24) and (25), Ocomp is obtained as follows:

Cutting tool

20

Axis position (mm)

Pz [mm]

∥Ocomp þ Oact ∥ ¼ 2ðOact  Oc Þ:

2.5

0

0.5

1

1.5

2

2.5

Time (s)

Fig. 8. Test tool path splines and axis position commands. (a) Tool tip position spline, (b) Tool orientation spline, (c) Translational axes and (d) Rotary axis.

6

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

tracking error prediction algorithm, the actual five-axis positions are predicted as q3 by considering the transfer functions of the drives (Eq. (9)). Applying the direct kinematics transformation (Eq. (5)) on q3 , the re-predicted tool tip position is obtained as P 3 . In Fig. 5(a), P ref is the reference tool tip position; P n is calculated from Eq. (18); and P act ¼ P 3 . The final pre-compensated tool tip position P comp is determined by adding the contour error vector P n P act to the reference P ref as follows: P comp ¼ P ref þ P n P act :

ð27Þ

Finally, by applying the inverse kinematics transformation (Eq. (4)) on the pre-compensated tool tip position P comp and tool orientation Ocomp , the pre-compensated five-axis position command qcomp is calculated. Theoretically, using the modified fiveaxis position command will lead the machine to follow the desired trajectory more closely since the effects of servo drive tracking errors on the path contouring are pre-compensated by modifying the original drive commands.

4. Simulation and experimental results A five-axis tool path shown in Fig. 7 is adopted in the experiments. The tool tip position and orientation splines fitted to five-axis tool path are shown in Fig. 8(a) and (b), respectively. The time stamped position commands are also generated as shown in Fig. 8(c) and (d). The sampling period of 1 ms with the maximum feedrate of 100 mm/s has been used in the experiments Workstation PC - Read NC-Program - Generate Trajectories - Uncompensated - Pre-compensated

conducted on a five-axis experimental machining center controlled by an in-house developed, modular and reconfigurable open CNC system running on dSPACE signal processing board, see Fig. 9. The drive parameters for each axis are identified and listed in Table 1. To validate its effectiveness in linear controllers, two PID controllers of varying bandwidth are used in the simulations and experiments. The PID controller parameters (Table 1) are tuned to match the dynamics of all drives. The friction as well as the equivalent inertia of each drive have been identified as outlined in [19]. The disturbance torque at each corresponding velocity is used to develop a friction model of the following form: T^ f



8 þ þ ω=Ωþ ω=Ωþ 1 þ T 2 Þ > < T stat e coul ð1e   ω=Ω1 ω=Ω 2 Þ ωÞ ¼ T  þ T coul ð1e stat e > : 0

if ω 4 0 ð28Þ

if ω o 0 otherwise

The values of the parameters in Eq. (28) are listed in Table 2. The model compensates the static friction at low velocities, while compensating the Coulomb friction at higher velocities. The viscous friction (B) is included in the linear, transfer function of the rigid drive. Since the proposed tracking error prediction model assumes a plant with only linear characteristics, the non-linear friction disturbances have been compensated by a feed-forward block as shown in Fig. 10. Typically, increasing the bandwidth of individual axis tracking control will reduce the contour errors as well. However, given the limited bandwidth of tracking control, contour errors will still exist unless the trajectory is pre-compensated. This can be seen in Fig. 11, where the simulated contouring results of a high performance, aggressive sliding mode controller [5] is compared to the simulated contouring results of a P–PI controller where both controllers had feed-forward friction control. Although the sliding mode controller has higher axis control bandwidth, it gives

PCI Bus (Reference Commands)

dSPACE Real Time Control (1kHz) - Position Control (PID), Feed-forward Friction Compensation - Read ADC - Quadrature Encoder decoding - Write DAC

Encoder Feedback Motor Current Command Fig. 9. Five-axis CNC test bed.

Table 2 Friction parameters for each axis. Parameters Tþ stat (N m) Ωþ 1 (rad/s)

Tþ coul (N m) Ωþ 2 (rad/s) T stat (N m) Ω 1 (rad/s) T coul (N m) Ω 2 (rad/s)

X-axis

Y-axis

Z-axis

A-axis

C-axis

3.1392 7.4600 2.3489

2.3553 4.8200 1.2032

1.9401 4.2000 1.7132

0.1033 7.3600 0.2172

0.1504 7.9000 0.2490

7.4600  2.0692  7.5700  1.5685  6.9200

4.8200  2.4599  6.6700  1.4039  6.6700

3.6000  1.4852  10.0000  1.1141  7.9000

10.0000  0.1908  7.6200  0.2795  10.0000

10.0000  0.1577  7.0000  0.2496  10.0000

Table 1 Drive and control parameters for each axis. Parameters Drive ka (A/V) kt (N m/A) rg (mm/rad) J (kg m2) B (kg m2/s)

X-axis

6.5723 0.4769 1.5915 7.00e  3 0.0236

Y-axis

6.2274 0.4769 1.5915 8.19e  3 0.0430

Z-axis

6.4841 0.4769 1.2732 7.67e  3 0.0323

A-axis

1.2700 0.3333 0.0056 8.11e  5 0.0011

High bandwidth (55.8 Hz) PID controller 30.000 kp (V/mm) ki (V/(mm s)) 650.00 kd (V/(mm/s)) 0.4000

37.031 802.34 0.4905

41.640 902.19 0.5535

736.64 15 960 9.4701

Low bandwidth (38.6 Hz) PID controller 10.000 kp (V/mm) ki (V/(mm s)) 50.000 kd (V/(mm/s)) 0.3000

12.344 61.718 0.3670

13.880 69.400 0.4147

245.55 1227.7 7.0147

C-axis

1.2700 0.3333 0.0111 2.21e  4 0.0029 1005.1 21 777 12.943 335.03 1675.1 9.5925

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

7

Fig. 10. Feed-forward friction compensation.

Fig. 12. Comparison of simulated contour errors with tracking error compensation and contour error compensation methods when the drives are controlled by PID controllers.

Fig. 11. Comparison of simulated contour errors with sliding mode controller, P–PI controller with contour error compensation, and P–PI controller without contour error compensation

similarly high contouring error as the P–PI controller. However, the P–PI controller with the contouring error compensation gives less path errors than the non-linear sliding mode controller applied to the individual drives with higher bandwidth. The example demonstrates the need to pre-compensate the trajectory commands either at the path contouring error or at the axis tracking error stage. Theoretically, it is possible to apply the same principles of compensation using axis tracking errors and achieve smaller contouring errors. As seen in Fig. 12, using contour error precompensation and tracking error pre-compensation for a PID controller reveal very similar contour error results. However, since the magnitude of tracking errors is higher, the potential for the tracking error compensation signals to excite the harmonic of drives is higher. As seen in Fig. 13, the tracking error compensation commands have higher magnitudes at a wider frequency range. As a result, the contour error pre-compensation can achieve the same results as tracking error compensation but with less lower frequency content to excite the natural modes of the machine during high speed-high feed machining. A set of simulations have been conducted to show the effectiveness of contour error compensation over axis tracking error compensation using cascaded P–PI controller used in commercial CNCs. The simulations revealed that the contouring performance decreased with tracking error pre-compensation. Fig. 14 shows the tracking errors at each axis for uncompensated, compensated axis tracking, and compensated contouring errors. Naturally the axis tracking errors are small when they are pre-compensated. However, the path contour errors even increase with the axis tracking error compensation, while the proposed path trajectory precompensation gives much smaller contouring errors as seen in Fig. 15. The P–PI structure of the controller tries to control both position and velocity of the drive, and it has to be more aggressive on drives even for a smaller contouring error. Also, any uncertainty in the drive dynamics will amplify the errors in axis tracking error predictions, which will make the control action more noisy and aggressive. However, since the reconstructed five-axis trajectory is

Fig. 13. Fast Fourier transform of tracking error compensation and contour error compensation reference commands when the drives are controlled by P–PI controllers.

smoothed, the contour error compensation has been observed to be more smooth. When changing the trajectory profile, the demands in velocity, acceleration, and jerk will also change. However, the trajectory modifications are small enough so that the limits should not change considerably, leaving the original feed profile valid. Furthermore, the preserved continuity of the proposed method prevents discontinuous changes in velocity, acceleration, and jerk. The comparison of the kinematic demands of the compensated and uncompensated trajectory in the x-axis for the low bandwidth PID controller can be seen in Fig. 16. Although very slight changes are experienced, should the new profile violates any assigned limit, a feed optimization algorithm [14,15] can be applied recursively, with a safety margin on the upper limit of acceleration and jerk. Then the pre-compensation can be performed again. In this respect, tracking error compensation may be limited. Since the contour error is the shortest distance that must be compensated, the trajectory will only result in minimal changes, as shown in Fig. 16. To validate the proposed method, several five-axis contouring experiments have been conducted. First, the prediction of tracking errors is validated. The tool path used in all experiments is shown in Fig. 7. Due to its high and varying curvature, this tool path is ideal in testing the effectiveness of the proposed algorithm. The predicted tracking errors for all five axes are compared favorably against actual measurements. Fig. 17 shows the comparison of

8

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

Fig. 16. Comparison of kinematic profiles of uncompensated and contour error precompensated reference commands using a low bandwidth PID controller.

Fig. 14. Comparison of simulated tracking errors with tracking error compensation and contour error compensation methods when the drives are controlled by cascaded P–PI controllers.

Fig. 15. Comparison of simulated contour errors with tracking error compensation and contour error compensation methods when the drives are controlled by cascaded P–PI controllers.

predicted and experimental results for the low bandwidth PID controller. The feed-forward friction compensation was active during the experiments, hence the prediction of tracking errors using the linear transfer functions of drives performed well.

It should be noted that the prediction of the tracking error is relatively less precise in the rotary axes. Due to the worm gear of the rotary axes, more non-linear effects exist. Since the presented tracking error method works best with systems that exhibit mostly linear characteristics, the rotary axes could not be modeled as accurately as the linear axes, resulting in more errors in the tracking error model. Regardless of this, the results still match fairly close, making the tracking error model acceptable for contour error prediction and compensation. To validate the proposed contour error compensation technique, experimental tests have been performed with two different PID controllers. The tests have been repeated with and without applying the proposed pre-compensation of contouring errors for the PID controllers. The tool tip and tool orientation contour errors, when using uncompensated and pre-compensated commands, are evaluated from Eqs. (19) and (21), compared in Fig. 18 for the high bandwidth controller, and Fig. 19 for the low bandwidth controller. The results are further summarized in Tables 3 and 4. The tool tip and tool orientation contour errors are reduced at an average of 56% and 51% for the high bandwidth controller, and 69% and 57% for the low bandwidth PID controller, thus validating the proposed method's effectiveness on linear controllers of various bandwidths. P–PI axis controller requires velocity feedback in its PI velocity control loop. Since the rotary drives were not equipped with the velocity sensors (i.e. tachogenerators), it was impossible to evaluate the velocity by digitally integrating the coarse encoders colored with friction in the worm gear. As a result, P–PI controller was not possible to be implemented on the experimental machine. However, the proposed method would work on the present rotary drives which are directly driven by the linear motors without needing the worm gear. Though the reduction in contour error is significant, some errors still exist. Since a continuous spline tool path is used as the reference, instead of discrete data, the contour errors should theoretically be calculated correctly, even in areas of high curvature. As a result, the source of errors is mainly from uncertainties in the drive models, which are reflected on the prediction of tracking and contouring errors. Mechanical interfaces create nonlinear phenomenon such as static friction or back lash that affect the accuracy of the model. In this work, the ball screws are considered to be rigid bodies which may not be entirely the case. By using more sophisticated modeling techniques, it is possible to improve the proposed method. However, complications may arise when attempting to preserve jerk continuity using the more sophisticated modeling techniques.

Predicted Experiment

0 −0.1 −0.2

0

0.5

1

1.5

2

2.5

Predicted Experiment

0.1 0 −0.1 −0.2

0

0.5

A−axis tracking error (rad)

Time (s) 4

x 10−3 Predicted Experiment

0 −2

0

1.5

2

2.5

9

0.2 Predicted Experiment

0.1 0 −0.1 −0.2

0

0.5

Time (s)

2

−4

1

0.5

1

1.5

2

2.5

C−axis tracking error (rad)

0.1

0.2

Z−axis tracking error (mm)

0.2

Y−axis tracking error (mm)

X−axis tracking error (mm)

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

6

1

1.5

2

2.5

Time (s)

x 10−3 Predicted Experiment

4 2 0 −2 −4

0

0.5

Time (s)

1

1.5

2

2.5

Time (s)

Fig. 17. Comparison of predicted and experimental tracking errors for the low bandwidth PID controller.

Table 3 Comparison of experimental contour errors for uncompensated and contour error pre-compensation for the high bandwidth PID controller. Trajectory Type

Uncompensated Pre-compensated

Tool tip contour error

Orientation contour error

Max. (μm)

Mean (μm)

Max. (mrad)

Mean (mrad)

43.10 20.00

12.60 5.50

0.56 0.27

0.10 0.05

Table 4 Comparison of experimental contour errors for uncompensated and contour error pre-compensation for the low bandwidth PID controller. Trajectory Type

Tool tip contour error (mm)

Fig. 18. Comparison of experimental contour errors for uncompensated and contour error pre-compensation for the high bandwidth PID controller.

Orientation contour error

Max. (μm)

Mean (μm)

Max. (mrad)

Mean (mrad)

265.26 81.37

85.62 26.72

2.30 0.99

0.70 0.30

5. Conclusions 0.4 Uncompensated Pre−compensated

0.3 0.2 0.1 0 0

0.5

1

1.5

2

2.5

Time (s) Orientation contour error (rad)

Uncompensated Pre-compensated

Tool tip contour error

2.5

x 10−3 Uncompensated Pre−compensated

2 1.5 1 0.5 0 0

0.5

1

1.5

2

2.5

Time (s)

Fig. 19. Comparison of experimental contour errors for uncompensated and contour error pre-compensation for the low bandwidth PID controller.

The tracking errors in five-axis contour machining may lead to significant contouring errors which may violate the tolerance of parts, unless low feed rates are used. The contouring errors are especially the main obstacle in high speed machining when sculptured surfaces with sharp curvatures are machined. The compensation of contouring errors using control laws is limited, since the axis controller usually acts after the error is imprinted on the surface of the machined part. Also, feed-forward precompensation of axis tracking errors requires highly accurate model of the drive's transfer function which has always uncertainties. This paper presents a pre-compensation of contouring errors by estimating them ahead of positioning the tool. The tracking errors of each drive are estimated analytically by using the transfer function of drives and fifth order polynomial model of reference axis commands. The jerk, acceleration and velocities are kept continuous along the curved path, allowing a smooth motion profile which could be analytically differentiated to estimate the

10

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

tracking errors accurately. The contouring errors are estimated; trajectory is regenerated smoothly with contour error compensation and decoupled into axis components using the kinematic model of the machine. It is shown that the contouring errors can be reduced significantly provided that the transfer function of the drives is kept linear by feed-forward compensation of friction disturbances.

which is denoted as L0 β ¼ ξ0 :

The velocity, acceleration and jerk continuity constraints at the point ðt i ; xi Þ can be written as 8 S_ ðh ÞS_ iþ1 ð0Þ ¼ 0; > > < i i S€ i ðhi ÞS€ iþ1 ð0Þ ¼ 0; ðA:8Þ >  > :  S ðh ÞS ð0Þ ¼ 0; i

Acknowledgments

ðA:7Þ

i

iþ1

which can be expanded as The study is supported by the National Key Basic Research Project of China (2011CB302400) and by National Science and Engineering Research Council of Canada under Strategic Canadian Network Grant in Machining Research (NSERC CANRIMT). Appendix A. Quintic spline fitting The overall jerk objective in Eq. (8) can be split into spline segments as Z h1  Z hn  J¼ S 1 ðtÞ2 dt þ … þ S n ðtÞ2 dt≔J1 þ … þ Jn : ðA:1Þ 0

0

The ith local objective can be derived as Z hi ð60b5;i t 2 þ 24b4;i t þ 6b3;i Þ2 dt Ji ¼ 0   2 5 4 2 3 ¼ 720b5;i hi þ 720b5;i b4;i hi þ 240b5;i b3;i þ 192b3;i hi

2

2hi

1

0

⋮

3hi

2 12hi

0

6hi

2

0

0

⋮

0

24hi

6

0

0

0

⋮

0

0

0

1

0

0

2

0

0

6

0

0

0

3 2 3 0 " # 0 7 βi 6 7 07 5 βiþ1 ¼ 4 0 5: 0 0

ðA:9Þ Denoting the matrix on the left hand side of Eq. (A.9) as derivative boundary condition are grouped together as 2 1 32 3 β1 L1 Ln 6 76 7 " 0 # 1 n β 6 76 2 7 L2 L 6 76 7 ¼ ⋮ ; 6 76 ⋮ 7 ⋱ ⋱ 4 54 5 0 3ðn1Þ1 βn L1n1 Ln

½L1i Ln ,

the

ðA:10Þ

ðA:11Þ

Setting the initial velocity, acceleration and jerk to be zero, the following conditions need to hold: 

S_ 1 ð0Þ ¼ S€ 1 ð0Þ ¼ S 1 ð0Þ ¼ 0; which can 2 0 0 0 60 0 0 4 0 0 6

ðA:12Þ

be expanded as 0

1

0

⋮

2

0

0

⋮

0

0

0

⋮

32

3

2 3 0 6 7 6 7 036ðn1Þ 7 54 ⋮ 5 ¼ 4 0 5: βn 0 β1

ðA:13Þ

Denote Eq. (A.13) as ðA:2Þ

The vector β contains the coefficients of all spline segments S(t) (Eq. (7)). The quintic spline has to pass through the reference points used in the fitting process ðt 0 ; x0 Þ; ðt 1 ; x1 Þ; …; ðt n ; xn Þ. For the ith segment, the following position boundary conditions need to hold, i.e., for x-axis:

which can be written as " # " # 0 0 0 0 0 1 xi1 ¼ : β 5 4 3 2 hi hi hi hi hi 1 i xi

3

4hi

L1 β ¼ ξ1 :

From Eqs. (A.(1) and A.2), the overall jerk objective in entire reference trajectory (Eq. (8)) can be evaluated by assembling local objective functions of all spline segments as 2 32 3 β1 i H1 1h T 76 ⋮ 7 1 T T 6 ⋱ β1 …βn 4 ðA:3Þ J¼ 54 5≔ β Hβ: 2 2 βn Hn

Si ð0Þ ¼ xi1 ; Si ðhi Þ ¼ xi ;

4

5hi 6 6 20h3i 4 2 60hi

which is denoted as

 1 2 2 b b b b b b þ144b4;i b3;i hi þ 36b3;i hi ¼ 2 5;i 4;i 3;i 2;i 1;i 0;i 3 2 32 5 4 3 b5;i 1440hi 720hi 240hi 0 0 0 6 7 6 7 3 2 6b 7 6 720h4i 384hi 144hi 0 0 0 7 6 76 4;i 7 6 76 b3;i 7 3 2 6 7 6 144hi 72hi 0 0 07 6 240hi 7 76 6b 7 6 0 0 0 0 0 07 6 76 2;i 7 7 6 76 4 0 0 0 0 0 0 54 b1;i 5 b 0;i 0 0 0 0 0 0 1 ≔ βTi H i βi : 2

2

ðA:4Þ

ðA:5Þ

Denoting the matrix on the left hand side of Eq. (A.5) as L0i , and the right hand side as ξ0i , the knot position constraints from all segments can be written as 2 32 3 2 3 β1 ξ0 L01 6 76 7 6 1 7 6 74 ⋮ 5 ¼ 4 ⋮ 5; ⋱ ðA:6Þ 4 5 βn ξ0n L0n

L2 β ¼ ξ2 :

ðA:14Þ

Similarly, by setting the final velocity, acceleration and jerk to be zero, the following conditions need to hold: 

S_ n ðhn Þ ¼ S€ n ðhn Þ ¼ S n ðhn Þ ¼ 0; which can be expanded as 2 4 3 ⋮ 5hn 4hn 6 3 6 036ðn1Þ ⋮ 20hn 12h2n 4 2 ⋮ 60hn 24hn

ðA:15Þ

2

3hn

2hn

1

6hn

2

0

6

0

0

32

3 2 3 β1 0 76 7 6 7 ⋮ 5 ¼ 4 0 5: 07 4 5 βn 0 0 0

ðA:16Þ Denote Eq. (A.16) as ðA:17Þ L3 β ¼ ξ3 : By combining all the constrains expressed in Eqs. (A.(7), A.11), (A.14), and (A.17) together in the following form: Lβ ¼ ξ

ðA:18Þ

where L ¼ ½L0 L1 L2 L3 T ; ξ ¼ ½ξ0 ξ1 ξ2 ξ3 T ; the minimum jerk optimization problem can now be expressed as min J ¼ β

1 T β Hβ 2

subject to : Lβξ ¼ 0:

ðA:19Þ

The above constrained quadratic minimization problem is solved using Lagrange Multiplier method. Introducing the vector of

K. Zhang et al. / International Journal of Machine Tools & Manufacture 74 (2013) 1–11

Lagrange Multiplier Λ, the augmented objective function is written as J′ðβ; ΛÞ ¼

1 T β Hβ þ ΛT ðLβξÞ: 2

ðA:20Þ

Equating ∂J′=∂β ¼ 0 and ∂J′=∂Λ ¼ 0 yields the linear equation system " #  " # 0 β H LT ¼ : ðA:21Þ ξ Λ L 0 The solution of linear equations (Eq. (A .21)) yields the coefficients β of S(t), which has a third order continuity and minimized jerk. References [1] Y. Altintas, A. Verl, C. Brecher, L. Uriarte, G. Pritschow, Machine tool feed drives, CIRP Annals – Manufacturing Technology 60 (2) (2011) 779–796. [2] M. Tomizuka, Zero phase error tracking algorithm for digital control, ASME Journal of Dynamic Systems, Measurement, and Control 109 (1) (1987) 65–68. [3] M. Weck, G. Ye, Sharp corner tracking using the IKF control strategy, CIRP Annals – Manufacturing Technology 39 (1) (1990) 437–441. [4] D. Renton, M. Elbestawi, High speed servo control of multi-axis machine tools, International Journal of Machine Tools and Manufacture 40 (4) (2000) 539–559. [5] Y. Altintas, K. Erkorkmaz, W. Zhu, Sliding mode controller design for high speed feed drives, CIRP Annals – Manufacturing Technology 49 (1) (2000) 265–270. [6] K. Erkorkmaz, Y. Altintas, High speed CNC system design. Part III: high speed tracking and contouring control of feed drives, International Journal of Machine Tools and Manufacture 41 (11) (2001) 1637–1658. [7] A. Kamalzadeh, K. Erkorkmaz, Accurate tracking controller design for highspeed drives, International Journal of Machine Tools and Manufacture 46 (10) (2006) 1124–1138.

11

[8] C. Okwudire, Y. Altintas, Minimum tracking error control of flexible ball screw drives using a discrete-time sliding mode controller, ASME Journal of Dynamic Systems, Measurement, and Control 131 (5/051006) (2009) 1–12. [9] Y. Koren, Cross-coupled biaxial computer control for manufacturing systems, ASME Journal of Dynamic Systems, Measurement, and Control 102 (4) (1980) 265–272. [10] Y. Koren, C. Lo, Variable-gain cross-coupling controller for contouring, CIRP Annals – Manufacturing Technology 40 (1) (1991) 371–374. [11] Y. Altintas, B. Sencer, High speed contouring control strategy for five-axis machine tools, CIRP Annals – Manufacturing Technology 59 (1) (2010) 417–420. [12] C. Ernesto, R. Farouki, Solution of inverse dynamics problems for contour error minimization in CNC machines, International Journal of Advanced Manufacturing Technology 49 (5) (2010) 589–604. [13] Y. Altintas, M. Khoshdarregi, Contour error control of CNC machine tools with vibration avoidance, CIRP Annals – Manufacturing Technology 61 (1) (2012) 335–338. [14] Y. Altintas, K. Erkorkmaz, Feedrate optimization for spline interpolation in high speed machine tools, CIRP Annals – Manufacturing Technology 52 (1) (2003) 297–302. [15] B. Sencer, Y. Altintas, E. Croft, Feed optimization for five-axis CNC machine tools with drive constraints, International Journal of Machine Tools and Manufacture 48 (7) (2008) 733–745. [16] G. Birkhoff, G. Rota, Ordinary Differential Equations, John Wiley, New York, 1969. [17] A. Yuen, K. Zhang, Y. Altintas, Smooth trajectory generation for five-axis machine tools, International Journal of Machine Tools and Manufacture 71 (2013) 11–19. [18] B. Sencer, Y. Altintas, Modeling and control of contouring errors for five-axis machine tools. Part I: modeling, ASME Journal of Manufacturing Science and Engineering 131 (3–031006) (2009) 1–8. [19] K. Erkorkmaz, Y. Altintas, High speed CNC system design. Part II: modeling and identification of feed drives, International Journal of Machine Tools and Manufacture 41 (10) (2001) 1487–1509.