the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An SQP method for Chebyshev and holepattern fitting with geometrical elements
Frank Härtig
Markus Schmidt
A customized sequential quadratic program (SQP) method for the solution of minimaxtype fitting applications in coordinate metrology is presented. This area increasingly requires highly efficient and accurate algorithms, as modern threedimensional 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 fourCPU parallelization, the associated calculation time has been less than 5 s.
 Article
(3082 KB)  Fulltext XML
 BibTeX
 EndNote
Threedimensional 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 nonpermissible 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 (PhysikalischTechnische Bundesanstalt) together with the German CCM manufacturer Werth Messtechnik GmbH have developed a sequential quadratic program (SQP) for the calculation of Chebyshev and holepattern 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, minimumcircumscribed and maximuminscribed geometrical elements (commonly all denoted by Chebyshev fitting) as well as holepattern 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 directionbased 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 approximation. 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 holepattern fitting is a point data set:
where ${P}_{i}={\left({P}_{ix},\phantom{\rule{0.125em}{0ex}}{P}_{iy},{P}_{iz}\right)}^{\mathrm{T}}\in {\mathbb{R}}^{\mathrm{3}}\left(i=\mathrm{1},\mathrm{\dots}m\right)$ 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 holepattern 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
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 nonlinear constrained form.
Here, $\mathit{v}:={\left({\mathit{a}}^{\mathrm{T}},s\right)}^{\mathrm{T}}\in {\mathbb{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}\left(\mathit{v}\right):={f}_{i}\left(\mathit{a}\right)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 secondorder 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.
For the application of an SQP method, the Lagrangian of program (2) is considered. For the Lagrange multipliers $\mathit{\lambda}={\left({\mathit{\lambda}}_{\mathrm{1}},\mathrm{\dots},{\mathit{\lambda}}_{m}\right)}^{\mathrm{T}}\in {\mathbb{R}}^{m}$ and $\mathit{\mu}={\left({\mathit{\mu}}_{\mathrm{1}},\mathrm{\dots},{\mathit{\mu}}_{q}\right)}^{\mathrm{T}}\in {\mathbb{R}}^{q}$, it is
Under the assumption of a regularity constraint qualification (e.g. the linear inequality constraint qualification; Geiger and Kanzow, 2002), a local minimum of (2) and its Lagrangian have multipliers for which the Karush–Kuhn–Tucker (KKT) condition
is satisfied. The SQP method is derived from this necessary condition by applying a Lagrange–Newton method to ${\mathrm{\nabla}}_{\mathit{v}}L\left(\mathit{v},\mathit{\lambda},\mathit{\mu}\right)=\mathit{v}$, ${h}_{j}\left(\mathit{v}\right)=\mathbf{0}\phantom{\rule{0.125em}{0ex}}\left(j=\mathrm{1},\mathrm{\dots},q\right)$ and ${g}_{i}\left(\mathit{v}\right)\le \mathrm{0}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\left(i=\mathrm{1},\mathrm{\dots},m\right)$.
Algorithm 1 (globalized SQP method)
(S0) Set initial values for v^{0}, λ^{0}, μ^{0}, α>0, $\mathit{\beta}\in \left(\mathrm{0},\mathrm{1}\right)$, $\mathit{\sigma}\in \left(\mathrm{0},\mathrm{1}\right)$ and ${\mathbf{H}}_{\mathbf{0}}={\mathbf{I}}_{n}\in {\mathbb{R}}^{n\times n}$, k=0.
(S1) If $\left({\mathit{v}}^{k},\phantom{\rule{0.125em}{0ex}}{\mathit{\lambda}}^{k},\phantom{\rule{0.125em}{0ex}}{\mathit{\mu}}^{k}\phantom{\rule{0.125em}{0ex}}\right)$ satisfy KKT condition (3), stop.
(S2) Calculate a decent direction Δv∈ℝ^{n} and $\left({\mathit{\lambda}}^{k+\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{\mathit{\mu}}^{k+\mathrm{1}}\right)$ from the quadratic program
If ∥Δv∥=0, stop with solution $\left({\mathit{v}}^{k},\phantom{\rule{0.125em}{0ex}}{\mathit{\lambda}}^{k+\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{\mathit{\mu}}^{k+\mathrm{1}}\right).$
(S3) Calculate a step width ${t}_{k}\ge \mathrm{0}\phantom{\rule{0.25em}{0ex}}({t}_{k}={\mathit{\beta}}^{l},\phantom{\rule{0.125em}{0ex}}l=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots})$ that satisfies
(S4) Update
${\mathit{v}}^{k+\mathrm{1}}={\mathit{v}}^{k}+{t}_{k}\mathbf{\Delta}\mathit{v}$, BFGS update of
${\mathbf{H}}_{k+\mathrm{1}},k=k+\mathrm{1}$ and go to step (S1).
In Algorithm 1, the matrix H_{k} is a positive semidefinite approximation of the Lagrange function's Hessian. It guarantees that the quadratic program (4) in step (S2) has a solution. Its update is made with the Broyden–Fletcher–Goldfarb–Shanno method (BFGS, Geiger and Kanzow, 2002).
The initial values for the parameter vector shall be a feasible point of the minimization program. Such a point satisfies the constraints ${g}_{i}\left({\mathit{v}}^{\mathrm{0}}\right)\le \mathrm{0}\phantom{\rule{0.125em}{0ex}}(i=\mathrm{1},\mathrm{\dots},m)$ and ${h}_{j}\left({\mathit{v}}^{\mathrm{0}}\right)=\mathrm{0}\phantom{\rule{0.125em}{0ex}}\left(j=\mathrm{1},\mathrm{\dots},q\right)$. 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 Armijotype step width is calculated from the exact l_{1}barrierpenalty function.
The term ${P}_{\mathrm{1}}^{\prime}({\mathit{v}}^{k},\mathbf{\Delta}\mathit{v},\mathit{\alpha})$ is its directional derivative in v^{k} towards Δv. Step (S3) is a commonly used globalization technique for enabling the superlinear 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.
Solving the QPs (4) in step (S2) of the SQP method is the most timeconsuming part when considering fitting large data sets. The applied algorithm is an activeset 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.
Through g_{A}∈ℝ^{A}, a vector is introduced whose components are g_{i}(v^{k}) with i∈A. The matrix ${\mathbf{G}}_{A}\in {\mathbb{R}}^{\leftA\right\times n}$ has the rows ${\mathrm{\nabla}}_{\mathit{v}}{g}_{i}^{\mathrm{T}}\left({\mathit{v}}^{k}\right)$ for all i∈A. Finally, h denotes a vector whose components are all h_{j}(v^{k}) and H is the matrix with the columns ${\mathrm{\nabla}}_{\mathit{v}}{h}_{j}^{\mathrm{T}}\left({\mathit{v}}^{k}\right)\left(j=\mathrm{0},\mathrm{\dots},q\right)$. The KKT condition for a local solution of program (4) introduces the Lagrange multipliers λ_{A}∈ℝ^{A} for all the active constraints and μ∈ℝ^{q} for the equality constraints.
With these preliminary definitions for linearly independent ${\mathrm{\nabla}}_{\mathit{v}}{g}_{i}^{\mathrm{T}}\left({\mathit{v}}^{k}\right)\phantom{\rule{0.125em}{0ex}}\left(i\in A\right)$ and ${\mathrm{\nabla}}_{\mathit{v}}{h}_{j}^{\mathrm{T}}\left({\mathit{v}}^{k}\right)\left(j=\mathrm{0},\mathrm{\dots},q\right)$, the KKT condition for QP (4) is
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=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots}$ that are
Algorithm 2 (customized activeset method for QP)
(S0) For l=0, calculate initial values $\mathbf{\Delta}{\mathit{v}}^{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}{A}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}{\mathit{\lambda}}_{{A}_{\mathrm{0}}}$ and μ by solving
and
(S1) If $\left(\mathbf{\Delta}{\mathit{v}}^{l},{\mathit{\lambda}}_{{A}_{l}},\mathit{\mu}\right)$ satisfies KKT condition (6), stop.
(S2) Calculate an update direction d by solving the linear system
(S3) If ∥d∥=0

(S3.1) If λ_{i}≥0 for all i∈A_{l}, stop with the solution $\left(\mathbf{\Delta}{\mathit{v}}^{l},{\mathit{\lambda}}_{{A}_{l}},\mathit{\mu}\right)$ from (8).

(S3.2) Otherwise create A_{l+1} from A_{l} by removing the index with the smallest λ_{i}<0 value and set $\mathbf{\Delta}{\mathit{v}}^{l+\mathrm{1}}=\mathbf{\Delta}{\mathit{v}}^{l}.$
(S4) If ∥d∥≠0

(S4.1) If ${g}_{i}\left({\mathit{v}}^{k}\right)+{\mathrm{\nabla}}_{\mathit{v}}{g}_{i}{\left({\mathit{v}}^{k}\right)}^{\mathrm{T}}(\mathbf{\Delta}{\mathit{v}}^{l}+\mathit{d})\le \mathrm{0}\phantom{\rule{0.125em}{0ex}}\left(i=\mathrm{1},\mathrm{\dots},m\right)$, set $\mathbf{\Delta}{\mathit{v}}^{l+\mathrm{1}}=\mathbf{\Delta}{\mathit{v}}^{l}+\mathit{d}$ and ${A}_{l+\mathrm{1}}={A}_{l}.$

(S4.2) Otherwise calculate
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{K}_{i}:={\displaystyle \frac{{g}_{i}\left({\mathit{v}}^{k}\right){\mathrm{\nabla}}_{\mathit{v}}{g}_{i}({\mathit{v}}^{k}{)}^{T}\mathbf{\Delta}{\mathit{v}}^{l}}{{\mathrm{\nabla}}_{\mathit{v}{g}_{i}({\mathit{v}}^{k}{)}^{T}\phantom{\rule{0.125em}{0ex}}\mathit{d}}}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}(i=\mathrm{1},\mathrm{\dots},m),\\ {\displaystyle}& {\displaystyle}t={\mathrm{min}}_{i\notin {A}_{l}}\left\{{K}_{i}\left(\right)open="">{\mathrm{\nabla}}_{\mathit{v}}{g}_{i}({\mathit{v}}^{k}{)}^{T}\mathbf{\Delta}{\mathit{v}}^{l}\mathrm{0}\right\}\end{array}$$and set $\mathbf{\Delta}{\mathit{v}}^{l+\mathrm{1}}=\mathbf{\Delta}{\mathit{v}}^{l}+t\mathit{d},{A}_{l+\mathrm{1}}={A}_{l}\cup \left\{j\right\}.$
(S5) Set $l=l+\mathrm{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 $\mathit{\rho}\mathbf{\Lambda}\left({\mathit{g}}_{{A}_{l}}\right)$, which is a diagonal matrix with the components of ${\mathit{g}}_{{A}_{l}}$ weighted with $\mathrm{0}<\mathit{\rho}\ll \mathrm{1}$. The equation system that is solved instead of Eq. (8) is as follows.
If in subsequent iterations of Algorithm 2 (8) is degenerated, the calculation is stopped and a nonoptimal 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 nonoptimal solution is also returned to the SQP.
Note that a solution $\left(\mathbf{\Delta}{\mathit{v}}^{l},{\mathit{\lambda}}_{{A}_{l}},\mathit{\mu}\right)$ from Algorithm 2 provides the Lagrange multipliers for the next SQP iteration ${\mathit{\mu}}^{k+\mathrm{1}}=\mathit{\mu}$, ${\mathit{\lambda}}_{i}^{k+\mathrm{1}}=\mathrm{0}\phantom{\rule{0.125em}{0ex}}\left(i\notin {A}_{l}\right)$, and ${\mathit{\lambda}}_{i}^{k+\mathrm{1}}$ is the corresponding value from ${\mathit{\lambda}}_{{A}_{l}}$ for all i∈A_{l}.
The aim of applying a secondorder 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 (secondorder sufficiency condition)
Let $\left(\mathit{v},\mathit{\lambda},\mathit{\mu}\right)$ be a solution of (2) that satisfies KKT condition (3). If
holds for all $\mathit{z}\in {\mathbb{R}}^{n},\mathit{z}\ne \mathrm{0}$ with
and
then the solution is a strict minimum of the fitting application.
Computationally, condition (10) can be treated by a quadratic program and by minimizing ${\mathit{z}}^{\mathrm{T}}{\mathrm{\nabla}}_{\mathit{v}\mathit{v}}L\left(\mathit{v},\mathit{\lambda},\mathit{\mu}\right)\mathit{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 ${\mathrm{\nabla}}_{\mathit{v}\mathit{v}}L\left(\mathit{v},\mathit{\lambda},\mathit{\mu}\right)$ of the Lagrangian is positive definite on the orthogonal complement of the affine subspace in ℝ^{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.
The performance of the software is evaluated for holepattern 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 (minimumcircumscribed 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}_{\mathrm{1}}^{P},\mathrm{\dots},{P}_{{m}_{p}}^{P}$. Points ${P}_{\mathrm{1}}^{C},\mathrm{\dots},{P}_{{m}_{C}}^{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}_{\mathrm{1}}^{k},\mathrm{\dots},{P}_{{m}_{k}}^{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
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 $\overrightarrow{\mathit{n}}$; ${f}_{i}({C}_{P},\overrightarrow{\mathit{n}})$ is the orthogonal distance between the point ${P}_{i}^{P}$ and the plane associated with ${C}_{P},\overrightarrow{\mathit{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=\mathrm{1}/{m}_{P}\phantom{\rule{0.125em}{0ex}}\left({\sum}_{i=\mathrm{1}}^{{m}_{P}}{P}_{i}^{P}\right)\phantom{\rule{0.125em}{0ex}}$ on the associated plane. In the SQP method, it is implemented by constructing two nonparallel vectors ${\overrightarrow{\mathit{v}}}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{\overrightarrow{\mathit{v}}}_{\mathrm{2}}$ that are orthogonal to $\overrightarrow{\mathit{n}}$ and by considering the equality constraints ${\left(G{C}_{P}\right)}^{\mathrm{T}}{\overrightarrow{\mathit{v}}}_{\mathrm{1}}=\mathrm{0}$ and ${\left(G{C}_{P}\right)}^{\mathrm{T}}{\overrightarrow{\mathit{v}}}_{\mathrm{2}}=\mathrm{0}$.
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 orthogonal distance between the cylinder axis and the point ${P}_{i}^{C}$ 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 $\overrightarrow{\mathit{n}}$ are used to specify a Cartesian workpiece coordinate system $\left({x}_{W},{y}_{W},{z}_{W}\right)$. The direction of the z axis is defined by $\overrightarrow{\mathit{n}}$. All the Cartesian axes intersect in C_{C}.
In the last calculation step, the holepattern 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 $\left({x}_{G},{y}_{G},{z}_{G}\right)$, which is denoted as a gauge coordinate system. This coordinate system is aligned with the workpiece coordinate 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
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}_{i}^{k}$ 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 stateoftheart 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},\phantom{\rule{0.125em}{0ex}}{\mathit{x}}_{G}^{R}{\mathit{y}}_{G}^{R}{\mathit{z}}_{G}^{R}$ 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 CTsensorbased measurements. The gauge coordinate system ${C}_{R},\phantom{\rule{0.125em}{0ex}}{x}_{G}^{R}{y}_{G}^{R}{z}_{G}^{R}$ 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 holepattern 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 nonoptimal decent direction which has no significant influence on the overall convergence of the algorithm. The stepwidth 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 twodimensional circle, twodimensional 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 nonminimaxtype 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 holepattern 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}_{\mathrm{1}}^{P}$ of the datum plane is shifted stepwise by setting ${P}_{\mathrm{1}}^{P}={P}_{\mathrm{1}}^{P}+a\cdot \overrightarrow{\mathit{n}}$ with $a=\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{0.1},\phantom{\rule{0.125em}{0ex}}\mathrm{0.2},\mathrm{\dots},\mathrm{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}_{\mathrm{1}}^{P}$. For a≥0.3 mm the deviation of the fitting parameter s of the holepattern 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 holepattern fit a more sophisticated method is required to determine an initial rotation φ^{0} of the virtual gauge. An example is the Gaussian centrepoint 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 holepattern fit is shown in Fig. 8.
The initial ankle ${\mathit{\phi}}^{\mathrm{0}}=\mathrm{0.034}$ rad is varied by setting $\mathit{\phi}\left(d\right)={\mathit{\phi}}^{\mathrm{0}}+d\frac{\mathit{\pi}}{\mathrm{45}}$ with $d=\mathrm{1}+i\cdot \mathrm{0.2}\phantom{\rule{0.125em}{0ex}}(i=\mathrm{0},\mathrm{\dots},\mathrm{10})$. 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 $\mathrm{0.2}\le d\le \mathrm{0.6}$. For $d=\mathrm{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=\mathrm{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.
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 activeset method. Due to the iterative structure of the SQP method and the QP activeset 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 mediumsized 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.
Minimaxtype fitting applications in coordinate metrology such as Chebyshev, minimumcircumscribed, maximuminscribed and holepattern fitting can be transformed into a convenient general nonlinear 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 fourcore 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 multisensor CCM manufacturer Werth (Werth Messtechnik GmbH, 2017). An implementation is provided for different threedimensional holepattern fitting applications such as the calculation for the flange ring. Further details are given in the PTB holepattern fitting guideline (Hutzschenreuter et al., 2017).
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).
The rotation angles R_{x}, R_{y} and R_{z} in Table A2 refer to a rotation of the gauge (virtual counterpart of holepattern fit) around the coordinate system axes. Each rotation is clockwise. First the gauge is rotated around the z axis. A rotation around the y axis and finally around the x axis follow.
∥v∥  Euclidian norm of vector v 
I_{n}  Unit matrix of dimension n 
∇_{v}  Gradient operator for variables v 
(⋅)^{T}  Transposed vector or matrix 
This paper and its content, numbers and figures were authorized for publication by M. Schmidt, who is coauthor and leader of the software development department of Werth Messtechnik GmbH.
This research is funded by the German Federal Ministry of Economic Affairs and Energy and Werth Messtechnik GmbH within the scope of the MNPQTransfer project “3DLochbildeinpassung für Koordinatenmesssysteme mit taktilen und optischen Sensoren sowie CTSystemen”.
Thanks go to Matthias Franke from PTB for the support provided with
numerical calculations and Harald Löwe from the Technical
University of Braunschweig for his many helpful remarks which improved the
presentation of the algorithms. Finally, the authors like to thank the
anonymous reviewers for their many suggestions that helped to improve the
understanding of the SQP method and its application.
Edited by: Rainer Tutsch
Reviewed by: two anonymous referees
Anthony, G. T., Anthony, H. M., Bittner, B., Butler, B. P., Cox, M. G., Drieschner, R., Elligsen, R., Forbes, A. B., Groß, H., Hannaby, S. A., Harris, P. M., and Kok, J.: Chebychev bestfit geometric elements, Centrum voor Wiskunde en Informatica, Report NMR9317, 1993.
DIN EN ISO 2692: Geometrical product specifications (GPS) – Geometrical tolerancing – Maximum material requirement (MMR), least material requirement (LMR) and reciprocity requirement (RPR) (ISO 2692:2006), German version EN ISO 2692:200704, 2007.
Forbes, A. B. and Minh, H. D.: Generation of numerical artefacts for geometric form and tolerance assessment, Int. J. Metrol. Qual. Eng., 3, 145–150, https://doi.org/10.1051/ijmqe/2012022, 2012.
Geiger, C. and Kanzow, C.: Theorie und Numerik restringierter Optimierungsaufgaben, Springer, https://doi.org/10.1007/9783642560040, 2002.
Hutzschenreuter, D.: Benchmarking data sets for flange hole pattern fitting, PhysikalischTechnische Bundesanstalt (PTB), https://doi.org/10.7795/710.20170606, 2017.
Hutzschenreuter, D., Härtig, F., Wendt, K., Lunze, U., and Löwe, H.: Online validation of Chebyshev geometrical element algorithms using the TraCIMsystem, Journal of Mechanical Engineering and Automation, 5, 94–111, https://doi.org/10.5923/j.jmea.20150503.02, 2015.
Hutzschenreuter, D., Härtig, F., and Schmidt, M.: LoBEK 3D – User Guide to 3 D pattern fitting with coordinate measuring machines, 2017 – Version 1, https://doi.org/10.5923.j.jmea.20150503.02, 2017.
IEEE 7542008: Standard for FloatingPoint Arithmetic, IEEE Standards Association, https://doi.org/10.1109/IEEESTD.2008.4610935, 2008.
JCGM: 101:2008, Evaluation of measurement data – Supplement 1 to the “Guide to the expression of uncertainty in measurement” – Propagation of distributions using a Monte Carlo method, available at: http://www.bipm.org/utils/common/documents/jcgm/JCGM_101_2008_E.pdf, (last access: January 2018), 2008.
Mircrosoft Visual Studio Express 2013: available at: https://www.microsoft.com/enus/download/details.aspx?id=40769, (last access: January 2018), 2013.
PTB (PhysikalischTechnische Bundesanstalt): The TraCIM service web pages, available at: http://www.tracim.ptb.de (last access: June 2017), 2015.
Werth Messtechnik GmbH: Companies web pages, available at: http://www.werth.de (last access: June 2017), 2017.
 Abstract
 Introduction
 The SQP method
 Customized activeset QP method
 Secondorder sufficiency condition
 Test calculations
 Parallelization
 Conclusions
 Data availability
 Appendix A: Test data set properties
 Appendix B: SQP solver parameters
 Appendix C: Nomenclature
 Competing interests
 Acknowledgements
 References
 Abstract
 Introduction
 The SQP method
 Customized activeset QP method
 Secondorder sufficiency condition
 Test calculations
 Parallelization
 Conclusions
 Data availability
 Appendix A: Test data set properties
 Appendix B: SQP solver parameters
 Appendix C: Nomenclature
 Competing interests
 Acknowledgements
 References