Minimum-Cost Ordering for Selective Assembly

Abstract

We consider the minimization of input cost for a selective assembly system which features two random inputs and a finite number of matching classes. This setup frequently arises in high-precision manufacturing when input tolerances are not tight enough for the required output precision. We determine optimality conditions for the cost-optimal input portfolio given an expected-output target, first using a normal approximation of the multinomial binning distribution, and second employing a simple upper envelope of the output objective. We show that the relative error tends to zero as the production scale becomes sufficiently large. The envelope optimization problem also yields a closed-form solution for the cost-minimizing inputs as well as total costs, which are easy to understand for managers. A numerical study tests the practicality of the envelope solution. The latter can be used as seed for a numerical solution of the exact problem, as well as a closed-form approximate solution which allows for an analysis of structural properties.

Keywords:

• Cost minimization
• selective assembly
• stochastic production

Disclaimer

As a service to authors and researchers we are providing this version of an accepted manuscript (AM). Copyediting, typesetting, and review of the resulting proofs will be undertaken on this manuscript before final publication of the Version of Record (VoR). During production and pre-press, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal relate to these versions also.

1 Introduction

By sorting low-precision inputs into predefined tolerance classes, so as to obtain a grade-specific matching for components of different types, it is possible to produce relatively high-precision output. This process of “selective assembly” has been known for more than 100 years—as part of interchangeable manufacturing. While much of the literature thus far has been concerned with optimizing the design of the matching bins, the available work on cost-effectively procuring the inputs for a given system of matching bins with arbitrary distributions of part characteristics is rather sparse. This paper provides such an analysis, characterizing the optimal order quantities for two different part types, together with a closed-form approximation of the optimal input quantities. The latter allows for structural insights and suggests a straightforward generalization to an arbitrary number of part types.

In the batch production of movements for mechanical watches, the oscillator—assembled from a hairspring and a balance wheel—must exhibit a well-defined resonance frequency within a small tolerance. 1 The accuracy of this output characteristic is key for the quality of the watch, and this is why many high-end horology companies have adopted a system of selective assembly, sorting hairsprings and balance wheels into different grades, so as to guarantee a sufficiently high clocking precision after joining any two components with matching grades. Naturally, when using selective assembly the probability that, say, of 100 hairsprings not all can be matched with 100 balance wheels is generally positive, and this mismatch likelihood increases in the number of matching classes. That is, a greater output precision comes at the price of a bigger loss of input components. Finding the cost-optimal input quantities needs to trade off the cost of the inputs against the benefits of the outputs, which would vary across the different assembly grades. Given two multinomial distributions for the input parts being in the various matching classes, we characterize the cost-minimizing input quantities so as to achieve a given expected output quantity. We also show that a good approximation of the output objective (with arbitrarily small relative error for sufficiently large production lots) leads to a constant ratio between inputs and outputs, taking the production process to be a weighted sum of perfectly complementary systems. The latter may then allow for the separation of input parts into “principal components” and “buffer components.” In mechanical watchmaking, for example, hairsprings turn out to be buffer components for the significantly more expensive balance wheels; cf. Sec. 4.

1.1 Literature

The idea of sorting inputs of a given type into different bins so as to be able to interchangeably match components across corresponding bins for different part types, also termed “selective assembly,” was discussed more than a century ago by Buckingham (1921) who notes: “[t]he chief purpose of manufacturing, by selective assembly (…), is the production of large quantities of duplicate parts as economically as possible, within such limits that they may be assembled without further machining” (p. 224). Our concern here is related to the requirement of “as economically as possible,” by specifically considering the tradeoff of input costs against output benefits. The latter may naturally depend on the specific “grade” of the output, that is, on the matching class from which inputs were taken to assemble a given unit of a finished product. Mansoor (1961) points out that in practice assembly processes with more than 20 matching classes can be observed, but that one of the main problems of selective assembly is evidently the possibility for mismatch. In contrast to our input-output view, which is driven by the input-cost minimization problem and an interest in the dependence of input cost on required output precision and the number of matching classes, Caputo and Di Salvo (2019) are interested in different manufacturing methods for which they disaggregate production-specific cost parameters including processing times, sorting costs, and so forth, which are essentially invariant once the production process has been fixed, leaving open the question about the optimal inputs.

There are numerous practical applications for selective assembly systems, such as hole-and-shaft combinations (Asha et al., 2008; Kannan and Pandian, 2021), sheet metal assembly (Rezaei Aderiani et al., 2019), flat-panel display manufacturing (Duenyas et al., 1997), camshafts in the automobile industry (Mease et al., 2004), as well as bimetal thermostats, Fortini’s clutch, and knuckle-joint assembly (Tan and Wu, 2012). Mease et al. (2004) note that cost imbalances may put into question an otherwise yield optimized selective assembly system. In the case where input prices differ substantially, as for hairsprings versus the at least 10-times more expensive balance wheels, the less costly parts should be ordered in excess so that the expensive parts are paired more often.

One main concern of the literature thus far has been the design of these matching classes, in the form of bins in the space of part characteristics which are often assumed to be normally distributed. For example, Pugh (1986) considers an equidistant partition of scalar part characteristics. Kwon et al. (1999) improve on the equidistant partition by determining matching classes so as to minimize a square deviation of the difference of the (assumed to be) commensurate parts from some target. While Matsuura (2011) allows for more general convex loss functions, the general problem with the loss-function approach is that it aggregates penalizations (by taking the expectation) irrespective of whether parts fit or not. In other words, under a set of matching classes which was designed using the penalization approach it is generally possible that parts do not match, even though they were taken from corresponding bins and therefore should match. We therefore follow here the robust binning approach by Weber (2021) which guarantees that parts across corresponding bins always match. This allows us to concentrate on the design of the minimum-cost input portfolio. We also mention in passing that in addition to the fixed-bin selective assembly discussed thus far there is the notion of a “direct” selective assembly which requires to keep track of every single part so as to be able to match parts one-to-one by solving a matching problem on a bipartite graph (Coullard et al., 1998). Since tracking very similar individual parts would often prove impractical and very costly in production settings with significant output volumes, we restrict attention to a selective assembly where bins are given, so that for each part type it is possible to specify a discrete probability distribution of parts being associated with a given bin. As noted earlier, we also assume that parts of different types can always be matched when taken from corresponding bins.

In Economics, the perfect complementarity of inputs for production, with transformation to outputs in constant proportions, goes back to Leontief (1941) who used it as a simplifying assumption to analyze a large economy. In assembly operations, constant proportions in the conversion of inputs to outputs arise because the required amount of input for a unit of output is typically fixed. Without any significant loss of generality, the different part types are assumed to be matched one-to-one, i.e., at equal proportions. For example, if each finished ball bearing requires 8 balls and 1 shaft, then we can think of the balls as being provided in lots of 8, so as to maintain equal proportions of inputs and outputs. Given that for any given matching class, the number of components of each type is ex-ante uncertain, the deterministic production function mapping inputs to expected outputs no longer exhibits perfect complementarity but rather a certain input substitutability, since additional quantities of type-1 parts (e.g., ‘balls’) provides an insurance against being long with type-2 inputs (e.g., ‘shafts’) in a given matching class. Thus, increasing the quantity of one input type must always lead to a nonzero increase in the expected output, despite the perfect complementarity applied to any realization of the “sorting product” (i.e., the vector of measured parts in the different matching classes). Our paper extends the work by Weber (2022) which focuses on the special case of binary random part matching, while we here allow for a selective assembly system with an arbitrary number of matching classes and also examine the scaling behavior.

2 Model

2.1 Matching Classes and Sorting Products

The selective assembly problem considered here involves components of two types (“1” and “2”) and $M\ge 2$ matching classes. A component of type $i\in \mathcal{I}=\{1,2\}$ features a (random) characteristic $S_i$, with realizations $s_i$ in the (measurable) space ${\mathcal{S}}_i$ which has at least M elements; see Sec. 3.7. The “matching classes” ${\mathcal{C}}_{i,m}$, for $m\in \mathcal{M}=\{1,\dots ,M\}$, are nonempty subsets which together form a partition of ${\mathcal{S}}_i$, in the sense that$$(m,\widehat{m}\in \mathcal{M},m\ne \widehat{m}\Rightarrow {\mathcal{C}}_{i,m}\cap {\mathcal{C}}_{i,\widehat{m}}=\varnothing )\mathit{and}\underset{m\in \mathcal{M}}{\cup }{\mathcal{C}}_{i,m}={\mathcal{S}}_i,$$

for all $i\in \mathcal{I}$. Let(1) $$p_{i,m}=\mathrm{\mathbb{P}}(S_i\in {\mathcal{C}}_{i,m})\in (0,1)$$(1)

be a given probability that a given component of type i has characteristics that lie in the matching class m. To avoid trivialities, we exclude the case where $p_{i,m}\in \{0,1\}$. Indeed, for $p_{i,m}=0$ no successful grade-m matching is possible and the class should be suppressed. Conversely if $p_{i,m}=1$, then no other matching class can have successful matches, contradicting the assumption of having at least two functional matching classes, the minimum number of classes for a meaningful selective assembly.

Two parts of characteristics $s_1$ and $s_2$ a referred to as “matched” (across types) if the vector $s=(s_1,s_2)\in \mathcal{S}={\mathcal{S}}_1\times {\mathcal{S}}_2$ is an element of the joint matching class ${\mathcal{C}}_m={\mathcal{C}}_{1,m}\times {\mathcal{C}}_{2,m}$, for some $m\in \mathcal{M}$. Since by assumption the set of matching classes $\{{\mathcal{C}}_{i,1},\dots ,{\mathcal{C}}_{i,M}\}$ forms a partition of the space of characteristics ${\mathcal{S}}_i$, all incoming parts can be sorted into matching classes; thus, taking two parts from corresponding matching classes ${\mathcal{C}}_{1,m}$ and ${\mathcal{C}}_{2,m}$ always produces a “matched pair;” see Fig. 1. In practice, the characteristics of different types of parts (e.g., hairsprings and balance wheels for mechanical watches) can vary substantially and be measured in different units (cf. footnote 1). When matching parts of corresponding grades, the number of feasible assemblies for any given matching class (referred to as “grade”) corresponds to the minimum number of parts available from either type. To clarify this strict complementarity, Fig. 1 depicts a situation where the smallest matching class (for each component type) produces zero output: indeed, the two type-1 components in ${\mathcal{C}}_{1,1}$ cannot be matched with any type-2 components, as ${\mathcal{C}}_{2,1}$ is empty.

Let $x_i\ge 0$ be the chosen input quantity for type-i components. Sorting them into matching classes yields the vector ${\mathit{x}}_i=(x_{i,1},\dots ,x_{i,M})$ (henceforth referred to as “sorting product”), where $x_{i,m}\ge 0$ denotes the number of components of type i in matching class m, so that$$x_i=\sum _{m\in \mathcal{M}}{x}_{i,m},i\in \mathcal{I}.$$

The multinomial probability that a size-$x_i$ ordering of type-i components into M matching classes produces a particular sorting product ${\mathit{x}}_i\in {\mathbb{N}}^{x_i}$ is$$\begin{array}{c}\mathrm{\mathbb{P}}({\mathbf{X}}_i={\mathit{x}}_i|{\mathbf{p}}_i,x_i)=\{\begin{array}{cc}\left(\prod _{m\in \mathcal{M}}\frac{p_{i,m}^{x_{i,m}}}{x_{i,m}!}\right)(x_i!),\hfill & \text{if}\mathrm{ }{x}_i=\sum _{m\in \mathcal{M}}{x}_{i,m},\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}\hfill \end{array}$$

where ${\mathbf{p}}_i=(p_{i,1},\dots ,p_{i,M})\in {(0,1)}^M$ is a given vector of positive likelihoods $p_{i,m}$ such that(2) $$\sum _{m\in \mathcal{M}}{p}_{i,m}=1,i\in \mathcal{I}.$$(2)

The preceding normalization constraint may be relaxed to include the possibility of wasteful off-spec components (cf. Remark 5).

2.2 Normal Approximation

We use the standard normal approximation of the multinomial distribution of the random vector ${\mathbf{X}}_i=(X_{i,1},\dots ,X_{i,M})$ with the mean ${\mathbf{\mu }}_i=({\mu }_{i,1},\dots ,{\mu }_{i,M})$, where ${\mu }_{i,m}=p_{i,m}{x}_i$, and the symmetric positive-definite covariance matrix ${\Sigma }_i\in {ℝ}^{M\times M}$. The latter features the diagonal elements ${\sigma }_{i,m}^2=\text{var}(X_{i,m})=x_i{p}_{i,m}(1-p_{i,m})$ and the off-diagonal elements $\text{cov}(X_{i,m},X_{i,\widehat{m}})=-x_i{p}_{i,m}{p}_{i,\widehat{m}}$ for $m,\widehat{m}\in \mathcal{M}$ with $m\ne \widehat{m}$. In what follows, we usually suppress the dependence on ${\mathbf{p}}_i$ for simplicity, as it is not a decision variable. Under the proposed normal approximation,$$\mathrm{\mathbb{P}}(X_{i,m}\le \xi |x_i)\approx G_{i,m}(\xi |x_i)=\Phi \left(\frac{\xi -{\mu }_{i,m}}{{\sigma }_{i,m}}\right),\xi \in [0,x_i],$$

denotes the cumulative distribution function (cdf) for the random number $X_{i,m}$ of type-i components in the m-th matching class, where the cdf of the standard normal distribution is$$\Phi (\xi )={\int }_{-\infty }^{\xi }\varphi (\zeta )d\zeta =\frac{1}{2}\left[1+\text{erf}\mathrm{ }(\xi /\sqrt{2})\right],\xi \in ℝ,$$

with probability density function (pdf)$$\varphi (\xi )=\frac{\hspace{0.17em}\exp \hspace{0.17em}\left(-{\xi }^2/2\right)}{\sqrt{2\pi }},\xi \in ℝ.$$

Clearly, the functions $G_{i,m}(·|x_i)$, for $m\in \mathcal{M}$, constitute marginal distributions of the joint cdf that describes the random sorting product of type-i components under the normal approximation, but it turns out that they are all that is required for our purposes.

Remark 1 (Approximation Error).

By the Berry-Esseen theorem (Berry, 1941; Esseen, 1942) the absolute deviation of the marginal distribution $G_{i,m}$ from the underlying binomial distribution can be bounded as follows:$$\underset{\xi \in [0,x_i]}{\sup }|\mathrm{\mathbb{P}}\left(X_{i,m}\le \xi |x_i\right)-G_{i,m}(\xi |x_i)|\le \left(\frac{p_{i,m}^2+{(1-p_{i,m})}^2}{{\sigma }_{i,m}}\right)\kappa ,$$

for all $(i,m)\in \mathcal{I}\times \mathcal{M}$, with $\kappa =0.4748$ (Shevtsova, 2011), which lies above the theoretical lower bound for $\kappa$ of $(\sqrt{10}+3)/\left(6\sqrt{2\pi }\right)\approx 0.4097$ (Esseen, 1956). Hence, as long as$$x_i>\left(\frac{\max \{p_{i,m},1-p_{i,m}\}}{\min \{p_{i,m},1-p_{i,m}\}}\right){\nu }^2,\nu \ge 0,$$

the approximation error corresponds to the probability of ‘extreme events,’ further than $\nu$ standard deviations away from the mean.

Remark 2 (Continuous Input Values).

The normal approximation allows for nonnegative input quantities $x_i$, lifting the burden of the integer constraint for the purposes of optimization. A practically implementable input vector can then usually be found ex post by rounding to the nearest integer or the nearest feasible batch size.

2.3 Matching Output

Given its type-1 and type-2 inputs, the firm produces the (sorted) output vector $\mathbf{Y}=(Y_1,\dots ,Y_M)$. The m-th random matching output $Y_m$, generated by using components with joint characteristics in matching class ${\mathcal{C}}_m$, is(3) $$Y_m=\min \{X_{1,m},X_{2,m}\},m\in \mathcal{M}.$$(3)

It is evident that the intrinsic value of the firm’s output may not be the same for all matching classes. For example, the matching-class index m can denote different quality grades in which case high-quality output would yield a higher price than low-quality output, without any difference in the unit cost for each random input component. To capture the differentiation of benefits across matching classes, let(4) $$Q=\mathrm{\eta }·\mathbf{Y}=\sum _{m\in \mathcal{M}}{\eta }_m{Y}_m$$(4)

denote the firm’s total (weighted) output, where $\mathrm{\eta }=({\eta }_1,\dots ,{\eta }_M)$ is a vector which contains the net unit value ${\eta }_m>0$ for the m-th matching output; see Fig. 2. The net unit value can be viewed as an absolute or a relative measure of benefits. When absolute, it may incorporate the revenues from selling the assembled grade-m product minus production costs, whereas when relative it serves to quantify a normalized net benefit. In the latter case, the normalizations ${\eta }_1+\cdots +{\eta }_m\in \{1,M\}$ are commonly used; for instance, by setting ${\eta }_m\equiv 1$ one obtains a simple physical output count, $Q=Y_1+\cdots +Y_m$, relevant for many (if not most) manufacturing applications. The assumption that ${\eta }_m$ is strictly positive reflects the firm’s concern for all units and its option of “free disposal.” If a certain matching class (e.g., $m=1$) represents “junk” components, then it is possible to reduce that weight reflecting any possible scrap value.

2.4 Expected Output

The need for grade-specific component matching renders the firm’s production stochastic, allowing generically for the possibility of even zero production output (with positive probability), irrespective of the chosen input quantities. While this may be true in principle, of course the likelihood of such poor output realizations becomes negligible in any realistic setting. For example, the likelihood of observing a type-i quantity of two standard deviations below the expected quantity in any given matching class is below 5%; for three standard deviations it is below $0.3\%$.

Since the firm is risk-neutral, its focus lies on its expected total (weighted) output,(5) $$F(\mathit{x})=\mathbb{E}\left[Q|\mathit{x}\right]=\sum _{m\in \mathcal{M}}{\eta }_m\mathbb{E}\left[Y_m|\mathit{x}\right]=\sum _{m\in \mathcal{M}}{\eta }_m{F}_m(\mathit{x}),\mathit{x}\in {ℝ}_{+}^2,$$(5)

which is a deterministic function of the corresponding input $\mathit{x}=(x_1,x_2)$, where the total output $F(·)$ is obtained as the weighted sum of the expected grade-m outputs $F_m(·)$. For any given m, the expected grade-m output depends on the statistical properties of the difference, ${\Delta }_m=X_{2,m}-X_{1,m}$, referred to as the (grade-m) “parts defect,” which has the mean$${\delta }_m={\mu }_{2,m}-{\mu }_{1,m},$$

and the standard deviation$${\sigma }_m=\sqrt{{\sigma }_{1,m}^2+{\sigma }_{2,m}^2}.$$

The expected output in any given matching class can be expressed as a convex combination of the expected type-1 and type-2 quantities, minus a positive loss term that increases superlinearly in the aforementioned standard deviation of the parts defect.

Lemma 1 (Grade-m Output).

Given any $m\in \mathcal{M}$, the expected grade-m output is given by (6) $$F_m(\mathit{x})={\Phi }_m{\mu }_{1,m}+(1-{\Phi }_m){\mu }_{2,m}-{\varphi }_m{\sigma }_m,$$(6)

for all $\mathit{x}\in {ℝ}_{+}^2$, where ${\Phi }_m=\Phi ({\delta }_m/{\sigma }_m)$ denotes the probability that the grade-m part defect reaches one standard deviation, and ${\varphi }_m=\varphi ({\delta }_m/{\sigma }_m)$ is the probability density at that point.

When the expected grade-m parts defect ${\delta }_m$ is positive (resp., negative), then the weight ${\Phi }_m\in [0,1]$ placed on the expected number ${\mu }_{1,m}$ of type-1 parts decreases (resp., increases) in ${\sigma }_m$, since with an increased coefficient of variation for parts defects the expected number of type-2 parts increases (resp., decreases) in relative importance. In other words, perhaps counterintuitively, the nominally expected matches (as measured by the convex combination of the expected type-specific yields ${\mu }_{1,m}$ and ${\mu }_{2,m}$) goes up in the dispersion of the grade-m parts defect. However, this effect is offset by an increased “matching loss” ${\sigma }_m{\varphi }_m$: independent of the sign of ${\delta }_m$, the latter grows in ${\sigma }_m$ superlinearly, for ${\varphi }_m>0$ also grows in ${\sigma }_m$.

Lemma 2 (Monotonicity).

Given any $m\in \mathcal{M}$, the expected grade-m output is increasing, so $$(\mathit{x}\ne \widehat{\mathit{x}}\mathit{and}\mathit{x}\le \widehat{\mathit{x}})\Rightarrow F_m(\mathit{x})<F_m(\widehat{\mathit{x}}),$$

for any two input vectors $\mathit{x},\widehat{\mathit{x}}\in {ℝ}_{+}^2$.

By increasing the amount of any input, the expected output increases over all matching classes. By increasing the type-i input from $x_i$ to ${\widehat{x}}_i>x_i$ the grade-m distribution $G_{i,m}(·|{\widehat{x}}_i)$ first-order stochastically dominates the grade-m distribution $G_{i,m}(·|x_i)$. Since first-order stochastic dominance is preserved by taking the minimum across the component types, the grade-m output ${\widehat{Y}}_m$ (with $x_i$ in $\mathit{x}$) first-order stochastically dominates the output $Y_m$ (with ${\widehat{x}}_i$ in $\widehat{\mathit{x}}$), all else equal, which in turn implies the (strict) monotonicity of the firm’s total output in its inputs.

3 Minimum-Cost Selective Matching

3.1 Problem Statement

Each unit of a type-i input is procured at the cost $c_i>0$, which implies that the firm is interested in minimizing the cost of generating a desired (expected) total output $q>0$. Thus, given the cost vector $\mathbf{c}=(c_1,c_2)$, the firm’s (minimum-cost) selective matching problem is to determine the optimal input ${\mathit{x}}^{*}(q)=(x_1^{*}(q),x_2^{*}(q))$ (also referred to as its “output-driven demand”) which achieves its minimum cost,(*) $$C(q)=\underset{\mathit{x}\in {ℝ}_{+}^2}{\min }\left\{\mathbf{c}·\mathit{x}\right\},\text{subject}\mathrm{ }\text{to:}F(\mathit{x})\ge q,$$(*)

for any output $q>0$. Clearly, the solution to the selective matching problem does not depend on the individual input costs, but only on their ratio $\gamma =c_1/c_2$.

3.2 Optimal Input

The firm’s output driven demand (or optimal input) ${\mathit{x}}^{*}(q)$ minimizes the cost of its inputs to achieve a given output objective $q>0$. It will come as no surprise that at the optimum the expected total output must exactly equal the output objective. This “output efficiency” is a consequence of the strict monotonicity of expected output in its inputs and the fact that inputs are costly at the margin. In addition, the technical rate of substitution between the two inputs must, at the optimum, be equal to the cost ratio $\gamma$. This provides optimality conditions for the otherwise nonconvex optimization problem ( $\ast$), summarized by the following result.

Theorem 1 .

For any expected total output $q>0$, a solution ${\mathit{x}}^{*}(q)$ to the minimum-cost selective matching problem ( $\ast$) satisfies output efficiency, (7) $$F({\mathit{x}}^{*}(q))=q;$$(7)

in addition, at $\mathit{x}={\mathit{x}}^{*}(q)$ the marginal rate of technical substitution between type-1 and type-2 components is equal to the cost ratio $\gamma =c_1/c_2$, so (8) $${\frac{\sum _{m\in \mathcal{M}}{\eta }_m{p}_{1,m}\left({\Phi }_m-\frac{1-p_{1,m}}{2{\sigma }_m}{\varphi }_m\right)}{\sum _{m\in \mathcal{M}}{\eta }_m{p}_{2,m}\left(1-{\Phi }_m-\frac{1-p_{2,m}}{2{\sigma }_m}{\varphi }_m\right)}|}_{\mathit{x}={\mathit{x}}^{*}(q)}=\gamma ,$$(8)

where ${\Phi }_m$ and ${\varphi }_m$, for all $m\in \mathcal{M}$, are specified in La. 1.

Remark 3 (Nonconvexity).

Generally, the output objective $F(\mathit{x})$ is can be nonconcave in the vector of inputs $\mathit{x}$. Indeed, in the simple case of two grades, where the firm only cares about grade-1 output, it is $M=2$ and $\mathit{\eta }=({\eta }_1,{\eta }_2)=(1,0)$. The corresponding Hessian of F is:$$\begin{array}{c}\text{det}\hspace{0.17em}{D}^{2}F=\text{det}\hspace{0.17em}\left[\begin{array}{cc}{\partial }_{x_1{x}_1}^{2}F& {\partial }_{x_1{x}_2}^{2}F\\ {\partial }_{x_1{x}_2}^{2}F& {\partial }_{x_2{x}_2}^{2}F\end{array}\right]=-\frac{p_{1,1}^2{p}_{2,1}^2{(2-p_{1,1}-p_{2,1})}^2}{8{\sigma }_1^4(\mathit{x})}{\varphi }^2\left(\frac{p_{2,1}{x}_2-p_{1,1}{x}_1}{{\sigma }_1(\mathit{x})}\right)<0,\hfill \end{array}$$

for all $\mathit{x}=(x_1,x_2)\in {ℝ}_{++}^2$. The preceding determinant of the Hessian is equal to the product of its eigenvalues, which must have therefore different signs. As a consequence the output objective is nonconvex in this case. Thus, the solution to the optimality conditions (7)–(8) may not be unique. As shown in Sec. 3.4, a simplified approximate output objective may yield a reasonable approximate solution (in closed-form).

3.3 Envelope Approximation

Instead of the normal approximation of the expected output F we now consider the envelope output,(9) $$\overline{F}(\mathit{x})=\sum _{m\in \mathcal{M}}{\eta }_m\min \left\{{\mu }_{i,m}:i\in \mathcal{I}\right\}=\sum _{m\in \mathcal{M}}{\eta }_m\min \left\{p_{i,m}{x}_i:i\in \mathcal{I}\right\},$$(9)

which is well-defined for all $\mathit{x}\in {ℝ}_{+}^2$. The envelope output $\overline{F}$ is an upper bound for F.

Lemma 3 .

$F(\mathit{x})\le \overline{F}(\mathit{x})$, for all $\mathit{x}\in {ℝ}_{+}^2$.

As the positive linear combination of minimum functions is concave, Jensen’s inequality implies that the expectation of the minimum is weakly smaller that the minimum of the expectations. An important question that arises of course is how far, for a given input $\mathit{x}$, the envelope $\overline{F}(\mathit{x})$ output may deviate from the reference output $F(\mathit{x})$? We thus consider the absolute deviation,$$R(\mathit{x})\triangleq |F(\mathit{x})-\overline{F}(\mathit{x})|=\overline{F}(\mathit{x})-F(\mathit{x})\ge 0,\mathit{x}\in {ℝ}_{+}^2,$$referred to as the envelope approximation error. The next result shows that this error increases only sublinearly with the inputs.

Lemma 4 (Approximation Error).

For any $\mathit{x}\in {ℝ}_{++}^2$ it is $$\frac{R(\mathit{x})}{\sigma (\mathit{x})}\le \frac{1}{\sqrt{2\pi }},$$

where $\sigma (\mathit{x})=\sum _{m\in \mathcal{M}}{\eta }_m{\sigma }_m(\mathit{x})$.

In other words, the envelope approximation error is at most proportional to the weighted sum of the grade-specific standard deviations ${\sigma }_m(\mathit{x})=\sqrt{p_{1,m}(1-p_{1,m})x_1+p_{2,m}(1-p_{2,m})x_2}$, using the given output weight ${\eta }_m$, for each grade $m\in \mathcal{M}$. Note also that the preceding error bound is tight when ${\eta }_m\equiv \text{const}.$, $p_{i,m}\equiv 1/2$, and $x_{1,m}\equiv x_{2,m}$. In a heterogeneous setting, the derived bound may be conservative as it presumes that $p_{1,m}{x}_1=p_{2,m}{x}_2$ holds for all $m\in \mathcal{M}$; cf. Weber (2022). While the absolute deviation between F and $\overline{F}$ may grow with the size of the inputs, in terms of the relative error,$$r(\mathit{x})=\frac{R(\mathit{x})}{F(\mathit{x})},\mathit{x}\in {ℝ}_{++}^2,$$

we obtain the following strong global approximation result.

Lemma 5 (Relative Approximation).

For any $\epsilon \in (0,1)$, there exists a finite minimal envelope output $\underset{¯}{q}(\epsilon )$, which is such that: 2$$\overline{F}(\mathit{x})\ge \underset{¯}{q}(\epsilon )\Rightarrow r(\mathit{x})\le \epsilon .$$

For small $\epsilon >0$, the minimal envelope output $\underset{¯}{q}(\epsilon )$ is $O({\epsilon }^{-2})$.

Instead of a minimal envelope output $\underset{¯}{q}(\epsilon )$ it is possible to specify a minimal input $\underset{¯}{x}(\epsilon )$, which is also $O({\epsilon }^{-2})$), since by the homogeneity of $\overline{F}$ it is proportional to $\underset{¯}{q}(\epsilon )$; cf. footnote 2 . Then $\min \{x_1,x_2\}\ge \underset{¯}{x}(\epsilon )$ also implies $r(\mathit{x})\le \epsilon$. A key driver of the relative approximation result in La. 5 is that, as a consequence of La. 4, the given upper bound of the relative error decreases with the inverse square-root of the minimal input (i.e., with ${\left(\min \{x_1,x_2\}\right)}^{-1/2}$). The given lower bounds for input or output are fairly conservative. In practice, it is preferable to perform an error correction based on the relative error at the solution to an envelope optimization problem.

3.4 Envelope Optimization

Based on the upper bound $\overline{F}$ for the production function introduced in Sec. 3.3, we now consider the envelope optimization problem,(**) $$\overline{C}(q)=\underset{\mathit{x}\in {ℝ}_{+}^2}{\min }\left\{\mathbf{c}·\mathit{x}\right\},\text{subject}\mathrm{ }\text{to:}\overline{F}(\mathit{x})\ge q,$$(**)

for any given target output $q>0$. By construction, this optimization problem is convex, and in order to find a solution to ( $\ast \ast$), which will be referred to as an optimal envelope input ${\overline{\mathit{x}}}^{*}$, it is useful to combine all grades which have the same ratio of likelihoods (across component types),(10) $$\{{\rho }_{\ell }:\ell \in \{1,\dots ,L\},{\rho }_1<\cdots <{\rho }_L\}=\left\{\frac{p_{1,m}}{p_{2,m}}:m\in \mathcal{M}\right\},$$(10)

to obtain $L\le M$ combined grades, where each such combined grade is characterized by a distinct likelihood ratio ${\rho }_{\ell }$. We can then rewrite the envelope output in Eq. (9) as follows:(11) $$\overline{F}(\mathit{x})=\sum _{\ell =1}^L{\widehat{\eta }}_{\ell }\min \left\{{\rho }_{\ell }{x}_1,x_2\right\},$$(11)

where$${\widehat{\eta }}_{\ell }=\sum _{m\in \mathcal{M}}{\mathbf{1}}_{\{{\rho }_{\ell }=p_{1,m}/p_{2,m}\}}{\eta }_m{p}_{2,m},\ell \in \{1,\dots ,L\},$$

are the corresponding combined weights. The alternative representation of the envelope output in Eq. (11) allows for a closed-form representation of the solution to the envelope optimization problem ( $\ast \ast$).

Theorem 2 .

For any $q>0$, a solution to the envelope optimization problem ( $\ast \ast$) is (12) $${\overline{\mathit{x}}}^{*}(q)={\overline{\mathit{\xi }}}^{{\ell }^{*}}q,$$(12)

where each ${\overline{\mathit{\xi }}}^{\ell }=({\overline{\xi }}_1^{\ell },{\overline{\xi }}_2^{\ell })$, for $\ell \in \{1,\dots ,L\}$, denotes a vertex of the iso-output contour for a unit output, with $${\overline{\xi }}_1^{\ell }={\left(\sum _{k=1}^L{\widehat{\eta }}_k\min \left\{{\rho }_k,{\rho }_{\ell }\right\}\right)}^{-1}\text{and}{\overline{\xi }}_2^{\ell }={\rho }_{\ell }{\overline{\xi }}_1^{\ell }.$$

The optimal combined grade, ${\ell }^{*}=\max \{\ell \in \{1,\dots ,L\}:s_{\ell }<\gamma \}$, determines the critical vertex ${\overline{\mathit{\xi }}}^{{\ell }^{*}}$, where the iso-output contour slopes $s_1,\dots ,s_L$ are such that $$0=s_1<s_{\ell }=\left({\overline{\xi }}_2^{\ell }-{\overline{\xi }}_2^{\ell -1}\right)/\left({\overline{\xi }}_1^{\ell -1}-{\overline{\xi }}_1^{\ell }\right),$$

for all $\ell \in \{2,\dots ,L\}$.

The optimal envelope solution ${\overline{\mathit{x}}}^{*}$ is such that the slope of the iso-cost curves lies in the subgradient, i.e., $-\gamma \in \partial \overline{F}({\overline{\mathit{x}}}^{*})$, and at the same time envelope output efficiency holds, in the sense that $\overline{F}({\overline{\mathit{x}}}^{*})=q$. The following example illustrates the result.

Example 1 .

Consider a selective-assembly system with $M=5$ matching classes featuring the matching-probability vectors ${\mathbf{p}}_1=(0.4,0.2,0.1,0.1,0.2)$ for parts of type 1 and ${\mathbf{p}}_2=(0.2,0.1,0.1,0.2,0.4)$ for parts of type 2. The unit-cost vector is $\mathbf{c}=(3,1)$, so $\gamma =3$. Since the set of likelihood ratios, $\{p_{1,m}/p_{2,m}:m\in \mathcal{M}\}=\{1/2,1,2\}$, features only $L=3$ distinct elements, by Eq. (10) the relevant likelihood ratios (in increasing order) are ${\rho }_1=1/2$, ${\rho }_2=1$, and ${\rho }_3=2$. Given a uniform output-weight vector $\mathit{\eta }=(1,1,1,1,1)$, the envelope output can be written as in Eq. (11), with the combined weights ${\widehat{\eta }}_1=0.6$ (for ${\rho }_1=p_{4,1}/p_{4,2}=p_{5,1}/p_{5,2}$), ${\widehat{\eta }}_2=0.1$ (for ${\rho }_2=p_{3,1}/p_{3,2}$), and ${\widehat{\eta }}_3=0.3$ (for ${\rho }_3=p_{1,1}/p_{1,2}=p_{2,1}/p_{2,2}$):$$\overline{F}(x_1,x_2)=0.6\min \{x_1/2,x_2\}+0.1\min \{x_1,x_2\}+0.3\min \{2x_1,x_2\}.$$

The corresponding unit-output vertices are ${\overline{\mathit{\xi }}}^1=(2,1)$, ${\overline{\mathit{\xi }}}^2=(10/7,10/7)$, and ${\overline{\mathit{\xi }}}^3=(1,2)$, resulting in the iso-output contour slopes $s_1=0<s_2=3/4<s_3=4/3$. Comparing the latter with the cost ratio $\gamma =3$, we find ${\ell }^{*}=\max \{\ell \in \{1,2,3\}:s_{\ell }<\gamma \}=3$, so that ${\overline{\mathit{\xi }}}^3$ becomes the critical vertex. By Thm. 2, the solution to the envelope optimization problem ( $\ast \ast$) is therefore ${\overline{\mathit{x}}}^{*}(q)=(1,2)q=(q,2q)$, for any $q>0$. Fig. 3 depicts the corresponding situation comparing, for $q=100$, the iso-F contour with the iso-$\overline{F}$ contour. We also note that to ${\ell }^{*}=3$ there correspond two critical grades, $m^{*}\in \{1,2\}$. The solution to the envelope optimization problem is consistent with minimizing input cost only with respect to the critical grades.

Remark 4 (Alternative Representation).

From a computational viewpoint, instead of determining the set $\{{\rho }_1,\dots ,{\rho }_L\}$ from the likelihood ratios $p_{1,m}/p_{2,m}$, it is usually more convenient to compute $\overline{F}({\overline{\mathit{x}}}^m(q))$ for all M candidate solutions ${\overline{\mathit{x}}}^m(q)=({\overline{x}}_1^m(q),{\overline{x}}_2^m(q))$, with(13) $${\overline{x}}_i^m(q)={\left(\sum _{k\in \mathcal{M}}{\eta }_k\min \left\{\frac{p_{j,k}}{p_{j,m}}:j\in \mathcal{I}\right\}\right)}^{-1}\frac{q}{p_{i,m}},(i,m)\in \mathcal{I}\times \mathcal{M},$$(13)

instead of using Eq. (12), and to then select the optimal envelope input as the cheapest of these, so ${\overline{\mathit{x}}}^{*}(q)={\overline{\mathit{x}}}^{m^{*}}(q)$, with(14) $$m^{*}\in \arg \underset{m\in \mathcal{M}}{\min }\left\{\mathbf{c}·{\overline{\mathit{x}}}^m(1)\right\}.$$(14)

The critical grade $m^{*}$ does thereby not depend on the magnitude of the output, since both the cost and the envelope output function are homogeneous (of degree 1).

Example 2 .

In the setting of Ex. 1, we find ${\overline{\mathit{x}}}^1(1)={\overline{\mathit{x}}}^2(1)=(1,2)$, ${\overline{\mathit{x}}}^3(1)=(10/7,10/7)$, and ${\overline{\mathit{x}}}^4(1)={\overline{\mathit{x}}}^5(1)=(2,1)$. Since $\mathbf{c}·{\overline{\mathit{x}}}^1(1)=5<\mathbf{c}·{\overline{\mathit{x}}}^3(1)=40/7<\mathbf{c}·{\overline{\mathit{x}}}^5(1)=7$, this implies the critical grades $m^{*}\in \arg {\min }_{m\in \mathcal{M}}\{c·{\overline{\mathit{x}}}^m(1)\}=\{1,2\}$, and thus the solution to the envelope optimization problem, ${\overline{\mathit{x}}}^{*}(q)=q{\overline{\mathit{x}}}^{m^{*}}(1)=(q,2q)$, for all $q>0$, as before.

3.5 Approximate Solution

By La. 3 the envelope output $\overline{F}$ exceeds the expected output F. For any target output $q>0$, output efficiency (i.e., the fact that $\overline{F}({\overline{\mathit{x}}}^{*}(q))=q$) means the solution ${\overline{\mathit{x}}}^{*}(q)$ to the envelope optimization problem ( $\ast \ast$) is not a feasible point of the original selective matching problem ( $\ast$). However, since the envelope objective is homogeneous (of degree 1), it is enough to scale the optimal envelope input to achieve feasibility, leading to an approximate solution to the selective matching problem, 3(15) $${\widehat{\mathit{x}}}^{*}(q)=\left[\frac{q}{F({\overline{\mathit{x}}}^{*}(q))}\right]{\overline{\mathit{x}}}^{*}(q)=\frac{{\overline{\mathit{x}}}^{*}(1)q^2}{F({\overline{\mathit{x}}}^{*}(q))},$$(15)

which in itself is sublinear in q (as F is convex in q due to a risk-pooling effect). Consequently, the approximate cost function,(16) $$\widehat{C}(q)={\left(\sum _{k\in \mathcal{M}}{\eta }_k\underset{i\in \mathcal{I}}{\min }\left\{\frac{p_{i,k}}{p_{i,m^{*}}}\right\}\right)}^{-1}\left(\sum _{i\in \mathcal{I}}\frac{c_i}{p_{i,m^{*}}}\right)\frac{q^2}{F({\overline{\mathit{x}}}^{*}(q))},$$(16)

is an upper bound for $C(q)$ in ( $\ast$), where the critical grade $m^{*}$ is determined by Eq. (14). On the other hand, the optimal envelope cost, $\overline{C}(q)=\mathbf{c}·{\overline{\mathit{x}}}^{*}(q)$, must constitute a lower bound for the optimal cost, so(17) $$\underset{\mathbf{c}·{\overline{\mathit{x}}}^{*}(q)}{\underbrace{\overline{C}(q)}}\le \underset{\mathbf{c}·{\mathit{x}}^{*}(q)}{\underbrace{C(q)}}\le \underset{\mathbf{c}·{\widehat{\mathit{x}}}^{*}(q)}{\underbrace{\widehat{C}(q)}},$$(17)

for all $q>0$. At the approximate solution ${\widehat{\mathit{x}}}^{*}(q)$, the relative output error becomes(18) $$\epsilon (q)=\frac{\overline{F}({\widehat{\mathit{x}}}^{*}(q))-F({\widehat{\mathit{x}}}^{*}(q))}{F({\widehat{\mathit{x}}}^{*}(q))}=\frac{q^2}{F({\overline{\mathit{x}}}^{*}(q))·F({\widehat{\mathit{x}}}^{*}(q))}-1\le \frac{q}{F({\overline{\mathit{x}}}^{*}(q))}-1,$$(18)

since, by construction, $F({\widehat{\mathit{x}}}^{*}(q))\ge q$. We note that the term on the right-hand side of the preceding inequality (18) also bounds the relative cost overage,(19) $$\frac{\widehat{C}(q)-C(q)}{C(q)}=\frac{\mathbf{c}·{\widehat{\mathit{x}}}^{*}(q)}{\mathbf{c}·{\mathit{x}}^{*}(q)}-1\le \frac{q}{F({\overline{\mathit{x}}}^{*}(q))}-1,$$(19)

where we have used Eqs. (15) and (17). Based on the preceding inequality, it is now possible to provide an a priori performance guarantee for the relative cost overage which depends only on the desired level of production output.

Lemma 6 (Relative Cost Overage).

For any target output $q>0$, the relative cost overage is such that (20) $$\frac{\widehat{C}(q)-C(q)}{C(q)}\le \left(2\sqrt{\frac{1-p_{\min }}{\pi p_{\min }}}\right)\frac{1}{\sqrt{q}},$$(20)

where $p_{\min }={\min }_{(i,m)\in \mathcal{I}\times \mathcal{M}}\left\{p_{i,m}\right\}>0$.

In particular, the relative cost overage is $O(q^{-1/2})$, in the sense that it is inversely proportional to the square root of the target output. For example, if $p_{\min }=10\%$, the upper bound for the relative cost overage (corresponding to the right-hand side in Eq. (20)) becomes $(6/\sqrt{\pi })/\sqrt{q}\approx (3.3851)/\sqrt{q}$; this conservative bound drops below 10% for target outputs of 1,146 units and beyond. The latter value corresponds to the input bound implied by La. 5; cf. footnote 2 .

Example 3 .

In the setting of Exs. 1 and 2, the solution to the envelope optimization problem, ${\overline{\mathit{x}}}^{*}(q)=(q,2q)$, produces the expected output $F({\overline{\mathit{x}}}^{*}(q))$ of about 94.6345 for $q=100$ (resp., an expected output of about 983.2032 for $q=1000$). Accordingly, the approximate solution of the selective matching problem is ${\widehat{\mathit{x}}}^{*}(100)=[100/94.6345](100,200)\approx (105.67,211.34)$ with relative output error $\epsilon (100)\approx 5.499\%\le 5.67\%$ (resp., ${\widehat{\mathit{x}}}^{*}(1000)\approx (1017.1,2034.2)$ with $\epsilon (1000)\approx 1.692\%\le 1.71\%$). By Eq. (19) the upper bound of 5.67% (resp., 1.71%) also majorizes the relative cost overage incurred by the approximate solution ${\widehat{\mathit{x}}}^{*}(q)$ to the selective matching problem compared to the optimal solution ${\mathit{x}}^{*}(q)$, for $q=100$ (resp., for $q=1000$).

3.6 Cost of Selective Matching

To gauge the cost of selective matching as a function of the grade cardinality $M=\left|\mathcal{M}\right|$, we consider a symmetric setting with maximum entropy, where all grades are attained with equal probability (i.e., ${\mathbf{p}}_i=(1/M,\dots ,1/M)$ and $\gamma =1$), and where the output weights and the unit costs are all equal (i.e., $(\eta ,\gamma )=\mathbf{1}$). Assuming a symmetric input $\mathit{x}=(N,N)$ for some $N\ge M$, La. 1 yields the expected grade-m output,$$F_m(N,N)=\frac{N}{M}-\sqrt{\frac{N}{\pi M}\left(1-\frac{1}{M}\right)},m\in \mathcal{M}.$$

Summing over all matching grades $m\in \mathcal{M}$, with $F(N,N)=M{F}_m(N,N)$) yields the (symmetric) relative output loss for selective-matching,$$\mathrm{L}(N,M)=\frac{N-F(N,N)}{N}=\sqrt{\frac{M-1}{\pi N}}\le \frac{1}{\sqrt{\pi }},$$

which is increasing (sublinearly) in $M\ge 2$ and decreasing in $N\ge M$. Since by symmetry Eq. (8) in Thm. 1 is automatically satisfied, the cost-minimizing input quantity $N^{*}(q)$ for both type-i inputs, given a desired output $q\ge M$, is entirely determined by the output-efficiency condition (7),$$N^{*}(q)=q+\frac{M-1}{2\pi }\left(1+\sqrt{1+\frac{4\pi q}{M-1}}\right),q\ge M.$$

This output-driven demand leads to a relative input overage,$$\mathrm{O}(q)=\frac{N^{*}(q)-q}{q}=\frac{M-1}{2\pi q}\left(1+\sqrt{1+\frac{4\pi q}{M-1}}\right),$$

which is decreasing in the target output quantity q, due to risk pooling across matching classes. The value of $\mathrm{O}(q)$ is also equal to the relative cost overage (i.e., $\mathrm{O}(q)=(C(q)-2cq)/(2cq)$, for any $\mathbf{c}=(c,c)$ with $c>0$) compared with perfect matching. 4 However, this relative overage increases almost linearly in the number of excess matching classes $M-1$.

Example 4 .

When transitioning from $M=20$ to $M=40$ matching classes, the relative overage in input or cost increases by a factor of at least 1.9279 (for $q\ge 1$), all else equal. For $q\to \infty$, the factor approximates $(40-1)/(20-1)\approx 2.0526$. That is, for small target outputs q we observe a slight sublinear increase of the relative overage in the number of matching classes, which becomes somewhat superlinear when target outputs are relatively large.

3.7 Multi-Type Matching

When there are more than two part types (i.e., for $\mathcal{I}=\{1,\dots ,I\}$ with $I>2$), an explicit formula for the expected output under the normal approximation is no longer available. This in turn makes it more difficult to derive solid approximation results. The corresponding I-type envelope optimization problem can still be solved explicitly. Eqs. (13) and (14) remain valid, and they determine an optimal solution ${\overline{\mathit{x}}}^{*}(q)$ of the I-type envelope optimization problem (**). The re-scaling in Eq. (15) still applies and so does Eq. (16) for the approximate cost. To illustrate the behavior of the approximation error $r({\overline{x}}^{*}(q))$ for $I=4$ compared to $I=2$, we can “double” Ex. 1 by setting ${\mathbf{p}}_3={\mathbf{p}}_1$ and ${\mathbf{p}}_4={\mathbf{p}}_2$ for the additional component types $i=3$ and $i=4$. This ensures a fair comparison when the number of component types increases by a factor of 2. Accordingly, the optimal solution to the envelope optimization problem for $I=4$ becomes ${\overline{x}}^{*}(q)=(q,2q,q,2q)$ (compared to ${\overline{x}}^{*}(q)=(q,2q)$ for $I=2$; cf. Ex. 2). Fig. 4(a) shows the behavior of the relative approximation error, $r({\overline{x}}^{*}(q))$ (defined after La. 4), as a function of $q\in [500,10000]$. The (true) expected output is obtained (approximately) as the sample-average using the corresponding normal distribution (for $N=10000$ joint samples) for each matching class, according to Eq. (9). Indeed, the law of large numbers guarantees a consistent approximation of the population mean by the sample mean.

The relative error, while increasing on average by a factor of about 2.94 on the given range (decreasing in q; see Fig. 4(b)) drops below 5% shortly after passing $q=1000$, so that for industrial-size problems the provided approximation continues to be reasonable. A more precise treatment of multi-type matching remains an interesting topic for further research.

4 Application

Consider the production of the oscillator for a mechanical watch movement where hairsprings (i.e., type-1 components) and balance wheels (i.e., type-2 components) are mated across $M\ge 2$ matching classes, where M is an even number. Modern hairsprings are often made in silicon for its antimagnetic properties and are delivered in silicon wafers, with substantial variations of their stiffness across any sample. We assume that the spring coefficient $S_1$ [in Nm/rad] is random with realizations in the interval $[a_1,b_1]$. Balance wheels are usually made from metal alloys (e.g., containing beryllium or gold among other constituents) and are characterized by the moment of inertia $S_2$ [in ${\text{kg}\mathrm{ }\mathrm{m}}^2$]. The latter can be considered a random variable with realizations in the interval $[a_2,b_2]$. The prespecified constants $a_i$ and $b_i$ are such that $0<a_i<b_i<\infty$. As already pointed out in footnote 1, the key characteristic of the assembled oscillator is the oscillation period,$$T=2\pi \sqrt{S_2/S_1},$$

which, in itself a random variable, should have realizations close to the target period $t_0$ for the particular movement under consideration. For example, $1/t_0\in \{2.5,3,3.5,4,5\}\hspace{0.17em}\text{Hz}$ corresponds to today’s most common beat rates, $\text{BR}=2(3600\mathrm{s})/t_0$ in $\{18000,21600,25200,28800,36000\}$, which are measured in half-periods (“ticks”) per hour.

4.1 Model Identification

The four sets of primitives for the selective-assembly model, (i)–(iv), are now identified in turn, first the matching classes (i), then the parameters of the multinomial distribution for the sorting products (ii), the output values (iii), and finally the unit input costs (iv).

(i) Matching Classes.

The M type-i matching classes ${\mathcal{C}}_{i,m}=[s_{i,m-1},s_{i,m})$, for $m\in \mathcal{M}\setminus \left\{M\right\}$, and ${\mathcal{C}}_{i,M}=[s_{i,M-1},s_{i,M}]$ for $m=M$ feature $M+1$ breakpoints $s_{i,0},\dots ,s_{i,M}$ such that$$a_i=s_{i,0}<s_{i,1}<\cdots <s_{i,M-1}<s_{i,M}=b_i,$$

for $i\in \mathcal{I}$. According to Weber (2021), the optimal set of breakpoints, which minimize the maximum absolute deviation of T from $t_0$, in this setup are such that the maximum matching error e is equal over all matching classes, so$$e=\underset{(s_1,s_2)\in {\mathcal{C}}_m}{\max }\left|2\pi \sqrt{s_2/s_1}-t_0\right|,m\in \mathcal{M},$$

where ${\mathcal{C}}_m={\mathcal{C}}_{1,m}\times {\mathcal{C}}_{2,m}$ denotes the m-th joint matching class. This approach requires a consistency condition, namely that the ratio of the endpoints $\beta =b_i/a_i$ is constant across the component types $i\in \mathcal{I}$, so$$\frac{e}{t_0}=\frac{{\beta }^{1/M}-1}{{\beta }^{1/M}+1}.$$

The latter can always be achieved by increasing (resp., decreasing) one of the four endpoints by extending the range of one characteristic (or by adding a zero-value matching class). The optimal set of breakpoints is then given by$$\begin{array}{c}{s}_{i,2k}={\left(\frac{t_0+e}{t_0-e}\right)}^{2k}{a}_i,s_{i,2k+1}=\{\begin{array}{cc}{s}_{2,2k}{\left[(2\pi )/(t_0-e)\right]}^2,\hfill & \text{if}\mathrm{ }i=1,\hfill \\ s_{1,2k}{\left[(t_0+e)/(2\pi )\right]}^2,\hfill & \text{if}\mathrm{ }i=2,\hfill \end{array}\hfill \end{array}$$

for $k\in \{0,\dots ,M/2\}$ and $i\in \mathcal{I}$, where we limit attention to the case where M is even.

(ii) Distributional Parameters.

Given a distribution of the random part characteristics, the multinomial probabilities $p_{i,m}$, for $(i,m)\in \mathcal{I}\times \mathcal{M}$, follow directly from Eq. (1) in Sec. 2.1. Alternatively, it is possible to avoid distributional assumptions about the underlying distribution of characteristics and proceed in a fully data-driven way, based on repeated observation of realizations ${\mathit{x}}_i^{(k)}=(x_{i,1}^{(k)},\dots ,x_{i,M}^{(k)})$ for the random sorting product ${\mathbf{X}}_i$. This yields the maximum-likelihood estimates ${\widehat{\mathbf{p}}}_i=({\widehat{p}}_{i,1},\dots ,{\widehat{p}}_{i,M})$ of the unknown parameters ${\mathbf{p}}_i=(p_{i,1},\dots ,p_{i,M})$ of the multinomial distribution,(21) $${\widehat{p}}_{i,m}=\frac{\sum _{k=1}^K{x}_{i,m}^{(k)}}{\sum _{m\in \mathcal{M}}\sum _{k=1}^K{x}_{i,m}^{(k)}},$$(21)

where $k\in \{1,\dots ,K\}$ is the batch number, $x_{i,m}^{(k)}$ denotes the number of grade-m components in batch k, and K the number of observed batches. That is, the maximum-likelihood estimates for the component likelihoods in their respective type-grade combinations correspond simply to their respective frequencies over all observed samples.

Remark 5 (Off-Spec Components).

Because of the empirical regularity that a received batch of input parts may well contain items that do not fit into any of the matching classes, the following question arises: What happens when relaxing the normalization in Eq. (2) to account for “off-spec” components? Specifically, if $x_i^{(k)}$ denotes the size of batch k, then the expression for ${\widehat{p}}_{i,m}$ in Eq. (21) implicitly allows for off-spec components (which do not fit in any ${\mathcal{C}}_{i,m}$), since generally $\sum _{m\in \mathcal{M}}{\sum }_{k=1}^K{x}_{i,m}^{(k)}\le x_i^{(k)}$, while at the same time ${\widehat{p}}_{i,1}+\cdots +{\widehat{p}}_{i,M}=1$ satisfying Eq. (2). Yet, to explicitly account for the possibility of “off-spec” components let$${\widehat{p}}_{i,0}=1-\frac{\sum _{m\in \mathcal{M}}\sum _{k=1}^K{x}_{i,m}^{(k)}}{\sum _{k=1}^K{x}_i^{(k)}},i\in \mathcal{I},$$

denote the observed frequency of type-i reject components, as the maximum-likelihood estimate of the true underlying reject probability $p_{0,i}$. Then by adding a matching class 0 with output weight ${\eta }_0$, by Eq. (13) one can simply rescale the approximate solution to the selective matching problem in Eq. (15) to account for the reject probability:(22) $${\widehat{x}}_i^{*}(q)=\frac{1}{1-p_{i,0}}\left[\frac{q}{F({\overline{\mathit{x}}}^{*}(q))}\right]{\overline{x}}_i^{*}(q),i\in \mathcal{I},$$(22)

where $q>0$ is any given output target and ${\overline{x}}_i^{*}(q)$ is given in Thm. 2 (for $p_{i,m}\equiv {\widehat{p}}_{i,m}$). Hence, we can decouple the consideration of off-spec components from the main analysis (which is purely “on-spec”) and then rescale the approximate solution ex post.

(iii) Net Unit Values.

The weighted total output in Eq. (4) requires a vector of positive weights $\mathit{\eta }=({\eta }_1,\dots ,{\eta }_M)$. Since the quality of the assembled mechanical oscillators does not depend on the index m of the matching grade, it is most natural for the firm to set all ${\eta }_m=1$ and consider standard output, i.e., $Q=Y_1+\cdots +Y_M$. In standard selective-manufacturing applications, no grade is privileged over others (except possibly for a designated “reject” grade of negligible net unit value); in that case, the grade sorting satisfies the sole purpose of sharpening the precision of the assembly output given the quality imprecisions of the component inputs.

(iv) Unit Input Costs.

The cost function in the selective matching problem ( $\ast$) depends on the vector $\mathbf{c}=(c_1,c_2)$ of the unit costs for the two component types. However, as evident from Thm. 1, for any given output quantity $q>0$ the optimal input vector ${\mathit{x}}^{*}(q)$ depends only on the ratio $\gamma =c_1/c_2$ of the unit costs, so that there is no need to know their absolute values to determine the composition of the firm’s optimal input portfolio. We assume that the unit cost of a balance wheel is about 10 to 20 times the unit cost of a silicon hairspring, so $\gamma \in [0.05,0.1]$. In the numerical example, we set $\mathbf{c}=(1,10)$, so $\gamma =0.1$; the results for $\gamma =0.05$ are qualitatively very similar.

4.2 Data Generation

To generate synthetic samples representing the observed stochastic supply characteristics, assume that the true distribution of the normalized type-i component characteristic,$$Z_i=\frac{S_i-{\mu }_{S_i}}{{\sigma }_{S_i}},$$

follows a standard normal distribution with cdf $\Phi (·)$, where ${\mu }_{S_i}\in (a_i,b_i)$ is the mean and ${\sigma }_{S_i}>0$ the standard deviation of the (stationary) component population. Given the distributional assumption, the multinomial probabilities (conditional on parts being on-spec) can be obtained by adjusting Eq. (1) as follows:$$p_{i,m}=\frac{1}{1-p_{i,0}}\left(\Phi \left(\frac{s_{i,m}-{\mu }_{S_i}}{{\sigma }_{S_i}}\right)-\Phi \left(\frac{s_{i,m-1}-{\mu }_{S_i}}{{\sigma }_{S_i}}\right)\right),(i,m)\in \mathcal{I}\times \mathcal{M}.$$

Note that with positive probability some realizations $s_i$ of the random component characteristic $S_i$ may be “off-spec,” in which case we use input scaling as noted in Sec. 4.1 (ii). The off-spec probability in the component population is $p_{i,0}=1-(\Phi ((b_i-{\mu }_{S_i})/{\sigma }_{S_i})-\Phi ((a_i-{\mu }_{S_i})/{\sigma }_{S_i}))$. As explained in Remark 5 on off-spec components, the renormalization of the multinomial probabilities has no effect on the likelihood ratios in Eq. (13), so that it is enough to rescale the approximate solution as in Eq. (22), an approach which is implemented here. Finally, we assume that both distributions are centered, so ${\mu }_{S_i}=(a_i+b_i)/2$ for $i\in \mathcal{I}$, but that the coefficient of variation is twice as large for hairsprings than for balance wheels, with ${\sigma }_{S_1}/{\mu }_{S_1}=2\%$ and ${\sigma }_{S_2}/{\mu }_{S_2}=1\%$.

4.3 Performance Analysis

Let $(a_1,b_1)=(2.94,3.06)\times {10}^{-7}$ [in Nm/rad] be the boundaries of the hairspring stiffness and $(a_2,b_2)=(4.655,4.845)\times {10}^{-10}$ [in $\text{kg}\hspace{0.17em}{\mathrm{m}}^2$] the boundaries of the balance-wheel inertia (in either case corresponding to a spread of 2% around a central value), so $\beta =a_2/a_1=b_2/b_1$, as required in Sec. 4.1 (i) where we also provide the optimal set of breakpoints for the movement frequency of $1/t_0=4\hspace{0.17em}\text{Hz}$, and thus the optimal matching classes ${\mathcal{C}}_{i,m}$. Clearly, the maximum matching error e decreases in the number of matching classes $M\ge 2$, which is a choice variable. For practical purposes it is useful to interpret $e/t_0$ in terms of “deviation per day” (dpd) [in seconds/day, i.e., s/d] by multiplying numerator and denominator by 345,600 ( $=24\times 3,600\times 4$). Thus, since $4\times t_0=1\mathrm{s}$, we obtain that $e/t_0=\text{dpd}=(345,600\times e)/(1\mathrm{d})$, so(23) $$M(\text{dpd})=\lceil {\left(\hspace{0.17em}\log \hspace{0.17em}\left(\frac{1+\text{dpd}}{1-\text{dpd}}\right)\right)}^{-1}\hspace{0.17em}\log \hspace{0.17em}(\beta )\rceil ,$$(23)

where $\beta \approx 1.040816327$. That is, Eq. (23) provides the number of matching classes necessary to guarantee a given dpd, as shown in Fig. 5. For example, it is $M(60)=29$, $M(30)=58$, and $M(15)=116$. Finally, for a rather fine partition of part characteristics into $M(5)=346$ matching classes it would be possible (at least theoretically) to achieve the COSC standard of $\pm$5s/d without any watchmaker intervention after assembly. Indeed, the Contrôle Officiel Suisse des Chronomètres (COSC) requires mechanical watch movements to run within –4s and +6s per day to issue a “chronometer”' certification, which is equivalent to a dpd of 5. For practical purposes, and to stay exactly in the framework of Weber 2021 used in Sec. 4.1 (i), we round up to the smallest even number of matching classes, thus considering $2\lceil M(\text{dpd})/2\rceil$ instead of $M(\text{dpd})$ in Eq. (23), as indicated in Fig. 5.

To guarantee practical implementability and feasibility, the approximate solution ${\widehat{\mathit{x}}}^{*}(q)$ in Eq. (15) is rounded to an input vector $\lceil {\widehat{\mathit{x}}}^{*}(q)\rceil$, with the next-highest integers as components (i.e., $\lceil {\widehat{x}}_i^{*}(q)\rceil$, for $i\in \mathcal{I}$). An integer discretization (determined by optimality) is also applied to the solution ${\mathit{x}}^{*}(q)$ of the original selective matching problem ( $\ast$). This solution discretization accounts for some of the fluctuations in Fig. 6, which depicts the optimal costs (in terms of C and $\widehat{C}$; cf. Eq. (17)) and optimal inputs. To ensure comparability, the costs are normalized to the input cost $C_0(q)=(c_1+c_2)q$ without matching classes (where $M=1$), and the input solutions are normalized to the “overhead-free” input q for all types $i\in \mathcal{I}$ in that case. As shown in Fig. 6(a), the normalized cost decreases in the target quantity q, which reflects a risk-pooling effect implied by the sublinear increase of the standard deviation of the part characteristics in the amount ordered. It is also apparent that the deviation between the optimal cost C (for optimal integer inputs determined via grid search) and the cost approximation $\widehat{C}$ is quite negligible for realistic output targets. It is also apparent that as the dpd-precision requirement becomes more demanding (i.e., smaller), the normalized cost increases very fast—well exceeding 4 for COSC chronometer specifications. Fig. 6(b) illustrates (for $\gamma =0.1$) that these cost increases are driven by a massive over-ordering of the 10-times cheaper hairsprings as “buffer components.” The over-ordering increases both in the production target and the precision requirement.

5 Conclusion

The envelope approximation of the firm’s output objective proved effective in determining a simple closed-form solution for the cost-optimal input portfolio in Eqs. (13) and (22), which includes the possibility of scrap parts. In that solution, somewhat intentionally, the notation of the index set $\mathcal{I}$ for the different part types (e.g., hairsprings and balance wheels in the watchmaking application) has been used for displaying the results. Indeed, an extension to $\mathcal{I}=\{1,\dots ,I\}$ with $I>2$ part types (e.g., balls, inner/outer rings, and cages, in the case of a ball-bearing assembly) is quite straightforward, and requires no adjustment of these formulas. While this naturally applies to the envelope output in Eq. (9), the same cannot be said for the expected output in Eq. (6), for which a simple closed-form expression may not be within easy reach for more than two input part types.

From a practical viewpoint, a doubling of the number of matching classes tends to roughly double the precision of output parts. 5 As noted in Sec. 3.6, it also leads to almost a doubling of the relative input overage, which is the number of extra units required compared with perfect matching, where the latter corresponds to random matching with equal inputs without quality concerns—equivalent to using a single matching class. That is, all else equal, improving the output tolerance by a certain factor basically multiplies the overage cost by that factor. In addition, one should not overlook the fact that unused parts appear in matching classes for which the counterpart of the other type is stocked-out and on backorder. Yet, the residual inventory vector serves as a buffer for the next order and thus increases expected output relative to an empty stock. A formulation of the corresponding dynamic reordering policy, taking into account parts residuals, is an interesting topic for future research.

Appendix: Proofs

Proof of Lemma 1. Let $m\in \mathcal{M}$ be a given matching grade and let $\mathit{x}\in {ℝ}_{+}^2$ be an arbitrary input vector. By Eq. (3) the firm’s random grade-m output (conditional on $\mathit{x}$),$$Y_m=\min \{X_{1,m},X_{2,m}\}=X_{1,m}+\min \{0,X_{2,m}-X_{1,m}\}=X_{1,m}+\min \{0,{\Delta }_m\},$$

can be viewed as a sum of the grade-m type-1 output and the negative part of the grade-m parts defect, ${\Delta }_m=X_{2,m}-X_{1,m}$. As the difference of two correlated and normally distributed random variables, the latter follows a normal distribution with mean ${\delta }_m={\mu }_{2,m}-{\mu }_{1,m}$, and variance ${\sigma }_m^2={\sigma }_{1,m}^2+{\sigma }_{2,m}^2$. Thus, the expected grade-m output is $F_m(\mathit{x})=\mathbb{E}\left[Y_m|\mathit{x}\right]={\mu }_{1,m}+\mathbb{E}\left[\min \{0,{\Delta }_m\}|\mathit{x}\right]$, where$$\mathbb{E}\left[\min \{0,{\Delta }_m\}|\mathit{x}\right]={\delta }_m+{\sigma }_m\mathbb{E}\left[\min \left\{-\frac{{\delta }_m}{{\sigma }_m},\frac{{\Delta }_m-{\delta }_m}{{\sigma }_m}\right\}|\mathit{x}\right].$$

On the other hand, with $({\Delta }_m-{\delta }_m)/{\sigma }_m$ following a standard normal distribution,$$\mathbb{E}\left[\min \left\{-\frac{{\delta }_m}{{\sigma }_m},\frac{{\Delta }_m-{\delta }_m}{{\sigma }_m}\right\}|\mathit{x}\right]=-\frac{{\delta }_m}{{\sigma }_m}\left(1-\Phi \left(-\frac{{\delta }_m}{{\sigma }_m}\right)\right)+{\int }_{-\infty }^{-{\delta }_m/{\sigma }_m}\xi d\Phi (\xi ).$$

Taking into account the symmetry of $\varphi (·)$ the last term can be computed explicitly,$${\int }_{-\infty }^{-\delta /{\sigma }_m}\xi d\Phi (\xi )={\left[-\frac{\hspace{0.17em}\exp \hspace{0.17em}(-{\xi }^2/2)}{\sqrt{2\pi }}\right]}_{\xi \to -\infty }^{\xi =-{\delta }_m/{\sigma }_m}=-\varphi \left(-\frac{{\delta }_m}{{\sigma }_m}\right)=-\varphi \left(\frac{{\delta }_m}{{\sigma }_m}\right).$$

Hence, the expected grade-m output becomes$$F_m(\mathit{x})={\mu }_{1,m}+{\delta }_m-{\delta }_m\left(1-\Phi \left(-\frac{{\delta }_m}{{\sigma }_m}\right)\right)-{\sigma }_m\varphi \left(\frac{{\delta }_m}{{\sigma }_m}\right).$$

By virtue of the fact that $\Phi (-\xi )=1-\Phi (\xi )$, for all $\xi \in ℝ$, the preceding expression simplifies,$$F_m(\mathit{x})={\Phi }_m{\mu }_{1,m}+(1-{\Phi }_m){\mu }_{2,m}-{\varphi }_m{\sigma }_m,$$

where we have set ${\Phi }_m=\Phi ({\delta }_m/{\sigma }_m)$ and ${\varphi }_m=\varphi ({\delta }_m/{\sigma }_m)$. This establishes the claim. $\square$

Proof of Lemma 2. Let $m\in \mathcal{M}$ be a given matching grade and let $\mathit{x}=(x_1,x_2)\in {ℝ}_{+}^2$ be an arbitrary input vector. For any $i\in \mathcal{I}$, the number of type-i components in grade-m follows a normal distribution with cdf $G_{i,m}(·|x_i)$. The cdf of the random grade-m output $Y_m=\min \{X_{1,m},X_{2,m}\}$ becomes$$\mathrm{\mathbb{P}}(Y_m\le \xi |\theta )=G_{1,m}(\xi |x_1)+G_{2,m}(\xi |x_2)-G_{1,m}(\xi |x_1)G_{2,m}(\xi |x_2),\xi \ge 0.$$

Consider now an increase of type-i input from $x_i$ to ${\widehat{x}}_i>x_i$, and note that$$\frac{\partial G_{i,m}(\xi |x_i)}{\partial x_i}=-\frac{p_{i,m}+(\xi /x_i)}{2{\sigma }_{i,m}}\varphi \left(\frac{\xi -{\mu }_{i,m}}{{\sigma }_{i,m}}\right)<0,x_i\ge 0.$$

This implies first-order stochastic dominance (FOSD) of the type-i quantity distribution after the quantity increase, $G_{i,m}(\xi |{\widehat{x}}_i)<G_{i,m}(\xi |x_i),\xi \ge 0$. Thus, if we denote by $({\widehat{X}}_{1,m},{\widehat{X}}_{2,m})$ the grade-m yields under $\widehat{\mathit{x}}$, then $(X_{1,m},X_{2,m}){\prec }_{\text{FOSD}}({\widehat{X}}_{1,m},{\widehat{X}}_{2,m})$, where ${\prec }_{\text{FOSD}}$ marks the (strict) FOSD-order. But this implies that $Y_m=\min \{X_{1,m},X_{2,m}\}{\prec }_{\text{FOSD}}\min \{{\widehat{X}}_{1,m},{\widehat{X}}_{2,m}\}={\widehat{Y}}_m$, so that we can now compare the expected values, $F_m(\mathit{x})=\mathbb{E}[Y_m|\mathit{x}]<\mathbb{E}[{\widehat{Y}}_m|\widehat{\mathit{x}}]=F_m(\widehat{\mathit{x}})$, establishing the claimed monotonicity. $\square$

Proof of Theorem 1. Let $(q,c)\in {ℝ}_{++}\times {ℝ}_{++}^2$. The Lagrangian associated with the cost-minimization problem ( $\ast$) is $L(\mathit{x};\lambda )=c_1{x}_1+c_2{x}_2-\lambda (F(\mathit{x})-q)$, where the adjoint variable $\lambda \ge 0$ measures the “shadow cost” of relaxing the output-attainment constraint. The first-order necessary optimality conditions are therefore(24) $$\frac{\partial L(\mathit{x};\lambda )}{\partial x_i}=c_i-\lambda \frac{\partial F(\mathit{x})}{\partial x_i}=0,$$(24)

for all $i\in \mathcal{I}$ (Zorich, 2004, p. 527ff). Since by La. 2 it is $\partial F(\mathit{x})/\partial x_i>0$, which in particular applies to input vectors $\mathit{x}$ such that $F(\mathit{x})=q$, the output-attainment constraint $F(\mathit{x})\ge q$ must be binding, so that (by complementary slackness) necessarily $\lambda >0$ and Eq. (7) must be satisfied if $\mathit{x}={\mathit{x}}^{*}(q)$ is to be a solution candidate for the selective matching problem ( $\ast$). Substitution of Eqs. (5) and (6) (in La. 1) into Eq. (24) for $i\in \{1,2\}$ yields the marginal rate of technical substitution,$$\frac{\partial F(\mathit{x})/\partial x_1}{\partial F(\mathit{x})/\partial x_2}=\frac{\sum _{m\in \mathcal{M}}{\eta }_m\partial F_m(\mathit{x})/\partial x_1}{\sum _{m\in \mathcal{M}}{\eta }_m\partial F_m(\mathit{x})/\partial x_2}=\frac{c_1}{c_2}=\gamma ,$$

where$$\begin{array}{cc}\frac{\partial F_m(\mathit{x})}{\partial x_1}\hfill & =p_{1,m}{\Phi }_m-{\delta }_m\frac{\partial {\Phi }_m}{\partial x_1}-{\sigma }_m\frac{\partial {\varphi }_m}{\partial x_1}-\frac{p_{1,m}(1-p_{1,m})}{2{\sigma }_m}{\varphi }_m,\hfill \\ \frac{\partial F_m(\mathit{x})}{\partial x_2}\hfill & =p_{2,m}(1-{\Phi }_m)-{\delta }_m\frac{\partial {\Phi }_m}{\partial x_2}-{\sigma }_m\frac{\partial {\varphi }_m}{\partial x_2}-\frac{p_{2,m}(1-p_{2,m})}{2{\sigma }_m}{\varphi }_m,\hfill \end{array}$$

and$$-{\delta }_m\frac{\partial {\Phi }_m}{\partial x_i}={\sigma }_m\frac{\partial {\varphi }_m}{\partial x_i}=\frac{{\delta }_m}{{\sigma }_m}\left({(-1)}^{i+1}{p}_{i,m}+\frac{{\delta }_m}{{\sigma }_m}\frac{p_{i,m}(1-p_{i,m})-2C_{i,m}}{2{\sigma }_m}\right){\varphi }_m,i\in \mathcal{I}.$$

As a result,$$\frac{\partial F_m(\mathit{x})}{\partial x_1}=p_{1,m}\left({\Phi }_m-\frac{1-p_{1,m}}{2{\sigma }_m}{\varphi }_m\right)\text{and}\hspace{0.17em}\frac{\partial F_m(\mathit{x})}{\partial x_2}=p_{2,m}\left(1-{\Phi }_m-\frac{1-p_{2,m}}{2{\sigma }_m}{\varphi }_m\right),$$

which yields Eq. (8) for $\mathit{x}={\mathit{x}}^{*}(q)$. This concludes our proof. $\square$

Proof of Lemma 3. Fix any input bundle $\mathit{x}=(x_1,x_2)\in {ℝ}_{+}^2$. By Jensen’s inequality, it is$$\mathbb{E}\left[Y_m\right|\mathit{x}]=\mathbb{E}\left[\min \left\{X_{1,m},X_{2,m}\right\}|\mathit{x}\right]\le \min \left\{\mathbb{E}\left[X_{1,m}\right|x_1],\mathbb{E}[X_{2,m}\left|x_2\right]\right\},m\in \mathcal{M}.$$

Hence, using Eq. (5) and the linearity of the expectation operator, one obtains$$F(\mathit{x})=\mathbb{E}\left[Q\right|\mathit{x}]=\sum _{m\in \mathcal{M}}{\eta }_m\mathbb{E}[Y_m\left|\mathit{x}\right]\le \sum _{m\in \mathcal{M}}{\eta }_m\min \left\{p_{1,m}{x}_1,p_{2,m}{x}_2\right\}=\overline{F}(\mathit{x}),$$

which concludes our proof. $\square$

Proof of Lemma 4. Let $\mathit{x}=(x_1,x_2)\in {ℝ}_{++}^2$. Consider first the case where ${\delta }_m(\mathit{x})=p_{2,m}{x}_2-p_{1,m}{x}_1\ge 0$, for some $m\in \mathcal{M}$. Then by La. 1 it is$$0\le \frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}-F_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\le \varphi \left(\frac{{\delta }_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\right)-\frac{{\delta }_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\left(1-\Phi \left(\frac{{\delta }_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\right)\right).$$

On the other hand, since $(d/d\xi )\left(\varphi (\xi )-\xi \left(1-\Phi (\xi )\right)\right)=-(1-\Phi (\xi ))<0$, for all $\xi \in ℝ$, the expression on the right-hand side of the preceding equality is downward-sloping in ${\delta }_m/{\sigma }_m>0$, so(25) $${\delta }_m(\mathit{x})\ge 0\Rightarrow 0\le \frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}-F_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\le \frac{1}{\sqrt{2\pi }}=\varphi (0).$$(25)

A similar reasoning applies for ${\delta }_m(\mathit{x})\le 0$, in which case La. 1 yields$$0\le \frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}-F_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\le \varphi \left(\frac{{\delta }_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\right)+\frac{{\delta }_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\Phi \left(\frac{{\delta }_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\right).$$

Thus, taking account of the fact that $(d/d\xi )\left(\varphi (\xi )+\xi \Phi (\xi )\right)=\Phi (\xi )>0$, for all $\xi \in ℝ$, the expression on the right-hand side of the preceding inequality must be upward-sloping in ${\delta }_m/{\sigma }_m<0$, whence(26) $${\delta }_m(\mathit{x})\le 0\Rightarrow 0\le \frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}-F_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\le \frac{1}{\sqrt{2\pi }}=\varphi (0).$$(26)

Combining Eqs. (25) and (26), and repeating this exercise for all grades, yields that in fact(27) $$0\le \frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}-F_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\le \frac{1}{\sqrt{2\pi }},m\in \mathcal{M}.$$(27)

Using the definitions of the weighted output F and the (equally weighted) envelope output $\overline{F}$, one obtains by mere summation that $0\le \overline{F}(\mathit{x})-F(\mathit{x})\le \sigma (\mathit{x})/\sqrt{2\pi }$, where $\sigma (\mathit{x})=\sum _{m\in \mathcal{M}}{\sigma }_m(\mathit{x})$. The (maximum) approximation error is therefore at most linear in the aggregate standard deviation $\sigma (\mathit{x})$, so $\frac{R(\mathit{x})}{\sigma (\mathit{x})}\le 1/\sqrt{2\pi }$, which concludes our proof. $\square$

Proof of Lemma 5. To exclude trivialities, we restrict attention to positive inputs, so $\mathit{x}=(x_1,x_2)\in {ℝ}_{++}^2$. By Eq. (27) we have$$0\le \frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}-F_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\le \frac{1}{\sqrt{2\pi }},m\in \mathcal{M},$$

so that$$\frac{F_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\ge \frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}}{{\sigma }_m(\mathit{x})}-\frac{1}{\sqrt{2\pi }}\ge \sqrt{\frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}}{2-p_{1,m}-p_{2,m}}}-\frac{1}{\sqrt{2\pi }},m\in \mathcal{M}.$$

On the other hand, if we set $p_{\min }={\min }_{(i,m)\in \mathcal{I}\times \mathcal{M}}\left\{p_{i,m}\right\}>0$, then$$\sqrt{\frac{\min \{p_{1,m}{x}_1,p_{2,m}{x}_2\}}{2-p_{1,m}-p_{2,m}}}\ge \sqrt{\frac{p_{\min }/2}{1-p_{\min }}}\sqrt{\min \{x_1,x_2\}},m\in \mathcal{M},$$

where the right-hand side of the preceding inequality does not depend on the grade m. Hence, it is easy to verify that$$\min \{x_1,x_2\}\ge \frac{4}{\pi }\left(\frac{1-p_{\min }}{p_{\min }}\right)\Rightarrow \frac{F_m(\mathit{x})}{{\sigma }_m(\mathit{x})}\ge \sqrt{\frac{p_{\min }/8}{1-p_{\min }}}\sqrt{\min \{x_1,x_2\}},m\in \mathcal{M}.$$

The latter implies:$$\min \{x_1,x_2\}\ge \frac{4}{\pi }\left(\frac{1-p_{\min }}{p_{\min }}\right)\Rightarrow \frac{F(\mathit{x})}{\sigma (\mathit{x})}\ge \sqrt{\frac{p_{\min }/8}{1-p_{\min }}}\sqrt{\min \{x_1,x_2\}}.$$

Combining this with La. 4, the relative error decreases in the smallest input,$$r(\mathit{x})=\frac{R(\mathit{x})/\sigma (\mathit{x})}{F(\mathit{x})/\sigma (\mathit{x})}\le \sqrt{\frac{1-p_{\min }}{p_{\min }}}\frac{2/\sqrt{\pi }}{\sqrt{\min \{x_1,x_2\}}},$$

as long as the scale of the application is sufficiently large, so$$\min \{x_1,x_2\}\ge \frac{4}{\pi }\left(\frac{1-p_{\min }}{p_{\min }}\right).$$

Thus, for any $\epsilon \in (0,1)$:$$\min \{x_1,x_2\}\ge \underset{¯}{x}(\epsilon )=\frac{4}{\pi }\left(\frac{1-p_{\min }}{p_{\min }}\right)\max \left\{1,\frac{1}{{\epsilon }^2}\right\}=\frac{4}{\pi {\epsilon }^2}\left(\frac{1-p_{\min }}{p_{\min }}\right)\Rightarrow r(\mathit{x})\le \epsilon ,$$

as claimed. But this means that$$\min \{x_1,x_2\}\ge \underset{¯}{x}(\epsilon )\Rightarrow \frac{\overline{F}(x_1,x_2)}{1+\epsilon }\le F(x_1,x_2),$$

which in turn implies the minimal output objective,$$\underset{¯}{q}(\epsilon )=\frac{\overline{F}(\underset{¯}{x}(\epsilon ),\underset{¯}{x}(\epsilon ))}{1+\epsilon }=\left(\sum _{\ell =1}^L{\widehat{\eta }}_{\ell }\min \{{\rho }_{\ell },1\}\right)\frac{\underset{¯}{x}(\epsilon )}{1+\epsilon }.$$

Finally, we note that $\underset{¯}{x}(\epsilon )$ is $O({\epsilon }^{-2})$, so that for small $\epsilon >0$ the minimum envelope input $\underset{¯}{q}(\epsilon )$ is also $O({\epsilon }^{-2})$, as claimed. This concludes the proof. $\square$

Proof of Theorem 2. As a positive linear combination of concave functions, the envelope output,$$\overline{F}(\mathit{x})=\sum _{\ell =1}^L{\widehat{\eta }}_{\ell }\min \{{\rho }_{\ell }{x}_1,x_2\},\mathit{x}=(x_1,x_2)\in {ℝ}_{+}^2,$$

is also a concave function. Thus, for any $q>0$ the upper contour set, $\mathcal{U}(q)=\left\{\mathit{x}\in {ℝ}_{+}^2:\overline{F}(\mathit{x})\ge q\right\}$, is convex. In addition, it is also possible to achieve at least as much envelope output with more inputs, so$$\mathit{x}\in \mathcal{U}(q)\Rightarrow \mathit{x}+{ℝ}_{+}^2\subset \mathcal{U}(q).$$

As a result, the upper contour set $\mathcal{U}(q)$ is a polygonal chain (i.e., a connected series of line segments),$$\mathcal{P}(q)=\left\{(x_1,x_2)\in {ℝ}_{+}^2:x_1\in [q{\overline{x}}_1^L,q{\overline{x}}_1^1],x_2=q\underset{\ell \in \{1,\dots ,L\}}{\max }\left\{{\overline{x}}_2^{\ell }-s_{\ell }\left(x_1-{\overline{x}}_1^{\ell }\right)\right\}\right\},$$

in the sense that $\mathcal{U}(q)=\mathcal{P}(q)+{ℝ}_{+}^2$, where the vertices ${\overline{\mathit{x}}}^{\ell }=(x_1^{\ell },x_2^{\ell })$ of the polygonal chain $\mathcal{P}(1)$ are implied by the conditions(28) $${\rho }_{\ell }{\overline{x}}_1^{\ell }={\overline{x}}_2^{\ell }\text{and}\overline{F}({\overline{\mathit{x}}}^{\ell })=1,$$(28)

for $\ell \in \{1,\dots ,L\}$; the edges of $\mathcal{P}(1)$ have the slopes of absolute values $s_1=0$, and $s_{\ell }=({\overline{x}}_2^{\ell }-{\overline{x}}_1^{\ell -1})/({\overline{x}}_1^{\ell -1}-{\overline{x}}_1^{\ell })$, for all $\ell \in \{2,\dots ,L\}$. The vertex conditions in Eq. (28) imply$$\sum _{k=1}^{\ell }{\widehat{\eta }}_k{\rho }_k{\overline{x}}_1^{\ell }+\sum _{k=\ell +1}^L{\widehat{\eta }}_k{\overline{x}}_2^{\ell }={\overline{x}}_1^{\ell }\left(\sum _{k=1}^L{\widehat{\eta }}_k\min \{{\rho }_k,{\rho }_{\ell }\}\right)=1,\ell \in \{1,\dots ,L\},$$

taking account of the fact that the likelihood ratios ${\rho }_{\ell }$ are size-ordered (i.e., ${\rho }_1<\cdots <{\rho }_L$) by construction. Hence, the coordinates of the vertices of $\mathcal{P}(1)$ are fully described by$${\overline{x}}_1^{\ell }={\left(\sum _{k=1}^L{\widehat{\eta }}_k\min \{{\rho }_k,{\rho }_{\ell }\}\right)}^{-1}\text{and}{\overline{x}}_2^{\ell }={\rho }_{\ell }{\overline{x}}_1^{\ell },$$

for all $\ell \in \{1,\dots ,L\}$. By the Hahn-Banach theorem (Berge, 1963, p. 157) there exists a hyperplane separating the convex set $\{\mathit{x}=(x_1,x_2)\in {ℝ}_{+}^2:c_1{x}_1+c_2{x}_2<C\}$, for some positive constant C, from the convex set $\mathcal{U}(q)=q\mathcal{U}(1)$ as long as they are disjoint, i.e., as long as$$c_1{x}_1+c_2{x}_2<C\Rightarrow \overline{F}(x_1,x_2)<q.$$

In particular, if $c_1{\overline{x}}_1^{\ell }+c_2{\overline{x}}_2^{\ell }=\overline{C}(1)$, then necessarily $\overline{F}({\overline{\mathit{x}}}^{\ell })=1$. The aforementioned separation for $C=\overline{C}(1)$ occurs at the vertex ${\overline{\mathit{x}}}^{{\ell }^{*}}$, with$${\ell }^{*}=\max \left\{\ell \in \{1,\dots ,L\}:\gamma =s_{\ell }<\frac{c_2}{c_1}\right\},$$

since that avoids an intersection of the separating hyperplane with the budget set $\{\mathit{x}=(x_1,x_2)\in {ℝ}_{+}^2:c_1{x}_1+c_2{x}_2<\overline{C}(1)\}$. Using the fact that the cost objective and the envelope output are each homogeneous of degree 1, we obtain that the input vector ${\overline{\mathit{x}}}^{*}=q{\overline{\mathit{x}}}^{{\ell }^{*}}$ solves the envelope optimization problem ( $\ast \ast$). $\square$

Proof of Lemma 6. Consider any fixed target output $q>0$ and let $p_{\min }={\min }_{(i,m)\in \mathcal{I}\times \mathcal{M}}\left\{p_{i,m}\right\}$ the minimal matching probability (which is positive). By La. 5, the right-hand side of Eq. (19) is such that$$q\ge \underset{¯}{q}(\epsilon )=\frac{4}{\pi {\epsilon }^2}\left(\frac{1-p_{\min }}{p_{\min }}\right)\Rightarrow \frac{q}{F({\overline{\mathit{x}}}^{*}(q))}-1\le \frac{q(1+\epsilon )}{\overline{F}({\overline{\mathit{x}}}^{*}(q))}-1=\epsilon ,$$

where we have used the fact that $\overline{F}({\overline{\mathit{x}}}^{*}(q))=q$, as pointed out in footnote 3. Thus, by eliminating $\epsilon$, the relative cost overage is such that$$\frac{\widehat{C}(q)-C(q)}{C(q)}\le \left(\sqrt{\frac{1-p_{\min }}{(\pi /4)p_{\min }}}\right)\frac{1}{\sqrt{q}},$$

which corresponds to Eq. (20). This concludes our proof. $\square$

Acknowledgements

The author would like to thank several anonymous referees, as well as participants of the 2022 CORS/INFORMS International Conference in Vancouver, Canada, and the 2022 EURO Conference in Espoo, Finland, for helpful suggestions.

Notes on Contributor

Thomas A. Weber is full professor of Operations, Economics and Strategy at the Swiss Federal Institute of Technology in Lausanne (EPFL). Earlier he was a faculty member in the Economics and Finance Group of the Department of Management Science and Engineering at Stanford University. Prof. Weber is an Ingénieur des Arts et Manufactures (Ecole Centrale Paris) and a Diplom-Ingenieur in Electrical Engineering (Technical University Aachen). He holds master’s degrees in Technology and Policy and Electrical Engineering and Computer Science from MIT, and a PhD in Applied Economics and Managerial Science from the Wharton School of the University of Pennsylvania. He was a visiting faculty in Economics at Cambridge University and in Mathematics at Moscow State University. Between 1998 and 2002, he was a senior consultant with the Boston Consulting Group. His current research interests include the economics of information and uncertainty, robust optimization, and strategy. Prof. Weber’s more than 100 scholarly publications have appeared, for example, in American Economic Journal: Microeconomics, Decision Support Systems, Economics Letters, Economic Theory, Health Care Management Science, Information Systems Research, Journal of Economic Dynamics and Control, Journal of Environmental Economics and Management, Journal of Management Information Systems, Journal of Mathematical Economics, Journal of Optimization Theory and Applications, Journal of Regulatory Economics, Management Science, Mathematical Reviews, Medical Decision Making, Operations Research, Operations Research Letters, Optimal Control: Applications and Methods, Physical Review E, Sustainable Production and Consumption, as well as Theory and Decision. He has been associate editor for Management Science and is the author of Optimal Control Theory with Applications in Economics (MIT Press, 2011).

Notes

¹ The hairspring features a stiffness (described by $\mathrm{\varkappa }$ [in Nm/rad]) while the balance wheel is characterized by a moment of inertia (given by I [in kg ${\mathrm{m}}^2$]). Their ratio determines the oscillation period, $T=2\pi \sqrt{I/\mathrm{\varkappa }}$.

2 As established in the proof of La. 5, a (very conservative) minimal envelope output is obtained by setting $\underset{¯}{q}(\epsilon )=\left({\sum }_{\ell =1}^L{\widehat{\eta }}_{\ell }\min \{{\rho }_{\ell },1\}\right)\underset{¯}{x}(\epsilon )/(1+\epsilon )$, where $\underset{¯}{x}(\epsilon )=(4/\pi )(1/{\epsilon }^2)(1-p_{\min })/p_{\min }$ is the corresponding minimal input, with $p_{\min }={\min }_{(i,m)\in \mathcal{I}\times \mathcal{M}}\left\{p_{i,m}\right\}>0$ being the smallest grade-achievement probability across all part types. For example, given $p_{\min }=10\%$ and $\epsilon =10\%$, the minimum input becomes $\underset{¯}{x}(\epsilon )={(60)}^2/\pi \approx 1,146$ units.

³ That ${\widehat{\mathit{x}}}^{*}(q)=\left[\frac{q}{F({\overline{\mathit{x}}}^{*}(q))}\right]{\overline{\mathit{x}}}^{*}(q)$ is feasible in (*) can be seen as follows. By Eq. (5) the production function F is homogeneous of degree 1 (i.e., $F(\alpha \mathit{x})\equiv \alpha F(\mathit{x})$, for any $\alpha >0$). As a consequence, $F({\widehat{\mathit{x}}}^{*}(q))=F\left(\left[\frac{q}{F({\overline{\mathit{x}}}^{*}(q))}\right]{\overline{\mathit{x}}}^{*}(q)\right)=\left[\frac{q}{F({\overline{\mathit{x}}}^{*}(q))}\right]F\left({\overline{\mathit{x}}}^{*}(q)\right)=q.$ Thus, not only is ${\widehat{\mathit{x}}}^{*}(q)$ feasible in (*), but it also satisfies output efficiency (since $F({\widehat{\mathit{x}}}^{*}(q))\le q$ is in fact binding).

⁴ Perfect matching happens when there is only a single matching class (i.e., for $M=1$); then all parts are matched. Naturally, this indiscriminate approach propagates lax input tolerances into the assembly which drives up other costs (e.g., in terms of rework or rejected output parts for a failure to satisfy tight output tolerances).

⁵ See, e.g., the deviation-per-day (dpd) precision of a mechanical watch movement in Fig. 5, which increases by a factor of 2 (from 35s to 17.5s) when doubling the matching classes (from $M=50$ to $M=100$).

Figure 1: Sorting $x_1=10$ type-1 and $x_2=11$ type-2 parts with characteristics $(s_1,s_2)$ in $\mathcal{S}={\mathcal{S}}_1\times {\mathcal{S}}_2$; random matching within $M=5$ matching classes ${\mathcal{C}}_m={\mathcal{C}}_{1,m}\times {\mathcal{C}}_{2,m}$, for m in $\mathcal{M}=\{1,\dots ,M\}$; three type-1 components and four type-2 components remain unmatched.

Figure 2: Selective assembly system, featuring the sorting of two types of inputs into $M\ge 2$ matching classes according to the part characteristics, random matching across corresponding classes, and subsequent assembly. The scalar output Q is a random variable.

Figure 3: The solution ${\mathit{x}}^{*}(q)$ to the selective matching problem ( $\ast$) can be approximated by scaling the solution ${\overline{\mathit{x}}}^{*}(q)$ to the envelope optimization problem ( $\ast \ast$). This yields the “approximate selective matching solution,” ${\widehat{\mathit{x}}}^{*}(q)=[q/F({\overline{\mathit{x}}}^{*}(q))]{\overline{\mathit{x}}}^{*}(q)$, feasible in ( $\ast$); cf. Exs. 1–3.

Figure 4: (a) Relative error $r({\overline{x}}^{*}(q))$, for $I=2$ and $I=4$, as a function of $q\in [500,10000]$; (b) Ratio of relative errors, for $I=4$ to $I=2$, as a function of $q\in [500,10000]$.

Figure 5: Smallest (even) number of matching classes to guarantee a maximum daily deviation, $\text{dpd}\in \{1,\dots ,60\}$ [s/d]. COSC chronometer standard is achieved for $\text{dpd}\le 5$ [s/d].

Figure 6: Optimal and approximately optimal cost (a) and input (b) as a solution to ( $\ast$) and ( $\ast \ast$), for a target quantity $q\in \{500,1000,2000\}$ units and for an output tolerance $\text{dpd}\le 60$ [s/d].