Evaluating the form errors of spheres,cylinders and cones using the primal–dual interior point method

Abstract

In precision metrology, the form errors between the measured data and reference nominal surfaces are usually evaluated in four approaches: least squares elements, minimum zone elements, maximum inscribed elements and minimum circumscribed elements. The calculation of minimum zone element, maximum inscribed element and minimum circumscribed element is not smoothly differentiable, thus very difficult to be solved. In this article, a unified method is presented to evaluate the form errors of spheres, cylinders and cones in the sense of minimum zone element, maximum inscribed element and minimum circumscribed element. The primal–dual interior point method is adopted to solve this non-linearly constrained optimisation problem. The solution is recursively updated by arc search until the Karush–Kuhn–Tucker conditions are satisfied. Some benchmark data are employed to demonstrate the validity and superiority of this method. Numerical experiments show that this optimisation algorithm is computationally efficient and its global convergence can be guaranteed.

Introduction

In precision engineering, form errors are crucial for their functionalities, for example, vibration of shafts and wear of contact elements. Form errors are defined as the relative deviation between the measured data and the reference nominal surfaces. The reference surfaces are usually fitted from the discrete data points in the sense of least squares elements (LSEs), minimum zone elements (MZEs), maximum inscribed elements (MIEs) and minimum circumscribed elements (MCEs). ¹ The MZE encloses all the measured data points between two coincident geometry elements. Among the four fitting criteria, the MZE approach can obtain the smallest form error, and it complies with the definition of form error in ISO 1101, ² which is consistent with the tolerance design. In contrast, the MCE and MIE are realised by one-sided fitting, which constraints all the data points to be outside the inscribed geometry elements or inside the circumscribed elements. MCE is suitable for the convex geometry with the maximum material condition, while the MIE is suitable for the geometry with the maximum material condition of a concave shape. The calculation of MZE, MIE and MCE is a constrained non-linear optimisation problem and very difficult to be solved. As a consequence, most commercial precision instruments conduct least squares fitting, and the form errors in the sense of MZE are evaluated by subtracting the maximum and minimum deviations from the data points to the fitted surface. But this will lead to overestimation, and unnecessary rejections of qualified parts may occur. Therefore, appropriate objective functions should be employed for the evaluation of form errors, and proper solving algorithms need to be developed accordingly.

Literature review

Spheres, cylinders and cones are ubiquitous in the shapes of industrial components. In fact, more than 70% of the functional parts in precision machines can be modelled using these simple shapes. Researchers have developed various methods to evaluate the form qualities of simple geometries in the sense of MZE.

The minimum zone fitting of circles and spheres are quadratic programming problems. Chetwynd ³ and Kanada ⁴ iteratively linearised the formulation and solved it using the simplex technique. Goch and Lübke ⁵ approximated the Chebyshev norm fitting by lp norm fitting with 50 ≤ p ≤ 100. Obviously, it is equivalent to Chebyshev norm fitting when p →∞, but when p grows larger, the observation matrix will be seriously ill conditioning, and the solution becomes very unstable. Al-Subaihi and Watson ⁶ also recursively approximated the minimax problem by linearisation, and the new optimisation problem is solved using the Gauss–Newton algorithm. All these approximation techniques suffer from slow convergence, especially when the solutions approach the regions of optima.

Rather than approximation, much effort has been made to directly solve the minimum zone fitting problems. Computational geometry–based methods have drawn extensive attention, such as vertex points, ⁷ Voronoi diagrams, ⁸ convex hulls, ⁹ α-hulls ¹⁰ and so on. However, these methods have serious limitations in restricted applicability. They require analysing the specific geometry conditions in each individual case for different shapes and error metrics, thus very inconvenient to apply. Moreover, some heuristic search methods, such as differential evolution ¹¹ and genetic algorithm, ¹² can also be adopted. These techniques are easy to implement and able to solve very complex optimisation problems, but they suffer from high computational complexity. In addition, the solutions are not deterministic, that is, the solution obtained at each run varies due to the randomness introduced in the evolution process.

Concerning the MIE and MCE fitting of spheres, cylinders and cones, little work can be found in literature, especially for cones. Goch and Lübke ⁵ first carried out Chebyshev norm fitting by reducing one degree of freedom, say, radius of the cylinder, and then the radius was adjusted to ensure that all the data points are inside the circumscribed element or outside the inscribed element. Similar to minimum zone fitting, the solutions attained by this approach are approximation of the real one, and the convergence rate is very slow. In fact, heuristic searching algorithms can be directly applied here, but they suffer from the same problems with the minimum zone fitting mentioned earlier.

In this article, the fitting problems are converted into a differentiable optimisation problems with non-linear constraints. Then, a unified solving algorithm, called the primal–dual interior point (PDIP) method, can be employed. As a consequence, the fitting of different shapes with different optimisation criteria can be accomplished in the same procedure, just with some extra operations of selecting appropriate objective functions and constraints.

The PDIP algorithm

For an optimisation problem with inequality constraints

$$\underset{\mathbf{x}}{min}f(\mathbf{x})\text{s}.\text{t}.d_j(\mathbf{x})\ge 0,j=1,2,\dots ,N$$

(1)

where f: R ⁿ → R and {d_j: R ⁿ → R, j = 1, …, N} are twice continuously differentiable functions. The generalised Lagrangian function of this constrained optimisation problem is

$$L(\mathbf{x})=f(\mathbf{x})-\sum _{j=1}^N{\omega }_j{d}_j(\mathbf{x})$$

(2)

where {ω_j≥ 0} are generalised Lagrangian multipliers. The first-order necessary optimality conditions of the constrained optimisation problem are ¹³

$$\begin{array}{cc}\mathbf{g}(\mathbf{x})-\mathbf{B}{(\mathbf{x})}^T\mathbf{\omega }=0\hfill & \text{Stationarity}\hfill \\ \multicolumn{1}{c}{\mathbf{d}(\mathbf{x})\ge 0}& \text{Primal feasibility}\hfill \\ \multicolumn{1}{c}{\mathbf{\omega }\ge 0}& \text{Dual feasibility}\hfill \\ \multicolumn{1}{c}{\mathbf{D}(\mathbf{x})\mathbf{\omega }\ge 0}& \text{Complementary slackness}\hfill \end{array}$$

where g(x) := ∂f(x)/∂x, B(x) := ∂d(x)/∂x and D(x) := diag(dj(x)).

These optimality conditions are also called the Karush–Kuhn–Tucker (KKT) conditions. The solution x* is called stationary if there exists a weighting vector ω, such that the KKT conditions hold. The active set is defined as the set of the indices of the active inequality constraints

$$\text{A}(\mathbf{x})=\{j|d_j(\mathbf{x})=0,j=1,\dots ,N\}$$

The multipliers ω measure the sensitivity of the optimal objective function f(x) against constraints {d_j(x)}. At optimum x*, for active constraints $d_j({\mathbf{x}}^{*})=0,\forall j\in \text{A}$ ; while for inactive indices, the complementary slackness condition suggests that ${\omega }_j^{*}=0,\forall j\notin \text{A}$ , which indicates that f(x) is not affected by these constraints.

A PDIP method can be applied to find the optimal solution x*. ¹⁴ This method needs an initial feasible solution x₀, that is, an interior point (IP) within the feasible domain restricted by the constraints. The objective function f(x) is decreased by iteratively solving the Lagrangian function of the optimisation problem. A correct solution is thought to be found when the KKT conditions are satisfied. In this method, x are called primal variables and ω are called dual variables. Sometimes the non-negativity of ω does not hold, but the convergence of the optimisation problem will not be influenced. The procedures of iteratively updating the solution x and vector ω are shown in Algorithm 1.

Algorithm 1. Procedure of the PDIP algorithm

• Step 1. Initialise parameters

Set $\xi \in (0,1/2),\eta \in (0,1),\nu >2,\theta \in (0,1),\tau \in (2,3),\kappa \in (0,1)$ and thresholds ${\omega }_{min}$ , ${\omega }_{max}$ and ε. Obtain an initial solution x₀ and then solve ω₀ from $\underset{\omega \prime }{\text{min}}\Vert \mathbf{\text{g}}({\mathbf{x}}_0)-\mathbf{\text{B}}({\mathbf{x}}_0)\mathbf{\omega }\prime \Vert$ and ${\omega }_0=min(0.1,\omega \prime )$ .

• Step 2. Compute incremental direction

1. Compute (Δx, Δω) by solving the Lagrangian function L(x, ω, W, µ)

   $$\left[\begin{array}{cc}\hfill -\mathbf{W}\hfill & {\mathbf{B}}^T\\ \multicolumn{1}{c}{\text{diag}(\mathbf{\omega })\mathbf{B}}& \mathbf{D}\end{array}\right]\left[\begin{array}{c}\hfill \Delta \mathbf{x}\hfill \\ \multicolumn{1}{c}{\Delta \mathbf{\omega }}\end{array}\right]=\left[\begin{array}{c}\hfill \mathbf{g}(\mathbf{x})-{\mathbf{B}}^T\mathbf{\omega }\hfill \\ \multicolumn{1}{c}{\mathbf{\mu }-\mathbf{D}\mathbf{\omega }}\end{array}\right]$$

   (3)

   with µ = 0 and $\mathbf{W}\approx {\nabla }^{2}f(\mathbf{x})-\sum _{j=1}^N{\omega }_j{\nabla }^2{d}_j(\mathbf{x})$ .

   If $\Vert \Delta \mathbf{x}\Vert \le \epsilon$ and $\mathbf{\omega }+\mathbf{\Delta }\mathbf{\omega }\ge 0$ , an optimum is found, then stop.

2. Set $\phi :=min\{max[0,-\mathbf{\omega }-\mathbf{\Delta }\mathbf{\omega }-{10}^3\mathbf{d}],1\}$ and $\gamma :=〈\mathbf{g},\Delta \mathbf{x}〉+{\displaystyle {\sum }_{j=1}^N\frac{{\omega }_j+\Delta {\omega }_j}{{\omega }_j}{\varphi }_j}$

   Set $\mathbf{\mu }=(1-\psi )\phi +\rho \psi \mathbf{\omega }$ with $\rho :=\Vert \Delta \mathbf{x}{\Vert }^{\nu }+\Vert \phi \Vert$ and

   $$\psi :=\{\begin{array}{cc}\hfill 1\hfill & \sum _{j=1}^N\frac{{\omega }_j+\Delta {\omega }_j}{{\omega }_j}(\rho {\omega }_j-{\phi }_j)\le 0\\ \multicolumn{1}{c}{min\{1,\frac{(1-\theta )|\gamma |}{\sum _{j=1}^N\frac{{\omega }_j+\Delta {\omega }_j}{{\omega }_j}(\rho {\omega }_j-{\phi }_j)}\}}& \text{otherwise}\end{array}$$

3. Compute (Δx, Δω) by solving L(x, ω, W, µ) in equation (3).

4. Set $I:=\{j|d_j(\mathbf{x})\le {\omega }_j\}$ .

5. Set δx be the solution of $min〈\delta \mathbf{x},\mathbf{W}\delta \mathbf{x}〉,\text{s}.\text{t}.d_j(\mathbf{x}+\Delta \mathbf{x})+〈\nabla d_j(\mathbf{x}),\delta \mathbf{x}〉=max(\Vert \Delta \mathbf{x}{\Vert }^{\tau },\underset{j\in I}{max}|{\rho }_j{|}^{\kappa }\Vert \Delta \mathbf{x}{\Vert }^2),\forall j\in I$

If this is infeasible or unbounded or $\Vert \delta \mathbf{x}\Vert >\Vert \Delta \mathbf{x}\Vert$ , then set δx = 0.

• Step 3. Arc search

Compute α, the first in the sequence {1, η, η², …} satisfying $f(\mathbf{x}+\alpha \Delta \mathbf{x}+{\alpha }^2\delta \mathbf{x})\le f(\mathbf{x})+\xi \alpha 〈\mathbf{g}(\mathbf{x}),\alpha \Delta \mathbf{x}〉,\forall d_j(\mathbf{x}+\alpha \Delta \mathbf{x}+{\alpha }^2\delta \mathbf{x})>0,j=1,\dots ,N$ .

• Step 4. Update

$\mathbf{x}\leftarrow \mathbf{x}+\alpha \Delta \mathbf{x}+{\alpha }^2\delta \mathbf{x}$ and $\mathbf{\omega }\leftarrow min(max(max(\Vert \mathbf{\Delta }\mathbf{x}{\Vert }^2,\mathbf{\omega }+\mathbf{\Delta }\mathbf{\omega }),{\mathbf{\omega }}_{min}),{\mathbf{\omega }}_{max})$

Go to Step 2.

The global convergence condition for the sequence of solution $\{{\mathbf{x}}_k\}$ is: given any index set K such that the sequence $\{{\mathbf{x}}_k|k\in K\}$ is bounded, there exist σ₁, σ₂ > 0 such that for all $k\in K$ , k≥ 1 ¹⁵

$$〈\mathbf{u},({\mathbf{W}}_k+\sum _j\frac{{\omega }_k^j}{d_j({\mathbf{x}}_k)}\nabla d_j({\mathbf{x}}_k)\nabla d_j{({\mathbf{x}}_k)}^T)\mathbf{u}〉\ge {\sigma }_1\Vert \mathbf{u}{\Vert }^2,\forall \mathbf{u}$$

and

$$\Vert {\mathbf{W}}_k\Vert \le {\sigma }_2$$

Therefore, at each iterate, it needs to check whether the matrix W+B ^T diag{ω_j/d_j}B is positive definite. If not, an appropriate diagonal matrix will be added onto it to guarantee the convergence of the solution.

One major superiority of this algorithm over other gradient-based iterative methods like sequential linear/quadratic programming is, the configuration parameters can be flexibly adapted according to the numerical properties of specific optimisation problems. Global optima can be obtained by relatively low computational cost.

Mathematical formulation of form error evaluation

In the field of metrology, the representations of ‘reference surfaces’ are conventionally adjusted to fit the measured data. This approach is convenient for planes and spheres since their mathematical functions keep similar forms after translation and rotation. But for cylinders and cones, the representations will become unbearably complicated when the surfaces are not at the ‘standard positions’. It has been strictly proved that when the form deviations are measured in the normal directions of the reference surfaces, moving the data is equivalent to transforming the fitted surfaces, ¹⁶ and the evaluated results of the form errors are the same. As a consequence, transformations are always performed on the measured data for the sake of numerical simplicity.

For the minimum zone fitting problem, symmetric tolerances are considered here only, thus MZE fitting is a Chebyshev norm fitting problem. It can be converted into a constrained optimisation problem very straightforwardly

$$\underset{\mathbf{x}}{min}\underset{j}{max}|e_j(\mathbf{x})|\iff \underset{\mathbf{x}}{min}D\text{s}.\text{t}.D\ge |e_j(\mathbf{x})|,j=1,2,\dots ,N$$

with e_i denoting the Euclidean distance from the data point p _i to the fitted surface.

The MIE and MCE fitting originate from the maximum material requirement on mating components in assembly of parts. ¹⁷ For spheres and cylinders, the MIE fitting is to maximise the radii of inscribed spheres/cylinders and the MCE fitting is to minimise the radii of circumscribed spheres/cylinders. But there is no such a simple intrinsic characteristic (shape parameter) to be taken as an objective function to indicate the sizes of cones. Therefore, an alternative approach is needed to perform the MIE/MCE fitting of cones.

We return to the MIE fitting of cylinders

$$\underset{\mathbf{x}}{min}R\text{s}.\text{t}.\sqrt{x_j^2+y_j^2}\ge R,j=1,2,\dots ,N$$

Obviously, it is equivalent to

$$\underset{\mathbf{x}}{min}D\text{s}.\text{t}.0\le \sqrt{x_j^2+y_j^2}-R\le D,j=1,2,\dots ,N$$

Here, $e_j=\sqrt{x_j^2+y_j^2}-R$ is the signed orthogonal distance from data point p _i to the fitted surface. So that the MIE fitting is converted to minimising the largest signed distance from the data points to the fitted surface while keeping all the distances positive. Similarly, the MCE fitting is converted to maximising the smallest signed distance from the data points to the fitted surface while keeping all the distances negative. This approach is generalised to the cases of circles, spheres and cones, by which the concept of ‘the size of a feature’ can be avoided.

An extra variable D is employed as the objective function. The variables and constraints of MZE, MIE and MCE are listed in Table 1 for spheres, cylinders and cones, respectively. Here, $l=\sqrt{x^2+y^2+z^2}$ , $h=\sqrt{x^2+y^2}$ , $c=\cos \varphi$ and $s=\sin \varphi$ .

Table 1.

Mathematical formulation of form error evaluation of the IP method.

Shape	X	d of MZE	d of MIE	d of MCE
Sphere	{D, r, tx, ty, tz }	D≥ (l−r)2	D+l≥ 0	D−l≥ 0
Cylinder	{D, r, θx, θy, tx, ty }	D≥ (h−r)2	D+h≥ 0	D−h≥ 0
Cone	{D, ϕ, θx, θy, tx, ty, tz}	D≥ (ch−sz)2	D≥ch−sz≥ 0	D≥sz−ch

MZE: minimum zone element; MCE: minimum circumscribed element; MIE: maximum inscribed element.

Numerical validation

The effectiveness of this algorithm is demonstrated by several benchmark data sets seen in literature. Here two spheres, two cylinders and two cones are employed. The best fitting results found in literature are also presented for comparison. The form evaluation programs using the PDIP algorithm are coded with MATLAB and run in a personal computer (PC).

The evaluation results of the spheres and cylinders are shown in Table 2. Chen and Liu ¹⁸ evaluated the form errors of spheres by identifying the distribution of the vertex points and then solving simultaneous linear equations. The cylindricity problem was solved by genetic algorithms in the study by Lai et al. ¹⁹ and by sequential linear programming in the study by Carr and Ferreira. ²⁰ It can be seen that the proposed algorithm can always obtain better or equivalently good results with the techniques in literature. It is worth noting that for the maximum inscribed evaluation of Sphere 1 given in the study by Chen and Liu, ¹⁸ the obtained radius of the IP method is greater than the reference, which implies that the PDIP algorithm overcomes the local minimum problem of Chen et al.’s algorithm.

Table 2.

Comparison on the evaluation results of spheres and cylinders

Data set	MZE form error (µm)		MCE radius (mm)		MCE form error (µm)	MIE radius (mm)		MIE form error (µm)
	Literature	IP method	Literature	IP method		Literature	IP method
Sphere 1 in Chen and Liu 18	8.327	8.327	1.0035	1.0035	8.346	0.9959	0.9960	8.365
Sphere 2 in Chen and Liu 18	9.67	9.67	50.0379	50.0379	12.59	50.0302	50.0302	9.93
Cylinder 1 in Lai et al. 19	2.788	2.788	12.0002	12.0002	2.826	12.0016	12.0016	2.788
Cylinder 3 in Carr and Ferreira 20	9.931	9.931	50.004	50.004	10.652	49.996	49.996	9.931

MZE: minimum zone element; MCE: minimum circumscribed element; MIE: maximum inscribed element; IP: interior point.

Table 3 presents the results of cone fitting. Chatterjee and Roth ⁷ adopted a simplex search method for the vertices to achieve an optimum for Chebyshev fitting. Huang and Lee ²¹ implemented the form error evaluation using the minimum potential energy algorithm, which is stimulated by the techniques in mechanics. The calculated results of MZE by the PDIP method is equivalent to the results obtained by Chatterjee and Roth ⁷ and can obtain a little better result for Cone 2 in the study by Huang and Lee. ²¹ As far as the authors know, there are no publications concerning the quantitative evaluation of form errors of cones in the sense of MCE and MIE yet, thus comparison cannot be undertaken for them. The results shown in Table 3 are to be verified and compared for the readers in the future.

Table 3.

Comparison on the evaluation results of cones.

Data set	MZE in literature	IP method
		MZE	MCE	MIE
Cone in Chatterjee and Roth 7	0.0032 in	0.0032 in	0.0042 in	0.0038 in
Cone 2 in Huang and Lee 21	4.590 µm	4.589 µm	6.264 µm	4.964 µm

MZE: minimum zone element; MCE: minimum circumscribed element; MIE: maximum inscribed element; IP: interior point.

It has to be stressed that the solutions of MZE, MCE and MIE fitting of spheres/cylinders/cones are not unique, ¹ which means that the locations and orientations of the fitted elements may differ even for the same data set. The initial solution x₀ of the PDIP optimisation program determines which optimum x* can be obtained. However, the form error, which is related to the objective function, should not be affected, as long as the global convergence can be guaranteed. Here, the initial position and shape parameters are calculated by the least squares orthogonal distance fitting, ²² and the initial value D in Table 1 is set as the smallest value satisfying the constraints for all the data points. The optimisation programs normally converge within 40 iterates, and the time taken is less than 0.5 s. Thus, this method is very computationally efficient and easy to implement.

Conclusions

In this article, form error evaluations based on MIE, MCE and MZE are discussed. These optimisation problems are solved using the PDIP method. During the iteration process, the Lagrangian functions are solved to work out the incremental directions of the solution and the auxiliary weighting vector. Arc search is employed to determine the step lengths. The Hessian matrix needs to keep positive definite to guarantee global convergence. Numerical experiments proved that the present algorithm works well for spheres, cylinders and cones. The obtained results are no worse than any techniques found in literature. The present method can obtain global optima by much less computational effort than sequential linear/quadratic programming techniques and heuristic searching algorithms. Compared to computational geometry methods, this algorithm is versatile and can be employed to the form error evaluation of various shapes by different error metrics, as long as the fitting problem can be converted into a differentiable optimisation problem with inequality constraints. Therefore, the PDIP method is very promising for evaluating the form errors of complex surfaces in the sense of MZE, MIE and MCE.