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

An SQP method for Chebyshev and hole-pattern fitting

An SQP method for Chebyshev and hole-pattern fitting with geometrical elementsAn SQP method for Chebyshev and hole-pattern fittingDaniel Hutzschenreuter et al.

Daniel Hutzschenreuter^{1},Frank Härtig^{1},and Markus Schmidt^{2}Daniel Hutzschenreuter et al. Daniel Hutzschenreuter^{1},Frank Härtig^{1},and Markus Schmidt^{2}

^{1}Abteilung Mechanik und Akustik, Physikalisch-Technische Bundesanstalt
(PTB), Braunschweig, 38116, Germany

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.

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:

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

$$}{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

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, $\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.

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

(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 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}\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 Armijo-type step width is calculated from the
exact l_{1}-barrier-penalty 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 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.

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.

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

(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}.$

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 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}\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 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

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

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 non-parallel 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 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

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 3Associated 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 4Hole-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

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

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},\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 CT-sensor-based 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.

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

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.

Figure 6Deviation 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^{−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}_{\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},\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.

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

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^{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 ${\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
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 9Calculation 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^{−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.

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^{−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).

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 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 A1Amount of measurement points per test data set.

This paper and its content, numbers and figures were
authorized for publication by M. Schmidt, who is co-author 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 MNPQ-Transfer
project “3D-Lochbildeinpassung für Koordinatenmesssysteme mit taktilen
und optischen Sensoren sowie CT-Systemen”.

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 best-fit geometric elements, Centrum
voor Wiskunde en Informatica, Report NM-R9317, 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:2007-04, 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.

Hutzschenreuter, D.: Benchmarking data sets for flange hole pattern fitting,
Physikalisch-Technische 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
TraCIM-system, 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.

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.