Penalty Function Optimization in Dual Response Surfaces Based on Decision Maker’s Preference and Its Application to Real Data

: The dual response surface methodology is a widely used technique in industrial engineering for simultaneously optimizing both the process mean and process standard deviation functions of the response variables. Many optimization techniques have been proposed to optimize the two ﬁtted response surface functions that include the penalty function method (PM). The PM method has been shown to be more efﬁcient than some existing methods. However, the drawback of the PM method is that it does not have a speciﬁc rule for determining the penalty constant; thus, in practice, practitioners will ﬁnd this method difﬁcult since it depends on subjective judgments. Moreover, in most dual response optimization methods, the sample mean and sample standard deviation of the response often use non-outlier-resistant estimators. The ordinary least squares (OLS) method is also usually used to estimate the parameters of the process mean and process standard deviation functions. Nevertheless, not many statistics practitioners are aware that the OLS procedure and the classical sample mean and sample standard deviation are easily inﬂuenced by the presence of outliers. Alternatively, instead of using those classical methods, we propose using a high breakdown and highly efﬁcient robust MM-mean, robust MM-standard deviation, and robust MM regression estimators to overcome these shortcomings. We also propose a new optimization technique that incorporates a systematic method to determine the penalty constant. We call this method the penalty function method based on the decision maker’s (DM) preference structure in obtaining the penalty constant, denoted as PMDM. The performance of our proposed method is investigated by a Monte Carlo simulation study and real examples that employ symmetrical factorial design of experiments (DOE). The results signify that our proposed PMDM method is the most efﬁcient method compared to the other commonly used methods in this study.


Introduction
Response surface methodology (RSM) was first developed by Box and Draper [1].It is an important tool to find the optimal factor settings of the design point, which can either maximize or minimize the given response function.RSM is widely used in many disciplines, such as in manufacturing industries, engineering, and agricultural sciences.The food industry, in particular, has been a prime user of RSM since the early 1970s [2].For example, Dey and Dora [3] studied the effects of temperature, pH, enzyme concentration/substrate concentration (E/S) ratio on the response, i.e., degree of hydrolysis (DH) for marine shrimp.They noted that RSM was successfully applied to determine the optimum operating conditions on the control variables for maximum DH value.The traditional RSM Symmetry 2022, 14, 601 2 of 20 emphasizes locating the optimal process parameters to achieve the target mean value where it assumes a homogeneous variance.However, in a real situation, this assumption may not be achieved.In this situation, both the mean and the standard deviation of the responses should be considered when determining the optimum conditions for the input variables.
The dual response surface methodology is the most efficient method for simultaneously optimizing the mean and the standard deviation functions of the responses to achieve the desired target while keeping the standard deviation small.It is common to employ both symmetrical and asymmetrical designs along with the regression method to model the mean and the standard deviations of the responses.Asymmetric design refers to a situation where all the factors have the same levels of variation while asymmetrical design is when not all the factors have the same levels of variation.However, symmetrical 2 k , 3 k and 4 k factorial designs of experiments are the commonly used designs in industrial engineering because a large amount of information can be acquired from a relatively small number of experiments.In dual RSM, the mean and the standard deviation of the responses at each design point are first calculated.Then, statistical models for the mean and standard deviation of the responses are established by fitting the response surface models (usually using a second-degree polynomial model) to the experimental data.The parameters of the models are then estimated prior to simultaneously optimizing both the estimated mean and the estimated standard deviation of the responses to obtain an optimal setting for the input/control variables.
Various optimization methods of dual response surface optimization (DRSO), which attempt to simultaneously optimize both the mean and the standard deviation of the responses, were proposed in the literature [4][5][6][7].Vining and Myers [8] were the first to introduce dual response surface optimization (DRSO) using response surface methodology (RSM).We denote this method as the VM method.However, their method received criticism from Del Castillo and Montgomery [4] due to using Lagrange multipliers to simultaneously optimize the process mean and process standard deviation.They noted that this approach may be impractical and imprecise for a global best solution due to the restriction of the process mean which was made to be equal to a specific target value.Thus, they proposed the generalized reduced gradient (GRD) technique with inequality constraints to rectify the drawback of the VM optimization method.Similarly, Lin and Tu [5] also pointed out the weakness of the VM method in which it did not guarantee a global optimal due to its constraint to a specific value.Hence, they proposed an optimization scheme which we refer to as the LT optimization scheme.It is based on the mean squared error (MSE) objective function that allows for a small bias.Kim and Lin [9] pointed out that the objective function of LT does not consider the measure of the violation of constraint, and it does not clearly specify how far the estimated mean response is from the desired target value.In an attempt to reduce the effect of bias and variance in obtaining the optimal setting conditions for the input variables, Ding et al. [7] modified the LT method by imposing relative weights, w on the bias and the variance terms of the MSE optimization scheme.Nonetheless, this method, which is referred to as WMSE, received criticism from Jeong et al. [10] due to the lack of consideration of the decision maker's interest in determining the weight functions.Therefore, [9][10][11] then introduced a proper procedure to determine the value of w by considering the decision maker's preference.Nevertheless, most of the existing methods do not show much improvement in terms of acquiring an estimated mean response that is closer to the target value with small variation.
Despite the shortcoming of the LT method, this optimization scheme is a widely used method to solve dual response and multiple response problems [12][13][14][15][16]. Baba et al. [17] developed a new optimization scheme by employing the penalty function (PM) approach which was found to be more efficient compared to other approaches discussed earlier.However, the weakness of this method is that it has neither a clear rule nor a systematic approach to determine the value of the penalty constant, ξ.
It is important to highlight that in most optimization methods of dual RSM including the PM, the ordinary least squares (OLS) is the commonly used method for estimating the models' parameters.The classical sample mean and the classical sample standard deviation are also popular used measures of location/central tendency and measures of spread, respectively.Unfortunately, many statistics practitioners are not aware that the classical mean, classical standard deviation, and OLS method are not resistant to outliers [18].Outlier(s) is a single or group of observations that are markedly different among the bulk of data or pattern set by most observations.It is now evident that outliers may have an unduly effect on the sample mean, sample standard deviation, and OLS estimates [19][20][21][22].As a result, using those classical approaches will give inaccurate estimates of the mean response and standard deviation of the response.
Therefore, to address the problem of outliers on the parameter estimates, robust methods which are not easily affected by outliers are put forward.Robustness refers to the ability of a method to remain unaffected or less affected by outliers [23].Many robust estimation methods such as the M, MM, LMS, and LTS can be found in the literature [24][25][26].Rousseeuw and Leroy [25] proved that the robust MM estimator is highly efficient and has a high breakdown point.
The work of [17] has inspired us to propose a new method called the penalty function method based on the decision maker's (DM) preference structure in obtaining the penalty constant, denoted as PMDM.The PMDM is the extension work of the PM optimization method, which is based on our newly developed method of determining penalty constant, ξ.To reduce the effect of outliers, robust MM estimators are incorporated in the algorithm of PMDM to estimate the mean and the standard deviation of the responses as well as to obtain the fitted models of the mean and the standard deviations.
The objectives of this study are: (1) to establish a new optimization method for simultaneously optimizing the fitted models of the mean and the standard deviation of the responses by incorporating a penalty function based on our newly proposed penalty constant, ξ; (2) to evaluate the performance of the proposed method compared to VM, LT, WMSE, and PM by using simulation study; (3) to apply the proposed method on two real datasets, namely the catapult study data and printing process study data.The significance of this study is that it can contribute to the development of a very powerful tool for simultaneously optimizing the fitted models of the mean and the standard deviation to determine optimum conditions on the models' control variables which could give an optimum response with the least bias and least variability.Researchers will find our proposed PMDM method very useful, especially to those who wish to determine the optimum combination of factor settings that result in the estimated mean response closer to the target value with the smallest bias and smallest variability, which is desirable.However, it is worth mentioning that the accuracy of the estimators relies on the selected data (i.e., experimental design).Hence, it is crucial to choose a good design for an experiment to be able to capture maximum information about its behavior.A poor design cannot produce valuable knowledge about the experiment, even if effective estimation methods are used [27,28].
The rest of the paper is structured as follows.Section 2 presents a brief review of the dual response surface optimization.The proposed Penalty Function Optimization Method is discussed in Section 3. Section 4 presents the results of a simulation study and two real datasets.Finally, concluding remarks are described in Section 5.

The Dual Response Surface Optimization
The dual response surface method aims to find the optimum setting condition of the controllable factors to diminish the variability and deviation from the desired target of the decision-maker.Three basic strategies are considered to achieve such aim: experimental design, regression fitting, and optimization aspect.Let us say an experiment was conducted to analyze the effect of some control factors for a given experimental design.Suppose the characteristics of interest, i.e., the response variable depends on a set of control variables (coded), x = x 1 , x 2 , . . ., x p .Let y ij represent the response variable at the ith design point (i = 1, 2, . . ., n) and the jth replicate (j = 1, 2, . . ., m).Table 1 presents the summary of the observed data where x, y i , s 2 i represent a set of controlled variables or factors, sample mean, and sample variance of the response, respectively.At each of the design points, the sample mean and sample standard deviation of the responses are usually calculated as in Equation (1).
Following Myers et al. [29], the control variables are firstly translated from the natural units (original) into coded as in Equation (2).
where δ is the original value of x.Max (.) and min (.) refer to the maximum and minimum values of the original values of x, respectively.
Subsequently, the second-order polynomial model response functions for the mean and standard deviation of the responses are formulated.The ordinary least squares (OLS) method is usually used to estimate the parameters of the models.The fitted dual response functions for the sample mean and sample standard deviation at each design point, denoted as ωµ and ωσ , respectively, are written as follows: and where β and γ are the estimates of the model parameters.The control variables, x i in Equations ( 3) and (4) are firstly translated from the natural units into coded units.Then, an optimization method is used to simultaneously optimize both the mean and standard deviation functions of Equations ( 3) and (4) to obtain the estimated optimal factor settings.Then, the estimated mean response and estimated standard deviation of the responses can be determined.
In this section, we briefly summarized some of the commonly used optimization methods.Prior to formulating the optimization scheme, the processes mean and process standard deviations' models need to be established at the outset.The VM method simultaneously optimizes the objective function of the process standard deviation in which it Symmetry 2022, 14, 601 5 of 20 restricts the process mean towards the desired target.Using the second-order polynomial model in Equations ( 3) and (4), the VM method of optimization scheme is to minimize ωσ subject to ωµ = τ (5) where τ is the target value.Lin and Tu [5] proposed an optimization scheme based on the squared-loss model which is known as the mean square error (MSE).This model deals with square bias and minimizes the variance component of the factor settings.The Lin and Tu (LT) method is to Ding et al. [7] proposed a natural extension of the LT method by imposing relative weights on the bias and the variance terms.Their method which is called weighted MSE where w = relative weight factor, (0 ≤ w ≤ 1).The data-driven approach is used to determine the value of w.
where ξ is the penalty constant, 0 ≤ ξ ≤ ∞.Nonetheless, this method concentrates more on the bias and forces the estimated mean response to be close to the target value, which is not adequately efficient.Moreover, no specific penalty constant value, ξ was given to Equation (8).Hence, practitioners will find these approaches difficult since they must do trial and error subjective judgment.

The Proposed Penalty Function Optimization in Dual Response Based on Decision Maker's Preference
In the preceding sections, we have reviewed some of the optimization schemes of the dual response surface methodology, and their weaknesses are noted.The penalty functionbased approach (PM) seems to be better than the other methods but still has shortcomings, whereby the determination of the penalty constant ξ is not clearly specified.Hence, we propose a new method which is called the penalty method based on the decision maker's (PMDM) preference, where the penalty constant ξ is determined by considering the decision maker's judgments.The proposed method is less complicated and can be implemented by practitioners.Prior to the development of the method, the sample mean, sample standard deviations, and the fitted response models of Equations ( 3) and (4) are estimated using the MM estimator which is highly robust and highly efficient with a high breakdown point.The resultant MM estimates are highly robust and less affected by outliers.

Robust Measures of Location and Spread
Mean and standard deviation are popularly used measures of location/central tendency and measure spread of data, respectively.However, there is evidence that they may perform poorly when outliers occur in data [23].As already mentioned, the robust method which is not easily affected by outliers is recommended.This section discusses the MM robust estimator as an alternative to the classical methods presented in Equation (1) for estimating the mean and standard deviation of the response values at each design point.This estimator, known as the MM estimator, was proposed by Yohai [26].It is highly efficient with high breakdown points and is not easily affected by outliers.Consider the following location-scale model where x 1 , x 2 , . . ., x n are n observations and ε i , i = 1, 2, . . ., n, are independent and identically distributed (i.i.d) observations with variance equals to 1.The following algorithm (see Maronna et al. [23]) summarized the MM algorithm for estimating µ, and the scale σ.
Step 1: The initial consistent estimator of the location µ 0 and scale σ 0 were computed using high breakdown point S-estimator [26].
Step 2: From Step 1, compute the residuals e i and subsequently compute the M estimate of standard deviation parameter (scale), σ where σ is the solution to where t = e i / σ and ψ = ρ 0 (t/c 0 ) must be a redescending ρ function such as Hampel, Tukey's biweight, and Tanh functions.Tukey's biweight function is employed in this paper.
Step 3: Compute the M estimate of μ using ρ 1 .Yohai [26] noted that for Tukey's bisquare weight function, employing, c 0 = 0.4685 and c 1 = 4.68, results in high breakdown and high-efficiency estimator, respectively.μ is a solution to where ψ = ρ 1 (t/c 1 ).Upon convergence, μ and σ are the estimates for the mean (MM- mean) and standard deviation (MM-standard deviation) based on the MM estimator.In the subsequent sections, they are referred to as MM l and MM s , respectively.The same procedures were applied to estimate the parameters of Equations ( 3) and ( 4) using the MM estimator with a slight change, where the polynomial regression model was considered instead of the location model.

The New Optimization Approach
The new optimization approach (PMDM) consists of two stages whereby in the first stage, penalty constant ξ is determined.Subsequently, in the second stage, Equation ( 13) is optimized to obtain the estimated optimal factor settings.

Stage 1: Determining the value ξ
Determining the value ξ is critical and challenging because of its large interval [11].The optimal penalty constant denoted as ξ * is defined as the value that can provide an estimated mean response with reasonably close to zero bias and at the same time minimizes the estimated standard deviation of the response.Here, we illustrate how the value of penalty value ξ is determined, by using, the roman catapult problem data.This data set was also used by Kim and Lin [9] to assess their proposed method.The aim of the experiment is to identify a combination of factor settings for which a variety of projectile types can be delivered to a target with minimal variability.The response variable is the distance of the projectile from the base of the catapult.The distance was measured in inches.The input variables are x 1 the arm length (0.32-3.68 inches), x 2 stop angle (30-90 • ) x 3 and pivot height (2.5-5.5 inches).The target distance of the projectile from the base is set to 80. Figure 1 depicts the corresponding estimated mean and estimated variance of responses for Catapult study data, for some values of penalty constant, ξ, in the interval of 1-200, i.e., [1,200].
From the figure, it can be seen that the estimated variance of the response increases as the ξ value increases.On the contrary, the estimated mean response decreases when the ξ value increases.However, when the ξ value reaches a certain point, the performance of both estimated responses becomes stable.In such circumstances, it is essential to diminish the unnecessary value of ξ.
response decreases when the  value increases.However, when the  value reaches a certain point, the performance of both estimated responses becomes stable.In such circumstances, it is essential to diminish the unnecessary value of .Jeong et al. [10] stated that in such a case, it is necessary to narrow down the interval to make a substantial choice of .Following the idea of [10], we firstly defined a vector  =  ,  where  = ( − ) (bias) and  =  (variance).The penalty function (PF) denoted as PF for each , 0 ≤  ≤ ∞ can be expressed as follows: The idea behind this basic premise is that the  value should be congruent with the DM's preference structure whereby the ranking of the  given by the DM should be consistent with the ordering of the corresponding PF values.For example, if the DM prefers  to  (denoted as  ≻  ), PF should be less than PF .We employed the pairwise ranking scheme that compares only one pair of vectors at each time.The maximum number of pairwise ranks is  .If the expression PF < PF is true, then the possible interval of  is obtained, otherwise, it will not be considered as the interval of .This procedure is repeated for all pairwise vectors.Next, the set of all interval  , that satisfies inequalities PF < PF is obtained by intersecting the  , for all possible pairs of vectors.
Owing to this, we propose a systematic algorithm to find the optimal penalty constant  * .Hence, the proposed algorithm of determining the penalty constant  method can be summarized as follows Step 1: Firstly compute  and  as in Equation (1).
Step 3: Calculate criteria value, i.e., CV = | − | +  .We suggest a simple algorithm of reducing the number of pairwise comparisons of  by computing △ where △ is an outlier-resistant estimator, i.e., median of CV .If CV value exceeds the △ value, then vectors of  =  ,  will not be considered in the ranking of  .
To clarify the idea, the process mean,  , process standard deviation,  , criteria values, CV along with the value  , are presented in Table 2 for Catapult data.In this example, the design point  = 1,3,6,8,9,10,11,12,13,14 were removed from the analysis since their corresponding CV exceed the median CV = 27.72,Consequently, only design point  = 2,4,5,7,15,16,17,18,19,20 are remained.After that, the remaining design Jeong et al. [10] stated that in such a case, it is necessary to narrow down the interval to make a substantial choice of ξ.Following the idea of [10], we firstly defined a vector where z 1 = (y i − τ) 2 (bias) and z 2 = s 2 i (variance).The penalty function (PF) denoted as PF ξ for each ξ, 0 ≤ ξ ≤ ∞ can be expressed as follows: The idea behind this basic premise is that the ξ value should be congruent with the DM's preference structure whereby the ranking of the z i given by the DM should be consistent with the ordering of the corresponding PF ξ values.For example, if the DM prefers z i to z j (denoted as z i z j ), PF z i ξ should be less than PF z j ξ .We employed the pairwise ranking scheme that compares only one pair of vectors at each time.The maximum number of pairwise ranks is C n 2 .If the expression PF z i ξ < PF z j ξ is true, then the possible interval of ξ is obtained, otherwise, it will not be considered as the interval of ξ.This procedure is repeated for all pairwise vectors.Next, the set of all interval ξ n i,j that satisfies inequalities PF z i ξ < PF z j ξ is obtained by intersecting the ξ n i,j for all possible pairs of vectors.Owing to this, we propose a systematic algorithm to find the optimal penalty constant ξ * .Hence, the proposed algorithm of determining the penalty constant ξ method can be summarized as follows Step 1: Firstly compute y i and s i as in Equation (1).
Step 3: Calculate criteria value, i.e., CV i = |y i − τ| + s i .We suggest a simple algorithm of reducing the number of pairwise comparisons of z i by computing where is an outlier-resistant estimator, i.e., median of CV i .If CV i value exceeds the value, then vectors of will not be considered in the ranking of z i .
To clarify the idea, the process mean, y i , process standard deviation, s i , criteria values, CV i along with the value z i , are presented in Table 2 for Catapult data.In this example, the design point i = 1, 3, 6, 8, 9, 10, 11, 12, 13, 14 were removed from the analysis since their corresponding CV i exceed the median CV i = 27.72,Consequently, only design point i = 2, 4, 5, 7, 15, 16, 17, 18, 19, 20 are remained.After that, the remaining design points were ranked in ascending order, i.e., z 1 , . . ., z 10 based on the value of CV i .The number of pairwise comparisons to be made is reduced to C k 2 = 45, i.e., z 1 , z 2 , z 1 , z 3 , z 1 , z 4 , . . ., z 9 , z 10 , where k is the remaining number of design points.Step 4: Conduct C k 2 pairwise comparisons.For each pair of z i , z j , where z i z j ), obtain the values for PF z i ξ and PF z j ξ for some values of ξ in an interval, i.e., 0 ≤ ξ ≤ a that satisfy the condition PF z i ξ < PF z j ξ .Denote the interval of the corresponding ξ for each pair in term of sets, i.e., S 1 S 2 . . .S C k 2 .Then, the intersection of these sets, i.e., S 1 S 2 . . .S C k 2 will be considered as the final set of the optimal penalty constant, ξ.

Stage 2: Optimization stage
In this stage, optimization models are formulated and solved.This yield estimated optimum values of process mean and process standard deviation of the responses.As mentioned earlier, the new procedure introduced in this paper is called a penalty method based on the decision maker's (PMDM) preference which simultaneously optimizes the process mean and process standard deviation functions of the response variables.The PMDM is defined as: where ξ * is the optimal penalty constant determined in Step 4 of Stage 1.The optimal setting (x * ) with the corresponding ξ value within the interval that gives the smallest RMSE and bias will be considered as the final optimal setting.The ωµ and ωσ from Equations ( 3) and ( 4) are the fitted values based on the MM regression estimator.It is also noted that the MM estimator is used to estimate the mean and the standard deviation of the response variables.The PMDM can also be used if it is based on classical estimators.However, we anticipate that the results based on a robust estimator are more efficient than on a classical estimator.
The PMDM in Equation ( 13) is optimized by using any unconstrained optimization method such as Newton's method, BFGS method, conjugate gradient method, and steepest ascent (descent) method.In addition, nonlinear optimization software can be used to find the optimal factor settings for dual response surface problems.In this study, the genetic algorithm (GA) optimization was implemented to obtain the estimated optimal factor settings.The R software is utilized to analyze this problem.The resultant estimated optimal factor settings are then substituted in Equations ( 3) and ( 4) to obtain the estimated optimum mean response and standard deviation of the responses.We expect that the PMDM could give an optimum response with the least bias and least variability, compared to other methods.

Results
Monte Carlo simulation study and two real datasets are presented in this section to assess the performance of the proposed PMDM compared to VM, LT, WMSE, and PM methods.All analyses were computed using the R program.The robust package in R, i.e., lmrob, rrcov, and Rsolnp are used to obtain the robust MM estimates of model parameters, process mean, and process standard deviation, respectively.

Monte Carlo Simulation Study
Monte Carlo simulation study is designed to assess the performance of our new proposed PMDM method and to compare its performance with other methods namely VM, LT, WMSE and PM.Following Park and Cho [12] and Goethals and Cho [13], a symmetric 3 3 factorial design of experiments is considered, whereby at each control factor settings, x i = (x i1 , x i2 , x i3 ), i = 1, 2, . . ., 27, five responses (y i1 , y i2 , . . . ,y i5 ) are randomly generated from normal distribution with mean µ(x i ) and standard deviation σ(x i ), as follows: where the desired target is assumed to be 500, i.e., τ = 500.We refer to this generated data as 'clean' data.The performance of the proposed method when compared to other methods is evaluated under three different frameworks as shown in Table 3.In order to see the effect of outliers on the different methods, three observations of the original responses are randomly replaced with contaminated values which are drawn from a normal distribution, i.e., N 150, 10 2 .The five methods are then applied to each set of the generated data.In each simulation run, 500 and 1000 iterations are considered.The performance of the proposed method is evaluated based on bias, standard error (SE) and root mean squares errors (RMSE) of the estimated optimal mean response, μ.The estimated mean and variance of the estimated optimal mean response, computed over k iterations are given by μ = ∑ k i=1 ( μ/k) and ∑ k i=1 μ − μ /k, respectively.The bias is computed as bias = μ − 500 and the root mean square error (RMSE) = (bias) 2 + var( μ).The MSE consists of two components; one measures the variability (precision) and the other measures its bias (accuracy).A good method is one that has the least values of bias, standard error (SE), and RMSE.The bias, standard error (SE), and RMSE of the estimated optimal mean response for models A, B, and C are exhibited in Tables 4-6, respectively.First, the focus is on the results of Table 4 for clean data where the process mean estimates, process standard deviation estimates, and regression estimates were calculated based on the classical methods.It can be observed that our proposed PMDM-optimization method outperforms the other four methods evident by having the smallest bias and RMSE for both the simulation iterations.
Similar results are obtained in Table 5 for the data having three outliers with the classical framework as described in Table 3 (model B).It can be observed from Table 5 that the biases and RMSE of all methods have markedly increased.This is due to the use of the classical methods to compute the estimates of mean and standard deviation and also the use of the OLS method in obtaining the estimates of the predicted models in Equations ( 3) and (4).However, our proposed PMDM method is still the most efficient method among the five optimization schemes because it has the smallest value of RMSE, followed by the PM, LT, WMSE, and VM methods.
The bias, SE, and RMSE of the estimated optimal mean responses for contaminated data based on robust estimators (model C) are presented in Table 6.It is interesting to observe that the bias, SE, and RMSE for the five estimates have significantly reduced.These attractive results are due to the usage of the robust MM estimator, which is resistant to outliers.The results of our proposed PMDM are very appealing.Similar to models A and B, the proposed PMDM emerges to be conspicuously the most efficient method.It has the smallest bias, SE, and RMSE among the five methods.Again, the performance of PM comes second and is followed by the LT, WMSE, and VM methods.Our proposed PMDM optimization method based on the robust MM-fitted response surface function has substantially improved the estimated optimal mean responses in terms of having the smallest bias, SE, and RMSE.Smaller RMSE indicates that the combination of optimum factor settings results in the estimated mean response being closer to the target value with less bias and less variability, which is desirable.

The Catapult Study Data
The roman style catapult problem data is the first real example to show the merit of our proposed PMDM method when compared to the other four methods.This data set was also used by Kim and Lin [9] and Baba et al. [17] to assess their proposed methods.The central composite design (CCD) with three replicates was evaluated at each design point.Following Kim and Lin [9], the target of interest for distance point is assumed to be 80, i.e., τ = 80.
In order to see the effect of an outlier on the estimated optimal mean response, we deliberately changed one data point, that is the y 3,20 observation corresponding to y value equals 79 (in bold) with higher values of 790. Figure 2a,b shows the boxplot of the original and modified data sets with one outlier.It can be observed from the boxplot that the distribution of the response variable of the original data set is fairly close to a normal distribution.However, the response variables for the modified data are skewed to the right with one outlier detected.The data set is shown in Table 7 (the control variables, x i are coded and the original values are in parenthesis) while in Table 8, the mean, standard deviation, MM-mean, MM-standard deviation, and criteria value (CV) for model A, model B, and model C are presented.What is immediately clear from Table 8 is that the process mean and process standard deviation based on the sample mean and sample standard deviation are very sensitive to the contaminated data (shown in bold) when compared with the proposed estimation approach using MM-mean and MM-standard deviation.

The Catapult Study Data
The roman style catapult problem data is the first real example to show the mer our proposed PMDM method when compared to the other four methods.This dat was also used by Kim and Lin [9] and Baba et al. [17] to assess their proposed meth The central composite design (CCD) with three replicates was evaluated at each de point.Following Kim and Lin [9], the target of interest for distance point is assumed t 80, i.e.,  = 80.
In order to see the effect of an outlier on the estimated optimal mean response deliberately changed one data point, that is the  , observation corresponding value equals 79 (in bold) with higher values of 790. Figure 2a,b shows the boxplot o original and modified data sets with one outlier.It can be observed from the boxplot the distribution of the response variable of the original data set is fairly close to a nor distribution.However, the response variables for the modified data are skewed to right with one outlier detected.The data set is shown in Table 7 (the control variable are coded and the original values are in parenthesis) while in Table 8, the mean, stand deviation, MM-mean, MM-standard deviation, and criteria value (CV) for mode model B, and model C are presented.What is immediately clear from Table 8 is tha process mean and process standard deviation based on the sample mean and sam standard deviation are very sensitive to the contaminated data (shown in bold) when c pared with the proposed estimation approach using MM-mean and MM-standard de tion.If all of the 20 alternative vectors were to be compared, the total number of pairwise comparisons would be 190.In addition, the number of pairwise comparisons can be reduced by screening out the vectors in which their CVs exceed the median of CV.Similar to the simulation experiment, we consider three model scenarios as depicted in Table 3; model A (the classical method is used to estimate the mean and the standard deviation of the responses for clean data), model B (the classical method is used to estimate the mean and the standard deviation of the responses for contaminated data), and model C (the robust MM method is used to estimate the mean and the standard deviation of the responses for contaminated data).From Table 8, it can be observed that only 10 remaining vectors for data with model A and model C, and 9 vectors for data with model B were considered out of 20 vectors.Pairwise comparisons among the vectors are conducted.For models A and B, they consist of C 10 2 while for model C, C 9 2 pairwise comparisons.Figure 3 shows that the proposed method has established a set value of ξ as a result of the intersection of the I(i, j).First, in Figure 3a, the resultant intersection plot is determined by I(3, 5) and I(6, 7) which corresponds to the lower and upper limits of ξ, i.e., ξ = 5 and ξ = 6.Then, the final set of optimal penalty constant ξ * to be used for model A is ξ * A [5, 6].Secondly, for model B, the lower and upper limits are ξ = 2 and ξ = 3 which are determined by I(2, 4) and I(4, 5).Hence, the optimal penalty constant is ξ * B [2, 3] as shown in Figure 3b.Finally, Figure 3c shows the final penalty constant for model C, i.e., ξ * C [1, 5].Subsequently, the optimal setting (x * ) for each model is obtained by minimizing the PMDM in Equation ( 13) using the resultant ξ * values.The optimal setting (x * ) with the corresponding ξ value within the interval that gives the smallest RMSE and bias will be considered as the final optimal setting.duced by screening out the vectors in which their CVs exceed the median of CV.Similar to the simulation experiment, we consider three model scenarios as depicted in Table 3; model A (the classical method is used to estimate the mean and the standard deviation of the responses for clean data), model B (the classical method is used to estimate the mean and the standard deviation of the responses for contaminated data), and model C (the robust MM method is used to estimate the mean and the standard deviation of the responses for contaminated data).From Table 8, it can be observed that only 10 remaining vectors for data with model A and model C, and 9 vectors for data with model B were considered out of 20 vectors.Pairwise comparisons among the vectors are conducted.For models A and B, they consist of  while for model C,  pairwise comparisons.
Figure 3 shows that the proposed method has established a set value of  as a result of the intersection of the (, ).First, in Figure 3a, the resultant intersection plot is determined by  (3,5) and (6,7) which corresponds to the lower and upper limits of , i.e.,  = 5 and  = 6.Then, the final set of optimal penalty constant  * to be used for model A is  *  5, 6 .Secondly, for model B, the lower and upper limits are  = 2 and  = 3 which are determined by (2,4) and (4,5) .Hence, the optimal penalty constant is  *  2, 3 as shown in Figure 3b.Finally, Figure 3c shows the final penalty constant for model C, i.e.,  *  1, 5 .Subsequently, the optimal setting ( * ) for each model is obtained by minimizing the PMDM in Equation (13) using the resultant  * values.The optimal setting ( * ) with the corresponding  value within the interval that gives the smallest RMSE and bias will be considered as the final optimal setting.The VM, LT, WMSE, and PM methods are then applied to the data and their performances are compared with the PMDM.The estimated optimal factor setting, estimated optimal mean response, estimated standard deviation, and RMSE for models A, B, and C are exhibited in Tables 9-11, respectively.It is interesting to observe the results of our proposed PMDM in Tables 9-11 for clean and contaminated data.Let us first focus on Table 9 for clean data.The PMDM method shows a slightly better result when compared to other methods for clean data.It can be observed from Table 9 that the estimated mean optimal response of our proposed PMDM is the closest to the target value and has the least standard deviation and RMSE compared to other methods.Similar results are obtained for model B (using classical methods to obtain various estimates for contaminated data) in which our PMDM method has the smallest standard deviations and RMSE values among the five methods.As displayed in Table 11 for contaminated data (model C) that use robust methods to obtain various estimates, our proposed PMDM method has also outperformed the other methods evident by having the smallest value of standard deviation and RMSE.Moreover, the estimated mean optimal response of the PMDM is very close to the target value, i.e., 80 with the least standard deviation, which is very desirable.The results of this data set agree reasonably well with the results of the simulation study.Here, we present the interpretation of the results, in particular for the PMDM method, model C as in Table 11 after the coded units of the estimated factor settings are translated back to natural units using Equation (2).The results indicate that the optimum combination of factors where arm length at 1.7628 inches, stop angle at 65.0340, and pivot height at 3.603 inches would give a maximum value of 79.9856 inches of the response, i.e., distance point where a projectile lands from the base of the catapult.The VM, LT, WMSE, and PM methods are then applied to the data and their performances are compared with the PMDM.The estimated optimal factor setting, estimated optimal mean response, estimated standard deviation, and RMSE for models A, B, and C are exhibited in Tables 9-11, respectively.It is interesting to observe the results of our proposed PMDM in Tables 9-11 for clean and contaminated data.Let us first focus on Table 9 for clean data.The PMDM method shows a slightly better result when compared to other methods for clean data.It can be observed from Table 9 that the estimated mean optimal response of our proposed PMDM is the closest to the target value and has the least standard deviation and RMSE compared to other methods.Similar results are obtained for model B (using classical methods to obtain various estimates for contaminated data) in which our PMDM method has the smallest standard deviations and RMSE values among the five methods.As displayed in Table 11 for contaminated data (model C) that use robust methods to obtain various estimates, our proposed PMDM method has also outperformed the other methods evident by having the smallest value of standard deviation and RMSE.Moreover, the estimated mean optimal response of the PMDM is very close to the target value, i.e., 80 with the least standard deviation, which is very desirable.The results of this data set agree reasonably well with the results of the simulation study.Here, we present the interpretation of the results, in particular for the PMDM method, model C as in Table 11 after the coded units of the estimated factor settings are translated back to natural units using Equation (2).The results indicate that the optimum combination of factors where arm length at 1.7628 inches, stop angle at 65.0340, and pivot height at 3.603 inches would give a maximum value of 79.9856 inches of the response, i.e., distance point where a projectile lands from the base of the catapult.

Printing Process Study Data
The second real example in this paper is the printing process study data [1].The experiment introduced by [1] was conducted to determine the effect of the three control variables: x 1 (speed), x 2 (pressure), and x 3 (distance) on the characteristics of the printing process, which include the machine's index to apply colored inks to package labels.A symmetric 3 3 factorial design of experiments with three replicates at each of the 27 design points is considered, assuming that the target of interest for distance point is τ = 500.
In order to see the effect of an outlier on the estimated mean optimal response, we deliberately changed one data point, that was the y 1,8 observation corresponding to y value equaling 259 (in bold) with higher values of 9259. Figure 4a,b displays the boxplot of the original and the modified data sets with an outlier.From the boxplot, it suggests that the response variable of the original data is reasonably close to normal distribution.On the other hand, the modified data of the response variables are skewed to the right with one outlier detected.The data set is shown in Table 12 (the control variable, x i are coded) while in Table 13, the mean, standard deviation, MM-mean, MM-standard deviation, and criteria value (CV) for model A, model B, and model C are presented.

Printing Process Study Data
The second real example in this paper is the printing process study data [1].The experiment introduced by [1] was conducted to determine the effect of the three control variables:  (speed),  (pressure), and  (distance) on the characteristics of the printing process, which include the machine's index to apply colored inks to package labels.A symmetric 3 3 factorial design of experiments with three replicates at each of the 27 design points is considered, assuming that the target of interest for distance point is  = 500.
In order to see the effect of an outlier on the estimated mean optimal response, we deliberately changed one data point, that was the  , observation corresponding to  value equaling 259 (in bold) with higher values of 9259. Figure 4a,b displays the boxplot of the original and the modified data sets with an outlier.From the boxplot, it suggests that the response variable of the original data is reasonably close to normal distribution.On the other hand, the modified data of the response variables are skewed to the right with one outlier detected.The data set is shown in Table 12 (the control variable,  are coded) while in Table 13,    If all 27 alternative vectors were to be compared, the total number of pairwise comparisons would be 351.In addition, the number of pairwise comparisons can be reduced by screening out the vectors in which their CVs exceed the median of CV.Similar to the simulation experiment and Catapult data example, we consider three model scenarios, i.e., models A, B, and C as depicted in Table 2. From Table 13, it can be observed that there are only 13 remaining vectors of data that include all three models.Thus, C 13 2 pairwise comparisons among the vectors are considered.
Figure 5 shows that the proposed method has established a set value of ξ as a result of the intersection of the I(i, j).First, in Figure 5a, the resultant intersection plot is determined by I(3, 4) which corresponds to the lower and upper limits of ξ, i.e., ξ = 3 and ξ = 4. Then the final set of optimal penalty constant ξ * to be used for models A and B is ξ * A,B [3, 4].Finally, Figure 5b   If all 27 alternative vectors were to be compared, the total number of pairwise comparisons would be 351.In addition, the number of pairwise comparisons can be reduced by screening out the vectors in which their CVs exceed the median of CV.Similar to the simulation experiment and Catapult data example, we consider three model scenarios, i.e., models A, B, and C as depicted in Table 2. From Table 13, it can be observed that there are only 13 remaining vectors of data that include all three models.Thus,  pairwise comparisons among the vectors are considered.
Figure 5 shows that the proposed method has established a set value of  as a result of the intersection of the (, ).First, in Figure 5a, the resultant intersection plot is determined by  (3,4) which corresponds to the lower and upper limits of , i.e.,  = 3 and  = 4. Then the final set of optimal penalty constant  * to be used for models A and B is  , *  3, 4 .Finally, Figure 5b shows the final penalty constant for model C, i.e.,  *  2, 4 .All five methods were used to determine the optimal factor setting, optimal mean response, standard deviation, and RMSE, which are exhibited in Tables 14-16, respectively.Similar to the Catapult data, it can be observed from Table 14 that the PMDM shows a slightly better result when compared to other methods for clean data.It is also interesting to observe that the PMDM consistently provides the smallest standard deviation and RMSE for models B and C when outliers are present in the data.The results of the real examples were consistent with the results of the simulation study.It is very interesting to note that the results that show the estimated mean optimal response of the PMDM is very close to the target value, i.e., 500 with the least standard deviation, which is very desirable.

Conclusions
The main aim of this study is to develop a more reliable alternative scheme (PMDM) that simultaneously optimizes both the mean and standard deviation functions of the response variables in the response surface model.Robust MM-mean, robust MM-standard deviation estimator, and MM regression estimator which are highly efficient and have high breakdown points were employed instead of using the classical non-robust methods.We also establish a systematic method to determine the penalty constant ξ which is based on the decision maker's (DM) preference structure and integrate this value in the PMDM optimization scheme.The existing optimization scheme is not effective enough as they fail to obtain an estimated mean response that is closer to the target value with a small variation.Moreover, those methods are easily affected by outliers due to the use of nonrobust methods to estimate the mean and the standard deviation of the response variables, and the non-robust least squares (OLS) method to estimate the parameters of the mean and standard deviation functions.The results obtained from the real data set and the Monte Carlo simulation study show that our proposed PMDM method is more efficient than the existing VM, LT, WMSE, and PM methods, irrespective of whether or not the outlier exists in the data set.Nonetheless, our proposed PMDM method has its own shortcoming as its computation running time is longer than that of the other methods in this study.A slightly longer running time compared to the other methods is a trade-off one has to consider when using our PMDM method as the method is very efficient in determining the optimum combination of factor settings that result in the estimated mean response closer to the target value with the smallest bias and smallest variability, which is desirable.In future work, we Symmetry 2022, 14, 601 19 of 20 look forward to developing a reliable method of estimating the parameters of the mean and standard deviations functions of the responses for an unbalanced design which is prone to produce heteroscedastic errors.After that, our proposed PMDM method will be employed to simultaneously optimize both functions to obtain more efficient optimal factor settings, hence improving processes.

Figure 1 .
Figure 1.The estimated mean and estimated variance versus penalty constant ().

Figure 1 .
Figure 1.The estimated mean and estimated variance versus penalty constant (ξ).

Figure 3 .
Figure 3.The individual set of ξ for each model for Catapult study data set.(a) ξ * A for model A, (b) ξ * B for model B and (c) ξ * C for model C.
the mean, standard deviation, MM-mean, MM-standard deviation, and criteria value (CV) for model A, model B, and model C are presented.(a) (b)

Figure 4 .
Figure 4. Boxplot for (a) original (no outlier) data set and (b) modified data set with outlier for Printing process data.

Figure 5 .
Figure 5.The individual set of ξ for each model for printing study data set (a) ξ * A and ξ * B for model A and B and (b) ξ * C for model C.

Table 1 .
The presentation of the observed data for a given experimental design.

Table 2 .
Summary of ranking process.

Table 3 .
Framework for the purpose of performance comparison.

Table 4 .
Estimated bias, SE, and RMSE of the optimal mean response for model A (classical method were used, data without outliers).

Table 5 .
Estimated bias, SE, and RMSE of the optimal mean response for model B (classical method were used, data with outliers).

Table 6 .
Estimated bias, SE, and RMSE of the optimal mean response for model C (robust method were used, data with outliers).

Table 7 .
The Catapult study data.

Table 8 .
Process mean and process standard deviation with the criteria value (CV) and ranking code for Catapult study data set.

Table 9 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response and RMSE for Catapult study data set, (model A, classical methods were used, data without outlier).

Table 9 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response and RMSE for Catapult study data set, (model A, classical methods were used, data without outlier).

Table 10 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response and RMSE for Catapult study data set, (model B, classical methods were used, data with outlier).

Table 11 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response and RMSE for Catapult study data set, (model C, robust methods were used, data with outlier).

Table 10 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response and RMSE for Catapult study data set, (model B, classical methods were used, data with outlier).

Table 11 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response and RMSE for Catapult study data set, (model C, robust methods were used, data with outlier).

Table 12 .
The printing process study data.

Table 13 .
Process mean and process standard deviation with the criteria value (CV) and ranking code for printing study data set.

Table 14 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response, and RMSE for printing study data set, (model A, classical methods were used, data without outlier).

Table 15 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response, and RMSE for printing study data set, (model B, classical methods were used, data with outlier).

Table 16 .
Estimated optimal settings, estimated mean response, estimated standard deviation of response, and RMSE for printing study data set, (model C, robust methods were used, data with outlier).