Metodologías de Superficie de Multirrespuesta

ROMÁN DE LA VARA SALAZAR

JORGE DOMÍNGUEZ DOMÍNGUEZ CIMAT, A. C. Apdo. Postal 402 Guanajuato, Gto., 36000

Abril, 1998

Resumen

Es común que la calidad de un producto dependa de los valores que toman más de dos de sus propiedades o características. De aquí la relevancia de conocer técnicas que sirvan para optimizar varias respuestas simultáneamente. En este artículo se presentan algunos métodos que se han propuesto hasta ahora para resolver dicho problema. Se estudia el desempeño de cuatro de los métodos al aplicarlos a los mismos datos experimentales. Se discuten las ventajas y desventajas de cada método y se proponen mejoras en dos de ellos en la manera de ponderar las variables.

1. Introducción

En la literatura de diseño de experimentos y de optimización de procesos se hace mucho énfasis en ejemplos y problemas que consideran el estudio de una sola propiedad o característica de calidad del producto o proceso, cuando lo típico es que la calidad del producto dependa del valor que asumen más de una de sus propiedades (ver por ejemplo Myers & Montgomery (1995) y Box et al. (1988)). Por ejemplo, un alimento tiene varias propiedades como: digestibilidad, textura, pH, sabor, aspecto, etc., y todas tienen su importancia para que el alimento sea bien aceptado por los consumidores. Muchas veces el optimizar solo una de las características del producto hace que las otras propiedades se vean afectadas y que el resultado sea un producto con peor calidad global que antes. De aquí la importancia de contar con técnicas que sirvan para optimizar de manera simultánea a todas las respuestas de interés, o al menos localizar al punto de operación del proceso en el que todas las variables tienen el mejor desempeño posible.

Sin embargo, el problema de optimización simultánea de varias respuestas puede no tener una solución completamente satisfactoria desde la perspectiva de cada respuesta individual, en el sentido de que generalmente los óptimos individuales no coinciden con el óptimo global; aunque siempre podemos decir que existe el punto de operación del proceso en el que todas las variables cumplen con los requerimientos “de la mejor manera posible”. Así pues, típicamente se obtiene una solución de compromiso, que busca balacear de manera adecuada la importancia relativa de las respuestas.

En este trabajo se describen algunos de los métodos de optimización simultánea de más de dos respuestas que se han propuesto hasta ahora por diferentes autores (ver Derringer y Suich (1980, 1994), Khuri y Conlon (1981), Ames, et al (1997), Pignatiello (1993) y Del Castillo (1996)). Se hace un estudio comparativo de cuatro de estos métodos vía su desempeño en un caso estudio y se discuten las ventajas y desventajas en cada uno de ellos. Se proponen mejoras a algunos las técnicas en lo concerniente a lograr un mejor balance de la importancia relativa de cada respuesta. Entre los métodos a comparar incluimos el método gráfico el cual proporciona una idea intuitiva del problema de multirrespuesta y puede dar buenos resultados si el número de factores controlados no es mayor que cuatro; consideramos que el método gráfico no ha sido valorado adecuadamente. Proponemos el método gráfico como un primer paso exploratorio del problema de optimización simultánea, para minimizar la posibilidad de una solución incorrecta o una solución que es óptimo local, y por tanto, no es el óptimo global que se busca. El análisis de un solo caso estudio no puede ser concluyente, pero basta un ejemplo para ilustrar la mayoría de las ventajas y desventajas de los métodos.

Todos los métodos de multirrespuesta reducen de alguna manera el problema multivariado al problema de optimizar una sola variable. Asimismo, todos los métodos tienen algún mecanismo que permite ponderar las variables de interés de acuerdo a su importancia relativa o de acuerdo a su variabilidad, aunque no es claro como elegir dichos pesos. Todos los métodos suponen que las k respuestas (Y_{1}, Y_{2}, ..., Y_{k}) a optimizar se pueden modelar adecuadamente por un modelo polinomial de segundo orden en términos de p factores de control (X_{1}, X_{2}, ..., X_{p}), es decir, se tienen los modelos estimados dados por

Ŷ_{1} = β̂_{01} + ∑_{i=1}^{p} β̂_{i1}X_{i} + ∑_{i=1}^{p} β̂_{ii1}X_{i}^{2} + ∑ ∑_{i<j}^{p} β̂_{ij1}X_{i}X_{j}  (1.1)

Ŷ_{2} = β̂_{02} + ∑_{i=1}^{p} β̂_{i2}X_{i} + ∑_{i=1}^{p} β̂_{ii2}X_{i}^{2} + ∑ ∑_{i<j}^{p} β̂_{ij2}X_{i}X_{j}
⋮
Ŷ_{k} = β̂_{0k} + ∑_{i=1}^{p} β̂_{ik}X_{i} + ∑_{i=1}^{p} β̂_{iik}X_{i}^{2} + ∑ ∑_{i<j}^{p} β̂_{ijk}X_{i}X_{j}.

Se supone que se cumplen los supuestos tradicionales de normalidad, varianza constante e independencia de los residuos, y que el coeficiente de determinación R_{aj}^{2} de cada modelo es al menos 0.70. Algunas de las técnicas (ver por ejemplo, Khuri y Conlon (1981)) suponen que las k variables de respuesta son linealmente independientes y en ese caso se deberán eliminar tales dependencias; algunas consideran la aleatoriedad de los valores predichos y otras no. Algunos de los métodos requieren el conocimiento de límites de especificación para cada una de las variables, y todos ellos necesitan saber que es lo que se busca con cada variable: si se busca un máximo, un mínimo o un cierto valor objetivo. Una vez definida la función objetivo global de acuerdo al método deseado, se utiliza un algoritmo de optimización para determinar el punto que optimiza dicha función. En este trabajo utilizamos el algoritmo de Nelder-Mead (Nelder y Mead (1965)), que no requiere el cálculo de derivadas, y el de Broyden Fletcher, Goldfarb, y Shanno (ver, Luenberger (1989)) que si las requiere, programados ambos en Gauss, V3.14.

Cabe agregar que el caso particular de optimizar dos respuestas, llamado el “problema de respuesta dual”, ha recibido bastante atención en la literatura reciente y en principio las técnicas de optimización que aquí discutimos son directamente aplicables a este caso, pero también se han propuesto varias técnicas específicamente para abordar este problema (ver por ejemplo, Copeland y Nelson (1996)).

2. Método gráfico

Quizás lo primero que se ocurre al optimizar varias respuestas es superponer todas las superficies sobre la región experimental y localizar adentro de ella subregiones en las cuales todos los modelos predicen valores aceptables para las respuestas.

Cuando se tienen solo dos factores de control (X_{1} y X_{2}) es bastante fácil superponer las superficies de respuesta en su forma de contornos sobre la región experimental, que en este caso se puede dibujar como un cuadrado o círculo centrado en el origen. El problema es que cuando se tienen más de dos factores las superficies no se pueden dibujar de una sola vez sobre toda la región experimental; aunque se pueden estudiar y superponer tomando dos factores a la vez y fijando a los restantes. Por ejemplo, si son tres los factores de control, la región experimental es un cubo (o esfera) centrado en el origen y se puede dibujar una superposición de las k superficies sobre cada corte del cubo, lo que implica fijar uno de los factores de control (ver Figura 2.1). En esta figura se fija el factor X_{2} en x_{20} y los contornos que vemos sobre esta lámina son en realidad un corte de los contornos originales que están en tres dimensiones. Con más de tres factores la situación se complica, pero se preserva la misma idea.

Figura 2.1: Aspecto de los contornos de la hipersuperficie en X_{2} = X_{20}

El método gráfico de optimización simultánea se puede aplicar en el paquete Design Expert, el cual a partir de los datos experimentales, ajusta los k modelos y es posible sobreponer las superficies sobre diferentes cortes bidimensionales de un cubo o hipercubo. Por ejemplo, en el caso de tres factores, el cubo se puede rebanar en el sentido de los planos (X_{1}, X_{2}), (X_{1}, X_{3}) y (X_{2}, X_{3}). En cada una de estas formas de rebanar al cubo, la ubicación específica de un corte está dada por el valor en que se fija el tercer factor. Así, por ejemplo, en la Figura 2.1 se dibuja la lámina paralela al plano (X_{1}, X_{3}) que pasa por el punto X_{2} = x_{20}, y resulta un corte diferente en cada valor de X_{2}; esto mismo pasa con los cortes en los otros dos sentidos.

Ahora imaginemos el punto óptimo global como una pequeña bola tridimensional en algún lugar del cubo. El problema de optimizar con el método gráfico consiste en encontrar las tres láminas que cortan dicha bola por su centro, en las tres direcciones posibles. El centro de la bola es el punto óptimo simultáneo. Note que en cada hoja se obtiene una circunferencia al rebanar la bola y la lámina que pasa por el centro tiene la circunferencia más grande.

Para el caso de k = 3 factores el método es como sigue:

1. para cada variable de respuesta y en cada posible corte del cubo, se dibujan las dos curvas de nivel que corresponden a las especificaciones. El software señala (dejando sin sombrear) la subregión donde todas las respuestas predicen valores factibles, es decir, adentro de las especificaciones. Existen varias posibilidades: si casi en cualquier corte existe región factible, es indicio de que las k restricciones se pueden cumplir fácilmente y se procede a estrechar las especificaciones. Si por el contrario, no se puede encontrar ningún corte donde exista región factible, entonces el problema no tiene solución para esas restricciones, y significa que al menos una de las variables tiene que tomar valores afuera de especificaciones para que las restantes sí puedan cumplir las suyas. Una tercera posibilidad es que se descubran varias regiones factibles separadas, en cuyo caso es probable que la más grande de ellas nos arroje los mejores resultados, pero también se pueden hacer consideraciones de costos de operación y de variabilidad para elegir una región.

2. Una vez encontrada una región factible, se localizan los planos que la cortan por su centro, que son aquellos donde la región sin sombrear es más grande.

3. Si ya se ubicó una región donde todas las variables cumplen los requerimientos, lo que procede es ir estrechando las especificaciones, comenzando por las variables más importantes, y para cada conjunto nuevo de valores se ubican otra vez los mejores cortes posibles por el centro de la región factible, que ahora es más pequeña. Así se sigue, hasta tener una región factible suficientemente pequeña cuyo centro es el óptimo simultáneo buscado. En todo este proceso se trabaja con las tres figuras en dos dimensiones, nunca con el cubo completo.

Muchos autores (ver por ejemplo, Del Castillo (1996)) han criticado la dificultad de aplicar el método gráfico cuando se tienen más de dos factores de control. En este trabajo se muestra que el método es bastante fácil de aplicar hasta con tres factores; con cuatro es el doble de difícil pero todavía es posible. Además, la gran mayoría de los casos prácticos de optimización se hacen con cuatro o menos factores de control, siendo el caso de tres el más frecuente. Por otra parte, con el método gráfico es posible detectar la presencia la presencia de óptimos locales o de varias regiones con la misma factibilidad, lo cual puede pasar desapercibido con los otros métodos.

3. Método de la función de deseabilidad

Este método fue propuesto originalmente por Harrington (1965) y posteriormente mejorado por Derringer y Suich (1980) y Derringer (1994). Consiste en definir una función en el espacio de factores que mide la deseabilidad global (DG) del producto predicha en cada punto, convirtiendo de esta forma el problema de optimización multivariado en un problema de optimización univariado. Basta maximizar dicha deseabilidad global para obtener el punto óptimo que buscamos. La deseabilidad global DG en un punto x′ = (x_{1}, x_{2}, ..., x_{p}) se define como la media geométrica de las deseabilidades individuales (d_{1}, d_{2},..., d_{k}), que a su vez son transformaciones que convierten los valores predichos de cada respuesta Ŷ_{1}, Ŷ_{2},..., Ŷ_{k} a números en el intervalo [0, 1]. Si d_{i} es igual a 1 significa que la correspondiente respuesta predicha Ŷ_{i} toma su valor máximo deseable; si d_{i} = 0 la respuesta toma un valor inaceptable y la deseabilidad global DG = 0, lo que significa que el todo el producto es inaceptable independientemente de las valores de las respuestas restantes. Esto último explica el uso de la media geométrica en la definición de la deseabilidad global DG.

Sean los k modelos ajustados dados por (1.1). Estos modelos se pueden evaluar en cualquier punto x′ = (x_{1}, x_{2}, ..., x_{p}) de la región experimental y resultan los k predichos Ŷ_{1}(x), Ŷ_{2}(x),..., Ŷ_{k}(x). Cada Ŷ_{i}(x) se transforma en un valor d_{i}(x) que cae en el intervalo [0, 1] y que mide el grado de deseabilidad de la respuesta en relación a las especificaciones y valor objetivo. Suponiendo que la variable Y_{i} tiene por especificaciones inferior y superior a LIE_{i} y LSE_{i}, y su valor objetivo o nominal es T_{i}, se define la transformación d_{i} como

d_{i}(x) = \{ \begin{matrix}[\frac{Ŷ_{i}(x)−LIE_{i}}{T_{i}−LIE_{i}} ]^{s}  LIE_{i} ≤ Ŷ_{i}(x) ≤ T_{i} \\ [\frac{Ŷ_{i}(x)−LSE_{i}}{T_{i}−LSE_{i}} ]^{t}  T_{i} < Ŷ_{i}(x) ≤ LSE_{i} \\ 0  Ŷ_{i}(x) < LIE_{i} ó Ŷ_{i}(x) > LSE_{i} \end{matrix}  (3.1)

donde s y t son exponentes que sirven para elegir la forma deseada de la transformación y con ello reflejar los deseos del experimentador: si se toman grandes (digamos s, t ≥ 10) significa que la deseabilidad d_{i} sólo toma valores grandes cuando Ŷ_{i} cae cerca de su valor objetivo; si se toman valores pequeños para s y t (s, t ≤ 0.1) significa que cualquier valor de Ŷ_{i} adentro del intervalo [LIE_{i}, LSE_{i}] es igualmente deseable; cuando no se tiene idea de grados de deseabilidad, se recomienda asignar el valor 1 a los exponentes, lo que sugiere un incremento lineal de la deseabilidad hacia el valor objetivo (ver Figura 3.1). Si la respuesta tiene especificaciones de un solo lado, lo que se hace es tomar el valor objetivo (T_{i}) igual al valor a partir del cual se considera que no hay ganancia adicional en la calidad de la respuesta. Es decir, en la función dada por (3.1) desaparece una de las restricciones y la Figura 3.1 se reduce a uno de los lados respecto al valor objetivo T_{i}.

Figura 3.1: Función de deseabilidad

Una vez calculadas las k deseabilidades individuales, se obtiene la deseabilidad global DG como la media geométrica

DG = (d_{1}^{W_{1}} × d_{2}^{W_{2}} × ⋯ × d_{k}^{W_{k}})^{1/ ∑ W_{i}}  (3.2)

donde los pesos W_{i} son constantes que permiten balancear la importancia relativa de cada variable; mientras más grande es el peso dado a una variable en relación a las restantes, más grande será su deseabilidad sobre el punto óptimo global. Si todas son igualmente importantes, W_{i} = 1 para i = 1, 2, ..., k. Note que los exponentes s y t definidos en la transformación (3.1) se pueden introducir como parte de los pesos W_{i}. El punto óptimo simultáneo es el punto x′_{0} = (X_{10}, X_{20}, ..., X_{p0}) sobre el cual la función DG es máxima (ver Derringer (1994)). Cabe decir que este método no toma en cuenta la aleatoriedad de Ŷ_{i}(x), ni la calidad de los modelos, ni la habilidad de las respuestas para cumplir con las tolerancias. Esta información puede servir para elegir pesos W_{i} más adecuados, como se propone más adelante.

En la maximización de la función dada por (3.2) puede haber problemas al utilizar un algoritmo de optimización que calcula derivadas, puesto que las funciones d_{i} que la componen (ver Figura 3.1) no son diferenciables en el valor objetivo T_{i}. Del Castillo et al (1996) proponen una modificación a las funciones de deseabilidad individuales dadas por (3.2) para hacerlas diferenciables en el punto T_{i} y así tener la opción de aplicar algoritmos basados en derivadas.

4. Método de la función de pérdida

En lugar de minimizar una función de deseabilidad, Ames et al (1997) proponen minimizar una función de pérdida global definida en cada punto x como la suma ponderada de las desviaciones cuadráticas de los predichos respecto al valor objetivo. Supongamos que se tienen los modelos ajustados para k variables de interés dados por (1.1), la función de pérdida a minimizar se define por

QLP = ∑_{i=1}^{k} W_{i}(Ŷ_{i}(x) − T_{i})^{2}  (4.1)

donde W_{i} es el peso que pondera la importancia relativa de la variable Y_{i} en la función; T_{i} es el valor objetivo de la variable y Ŷ_{i}(x) es el valor predicho por el modelo en el punto x′ = (x_{1}, x_{2}, ..., x_{p}). Los valores de las constantes W_{i} se eligen en primera instancia de tal manera que se elimine el efecto de las diferentes unidades de medición que tienen las variables, además de que se puede incluir allí mismo un peso adicional que balancee su importancia relativa.

Los autores sugieren que los valores nominales T_{i} en la función objetivo (4.1) se pueden mover arbitrariamente para lograr mejores resultados de su método o soluciones alternativas, en el caso de que todas las respuestas buscan máximo o mínimo, pero no explican cómo elegir dichos valores objetivo de tal manera que el punto óptimo resultante sea el más razonable. Además, cuando en algunas respuestas los valores nominales T_{i} están dados como el punto medio de las especificaciones, no es claro como conservar el balance en la importancia relativa de las variables. Más adelante mostramos que los resultados que obtienen para uno de los ejemplos que presentan, con respuestas de ambos tipos, son demasiado inconsistentes por no fijar los valores objetivo de manera adecuada.

5. Método de la función de distancia

Este método es diferente a los anteriores en al menos cuatro aspectos: es el más formal de todos, involucra menos el sentido común del experimentador, exige que las respuestas sean linealmente independientes y puede tomar en cuenta la aleatoriedad de las respuestas predichas por cada modelo. El hecho de que suponga que las variables de respuesta tienen que ser independientes es bastante fuerte, puesto que frecuentemente existe dependencia entre ellas. Al reducir el conjunto de variables a optimizar, el resultado ignora la información de las variables eliminadas y puede sesgarse fuertemente; el hecho de conocer la correlación de una variable eliminada con alguna optimizada no implica que la primera pueda también optimizarse.

Propuesto originalmente por Khuri y Conlon (1981), el método consiste en minimizar una función de distancia entre el vector de valores óptimos individuales y el vector de predichos. En un primer paso se determina cuáles de las variables son estadísticamente independientes para trabajar solo ese grupo. Cada variable seleccionada se optimiza individualmente utilizando técnicas de superficie de respuesta univariada (ver por ejemplo Khuri y Cornell (1997) y Cornell (1990)). Sea ϕ′ = (ϕ_{1}, ϕ_{2}, ..., ϕ_{k}) el vector de óptimos individuales de las k variables independientes y sea otra vez x′ = (X_{1}, X_{2}, ..., X_{p}) un punto cualquiera en el espacio de factores de control. Los modelos dados por (1.1) se pueden escribir en forma matricial como

Ŷ_{i}(x) = z′(x)β̂_{i} (i = 1, 2, ..., k)

donde z′(x)_{(1×r)} = (1, X_{1}, X_{2},..., X_{p}, X_{1}X_{2}, X_{1}X_{3},..., X_{p−1}X_{p}, X_{1}^{2}, X_{2}^{2},..., X_{p}^{2}) es el vector renglón de variables en el modelo, que tiene dimensión r = \frac{p(p+3)}{2} + 1, o sea el número de términos en cada modelo y β̂_{i} es el vector columna de parámetros estimados. En cada punto del espacio de factores se puede calcular el vector de predichos de las k respuestas dado por Ŷ(x) = (Ŷ_{1}(x), Ŷ_{2}(x),..., Ŷ_{k}(x))′. Se define la función de distancia entre Ŷ(x) y ϕ como

ρ_{1}[Ŷ(x), ϕ] = [\frac{(Ŷ(x) − ϕ)′ˆ∑^{−1}(Ŷ(x) − ϕ)}{z′(x)(X′X)^{−1}z(x)}]^{1/2} ,  (5.1)

donde ˆ∑ es la matriz de varianzas y covarianzas estimada de las respuestas y X_{(n×r)} es la matriz de diseño aumentada con una primera columna de unos y r − (p + 1) columnas que corresponden a los términos de segundo orden; es decir, esta matriz tiene tantas columnas como términos tiene el modelo ajustado y n renglones que son número de observaciones en cada respuesta. Nótese que en el caso de una variable de respuesta esta distancia se reduce a la diferencia (distancia) entre el predicho y el valor óptimo dividida entre el error estándar del predicho. Un estimador de la matriz de varianzas y covarianzas de los predichos está dado por

ˆ∑ = \frac{Y′[I_{n} − X(X′X)^{−1}X′]Y}{n − r}

donde I_{n} es la matriz identidad de dimensión n y Y es la matriz cuyas columnas son los valores observados para cada respuesta.

Otra posible función de distancia global sugerida por Khuri y Conlon (1981) está dada por

ρ_{2}[Ŷ(x), ϕ] = [∑_{i=1}^{k} \frac{(Ŷ_{i}(x) − ϕ_{i})^{2}}{ϕ_{i}^{2}} ]^{1/2} .  (5.2)

6. Otros métodos de optimización simultánea

Además de los métodos arriba mencionados, se han propuesto otros métodos para el problema de optimizar simultáneamente más de dos respuestas. Por ejemplo, Pignatiello (1993) propone un método basado también en la función de pérdida generalizada, que a diferencia de Ames et al (1997), si minimiza una auténtica función de pérdida en el sentido de Taguchi (ver Phadke (1989)). Esto es, lo que se minimiza es una pérdida en dinero que depende de las desviaciones respecto al valor objetivo; precisamente el problema para definir dicha función a minimizar es que se debe conocer a priori cual sería la pérdida incurrida en ciertos puntos.

Del Castillo (1996) propone un método de optimización basado en técnicas de programación no lineal. En un primer paso se obtiene el óptimo de cada variable con su correspondiente región de confianza; posteriormente se elige una respuesta primaria y se optimiza sujeta a la intersección de las k regiones de confianza y a la región experimental. Así pues, el punto óptimo global optimiza a la respuesta primaria, se encuentra en la intersección de las regiones de confianza y adentro de la región experimental.

6.1. El caso de dos respuestas

El problema de optimizar al mismo tiempo dos respuestas ha recibido bastante atención en los últimos años y se han propuesto diferentes métodos para resolverlo, aunque la mayoría de las propuestas caen dentro del contexto del “problema de respuesta dual” introducido por Myers y Carter (1973). El problema de respuesta dual consiste en optimizar una variable primaria sujeta a la restricción de que una variable secundaria toma cierto valor fijo, donde ambas respuestas son explicadas por modelos de segundo orden. Vining y Myers (1990) aplican las ideas de respuesta dual al problema señalado por Taguchi sobre la necesidad de considerar, además de la media, la variabilidad de una respuesta de interés. Suponiendo que la media y la variabilidad se pueden describir en la región de interés por modelos de segundo orden, el problema es ahora optimizar la media (respuesta primaria) sujeta a que la variabilidad (respuesta secundaria) asume un valor específico. Este problema se ha resuelto considerando los tres escenarios básicos dados por Taguchi (ver Phadke (1989)) en relación al tipo de respuesta que se puede tener: “mientras más grande es mejor”, “mientras más pequeña es mejor” y “su valor objetivo es lo mejor” (ver también Del Castillo y Montgomery (1993), Lin y Tu (1995) y Copeland y Nelson (1996)).

7. Ejemplo: Optimización de neumáticos

En esta sección se analiza, utilizando los cuatro métodos presentados arriba, el ejemplo de las llantas discutido por Derringer y Suich (1980), que consiste en encontrar la combinación óptima simultánea de tres ingredientes de un compuesto de las bandas de llanta (sílice (X_{1}), silenio (X_{2}) y sulfuro (X_{3})), para cuatro variables de interés cuyos nombres y especificaciones son:

Y_{1} > 120, Índice de abrasión
Y_{2} > 1000, Módulo 200%
400 < Y_{3} < 600, Elongación

60 < Y_{4} < 75, Dureza.

Si bien las dos primeras variables no tienen límite de especificación superior, desde el punto de vista práctico se considera que no hay ninguna ganancia adicional en estas propiedades si Y_{1} y Y_{2} toman valores mayores a 170 y 1300, respectivamente. En este sentido se pueden considerar estos últimos números como los valores objetivo de dichas variables. Por su parte las variables Y_{3} y Y_{4} toman como valores objetivo el punto central de su rango de tolerancia. Se supone que las cuatro variables tienen la misma importancia desde la perspectiva de los clientes, así que en todos los métodos se les dará (o tratará de dar) el mismo peso, aproximadamente. La idea es observar el desempeño de los diferentes métodos en cuanto a su habilidad para atrapar la combinación óptima de los ingredientes.

7.1. Método Gráfico

Utilizando el paquete Design Expert se localiza la región que mejor cumple con los requerimientos de las cuatro variables de respuesta y los resultados se muestran en la Figura 7.1. Esta región se logra dibujando dos curvas de nivel de la superficie descrita por el modelo ajustado en cada respuesta, cuyos valores son

Y_{1} → 131 − 170
Y_{2} → 1192 − 1300
Y_{3} → 464 − 536
Y_{4} → 64.8 − 70.2.

Esto equivale a una reducción de 64% del tamaño del rango de especificación original, excepto en la variable Y_{1} donde solo se reduce 22%. Esto es, como las variables son igualmente importantes se trataba de reducir sus rangos en la misma proporción para sacrificar a todas por igual. Sin embargo, no fue posible reducir más el rango de la primera variable sin perjudicar fuertemente el desempeño de las otras tres, en particular el de la tercera. En otras palabras, en este ejemplo es necesario sacrificar a la variable Y_{1} en aras de un mejor desempeño de las demás variables; este hecho también lo registran los otros métodos. Los resultados de una solución gráfica óptima bajo las condiciones arriba descritas se observan en la Figura 7.1 y se resumen en la Tabla 7.1. El punto óptimo simultáneo gráfico es la combinación de niveles dada por (x_{10}=-0.217, x_{20}=0.400, x_{30}=-0.725). En la Tabla 7.1 también se muestran los valores predichos para cada respuesta sobre el punto óptimo.

Figura 7.1: Optimización simultánea por método gráfico 13 Tabla 7.1. Resultados de optimización por el método gráfico

\begin{matrix}Ŷ_{1}(x_{0}) & Ŷ_{2}(x_{0}) & Ŷ_{3}(x_{0}) & Ŷ_{4}(x_{0}) & x_{10} & x_{20} & x_{30} \\ 131.61 & 1257.33 & 466.28 & 69.91 & –0.217 & 0.400 & –0.725 \end{matrix}

7.2. Método de la función de deseabilidad

Como se ha mencionado, Derringer y Suich (1980) (de aquí en adelante simplemente DS) ilustran su método precisamente con este ejemplo de los neumáticos. Ellos reportan como solución óptima el punto x′_{0}=(-0.050, 0.145, -0.868), sobre el que se predicen los valores Ŷ(x_{0})=(129.5, 1300, 465.5, 68.5)′, cuyas deseabilidades individuales son, respectivamente d̂(Ŷ(x_{0})) =(0.189, 1.000, 0.656, 0.932)′. Para obtener estos resultados definen las cuatro funciones de deseabilidad con exponentes unitarios y pesos iguales. Nótese que la variable más sacrificada es Y_{1}, al encontrarse a un nivel de 18.9% de deseabilidad, contra la Y_{2} que se predice en un 100% de su deseabilidad.

Este resultado, reportado como la solución DS en la Tabla 7.2, considera pesos iguales para las cuatro respuestas. Si se ponderan de otra manera las variables es posible obtener un mejor valor para la primera variable, manteniendo casi la misma deseabilidad de las otras tres variables. Por ejemplo, ponderando el doble las tres primeras variables en relación a la cuarta, la solución resultante predice los valores Ŷ(x_{0})=(131.13, 1300.00, 465.68, 69.44)′, donde se obtiene un mejor valor para la primera variable, a costa de empeorar un poco el valor de la tercera. Sin embargo, sabemos que es impráctico que el experimentador pruebe una serie de ponderaciones arbitrarias para quedarse con la que considera mejor. Así pues, se sugiere ponderar las variables con los pesos W_{i} = 1/C_{pmi} (ver Apéndice), que consideran la desviación de la respuesta al valor objetivo, así como su variabilidad. Estos pesos arrojan los resultados de la solución DS2 de la Tabla 7.2. Se observa una mejora en la variable Y_{1} a expensas de la variable Y_{4}.

Tabla 7.2. Resultados del método de la función de deseabilidad

Solución Ŷ_{1}(x_{0}) Ŷ_{2}(x_{0}) Ŷ_{3}(x_{0}) Ŷ_{4}(x_{0}) x_{10} x_{20} x_{30}
DS  129.50 1300.00 465.50 68.50 -0.050 0.145 -0.868
DS1 131.13 1300.00 465.68 69.44 -0.126 0.404 -0.843
DS2 130.51 1300.00 469.95 69.55 -0.151 0.424 -0.871

En resumen, si a priori las variables tienen la misma importancia relativa, parecería razonable ponderarlas parejo a todas. Pero lo cierto es que cuando las variables tienen diferente variabilidad y distinta desviación al valor objetivo, la ponderación unitaria no es la más adecuada.

7.3. Método de la función de pérdida

Ames et al (1997) también ilustran su método utilizando el ejemplo de los neumáticos de DS, sin embargo el punto óptimo que reportan no es el mejor que se puede encontrar, debido en nuestra opinión a dos aspectos: primero, la función no está bien definida, en el sentido de que toma como valores objetivo de las dos primeras variables el mínimo de sus valores deseables, cuando los autores originales comentan claramente que lo más deseable es tener valores máximos de estas dos variables (DS (1980)). Segundo, la manera de ponderar las variables no es la más adecuada. Ellos definen la función global a minimizar como

QLP = 0.0318U[120, Ŷ_{1}(x)]{Ŷ_{1}(x) − 120}^{2}  (7.1)
+0.00000925U[1000, Ŷ_{2}(x)]{Ŷ_{2}(x) − 1000}^{2}
+0.00237{Ŷ_{3}(x) − 500}^{2} + 0.62{Ŷ_{4}(x) − 67.5}^{2},

donde U[a, b] = 1 si b ≥ a; U[a, b] = 0 si b < a y las constantes son pesos iguales al inverso del cuadrado medio del error (CM_{E}) del modelo correspondiente. Es debido a la función escalón U que cualquier punto que prediga valores de las primeras dos variables arriba de 120 y 1000, respectivamente, y las dos últimas variables en sus valores objetivos, es un punto óptimo que tiene una pérdida mínima QLP ≅ 0. El problema es que hay una infinidad de puntos que pueden predecir estas condiciones, los cuales se representan con la línea oscura de la Figura 7.2, y difieren bastante en los valores que predicen para las primeras respuestas Y_{1} y Y_{2}. Así pues, todos los puntos señalados en la figura minimizan la función de pérdida dada por (7.1). Da la impresión que los autores (Ames et al) no se percatan de este hecho y lo confunden con el problema de tener muchos mínimos locales, que no es lo mismo.

Figura 7.2: Familia de puntos que hacen a QLP=0

El punto óptimo que ellos reportan (ver Figura 7.2) está dado por x′_{0}=(-0.461, -0.283, -0.528) y predice como valores de las respuestas a Ŷ(x_{0})=(122.7, 1069.0, 500.0, 67.0)′. Inmediatamente se observa que los valores predichos de las dos primeras respuestas son bastante bajos, lo cual se debe, entre otras cosas, a que el algoritmo de optimización utilizado convergió “arbitrariamente” a este punto existiendo muchos otros puntos que pueden predecir mejores valores para Y_{1} y Y_{2}. Es fácil mostrar que con su función objetivo (incorrecta) dada por (7.1) se pueden obtener mejores resultados simplemente cambiando el punto inicial del algoritmo de optimización. Por ejemplo, con el algoritmo de Nelder-Mead y debido a la manera de definir la función de pérdida global, se pueden encontrar una infinidad de puntos donde la pérdida toma un valor mínimo cercano a cero. Así, tomando como punto inicial el centro de la región experimental se converge al punto óptimo x′_{0}=(-0.100, 0.213, -1.278) sobre el que se predicen las respuestas Ŷ(x_{0})=(123.23, 1409.41, 500.00, 67.5)′. Este resultado es mucho mejor en las primeras dos variables que el reportado por Ames y colaboradores.

También se puede encontrar una solución que se ajuste más a la información original de Derringer sobre las respuestas y se obtiene, por ejemplo, que considerando como punto inicial del algoritmo al vector (0.0, 0.0, -1.0) éste converge al punto óptimo simultáneo x′_{0}=(-0.211, 0.057, -1.059), donde los predichos son Ŷ(x_{0})=(124.21, 1303.95, 500.00, 67.5)′. Claramente esta última solución también es mejor que la de Ames et al (1997), ya que las tres últimas variables toman sus valores objetivo, y la primera se sacrifica todavía excesivamente, pero sigue prediciendo un mejor resultado que el de Ames.

En resumen, los puntos más arriba que la solución de Ames en el sentido de la variable X_{3} (ver Figura 7.2) son peores soluciones que ésta; del punto de Ames hacia abajo, hasta un valor de alrededor de x_{3} = −1.28, los puntos son mejores soluciones.

Función de pérdida corregida y con ponderación alternativa

En la discusión del ejemplo de los neumáticos Ames et al (1997) argumentan que DS (1980) no proporciona mayor información para ponderar las variables, sin embargo es claro que dichos autores dan el mismo peso a cada variable al tomar coeficientes r=s=1 en las funciones de deseabilidad. Así pues, en principio se trata de dar el mismo peso a las cuatro variables, cosa que Ames et al (1997) de hecho asumen al ponderar cada variable por el inverso del cuadrado medio del error (CM_{E}) de cada modelo, lo que se supone estaría reduciendo el efecto de las diferentes unidades de medición. Desgraciadamente ellos “definen incorrectamente” la función de pérdida objetivo, la cual, considerando sus ponderaciones (ver 7.1) debería ser

QLP = 0.0318U[170, Ŷ_{1}(x)]{Ŷ_{1}(x) − 170}^{2}  (7.2)
+0.00000925U[1300, Ŷ_{2}(x)]{Ŷ_{2}(x) − 1300}^{2}
+0.00237{Ŷ_{3}(x) − 500}^{2} + 0.62{Ŷ_{4}(x) − 67.5}^{2},

al colocar los objetivos correctos en las primeras dos variables. Esta función alcanza su valor mínimo en el punto x′_{0}=(0.060, 0.536, -0.545), sobre el que se predicen las respuestas Ŷ_{0}=(139.90, 1325.02, 422.72, 70.19)′. Otra vez, se observa que el inverso del CM_{E} no es una ponderación suficiente, puesto que en este caso castiga demasiado a la variable Y_{3}. Esto es, no basta eliminar el efecto de las diferentes unidades, sino que la ponderación debería dar un poco más de peso a las variables que más dificultad tienen en cumplir los requerimientos individuales.

Una ponderación que parece dar mejores resultados consiste en utilizar los pesos W_{i}, inspirados en un índice de capacidad de procesos (Spiring (1997)), que se definen en el Apéndice. Con ellos se obtiene por punto óptimo x′_{0}=(-0.037, 0.188, -0.710) cuyos predichos son Ŷ_{0}=(132.31, 1282.34, 451.45, 68.51)′. Este punto es mejor solución que las dos anteriores (ver solución AMES2 en la Tabla 7.3). En la Tabla 7.3 se resumen los resultados obtenidos por el método de la función de pérdida: AMES es la solución original, AMES1 es la solución obtenida con la función corregida y AMES2 es la que resulta de la función corregida y ponderada por los pesos dados en el Apéndice.

Tabla 7.3. Resultados del método de la función de pérdida

solución Ŷ_{1} Ŷ_{2} Ŷ_{3} Ŷ_{4} x_{10} x_{20} x_{30}
AMES 122.70 1069.00 500.00 67.00 -0.461 -0.283 -0.528
AMES1 139.90 1325.02 422.72 70.19 0.060 0.536 -0.545
AMES2 132.31 1282.34 451.45 68.51 -0.037 0.188 -0.710

7.4. Método de la función de distancia

Después de elegir el grupo de variables de respuesta independientes que se van a optimizar, el siguiente paso del método de la función de distancia es optimizar cada variable individualmente utilizando técnicas de superficie de respuesta univariada. Del análisis canónico se encuentra que los cuatro puntos estacionarios de las respuestas son puntos silla, por lo que es necesario utilizar la técnica de análisis de cordillera con cada respuesta para encontrar los mejores puntos “óptimos” dentro de la región experimental (ver Myers y Montgomery (1995)). Haciendo los cálculos con el procedimiento RSREG de SAS, se obtienen los cuatro óptimos individuales (x_{0i}; i = 1, 2, 3, 4) que se muestran en la Tabla 7.4, donde también se reportan los correspondientes predichos (ϕ_{i}; i = 1, 2, 3, 4) sobre ellos. Estos puntos se ubican sobre la esfera de radio 1.7

Es importante hacer notar que el punto óptimo de la variable Y_{1} está en el lado opuesto del punto óptimo de la variable Y_{3}, es decir, la dirección óptima para una de estas variables es prácticamente la peor dirección para la otra. Esta contraposición entre las variables Y_{1} y Y_{3} afecta a cualquier método de optimización simultánea, por lo que el logro de mejores valores predichos de una de estas variables necesariamente se hace a expensas de peores valores de la otra.

Aplicando dos de las funciones de distancia sugeridas por Khuri y Conlon (1981), se obtienen los resultados de la Tabla 7.5. La solución KCd1 se refiere a la función de distancia ρ_{1} definida en la relación (5.1), mientras que KCd2 es la distancia ρ_{2} dada en (5.2), restringiendo en ambos casos la región a las condiciones dadas por DS (1980). La solución KCd1 castiga demasiado a la variable Ŷ_{3}, mientras que la KCd2 reporta un valor muy bajo para Ŷ_{1}. Así, las dos soluciones son malas.

Tabla 7.4. Puntos óptimos individuales

predicho (ϕ_{i})  punto óptimo
respuesta  Ŷ (x_{0i})  x_{10}  x_{20}  x_{30}
Y_{1}  194.32  0.933 1.024 0.986
Y_{2}  2334.51  0.607 0.543 1.492
Y_{3}  499.75  −0.571 −0.226 −0.394
Y_{4}  67.49  0.053 −0.272 −0.115

Parte del problema de las primeras dos soluciones es el valor ϕ_{2} = 2334.51 (ver Tabla 7.4) tan grande, que al utilizarlo como el valor objetivo de la variable Ŷ_{2}, acaba dándole más peso a esta variable en perjuicio de la respuesta Y_{3}. Buscando resultados más razonables, decidimos utilizar los valores nominales sugeridos por DS para las variables Y_{1} y Y_{2}, y es así como se llega a las soluciones KCa1 y KCa2 listadas en la Tabla 7.5. Efectivamente, ahora la variable Y_{2} tiene menos peso que antes, y se logra un mejor balance entre las cuatro respuestas, en particular con la solución KCa1 que es la mejor.

Tabla 7.5 Resultados con las distancias 1 y 2

solución Ŷ_{1} Ŷ_{2} Ŷ_{3} Ŷ_{4} x_{10} x_{20} x_{30}
KCd1 134.80 1460.01 426.35 69.65 0.452 0.864 -1.137
KCd2 133.08 1425.99 401.67 67.17 0.693 0.143 -1.000
KCa1 127.41 1334.91 487.55 71.24 -0.179 0.864 -1.137
KCa2 139.00 1250.00 419.22 70.31 -0.181 0.279 -0.121

7.5. Discusión

En la Tabla 7.6 se resumen los mejores resultados obtenidos con cada uno de los cuatro métodos presentados arriba para el problema de los neumáticos, aplicándolos tal como se describen originalmente. El método DS es el de la función de deseabilidad de Derringer y Suich (1980); las etiquetas KCa1 y KCa2 se refieren al método de la función de distancia considerando las dos funciones de distancia sugeridas por Khuri y Conlon (1981) y restrigiendo la región a las condiciones de DS; AMES es la solución dada por Ames et al (1997) y AMES1 es la solución que resulta al corregir los valores nominales en la función de pérdida.

Las mejores soluciones son las obtenidas con los métodos de la función de deseabilidad (DS) y el gráfico (Tabla 7.6), quedando muy parejas en el porcentaje de mediciones afuera de especificaciones que se espera tener en el futuro sobre el punto óptimo correspondiente. El método gráfico tiene la ventaja de ser más flexible puesto que se trabaja con las superficies de respuesta originales y en la escala original. Las peores soluciones son las llamadas KCa2 y AMES, puesto que reportan predichos demasiado alejados de los valores objetivos en las variables Y_{3} y Y_{2}, respectivamente (Tabla 7.6).

Ninguno de los métodos considera, además de la variabilidad y la media de Ŷ_{i}, su habilidad para cumplir con las tolerancias, de aquí que todas las soluciones dadas en la Tabla 7.6 sean susceptibles de mejora si se encuentra la manera de involucrar los tres aspectos. Es decir, algunos de los métodos utilizan solo la media, otros solo la media y la variabilidad, alguno más solo la media y las especificaciones, pero ninguno las tres cosas. Como resultado de esto, en los puntos óptimos simultáneos reportados se predice un excesivo número artículos defectuosos en el futuro, en términos de algunas de las variables. La idea entonces es determinar puntos óptimos que predigan un mínimo de defectuosos en las mediciones futuras de todas las respuestas.

Tabla 7.6. Resultados con los métodos originales

punto óptimo  predichos
método x_{10} x_{20} x_{30} Ŷ_{1} Ŷ_{2} Ŷ_{3} Ŷ_{4}
Gráfico -0.217 0.400 -0.725 131.61 1259.73 466.28 69.91
DS  -0.540 0.150 -0.869 129.43 1300.00 466.05 68.03
AMES -0.460 -0.284 -0.525 122.73 1069.10 500.00 67.50
AMES1 0.060 0.536 -0.545 139.90 1325.02 422.72 70.19
KCa1 -0.179 0.864 -1.137 127.41 1334.91 487.55 71.24
KCa2 -0.187 0.279 -0.121 139.00 1250.00 419.21 70.31

En otras palabras, no basta cumplir los requerimientos, sino que para evitar artículos defectuosos idealmente se debe cumplir que el intervalo de confianza para la predicción futura sobre el punto óptimo se encuentre completamente adentro de las especificaciones y no solo parcialmente como ocurre en algunas variables. Por ejemplo, la solución DS predice un valor de 129.43 para la primera variable, pero dado su error estándar, el límite inferior del intervalo del pronóstico es de alrededor de 113, es decir, se predice un porcentaje importante de defectivos. Lo mismo pasa con la segunda variable en esta solución, que aunque toma su valor óptimo de 1300 en promedio, el error estándar de las futuras mediciones es de alrededor de 420, por lo que también se tendrán defectuosos por este concepto. La pregunta es: ?‘cómo obtener el punto óptimo simultáneo con cada método, que además tenga el mejor centrado del intervalo de predicción adentro de las especificaciones?.

Soluciones alternativas

Buscando tomar en cuenta lo antes dicho, en los métodos de la función de deseabilidad y en el de la función pérdida, se proponen otras maneras de ponderar las funciones respectivas (Secciones 7.2 y 7.3). Las soluciones alternativas que resultan se resumen en la Tabla 7.7. La solución gráfica alternativa se obtiene utilizando mayor resolución (aumentos) en las gráficas y aplicando el criterio adicional de considerar el error estándar de los predichos. En nuestra opinión es la mejor solución, puesto que predice menos defectivos en las mediciones futuras de las cuatro variables. El hecho de permitir a la variable Ŷ_{2} tomar valores por arriba de 1300 hace posible que la solución gráfica sea la mejor; si se limita a la cota de 1300, los resultados serían más parecidos a las otras dos soluciones alternativas (Tabla 7.7).

Una cualidad importante del método gráfico es que se explora la hipersuperficie y se detecta la presencia de varias regiones factibles cuando estas existen, habilidad que no tienen los otros métodos.

El método alternativo DS2 consiste en ponderar las deseabilidades con los pesos W_{i}=1/C_{pmi}, lo cual mejora un poco los resultados originales (ver Tablas 7.6 y 7.7). La solución AMES2 alternativa consiste en una modificación a los pesos sugeridos por los autores del método y con ello (ver Tabla 7.7) se logra un resultado que compite con el DS2. Pensamos que se pueden mejorar un poco más estas dos soluciones si las especificaciones dadas por DS para las primeras dos variables se modifican considerando el error estándar del pronóstico de dichas variables.

Tabla 7.7. Soluciones alternativas

punto óptimo  predichos
método x_{1} x_{2} x_{3} Ŷ_{1} Ŷ_{2} Ŷ_{3} Ŷ_{4}
Gráfico 0.010 0.390 -0.890 132.00 1337.79 455.45 68.96
DS2 -0.151 0.424 -0.871 130.51 1300.00 469.95 69.55
AMES2 -0.037 0.188 -0.710 132.31 1282.34 451.45 68.51

En conclusión, las tres alternativas son mejores soluciones que las originales en cuanto a la proporción esperada de producto afuera de especificaciones. Las modificaciones que proponemos a los métodos de las funciones de deseabilidad y de pérdida deben probarse considerando otros escenarios, no solo el de este ejemplo, antes de recomendar su aplicación a cualquier problema de multirrespuesta.

8. Apéndice: pesos basados en capacidad

El índice de capacidad de procesos C_{pm} se define, de acuerdo a Spiring (1997), como

C_{pm} = \frac{T − LIE}{3\sqrt{σ^{2} + (μ − T)^{2}}} ,

donde μ es la media del proceso y T es el punto medio de las especificaciones, en el caso de las variables con doble especificación; en las variables con especificación de un solo lado, T es el punto donde la pérdida esperada es igual a cero. Este índice mide la habilidad que tiene una variable para cumplir con los requerimientos, considerando su dispersión, su sesgo respecto del valor objetivo y la amplitud de las tolerancias. La cuestión es ?‘como trasladar estas ideas a superficie de respuesta?.

Ponderación de las funciones de deseabilidad

A partir de los datos experimentales de cada variable, y del modelo ajustado, es posible calcular un número que mida el potencial de cada modelo para predecir valores adentro de las especificaciones. Considerando como estimador de la varianza al cuadrado medio del error (CM_{E}) y como estimador de la media la depth mean, se obtiene la cantidad

C_{pm_{i}} = \frac{T_{i} − LIE_{i}}{3\sqrt{CM_{E_{i}} + (μ̂_{i} − T_{i})^{2}}}  (8.1)

para cada una de las respuestas, donde T_{i} y LIE_{i} se toman de acuerdo a Derringer y Suich (1980). Como la idea es apoyar a las variables con mayor dificultad para cumplir los requerimientos, los pesos que se sugieren en el método de la función de deseabilidad están dados por

W_{i} = \frac{1}{C_{pmi}}; i = 1, 2, 3, 4.  (8.2)

Ponderación de la función de pérdida

Recordemos que Ames et al (1997) utilizan los pesos W_{i} = 1/CM_{E_{i}} en el ejemplo de los neumáticos, que en principio busca eliminar el efecto de las diferentes unidades de las respuestas, y eso pasa cuando la varianza es proporcional a las magnitudes de los datos, lo cual de hecho no ocurre en este ejemplo. Además estos pesos tampoco toman en cuenta el sesgo de los datos respecto a los valores nominales. Buscando subsanar estas limitantes, los pesos propuestos para el método de la función de pérdida son los valores

W_{i} = \frac{|T_{i} − μ̂_{i}| /\sqrt{CM_{E_{i}}}}{CM_{E_{i}} + (μ̂_{i} − T_{i})^{2}} ; i = 1, 2, 3, 4.  (8.3)

9. Referencias

AMES, A. E., MATTUCCI, M., STEPHEN, M., SZONYI, G. and HAWKINS, D. M. (1997). “Quality Loss Functions for Optimization Across Multiple Response Surfaces”. Journal of Quality Technology 29, pp. 339-346.

BOX, G. E. P., HUNTER, J. S. y HUNTER, W. G. (1988). “Estadística para Investigadores”. Editorial Reverté, Barcelona, España.

COPELAND, K. A. F. and NELSON, P. R. (1996). “Dual Response Optimization via Direct Function Minimization”. Journal of Quality Technology 28, pp. 331-336.

CORNELL, J. A. (1990). How to Apply Response Surface Methodology. ASQC Press, Milwaukee, Wisconsin.

DEL CASTILLO, E. (1996). “Multiresponse Process Optimization via Constrained Confidence Regions”. Journal of Quality Technology 28, pp. 61-70.

DEL CASTILLO, E. and MONTGOMERY, D. C. (1993). “A Nonlinear Programming Solution to the Dual Response Problem”. Journal of Quality Technology 25, pp. 199-204.

DEL CASTILLO, E., MONTGOMERY, D. C. and McCARVILLE, D. R. (1996). “Modified Desirability Functions for Multiple Response Optimization”. Journal of Quality Technology 28, pp. 337-345.

DERRINGER, G. and SUICH, R. (1980). “Simultaneous Optimization of Several Response Variables”. Journal of Quality Technology 12, pp. 214-219.

DERRINGER, G. (1994). “A Balancing Act: Optimizing a Product’s Properties”. Quality Progress pp. 51-58.

DESIG-EXPERT (1991). Stat-Ease, Version 2.02, Minneapolis, MN.

HARRINGTON, E. C. JR. (1965). “The Desirability Function”. Industrial Quality Control 21, pp. 494-498.

KHURI, A. and CONLON, M. (1981). “Simultaneus Optimization of Multiple Responses Represented by Polynomial Regression Functions”. Technometrics 23, pp.363-375.

KHURI, A, and CORNELL, J. A. (1997). Response Surfaces: Designs and Analyses. Second Edition. Marcel Dekker, New York.

LIN, D. K. J. and TU, W. (1995). “Dual Response Surface Optimization”. Journal of Quality Technology 27, pp. 34-39.

LUENBERGER, D.G. (1984). Linear and Nonlinear Programming. Reading, MA: Addison - Wesley.

MYERS, R. H. and CARTER, W. H. (1973). “Response Surface Techniques for Dual Response Systems”. Technometrics 15, pp. 301-317.

MYERS, R. H. and MONTGOMERY, D. C. (1995). Response Surface Metodology: Process and Product Optimization Using Designed Experiments. John Wiley & Sons, New York.

NELDER, J. A. and MEAD, R. (1965). “A Simplex Method for Function Minimization”. The Computer Journal 7, pp. 308-313.

PIGNATIELLO, J. J. JR. (1993). “Strategies for Robust Multiresponse Quality Engineering”. IIE Transactions 25, pp. 5-15.

SAS INSTITUTE (1990). SAS Versión 6.0, SAS Institute, Cary, NC.

SPIRING, F. A. (1997). “A Unifying Approach to Process Capability Indices”. Journal of Quality Technology 29, pp. 49-58.

VINING, G. G. and MYERS, R. H. (1990). “Combining Taguchi and Response Surface Philosophies: A Dual Response Approach”. Journal of Quality Technology 22, pp.38-45.