Stochastic approximation algorithm for industrial process optimisation

Algoritmo de aproximaciones estocásticas para la optimización de procesos industriales

Stochastic approximation algorithms are alternative linear search methods for optimising control systems where the functional relationship between the response variable and the controllable factors in a process and its analytical model remain unknown. These algorithms have no criteria for selecting succession measurements ensuring convergence, meaning that, when implemented in practice, they may diverge with consequent waste of resources. The objective of this research was to determine industrial processes' optimum operating conditions by using a modified stochastic approximation algorithm, where its succession measurements were validated by obtaining response variable values for each iteration through simulation. The algorithm is presented in nine stages; its first six describe which are process independent and dependent variables, the type of experimental design selected, the experiments assigned and developed and the second order models obtained. The last three stages describe how the algorithm was developed, and the optimal values of the independent variables obtained. The algorithm was validated in 3 industrial processes which it was shown to be efficient for determining independent variables' optimum operating conditions (temperature and time): the first three iterations were obtained at 66°C in 3 hours 42 minutes for process 1, unlike processes 2 and 3 where the first iteration was obtained at 66°C in 6 hours 06 minutes and 80°C in 5 hours 06 minutes, respectively.

Los algoritmos de aproximaciones estocásticas son métodos alternativos de búsqueda lineal para optimizar o controlar sistemas donde la relación funcional entre la variable de respuesta y los factores controlables de un proceso y su modelo analítico son desconocidos. En estos algoritmos no existe un criterio en la selección de sus medidas de sucesión que garanticen la convergencia, lo cual puede llevar a que al implementarlos en la práctica diverjan, con el consecuente desperdicio de recursos. El objetivo de la investigación es determinar las condiciones óptimas de operación de procesos industriales mediante un algoritmo de aproximaciones estocásticas modificado, donde sus medidas de sucesión son validadas al obtener valores de la variable de respuesta de cada iteración mediante simulación. El algoritmo es presentado en nueve etapas. En sus primeras seis se describen cuáles son las variables independientes y dependientes del proceso, se selecciona la clase del diseño experimental, se asignan y desarrollan los experimentos y se obtienen los modelos de segundo orden; en las últimas tres etapas se desarrolla el algoritmo y se obtienen los valores óptimos de las variables independientes. El algoritmo se validó en tres procesos industriales, demostrándose que es eficiente para determinar las condiciones óptimas de operación de las variables independientes (temperatura y tiempo); en el proceso 1 se obtienen en las primeras tres iteraciones en 66 °C y 3 h 42 min, a diferencia de los procesos 2 y 3, que se obtienen en la primera iteración con 66 °C y 6 h 06 min y 80 ° C y 5 h 06 min, respectivamente.

• HTML (es)
• HTML (en)

Jesús Everardo Olguín Tiznado¹, Rafael García Martínez², Claudia Camargo Wilson³, Juan Andrés López Barreras⁴

¹ Ingeniero Industrial, Instituto Tecnológico de Huatabampo. Maestro en Ciencias en Ingeniería Industrial, Instituto Tecnológico de Hermosillo. Profesor e Investigador, Universidad Autónoma de Baja California. México. jeol79@uabc.edu.mx

² Licenciado en Matemáticas, Universidad Autónoma de Nuevo León. Maestro y Doctor en Ciencias en Ingeniería Industrial, Instituto Tecnológico de Ciudad Juarez. Director del Instituto Tecnológico del Valle del Yaquí. México. ra_garcia@ith.mx

³ Ingeniero Industrial, Instituto Tecnológico de Los Mochis. Maestro en Ciencias en Ingeniería Industrial, Instituto Tecnológico de Hermosillo. Profesor e investigador, Universidad Autónoma de Baja California. México. ccamargo@uabc.edu,mx"

⁴ Ingeniero Industrial, Instituto Tecnológico de Huatabampo. Maestro en Ciencias en ingeniería Industrial, Instituto Tecnológico de Tijuana. Profesor e investigador, Universidad Autónoma de Baja California. jalopez@uabc.edu.mx

---

RESUMEN

Keywords: algoritmos de aproximaciones estocásticas, variables independientes, variables dependientes, proceso iterativo, simulación.

---

ABSTRACT

Palabras clave: stochastic approximation algorithm, dependent variable, independent variable, iterative process, simulation

---

Recibido: octubre 22 de 2010 Aceptado: noviembre 15 de 2011

---

Introducción

El método de aproximaciones estocásticas presentado por Robbins y Monro (1951) es un método de búsqueda lineal de la raíz de la función desconocida    que representa al valor esperado de una variable aleatoria. Kiefer y Wolfowitz (1952) lo modifican para que pueda ser usado en la determinación del óptimo de  . Blum (1954) extiende los resultados de los autores anteriores a espacios cartesianos de dimensión mayor que 1.

A partir del trabajo presentado por Blum (1954) se da un incremento en la cantidad de métodos de aproximaciones estocásticas (Kushner y Clark, 1978; Polyak, 1991; Polyak y Juditsky, 1992; Andradóttir, 1995 (i, ii); Delyon, 1996; Kulkarni y Horn, 1996; Maeda, 1996). Pero Andradóttir (1996) asegura que todos estos métodos son procedimientos sin un criterio teórico de determinación, usados para determina ^X* en  , de tal forma que  , donde   es la función que corresponde al vector gradiente de la función  , de la cual se desconoce su expresión analítica, pero es posible cuantificar su valor para una combinación específica de valores o niveles de los factores controlables, medición que está sujeta a un error experimental del que no se establece ningún supuesto en cuanto a su distribución de probabilidad.

Chin (1997) clasifica los procedimientos de aproximaciones estocásticas en dos tipos: el de Robbins-Monro y el tipo KieferWolfowitz. Los primeros se caracterizan por requerir las observaciones directas de ^h, dentro de los cuales se encuentran: los métodos de Robbins-Monro, pasos ascendentes, NewtonRaphson, análisis de perturbación y tasa de verosimilitud, mientras que los segundos requieren estimaciones o aproximaciones ^h como lo son: Kiefer-Wolfowitz, diferencias finitas, método de direcciones aleatorias, el método escalado y el algoritmo estocástico de perturbación simultánea). Estos últimos los considera más útiles, dado que no requieren un conocimiento profundo del sistema a optimizar, es decir, son aplicables en situaciones en las cuales se desconoce la relación funcional entre la variable de respuesta denotada como ^yi y los ^d factores controlables denotados por el vector   ( representa el espacio cartesiano de dimensión ^d ), situación que se presenta con mayor frecuencia en la práctica. Fu y Ho (1988) y Chin (1997) señalan al algoritmo estocástico de perturbaciones simultáneas como el más eficiente, tanto teórica como prácticamente, ya que presenta mayor tasa de convergencia y requiere de un menor número de observaciones en cada iteración; esta última resulta de gran interés pues de ella depende en forma directamente proporcional el costo económico y la sencillez del trabajo experimental.

El presente trabajo de investigación propone un algoritmo de aproximaciones estocásticas modificado en el que las sucesiones   se validan con simulación y aplicación real mediante modelos de segundo orden al obtener los valores óptimos de operación en las variables de respuesta que intervienen en tres procesos industriales similares, esto debido a un pobre control en las condiciones de operación en dichos procesos, generando con esto desperdicios de productos. Además, demostrar que los algoritmos de aproximaciones estocásticas con perturbación simultánea (AAEPS) son eficientes al trabajar con modelos de segundo orden, como el que se muestra en la ecuación 1:

Desarrollo experimental

Los materiales utilizados para el desarrollo y validación de este proyecto de investigación son: una computadora Pentium PC, portege R200, M, procesador 1,2 GHz y 598 MHz en RAM. Los software para el análisis estadístico de los datos son: Statistica^® y Matlab^®.

A continuación describimos el método utilizado en la obtención de la información requerida para el análisis experimental a los fines de determinar las condiciones óptimas de operación de tres procesos industriales analizados, de los cuales se ejemplifica con los datos del proceso 1.

Primero, se elabora una lista de las variables independientes significativas o factores controlables, incluyendo sus rangos; denotados por el vector X=(X1,...Xd) ε R^d donde R^d representa el espacio cartesiano de dimensión d. En esta investigación las variables significativas utilizadas son: X1 representa la temperatura, con un rango de inicio de 60 a 70 °C, y X2 representa el tiempo, con un rango de inicio de 4 a 5 horas en los tres procesos evaluados. En el análisis de los procesos se tiene el inicio para las variables independientes en: X1= 65°C y para X2 = 4 h 30 min de proceso.

Segundo, se elabora una lista de las variables dependientes (respuestas) y sus unidades, denotadas por el vector (nx1) observaciones Y=(Y1,...,Yd) ε R^d. Se listaron las variables de respuestas que darán sustento a la investigación, siendo estas: ^y1 representa la respuesta 1 de los procesos; su unidad está dada en porcentaje; ^y2 representa la respuesta 2 de los procesos y su unidad está dada en grados. Los valores nominales o meta que se persigue obtener para las respuestas son de 5% (Target) de humedad final para la variable ^y1 y 80 (Target) en escala de color de la Hunter Laboratories para la variable ^y2. Esto con la finalidad de que los procesos cumplan con los requerimientos que solicita el cliente ( de humedad final en el producto y escala de color).

Tercero, se seleccionó una clase de diseño experimental. En este caso el diseño es generado de un factorial 3^k, que significa k factores a tres niveles de análisis experimental (Montgomery, 2009); para este trabajo será de dos factores (temperatura y tiempo) a tres niveles (60, 65 y 70 °C en la variable de temperatura, y 4 h, 4 h 30 min y 5 h para la variable del tiempo). En el paso 1 se mencionó cómo se trabajará en este experimento con los valores iniciales de estas variables por cada proceso.

Cuarto, se asignaron los experimentos aleatoriamente. En cada una de las etapas experimentales las corridas se hicieron aleatoriamente.

Quinto, desarrollo de los experimentos y recopilación de los datos. Se realizaron cinco réplicas en el experimento, bajo las condiciones de los valores iniciales de las variables independientes de cada proceso con respecto a sus valores meta de cada variable de respuesta. Esto, con el fin de obtener la media (µ) y su respectiva desviación estándar (σ).

Sexto, ya recabados los datos se procede a obtener los modelos de segundo orden en las variables de respuestas y₁ y y₂ para su medida (µ) y su desviación estandar (σ) por cada proceso industrial. Por ejemplo, las ecuaciones de regresión de segundo orden (  ) y (  ) en el proceso industrial 1, son:

Séptimo, el algoritmo estocástico de perturbación simultánea (AEPS) se calcula de acuerdo con los siguientes pasos, según Spall J. C. (1998):

Paso 1. Inicialización y coeficiente de selección. Seleccione el índice con tador k=1. Tome un valor supuesto del vector gradiente ^θ0 y los coeficientes de no negatividad , ^c, ^A, y . Delyon (1996), Spall (2003) y Chien y Luo (2008) establecen que el valor que asume típicamente para y para cuando el vector gradiente es igual a la media aritmética de ^m estimaciones. Valores en la práctica efectivos y teóricamente válidos para y son ^0.602 y ^0.101 respectivamente (los valores óptimos asintóticos de ^1.0 y ^1/6 pueden ser usados también); los valores de , ^c, ^A pueden ser determinados como se mostrará más adelante. Una guía útil al seleccionar ^A es hacerlo como si fuera mucho menor que el máximo número de iteraciones permitidas o esperadas, es por ello que se seleccionó ^A=100, ^=0.16, ^c=1. Los resultados obtenidos con base en los datos mostrados son, para    y  .

Paso 2. Generación del vector de perturbación simultánea. Generado por el método de Montecarlo, un vector de perturbación aleatorio p-dimensional , donde cada uno de los ^p componentes de son generados independientemente de una distribución de probabilidad con una media cero. Una simple (y teóricamente valida) opción para cada uno de los componentes usar una distribución Bernoulli ^±1 con probabilidad de ^1/2 para cada resultado ^±1 . Nótese que variables uniformes y normales aleatorias no son permitidas para los elementos del por las condiciones regulares del AEPS (Brooks O., 2007; Maryak y Chin, 2008). Por lo tanto, los valores del vector en los tres procesos son: ^{+ =3°C} y el de ^{- = -0.3} horas, valores dados por el experimentador para los tres procesos industriales con base en sus características.

Paso 3. Evaluaciones de la función a minimizar (3a). Se seleccionan valores iniciales para   en los que realizamos la simulación partiendo de las corridas en las cuales se trabajó, de 60, 65 y 70 °C, e incluso de otras que no se efectuaron, de 50, 55, 75 y 80 °C, con 4 h, 4h 30 min, y 5 h, además de 3 h, 3 h 30 min, 5 h 30 min y 6 h, respectivamente. En los procesos se tiene que sus valores iniciales de las variables independientes son: .

Paso (3b). Después de seleccionar se procede a sustituir los valores correspondientes a las distintas variables independientes en las ecuaciones de regresió ^yµ  y  ^yσ obtenidas anteriormente con los datos del diseño experimental ^3k mencionado en el paso tercero. Las ecuaciones de regresión de segundo orden (  )  y  (   ) representan el proceso industrial 1, como sigue:

Las ecuaciones de regresión de segundo orden (   ) y (   ) representan el proceso industrial 2, como sigue:

Las ecuaciones de regresión de tercer orden (   ) y (  ) representan el proceso industrial 3, como sigue:

Sustituyendo los valores iniciales de   y ,  en las ecuaciones de regresión de segundo orden (    ) y (   ) para el análisis del proceso industrial 1 se tiene que los valores de las ecuaciones de regresión de segundo orden son    ^=5.011 y para    ^=-1.050.

Paso (3c). Obtenidos los valores de   ,   ,    y   ,   ,    para    se procede a sustituir los valores en la ecuación (2) a los fines de calcular el error cuadrático medio (ECM) de cada proceso industrial, tal como sigue:

donde ^yµi representa la variable de respuesta para su media en el proceso i; i =1, 2, 3; T representa el valor meta u objetivo del proceso, llamado también target, y ^yσi representa la variable de respuesta para su variación en el proceso i, donde i = 1, 2, 3.

Por lo tanto, el ECM en el proceso industrial 1 es 1,103, es de, cir  ^{ECM = (5.011 - 5)2 + (-1.050 )2}

Paso (3d). Calculado el valor del ECM de    y   para  se procede a obtener dos medidas de la función a minimizar basadas en la perturbación simultánea a partir del valor actual , con las ^Ck y   de los pasos 1 y 2, utilizando las siguientes ecuaciones para obtener :

Los resultados obtenidos al sustituir los valores iniciales de las variables independientes del proceso 1 en las ecuaciones 3 y 4 son: ^{X1+ = 67.8°C}  y  ^{X2+ = 4 h 13 min}

Paso (3e). Después de calcular  se procede a sustituir los valores correspondientes a cada proceso industrial.

En el proceso industrial 1 se utilizan las ecuaciones ( ^- ) y (^- ); en el proceso industrial 2 se utilizan las ecuaciones ( ^- ) y ( ^- ); y en el proceso industrial 3 se utilizan las ecuaciones ( ^- ) y ( ^- ). Para el análisis del proceso 1 con los valores de   se sustituyen en la ecuación de regresión de segundo orden, dando como resultado ^{yµ1+ = 4.417} y para ^{+ = -0.902}.

Paso (3f). Obtenidos los valor es de ⁺, ⁺, ⁺  y  ⁺, ⁺, ⁺ para se procede a sustituir los valores en la ecuación 2 con la finalidad de calcular el error cuadrático medio (ECM) de . El ECM para el proceso uno sería de 1.153, es decir, .

Paso (3g). Se procede a obtener la otra medida de la función a minimizar basada en la perturbación simultánea a partir del valor actual , con las  ^Ck y   de los pasos 1 y 2, utilizando las siguientes ecuaciones para obtener  :

Los resultados obtenidos al sustituir los valores iniciales de las variables independientes del proceso 1 en las ecuaciones 5 y 6  ^{X1 =65°C} y ^{X2 = 40 h 30 min}, la sucesión de números reales ^{Ck = 0.932} y los vectores de perturbación simultánea ^{+ =3} y  ^{- = -0.3}; se obtienen los valores ^{X1- = 62.2°C}  y  ^{X2- = 4 h 47 min}

Paso (3h). Después de calcular   se procede a sustituir los valores correspondientes a los procesos industriales.

En el proceso industrial 1 se utilizan las ecuaciones (  ^-  ) y ( ), en el proceso industrial 2 se emplean las ecuaciones (  ^-  ) ( ); y en el proceso industrial 3 las ecuaciones ( ^- ) y (^- ). El análisis del proceso uno con los valores de   se sustituyen en la ecuación de regresión de segundo orden, dando como resultando  ^{- =5.616} y para  ^{- = -0.883}

Paso (3).i Obtenidos los valores de  ^{-,  -,   -,}  y  ^{-  -  -} para   se procede a sustituir los valores en la ecuación 2 con el objetivo de calcular el error cuadrático medio (ECM) de  .

El ECM para el proceso 1 sería de 1.160, es decir,  .

Paso 4. Aproximación del gradiente. Generar la aproximación de perturbación simultánea del gradiente desconocido como sigue:

Para obtener los resultados de las ecuaciones 7 y 8 en el proceso industrial 1 se sustituyen los valores del   y  , la ecuación de sucesión   ^{= 0.932} y la perturbación simultánea  ^{+ =3},y  ^{- =-0.3} por lo que la solución es  ^{+ = -0.001} y  ^{- =0.0012}.

Paso 5. Actualizando el valor de   estimado. Actualizar el valor de   a un nuevo valor   se hace utilizando las fórmulas estándar del algoritmo estocástico como sigue:

Aplicando los resultados obtenidos en los pasos anteriores, se sustituyen los valores en las ecuaciones 9 y 10 obteniendo los nuevos valores para   y para  , resultados que son para iniciar la siguiente iteración del proceso industrial 1 dado en  .

Paso 6. Iteración o terminación. Regresar al paso 2 con ^k+1 reemplazando ^k . Terminar el algoritmo si hay un pequeño cambio en el ECM en varias iteraciones sucesivas o el número máximo de iteraciones ha sido rechazado.

Octavo, con la aplicación del procedimiento experimental mencionado en este ejemplo se encontraron los valores óptimos en las variables independientes   y    de los tres procesos industriales evaluados, en función del valor mínimo obtenido del ECM. Las simulaciones del AEPS fueron realizadas mediante el software Matlab®, obteniéndose buenos resultados para su validación.

Noveno, obtenidos esos valores óptimos se realizaron corridas de verificación y validación a la réplica 34 en cada proceso industrial.

Resultados

Los resultados obtenidos durante el desarrollo de esta investigación se presentan con base en el procedimiento planteado anteriormente para los tres procesos industriales evaluados. A estos efectos, se obtuvieron los siguientes resultados, mostrados en las tablas 1, 2 y 3.

Los resultados muestran que en el proceso 1 se logra obtener un ECM mínimo en la tercera iteración, llamada  ₂, ya que en la iteración 4, denominada  ₃, hubo un aumento del ECM y por tal motivo se detiene el algoritmo, tal y como se menciona en el paso 6, lo cual nos indica que la mejor alternativa para hacer las pruebas experimentales de validación es la obtenida en la tercera iteración con los valores para la temperatura de  ^{= 66°C} y para el tiempo de  ^{= 3h 42 min}, con un ^{ECM = 0.009} como se muestra en la tabla 1.

En la tabla 2 se ofrecen los resultados del proceso 2, en el cual se logra obtener un ECM mínimo en la iteración inicial llamada   ya que en las iteraciones 2 y 3 hubo aumentos considerables del ECM y por tal motivo se detuvo el algoritmo, tal y como se menciona en el paso 6, lo cual nos indica que la mejor alternativa para hacer las pruebas experimentales de validación es la obtenida en la iteración inicial con los valores óptimos para la variable de temperatura de  ^{= 65°C} y para la variable de tiempo de  ^{= 6h 06 min}, con un ^{ECM = 0.025}. Los resultados señalan que para el proceso 3 se logra obtener un ECM mínimo en la iteración inicial, llamada  , ya que en las iteraciones 2 y 3 hubo aumentos considerables del ECM, al igual que en el proceso 2 y por tal motivo se detuvo el algoritmo, tal y como se menciona en el paso 6, lo cual nos indica que la mejor alternativa para hacer las pruebas experimentales de validación es la obtenida en la iteración inicial con los valores óptimos para la variable de temperatura de  ^{= 80°C} y para la variable detiempo de   ^{=5h 06 min}, como un ^{EMC = 0.004}, como se indica en la tabla 3.

Con base en los resultados obtenidos mediante la aplicación del AEPS a la simulación, se llevaron a cabo pruebas de verificación y validación a la réplica 34 en cada uno de los procesos evaluados. Los resultados promedios de estos experimentos de validación en las variables independientes de temperatura y tiempo, y las variables de respuesta 1 y 2, se exponen en la tabla 4.

Los resultados de la tabla 4 muestran que el proceso 1 trabajará con una temperatura de 66 °C y un tiempo de proceso de 3 h 42 min para obtener un valor promedio de humedad final de 5,3% en el producto terminado y con un valor de 80° en la escala de color. El proceso 2 trabajará con una temperatura de 65 °C y un tiempo de proceso de 6 h 6 min para obtener un valor promedio de humedad final de 4,9% en el producto terminado y con un valor de 79° en la escala de color. El proceso 3 trabajará con una temperatura de 80 °C y un tiempo de proceso de 5 h 06 min para obtener un valor promedio de humedad final de 3,5% en el producto terminado y con un valor de 78° en la escala de color.

Como se puede observar, este proceso no logra establecer condiciones óptimas favorables con respecto al valor meta de la variable de respuesta 1 (T = 5%), ya que nos resulta un valor de humedad final menor al 5% y por debajo del intervalo mencionado al inicio de este artículo (4 al 6%). En la variable de respuesta 2 (color) sí logra obtener un valor dentro del intervalo de 75 a 85° en la escala de color, como se muestra en la tabla 4.

Conclusiones y trabajos futuros

En este trabajo de investigación se propone un algoritmo de aproximaciones estocásticas modificado, el cual trabaja con modelos de segundo orden para determinar el valor óptimo de las variables que intervienen en tres procesos industriales. Los resultados obtenidos son: para el proceso 1 se operará con una temperatura de 66 °C con un tiempo de proceso de 3 h 42 min con la finalidad de obtener un valor promedio de humedad final de 5,3% y con un valor de 80^o en la escala de color en el producto terminado. Para el proceso 2, con una temperatura de 65 °C y un tiempo de proceso de 6 h 06 min para obtener un valor promedio de humedad final de 4,9% y con un valor de 79° en la escala de color en el producto terminado. Para el proceso 3, con una temperatura de 80 °C y un tiempo de proceso de 5 h 06 min para obtener un valor promedio de humedad final de 3,5% y con un valor de 78° en la escala de color en el producto terminado. Se concluye que al trabajar con un modelo de segundo orden dicho algoritmo nos proporciona resultados satisfactorios en los procesos 1 y 2 evaluados, en los cuales las variables de respuestas analizadas se encuentran dentro de los parámetros satisfactorios para el cumplimiento de calidad del producto en cuanto a humedad final y color, ya que el proceso 3 en su variable de respuesta 1 está por debajo del parámetro establecido de humedad final, pero sí es satisfactoria la variable de respuesta 2, "color". Además se concluye que es un algoritmo simple dado que no requiere un conocimiento profundo sobre el proceso, ni de la verdadera relación funcional entre la variable de respuesta y los factores controlables, además de ser fácil de usar ya que no requiere ser operado por personal altamente calificado.

Los trabajos futuros de investigación evaluarán el índice de capacidad de los procesos   para medir la capacidad o aptitud de los procesos, así como el AAEPS mediante los diseños compuestos centrales (DCC) para analizar si existe una mejor eficiencia con respecto a los diseños 3^k.

---

Referencias

Andradóttir, S., A stochastic Approximation Algorithm with Varying Bounds., Operations Research, Vol. 43, número 6, 1995i, pp 1037-1048.

Andradóttir, S., A method for Discrete Stochastic Approximation., Management Science, Vol. 41, número 12, 1995ii, pp 1946 -1961.

Andradóttir, S., A Scaled Stochastic Approximation Algorithm., Management Science, Vol. 42, número 4, 1996, pp 475-498.

Blum, J.R., Multidimensional Stochastic Approximation Methods., Annals of Mathematical Statistic, Vol. 25, 1954, pp 737-744.

Brooks, O., Solving Discrete Resource Allocation Problems using the Simultaneous Perturbation Stochastic Approximation (SPSA) Algorithm, Proceedings of the Spring Simulation Multi-conference, 25-29 March 2007, Norfolk, VA, USA, pp. 55- 62.

Chien, S.I., Luo, J., Optimization of Dynamic Ramp Metering Control with Simultaneous Perturbation Stochastic Approximation., Control and Intelligent Systems, scheduled for fall 2008 issue, pp 8-10.

Chin, D.C., Comparative Study of Stochastic Algorithms for System Optimization Based on Gradient Approximation., IEEE Transaction on Systems, Man, and Cybernetics-partb:Cybernetics, Vol. 27, número 2, 1997, pp 244-249.

Delyon, B., General Results on the Convergence of Stochastic Algoritms., IEEE Transaction on Automatic Control, Vol. 41, número 3, 1996, pp 1245-1255.

Fu, M.C., Ho, Y.C., Using perturbation analysis for gradient estimation, averaging and updating in a stochastic approximation algorithm., Winter Simulation Conference Proceedings of the 20th conference on Winter simulation, 1988, pp 509-517.

Kiefer, J., Wolfowitz, J., Stochastic Estimation of the Maximum of a Regression Function., Annals of Mathematical Statistic, Vol. 23, número 3, 1952, pp. 462-466.

Kulkarni, S.R., Horn, C.S., An Alternative Proof for Convergence of Stochastic Approximation Algorithms., IEEE Transactions on Automatic Control, Vol. 41, número 3, 1996, pp 419-424.

Kushner, H.J., Clark, D.J., Stochastic Approximation Methods for Constrained and Unconstrained Systems., New York, Springer-Verlag, 1978.

Maeda, Y., Time difference Simultaneous Perturbation Method., Electronic Letters, Vol. 32, número 11, 1996, pp 1016-1017.

Maryak, J.L., Chin, D.C., Global Random Optimization by Simultaneous Perturbation Stochastic Approximation., IEEE Transactions on Automatic Control, vol. 53, número 3, 2008, pp. 780-783.

Montgomery, D.C., Desing and Analysis of Experiments, Seventh ed., NJ, John Wiley & Sons, 2009. pp 360-368

Polyak, B.T., New Method of Stochastic Approximation Type Procedures., Automatica I telemekhanika, Vol. 51 (1990) pp 98-107 en Ruso, trasladado al Inglés en Automatica Remote Control, Vol. 51, 1991, pp 937-945.

Polyak, B.T., Juditsky, A.B., Acceleration of Stochastic Approximation by Avering., SIAM Journal on Control and Optimization, Vol. 30, número 4, 1992, pp 838-855.

Robbins, H., Monro, S., A Stochastic Approximation Method., Annals of Mathematical Statistic, Vol.22, 1951, pp 400-407.

Spall, J.C., Implementation of the Simultaneous Perturbation Algorithm Stochastic Optimization., IEEE Transactions on Aerospase end Electronic Systems, Vol. 34, Número 3, 1998, pp 817-823.

Spall, J.C., Introduction to Stochastic Search and Optimization Estimation., Simulation and Control, NJ, John Wiley & Sons, Hoboken, NJ. 2003.

Jesús Everardo Olguín Tiznado¹, Rafael García Martínez², Claudia Camargo Wilson³, Juan Andrés López Barreras⁴

¹ Industrial Engineering, Instituto Tecnológico de Huatabampo. Master of Science inIndustrial Engineering , Instituto Tecnológico de Hermosillo. Researcher and Professor, Universidad Autónoma de Baja California. Mexico. jeol79@uabc.edu.mx

² B.Sc. in Mathematics, Universidad Autónoma de Nuevo León. Maestro Master andDoctor of Science in Industrial Engineering, Instituto Tecnológico de Ciudad Juarez. Director of the Instituto Tecnológico del Valle del Yaquí. Mexico. ra_garcia@ith.mx

³ Industrial Engineering, Instituto Tecnológico de Los Mochis. Master of Science inIndustrial Engineering, Instituto Tecnológico de Hermosillo. Professor and researcher, Universidad Autónoma de Baja California. Mexico. ccamargo@uabc.edu,mx

⁴ Industrial Engineer, Instituto Tecnológico de Huatabampo. Master of Science inIndustrial Engineering , Instituto Tecnológico de Tijuana. Profesor Professor and researcher, Universidad Autónoma de Baja California. jalopez@uabc.edu.mx

---

ABSTRACT

Palabras clave: stochastic approximation algorithm, dependent variable, independent variable, iterative process, simulation

---

Received: October 22th 2010 Accepted: November 15th 2011

---

Introduction

The stochastic approximation method presented by Robbins and Monro (1951) was a linear search method from the root of the unknown function    representing a random variable's expected value. Kiefer and Wolfowitz (1952) modified it so that it could be used in determining optimum  . Blum (1954) extended the findings of previous authors to Cartesian spaces of dimension greater than 1.

From the work presented by Blum (1954), there was an increase in the number of stochastic approximation methods (Kushner and Clark, 1978; Polyak, 1991, Polyak and Juditsky, 1992; An-dradóttir, 1995 (i, ii); Delyon, 1996; Kulkarni and Horn, 1996; Maeda, 1996). Howver, Andradóttir (1996) has stated that all these methods are theoretical procedures lacking termination criterion used to determine ^X* in  , so that  , where   would be the function corresponding to the gradient vector of function  . Its analytical expression is unknown, but it is possible to quantify its value to a specific combination of values or levels of controllable factors; such measurement is subject to experimental error, which of course is not set in terms of its probability distribution.

Chin (1997) has classified stochastic approach procedures into two types: Robbins-Monro and Kiefer-Wolfowitz types. The former are characterised by requiring direct observations of ^h, including the Robbins-Monro ascending steps methods and the Newton-Raphson perturbation analysis and likelihood rate, the latter requiring estimates or approximations of ^h, such as Kiefer-Wolfowitz finite difference, the random directions method, the scaling method and the simultaneous perturbation stochastic approximation algorithm. The latter are considered more useful, since they do not require in-depth knowledge of the system to be optimised, i.e. they are applicable in situations where the functional relationship between the response variable denoted  ^yi and controllable factors ^d denoted by vector   are unknown.  represents the Cartesian space of dimension ^d, a situation that occurs most frequently in practice. Fu and Ho (1988) and Chin (1997) have pointed to the simultaneous perturbation stochastic algorithm as being the most efficient (as much as is theoretically practical) because it has a higher convergence rate and requires a smaller number of observations in each iteration. The latter is of great interest since it is directly proportional to the economic cost and simplicity of the experimental work.

This research paper proposes a modified stochastic approximation algorithm in which sequences   are validated

by simulation and real implementation, through a second-order model on obtaining the optimum operation values in the response variables involved in three similar industrial processes. This was due to poor operating condition control in such processes generating this waste product. It also demonstrates that the simultaneous perturbation stochastic approximation algorithm (SPSAA) is efficient to work with second order models, as shown in equation (1):

Experimental development

The materials used for developing and validating this research project were a Pentium PC, Protégé R200, M, 1.2 GHz processor and 598 MHz RAM and STATISTICA and MATLAB software for statistical analysis of the data.

The method used to collect the information required in experimental analysis for determining optimum operating conditions for three industrial processes analysed is described.

A list of significant independent variables or controllable factors was drawn up, including their ranges, denoted by the vector X=(X1,...Xd) ε R^d where R^d represented the Cartesian space for dimension d. In this research, the significant variables used were X1 representing temperature ranging from 60°C to 70°C and X2 representing the time of onset ranging from 4 to 5 hours within the three processes. Independent variables' starting point was X1= 65°C and X2 = 4 h 30 processing.

A list of the dependent variables (responses) and their units was then drawn up, denoted by vector (nx1) observations Y=(Y1,...,,Yd) ε R^d. The response variables supporting the research were listed, these being ^y1 representing response 1 of the processes, its unit being given as a percentage, ^y2 representing response 2, its unit given in degrees. The nominal values or intended goal obtained for the answers were 5% (target) of final moisture for ^y1 and 80° (target) on Hunter Laboratories' colour scale for ^y2 . This was so that the processes met the requirements requested by the customer ( ^{y1 =4% to 6%} final moisture in the product and ^{y2 =75° to 85°} colour scale).

An experimental type of design was selected. In this case the design was generated from a 3^k factorial, which meant k factors at three levels of experimental analysis (Montgomery, 2009). It will be two factor (temperature and time) at three levels (60°C, 65°C and 70°C for temperature, and 4 h, 4 h 30 min, and 5 h for time) for this work. Step one mentioned how the work was done for this experiment with these variables' initial values for each process.

The experiments were randomly assigned. The runs were made at random in each experimental stage.

Experiments and data collection involved five replicates in the experiment, regarding conditions for the independent variables' initial values for each process regarding target values for each response variable. This was done to obtain the mean (µ) and the corresponding standard deviation (σ).

Once the data had been collected, second-order models were obtained for response variables y₁ and y₂for mean (µ) and standard deviation (σ) for each industrial process. For example, the second-order regression equations (   ) and (   ) for industrial process one were:

The simultaneous perturbation stochastic approximation algorithm (SPSAA) was calculated according to the following steps, as in Spall J.C. (1998):

Step 1: initialisation and selection coefficient. Index counter k=1 was selected. An assumed value of the initial gradient vector ^θ0 and non-negativity coefficients ^c, ^A, and were taken. Delyon (1996), Spall (2003) and Chien and Luo (2008) have stated that the value typically assumed for when the gradient vector is equal to the arithmetic average of ^m estimates. Values were practically effective and theoretically valid for and , being ^0.602 and ^0.101 respectively (asymptotic optimal values of ^1.0 and ^1/6 could be used as well) and values for , ^c, and ^A, and could be determined as shown below. This was a useful guide for selecting ^A as it was much less than the maximum number of iterations allowed or expected, which is why ^{A=100 =0.16}, ^c=1 were selected. The results based on the data were   and  .

Step 2: generating a simultaneous perturbation vector. A random perturbation vector p-dimensional, was generated by the Monte Carlo simulation method where each ^p component of was generated independently from a mean zero probability distribution. A simple (and theoretically valid) choice for each component was to use a Bernoulli distribution ^±1 with ^1/2 probability or each ^±1Outcome. It should be noted that uniform variables and random normal were not allowed for items of   for regular SPSA conditions (Brooks O. 2007; Maryak and Chin, 2008). Therefore, the vector values for all three processes were:  ^{+ =3°C} and  ^{- = -0.3} hours; these were values assigned by the experimenter for the three industrial processes based on their characteristics.

Step 3: evaluating the functions to be minimised. values for   were selected for simulations based on the runs involving 60°C, 65°C and 70°C, and even for others which were not performed 50°C, 55°C, 75°C and 80°C with 4 h, 4 h 30 min and 5 h, plus 3 h, 3 h 30 min, 5 h 30 min, and 6 h, respectively. The independent variables' initial values for processes were  .

3.b. After   had been selected, values were replaced for various independent variables in regression equations  ^yµ and ^yσ previously obtained with experimental design data ^3k referred to in the third step. The second-order regression equations (    ) and (   ) represented the industrial process, as follows:

Second-order regression equations (   ) and (    ) represented the industrial process, as follows:

Second-order regression equations (   ) and (   ) represented the industrial process, as follows:

Substituting the initial values of   and , , in second-order regression equations (  ) and (   ) for the industrial process analysis one led to regression equation values    ^=5.011 and    ^=-1.050.

3.c. Once   ,   ,    and   ,   ,    values had been obtained for     the values in equation (2) were replaced to calculate the mean square error (MSE) for each industrial process, as follows:

where: ^yµi represented the response variable for the average in process i, where i= 1, 2, 3. T represented the target value for the process, representing the response variable for variation in process i where i = 1, 2, 3.

The MSE for industrial process one was thus 1.103, i.e.^{ECM = (5.011 - 5)2 + (-1.050 )2}

3.d. Once the MSE had been calculated for   and   for   then two measurements of the function to be minimised were obtained, based on simultaneous perturbation from current value  , with ^Ck and   from steps 1 and 2; using the following equations to obtain :

The results obtained by replacing the independent variables' initial values for a process in equations (3) and (4) were: ^{X1+ = 67.8°C} and ^{X2+ = 4 h 13 min}.

3.e. After calculating , then the corresponding values for each industrial process were replaced. Equations ( ^-) and (^- ) were used for industrial process one, equations ( ^- ) and ( ^- ) for industrial process two and ( ^- ) and ( ^-) for industrial process three. Values of   were substituted into second-order regression equations for process one analysis: ^{yµ1+ = 4.417} and ^{+ = -0.902}.

3.f. Once ⁺, ⁺, ⁺   and ⁺, ⁺, ⁺ had been calculate for   then values were replaced in equation (2) to calculate the MSE for . The MSE for industrial process one would be 1.153, i.e. .

3.g. The other measurement of the function to be minimised was obtained by simultaneous perturbation from current value ;, with ^Ck and from steps 1 and 2, using the following equations to obtain :

The results obtained by replacing the independent variables' initial values for process one in equations (5) and (6) were:  ^{X1 =65°C} and ^{X2 = 40 h 30 min} . The sequence of real numbers ^{Ck = 0.932} and simultaneous perturbation vectors ^{+ =3} and  ^{- = -0.3} led to obtaining ^{X1- = 62.2°C}  and   ^{X2- = 4 h 47 min}.

3.h. After calculating   m then the industrial processes' corresponding values were replaced. Equations (  ^-) and ( ^-) were used for industrial process one, (  ^-) and ( ^-) for industrial process three. Analysis of process one having values was substituted into the second-order regression equation: ^{- =5.616} and^{- = -0.883}.

3.i. Once values had been obtained for  ^{-,  -,   -,} and ^{-  -  -} for then values were replaced in equation 2 to calculate MSE.

MSE for process one would be 1.160, i.e.  .

Step 4. Gradient approximation. The simultaneous perturbation approximation of unknown gradient was as follows:

The values of   and   industrial process one were replaced to obtain results in equations (7) and (8) giving   ^{= 0.932} succession equation and  ^{+ =3} and ^{- =-0.3} simultaneous perturbation. The solution was thus  ^{+ = -0.001} and  ^{- =0.0012}.

Step 5. Updating the estimated value of   . Updating the value of   to a new value of   was done by using stochastic algorithm standard formulae, as follows:

Applying the results obtained in the previous steps, then values in equations (9) and (10) were replaced for obtaining new values:   and   . The results were used to start the next iteration of the industrial process given in  .

Step 6. Iteration or termination. Step 2 was repeated with ^k+1 replacing ^k. The algorithm was terminated if there was a small change in MSE in several successive iterations or the maximum number of iterations had been rejected.

When applying the experimental procedure mentioned in this example then optimal values were found for the three industrial processes' independent variables   and   based on the minimum value obtained from the MSE. SPSAA simulations were performed using MATLAB software with good validation results.

Once these optimal values had been obtained, verification and validation runs were performed for each manufacturing process.

Results

Tables 1, 2 and 3 show the results obtained for the three industrial processes, respectively.

The results show that a minimum MSE in the third iteration called  ₂ was obtained in process one since MSE increased in the fourth iteration named  ₃ and thus the algorithm was stopped, as mentioned in step 6. This indicated that the best alternative for testing experimental validation was obtained in the third iteration: w  ^{= 66°C} and   ^{= 3h 42 min} hours ^{ECM = 0.009} ( Table 1).

Table 2 shows the results for process two where a minimum MSE was obtained in the initial iteration called   since there were significant increases in MSE in iterations two and three and thus the algorithm was stopped. This indicated that the best alternative for testing experimental validation was obtained in initial iteration:  ^{= 65°C} and  ^{= 6h 06 min}, ^{ECM = 0.025}.

The results showed that a minimum MSE was obtained in the initial iteration called   for process three since there were significant increases in MSE in iterations two and three for process two, and thus the algorithm was stopped, indicating that the best alternative for testing experimental validation was obtained in the optimal values' initial iteration:  ^{= 80°C},   ^{=5h 06 min}, ^{EMC = 0.004} (as shown in Table 3).

Based on the results obtained by applying SPSAA, simulation, verification and validation tests were carried out on the 34 replicas in each process evaluated. Table 4 shows the average results for these experiments for validating the independent variables of temperature and time, and response variables 1 and 2.

Table 4 shows that process one would work at 66°C and 3 h 42 min processing time to obtain 5.3% average final moisture content in the finished product and 80° on the colour scale. Process two would work at 65°C and 6 h 06 min processing time to obtain 4.9% average final moisture content in the finished product and 79° on the colour scale. Process three would work at 80° C having a 5 h 06 min processing time to obtain 3.5% average final moisture content in the finished product and 78° on the colour scale.

It can be seen that this process failed to establish optimal favourable conditions regarding variable one target value response (T = 5%) and final moisture value was below 5% and below the range mentioned at the beginning of this article (4%-6%). A value was obtained within the 75° to 85° colour scale range in response variable two (colour) (Table 4).

Conclusions and future work

A modified stochastic approximation algorithm was proposed in this research project, working with second-order models to determine the optimal value of the variables involved in 3 industrial processes. Table 4 shows the results obtained. It was thus concluded that this algorithm provided satisfactory results regarding processes 1 and 2 when working with a second-order model in which the response variables analysed came within satisfactory parameters for achieving product quality in terms of specified final moisture content and colour. Variable response 1 for process 3 was below the established final moisture parameter while response variable 2 "colour" was certainly successful It was also concluded that this is a simple algorithm to apply, given that it did not require a deep understanding of a particular process or of the true functional relationship between the response variable and controllable factors and it was easy to use as it did not need not be operated by highly qualified personnel.

Future research will be aimed at evaluating   process capability index to measure the process' ability or aptitude. The simultaneous perturbation stochastic approximation algorithm (SPSAA) will be evaluated by using central composite design (CCD) to analyse whether better efficiency can be achieved with regarding 3k designs.

---

References

Andradóttir, S., A stochastic Approximation Algorithm with Varying Bounds., Operations Research, Vol. 43, número 6, 1995i, pp 1037-1048.

Andradóttir, S., A method for Discrete Stochastic Approximation., Management Science, Vol. 41, número 12, 1995ii, pp 1946 -1961.

Andradóttir, S., A Scaled Stochastic Approximation Algorithm., Management Science, Vol. 42, número 4, 1996, pp 475-498.

Blum, J.R., Multidimensional Stochastic Approximation Methods., Annals of Mathematical Statistic, Vol. 25, 1954, pp 737-744.

Brooks, O., Solving Discrete Resource Allocation Problems using the Simultaneous Perturbation Stochastic Approximation (SPSA) Algorithm, Proceedings of the Spring Simulation Multi-conference, 25-29 March 2007, Norfolk, VA, USA, pp. 55- 62.

Chien, S.I., Luo, J., Optimization of Dynamic Ramp Metering Control with Simultaneous Perturbation Stochastic Approximation., Control and Intelligent Systems, scheduled for fall 2008 issue, pp 8-10.

Chin, D.C., Comparative Study of Stochastic Algorithms for System Optimization Based on Gradient Approximation., IEEE Transaction on Systems, Man, and Cybernetics-partb:Cybernetics, Vol. 27, número 2, 1997, pp 244-249.

Delyon, B., General Results on the Convergence of Stochastic Algoritms., IEEE Transaction on Automatic Control, Vol. 41, número 3, 1996, pp 1245-1255.

Fu, M.C., Ho, Y.C., Using perturbation analysis for gradient estimation, averaging and updating in a stochastic approximation algorithm., Winter Simulation Conference Proceedings of the 20th conference on Winter simulation, 1988, pp 509-517.

Kiefer, J., Wolfowitz, J., Stochastic Estimation of the Maximum of a Regression Function., Annals of Mathematical Statistic, Vol. 23, número 3, 1952, pp. 462-466.

Kulkarni, S.R., Horn, C.S., An Alternative Proof for Convergence of Stochastic Approximation Algorithms., IEEE Transactions on Automatic Control, Vol. 41, número 3, 1996, pp 419-424.

Kushner, H.J., Clark, D.J., Stochastic Approximation Methods for Constrained and Unconstrained Systems., New York, Springer-Verlag, 1978.

Maeda, Y., Time difference Simultaneous Perturbation Method., Electronic Letters, Vol. 32, número 11, 1996, pp 1016-1017.

Maryak, J.L., Chin, D.C., Global Random Optimization by Simultaneous Perturbation Stochastic Approximation., IEEE Transactions on Automatic Control, vol. 53, número 3, 2008, pp. 780-783.

Montgomery, D.C., Desing and Analysis of Experiments, Seventh ed., NJ, John Wiley & Sons, 2009. pp 360-368

Polyak, B.T., New Method of Stochastic Approximation Type Procedures., Automatica I telemekhanika, Vol. 51 (1990) pp 98-107 en Ruso, trasladado al Inglés en Automatica Remote Control, Vol. 51, 1991, pp 937-945.

Polyak, B.T., Juditsky, A.B., Acceleration of Stochastic Approximation by Avering., SIAM Journal on Control and Optimization, Vol. 30, número 4, 1992, pp 838-855.

Robbins, H., Monro, S., A Stochastic Approximation Method., Annals of Mathematical Statistic, Vol.22, 1951, pp 400-407.

Spall, J.C., Implementation of the Simultaneous Perturbation Algorithm Stochastic Optimization., IEEE Transactions on Aerospase end Electronic Systems, Vol. 34, Número 3, 1998, pp 817-823.

Spall, J.C., Introduction to Stochastic Search and Optimization Estimation., Simulation and Control, NJ, John Wiley & Sons, Hoboken, NJ. 2003.