Numerical Aspects of Multilateration for Volumetric Error Calculation

Citation: Holub, M.; Zatoˇcilová, J.; Marek, T.; Blecha, P.; Heinrich, P. Numerical Aspects of Multilateration for Volumetric Error Calculation. Machines 2022,10, 833. https:// doi.org/10.3390/machines10100833 Academic Editor: Feng Gao Received: 1 August 2022 Accepted: 15 September 2022 Published: 21 September 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). machines Article Numerical Aspects of Multilateration for Volumetric Error Calculation Michal Holub 1, Jitka Zatoˇcilová 1, Tomáš Marek 1, Petr Blecha 1,* and Petr Heinrich 2 1Faculty of Mechanical Engineering, Brno University of Technology, 616 69 Brno, Czech Republic 2KOVOSVIT MAS Machine Tools, 391 02 Sezimovo Ústí, Czech Republic *Correspondence: [email protected] Abstract: We discuss two approaches to multilateration with a particular focus on numerical aspects for a given dataset. More precisely, we calculate the volumetric errors of the MCV 754 milling machine from the data measured by a LaserTRACER and provide criteria showing which numerical method is appropriate for the solution of the multilateration problem with respect to a given measurement dataset. In the first case, we consider two cost functions of the optimisation problem; in the second, we use the extremal principle method. We discuss the calculation accuracy depending on the matrix condition number. Our results are compared to the reference results produced by the Trac-Cal software, which is a standard used by most producers for error compensations. Keywords: multilateration; numerical model; volumetric errors; machine tool error 1. Introduction In this paper, we discuss the numerical aspects of methods for multilateration and the calculation of the volumetric errors of a machine tool (MT). More precisely, we compare two optimisation methods in several modifications regarding their accuracy and stability. We claim that some authors are using mathematical methods improperly, meaning that their method is sensitive to the input dataset and, therefore, unsuitable for statistical refinement, see [ 1 , 2 ]. In general, for multilateration problems, various numerical methods may be used, see [ 3 ] for an overview with the description of advantages and disadvantages. In our paper, we consider two particular methods and compare them. We compare the standard optimisation by the trust region method with two different cost functions and an extremal principle solution. Although the methods are known, we specify the better version with respect to the given dataset, more precisely with respect to the condition number of the appropriate system matrix. We stress that the choice of a correct numerical method is vital because the results of the multilateration problem are used for volumetric error calculation, determining the volumetric accuracy model and, consequently, the geometric error calculation [ 4 ]. We are aware that the proportion of quasi-static errors in three-axis machine tools is 60–70% of the total working accuracy. In [ 5 , 6 ], the authors have extended this estimate further in tests carried out on five-axis machining centres, where they claim that the proportion of quasi-static errors is even higher, up to 80%. Quasi-static errors are defined as errors in the relative position of the tool centre point and the workpiece, with the errors changing slowly over time. These errors are directly related to the structure of the machine tool itself and can be broken down into geometric, kinematic and thermal errors. For the measurement itself, we use the LaserTRACER [ 7 – 11 ], but the data may be obtained by other tools such as the Ballbar-type apparatus [ 12 ], laser tracker [ 13 , 14 ] or laser interferometer [ 15 – 17 ]. We also describe the way of obtaining an appropriate kinematic chain because various versions for MTs are known, see, e.g., [16,18–21]. The paper is structured as follows. First, we present a verification of the MT kinematic chain in Section 2. Consequently, in Section 3.1, we consider two optimisation approaches to multilateration and volumetric error calculation. In the first approach, referred to as Case Machines 2022,10, 833. https://doi.org/10.3390/machines10100833 https://www.mdpi.com/journal/machines Machines 2022,10, 833 2 of 15 A, we compare two different cost functions for the trust region method and discuss their numerical properties. The second approach, Case B, uses a Lagrange multipliers optimisation method. Once all the methods are applied to the data measured by the LaserTRACER, they are verified and compared with the results obtained from the Trac-Cal software. This is a software tool widely used by producers of MTs and CMMs (coordinate measuring machines) for geometric error compensations and, thus, we consider it a good benchmark for comparisons. When appropriate numerical methods are applied, faster computation of geometric errors may be achieved and the uncertainty of the whole measuring system may drop. Therefore, the software for error compensations of CNC machine tools may be improved. The influence of such software compensations on the resulting shape accuracy is described in [ 22 – 24 ]. Finally, in Section 4, we give our recommendations for the methods studied. 2. Three-Axis Machine Verification We use the demonstrator MCV 754 QUICK (Figure 1left) for the measurement in order to validate and verify our results. The working space (WS) of the demonstrator is defined by travels of the individual axes with the rank of 754 × 500 × 550 mm, see Figure 1, right. The machine specifications are given in Table 1. Z = 550 mm Y = 500 mm Z-axis X-axis Y-axis Figure 1. Demonstrator MCV 754 QUICK. Table 1. CNC machine tools technical parameters. Item Travel of Xaxis/Measurement range of Xaxis 754/750 mm Travel of Yaxis/Measurement range of Yaxis 500/500 mm Travel of Zaxis/Measurement range of Zaxis 550/550 mm Bi-directional systematic positioning error of an axis ISO 230-2:2014 0.008 mm Temperature sensors uncertainty, Pt100, class A 0.25 °C In Table 2, we provide the fundamental technical specification of the LaserTRACER measurement system. Note that measurement uncertainty is one of the crucial parameters that is affected by a bundle of factors, among others by the calculation of deviations. Table 2. LaserTRACER specifications. Item Resolution 0.001 µm Measuring uncertainty (k=2) (0.2 + 0.3 ×L) µm Measuring range 200–18,000 mm Machines 2022,10, 833 3 of 15 Figure 2shows the error vector map of the MCV working space. The maximal deviation within the measured area is 64 µ m. The magnitude of the errors is tripled in order to visualise the error vector better. The multilateration-based calculation affects the resulting magnitude and direction of the error vectors. These are consequently imported to the CNC milling machine as error corrections. Figure 2. Error vector map (TRAC−CAL), the axes’ units are mm. Note that the Cartesian coordinate system X , Y , Z is assigned to the WS in the kinematic chain WXYZT (workpiece tool, W—workpiece, T—tool). In order to recall the notation for geometric errors used in the kinematic chain, we refer to Tables 3and 4. Each translational axis has six geometrical deviations. According to ISO 230-1 [ 25 ], the deviations of the X -axis are designated as EXX, EYX, EZX, EAX, EBX and ECX. Note also that the squareness deviations between two axes are mentioned, e.g., between X and Y —C0Y. This makes a total of 6 × 3 errors for the translational axes and three squareness errors of the appropriate axis pairs. Table 3. Geometric errors for three-axis MT. Axis Axis Error (ISO/Paper Symbol) X-axis EXX/δxx EYX/δyx EZX/δzx EAX/εxx EBX/εyx ECX/εzx Y-axis EXY/δxy EYY/δyy EZY/δzy EAY/εxy EBY/εyy ECY/εzy Z-axis EXZ/δxz EYZ/δyz EZZ/δzz EAZ/εxz EBZ/εyz ECZ/εzz Table 4. Squareness errors for three-axis MT. Axis Squareness Error (ISO/Paper Symbol) X-axis B0Z/Sxy C0Y/Sxz Y-axis A0Z/Syz Z-axis C0Y/Sxz The measured data were obtained on the test demonstrator with the “tracking laser interferometer” (LaserTRACER), Figure 3. Machines 2022,10, 833 4 of 15 Spindle Reflector LaserTRACER Temperature sensor Z-axis under machine covers Worktable Temperature sensor X-axis under machine covers Ambient sensor Figure 3. Configurations of tracking the laser interferometer. If we denote X=(x,y,z)T a vector of nominal coordinates, Xa=(xa,ya,za)T a vector of actual coordinates and Xp=xp,yp,zpT the probe offset coordinates, their relation may be expressed by the kinematic chain   xa ya za =Tx Ty Tz  xp yp zp +Lz +Ly +Lx, (1) where Tx , Ty , Tz are modified rotational matrices and Lx , Ly , Lz are translational vectors in the form Tx=  1−εzxι εyxι εzxι1−εxxι −εyxι εxxι1 ,Ty=  1−εzyι εyyι εzyι1−εxyι −εyyι εxyι1 , Tz=  1−εzzι εyzι εzzι1−εxzι −εyzι εxzι1 , Lx=  x+δxxι δyxι δzxι ,Ly=  δxyι−Sxyy y+δyyι δzyι ,Lz=  δxzι−Sxzz δyzι−Syzz z+δzzι , where we use the notation from Table 3for the geometric errors and ι is an auxiliary variable such that ι2=0, also known as a dual number unit, see [26,27]. Therefore, the volumetric error (∆x,∆y,∆z)Tcan be calculated from (1) as ∆x=xa−(x+xp) = δxx +δxy +δxz +y(−εzx −Sxy) + z(Sxz +εyx +εyy) +yp(−εzx −εzy −εzz) + zp(εyx +εyy +εyz), (2) ∆y=ya−(y+yp) = δyx +δyy +δyz +z(−Syz −εxx −εxy) +xp(εzx +εzy +εzz) + zp(−εxx −εxy −εxz), (3) ∆z=za−(z+zp) = δzx +δzy +δzz +yεxx +xp(−εyx −εyy −εyz) +yp(εxx +εxy +εxz)(4) Machines 2022,10, 833 5 of 15 after neglecting the terms with ι2 . Note that we chose the notation of [ 26 , 27 ] with dual number unit ιto neglect the terms with ι2algebraically, i.e., directly in the computation. The correctness of this kinematic chain has been verified on the Trac-Cal data. Indeed, we substitute the rotational and translational errors from Trac-Cal to the kinematic chain (1) and calculate the volumetric errors. Consequently, we consider the error evolution in the directions of the axes (red data in Figure 4) and no evolution (green data in Figure 4). By error evolution, we understand that, for instance, if we move on the y -axis, the ∆x calculation is affected by the change in the y coordinates, i.e., to the kinematic chains (2) – (4) we, respectively, substitute zeros for x and z , and for y we substitute the position on the y axis, while in the case with no evolution we also set y= 0. Thus, squareness and rotational errors only are added to ∆x . In this sense, when calculating ∆x on the x -axis, the remaining coordinates are 0 and, therefore, the terms containing y and z vanish in the kinematic chain (2) – (4) . We compare these errors to the Trac-Cal volumetric errors (blue isolated points in Figure 4). We conclude that it is necessary to take into account the evolution of the error with respect to the shift in the direction of the axis. Otherwise, the deviation is unacceptable. 0 200 400 600 800 x coordinate 0 0.01 0.02 0.03 x 0 200 400 600 800 x coordinate -2 -1 0 1 2 y 10-3 0 200 400 600 800 x coordinate -15 -10 -5 0 5 z 10-4 0 100 200 300 400 500 y coordinate -15 -10 -5 0 x 10-3 0 100 200 300 400 500 y coordinate -0.02 -0.015 -0.01 -0.005 0 y 0 100 200 300 400 500 y coordinate -2 0 2 4 z 10-3 -600 -400 -200 0 z coordinate -0.02 0 0.02 0.04 0.06 x -600 -400 -200 0 z coordinate 0 0.005 0.01 0.015 0.02 y -600 -400 -200 0 z coordinate -2 0 2 4 6 8 z 10-3 Figure 4. Difference between the calculated volumetric errors without evolution (green) and those with evolution (red) compared with the Trac-Cal data (blue). Units: mm. 3. Numerical Methods for Multilateration The general principle of multilateration lies in determining a point’s precise position from multiple distance measurements. As a measuring device, we use LaserTRACER (LT). Given the distances between the ideal (also referred to as nominal) positions and the LaserTRACER positions, we calculate the actual (inaccurate) positions of the tool. The accuracy of the location also depends on the number of LTs and their location [ 28 ]. We use four or six LTs. As the first step, we determine the exact positions of LTs and, consequently, we measure the actual positions. Machines 2022,10, 833 6 of 15 3.1. LTs’ Positions Calculation Now, we calculate the LaserTRACERs’ positions. Let us denote the nominal positions as Xi= (xi , yi , zi) for i= 1, . . . , N , where N is the number of the measured points and LT(k)= (x(k) LT , y(k) LT , z(k) LT ) with k= 1, . . . , NLT , the LaserTRACER positions, where NLT stands for the number of LTs. The relative distance between the LaserTRACER position LT(k) and the measured points Xican be expressed in two ways: 1. In the form of the following system of non-linear equations rx1−x(k) LT 2+y1−y(k) LT 2+z1−z(k) LT 2=Lk+P(k) 1 rx2−x(k) LT 2+y2−y(k) LT 2+z2−z(k) LT 2=Lk+P(k) 2 . . . (5) rxN−x(k) LT 2+yN−y(k) LT 2+zN−z(k) LT 2=Lk+P(k) N, where k= 1, . . . , NLT , Lk is the residual error of the k -th LT and P(k) i is the measured distance between LT(k)and the nominal positions Xi,i=1, . . . , N. 2. If particular non-linear equations of (5) are squared, we receive a system x1−x(k) LT 2+y1−y(k) LT 2+z1−z(k) LT 2=Lk+P(k) 12 x2−x(k) LT 2+y2−y(k) LT 2+z2−z(k) LT 2=Lk+P(k) 22 . . . (6) xN−x(k) LT 2+yN−y(k) LT 2+zN−z(k) LT 2=Lk+P(k) N2. Consequently, the positions of a particular LT(k) are determined from system (5) or (6) , respectively, by the least square minimisation. Yet, even this optimisation technique may be treated in various ways. In the sequel, we discuss the numerical properties and the optimisation methods applied to both systems (5) and (6). 3.1.1. Case A In this case, we transform system (5) into a non-linear function of four unknowns Fk(x(k) LT ,y(k) LT ,z(k) LT ,Lk) = N ∑ i=1 rxi−x(k) LT 2+yi−y(k) LT 2+zi−z(k) LT 2−Lk+P(k) i!(7) and apply minimisation by the trust-region method with respect to particular positions LT(k) with k= 1, . . . , NLT and Lk (in the sense of the least squares). Indeed, we are searching for min x(k) LT ,y(k) LT ,z(k) LT ,Lk kFk(x(k) LT ,y(k) LT ,z(k) LT ,Lk)k2 2, where k·k2 is a Euclidean norm. The positions calculated by this method are consequently compared with those acquired from the Trac-Cal software and the results are summarised in Tables 5and 6. Clearly, the difference between the computed values and the reference values is negligible for this numerical model. To justify this conclusion, let us define that the difference is considered tolerable if it falls into the statistical deviation interval of the Trac-Cal data, which is 5% for our measurements. Therefore, our deviations are negligible. Machines 2022,10, 833 7 of 15 Table 5. Calculated LT positions for Case A with function (7). Calculated Values (mm) Trac-Cal (mm) LT(1)[−154.799406, 125.978808, −471.460113] [−154.799400, 125.978800, −471.460100] LT(2)[−143.324882, 328.559001, −471.386505] [−143.324800, 328.559000, −471.386500] LT(3)[862.105512, 291.242899, −340.459834] [862.105400, 291.242700, −340.460000] LT(4)[833.3820, 103.1588, −340.4466] [833.382200, 103.158700, −340.446600] LT(5)[891.307748, 152.634701, −507.868235] [891.308200, 152.634600, −507.868500] LT(6)[885.579601, 154.917244, −507.865002] [885.581100, 154.917200, −507.865700] Table 6. Difference values from Table 5. Difference (mm) LT(1)[0.000006, 0.000008, 0.000013] LT(2)[0.000082, 0.000001, 0.000005] LT(3)[0.000112, 0.000199, 0.000166] LT(4)[0.000220, 0.000129, 0.000001] LT(5)[0.000452, 0.000101, 0.000265] LT(6)[0.001499, 0.000044, 0.000698] If the same approach is used for the non-linear system (6) , we apply the trust-region method, i.e., we are searching for a minimum (in the sense of the least square method) min x(k) LT ,y(k) LT ,z(k) LT ,Lk kFk(x(k) LT ,y(k) LT ,z(k) LT ,Lk)k2 2 of a non-linear function of four unknowns Fk(x(k) LT ,y(k) LT ,z(k) LT ,Lk) = N ∑ i=1xi−x(k) LT 2+yi−y(k) LT 2+zi−z(k) LT 2−Lk+P(k) i2. (8) The calculated positions of LTs are again compared with the reference data acquired from the Trac-Cal software and displayed in Table 7. Table 7. Calculated LT positions for Case A with function (8). Calculated Values (mm) Difference (mm) LT(1)[−154.794962, 125.982784, −471.450112] [0.004438, 0.003984, 0.009988] LT(2)[−143.318839, 328.556174, −471.377919] [0.005961, 0.002826, 0.008581] LT(3)[862.109478, 291.237421, −340.474264] [0.004078, 0.005279, 0.014264] LT(4)[833.386275, 103.153101, −340.458882] [0.004075, 0.005599, 0.012282] LT(5)[891.340682, 152.629665, −507.882949] [0.032482, 0.004935, 0.014449] LT(6)[885.615803, 154.916576, −507.878724] [0.034703, 0.000624, 0.013024] Even though the deviation may still be acceptable, clearly the results are worse than in the first case. Let us note that this may be partially caused even by the computer representation of decimal numbers. Indeed, for a large number, more digits after the decimal points may be neglected in the representation. Clearly, function (8) reaches higher values by definition. This effect may be avoided by a change in the data type but this may increase the computational time enormously. Machines 2022,10, 833 8 of 15 3.1.2. Case B In the second case, see e.g., ref. [ 1 ], the non-linear system (6) is rewritten into the following form x2 1+y2 1+z2 1−2x1x(k) LT −2y1y(k) LT −2z1z(k) LT −2LkP(k) 1−(P(k) 1)2+Ck=0 x2 2+y2 2+z2 2−2x2x(k) LT −2y2y(k) LT −2z2z(k) LT −2LkP(k) 2−(P(k) 2)2+Ck=0 . . . x2 N+y2 N+z2 N−2xNx(k) LT −2yNy(k) LT −2zNz(k) LT −2LkP(k) N−(P(k) N)2+Ck=0, where Ck= (x(k) LT )2+ (y(k) LT )2+ (z(k) LT )2−(Lk)2, and we are searching for a minimum (in the sense of the least squares) of a function Fk(x(k) LT ,y(k) LT ,z(k) LT ,Lk,Ck) = N ∑ i=1x2 i+y2 i+z2 i−2xix(k) LT −2yiy(k) LT −2ziz(k) LT −2LkP(k) i−(P(k) i)2+Ck2. (9) According to the extremal principle, the first-order partial derivative of the objective function should be equal to zero, while the second-order partial derivative should be greater than zero. The first-order partial derivatives are expressed as ∂Fk ∂x(k) LT = N ∑ i=1 4xi−x2 i−y2 i−z2 i+ (P(k) i)2+2xix(k) LT +2yiy(k) LT +2ziz(k) LT +2LkP(k) i−Ck ∂Fk ∂y(k) LT = N ∑ i=1 4yi−x2 i−y2 i−z2 i+ (P(k) i)2+2xix(k) LT +2yiy(k) LT +2ziz(k) LT +2LkP(k) i−Ck ∂Fk ∂z(k) LT = N ∑ i=1 4zi−x2 i−y2 i−z2 i+ (P(k) i)2+2xix(k) LT +2yiy(k) LT +2ziz(k) LT +2LkP(k) i−Ck(10) ∂Fk ∂Lk = N ∑ i=1 4Lk−x2 i−y2 i−z2 i+ (P(k) i)2+2xix(k) LT +2yiy(k) LT +2ziz(k) LT +2LkP(k) i−Ck ∂Fk ∂Ck = N ∑ i=1 2x2 i+y2 i+z2 i−(P(k) i)2−2xix(k) LT −2yiy(k) LT −2ziz(k) LT −2LkP(k) i+Ck It is found that the second-order partial derivative is always positive, and thus the second condition is satisfied. For unknowns (x(k) LT , y(k) LT , z(k) LT , Lk , Ck) , k= 1,..., NLT , the system of five equations (10) can be rewritten in the matrix form as A·x(k) LT ,y(k) LT ,z(k) LT ,Lk,CkT=b, (11) where the 5 ×5 matrix Ais of the form A=                  2N ∑ i=1x2 i2N ∑ i=1xiyi2N ∑ i=1xizi2N ∑ i=1xiP(k) i−N ∑ i=1xi 2N ∑ i=1xiyi2N ∑ i=1y2 i2N ∑ i=1yizi2N ∑ i=1yiP(k) i−N ∑ i=1yi 2N ∑ i=1xizi2N ∑ i=1yizi2N ∑ i=1z2 i2N ∑ i=1ziP(k) i−N ∑ i=1zi 2N ∑ i=1xiP(k) i2N ∑ i=1yiP(k) i2N ∑ i=1ziP(k) i2N ∑ i=1(P(k) i)2−N ∑ i=1P(k) i −N ∑ i=1xi−N ∑ i=1yi−N ∑ i=1zi−N ∑ i=1P(k) iN 2                  Machines 2022,10, 833 9 of 15 and the right-hand side vector bof dimension 5 is b=                  N ∑ i=1xix2 i+y2 i+z2 i−(P(k) i)2 N ∑ i=1yix2 i+y2 i+z2 i−(P(k) i)2 N ∑ i=1zix2 i+y2 i+z2 i−(P(k) i)2 N ∑ i=1P(k) ix2 i+y2 i+z2 i−(P(k) i)2 N ∑ i=1−1 2x2 i+y2 i+z2 i−(P(k) i)2                  . The system matrix A is always regular but if its condition number κ(A) is calculated, it is clear that A is ill-conditioned. For a system Ax =b , the matrix A condition number reads how sensitive the system solution x is with respect to a small change of the input data, i.e., how much the solution x deviates when a small numerical change is made in the system matrix Aor in the right-hand side vector b. For a general square regular matrix A , the condition number κ(A) may be determined as κ(A) = kAk · kA−1k. Consequently, if δA , δb denotes perturbations in matrix A and vector b , respectively, and if for any associated matrix norm it holds that kδAk<1 kA−1k, then [ 29 ], the solution ˜ x=x∗+δx to the system (A+δA)(x+δx) = b+δb approximates the solution x∗to the system Ax =bwith an error kδx∗k kx∗k≤κ(A) 1−κ(A)kδAk kAkkδbk kbk+kδAk kAk. Therefore, if the condition number κ(A) is greater than 1, then the relative error of the result may be κ(A) times greater than the sum of the relative errors of A and b . However, the calculation of the inverse matrix is too demanding and, therefore, we use the following assertion: κ(A) = σmax(A) σmin(A), where σmax(A) and σmin(A) are the greatest and the least singular values of A , respectively. To determine the singular values of A, we use the following theorem [30]. Theorem 1. Let ai· denote the i -th row of matrix An,n and a·j denote the j -th column of matrix An,n. The lower bound for the largest singular value is σmax(A) = maxmax 1≤i≤n{kai·k2}, max 1≤j≤nka·jk2 and the upper bound for the smallest singular value is σmin(A) = minmin 1≤i≤n{kai·k2}, min 1≤j≤nka·jk2. Because in the submatrix formed by the first four rows and columns of the matrix A there are numbers of order 10 7 in general, while in the last row and column of the matrix A