Application of Solution Strategies for Numerical Estimation of Thermodynamic Equilibrium Parameters for an Acetone–Butanol Mixture

: The development of reliable numerical estimation of thermodynamic parameters is a crucial aspect in the ongoing research about process engineering and design. The consideration of these concepts lets to design more precise processing units and separations stages based on the predicted nature of substances. Therefore, this study presents an application of di ﬀ erent solution methods for the estimation of thermodynamic equilibrium parameters of an acetone–butanol mixture. This dissolution is a non-ideal system, so, the non-ideal Raoult’s Law and Wilson’s equation were used to model the liquid–vapor equilibrium. Otherwise, the solution of this system required the application of nonlinear least squares (NLS) for determination of adjustable parameters. As the above step transformed Wilson’s equation into a system of nonlinear equations, solution algorithms such as; Newton–Raphson method (NRM), Broyden’s method (BM) and Levenberg–Marquardt method (LMM) were applied. All algorithms converged towards the same solution ( Λ 12 = 0.689 and Λ 21 = 0.798), but Newton’s and Broyden’s methods employed fewer computational time and number of iterations compared to performance showed by the Levenberg–Marquardt algorithm.


Introduction
Nowadays, there is a rising interest in modeling different environments and systems in terms of sustainable production regarding the use of renewable feedstock for transformation of these sources into value-added chemicals [1]. The sustainable design approach arises, considering the described motivation as a baseline framework for generating sustainable alternatives based on computer-aided process engineering (CAPE) tools. In this sense, the development of the principles of this approach requires performing reliable estimations of physicochemical parameters in order to obtain more accurate and precise calculations considering the prediction of substance properties [2]. The above consideration is essential because many transformation processes usually handle a variety of chemical components, including solvents and biofuels such as ethanol, butanol, acetone-among others-that present very complex behaviors. Thus, the modeling of mixture equilibria becomes a key aspect in the design of downstream units for this type of plant [3]. Carlon [4] mentioned that an essential task on process simulation is to find the adequate property parameters for equilibria modeling, taking into account that missing calculations or mistakes in the selection of estimation method lead to generating not trustable or reliable data. In this regard, there are several methods and equations reported in the literature for the estimation of thermodynamic properties of substances. These include Raoult's law, Margules', Wilson's or Van Laar equations, among others [5]. Regarding the beforementioned point, the Redlich-Kwon equation of state was used to model the solubility of methane in water under different pressure and temperature conditions [6]. Bader et al. [7] studied the influence of non-ideal vapor-liquid equilibrium on the evaporation operation of ethanol/iso-octane blends considering the importance of measuring thermal equilibrium in pretreatment systems of second. Galeotti et al. [8] estimated a vapor-liquid equilibrium in the ternary system formed of acetic acid, water and glucose or xylose. It is an interesting contribution -generation biorefineries. Kobuchi et al. [9] described mathematical expressions for correlation of vapor-liquid equilibria using the Wilson equation for calculating solubility and molar volumes. Other studies have shown interest in studying multi-component solvents (such as butanol, acetone, ethanol and water) mixtures under no-ideal thermodynamic models. In this regard, Kaminski et al. [10] modeled liquid-liquid equilibrium in a multi-component system comprised of acetone, butanol, ethanol and ionic liquid. Another issue that is needed to be considered is that the formulation of vapor-liquid equilibria generates systems of nonlinear equations that require the programming of solution strategies. In this sense, the Newton-Raphson, Regula-Falsi, Broyden or Levenberg-Marquardt methods are the most used algorithm for solving frequent problems about a system of nonlinear equations. Abderafi and Bounahmidi [11] measured the vapor-liquid equilibrium for industrial sugar juice using the Peng-Robinson equation of state using the mean square error and Davidson's methods. The last one is a quasi-Newton algorithm. A continuum framework for modeling liquid-vapor interfaces was proposed by Chandra et al. [12], using a locally discontinuous arbitrary Lagrangian-Eulerian finite element formulation. It is worth mentioning that the selection of an appropriate strategy for analyzing and solving problems of vapor-liquid equilibrium is crucial in the design of chemical processes through computer-aided tools, emphasizing the relevance to obtain accurate solutions, short computational times and high converge rates. This fact gains attention highlighting that the complex problems of sustainable design under CAPE tools just not involve the estimation of thermodynamic properties but includes calculation of chemical reactions, mass/energy balances, downstream processing, among others [13].
Under the beforementioned, this study presents the application of different solution strategies based on nonlinear modeling for the numerical estimation of thermodynamic equilibrium parameters and variables of the acetone-butanol mixture. The novelty of this work is the programming and comparison of three methods for the solution of the system of nonlinear equations derived from the formulation of vapor-liquid equilibria of the acetone-butanol system under the Wilson model. It will provide insights about the utilization of the available solution strategies considering computational times and converge rates for future application of more complex systems under computer-aided process engineering problems.

Materials and Methods
The development of this study used reported thermodynamic data in the literature of the acetone-butanol mixture (12). Table 1 shows data for acetone-butanol mixture regarding different pressures. x 1 and y 1 represent the molar concentration in liquid and vapor phases for component 1 (Acetone), respectively. Considering the physical-chemical properties and chemical affinity of the mixture, the Wilson model used to estimate the thermodynamic parameters as follows in Equations (1) and (2) [14]: Λ 12 is the interaction parameter of component 1 to component 2, x 2 is the mole fraction of component 2 in the liquid phase, y 2 is the molar fraction of component 2 in the vapor phase and Λ 21 is the interaction parameter of component 2 to component 1. Finally, γ 1 and γ 2 are the activity coefficients of components 1 and 2, respectively. Considering the Wilson model showed in Equations (1) and (2), the following step was to estimate the activity coefficients (γ 1 and γ 2 ) of mixture component. For this, the Modified Raoult's Law was used. This equation [see Equation (3)] is applied for binary non-ideal vapor-liquid equilibrium systems, as follows [7]: y i and x i are the concentration of component i in the vapor and the liquid phase, respectively. F i is the fugacity coefficient, P and P sat i are the system and saturation pressure, respectively. Next, the modeling of the binary system required the estimation of the fugacity coefficient for both components, according to Equation (3). This estimation was developed by using the virial equation of state described by Equations (4)-(10) as follows [5]: T r ij , w ij , Bij 0 , Bij 1 , Bij are the crossed reduced temperature, acentric factor, zero-order virial coefficient, 1st order virial coefficient and virial coefficient for a binary mixture, respectively. w i , T i and w j , T j are the acentric factor and reduced temperature of component i and component j, respectively. Once all required thermodynamic data were estimated, the solution of the nonlinear system composed of Equations (1) and (2), is formulated from the application of the nonlinear least squares (NLS) [15], taking into consideration that Λ 12 and Λ 21 are the adjustable parameter of the nonlinear system. Equation (11) describes the formulation of the NLS generatrix equation as follows: G is the generatrix function of NLS, F ex is the evaluated function, while F theo the theoretical nonlinear function and N represents the number of points. In the case of Wilson's equations, two generatrix functions were formulated as described by Equations (12) and (13): Lnγ 1ex and Lnγ 2ex represent the experimental data for the natural logarithm of the activity coefficient of component 1 and 2, respectively. G 1 and G 2 are the generatrix function for components 1 and 2, respectively. The Λ 12 and Λ 21 are the adjustable model parameters for component 1 and component 2, respectively. Next, the method involved the application of partial derivatives for each generatrix function to each adjustable parameter. In this case, such parameters are the unknows of the equation system. Regarding Equations (12) and (13), each generatrix function was derived from Wilson's equation parameters, so it gives a nonlinear system of four equations and two unknowns. Equation (14) presents this system of equations as follows: F represents the derivation matrix (with four derivates) and X is the adjustable parameter vector (with two parameters). In the literature are reported several methods for solving nonlinear equation systems, these solution strategies include; Newton´s method [16], Differential homotopy analysis [17], Broyden's method [18], Regula-Falsi method [19], tangent method [20], among other. This study used Newton's, Broyden's and Levenberg-Marquardt methods were to solve the system of nonlinear equations derived from the formulation of the NLS.

Newton's Method
Newton's method, which is also called the Newton-Raphson method (NRM), is an iterative method or numerical approach for attempting to find the roots of a function or system of equations [21]. NRM is formulated as an expansion of Taylor's series [22], as follows: f(x) and f (x) represent the function and the first order derivate evaluated in x o . If the Equation (15) is truncated in the first derivate and also, the system is extended to a matrix of nonlinear equations and an unknown vector is considered, then Equation (15) is formulated as follows: F (x) and F (x o ) are the matrices of the equation evaluated at x and x o vector, respectively. J(x o ) is the Jacobian matrix evaluated at x o vector. Rearranging Equation (16), and establishing an iteration system, the algorithm gets the calculation of vector x at the current k iteration as follows: where x k and x k−1 are the vector evaluated at k and k-1 point, respectively. Equation (17) represents the generalization of Newton-Raphson´s method in an iterative form.

Broyden's Method
Broyden's method (BM) is proposed based on NRM, which is generalized from Taylor's series. Moreover, BM comes from the incorporation of a damping loop in order to avoid losses of computational time (CT) evaluating inverse Jacobian iteration [23], as follows: Equation (18) is incorporated in the iteration system shown in Equation (16) and the evaluation of the actual iteration based on Broyden's modification is stated as follows:

Levenberg-Marquardt Method
As developed by Broyden´s (and Newton´s) method, the Levenberg-Marquardt method (LMM) is formulated from Taylor's series and is a combination of the gradient descendent rule and Newton's method [24]. The LMM employs a parameter for controlling the step-size, which takes enormous values in the first iterations (equivalent with the gradient descent algorithm) and small values in the latest ones (equivalent with the Gauss-Newton method). Therefore, the LMM function generalizes the system of equations as follows in Equation (20): I is the identity matrix, J T (x o ) is the transpose of the Jacobian and λ is the damping parameter for the LMM function. Moreover, Equation (21) defines the expression J(x o ) J T (x o ) as a matrix called "Hessian" as follows [25]: Therefore, the LMM iteration routine solves the system as follows in Equation (22): This method involves a local optimization to find the optimal value of λ k−1 that approximates x k−1 to x k , thereby achieving a faster conversion in the numerical method. A primary error E x k−1 = x o is evaluated, so, an initial value of λ k−1 = λ o = 0.01 is estimated. Then, the system is solved according to the LMM to obtain a new solution vector x k , and for this iteration the error E LM x k is re-evaluated, therefore; it states a local search as follows: Finally, a new solution vector is assigned to be x k−1 = x k and λ k−1 = λ k , likewise, the cycle is restarted until the allowed tolerance is reached. Such condition is established according to the discrepancy error set (for all used algorithms) to be 10 −7 , Equation (23) displays the error function (E) this as follows:

Results
Based on thermodynamic and state equations, the required parameters for estimating the liquid-vapor equilibrium for the mixture were calculated by using Equations (4)-(10). This system considered a temperature of 353.15 K, along with pressure conditions that vary between saturation pressures of both components. Table 2 shows the thermodynamic properties and parameters of pure components of the mixture. Such values were needed to establish the system for Wilson's equations. The data generated from the application of this analysis can be mostly used (in process design) to design distillation towers and separation systems. For such purpose, the system can deal with the intrinsic measurement and experimental errors, considering further process optimization (or intensification) and through the use of safety design factors [26].  (24) and (25).
where λ 12 − λ 11 and λ 21 − λ 22 are the interaction Wilson parameter for components 1 and 2, respectively, while v 1 and v 2 are the pure component molar volume of components 1 and 2 [27].

Estimation of Thermodynamic Properties and Parameters
By using the information provided by Table 1, Table 2 and the modified Raoult's Law, it was constructed the data for fugacity and activity coefficients for acetone and ethanol, as reported in Table 3. Moreover, the above step allows the definition of the NLS system for determining the adjustable coefficients (Λ 12 and Λ 12 ) of Wilson's equation. Such nonlinear equations obtained from the application of the NLS are solved by using the methods described in Section 3.1, Section 3.2, Section 3.3 and the results regarding such solutions are showed and discussed in the following sections. Table 3. Estimation of fugacity and activity coefficient for acetone-butanol mixture.

Solution by Newton's Method
The Newton-Raphson algorithm was implemented to determine the adjustable parameters of Wilson's equation for the binary system presented in this study. In addition, the convergence speed of the model is verified based on an estimation of the discrepancy error, which is commonly used for Newton's method [28]. This discrepancy function can be estimated as the Euclidean norm (E(x) p ) of the error vector as follows: p represents the maximum absolute row sum of the unknows vector. Simulations were developed, taking into consideration equations derived from the problem presented in this study. Moreover, Equation (27) helped to estimate the iteration error. It was also applied for both Broyden's and Levenberg-Marquardt methods. This study performed a step by step algorithm based on NRM to solve the described problem. Figure 1 shows the algorithm for attempting to find the adjustable binary parameters of acetone-butanol non-ideal mixture based on NRM. Otherwise, all performed numerical approaches and algorithms for simulations were developed in an Intel ® Core ™ i7-7500 CPU @ 2.70 GHz with 8 Gb of RAM. As shown in Figure 1, an initial value vector has to be set for solving the nonlinear system. Such consideration is a key-aspect to guarantee convergence in the solution and low Appl. Sci. 2020, 10, 3136 8 of 16 computational times. Therefore, the values established for the initial vector are assigned depending on the phenomenon's nature. This fact may be relevant is due to iterative numerical methods are further sensible respect to the initial vector. Table 4 shows the results for simulations based on the Newton-Raphson method. Therefore, Table 4 (and the study presented here) just analyzed such runs that obtained reliable results. Based on results showed in Table 4; the developed Newton-Raphson algorithm converged towards the same solution (x = [0.688; 0.798]) for simulation 2 and 3 while simulation 1 diverged. This outcome means that the nonlinear system analyzed in this study is highly sensible concerning the initial vector. Moreover, very deviated values of initial value compared to the roots of the equations could cause a destabilization of the algorithm. Figures 2 and 3 show the behavior of the discrepancy error and the solution vector for simulations 2-3, respectively.  Therefore, Table 4 (and the study presented here) just analyzed such runs that obtained reliable results. Based on results showed in Table 4; the developed Newton-Raphson algorithm converged towards the same solution (x = [0.688; 0.798]) for simulation 2 and 3 while simulation 1 diverged. This outcome means that the nonlinear system analyzed in this study is highly sensible concerning the initial vector. Moreover, very deviated values of initial value compared to the roots of the equations could cause a destabilization of the algorithm. Figures 2 and 3 show the behavior of the discrepancy error and the solution vector for simulations 2-3, respectively.

Solution by Broyden's Method
The modified Newton-Raphson algorithm with Broyden's constraint was implemented to determine the adjustable parameters of Wilson's equation for the binary system presented in this study. For this purpose, an algorithm was formulated to solve the NLS along with the nonlinear system generated from it. Figure 4 shows the proposed algorithm based on the described Broyden's method. As shown in Figure 4, an initial value vector has to be set for solving the nonlinear system. Such consideration as a key-aspect to guarantee convergence in the solution and low computational time. Therefore, the values set for the initial vector are assigned depending on the phenomenon. Table  5 shows the results for simulations based on Broyden's method.

Solution by Broyden's Method
The modified Newton-Raphson algorithm with Broyden's constraint was implemented to determine the adjustable parameters of Wilson's equation for the binary system presented in this study. For this purpose, an algorithm was formulated to solve the NLS along with the nonlinear system generated from it. Figure 4 shows the proposed algorithm based on the described Broyden's method. As shown in Figure 4, an initial value vector has to be set for solving the nonlinear system. Such consideration as a key-aspect to guarantee convergence in the solution and low computational time. Therefore, the values set for the initial vector are assigned depending on the phenomenon. Table  5 shows the results for simulations based on Broyden's method.

Solution by Broyden's Method
The modified Newton-Raphson algorithm with Broyden's constraint was implemented to determine the adjustable parameters of Wilson's equation for the binary system presented in this study. For this purpose, an algorithm was formulated to solve the NLS along with the nonlinear system generated from it. Figure 4 shows the proposed algorithm based on the described Broyden's method. As shown in Figure 4, an initial value vector has to be set for solving the nonlinear system. Such consideration as a key-aspect to guarantee convergence in the solution and low computational time. Therefore, the values set for the initial vector are assigned depending on the phenomenon. Table 5 shows the results for simulations based on Broyden's method.
As mentioned, due to Broyden's method is an iterative solution strategy that requires assuming the parameter of the initial solution vector, several simulations were developed in order to search the best performance. Based on results showed in Table 5, it can be observed that the developed Broyden's algorithm is stable and also it always converged towards the same solution (x = [0.688; 0.798]), which were the same for Newton's method (in simulation 2 and 3). In addition, these results indicate that when the values of the initial vector are decreased around the range 0.5 to 1 (for both parameters), the speed of convergence rises, obtaining better computational times. The above is supported by considering the changes generated in simulations 2 and 3 respect to the first one. Globally, this algorithm is stable and achieved convergence employing a small number of iterations. Figures 5-7 shows the behavior of the discrepancy error and the solution vector for simulations 1-3, respectively.  As mentioned, due to Broyden's method is an iterative solution strategy that requires assuming the parameter of the initial solution vector, several simulations were developed in order to search the best performance. Based on results showed in Table 5, it can be observed that the developed Broyden's algorithm is stable and also it always converged towards the same solution (x = [0.688; 0.798]), which were the same for Newton's method (in simulation 2 and 3). In addition, these results Y i , X 1 i , X 2 i , … , X m , E (y) = f(x 1 , x 2 , … , x m ; β 1 , β 2 , … , β k )   indicate that when the values of the initial vector are decreased around the range 0.5 to 1 (for both parameters), the speed of convergence rises, obtaining better computational times. The above is supported by considering the changes generated in simulations 2 and 3 respect to the first one. Globally, this algorithm is stable and achieved convergence employing a small number of iterations. Figures 5-7 shows the behavior of the discrepancy error and the solution vector for simulations 1-3, respectively.
(a) (b)   indicate that when the values of the initial vector are decreased around the range 0.5 to 1 (for both parameters), the speed of convergence rises, obtaining better computational times. The above is supported by considering the changes generated in simulations 2 and 3 respect to the first one. Globally, this algorithm is stable and achieved convergence employing a small number of iterations. Figures 5-7 shows the behavior of the discrepancy error and the solution vector for simulations 1-3, respectively. indicate that when the values of the initial vector are decreased around the range 0.5 to 1 (for both parameters), the speed of convergence rises, obtaining better computational times. The above is supported by considering the changes generated in simulations 2 and 3 respect to the first one. Globally, this algorithm is stable and achieved convergence employing a small number of iterations. Figures 5-7 shows the behavior of the discrepancy error and the solution vector for simulations 1-3, respectively.

Solution by Levenberg-Marquardt Method
As was developed by the above section, the problem regarding NLS was solved by using the Levenberg-Marquardt method. Figure 8 displays the implemented algorithm based on the Levenberg-Marquardt method. Simulations were developed considering the same initial vectors for the previous solution method. Table 6 shows the results for the three main simulations developed by using the Levenberg-Marquardt method.
As was developed by the above section, the problem regarding NLS was solved by using the Levenberg-Marquardt method. Figure 8 displays the implemented algorithm based on the Levenberg-Marquardt method. Simulations were developed considering the same initial vectors for the previous solution method. Table 6 shows the results for the three main simulations developed by using the Levenberg-Marquardt method. J(x o ); F(x 0 ) λ 0 = 10 −2 , x o , v = 10  As same for Broyden's method, the Levenberg-Marquardt algorithm is an iterative solution method that requires the assumption of an initial x o vector. Therefore, many simulations were developed for achieving the best performance and obtain results comparable with the application of Broyden's method. Based on the results presented in Table 6, the LMM applied for the solution of the NLS problem presented in this work shows moderate positive performance. In simulation 2 and 3, this algorithm reached converged towards the same solution (Λ 12 = 0.688; Λ 21 = 0.798) that was obtained by employing Broyden's method. Taking into account results for simulation 1, LLM also shows that if the initial solution vector is far from the final values of it, the algorithm could diverge or employ several iterations (and computational time) to achieve the real solution. Figures 9 and 10 show curves for discrepancy error and solution vector respect to iterations for simulations 2 and 3 in the LMM, respectively. As same for Broyden's method, the Levenberg-Marquardt algorithm is an iterative solution method that requires the assumption of an initial x o vector. Therefore, many simulations were developed for achieving the best performance and obtain results comparable with the application of Broyden's method. Based on the results presented in Table 6, the LMM applied for the solution of the NLS problem presented in this work shows moderate positive performance. In simulation 2 and 3, this algorithm reached converged towards the same solution (Λ = 0.688; Λ = 0.798) that was obtained by employing Broyden's method. Taking into account results for simulation 1, LLM also shows that if the initial solution vector is far from the final values of it, the algorithm could diverge or employ several iterations (and computational time) to achieve the real solution. Figures 8 and 9 show curves for discrepancy error and solution vector respect to iterations for simulations 2 and 3 in the LMM, respectively.  As same for Broyden's method, the Levenberg-Marquardt algorithm is an iterative solution method that requires the assumption of an initial x o vector. Therefore, many simulations were developed for achieving the best performance and obtain results comparable with the application of Broyden's method. Based on the results presented in Table 6, the LMM applied for the solution of the NLS problem presented in this work shows moderate positive performance. In simulation 2 and 3, this algorithm reached converged towards the same solution (Λ = 0.688; Λ = 0.798) that was obtained by employing Broyden's method. Taking into account results for simulation 1, LLM also shows that if the initial solution vector is far from the final values of it, the algorithm could diverge or employ several iterations (and computational time) to achieve the real solution. Figures 8 and 9 show curves for discrepancy error and solution vector respect to iterations for simulations 2 and 3 in the LMM, respectively.

Discussion
Regarding the development of the three numerical methods performed to solve the presented system, all of them showed reliable numerical stability and quick converged time despite that the topology of the derived system of equations from the application of NLS is highly nonlinear. In the case of Newton's method, the best performance was obtained by simulation 3 (6 iterations), getting a high rate of converged. In this case, the selection of the initial vector is a crucial aspect of obtaining reliable outcomes. Broyden's method showed excellent results for all three simulations. Analyzing in detail, simulations 2 and 3 posted the best performance solving with a similar computational time (4.950 and 4.263 s) and the number of iterations (6 iterations). Otherwise, LMM showed higher computational times for simulations 2 and 3 (94.312 and 4.263 s), indicating that its converged rate was lower than the other two algorithms. In this case, both Broyden and Newton methods presented high conversion rates, which is an expected result considering that these algorithms are very similar in their mathematical structure. Otherwise, the damping parameter inserted in the Broyden cycle helps to improve the converge rate for those initial vectors that are far from the real solution. This fact is evidenced by comparing the results for simulation 1 for both methods, which ran with the same initial vector, but the obtained outcomes for these were quite different (divergence for Newton's method, while convergence in 21 iterations for Broyden's method). Once the initial vector is close enough to the final solution, these methods behave similarly. Finally, we implemented the Levenberg-Marquardt method for the problem addressed in this study. It showed very different behavior compared with the simulations made by Newton and Broyden methods. In this case, this algorithm solves simulations 2 and 3, with a low rate of convergence and it took several iterations (1370 and 1667, respectively) to get a final vector under the minimum acceptable discrepancy error. The above means that for this case study, it is not recommended the use of the Levenberg-Marquardt algorithm and through the implementation of Newton-based methods, we can get better convergence rates and more stability in the simulations. It is worth mentioning that this is an unexpected result, taking into account that the LMM is more robust than the NRM, but in some cases, the Levenberg-Marquardt can be a little bit slower than the Newton algorithm. For this problem, the lower performance may be related to the selection and optimization of its damping value [29], which in the case of Broyden's method worked very well, as the simulations demonstrated. The obtained results can be useful for further development of calculation in the design of separation processes. The above is relevant, considering that there were efforts to model this system or similar ones, but applying different equation models such as NRTL [30]. One of the main outcomes of these findings is the further use of the data for modeling separation systems, which involve the acetone-butanol mixture [31]. An example of this is the design of distillation columns for the purification stage of the acetone-butanol-ethanol (ABE) fermentation. That reaction generates a complex mixture formed by ABE components, water and other substances [32].

Conclusions
Different optimization and nonlinear equations solution algorithms were used to numerical estimate thermodynamic equilibrium parameters for an acetone-butanol mixture. This process was first modeled using Wilson's equation and the nonlinear least square method. Hence, a system of nonlinear equations was generated from the above step. The application of numerical methods studied the problem addressed from the above. This includes Newton, Broyden and Levenberg-Marquardt methods. All of the performed algorithms reached the same solution vector (Λ 12 = 0.689 and Λ 21 = 0.798), showing different behavior. Newton and Broyden methods presented a quite similar performance, posting very similar convergence rates and computational times. Otherwise, the Levenberg-Marquardt method was the algorithm with the most inferior performance, reaching the solution vector after more than 1000 iterations (for simulations 2 and 3). Likewise, the first simulation using this method diverged, indicating that this strategy is not good enough stable for the nonlinear system presents in this study. In the case of Newton-based algorithms, Broyden's method presented the best performance due to simulation 1 demonstrated that this algorithm could compensate that the initial vector is very different through the damping factor. For future work, the consideration of dealing with measurement and experimental errors in order to decrease the uncertainty in the measure and obtained more exact calculations. Moreover, the inclusion of parallel numerical approaches in the solution of this type of system may generate insights for the solution of more complex problems. Further studies also can develop computational tools to assess mixtures and thermodynamic equilibrium systems incorporating these into a decision-making strategy for the selection of the most optimum and realistic models using foundations of machine learning applied to process engineering and design.