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⁻¹⁴ is observed. With four-CPU parallelization, the associated calculation time has been less than 5 s.

1 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 three-dimensional 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 maximum-inscribed 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 active-set 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 hole-pattern fitting is a point data set:

$${\displaystyle }P=\left\{P_{\mathrm{1}},\mathrm{\dots },P_m\right\},$$

where $P_i={\left(P_{ix},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:

$${\displaystyle }\mathit{a}={\left(a_{\mathrm{1}},\mathrm{\dots },a_{n-\mathrm{1}}\right)}^{\mathrm{T}}\in {\mathbb{R}}^{n-\mathrm{1}}.$$

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,

$${\displaystyle }{f}_i\left(\mathit{a}\right),$$

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

$$\begin{array}{}\text{(1)}& {\text{min}}_{\mathit{a}\in {\mathbb{R}}^{n-\mathrm{1}}}{\text{max}}_{i=\mathrm{1},\mathrm{\dots },m}{f}_i\left(\mathit{a}\right).\end{array}$$

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.

$$\begin{array}{ll}\text{(2)}& {\displaystyle }& {\displaystyle }{\text{min}}_{\mathit{v}\in {\mathbb{R}}^n}f\left(\mathit{v}\right){\displaystyle }& {\displaystyle }\text{subject to}{g}_i\left(\mathit{v}\right)\le \mathrm{0}\left(i=\mathrm{1},\mathrm{\dots },m\right)\\ {\displaystyle }& {\displaystyle }{h}_j\left(\mathit{v}\right)=\mathrm{0}\left(j=\mathrm{1},\mathrm{\dots },q\right)\end{array}$$

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 non-regular 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.

2 The SQP method

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

$${\displaystyle }L\left(\mathit{v},\mathit{\lambda },\mathit{\mu }\right):=f\left(\mathit{v}\right)+{\sum }_{i=\mathrm{1}}^m{\mathit{\lambda }}_i{g}_i\left(\mathit{v}\right)+{\sum }_{j=\mathrm{1}}^q{\mathit{\mu }}_j{h}_j\left(\mathit{v}\right).$$

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

$$\begin{array}{ll}\text{(3)}& {\displaystyle }& {\displaystyle }{\mathrm{\nabla }}_{\mathit{v}}L\left(\mathit{v},\mathit{\lambda },\mathit{\mu }\right)=\mathbf{0}{\displaystyle }& {\displaystyle }{h}_j\left(\mathit{v}\right)=\mathrm{0}\left(j=\mathrm{1},\mathrm{\dots },q\right)\\ {\displaystyle }& {\displaystyle }{g}_i\left(\mathit{v}\right)\le \mathrm{0},{\mathit{\lambda }}_i\ge \mathrm{0},{\mathit{\lambda }}_i{g}_i\left(\mathit{v}\right)=\mathrm{0}\left(i=\mathrm{1},\mathrm{\dots },m\right)\end{array}$$

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}\left(j=\mathrm{1},\mathrm{\dots },q\right)$ and $g_i\left(\mathit{v}\right)\le \mathrm{0}\left(i=\mathrm{1},\mathrm{\dots },m\right)$.

Algorithm 1 (globalized SQP method)

(S0) Set initial values for v⁰, λ⁰, μ⁰, α>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,{\mathit{\lambda }}^k,{\mathit{\mu }}^k\right)$ satisfy KKT condition (3), stop.

(S2) Calculate a decent direction Δv∈ℝⁿ and $\left({\mathit{\lambda }}^{k+\mathrm{1}},{\mathit{\mu }}^{k+\mathrm{1}}\right)$ from the quadratic program

$$\begin{array}{ll}\text{(4)}& {\displaystyle }& {\displaystyle }{\text{min}}_{\mathbf{\Delta }\mathit{v}}{\mathrm{\nabla }}_{\mathit{v}}f{\left({\mathit{v}}^k\right)}^{\mathrm{T}}\mathbf{\Delta }\mathit{v}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{\Delta }{\mathit{v}}^{\mathrm{T}}{\mathrm{H}}_k\mathbf{\Delta }\mathit{v}\text{subject to}{\displaystyle }& {\displaystyle }{g}_i\left({\mathit{v}}^k\right)+{\mathrm{\nabla }}_{\mathit{v}}{g}_i{\left({\mathit{v}}^k\right)}^{\mathrm{T}}\mathbf{\Delta }\mathit{v}\le \mathrm{0}(i=\mathrm{1},\mathrm{\dots },m)\\ {\displaystyle }& {\displaystyle }{h}_j\left({\mathit{v}}^k\right)+{\mathrm{\nabla }}_{\mathit{v}}{h}_j{\left({\mathit{v}}^k\right)}^{\mathrm{T}}\mathbf{\Delta }\mathit{v}=\mathrm{0}\left(j=\mathrm{1},\mathrm{\dots },q\right).\end{array}$$

If ∥Δv∥=0, stop with solution $\left({\mathit{v}}^k,{\mathit{\lambda }}^{k+\mathrm{1}},{\mathit{\mu }}^{k+\mathrm{1}}\right).$

(S3) Calculate a step width $t_k\ge \mathrm{0}(t_k={\mathit{\beta }}^l,l=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots })$ that satisfies

$$\begin{array}{}\text{(5)}& P_{\mathrm{1}}\left({\mathit{v}}^k+t_k\mathbf{\Delta }\mathit{v},\mathit{\alpha }\right)\le P_{\mathrm{1}}\left({\mathit{v}}^k,\mathit{\alpha }\right)-t_k\mathit{\sigma }{P}_{\mathrm{1}}^{\prime }\left({\mathit{v}}^k,\mathbf{\Delta }\mathit{v},\mathit{\alpha }\right).\end{array}$$

(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 semi-definite 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}(i=\mathrm{1},\mathrm{\dots },m)$ and $h_j\left({\mathit{v}}^{\mathrm{0}}\right)=\mathrm{0}\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 Armijo-type step width is calculated from the exact l₁-barrier-penalty function.

$$\begin{array}{ll}{\displaystyle }{P}_{\mathrm{1}}\left(\mathit{v},\mathit{\alpha }\right)& {\displaystyle }:=f\left(\mathit{v}\right)+\mathit{\alpha }\left\{{\sum }_{i=\mathrm{1}}^m\text{max}\left\{\mathrm{0},g_i\left(\mathit{v}\right)\right\}\right.\\ {\displaystyle }& {\displaystyle }\left.+{\sum }_{j=\mathrm{1}}^q\left|h_j\left(\mathit{v}\right)\right|\right\}\end{array}$$

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

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

$${\displaystyle }A:=\left\{i\in \left\{\mathrm{1},\mathrm{\dots },m\right\}\mathrm{|}{g}_i\left(\mathit{v}\right)+{\mathrm{\nabla }}_v{g}_i{\left(\mathit{v}\right)}^{\mathrm{T}}\mathbf{\Delta }\mathit{v}=\mathrm{0}\right\}$$

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}}^{\left|A\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)\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

$$\begin{array}{ll}\text{(6)}& {\displaystyle }& {\displaystyle }\left(\begin{array}{ccc}{\mathbf{H}}_k& {\mathbf{G}}_A^{\text{T}}& {\mathbf{H}}^{\text{T}}\\ {\mathbf{G}}_A& \mathbf{0}& \mathbf{0}\\ \mathbf{H}& \mathbf{0}& \mathbf{0}\end{array}\right)\left(\begin{array}{c}\mathbf{\Delta }\mathit{v}\\ {\mathit{\lambda }}_A\\ \mathit{\mu }\end{array}\right)=\left(\begin{array}{c}-{\mathrm{\nabla }}_{\mathit{v}}f\left({\mathit{v}}^k\right)\\ -{\mathit{g}}_A\\ -\mathit{h}\end{array}\right){\displaystyle }& {\displaystyle }{g}_i\left({\mathit{v}}^k\right)+{\mathrm{\nabla }}_{\mathit{v}}{g}_i{\left({\mathit{v}}^k\right)}^{\text{T}}\mathrm{\Delta }v<\mathrm{0}\left(i\notin A\right)\\ {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_i\ge \mathrm{0}\left(i\in A\right),{\mathit{\lambda }}_i=\mathrm{0}\left(i\notin A\right)\end{array}$$

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

$${\displaystyle }\mathbf{\Delta }{\mathit{v}}^l:={\left({\left(\mathbf{\Delta }{\mathit{a}}^l\right)}^{\mathrm{T}},\mathrm{\Delta }{s}^l\right)}^{\mathrm{T}}={\left(\mathrm{\Delta }{a}_{\mathrm{1}}^l,\mathrm{\dots },\mathrm{\Delta }{a}_{n-\mathrm{1}}^l,\mathrm{\Delta }{s}^l\right)}^{\mathrm{T}}.$$

Algorithm 2 (customized active-set method for QP)

(S0) For l=0, calculate initial values $\mathbf{\Delta }{\mathit{v}}^{\mathrm{0}},A_{\mathrm{0}},{\mathit{\lambda }}_{A_{\mathrm{0}}}$ and μ by solving

$$\begin{array}{ll}\text{(7)}& {\displaystyle }& {\displaystyle }\left(\begin{array}{cc}{\mathrm{I}}_n& {\mathrm{H}}^{\text{T}}\\ \mathrm{H}& \mathrm{0}\end{array}\right)\left(\begin{array}{c}\mathbf{\Delta }{\mathit{v}}^{\ast }\\ {\mathit{\mu }}_{\mathrm{0}}\end{array}\right)=\left(\begin{array}{c}\mathbf{0}\\ -\mathit{h}\end{array}\right){\displaystyle }& {\displaystyle }\mathbf{\Delta }{\mathit{a}}^{\mathrm{0}}={\left(\mathrm{\Delta }{v}_{\mathrm{1}}^{\ast },\mathrm{\dots },\mathrm{\Delta }{v}_{n-\mathrm{1}}^{\ast }\right)}^{\mathrm{T}}\\ {\displaystyle }& {\displaystyle }\mathrm{\Delta }{s}^{\mathrm{0}}=\underset{i=\mathrm{1},\mathrm{\dots },m}{max}\left\{g_i\left({\mathit{v}}^k\right)+{\mathrm{\nabla }}_{\mathit{v}}{g}_i{\left({\mathit{v}}^k\right)}^{\mathrm{T}}\mathbf{\Delta }{\mathit{v}}^{\ast }\right\}\\ {\displaystyle }& {\displaystyle }\mathbf{\Delta }{\mathit{v}}^{\mathrm{0}}={\left(\mathbf{\Delta }{\mathit{a}}^{\mathrm{0}},\mathrm{\Delta }{s}^{\mathrm{0}}\right)}^{\mathrm{T}}\\ {\displaystyle }& {\displaystyle }{A}_{\mathrm{0}}=\left\{i\mathrm{|}{g}_i\left({\mathit{v}}^k\right)+{\mathrm{\nabla }}_{\mathit{v}}{g}_i{\left({\mathit{v}}^k\right)}^{\mathrm{T}}\mathbf{\Delta }{\mathit{v}}^{\mathrm{0}}=\mathrm{0}\right\}\\ {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_i\gg \mathrm{1}(i\in A_{\mathrm{0}}),\\ {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_i=\mathrm{0}\left(i\notin A_{\mathrm{0}}\right)\end{array}$$

and

$${\displaystyle }\mathit{\mu }={\mathit{\mu }}_{\mathrm{0}}.$$

(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

$$\begin{array}{ll}\text{(8)}& {\displaystyle }& {\displaystyle }\left(\begin{array}{ccc}{\mathbf{H}}_k& {\mathbf{G}}_{A_l}^{\text{T}}& {\mathbf{H}}^{\text{T}}\\ {\mathbf{G}}_{A_l}& \mathbf{0}& \mathbf{0}\\ \mathbf{H}& \mathbf{0}& \mathbf{0}\end{array}\right)\left(\begin{array}{c}\mathit{d}\\ {\mathit{\lambda }}_{A_l}\\ \mathit{\mu }\end{array}\right){\displaystyle }& {\displaystyle }=\left(\begin{array}{c}-{\mathrm{\nabla }}_{\mathit{v}}f\left({\mathit{v}}^k\right)-{\mathbf{H}}_k\mathbf{\Delta }{\mathit{v}}^l\\ \mathbf{0}\\ \mathbf{0}\end{array}\right)\end{array}$$

(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}\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\mathit{d}}}}(i=\mathrm{1},\mathrm{\dots },m),\\ {\displaystyle }& {\displaystyle }t={\min }_{i\notin A_l}\left\{K_i\left|{\mathrm{\nabla }}_{\mathit{v}}{g}_i({\mathit{v}}^k{)}^T\mathbf{\Delta }{\mathit{v}}^l>\mathrm{0}\right.\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⁰ and Δs⁰ 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.

$$\begin{array}{ll}\text{(9)}& {\displaystyle }& {\displaystyle }\left(\begin{array}{ccc}{\mathbf{H}}_k& {\mathbf{G}}_{A_l}^{\text{T}}& {\mathbf{H}}^{\text{T}}\\ {\mathbf{G}}_{A_l}& -\mathit{\rho }\mathbf{\Lambda }\left({\mathit{g}}_{A_l}\right)& \mathrm{0}\\ \mathbf{H}& \mathrm{0}& \mathrm{0}\end{array}\right)\left(\begin{array}{c}\mathit{d}\\ {\mathit{\lambda }}_{A_l}\\ \mathit{\mu }\end{array}\right){\displaystyle }& {\displaystyle }=\left(\begin{array}{c}-{\mathrm{\nabla }}_{\mathit{v}}f\left({\mathit{v}}_k\right)-{\mathbf{H}}_k\mathbf{\Delta }{\mathit{v}}^l\\ \mathbf{0}\\ \mathbf{0}\end{array}\right)\end{array}$$

If in subsequent iterations of Algorithm 2 (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 $\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}\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.

4 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 $\left(\mathit{v},\mathit{\lambda },\mathit{\mu }\right)$ be a solution of (2) that satisfies KKT condition (3). If

$$\begin{array}{}\text{(10)}& {\mathit{z}}^{\mathrm{T}}{\mathrm{\nabla }}_{\mathit{v}\mathit{v}}L\left(\mathit{v},\mathit{\lambda },\mathit{\mu }\right)\mathit{z}>\mathrm{0}\end{array}$$

holds for all $\mathit{z}\in {\mathbb{R}}^n,\mathit{z}\ne \mathrm{0}$ with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{\nabla }}_{\mathit{v}}{g}_i{\left(\mathit{v}\right)}^{\mathrm{T}}\mathit{z}=\mathrm{0}\text{for all}i\in A\text{with}{\mathit{\lambda }}_i>\mathrm{0}\\ {\displaystyle }& {\displaystyle }{\mathrm{\nabla }}_{\mathit{v}}{g}_i{\left(\mathit{v}\right)}^{\mathrm{T}}\mathit{z}\le \mathrm{0}\text{for all}i\in A\text{with}{\mathit{\lambda }}_i=\mathrm{0}\end{array}$$

and

$${\displaystyle }{\mathrm{\nabla }}_{\mathit{v}}{h}_j{\left(\mathit{v}\right)}^{\mathrm{T}}\mathit{z}=\mathrm{0}\text{for all}j=\mathrm{1},\mathrm{\dots },q$$

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 ℝⁿ, 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.

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

Figure 1 Nominal flange model and measured (extracted) features.

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.

Figure 2 Extracted flange.

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

$$\begin{array}{ll}\text{(11)}& {\displaystyle }& {\displaystyle }\text{min}{\text{max}}_{i=\mathrm{1},\mathrm{\dots },m_P}\left|f_i\left(C_P,\overrightarrow{\mathit{n}}\right)\right|{\displaystyle }& {\displaystyle }\text{subject to}∥\overrightarrow{\mathit{n}}∥=\mathrm{1},{\left(G-C_P\right)}^{\mathrm{T}}\overrightarrow{\mathit{v}}=\mathrm{0}\forall \overrightarrow{\mathit{v}}\perp \overrightarrow{\mathit{n}}.\end{array}$$

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\left({\sum }_{i=\mathrm{1}}^{m_P}{P}_i^P\right)$ on the associated plane. In the SQP method, it is implemented by constructing two non-parallel vectors ${\overrightarrow{\mathit{v}}}_{\mathrm{1}},{\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 minimum-circumscribed 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

$$\begin{array}{ll}\text{(12)}& {\displaystyle }& {\displaystyle }\text{min}{\text{max}}_{i=\mathrm{1},\mathrm{\dots },m_C}{f}_i\left(C_C\right){\displaystyle }& {\displaystyle }\text{subject to}{\left(C_P-C_C\right)}^{\mathrm{T}}\overrightarrow{\mathit{n}}=\mathrm{0}.\end{array}$$

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 Associated Chebyshev plane (element with dashed edges) and minimum-circumscribed cylinder (element with dash-dotted edges) for the flange.

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

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.

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

$$\begin{array}{}\text{(13)}& \text{min}{\text{max}}_{\begin{array}{c}k=\mathrm{1},\mathrm{\dots }\mathrm{5}\\ i=\mathrm{1},\mathrm{\dots },m_k\end{array}}{f}_{ki}\left(\mathit{\phi }\right).\end{array}$$

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

$$\begin{array}{}\text{(14)}& {\displaystyle }s:={\text{max}}_{\begin{array}{c}k=\mathrm{1},\mathrm{\dots }\mathrm{5}\\ i=\mathrm{1},\mathrm{\dots },m_k\end{array}}{f}_{ki}\left(\mathit{\phi }\right)\le \mathrm{0}\end{array}$$

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 state-of-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,{\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 CT-sensor-based measurements. The gauge coordinate system $C_R,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.

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

The deviations between both coordinate systems are calculated by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }ds:=\left|s^R-s^C\right|\\ {\displaystyle }& {\displaystyle }dC:=∥C_R-C_C∥\\ {\displaystyle }& {\displaystyle }\mathit{\alpha }:=∥x_G^R\times x_G^C∥\\ {\displaystyle }& {\displaystyle }\mathit{\beta }:=∥y_G^R\times y_G^C∥\\ {\displaystyle }& {\displaystyle }\mathit{\gamma }:=∥z_G^R\times z_G^C.∥\end{array}$$

Calculations for the six data sets were made with an absolute convergence tolerance for the SQP steps of 1 × 10⁻¹⁴. 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.

Figure 6 Deviation of the calculated flange 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⁻¹³ and 1 × 10⁻¹⁵ mm, which corresponds to the convergence tolerance. Similarly, the deviation of all coordinate system axes is of the size 1 × 10⁻¹⁶ 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⁻¹⁴ 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_{\mathrm{1}}^P$ of the datum plane is shifted step-wise by setting $P_{\mathrm{1}}^P=P_{\mathrm{1}}^P+a\cdot \overrightarrow{\mathit{n}}$ with $a=\mathrm{0},\mathrm{0.1},\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.

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 Effect of a variation of the initial value for the flange hole-pattern fit on the calculated SQP solution and its performance with DS1.

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 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⁰ 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 φ⁰ 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 ${\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}(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.

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

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

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⁻¹⁴. 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.

7 Conclusions

Minimax-type fitting applications in coordinate metrology such as Chebyshev, minimum-circumscribed, maximum-inscribed 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⁻¹⁴ 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).

Appendix A: Test data set properties

The rotation angles R_x, R_y and R_z in Table A2 refer to a rotation of the gauge (virtual counterpart of hole-pattern 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.

Table A1 Amount of measurement points per test data set.

Table A2 Reference solutions for hole-pattern fit of test data sets.

Appendix B: SQP solver parameters

Table B1 SQP method parameters for test calculations.

Appendix C: Nomenclature

∥v∥	Euclidian norm of vector v
In	Unit matrix of dimension n
∇v	Gradient operator for variables v
(⋅)T	Transposed vector or matrix