The performance of simulated annealing in parameter estimation for vapor-liquid equilibrium modeling

THERMODYNAMICS AND SEPARATION PROCESSES

The performance of simulated annealing in parameter estimation for vapor-liquid equilibrium modeling

A. Bonilla-PetricioletI, ^*; U. I. Bravo-Sánchez^I; F. Castillo-Borja^I; J. G. Zapiain-Salinas^I; J. J. Soto-Bernal^II

^IDepartamento de Ingeniería Química, Instituto Tecnológico de Aguascalientes, Phone and Fax: [phone] ext. 127, CP: 20256, Aguascalientes, Aguascalientes, México. Email: [e-mail]

^IIDepartamento de Sistemas y Computación Instituto Tecnológico de Aguascalientes, CP: 20256, Aguascalientes, Aguascalientes, México

ABSTRACT

In this paper we report the application and evaluation of the simulated annealing (SA) optimization method in parameter estimation for vapor-liquid equilibrium (VLE) modeling. We tested this optimization method using the classical least squares and error-in-variable approaches. The reliability and efficiency of the data-fitting procedure are also considered using different values for algorithm parameters of the SA method. Our results indicate that this method, when properly implemented, is a robust procedure for nonlinear parameter estimation in thermodynamic models. However, in difficult problems it still can converge to local optimums of the objective function.

Keywords: Vapor-liquid equilibrium; Simulated annealing; Nonlinear parameter estimation; Error-in-variable method; Global optimization.

INTRODUCTION

In many important areas of chemical, biochemical and petroleum engineering, mathematical models form the basis for the design, optimization and control of process systems (Esposito and Floudas, 1998; Englezos and Kalogerakis, 2001). Parameter estimation is a common problem in many engineering applications and is one of the steps involved in the formulation and validation of a mathematical model that describes a process of interest. Specifically, it refers to the process of obtaining values of the parameters by matching the model-based calculated values with the set of measurements (Englezos and Kalogerakis, 2001). If we use linear model equations, the problem is called linear estimation, while nonlinear estimation refers to the more general problem where the model equations are nonlinear functions of the parameters. As indicated by Esposito and Floudas (1998), the use of nonlinear models introduces an added level of complexity into the numerical estimation of model parameters.

Parameter estimation is generally based on the classical least squares or maximum likelihood approaches (Esposito and Floudas, 1998; Gau et al., 2000; Gau and Stadtherr, 2002). In the classical least squares approach, it is assumed that there is a set of independent variables not subject to measurement error, while errors in all measured variables are accounted for in the maximum likelihood approach. Both methods involve the minimization of an objective function subject to constraints representing the model equations. Depending on the system characteristics, we can solve a constrained or unconstrained optimization problem, which in general is nonlinear and potentially nonconvex.

Several local and global optimization methods have been used for parameter estimation in both linear and nonlinear models (Bard, 1974; Britt and Luecke, 1973; Fabries and Renon, 1975; Gmehling et al., 1977-1990; Anderson et al., 1978; Valko and Vajda, 1987; Tjoa and Biegler, 1992; Vamos and Hass, 1994; Esposito and Floudas, 1998; Gau et al., 2000; Gau and Stadtherr, 2002). Local optimization methods are not reliable for finding the global minimum of the problem; however, these methods are the most widely used. In the other hand, some applications of global optimization methods in parameter estimation have been reported by Esposito and Floudas (1998), Kleiber and Axmann (1998), Park and Froment (1998), Gau et al. (2000), Costa et al. (2000), Gau and Stadtherr (2002) and Dominguez et al. (2002).

An attractive alternative for the reliable solution of nonlinear parameter estimation is the use of stochastic optimization methods such as simulated annealing (SA) or genetic algorithm (GA). Even though these methods provide no formal guarantee for global optimization, they are reliable strategies and offer a reasonable computational effort in the optimization of multivariable functions (Michalewicz and Fogel, 1999). In fact, some papers have reported the application of SA and GA in parameter estimation using chemical engineering models (Kleiber and Axmann, 1998; Park and Froment, 1998; Costa et al., 2000).

In the context of thermodynamics, parameter estimation may also encounter computational difficulties due to the possibility of several local optimums in the objective function used as optimization criterion. As indicated by Dominguez et al. (2002), failing to identify the global optimum in parameter estimation may cause errors and uncertainties in equipment design and erroneous conclusions about model performance. In view of this, it is important to test different optimization methods to identify a reliable and efficient technique for this purpose.

In this paper, we report the application and evaluation of the SA method in parameter estimation for vapor-liquid equilibrium (VLE) modeling. We have tested the numerical performance of SA using the classical least squares and maximum likelihood approaches. Reliability and efficiency in the data-fitting procedure are also considered using different values for the algorithm parameters of the SA method.

METHODOLOGY

Problem Formulation

In accordance with Gau et al. (2000), consider that a set of observations y_ij of i= l, ...,m dependent response variables from j=1,..., ndat experiments is available, where the responses will be adjusted to an explicit model with independent variables and parameters . Measurement errors in can either be treated or neglected, and depending on the choice, we can have a least squares or maximum likelihood formulation. Different objective functions can be used to obtain the parameter values that provide the best fit. For the case of the classical least squares (LS) criterion, we use the following objective function

This function is optimized with respect to the model parameters and Gmehling et al. (1977-1990), Gau et al. (2000) and Dominguez et al. (2002) have used this formulation in data fitting for VLE modeling.

If we assume that there are measurement errors in the state variables z_ij for the experiments in the system to be modeled, the error-in-variable (EIV) formulation has the form

subject to

where is a vector of np model functions, nest is the number of state variables, represents the unknown "true" values of state variables for each measurement and s_i represents the standard deviation associated with the measurement of state variable i. The optimization variables of this problem are the set of and the model parameters . In this case, there is a substantial increase in the dimensionality of the optimization problem, which depends on the number of experiments. Esposito and Floudas (1998) and Gau and Stadtherr (2002) have applied their deterministic optimization approaches to solving Eqs. [2] and [3] with different chemical engineering models including vapor-liquid equilibrium equations.

The last functions can be minimized either as a constrained or as an unconstrained optimization problem. In all our examples, we will consider only the unconstrained formulation. We have studied the numerical performance of the SA method in the global minimization of these functions using vapor-liquid equilibrium models. In the next section, we describe the algorithm used for the SA method.

Simulated Annealing (SA)

Simulated annealing is a stochastic optimization technique inspired by the thermodynamic process of cooling of molten metals to attain the lowest free energy state (Kirkpatrick et al., 1983). Starting with an initial solution and armed with adequate perturbation and evaluation functions, the algorithm performs a stochastic partial search of the state space. In minimization problems, uphill moves are occasionally accepted with a probability controlled by a parameter called annealing temperature, T_SA. The probability of acceptance of uphill moves decreases as T_SA decreases. At high temperatures, the search is almost random, while at low temperatures the search becomes almost greedy. At zero temperature, the search becomes totally greedy, i.e., only good moves are accepted (Kirkpatrick et al., 1983). The core of the algorithm is the Metropolis procedure, which simulates the annealing process at a given T_SA (Metropolis et al., 1953). The Metropolis criterion is used to accept or reject the uphill moves. Several algorithms have been proposed for the SA method. We have used the algorithm proposed by Corana et al. (1987) because previous papers have reported its reliability and efficiency in thermodynamic calculations (Henderson et al., 2001; Rangaiah, 2001).

In the algorithm proposed by Corana et al. (1987), a trial point is randomly chosen within the step length VM (a vector of length n variables) of a starting point. The function is evaluated at this trial point and its value is compared to its value at the initial point, where the Metropolis criterion is used to accept or reject the trial point. If the trial point is accepted, the algorithm moves on from that point. If it is rejected, another point is chosen instead for a trial evaluation. Each element of VM is periodically adjusted so that half of all function evaluations in that direction are accepted. A fall in T_SA, after NT iterations, is imposed upon the system with the RT variable by

where j is the iteration counter and RT is the temperature reduction factor. For our examples, we have defined this parameter equal to 0.85 which is the value suggested by Goffe et al. (1994). Thus, as T_SA declines, downhill moves are less likely to be accepted and the percentage of rejections rises. Given the scheme for the selection of VM, VM falls. Thus, as T_SA declines, VM falls and SA focuses upon the most promising area for optimization. A full description of this algorithm is found in Corana et al. (1987) and we have used the subroutine developed by Goffe et al. (1994). Figure 1 shows a flow diagram of this optimization method. The choice of cooling schedule is a crucial aspect in the implementation of SA because it affects the numerical performance of the optimization procedure. Considering this fact, we have tested the effect of parameters T_SA and NT on the numerical performance of SA to find the global optimum of the objective function used in the data-fitting procedure. Several calculations are done to perform this evaluation and a statistical analysis is used to establish their relative effect on the reliability of the SA method.

RESULTS AND DISCUSSION

We test the SA method in the data-fitting procedure using several experimental VLE data and different thermodynamic models. All systems used in this paper have an objective function with multiple optimums and they have been applied for testing global deterministic optimization methods (Esposito and Floudas, 1998; Gau et al., 2000; Gau and Stadtherr, 2002). The reliability and efficiency of SA is evaluated based on the following criteria: a) success rate in finding the global minimum SR and b) average number of objective function evaluations during the optimization procedure NFEV. All examples are solved 100 times (each time with a different random initial value and random number seed) and the results obtained are divided into four groups (each group with 25 calculations). For each group of results, SR and NFEV are calculated and a variance analysis is performed to establish the effect of parameters T_SA and NT on the reliability of the data-fitting procedure. We consider the statistical effect of these parameters significant if the value of the p-level obtained from the statistical analysis is lower than 0.05. In all our examples, we have used a tolerance of 1.0E-09 for the convergence of the SA algorithm. The experimental data are taken from DECHEMA and Reid et al. (1987). In all examples, we have assumed gas ideal behavior and the vapor pressures of pure components are calculated using the Antoine equation. For problems 1 – 4, we use the constants reported in DECHEMA for the Antoine equation.

Problem 1. Tert Butanol-1 Butanol

Our first example is the VLE for the binary system tert butanol-1 butanol. This system was studied by Gau et al. (2000) using an interval analysis approach and classical least square formulation. We use three sets of experimental data under different isobaric conditions. The Wilson equation is used to calculate the liquid-phase activity coefficients, which are defined by

The binary parameters L₁₂ and L₂₁are given by

where v₁ and v₂ are the pure component liquid molar volumes, T is the system temperature and q₁ and q₂ are the energy parameters. Assuming an ideal vapor phase, the experimental values for the activity coefficients are calculated using the next equation

where is the vapor pressure of pure component i at the system temperature T and c is the number of components. The vapor pressure is calculated using

where T is in °C and, in mmHg. In accordance with Gmehling (1977-1990) and Gau et al. (2000), we use a relative least squares formulation to fit the data

We optimize this function with respect to the Wilson model parameters inside the intervals 0₁Î (–8500, 320000) and 0₂Î (–8500, 320000). The initial values for each calculation are randomly generated within these intervals. The performance of SA is tested using the next arbitrary values for its parameters: T_SA (10, 100, 1000) and NT(10, 20). The results obtained by other thermodynamic calculations suggest that our proposed values for T_SA and NT are suitable choices and favor the numerical performance of SA. We don't use higher values of NT because the computational time of the SA method is significantly increased. The results obtained for all combinations of T_SA – NT are reported in Table 1 and the statistical effect of these parameters on the reliability of the data-fitting procedure is reported in Table 2. For all pressure conditions, SA fails several times to find the global minimum of the objective function. In all data sets, the parameter T_SA shows a significant statistical effect on the reliability of the parameter estimation procedure (p-level £ 0.05). In general, an increase in T_SA improves the reliability of the SA method, of course at the expense of a major computational effort. Only for the experimental data at 700 mmHg, the parameter NT and the interaction T_SA – NT show a significant statistical effect on the reliability of the data fitting. The global optimums reported in Table 1 are consistent with the results of Gau et al. (2000).

Problem 2. Water-1,2 Ethanediol

We use one data set for the system water-1,2 ethanediol at 430 mmHg. This system was also analyzed by Gau et al. (2000) using the least squares formulation and they reported several local optimums for the objective function. Also, they indicated that global optimal parameters provide a better fit to the experimental data. We optimize Eq. 11 with respect to the energy parameters of the Wilson, NRTL and UNIQUAC models. The energy parameters for the NRTL and UNIQUAC equations are defined as in DECHEMA. We use the same initial intervals employed in the last example for the Wilson model, while the feasible domain for the NRTL and UNIQUAC parameters are 0₁,0₂Î (-2000, 5000) and a12 Î (0.01, 10) and 0₁,0₂Î (-5000, 20000), respectively. These domains include the parameters reported in DECHEMA. Our results are presented in Table 3 and the statistical analysis is reported in Table 2. For the Wilson and UNIQUAC models, the parameter T_SA and the interaction T_SA – NT show a significant statistical effect on the reliability of the data-fitting procedure. Also, SA shows some failures in the global minimization of the objective function. The global minimum reported for the Wilson model is in agreement with that published by Gau et al. (2000). The energy parameters reported in DECHEMA correspond to local optimums for the three models. Similar findings with other systems have been reported by Dominguez et al. (2002) using the Wilson and UNIQUAC equations.

Problem 3. Benzene-Hexafluorobenzene

(LS Formulation)

The parameter estimation of this problem is performed assuming VLE data under isothermal and isobaric conditions. This system has also been analyzed by Gau et al. (2000) using interval analysis. Four data sets are used to test the numerical performance of the SA method. We use a least squares formulation in the parameter estimation and the same initial intervals as those for the energy parameters of the Wilson model are applied. The results of parameter estimation are reported in Table 4, while the statistical analysis appears in Table 2. In this example, SA has a better numerical performance; however, it still shows some failures in the global minimum of the objective function. Only for the data at 760 mmHg, SA finds the global minimum with a reliability of 100% and its numerical performance is independent of the values for parameters T_SA and NT. We find a significant statistical effect of parameter NT in the reliability of the parameter estimation procedure for the data at 30 °C, 50 °C and 300 mmHg. For the case of T_SA, it affects the reliability of SA for data at 50 °C. As indicated by Gau et al. (2000), a failure in the global minimization in parameter estimation affects the capability of the thermodynamic model in prediction of homogeneous azeotropes in this system. That failure in the model's performance in the calculation of azeotropes could result in design problems in separation processes. The global optimums found with SA are in agreement with those reported by Gau et al. (2000).

Problem 4. Benzene-Hexafluorobenzene (EIV Formulation)

In this problem, we consider the modeling of vapor-liquid equilibrium of the binary system benzene-hexafluorobenzene using the error-in-variable approach. This system was studied by Gau and Stadtherr (2002) using interval analysis. Four data sets are used and parameter estimation is considered using the Wilson equation for liquid-phase activity coefficients. The state variables are (x₁,y₁,P,T) where P is the system pressure in mmHg, T is the system temperature in °C and x₁ and y₁ are the liquid and vapor mole fractions of benzene. In accordance with Gau and Stadtherr (2002), a standard deviations of 0.001, 0.01, 0.75 and 0.1 are assumed for these variables. The vapor-liquid equilibrium can be described by the following equations

We formulate the data-fitting problem as unconstrained optimization using Eqs. [12] and [13] to eliminate P and y₁ from the objective function. For the unconstrained problem, the independent state variables are the set of for all measurements, while the optimization variables are and the set of . The objective function is defined as

This function is optimized with respect to 2+2ndat variables where ndat is the number of experimental data used in parameter estimation. The initial intervals of the independent state variables are set using plus and minus three standard deviations, while the intervals for the Wilson model parameters are defined as 0₁Î (–8500, 320000) and 0₂Î (–8500,320000) . The results of parameter estimation are reported in Table 5 and the statistical analysis appears in Table 2. Due to the increase in problem dimensionality, SA fails several times in finding the global minimum for all data sets. For only two cases, it shows 100% reliability in the global minimization of the objective function. Also, there is a significant increase in the computational effort of the data-fitting procedure, also caused by problem dimensionality. Of course, we can reduce the numerical effort of the SA method by decreasing the tolerance value and using a local search method to refine the solution obtained by SA. On the other hand, the variance analysis indicates that parameter T_SA affects the reliability of parameter estimation for all data sets. For data at 760 mmHg, NT has a significant effect on the success rate for finding the global optimum. An increase in both T_SA and NT improves the reliability of SA. The global optimums reported for all data are consistent with the results of Gau and Stadtherr (2002).

Problem 5. Methanol-1,2 Dichloroethane

This system was studied by Esposito and Floudas (1998) using a deterministic global optimization method. We consider the EIV parameter estimation using the van Laar model for liquid phase. The experimental data consists of five points for the state variables: pressure in mmHg, temperature in Kelvin and liquid and vapor mole fractions of methanol. The vapor pressure of pure components is calculated using

where is given in mmHg, T is given in Kelvin and the constants of the Antoine equation are taken from Esposito and Floudas (1998). The activity coefficients are defined by

We use an unconstrained optimization formulation in the data-fitting procedure. In accordance with Esposito and Floudas (1998) and Gau and Stadtherr (2002), the temperature is scaled by a reference temperature T_r = 323.15 K and the van Laar model parameters are defined considering the reference temperature as . The state variables of the system are and their standard deviations are defined as 0.005, 0.015, 0.75 and 0.000309. The objective function is given by

where TS= T/T^r. Equations [12] and [13] are used to eliminate P and y₁ in the objective function and we have twelve optimization variables where the initial intervals for the independent state variables are set using plus and minus three standard deviations. For the case of the van Laar model parameters, we use and . The results of parameter estimation are reported in Table 6. For all calculations performed, SA finds the global minimum of the objective function with 100% reliability. In this case, the parameters T_SA and NT affect only the efficiency of the data-fitting procedure. The global minimum reported is consistent with the results of Gau and Stadtherr (2002).

We have studied the numerical performance of SA with other examples reported by Gau et al. (2000) and Gau and Stadtherr (2002). The results obtained in those calculations, not reported here, indicate that SA is generally a reliable method for parameter estimation in VLE modeling.

Based on our experience with other SA algorithms, such as very fast simulated annealing (Sharma and Kaikkonen, 1999) and direct search simulated annealing (Ali et al., 2002), we can expect that the algorithm of Corana et al. (1987) has a better numerical performance and we consider it a suitable choice for data fitting in thermodynamic models.

CONCLUSIONS

The numerical performance of the simulated annealing method has been tested in parameter estimation for VLE modeling. The reliability and efficiency of the data-fitting procedure was evaluated with respect to the SA parameters T_SA and NT, using classical least squares and error-in-variable formulations. Our results indicate that parameter T_SA significantly affects the reliability of the data-fitting procedure. On the other hand, both parameters affect the efficiency of the SA method. The SA method has a better numerical performance using the classical least squares instead of the error-in-variables formulation. In general, SA is a robust optimization procedure, when properly implemented for parameter estimation using both formulations. However, it can fail in the global minimization of difficult systems.

ACKNOWLEDGMENTS

The authors acknowledge the financial support received from CONACyT and Instituto Tecnológico de Aguascalientes.

NOMENCLATURE

A, B parameters of the van Laar model a₁,b₁,c₁ parameters of the Antoine equation NT iteration number before temperature reduction P Pressure vapor pressure of pure component T temperature T_SA annealing temperature T_r reference temperature x₁, y₁ liquid and vapor mole fraction y_ij dependent response variable unknown "true" value of state variable s_i standard deviation of state variable q₁, q₂ energy parameters g activity coefficient

(Received: May 30, 2005 ; Accepted: December 4, 2006)

• Ali, M.M., Törn, A., and Viitanen, S., A direct search variant of the simulated annealing algorithm for optimization involving continuous variables. Computers & Operations Research, 29 (1), 87-102 (2002).
• Anderson, T.F., Abrams, D.S., and Grens, E.A. II., Evaluation of parameters for nonlinear thermodynamics models. AIChE Journal, 24 (1), 20-29 (1978).
• Bard, Y., Nonlinear parameter estimation, Academic Press, New York (1974).
• Britt, H.I., and Luecke, R.H., The estimation of parameters in nonlinear implicit models. Technometrics, 15, 233-247 (1973).
• Corana, A., Marchesi, M., Martini, C., and Ridella, S., Minimizing multimodal functions of continuous variables with the simulated annealing algorithm. ACM Transactions on Mathematical Software, 13 (3), 262-280 (1987).
• Costa, A.L.H., da Silva, F.P.T., and Pessoa, F.L.P., Parameter estimation of thermodynamic models for high-pressure systems employing a stochastic method of global optimization. Brazilian Journal of Chemical Engineering, 17 (3), 349-353 (2000).
• Dominguez, A., Tojo, J., and Castier, M., Automatic implementation of thermodynamic models for reliable parameter estimation using computer algebra. Computers and Chemical Engineering, 26 (10), 1473-1479 (2002).
• Englezos, P., and Kalogerakis, N., Applied parameter estimation for chemical engineers, Marcel-Dekker, New York (2001).
• Esposito, W.R., and Floudas, C.A., Global optimization in parameter estimation of nonlinear algebraic models via the error-in-variables approach. Ing. Eng. Chem. Res., 37 (5), 1841-1858 (1998).
• Fabries, J., and Renon, H., Method of evaluation and reduction of vapor-liquid equilibrium data of binary mixtures. AIChE Journal, 21 (4), 735-743 (1975).
• Gau C.Y., Brennecke J.F., Stadtherr M.A., Reliable nonlinear parameter estimation in VLE modeling. Fluid Phase Equilibria, 168 (1), 1-18 (2000).
• Gau, C.Y., and Stadtherr, M.A., Deterministic global optimization for error-in-variables parameter estimation. AIChE Journal, 48 (6), 1192-1197 (2002).
• Gmehling, J., Onken, U., and Arlt, W., Vapor-liquid Equilibrium Data Collection, Chemistry Data Series, Vol. I, Parts 1-8, DECHEMA (1977-1990).
• Goffe, B., Ferrier, G., and Rogers, J., Global optimization of statistical functions with simulated annealing. Journal of Econometrics, 60 (1-2), 65-99 (1994).
• Henderson, N., de Oliveria, J.R., Amaral Souto, H.P., and Pitanga Marques, R., Modeling and analysis of the isothermal flash problem and its calculation with the simulated annealing algorithm. Ind. Eng. Chem. Res., 40 (25), 6028-6038 (2001).
• Kirkpatrick, S., Gelatt, C., and Vecchi, M.P., Optimization by simulated annealing. Science, 220 (4598), 671-680 (1983).
• Kleiber, M., and Axmann, J.K., Evolutionary algorithms for the optimization of modified UNIFAC parameters. Computers and Chemical Engineering, 23 (1), 63-82 (1998).
• Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller A., and Teller, E., Equation of state calculations by fast computing machines. J. Chem. Phys., 21 (6), 1087-1092 (1953).
• Michalewicz Z., Fogel D.B., How to solve it? Modern heuristics, Springer (1999).
• Park, T.Y., and Froment, G.F., A hybrid genetic algorithm for the estimation of parameters in detailed kinetic models. Computers and Chemical Engineering, 22 (S1), S103-S110 (1998).
• Rangaiah, G.P., Evaluation of genetic algorithms and simulated annealing for phase equilibrium and stability problems. Fluid Phase Equilibria, 187-188, 83-109 (2001).
• Reid, R.C., Prausnitz, J.M., and Sherwood, T.K., The properties of gases and liquids, 3rd Ed., McGraw-Hill, New York (1987).
• Sharma, S.P., and Kaikkonen, P., Global optimisation of time domain electromagnetic data using very fast simulated annealing. Pure and Applied Geophysics, 155 (1), 149-168 (1999).
• Tjoa, T.B., and Biegler, L.T., Reduced successive quadratic programming strategy for errors-in-variables estimation. Computers Chemical Engineering, 16 (6), 523-533 (1992).
• Valko, P., and Vajda, S., An extended marquardt-type procedure for fitting error-in-variables methods. Computers Chemical Engineering, 11 (1), 37-43 (1987).
• Vamos, R.J., and Hass C.N., Reduction of ion-echange equilibria data using an error in variable approach. AIChE J., 40 (3), 556 569 (1994).

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  31 May 2007
• Date of issue
  Mar 2007

History

• Received
  30 May 2005
• Accepted
  04 Dec 2006