On the Application of a Design of Experiments along with an ANFIS and a Desirability Function to Model Response Variables

Abstract

In manufacturing engineering, it is common to use both symmetrical and asymmetrical factorial designs along with regression techniques to model technological response variables, since the in-advance prediction of their behavior is of great importance to determine the levels of variation that lead to optimal response values to be obtained. For this purpose, regression techniques based on the response surface method combined with a desirability function for multi-objective optimization are commonly employed, since it is usual to find manufacturing processes that require simultaneous optimization of several variables, which exhibit in many cases an opposite behavior. However, these regression models are sometimes not accurate enough to predict the behavior of these response variables, especially when they have significant non-linearities. To deal with this drawback, soft computing techniques are very effective in overcoming the limitations of conventional regression models. This present study is focused on the employment of a symmetrical design of experiments along with a new desirability function, which is proposed in this study, and with soft computing techniques based on fuzzy logic. It will be shown that more accurate results than those obtained from regression techniques are obtained. Moreover, this new desirability function is analyzed in this study.

1. Introduction

In manufacturing engineering, it is necessary to model the behavior of response variables as a function of several input variables in order to be able to select the most suitable operating conditions in advance. However, this selection of optimal operating conditions sometimes conflicts with the desired levels of variation for the response variables, since the levels of the independent variables that can optimize a given output (for example, the surface roughness or the wear of a tool) can conflict with those that optimize another different one (for example, the material removal rate). Therefore, the optimal selection of operating conditions is not an easy task and it is necessary to use statistical tools for modeling of the response variables versus variations in the process parameters. This can sometimes lead to a high number of experiments being carried out, and therefore it is common practice to employ design of experiments techniques that allow, with a relatively small number of experiments, some very useful information on the behavior of the variables under study to be obtained. In this way, response surface method (RSM), artificial neuronal networks (ANN), fuzzy inference systems (FIS) [1,2], and adaptive neural fuzzy inference system (ANFIS) [3] have been widely employed. With regard to fuzzy inference systems, Takagi-Sugeno [1] and Mamdani [2] are the most commonly employed. Several studies can be found in the literature dealing with the application of these aforementioned techniques [4,5,6].

Fuzzy inference systems (FIS) have been employed in several scientific fields, for example, dealing with modeling of material removal in abrasive belt grinding processes, as shown in Pandiyan et al. [7], classification of defects, as shown in Versaci [8] and in Burrascano et al. [9], to improve talc pellet manufacturing processes [10] and to predict the compressive strength of cement [11]. Further examples can be found in [12,13,14], among many others. In order to use the above-mentioned models, it is quite common to employ symmetric designs, in which all the factors have the same levels of variation as well as asymmetric ones, where not all the factors have the same levels of variation [15,16]. In manufacturing engineering 2^k, 3^k, and 4^k symmetric designs are commonly used, because a large amount of information is available from a relatively small number of experiments. As a consequence of the technological interest of modeling output variables in manufacturing processes, a large number of models have been developed in recent years based both on conventional regression techniques and on the use of soft computing, both through fuzzy inference systems and through neural networks, generally by using supervised learning and feed forward networks. It should be underlined that when regression models have high values in the coefficients of determination, they may approach the behavior of response variables with relatively high accuracy. However, when the response variables exhibit high non-linearities, conventional regression models are not capable of modeling the behavior of output variables adequately and hence some other types of techniques need to be employed.

Since it was first proposed by Harrington [17], and later modified by the classification function of Derringer and Suich [18], the methodology of multiple optimization based on the use of a desirability function, has been used in a large number of research studies for the simultaneous optimization of response variables. In this present study a symmetrical 4³ factorial design of experiments combined with an ANFIS and a new desirability function is employed. The proposed function to transform the output responses is of an arctangent type, which is analyzed in this present study. The ANFIS along with this new desirability function will be used to determine the values that can approach the optimum value of the independent variables that simultaneously optimize the response outputs.

As can be observed in Figure 1, three independent variables x₁, x₂, and x₃ are varied in order to show the values that optimize a desirability function by using a symmetric design of experiments and an ANFIS, following the methodology shown in this present study.

However, as previously mentioned if regression models are not capable of adequately predicting the actual behavior of output variables some other techniques should be used. Moreover, results obtained in this present study will be compared with those obtained with the transformation proposed by Derringer and Suich [18], which is commonly employed by several research studies found in the literature, as will be shown in the Section 3.

2. State of the Art

The research study of Harrington [17] showed how to optimize several variables simultaneously by using a desirability function combined with an exponential classification function. Subsequently, Derringer and Suich [18] proposed a transformation for the optimization of multiple variables simultaneously using a desirability function and assuming that the estimated values are a continuous function of the dependent variables and therefore the desirability function, which they obtain from the geometric mean of the transformed variables in the interval [0,1], is a continuous function of the independent variables [18]. As is stated in the research study of Del Castillo et al. [19], the desirability function approach consists in using response surface methods to fit polynomials for each response variable in order to have a single optimization problem. In their study Del Castillo et al. [19] present modified desirability functions that are differentiable everywhere. This is important because it is quite common for the obtained overall desirability functions obtained from RSM to be non-differentiable, especially when there are several responses and design variables [20]. In the study of He and Zhu [20], the authors propose a hybrid approach, which combines the genetic algorithm with the pattern search method, to deal with this point. In Costa el al. [21] a review on different desirability methods is carried out. Some other authors have proposed several transformation functions. For example, an exponential transformation is suggested by Kim and Lin [22]. In another study, Kros and Mastrangello [23] study the relationship between response types when they are mixed (the larger the better, the smaller the better, and nominal the best). On the other hand, Das et al. [24] employ a genetic algorithm-based optimization to determine the tensile properties of a high-strength low-alloy steel, where multi-layer perceptron ANNs are developed for the output responses and then a desirability function based on employing a sigmoidal function is obtained [24] and Zong et al. [25] combine variations due to noise factors and controllable factors. In another study developed by He et al. [26], a robust desirability function approach to simultaneously optimize multiple responses is proposed. Some other studies are those of Ribardo and Allen [27] that propose a method, which is expressed as a function of the mean and standard deviation of the associated criteria or quality characteristics; that of Wu and Chyu [28] which present an approach to optimizing correlated multiple quality characteristics with asymmetric loss function, that of Ortiz et al. [29] which proposes a multiple-response solution technique using a genetic algorithm in conjunction with an unconstrained desirability function and that of Pasandideh and Niaki [30] which models a multiple variable optimization problem by using the desirability function approach and then they employ four genetic algorithm methods to solve the problem.

On the other hand, Das [31] make use of the property of desirability function in the neural network architecture and evaluate their performances. Three combinations of transfer/desirability functions are analyzed in this study [31]. Lee and Kim [32] proposed a desirability function defined as the average of the conventional desirability values based on the probability distribution of the predicted response variable and in He et al. [33] a robust methodology for response surface analysis, when there are several responses, is studied. Moreover, in their study the authors propose a measure of robustness for second order response surface models [33]. In Lee et al. [34], a desirability function method to simultaneously optimize both the mean and variability of multiple responses is proposed.

Although the desirability function is a very interesting method to simultaneously optimize several variables, the transformation into the desirability space should take the preferences of the decision maker into account in order to obtain the optimal parameters as stated by Fuller and Scherer [35].

Derringer and Suich [18] transformation has been widely employed over the last few years. For example, the use of models based on the use of design of experiments along with the response surface methodology using the Derringer and Suich methodology has been employed in several research studies as shown in Vera et al. [36]. In the research study of Padilla et al. [37] the authors carry out an analysis on the optimization of a spark ignition engine. A 3³ design of experiments (DOE) is conducted where the independent variables are the revolutions, the load produced by a dynamometer, and an electrolyte concentration and the response variables are torque, hydrocarbon emissions, and power. In another study, Hur et al. [38] present the results of surface response optimization using the desirability function with a central composite design to determine the design variables to control the distribution of flow rate in a cooling system employed in an electric vehicle. On the other hand, Saleem et al. [39] consider the deflection angle, input power, and micromirror temperature rise from the ambient as response variables to be simultaneously optimized by using response surface models and a desirability function. The application of design of experiments, response surface methodology and a desirability function is also discussed by Akçay and Anagün [40] to model and optimize a manufacturing process. The influence of cutting speed, feed rate, and depth of cut on the machinability characteristics of Si optical lenses is analyzed by Jumare et al. [41]. These authors develop predictive models in order to obtain minimum arithmetical mean roughness and tool wear as well as maximum material removal rate using the desirability function approach [41].

Other studies such as that of Chahal et al. [42] show that RSM along with a desirability function is effective in order to optimize surface roughness in end milling; Zhao et al. [43] analyze the effect of welding parameters on the mechanical properties of welding bead and welding heat input using a DOE and a desirability function and Kribes et al. [44] study the effect of machining conditions on the obtained surface roughness in hard turning by using response surface methodology and a desirability function. Some other studies which employ the response surface methodology and a desirability function are Ahmad et al. [45] that optimizes the recycling conditions of aluminum (AA6061) chips, Qazi et al. [46] that analyzes the machinability of AA5005H34, and Osman et al. [47] that studies continuous multi-pass friction stir welding to clad dissimilar materials. Cutting speed, feed rate, and depth of cut are used as design parameters along with RSM and desirability function in Laghari et al. [48] for modeling cutting forces in the machining of a composite material; Pradhan et al. [49] analyze surface roughness in machining of Al/SiCp metal matrix composite; optimization of friction surface deposition of stainless steel over medium carbon steel is shown in Sahoo et al. [50] and in Du et al. [51] a multi-objective optimization is carried out in the dry machining of a 304 stainless steel, among many others.

The number of studies that exist in the literature related to the combined application of soft computing techniques and desirability functions is much lower. Among these studies, it is worth mentioning that of Mostafaei [52] that employs ANFIS models for prediction of biodiesel fuels cetane number using a desirability function, which employs a linear transformation and a comparison of the developed models is performed by statistical criteria. Another study is that of Labidi et al. [53] where ANN, RSM, and desirability function are used to determine optimal machining conditions in the turning of a hardened steel. On the other hand, Sengottuvel et al. [54] study the effects of various electrical discharge machining input parameters as well as the influence of different tool geometry on material removal rate, tool wear rate, and surface roughness. These authors also develop a fuzzy inference system which is validated with experimental results. Another study worth mentioning is that of Singh et al. [55] where the application of Taguchi’s robust design coupled with fuzzy based desirability function approach for optimizing multiple bead geometry parameters of submerged arc weldment is developed. Response surface methodology and the desirability function results were also employed for modeling Ti-6Al-4V milling under different lubrication conditions in the study of Paschoalinoto et al. [56]. These authors also compared the statistical results with those obtained from using supervised artificial neural networks and they found that both artificial neural network and the experimental design predict similar results [56]. Tank et al. [57] use a desirability function for each of the surface roughness parameters analyzed in the turning of glass fiber reinforced with plastic and then employ a fuzzy inference system to determine values of a characteristic index which then are used to find the signal to noise ratio. In the research study of Salmasnia et al. [58] an approach for optimization of correlated multi-response variables is analyzed. In their study an ANFIS is used as a tool for predicting system behavior. On the other hand, the study of Gajera et al. [59] deals with the application of response surface methodology with fuzzy-based desirability function approach in order to optimize multiple process parameters in direct metal laser sintering, among others.

3. Methodology

As shown in the review of the state of the art, statistical tools based on the design of experiments and the response surface method have been widely used to model the behavior of different response variables and it has also been shown that the use of a desirability function based on RSM and the model proposed by Derringer and Suich [18] have been commonly used in order to simultaneously optimize response variables, since it is usual to find manufacturing processes that require simultaneous optimization of several variables, which exhibit, in many cases, an opposite behavior. In manufacturing processes, is it quite common for certain technological variables to have a more regular response behavior and, therefore, they can be adequately modeled by methodologies such as RSM, and, on the other hand, there also exist some other technological variables that cannot be adequately modeled using conventional regression. These regression models are sometimes not accurate enough to predict the behavior of these response variables, especially when they have significant non-linearities.

In order to address these problems, in this present research study, the use of a new desirability function which is based on the inverse tangent (arctangent) is first proposed to classify the output responses of a design of experiments, so that for each of the response functions $f_1,\dots , f_n$ a transformed function that lies between [0,1] will be obtained from Equation (1).

$$y_j\stackrel{transformed}{\to }{y}_{tj}=c_1+c_2\ast ta{n}^{-1}\left(a_j\ast y_j+b_j\right)$$

(1)

where $c_1$ and $c_2$ are first adjusted so that the arctangent function lies between 0 and 1. This is easily obtained from simultaneously solving for $c_1$ and $c_2$ the equations $\left[1=c_1+c_2\frac{\pi }{2}, 0=c_1-c_2\frac{\pi }{2}\right]$ which yields to $\left[c_1=\frac{1}{2}, c_2=\frac{1}{\pi }\right]$. Therefore, the transformation to be used in this present study is given by Equation (2).

$$y_j\stackrel{transformed}{\to }{y}_{tj}=\frac{1}{2}+\frac{1}{\pi }\ast ta{n}^{-1}\left(a_j\ast y_j+b_j\right)$$

(2)

The values of constants $a_j$ and $b_j$ will be selected depending on the desired level for each variable (maximize, minimize, or lying within a range of values).

In the event that the desired variation of an output response ($y_j)$ is its maximum value, the constants of Equation (2) will be given by Equations (3) and (4).

$$a_j=\frac{1}{M_j-c_{M_j}}\left(H_{M_j}-tan\left(\pi (t_{c_{M_j}}-\frac{1}{2} \right)\right)$$

(3)

$$b_j=\frac{-1}{M_j-c_{M_j}}\left(H_{M_j}\ast c_{M_j}-M_j\ast tan\left(\pi (t_{c_{M_j}}-\frac{1}{2} \right)\right)$$

(4)

where $c_{M_j}$ for the case of a variable whose preferred level is the maximum ($M_j$) is defined from Equation (5), where $k_{c_{M_j}}$ is a constant which can be used to set $c_{M_j}$at a specific point (default value $k_{c_{M_j}}=1)$ and $t_{c_{M_j}}$ is defined to fix the value of the arctangent at this point. Moreover, the value of the constant $H_{M_j}$ is selected, for simplicity, to be a multiple of the difference between the maximum $\left(M_j\right)$ and the minimum $(m_j)$. Then $H_{M_j}$ is set by using the constant $k_{H_{M_j}}$, as Equation (6) shows (however, this value could be selected in whatever different way) and it is selected in order to set the transition of the output variable $(y_{tj})$ to the level of “1”. As previously mentioned, $M_j and m_j$ are the maximum and the minimum levels of each variable. In this present study these values will be selected from the DOE, however they could be whatever.

$$c_{M_j}=\frac{m_j+k_{c_{M_j}}\ast M_j}{1+k_{c_{M_j}}}$$

(5)

$$H_{M_j}=k_{H_{M_j}}\ast \left(M_j-m_j\right)$$

(6)

Analogously, if the preferred level of an output variable$\left(y_j\right)$ is the minimum, the $a_j$, $b_j$ and $c_j$ constants as well as $H_{m_j}$ are obtained from Equations (7)–(10) respectively, so that if the output variable is close to its minimum value ($m_j$), then the transformed function $(y_{tj})$ will approach to “1” and, on the contrary, if the value of this variable approach the maximum, in this case non-desirable value, then the transformed output will correspond to a value close to “0”, where the sub-index “t” in the transformed variable $(y_{tj})$ refers to an output variable transformed by using the arctangent function, as previously shown.

$$a_j=\frac{1}{m_j-c_{m_j}}\left(H_{m_j}-tan\left(\pi (t_{c_{m_j}}-\frac{1}{2} \right)\right)$$

(7)

$$b_j=\frac{-1}{m_j-c_{m_j}}\left(H_{m_j}\ast c_{m_j}-m_j\ast tan\left(\pi (t_{c_{m_j}}-\frac{1}{2} \right)\right)$$

(8)

$$c_{m_j}=\frac{k_{c_{m_j}}\ast m_j+M_j}{k_{c_{m_j}}+1}$$

(9)

$$H_{m_j}=k_{H_{m_j}}\ast \left(M_j-m_j\right)$$

(10)

As previously mentioned, ($c_{M_j} and c_{m_j})$, are defined from Equation (5), when the preferred level of a variable is its maximum, and from Equation (9), when the preferred level is its minimum. Moreover, a weighted average can be used to set ($c_{M_j} or c_{m_j})$ at a specific point by using $(k_{c_{M_j}} or k_{c_{m_j}})$ constants, respectively. The influence of these constants is shown in Figure 2. The $H$ parameter $\left(H_{M_j} or H_{m_j}\right)$ is used to select the transition to the “1” value of the transformed variable. For a given value of $k_{H_{M_j}} or k_{H_{m_j}}$, the values of $t_{c_{M_j}} or t_{c_{M_j}}$ will be selected within the interval given by Equation (11) and by Equation (12), respectively. The influence of $t_{c_{M_j}}$ and $t_{c_{m_j}}$ is shown in Figure 3.

$$0\le t_{c_{M_j}}<\frac{1}{2}+\frac{1}{\pi }\ast ta{n}^{-1}\left(k_{H_{M_j}}\ast \left(M_j-m_j\right)\right)$$

(11)

$$0\le t_{c_{m_j}}<\frac{1}{2}+\frac{1}{\pi }\ast ta{n}^{-1}\left(k_{H_{m_j}}\ast \left(M_j-m_j\right)\right)$$

(12)

$M_j and m_j$ are the maximum and the minimum levels of each variable, as previously mentioned. If values of $t_{c_{M_j}}$ and $t_{c_{m_j}}$ are selected within the range given by

$$\frac{1}{2}+\frac{1}{\pi }\ast ta{n}^{-1}\left(k_{H_{M_j}}\ast \left(M_j-m_j\right)\right)\le \{t_{c_{M_j}}or t_{c_{m_j}}\}\le 1$$

The shape of the function will be inverted. It could have been possible to select values higher than 1 but it makes no sense for the present parameter because its variation is desired to be within the interval [0, 1]. Finally, the influence of the $H$-parameter $\left(H_{M_j} or H_{m_j}\right)$ is shown in Figure 4. If a given value of the t-parameter ($t_{c_{M_j}} or t_{c_{M_j}})$is given, then the value of $k_{H_{M_j}} or k_{H_{m_j}}$, which will define the value of $H_{M_j}$ or $H_{m_j}$, respectively should be selected within the interval given by Equation (13) and by Equation (14), respectively.

$$k_{H_{M_j}}>\frac{1}{M_j-m_j}tan\left(\pi \left(t_{c_{M_j}}-\frac{1}{2}\right)\right)$$

(13)

$$k_{H_{m_j}}>\frac{1}{M_j-m_j}tan\left(\pi \left(t_{c_{m_j}}-\frac{1}{2}\right)\right)$$

(14)

Therefore, the shape of the arctangent transformation can be adjusted by varying ($H_{M_j}$,$H_{m_j})$, ($t_{c_{M_j}},t_{c_{m_j}})$and ($k_{H_{M_j}},k_{H_{m_j}})$ parameters whose influence is respectively shown in Figure 2a, Figure 3a, and Figure 4a, in the case where the preferred level of a response variable is the maximum, and in Figure 2b, Figure 3b, and Figure 4c, in the case where the preferred level of a response variable is the minimum.

If the output response has to be within a specified range [LB,UB], then it is possible to combine two transformation functions similar to those shown in Figure 5a,b, so that the result attained is shown in Figure 6, where LB and UB are both the lower bound and the upper bound selected.

The shape of the transformed function shown in Figure 6 is obtained from Equation (15), where the midpoint (LB + UB)/2 has been considered. However, this point could be whatever within the specified range [LB,UB]

$$y_{tj}=\left|1-\left(y_{t1}+y_{t2}\right)\right|$$

(15)

However, if the output response has to be within the range defined by [LB, UB] and the external values are not acceptable, it is possible to set a value of $t_{c_{M_j}}and t_{c_{m_k}}$ close to zero, so that results shown Figure 7 will be obtained and then by applying Equation (15) the output shown Figure 8, is obtained.

As can be observed in Figure 8, the variable remains within the specified range. As in the previous case, the shape of the transformed function shown in Figure 8 is obtained from Equation (15).

Note that it could be possible to employ this arctangent function by selecting the levels of the transition to “1” and “0” instead of the method shown below by employing two H-parameter (one to adjust the transition to the “1” level and another to the “0” level). In this case, the obtained values will be those shown in Figure 9. However, parametrization shown in Figure 9 will not be used in this present study and the methodology previously shown will be used instead.

Once the response variables have been transformed, a mean is selected to be maximized. As was previously shown in the state-of-the-art Section, most of the previously published studies employ the geometric mean, shown by Equation (16). This mean was proposed to be employed in the research study of Harrington [17] and later was used by Derringer and Suich [18], among many others. Although other means could be used in this present study such as the harmonic one, in this present study the geometric mean will be used, without weighting the functions.

$$D_{gt}={\left(f_{1t}\ast f_{2t}\ast f_{3t}\right)}^{\frac{1}{3}}$$

(16)

A symmetric 4³ DOE with three output responses, which will be simultaneously optimized has been considered. However, it should be mentioned that this DOE as well as the output responses could be whatever. From the above, in this present study, one of the response variables is preferred to be at its minimum level, another response variable is preferred to be at its maximum and the third output response has to be within a certain range. With the proposed methodology, it is possible to approximate the optimum that simultaneously optimizes the desired levels of variation of these functions within the range of study under consideration. As shown in Section 2, this is usually done through the use of models based on RSM and most of these studies employ the function proposed by Derringer and Suich [18].

In the present study, the levels of variation of the input variables (process parameters) have been varied into the range $\left[-1,1\right]$ as

Table 1. Factors and levels of the symmetrical DOE.

Design Factors	Levels and Values
$x_1$	−1	−1/3	1/3	1
$x_2$	−1	−1/3	1/3	1
$x_3$	−1	−1/3	1/3	1

shows. By using the DOE shown in Table 1, the results shown in

Table 2. Values of the dependent variables using a symmetrical DOE.

Exp.	${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$	Exp.	${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$	Exp.	${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$	Exp.	${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$
1	3.27	3.10	7.72	17	3.18	3.32	5.72	33	3.12	5.21	5.19	49	3.13	19.74	5.05
2	3.09	2.96	6.40	18	3.04	3.19	5.37	34	3.01	5.08	5.10	50	3.02	19.60	5.03
3	3.04	8.34	5.72	19	3.01	8.57	5.19	35	2.99	10.46	5.05	51	3.00	24.98	5.01
4	3.02	19.93	5.37	20	3.01	20.15	5.10	36	3.00	22.04	5.03	52	3.00	36.57	5.01
5	7.09	8.14	21.08	21	6.11	8.27	9.24	37	5.91	9.33	6.12	53	6.50	17.48	5.29
6	4.59	8.01	13.26	22	4.43	8.13	7.18	38	4.23	9.20	5.57	54	4.15	17.35	5.15
7	3.51	13.39	9.24	23	3.45	13.51	6.12	39	3.32	14.58	5.29	55	3.26	22.73	5.08
8	3.17	24.97	7.18	24	3.13	25.10	5.57	40	3.06	26.16	5.15	56	3.04	34.31	5.04
9	10.65	7.31	21.08	25	7.40	7.44	9.24	41	7.71	8.50	6.12	57	11.94	16.65	5.29
10	8.29	7.18	13.26	26	6.39	7.30	7.18	42	6.95	8.36	5.57	58	10.31	16.51	5.15
11	5.27	12.56	9.24	27	4.89	12.68	6.12	43	5.68	13.74	5.29	59	7.51	21.89	5.08
12	3.47	24.14	7.18	28	3.61	24.26	5.57	44	4.37	25.33	5.15	60	5.14	33.48	5.04
13	5.78	8.10	7.72	29	4.11	8.32	5.72	45	4.65	10.21	5.19	61	8.70	24.74	5.05
14	4.23	7.96	6.40	30	3.38	8.19	5.37	46	4.43	10.08	5.10	62	9.82	24.60	5.03
15	2.39	13.34	5.72	31	2.41	13.57	5.19	47	4.39	15.46	5.05	63	12.57	29.98	5.01
16	2.24	24.93	5.37	32	1.69	25.15	5.10	48	5.19	27.04	5.03	64	18.14	41.57	5.01

are obtained for three output responses that have been named $f_1$, $f_2$ y $f_3$, respectively (these values have been generated by using three analytical functions). It should be underlined that either an actual study case or any other function could have been used to carry out this study. As can be observed in Table 2, 4³ experiments will be available. The values shown in Table 2 have been rounded to two decimal places.

Therefore, from the above, the proposed methodology shown in Equations (2)–(15) is to be employed in order to simultaneously minimize the output response $f_1$, maximize the output response $f_2$, and to keep the output response $f_3$ within a given range. In order to analyze a case study this range is considered to be [LB,UB] = [10,15], although any other could be selected. The actual shape of the functions: $f_1$, $f_2$, and $f_3$, used in this study to generate the values of Table 2 is shown in Figure 10, Figure 11, Figure 12, respectively.

From the symmetric DOE shown in Table 2, three regression models and three ANFIS will be developed, and they will be used to determine the values that simultaneously optimize $f_1$, $f_2$, and $f_3$ functions. As will be shown later, regression is not able to adequately model all the response variables. Therefore, in these cases the fact of optimizing surface models which are not accurate enough may lead to erroneous values being obtained. In order to overcome this drawback, soft computing techniques are a powerful tool for modeling output responses. In this present study an ANFIS modeling is employed which can predict the actual behavior of the output responses with higher accuracy than that obtained by using conventional regression.

As previously mentioned, in this present study $f_1$ is preferred to be set at its minimum level, $f_2$ at its maximum, and $f_3$ has to be kept within a range given by a lower and upper bound $\left[LB, UB\right]$, any other combination could have been used, but the one analyzed is considered general enough to show the proposed methodology. Moreover, results obtained with the new proposed desirability function will be compared with those obtained with that proposed by Derringer and Suich [18]. Additional data could be employed to validate the models as suggested in [60]. However, this will not be considered in this present study and all data from the symmetric 4³ factorial DOE will be used to obtain both the regression and the ANFIS models. As previously mentioned, if the regression model can provide a highly-adjusted R-squared value, these models may be used. However, if more accuracy is needed some other techniques such as those based on soft computing should be considered [60].

Summary of the Proposed Desirability Function

In order to summarize the proposed desirability function, this section is included.

The minimum and maximum ranges of variation for each variable will be first determined. In the event where the preferred variation of a variable is the maximum, the procedure shown in Equation (17) will be followed; in the case where the preferred variation is the minimum, that indicated in Equation (18) and if the variable has to be within a specific interval [LB,UB] then the procedure shown in Equation (19) should be followed. Once the variables have been transformed, then a mean is to be maximized, for example, similar to that shown by Equation (16).

$$\begin{gathered} If max (y_j)\stackrel{transform to}{\to }{y}_{tj} \\ =\frac{1}{2}+\frac{1}{\pi }\ast ta{n}^{-1}\left(a_j\ast y_j+b_j\right) where \left\{\begin{array}{c}{a}_j=\frac{1}{M_j-c_{M_j}}\left(H_{M_j}-tan\left(\pi (t_{c_{M_j}}-\frac{1}{2} \right)\right)\\ b_j=\frac{-1}{M_j-c_{M_j}}\left(H_{M_j}\ast c_{M_j}-M_j\ast tan\left(\pi (t_{c_{M_j}}-\frac{1}{2} \right)\right)\\ c_{M_j}=\frac{m_j+k_{c_{M_j}}\ast M_j}{1+k_{c_{M_j}}}\\ ---\\ If a value of t_{c_{M_j}}is to be selected, then \\ Select k_{H_{M_j}} so that\\ H_{M_j}=k_{H_{M_j}}\ast \left(M_j-m_j\right)\\ where\\ k_{H_{M_j}}> \frac{1}{M_j -m_j}tan\left(\pi \left(t_{c_{M_j}} -\frac{1}{2}\right)\right)\\ ---\\ or if a value of k_{H_{M_j}}is to be selected, then\\ Select t_{c_{M_j}} so that\\ 0\le t_{c_{m_j}}<\frac{1}{2}+\frac{1}{\pi }\ast ta{n}^{-1}\left(k_{H_{m_j}}\ast \left(M_j -m_j\right)\right)\\ H_{M_j}=k_{H_{M_j}}\ast \left(M_j-m_j\right)\end{array}\right. \end{gathered}$$

(17)

$$\begin{gathered} If min (y_j)\stackrel{transform to}{\to }{y}_{tj} \\ =\frac{1}{2}+\frac{1}{\pi }\ast ta{n}^{-1}\left(a_j\ast y_j+b_j\right) where \left\{\begin{array}{c}{a}_j=\frac{1}{m_j-c_{m_j}}\left(H_{m_j}-tan\left(\pi (t_{c_{m_j}}-\frac{1}{2} \right)\right)\\ b_j=\frac{-1}{m_j-c_{m_j}}\left(H_{m_j}\ast c_{m_j}-m_j\ast tan\left(\pi (t_{c_{m_j}}-\frac{1}{2} \right)\right)\\ c_{m_j}=\frac{k_{c_{m_j}}\ast m_j+M_j}{k_{c_{m_j}}+1}\\ ---\\ If a value of t_{c_{m_j}}is to be selected, then \\ Select k_{H_{m_j}} so that\\ H_{m_j}=k_{H_{m_j}}\ast \left(M_j-m_j\right)\\ where\\ k_{H_{m_j}}> \frac{1}{M_j -m_j}tan\left(\pi \left(t_{c_{m_j}} -\frac{1}{2}\right)\right)\\ ---\\ or if a value of k_{H_{m_j}}is to be selected, then\\ Select t_{c_{m_j}} so that\\ 0\le t_{c_{m_j}}<\frac{1}{2}+\frac{1}{\pi }\ast ta{n}^{-1}\left(k_{H_{m_j}}\ast \left(M_j -m_j\right)\right)\\ H_{m_j}=k_{H_{m_j}}\ast \left(M_j-m_j\right)\end{array}\right. \end{gathered}$$

(18)

$$\left\{\begin{array}{c}If LB\le y_j\le UB \stackrel{transform to}{\to }{y}_{tj}=\left|1-\left(y_{t1}+y_{t2}\right)\right| \\ \\ where y_{t1} is a variable that has been transformed as indicated in \mathrm{Equation} \left(17\right)\\ \left[m_1,M_1\right]=\left[\frac{LB+UB}{2},UB\right]\\ \\ and y_{t2} is a variable that has been transformed as indicated in \mathrm{Equation} \left(18\right)\\ \\ \left[m_2,M_2\right]=\left[LB,\frac{LB+UB}{2}\right]\end{array}\right.$$

(19)

Although $m_1$ and $M_2$ values that are shown in Equation (19) could be selected in a different way, in this study these values have been selected as they are shown.

Finally, once some or all of the transformation shown in Equation (17), Equation (18), and Equation (19) have been carried out over the response variables, previously obtained with the ANFIS, then an average with these transformed variables (geometric, harmonic, etc) should be employed and it should be maximized. In this present study, the non-weighted geometric mean is to be used.

In the case where, when maximizing or minimizing a variable, it is not possible to obtain response values of this variable higher or lower than a given one, it is possible to use the arctangent transformations which lead to shapes as shown in Figure 5, Figure 6, Figure 7, Figure 8. The procedure to follow in these cases, consists in filtering said output responses, as follows. Let us suppose that the desired level of a variable is at its maximum. In this case, first a transformed variable $\left(f_{tj}{}_1\right)$ is used to maximize it, according to Equation (17), and second a transformed variable $\left(f_{tj}{}_2\right)$ is used to limit the output response, according to Equation (18); then, the resulting data will be filtered above the maximum (i.e, $f_{tj}=f_{tj}{}_1\ast f_{tj}{}_2$). Similarly, if the output response has to be minimized, first a transformed variable $\left(f_{tk}{}_1\right)$ is used, according to Equation (18) and the output response is limited by using a transformed varible $\left(f_{tk}{}_2\right)$, according to Equation (17), then the resulting data will be filtered below the minimum (i.e., $f_{tk}=f_{tk}{}_1\ast f_{tk}{}_2$). On the other hand, in the case where a variable $\left(f_j\right)$ has to be maximized or minimized and at the same time has to be within a certain range. Then, first a transformed variable $\left(f_{ti}{}_1\right)$, according to Equation (17) or Equation (18) is used, and second a transformed variable $\left(f_{ti}{}_2\right)$ is used to keep the output response in the pre-set interval, according to Equation (19).

The geometric mean is calculated according to Equation (20) $\left(l=1..m\right)$, where $m$ is the number of said variables and $n$ is the total number of variables).

$$D_{gt}={\left(\left({\displaystyle \prod }_{l=1}^m{f}_t{}_l\right)f_t{}_{m+1}\ast f_t{}_{m+2}\ast \dots \ast f_t{}_n\right)}^{\frac{1}{n}}$$

(20)

This process can be done if there are $m$ variables that have to be maximized/minimized and other restrictions also apply to them.

In order to simplify this present study, unconstrained restrictions will be considered for the variables. That is to say, the output responses are allowed to take values higher or lower, depending on the case, than those provided by the desired min and desired max levels, which have been taken from the DOE, as Table 2 shows. However, they could be restricted following the above-mentioned procedure.

4. Results and Discussion

This section presents the results obtained by modeling the response variables, first, by means of using conventional regression. Given the limitations that such modeling presents, an ANFIS will be used later to overcome these stated drawbacks. As shown Section 2, it is very common to employ Derringer and Suich [18] desirability function for multiple optimization of output responses. In this present study results obtained with this desirability function and those obtained with the proposed one, which is based on the arctangent function, will be compared in order to simultaneously optimize three output responses.

The use of models based on response surfaces, together with a desirability function, is commonly used to model response variables in manufacturing engineering and to solve a multiple optimization problem because it is possible to find a solution using optimization methods. The proposed desirability function, based on a transformation using the arctangent function, as previously shown, could also be employed with RSM and optimization methods. In this case, t has to be higher than zero $(t>0)$. If necessary, a value of t close to zero could be employed, for example, $t=\epsilon \left({10}^{-3}\dots {10}^{-5}\right)$ would be enough although, if necessary, it could be reduced. However, in the event that the regression does not provide high values of the coefficient of determination, it is very likely that these regression models will not provide accurate results and therefore it would be better to employ other models for the output responses for example based on soft computing, because, in this way, more accurate results can be obtained.

4.1. Regresion Model

Equation (21) shows the regression model that will be first employed to predict the behavior of the output responses obtained from the DOE shown in Table 2. The full regression models will be used instead of those obtained from the higher adj-R², since the fact of using the model adjusted to the degrees of freedom could eliminate some of the independent variables.

$$\begin{gathered} y ~ (a_0+a_1\ast x_1+a_2\ast x_2+a_3\ast x_3+a_4\ast x_1\ast x_2+a_5\ast x_1\ast x_3+ \\ {a}_6\ast x_2\ast x_3+a_7\ast x_1^2+a_8\ast x_2^2+a_9\ast x_3^2) \end{gathered}$$

(21)

Equation (22) shows that the validation metrics employed in this study are: mean squared error ($MSE$), root mean squared error ($RMSE$), and mean absolute error ($MAE$), where $y_j$ are the values provided by the DOE and ${\widehat{y}}_j$ are the predicted values. As previously mentioned, the model which exhibit the higher coefficient of determination R² is to be used.

$$MSE=\frac{1}{n}{\displaystyle \sum }_{j=1}^n{\left(y_j-{\widehat{y}}_j\right)}^2, RMSE=\sqrt{MSE} and MAE=\frac{1}{n}{\displaystyle \sum }_{j=1}^n\left|y_j-{\widehat{y}}_j\right|$$

(22)

Table 3. Regression results for the first output response.

${\mathit{f}}_1$	${\mathit{R}}^2$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{A}\mathit{E}$
Regression	0.6273	1.7893	3.2015	1.2549
$$\begin{gathered} f_1=4.4567+1.2109\ast x_1+1.6519\ast x_2-0.98981\ast x_3+2.0487\ast x_1\ast x_2+0.85815 \\ \ast x_1\ast x_3+0.07785\ast x_2\ast x_3+1.7958\ast x_1^2-1.2473\ast x_2^2+0.47391 \\ \ast x_3^2 \end{gathered}$$

,

Table 4. Regression results for the second output response.

${\mathit{f}}_2$	${\mathit{R}}^2$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{A}\mathit{E}$
Regression	0.9523	1.9817	3.9273	1.6819
$$\begin{gathered} f_2=8.7298+6.0665\ast x_1+2.125\ast x_2+8.3802\ast x_3-0.00039405\ast x_1\ast x_2 \\ +0.00017002\ast x_1\ast x_3-0.00016849\ast x_2\ast x_3+6.2793\ast x_1^2 \\ -0.053086\ast x_2^2+6.5921\ast x_3^2 \end{gathered}$$

,

Table 5. Regression results for the third output response.

${\mathit{f}}_3$	${\mathit{R}}^2$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{A}\mathit{E}$
Regression	0.6601	1.8221	3.3201	1.2238
$$\begin{gathered} f_3=6.5491-2.1172\ast x_1+5.8144\mathrm{e}-06\ast x_2-1.3545\ast x_3+2.75\mathrm{e}-06\ast x_1\ast x_2 \\ +1.8884\ast x_1\ast x_3+7.7464\mathrm{e}-07\ast x_2\ast x_3+1.7325\ast x_1^2-2.4258\ast x_2^2 \\ +0.63984\ast x_3^2 \end{gathered}$$

show the obtained results for $f_1$, $f_2$, and $f_3$ functions, respectively. As can be observed, regression results are not good enough for $f_1$ and $f_3$ functions and hence these models will not be accurate to predict the behavior of these response variables. On the other hand, the model obtained for $f_2$ function could be used to approximately predict the behavior of this output response which is a relatively simple equation. However, if more accuracy is required, then other models should be used, for example those based on soft computing.

Figure 13 shows the response surface obtained with the regression model of $f_1$. Because the coefficient of determination of the regression is not very high, the results will not be very reliable with this modeling. Therefore, the fact of considering an optimization using a desirability function with this kind of functions may lead to values with a lack of precision obtained. On the other hand, if it is compared with the actual shape of the real surface of $f_1$ function, which has been shown in Figure 10, the aforementioned can be observed.

Table 4 shows the results obtained for the case of the second function considered in this study. Figure 14 shows the response surface predicted by the regression model. In this case, the function could be used to address the proposed problem since, although the coefficient is not very high, it could give results close to the real optimum of the variable. It may be seen that the shape of the response surface given by the regression model is more similar to the actual shape of the function in this case than in the previous one.

Table 5 shows the results obtained for the case of the third function considered in this study and Figure 15 shows the response surface predicted by the regression model. In this case, as with the first function, the regression model is not able to adequately predict the behavior of the function, which can be observed when comparing the shape of the response surface obtained with the regression model with the actual form of the function, which is shown in Figure 12.

Based on the above, although there are three polynomial models available with which it is possible to use a desirability function and then obtain an optimum, this value will not be precise and it may also occur that this optimum falls outside the limits established in the optimization for some of the variables, as a consequence of the models thus developed; therefore, not having high values of the regression determination coefficients, they will not adequately predict the actual behavior of the response variables.

4.2. Adaptive Neuro Fuzzy Inference System (ANFIS) Modeling

In order to avoid the above-mentioned drawbacks of the regression models, soft computing modeling based for example on artificial neuronal networks (ANNs) could be used. In [61], different types of ANNs are shown. However, in this present study a FIS which is later tuned by using an ANFIS is considered [3]. Specifically, a zero-order Sugeno FIS is employed by using the Fuzzy Logic Toolbox™ of MATLAB^TM2020a [62], because the de-fuzzification process for a Sugeno system is computationally more efficient compared to that of a Mamdani system [62,63]. However, it should be stated that a Mamdani FIS could have been used. Furthermore, the membership functions for fuzzification of the independent variables are Gaussian as shown in Equation (23). These membership functions have been developed from the methodology shown in [60], although some other membership functions could have been used.

As previously mentioned, in the present study, no additional points will be used to optimize the FIS models as suggested in [60]. It should be noted that if a greater number of data is available to this end, then the ANFIS models obtained will be more reliable. In any case, it can be observed in this section that the models that the ANFIS provide much better results than those obtained with the regression.

$${\mu }_x=e^{\frac{-{\left(x-c\right)}^2}{2{\sigma }^2}}$$

(23)

The aggregation method is the sum of the fuzzy sets and the aggregated output is obtained from the weighted average of all output rules, where product implication method is used in Sugeno systems [62]. For the ${\mathrm{i}}_{\mathrm{th}}$ rule, the implication method is obtained from Equations (24) and (25) shows the output of the FIS.

$${\lambda }_j\left(x\right)=AndMethod\left\{{\mu }_{j1}\left(x_1\right),\dots ,{\mu }_{jn}\left(x_n\right)\right\}$$

(24)

$$\left\{f_1,f_2,f_3\right\}=\frac{{{\displaystyle \sum }}_{j=1}^{Number of rules}{\lambda }_j\ast z_j}{{{\displaystyle \sum }}_{j=1}^{Number of rules}{\lambda }_j}$$

(25)

The created FIS will have a set of $“l”$ rules of the form

$$If \left(x_1 is x_{1,i}\right) and \left(x_2 is x_{2,j}\right) and \left(x_3 is x_{3,k}\right) then \left(\left[f\right] is outpu{t}_l\right)$$

Figure 16 shows the structure of the ANFIS employed in this present study and Figure 17 shows the response surfaces obtained with the ANFIS developed for the first response variable and the values of the statistics employed are shown in

Table 6. ANFIS results for the first output response.

${\mathit{f}}_1$	${\mathit{R}}^2$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{A}\mathit{E}$
ANFIS	1.0000	0.0050	2.4881 × 10−5	0.0039

. It can be noted that they are much better than those obtained in the case of the regression for this variable, which were shown in Table 3. Likewise, it can be observed that the ANFIS is capable of adequately predicting the real form of the function despite not having optimized the FIS with additional points.

A similar result is obtained for ${\mathrm{f}}_2$ function as can be observed both in

Table 7. ANFIS results for the second output response.

${\mathit{f}}_2$	${\mathit{R}}^2$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{A}\mathit{E}$
ANFIS	1.0000	7.9133 × 10−4	6.2621 × 10−7	6.5990 × 10−4

and in Figure 18. Analogously to the previous case, results obtained with the ANFIS improve very much on those obtained with the regression.

Likewise,

Table 8. ANFIS results for the third output response.

${\mathit{f}}_3$	${\mathit{R}}^2$	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{S}\mathit{E}$	$\mathit{M}\mathit{A}\mathit{E}$
ANFIS	1.0000	0.0038	1.4506 × 10−5	0.0033

and Figure 19 show the obtained results with the ANFIS for the case of $f_3$ function. It is observed that this function is now adequately modeled, contrary to that observed with the regression model.

As can be seen in the results obtained in this Section, the precision obtained with the ANFIS models is much higher than that obtained with the regression. Therefore, the prediction developed using ANFIS will be more accurate, especially for regression models not having highly R-squared coefficients. It is worth mentioning the fact that the availability of analytical models, such as those provided by the regression, makes it possible to employ them with optimization algorithms when they are combined with a desirability function. In that way, the desirability function based on the arctangent, which has been shown in Section 3, could be employed to obtain a value that simultaneously optimizes the required conditions by using a desirability function. As previously mentioned, in this case t has to be higher than zero $(t>0)$. A value of $t=\epsilon {10}^{-3}\dots {10}^{-5}$ may be enough, if necessary, it could be reduced. However, if regression models are not accurate it is preferrable either to employ, as an approach, the values obtained directly from the DOE or to employ soft computing techniques. Because, if the regression models are not accurate, it is very likely that the optimization will not be so either. In order to overcome this drawback, the ANFIS models will be employed in order to approach the optimum of the functions by using a desirability function.

4.3. Analysis of the Case Study

In order to compare the results attained, the values obtained with the proposed methodology will be compared with those obtained by using the transformation suggested by Derringer and Suich [18]. The symmetrical DOE shown in Table 2 with three output variables $f_1$, $f_2$, and $f_3$ will be used. As previously mentioned, the preferred level of $f_1$ is its minimum, the preferred level of $f_2$ is its maximum and $f_3$ has to be within a set range, defined by a lower bound and by an upper bound. For example, $\left[LB,UB\right]=\left[10,15\right]$, although this range could be whatever. In order to set the maximum and minimum level of each variable, those provided by the DOE are considered, as shown in

Table 9. Ranges of the output responses.

Function	Min	Max
$f_1$	1.69	18.14
$f_2$	2.96	41.57
$f_3$	5.01	21.08

, although they could be any other value.

Figure 20b₁–b₃ show how the data of the output variables obtained from the DOE are classified by using the Derringer and Suich [18] transformation, where it has been considered (r = 3; for one side transformation, s = t = 3, for two-side transformation and c = (LB + UB)/2). Figure 20a₁–a₃ show these data classified by using the proposed arctangent transformation, in the case of minimizing $f_1$, maximizing $f_2$, and in order to keep $f_3$ within a certain range.

Table 10. Values of the constants to transform the $f_1$ and $f_2$ outputs.

Function	Desired Level	${\mathit{k}}_{{\mathit{c}}_{{\mathit{M}}_{\mathbf{1}}}}$	${\mathit{t}}_{{\mathit{c}}_{{\mathit{M}}_{\mathbf{1}}}}$	${\mathit{k}}_{{\mathit{H}}_{{\mathit{M}}_{\mathbf{1}}}}$	${\mathit{m}}_{\mathbf{1}}$	${\mathit{M}}_{\mathbf{1}}$
${\mathrm{f}}_1$	Max.	2	0.2	1	1.69	18.14
		${\mathit{k}}_{{\mathit{c}}_{{\mathit{m}}_{\mathbf{2}}}}$	${\mathit{t}}_{{\mathit{c}}_{{\mathit{m}}_{\mathbf{2}}}}$	${\mathit{k}}_{{\mathit{H}}_{{\mathit{m}}_{\mathbf{2}}}}$	${\mathit{m}}_{\mathbf{2}}$	${\mathit{M}}_{\mathbf{2}}$
${\mathrm{f}}_2$	Min.	1	0.5	1	2.96	41.57

shows the values of the constants to transform the $f_1$ and $f_2$ outputs using the arctangent transformation. Where it was desired to weight the values of $f_1$. Hence the values shown in this table. On the other hand, to keep $f_3$ within the desired range [LB,UB] = [10,15], where $k_H=5, k_c=1, t_c={10}^{-5}$. However, other values could have been used. As can be observed in Table 2 the range of values specified for $f_3$ to stay within only has one output value of this variable although it is obtained in two combinations of the independent variables.

Table 11. Desirability values obtained from the DOE.

Exp.	${\mathit{D}}_{\mathit{g}\mathit{t}}$	${\mathit{x}}_{\mathbf{1}}$	${\mathit{x}}_{\mathbf{2}}$	${\mathit{x}}_{\mathbf{3}}$	${\mathit{f}}_{\mathbf{1}}$	${\mathit{f}}_{\mathbf{2}}$	${\mathit{f}}_{\mathbf{3}}$
6	0.2201	−1	−1/3	−1/3	4.59	8.01	13.26
Exp.	${\mathit{D}}_{\mathit{g}\mathit{d}}$	${\mathit{x}}_{\mathbf{1}}$	${\mathit{x}}_{\mathbf{2}}$	${\mathit{x}}_{\mathbf{3}}$	${\mathit{f}}_{\mathbf{1}}$	${\mathit{f}}_{\mathbf{2}}$	${\mathit{f}}_{\mathbf{3}}$
6	0.0750	−1	−1/3	−1/3	4.59	8.01	13.26

shows the optimum values found directly with the data provided by the symmetrical DOE, which is shown in Table 2, where $D_{gt}$ refers to the desirability value obtained with the arctangent transformation and $D_{gd}$ that obtained by using Derringer and Suich [18], both results are obtained by using the non-weighted geometric average. However, it may occur that the DOE does not provide values within a certain range for example, in this case study, within the range [15,20] or that we want to explore the possibility of finding improved values, then if regression results are not accurate, as is the present case, an ANFIS is to be used to overcome this drawback.

It could be possible to define a relative index to compare both desirability functions $(D_{gt}$, $D_{gd}$) however this is not carried out in this study. As can be observed in Table 11 the optimal values that both functions find, using only the data provided by the DOE, are the same which is logical as the number of points is not high enough to allow significant differences to exist. In order to employ the ANFIS models, whenever the regression is not accurate enough, first the initial DOE which was employed to develop these ANFIS models, is discretized by using N points. In this study the same number of points has been employed for each of the three independent variables, in order to have an “augmented-symmetrical N³ design of experiments”, although the number of points to discretize each independent variable could be different. From the above, the DOE is discretized as Equation (26) shows. That is, N points are obtained within the range $\left[\min \left(x_j\right),\max \left(x_j\right)\right]$.

$$x_j=\min \left(x_j\right) : \frac{\max \left(x_j\right)-\min \left(x_j\right)}{N-1} : \max \left(x_j\right)$$

(26)

First of all, all combinations of the data provided by Equation (26) are evaluated making use of the previously developed ANFIS models, and hence N³ outputs will be obtained for the response variables, which can be then combined to determine the optimum value of a previously defined parameter, which is, in this present study the geometric mean. This optimum value will be linked to a point, in this case, ($x_1$,$x_2$,$x_3$) withing the selected range and this will be the value that approximate to the optimum of the function. The higher the number of the selected levels, the higher the accuracy of the approach. If N is high enough the approach to the optimum will be accurate. A convergence criterion could be established based, for example, on the convergence band obtained from the optimal values found in each iteration with the desirability functions depending on the number of points (N), among others. However, this is not considered in this study and the optimal values are determined after considering a value of N high enough to observe the convergence of the optimal values that provide the desirability functions in each iteration. Figure 21a,b show the optimal values of the desirability function versus the number of points (N) employed in the ANFIS (in this case $N^3$ combinations of the variables, because there are three independent variables) where $D_{gt}$ is the desirability function obtained by using the geometric mean and the proposed transformation and $D_{gd}$ the desirability values obtained from the geometric mean and the methodology of Derringer and Suich [18]. Moreover, a comparison with the value provided by using the initial data of the DOE can be observed in this figure for both desirability functions when using the ANFIS to generate the values of the output responses. As can be observed, results attained improve those obtained directly with the initial DOE.

Figure 22 show the transformed function by using the ANFIS models both with the proposed desirability function and with that of Derringer and Suich [18].

Table 12. Results attained with the ANFIS models and the desirability functions.

Arctangent	${\mathit{D}}_{\mathit{g}\mathit{t}}$	${\mathit{x}}_{\mathbf{1}}$	${\mathit{x}}_{\mathbf{2}}$	${\mathit{x}}_{\mathbf{3}}$	${\mathit{f}}_{\mathbf{1}}$	${\mathit{f}}_{\mathbf{2}}$	${\mathit{f}}_{\mathbf{3}}$
Num. Points: 55	0.2535	−1.0000	−0.0526	0.2982	4.5892	12.9599	10.0063
Actual values (i.e., obtained with the actual functions)					4.2857	11.6280	10.4535
Derringer & Suich	${\mathit{D}}_{\mathit{g}\mathit{d}}$	${\mathit{x}}_{\mathbf{1}}$	${\mathit{x}}_{\mathbf{2}}$	${\mathit{x}}_{\mathbf{3}}$	${\mathit{f}}_{\mathbf{1}}$	${\mathit{f}}_{\mathbf{2}}$	${\mathit{f}}_{\mathbf{3}}$
Num. Points: 97	0.1202	−1.0000	−0.2121	−0.1717	4.9343	8.7416	12.5005
Actual values (i.e., obtained with the actual functions)					4.8071	8.1795	13.0183

presents the optimal values found with the ANFIS. As can be observed in Table 12, the optimum value with the proposed desirability function is obtained after using N = 55 points, to be evaluated with the ANFIS models, that is (55³ combinations of the variables) and it leads to the values of ${\left[f_1,f_2,f_3\right]}_{D_{gt}}=\left[4.5892, 12.9599, 10.0063\right]$ which, in the case of this present study, are improved values, of the response variables, compared with those obtained with Derringer and Suich [18] desirability function (after N = 97 iterations, 97³ combinations of the variables evaluated with the ANFIS models),which leads to ${\left[f_1,f_2,f_3\right]}_{D_{gd}}=\left[4.9343, 8.7416, 12.5005\right].$ That is, with the proposed desirability function $f_1$ is more minimized and $f_2$ is more maximized while $f_3$ is kept within the desired range, which may be a consequence of the weighting values employed in the Derringer and Suich [18] transformation, which were fixed at (r = s = t = 3), as previously mentioned. It may occur that other values of this transformation could lead to improved results. Moreover, it is also likely that the attained results with the arctangent may also vary if different constants of those shown in Table 10 were selected. However, this is not analyzed in this present study.

Results obtained with the regression models are shown in

Table 13. Results attained with the regression models and the desirability functions.

Arctangent	${\mathit{x}}_{\mathbf{1}}$	${\mathit{x}}_{\mathbf{2}}$	${\mathit{x}}_{\mathbf{3}}$	${\mathit{f}}_{\mathbf{1}}$	${\mathit{f}}_{\mathbf{2}}$	${\mathit{f}}_{\mathbf{3}}$
Regression Results	−1.0000	0.1489	0.1064	4.7649	10.2241	10.0073
Actual values (i.e., obtained with the actual functions)				5.7264	9.2598	11.3551
Derringer & Suich	${\mathit{x}}_{\mathbf{1}}$	${\mathit{x}}_{\mathbf{2}}$	${\mathit{x}}_{\mathbf{3}}$	${\mathit{f}}_{\mathbf{1}}$	${\mathit{f}}_{\mathbf{2}}$	${\mathit{f}}_{\mathbf{3}}$
Regression Results	−1.0000	0.8571	−1.0000	6.0402	8.9376	12.4994
Actual values (i.e., obtained with the actual functions)				7.0248	8.1123	9.6217

. As can be observed, the results provided by the regression do not fit very well with the actual values of the function and in the case of using the Derringer and Suich [18] desirability function are outside specified tolerances for $f_3$ since the tolerance range specified for this response variable was $10\le f_3\le 15$.

Therefore, as shown in Table 12 and Table 13 results obtained with the ANFIS are more accurate than those provided by using the regression models. Moreover, in the case of using regression models with lack of precision non-desirable results could be obtained as Table 13 shows.

Table 14. Comparison of the results obtained.

		${\mathit{f}}_1$	${\mathit{f}}_2$	${\mathit{f}}_3$
Derringer & Suich	ANFIS	4.9343	8.7416	12.5005
	Regression	6.0402	8.9376	12.4994
Arctangent	ANFIS	4.5892	12.9599	10.0063
	Regression	4.7649	10.2241	10.0073

presents a summary of the results obtained with the Derringer&Suich desirability function [18] and the new proposed desirability function and

Table 15. Actual values of $f_1, f_2$ and $f_3$ functions.

$f_1=3+\frac{{\left(x_1^2-x_2^2+x_1\ast x_2\ast \exp \left(x_3\right)+2\right)}^2-\sin {\left(x_2\right)}^2-sin\left(x3\right)}{1+exp\left(x_1^2-4\ast x_2+2\ast x_3\right)}$

$$\begin{gathered} f_2=10+{\left(0.65+x_1\right)}^3\ast exp\left(0.65\ast (x_1^2+x_2^2\right))-{\left(x_2+1.25\ast \cos \left(\pi \ast x_2\right)\right)}^2+2 \\ \ast {\left(2+x_3\right)}^2\ast \sin \left(x_3\right) \end{gathered}$$

$f_3=5+\frac{1}{\exp \left(2\ast x_1+2\ast x_2^2+x_3\right)}$

shows the actual values of the functions used to generate the data shown in Table 2.

It should be mentioned that in this case study the desirability functions do not provide high values due to the fact that $f_3$ has to be kept within a certain range. Some other combinations could lead to different values being found. For example, Figure 1 shows the values of the proposed arctangent desirability function for a multiple optimization problem where the preferred level of variation of $f_3$ was selected to be the minimum value instead of being fixed within a tolerance range and the maximum was the preferred level for both $f_1$ and $f_2$. However, it may be that a technological interest exists in maintaining an output response within a certain range while other variables are being simultaneously optimized, and as has been shown in this present study, a design of experiments combined with an ANFIS and a desirability function may obtain an accurate approach to the optimum value even in the event that the DOE does not provide values in a certain range which has to be selected for a variable to be kept within. On the other hand, as can be seen for example in the research study of Derringer [64] it is often observed that when one property improves, it is usually at the expense of others. By using the proposed methodology, it is possible to have all available data attained by the ANFIS, that is, the combination of both independent and dependent variables which lead to a specific desirability value and then allow the user to select the most appropriate combination of the independent variables (process parameters) which lead to reaching a certain value in the response variables. Therefore, it could be possible to establish classes of desirability and several solutions can be obtained instead of a single optimum since it may happen that several combinations of the independent variables may lead to very close values and the fact of selecting the better result could avoid some others being taken into account, which may, in turn, be of technological interest.

Figure 23a–d and Figure 24a–d show the desirability functions versus the number of variables in both cases examined in this present study as well as the values of the response variables which correspond to a selected level of the desirability function (these figures show the combination of the independent variables that lead to desirability function values being obtained that are higher than 80% of the maximum desirability value of each desirability function respectively). Therefore, it is possible to see the desired levels of the process parameters as well as their influence on the desirability function.

5. Conclusions

This present research study shows a methodology based on using a design of experiments along with an ANFIS and a desirability function for simultaneous multiple response optimization which can be applied to overcome the drawbacks of the regression techniques.

A new desirability function has been proposed and analyzed in this research study. The proposed desirability function is based on the arctangent function and it may adopt several shapes as a function of its parameters.

By using the proposed methodology, it could be possible to establish classes of desirability along with the values of the process parameters and the responses variables and then select the most appropriate ones with higher accuracy than that obtained by conventional regression techniques.