An SQP method for Chebyshev and hole-pattern fitting with geometrical elements

A customized sequential quadratic program (SQP) method for the solution of minimax-type fitting applications in coordinate metrology is presented. This area increasingly requires highly efficient and accurate algorithms, as modern three-dimensional geometry measurement systems provide large and computationally intensive data sets for fitting calculations. In order to meet these aspects, approaches for an optimization and parallelization of the SQP method are provided. The implementation is verified with medium (500 thousand points) and large (up to 13 million points) test data sets. A relative accuracy of the results in the range of 1× 10−14 is observed. With four-CPU parallelization, the associated calculation time has been less than 5 s.


Introduction
Three-dimensional coordinate metrology is an essential element of modern economic production.The application of coordinate measuring machines (CMMs) offers an efficient way to inspect geometrical properties of workpieces.The manufactured workpieces are measured in a first step.Then, through the mathematical fitting of ideal geometrical elements to the measured features, an evaluation of workpiece deviations from the nominal shape in a technical drawing is possible.The information generated on the deviations is used to classify workpieces as, for example, scrap parts and to adjust the manufacturing parameters in order to reduce the percentage of non-permissible parts in the production process.
Optical and computer tomography (CT) measurement systems have gained increasing importance for threedimensional coordinate metrology in recent years.The rapid development of new applications and measurement capabilities poses a great challenge for developers of CMM software as well as for metrological institutes to provide traceability for different measurands.Fitting software has to keep up with an increasing amount of measurement data and simultaneously the accuracy of calculation results needs to meet high requirements.Under these conditions, PTB (Physikalisch-Technische Bundesanstalt) together with the German CCM manufacturer Werth Messtechnik GmbH have developed a sequential quadratic program (SQP) for the calculation of Chebyshev and hole-pattern fitting applications with different geometrical elements.Through very simple modifications and computational parallelization, this method is capable of highly accurate and efficient calculations with large data sets as necessary for optical and CT measurements of workpieces in manufacturing.
The fitting algorithm in this work is subject to applications with different types of minimax calculations.These are Chebyshev, minimum-circumscribed and maximuminscribed geometrical elements (commonly all denoted by Chebyshev fitting) as well as hole-pattern fitting of structures with multiple geometrical elements.Modern algorithms rely on a combination of methods for calculating an initial course fit with a subsequent refinement by a decent direction-based solver.Linear methods are used to calculate decent directions as well as approximations of the fitting model gradients.The following algorithm presents a decent direction method.Its implementation is based on the general SQP and the activeset quadratic program (QP) for decent direction calculation from Geiger and Kanzow (2002).In comparison to the aforementioned linear approaches it calculates with exact model gradients.This requires more calculation time than the ap-D.Hutzschenreuter et al.: An SQP method for Chebyshev and hole-pattern fitting proximation.However, a significant higher accuracy of fitting results can be achieved.Furthermore, the QP decent directions provide quadratic or at least super linear convergence speed of the method, which is an advantage compared to a linear approach.
The input for the Chebyshev and hole-pattern fitting is a point data set: where P i = P ix , P iy , P iz T ∈ R 3 (i = 1, . ..m) are Cartesian coordinates of points on the surface of a real workpiece.The ideal geometrical element for the fit is defined by a number of parameters: These give either the geometrical element parameters in the case of Chebyshev applications, or transformation parameters in the case of hole-pattern fitting.In general, the parameter values are overestimated as the amount of measurement points m is much larger than the number n − 1 of unknown parameters.Hence, the fit considers the orthogonal distances, between the points P i ∈ P and the ideal geometry with parameters a.The model for the calculation of the geometrical parameters is then the minimax program min a∈R n−1 max i=1,...,m f i (a). (1) The objective is to find parameter values for a that minimize the maximum orthogonal distance between all points and the geometrical element.For the development of SQP methods, program (1) is equivalently transformed into an ordinary non-linear constrained form.
Here, v := a T , s T ∈ R n denotes the extended parameter vector.The new parameter s is an upper bound for the maximum over all f i (a) in program (1).By introducing f (v) := s as an objective and the set of inequality constraints g i (v) := f i (a) − s, it equivalently replaces the minimization of the maximum term.The equality constraints h j (v) complete the model and implement additional conditions related to the definition of the parameter vector a.
Section 2 gives the SQP core algorithm.It subsequently requires the application of a solver for quadratic programs.An optimized quadratic program solver for the fitting applications of interest is given in Sect.3.For specific nonregular solutions of program (2), a second-order criterion for a solution is presented in Sect. 4. The numerical precision of calculations by the customized SQP method is then investigated in Sect. 5 using the example of flange ring inspections (Hutzschenreuter et al., 2017) that require different types of minimax fits to be calculated.Finally, the effect of the partial parallelization of the algorithm and the resulting runtime improvements are investigated in Sect.6.
The initial values for the parameter vector shall be a feasible point of the minimization program.Such a point satisfies the constraints g i v 0 ≤ 0 (i = 1, . .., m) and h j v 0 = 0 (j = 1, . .., q).The initial Lagrange multiplier values are selected as zero values.In each iteration of the algorithm, the update of the Lagrange multipliers is calculated by solving the QP in step (S2).Thereby, λ k+1 and µ k+1 are the Lagrange multipliers associated with program (4).Details on their definition are presented in the following section.
In step (S3) formula ( 5), the Armijo-type step width is calculated from the exact l 1 -barrier-penalty function.
Step (S3) is a commonly used globalization technique for enabling the super-linear convergence up to the quadratic convergence of the SQP method towards a solution of the KKT system (3).Parameter α controls the influence of invalid constraints on the step width; β is the basis for the Armijo steps and σ adds some additional damping.

Customized active-set QP method
Solving the QPs (4) in step (S2) of the SQP method is the most time-consuming part when considering fitting large data sets.The applied algorithm is an active-set method, which overcomes the efficiency problem by introducing simple modifications.A minimum of program (4) defines the active set of inequality constraints that are equal to zero.
Finally, h denotes a vector whose components are all h j (v k ) and H is the matrix with the columns ∇ v h T j (v k ) (j = 0, . .., q).The KKT condition for a local solution of program (4) introduces the Lagrange multipliers λ A ∈ R |A| for all the active constraints and µ ∈ R q for the equality constraints.
With these preliminary definitions for linearly independent The following algorithm uses approximations A l of A and the KKT condition ( 6) to calculate a solution of program (4).It uses iteratively calculated approximations of v with the index l = 0, 1, 2, . . .that are Algorithm 2 (customized active-set method for QP) (S0) For l = 0, calculate initial values v 0 , A 0 , λ A 0 and µ by solving (S3.2) Otherwise create A l+1 from A l by removing the index with the smallest λ i < 0 value and set (S4.2) Otherwise calculate (S5) Set l = l + 1 and go to step (S1).
The initialization step (S0) is a specific customization of the method that speeds up the calculations significantly.Both vector a 0 and s 0 are derived from the solution v * of the linear equation system ( 7).These values satisfy the constraints of QP (4).It is only applicable for the minimax fitting models (1) and respectively for its transformation (2).
For large data sets, the degeneration of the linear equation system (8) may occur.Then some of the active constraint columns become linearly dependent and a unique solution of the update direction d is not available.Through a heuristic approach, the system matrix in Eq. ( 8) is extended by the regularization term −ρ g A l , which is a diagonal matrix with the components of g A l weighted with 0 < ρ 1.The equation system that is solved instead of Eq. ( 8) is as follows.8) is degenerated, the calculation is stopped and a non-optimal solution for v = v l is returned to the SQP method.Furthermore, a maximum for the iteration index l should be set to guarantee the finite termination of the method.If the maximum index is reached, a non-optimal solution is also returned to the SQP.
Note that a solution v l , λ A l , µ from Algorithm 2 provides the Lagrange multipliers for the next SQP iteration µ k+1 = µ, λ k+1 i = 0 (i ∈ A l ), and λ k+1 i is the corresponding value from λ A l for all i ∈ A l .

Second-order sufficiency condition
The aim of applying a second-order sufficiency condition is to verify that a solution of the SQP algorithms is a valid minimum of the fitting program (2).For this purpose, the following theorem form (1) is utilized.Theorem (second-order sufficiency condition) Let (v, λ, µ) be a solution of (2) that satisfies KKT condition (3).If then the solution is a strict minimum of the fitting application.
Computationally, condition (10) can be treated by a quadratic program and by minimizing z T ∇ vv L (v, λ, µ) z.The application of Algorithm 2 is valid.In the implementation, the Hessian is approximated by central difference quotients.The active set A and the Lagrange multipliers λ, µ used in the theorem above are those from the last QP iteration in Algorithm 1.
Condition (10) implies that the Hessian matrix ∇ vv L (v, λ, µ) of the Lagrangian is positive definite on the orthogonal complement of the affine subspace in R n , which is spanned by the gradients of active constraints with positive Lagrange multiplier values and the gradients of the equality constraints.Its application is necessary for special data sets where active constraints have λ i = 0 or where the number of linearly independent active constraint gradients and equality constraint gradients is less than n.

Test calculations
The performance of the software is evaluated for hole-pattern fitting with multiple geometrical elements which are located inside flange rings.This is a typical application for the inspection of workpieces by means of the geometrical product specification (GPS) standard ISO 2692 (2007) using a CMM.
The position and size of the five clearance holes of the flange in Fig. 1 are inspected by fitting a virtual counterpart consisting of five pins with an ideal geometrical shape into the holes.The alignment of the virtual counterpart to the holes is constrained by the datum elements A (Chebyshev plane) and B (minimum-circumscribed cylinder).A CMM measures points on the inner surface of all clearance holes and on the surfaces A and B. The data are denoted as the extracted workpiece and outlined in Fig. 2.
The extracted geometry of plane A is outlined by large round dots.These give the measurement points P P 1 , . .., P P m p .Points P C 1 , . .., P C m C are the extracted cylinder surface B which is marked by rectangular dots in Fig. 2. Finally, the small round dots outline the geometry measured in the holes.These have the point sets P k 1 , . .., P k m k , where k is a unique index for each hole.
In the first calculation step, datum plane A is fitted to the extracted flange.The fitting model is min max i=1,...,m P f i C P , − → n (11) The amount of points for the plane is m P .The geometrical parameters are the point C P on the plane and its normal vector − → n ; f i (C P , − → n ) is the orthogonal distance between the point P P i and the plane associated with C P , − → n .Furthermore, the geometrical parameters have two constraints.The length of the normal vector must be one.The point on the plane is defined as the orthogonal projection of the extracted plane centroid G = 1/m P m P i=1 P P i on the associated plane.In the SQP method, it is implemented by constructing two nonparallel vectors − → v 1 , − → v 2 that are orthogonal to − → n and by considering the equality constraints (G − In the second calculation step, the minimumcircumscribed cylinder is fitted to the extracted outer surface B. The axis direction of the associated cylinder is constrained to be parallel to the plane normal vector that was calculated before.Then the fitting model is The geometrical parameters C C for the fit are the coordinates of a point on the cylinder axis.Here, f i (C C ) gives the or-  thogonal distance between the cylinder axis and the point P C i on the extracted surface B. Furthermore, the position of the cylinder axis is constrained to be an element of the associated plane A.
Figure 3 shows an outline of the associated plane and cylinder.The point C C and normal vector − → n are used to specify a Cartesian workpiece coordinate system (x W , y W , z W ). The direction of the z axis is defined by − → n .All the Cartesian axes intersect in C C .
In the last calculation step, the hole-pattern fit of the clearance holes is calculated.Figure 4 shows the fitting of the gauge and the extracted geometry in parallel projection onto the x-y plane of the workpiece coordinate system.
The virtual counterpart for the fit represents five cylinders that are the solid circular elements in Fig. 4. Their axes are parallel to the z axis of the Cartesian coordinate system (x G , y G , z G ), which is denoted as a gauge coordinate system.This coordinate system is aligned with the workpiece coor-dinate system in the centre point C C and with z G parallel to z W .At fitting, the whole gauge coordinate system and the counterpart elements are rotated around the z axis by the angle ϕ.If there is one position where all cylinder elements fit into the extracted holes and no point overlaps with the inside of the cylinders, the flange ring is within its tolerance for the maximum permissible shape deviations.The minimax program for the calculation of such an angle is min max k=1,...5 i=1,...,m k f ki (ϕ) . ( Here k is an index that uniquely identifies each clearance hole and its associated element of the counterpart.Number m k gives the amount of measurement points for each hole.The orthogonal distance f ki (ϕ) between the outer surface of the cylinder with the index k and the point P k i is positive if the point is inside the cylinder (overlapping).Otherwise the distance has a negative value.All cylinder elements and hence the virtual counterpart fit into the clearance holes at the same time if at a solution of program ( 13).
In order to investigate the precision of the customized SQP method for the calculations ( 11), ( 12) and ( 13) of the flange ring inspection, test data sets DS1, . . ., DS6 were created by an inverse data generator.It relies on the stateof-the-art methods in Anthony et al. (1993), Forbes and Minh (2012) and Hutzschenreuter et al. (2015).Each data set refers to a nominal flange with ten clearance holes of radius 5 mm, an outer ring radius of 120 mm and a height of 20 mm.A reference solution of the flange fit is input for the data generator.It is given by the gauge coordinate system C R , x R G y R G z R G and the maximum distance s R for the solution of program ( 12).The generator then computes positions of measurement points that give a test data set for which the reference parameters represent the unambiguous global minimum of fitting.The points differ from the nominal shape of the flange by a random uniform form deviation.All points are equally distributed along rectangular surficial grids on the nominal flange surface.It simulates the dense probing from optical or likewise CT-sensor-based measurements.The gauge coordinate system C R , x R G y R G z R G and s C from the calculation with the SQP method is compared to the reference coordinate system according to Fig. 5.
The deviations between both coordinate systems are calculated by  γ Calculations for the six data sets were made with an absolute convergence tolerance for the SQP steps of 1 × 10 −14 .The algorithm is implemented in C++ (full compiler optimization, Microsoft Visual Studio, 2013), double precision floating point number format (IEEE 754, 2008).Figure 6 presents the deviations of the computed fit to the reference solution.For all data sets, the SQP method converged to a solution of the flange ring fitting.The maximum distance deviation ds and the gauge coordinate system position deviation dC are between 1 × 10 −13 and 1 × 10 −15 mm, which corresponds to the convergence tolerance.Similarly, the deviation of all coordinate system axes is of the size 1 × 10 −16 rad.This deviation is in the size of the relative machine precision of real number representation in floating point arithmetic (IEEE 754, 2008).In conclusion, the SQP method was able to calculate highly numerically precise solutions for the flange ring fitting examples.
Details on each test data set and SQP method parameters are summarized in the Appendices A and B. For all calculations the tolerance ε = 1 × 10 −14 was used to stop iterations when the decent direction converges with v ≤ ε.The convergence of the QP method has to be tested with a tolerance threshold smaller or equal to ε, or else the SQP method may not converge to a solution.Additional bounds for the maximum of SQP iterations and QP iterations are set to ensure the finite termination of the algorithm.For the flange these are 100 SQP iterations.The number of QP iterations is, in comparison, set to a maximum of 10 for the Chebyshev plane and hole-pattern fit as well as 20 for the MC cylinder.If the QP method does not converge to a solution within the given number of iterations, then the SQP method continues with a non-optimal decent direction which has no significant influence on the overall convergence of the algorithm.The step-width control parameters are set to α = 1 for all fitting applications.It allows large decent steps which violate boundary constraints.All violations are corrected by the initial value selection of each subsequent QP method call.This approach bypasses slow convergence properties like the Maratos effect (Geiger and Kanzow, 2002).In addition, the step scaling parameter was set to β ≥ 0.7 and the maximum number of Armijo steps is 2, which forces large steps in each SQP iteration.
Further verification of the customized SQP method was made by using the public test for Chebyshev geometrical elements of PTB (Hutzschenreuter et al., 2015).Calculations were made for 50 different geometrical element data sets.The elements covered by the test are a two-dimensional circle, two-dimensional straight line, plane sphere and cylinder.In comparison to the previously used flange ring data sets the Chebyshev test data sets also simulate the effect of systematic form deviations such as harmonic deviations and convexity as well as full and partitial features.In these cases it is more likely that an insufficient implementation of an fitting algorithm will end up in an local minimum which is not the required fitting solution.For all elements in the test, the initial values for the Chebyshev fitting were calculated by Gaussian fitting.Conveniently, the implementation of the Gaussian fitting is possible with the SQP method.Only the convergence speed towards a solution is very slow for this non-minimax-type fitting.With these preliminaries the public Chebyshev element test was passed for all 50 test data sets with an maximum permissible error of 0.1 µm for position and size parameters, 0.01 µm for form deviation and 0.1 µrad for the orientation of the geometrical element parameters.
Finally, the following additional test calculations point out the effect of outliers in the given measurement data sets and the influence of different selections of initial values for the SQP method.
Outliers denote perturbations in the measurement data of the extracted geometry that can occur when dirt particles such as dust stick to the workpiece surface at the measurement.A reliable evaluation of the geometrical measurands requires removal of these outlier points form the data before applying Chebyshev and hole-pattern fitting as they have a considerable influence on the solution of fitting program (1).The situation is outlined for the flange data set DS1. Measurement point P P 1 of the datum plane is shifted step-wise by setting P P 1 = P P 1 +a• − → n with a = 0, 0.1, 0.2, . .., 1. Thereby, the normal n is the direction given by the reference plane fit of the test data set.The resulting deviation between reference flange fit and the fit with the SQP method is shown in Fig. 7.
The deviation of the fitting solution is strongly correlated with the shift of P P 1 .For a ≥ 0.3 mm the deviation of the fitting parameter s of the hole-pattern fit becomes larger than 0.01 mm.In this case the solution of the SQP method changes from the virtual gauge fitting into all holes to an overlapping.The workpiece would be identified as out of its tolerance only because of the presence of one outlier.
The choice of an initial solution v 0 affects the calculation time of the SQP method and whether the algorithm converges to the global fitting minimum or to some inadequate local solution.For the Chebyshev plane fit in program (11) it is sufficient to select the gravity centre point of the given measurement data as initial position.An initial normal can be drawn from a list of candidate direction vectors as the normal direction that provides the smallest form deviation with the initial centre point.Similarly the initial approximation for the MC cylinder in program ( 12) can be the fitting result of the Chebyshev plane.For the hole-pattern fit a more sophisticated method is required to determine an initial rotation ϕ 0 of the virtual gauge.An example is the Gaussian centre-point fit presented in Hutzschenreuter et al. (2017).It provides fast convergence for the flange fit examples.The effect of choosing different initial values for the hole-pattern fit is shown in Fig. 8.
The initial ankle ϕ 0 = −0.034rad is varied by setting ϕ ).In the upper half of the figure the asynchronous behaviour of the calculation time can be seen.Here, t is the total calculation time, tA is the time for assembling distance values and their gradients, tQP is the time for solving the customized active set method, tL1 is the time for the line search and tKKT gives the time for checking the Karush-Kuhn-Tucker condition for an optimal solution.The area for the best convergence of the SQP method towards the global minimum is observed for −0.2 ≤ d ≤ 0.6.For d = −0.6 and d = 0.8 the amount of QP iterations increases significantly and results in a calculation time about twice as long as for the previous cases.For d = −1, −0.8, −0.4 and d = 1 the SQP method converges to different local solutions that give a wrong fitting result.These calculations have up to 50 SQP iterations and more than 350 QP iterations.Moreover, all calculations show that the QP method has the greatest share of the total calculation time.In the best case with d = 0, the SQP method requires only four iterations to converge to the global minimum with relative precision, as shown in Fig. 6.The amount of QP iterations is 10.In the worst case -with a correct solution -28 SQP iterations and 214 QP iterations have been counted at d = 0.8.Especially for data sets with uneven size ratios such as long cylinder elements or partitial features, the area for fast convergence of the SQP method becomes more narrow than in the test examples DS1-DS6.

Parallelization
The efficiency of the calculation is improved by the partial parallelization of the Algorithms 1 and 2. What is suitable for parallelization is the computation of orthogonal distance values and their gradients as well as expensive vector operations that frequently have to be made within the QP active-set method.Due to the iterative structure of the SQP method and the QP active-set method, the basic part of the algorithms remains serial.Hence, it is expected that the speedup will not be proportional to the amount of available threads for parallel calculation (e.g.available CPU cores).
The times in Fig. 9 were measured for calculations of the flange ring fit to data sets DS1, . . ., DS6 of the previous section at a workstation with an Intel(R) W3550 CPU and absolute convergence tolerance 1 × 10 −14 .For all data sets, the speedup for four threads is approximately between factor 3 and factor 4 compared to the serial computation with 1 thread.The fit for the small and medium-sized data sets DS1, . . ., DS4 (up to 6 million points) is less than 1 s for four thread calculations.The speedup decreases with an increasing amount of threads.

Conclusions
Minimax-type fitting applications in coordinate metrology such as Chebyshev, minimum-circumscribed, maximuminscribed and hole-pattern fitting can be transformed into a convenient general non-linear constrained form.For this formulation of the applications, a very efficient solver was developed from a basic SQP method by means of simple modifications.For flange ring fitting test data sets, the software was able to compute highly numerically precise solutions with a relative deviation of 1 × 10 −14 of the numerical values from the reference solutions.In addition, a verification with TraCIM Chebyshev element test data sets was made.
The parallelization of the algorithms showed a significant speedup of the calculation time even for a small amount of parallel threads used.It allows the fitting with a huge amount of data (13 million points in the test data sets) in less than 5 s on a four-core CPU that is increasingly becoming the standard even for office computers.Subsequently, the performance capacity of modern multicore workstation CPUs (8, 16 and more cores) can almost fully be exploited for further reduction of the calculation time.
In further metrological applications, the SQP method is suitable for the simulation of numerical uncertainties of geometrical element and fitting calculation results by Monte Carlo methods (JCGM, 2008).Moreover, due to its accuracy, it can be used for the validation of data sets from industrial test services -for example the TraCIM validation service (PTB, 2015).
Finally, one aim is to provide the software to users of CMM as part of the WinWerth evaluation software offered by the multi-sensor CCM manufacturer Werth (Werth Messtechnik GmbH, 2017).An implementation is provided for different three-dimensional hole-pattern fitting applications such as the calculation for the flange ring.Further details are given in the PTB hole-pattern fitting guideline (Hutzschenreuter et al., 2017).
Data availability.The authors' original test data are available at PTBs public repository (Hutzschenreuter, 2017).It comprises several ASCII formatted files for each test case in Sects.5 and 6.Supplemental plot data from Figs. 6 and 7 are also located within this repository.The public Chebyshev test data sets are available at PTB's TraCIM service website (PTB, 2015).

Figure 3 .
Figure 3. Associated Chebyshev plane (element with dashed edges) and minimum-circumscribed cylinder (element with dash-dotted edges) for the flange.

Figure 4 .
Figure 4. Hole-pattern fit of the flange.In the actual alignment, the gauge overlaps with the extracted workpiece which is outlined by the black dots intersecting with the circular gauge elements.

Figure 5 .
Figure 5. Principle of the comparison between reference and calculated gauge coordinate systems.

Figure 6 .
Figure 6.Deviation of the calculated flange fit to the reference solution.

Figure 7 .
Figure 7. Deviation between the calculated flange fit and the reference solution for a simulated outlier point in the extracted plane data set.

Figure 8 .
Figure 8.Effect of a variation of the initial value for the flange hole-pattern fit on the calculated SQP solution and its performance with DS1.

Figure 9 .
Figure 9. Calculation times for flange ring fitting with different amounts of parallel threads.