Interpolation and Extrapolation of Optimally Fitted Kinematic Error Model for Five-Axis Machine Tools

Machine tool accuracy is critically important in many industrial applications such as dies and molds [1], large mechanical parts [2], and flank milling for free-form surfaces [3]. A machining process includes many error sources that lead to the inaccuracy of the final parts. Those error sources include geometric errors, kinematic errors, thermal errors, and dynamic errors and are categorized and discussed in Ref. [4]. Substantial research has been conducted in the improvement of machine tool and machining process accuracy. In Ref. [5], a general mathematical model was proposed to consider both geometric errors and thermal effects. A calibration procedure was developed to compensate those errors in real time. A new method was proposed in Ref. [6] to measure the geometric errors of the rotary axes of five-axis machine tools more efficiently and accurately. Regarding model validation, Hashemiboroujeni et al. [7] used a scale and master balls artefact (SAMBA) calibration method, which identified eight-axis location errors, two spindle offsets, and three linear gain errors, on a five-axis machine tool and validated the method with a ball dome artefact instead of a coordinate measuring machine. The proposed method improves both accuracy and efficiency. For machining process, a method based on mirror compensation was proposed in Ref. [8] to compensate errors of parts having free-form surfaces. A virtual machine tool also requires the true machine to be accurate enough to be consistent with the simulated machining process [9], which requires accurate calibration.

In a general machine tool calibration process, accuracy is achieved by using various metrology instruments to measure link lengths variations, assembly offsets, and alignment errors that generate kinematic errors, and then correcting the errors by mechanical adjustments or constructing kinematic error models used to generate compensating position command algorithms [10–12]. While kinematic offsets and alignments can be directly measured, it is often time-consuming to measure each source of error independently. Indirect measurement methods measure and analyze tool tip volumetric inaccuracy, which is treated as a contribution of all geometric error sources. Models are constructed that are used to correct all errors simultaneously. For five-axis machine tools, typical metrology instruments for indirect methods include ball bars [13,14], R-test [15,16], touch-trigger probes [17,18], laser trackers [19,20], and machining tests [21]. Indirect measurement methods are reviewed and discussed in Ref. [22]. An example of a kinematic error modeling method using indirect measurement methods is given in Ref. [2]. In this study, a laser tracker is used to quickly gather error measurements throughout the axis space, and a maximum likelihood method is used to construct kinematic error models with high-order, position-dependent translational and rotational errors for each axis. The method was applied to an industrial five-axis machine tool and was able to reduce volumetric error by an order of magnitude.

Curve fitting is an essential element in indirect measurement methods for geometric error compensation, although not a topic well addressed in the literature. In curve fitting, the finite measurements are extended across the entire working volume by the fitted curve, interpolating to regions surrounded by measurements and extrapolating to regions outside of the measurements. While most curve-fitting methods focus on the accuracy of the fitted curve at the measurement locations, accuracy in interpolated and extrapolated regions is highly dependent on the curve-fitting method. In this paper, a method of controlling the fitting process to ensure the models generate realistic solutions over the interpolation and extrapolation spaces is proposed.

Extrapolation is used for machine tool compensation in situations where geometric errors can only be identified over a limited range of the machine tool’s workspace due to the design and size of the metrology instruments. For devices such as ball bars, R-tests, and touch-trigger probes, which must maintain contact with the machine tool’s spindle and table, orientation errors may be difficult to measure near the table boundaries. For other metrology instruments such as laser interferometers and laser trackers, no contact between the spindle and the instrument is needed and, thus, measurements can be collected over a much larger volume. However, limitations still exist. Figures 1 and 2 give two examples. In Fig. 1, the minimum commanded Z axis when machining is lower than the minimum commanded Z axis when collecting data. Thus, all of the compensation tables having Z axis values less than the minimum commanded Z axis when collecting data cannot be populated unless the geometric error functions are extrapolated. The minimum commanded Z axis can be decreased, as shown in Fig. 2, if the machine tool spindle is attached to an A or B rotary axis. However, the cutting tool will be in a different orientation, illustrating that machine tools with more than three axes have unmeasured spaces in the joint space as opposed to the workspace. While it is good practice to use retroreflectors that are smaller than the tooling that will be used, state-of-the-art retroreflectors, which are automated to follow the laser tracker for improved efficiency, have form factors larger than many cutting tools. In addition to collision avoidance, the measured region is also limited due to a line of sight constraint when a laser tracker or an interferometer is used. Figure 3 shows the example from Ref. [2] that the measurements of a five-axis machine tool in some joint spaces are not collected due to collision avoidance and line of sight constraints. As shown in Fig. 3, the top plot illustrates the consequence of Fig. 2 and the bottom plot illustrates the consequence when the laser beam is blocked by the machine tool spindle. While the laser tracker can be repositioned to measure data in regions blocked by line of sight constraints, this will increase measurement time by requiring multiple setups and increase errors by creating data sets in multiple measurement frames, each of which have different uncertainties.

Compensation values over the unmeasured spaces are often set to the last compensation value. In this case, the geometric errors that exist at these points will be uncompensated, resulting in part geometry errors. Another method is to extrapolate the geometric error model, which are typically constructed from polynomials. In Ref. [23], the angular errors are modeled with second-order polynomials for a three-axis machine tool, while in Ref. [24] geometric errors are fitted with third-order polynomials as a function of axis position. In the two proposed models described in Ref. [2], error models with sixth- and eighth-order polynomials provided the best performance. In selecting a polynomial order, it should be noted that as noise in the data increases the polynomial order providing the least residual will increase. However, these high-order polynomials do not properly describe the characteristics of the machine’s kinematic error. This is especially evident when attempting to extrapolate the model. Polynomial curve fitting possesses provide good interpolation characteristics; however, they often have poor extrapolation properties [25]. The poor extrapolation phenomenon is especially obvious when using high-order polynomials. Although polynomial curve fitting provides small error approximation to the data, the fitted curve can “take off” sharply outside the data domain [26] due to the increase of the magnitude of the fitted curve slope. Thus, the fitting accuracy deteriorates rapidly outside the range of the data used to construct the model. An illustrative example is given in Fig. 4 in which 13 simulated measurements (with measurement errors) are fitted with a continuous modeled function. The slope of the modeled function changes rapidly and unrealistically when extrapolated, causing inaccuracy in predicting the end points and a continuous increase in the inaccuracy as the modeled function is further extrapolated. Therefore, a method of controlling the poor extrapolating behavior and, thus, improving extrapolated model accuracy is needed for machine tool geometric error models.

In addition to the unmeasured space limitation, curve fitting also suffers an interpolation process since indirect measurement methods use continuous polynomials to fit discrete measurements. Low-order polynomials will make the model less accurate on the collected data while high-order polynomials may overfit the measurements and, thus, lose accuracy. An illustrative example of using low- and high-order functions fitting 13 simulated measurements is shown in Fig. 5. Compared with the low-order function, the high-order function fits the actual measurements better. However, for the regions between the actual measurements, the fitted function has an oscillation, which is known as Runge’s phenomenon, when the high-order function is used [27]. Thus, the interpolation of machine tool geometric error models between actual measurements must be carefully considered and treated.

One may consider using the analysis tool error budgeting to predict and control a machine’s errors [28]. This method is commonly used at the design stage of machines [29]. Error budgeting is typically applied to account for economic factors. As examples, an improved value analysis method, which considers both machining accuracy and manufacturing cost, was used in Ref. [30] and the total cost was optimized in Ref. [31]. Error budgeting was used in Ref. [32] to consider both thermal and tool wear errors. Error budgeting with constraints naturally accounts for the stochastic nature of machine errors. An example is Ref. [33] that utilized Monte-Carlo simulation at the design stage. The method proposed in this paper focuses on machine calibration for geometric errors, as opposed to machine design, and computes a deterministic compensation while accounting for the stochastic nature of the machine and measurement device.

This paper proposes a method to construct machine tool geometric error models that interpolate and extrapolate well throughout the entire machine tool workspace. Based on an error model proposed in Ref. [2], an analytical form of the tool tip modeled error slope is described using all of the single-axis error slopes. By using a uniform constraint to bound the magnitude of the tool tip modeled error slope, all single-axis error slopes and, thus, magnitudes will be constrained. A systematic methodology is given to determine the constraint value. Compared with the unconstrained model, the proposed constrained model provides a more consistent error description over the entire machine tool workspace. The rest of this paper is organized as follows. Section 2 gives the geometric error compensation methodology used in this paper. Section 3 proposes a constraint design used to construct constrained models. Section 4 gives the experimental results implemented on an industrial five-axis machine tool, and a comparison between the unconstrained and constrained models is described and analyzed. The paper is summarized and conclusions are drawn in Sec. 5.

While there are many methods to describe and compensate machine tool geometric errors, this paper utilizes a table-based compensation methodology given in Ref. [2]. Most machine tool controllers offer a set of compensation tables that map single-axis commands into small corrections to a single-axis position in real time. This section briefly describes the model used for table-based compensation, denoted the axis perturbation model. The description of actual measurements and the identification of model parameters are also given.

An optimization algorithm based on the implicit loop method [36] is used to determine the error model parameters. An advantage of this method is that axis positioning errors and measurement errors can be directly incorporated. They are assumed to be independent and follow normal distributions. These errors are identified in the optimization process and, thus, are separate from the machine tool’s geometric errors.

As discussed previously, large slopes in the polynomials at the ends of the measurement data used to construct the polynomials substantially reduce extrapolation accuracy. In Ref. [26], extrapolation accuracy is improved by using the asymptotic behavior (the function slope approaches zero) of rational functions. However, rational functions will introduce more polynomial coefficients, bringing challenges to polynomial order selection and model identification. Also, zero slope will tend to over constrain the error model in the extrapolation region. Thus, while high-order polynomial functions are still used to describe the collected measurements, a constraint on the modeled function slope is added to improve the accuracy of the model function over the entire workspace. This will allow the methodology to construct models that not only fit the measured data but also capture the characteristics of the slope of the data, even outside of the measurement space. This allows the methodology to avoid large slopes at the ends of measurement space and provide better prediction in the extrapolated space. In this section, a constraint that aims to control all perturbation function slopes over the entire workspace is developed. A procedure to utilize this constraint for geometric error model construction will also be given.

J_ω,nom,i and J_ω,AP,i are the ith columns of the angular velocity portions of the nominal and modeled Jacobian matrices, respectively, and J_v,nom,i and J_v,AP,i are the ith columns of the linear velocity portion of the nominal and modeled Jacobian matrices, respectively. The machine tool Jacobian matrix is a matrix relating the small changes between the axes’ positions in joint coordinates and the end effector position in Cartesian coordinates [38].

To construct a geometric error model with consistent behavior over the entire workspace, one must determine the appropriate constraint value, i.e., the value c in Eq. (31). A large value of c will fail to constrain unrealistic model behavior when extrapolating, while a small value of c will over constrain the perturbation function slopes, reducing the model accuracy over the entire space. Boundary points, which by definition are the points located around the edges of the measurement data [39], are distributed at the intersection of the interpolation and extrapolation spaces. The slopes of the identified kinematic errors functions are very sensitive at the boundary points. Thus, the boundary points are used to determine the appropriate constraint value, c, which will give a balance between unconstrained and overconstrained models. To detect boundary points in a two-axis space, a database technique named BORDER [39] is used. Typically BORDER uses three steps to determine the boundary points. The first step is to find the k-nearest neighbors (kNN) for each point in the data set, where k is a tuned integer. The second step is to count the number of reverse k-nearest neighbors (RkNN) and the last step is to sort the points according to the RkNN number. As boundary space points tend to have fewer RkNN, a bound on the RkNN values is selected to determine which points will be considered to be boundary points.

To create a set of measurements in the axis space with low discrepancy, a quasi random sequence is used to generate the axis commands [40]. For both tool lengths, the same 295 quasi random measurements are collected. Table 1 gives the minimum and maximum axis commands and Fig. 9 gives the distributions of the measurements projected in the two-axis spaces.

The 295 measurements are divided into three sets, shown in the BZ axis space in Fig. 10. To analyze the ability of the methodology to extrapolate geometric error models, 25 measurements at the bottom of the BZ space are taken to be the extrapolation validation set. Note that while these points are measured and, thus, the kinematic error is known, for the purpose of this study they are assumed to reside in the unmeasurable space and extrapolation must be used to predict them. They will be used to validate the extrapolated model performance. Over the interpolation space, another 25 points are randomly selected to be the interpolation validation set. They will be used to validate the interpolated model performance. The remaining 245 measurements are used as the identification set for model construction.

As a baseline, a model will be constructed without imposing the slope constraint. This model will be referred to as the unconstrained model. According to Ref. [2], perturbation functions described by sixth-order polynomials are appropriate for this specific machine tool. Numerical optimization is used to minimize (11) for model identification. Here, the matlab optimization solver fmincon is used as it is capable of including nonlinear constraints during the optimization process.

The unconstrained model is constructed with the 245 identification points. Table 2 lists the mean and maximum residuals for the identification and validation sets of the nominal (i.e., uncompensated) machine and the unconstrained model. Figure 11 shows the compensation table functions, which are generated from the geometric error models, over the interpolation and extrapolation spaces. The model reduces the mean residual from 0.307 to 0.038 mm for the identification data set and 0.261 to 0.044 mm for the interpolation validation data set, providing 87.6% and 83.1% reductions, respectively. However, for the extrapolation validation data set, the constrained model increases the mean residual from 0.321 to 0.451 mm and the maximum residual from 0.525 to 1.909 mm. The reason for the poor model performance for the extrapolation validation set can be seen in Fig. 11. For the perturbation functions that are dependent on the Z axis position, the behavior in the extrapolation space for small values of Z is not consistent with the behavior in the interpolation space. In this case, the function slopes are much larger in the extrapolation space. This results in unrealistic error magnitudes in the extrapolation space and, thus, poor accuracy of the extrapolated model.

In Fig. 11, the perturbation functions dependent on the Z axis position have unrealistically large slopes for the small values of Z since no data in this range was used to construct the models. Thus, boundary space points located in the BZ space are used to determine the constraint value c. BORDER is applied to the 245 identification points. Here, k is selected as 30. As it is preferred to keep most measurements as identification points, only the first ten boundary space points identified by BORDER, denoted BORDER points, are used. The remaining 235 points are the identification subset points. Figure 12 gives the distribution of the BORDER points and the identification subset points.

Using the identification subset, constrained models are constructed with different constraint values. Figure 13 shows the performance of the identification subset points and the BORDER points. As the constraint goes from infinity (i.e., unconstrained model) to 0.4, the maximum residual of the BORDER points decreases while the mean residual of the identification subset remains nearly constant. This indicates that the interpolation space perturbation functions are refitted to best fit the identification subset points and the perturbation function slopes over the extrapolation space are being constrained. When the constraint is set smaller than 0.4, the perturbation function slopes are overconstrained. As a result, the model performance becomes worse for both the identification subset points and the BORDER points. From the trend of the maximum residual of the BORDER points, c = 0.4 mm/half axis motion is selected for constrained model construction. Note that for other machine tools or other tool lengths the constraint value may need to be modified.

The final constrained model is constructed with the original 245 identification points and c = 0.4. The model performance and the identified error curves are compared between the unconstrained and constrained models. Figure 14 shows the compensation table functions generated from the unconstrained and constrained models. Figure 15 shows the modeled residuals of the unconstrained and constrained models on the identification set and the extrapolation validation set with respect to the position of Z axis. Table 3 lists the mean and maximum residuals of the nominal machine and the unconstrained and constrained models for the identification and validation, both interpolation and extrapolation, data sets and Fig. 16 gives the mean and maximum residuals of the nominal machine and the constrained and unconstrained models for both validation data sets.

As shown in Fig. 14, the perturbation functions dependent on the Z axis position now have slopes at small values of Z that are consistent with the slopes of the perturbation functions in the nonextrapolated region. Of note are the sinusoidal errors seen in the X and Y axes as the C axis rotates. These errors are the result of a mechanical offset of the C axis and illustrate complex errors that cannot be properly captured by simple models. In Fig. 15, the unconstrained model residuals increase rapidly as the Z position decreases, indicating poor performance over the extrapolation space. This is consistent with the third column of the compensation table functions over the extrapolation space shown in Fig. 11. While the constrained model residuals increase as the Z position decreases, the magnitude and rate of increase of the constrained model residuals are significantly less than those of the unconstrained model. Further, the constrained model maintains consistent performance over the measured space. In Table 3 and Fig. 16, it can be seen that the constrained model performs much better than the unconstrained model for the extrapolation validation data. The constrained model reduces the mean residual by 40%, from 0.321 mm to 0.191 mm, and the maximum residual by 16%, from 0.525 mm to 0.443 mm, for the extrapolation validation data. For the interpolation validation data, the unconstrained model reduces the mean and maximum residuals from the nominal machine by 83%. This performance is the same as the calibration method proposed in Ref. [7] which reduces the maximum and average volumetric errors by 83% and 82%, respectively. The constrained model reduces the residuals by 79% for the interpolation validation data with only 4% loss comparing with the unconstrained model. Although the mean and maximum residuals for the interpolation validation data increase by 0.012 and 0.020 mm, respectively, from the unconstrained to constrained models, these increases are one and two orders of magnitude, respectively, less as compared with the decreases in the mean and maximum residuals, i.e., 0.26 and 1.466 mm, respectively, for the extrapolation validation data. The constrained model provides a more consistent description of the machine tool geometric errors over the entire machine tool space (i.e., interpolated and extrapolated).

Machine tool geometric errors are frequently corrected by constructing geometric error models using measurements from external metrology equipment. Due to measurement device limitations, extrapolation of error models often yields poor results. A new method of interpolating and extrapolating geometric error models is proposed in this paper. Based on the axis perturbation model, the proposed method uses the tool tip modeled error slopes at the boundary of the collected data as a constraint to ensure slopes in extrapolated regions are consistent with the slopes in interpolated regions. By adding the constraint during the model identification process, geometric errors in the interpolation regions are still well identified and geometric errors in the extrapolation regions behave more consistently.

The method was validated in an experimental study on an industrial five-axis machine tool. For the unconstrained model, although the maximum residual of the interpolation validation set for the nominal machine is reduced 83%, the maximum residual of the extrapolation validation set for the nominal machine increased 264% since this model fails to properly describe the geometric errors in the extrapolation regions. A constrained model is constructed and the slopes of the perturbation functions are shown to be dramatically more consistent over the entire machine tool workspace. While the constrained model only reduces the maximum residual of the interpolation validation set for the nominal machine by 79%, it reduced the maximum residual of the extrapolation validation set for the nominal machine by 16%. The new geometric error modeling method developed in this paper provides a systematic means to construct models that are valid when extrapolated, creating an alternative to ad hoc methodologies to populating machine tool compensation tables in extrapolation regions.

The authors gratefully acknowledge the financial support for this work from the National Science Foundation (Grant No. CMMI-1335340) and the Center for Aerospace Manufacturing Technologies at the Missouri University of Science and Technology.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.

c =

constraint value for constrained model construction

b =

vector of error parameters to be identified

q =

nominal axis command vector

S =

skew symmetric matrix

l_t =

measured tool length

$q¯i$ =

ith axis command scaled to range [−1 1]

e_t =

tool tip modeled error

p_a =

actual measurement

F_nom =

nominal kinematics

F_AP =

axis perturbation model kinematics

J_v =

linear velocity portion of machine tool Jacobian matrix

J_ω =

angular velocity portion of machine tool Jacobian matrix

T_i =

linear homogeneous transformation for ith axis

$h¯i$ =

tool tip modeled error change per half range motion of ith axis

E_ij(q_j) =

axis perturbation function mapping axis command q_j onto an adjustment for axis command q_i

δq =

axis command perturbation vector

ν =

axis positioning error vector

ξ =

measurement error vector

σ =

standard deviation

Σ_ν =

covariance matrix for axis positioning error

Σ_ξ =

covariance matrix for measurement error