Dynamic Robust Parameter Design Using Response Surface Methodology based on Generalized Linear Model

Purpose: When designing an input-output system susceptible to noise, engineers assume a functional relation between the input and the output. The Taguchi method, which uses a dynamic, robust parameter design (RPD) to evaluate the robustness of the input-output relation against noise, is employed. This study aims to address extending the scope of use of a dynamic RPD. Methodology/Approach: A target system in a typical dynamic RPD can be interpreted as one in which the relation between the input and the output is a linear model, and the output error follows a normal distribution. However, an actual system often does not conform to this premise. Therefore, we propose a new analysis approach that can realize a more flexible system design by applying a response surface methodology (RSM) based on a generalized linear model (GLM) to dynamic RPD. Findings: The results demonstrate that 1) a robust solution can be obtained using the proposed method even for a typical dynamic RPD system or an actual system, and 2) the target function can be evaluated using an adjustment parameter. Research Limitation/implication: Further analysis is required to determine which factor(s) in the estimated process model largely contribute(s) to changes in the adjustment parameter. Originality/Value of paper: The applicability of typical dynamic RPD is limited. Hence, this study’s analytical process provides engineers with greater design flexibility and deeper insights into dynamic systems across various contexts.


INTRODUCTION
When designing an input-output system, engineers assume a functional relation between the input and the output (Tatebayashi, 2004).However, in practice, this functional relation does not consistently hold because of noise, which results in critical quality defects.Therefore, the Taguchi method uses a dynamic, robust parameter design (RPD), an indispensable design method, for evaluating the robustness of the input-output relation against noise.The method is centred around reducing variations in a process, and it accounts for noise at the design stage.
In a typical dynamic RPD, the relation between the input  and output  of a system is assumed to be linear, while the output error  is assumed to follow a normal distribution and possess equal variability.Based on Nagata's (2009) explanation, this assumption can be defined using the following additive model:  =  + , ~(0,  2 ). (1) Hence, methods derived from the RPD (e.g., Kawamura and Takahashi, 2013) also assume that the data follows these assumptions.However, actual systems often do not conform to these assumptions.For example, Mikami and Yano (2004) employed a typical dynamic signal-to-noise ratio analysis with the number of thermotolerant bacteria as the output and the incubation time as the input.
However, the number of bacteria per unit of time generally followed a Poisson distribution.Moreover, analyses such as the growth rate of bean sprouts (Yoshino, 1995) assume a growth curve where the relation between the input and output of the system is nonlinear.
To realize a more flexible system design, we focus on a generalized linear model (GLM), a nonlinear model with a nonlinear input-output relation.It is useful for analyzing experimental data because it can handle output errors that follow a nonnormal distribution with unequal variability (Lee and Nelder, 1998).If the typical dynamic RPD is reconsidered in the GLM context, the linear predictor then becomes  , the link function is absent, and the error structure is normally distributed.
Therefore, we propose a new analytical approach that can realize a more flexible system design by applying response surface methodology (RSM) based on GLM to dynamic RPD.The formulation of the proposed method is based on the framework by Myers et al. (2005), in which an approach to RPD was developed to verify its usefulness.In this study, we demonstrated that our proposed method can enable engineers to attain greater design freedom and gain insights into the experimental data of dynamic systems with various backgrounds.
The remainder of this paper is organized as follows.Section 2 explains the approach RPD proposed by Myers et al. (2005).Section 3 proposes a new approach to RPD for dynamic systems.Section 4 analyzes actual data to verify the performance of the proposed method.Section 5 presents the verification results.
Section 6 summarizes the study and provides future recommendations.
2 STATIC RPD USING RSM BASED ON GLM Myers et al. (2005) proposed an approach for static RPD using RSM based on GLM for a static system-a system with fixed inputs.The method aimed to design the system to ensure the output is constantly maintained at the target value.They considered the conditional population mean   of response   , as follows: ) ,   ~(0,   2   ). (2) Provided the control factors   and noise factors   use the GLM framework,  =  and  = (, ) is the linear predictor.The error structure is chosen to fit the data and the link function (•) is chosen accordingly.The index  = 1,2, … ,  denotes the combination of control factors, and  = 1,2, … ,  denotes the combination of noise factors. 0 is the intercept parameter.The  dimensional vector  denotes the vector of coefficients for the control factors, while the  dimensional vector  represents the vector of coefficients for the noise factors.The  ×  matrix  denotes the matrix of control by noise interaction coefficients, and the  ×  matrix  represents the matrix of second order effect coefficients of the control factors.Each level of the noise factor follows a normal distribution (0,   2 ).
Subsequently, the process mean   (  ) and process variance (  ) were derived.Using the target value  of the response, the following evaluation function was defined as In the RPD proposed by Myers et al. (2005), the design solution is a combination of the control factors that minimize the evaluation function.

Proposal of a new method
In this section, we formulate an approach for dynamic systems based on Section 2. This method comprises four steps.In step 1, the RSM is estimated using a GLM or double generalized linear model (DGLM).The GLM was proposed by Nelder and Wedderburn (1972) and holds the dispersion parameter as a constant.The DGLM, an extension of the GLM proposed by Smyth (1989), also models the dispersion parameter by enabling additional flexible modelling.
We consider the conditional population mean   of the response   , as follows: where  is the inverse of the link function .
Second, using a second order Taylor series approximation around the mean of the linear predictor  0 =   (), the process mean is given as: The variance of the linear predictor   () is defined as In step 3, we derive the process variance (  ) using Lee and Nelder's (2003) framework.(  ) is given as: The first term shows the variation due to noise factors, while the second term expresses the variation independent of the noise factors.  (dispersion parameter) represents the variation in population mean   , and the variance in independent distribution, indicating that it is unaffected by the increase or decrease in the process mean.ior simplicity, we consider the dispersion parameter as   , which remains constant regardless of noise factors.(  ) denotes the change in variance relative to the population mean   of the distribution.This indicates that the variation is affected by the increase or decrease in the process mean.
When estimating the RSM with GLM,   is fixed to a constant.When the RSM is estimated using DGLM, we consider the dispersion parameter for the conditional population mean   of the response   , as follows: given   and   .The index  = 1,2, … ,  denotes the combination of control factors, and  = 1,2, … ,  represents the level of the signal factor.Based on Smyth and Verbyla's (1999) explanation,   refers to (, ) obtained by transforming the exponential family distribution as shown in is a distance measure between  and . 0 is the intercept parameter,   denotes the control factor and   represents the function of the signal factor.The  dimensional vector  represents the vector of coefficients for the control factors while the  dimensional vector  denotes that for functions of the signal factor.
The  ×  matrix  denotes the matrix of the control by function of the signal interaction coefficients, and the link function is ℎ(•) .Notably, the variables, vectors, and matrices of the coefficient that constitute   differ from those of   .We adopt the estimation method proposed by Smyth and Verbyla (1999).
iirst, we derive the first term of the process variance, represented in This equation can be approximated as: using the delta method.The latter can be represented as: using the chain rule.Ψ =   ⁄ and  are location parameters.
Second, we derive the second term of the process variance using a second order Taylor series approximation around  0 , given as: Finally, (  ) is defined as In step 4, the parameters are optimized.Hence, we define the necessary evaluation function as the adjusted sum of square errors ( ) for optimization, as follows: M is the set of signal factor levels. denotes the target function ( 0 * +  * ′ ) and  represents the adjustment parameter (0 ≤ ) .The first term measures the deviation between  and   () while the second measures the process variation. is used to make adjustments between the target function and the process mean.Hence, the  2 scale of variance is omitted.This adjustment allows for the process mean to approach the target function, which is the design concept.If the target function is less than the estimated process mean, then  is less than 1; otherwise,  is greater than 1.Thus, the target function can be evaluated.
Next, we explain the optimization method.iirst, we minimize the   with an adjustment parameter of one to derive the design solution.We use the design solution and  = 1 to calculate the first term of the  , which is the evaluation value.Second, we derive the design solution by generating the adjustment parameter and minimizing the   to which it is assigned.We use the design solution and the generated adjustment parameter to calculate the first term of the  .If the first term is less than the evaluation value, the generated adjustment parameter is adopted and the evaluation value is updated.
Otherwise, a new adjustment parameter is generated.This series of processes is repeated until the adjustment parameter yielding the lowest evaluation value is derived.iinally, we find the combination of control factors that minimize the   substitution using the adjustment parameter identified in the previous step.The   is minimized by determining an adjustment parameter using the simulated annealing method and deriving the design solution using the quasi-Newton method.
The typical dynamic RPD is a two stage design that first decreases variance and subsequently enables the mean to approach the target function (e.g., dynamic signal-to-noise ratio analysis).Contrastingly, our method first draws the mean closer to the target function and subsequently decreases the variance.This departure arises because the system under consideration entails a tradeoff between mean and variance, a complexity beyond the scope of conventional two stage design.

Related research
Kume and Nagata (2013) designed the parameters by defining and minimizing the following equation: where  denotes an arbitrary weight to balance the first and second terms (0 ≤  ≤ 1).If   (  ) needs to be adjusted to ,  should be set to a higher value; otherwise, it should be set to a lower value.However, a robust solution cannot be obtained without setting an appropriate target function for the estimate process model.This is particularly true when the difference between the process mean and target function is significant.Additionally, the arbitrary nature of the decision of  presents another challenge.

Overview of data
In this section, we verify the design solution using   and actual data.The experimental data of the high-speed response valve (also known as a highspeed on/off solenoid valve) employed by Enkawa and Miyakawa (1992) is used along with dynamic systems.In a high-speed response valve, the input-output relation for the generic function is given by the following equation: The input (signal factor) is the square-root transformation of the pressure difference, the output (response ) represents the flow rate, and the duty ratio denotes the ratio of the time the valve is on to the times it is on and off.In this experiment, four control factors are assigned to an inner  8 orthogonal array. is the stroke,  is the spring-attachment load,  is the pressure balance, and  is the oil passage area.The noise factor  and signal factor  are arranged in a one-by-one outer array.In the following, the notation of the noise factor is , and that of the signal factor is .The noise factor is the input voltage.The first level of the control factor and noise factor is 1, and the second level is −1.The signal factor is 4 for the first level, 8 for the second level, and 12 for the third level.
These data follow the typical dynamic RPD assumptions.Therefore, it is ideal that the design solution using   is consistent with their conclusions-the good design solution is to set (, ) = (1, −1) to reduce variation in the inputoutput relation and  to adjust the output.

Data analysis
First, we estimate the response surface using  of the link function and DGLM in the normal distribution as follows: Nair (1992) demonstrated that when the   error structure is normally distributed, the dispersion parameter model follows a gamma distribution with log() of the link function.Therefore, log() is used as the link function in the dispersion parameter model, which assumes a gamma distribution.Our analysis uses the mean square error calculated using leave-one-out cross-validation.Because control factor  was not selected, setting the level is unnecessary.
In the population mean model, the control factors  and  interact with the noise factors and the control factors  and  are selected in the dispersion parameter model.Therefore, a tradeoff occurs because the control factors decrease the variation while simultaneously adjusting to the target function.
We derive , (), and canonical link functions in the normal distribution.The normal distribution (,  2 ) can be transformed following Myers and Montgomery (1997) to obtain the following equation: Hence, in the normal distribution,  =  , ℎ() =  2 , and () = 1 .The canonical link function is .
Second, we derive the process mean   (  ).Considering that the link function is  and  ′′ [ 0 ] = 0, we can derive the following equation: When the population mean model of Equation ( 19) is substituted, the estimation model of the process mean is given as: Third, the process variance, (  ), is derived.Notably, Ψ = 1 because we use the established link functions and () = 1 because we use normal distributions.Therefore, the following equation is derived: (  ) ≈ ( +     +    )  (  )( +     +    )  +   . ( When the population mean model and dispersion parameter model of Equation ( 19) are substituted, the estimation model of the process variance is given as: ̂(  ) = (−10.92+ 4.50 − 4.23) 2   2 + exp(4.25 − 0.73 − 0.92). (24) Fourth, we define the evaluation function that requires parameters that are arbitrarily determined by the designer based on the purpose and state of the process.The setting used for this analysis is as follows: = 1.0,  = (4.0,6.0,8.0,10.0,12.0), where  denotes the target function,  denotes the standard set of signal factors, and   2 denotes the variance of the noise factors.Because an appropriate target function was unknown in the original experiment (Enkawa and Miyakawa, 1992), we define a new one.We define each evaluation function based on the setting and optimize in the range of −1 ≤  ≤ 1 for each evaluation function.The equation for the evaluation function is omitted.
The details of each design solution are listed in Table 2.  is the evaluation function with the normal RSM approach when the   is fixed at  = 1.The  is minimized using the quasi-Newton method.
As presented in Table 2, the design solution using the   is consistent with the conclusions of Enkawa and Miyakawa (1992).The first term of the   is diminutive owing to the adjustment parameter.Therefore, we calculate and compare the first term of the  using the design solution of the  .It can be observed that the  is smaller, indicating that it derives the design parameters by drawing the process mean closer to the target function.Moreover, because the adjustment parameter is 0.17, the target function is less than the estimated process model.The second term allows for a simple comparison because no adjustment exists in either evaluation function.By comparing the values of the second term, it is evident that the   derives solutions that have a more attenuated process variation than the .Next, we visually evaluate the quality of each design solution.We show the relation plot between the level of the signal factor and the response in the obtained design solution in iigure 1.The vertical axis of the graph shows the response  and the horizontal axis shows the level of signal factor .The target function t is indicated by a solid line; the process mean   (  ) is indicated by a dotted line, and the population mean   of the response at each level of the noise factor  is indicated by four white points.iurthermore, the level of the signal factor  is plotted in 1.0 increments 4.0 of 12.0.The noise factor  has four levels (−1.0, −0.5, 0.5, 1.0).
As illustrated in Figure 1, the solution designed using  gets closer to the target function by increasing the process variation.Conversely, the solution designed using   decreases the process variation while drawing the process mean sufficiently closer to the target function.

Simulation settings
This section describes the simulations performed to verify  .The differences and variations between the ideal and design solutions in the simulation model are verified as evaluation criteria to clarify whether our solution is robust.
We set the simulation model for response , as follows: The error  follows a Poisson distribution () and  = .Hence, this simulation does not conform to the typical dynamic RPD assumption.The conditional mean and variance are expressed as (  |  ,   ,   ) =   .
Moreover, our simulation does not assume over-or under-dispersion.Therefore, we use the GLM to estimate the RSM.
The solution (, , , ) = (1.0,1.0, −0.5, −0.5) is considered ideal because it maximizes the intercept and slope while nullifying the effect of the noise factor.However, extrapolated solutions are not assumed.Substituting the ideal solution into Equation (26) yields the following: Therefore, we set the target function  as in log() = 2.0 + 1.0.
The slope of this target function exceeds that in Equation ( 29).If the design solution increases the process variation, the process mean can be brought closer to the target function.However, this design solution is not robust.Therefore, in this simulation, the average of the design solutions is close to (, , , ) = (1.0,1.0, −0.5, −0.5); the smaller the variance in the design solutions, the better the solution.
In the framework in GLM, where neither overdispersion nor underdispersion is assumed, the Poisson distribution () can be transformed following Myers and Montgomery (1997) Each parameter used for the design is given in the following equation:  = (1.0,1.5,2.0,2.5,3.0), The equation of the evaluation function is omitted.
The data used for the simulation were generated using the orthogonal array  16 .Factors , , , and  are assigned to 1 st , 2 nd , 4 th , and 8 th columns, respectively.The noise factor  and signal factor  are arranged in a one-by-one outer array.
The first level of the noise factor is 1, while the second is -1.The first, second, and third levels of the signal factor are 1, 2, and 3, respectively.The number of simulations is 10,000.In the real data analysis, the  -a normal RSM approach-is also used for comparison.

Simulation results
Figure 2 shows the simulation results.The red cross symbol in the box plot represents the ideal solution level, and the blue point is the average of each design solution.
As indicated in Figure 2, the control factors  and  optimized using  are closer to the ideal solution, whereas the control factors  and  are significantly far from the ideal solution.In the simulation model, the design solutions other than  are related to the noise factor.These results imply that the process variation increases.Therefore, it can be inferred that the process mean can approach the target function by increasing the process variation without decreasing the influence of the noise factor.
Conversely, the control factor optimized using the   is very close to the ideal solution, and the variability of all control factors is diminutive, indicating that the variation is reduced while achieving appropriate input-output relations.
Based on the above, in the simulation of a realistic system, the   can obtain a more robust stable design solution against changes in the noise factor and draw the process mean closer to the target function than the .

CONCLUSION AND FUTURE ISSUES
The dynamic RPD is indispensable for evaluating the robustness of the inputoutput relation.However, realistic systems often deviate from typical dynamic RPD assumptions.Therefore, this study proposed a novel approach to RPD using RSM based on GLM for dynamic systems based on the method proposed by Myers et al. (2005).Additionally, we demonstrated the effectiveness of our method using dynamic experimental data and simulations.
The actual data analysis reveals that a system design using   is found to derive a more robust design solution against the fluctuation of the noise factor than the one using .Moreover, the target function could be evaluated using the adjustment parameter.In the simulation, when the system design using   sets a free target function, we could derive a design solution close to the assumed ideal solution.This enables the designer to freely set the target function, design the system with various backgrounds, and subsequently evaluate it.
Future research should analyze the factor(s) in the estimate process model that significantly contributes to the adjustment parameter's variation.We propose introducing the contribution rate as an index to identify these factors.Thus, significant improvements are expected in the process adjustments.

Figure 1 -
Figure 1 -Relation plots between signal factor levels and responses in each design solution

Figure 2 -
Figure 2 -Boxplot of each design solution derived in simulation denotes the level of the signal factor. 0 is the intercept parameter,   denotes the control factors,   represents the noise factors and   denotes the functions of the signal factor.The  dimensional vector  signifies the vector of coefficients for the control factors, the  dimensional vector  denotes that for the noise factors, and the  dimensional vector  denotes that for functions of the signal factor.The  ×  matrix  represents the matrix of control by noise interaction coefficients, the  ×  matrix  is the matrix of the control by function of the signal interaction coefficients, and the  ×  matrix  represents the matrix of the noise by function of the signal interaction coefficients.