Space-Filling Designs for Robustness Experiments
V. Roshan Joseph Li Gu Shan Ba William R. Myers
stat.ME 62K25





Space-Filling Designs for Robustness Experiments



V. Roshan Joseph

H. Milton Stewart School of Industrial and Systems
Engineering, 
 Georgia Institute of
Technology,
Atlanta, GA 30332-0205



Li Gu

    Wells Fargo & Company, 

    Charlotte, NC 28202 

    


Shan Ba and William R. Myers

Quantitative Sciences
 The Procter & Gamble Company,
Mason, OH 45040








To identify the robust settings of the control factors, it is very important to understand how they interact with the noise factors. In this article, we propose space-filling designs for computer experiments that are more capable of accurately estimating the control-by-noise interactions. Moreover, the existing space-filling designs focus on uniformly distributing the points in the design space, which are not suitable for noise factors because they usually follow non-uniform distributions such as normal distribution. This would suggest placing more points in the regions with high probability mass. However, noise factors also tend to have a smooth relationship with the response and therefore, placing more points towards the tails of the distribution is also useful for accurately estimating the relationship. These two opposing effects make the experimental design methodology a challenging problem. We propose optimal and computationally efficient solutions to this problem and demonstrate their advantages using simulated examples and a real industry example involving a manufacturing packing line.




KEY WORDS: Computer experiments, Experimental design, Gaussian process, Optimal designs, Quality improvement, Robust parameter design.




§ INTRODUCTION



Robust parameter design is a cost-efficient technique for quality improvement. Originally proposed by <cit.>, the technique has been widely adopted in industries for system (product or process) optimization. The core idea is to first divide the factors in the system into two groups: control factors and noise factors. Control factors are those factors in the system than can be cost-effectively controlled. On the other hand, noise factors are those factors which are either impossible or too expensive to control. For example, in product design of a razor for shaving, blade thickness, gap between the blades, angle of the blades, etc. are control factors, whereas consumer attributes like the skin type, hair length and density, and product usage attributes like the pressure applied on the skin, handle angle, etc. are noise factors. Since the noise factors are uncontrollable, they introduce variability in the performance of the product. Robust parameter design is a technique to find a setting of the control factors (also known as parameter design) that will make the system robust or insensitive to the noise factors. Thus, under a robust parameter design, the output becomes less affected by the noise variability even when the noise factors are left uncontrolled. This is why the approach using robust parameter design is less costly than the other quality improvement techniques, which try to directly control the noise factors in the system.

The key to a successful robust parameter design is in identifying important control-by-noise interactions of the system. Only when such interactions exist we can use the control factors to reduce the sensitivity of the noise factors. These interactions are usually unknown in practice and their existence need to be investigated through experimentation. Thus designing good experiments is a crucial step in robustness studies. Many efficient experimental design techniques are proposed in the literature such as cross arrays <cit.> and single arrays <cit.>. A thorough discussion of these techniques can be found in the books by <cit.> or <cit.>.

The aforementioned experimental design techniques are mainly proposed for physical experimentation except for the work of <cit.>. Recently computer experiments have become very common in industry. That is, if a computer model is available that can simulate the physical system, then the experiments can be performed in computers instead of the physical system. It is becoming very common for industry to develop a computer model for product design like simulating the performance of a razor. One particular example from Procter & Gamble involving the development of a computer model that simulates a critical transformation of a packing line, that involves both control and noise factors, will be discussed in more detail later in the paper. Computer experiments can bring in tremendous cost savings because direct experimentation with the real physical system is always more expensive than investing on some computer time. However, there are several aspects of computer experiments that necessitate the use of a different experimental design technique or philosophy compared to those of physical experiments <cit.>. Since most computer models are deterministic in nature, randomization and replications are not needed. Fractional factorial and orthogonal array-based design techniques that are prevalent in physical experiments lead to replications when projected onto subspace of factors and thus are unsuitable for computer experiments. Split-plot designs that are considered to be useful in robustness studies <cit.> become unnecessary as run orders and restrictions on randomization will not affect the computer model outputs. This led to the development of space-filling designs in computer experiments.

The existing work on robust parameter design using space-filling designs do not make any distinction between control and noise factors. A distinction is made only at the analysis stage <cit.>. Sequential designs that directly attempt to find robust settings of control factors using expected improvement-type algorithms are proposed in the literature <cit.>, but we are not aware of any work on space-filling designs. It is important to develop space-filling designs that distinguish control and noise factors because their distributional properties are entirely different. Control factors are assumed to follow a uniform distribution, whereas noise factors typically follow non-uniform distributions such as normal distribution. Their nature of randomness is also different. Noise factors are intrinsically random and can vary over time and space. On the other hand, control factors remain fixed once their levels are chosen. A uniform distribution is imposed on the control factors only to represent our indifference on the choice of level given the range of possible values for each control factor.  Thus, unlike the control factors, most of the “action” in the noise factor space takes place in the regions of high probability mass. Therefore, space-filling designs that uniformly spread out points in the experimental region are not adequate for robust parameter design experiments. Moreover, the existing space-filling designs are not designed for precise estimation of control-by-noise interactions. In this article we propose a new version of space-filling design that is capable of estimating the control-by-noise interactions more precisely and puts more points in regions that matters the most.

The article is organized as follows. In Section 2, we propose a model-based optimal experimental design for robustness studies. Because of certain computational and practical difficulties associated with this approach, in Section 3, we propose a modified space-filling design as an alternative. In Section 4, we investigate the optimal choice of noise levels for the experiment. Extension of designs to deal with internal noise factors is proposed in Section 5. The proposed methodology is applied in a simulated example and the packing line computer experiment from Procter & Gamble in Section 6. We conclude with some remarks in Section 7.




§ MODEL-BASED OPTIMAL DESIGNS



Let x=(x_1,…,x_p)' be the set of control factors and z=(z_1,…,z_q)' the set of (external) noise factors.  The case of internal noise factors and the differences between the two types will be discussed later. We assume that x ∈𝒳=[0,1]^p and z ∈𝒵, the support of the distribution of z which could be ℝ^q. The response y is a deterministic function of both control and noise factors given by y=g( x, z). Depending on the type of characteristic such as smaller-the-better, larger-the-better, or nominal-the-best, we can impose a quality loss function on y. Let L(y) be such a loss function. Then, the objective of robust parameter design is to find the setting of control factors that minimizes the expected loss, where the expectation is taken with respect to the distribution of noise factors. Let f( z) denote the probability density function of z. Then, the robust parameter design can be obtained by

    min_ x ∈𝒳∫_𝒵 L{g( x, z)}f( z)d z.

Since the function g(·,·) is available only as a computer code, an experiment will be conducted to estimate it. Let D={ x_1,…, x_p,  z_1,…, z_q} be the experimental design with n runs, where x_j=(x_1j,…,x_nj)' and z_k=(z_1k,…,z_nk)' denote the settings of the jth and kth control and noise factors, respectively. Let y_i  be the ith output from the computer model, i=1,…,n.

A Gaussian process model or kriging <cit.> is commonly used for estimating g(·,·). So assume

    g( x, z)∼ GP(μ,τ^2 R(·,·)),

where μ and τ^2 are the unknown mean and variance parameters, and R(·,·) is the correlation function. A commonly used correlation function is the Gaussian correlation function given by

    R( x_i- x_j,  z_i- z_j)=exp{-∑_l=1^pθ^x_l(x_il-x_jl)^2-∑_l=1^qθ^z_l(z_il-z_jl)^2},

where θ^x=(θ^x_1,…,θ^x_p)' and θ^z=(θ^z_1,…,θ^z_q)' are the unknown correlation parameters of the control and noise factors, respectively. Let θ be a p+q column vector containing θ^x and θ^z. We will use this correlation function throughout this article, but other correlation functions are also allowed as long as they can produce smooth realizations of the response in the noise factor space. The smoothness assumption with respect to noise factors is critical for our methodology and we will exploit it for developing the experimental designs. The Gaussian process can be viewed as a prior on the unknown function and therefore, we can obtain its posterior distribution using Bayes theorem <cit.>:

    g( x, z)| y∼ N(ĝ( x, z), τ^2MSE( x, z;  D, θ)),

where ĝ( x, z)=μ+r( x,  z)' R^-1( y-μ 1) and

    MSE( x, z;  D, θ)=1- r( x,  z)' R^-1 r( x,  z),

where y=(y_1,…,y_n)', r( x,  z) is an n× 1 vector with ith element R( x- x_i,  z- z_i), R is an n× n matrix with ijth element R( x_i- x_j,  z_i- z_j), and 1 is a vector of 1's having length n. For simplicity, we chose to ignore the extra variability due to the estimation of μ and θ. See <cit.> and <cit.> for a discussion on the effect of this extra variability on experimental design.

The posterior mean ĝ( x, z) can be used as an estimate of the response function from the experiment (also known as metamodel, surrogate model, or emulator). Then, the optimization in (<ref>) can be simplified as

    min_ x ∈𝒳∫_𝒵 L{ĝ( x, z)}f( z)d z.

It is also possible to incorporate the uncertainties in the estimation of g(·,·) in the optimization as in <cit.> and <cit.>, but it will not be considered here for the sake of simplicity. The problem we are trying to solve is how to design the experiment D so that we can accurately estimate the solution to the optimization problem in (<ref>).

Clearly, the optimization in (<ref>) will give the true robust setting if ĝ( x, z) is the true response surface, that is, if MSE( x, z;  D, θ)=0 for all x∈𝒳 and z∈𝒵. Thus we should design the experiment so that MSE( x, z;  D, θ) is as small as possible. Furthermore, a careful examination of (<ref>) reveals an important insight on the experimental design problem. We need an accurate g(·,·) only in the regions of z where f( z) is large. In other words, if f( z) is small in some regions, then the inaccuracies in the estimation of g(·,·) in those regions will not affect the robust settings. This makes the experimental design problem for robustness different from that of a usual computer experiment. In fact, it makes sense to focus on the estimated solution to the optimization problem in (<ref>) as proposed in <cit.> rather than the estimation of g(·,·). However, their approach works only for linear models fitted to physical experimental data. In contrast, the models considered in computer experiments are highly nonlinear and thus, finding an explicit solution to (<ref>) is not feasible. So in this work we will focus on the estimation of g(·,·).  Although this approach may not look ideal for the robustness objective, it does have certain advantages. The loss functions are many times loosely defined and one may want to investigate solutions to different possible choices of loss function <cit.>. Moreover, in real problems, there can be multiple quality characteristics and thus one may need to be satisfied with a compromise solution, which can be different from the optimal solution in (<ref>). Thus an accurate g(·,·) in the region of interest can be more beneficial than an accurate solution to (<ref>) obtained for a specific choice of loss function and quality characteristic.

Thus, our aim is to find D such that MSE( x, z;  D, θ) is small. However, since MSE( x, z;  D, θ) is a function of x and z, it is not possible to find such a design over the entire experimental region. Instead, a feasible approach is to minimize the average of MSE( x, z;  D, θ), that is

    min_ D∫_𝒳∫_𝒵 MSE( x, z;  D, θ) f( z) d z d x.

This design criterion is the same as the  integrated mean squared error criterion in the literature <cit.> except that we use the density of z as a weight function. This is quite a natural modification of the existing criterion and agrees with our intuition that we should give more weights for regions where f( z) is large. Surprisingly, we found that the solution to (<ref>) places points in extremely low probability regions which are not very useful for finding the robust setting.  This problem can be alleviated if we use  the root-mean squared prediction error, which directly corresponds to confidence intervals of the prediction. Thus, consider a modified criterion

    min_ D IRMSE ( D, θ)=min_ D∫_𝒳∫_𝒵√(MSE( x, z;  D, θ)) f( z) d z d x.

This is a more meaningful criterion as it tries to minimize the expected volume of the confidence region of the predictions. We may generalize this criterion to

    min_ D IRMSE_k ( D, θ)=min_ D[∫_𝒳∫_𝒵{√(MSE( x, z;  D, θ)) f( z)}^k/C_k d z d x]^1/k

for k>0 and C_k=∫_𝒵f^k( z) d z. The special case of k=2 is of great interest as it is analytically tractable in some situations. Let IMSE=IRMSE_2^2. Thus,

    min_ D IMSE ( D, θ)=min_ D∫_𝒳∫_𝒵 MSE( x, z;  D, θ) f^2( z)/C_2 d z d x.

Interestingly, this is the same as the integrated mean squared error criterion in the literature but with a weight function f^2( z). In the case of uniform distributions,
f(·) or f^2(·) doesn't make any difference, but for non-uniform distribution this does
 make a big difference. We will see later that f^2(·) gives the right scaling and provides meaningful solutions to the robust parameter design problem.


A major challenge of using the foregoing criteria is that they are functions of the unknown correlation parameters θ. One can minimize IRMSE (or IMSE) for a guessed value of θ, but the optimal design may not work well for another value of θ. A potential fix to overcome this problem is to first average the IRMSE over a prior distribution of θ and find the design using

    min_ D∫ IRMSE( D,θ)  p(θ)dθ.

This is a computationally intensive problem because R is a function of θ and thus, inverting R and then integrating the IRMSE is time consuming. Moreover, this criterion doesn't work well in practice because MSE( x, z; D,θ) increases with θ and therefore, the criterion is dominated by the large values of θ. <cit.> proposed to overcome this problem by standardizing the criteria with respect to the optimal design obtained for a given value of θ. See <cit.> for a Bayesian version of this approach. Let Θ be a compact set containing the possible values of θ. Then, Sacks et al.'s approach is to find the design to maximize the minimum efficiency:

    max_ Dmin_θ∈ΘIRMSE( D^*(θ),θ) /IRMSE( D, θ),

where D^*(θ)=argmin_ DIRMSE( D, θ). This criterion is extremely computationally intensive because one needs to find the optimal design for every possible value of θ∈Θ, which is difficult in high dimensions. <cit.> tried to circumvent this problem by letting θ_i^x=θ_j^z=θ_0 for all i=1,…,p and j=1,…,q and then taking a few discrete values of θ_0. In our experience, this simplification results in designs having poor projections in subspaces of the factors, which is undesirable. One possible approach to improve the projections is to assign independent prior distributions for each of the unknown correlation parameters as in <cit.>. However, this would make the computation of (<ref>) very expensive. For example, even if we choose only three discrete values for each correlation parameter, the total number of possible values in Θ would become 3^p+q which can be quite high when p and/or q is large.

Because of the foregoing difficulties, in the next section we will develop space-filling designs which are easy to compute and are model-robust. We will modify them so that they will perform well according to the criterion in (<ref>).



§ SPACE-FILLING DESIGNS

Space-filling designs aim at filling the experimental region evenly with as few gaps as possible. These designs are robust to modeling choices and thus, are widely used as designs for computer experiments.  Popular space-filling designs include Latin hypercube designs <cit.>, distance-based designs such as maximin and minimax <cit.>, uniform designs <cit.>, and several useful variants of them such as maximin Latin hypercube designs <cit.> and maximum projection designs <cit.>.  See <cit.> for a recent review of space-filling designs. However, these designs are developed for general purpose applications such as function approximation and not specifically for robustness experiments. As mentioned in the introduction, control-by-noise interactions are especially important for identifying robust settings. Therefore, we may hope to improve the performance of space-filling designs by improving their ability to estimate the control-by-noise interactions, possibly by sacrificing other not so important effects.

In the physical experiments' literature, there are mainly two classes of designs suitable for robustness experiments: cross arrays <cit.> and single arrays <cit.>. To develop single arrays, one needs to first quantify the importance of each effect <cit.>. This was not too difficult with the fractional factorial experiments because such designs usually have only two or three levels for each factor. In contrast, computer experiments have large number of levels for each factor and therefore, numerous effects are involved in the modeling. This makes the effect ordering a difficult task. Thus, cross arrays seems to be an easier and straightforward approach for computer experiments.

To develop cross arrays, we first need to choose an n_1-run design D_x for the control factors (known as control array) and another n_2-run design D_z for the noise factors (known as noise array). Cross array can then be obtained by repeating the noise array for each run of the control array, which will have a total of n_1n_2 number of runs. We will denote the cross array by D= D_x × D_z. The suitability of a cross array for robustness experiments should be evident from its construction. We can estimate the effect of noise factors under each settings of the control factors in D_x, which enables one to choose the control factor setting that makes the noise factors' effect on the response as small as possible. As shown in <cit.> (see Theorem 11.1), if control factor effects are estimable from D_x and noise factor effects from D_z, then the two-factor interactions between control and noise factors are estimable and clear in D.

The following result shows how to construct an optimal cross array that minimizes the IMSE criterion in (<ref>). The proof is given in the Appendix.

If a product correlation is used between control and noise factors, then an IMSE-optimal cross array can be obtained by crossing an IMSE-optimal control array and an IMSE-optimal noise array.

The IMSE optimal designs are closely related to the space-filling designs. For example, when f( z)=1, as k→∞ in (<ref>),

    IRMSE_∞( D)=max_𝒳max_𝒵MSE( x, z; D),

which is minimized by a minimax distance design when the correlations are small (Johnson et al. 1990). The assumption of small correlations is justifiable in the control factor space, but not in the noise factor space. Moreover, Noise factors have nonuniform distributions. We will discuss on how to modify the space-filling designs for the noise array in the next section.

A simple example can be used to illustrate why the cross arrays are useful in estimating control-by-noise interactions. Consider two control factors (x_1 and x_2) and two noise factors (z_1 and z_2). Suppose we choose a Maximin Latin hypercube design (MmLHD) with four runs for the control array and another MmLHD with five runs for the noise array. Then the cross array will have 20 runs. Their two-factor projections are shown in Figure <ref>. We can see that the projections are very poor in the control factor space and the noise factor space, but the projections are excellent on the control-by-noise factor space. Thus, we can obtain a good estimation of the control-by-noise interactions using this design. However, the estimation of the pure control or noise factor effects can be poor compared to using a 20-run MmLHD for the four factors.

./Figures/


Although the cross array looks promising for robustness experiments, it has certain disadvantages for using in computer experiments. Each control factor level is replicated n_2 number of times and noise factor level n_1 number of times. Thus, the number of levels is much smaller than that of a comparable single array, which has n_1n_2 number of levels for each factor. This can lead to poor estimation of nonlinear effects and higher order interactions in the control and noise factor spaces. This disadvantage of the cross array is amplified if there are only a few factors that are active. We propose an idea to overcome this disadvantage.

We can jitter each point in the cross array to increase the number of levels for each factor. The resulting design will still posses the good estimation ability of control-by-noise interactions because the response values observed over adjacent points are expected to be highly correlated. However, if the jittering radius is small, then there will be no improvement in the estimation of the pure control and noise factor effects. On the other hand, if the jittering radius is large, then the resulting design can lose its ability to efficiently estimate the control-by-noise interactions. So how much to jitter is a critical question and we have an intuitive solution for this. Let r be the fill distance or covering radius of the cross array, which is defined as

    r=max_ u ∈ [0,1]^p+q || u-Q( u, D)||,

where Q( u, D) is the closest point in D to u and ||·|| denotes the Euclidean distance. Let B_i denote the ball with center at the ith design point D_i and radius r and let C_i be the hypercube inscribed in the ball. Now we can jitter the ith point within C_i. We chose C_i instead of B_i because in high dimensions the projected points inside the ball can be far away from the center. We will call the resulting design a Jittered Cross Array (JCA).

We can do better than a random jittering. Since our aim is to overcome the issue with projections, we can choose the points in C_i that will ensure good projections. We propose a sequential algorithm for doing this. Start from the center point (.5,…,.5) after rearranging the rows in the cross array so that the first run is closest to the center point. Now add one point from C_i sequentially using the maximum projection (MaxPro) criterion <cit.>:

    D_i=min_ u ∈𝒞_i∑_j=1^i-11/∏_l=1^p|u_l- D_jl|^s,

for i=2,…,n and s=2. This algorithm adds points sequentially in a greedy manner such that the ith point is as far as possible from the previously chosen points under the MaxPro criterion. The algorithm is very easy to implement, but can converge to a local optimum. To improve the performance, we repeat this procedure many times by randomizing the order in which the C_i's are chosen. Although the MaxPro criterion ensures that no two levels can be the same, the levels may not be equally spaced. So at the end of each iteration, we force them to be equally spaced, which can be easily done by ordering the levels. So the final design is like a Latin hypercube design (LHD), but with some clustering in the control and noise factor spaces. JCA shouldn't be confused with a cascading LHD <cit.>, which has clusters in the full-dimensional space and doesn't have a crossed array structure.

The two-dimensional projections of the JCA for the previous example is shown in Figure <ref>. We can now see the 20 points in the x_1× x_2 and z_1× z_2 projections as opposed to only four and five points in the cross array. The projections in the control-by-noise spaces are still very good. We can also observe the clusters in the control (symbols) and noise (color) spaces, which shows that the cross array structure is approximately maintained. We will study the performance of these designs using simulated and real examples in a later section after deciding the optimal choice of noise array.
./Figures/




§ NOISE ARRAY


In this section we will discuss three possibilities for the choice of noise array. The first one is the most intuitive choice, but the latter two are better for robustness experiments.



 §.§ Transformed Design

As mentioned earlier, the distribution of the noise factors are usually nonuniform. On the other hand, most space-filling designs are closely related to a uniform distribution as they try to spread out the points evenly in the experimental region. This suggests that we can possibly use the inverse probability transform method to transform the noise factor columns in a space-filling design to have the right distribution. More specifically, assume that the noise factors are independent. Let F_l(z_l) be the distribution function of z_l, l=1,…,q. Then it can be shown that if the design { z_1,…, z_q} minimizes the discrepancy from a uniform distribution, then {F_1^-1( z_1),…,F_q^-1( z_q)} minimizes the F-discrepancy <cit.>, where the transformation is applied element-wise. We will call this design a transformed design.

The independence assumption is crucial for the above simplification. This assumption can be easily relaxed when the noise factors follows a multivariate normal distribution: N( 0,Σ_z). We can first find the space-filling design assuming independence and then transform using Σ^1/2Φ^-1( z^i), where Φ(·) is the distribution function of a standard normal variable and z^i is the ith row of the space-filling design. If the distribution is not normal, then the design can be found using the idea of support points <cit.>.

For illustration, let z be a normal distribution with mean 0.5 and standard deviation σ=1/6. These values are chosen so that the ± 3σ limits of z  coincide with [0,1]. Suppose we use n=10 and a Gaussian correlation function: R(z_i-z_j)=exp{-θ (z_i-z_j)^2}. As shown in <cit.>, d_0={.5/n, 1.5/n,…,(n-.5)/n} minimizes the discrepancy from a uniform distribution. So the desired design can be obtained as F^-1( d_0). Figure <ref> plots the weighted root mean squared error

    WRMSE( x, z)=√(MSE( x, z))f( z),

for θ=10 and θ=1000. The design points are shown as crosses in the same plots. We can see that the points are pulled towards the center as one would expect for a normal distribution. However, although the transformation seems to balance the WRMSE well throughout the region when θ=1000, it seems to be too high in the tail regions when θ=10. Thus, the benefit of using a transformed design seems to depend on the smoothness of the underlying response function. If the function is wiggly, then transformation will work great, but if the function is smooth, the transformation may do more harm than good. From our experience, the external noise factors usually have  a smooth relationship with the response. This could be because the realistic ranges of noise factors are much smaller than the possible ranges of control factors and therefore, the noise-response relationship can be adequately modeled using a smooth Gaussian process. Thus, θ is expected to be small in the Gaussian correlation function. In summary, a transformed design does not seem to be a good choice for the noise array and we need to look for other alternatives.

./Figures/



 §.§ Hybrid Design

We can find the design that minimizes IRMSE in (<ref>) for a given value of θ or even better by maximizing the efficiency in (<ref>). However, as mentioned before, this optimization is hard to perform, especially in high dimensions. We propose a simple idea to overcome this problem. We will find the  optimal design for one factor, which is easier. The optimal design can be viewed as a transformation of a uniform design. Now we use this  optimal transformation on each column of the noise factors of a space-filling design. Since the final design is obtained using a combination of space-filling  and model-based optimal design criteria, we will call the design a hybrid design.

To fix the idea, consider a single noise factor z. Let d^*={z_1^*,…,z_n^*} be the optimal design obtained using the model-based criterion

    max_ dmin_θ∈ΘIRMSE( d^*(θ),θ) /IRMSE( d, θ),

where

    d^*(θ)=argmin_ d∫_𝒵√(MSE(z;  d, θ)) f(z) dz.

We want to emphasize that these optimizations are computationally much simpler than those of (<ref>) and (<ref>) which use the full factor space.
Let {u_i^*=(i-.5)/n, i=1,…,n} be the uniform design points in [0,1]. It is easy to show that the optimal design points are distinct, that is z_i^* z_j^* for i  j. Therefore, there exists a one-to-one transformation:

    z_i^*=T(u_i^*),  for  i=1,…,n.

For the case of multiple factors, let U_z be a space-filling design with levels {u_1^*,…,u_n^*} for each of the noise factor. Then, the hybrid design can be obtained as

    D^*_z=T( U_z),

where the transformation T(·) is applied on each element of U_z.

Consider again the same example of the previous subsection with z∼ N(.5,1/6) and θ=10. The WRMSE for the 10-run optimal design is plotted in Figure <ref>. We can see that the optimal design performs much better than the  transformed design. However, this improvement is not realizable in practice because we never know the true value of θ before the experiment. We need to choose a robust value of θ using (<ref>).

./Figures/


Consider a set Θ={5,10,20,30} with z∼ N(.5,1/6) and n=50. We first find the optimal designs for each of the four values of θ∈Θ. A histogram of the points for one case is shown in Figure <ref>. We can see that the IRMSE-optimal design points are slightly more dispersed than the noise distribution. The efficiencies of the four designs are computed using (<ref>) and are plotted over θ in Figure <ref>. We can see that the optimal design found based on smaller values of θ perform poorly for larger values of θ. On the other hand, the optimal designs found using larger values of θ perform not so poorly for smaller values of θ. This suggests that we should find the set of probable values of θ and use the largest value in that set to generate the optimal design.

./Figures/


./Figures/


Guessing the largest possible value of θ can be challenging in a practical problem. Moreover, the numerical inaccuracies in computing IRMSE increases with n and the optimization becomes harder and unstable. Due to these difficulties, in the next subsection, we will try to identify an approximate optimal design that is easier to use in practical applications.



 §.§ Double Transformed Design

Let d^*={z_1^*,…,z_n^*} be the optimal set of n points that minimizes the IRMSE in (<ref>) for a given value of θ. The empirical distribution function of this point set is given by F_n(z)=1/n∑_i=1^n I(z_i^*<z), where I(·) is the indicator function. Our aim is to understand the limiting distribution of F_n(z) as n→∞, which we denote by F̃(z). This is not an easy problem because there is no explicit solution for the optimal design. Moreover, IRMSE is a complex function of the design points. Therefore, we will make use of an existing result on optimal designs for uniform distributions to get an idea of the limiting distribution.

For uniformly distributed variables, <cit.> showed that a beta distribution with density b(z; α,α) for α∈ [0.5,1] is optimal for a reciprocal distance criterion which can be viewed as a surrogate for (<ref>). See also <cit.> for a rigorous justification of this result. We know that if F(·) is the distribution function of z, then F(z)∼ U(0,1). Thus, by using change of variables, the optimal density of the design can be obtained as

    f̃(z)   =   b(F(z); α,α)f(z)
       =   Γ(2α)/Γ^2(α)f(z)/{F(z)[1-F(z)]}^1-α.

Let B_α(z)=∫_0^z b(u;α,α) du be the distribution function of the beta distribution. Then, the asymptotic distribution function of the optimal design is given by

    F̃(z)=B_α(F(z)).

Thus, if d_0 denotes the uniform design, then an approximation to the optimal design can be obtained as

    d^*=F^-1(B_α^-1( d_0)).

To distinguish from the previous transformed design F^-1( d_0), we will call this the double transformed design.

What should be the value of α? If θ is large in the Gaussian correlation function, we should use α=1, which leads to the transformed design. We have seen in Figure <ref> that the transformed design indeed works well with large θ. However, as mentioned earlier, we are more interested in small values of θ. <cit.> recommended using α=1/2, which is the limiting distribution of a D-optimal design for large degree polynomial regression. Theorem 10.1 in <cit.> shows that as the correlation parameter θ→ 0, the Gaussian process predictor tends to a high degree polynomial interpolator and therefore, it makes sense to use an α value close to 1/2. The design points, in this case, are the same as Chebyshev nodes, which possess minimax optimality properties for polynomial interpolation <cit.>. However, as discussed towards the end of previous subsection, the optimal design based on a small value of θ may work poorly for large values of θ. Thus, intuitively, a value of α slightly larger than 1/2 such as 2/3 or 3/4 might be a more robust choice. We investigate this more carefully below for the case of a normal distribution.

It is easy to show that MSE(z; d)=1- r(z)' R^-1 r(z)≤ 1-R^2(z-Q(z, d)), where Q(z, d) is the closest point in d to z. Let

    IRMSE( d)=∫_𝒵√(1-R^2(z-Q(z, d)))f(z)dz,

which is an upper bound of IRMSE( d). For Gaussian correlation function,

    IRMSE( d)   =   ∫_𝒵√(1-exp{-2θ(z-Q(z, d))^2})f(z)dz 
       ≈   √(2θ)∫_𝒵|z-Q(z, d)|f(z)dz,

where the approximation is valid for large n. <cit.> has shown that the design that minimizes ∫ |z-Q(z, d)|^kf(z)dz has an asymptotic distribution proportional to f^1/(1+k)(z). This implies that the asymptotic distribution of d that minimizes IRMSE( d) should be proportional to √(f(z)). Now consider the case of a normal distribution f(z)=ϕ(z;.5,σ). Based on the simple approximation given by <cit.>, F(z)[1-F(z)]≈ f^4/π(z). Substituting this approximation in (<ref>), we obtain f̃(z)≈Γ(2α)/Γ^2(α) {f(z)}^1-4(1-α)/π. This will be proportional to √(f(z)) if α=1-π/8≈ 0.607. Similar exercise using IMSE gives α=1-π/12≈ 0.738. Based on these values, we choose α=2/3 ≈ (.607+.738)/2. Interestingly, in a totally different problem setting of searching for the maximum of a continuous function using non-adaptive algorithms, <cit.> showed that the same beta density with α=2/3 is optimal.

./Figures/

Consider the previous example with z∼ N(.5,1/6).  The WRMSE for the double transformed design using (<ref>) with α=2/3 is shown in Figure <ref>. We can see that WRMSE for the double transformed design is much smaller than that using the transformed design, but not as good as the optimal design in Figure <ref>.

The left panel of Figure <ref> compares the noise distribution with the asymptotic optimal density in (<ref>). We can see that the density for the optimal design is more dispersed than the original noise distribution.  The right panel of Figure <ref> shows the optimal density when the noise distribution is truncated to [0,1]. This density has three modes: one at the center and two at the boundaries. This is an interesting result, because Taguchi (1987) has recommended using three levels for the noise factor, one at the mean and two at the extremes. The result in Figure <ref> can be viewed as an extension of this three-level design for physical experiment to an n-level design for computer experiment, where projections are important.

./Figures/


A similar investigation using the initial IMSE criterion in (<ref>) shows that α≈ 0.476 is the optimal choice. However, this makes the optimal density quite dispersed and places points in very low probability regions. For example, when n=100, the points can be as far as ± 3.95σ from the center. On the other hand, α=2/3 places points within ± 3.25σ from the center, which looks more reasonable. This is why we feel the IRMSE criterion in (<ref>) or the IMSE criterion in (<ref>) is more meaningful than the IMSE criterion in (<ref>). This was also verified using the prediction performance on some test cases.




§ FACTORS WITH INTERNAL NOISE

There are some factors whose nominal values can be controlled, but they can vary around their nominal values. Such factors are said to have internal noise. Examples include, part-to-part variability within their manufacturing tolerances and process parameter variability around its target. On the other hand, external noise factors are completely uncontrollable including their nominal values. Examples of external noise factors include user conditions, incoming raw material properties,  etc. In this section we will propose methods for designing experiments with internal noise factors. One may wonder why we need to consider these factors differently from the external noise factors. Why not just merge them with the external noise factors and use the techniques described in the previous two sections? The reason is that we don't need to vary the internal noise factors in the experiment! They can be easily introduced at the modeling stage. Thus, internal noise factors can be ignored at the design stage although this may not be the “optimal” approach. This topic has received scant attention in the literature except possibly for the work of <cit.> for the case of physical experiments.

A factor with internal noise can be represented as X=x+e, where the nominal value x is controllable and the internal noise e is uncontrollable. Here we have used additive noise, but the case of multiplicative noise can be handled similarly. Suppose we have estimated the relationship with the nominal values of the factors (x) and external noise factors (z): ĝ( x, z), then we can easily obtain the relationship with the internal noise as ĝ( x+ e, z). This is why we don't need to vary e in the experiment.

Since a factor with internal noise is both a “control”  and a “noise” factor, it make sense to cross this factor with the other factors in the experiment. Thus, if D_X denotes the design for the factors with internal noise, then we can obtain the cross array using D_x × D_z × D_X. From this we can obtain the jittered cross array using the same algorithm discussed earlier. Now, we only need to decide on how to choose the levels for a factor with internal noise.

As before let's assume x to follow a uniform distribution in [0,1] and e to have a noise distribution with density f_e(e). Because the internal noise factor is the result of not controlling the process well, it is mostly going to have a normal distribution.  So let e∼ N(0,σ_e). The optimal design d^*={x_1^*,…,x_n^*} can be obtained by minimizing

    IRMSE( d,θ)=∫_0^1∫_-∞^∞√(1- r(x+e)' R^-1 r(x+e))ϕ(e;0,σ_e)de dx,

where ϕ (e; 0,σ_e) denote the density of a normal distribution with mean 0 and standard deviation σ_e. Here it is better to consider the IMSE criterion in (<ref>) because an explicit expression for the integral can be obtained under a Gaussian correlation function, R(h)=exp(-θ h^2). Thus,

    IMSE( d,θ)   =   ∫_0^1∫_-∞^∞{1- r(x+e)' R^-1 r(x+e)}ϕ^2(e;0,σ_e)/C_2de dx
       =   1/2σ_e√(π)C_2[1-∫_0^1 tr{ R^-1 A(x)}dx],

where the ijth element of A(x) is given by

    A_ij(x)=1/√(1+2θσ_e^2)exp{-2θ/1+2θσ_e^2(x-x_i+x_j/2)^2}exp{-θ/2(x_i-x_j)^2},

for i,j=1,…,n. The integration with respect to x can also be done explicitly to obtain

    IMSE( d,θ)=1/2σ_e√(π)C_2[1- tr{ R^-1 A}],

where

    A_ij=√(π)/√(2θ){Φ( √(θ)(2-x_i-x_j)/√(1+2θσ^2_e))-Φ( -√(θ)(x_i+x_j)/√(1+2θσ^2_e))}exp{-θ/2(x_i-x_j)^2},

where Φ(·) is the standard normal distribution function. Similar explicit expressions could have been obtained in the previous section as well for the external noise factors, but we didn't do it because of numerical issues. We found the formula in (<ref>) to be quite vulnerable to numerical issue when R is nearly singular, which happens when θ is small. For the external noise factors, we need to consider small values of θ because they are expected to have a smooth relationship with y. On the other hand, the relationship with x can be quite rough and therefore, here we should use large values of θ which doesn't lead to numerical problems.


As an example, consider a factor with internal noise distribution e∼ N(0,σ_e). Let σ_e=1/12 and θ=50 in the Gaussian correlation function. A 10-point optimal design is obtained by numerically minimizing (<ref>). Figure <ref> plots the expected mean squared error for the 10-point uniform design and the optimal design. We can see that the optimal design is almost equally spaced but  with points placed at the boundaries. Our simulations show that {0,1/(n-1),2/(n-1),…,1} is close to optimal, which is only a slight change from the uniform design.

./Figures/




§ EXAMPLES



 §.§ A Simulated Example

Consider a simple example with one control and four external noise factors. Let

    y=∑_i=1^4β_i(x-γ_i)z_i^2e^-(x-γ_5)^2,

where z_1,z_2∼^iid N(0,1). We constructed a cross array using n_1=6 equally spaced levels for the control factor and an MmLHD with n_2=9 runs for the noise factors. From this a JCA is obtained using the sequential MaxPro algorithm described in Section 3. For comparison, we also constructed a MaxProLHD in 54 runs. We consider two versions of these two designs using transformed noise array and double transformed noise array. The resulting four designs are denoted as TrMaxProLHD, DTMaxProLHD, TrJCA, and DTJCA. Two-dimensional projections of DTJCA is shown in Figure <ref>. We simulated 200 cases by randomly sampling β_i's and γ_i's from U(0,1) and kriging models were fitted using the data generated by each of the three designs. 100 test samples were generated using a scrambled Sobol sequence with the noise factor columns transformed using the inverse distribution function of the noise distribution. The root-mean squared prediction errors from the kriging models are plotted in the left panel of Figure <ref>. We can see that the the double transformation on the MaxproLHD and JCA has helped to improve the prediction errors.

./Figures/


./Figures/


Suppose our aim is to minimize the variance of the response due to the noise factors. We have computed the true robust setting of the control factor (x^*) by minimizing the variance and also the robust settings obtained from the fitted models based on the four designs (x_1^*,x_2^*,x_3^*, and x_4^*). The right panel of  Figure <ref> shows the density plots of the errors x_i^*-x^*, for i=1,2,3,4. We can see that the DTJCA gives the best performance followed by TrJCA and DTMaxProLHD. Thus, although JCA didn't help improve the prediction errors over the existing TrMaxProLHD, it does seems to improve the identification of the robust setting.




 §.§ A Real Example

Computer experiments with noise factors are quite common for the simulations conducted at The Procter & Gamble Company. The specific example we will use in our study involves a manufacturing packing line.  The example has been slightly modified for the benefit of simplicity and to prevent disclosure of any potential sensitive information.  A computer simulator was developed for one critical transformation of the packing line involving transport of the package for product fill. A computer experiment with nine input factors was performed and an emulator was built. In this study, we use this emulator for investigating the robustness. Variables x_1, x_2, x_3, x_4, x_5, and  x_6 are process variables such as speed of the line and dimensions of the puck that transports the package that will remain fixed or easy to control once they are chosen and therefore are defined as control factors. Variables z_1, z_2, and z_3 are material properties  of the packaging component such as density and modulus which are defined as noise factors given that there is variation in normal production of the material supplier.  The output response from the computer simulation measures the deflection (deviation from a vertical orientation) after the package holder stops on the packing line, which impacts the quality of the given packing line transformation.  Figure <ref> illustrates the “passing” and “failing” scenarios of the package holder. The objective of this study is to find the settings for the six control factors that are robust to the variation of the three noise factors. From historical data, the noise factors are found to be approximately normally distributed. After re-scaling, we let z_i∼^iid N(.5,σ) for i=1,2,3 with σ=1/6.

./Figures/




First we generated an MmLHD with n_1=13 runs for the six control factors and an MmLHD with n_2=7 runs for the noise factors. The JCA in 91 runs is thus obtained using the sequential MaxPro algorithm and then performed the double transformation on the noise factor columns using (<ref>). We also generated a MaxProLHD in 91 runs for comparison and transformed using the noise distribution. Kriging models were fitted to the data generated from the two designs. We found the prediction errors from the two fitted models to be close, but there was some major differences in the estimation of control-by-noise interactions. Figure <ref> shows the interaction between x_2 and z_3, which is the most significant interaction in this experiment. Clearly, the new double transformed JCA did a much better job in accurately estimating the interaction than the existing transformed MaxProLHD.

./Figures/




§ CONCLUSIONS

In this paper we have proposed space-filling designs that are suitable for identifying robust settings using computer experiments. The key idea was to modify the well-known cross array designs using a space-filling criterion such as the maximum projection criterion. We have also proposed how to optimally choose the noise array. The most intuitive way to construct the noise array is to transform the columns of the array using the inverse cumulative distribution function of the noise factors. This will pull the design points to the high probability region of the noise distribution. However, we found this “pulling effect” to be too extreme. This was mainly because the response is usually a smooth function of the noise factors and to precisely estimate a smooth function it is desirable to push the points outward from the center. The optimal design balances this “pulling” and “pushing” effects. We found that pushing the points uniformly distributed in the unit interval using a Beta(2/3,2/3) distribution before applying the inverse probability transform to be close to asymptotically optimal for normally distributed noise variables. We have also proposed model-based methods to obtain the optimal transformation for any noise distribution, but it requires specification of the upper bound of certain correlation parameters.




    ACKNOWLEDGMENTS


This research is supported by a U.S. National Science Foundation grant DMS-1712642 and a U.S. Army Research Office grant W911NF-17-1-0007.



APPENDIX: Proof of Theorem 1


Consider a production correlation R( x_i- x_j, z_i- z_j)=R_x( x_i- x_j)R_z( z_i- z_j). Now, if D= D_x × D_z, then R= R_x ⊗ R_z and r( x, z)= r_x( x) ⊗ r_z( z), where ⊗ denotes Kronecker product (see, for example, Hung et al. 2015). Then, by using the properties of Kronecker products

    r( x, z)' R^-1 r( x, z)    =    ( r_x( x) ⊗ r_z( z))'( R_x ⊗ R_z)^-1 r_x( x) ⊗ r_z( z)
       =   ( r_x( x)' ⊗ r_z( z)')( R_x^-1⊗ R_z^-1) r_x( x) ⊗ r_z( z)
       =   ( r_x( x)' R_x^-1⊗ r_z( z)' R_z^-1) r_x( x) ⊗ r_z( z)
       =   ( r_x( x)' R_x^-1 r_x( x))⊗ ( r_z( z)' R_z^-1 r_z( z))
       =    r_x( x)' R_x^-1 r_x( x) r_z( z)' R_z^-1 r_z( z).

Thus,

    IMSE( D)    =   ∫_𝒳∫_𝒵{1- r_x( x)' R_x^-1 r_x( x) r_z( z)' R_z^-1 r_z( z)} f^2( z)/C_2 d z d x 
       =    1-∫_𝒳 r_x( x)' R_x^-1 r_x( x)d x ∫_𝒵 r_z( z)' R_z^-1 r_z( z)f^2( z)/C_2 d z
       =    1-(1-IMSE( D_x))(1-IMSE( D_z)).

Thus, min_ D IMSE( D)=1-(1-min_ Dx IMSE( D_x))(1-min_ D_zIMSE( D_z)).
ECA_jasa


