Evaluating the form errors of spheres,cylinders and cones using the primal

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.

Keywords

Form error interior point method minimum zone elements maximum inscribed elements minimum circumscribed elements

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

min_{x} f (x) s . t . d_{j} (x) \geq 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 (x) = f (x) - \sum_{j = 1}^{N} ω_{j} d_{j} (x)

(2)

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

\begin{matrix} g (x) - B {(x)}^{T} ω = 0 & Stationarity \\ d (x) \geq 0 & Primal feasibility \\ ω \geq 0 & Dual feasibility \\ D (x) ω \geq 0 & Complementary slackness \end{matrix}

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

A (x) = {j | d_{j} (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} (x^{*}) = 0, \forall j \in A$ ; while for inactive indices, the complementary slackness condition suggests that $ω_{j}^{*} = 0, \forall j \notin 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 $ξ \in (0, 1 / 2), η \in (0, 1), ν > 2, θ \in (0, 1), τ \in (2, 3), κ \in (0, 1)$ and thresholds $ω_{min}$ , $ω_{max}$ and ε. Obtain an initial solution x₀ and then solve ω₀ from $min_{ω'} ‖ g (x_{0}) - B (x_{0}) ω' ‖$ and $ω_{0} = min (0.1, ω')$ .

Step 2. Compute incremental direction

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

[\begin{matrix} - W & B^{T} \\ diag (ω) B & D \end{matrix}] [\begin{matrix} Δ x \\ Δ ω \end{matrix}] = [\begin{matrix} g (x) - B^{T} ω \\ μ - D ω \end{matrix}]

(3)

with µ = 0 and $W \approx \nabla^{2} f (x) - \sum_{j = 1}^{N} ω_{j} \nabla^{2} d_{j} (x)$ .

If $‖ Δ x ‖ \leq ε$ and $ω + Δ ω \geq 0$ , an optimum is found, then stop.

Set $φ : = min {max [0, - ω - Δ ω - 10^{3} d], 1}$ and $γ : = 〈 g, Δ x 〉 + \sum_{j = 1}^{N} \frac{ω_{j} + Δ ω_{j}}{ω_{j}} ϕ_{j}$

Set $μ = (1 - ψ) φ + ρ ψ ω$ with $ρ : = ‖ Δ x ‖^{ν} + ‖ φ ‖$ and

ψ : = {\begin{matrix} 1 & \sum_{j = 1}^{N} \frac{ω_{j} + Δ ω_{j}}{ω_{j}} (ρ ω_{j} - φ_{j}) \leq 0 \\ min {1, \frac{(1 - θ) | γ |}{\sum_{j = 1}^{N} \frac{ω_{j} + Δ ω_{j}}{ω_{j}} (ρ ω_{j} - φ_{j})}} & otherwise \end{matrix}

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

Set $I : = {j | d_{j} (x) \leq ω_{j}}$ .

Set δx be the solution of $min 〈 δ x, W δ x 〉, s . t . d_{j} (x + Δ x) + 〈 \nabla d_{j} (x), δ x 〉 = max (‖ Δ x ‖^{τ}, max_{j \in I} | ρ_{j} |^{κ} ‖ Δ x ‖^{2}), \forall j \in I$

If this is infeasible or unbounded or $‖ δ x ‖ > ‖ Δ x ‖$ , then set δx = 0.

Step 3. Arc search

Compute α, the first in the sequence {1, η, η², …} satisfying $f (x + α Δ x + α^{2} δ x) \leq f (x) + ξ α 〈 g (x), α Δ x 〉, \forall d_{j} (x + α Δ x + α^{2} δ x) > 0, j = 1, \dots, N$ .

Step 4. Update

$x \leftarrow x + α Δ x + α^{2} δ x$ and $ω \leftarrow min (max (max (‖ Δ x ‖^{2}, ω + Δ ω), ω_{min}), ω_{max})$

Go to Step 2.

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

〈 u, (W_{k} + \sum_{j} \frac{ω_{k}^{j}}{d_{j} (x_{k})} \nabla d_{j} (x_{k}) \nabla d_{j} {(x_{k})}^{T}) u 〉 \geq σ_{1} ‖ u ‖^{2}, \forall u

and

‖ W_{k} ‖ \leq σ_{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

min_{x} max_{j} | e_{j} (x) | \Leftrightarrow min_{x} D s . t . D \geq | e_{j} (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

min_{x} R s . t . \sqrt{x_{j}^{2} + y_{j}^{2}} \geq R, j = 1, 2, \dots, N

Obviously, it is equivalent to

min_{x} D s . t . 0 \leq \sqrt{x_{j}^{2} + y_{j}^{2}} - R \leq 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 ϕ$ and $s = \sin ϕ$ .

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, t_x, t_y, t_z }	D≥ (l−r)²	D+l≥ 0	D−l≥ 0
Cylinder	{D, r, θ_x, θ_y, t_x, t_y }	D≥ (h−r)²	D+h≥ 0	D−h≥ 0
Cone	{D, ϕ, θ_x, θ_y, t_x, t_y, t_z}	D≥ (ch−sz)²	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¹⁸	8.327	8.327	1.0035	1.0035	8.346	0.9959	0.9960	8.365
Sphere 2 in Chen and Liu¹⁸	9.67	9.67	50.0379	50.0379	12.59	50.0302	50.0302	9.93
Cylinder 1 in Lai et al.¹⁹	2.788	2.788	12.0002	12.0002	2.826	12.0016	12.0016	2.788
Cylinder 3 in Carr and Ferreira²⁰	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⁷	0.0032 in	0.0032 in	0.0042 in	0.0038 in
Cone 2 in Huang and Lee²¹	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.

Footnotes

Appendix 1 Funding

This study was supported by the European Research Council through its program ERC-2008-AdG 228117-Surfund and National Natural Science Foundation of China through its grant 51205064.

References

Whitehouse

. Surfaces and their measurement. London: Hermes Penton Science, 2002.

ISO 1101:2004. Geometrical product specifications – geometrical tolerancing – tolerances of form, orientation, location and run-out.

Chetwynd

. Applications of linear programming to engineering metrology. Proc IMechE, Part B: J Engineering Manufacture 1985; 199(2): 93–100.

Kanada

. Evaluation of spherical form errors: computation of sphericity by means of minimum zone method and some examinations with using simulated data. Precis Eng 1995; 17(4): 281–289.

Goch

Lübke

. Tschebyscheff approximation for the calculation of maximum inscribed/minimum circumscribed geometry elements and form deviations. CIRP Ann: Manuf Techn 2008; 57: 517–520.

Al-Subaihi

Watson

. Fitting parametric curves and surfaces by l∞ distance regression. BIT 2005; 45: 443–461.

Chatterjee

Roth

. Chebyshev approximation methods for evaluating conicity. Measurement 1998; 23(3): 63–76.

Huang

. An exact minimum zone solution for sphericity evaluation. Comput Aided Design 1999; 31: 845–853.

Venkaiah

Shunmugam

. Evaluation of form data using computational geometric techniques – part II: cylindricity error. Int J Mach Tool Manu 2007; 47: 1237–1245.

10.

Devillers

Preparata

. Culling a set of points for roundness and cylindricity evaluation. Int J Comput Geom Ap 2003; 13: 231–240.

11.

Zhang

Jiang

Scott

. A reliable method of minimum zone evaluation of cylindricity and conicity from coordinate measurement data. Precis Eng 2011; 35: 484–489.

12.

Sharma

Rajagopal

Anand

. A genetic algorithm based approach for robust evaluation of form tolerances. J Manuf Syst 2000; 19(1): 46–57.

13.

Nocedal

Wright

. Numerical optimization. New York: Springer, 1999.

14.

Bakhtiari

Tits

. A simple primal-dual feasible interior-point method for nonlinear programming with monotone descent. Comput Optim Appl 2003; 25: 17–38.

15.

Tits

Wächter

Bakhtiari

. A primal-dual interior-point method for nonlinear programming with strong global and local convergence properties. ISR Tech Report TR 2002-29, 19 July, 2002. College Park, MD: Institute for Systems Research, University of Maryland.

16.

Atieg

Watson

. A class of methods for fitting a curve or surface to data by minimizing the sum of squares of orthogonal distances. J Comput Appl Math 2003; 158(2): 277–296.

17.

BS EN ISO 2692:2006. Geometrical product specifications – geometrical tolerancing – maximum material requirement, least material requirement and reciprocity requirement.

18.

Chen

Liu

. A study on analyzing the problem of the spherical form error. Precis Eng 2000; 24: 119–126.

19.

Lai

Jywe

Chen

. Precision modeling of form errors for cylindricity evaluation using genetic algorithms. Precis Eng 2000; 24(4): 310–319.

20.

Carr

Ferreira

. Verification of form tolerances – part II: cylindricity and straightness of a median line. Precis Eng 1995; 17(2): 144–156.

21.