# An Alternative Approach of Dual Response Surface Optimization

Authors: Ishaq Baba; Habshah Midi; Sohel Rana; Gafurjan Ibragimov
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The dual response surface for simultaneously optimizing the mean and variance models as separate functions suffers some deficiencies in handling the tradeoffs between bias and variance components of mean squared error (MSE). In this paper, the accuracy of the predicted response is given a serious attention in the determination of the optimum setting conditions. We consider four different objective functions for the dual response surface optimization approach. The essence of the proposed method is to reduce the influence of variance of the predicted response by minimizing the variability relative to the quality characteristics of interest and at the same time achieving the specific target output. The basic idea is to convert the constraint optimization function into an unconstraint problem by adding the constraint to the original objective function. Numerical examples and simulations study are carried out to compare performance of the proposed method with some existing procedures. Numerical results show that the performance of the proposed method is encouraging and has exhibited clear improvement over the existing approaches.


---

## Body

## 1. Introduction

Response surface methodology (RSM) is a design of experimental technique which shows relationship between several designs and response variables. The goal of the experimenter is to determine the optimal settings for the design variables that minimize or maximize the fitted response. For more explanation on response surface techniques see [1–3]. Most of the early work in RSM is centered on a single response problem. This methodology works effectively under the assumption of the homogeneous variance of the response. However, such an assumption may not hold in solving real-life applications. Myers and Carter [4] suggested the need for developing statistical methodology known as dual response surface methodology, which can simultaneously optimize the mean and the variance function as to achieve the desired target while keeping the variance small. Generally, they defined the two responses as primary and secondary. The objective is to find the condition x on the design factors that minimize or maximize the primary response function yp(x) subject to the secondary response ys(x). In order to achieve this, three basic strategies are involved: experimental design, regression fitting, and optimization aspect. For the regression fitting, the method of the least squares is usually used to obtain the adequate response functions for the processes mean and variance by assuming that the collected data comes from a normal distribution. While, in the optimization stage, the interest is on what to optimize (i.e., determination of the objective function) and how to optimize (the optimization algorithm), in this paper, we propose a new optimization technique in dual response surface methodology based on the penalty function method for simultaneously optimizing both the location and scale functions. The usefulness of our newly proposed method for estimating the mean and variance of the optimal mean response is studied by some well known data sets and a simulation study. The outline of this paper is organized as follows. In the next section, we discussed some basic concept of the dual response surface followed by the description of the proposed method in Section 3. Numerical examples and simulation study are given in Sections 4 and 5, respectively. Finally, the conclusion is given in Section 6.


## 2. Dual Response Surface Review

Dual response surface technique consists of finding the optimum setting condition of the controllable factors in order to diminish the performance variability and deviation from the desired target of the decision maker. This method is an extension of the standard ridge analysis procedure which was introduced by Myers and Carter [4]. Ridge analysis has been used by researchers in searching the optimum setting condition for a single response problem [1, 2, 5]. The dual response used mean and variance as separate functions for the system under examination. Then, these functions are optimized based on the chosen optimization technique to determine the optimum operating conditions of the system. Following the strategy of Vining and Myers [6], the mean and the standard deviation fitted response surfaces can be written as(1)ω^μ=b0+x′b+x′Bx,where (2)x=x1x2⋮xk,b=β^1β^2⋮β^k,B=122β^11β^12⋯β^1kβ^122β^22⋯β^2k⋮⋮⋱⋮β^1kβ^2k⋯2β^kk,(3)ω^σ=c0+x′c+x′Cx,where(4)x=x1x2⋮xk,c=α^1α^2⋮α^k,C=122α^11α^12⋯α^1kα^122α^22⋯α^2k⋮⋮⋱⋮α^1kα^2k⋯2α^kk,where b0, c0, b, c, B, and C are estimated vectors and matrices of coefficients obtained from the least squares method. The simultaneous optimization of (1) and (3) using the Lagrangian multipliers method was proposed by Vining and Myers [6]. Lin and Tu [7] noted that the Vining and Myers approach does not always guarantee global optimum solutions due to the restriction of the optimization to equality constraints. Based on this, they proposed the minimization of mean squared error model (MSE) by introducing a slight bias in order to minimize the variability in the responses. This method includes two major parameters (the bias and the variance). Copeland and Nelson [8] observed that minimizing MSE function does not specify how large the estimated mean ω^μ might be from the specified target value T. Instead, they modified VM model by placing some restriction on ω^μ that minimizes ω^σ subject to (ω^μ-T)2≤Δ2. Kim and Lin [9] introduced a fuzzy modeling methodology using the idea of desirability function method. Several other techniques for solving the dual response surface problem have been presented, for example [10–12] proposed modification of mean squared error model. Furthermore, [13–16] presented a robust design for contaminated and nonnormal data using squared loss optimization scheme, a highly efficient and outlier resistant robust design estimator, a dual response approach to multiple response robust design problem, a multivariate robust design using MSE and dual response modeling, and robust parameter design, respectively. Recently, a biobjective robust design model has been developed in [17] and a robust cutting parameter design using computer simulation has been studied in [18]. In many real world situations, the experimenter or the decision maker often needs to keep a balance between the process mean and the process variance to achieve the desired target. It is known that getting all efficient solutions with the class of LT methods is challenging due to the large resulting process variance. However, most of the preceding optimization schemes are derived from the LT model, except the Vining and Myers model which is basically based on the Lagrangian multipliers approach. In this paper, an alternative objective function is considered based on the penalty function method.


## 3. Proposed Optimization Scheme for Dual Response Surface

In the present study, we present a new optimization technique for dual response surface methodology based on the penalty function method. The penalty function approach swaps a constrained optimization problem by a sequence of unconstrained optimization problems whose approximate solution ideally converges to a true solution of the original constrained problem. The unconstrained problem is formulated by adding a penalty term to the original objective function which consists of the penalty parameter multiplied by a measure of violation of the constraints [19, 20]. Consider the general formulation of the constrained optimization problem given below:(5)Minimize:fxSubject  to:gjx≤0,j=1,2,…m,hix=0,i=1,2,…,n.By applying the penalty function method, we can obtain the solution of (5) using the modified objective function:(6)Fx=fx+∑j=1mμgjx+∑i=1nμhix,where f(x) is the original objective function to be minimized and gj(x) and hi(x) are set of inequality and equality constraints, respectively. This paper specifically considers a quadratic penalty function of the form(7)Fx,μ=fx+μ2∑i=1nhix2,where μ is called the penalty constant which penalizes the equality constraints when the constraints relations are not satisfied. For the purpose of clarity, we replaced f(x) and hi(x) in (7) with ω^σ and (ω^μ-T) and then write the following quadratic unconstrained minimization problem as(8)min⁡⁡Fx,μ=ω^σ+μ2ω^μ-T2,where ω^μ is the fitted response surface for mean, ω^σ is the fitted response surface for the standard deviation function, and T is the target value (usually specified by experimenter). If μ=∞, the method gives exact solution. Since in this case it is necessary that ω^μ=T which implies the bigger penalty parameter μ, thus the more exact solution is achieved. For (8), we can apply any unconstrained optimization method such as Newton’s method, BFGS method, Conjugate gradient method, and Steepest ascent (descent) method. Moreover, any nonlinear optimization software may be used to find the optimal design settings for the dual response surface problem. We used the package Rsolnp introduced in [21, 22] in R language, which is open source statistical software to perform the numerical computations and analysis. Our aim is to find an optimum solution such that the estimated mean value will be very close or equal to the target value T, while the variance is kept small. The proposed approach has some advantages over some existing methods. Firstly, the proposed method takes into consideration the measure of violation of the constraint, whereas the VM [6] and the class of LT methods are minimized without regard to the relative magnitude of violation of constraint. Secondly, the penalty parameter in (8) forced the (ω^μ-T)2 to be close to zero or equal to zero so as to achieve the target of the experimenter (decision maker). Therefore, we anticipate the optimal setting condition obtained by the proposed method would be more efficient in terms of the contribution of both bias and variance components of of the estimated mean response compared with other existing methods.


## 4. Simulation Study and Results

In this section, a simulation study is conducted to assess the performance of the newly proposed method and compare it with the commonly used methods, such as VM, LT, and WMSE. Following [13, 23] the five responses (yi1,…,yi5) are randomly generated from a normal distribution with mean ωμ and standard deviation ωσ at each control factor of settings xi=(xi1,xi2,xi3), i=1,…,27. The mean ωμ and the standard deviation ωσ are given as(9)ωμ=500+x1+x2+x32+x1+x2+x3,ωσ=100+x1+x2+x32+x1+x2+x3,where T=500. All the four methods were then applied to the data. The total of 500, 1000, and 2000 iterations is considered. Some summary values, such as the estimated mean of the optimal mean response, computed over m iterations are defined by μ^¯=∑i=1m(μ^/m), bias = μ^¯-500, and var(μ^)=∑i=1m(μ^-μ^¯)2/m. The mean squared error, denoted by (MSE), is written as MSE(μ^)=(Bias)2+var(μ^). Hence, the root mean squared error (RMSE) is given by [MSE(μ^)]1/2. The bias, standard error, and root mean squared error (RMSE) of the estimates of the optimal mean response are exhibited in Table 1.Table 1
Estimated bias, standard error (SE), and RMSE of the optimal mean response.
Methods
500 iterations
1000 iterations
2000 iterations
Bias
SE
RMSE
Bias
SE
RMSE
Bias
SE
RMSE
VM
0.84
4.38
4.46
1.30
4.41
4.52
0.96
4.53
4.64
LT
1.03
3.75
3.97
1.46
3.73
4.01
1.48
3.82
4.10
WMSE
1.32
3.82
4.04
1.47
3.78
4.06
1.50
3.87
4.14
PM
0.07
0.22
0.22
0.07
0.23
0.24
0.07
0.22
0.24Figures1 and 2 show the estimated bias and mean squared error based on the total number iteration for the various methods. It can be observed that the bias of the VM estimate is smaller than the LT and WMSE estimates. However, its RMSE is the largest among the three estimates since the variance of the VM estimate makes up most of the MSE. It is interesting to see that our proposed method is the best in terms of the smallest bias and RMSE values. Due to space constraint, Figure 3 presents kernel density estimates of VM, LT, WMSE, and PM for 1000 iterations only. The plotted results indicate that the proposed approach has a good behavior in which it is very close to the desired target. Therefore, one can say that the behavior of the constructed objective function based on the penalty function technique is more efficient and robust than other existing methods for solving dual response surface optimization problem.Figure 1
Estimated bias for the various methods.Figure 2
Estimated MSE for the various methods.Figure 3
Kernel density estimates for various methods.


## 5. Numerical Examples

### 5.1. Printing Process Study Data

To show a clear comparison, we consider the data set used by Vining and Myers [6] and Lin and Tu [7] which is given in Table 2. The experiment was conducted to determine the effect of the three variables x1 (speed), x2 (pressure), and x3 (distance) on the quality of the printing process, that is, on the machines ability to apply colored inks to package labels. The experiment is a 33 factorial design with 3 replicates at each point. Firstly, the average and variance of the 3 responses at each design point are computed, respectively. Vining and Myers [6] used the least squares method to fit a quadratic response surface model for mean and standard deviation as follows:(10)ω^μ=327.6+177.0x1+109.4x2+131.5x3+32.0x12-22.4x22-29.1x32+66.0x1x2+75.5x1x3+43.6x2x3,ω^σ=34.9+11.5x1+15.3x2+29.2x3+4.2x12-1.3x22-16.8x32+7.7x1x2+5.1x1x3+14.1x2x3.Based on the models in (10), Table 3 gives the summary of the results for the four different approaches obtained using the cuboidal region -1<xi<1, i=1,2,3. The optimal setting, estimated mean response, estimated standard deviation, and RMSE are presented in Table 3. The RMSE is calculated using the formula RMSE=(ω^μ-T)2+ω^σ21/2, where T=500. It can be seen that the VM approach leads to the optimum setting (x1,x2,x3)=(0.62,0.23,0.1) which resulted in an estimated mean response of (501.57) and root mean squared error of (51.94). This optimal setting is very close to the target mean response but with larger RMSE. The second approach LT produces a target value of (494.69) and RMSE of (44.78) which is smaller than the RMSE of the VM. The third approach is the WMSE which gives the estimated mean response of (496.44) and RMSE of (44.81) which indicate a slight increase in RMSE and little improvement in mean response. The overall performance of the proposed approach is better than those three approaches mentioned with approximate mean value of (500.00) and RMSE of (44.75). Similar procedure is repeated in the next example in order to demonstrate a clear advantage of using the proposed method in terms of closeness to target mean response and the smallest RMSE.Table 2
The printing process study data.
Index
x
1
x
2
x
3
y
u
1
y
u
2
y
u
3
y
-
u
s
u
1
−1
−1
−1
34
10
28
24.0
12.49
2
0
−1
−1
115
116
130
120.4
8.39
3
1
−1
−1
192
186
263
213.7
42.80
4
−1
0
−1
82
88
88
86
3.46
5
0
0
−1
44
178
188
136.7
80.41
6
1
0
−1
322
350
350
340.7
16.17
7
−1
1
−1
141
110
86
112.3
27.57
8
0
1
−1
259
251
259
256.3
4.62
9
1
1
−1
290
280
245
271.7
23.63
10
−1
−1
0
81
81
81
81.0
0.00
11
0
−1
0
90
122
93
101.7
17.67
12
1
−1
0
319
376
376
357.0
32.91
13
−1
0
0
180
180
154
171.3
15.01
14
0
0
0
372
372
372
372.0
0.00
15
1
0
0
541
568
396
501.7
92.5
16
−1
1
0
288
192
312
264.0
63.50
17
0
1
0
432
336
513
427.0
88.61
18
1
1
0
713
725
754
730.7
21.08
19
−1
−1
1
364
99
199
220.7
133.80
20
0
−1
1
232
221
266
239.7
23.46
21
1
−1
1
408
415
443
422.0
18.52
22
−1
0
1
182
233
182
199.0
29.45
23
0
0
1
507
515
434
485.3
44.64
24
1
0
1
846
535
640
673.7
158.20
25
−1
1
1
236
126
168
176.7
55.51
26
0
1
1
660
440
403
501.0
138.90
27
1
1
1
878
991
1161
1010.0
142.50Table 3
Comparison with other methods using printing process study data.
Method
Optimal settings
ω
^
μ
ω
^
σ
RMSE
VM
(0.62, 0.2300, 0.10)
501.57
51.92
51.94
LT
(1.00, 0.07, −0.25)
494.69
44.46
44.78
WMSE
(1.00, 0.08,  −0.25)
496.44
44.67
44.81
PM
(1.04, 0.02,  −0.23)
500.00
44.75
44.75


### 5.2.  The Catapult Study Data

This example will consider the data used by Luner [24] and Kim and Lin [9]. Three variables, x1 (arm length), x2 (stop angle), and x3 (pivot height), are under consideration to predict the distance to the point where a projectile landed from the base of the roman style catapult. The experiment is a central composite design with three replicates as given in Table 5. In [9], the fitted second order polynomial regression models for the mean and standard deviation functions are given by(11)ω^μ=84.88+15.29x1+0.24x2+18.80x3-0.52x12-11.80x22+0.39x32+0.22x1x2+3.60x1x3-4.42x2x3,ω^σ=4.53+1.84x1+4.2x2+3.73x3+1.16x12+4.40x22+0.94x32+1.20x1x2+0.73x1x3+3.49x2x3.Here, the assumed target mean value is T=80. Table 4 shows the estimated mean and RMSE of the estimated optimal mean response. It is evident from Table 4 that the proposed method outperformed the other existing procedures.Table 4
Comparison with other methods using catapult study data.
Method
Optimal settings
ω
^
μ
ω
^
σ
RMSE
VM
(0.09, 0.65, 0.00)
81.13
2.79
3.00
LT
(−0.12, 1.00, 0.39)
79.90
2.59
2.59
WMSE
(−0.14, −1.00, 0.35)
78.68
2.54
2.87
PM
(−0.12,  −1.00,  0.39)
80.00
2.59
2.59Table 5
The catapult study data.
Index
x
1
x
2
x
3
y
u
1
y
u
2
y
u
3
y
-
u
s
u
1
−1
−1
−1
39
42
42
38.3
4.0
2
−1
−1
1
80
91
71
80.7
10
3
−1
1
−1
52
45
44
47
4
4
−1
−1
1
97
60
68
75
19.5
5
1
−1
−1
60
68
53
60.3
7.5
6
1
1
1
113
127
104
114.7
11.6
7
1
1
−1
78
65
64
69
7.8
8
−1.682
0
1
130
75
79
94
30.7
9
1.682
0
0
59
60
51
56.7
4.9
10
0
−1.682
0
115
117
102
111.3
8.1
11
0
1.682
0
50
57
43
50.0
7.0
12
0
0
−1.682
88
43
49
60.0
24.4
13
0
0
1.682
54
60
50
54.7
5.0
14
0
0
0
122
119
109
116.7
6.8
15
0
0
0
87
89
78
84.7
5.9
16
0
0
0
86
85
79
83.3
3.8
17
0
0
0
88
87
81
85.3
3.8
18
0
0
0
89
87
82
86.0
3.6
19
0
0
0
86
88
79
84.3
4.7
20
0
0
0
88
90
79
85.7
5.9


## 5.1. Printing Process Study Data

To show a clear comparison, we consider the data set used by Vining and Myers [6] and Lin and Tu [7] which is given in Table 2. The experiment was conducted to determine the effect of the three variables x1 (speed), x2 (pressure), and x3 (distance) on the quality of the printing process, that is, on the machines ability to apply colored inks to package labels. The experiment is a 33 factorial design with 3 replicates at each point. Firstly, the average and variance of the 3 responses at each design point are computed, respectively. Vining and Myers [6] used the least squares method to fit a quadratic response surface model for mean and standard deviation as follows:(10)ω^μ=327.6+177.0x1+109.4x2+131.5x3+32.0x12-22.4x22-29.1x32+66.0x1x2+75.5x1x3+43.6x2x3,ω^σ=34.9+11.5x1+15.3x2+29.2x3+4.2x12-1.3x22-16.8x32+7.7x1x2+5.1x1x3+14.1x2x3.Based on the models in (10), Table 3 gives the summary of the results for the four different approaches obtained using the cuboidal region -1<xi<1, i=1,2,3. The optimal setting, estimated mean response, estimated standard deviation, and RMSE are presented in Table 3. The RMSE is calculated using the formula RMSE=(ω^μ-T)2+ω^σ21/2, where T=500. It can be seen that the VM approach leads to the optimum setting (x1,x2,x3)=(0.62,0.23,0.1) which resulted in an estimated mean response of (501.57) and root mean squared error of (51.94). This optimal setting is very close to the target mean response but with larger RMSE. The second approach LT produces a target value of (494.69) and RMSE of (44.78) which is smaller than the RMSE of the VM. The third approach is the WMSE which gives the estimated mean response of (496.44) and RMSE of (44.81) which indicate a slight increase in RMSE and little improvement in mean response. The overall performance of the proposed approach is better than those three approaches mentioned with approximate mean value of (500.00) and RMSE of (44.75). Similar procedure is repeated in the next example in order to demonstrate a clear advantage of using the proposed method in terms of closeness to target mean response and the smallest RMSE.Table 2
The printing process study data.
Index
x
1
x
2
x
3
y
u
1
y
u
2
y
u
3
y
-
u
s
u
1
−1
−1
−1
34
10
28
24.0
12.49
2
0
−1
−1
115
116
130
120.4
8.39
3
1
−1
−1
192
186
263
213.7
42.80
4
−1
0
−1
82
88
88
86
3.46
5
0
0
−1
44
178
188
136.7
80.41
6
1
0
−1
322
350
350
340.7
16.17
7
−1
1
−1
141
110
86
112.3
27.57
8
0
1
−1
259
251
259
256.3
4.62
9
1
1
−1
290
280
245
271.7
23.63
10
−1
−1
0
81
81
81
81.0
0.00
11
0
−1
0
90
122
93
101.7
17.67
12
1
−1
0
319
376
376
357.0
32.91
13
−1
0
0
180
180
154
171.3
15.01
14
0
0
0
372
372
372
372.0
0.00
15
1
0
0
541
568
396
501.7
92.5
16
−1
1
0
288
192
312
264.0
63.50
17
0
1
0
432
336
513
427.0
88.61
18
1
1
0
713
725
754
730.7
21.08
19
−1
−1
1
364
99
199
220.7
133.80
20
0
−1
1
232
221
266
239.7
23.46
21
1
−1
1
408
415
443
422.0
18.52
22
−1
0
1
182
233
182
199.0
29.45
23
0
0
1
507
515
434
485.3
44.64
24
1
0
1
846
535
640
673.7
158.20
25
−1
1
1
236
126
168
176.7
55.51
26
0
1
1
660
440
403
501.0
138.90
27
1
1
1
878
991
1161
1010.0
142.50Table 3
Comparison with other methods using printing process study data.
Method
Optimal settings
ω
^
μ
ω
^
σ
RMSE
VM
(0.62, 0.2300, 0.10)
501.57
51.92
51.94
LT
(1.00, 0.07, −0.25)
494.69
44.46
44.78
WMSE
(1.00, 0.08,  −0.25)
496.44
44.67
44.81
PM
(1.04, 0.02,  −0.23)
500.00
44.75
44.75


## 5.2.  The Catapult Study Data

This example will consider the data used by Luner [24] and Kim and Lin [9]. Three variables, x1 (arm length), x2 (stop angle), and x3 (pivot height), are under consideration to predict the distance to the point where a projectile landed from the base of the roman style catapult. The experiment is a central composite design with three replicates as given in Table 5. In [9], the fitted second order polynomial regression models for the mean and standard deviation functions are given by(11)ω^μ=84.88+15.29x1+0.24x2+18.80x3-0.52x12-11.80x22+0.39x32+0.22x1x2+3.60x1x3-4.42x2x3,ω^σ=4.53+1.84x1+4.2x2+3.73x3+1.16x12+4.40x22+0.94x32+1.20x1x2+0.73x1x3+3.49x2x3.Here, the assumed target mean value is T=80. Table 4 shows the estimated mean and RMSE of the estimated optimal mean response. It is evident from Table 4 that the proposed method outperformed the other existing procedures.Table 4
Comparison with other methods using catapult study data.
Method
Optimal settings
ω
^
μ
ω
^
σ
RMSE
VM
(0.09, 0.65, 0.00)
81.13
2.79
3.00
LT
(−0.12, 1.00, 0.39)
79.90
2.59
2.59
WMSE
(−0.14, −1.00, 0.35)
78.68
2.54
2.87
PM
(−0.12,  −1.00,  0.39)
80.00
2.59
2.59Table 5
The catapult study data.
Index
x
1
x
2
x
3
y
u
1
y
u
2
y
u
3
y
-
u
s
u
1
−1
−1
−1
39
42
42
38.3
4.0
2
−1
−1
1
80
91
71
80.7
10
3
−1
1
−1
52
45
44
47
4
4
−1
−1
1
97
60
68
75
19.5
5
1
−1
−1
60
68
53
60.3
7.5
6
1
1
1
113
127
104
114.7
11.6
7
1
1
−1
78
65
64
69
7.8
8
−1.682
0
1
130
75
79
94
30.7
9
1.682
0
0
59
60
51
56.7
4.9
10
0
−1.682
0
115
117
102
111.3
8.1
11
0
1.682
0
50
57
43
50.0
7.0
12
0
0
−1.682
88
43
49
60.0
24.4
13
0
0
1.682
54
60
50
54.7
5.0
14
0
0
0
122
119
109
116.7
6.8
15
0
0
0
87
89
78
84.7
5.9
16
0
0
0
86
85
79
83.3
3.8
17
0
0
0
88
87
81
85.3
3.8
18
0
0
0
89
87
82
86.0
3.6
19
0
0
0
86
88
79
84.3
4.7
20
0
0
0
88
90
79
85.7
5.9


## 6. Conclusion

Numerous procedures have been developed in the literature to obtain an optimal setting condition for the dual response methodology. This paper discusses four different objective functions for dual response optimization approach, based on the mean and variance models as separate response functions. The proposed objective function is based on the penalty function method. We have proposed a new objective function which is more efficient compared with the other existing methods. Numerical examples and simulations study are carried out to compare the performance of the newly proposed method with the frequently used methods. The numerical results clearly show an improvement of the proposed method over the existing methods in terms of having the smallest bias and RMSE. Moreover, the proposed approach can be applied to ridge analysis method and robust parameter design optimization.



---

*Source: 450131-2015-05-12.xml*