Bayesian modeling and optimization for multi-response surfaces

Computers & Industrial Engineering 142 (2020) 106357

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

Bayesian modeling and optimization for multi-response surfaces a

a

b

c,⁎

Jianjun Wang , Yizhong Ma , Linhan Ouyang , Yiliu Tu

T

a

School of Economics and Management, Department of Management Science and Engineering, Nanjing University of Science and Technology, Nanjing, Jiangsu 210094, China College of Economics and Management, Department of Industrial Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing, Jiangsu 210094, China c SCHULICH School of Engineering, Department of Mechanical and Manufacturing Engineering, University of Calgary, University Drive 2500, N.W., Calgary, Alberta T2N 1N4, Canada b

ARTICLE INFO

ABSTRACT

Keywords: Multi-response optimization Quality loss function Bayesian analysis Conformance probability Model parameter uncertainty Seemingly unrelated regression

This paper proposes a new Bayesian modeling and optimization method for multi-response surfaces (MRS). The proposed approach not only measures the conformance probability (i.e., the probability of all responses simultaneously falling within their corresponding specification limits) through the posterior predictive function but also takes into account the expected loss (i.e., bias, quality of predictions and robustness) with the expected loss function. Also, it is shown that the Bayesian SUR models can provide more flexible and accurate modeling than the standard multivariate regression (SMR) models with the same covariate structure across different responses. Besides, the proposed approach also takes into account the correlation structure of the response data, the variability of the process distribution, and the uncertainty of model parameters as well as the prediction performance of the response model. A Polymer experiment and a simulation experiment are used to demonstrate the effectiveness of the proposed approach. The comparison results show that the Bayesian SUR model has higher conformance probability and lower expected loss than the Bayesian SMR model.

1. Introduction Response surface methodology (RSM) is a collection of statistical design, empirical modeling methodologies, and numerical optimization techniques used to optimize processes and product design (Myers, Montgomery, & Anderson-Cook, 2009). RSM, broadly understood, has become the core of industrial experimentation (Myers, Montgomery, Vining, Borror, & Kowalski, 2004). In industrial experiments, the RSM can establish a functional relationship between the input factors and the output responses. Moreover, the constructed response surfaces can analyze how the input factors affect the output responses and can obtain the optimal parameter settings. Therefore, the RSM method has become an increasingly important role for continuous quality improvement activities. Numerous industrial experiments often involve more than one quality characteristic of interest. A common problem in quality design is how to optimize multiple response surfaces simultaneously, which is called a multi-response surfaces (MRS) optimization problems (Kim & Lin, 2006). In MRS optimization problems, there exists a series of research issues such as correlation among multiple responses (Chiao & Hamada, 2001; Zhang, Ma, Ouyang, & Liu, 2016), robustness measurement of multivariate process (He, Zhu, & Park, 2012), making

⁎

trade-offs among multiple optimization goals (Shin, Samanlioglu, Cho, & Wiecek, 2011), uncertainty of model parameters or model form itself (Ng, 2010; Ouyang, Ma, Wang, & Tu, 2017; Wang, Ma, Ouyang, & Tu, 2016), prediction performance of process model (Ouyang, Ma, Byun, Wang, & Tu, 2016) and reliability assessment for optimization results (Peterson, 2004; Peterson, Miro-Quesada, & Del Castillo, 2009). Indeed, it is very intractable to consider these issues above simultaneously in a unified framework of modeling and optimization. In the last three decades, various innovative approaches have been proposed to solve the above MRS problems (Murphy, Tsui, & Allen, 2005; Myers, Brenneman, & Myers, 2005). Generally speaking, MRS problems include three different stages: performance index construction, response surface modeling, and model parameter optimization (Wang et al., 2016). In MRS optimization problems, one of the most popular approaches is to use the dimensionality reduction strategy. This strategy usually finds a simplified performance index that can transform an MRS optimization problem into a single-objective optimization problem. Several common simplified performance indices have been proposed for MRS optimization problems such as the desirability function (Del Castillo, Montgomery, & McCarville, 1996; Goethals & Cho, 2012; He et al., 2012; Jeong & Kim, 2009), the loss function (Ko, Kim, & Jun, 2005;

Corresponding author. E-mail addresses: [email protected] (J. Wang), [email protected] (Y. Ma), [email protected] (L. Ouyang), [email protected] (Y. Tu).

https://doi.org/10.1016/j.cie.2020.106357 Received 25 October 2018; Received in revised form 31 October 2019; Accepted 9 February 2020 Available online 14 February 2020 0360-8352/ © 2020 Elsevier Ltd. All rights reserved.

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

Ouyang et al., 2017; Ouyang, Ma, & Byun, 2015) and the posterior predictive function (Peterson et al., 2009; Peterson, 2004; Wang et al., 2016). The desirability function method effectively solves the conflict and makes a trade-off between multiple optimization goals according to the types of different quality characteristics such as the smaller the better, the larger the better and the nominal is the best. The desirability function method is widely used by some quality engineers because of its simplicity and ease of operation and is embedded in some typical quality professional analysis software such as MINITAB. However, there are some obvious deficiencies in the commonly used desirability function methods. For example, the variance-covariance relationship between responses, the correlation between multiple responses, and the predictive performance of responses are often not adequately considered in these standard desirability functions. Harrington (1965) firstly proposed a desirability function, which gives a numerical value between zero and one for quality characteristics of a product or process. Derringer and Suich (1980) proposed an improved desirability function to solve the conflict among multiple responses through the relative weights for multiple responses. However, the improved desirability function approach does not take into consideration the correlations among multiple responses, and the variance-covariance structure of the responses as well as the variability of the predicted responses. Wu (2005) proposed an approach to optimizing the multiple correlated responses based on the modified double-exponential desirability function. He, Wang, Oh, and Park (2010) proposed a robustness desirability function, which makes a trade-off between robustness and optimization for MRS optimization problems. Also, Goethals and Cho (2012) extended the above robustness desirability function to account for the variability among multiple responses by fitting high-order variance and covariance models. However, high-order response surface models can result in over-fitting problems (Zhou, Ma, Tu, & Feng, 2013). Vining and Bohn (1998) stressed that the process variance is a “noise” system, which results in poor fit or prediction results. Besides, the response variability associated with the mode parameter uncertainty can have an impact on the optimization results for the MRS problem. Chapman, Lu, and Anderson-Cook (2014) proposed an alternative for desirability function approaches to consider the impact of response variability and estimate uncertainty on solution selection, which can facilitate improved decision making for MRS problems. He et al. (2017) developed the robust fuzzy desirability sets and constructed the robust membership functions of the responses, considering the model uncertainty, location, and dispersion effects of the responses simultaneously. Another simplified approach to tackle MRS optimization problems is the use of loss function. The loss function was initially introduced into the field of quality engineering by the Japanese quality expert Genichi Taguchi. Taguchi has defined “loss to society” as those losses incurred when a product’s characteristics do not meet the customer’s requirements, which are related to the deviation and the variation of products or processes. Taguchi’s quality loss function was regarded as a breakthrough in describing the quality and helped fuel the continuous improvement that since has become known as lean manufacturing. Hence, Taguchi’s loss function can help quality engineers better understand the importance of designing for variation. Pignatiello (1993) extended Taguchi’s univariate loss function (Taguchi, Elsayed, & Hsiang, 1989) to consider the correlation problem among multiple responses. Vining (1998) extended Pignatiello’s method to further take into account the correlation among multiple responses, process economics, and the quality of model predictions. Furthermore, Ko et al. (2005) proposed a new loss function method combining the strengths of Pignatiello’s and Vining’s methods, which allows analysts to consider both the robustness and the quality of predictions as well as the bias for MRS optimization problems in a single framework of the loss function. The significant advantage of Ko et al.’s loss function approach is its ability to incorporate the quality of predictions and process economics as well as the variance-covariance structure of the responses (Ouyang et al., 2015). As noted by Ko et al. (2005), their loss function approach

neglects the robustness of process parameter fluctuations. Peterson (2004) pointed out that some of the quadratic loss functions do not take into consideration the model parameter uncertainty associated with the variance-covariance matrix of the error terms. Ouyang et al. (2015) proposed an integrative loss function approach to simultaneously consider the location and dispersion performances of squared error loss to optimize multiple correlated responses with model uncertainty. Besides, Ouyang et al. (2017) proposed a new loss function approach to consider the model parameter uncertainty and implementation errors. However, these loss functions only focus on the bias and robustness of the multivariate process, while neglecting the conformance probability of responses. In recent years, a posterior predictive function has attracted tremendous interest in tackling MRS optimization problems. The posterior predictive function can measure the conformance probability of future multiple responses that meet their specified quality conditions. Some quality experts also viewed the conformance probability as the reliability of an acceptable quality result (Peterson, 2004). Chiao and Hamada (2001) proposed a procedure to maximize the conformance probability that all responses simultaneously meet their respective specification limits. As noted by Peterson (2004), a significant drawback of their method is that it failed to take into account the uncertainty of model parameters. Peterson (2004) proposed a posterior predictive approach that considers the correlation among multiple responses, the variability of the process distribution, and the model parameter uncertainty as well as the reliability of optimization results. MiroQuesada, Del Castillo, and Peterson (2004) proposed a new Bayesian approach for MRS optimization problems in the presence of noise variables. Their approach not only includes noise variables in the integration of the predictive density but also performs the optimization for the numerical maximizations using nonlinear programming algorithms. As pointed out by Kazemzadeh, Bashiri, Atkinson, and Noorossana (2008), the posterior predictive function can help practitioners to control the conformance probability of the responses falling within their corresponding specification regions, but it fails to take into consideration the deviation from the targets. Peterson et al. (2009) proposed a Bayesian predictive approach to MRS optimization problems with seemingly unrelated regression models, which provides more flexible and accurate modeling than the standard multivariate regression (SMR) models with the same covariate structure across different responses. Besides, Del Castillo, Colosimo, and Alshraideh (2012) presented a Bayesian predictive modeling and optimization approach for functional response systems with a hierarchical two-stage mixedeffects model. Robinson, Pintar, Anderson-Cook, and Hamada (2012) proposed a Bayesian predictive approach to process optimization for a general class of robust parameter design models, including both normal and non-normal responses, in the split-plot context using an empirical approximation of the posterior distribution for an objective function of interest. That above posterior predictive function performs well since it provides a more feasible solution by measuring the conformance probability of all future responses falling within their corresponding specification limits. However, existing posterior predictive function approaches may pay too much attention to the conformance probability rather than the expected loss of product or process. The conformance probability can reflect some information on the production qualification rate required by customers or experimenters. However, the purpose of continuous quality improvement is not only higher conformance probability (i.e., reducing the defects of product or process) but also lower expected quality loss (i.e., decreasing deviations and improving robustness). Wang et al. (2016) proposed a new Bayesian approach to MRS optimization problems through integrating the loss function with the posterior probability function using a standard multivariate regression (SMR). As pointed out by Wang et al. (2016), it may be unreasonable to assume the same model forms for different responses in some cases. Motivated by the present literature, we propose a new method for 2

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

the MRS optimization problems in a unified framework of Bayesian seemingly unrelated regression (SUR) modeling and optimization. The proposed method attempts to make a trade-off between the conformance probability and the expected quality loss. The rest of this paper is organized as follows. Section 2 presents a Bayesian analysis of SUR models. The proposed method is proposed in Section 3. In Section 4, a polymer experiment from Box, Hunter, and Hunter (1978) is used to compare the proposed Bayesian SUR models with the Bayesian SMR models. Section 5 further illustrates the advantages of the proposed method through a simulation experiment with replicate response data. Finally, some conclusions and future research work are given in Section 6.

models. In recent years, some scholars have constructed a Bayesian inference method for SUR models. Bayesian prior distribution can be roughly divided into two categories: no information prior and conjugate prior, so Bayesian parameter estimation and statistical inference of the SUR model are often analyzed from the above two perspectives. The first category is the use of non-information priors as a parameter in the absence of knowledge of the parameter information. One of the most popular non-informative priors is Jeffreys’s invariant prior (Jeffrey, 1946). Assuming that parameters and are independent of each other, the prior joint distribution of parameters and is 1(

2.1. Overview of the SUR model The seemingly unrelated regression (SUR) models were originally proposed by Zellner (1962) as a modeling tool for response surfaces and have been widely used in many fields. The linear SUR models involve a set of regression equations with cross-equation parameter restrictions and correlated error terms with different variances (Zellner & Ando, 2010). The SUR model allows the experimenter to fit a different linear model for each response type, which will provide more flexible and accurate response surface modeling than one would obtain with the standard multivariate regression (SMR) models (Peterson et al., 2009). Here, Y = (Y1, Y2, , Ym) is a n × m matrix of m response-types and x is a k × 1 vector of factors. There are n observations for each response, and the SUR models are given as j

+

j

j = 1, 2,

,m

g1 ( , |D)

Cov ( i, j ) =

ij In

(1)

,m

i

1 (2 )nm 2 | |n

2

exp

1 tr{R 2

)

m+1 2

| |

(3)

(n + m + 1) 2 exp

| |

1 tr{R 2

1}

(4)

g1 ( | D, )

N( ,

)

(5)

g1 ( | D, )

IW (R , n)

(6)

= {ZT (

1

I) Z} 1ZT (

= (ZT (

1

I) Z)

1

(7)

I) Y

(8)

1

where represents the covariance matrix of (i.e., cov[ ]), IW (·,·) denotes the inverse Wishart distribution. R has the same form as in Eq. (2). It is worth noting that the conditional posterior density function of and depend upon each other. The second category uses an empirical Bayesian prior to model parameters based on past empirical knowledge (Zellner & Ando, 2010). The prior distribution of the parameters is assumed to have the same form as the posterior distribution, the prior joint distribution of parameters and is given as 2(

, )=

2(

)

2(

)

with

j

As shown in Eq. (1), the equations have different independent variables and variances. Also, the SUR models allow error terms in different equations to be correlated. In the matrix form, the SUR models in Eq. I) (1) are expressed as Y = Z + , Ñ(0, I) , where N (0 , denotes the normal distribution with mean 0 and covariance matrix I , is the tensor product, is the m × m matrix with the diagonal elements { 12, , m2 } , and the off-diagonal ijth elements are . The estimates of and in SUR models are obtained by maximizing the likelihood function.

f (D| , ) =

1(

with

jj In

i, j = 1, 2,

)

From the joint posterior density function given in Eq. (4), the conditional posteriors for and can be obtained as follows:

where Yj is an n × 1 vector of observations on the jth response, Zj (x) is a n × mj pre-determined matrix obtained by the specific model form for jth response, j is a mj × 1 vector of unknown parameters, and j is a random error vector associated with the jth response. Also, the error terms j in any individual model have constant variance, but error terms are correlated in different models (Shah, Montgomery, & Carlyle, 2004). Algebraically, this means that error terms j in Eq. (1) have the following properties:

E ( j ) = 0 Var ( j ) =

1(

which is proportional to the square root of the determinant of the Fisher information matrix. Using the likelihood function in Eq. (2) and the prior formation in Eq. (3), the joint posterior density function for parameters and is given as:

2. Bayesian analysis of SUR models

Yj = Zj (x)

, )=

1}

2(

) = N ( 0 , A - 1)

2(

) = IW (

0,

0)

According to the prior distributions of parameters and , it can be inferred that the conditional posterior distribution of parameters and is given as

g2 ( | D, )

N( ,

h2 ( | Y , X , )

IW (

) 0

+ R, k +

0)

with

(2)

where “tr” denotes the trace of a matrix. | | = det( ) denotes the value of the determinant of , the ijth element rij of m × m matrix R = (rij ) is rij = (Yi Zi (x) i )T (Yj Zj (x) j ) . D denotes the given experimental data. The maximum likelihood estimates of and can be obtained by using the iterative approach (Zellner, 1962).

= (X (

1

I) X + A ) 1 ( X (

= (X (

1

I) X + A )

1

I) X

+A

0)

1

where A denotes the variance-covariance matrix estimation for model parameters , and k = dim{ } . Currently, one of the most widely used Bayesian estimation methods for the SUR model is the use of the Markov chain Monte Carlo (MCMC) approach that is applied in the following Bayesian inference.

2.2. Prior and posterior analysis

2.3. Bayesian inference using Gibbs samplers

Since the Bayesian method can adequately consider a priori information or expert knowledge in the modeling process, it is widely used in parameter estimation and statistical inference of complex

In order to generate new response values Ynew at any given 3

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

Y (x) is given, we can define the estimated model as follows:

parameter setting x , we will simulate a large number of response values from the posterior distributions of and in the above SUR models via the Gibbs sampling procedure. First, the response surface model is fitted to each response in combination with the experimental data to determine the primary form of each response model, and the initial estimations (0) and (0) for each response are obtained by fitting the generalized linear square (GLS) method. Then, the coefficient vectors (k ) are updated by drawing a new value from the marginal posterior density ( | D, (k 1) ) in Eq. (5). Also, the variance-covariance matrices (k ) are updated by drawing a new value from the marginal posterior density ( | Y , Z, (k ) ) in Eq. (6). Finally, we repeat the above iterative sampling steps for a large number of times (k = 1, 2, , N ). The detailed Gibbs sampling algorithm (Ando, 2010) is summarized as follows:

Ynew (x) = Y (x) +

where new (x) is a random error term. With Eqs. (10) and (11), the expected loss function is given by

update X),1)

update end

(0)

and

+ trace [C

in Eq. (5) with the matlab code mvnrnd(inv(XT*X)*XT*Y, inv(XT*O*-

(k )

in Eq. (6) with the matlab code reshape(iwishrnd(R,n), 1,2*neq)

The Gibbs sampling produces samples by sequentially updating each model parameter conditional on the current values of the other parameters. Using the Gibbs sampling procedure as mentioned above, the parameter values of and can be estimated by a large number of posterior samples of (k ) and (k ) . It is worth noting that the initial samples (i.e., burn-in terms) should be discarded as being unrepresentative of the posterior distribution. The remaining samples are then employed for the posterior inference. Using the Gibbs samples of (k ) and (k ) , the simulated responses Y (k ) (x) based on the SUR models new can be generated as follows: (k )

(k ) Y new (x) = Y new +

(k )

=Z

(k )

+

(k )

k = B + 1, B + 2,

,N

(9)

ij

) y (x)]

(12)

(13)

=

T

i

j

n

i , j = 1, 2,

,m

where is the residual vector from the OLS estimation based on the posterior mean of (k ) . Furthermore, the prediction variance-covariance matrix y (x) (i.e., y (x) = [ yij (x)]) can be calculated as yij (x )

= zTi cov[ i ,

where cov[ i , for

j]

j ] zj

i, j = 1, 2,

,m

represents the block matrix of the covariance matrix

which can be calculated using Eq. (8). In another case with replicated response data, the samples of (k ) and (k ) can be obtained by using Gibbs sampling based on the sample mean Y of replicated response data instead of the responses Y in Eqs. (7) and (8). The mean vector E [Y (x)] also can be calculated using the mean of sampling values Z (k ) , which are the same as in the case with unreplicated response data. In this case, the variance-covariance matrix y (x) can be calculated using the variance-covariance matrix of residual terms, i.e., cov(Ei (x), Ej (x)) , where Ei = (ei1 (x), , eir (x)) represents the residual vector between the response values Yir and the predictive values Yir , and r is the replicated times of response Yi . It is noted that the predictive values Yir can be calculated by

From the perspectives of quality and reliability, a natural way to optimize a multi-response process is to integrate the loss function and the posterior predictive function in a single framework of Bayesian SUR modeling and optimization. In doing so, we can address a series of critical issues in a multi-response optimization process simultaneously. 3.1. Bayesian analysis of quality loss function Ko, Kim, and Jun (2005) proposed an improved quality loss function by integrating the strengths of Pignatiello (1993) and Vining’s methods (1998). The improved quality loss function was defined as

)

(x)] + trace [C

where EL (x) denotes E [L (Ynew (x), )]. The detailed proof of Eq. (13) can be found in Ko et al. (2005). In the framework of Bayesian SUR modeling, the loss function proposed by Ko et al. (2005) can be analyzed to consider the uncertainty of model parameters under two different cases, including unreplicated response data and replicated response data. In the case with unreplicated response data, the mean vector E [Y (x)] can be calculated by using the posterior samples of (k ) , namely, the mean of sampling values Zn (k ) . The matrix y (x) (i.e., y (x) = [ ij]) can be obtained from the OLS estimation of the individual responses, i.e.,

3. Multi-Response surface optimization using Bayesian SUR models

L (Ynew (x), ) = (Ynew (x)

y

EL (x) = ELbias + ELpred + ELrobust

I), and B denotes the number of where (k ) is sampled from N (0 , (k) initial samples discarded as burn-in terms. Simulation-based Bayesian inference requires using simulated draws to calculate any relevant quantities of interest. However, we need to treat the simulation draws with caution. It is necessary to judge whether the simulation draws has reached its stationary, or the desired posterior distribution through convergence diagnostics before proceeding to make any inference. Furthermore, it is necessary to check the convergence of all parameters, not just those of particular interest. See (Cowles & Carlin, 1996) for discussions about convergence diagnostics.

)T C (Ynew (x)

) C (E [Y (x)]

where y (x) represents the variance-covariance matrix of the predictive values and y (x) is the variance-covariance matrix of the responses. The expected loss function consists of three terms in its right-hand side. The T ) C (E [Y (x)] ) represents the bias due to the first term (E [Y (x)] deviation from the target. The second term, trace [C y (x)], evaluates the quality of predictions due to the uncertainty of predicted responses, which also reflects the variability of prediction responses because of model parameter uncertainty in the simulated modeling process. The third term, trace [C y (x)], measures the robustness due to the intrinsic variation of responses, which also reflect the variability of simulated responses because of the uncertainty of the model itself in the simulated modeling process. The three terms in Eq. (12) represent the bias (ELbias ), the quality of predictions (ELpred ), and the robustness (ELrobust ), respectively. Consequently, the expected quality loss can be rewritten as

(0)

(k )

T

E [L (Ynew (x), )] = (E [Y (x)]

Gibbs sampling algorithm Obtain the initial for j = 1 to N do

(11)

new (x)

Yir = Zr (x)(Z

Z) 1Z

Yir

Therefore, the variance-covariance matrix

(10)

log trace[C

where Ynew (x) is a m × 1 vector of responses at a new observation at x , is the m × 1 vector of the specified target values, and C is a m × m positive definite cost matrix which represents the losses incurred when Y (x) deviates from the target . Assuming the response prediction value

y (x)]

= log trace[C

y (x)

is calculated by

cov(Ei (x), Ej (x))]

(14)

Consequently, the regression model can be estimated based on the transformed response data log trace[C y (x)] and controllable factors x . On the other side, the predictive variance matrix yµ (x) (i.e., yµ (x) = [ yµij (x)]) can be calculated by 4

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al. yµij (x)

= zTi cov[

µi ,

µj ] zj

i, j = 1, 2,

,m

s. t P (Ynew (x )

(15)

i

(k ) .

On the other hand, cov[

presents the block matrix of the covariance matrix for

µ

µi ,

µj ]

re-

which can also

be calculated using Eq. (8). Here, it is worth noting that (i.e., (x) ) in Eq. (8) should be estimated by the residuals which are obtained from the sample mean of individual responses Yµ and the mean of predictive values for individual responses Yµ . Compared with Ko et al. ’s loss function, our proposed loss function extends their approach to take into account the uncertainty of model parameters to improve the robustness of process parameter fluctuations in the Bayesian SUR modeling framework. 3.2. Conformance probability using Bayesian SUR models As shown in the literature, the posterior predictive function has attracted some attentions to address the MRS optimization problems. The posterior predictive function can help experimental analysts to measure the conformance probability of future multiple responses that meet their specified quality conditions, which also reflects the reliability of an acceptable quality of products or processes. A Bayesian reliability approach initially proposed by Peterson (2004) to compute the conformance probability that a future multivariate response Ynew (x) will satisfy some specified quality conditions A for given experimental data in the framework of the SMR models, i.e., P (Ynew (x) A |data) . Here, the Bayesian reliability approach is extended to compute the conformance probability that a new simulated response values Ynew (x) obtained in Eq. (9) from the Bayesian SUR models. Hence, the conformance probability for given experimental data A can be calculated by the following equation:

P (Ynew (x)

A |data) =

A

P (Ynew (x)|data) d Ynew

(16)

A |data)

1 nsim

Since the proposed model in the Bayesian SUR modeling framework is highly nonlinear and heavily constrained, conventional optimization algorithms (e.g., gradient-based algorithms) and direct search methods (e.g., the Nelder-Mead simplex) can find only a local optimum or even fail to find a feasible solution (Ortiz, Simpson, Pignatiello, & HerediaLangner, 2004). In these situations, one alternative approach is to use a hybrid heuristic search procedure that can effectively integrate the advantages of global optimization algorithms (e.g., genetic algorithm (GA), or simulated annealing) with local search techniques (e.g., pattern search). Also, Jourdan, Basseur, and Talbi (2009) demonstrate that the hybrid optimization approach often performs better than one algorithm alone. Therefore, we utilize a hybrid genetic algorithm with a pattern search used by He et al. (2012) to optimize our proposed model. In the framework of Bayesian SUR modeling, the procedures of our proposed approach to multi-response optimization can be summarized as follows:

nsim k) I (Y (new

A)

s=1

(18)

3.4. Model optimization using hybrid genetic algorithm

Since the calculation of Eq. (16) with numerical integration is an intractable issue, one can approximate the conformance probability using a Gibbs sampling method as follows:

P (Ynew

p0

where p0 is the expected probability which the analysts or customers wish to be achieved. If the prior information for the expected probability p0 in practical applications lacks, an effective reference value can be obtained by maximizing the conformance probability in Eq. (17). The reference value will provide some information to help analysts to set a suitable expected probability. Assuming the specified quality conditions A to be product specification limits in the proposed approach, the constraint function in Eq. (18) focuses more on the conformance probability of the future responses or simulated responses falling within the corresponding specification limits. In this case, the expected probability p0 in Eq. (18) provides some information of the product qualification rate required by customers and experimenters. In doing so, the constraint function of the proposed approach can directly reflect the requirements of customers or experimenters, which also conforms to the basic definition of quality, namely, the characteristics of a product or service that bear on its ability to satisfy stated or implied needs. This also reflects Philip Crosby’s definition of quality “conformance to requirements.” On the other hand, product quality is also measured by the degree of conformance to predetermined specification limits or standards, and deviations from these standards (i.e., bias and robustness) can lead to poor quality and low reliability. Furthermore, the purpose of quality improvement is aimed at eliminating defects (i.e., products that are out of conformance), reducing deviation and improving process robustness to reduce overall product cost. The expected loss function in Eq. (18) is regarded as the objective function of the proposed approach to account for higher quality requirement such as process bias and robustness, which reflects another definition of product quality, namely, “the loss a product imposes on society after it is shipped” given by Genichi Taguchi from a more comprehensive view of the production system. Hence, it is crucial to take into account the conformance probability and quality loss in a single framework of Bayesian SUR modeling and optimization.

where yµi represents the predictive value of the mean for ith replicated response variables, namely, Zi (x) µi , and µi is the mean of the ith

parameter sampling values

A data)

(17)

where nsim is the number of simulated samples, and I ( ) is a 0–1 binary (k ) indicator function. If the sampling response values Y new from Bayesian SUR model falling within the region A , the binary indicator is equal to 1, and vice versa. The Bayesian posterior predictive approach not only allows analysts to take into account the correlation structure of the data, the variability of the process distribution, and the model parameter uncertainty, but also it can assess the conformance probability of future responses that meet the predetermined quality conditions. In addition, the Bayesian SUR models provide more flexible and accurate response surface models than the SMR models with the same covariate structure across different responses. Also, the proposed method allows an investigator to perform a comparative analysis between SUR models and SMR models to judge how different model forms for individual responses can modify the reliability of optimization results.

(1). Select a suitable model form for each response by analyzing the experimental data with the ordinary least square (OLS) regression method and the effect heredity principle (Wu & Hamada, 2000). (2). Build the response surface models between the controllable factors and responses with the Bayesian SUR modeling method. (3). Simulate the posterior samples of model parameters and responses using a Gibbs sampling method from Bayesian SUR models (4). Check the convergence of all parameters through employing the tools of convergence diagnostics before proceeding to make any inference. (5). Construct a quality loss function according to two different cases (i.e., un-replicate or replicate response data) by using posterior

3.3. Integrating quality loss with conformance probability Given that the advantages and drawbacks of the quality loss function and the posterior predictive function described above, a new multiresponse optimization model integrating the quality loss function with the posterior predictive function is constructed as follows:

min E [L (Ynew (x), )] = ELbias + ELpred + ELrobust 5

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

samples of model parameters in step 2, and consider it as the objective function of the proposed approach. (6). Establish a posterior predictive function through posterior samples of model responses to evaluate the reliability of predicted responses that satisfy the specified region. (7). Maximize conformance probability by using the posterior predictive function individually to obtain a reference value, and then view the inequation (when the conformance probability is more than some reference value) as the constrained function of the proposed approach. (8). Utilize the hybrid genetic algorithm with a pattern search method to optimize the proposed approach and find the most desirable setting of the controllable factors.

y2 = 59.95 + 3.59x1 + 2.23x3 + 0.82x12 The prediction models given above are similar to those models obtained by the SUR method in Ko et al. (Ko et al., 2005) except for adding a term x12 in the predicted model of response y2 . Therefore, the covariate vectors (e.g., model forms) in Eq. (1) for y1 and y2 can be assumed as follows:

z1 (x) = (1, x1, x2 , x3, x12 , x 22 , x32 , x1 x2 , x1 x3, x2 x3),

z2 (x) = (1, x1, x3, x12). 4.1. Bayesian inference and convergence diagnostics Here, we use the Gibbs sampling algorithm to obtain the posterior samples of parameters (e.g., and ) and responses in the Bayesian SUR modeling framework. First, we run additional 16,000 iterations to obtain nice convergence samples after discarding initial 4000 iterations as a burn-in term. For the polymer experimental data, summaries of the posterior draws for model parameters (i.e., and ) are given in Table 2.Note that the SUR model parameters for responses y1 and y2 are notated as 1 ~ 10 and 11 ~ 14 in Table 2, respectively. Moreover, each element of is set to ( 11, 12; 21, 22) in Table 2. Table 2 reports the posterior mean, standard deviations (SD), and 95% posterior credible intervals (CI). In addition, the potential scale reduction factor (PSRF) which is a convergence diagnostic test statistics of MCMC algorithm is also presented in Table 2. If the PSRF is close to 1, we can conclude that each of the model parameter draws has stabilized, and they are likely to reach the target distribution (Brooks & Gelman, 1998). Hence, the PSRF results in Table 2 indicate that posterior draws for each model parameter show good convergence behaviors. Moreover, some visual tools (e.g., trace, autocorrelation and density plots) can also be provided to help us to judge whether the Markov chains for all parameters in the SUR models have good convergence behaviors. Fig. 1 displays the diagnostic plots for the regression coefficient 1 in the polymer example. Also, the diagnostic plots for the parameter 12 (non-diagonal elements in ) are shown in Fig. 2. The diagnostic plots for the other parameters (not shown here) involved in this experiment have similar results. Besides, these chains in the trace plots traverse their posterior spaces rapidly, and they can jump from one remote region of the posterior to another in relatively few steps. Therefore, we can conclude that the mixing of these chains is quite good here. The autocorrelation plots in Figs. 1 and 2 have quick drop-offs, which suggest that the autocorrelation for these parameters is very weak and

The proposed approach is implemented by using MATLAB, and these codes are available under request from the authors of this paper. 4. A polymer experiment Box, Hunter, and Hunter (1978) present a polymer experiment, which is further analyzed by Vining (1998) and Ko et al. (2005). The purpose of the experiment is to find the optimal settings of controllable factors to maximize the conversion ( y1) of a polymer and achieve a target value of 57.5 for the thermal activity ( y2 ). The controllable factors chosen in the experiment are reaction time ( x1), reaction temperature ( x2 ), and amount of catalyst ( x3 ). A central composite design (CCD) is used by the analysts to perform the study. The experimental results are summarized in Table 1. The acceptable ranges for responses y1 and y2 are (80,100) and (55, 60), respectively. The response y1 is a larger-the-better (LTB-type) response. Hence, an appropriate target value for the conversion is set by 1 = 100 . The response y2 is a nominal-the-best (NTB-type) response, and hence a suitable target value for the thermal activity is set by 2 = 57.5. The cost matrix C , referring to Ko et al. (2005), is set by

C=

(0.100 0.025

0.025 0.500

)

The OLS regression method and effect heredity principle (Wu & Hamada, 2000) are utilized to fit response surface models for y1 and y2 , respectively. The predicted models for y1 and y2 are given as follows:

y1 = 81.09 + 1.03x1 + 4.04x2 + 6.20x3 + 2.95x 22

1.84x12

5.20x 32 + 2.13x1 x2 + 11.38x1 x3

3.88x2 x3

Table 1 The results of the polymer experiment. run

x1

x2

x3

y1

y2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

−1 1 −1 1 −1 1 −1 1 −1.68 1.68 0 0 0 0 0 0 0 0 0 0

−1 −1 1 1 −1 −1 1 1 0 0 −1.68 1.68 0 0 0 0 0 0 0 0

−1 −1 −1 −1 1 1 1 1 0 0 0 0 −1.68 1.68 0 0 0 0 0 0

74 51 88 70 71 90 66 97 76 79 85 97 55 81 81 75 76 83 80 91

53.2 62.9 53.4 62.6 57.3 67.9 59.8 67.8 59.1 65.9 60.0 60.7 57.4 63.2 59.2 60.4 59.1 60.6 60.8 58.9

Table 2 Summaries of posterior draws for model parameters. Parameter 1 2

3 4

5 6

7 8

9

10 11 12

13 14

11

12 22

6

Mean

SD

95%CI

81.076

1.346

81.047

81.106

0.999

4.033

0.892

4.014

4.053

1.000

1.034 6.211

−1.839 2.962

−5.202 2.122

11.382

−3.868 59.950 3.587 2.228 0.824

19.623 −0.811 2.390

0.904 0.901 0.883 0.872 0.871 1.181 1.190 1.169 0.377 0.352 0.354 0.345 7.943 1.884 0.933

1.014 6.192

−1.859 2.942

−5.222 2.096

11.356

−3.893 59.942 3.579 2.220 0.816

19.450 −0.853 2.369

PSRF 1.054 6.231

−1.820 2.981

−5.183 2.148

11.409

−3.842 59.959 3.595 2.235 0.832

19.797 −0.770 2.410

0.999 1.000 0.999 1.000 1.000 1.000 1.000 1.001 0.999 0.999 1.001 1.000 0.999 0.999 1.000

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

Diagnostics for the parameterβ 1

90

β1

85

80

75 0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1

0.5

0

-0.5

4

x 10 0.4

Posterior Density

ACF for the parameter β 1

iteration

0

5

10

15

0.3 0.2 0.1 0 74

20

76

78

80

82

Fig. 1. Diagnostic plots for the regression coefficient

can be ignored in the following inference. The posterior density plot in Fig. 1 is very smooth, and it approximates to a normal density curve. This result indicates that the prior assumptions for SUR model parameters seem reasonable. In summary, the convergence diagnostics for these parameters show that the draws from the posterior distributions of the model parameters can be used for prediction and assessment of various performance criteria in the following analysis.

84

β1

Lag

86

88

90

1.

4.2. Model optimization and comparative analysis Given that there is no prior information for p0 in Eq. (18), we can obtain a sufficient reference value by maximizing the conformance probability in Eq. (17) with the hybrid genetic algorithm (GA) with pattern search. The population size is set at 100, and the pattern search is selected as the hybrid function as well as other parameters are selected as default values in Matlab GA Optimization Toolbox. Assuming

Diagnostics for the parameterσ 10

12

5

σ 12

0 -5 -10 -15 -20 0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

iteration

x 10 0.4

0.3 Posterior Density

ACF for the parameter σ 12

1

4

0.5

0

-0.5

0

5

10

15

0.2

0.1

0 -25

20

-20

-15

Lag Fig. 2. Diagnostic plots for the regression coefficient 7

12 .

-10

σ

-5 12

0

5

10

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

Table 3 Optimization results for different methods. Different Methods

Posterior Probability Loss Function Ko et al’s Method Pignatiello’s Method Proposed Method

Setting

Components of Expected Loss

Marginal Probability

x1

x2

x3

ELbias

ELpred

ELrobust

EL

P1

P2

P

0.228 −0.313 −0.420 −0.450 −0.388

1.067 1.675 1.680 1.680 1.568

0.180 −0.422 −0.470 −0.510 −0.005

14.986 2.182 2.368 2.344 5.605

0.312 0.847 0.910 0.942 0.653

1.935 1.935 1.935 1.935 1.935

17.233 4.964 5.213 5.221 8.193

0.966 0.828 0.821 0.821 0.9128

0.997 0.967 0.940 0.930 0.9866

0.963 0.803 0.776 0.768 0.9030

that the fitness function is 1 P (Ynew A |data) , the conformance probability in which the responses satisfy the specific region A = {80 y1s 100, 55 y2s 60} can be optimized by using the hybrid genetic algorithm with pattern search. The values in bold in Table 3 indicate that the optimization criteria (e.g., the expected loss or the conformance probability) are utilized by the corresponding methods. The optimization results based on the posterior predictive function are given in the second row of Table 3. In this case, the optimal setting is (0.228, 1.067, 0.180) and the conformance probability (i.e., P ) is equal to 0.963. In the meantime, the expected quality loss (i.e., EL ) of the corresponding optimal setting is 17.233 based on Eq. (13). Furthermore, the quality loss function can be optimized by using Eq. (13) in the Bayesian SUR modeling framework without considering the constraint of conformance probability. Assuming that the maximum number of iterations GA performs is 100, the iterative optimization process using hybrid GA is displayed in Fig. 3. The optimization results based on quality loss function are shown in the third row of Table 3. In this case, the optimal setting is (−0.313, 1.675, −0.422), and the expected quality loss is equal to 4.964, where the conformance probability of the corresponding optimal setting is 0.803 based on Eq. (17). Referring to Ko et al. (2005), the expected quality loss and the Conformance probability can be calculated, respectively. The optimization results based on Ko et al.’s method are given in the fourth row of Table 3. In this case, the expected quality loss is 5.213, and the conformance probability is 0.776. Moreover, the optimization results based on Pignatiello’s method are also displayed in

the fifth row of Table 3, which are close to the results of Ko et al. (2005). Given the optimization results given above, we assume that the expected conformance probability (i.e., p0 ) which the analysts or customers are satisfied with is equal to 0.9 in our proposed method. Given the same assumption in the optimization process, the proposed model, namely, Eq. (18) is optimized by using the hybrid GA. The optimization results based on the proposed method are shown in the final row of Table 3. In this case, the optimal setting is (−0.388, 1.568, −0.005), the quality loss is 8.193, and the conformance probability is 0.903. The component ELrobust is the same in all five methods because of the equal robustness, i.e., y (x) = in the case with un-replicate response data. Also, It is noted that the marginal conformance probability (P1) of response y1 is slightly less than the marginal conformance probability (P2 ) of response y2 in Table 3. The results show that the conformance probability can be affected by model predictive ability because the prediction R2 of response y1 is less than that of response y2 . Compared to other methods in Table 3, the posterior probability method has the highest conformance probability and the largest expected quality loss. This result implies that this method does pay more attention to the reliability of optimization results but often ignores the quality loss. On the contrary, the Bayesian quality loss approach in Table 3 gives an optimal setting at which the quality loss is the best while the conformance probability is undesirable. The optimization result of the Bayesian loss function is similar to the result of Ko et al. (2005), which reflects the effectiveness of Bayesian SUR modeling and

Best: 4.96394 Mean: 5.03076

40

Best fitness Mean fitness

35

30

Fitness value

25

20

15

10

5

0

0

10

20

30

40

50

60

70

Generation Fig. 3. Optimization process using the Bayesian loss function. 8

80

90

100

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

optimization techniques for the multi-response problems. Compared with the posterior predictive function, the proposed method not only dramatically reduces the expected quality loss but also yields a very close conformance probability in the polymer experiment. On the other hand, compared with the Bayesian quality loss function or Ko et al. ’s method, the proposed method can significantly improve the conformance probability. In summary, the proposed method gives more consistent results than the existing approaches (e.g., posterior predictive function or quality loss function) when both conformance probability and quality loss are essential issues in dealing with multiresponse optimization problems.

provide more flexible and accurate modeling than the Bayesian SMR models do. The results also suggest that the predictive ability of the response surface models have a significant impact on the optimization results from the SMR and SUR models as presented in this experiment. 5. A simulation experiment 5.1. Simulated experimental data The polymer experiment gave above only involves un-replicate response data, but it is necessary to analyze this experiment with replicate response data using the proposed method in the paper. Referring to Ko et al. (2005), we generate the simulation data with five replicates using the following model:

4.3. Comparative analysis of Bayesian SMR and SUR modeling As noted before, the standard multivariate regression (SMR) model was used by Peterson (2004) to analyze multi-response optimization (k ) (x) (i.e., posterior problems. In this case, it is easy to simulate Y new draws of responses) from the corresponding predictive density function since the associated predictive density functions based on the SMR models have the multivariate t -distribution form. However, as noted by Peterson et al. (2009), a practical drawback of the methodology in Peterson (2004) is that the regression model is limited to the SMR models with the same covariate structure across response types. Hence, the covariate vectors in the SMR models for responses y1 and y2 should be assumed as:

z1 (x) = z2 (x) = (1,

x1, x2 , x3 , x12 ,

x 22 ,

x 32 ,

y1i (x) y (x) = 1 + y2i (x) y2 (x)

iÑ

x1 x2 , x1 x3 , x2 x3)

x1

x2

Posterior Probability Loss Function Pignatiello’s Method Ko et al’s Method Proposed Method

−0.46

(00),

SUR Models

x3

EL

P

EL

P

1.15

−0.48

20.9581

0.6478

16.5875

0.8851

−0.29 −0.48

1.68 1.68

−0.41 −0.56

11.5613 11.8398

0.4612 0.4810

4.9623 5.4594

0.8013 0.7504

−0.38

1.68

−0.49

11.6662

0.4806

5.1536

0.7892

−0.43

1.44

−0.49

14.6636

0.6028

9.4352

0.8823

11 (x)

12 (x)

21 (x)

22 (x)

11 (x)

= exp(3

x12

22 (x)

= exp(2

2x12

12 (x)

=

21 (x)

2x32 )

= 0.03

x 32), 11 (x) 22 (x).

The simulation data are displayed in Table 5. Using the simulated experimental data, two response surface models for the mean of replicate responses are estimated by using the OLS regression method as follows:

y1µ = 80.8 + 0.87x1 + 4.29x2 + 6.36x3

1.59x12

+ 3.32x 22 - 5.55x 32 + 2.04x1 x2 + 11.2x1 x3

4.01x2 x3

y2µ = 59.90 + 3.55x1 + 2.15x3 + 0.84x12 The discussion of Bayesian convergence diagnostics for model parameters in the polymer experiment is presented in Section 4.1, and no violations of assumptions are observed in this experiment. Given that the simulated experimental data are generated based on the polymer experiment, we do not display the results of the convergence diagnostics for model parameters in this simulation experiment. 5.2. Bayesian modeling and the comparative analysis Using Eq. (14) as described above, the regression model based on the response data log trace[C y (x)] and the controllable factors x can be estimated as follows:

log trace[C

SMR Models

,5

where

Table 4 Comparative analysis between Bayesian SMR and SUR models. Setting

i = 1, 2,

Here, y1 (x) and y2 (x) represent the original response data based on the polymer experiment in Table 1, and

Given that the SUR models with different model forms can provide more flexible and accurate modeling than the SMR models do, it is meaningful to investigate the impact of different modeling forms on the optimization results through comparing SMR with SUR models in the framework of Bayesian modeling and optimization. The results of a comparative analysis between Bayesian SMR and SUR models are displayed in Table 4. The values in bold given in Table 4 indicate that the optimization criteria (i.e., the expected quality loss EL or conformance probability P ) are considered by the corresponding methods in the Bayesian SMR modeling framework. Moreover, the optimal parameter settings in Table 4 are also obtained by their corresponding methods based on the Bayesian SUR models. Note that the parameter settings of two different loss function methods (i.e., Pignatiello’s method and Ko et al. ’s method) in Table 4 refer to the optimization results of the case (i.e., low quality of predictions and equal robustness) in Ko et al. (2005). Based on the same optimal parameter settings found by the corresponding methods with the Bayesian SMR modeling in Table 4, we can obtain the optimization results of the expected loss and the conformance probability using Bayesian SUR models. Compared with the optimization results based on the Bayesian SMR models, all methods with the Bayesian SUR models can significantly improve the corresponding posterior probabilities, and significantly reduce their expected losses at the same parameter settings. It is because the Bayesian SUR models can

Different Methods

i

y (x)]

= 1.36 + 0.29x1

0.24x3

1.31x12

0.65x 32

In addition, the predictive variance matrix yµ (x) is calculated by using Eq. (15). For the simulation experiment with replicate response data, the optimization process is performed based on the same cost matrix and target values given in the polymer experiment. The results are presented in Table 6. As noted before, the posterior predictive function can obtain the highest reliability of optimization results, while the quality loss at the same parameter settings is the worst. As far as quality loss is concerned, the quality loss function approaches (e.g., Bayesian loss function or Ko et al. ’s method) can find the desired parameter settings. However, the conformance probability for the same parameter settings may be unsatisfactory. Assuming that the reaction temperature (x2 ) is set at 1.501, a mesh plot of the conformance probability P versus the reaction time 9

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

Table 5 Simulated experimental data with replicate response data. x1

x2

x3

y11

y12

y13

y14

y15

y21

y22

y23

y24

y25

−1.00 1.00 −1.00 1.00 −1.00 1.00 −1.00 1.00 −1.68 1.68 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

−1.00 −1.00 1.00 1.00 −1.00 −1.00 1.00 1.00 0.00 0.00 −1.68 1.68 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

−1.00 −1.00 −1.00 −1.00 1.00 1.00 1.00 1.00 0.00 0.00 0.00 0.00 −1.68 1.68 0.00 0.00 0.00 0.00 0.00 0.00

73.6 51.4 88.7 70.1 70.7 89.8 66.0 96.9 77.5 79.2 85.1 97.0 50.3 82.7 85.6 65.4 79.9 84.7 82.8 90.2

73.9 51.7 88.7 70.5 71.8 89.5 66.2 96.9 76.0 79.1 81.5 102.8 54.6 80.6 79.7 76.0 81.5 90.0 80.8 90.3

73.7 50.7 88.4 70.7 71.0 90.3 66.6 96.8 76.1 78.9 84.9 96.9 49.3 80.6 83.7 67.5 77.0 80.1 74.5 82.8

74.0 50.9 89.3 69.7 72.2 90.7 67.7 96.5 76.1 79.1 87.0 97.3 55.1 81.6 84.0 74.1 77.6 84.2 74.4 89.1

74.5 51.2 88.8 70.7 71.5 88.9 65.5 96.9 75.9 78.1 85.0 97.1 53.5 75.9 76.7 75.3 72.0 83.6 81.5 96.9

54.6 63.8 52.6 61.4 57.9 67.7 61.4 67.4 59.1 65.8 60.6 62.4 57.2 63.0 57.2 58.3 55.9 60.7 62.2 61.8

53.2 62.4 54.0 63.7 57.0 67.6 60.7 67.5 57.4 66.8 59.9 60.7 55.3 62.6 64.2 60.2 67.7 61.1 64.5 52.6

53.6 63.1 54.1 61.7 56.9 67.1 60.5 68.4 59.2 65.9 62.4 59.0 56.7 63.5 59.0 59.6 66.0 61.8 58.7 61.4

52.5 62.1 53.0 63.0 57.3 67.7 60.2 68.2 59.0 65.3 60.0 60.6 61.4 60.1 57.9 54.2 62.1 59.0 55.4 61.8

52.7 63.2 53.0 62.6 57.6 67.9 60.0 68.0 59.1 65.8 57.1 54.7 57.0 62.5 55.5 59.7 60.2 62.9 57.4 60.8

( x1) and the amount of catalyst ( x3 ) is displayed in Fig. 4, which reflects the tendency of the expected conformance probability with respect to the controllable factors in this simulated experiment. Fig. 4 indicates that the further the parameter setting of the reaction time ( x1) and the amount of catalyst ( x3 ) deviates from the center of the specific region, the lower the conformance probability of responses will be achieved in this simulated experiment. Two mesh plots of a bias component (i.e., ELbias ) and a predictive quality component (i.e., ELpred ) versus the reaction time ( x1) and the amount of catalyst ( x3 ) are displayed in Fig. 5(a) and (b), respectively. Fig. 5(a) indicates that the bias component ELbias plays a dominant role when compared with other components in the total expected loss (i.e., EL ). Similarly, different methods in Table 6 provide different results of the predictive quality in this simulated experiment with replicated responses. The predictive component ELpred of the expected loss can be measured by the prediction variance-covariance matrix y (x) , which is relevant to the block matrix of the covariance matrix for the model parameter estimations and the variable vectors z(x ) in Eq. (15). Hence, the predictive component ELpred can reflect the variance of the predictive responses due to the model parameter uncertainty. It can be observed from Fig. 5(b) that the further the parameter setting of the reaction time (x1) and the amount of catalyst ( x3 ) closes to the center of the specific region, the smaller expected predictive loss ELpred will be achieved in this simulated experiment. Therefore, the predictive component ELpred of the posterior predictive function is the smallest when compared with other methods in Table 6 because their parameter setting for the reaction time ( x1) and the amount of catalyst ( x3 ) closes to the center of the specific region where the effect of model parameter uncertainty is small. This result shows that the movement of parameter setting can mitigate the adverse effect of model parameter uncertainty.

1 0.8

P

0.6 0.4 0.2 0 2

1

x3

0

-1

-2 -2

-1

0

2

1

x1

Fig. 4. Mesh plot of the conformance probability.

Unlike in the polymer experiment, different methods in Table 6 provide unequal robustness results, i.e., ELrobust . It is because the robustness component ELrobust of the expected loss can be measured by the variance-covariance matrix, y (x) which is relevant to the controllable factors x in this simulation experiment with replicate response data. A mesh plot of the robustness component ELrobust versus the reaction time ( x1) and the amount of catalyst (x3 ) is displayed in Fig. 5(c), which reflects the intrinsic variation of responses in the specified region. Fig. 5(c) indicates that the further the parameter setting of the reaction time ( x1) and the amount of catalyst (x3 ) deviates from the center, the better robustness will be achieved in this simulated experiment. The expected quality loss in Vining’s method (G G Vining, 1998) is minimized by using the sum of the bias component ELbias and the

Table 6 Optimization results for different methods. Different Methods

Posterior probability Loss Function Ko et al’s Method Vining’s method Pignatiello’s method Proposed Method

Setting

Components of Expected Loss

Marginal Probability

x1

x2

x3

ELbias

ELpred

ELrobust

EL

P1

P2

P

0.076 −0.681 −0.657 −0.400 −0.817 −0.372

1.201 1.680 1.680 1.679 1.679 1.501

−0.118 −0.536 −0.575 −0.386 −0.653 −0.325

11.846 1.532 1.484 1.390 2.254 4.040

0.390 1.288 1.296 1.007 1.526 0.720

4.344 2.171 2.205 3.369 1.560 3.525

16.580 4.991 4.998 5.766 5.340 8.285

0.961 0.821 0.816 0.831 0.799 0.922

0.996 0.856 0.855 0.966 0.743 0.980

0.958 0.707 0.705 0.804 0.604 0.905

Note that the values in bold are the best. 10

Computers & Industrial Engineering 142 (2020) 106357

6

5

5

4

ELrobust

400 350 300 250 200 150 100 50 0 2

ELpred

ELbias

J. Wang, et al.

4 3 2

1

0

x3

-1

-2 -2

-1

0

1

2

0 2 1

0

x3

x1

(a)

2 1

1 0 2

3

-1

-1

-2 -2

0

1

2

1

0

x3

x1

(b)

-1

-1

-2 -2

0

1

2

x1

(c)

Fig. 5. Mesh plot of different components for the expected loss.

predictive quality component ELpred . Hence, the results in Table 6 show that Vining’s method can obtain the smallest bias and the best predictive quality while robustness is the worst compared with other loss function methods (i.e., Bayesian loss function, Ko et al. ’s method, and Pignatiello’s method). Meanwhile, the conformance probability of the parameter settings found by Vining’s method is higher than other loss function methods. This result indicates that the quality of model predictions may have some influence on the results of the conformance probability. On the other hand, the expected quality loss in Pignatiello’s method is optimized by the sum of a bias component and a robustness component. Compared with other several loss function methods, the robustness component ELrobust is the smallest, while the predictive quality component ELpred is the worst in Pignatiello’s method. The result indicates that Pignatiello’s method pays more attention to the robustness but ignores the quality of predictions. Furthermore, the conformance probability of Pignatiello’s method for the corresponding parameter settings is the lowest in comparison with other loss function methods. It is worth noting that there is an inverse relationship between the component ELpred and the corresponding conformance probability, i.e., as the component ELpred increases, the corresponding conformance probability decreases, and vice versa. The result again shows that the conformance probability is closely related to the predictive ability of models. Generally speaking, utilizing either the conformance probability method or the loss function methods cannot find desirable parameter settings when quality loss (i.e., bias, quality of predictions and robustness) and conformance probability need to be considered simultaneously. Compared with the posterior predictive method, the proposed method significantly reduces the expected quality loss when the expected conformance probability p0 is 0.9, which satisfies the expectation of the analysts or customers. Moreover, compared with other loss function methods, the proposed method can greatly improve the reliability of optimization results. Just like the polymer experiment, the results in Table 6 demonstrate that the proposed method can find desirable parameter settings to make a trade-off between quality loss and conformance probability. It is very significant to display how to make the trade-offs between the two optimization criteria (i.e., conformance probability and quality loss) by using a standard Pareto front plot, which will help us to clarify the choices in Table 6. The expected loss function and the posterior probability function are simultaneously optimized by using the solver “Gamultiobj-Multiobjective Optimization using Genetic algorithm” in the MATLAB optimization tool. The Pareto front plot between the expected loss and conformance probability is displayed in Fig. 6.

University of Calgary, Canada. The laser machining center includes a femtosecond laser beam generator, a machining control computer, a micro-machining workstation, and measuring instruments (Wang, Ma, Tsung, Chang, & Tu, 2019). The whole experimental platform is shown in Fig. 7. The purpose of our experiment was to investigate a more accurate functional relationship between the input factors and the output responses during the laser micro-drilling process and find the optimal parameter setting. Three significant controllable factors are average power ( x1), pulse frequency (x2 ), and cutting speed (x3 ). Three key quality characteristics (diameter ( y1), roundness ( y2 ), and circularity ( y3 )) are chosen to reflect geometrical property and machining precision of micro-holes. Circularity tends to focus on whether the circular peripheral shape is closer to a standard circle. The formula of the circularity is defined as Circularity = 4 × Area [Perimeter ]2 . The formula of the circularity is related to the area and perimeter of micro-holes, which is different with the formula of the roundness. Roundness often reflects the sharpness or radius of curvature of the edges and the smoothness of the outer edge of the micro-hole. The formula of the roundness is shown as Roundness = 4 × Area Major axis 2 . The experimental data of Area , Perimeter , Major axis can be measured by analyzing the pictures of micro-holes with the software IMAGEJ. The photo software IMAGEJ and its Java source code are freely available in the following website: http://imagej.nih.gov/ij/. A CCD was used to study the micro-drilling process and find the optimal parameter settings. The experimental plan with actual values and observed experimental results are shown in Table 7. 6.2. Bayesian modeling and comparison analysis Here, the acceptable ranges for diameter, roundness and circularity are assumed to be [18, 24], [0.9, 1.0], and [0.6 1.0], respectively. Clearly, the response diameter ( y1) is nominal-the-best (NTB-type) response, while the response roundness ( y2 ) and circularity ( y3) are the larger-the-better (LTB-type) response. Hence, their target values are assumed to T1 = 21, T2 = 1, and T3 = 1, respectively. Moreover, referring to prior information, the cost matrix is assumed as

5 2.5 2.5 C = 2.5 3 2.5 2.5 2.5 3 According to the results of statistical modeling analysis, the predicted models for y1, y2 , and y3 are given as follows

1.642x2 + 2.047x12 + 1.855x 22 - 1.349x1 x2

6. Laser micro-drilling experiment

y1 = 12.778 + 3.212x1

6.1. Experimental platform and experimental data

y2 = 0.892

0.010x1

y3 = 0.673

0.098x1 + 0.083x2

A five-axis femtosecond laser machining center has been developed in the laboratory of laser micro/nano manufacturing from the

0.003x2

0.015x3 + 0.018x12 + 0.014x 22 + 0.014x 32

0.006x3 + 0.051x 32

Then, the model forms in Eq. (1) for y1, y2 , and y3 can be assumed as 11

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

Pareto front

18 16 14

EL

12 10 8 6 4 0.95

0.90

0.85

0.80

P

0.75

0.70

0.65

0.60

Fig. 6. Pareto front plot between the expected loss and conformance probability.

follows:

z1 (x) = (1, x1, x2 , x12 , x 22 , x1 x2) z2 (x) = (1, x1, x2 , x3 , x12 , x 22 , x 32) = (1,

It can be observed from Table 8 that the posterior probability method has the highest conformance probability and the largest expected loss. On the contrary, the expected loss method in Table 8 can obtain the smallest expected loss while the conformance probability is undesirable. Similarly, the results in Table 8 indicate that the posterior probability can significantly improve the conformance probability, but the loss function does pay more attention to the expected loss. Compared with the posterior probability method and the loss function method, the proposed method can make a trade-off between the expected loss and the conformance probability.

z3 (x)

x1, x2 , x3, x 32)

As shown in Section 4.1, Bayesian inference and convergence diagnostics for model parameters are analyzed by using the Gibbs sampling procedure. Likewise, the convergence diagnostics for the chains of all parameters in this experiment indicate that the posterior distributions of the model parameters in the SUR models can be further used for predictions and assessments of various performance indices (e.g., posterior probability and expected loss) in the following analysis. Here, we run additional 4000 iterations to obtain nice convergence samples after discarding initial 1000 iterations as a burn-in term. The population size is set at 100, and the maximum number of iterations is set at 100. Likewise, the expected conformance probability in Eq. (18) is assumed to be 0.9 in this experiment. The optimization results for several methods are displayed in Table 8.

6.3. Uncertainty analysis of optimization results It is a very insightful suggestion from an anonymous reviewer who inspires us to discuss how the Gibbs sampling method and the hybrid GA algorithm affect the performance criteria (i.e., expected loss and conformance probability) of our proposed method. Here, we carried out

Fig. 7. The experimental platform. 12

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

Table 7 Actual values and observed experimental results.

0.93

Responses

0.92

x2

x3

Diameter

Roundness

Circularity

50.0 150.0 50.0 150.0 50.0 150.0 50.0 150.0 15.9 184.1 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0

500.0 500.0 800.0 800.0 500.0 500.0 800.0 800.0 650.0 650.0 397.7 902.3 650.0 650.0 650.0 650.0 650.0 650.0 650.0 650.0

6.0 6.0 6.0 6.0 16.0 16.0 16.0 16.0 11.0 11.0 11.0 11.0 2.6 19.4 11.0 11.0 11.0 11.0 11.0 11.0

16.43 33.80 15.00 16.53 13.47 32.81 14.04 17.43 10.41 34.23 19.13 13.42 11.51 12.87 12.59 12.38 11.90 12.14 12.69 12.55

0.965 0.940 0.968 0.900 0.923 0.924 0.921 0.929 0.969 0.925 0.948 0.923 0.963 0.903 0.905 0.908 0.892 0.873 0.886 0.887

0.802 0.571 0.839 0.748 0.733 0.573 0.865 0.728 0.910 0.483 0.501 0.876 0.811 0.802 0.652 0.695 0.630 0.640 0.629 0.669

0.91

P

x1

0.9 0.89 0.88 0.87 0

100

200

300

400

500

600

Repeated runs

700

800

900

1000

(a) the posterior probability 28 27 26

EL

Controllable factors

0.94

25 24 23

1000 repeated runs at the optimal parameter settings obtained by the proposed method. In doing so, we can further analyze the influence of the uncertainty associated with the model parameter uncertainty and output response variability on the optimization results. The mean value of conformance probability is 0.8995 and its standard deviation is 0.0076, meanwhile the mean value of quality loss is 24.7162 and its standard deviation is 0.5761. The trace plots of 1000 repeated runs for two optimization criteria (i.e., posterior probability and expected loss) are show in Fig. 8.

22 21 0

100

200

300

400

500

600

700

800

900

1000

Repeated runs

(b) the expected loss Fig. 8. The results of 1000 times repeated runs for two optimization criteria.

and replicate response data. Results of the comparative analysis demonstrate that the Bayesian SUR models can provide more reliable and reasonable solutions than the Bayesian SMR models do. This research result also indicates that the quality of model predictions will have a significant impact on the optimization results. Furthermore, the impact of changed expected conformance probabilities on the optimization results of the proposed method is further discussed and analyzed by using the simulated experimental data. The results show that a suitable expected probability can effectively make a trade-off between quality loss and conformance probability. Also, the proposed approach takes into account the correlation structure of the data, the variability of the process distribution, and the uncertainty of the model parameters as well as the prediction performance of the process model. It is worth noting that the uncertainty of the model form itself is ignored in the proposed method (Apley & Kim, 2011). A natural extension is to consider the uncertainty of the model structure itself in the proposed method. In this aspect, the Bayesian model averaging (BMA) approach can be an alternative modeling method (Ng, 2010; Rajagopal & Del Castillo, 2005). Besides, more complex non-linear response surface models also are built to consider the case of robust parameter design for the proposed method through using Gaussian process models

7. Conclusions and further work The proposed method allows analysts to consider both conformance probability and quality loss in a single framework of Bayesian SUR modeling and optimization, which will significantly highlight the research contribution of this paper. Two different performance indices (i.e., conformance probability and quality loss) in multi-response optimization problems reflect two different implications in quality management. The posterior predictive function does pay more attention to the basic definition of product quality, namely, the collection of features and characteristics of a product that contributes to its ability to meet given requirements from the viewpoint of manufacturers. Nevertheless, the loss function usually focuses on reducing the quality loss incurred by the defective products with high deviation and low robustness, which reflects the philosophy of Taguchi’s quality loss from the viewpoint of customers. In such a case, the optimization results for multi-response optimization problems using a single performance index (i.e., quality loss or conformance probability) may be misleading when both quality loss and conformance probability all are important issues. Two examples are utilized to validate the effectiveness of the proposed method under two different cases, i.e., un-replicate response data Table 8 Optimization results for different methods. Different Methods

Posterior Probability Loss Function Proposed Method

Setting

Components of Expected Loss

Marginal Probability

x1

x2

x3

ELbias

ELpred

ELrobust

EL

P1

P2

P3

P

1.333 0.805 0.673

0.868 −0.137 −0.388

−1.380 −1.648 −1.603

4.537 0.187 4.420

8.065 2.195 3.131

16.794 16.794 16.794

29.396 19.806 24.349

0.914 0.853 0.914

0.996 0.997 0.998

0.981 0.988 0.991

0.916 0.844 0.906

Note that the optimization criteria are in bold. 13

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al.

(Tan & Ng, 2009; Tan & Wu, 2012). Apart from that, the proposed approach can be extended to consider non-normal multi-response optimization problems. It could be a useful tool for modeling non-normal responses using generalized linear models (GLM). A challenging work for further research is to simulate the posterior predictive distribution of the GLM through using the MCMC method (see Wang & Ma, 2013).

multiple response surface optimization in the presence of noise variables. Journal of Applied Statistics, 31(3), 251–270. Murphy, T. E., Tsui, K. L., & Allen, J. K. (2005). A review of robust design methods for multiple responses. Research in Engineering Design, 15(4), 201–215. Myers, R. H., Montgomery, D. C., & Anderson-Cook, C. M. (2009). Response surface methodology: Process and product optimization using designed experiments (3rd ed.). New York, NY: John Wiley & Sons Inc. Myers, R. H., Montgomery, D. C., Vining, G. G., Borror, C. M., & Kowalski, S. M. (2004). Response surface methodology: A retrospective and literature survey. Journal of Quality Technology, 36(1), 53–77. Myers, W. R., Brenneman, W. A., & Myers, R. H. (2005). A dual-response approach to robust parameter design for a generalized linear model. Journal of Quality Technology, 37(2), 130–138. Ng, S. H. (2010). A Bayesian model-averaging approach for multiple-response optimization. Journal of Quality Technology, 42(1), 52–68. Ortiz, F., Simpson, J. R., Pignatiello, J. J., & Heredia-Langner, A. (2004). A genetic algorithm approach to multiple-response optimization. Journal of Quality Technology, 36(4), 432–450. Ouyang, L., Ma, Y., & Byun, J. H. (2015). An integrative loss function approach to multiresponse optimization. Quality and Reliability Engineering International, 31(2), 193–204. Ouyang, L., Ma, Y., Byun, J. H., Wang, J., & Tu, Y. (2016). A prediction region-based approach to model uncertainty for multi-response optimization. Quality and Reliability Engineering International, 32(3), 783–794. Ouyang, L., Ma, Y., Wang, J., & Tu, Y. (2017). A new loss function for multi-response optimization with model parameter uncertainty and implementation errors. European Journal of Operational Research, 258(2), 552–563. Peterson, J., Miro-Quesada, G., & Del Castillo, E. (2009). A Bayesian reliability approach to multiple response optimization with seemingly unrelated regression models. Quality Technology and Quantitative Management, 6(4), 353–369. Peterson, J. J. (2004). A posterior predictive approach to multiple response surface optimization. Journal of Quality Technology, 36(2), 139–153. Pignatiello, J. J. (1993). Strategies for robust multiresponse quality engineering. IIE Transactions, 25(3), 5–15. Rajagopal, R., & Del Castillo, E. (2005). Model-robust process optimization using Bayesian model averaging. Technometrics, 47(2), 152–163. Robinson, T. J., Pintar, A. L., Anderson-Cook, C. M., & Hamada, M. S. (2012). A Bayesian Approach to the Analysis of Split-Plot Combined and Product Arrays and Optimization in Robust Parameter Design. Journal of Quality Technology, 44(4), 304–320. Shah, H. K., Montgomery, D. C., & Carlyle, W. M. (2004). Response surface modeling and optimization in multi-response experiments using seemingly unrelated regressions. Quality Engineering, 16(3), 387–397. Shin, S., Samanlioglu, F., Cho, B. R., & Wiecek, M. M. (2011). Computing trade-offs in robust design: Perspectives of the mean squared error. Computers & Industrial Engineering, 60(2), 248–255. Taguchi, G., Elsayed, E. A., & Hsiang, T. C. (1989). Quality engineering in production systems. New York: McGraw-Hill College. Tan, M. H., & Wu, C. J. (2012). Robust design optimization with quadratic loss derived from Gaussian process models. Technometrics, 54(1), 51–63. Tan, M. H. Y., & Ng, S. H. (2009). Estimation of the mean and variance response surfaces when the means and variances of the noise variables are unknown. IIE Transactions, 41(11), 942–956. Vining, G. G. (1998). A compromise approach to multiresponse optimization. Journal of Quality Technology, 30(4), 309–313. Vining, G. G., & Bohn, L. L. (1998). Response surfaces for the mean and variance using a nonparametric approach. Journal of Quality Technology, 30(3), 282–291. Wang, J. J., & Ma, Y. Z. (2013). Bayesian Analysis of Two-Level Fractional Factorial Experiments with Non-Normal Responses. Communications In Statistics-Simulation and Computation, 42(9), 1970–1988. Wang, J., Ma, Y., Ouyang, L., & Tu, Y. (2016). A new Bayesian approach to multi-response surface optimization integrating loss function with posterior probability. European Journal of Operational Research, 249(1), 231–237. Wang, J., Ma, Y., Tsung, F., Chang, G., & Tu, Y. (2019). Economic parameter design for ultra-fast laser micro-drilling process. International Journal of Production Research, 57(20), 6292–6314. Wu, C. F. J., & Hamada, M. (2000). Experiments: Planning, analysis, and parameter design optimization. New York: John Wiley & sons Inc. Wu, F.-C. (2005). Optimization of Correlated Multiple Quality Characteristics Using Desirability Function. Quality Engineering, 17(1), 119–126. Zellner, A. (1962). An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias. Journal of The American Statistical Association, 57(298), 348–368. Zellner, A., & Ando, T. (2010). A direct Monte Carlo approach for Bayesian analysis of the seemingly unrelated regression model. Journal of Econometrics, 159(1), 33–45. Zhang, L., Ma, Y., Ouyang, L., & Liu, J. (2016). Balancing the subjective and objective weights for correlated multiresponse optimization. Quality and Reliability Engineering International, 32(3), 817–835. Zhou, X. J., Ma, Y. Z., Tu, Y. L., & Feng, Y. (2013). Ensemble of Surrogates for Dual Response Surface Modeling in Robust Parameter Design. Quality and Reliability Engineering International, 29(2), 173–219.

Acknowledgements The authors are grateful to the editor and anonymous reviewers for their insightful comments and suggestions.This work is supported by the National Natural Science Foundation of China (Nos. 71771121, 71931066, and 71702072), the Fundamental Research Funds for the Central Universities (No.3091511102), and the Natural Science Foundation of Jiangsu Province (No. BK20170810) as well as NSERC (Natural Sciences and Engineering Research Council of Canada) Discovery Grant (No. RGPIN-2018-03862). Appendix A. Supplementary material Supplementary data to this article can be found online at https:// doi.org/10.1016/j.cie.2020.106357. References Ando, T. (2010). Bayesian model selection and statistical modeling. New York: Chapman & Hall. Apley, D. W., & Kim, J. (2011). A cautious approach to robust design with model parameter uncertainty. IIE Transactions, 43(7), 471–482. Box, G. E. P., Hunter, W. G., & Hunter, J. S. (1978). Statistics for experiments. In. New York, NY: John Wiley & Sons. Brooks, S. P., & Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7(4), 434–455. Chapman, J., Lu, L., & Anderson-Cook, C. (2014). Incorporating response variability and estimation uncertainty into Pareto front optimization. Computers & Industrial Engineering, 76, 253–267. Chiao, C. H., & Hamada, M. (2001). Analyzing experiments with correlated multiple responses. Journal of Quality Technology, 33(4), 451–465. Cowles, M. K., & Carlin, B. P. (1996). Markov chain Monte Carlo convergence diagnostics: A comparative review. Journal of The American Statistical Association, 91(434), 883–904. Del Castillo, E., Colosimo, B. M., & Alshraideh, H. (2012). Bayesian modeling and optimization of functional responses affected by noise factors. Journal of Quality Technology, 44(2), 117–135. Del Castillo, E., Montgomery, D. C., & McCarville, D. R. (1996). Modified desirability functions for multiple response optimization. Journal of Quality Technology, 28(3), 337. Derringer, G., & Suich, R. (1980). Simultaneous optimization of several response variables. Journal of Quality Technology, 12(4), 214–219. Goethals, P. L., & Cho, B. R. (2012). Extending the desirability function to account for variability measures in univariate and multivariate response experiments. Computers & Industrial Engineering, 62(2), 457–468. Harrington, E. C. (1965). The desirability function. Industrial Quality. Control, 21(10), 494–498. He, Y., He, Z., Lee, D.-H., Kim, K.-J., Zhang, L., & Yang, X. (2017). Robust fuzzy programming method for MRO problems considering location effect, dispersion effect and model uncertainty. Computers & Industrial Engineering, 105, 76–83. He, Z., Wang, J., Oh, J., & Park, S. H. (2010). Robust optimization for multiple responses using response surface methodology. Applied Stochastic Models In Business and Industry, 26(2), 157–171. He, Z., Zhu, P. F., & Park, S. H. (2012). A robust desirability function method for multiresponse surface optimization considering model uncertainty. European Journal Of Operational Research, 221(1), 241–247. Jeffrey, H. (1946). An invariant form for the prior pobability in etimation poblems. Proceedings Of The Royal Society of London Series A-Mathematical Physical A, 196, 453–461. Jeong, I. J., & Kim, K. J. (2009). An interactive desirability function method to multiresponse optimization. European Journal of Operational Research, 195(2), 412–426. Jourdan, L., Basseur, M., & Talbi, E.-G. (2009). Hybridizing exact methods and metaheuristics: A taxonomy. European Journal of Operational Research, 199(3), 620–629. Kazemzadeh, R., Bashiri, M., Atkinson, A., & Noorossana, R. (2008). A general framework for multiresponse optimization problems based on goal programming. European Journal Of Operational Research, 189(2), 421–429. Kim, K., & Lin, D. (2006). Optimization of multiple responses considering both location and dispersion effects. European Journal of Operational Research, 169(1), 133–145. Ko, Y. H., Kim, K. J., & Jun, C. H. (2005). A new loss function-based method for multiresponse optimization. Journal of Quality Technology, 37(1), 50–59. Miro-Quesada, G., Del Castillo, E., & Peterson, J. J. (2004). A Bayesian approach for

Jianjun Wang is an Associate Professor of the Department of Mangement Science and Engineering at Nanjing University of Science and Technology. He received his BSc in Applied Mathematics from Jishou University, China, and his MSc in Applied Statistics from Hunan University. He earned his Ph.D. in Quality Engineering from Nanjing

14

Computers & Industrial Engineering 142 (2020) 106357

J. Wang, et al. University of Science and Technology, China. He is a member of QSR and INFORMS, and a senior member of Chinese Society of Optimization, Overall Planning and Economical Mathemics. He has authored over 50 refereed journal pulications, and is a reviewer of several international journals such as JQT, EJOR, IJPR, CIE and QTQM. His research interests include parameter design and optimization, Bayesian statistics and modeling, industrial statistical and data analysis.

Linhan Ouyang is an assistant professor in the Department of Industrial Engineering, Nanjing University of Aeronautics and Astronautics, China. He earned his Ph.D. in Quality Engineering from Nanjing University of Science and Technology, China. His research interests include process modeling and design of experiments. Yiliu Tu is a professor at Department of Mechanical and Manufacturing Engineering, University of Calgary, Canada, and Zi Jin Scholar Chair Professor at School of Economics and Management, Nanjing University of Science and Technology, China. He received his BSc in electronic engineering and MSc in mechanical engineering, both from Huazhong University of Science and Technology (HUST), China; PhD in production engineering from Aalborg University (AU), Denmark. His present research interests are OKP (One-ofa-Kind Production) product design and manufacture, ultra-fast laser micro-machining technology, and modeling and simulation of networked critical infrastructure. He has published more than 170 research papers (SCI 78). He is a senior member of SME (Society of Manufacture Engineers) and a professional engineer of APEGGA (The Association of Professional Engineers, Geologists, and Geophysicists of Alberta).

Yizhong Ma is a professor in the Department of Management Science and Engineering at Nanjing University of Science and Technology. He received his BSc in Applied Mathematics from Huazhong Normal University, Wuhan, China, and his MSc in Quality Engeering and Ph.D. in Control Science from Northwestern Polytechnical University, China. He is an executive director of several associations such as Quality Society of China (QSC), Society of Management Science and Engineering of China (SMSEC) and Chinese Society of Optimization, Overall Planning and Economical Mathemics. His research interests include quality engineering and quality Management.

15