Ind. Eng. Chem. Process Des. Dev. 1985, 24, 494-498
494
that the.last PP rows of A are linearly independent. Thus, there exists an infinite number of solutions to (A4) which satisfy the last r.P equations. By carrying out a Gaussian elimination on the top and permuting the elements of U(L) if necessary, we obtain -
0
I 1
I
A HTo+i
I
. . . H,
U ( Z )+
Hi
-HTo+p . . . . . . . . . .
Q V ( Z - 1) (A5)
where the nonzero last rows of the ney matrix h’ are linearly independent. All solutions U ( k )which satisfy the nontrivial bottom half Of the eq A5 produce the decoupled response y i ( k ) = Si k 2
~ i +
+1
But this is the response of the perfect controller, and since
it is unique, it must be true that all these solutions are equal to the perfect control law.
Literature Cited Foss, A. S.; Edmunds. J. M.; Kouvaritakis, B. Mol. Eng. Chem. Fundam. 1980, 19, 109. Gantmacher, F. R. “The Theory of Matrices”; Chelsea: New York, 1959; voi. 1. Garcia, C. E. Ph.D. Thesis, University of Wisconsin-Madison, 1982. Garcla, C. E.; Morari, M. Ind. Eng. Chem. Process Des. D e v . 1982, 2 1 , 308. Hoit, 6.; Morari, M. Chem. Eng. Sci., in press. Kucera, V. “Dlscrete Linear Control”; Wiley: New York, 1979. Kwakernaak. H.: Sivan. R. “Linear ODtimai Control Systems”; Wiley-Interscience: New York, 1972. Penrose, R. proc, Cambridge phi/os, sot. 1955, 51, 406, ReM, J. G.; Mehra, R. K.; Kirkwood, E. Proc. IEEE Conf. Dec. Control 1979, 307. Wood, R. K,; Berry, M, W. Chem, Eng. sei. 1973, 28,1707,
Received for review April 25,1983 Revised manuscript received May 23, 1984 Accepted July 23, 1984
COMMUNICATIONS Improved Algorithm for Obtaining Liquid-Liquid EquMibrium Parameters by Use of a Penalty Term An algorithm is proposed to improve optimization on liquid-liquid equilibrium data to obtain liquid solution model parameters. The method increases the efflciency of the Marquardt technique coupled with a penalty function. The algorithm was tested with the UNIQUAC and NRTL liquid solution models.
Introduction Due to the increased interest in extraction and multiple-phase distillation, numerous papers concerning liquid-liquid equilibria have appeared in the recent literature. Several molecular thermodynamic models, for example NRTL (Renon and Prausnitz, 1968), UNIQUAC (Abrams and Prausnitz, 1975), and UNIFAC (Fredenslund, Jones, and Prausnitz, 1975), can be used for LLE calculations. The predictive ability of these semiempirical representations of the liquid state is strongly dependent on the method used to determine the model parameters from experimental data. Simonetty et al. (1982) have shown that regressing LLE tie line data alone yields the best model parameter set for type I liquid-liquid systems. LLE correlations with these models are quite nonlinear because of the exponential interaction energy term. Regression on experimental tie line data with a standard numerical technique will yield parameters that give a good fit to the experimental data included in the regression set. However, in subsequent calculations unstable liquid phases may be obtained. Ma,ny techniques to avoid this have been developed. Serrensen (1980) adds a penalty term to the objective fuhction to prevent the model parameters from assuming
large numerical values. This work proposes an algorithm to automatically adjust the penalty function to improve optimization of the model parameters.
Criterion of Equilibrium Three conditions must be met for an LLE system (Sarensen et al., 1979). First, the material balance must be preserved. Second, the chemical potentials of each compound in the two liquid phases must be equal. Finally, the Gibbs energy must be at a minimum at the system temperature and pressure. The last restriction is a necessary and sufficient condition, whereas the second restriction is only necessary but not sufficient. The first and second restrictions are used due to their simplicity and are normally written as
EX;’= 1 i
(1)
Since these are necessary but not sufficient conditions, computations based on eq 1 through 3 can lead to false
0 1985 American Chemical Society Q196-43Q5/85/1124-Q4945Q1.5Q/Q
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985
495
Eval. Eq. (11)
0
I
Local
No
Exceed
Check Meet
Pass
Accept P ( r + l ) P R M RESULT
+
I
Figure 1. Schematic diagram of the proposed algorithm.
solutions (Heidemann and Mandhane, 1973;De Fre and Verhoeye, 1976). Unstable phases such as a local maxima can be detected with a determinant stability test (Smensen, 1980). This test does not rule out the possibility of a false solution due to a local minima.
Objective Functions for LLE The current molecular thermodynamic models contain some imperfections, and the data are not without error, so parameters to do LLE predictions must be determined
with care. Two objective functions commonly used to regress LLE data are
Fa=C C [ aik' k i F, = C min k
+ aiklI
c C ( x i j k
J i
]
2
ai: -
- xijk)2
(4)
(5)
Equation 5 is considered the more attractive objective function because it fits experimental compositions, X,
496
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985
better. However, the latter method requires the construction of a binodal curve and the use of an interpolation procedure to determine which mole fractions on the binodal curve should be compared to the experimental data (Smensen, 1980). It is more convenient to start the correlation with eq 4 and then switch to the F, objective function. Of course, use of the F, objective function does not take into account other considerations such as regeneration of the plait point or solute distribution ratio. Minimization of the objective functions of eq 4 and 5 allows the model parameters to assume any arbitrary value. To prevent the model parameters from assuming large numerical values, Sarensen (1980) adds a penalty function to each of the objective functions in eq 4 and 5 to obtain the following objective functions
+ QaCPn2
(6)
G, = F, + Q,CP,2
(7)
Ga = Fa
The penalty term, Q, is assigned a numerical value based on experience. Unfortunately, it may be different for some LLE systems. In this work, the Marquardt (1963) method is used to minimize the objective functions according to the familiar equation (A(')+ X ( r ) n s ( r ) = g") (8) Addition of the penalty function transforms eq 8 into the form [A(')+
(X'd
+ Q ( r ) ) q g ( r ) = g ( r ) + Q(r)p(r)
(9)
The detailed mathematical expressions of each term in eq 9 are available from Yu (1983). If the value of Q is arbitrarily assigned, then one of the following four cases may result (A(') +
A(r)n$r)
= g(d;
(Q=
+ (A(') + A(')ns(') = g(') + Q(')p(r);
0)
(A(') Q(r)I)8(r) = g(') + Q @ ) P r ) ; ( Q >> A)
A(r)a(r)= g");
(10) (11)
(A >> Q) (12)
(A = 0, Q = 0)
(13)
When Q is very small, the standard Marquardt form is recovered and the model parameters may assume large values. If Q is too large, the penalty function controls the regression and the parameters may be forced to numerical values of zero. This is the case shown in eq 11. For the situation depicted by eq 12, Q has little influence on the increment or step size g"). This factor needs to be considered in constructing en algorithm to avoid a slow convergence rate or a divergence. In the last case, the mathematical formula behaves like the Gauss-Newton method. For the penalty term to properly direct the optimization, the numerical value of Q must be adjusted according to the value of the Marquardt parameter X and the correlation rate. There is an implicit assumption involved in using a penalty function to get LLE parameters. This is that large numerical values of the model parameters tend to give excessively complicated Gibbs energy curves which can lead to a false solution when the regression converges on a local minimum or local maximum. It is too general a statement to say that the minimum with the smallest set of parameters is tHe global minimum. However, in the majority of the cases examined in this study, the minimum with the smallest set of parameters had a lower objective
function value than any of the other local minima with a larger set of parameter values. Therefore, the basis of this work is that the best set of model parameters is obtained when the objective function, F,, is at its lowest value with the smallest possible set of parameters. No further attempt to justify this assumption is made in this work. Penalty Term Adjustment Algorithm Figure 1 is the proposed algorithm to continuously adjust the penalty term. The very fast method of Michelsen (1980) to construct the binodal curve is used to check the new model parameter matrix P1). The regression is allowed to proceed until the first binodal curve is produced. Any subsequent parameter set which gives a predicted binodal curve that either fails to pass the determinant stability test or is not a type I system will be rejected and either Q or X is increased to change both the increment size and direction. After the binodal curve test is executed, the improvement rates, RT1 and RT2, are computed. RT1 is the actual improvement in the objective function G divided by the theoretical improvement corresponding to a "parabolic" surface of the value G as a function of the parameters (Sarensen, 1980). RT2 is the relative improvement rate of the objective function F. The convergence criterion is that the relative improvement rate of F (RT2) must be less than 1% for a few consecutive iterations (the suggested number here is 7) in order to assure that both the penalty term Q and the Marquardt parameter X are well adjusted before the program is terminated. After a few computing iterations with the same penalty term Q, the correlation improvement rate slows down not only because the objective function value is closer to a minimum point, but also because the penalty term is large enough to make the parameter increment keep getting smaller. As the penalty term is reduced by a factor (the suggested value here is 3), the correlation rate speeds up again. The Marquardt parameter is adjusted according to the improvement rate RT1 (Sarensen, 1980). The new set of model parameters is accepted if either the G or F objective function is smaller. The regression proceeds with the new set of model parameters until the iteration limit is exceeded. This iteration limit was given a very high value so that regression continued until the convergence criterion was satisfied. This algorithm may be extended to systems of more than three components and to other nonlinear problems which have to be correlated with specific constraints. The major difficulty of correlating a system with more than three components will be the computing speed in constructing an LLE curve and the complete test mechanism to assure a satisfactory LLE curve. Results The UNIQUAC and NRTL liquid models were used to test this algorithm with nine ternary type I systems from the DECHEMA collection (Sarensen and Arlt, 1980). For calculations with the NRTL model, the initial values of the penalty term for all LLE systems were lo4 for Q, and lo4 for Q,. The UNIQUAC model used the starting penalty terms of for Q, and lo-' for Q, for all but 1 of the LLE systems. All of the optimization runs yield final penalty terms of 10-lo to when regressing on compositions. This penalty term is small enough so that the objective function G, is equal to the desired objective function F,. The results of this work are compared to nine of the systems with specific parameters listed by Sarensen and Arlt (1980). Each optimization with an adjustable penalty had a smaller mole fraction deviation than the correlation
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985
Table I. Experimental Systems in this Work" no. 1 2 3 4 5 6 7 8 9
left component water water water hexane pentane heptane l-hexanol l-heptanol l-octanol
"IB: f,'
upper component ethanol ethanol 2-propanol 2-propanone 2-propanone 2-propanone 2-propanone 2-propanone 2-propanone
right component heptane hexane heptane water water water water water water
> 0.01; 0.99 > f,'I > 0.86; IC:
2:
SYMBOLS.
type IC IC IC IC
IC IC IB IB IB
temp, "C 30 20 25 25 20 25 30 30 30
= (F,/6M)1/2 X 100%; M = number of experimental tie
lines.
Table 111. Individual Mean Deviation (mol %) between Experimental and Calculated Concentration for Nine Type I Ternary LLE Test Systems system
UN IQUA C fixed adjustable penalty penalty 2.48 1.56 0.92 0.32 0.82 0.64 0.50 0.53 0.58 1.81 0.59 0.84 0.25 0.30 0.26 0.49 0.44 0.55
Llne~
Left Component: Centen.
0 ?hied Penalty MOlhod 0 Adluoleble Cenelty Method
A
Upper Component: I-Croponone Ulghl Component: Water
< 0.01; f3u> 0.99.
Table 11. Overall Mean Deviation (mol % ) between Experimental and Calculated Concentration for Nine Type I Ternary LLE Test Systems model correlation method 6s" UNIQUAC with fixed penalty 0.97 UNIQUAC with adjustable penalty 0.57 NRTL with fixed penalty 0.88 NRTL with adjustable penalty 0.58 a$x
SYSTEM 6
Experlmantal Date
- Laperlmentel TI.
497
fixed penalty 2.04 1.02 0.59 0.54 0.69 1.80 0.33 0.40 0.47
NRTL adjustable penalty 1.21 0.52 0.39 0.37 0.69 1.14 0.26 0.24 0.40
with a fixed penalty term. Table I gives the nine test systems. Six of them are type IC systems which are difficult to fit. The overall comparison is summarized in Table I1 which is the mean deviation of concentration for the nine systems. The individual comparison of mean deviation of concentration for each system is shown in Table 111. Figure 2 is an example to show different binodal curves regenerated by the fixed penalty method and the adjustable penalty method. The stars in Figure 2 connected by solid lines are the experimental data points. Starting in the lower left hand corner, a circle and a square are drawn in every 15 calculated points by the plotting program in order to differentiate between the two binodal curves. One of the squares is left off of Figure 2 because it happens to coincide with the plait point. Parameter sets for the UNIQUAC and NRTL models for the nine test systems from the adjustable penalty method are available on microfilm as Tables IV-VII. As is well-known, determination of the parameters for one of the current liquid solution models is often more important than which model is selected for the LLE calculations. The algorithm described in this work was found to provide a stable regression that minimizes the desired objective function while maintaining the numerically smallest possible set of model parameters. The proposed algorithm is not sensitive to the initial penalty term value so a time consuming search for the proper values of Q is eliminated.
Figure 2. Regenerated binodal curves from the UNIQUAC model with the adjustable and fixed penalty method.
Acknowledgment We wish to thank Dr. Aa. Fredenslund of the Instituttet for Kemiteknik, Danmarks Tekniske Hajskole, Lyngby, Denmark, for supplying us with a computer program for calculating LLE. This work was supported by The University of Alabama School of Mines and Energy Development Grant 977R.
Nomenclature A = BTB,B is a matrix of partial derivatives, eq 8 a = yX = activity F = objective function G = objective function with penalty function g = BTf,f is the difference between experimental and calculated values I = unit matrix M = number of experimental tie lines P = model parameter to be correlated Q = penalty term RT1 = improvement rate of G RT2 = improvement rate of F X = liquid phase mole fraction Greek Letters
S = model parameter increment y = activity coefficients X = Marquardt parameter $ = average mean deviation between experimental and calculated values Superscripts A=
experimental value
r = iteration number
I = liquid phase on left of triangular diagram I1 = liquid phase on right of triangular diagram Subscripts a = activity i = component i j = phase j
k = tie line k n = model parameter n x = liquid mole fraction 1 = left component on triangular diagram 2 = upper component on triangular diagram 3 = right component on triangular diagram
Ind. Eng. Chem. Process Des. Dev. 1985, 24, 498-500
498
Literature Cited
Yu, C. Y. MS Thesis, University of Abbama, Tuscaloosa, AL, 1983.
Abrams, D. S.;Prausnitz, J. M. AICh.6 J. 1975, 27, 118. De Fre, R. M.; Verhoeye. L. A. J. Appl. Chem. Blotechnol. 1978, 26, 469. Fredenslund. Aa.; Jones, R. L.; Prausnh, J. M. AIChE J. 1975, 27, 1086. Heidemann, R. A.: Mandhane, J. M. Chem. Eng. Sci. 1973, 28. 1213. Marquardt, D. W. J. Soc.Ind. Appl. Math. 1963, 7 1 . 431. Micheisen, M. L. fluid Phase Equilb. 1980, 4 , 1. Renon, H.: Prausnltz, J. M. AIChE J. 1968, 14, 135. Simonetty, J.: Yee, D.; Tassios, D. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 174. S0rensen. J. M. Ph.D. dissertation, Instituttet for Kemiteknik, Lyngby, Den. . mark, 1980. Ssrensen, J. M.; Magnussen, T.: Rasmussen, P.: Fredensiund. Aa. FlUM Phase Equilib. 1979, 3 47. Sarensen, J. M.; Ark, W. "Liquid-Liquid Equilibrium Data Collection": Chemical Data Series, voi. v., Parts 2 and 3, DECHEMA, Schon & Wekel GmbH: Frankfurt, W. Germany, 1980.
Department of Chemical and Metallurgical Engineering The University of Alabama University, Alabama 35486
Chang-Yu Yu
David W . Arnold*
Received for review December 17, 1982 Revised manuscript received May 21, 1984 Accepted June 4, 1984
Supplementary Material Available: Tables IV, V, VI, VII. NRTL and UNIQUAC parameters from adjustable and fixed penalty method (4pages). Ordering information is given on any current masthead page.
Estimating the Molecular Welght of Petroleum Fractions A set of simple equations is presented derived from an earlier chart method for estimating the molecular weights of petroleum fractions boiling at different temperatures from a blend such as gasoline. The only data needed for the computations are the ASTM distillation curve and the specific gravity of the blend.
Introduction In the course of an investigation on fuellair mixture formation in spark ignition engines, it became necessary to estimate the molecular weight of the fuel vapor as a function of the extent of evaporation of the liquid petroleum fuel. A need of this nature may also arise in other situations of engineering interest. The only available data for the liquid fuel, in most cases, would usually be comprised of the ASTM distillation curve and the specific gravity. With a knowledge of these two alone, it is then desirable to devise a method of estimating the molecular weights of the vapor at different percentages of evaporation of the fuel. A survey of the relevant published literature revealed that the most expedient way of doing this is by using the U.O.P. (Universal Oil Products) characterization factor proposed by Watson and Nelson (1933) and now commonly referred to as the Watson K factor. The molecular weight of petroleum fractions cannot be satisfactorily correlated with any single property. However, it can be correlated with two properties with an accuracy sufficient for most purposes. This is found to be the case with other properties too. Watson and co-workers (1935) developed a series of charts to evaluate various properties of petroleum fractions. Each property is related to two parameters such as,for example, the Watson K factor and the A.P.I. gravity. Watson et al. state that the molecular weight of petroleum fractions can be estimated with errors rarely exceeding 5% with the help of their chart. Winn (1957) later converted the data in these charts to the form of a nomograph in which the values of any two properties determine the values of the rest. The U.O.P. characterization factor or the Watson K factor is defined by the relation (Tb,R)1'3 K=-
S
(1)
where Tb,R = molal average boiling point, OR,and s = specific gravity at 60/60OF. The value of K and the molal average boiling point uniquely determine the average molecular weight of any petroleum fraction and can be determined from the Watson chart or the Winn nomograph. It is stated that the Watson K factor is a measure of "paraffinicity" and has approximately the same value 01 96-43051a51 1i24-0498$01 .solo
for all fractions of a given stock (Hougen and Watson, 1947). The assumption of a constant K over the boiling range necessarily implies that the method can be applied with confidence only to straight-run cuts. Watson et al. provide corrections to the volumetric average boiling point determined from the distillation curve to estimate the molal average boiling point. If the distillation curve is symmetrical, the volumetric average boiling point is close to the 50% evaporation temperature. Further, unless the slope of the distillation curve is steep, the value of K is not significantly altered if the volumetric average boiling point is used instead of the molal average boiling point. Present Work It was felt desirable to develop a mathematical relationship among the properties-characterization factor, average boiling point, and molecular weight-in order to facilitate machine computations. Mair and Willingham (1936) and Lucy (1938) proposed that molecular weights of heavy hydrocarbons of low volatility could be estimated from the relation
M = (T1/35)2'276 (2) where M = molecular weight TI= uncorrected distilling temperature at 1 mmHg, K. If the vapor pressure/temperature relationship of hydrocarbons can be approximated by the Clapeyron equation as is usually done, then T I can be related to the normal boiling point. It was therefore decided to investigate whether the correlation of Watson et al. presented in chart form could be approximated by the following expression M = (Tb/A)B (3) where M = molecular weight, T b = boiling point, K, and A and B are functions of K , the characterization factor. Equation 1 can also be written as In M = B In Tb - B In A (4) A large number of sets of values of M , Tb, and K were read off from the chart of Watson et al. (1935). A linear regression analysis between values of h M and In T b showed that they indeed vary linearly, the correlation coefficient in all cases being in excess of 0.99. The analysis also 0 1985 American Chemical Society
