Stochastic approximation algorithm for industrial process optimisation

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

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.

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.