The Parametrization Problem in the Modeling of the Thermodynamic

Correlation Cite This: Ind. Eng. Chem. Res. 2019, 58, 12876−12893

pubs.acs.org/IECR

The Parametrization Problem in the Modeling of the Thermodynamic Behavior of Solutions. An Approach Based on Information Theory Fundamentals Adriel Sosa,† Luís Fernández,† Juan Ortega,*,† and Laureano Jiménez‡ †

Downloaded via NOTTINGHAM TRENT UNIV on August 13, 2019 at 08:21:12 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Grupo de Ingeniería Térmica e Instrumentación (IDeTIC), Parque Científico- Tecnológico, Universidad de Las Palmas de Gran Canaria, 35071-Canary Islands, Spain ‡ Departament d’Enginyeria Química, Universitat Rovira i Virgili, Av. Paisos Catalans 26, Tarragona, Spain S Supporting Information *

ABSTRACT: This work shows a new approach to the parameter-fitting problem useful in the solutions thermodynamic field, providing a more objective framework to obtain better empirical/semiempirical models for chemical engineering applications. A model based on the excess Gibbs energy function gE is used to represent the behavior of real solutions, together with its first derivative hE, using a combined modeling under the paradigm of multiobjective optimization. The problem is formulated as an MINLP methodology to simultaneously consider two aspects: the model complexity and the best parametrization to prevent the overfitting, controlling the trade-off between them by applying the Akaike Information Criterion to gE residuals. Two different solvers, one deterministic (SBB/ CONOPT) and another evolutionary (GA), are used, and their ability to solve the problem is analyzed. The designed methodology is applied to three highlighted VLE cases in chemical engineering, and the results obtained show the ability of the method to get the best model in each case. The proposed methodology proved useful for modulating the number of parameters considering the imposed requirements, which decrease as the accuracy requirements for hE are relaxed. The ef f icient-f ronts obtained show a small trade-off region, noting that the proposed framework provides the simplest models with the minimum completeness uncertainty. Gibbs function models, generating studies on the EOS-gE binomial. However, owing to the complexity of some of the models, an important goal is to achieve an optimal parametrization without an excessive computational cost. The use of correlative models implies, among other questions, arbitrarily deciding on the shape of the equations to improve the fit. The simplest and most common approach is to use a least-squares procedure to represent the properties involved. If the selected model does not reproduce the real information, a possible solution is to modify the number of parameters providing it with a greater flexibility. Thus, Ko et al.12 proposed modifications for the NRTL model by improving the modeling of the functional gE = gE(x,T) for liquid−liquid equilibria, Gmehling et al.13 modified the relationships between the interaction parameters and their temperature dependence for UNIFAC, and so on. This practice is also common for the EOS, where some temperature-dependent relationships are generated for the kij binary interaction parameter.14 These modifications, and other more specific ones, can be sufficient to improve the representation of

1. INTRODUCTION Precise knowledge of the thermodynamic properties of solutions is of great scientific and technical interest. On the one hand, scientists need to accurately represent the behavior of fluid matter knowing the bases that govern the mixing processes and equilibrium phenomena. Alternatively, engineers require exact modeling to develop optimum designs of processes. In both cases, thermodynamics plays an important role as a scientific tool, since it allows the real behavior of solutions to be combined with the formalism of said science and the subsequent applications. The modeling has important repercussions on the results obtained by the engineer in the analysis, design, and optimization of the processes involved. In this context, the literature contains information on many models, some with a scientific basis,1−6 which aim to predict the behavior of solutions without previous experimentation, and others of an empirical nature,7−9 generally used as correlative and nonpredictive to get a better representation of the experimental data. Occasionally, some of these models have been redirected as correlatives of experimental information; this is the case of group contribution methods, such as UNIFAC10 and ASOG.11 However, there is no single modeling that represents the behavior of all fluid matter. Therefore, researchers have taken a different approach to develop theories by combining equations of state (EOS) with © 2019 American Chemical Society

Received: Revised: Accepted: Published: 12876

March 17, 2019 June 10, 2019 June 18, 2019 June 18, 2019 DOI: 10.1021/acs.iecr.9b01493 Ind. Eng. Chem. Res. 2019, 58, 12876−12893

Correlation

Industrial & Engineering Chemistry Research fluid behavior but can also lead to problems of data reduction. The inconvenience of using this strategy is related to the overfitting problem. This allows the model to represent well the experimental data used in the modeling but decreases its generalization capacity.15 A more comprehensive way to address the modeling is through a multiproperty approach, where several related quantities are considered simultaneously in the optimization process. An alternative is to minimize the aggregated errors of the properties used to obtain the better parameters, such as the case of the group contribution models.10,11 The final model is a trade-off solution regarding each of the involved thermodynamic quantities. (Throughout the text the modeling results are indicated as “results” (in italics), thus avoiding confusion with words such as “solutions” or “binary dissolutions” used.) If the model definition is not sufficient to describe the whole thermodynamic behavior, imbalances occur, meaning that one of the properties is better represented than the others. In this case, a multiobjective approach is required to solve the problem, defining a set of possible results mathematically equivalent (nondominated Pareto-results);15−19 the advantage of this procedure is that the model search is not biased toward one property or another. Nevertheless, the specialized literature pays little heed to defining the models, in other words the search for the appropriate functional and its parametrization. This work proposes a methodology to improve model definition using objective criteria based on information theory. It provides a framework to compute the loss of information related to the model used to represent the experimental data, which is linked to the precision and complexity of the model, i.e., the number of free parameters. The main difference with conventional practices is that, with this approach, initially only a general form of the model is specified; and the final form is deduced during the fitting procedure. Consequently, two tasks are addressed simultaneously: (a) the fitting-procedure and (b) the selection of the model, considering, respectively, the minimization of the model residuals and complexity. Task (a) is addressed by two algorithms, one deterministic20 and another evolutionary,21 comparing their ability to generate a well-distributed set of nondominated results applied to three highlighted cases of vapor−liquid equilibria (VLE) of binary systems; the Akaikés criterion22 is used for task (b).

bound vector of linear constraints. Nonlinear constraints are denoted by g(θ); the matrix F(θ) stores the differences between the observed pairs {xj, yj} and the model estimates f(xj, θ).

ij f (x1, θ) − y yz 1 z jj zz jj jj f (x , θ) − y zzz jj 2 zz 2z F(θ) = jjj zz jj zz ∂ jj zz jj z jj f (x , θ) − y zzz N N (2) k { Alternatively, another popular error metric used in problem LS1 is the root of the mean square error, RMSE or s[F(θ)]. Its relation to the previous measure of error is expressed in the following equation: ÅÄÅ N ÑÉÑ1/2 2 1/2 ÅÅ ÑÑ jij F(θ) 2 zyz 2 Å z = ÅÅ∑ (f (x k , θ) − yk ) /N ÑÑÑ s(F(θ)) = jj j |F(θ)| zz ÅÅ ÑÑ k { ÅÇ k = 1 ÑÖ (3)

where N = |F(θ)| is the number of instances in the data set. Hereon, we will refer to this latter metric in the manuscript. The parameters are tuned during the metric minimization until a desired tolerance is achieved. The problem may include constraints, setting the upper and lower limits for each parameter. This task is not trivial because, a priori, the actual boundaries of the parameters are unknown. Alternatively, if the parameter space is not specified, this leads to problems in the fitting process, such as local minima stagnation, disturbing the mathematical procedure for some optimization algorithms. The previous approach refers to the representation of a single property. However, thermodynamic quantities that represent the phase equilibria involved in separation processes have a common origin and are included in the differential, eq 4, and derived relationships, eqs 5 and 6. E E E d( g RT ) = ( v RT ) dp − ( h RT 2 ) dT +

ÄÅ ÅÅ 1 ÅÅÅÅ E ln γi = Åg − RT ÅÅÅÅ ÅÅÇ

ij ∂g E yz zz hE = g E + Ts E = g E − T jjjj zz ∂ T k { p,x

2. THE MODELING PROBLEM 2.1. Classical Approach. The most extended form for modeling data of thermodynamic properties involves minimizing the distance between the value calculated by the chosen model and the experimental values. Frequently a L2 measure is used resulting in a least-squares (LS) procedure, which is described by problem LS1. Problem LS1. min θ

F ( θ)

i

(4)

(5)

(6)

These relationships, generated by the thermodynamic formalism for Gibbs function, ensure a theoretical consistency between the different properties of the dissolution studied. Particularly, for the excess Gibbs function, gE= gE(xi.p,T), several mathematical expressions have been postulated to represent the activity coefficients of the liquid phase γi, as a function of the canonical variables. 2.2. The Multiobjective Nature of the Problem. As we can see from the previous section, the modeling can be formulated as a mono-objective problem. However, for the purposes of model completeness and due to the uncertainty resulting from data acquisition and measurement, fitting the different properties becomes a multiobjective problem.16,23 This leads to a set of results in which each of them is no better than the others, constituting a Pareto set, 7*. Since the results

2 2

s. t. Aθ ≤ b g (θ) ≤ 0 θ∈P

ÉÑ ÑÑ ÑÑ ij ∂g E yz ÑÑ j z j z ∑ xk jj zz ÑÑ Ñ x ∂ k≠i k k { p , T , xj ≠i,k ÑÑÑ ÑÖ

∑ ln γi dxi

(1)

where P represents the parameter space of dimension m. The variables of the problem correspond to the parameters (θ) of the model. A and b are, respectively, the coefficient matrix and 12877

DOI: 10.1021/acs.iecr.9b01493 Ind. Eng. Chem. Res. 2019, 58, 12876−12893

Correlation

Industrial & Engineering Chemistry Research obtained are not guaranteed to be global optima, the final set of nondominated results, see Appendix A (Def 1), constitutes the efficient-results set.

models offer the same accuracy, one is preferred over the other if it contains fewer parameters, according to the parsimony principle, i.e.: Occam’s razor.29 The balance between precision and the complexity of the AIC results from the fact that the maximized log-likelihood, log(3(θ|y), is a biased estimator of expected relative Kullback−Leibler (K-L) divergence.30 The latter is asymptotically corrected by the number of parameters in the model (i.e. K = |θ| = m).

Problem LS2. min F(θ) = [s(F1(θ)), s(F2(θ)), ..., s(Fq(θ))] θ

s. t. Aθ ≤ b g (θ) ≤ 0 θ∈P

AIC = − 2log(3(θ|̂ y) + 2K

(7)

2.3. Mixed Integer Nonlinear Programming (MINLP) Approach for Parameter Selection. Although problem LS2 is sufficient to fulfill goal (a) in section 1, the best functional selection is required to fulfill goal (b) of said section, to achieve model development. To do this, the parameter space is extended to include a binary matrix a, which determines the selection of a particular model parameter. Therefore, the previous least-squares method is transformed from a nonlinearprogramming problem (NLP) to a mixed-integer-nonlinearprogramming (MINLP) one. Problem LS3.

For least-squares with Gaussian distribution of elements of matrix F(θ) (i.e., estimation residuals), eq 10 adopts the following expression: AIC = N (log(2π /N ) + log F(θ)

θ

s. t. Aθ ≤ b g (θ) ≤ 0

AICc = AIC +

θ∈P (8)

+ 1) + 2K

(11)

2K (K + 1) N−K−1

(12)

The implementation of AICc is modified in this work by including knowledge of the domain field to prevent the models from being built without physical meaning. This modification is shown in section 3. 2.4. Solution to Problems of the Form LS4. There are several alternatives to address LS4 or LS2 problems.31 A coarsegrain classification of available methods consists of two groups depending on the strategies used to obtain the efficient-front: (1) Stepwise construction of the Pareto-front transforming the original problem into several mono-objective subproblems; (2) Use of a multiobjective procedure. The ε-constraint algorithm,32 which belongs to the first group, is used in this work in its essential implementation, since it is a solverindependent methodology and is appropriate for generating spread-efficient convex and concave fronts. 2.4.1. The ε-Constraint Method. In the ε-constraint method, all the objectives but one become inequality constraints, setting limits to their εI values. The one that remains is used as the true objective. The procedure starts by determining extreme (best and nadir) values for each objective. Thereafter, a grid of values is constructed of εi over given ranges and one mono-objective problem is solved for each of its points. An advantage over other methods is that the diversity of ef ficient results is easily controlled by adjusting the size of the grid. The problem of LS4 is rewritten by setting the true objective function as the AIC metric of the model for one of the properties. The remaining squared-errors are converted into constraints that are added to the existing Big-M constraints.

The selection of variables in the LS3 problems is handled by manipulating the constraints, adapting it well for deployment for the final problem. The domain of each parameter in the global model is altered by the corresponding element aj, from matrix a. When a given element of a is set to zero, the associated parameters are excluded from the model; otherwise, they are included but limited to their upper and lower values. The combinatorial nature of this formulation is often addressed by relaxing integer or binary constraints, letting the program explore the search space continuously. In each iteration step, the result obtained is made to satisfy integer or binary constraints by applying strategies such as refinement or branching.24−26 In this work, the relaxation problem is treated using the Big-M formulation.27 However, in this approximation the disjunction impositions are eliminated, making it no longer necessary to fulfill complementarity constraints, i.e., ∑jaj = 1. Problem LS4: The Big-M method. min F(θ) θ

s. t. g(a , θ) ≤ 0 θ∈P −M1a1 ≤ θ1 ≤ M1a1 ∂ −M mam ≤ θm ≤ M mam a ∈ {0, 1}m

2

In the proposed methodology, the RMSE is replaced by the AIC to drive the search for the best model. If the quadratic errors of the model, or the number of parameters, decrease in eq 11, then the AIC also decreases; that is, for two models with the same accuracy, the simplest is preferred. Therefore, one objective is to find the parameters that minimize the AIC. Considering the finite size of the available data sets, in practice a modified version of the AIC is applied instead,28 noted as

min F(θ , a)

a ∈ {0, 1}m

(10)

(9)

2.3.1. Information Driven Search. The Akaike Information Criterion (AIC)22,28 is used for the model selection. When two 12878


Correlation

Industrial & Engineering Chemistry Research Problem LS5.

to have real values between their limits. This problem is then solved by an appropriate NLP solver. The global problem is partitioned into nodes that are arranged in a tree-like structure, each of which contains a new relaxed subproblem. When the result of the relaxed subproblem verifies the integer constraints, it is not further partitioned, and the result is stored to compare it with those obtained in other nodes. If a certain subproblem is not feasible, the node is not further divided, since any subproblem that has emerged from that node is unfeasible. The SBB method is implemented here in the modeling language GAMS (General Algebraic Modeling Language).35 The main parameters of the method are summarized in Table 1; the remaining ones maintain their default values.

min AICc (F1(θ)) θ

s. t. s(F2(θ)) ≤ ε2 ∂ s(Fq(θ)) ≤ εq −M1a1 ≤ θ1 ≤ M1a1 ∂ −M mam ≤ θm ≤ M mam a ∈ {0, 1}m

(13)

Table 1. Branch-and-Bound and Genetic Algorithm Settings

2.4.2. Results Set Refinement. When the optimization process ends, the set of results obtained (each of which corresponds to a given model) contains some dominated members depending upon the capability of the optimization solver used. Hence, some of the results coming from the procedure must be filtered out to finally have an efficient-set that guarantees the real efficiency of each result considered in the decision-making task, as defined in Appendix A (Def 3). 2.4.3. Quality Evaluation of Fronts. To study the quality of the achieved ef f icient-f ronts, the / -hypervolume metric, proposed by Zitzler and Thiele,33 is used. The formalism of the metric is described in terms of the Lebesgue measure, λ(EF),34 on the ef f icient-front. /(EF *, yref ) = λ( ∪ {y′|y ≺ y′ ≺ yref }) y∈E*

Parameters

Value

Branch-and bound algorithm NLP Solver CONOPT Integer tolerance 1 × 10−5 Multistart repetitions 1000 Maximum number of nodes ∞ Genetic algorithm Max. generations 2000 Population size 150  → Real/B → binary Coding Selection operator Linear scaling roulette Crossover operator SBX → Real/BUX → binary Crossover rate 0.95 Mutation operator UM → Real/BF → binary Mutation rate 0.25 Elite fraction 0.25 Number of generations without improvement 300 (early stopping)

(14)

In practice, the / -hypervolume metric gives the space covered by the nondominated results in relation to an arbitrarily selected reference point (yref), which must be dominated by every result in the ef f icient-front. The measure rewards both the departure of the reference point from the ef ficient-f ront (related to convergence) and the distribution of the results in the front (diversity). To check this last criterion, in this work the metric used is such that the reference point is selected as the nadir point of the ef f icient-f ronts achieved after scaling and normalizing the objective space. 2.5. Algorithm Notes. The resolution of each step of the ε-constraint algorithm requires an optimization method to obtain the solution of an LS5 problem. Previous works17 have shown that even the algorithms optimizing the parameterfitting problem (LS2) have difficulty obtaining the global optimum, due to the strong influence of the starting point. This means that several runs, with different initial values, must be carried out to ensure that the parametrization obtained is probably among the best possible. In this work we compare the performance of deterministic strategies20 and evolutionary algorithms21 to address thermodynamic model selection. However, this latter approach for parameter-fitting is not always justified because of the additional computational cost. As the influence of initial guesses becomes evident, it is possible to exploit the implicit parallelism of evolutionary heuristics to overcome that obstacle. 2.5.1. The Simple Branch-and-Bound Algorithm (SBB). The SBB algorithm is a method for Mixed Integer Programming (MIP) problems that, when combined with an NLP method, enables the optimization of LS3 problems. The algorithm (Figure S1 of Supporting Information) starts by relaxing the search space by allowing integer/binary variables

2.5.2. Real-Coded Genetic Algorithm (GA). A GA proposed by Deep36 is chosen as a rough basis to optimize the LS5 problem. This concept is suitable to deal with MINLP problems, using a phenotypic representation of the chromosomes that match with genotypic coding. Different operators have been used for the recombination and mutations of real (θ) and binary (a) parameters: The Simulated Binary Crossover (SBX)37 and the Uniform Mutation (UM)38 are used for real genes (variables), while the Binary Uniform Crossover (BUX) and Flip Mutation (BF) are considered for binary ones. Constraints are handled by a parameter-free penalty function as suggested by Deb.39 The iterative process depicted in Figure 1 was implemented in MATLAB while the parameters of the method are recorded in Table 1. The constraints related to the Big-M method (see LS3) were reformulated to increase the efficiency of the evolutionary algorithm. Due to the nature of the genetic operators, the BigM constraints will cause the crossing and mutation operators to produce new candidates that violate the constraints, since the real and binary variables are modified independently. Thus, the evolutionary algorithm is forced to continuously modify or penalize such individuals, thus increasing the randomness of the search. To overcome this problem, the genetic operators are adjusted so that the real variables are made dependent on the binary ones. Hence, binary variables are first crossed and/ or mutated giving new forms of the candidate model. After this, the search space for each real variable is extended or 12879


Correlation

Industrial & Engineering Chemistry Research n

Zn(T , p , x) =

∏ zi(T , p , x)

(18)

i=1

Hence, the particularization of the model, defined as eq 17, for a binary solution, which is the most common case and that we will use in the current work, is written as g2,EQ (p , T , x) = z1(T , p , x)z 2(T , p , x) × Q −2

∑ gj(p , T )z1j(T , p , x)

(19)

j=0

The functional gj(p,T) are linear-combinations of the influencecoef f icients generating the polynomic form of the model. In this work, the vE’s are not considered and, as gE varies weakly with pressure, this variable can be removed from eq 16 for the activef raction and following equations. So, the coefficients gj(p,T) are formulated as

Figure 1. Description of genetic algorithm main steps and operators for MINLP problems

g j(T ) = g j1 + g j2T + g j3 /T + g j4 T 2

reduced according to the values of its associated binary variables. According to this strategy, Big-M linear constraints are transformed into dynamic boundary-constraints. All calculations are made in a computer with an Intel processor i7-6700 @ 3.40 GHz and 16 GB of RAM.

(20)

Attending to the thermodynamic relations already presented in eqs 5 and 6, expressions for activity coefficients γi and hE, in the context of the present model, are RT ln γi(p , T , x) = z1z 2

k il(p , T )xi ςi(p , T )xi = ∑j ςj(p , T )xj ∑j k jl(p , T )xj

∑

+ ( −1)

(21)

(22)

where Q −2

Y = (z 2 − z1)

∑ j=0

Q −2

g jz1j + z1z 2

∑ jgjz1j − 1 j=1

(23)

In practice, the interaction of order Q and the number of terms in eq 20 are established, a priori, according to experience. If necessary, they are modified during the modeling process to guarantee a good representation of the properties being modeled. However, this criterion is not rigorous enough and, sometimes, leads to overfitting issues. In addition, the dependencies of functional gj on temperature do not have to be the same for each interaction order. It is important to know the influence of the model form on the modeling and how errors are propagated to the simulation of certain processes. This precedent motivates the automatic configuration of eqs 15, 19, and 20, which is addressed in this work by setting a sufficiently large value of Q to define a metamodel that will be optimized according to the proposed methodology. Section 2.3.1 explains how the model is selected with the AICc approach based on Information Theory. According to the Akaike criterion, if only the number of parameters were counted, all of them would have the same probability of being included, or considered, in the final set of model parameters. In other words, the resulting model would contain some parameters that do not contribute to the physical interpretation of the property under study. For example, if in the model expressed by eq 19 we consider a fourth-order interaction for a binary system, gE2,4(p,T,x), this consideration would integrate the previous orders but may have no physical sense, and those

(16)

1− i 2 − ... −i n − 1) gnE(i (p , T , x ) − 1, Q

i1,i 2 ... i n − 1∈ CR * (n , n − 1)

+ Zn(T , p , x) ·χQ (p , T , x)

i−1

Ä ÑÉÑ Q −2 Å ÅÅ i ∂z y jij ∂gi zyzÑÑÑ i Å E Å h (p , T , x) = z1z 2 ∑ ÅÅg1 − T jj zzÑÑz1 − TY jjj 1 zzz j ∂T zÑÑ ÅÅ k ∂T { k {ÑÖ i=0 Å Ç

in eq 15, ai1,i2,···,il are the influence-coef f icients related to the interaction of l-th order (i.e. involves l molecules), n is the number of components, and Q is the maximum order of interaction considered. Equation 15 may be expressed in a more useful form by expanding the summations and grouping the terms in a recursive expression. gnE, Q (p , T , x) =

g jz1j

j=0

where zi is formulated as z i(x , p , T ) =

∑

ij z yz x3 − ijjj 1 zzz kY jx z k 1{ 2

Q −2

3. THERMODYNAMIC MODEL The combined modeling of the vapor−liquid equilibrium (VLE) data and mixing enthalpies hE is carried out by applying a polynomial model previously used by the authors,40−43 for multiproperty correlations. The model is formulated for the excess Gibbs function, gE, of the dissolutions in terms of the socalled active-f raction, zi, of component i, which takes into account the contribution of each individual component to the mixing property. The general expression for the model, obtained from statistical principles on possible molecular interactions, is Ä ÑÉÑ ÅÅ Q Å ÑÑ ÅÅ E gn , Q (p , T , x) = ∑ ÅÅÅ a i1,i 2 ... ilz i1z i 2 ... z il ÑÑÑÑ ∑ Å ÑÑ l=2 Å ÅÇÅ i1,i 2 ... il ∈ CR * (n , l) ÑÖÑ (15)

(17)

E(i1−i2−...−in−1) gn−1,Q

where is the excess Gibbs energy of all subsystems due to interaction of (n − 1) components up to a maximum of Q molecules and χQ(p, T, x) is a polynomial whose degree depends on Q (in general it will be of degree Q − n). The general form of Zn(T, p, x) expressed in terms of activef ractions, is as follows 12880


Correlation

Industrial & Engineering Chemistry Research

Figure 2. Comparison between the eff icient-f ronts obtained with NLP and MINLP procedures (a) and the resulting gE/RT surfaces: (b) NLP model (hE = 92 J·mol−1); (c) MINLP model (hE = 97 J·mol−1). Both procedures are compared on the basis of the ethyl ethanoate (1) + ethanol system (2).

a bad modeling when the real cause is the poor or fair quality of data. 4.1. Nonideal Azeotropic Systems. The presence of azeotropes produces technical difficulties in the separation processes and significantly influences on their designs. Therefore, it is important to obtain a model able to accurately describe the conditions of that singular point. The binary ethanol + ethyl ethanoate is an example of this type of dissolution. In the literature, 23 publications dealing with the VLE of this system were found;46−68 six of them are iso-T, and the remaining ones are iso-p. VLE data by Perelygin and Suntsov65 were discarded for violating the condition (c) established in section 4; the other data series are depicted in Figure S2. Another 12 papers58,59,69−78 in which hE is measured in the interval [295−413] K are incorporated into the data set. The consistency test (see Table S1) declared 13 VLE data series as inconsistent, and of the others, four were used to define the model and five to test its extrapolation capability. 4.2. Nonideal Zeotropic Systems. As the difference between the boiling points of the components in the mixture increases, the possibility of azeotropic formation decreases; so, many nonideal systems are zeotropic. The binary chloroform + benzene presents a negative deviation from Raoult’s law, contrary to the case of section 4.1 (Figure S3). The Thermolitdatabase compiles 16 references79−94 of VLE data of the aforementioned system, measured under different conditions. Seven data series80,83,89−94 do not meet the requirement (c) and were discarded. The mixing effects for the binary chloroform + benzene are exothermic (hE < 0), and the hE data are presented in six articles.73,95−99 The extrapolation capacity of the models is used to reproduce the enthalpic effects of the mixture and is limited because the literature data are in the same p,T conditions. The results of the consistency test are summarized in Table S2. 4.3. High Relative Volatility Systems. Differences in the molecular interactions and/or in the molecular weights of the pure components in the mixture increase the boiling point interval and the relative volatility of the mixtures. In these cases, the experimental data have a greater uncertainty level, because of the low composition of the nonvolatile compound in the vapor phase, so careful modeling is required to avoid an

of a higher order (i.e., fourth) prevail over those of a lower order. For a simple correlative model, this is not a problem from a thermodynamic point of view, since as the order of the model increases, the number of parameters increases; then, the parameters associated with a greater interactional order would have less weight in the studied property calculation. To avoid this drawback, the K factor of eq 12 is modified according to the following relationship: Ñ |a| iÅ yz jÅÅ i ÑÑÑÑ K = ∑ jjjÅÅÅÅ + 1zzz·ai jÅÅ Q ·v ÑÑÑÑ z (24) Ö i = 1 kÇ { where ν is the number of parameters of the functional gj(T), eq 20, Q is the maximum order of the model, and ai is the indicator variable that establishes if a parameter is included or not. In this way, the indicator variables are weighted as in eq 24, and the model complexity is evaluated by multiplying its value by the order of the influence-coef ficient to which it is associated.

4. DATASET DESCRIPTION The ability of the proposed modeling methodology to represent thermodynamic properties, using empirical equations, is tested for several systems that exhibit interactions of different class and magnitude. Three VLE types are analyzed, which are classified as (a) nonideal azeotropic; (b) nonideal zeotropic; and (c) those of high relative volatility (α12 = y1·x2/ y2·x1). In the Thermolit database from NIST,44 systems suitable for each of the defined classes are selected, provided that sufficient experimental information is available to train and test the models. In this sense, the selected systems should meet the following requirements: (a) There are iso-p VLE data at different pressures, or at least, there are iso-p and iso-T VLE data. (b) There is at least one data series for hE. (c) The information on VLE must be complete, that is, with well-defined liquid and vapor compositions. Besides these requirements, data sets should be checked for their thermodynamic consistency, which is carried out using a rigorous methodology,45 see S1 in Supporting Information, by using the binomial {data+model}. This checking is important because, if it is not done, this can lead to false conclusions, i.e., 12881


Correlation


Table 2. Analysis of Importance of Model Parameters (gij) on gE/RT Description After NLP-simplex/ε-Constraint Runsa

Model with a maximum interaction order of 5 (16 temperature-dependent parameters). Δϵij(gE/RT) = | ϵ1̅∀ij − ϵ0ij| . Error in excess enthalpy increases as row number grows. a

Table 3. Analysis of Importance of Model Parameters (gij) on hE Description after NLP-Simplex/ε-Constraint Runsa

Model with a maximum interaction order of 5 (16 temperature-dependent parameters) Δϵij(hE) = | ϵ1̅∀ij − ϵ0ij|. Error in excess enthalpy increases as row number grows. a

proposed consistency test. Seven references containing hE data117−123 were included in the database for temperatures ranging from 298 to 338 K.

incorrect representation of the thermodynamic behavior. As an example, we will use the VLE data of the binary ethylene glycol (1) + water (2), as shown in Figure S4. In this case, the molar mass of water is almost one-third of that of ethylene glycol, which in addition to the greater number of −OH groups in alcohol, yields a difference in their boiling points of ΔToi ≈ 100 K. Although both compounds have hydrogen-bonding interactions, the cross-association increases the nonideality of the mixture causing the actual excess Gibbs energy to be lower than that of the ideal solution, gE/RT < 0. For these mixtures, the mixing process favors a better molecular aggregation with hE ≈ −800 J·mol−1 at intermediate composition. In the literature we found 18 references100−117 containing VLE data of ethylene glycol + water, however, requirement (c) eliminates ten of them, and the analysis of the remaining ones100,102,105,108,112−114,116 is shown in Table S3. These eight references comprise a total of 16 VLE at different conditions (11 iso-p, 5 iso-T), ten of which were validated by the

5. MODELING RESULTS 5.1. Results for Ethanol (1) + Ethyl Ethanoate (2) System. The results obtained with the methods indicated in the previous sections are discussed below. 5.1.1. Nonlinear vs Mixed Integer Nonlinear Programming (LS2 vs LS5). To demonstrate the usefulness of the proposed methodology in semiempirical models, eq 19 is fitted to the ethanol + ethyl ethanoate data series. With an interaction-order of Q = 5, the resulting models of the optimization process have, at constant temperature, four parameters of composition, while temperature-dependence is given by eq 20. In Figure 2(a), the ef f icient-f ront achieved with the NLP-simplex/ε-constraint algorithm is shown. The final set of results is enclosed within an interval of enthalpy errors, 12882


Correlation


Figure 3. Obtained eff icient-f ronts for the following systems: (a−c) ethanol + ethyl ethanoate, (d) chloroform + benzene, (e) ethylene glycol + water. (black diamonds) solution from the final population; (red diamonds) efficient solution; (blue diamonds) efficient solution selected for model comparison.

s(hE)/J·mol−1 ϵ [10−100]. The multiproperty combinedmodeling of {VLE + hE} set for binary solutions does not usually require a model with Q > 4 [see refs 41, 42], since it would be overparametrized. To confirm this, an importanceanalysis on model parameters was made in order to check the influence of each individual parameter in the description of the properties mentioned. This was done sequentially: cutting out the parameters from each model in the ef f icient-f ront without modifying the other parameters, observing the changes that appear in the absolute average relative deviation, AARD, eq 25, for hE and for gE is computed according to eq 26: ϵ[f (θ)] =

1 N

N

∑ k=1

Δϵij = |ϵ 1∀̅ ij − ϵ0ij|

for this property. The parameters of higher interactional-order are more dispersed, and many of them produce error changes below 10%. Since errors in the hE increase, Table 3, from top to bottom, a greater number of parameters is required to describe both properties simultaneously. The colored maps reflect a clear partial conclusion; however, there are several cases in which the error share of a given parameter is less than 20%, so we can ask ourselves: Would it be possible to amend the error by cutting-out certain parameters and by fine-tuning the others? The hE results for Table 3 show that for small model errors, see Figure 2(a−c), each parameter shares almost the same importance-degree, except for the gi2 coefficients, which do not influence the estimate of hE, as deduced from eq 6. However, as the tolerance in the description of hE is relaxed, many parameters become negligible (with relative error increments below 0.1%). Since these models do not provide much quality in the description of these properties by default; the weak influence of several of the parameters is justified. The results obtained for the parameter-fitting problem (LS2) suggest that the optimization algorithm ignores the unnecessary parameters of the model, although these are not completely cancelled-out, becoming a problem for extrapolation purposes since they introduce noise in the models. Proof of this is observed in Figure 2(b), in which the extrapolation of gE/RT at high temperatures is presented and compared with the results achieved with the MINLP approach,

|yk − f (x k , θ)| yk

(25) (26)

where ϵ0ij is the AARD when parameter j from model i is withdrawn from the latter. Conversely, ϵ1∀ij is the AARD when the complete model is considered. Tables 2 and 3 present two colored maps indicating the magnitude of the change in error. Those schemes show how the influence of several parameters in each individual model is negligible. Table 2 shows that the parameters of the temperature-functional for the second-order interaction parameters (g11 − g14) are, on average, those that produce the greatest error changes for gE/RT, hence the most relevant 12883


Correlation


Table 4. Root Mean Square Errors (RMSE) of VLE and hE of Efficient Models Obtained with the MINLP-SBB/ε-Constraint Approach for the System Ethanol + Ethyl Ethanoate Eff icient results (models) s(F)

1

2

3

4

γ1 γ2 gE/RT yiso‑P T yiso‑T p hE

0.652 0.665 0.192 0.119 6.068 0.100 16.269 11.620

0.469 0.578 0.127 0.079 4.161 0.074 11.401 12.992

0.256 0.144 0.012 0.012 0.359 0.012 1.068 14.795

0.237 0.126 0.010 0.011 0.307 0.009 0.877 22.867


1.099 0.695 0.206 0.145 5.470 0.121 14.167 77.597

0.972 0.623 0.153 0.099 3.843 0.093 9.615 75.145

0.852 0.253 0.101 0.020 0.333 0.037 0.607 85.107

0.849 0.246 0.101 0.021 0.341 0.038 0.570 76.525

5

6

Training Set 0.239 0.240 0.138 0.108 0.011 0.010 0.011 0.010 0.336 0.302 0.011 0.010 1.041 0.953 28.033 30.757 Checking Set 0.849 0.849 0.251 0.234 0.101 0.101 0.021 0.019 0.322 0.310 0.037 0.036 0.576 0.536 89.585 95.955

7

8

9

10

0.256 0.136 0.009 0.012 0.290 0.012 0.720 32.538

0.248 0.085 0.008 0.011 0.277 0.009 0.649 49.225

0.179 0.105 0.008 0.010 0.229 0.012 0.810 63.334

0.227 0.100 0.008 0.010 0.236 0.011 0.802 111.631

0.852 0.249 0.101 0.021 0.327 0.039 0.761 79.567

0.853 0.226 0.102 0.020 0.365 0.039 0.679 110.196

0.830 0.233 0.102 0.021 0.267 0.037 0.569 96.456

0.840 0.232 0.101 0.020 0.275 0.036 0.587 165.227

Table 5. Root Mean Square Errors (RMSE) of VLE and hE of Efficient Models Obtained with the MINLP-GA/ε-Constraint Approach for the System Ethanol + Ethyl Ethanoate Eff icient results (models) s(F)

1

2

3

4


0.290 0.568 0.071 0.065 1.994 0.078 4.401 18.110

0.404 0.213 0.019 0.030 0.682 0.029 1.583 18.996

0.227 0.076 0.007 0.010 0.212 0.010 0.392 20.672

0.197 0.060 0.006 0.009 0.177 0.008 0.454 38.800


0.825 0.632 0.118 0.064 1.932 0.091 2.975 85.401

0.896 0.302 0.101 0.032 0.601 0.051 1.972 81.784

0.842 0.223 0.102 0.020 0.315 0.036 0.764 88.367

0.835 0.220 0.102 0.019 0.295 0.035 0.608 78.983

5

6


7

8

9

10

0.193 0.056 0.006 0.009 0.177 0.008 0.470 53.197

0.182 0.053 0.006 0.008 0.175 0.008 0.475 67.982

0.176 0.062 0.006 0.009 0.174 0.009 0.505 73.000

0.184 0.054 0.006 0.009 0.177 0.008 0.495 95.970

0.834 0.220 0.102 0.019 0.291 0.035 0.591 76.673

0.832 0.221 0.102 0.019 0.301 0.035 0.570 98.287

0.830 0.222 0.102 0.019 0.294 0.035 0.589 98.306

0.833 0.222 0.102 0.019 0.304 0.035 0.566 124.046

under such conditions. However, the surface generated by the model with the NLP procedure, Figure 2(b), shows an unnatural fold in the rich ester as the temperature departs from the conditions at which the model was obtained. As a result gE/ RT becomes wavy at high temperatures, an unusual behavior for a binary system. These two effects indicate an overparametrization in eq 19. In contrast, transitions in the gE/RT surface obtained with the model generated with the MINLP procedure are smoother: the fold at low ethyl ethanoate composition disappears, and the variation with composition is more regular at high temperatures. 5.1.2. Comparison of Solvers. The training samples for the ethanol + ethyl ethanoate system were used to build some models with LS5 strategies: one with the GA/ε-constraint, Figure 3(a), and the other with the SBB/ε-constraint

Figure 2(c). Consequently, it is convenient to formulate the problem of obtaining the model in the form of LS5. Additionally, the MINLP procedure for LS5 shows a trend for the zero-valued parameters which generally coincides with the results of the LS2 problem. Parameters of the functional temperature, eq 20, related to second-order interactions in eq 19 (i.e.: g11 − g14), have nonzero values in most cases. Furthermore, the models become denser as the errors in the hE become lower, which is also observed for the NLP results. In contrast to the models achieved with the NLP-Simplex/εconstraint method, the nonimportant parameters are cut-out from the final model, preventing problems in the extrapolation. Since there are no experimental values at higher temperatures than those represented, it is not possible to decide whether a certain model better describes the behavior of the dissolution 12884


Correlation


Figure 4. Comparison of efficient model performances for the system ethanol (1) + ethyl ethanoate (2) over the training set (a−c): (△) ref 46, (○) ref,51 (◇) ref 55, (dotted circles) ref 63, (∗) ref 69, (▼) ref 75, (●) ref 59, (■) ref 77; and checking set (d−f): (red plus signs) ref 57, (dotted diamonds) ref 70, (∗) ref 71, (slashed circles) ref 72, (★) ref 73, (▲) ref 74; (a, d) y1 − x1 vs x1; (b, e) T vs y1, x1; (c, f) hE. Model 1 (blue lines), Model 5 (green lines), Model 6 (black lines).

constraint requires a smaller step for ε to obtain a high coverage of the ef ficient-f ront. In practice, this involved using 1000 optimization steps for the deterministic approach, in contrast to the 50 steps used for the GA/ε-constraint algorithm. From all the runs, 13.7 ± 1.57 ef f icient-results were found with the SBB/ε-constraint and 12.3 ± 1.25 with the GA/ε-constraint. Implicit parallelism of the evolutionary method favors convergence toward the ef ficient- front, so a significant reduction in the number of runs is achieved using the latter procedure to obtain a f ront formed by a similar number of results. Given the somewhat better quality of the ef f icient-f ront produced by the GA and its higher efficiency, we relied on this approach for the two remaining cases. 5.1.3. GA/ε-Constraint Efficient-Front Evaluation. Threeout-of-12 models were extracted from the ef f icient-f ront obtained with the GA, light-blue in Figure 3(a), and were compared in Figure 4. Results with unreasonable errors in both gE/RT and hE were not considered. Since the shape of the ef f icient-f ront is very steep at low AIC values, all models describe well the VLEs, see Table S7, even for the checking samples. The T − x, y diagram shows that the checking data series at 77 kPa is the worst estimated, Figure 4(e). However, since the remaining VLEs in the checking set are well described and their conditions are more extreme, it is possible that the equilibrium at 77 kPa is not coherent with the whole data set, which cannot be detected by the consistency test used for raw data evaluation. In practice, the main difference among these models is the description of the hEs. Even though the obvious choice would be to keep the set of parameters that

algorithm, Figure 3(b). Tables S4 and S5 show that the models achieved with the first procedure have fewer parameters than those achieved with the second one. Since errors in gE/RT are comparable, Tables 4 and 5, it follows that values for the AIC are better (less complex) for those models coming from the evolutionary-based methodology. Nevertheless, the SBB algorithm can produce results with lower hE errors (Table 4) than the minimum value observed with the GA. Hence, none of the ef f icient-fronts obtained with both algorithms completely dominates over the other, see Figure 3(c). Note that the differences in the representation of the hEs are not significant for models obtained from any solver as seen for two trade-off models (GA: Table S5, model 3; SBB: Table S4, model 3) depicted in Figure S5. Therefore, the greater error in hE for the GA model does not worsen the qualitative and quantitative description of the property compared to that of the SBB algorithm. Regarding the quality of the ef ficient-f ronts, those produced with the GA have a higher hypervolume metric, see Table S7, compared with that of the SBB\ε-constraint; the first approach also shows smaller standard deviation (0.008 vs 0.030) of the convergence metric, demonstrating the robustness of the solver. Differences in the hypervolume values are more related with the shape of the ef f icient-f ront than with the actual diversity of the results. This is caused by the size of the tradeoff region obtained with both algorithms, being sharper for the GA. Besides, the efficiency of the procedures, quantified as the number of ef f icient-results over the total number of screened εvalues, is the most significant existing difference. The SBB/ε12885


Correlation


Table 6. Root Mean Square Errors (RMSE) of VLE and hE of Efficient Models Obtained with the MINLP-GA/ε-Constraint Approach for the System Chloroform + Benzene Eff icient results (models) s(F)

1

2

3

4


8.889 7.350 0.686 0.390 30.626 0.395 72.720 10.333

0.021 0.061 0.013 0.005 0.472 0.006 0.382 11.800

0.021 0.062 0.013 0.005 0.466 0.006 0.380 15.400

0.020 0.065 0.013 0.006 0.458 0.005 0.379 44.177


10.101 9.016 0.673 0.396 32.189 0.381 111.352 30.975

0.023 0.031 0.010 0.005 0.269 0.006 0.573 35.248

0.023 0.031 0.010 0.005 0.267 0.005 0.591 37.104

0.025 0.032 0.010 0.005 0.253 0.004 0.633 57.332

5

6


7

8

9

10

0.020 0.062 0.013 0.005 0.459 0.005 0.343 76.600

0.020 0.062 0.013 0.005 0.458 0.005 0.345 78.400

0.020 0.062 0.013 0.005 0.458 0.005 0.346 80.200

0.019 0.064 0.013 0.005 0.456 0.005 0.358 91.000

0.027 0.027 0.007 0.005 0.210 0.005 0.419 87.603

0.027 0.027 0.008 0.005 0.210 0.005 0.437 89.486

0.027 0.028 0.008 0.005 0.210 0.005 0.455 91.382

0.027 0.029 0.009 0.005 0.215 0.005 0.564 102.776

Figure 5. Efficient model performance comparison for the system chloroform (1) + benzene (2) over the training set (a−c): (▷) ref 86, (◁) ref 87, (red circles) ref 88, (red squares) ref 89 (stars) ref 73, (○) ref 97, (□) ref 98; and checking set (d−f): (△) ref 81, (▽) ref 82, (+) ref 95, (open plus signs) ref 96, (×) ref 99; (a, d) y1 − x1 vs x1; (b, e) T vs y1, x1; (c, f) hE. Model 2 (blue lines), Model 5 (green lines), Model 8 (black lines).

constraint steps after applying the globality filter. The strict implementation of the filtering process eliminates several results from the final population, producing a gap in the front, see Figure 3(d). Table S6 summarizes the parameters for ten selected models of the 31 that proved to be efficient. As for the ethanol + ethyl ethanoate system, the models that yielded the most precise hE values contain a greater number of parameters. However, the first model in the table is useless because of its

produces the lowest error for the hE, a thorough analysis of error propagations to the simulation should be done, as indicated by the authors.17 However, the objective of the present work is not to provide a rigorous exploration of the ef ficient-f ront. 5.2. Results for System Chloroform (1) + Benzene (2). The efficient-front achieved with the MINLP-GA/ε-constraint approach for this system consists of 31 results, with 51 ε12886


Correlation


Figure 6. Efficient model performance comparison for the system ethylene glycol + water over the training set (a−c): (red triangles) ref 100, (▷) ref 105, (□) ref 116, (+) ref 119, (□) ref 117; and checking set (d−f): (▽) ref 102, (red left pointing triangles) ref 108, (○) ref 113, (◊) ref 118, (open plus signs) ref 120, (stars) ref 121, (○) ref 122, (×) ref 123; (a, d) y1 − x1 vs x1; (b, e) T vs y1, x1; (c, f) hE. Model 1 (blue lines), Model 5 (green lines), Model 9 (black lines).

AIC value. Removing this model from the efficient set leaves a quasi-vertical front composed of a wide variety of models in terms of hE performance but with close errors for the VLE data set (gE/RT and γi’s), see Table 6. Figure 5 shows three examples with this situation. Two temperature-dependent parameters and an interactionorder of Q = 2 in eq 19 are sufficient to describe both thermodynamic properties VLE (gE/RT) and hE. As constraints for hE are relaxed, VLE can be described with only one parameter. The activity coefficients of the two components in the solution are similar in magnitude and in the variation with composition, giving rise to a symmetric gE /RT. The representation of hE is slightly asymmetric toward pure benzene, justifying that the values of the parameter k21 of the best models for that property are somewhat greater than one. On the contrary, as the error tolerance for hE increases, k21

adopts values of around one, caused by the symmetry already mentioned for gE/RT. 5.3. Results for the System Ethylene Glycol (1) + Water (2). From the three cases studied, the efficient-front found in Figure 3(e) is the one with the largest trade-off region, ranging from 29.8 to 40 for s(hE), and with AIC values ranging from −13 to −31. On both sides of this region, the front extends vertically or horizontally. In other words, one criterion remains almost constant while the other changes. Thus, two results of the trade-off region were selected (models 1 and 5 in Table S7) along with another one with a large hE error (model 9) that is plotted in Figure 6 (highlighted in light blue). Model 1 outperforms models 5 and 9 in the description of hE, whose errors are listed in Table 7. However, as shown in Figure 6(i), even this model is not capable of describing the variation in this last property with temperature according to 12887


Correlation


Table 7. Root Mean Square Errors (RMSE) of VLE and hE of Efficient Models Obtained with the MINLP-GA/ε-Constraint Approach for the System Ethylene Glycol + Water Eff icient results (models) s(F)

1

2

3

4


0.141 0.080 0.068 0.020 1.689 0.004 1.720 29.800

0.113 0.073 0.055 0.022 0.762 0.004 0.710 31.600

0.139 0.139 0.048 0.046 2.629 0.004 1.078 33.397

0.107 0.068 0.056 0.020 0.632 0.004 0.710 35.200


0.234 0.070 0.085 0.007 2.183 0.008 3.148 31.009

0.210 0.068 0.079 0.005 1.603 0.011 3.962 35.353

0.230 0.175 0.086 0.028 2.776 0.042 7.443 34.391

0.213 0.067 0.079 0.005 1.598 0.010 3.878 35.532

5

6


the checking set. This data series119 exhibits a minimum close to 308 K, which does not agree with the series used for training the model, which shows a monotone growing with temperature. Also, the two references used as training data series117,119 do not match, so the models obtained only reflect an average behavior of the whole data set considered for this system. Model 9 has the best AIC value, presenting significant differences with model 2 when iso-p VLE cases are described. For models 5 and 9 the description of the VLE could be considered equivalent without further analysis, but model 5 gives a better description of the average behavior of the hE. These results reinforce the comments for the two previous cases, for which it was observed that more parameters are necessary to represent the hE fitting. As constraints on this property are relaxed, the use of more than two parameters for the present case is redundant and might lead to the problems already mentioned in the ethanol + ethyl ethanoate system.

7

8

9

10

0.109 0.070 0.055 0.021 0.658 0.004 0.716 65.800

0.109 0.071 0.055 0.021 0.668 0.004 0.722 71.200

0.109 0.072 0.055 0.021 0.689 0.004 0.720 85.600

0.109 0.072 0.055 0.022 0.703 0.004 0.715 96.400

0.215 0.069 0.079 0.006 1.611 0.011 3.810 65.228

0.215 0.070 0.079 0.006 1.624 0.012 3.803 69.458

0.215 0.071 0.079 0.006 1.622 0.012 3.790 84.790

0.215 0.071 0.079 0.006 1.610 0.012 3.782 97.080

and scenarios in terms of VLEs (i.e., azeotropy, high volatility ratios) and hE’s (endothermic and exothermic). Data quality was previously checked by a consistency method.45 A total of 37 data series were validated, and they were used in the subsequent modeling stage, both as training or checking sets. With the modeling of the azeotropic ethanol + ethyl ethanoate system, two optimization procedures were tested: a deterministic one, based on the SBB algorithm coupled with the CONOPT solver, and an evolutionary algorithm tailored to the current application. With both of these, the ε-constraint strategy was used to generate the ef f icient-f ronts. The two methodologies produced well-distributed fronts according to the hypervolume metric, but the evolutionary procedure gave better results. However, models from the MINLP-SBB/εconstraint method require more parameters and show worse AIC values than those obtained with the MINLP-GA/εconstraint for this system. The smallest errors for hE are attained with the deterministic algorithm. For the azeotropic ethanol + ethyl ethanoate system, the evolutionary approach outperformed the deterministic one, since the former produces ef f icient-results in (24.6 ± 2.5)% of the total number of ε-values selected in comparison with the determinist one,

The Parametrization Problem in the Modeling of the Thermodynamic

Recommend Documents