## 1 Introduction  
Among the many quality engineering methodologies, robust design (RD) based on statistical design and analysis methods and optimization methods have contributed significantly to the improvement of product/process quality for more than 20 years. The main objective of RD is to identify the optimal factor settings that can minimize both variability and bias (namely, the deviation between the process mean and desired target value) in a process/product. To demonstrate this RD principle, [PERSON] introduced a two-step model based on new design and analysis approaches, particularly orthogonal array and signal-to-noise ratios [1]. However, orthogonal arrays, statistical analysis, and signal-to-noise ratios are highly controversial [2, 3, 4, 5]. Consequently, [PERSON] and [PERSON] [6] improved [PERSON]'s model by proposing dual-response (DR) estimation and optimization approaches, in which the separate functions of the mean and variance responses are estimated using response surface methodology (RSM) based on the least-squares method (LSM). Their proposed optimization approach can calculate robust factor settings by minimizing process variability while maintaining the process mean at the target value. This results in the primary procedure of RD being identified in three sequential steps: design of experiment (DoE), response function estimation, and optimization. However, [PERSON] and [PERSON] [7] pointed out that process bias and variability are not considered simultaneously in the DR model and proposed a mean square error (MSE) model as an alternative approach providing more model flexibility when compared with that provided by the DR approach. These DR and MSE models were further extended by [PERSON] and [PERSON] [8], [PERSON] and [PERSON] [9], [PERSON] et al. [10], [PERSON] et al. [11], [PERSON] [12], [PERSON] and [PERSON] [13], [PERSON] and [PERSON] [14], [PERSON] et al. [15], [PERSON] and [PERSON] [16; 17; 18], [PERSON] et al. [19], [PERSON] [20], [PERSON] and [PERSON] [21], [PERSON] and [PERSON] [22; 23], [PERSON] et al. [24], [PERSON] et al. [25], and [PERSON] et al. [26]. In the response function estimation step of the RD procedure, RSM is used extensively to obtain the estimated function associated with an output response and its associated input factors based on the DoE results. RSM is applied extensively in real-world industrial situations, mainly where several input factors potentially affect a product or process's performance measures or quality characteristics [27]. The relationship between several input factors and one or more associated output responses can be estimated using RSM. Conventionally, LSM can estimate the unknown coefficients of a regression model with many specific error assumptions, such as the errors in the output responses should be followed by a normal distribution (i.e., zero mean and constant variance) identically and independently. However, these assumptions are often violated in many practical problems. Several alternative approaches, such as transformation, weighted least squares (WLS), maximum-likelihood estimation (MLE), Bayesian, and inverse problems, can be used to address this problem. For example, the WLS method has been applied to RSM to estimate the model parameters for unbalanced data [28; 29], and [PERSON] and [PERSON] [30] and [PERSON] et al. [31] have integrated an MLE method and an expectation-maximization algorithm for unknown parameter estimations using incomplete data. In addition, [PERSON] et al. [32] and [PERSON] and [PERSON] [33; 34] proposed a Bayesian approach to estimate the coefficients of response functions where the variances follow a log-normal distribution. [PERSON] and [PERSON] [22; 23] proposed a new estimation method that used an inverse problem approach based on Bayesian perspectives. Most of the existing estimation methods in the RSM and RD literature attempt to generate functional forms in input factors for output responses (namely, the process mean, process variance, and quality characteristic).