Approach to a reliable solution strategy for performing phase equilibrium calculations using MINLP optimization

I.J. Jerez, F. Muñoz and J.M. Gomez

Dto. Ing. Química, Universidad de los Andes, Bogotá, Colombia. E-mail: jorgomez@uniandes.edu.co

Abstract — The objective of this contribution is to propose a reliable strategy to solver the problem of phase equilibrium calculations for non-ideal systems, using the Gibbs free energy minimization. This type of problem, using the Gibbs free energy minimization, is usually formulated as a Mixed Integer Non-Linear Programming (MINLP) Optimization. This optimization problem allows the compositions to be associated with continuous variables, and the presence of phases in the equilibrium to be associated with the integer variables. The solution strategy proposes a bi-level approach. The first level combines a stochastic (Simulated Annealing - SA) and a local deterministic algorithm (Sequential Quadratic Programming - SQP), and solves a Non Linear Programming Problem (NLP). The continuous variables are considered at this level. The second level considers the integer variables. The advantage of this bi-level strategy lies in its easy implementation and in its proven efficiency to locate global optima with acceptable computational load. This article includes the study of the Water-Ethanol-Cyclohexane and Water-Ethanol-Glycerin systems. A comparative analysis was conducted using experimental data reported in published works and theoretical calculations by means of the Gamma-Phi classic method.

Keywords — Phase Equilibrium; Gibbs Free Energy; Non-Linear Optimization.

I. INTRODUCTION

In the chemical industry, most of the separation processes are based on phase equilibrium (Kojima et al., 1969; Seader and Henley, 2006), which streamlines the study and development of reliable theoretical proposals in order to estimate equilibrium compositions (Parodi and Campanella, 2009). From a thermodynamic point of view, with a given temperature and pressure, a multiple phase and component system is in equilibrium if it has the global minimum of the total Gibbs free energy.

Two different solution strategies have been proposed for this problem: the first one using the solution of simultaneous algebraic non-linear equations (Nichita et al., 2002) and the second one using numerical optimization. This last approach sets out an optimization problem where an objective function is formulated based on the total Gibbs free energy. The objective function is constrained by a subset of equations that are established by calculating the total Gibbs energy for each component in the different phases and are supported in the mass balance. The set of objective function and equality constrains constitutes the formulation of the optimization problem. The presence or absence of phases is determined by the value of the integer variables. The continuous variables correspond to the composition of the phases in the equilibrium system. The optimization problem is composed of a function and a group of non-linear equations. The optimization problem must be solved with suitable techniques.

It is important to emphasize that in both approaches, the necessary condition to find the solution is a minimization of the total Gibbs free energy and that the sufficient condition depends on finding the global minimum. Proving that the solution corresponds to the global minimum is a complex problem because it has to fulfill the second order optimality conditions and these conditions are not easy to be mathematically verified (Revollar et al., 2010).

Formulating and solving the phase equilibrium problems based on approximation by optimization has already been studied. The proposed strategies have been based on the use of local and global optimization algorithms and in some cases, like in this paper, propose a mix of two kinds of algorithms (Saber and Shaw, 2008; McDonald and Floudas, 1996; Inoue and Zhu, 2001). There are three papers that are particularly important to mention because they represent the different strategies used: Zhu and Xu (1999), propose a strategy with a stochastic algorithm, Simulated Annealing, combined with a Newton deterministic algorithm. The stochastic algorithm is in charge of finding local points and the deterministic algorithm calculates composition in points established beforehand. Afterwards, the minimums are compared in order to retain the lowest. Rangaiah (2001), solved the optimization problem separately with two stochastic algorithms, Simulated Annealing and Genetic Algorithms. The results of both algorithms were compared in their efficiency in order to find the global solution. Finally, Reneaume et al. (1996) used a variable binary integer vector in order to formulate a Non-Lineal Subproblem (NLP-subproblem). The fixed point homotopy method is used in order to locate different stationary points of the NLP - subproblems.

This article proposes a bilevel strategy that breaks up an MINLP optimization into a series of solutions of NLP optimization problems. The second of its levels includes a variable integer binary vector. Taking into account that the maximum number of phases of the equilibrium system is not large, it conducts a systematic search in all the possible different combinations. At the same time, in the first level, the NLP problem uses deterministic and stochastic algorithm coupling in order to avoid local solutions when they appear. The biggest problem in determining phase equilibrium is with establishing the global minimum, so the strategy must emphasize fulfilling this condition. If it ensures finding global solutions in the NLP problem, the MINLP problem will be resolved perfectly. While trying to find the global optimums and without being able to ensure this condition mathematically, stochastic and deterministic algorithms are coupled at the NLP problem level. The simultaneous use of these algorithms happens when the deterministic algorithm is insufficient.

It is important to specify that each type of algorithm has its strengths and weaknesses in the context of this thermodynamic problem. Stochastic algorithms are known in the search for global solutions, but if the search is done in systems with multiple phases, high computational times are needed. Meanwhile, deterministic algorithms based on derivatives present good search directions and take less computational time. However, deterministic algorithms frequently encounter local optima and show convergence problems because of their mathematical foundation. In this paper the Sequential Quadratic Programming SQP was chosen as a deterministic algorithm because of its capacity to solve non-linear optimization problems with very low solution time. The stochastic algorithm, the Simulated Annealing (SA) was implemented taking into consideration the physical coincidence between the theoretical base where this algorithm was developed and the Gibbs free energy minimization problem in multiphase systems. This algorithm is frequently used in situations where the search area is highly non-convex and requires leaving local solutions in order to look for global solutions.

This article is divided into three parts: the part included in numerals two and three describes the formulation of the problem and the solution strategy. Numerals four and five present: results obtained and their analysis in two case studies: Water-Ethanol-Cyclohexane and Water-Ethanol-Glycerin. Finally numeral six presents general conclusions.

II. OPTIMIZATION PROBLEM CONSTRUCTION

The formulation chosen for developing in this paper is the minimization of Gibbs free energy approximation at a given pressure and temperature. Equations representing Gibbs energy calculations, the mass balance and the objective function must be flexible so they can incorporate new continuous variables and binary integers. In that way, the mathematical construction will be extensible to any number of phases and components. Initially the problem is formulated as a Mixed Integer Non-Linear problem, where the presence or absence of a phase will be determined by a binary variable. The formulations corresponding to the objective function and the constraints are shown.

A. Objective function

(1)

The objective function represents the sum of Gibbs free energies of the i-components in the j-phases that constitute the system. In other words, the objective function consists of the sum of Gibbs free energy of the compounds present in each equilibrium phase. These equilibrium phases are determined by the value of the binary variable Y_j.

B. Constraints

Constraints are distributed in two different groups: the linear type and the non-linear type. The first type belongs to the mass balance that must be reached. It should be noted that some non linear equations change to be linear, because the solution strategy conducts a systematic search in all the possible different combinations of the number of phases, fixing the binary variables (Y_⊥j) as parameters at these expressions. The second type belongs to the constraints that must be used to express the Gibbs free energy in terms of the different compositions in the phases present in the equilibrium.

Mass balance constraints

(2)

These constraints correspond to the mass balance posed for each component (Eq. 2) and describe the sum of moles of the i-compound in each phase where it belongs and must be equal to the total number of moles of the i-compound specified before conducting the system solution. For example, if it is necessary to perform mass balance for a compound A present in the equilibrium of a vapor phase and three liquid phases, the expression obtained will be:

(3)

Gibbs free energy constraints

(4)

These Non-linear constraints are formed by the expressions to calculate the Gibbs free energy of an i-compound in a j-phase. In order to describe the chemical potential in terms of measurable parameters, Eq. 4 must be replaced. Equations 5 and 6 show functions replacing the term .

For a vapor phase:

(5)

Meanwhile for a liquid phase:

(6)

Due to the algebraic form of Eqs. 5 and 6, it is necessary to use of following bound constraints: . For all the cases: and . The main drawback, with this type of algebraic form of equations, is the impossibility to consider only one NLP problem with all the phases and obtain a global solution with some of them (phases) with component number of moles at zero.

It is worth mentioning that these expressions are valid taking into consideration that the vapor phase behaves in an ideal mixture, while the liquid phase can be deviated from ideal behavior. Finally, in Eq. 4 the following expressions are obtained by replacing Eqs. 5 and 6 and depending on the phase being described.

The constraint for the vapor phase will be:

(7)

The constraint for the liquid phase will be:

(8)

One should note that the number of constrains and variable is a function of the number of phases (F) and components (C). For this kind of MINLP problem, it has: F*C+C equality constraints, 2*F*F continuous variables and F binary variables.

III. SOLUTION STRATEGY

The solution to the Gibbs free energy minimization problem is reached using a two-level strategy. On the second level, and in a systematic search, different phases are proved in equilibrium. This level proposes, systematically, seven mono/multiphase systems: vapor, liquid, liquid-liquid, vapor-liquid, liquid-liquid-liquid, liquid-liquid-vapor and liquid-liquid-liquid-vapor. The algorithm activates or deactivates the phases according to the value of the variables in the vector of integers. The binary integer values are organized into a vector, which establishes the presence or absence of the phases, where 1 signifies the presence of a phase and 0 signifies its absence. The vector has four positions, with the first one belonging to the vapor phase. If necessary, it is possible to increase the number of phases, creating more positions in the vector. The aim of this level is transform the MINLP problem into NLP subproblems.

The second level, which solves the NLP subproblems, has three main subroutines: one in charge of randomly generating the initialization of the compositions of the phases, a deterministic, and a stochastic algorithm. Initially, 5 different continuous variable initializations are proven, solving the problem with a Sequential Quadratic Programming -SQP algorithm. If at least one solution is different from the others, the problem is considered to have multiple local optimums and goes on to operate the stochastic algorithm; Simulated Annealing -SA. In this situation, the SA is in charge of administrating the initialization subroutine and the acceptance or rejection of the objective function obtained by the deterministic algorithm. The accepting or rejecting is determined by the Metropolis criterion and it depends on the parameter called "temperature", which decreases periodically during the optimization. If multiple local optimums are detected, the process continues until some predetermined stopping criterion has been met: no more improvement in the objective function or maximum number of move attempts. A stopping criterion of a maximum of 2000 iterations, or 10 iterations without evolution of the optimal objective function, was established. These last two parameters were chosen taking into account that it was necessary to have enough iterations to ensure the global solution, without excessively increasing computation time. If there are not multiple local optima, the first solution of the SQP is enough. Fig. 1 presents a diagram of the proposed strategy.

Figure 1. Solution strategy diagram

It is important to say that Simulated Annealing algorithm needs to specify three parameters that play a determining role in its efficiency: initial temperature, the number of moves per temperature level and an exponential temperature schedule a. In order to determine the initial temperature, the proposal of Pibouleau et al. (2005) is used. Temperature is from the same order of magnitude as the initial objective function, being defined by:

T_o = f(s)/A (9)

where f(s) represents the initial objective function value, A is a parameter that depends on the initial acceptance probability wanted. An initial temperature value of 1.0 is proposed in this paper. For the moves per temperature a maximum number of 100 iterations was established. The exponential temperature has been chosen, with a coefficient a of reduction:

T_i-1 = αT_i. (10)

Standard values for this coefficient are placed in the interval [0.75-0.95] (Pibouleau et al., 2005). In this paper the proposed value is 0.8.

The mathematical implementation was done in MATLAB© 7.10.0. The Simulated Annealing algorithm was programmed directly in a mathematical environment, while for the SQP deterministic algorithm the tool provided in the software package was used. The programming included the creation of independent functions that generate constraints according to the number of integer variables. The program allows this mathematical structure to be worked with quite easily because of its matrix management. With this approach to the optimization problem it is easy to extend or modify the conditions of the phases present and obtain results for new systems under study.

A more detailed description of an SQP algorithm is included in Appendix A. The description of the Simulated Annealing algorithm is included in the Appendix B.

IV. CASE STUDY - RESULTS

This study includes the implementation of the Water-Ethanol-Cyclohexane and Water-Ethanol-Glycerin systems.

A. Water-Ethanol- Cyclohexane case

This example is extracted from the article by Reneaume et al. (1996). In the aforementioned document, the authors study seven possible phase configurations, and the stationary points found in each configuration. In order to find the stationary points, the Newton homotopy method is used, retaining the lowest value of the total Gibbs free energy as a solution.

In this paper, in order to prove the constraint equations of the optimization problem, a test was done using the values reported by the authors as optimum. In this test, the optimization problem became a problem with zero degrees of freedom, where the nature and composition of the phases was known. In order to establish the activity coefficients, the NRTL thermodynamic model and the binary interaction parameters reported by Reneaume et al. (1996) were taken. With these modifications, it is possible to compare the value of the objective function to the value of the objective function reported in the published works. It is also possible to verify if the constraints are properly formulated. Reneaume reports the Table 1 optimum solution, a three-phase system: one vapor and two liquids.

Table 1. Molar compositions and objective function optimization results (Reneaume et al., 996)

When the molar composition values are introduced in the formulated model, the equations are numerically equal to zero and a value for the objective function of -1.2290 is obtained, which implies an error of 0.02% in relation to what Reneaume et al. (1996) reported.

A. Optimization

Optimization was done taking into account the seven possible systems previously mentioned and using the algorithm explained in numeral three. The initialization values of the continuous variables are randomly generated. The lowest value of the objective function belongs to the system with one vapor phase and two liquid phases, as reported in Reneaume et al. (1996). However the minimum found is significantly lower [8%] (Table 2). This solution involves a composition of different phases.

Table 2. Compositions and objective function optimization results

Table 3 shows a comparison of the minimums obtained in this study and those reported in the article mentioned above.

Table 3. Comparison of optimized values of objective function

It is possible to observe, that even when the homotopy method tries to predict the inflection points, the method is not able to find every single one and in some sensitive cases, ignores the global minimum value, which is really to solve the problem.

B. Water-Ethanol-Glycerin Case

This case is especially interesting because of the presence of azeotrope in the Water-Ethanol equilibrium and the possibility of its elimination by extractive distillation using glycerin as an extracting agent (solvent). The constants of the Antoine equation, as well as the parameters for the NRTL model were taken from the Aspen Plus database (Aspen Tech®). The optimization has the same structure used in the previous system. It allows a maximum of three liquid phases in equilibrium with a vapor phase, therefore there are four binary variables and 24 continuous variables. Table 4 shows the results of each NLP problem for a temperature equal to 363.15 K at a pressure equal to 101.325 kPa. For initialization values when there is just one phase, molar fractions of 0.3, ethanol of 0.4 and glycerin of 0.3 were taken. When there is more than one phase, the initializations of the compositions are random and the only requirement is that it fulfills the global balance for each component.

Table 4. Solution Results NLP optimization problems with SQP

Table 4 has the results obtained. The minimum corresponds with the distribution of two phases in equilibrium: one vapor phase and one liquid phase. The optimal condition refers to the sum of the Lagrangian derivatives with regard to the variables and the Lagrange multiplier at the point of optimality. This condition is necessary in order to obtain a real solution.

The result is consistent with what was observed experimentally by Lee and Pahl (1985). As an additional test, the minimum obtained in the two-phase system is confirmed based on three different initializations. The first one starts at the value of the compositions that were obtained with the system solved by the traditional Gamma-Phi method using the same acceptance tolerances of the response that were used for the optimization. The other two initializations are chosen randomly because they only fulfill the global mass balance. Table 5 shows initializations used and the objective function value that each initialization reports without being optimized.

Table 5. Liquid-Vapor Subproblem Initializations

The same result shown in Table 4 is reached with all the initializations for the two-phase system, with objective function variations only in the sixth decimal digit.

Taking into account the interest of knowing the thermodynamic pseudo-binary behavior of the Water-Ethanol-Glycerin system, multiple optimizations were made at different temperatures and at the pressure at 1 atmosphere. In order to make a pseudo-binary diagram, a restriction was added. This restriction sets the ratio of glycerin (solvent) in respect to the other two components (water-ethanol). This new restriction reduces the degrees of freedom, but maintains the essence of the problem as an optimization problem. Figure 2 shows the results obtained with a fixed ratio of Solvent/(Water-Ethanol) and it compares them with the results obtained experimentally by Lee et al, and with the results calculated theoretically by the classic Gamma-Phi method. It is important to note that all the points in the figure are evaluated with a glycerin-free base.

Figure 2. Pseudo-Binary Water-Ethanol/Glycerin P=1 atm. diagram.

V. DISCUSSION OF RESULTS

As a computational tool a conventional PC with a 2GHz Pentium dual processor was used. The average computational times that each NLP subproblem require, taking into account the five basic initializations to verify the presence of local optima and without to stop operating the stochastic algorithm, are shown in Table 6.

Table 6. Solution Times for NLP Subproblems

Total computational time was 4512 seconds for the first case and only required the use of the stochastic method in two of the seven combinations of proposed phases. For the second case, the computational total time was 5420 seconds and also it required the stochastic algorithm intervention in two cases: in the Liquid-Liquid subproblem and in the subproblem Liquid-Liquid-Liquid-Vapor subproblem solution. The computational times are greater than in several of the papers conducted previously, but it is important to consider the elements incorporated that help to find the global optimum and the largest number of proven phases.

In the two cases studied, results were obtained that indicate, presumably, that the global minimum of the thermodynamic optimization problem proposed was reached. Because of the highly non-linear nature of the problem, it is necessary to take several precautions in order to increase the probability of finding the proper response. These precautions include the verification of the solution starting with various initializations and the combined use of algorithms that complement each other, as presented in this article. With this Water-Ethanol-Cyclohexane system, it was shown that it is not enough to use one algorithm that tries to predict global solutions.

The second case also proves, that even when the classic Gamma-Phi method delivers an apparently proper response, it is necessary to evaluate if the response corresponds to one of the multiple roots that solve the non-linear algebraic equations system. It is not excessive to recommend the use of different resolution strategies in order to obtain the result of a phase system in equilibrium.

The second case tests the solution capacity of the proposed strategy. The system is highly non-linear and the nature and composition of the phases is not known beforehand. The results are relatively reliable for finding a global optimum, if the experimental data obtained for this system is compared. It is important to emphasize that the alignment of the data as a product of the optimization, with respect to the experimental data, does not represent an advantage of the proposed strategy over the standard Gamma-Phi method. In order to make these kinds of statements, it is necessary to carry out a numerical analysis. This analysis must ensure that the binary interaction parameters are adjusted, starting with the experimental data, which the strategies will be compared against.

It is important to note that in the cases that were studied, many tests of the objective function were not required. The computational times were reasonable when compared to the ones reported by the authors (McDonald and Floudas, 1996; Zhu & Xu, 1999; Inoue and Zhu, 2001), taking into account the hardware conditions available in each case.

VI. CONCLUSIONS

A strategy is proposed for the calculation of the phase equilibrium and the phases for non-ideal systems, using the approximation of the Gibbs free energy minimization. Two examples illustrate the use and the efficiency of the method in order to find phases and compositions in the equilibrium. The strategy was implemented in systems where the liquid phases show a high deviation from the ideal behavior. It is important to mention that the method converged in all the cases and solutions were found with high probabilities of being global.

The strategy presents several elements that when combined look to obtain global solutions, minimizing the probability of falling into improper local solutions to solve the thermodynamic problem proposed, without showing convergence limitations.

It is feasible to propose optimization strategies for phase equilibrium calculations of non-ideal systems from the approximation of the Gibbs free energy minimization using software packages commercially available and easy to acquire. Few lines of code are required in order to produce a solid program, a program capable of calculating the equilibrium of multiple phases and components.

NOMENCLATURE

α: Reduction algorithm Simulated Annealing coefficient.

g_ij: i-compound Gibbs free energy in the j phase (dimmensionless).

n_ij : i-compound moles number in the j phase (mole).

n^T_i: i-compound total number of moles in the system (mole).

m_ij(nij): i-compound chemical potential in the j phase.

P_i^sat: i-compound saturation pressure (atm).

g_i: i-compound activity coefficient in the liquid phase (dimensionless)

C: Components total number

F: Phase total number

R: Gases universal constant (J/mol K).

T: Absolute system temperature (K).

P: System absolute pressure (atm).

Subscript

i: Subscript components.

j: Subscript phases.

Superscript

v: Superscript vapor phases.

l= Superscript liquid phases.

T: Matrix or vector transposition.

Abbreviations

MINLP: Mixed Integer Non Linear Programming.

NLP: Non Linear Programming.

SQP: Sequential Quadratic Programming.

s.t.: Subject to

min: Minimize value

Mathematical

(x,l): Objective Function transformed with Lagrange multiplier.

λ^T: Lagrange multiplier.

∇: Differential operator.

∇²_xx: Laplacian operator

APPENDIX A. SEQUENTIAL QUADRATIC PROGRAMMING (SQP)

This is one of the methods widely used in non-linear systems with continuous variable (Nocedal and Wright, 2006) resolution. Consider the following problem representing an optimization of a function with restrictions:

min f(x) (1a)

s.t. c(x) (2a)

where f(x) is a function of Rⁿ→R and as long as the objective function as well as the restrictions fulfill the condition of being differentiable functions. It is possible to approximate the model described by equations 1a an 1b to a quadratic programming problem (Nocedal and Wright, 2006); Although in order to do this approximation, it is necessary to apply optimality conditions of first order based on the Lagrange function (Eq. 3a).

(x,λ) = f(x) - λ^T c(x) (3a)

On the other hand, for the restrictions it is possible to obtain the Jacobian matrix:

J(x)^T = [∇c₁(x), ∇c₁(x), ..., ∇c_m(x),] (4a)

where c_i(x) is the i-th component of vector of restrictions (corresponding to the i-th restriction). Equations 3a and 4a can transform into a second problem where, under this same problem resolution, the fulfillment of the first order optimality conditions (Nocedal and Wright, 2006) is ensured. The new system will be:

(5a)

Therefore, the Jacobian of this new problem will be:

(6a)

A general iteration using Newton's method:

(7a)

Applying the iteration to equations 5a and 6a it will be:

(8a)

The resolution of the iterations of the system of 7 equations approaches a resolution by quadratic programming in the following manner:

	(9a)
	(10a)

Assuming that the problem has a unique solution (p_k,l_k) that fulfills these equations:

	(11a)
	(12a)

Vectors p_k and l_k will be identified as the solution for the system of Eq. 8a (It is good to remember that the solution to this problem is clearly algebraic). And so the convergence algorithm is:

a. Choose an initial value of (x₀,l₀)

b. Repeat until the convergence will be reached with the following steps:

b.1. Evaluate f_k, ∇f_k, ∇²_xxk, c_k and J_k.

b.2. Obtain the respective values of: λ_k+1 and x_k+1

b.3. Evaluate f_k+1, ∇f_k+1, ∇²_xxL_k+1, c_k+1 and J_k+1

Convergence is reached when the values λ_k+1 and x_k+1 are close to the values λ_k and x_k

APPENDIX B. SIMULATED ANNEALING (SA)

The Simulated Annealing algorithm is inspired by the process used in metallurgy. Metal is brought to a high temperature. The metal is subsequently cooled down. In this way the internal energy is minimized. If the cooling speed is rapid, the material is weak (crystalline state corresponding to a local minimum of internal energy). Cooling must be done in steps. Each stage is a step where the constant temperature is kept until the equilibrium is established. Once the equilibrium is established the temperature is lightly decreased. In this way the thermodynamic state changes and a new equilibrium is reached. The process continues until reaching a minimum internal energy. The best lowest internal energy value is obtained (in other words the global minimum). This value corresponds with a better quality material. For a given temperature (T), the probability of an internal energy system existence (U) is proportional where K_b is the Boltzmann constant. The probability of passing from a weak energy level to a stronger one is possible, but decreases with the temperature. With this procedure of "jumps" it is possible to avoid probable internal local energy points.

From the point of view of the algorithm, optimization variables are successively and randomly perturbed at each iteration (temperature step). For each perturbation, the objective function (energy) is evaluated. If the perturbation leads to an object function decrease, it is applied to an actual solution. Otherwise, it is accepted as a probability that depends on the system temperature. Metropolis criterion signals that at a high temperature, the system is almost free to move inside the solution space: the probability of accepting a perturbation that increases the objective function is high. At a low temperature, the probability of accepting a perturbation that increases the objective function is low. Only perturbations that lead to a decrease of objective function are accepted.

ACKNOWLEDGEMENT
Authors are grateful to Colciencias (Science administrative, Technology and Innovation Department) for its collaboration. Colombia. Research Project code: 1101-452-21113.

REFERENCES
1. Inoue, K. and Y. Zhu, "Calculation of chemical and phase equilibrium based on stability analysis by QBB algorithm: application to NRTL equation," Chemical Engineering Science, 56, 6915-6931 (2001).         [ Links ]
2. Kojima, K., K. Ochi and Y. Nakazawa, "Relationship Between Liquid Activity Coefficient and Com-position for Ternary Systems," International Chemical Engineering, 9, 342-347 (1969).         [ Links ]
3. Lee, F.-M. and R.H. Pahl, "Solvent Screening Study and Conceptual Extractive Distillation Process To Produce Anhydrous Ethanol from Fermentation Broth," Industrial & Engineering Chemistry Process Design and Development, 24, 168-172 (1985).         [ Links ]
4. McDonald, C. and C. Floudas, "Glopeq: a new computational tool for the phase and chemical equilibrium problem," Computers& Chemical Engineering, 21, 1-23 (1996).         [ Links ]
5. Nichita, D. V., Gomez, S., and Luna, E.. "Multi-phase equilibria calculation by direct minimization of Gibbs free energy with a global optimization method". Computers and Chemical Engineering, 26, 1703-1724(2002).         [ Links ]
6. Nocedal, J. and S. Wright, Numerical Optimization, 2^nd Ed., Springer (2006).         [ Links ]
7. Parodi, C.A. and E.A. Campanella, "Effect of vapor-liquid equilibrium data on the design of separation sequences by distillation," Latin American Applied Research, 39, 33-40 (2009).         [ Links ]
8. Pibouleau, L., S. Domenech, A. Davin and P. Azzaro, "Expérimentations Numériques sur les variantes et paramètres de la mèthode du recuit simulé". Chemical Engineering, 105, 117-130 (2005).         [ Links ]
9. Rangaiah, G.P., "Evaluation of genetic algo-rithms and simulated annealing for phase equilibrium and stability problems," Fluid Phase Equilibria, 187, 83-109 (2001).         [ Links ]
10. Reneaume, J.M., M. Meyer, J.J. Letouneau and X. Joulia, "A global MINLP approach for phase equilibrium calculations," Computers & Chemical Engineering , 20, 303-308 (1996).         [ Links ]
11. Revollar, S., M. Francisco, P. Vega and R. Laman, "Stochastic optimization for the simultaneous synthesis and control system design of an activated sludge process," Latin American Applied Research, 40, 137-146 (2010).         [ Links ]
12. Saber, N. and J.M. Shaw, "Rapid and robust phase behaviour stability analysis using global optimi-zation," Fluid Phase Equilibria, 264, 137-146 (2008).         [ Links ]
13. Seader, J.D. and E.J. Henley, Separation Process Principles, 2^nd Ed., Wiley (2006).         [ Links ]
14. Zhu, Y. and Z. Xu, "A reliable prediction of the global phase stability for liquid-liquid equilibrium through the simulated annealing algorithm: Application to NRTL and UNIQUAC equations," Fluid Phase Equilibria, 154, 55-69 (1999).         [ Links ]

Received: August 4, 2011
Accepted: April 19, 2013
Recommended by Subject Editor: Orlando Alfano