Abstract

Machine tool geometric errors are frequently corrected by populating compensation tables that contain position-dependent offsets to each commanded axis position. While each offset can be determined by directly measuring the individual geometric error at that location, it is often more efficient to compute the compensation using a volumetric error model derived from measurements across the entire axis space. However, interpolation and extrapolation of measurements, once explicit in direct measurement methods, become implicit and obfuscated in the curve-fitting process of volumetric error methods. The drive to maximize model accuracy while minimizing measurement sets can lead to significant model errors in workspace regions at or beyond the range of the metrology equipment. In this paper, a method of constructing machine tool volumetric error models is presented in which the characteristics of the interpolation and extrapolation errors are constrained. Using a typical five-axis machine tool compensation methodology, a constraint bounding the tool tip modeled error slope is added to the error model identification process. By including this constraint over the entire space, the geometric errors over the interpolation space are still well identified. Also, the model performance over the extrapolation space is consistent with the behavior of the geometric error model over the interpolation space. The methodology is applied to an industrial five-axis machine tool. In the experimental implementation, for measurements outside of the measured region, an unconstrained model increases the mean residual by 40% while the constrained model reduces the mean residual by 40%.

1 Introduction

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 [1012]. 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.

Fig. 1
Fig. 1
Close modal
Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal

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.

Fig. 4
Fig. 4
Close modal

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.

Fig. 5
Fig. 5
Close modal

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.

2 Geometric Error Compensation Background

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.

2.1 Geometric Error Model Construction.

Using the zero reference model in Ref. [34] to describe the linear homogeneous transformations [35] between the coordinate systems of two axes, the nominal kinematics of an n-axis machine tool is
$Fnom(q)=T1(q1)T2(q2)…Tn(qn)$
(1)
where Fnom describes the nominal orientation and position of the machine tool’s last frame with respect to the machine tool’s base frame, q = [q1q2qn]T is the nominal axis command vector, and Ti is the nominal linear homogeneous transformation from the (i−1)th axis coordinate system to the ith axis coordinate system. Machine tool compensation tables are look-up tables that depend on the input nominal axis commands and contain small adjustments to the nominal axis commands [2]. The axis perturbation model is a type of kinematic error model that can be used to efficiently generate the compensation tables. The axis perturbation kinematic model is given by
$FAP(q)=Fnom(q+δq(q))$
(2)
where FAP describes the uncompensated machine and δq(q) = [δq1(q) δq2(q) … δqn(q)]T is a perturbation in the axis command vector, where δqi(q) is the error in each axis at the axis command qi. Compensation tables can be directly generated as −δq(q) such that the machine recovers nominal kinematics using the compensated input $q^=q−δq(q)$ or $FAP(q^)=Fnom(q)$.
Based on the structure of common machine tool compensation tables, each error is described as the sum of n perturbation functions
$δqi(q)=Ei1(q1)+Ei2(q2)+⋯+Ein(qn)$
(3)
where Eij(qj) is a table function (if it exists) that maps the command of axis j onto a perturbation for the command position of axis i [11]. To capture both constant and complex errors, each function is mathematically described by a polynomial with sufficient order. Here, a basis set of functions described by Chebyshev polynomials given on a normalized scale is used. Given λ ∈ [−1, 1], a Chebyshev polynomial has the form
$C(λ)=a0g0(λ)+a1g1(λ)+a2g2(λ)+⋯+amgm(λ)$
(4)
where
$g0(λ)=1,g1(λ)=λ,g2(λ)=2λ2−1,…,gm+1(λ)=2λgm(λ)−gm−1(λ)$
(5)
m is the Chebyshev polynomial order and a0, a1, a2, …, am are the polynomial coefficients. Thus, each perturbation function can be represented by mth order Chebyshev polynomials as
$Eij(q¯j)=a0,ijg0(q¯j)+a1,ijg1(q¯j)+a2,ijg2(q¯j)+⋯+am,ijgm(q¯j)$
(6)
where
$q¯j=2(qj−qj,min)(qj,max−qj,min)−1$
(7)
is the jth linearly mapped axis command by scaling the axis range to the interval [−1 1] and qj,min and qj,max are the minimum and maximum jth axis commands, respectively.

2.2 Measurement.

A laser tracker coupled with an active retroreflector attached to the machine tool spindle is used to measure machine position. To describe the measurements, a complete closed kinematic loop between the laser tracker and retroflector is needed. Considering the measurement errors and the potential axis positioning errors, the measured machine position can be described as
$pa(q)=TmfEmfFAP(q+ν)ptl+ξ$
(8)
where Tmf is the nominal transformation from the laser tracker frame to the machine tool’s base frame, Emf is the correction of Tmf, ν is a stochastic axis positioning error vector, ptl is a tool length vector to the retroreflector, and ξ is the measurement error vector. Figure 6 gives a schematic description of the actual measurement with respect to the laser tracker frame. The base frame correction, Emf, is described with a fixed six degree of freedom error matrix
$Emf=[1−EC0EB0EX0EC01−EA0EY0−EB0EA01EZ00001]$
(9)
where EA0, EB0, and EC0 are small rotations about the X, Y, and Z axes, respectively, and EX0, EY0, and EZ0 are small translations along the X, Y, and Z axes, respectively. The tool length vector can be modeled as
$ptl=[00lt+δt1]T$
(10)
where lt is the measured tool length between the origin of the machine tool’s last frame and the retroreflector and δt is a correction to lt. The retroreflector can only determine 3D positional information. To obtain orientation error, two sets of measurements are taken, each time with the retroreflector mounted on a tool with a different length. Both groups of measurements use the same axis commands and, thus, lie on the same spindle axis orientation.
Fig. 6
Fig. 6
Close modal

2.3 Model Identification.

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.

Let Σν and Σξ be the covariance matrices for the axis positioning and measurement errors, respectively. Considering N measured machine tool poses, the kinematic errors are identified by minimizing the following function [37]
$[ν1*,…,νN*,b*]=argminν1,…,νN,b∑k=1N(νkTΣν−1νk+ξkTΣξ−1ξk)$
(11)
subject to the measurement constraints
$ξk=[ξk,1ξk,2]=[pa,k,1(qk)−TmfEmf(b)Fnom(qk+δqk(b,qk)+νk,1)ptl,1(b)pa,k,2(qk)−TmfEmf(b)Fnom(qk+δqk(b,qk)+νk,2)ptl,2(b)]T$
(12)
where $ξk=[ξk,1Tξk,2T]T$ is the measurement error vector and $νk=[νk,1Tνk,2T]T$ is the axis positioning error vector for the kth pose measured by the retroreflector mounted on two tools with different lengths, b is the error model parameter vector including the polynomial coefficients in the modeled machine tool kinematics, six static errors in Emf and the tool length errors δl,1 and δl,2 corresponding to the two tool length vectors ptl,1 and ptl,2, respectively.

3 Constrained Geometric Error Model

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.

3.1 Constraint Design.

Let the nominal and modeled transformations from the machine tool’s base frame to the machine tool’s last frame, respectively, be
$Fnom=[Rnompnom01×31]$
(13)
$FAP=[RAPpAP01×31]$
(14)
where Rnom and RAP are the nominal and modeled rotations from the machine tool’s base frame to the machine tool’s last frame, respectively, and pnom and pAP are nominal and modeled positions of the machine tool’s last frame with respect to the machine tool’s base frame, respectively. Let the measured tool length vector be
$p¯tl=[00lt]T$
(15)
To combine the rotation and position information into one representation, the nominal and modeled tool tip positions with respect to the machine tool’s base frame are used
$pt,nom=Rnomp¯tl+pnom$
(16)
$pt,AP=RAPp¯tl+pAP$
(17)
The gradients of Eqs. (16) and (17) with respect to the axis command vector q, respectively, are
$∇qpt,nom=[∂pt,nom∂qX∂pt,nom∂qY∂pt,nom∂qZ∂pt,nom∂qC∂pt,nom∂qB]$
(18)
$∇qpt,AP=[∂pt,AP∂qX∂pt,AP∂qY∂pt,AP∂qZ∂pt,AP∂qC∂pt,AP∂qB]$
(19)
For the ith axis
$∂pt,nom∂qi=∂Rnomp¯tl∂qi+∂pnom∂qi=S(Jω,nom,i)Rnomp¯tl+Jv,nom,i$
(20)
$∂pt,AP∂qi=∂RAPp¯tl∂qi+∂pAP∂qi=S(Jω,AP,i)RAPp¯tl+Jv,AP,i$
(21)
where the skew symmetric matrix operator for a vector a = [a1a2a3] is
$S(a)=[0−a3a2a30−a1−a2a10]$
(22)

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 Jv,nom,i and Jv,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].

The difference between Eqs. (20) and (21) is taken, giving the slope of the tool tip modeled error, et, with respect to the ith axis
$∂et∂qi=∂(pt,nom−pt,AP)∂qi=S(Jω,nom,i)Rnomp¯tl+Jv,nom,i−S(Jω,AP,i)RAPp¯tl−Jv,AP,i=(S(Jω,nom,i)Rnom−S(Jω,AP,i)RAP)p¯tl+(Jv,nom,i−Jv,AP,i)$
(23)
Assuming the difference between Rnom and RAP is small, such that
$Rnom≈RAP$
(24)
Equation (23) can be written as
$∂et∂qi=(S(Jω,nom,i)Rnom−S(Jω,AP,i)Rnom)p¯tl+(Jv,nom,i−Jv,AP,i)=S(Jω,nom,i−Jω,AP,i)Rnomp¯tl+(Jv,nom,i−Jv,AP,i)=S(ΔJω,i)Rnomp¯tl+ΔJv,i$
(25)
where ΔJω,i is the difference between the ith columns of the angular velocity portions of the nominal and modeled Jacobian matrices and ΔJv,i is the difference between the ith columns of the linear velocity portions of the nominal and modeled Jacobian matrices.

3.2 Jacobian Matrix.

The detailed derivation of Jω,AP, Jω,nom, Jv,AP, and Jv,nom are shown in the  Appendix. Using the derived formulations, the differences in the linear and angular Jacobian matrices, respectively, are
$ΔJv=−[dEXX(qX)dqXdEXY(qY)dqYdEXZ(qZ)dqZdEXC(qC)dqCdEXB(qB)dqBdEYX(qX)dqXdEYY(qY)dqYdEYZ(qZ)dqZdEYC(qC)dqCdEYB(qB)dqBdEZX(qX)dqXdEZY(qY)dqYdEZZ(qZ)dqZdEZC(qC)dqCdEZB(qB)dqB]$
(26)
$ΔJω=[000000000001−sin(qC)cos(qC)0]−A[dEBX(qX)dqXdEBY(qY)dqYdEBZ(qZ)dqZdEBC(qC)dqC1+dEBB(qB)dqBdEBX(qX)dqXdEBY(qY)dqYdEBZ(qZ)dqZdEBC(qC)dqC1+dEBB(qB)dqBdECX(qX)dqXdECY(qY)dqYdECZ(qZ)dqZ1+dECC(qC)dqCdECB(qB)dqB]$
(27)
where
$A=[−sin(qC+δqC(q))000cos(qC+δqC(q))0001]$
(28)

3.3 Constrained Model Construction.

As shown in Eqs. (26) and (27), the linear and angular Jacobian differences depend on the slope of each of the perturbation functions. Let Eq. (25) be rewritten as
$hi=∂et∂qi=S(ΔJω,i)Rnomp¯tl+ΔJv,i$
(29)
which represents the tool tip modeled error change per axis unit. Note that the unit of hi is mm/mm for translational axes and mm/deg for rotational axes. A unification of the unit is made
$h¯i=∂et∂q¯i=∂et∂qi(qi,max−qi,min)2=hi(qi,max−qi,min)2=(S(ΔJω,i)Rnomp¯tl+ΔJv,i)(qi,max−qi,min)2$
(30)
which represents the tool tip modeled error change per half range motion for the ith axis. By constraining the value of Eq. (30) during model identification, the magnitude of all perturbation function slopes will be regulated. Previously, Eqs. (11) and (12) are used for model parameter identification. To construct constrained geometric error models and control the error function slopes, the following constraint is added to Eq. (12) when optimizing (11) during the model parameter identification process
$‖h¯i‖2≤c(i=X,Y,Z,C,B)$
(31)
where c is the constraint applied to all machine tool axes throughout the entire workspace.

3.4 Constraint Value Determination.

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.

4 Experimental Results

4.1 Experimental Setup.

An industrial five-axis machine tool with axis sequence XYZCB and a Siemens 840D controller is used for the experimental studies. The machine tool moves along a gantry in the X direction and the X, Y, Z, C, and B axes have ranges of 6.1 m, 2.6 m, 1 m, 544 deg, and 222 deg, respectively. This machine tool is a “gantry with a two-rotary-axes head” type, one of the three types of common machine tool configurations described in Ref. [4]. While our extrapolation method will be demonstrated on this type of machine tool, our method is also applicable to the other two types of machine tools, namely, “transverse column with a two-rotary-axes tilting table” and “a machine with a swiveling spindle and rotating plate” [4] since these types of machine tools are also described by a kinematic chain and their errors are corrected using machine tool compensation tables. Figures 7 and 8, respectively, give a picture and structural schematic of the machine tool used in this paper. An Automated Precision Inc., T3 laser tracker, is located on the machine tool table and an Active Target (AT) retroreflector is mounted in the spindle. The laser tracker has a measurement error of 5 μm/m. The average distance between the laser tracker and the measured points is approximately 2.5 m, thus, the average measurement error from the laser tracker is 12.5 μm. The manufacturer specifications for the Active Target report a measurement error of ±12.5 μm, which is treated as deterministic. Before taking measurements, the laser tracker and the machine were fully warmed up using specified procedures and the measurement process was implemented in a relatively temperature-controlled environment which can be treated as a thermally isolated environment. In Fig. 8, lBs is the length between the B axis rotating center and the machine tool spindle surface. Nominally, lBs = 98 mm. The lengths between the spindle surface and the AT retroreflector, lst, for the two AT mountings are 304.9 mm and 402.9 mm for the short and long tool length, respectively. Thus, the total tool length is
$lt=lBs+lst$
(32)
and the short tool and long tool lengths are 402.9 mm and 500.9 mm, respectively. Each the short or long tool was attached to the spindle, the spindle was rotated and circle measurements were collected, and the data were fitted. Then a wobble plate connecting the Active Target and spindle was adjusted until the radius of the fitted circle was smaller than 25 μm. This process ensures the short and long tools are coaxial and the measurement location of the Active Target is independent of the spindle rotation.
Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal

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.

Fig. 9
Fig. 9
Close modal
Table 1

Minimum and maximum axis commands for collected 295 measurements

AxisMinimumMaximum
X (mm)83.26081.1
Y (mm)37.22557.1
Z (mm)7.4988.2
C (deg)−269.7269.9
B (deg)−109.5109.9
AxisMinimumMaximum
X (mm)83.26081.1
Y (mm)37.22557.1
Z (mm)7.4988.2
C (deg)−269.7269.9
B (deg)−109.5109.9

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.

Fig. 10
Fig. 10
Close modal

4.2 Unconstrained Model.

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.

Fig. 11
Fig. 11
Close modal
Table 2

Mean and maximum residuals for the identification and validation data sets (mm)

Identification setInterpolation validation setExtrapolation validation set
MeanMaxMeanMaxMeanMax
Nominal machine0.3070.7390.2610.5390.3210.525
Unconstrained model0.0380.1400.0440.0920.4511.909
Identification setInterpolation validation setExtrapolation validation set
MeanMaxMeanMaxMeanMax
Nominal machine0.3070.7390.2610.5390.3210.525
Unconstrained model0.0380.1400.0440.0920.4511.909

4.3 Constrained 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.

Fig. 12
Fig. 12
Close modal

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.

Fig. 13
Fig. 13
Close modal

4.4 Comparison and Analysis.

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.

Fig. 14
Fig. 14
Close modal
Fig. 15
Fig. 15
Close modal
Fig. 16
Fig. 16
Close modal
Table 3

Mean and maximum residuals of nominal machine and the unconstrained and constrained models for identification and validation sets (mm)

Identification setInterpolation validation setExtrapolation validation set
MeanMaxMeanMaxMeanMax
Nominal machine0.3070.7390.2610.5390.3210.525
Unconstrained model0.0380.1400.0440.0920.4511.909
Constrained model0.0540.2170.0560.1120.1910.443
Identification setInterpolation validation setExtrapolation validation set
MeanMaxMeanMaxMeanMax
Nominal machine0.3070.7390.2610.5390.3210.525
Unconstrained model0.0380.1400.0440.0920.4511.909
Constrained model0.0540.2170.0560.1120.1910.443

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).

5 Summary and Conclusions

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.

Acknowledgment

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.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

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.

Nomenclature

c =

constraint value for constrained model construction

b =

vector of error parameters to be identified

q =

nominal axis command vector

S =

skew symmetric matrix

lt =

measured tool length

$q¯i$ =

ith axis command scaled to range [−1 1]

et =

tool tip modeled error

pa =

actual measurement

Fnom =

nominal kinematics

FAP =

axis perturbation model kinematics

Jv =

linear velocity portion of machine tool Jacobian matrix

Jω =

angular velocity portion of machine tool Jacobian matrix

Ti =

linear homogeneous transformation for ith axis

$h¯i$ =

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

Eij(qj) =

axis perturbation function mapping axis command qj onto an adjustment for axis command qi

δ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

Appendix A: Linear Velocity Portion of Model Jacobian Matrix

From Ref. [37], the linear velocity portion of the model Jacobian matrix is
$Jv,AP=∇qpAP$
(A1)
For an industrial five-axis machine tool with axis sequence XYZCB, the modeled position of the machine tool’s last frame with respect to the machine tool’s base frame is
$pAP=[qX+δqX(q)qY+δqY(q)qZ+δqZ(q)]T$
(A2)
Thus, using Eq. (A1), the linear spatial rate portion of the model Jacobian matrix is
$Jv,AP=[∂pAP∂qX∂pAP∂qY∂pAP∂qZ∂pAP∂qC∂pAP∂qB]=[1+dEXX(qX)dqXdEXY(qY)dqYdEXZ(qZ)dqZdEXC(qC)dqCdEXB(qB)dqBdEYX(qX)dqX1+dEYY(qY)dqYdEYZ(qZ)dqZdEYC(qC)dqCdEYB(qB)dqBdEZX(qX)dqXdEZY(qY)dqY1+dEZZ(qZ)dqZdEZC(qC)dqCdEZB(qB)dqB]$
(A3)

Appendix B: Angular Velocity Portion of Model Jacobian Matrix

Continuing to consider a five-axis machine tool with the sequence XYZCB, the modeled rotation transformation from the machine tool’s base frame to the machine tool’s last frame is
$RAP=RCRB$
(B1)
where RC and RB are the modeled rotation transformations for the C and B axes, respectively
$RC=[cos(qC+δqC(q))−sin(qC+δqC(q))0sin(qC+δqC(q))cos(qC+δqC(q))0001]$
(B2)
$RB=[cos(qB+δqB(q))0sin(qB+δqB(q))010−sin(qB+δqB(q))0cos(qB+δqB(q))]$
(B3)
Taking the time derivative of RAP
$R˙AP=dRCdtRB+RCdRBdt=∑i=X,Y,Z,C,B(q˙i(qC+dECi(qi)dqi)S(k)RCRB+q˙i(qB+dEBi(qi)dqi)RCS(j)RB)=∑i=X,Y,Z,C,B(S(q˙i(qC+dECi(qi)dqi)k)R+S(q˙i(qB+dEBi(qi)dqi)RCj)R)=S(∑i=X,Y,Z,C,B(q˙i(qC+dECi(qi)dqi)k+q˙i(qB+dEBi(qi)dqi)RCj))R$
(B4)
where
$k=[001]T,j=[010]T$
(B5)
Thus, the modeled angular velocity is
$ωAP=∑i=X,Y,Z,C,BωAP,i=∑i=X,Y,Z,C,B(q˙i(qC+dECi(qi)dqi)k+q˙i(qB+dEBi(qi)dqi)RCj)$
(B6)
where
$ωAP,i=[−sin(qC+δqC(q))(dEBi(qi)dqi)cos(qC+δqC(q))(dEBi(qi)dqi)dECi(qi)dqi]q˙i,i=X,Y,Z$
(B7)
$ωAP,C=[−sin(qC+δqC(q))(dEBC(qC)dqC)cos(qC+δqC(q))(dEBC(qC)dqC)1+dECC(qC)dqC]q˙C$
(B8)
$ωAP,B=[−sin(qC+δqC(q))(1+dEBB(qB)dqB)cos(qC+δqC(q))(1+dEBB(qB)dqB)dECB(qB)dqB]q˙B$
(B9)
The angular velocity portion of the model Jacobian matrix is thus
$Jω,AP=[ωAP,Xq˙XωAP,Yq˙YωAP,Zq˙ZωAP,Cq˙CωAP,Bq˙B]$
(B10)

Appendix C: Linear and Angular Velocity Portion of Nominal Jacobian Matrix

The linear and angular velocity portions of the nominal Jacobian matrix are given when the perturbation functions and their slopes are zero in Eqs. (A3) and (B10)
$Jv,nom=[100000100000100]$
(C1)
$Jω,nom=[000000000001−sin(qC)cos(qC)0]$
(C2)

References

1.
Gu
,
J.
,
Agapiou
,
J. S.
, and
Kurgin
,
S.
,
2015
, “
CNC Machine Tool Work Offset Error Compensation Method
,”
J. Manuf. Syst.
,
37
, pp.
576
585
.
2.
Creamer
,
J.
,
Sammons
,
P. M.
,
Bristow
,
D. A.
,
Landers
,
R. G.
,
Freeman
,
P. L.
, and
Easley
,
S. J.
,
2017
, “
Table-Based Volumetric Error Compensation of Large Five-Axis Machine Tools
,”
ASME J. Manuf. Sci. Eng.
,
139
(
2
), p.
021011
.
3.
Calleja
,
A.
,
Bo
,
P.
,
González
,
H.
,
Bartoň
,
M.
, and
López de Lacalle
,
L. N.
,
2018
, “
Highly Accurate 5-Axis Flank CNC Machining with Conical Tools
,”
,
97
(
5
), pp.
1605
1615
.
4.
Lamikiz
,
A.
,
López de Lacalle
,
L. N.
,
Ocerin
,
O.
,
Díez
,
D.
, and
Maidagan
,
E.
,
2008
, “
The Denavit and Hartenberg Approach Applied to Evaluate the Consequences in the Tool Tip Position of Geometrical Errors in Five-Axis Milling Centres
,”
,
37
(
1–2
), pp.
122
139
.
5.
Donmez
,
M. A.
,
Blomquist
,
D. S.
,
Hocken
,
R. J.
,
Liu
,
C. R.
, and
Barash
,
M. M.
,
1986
, “
A General Methodology for Machine Tool Accuracy Enhancement by Error Compensation
,”
Precis. Eng.
,
8
(
4
), pp.
187
196
.
6.
Li
,
C.
,
Liu
,
X.
,
Li
,
R.
,
Wu
,
S.
, and
Song
,
H.
,
2020
, “
Geometric Error Identification and Analysis of Rotary Axes on Five-Axis Machine Tool Based on Precision Balls
,”
Appl. Sci.
,
10
(
1
), p.
100
.
7.
Hashemiboroujeni
,
H.
,
Marzdashti
,
S. E.
,
Xing
,
K.
, and
Mayer
,
J. R. R.
,
2019
, “
Five-Axis Machine Tool Coordinate Metrology Evaluation Using the Ball Dome Artefact Before and After Machine Calibration
,”
J. Manuf. Mater. Process.
,
3
(
1
), p.
20
.
8.
Lee
,
W.-C.
,
Lee
,
Y. T.
, and
Wei
,
C. C.
,
2019
, “
,”
Appl. Sci.
,
9
(
15
), p.
3073
.
9.
Altintas
,
Y.
,
Brecher
,
C.
,
Weck
,
M.
, and
Witt
,
S.
,
2005
, “
Virtual Machine Tool
,”
Ann. CIRP
,
54
(
2
), pp.
115
138
.
10.
Bringmann
,
B.
,
Besuchet
,
J. P.
, and
Rohr
,
L.
,
1994
, “
Systematic Evaluation of Calibration Methods
,”
CIRP Ann.
,
57
(
1
), pp.
529
532
.
11.
ISO
,
2012
, “
Test Code for Machine Tools Part I: Geometric Accuracy of Machine Tools Operating Under No-Load or Quasi-Static Conditions
,”
International Organization for Standardization
,
Geneva, Switzerland
, Standard No. ISO 230-1.
12.
Schwenke
,
H.
,
Knapp
,
W.
,
Haitjema
,
H.
,
Weckenmann
,
A.
,
Schmitt
,
R.
, and
Delbressine
,
F.
,
2008
, “
Geometric Error Measurement and Compensation of Machines—An Update
,”
CIRP Ann.
,
57
(
2
), pp.
660
675
.
13.
Tsutsumi
,
M.
, and
Saito
,
A.
,
2003
, “
Identification and Compensation of Systematic Deviations Particular to 5-Axis Machining Centers
,”
Int. J. Mach. Tools Manuf.
,
43
(
8
), pp.
771
780
.
14.
Lei
,
W. T.
,
Sung
,
M. P.
,
Liu
,
W. L.
, and
Chuang
,
Y. C.
,
2007
, “
Double Ballbar Test for the Rotary Axes of Five-Axis CNC Machine Tools
,”
Int. J. Mach. Tools Manuf.
,
47
(
2
), pp.
273
285
.
15.
Ibaraki
,
S.
,
Oyama
,
C.
, and
Otsubo
,
H.
,
2011
, “
Construction of an Error Map of Rotary Axes on a Five-Axis Machining Center by Static R-Test
,”
Int. J. Mach. Tools Manuf.
,
51
(
3
), pp.
190
200
.
16.
Bringmann
,
B.
, and
Knapp
,
W.
,
2006
, “
Model-Based ‘Chase-the-Ball’ Calibration of a 5-Axes Machining Center
,”
CIRP Ann.
,
55
(
1
), pp.
531
534
.
17.
Erkan
,
T.
, and
Mayer
,
J. R. R.
,
2010
, “
A Cluster Analysis Applied to Volumetric Errors of 5-Axis Machine Tools Obtained by Probing an Uncalibrated Artifact
,”
CIRP Ann.
,
59
(
1
), pp.
539
542
.
18.
Ibaraki
,
S.
,
Iritani
,
T.
, and
Matsushita
,
T.
,
2012
, “
Calibration of Location Errors of Rotary Axes on 5-Axis Machine Tools by On-the-Machine Measurement Using a Touch-Trigger Probe
,”
Int. J. Mach. Tools Manuf.
,
58
(
1
), pp.
44
53
.
19.
,
S.
,
Samper
,
D.
,
Santolaria
,
J.
, and
Aguilar
,
J. J.
,
2012
, “
Identification Strategy of Error Parameter in Volumetric Error Compensation of Machine Tool Based on Laser Tracker Measurements
,”
Int. J. Mach. Tools Manuf.
,
53
(
1
), pp.
160
169
.
20.
Wang
,
J.
,
Guo
,
J.
,
Zhang
,
G.
,
Guo
,
B. A.
, and
Wang
,
H.
,
2012
, “
The Technical Method of Geometric Error Measurement for Multi-Axis NC Machine Tool by Laser Tracker
,”
Meas. Sci. Technol.
,
23
(
4
), p.
045003
.
21.
Hong
,
C.
,
Ibaraki
,
S.
, and
Matsubara
,
A.
,
2001
, “
Influence of Position Dependent Error of Rotary Axes on a Machining Test of Cone Frustram by 5-Axis Machine Tools
,”
Precis. Eng.
,
35
(
1
), pp.
1
11
.
22.
Ibaraki
,
S.
, and
Knapp
,
W.
,
2012
, “
Indirect Measurement of Volumetric Accuracy for 3-Axis and 5-Axis Machine Tools: A Review
,”
Int. J. Autom. Technol.
,
6
(
2
), pp.
110
124
.
23.
Choi
,
J. P.
,
Min
,
B. K.
, and
Lee
,
S. J.
,
2004
, “
Reduction of Machining Errors of a Three-Axis Machine Tool by On-Machine Measurement and Error Compensation System
,”
J. Mater. Process. Technol.
,
155
, pp.
2056
2064
.
24.
Xiang
,
S.
, and
Altintas
,
Y.
,
2016
, “
Modeling and Compensation of Volumetric Errors for Five-Axis Machine Tools
,”
Int. J. Mach. Tools Manuf.
,
101
, pp.
65
78
.
25.
Potgieter
,
G.
, and
Engelbrecht
,
A. P.
,
2007
, “
Genetic Algorithms for the Structural Optimisation of Learned Polynomial Expressions
,”
Appl. Math. Comput.
,
186
(
2
), pp.
1441
1466
.
26.
Aaen
,
P.
,
Plá
,
J. A.
, and
Wood
,
J.
,
2007
,
Modeling and Characterization of RF and Microwave Power FETs
,
Cambridge University Press
,
Cambridge, UK
.
27.
Runge
,
C.
,
1901
, “
Über Empirische Funktionen und die Interpolation Zwischen Äquidistanten Ordinaten
,”
Z. Angew. Math. Phys.
,
46
, pp.
224
243
.
28.
Donaldson
,
R. R.
,
1980
, “
Error Budgets
,”
Technol. Mach. Tools
,
5
, pp.
1
14
.
29.
Shen
,
Y. L.
, and
Duffie
,
N. A.
,
1993
, “
Comparison of Combinatorial Rules for Machine Error Budgets
,”
CIRP Ann.
,
42
(
1
), pp.
619
622
.
30.
Cheng
,
Q.
,
Dong
,
L.
,
Liu
,
Z.
,
Li
,
J.
, and
Gu
,
P.
,
2018
, “
A New Geometric Error Budget Method of Multi-Axis Machine Tool Based on Improved Value Analysis
,”
Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci.
,
232
(
22
), pp.
4064
4083
.
31.
Dorndorf
,
U.
,
Kiridena
,
V. S. B.
, and
Ferreira
,
P. M.
,
1994
, “
Optimal Budgeting of Quasistatic Machine Tool Errors
,”
ASME J. Eng. Ind.
,
116
(
1
), pp.
42
53
.
32.
Zhang
,
Z.
,
Cai
,
L.
,
Cheng
,
Q.
,
Liu
,
Z.
, and
Gu
,
P.
,
2019
, “
A Geometric Error Budget Method to Improve Machining Accuracy Reliability of Multi-Axis Machine Tools
,”
J. Intell. Manuf.
,
30
(
2
), pp.
495
519
.
33.
Eisenbies
,
S. K.
,
2001
, “
Error Budget by Constraints
,” Master’s thesis,
University of North Carolina
,
Charlotte
.
34.
Gupta
,
K.
,
1986
, “
Kinematic Analysis of Manipulators Using the Zero Reference Position Description
,”
Int. J. Rob. Res.
,
5
(
2
), pp.
5
13
.
35.
Denavit
,
J.
, and
Hartenberg
,
R. S.
,
1955
, “
A Kinematic Notation for Lower-Pair Mechanisms Based on Matrices
,”
ASME J. Appl. Mech.
,
22
(
2
), pp.
215
221
.
36.
Wampler
,
C. W.
,
Hollerbach
,
J. M.
, and
Arai
,
T.
,
1995
, “
An Implicit Loop Method for Kinematic Calibration and Its Application to Closed-Chain Mechanisms
,”
IEEE Trans. Rob. Autom.
,
11
(
5
), pp.
710
724
.
37.
Ma
,
L.
,
Bazzoli
,
P.
,
Sammons
,
P. M.
,
Landers
,
R. G.
, and
Bristow
,
D. A.
,
2018
, “
Modeling and Calibration of High-Order Joint-Dependent Kinematic Errors for Industrial Robots
,”
Rob. Comput. Integr. Manuf.
,
50
, pp.
153
167
.
38.
Wenger
,
P.
, and
Chablat
,
D.
,
2000
, “
Kinematic Analysis of a New Parallel Machine Tool: the Orthoglide
,”
, pp.
305
314
.
39.
Xia
,
C.
,
Hsu
,
W.
,
Lee
,
M. L.
, and
Ooi
,
B. C.
,
2006
, “
BORDER: Efficient Computation of Boundary Points
,”
IEEE Trans. Knowl. Data Eng.
,
18
(
3
), pp.
289
303
.
40.
Niederreiter
,
H.
,
1988
, “
Low-Discrepancy and Low-Dispersion Sequences
,”
J. Number Theory
,
30
(
1
), pp.
51
70
.