Ind. Eng. Chem. Res. 1993,32, 1367-1371
1367
PROCESS DESIGN AND CONTROL Effect of Uncertainties in Thermodynamic Data and Model Parameters on Calculated Process Performance Wallace B. Whiting' and Ting-Man Tong Department of Chemical Engineering, West Virginia University, Morgantown, West Virginia 26506-6101
Michael E. Reed Morgantown Energy Technology Center, Morgantown, West Virginia 26507-0880
A technique is presented to use Monte Carlo simulation with Latin-hypercube sampling to determine the uncertainty of process operation to thermodynamic model parameters. Significant changes in process variables are shown for superfractionator and flash examples when reasonable uncertainties are applied to pure-component and binary-interaction parameters in an equation of state model. These uncertainties are quantified both as cumulative probability distributions and as sensitivity coefficients. Thus, the technique provides a quantitative alternative to design safety factors, and it can be used as an indicator of data deficiencies. The examples show the importance of using consistent error estimates for thermodynamic data and the importance of studying the uncertainty of operation of proposed units. Introduction Chemical processes are commonly designed using process simulators. These simulators depend on thermodynamic models for the basic properties. While process engineers understand that they must specify the best available thermodynamic model to get reasonable results, they need some quantitative measure of the design uncertainty, based on the fact that the chosen thermodynamic model is not perfect. Previous studies (for example, Zudkevitch, 1980;Larsen, 1986; Zeck, 1991)have demonstrated the consequences of bad thermodynamic data and inaccurate models. Our approach to design uncertainty is based on Monte Carlo simulation with assumed model parameter uncertainty (Reed and Whiting, 1992). This approach goes beyond first-order sensitivity analysis (Nelson et al., 1983; Macchietto et al., 1986) and provides a quantitative alternative to the application of safety factors in design. Two further problems are addressed in this paper. Once a design is chosen, the uncertainty that is important is the operationaluncertainty, i.e., how must the unit be operated to attain a given specification (e.g., product purity). The other problem is how to specify the uncertainty of model parameter inputs. In the examples, we used the Soave-Redlich-Kwong equation of state for equilibrium properties. The uncertain inputs were the critical pressures, critical temperatures, acentric factors, and binary-interaction parameters. Uncertainty of Model Parameters As pointed out by Zeck (19911, the accuracies of data and model parameters should be estimated and reported. For our examples, we use one of the few sources of such information, the DIPPR compilation (Daubert and Danner, 19891, for pure-component uncertainties. This is probably the most comprehensive source of pure-com-
* Author to whom correspondence should be addressed.
ponent thermophysical data quality judgments. For critical temperatures, critical pressures, and vapor pressures, quality codes are given for each entry. These quality codes are based on the commonly used average absolute relative deviation (AARD) (Daubert et al., 1990;Danner, 1992). To apply the Monte Carlo uncertainty technique, we assumed that each of the pure-component inputs to the model (criticalpressures, criticaltemperatures, acentric factors) were characterized by Gaussian distributions with means equal to the listed values in the DIPPR compilation and with standard deviations such that the Gaussian distribution has the AARD listed. Thus, we used a commonly available, critical compilation for pure-component parameters. For the acentric factars, we calculated the uncertainty from the propagation of the vapor pressure uncertainty, as the acentric factors were derived from the vapor pressure correlation. For the binary-interaction parameters, we used Gaussian distributions with means equal to the values used by Macchietto et al. (1986) in their first-order sensitivity investigation of the two examples that we use (for comparison) and with a standard deviation equivalent to an average absolute deviation of 0.015 for each binary-interaction parameter. We chose this value because reported binary-interaction parameters often very by about this much from one study to another for a given system. A more formal approach would be to use the uncertainty generated directly from the regression technique used to determine the binary-interaction parameter (Whiting et al., 1992). We did not use this approach here so as to simplify the interpretation of our results of the two examples.
The Monte Carlo Simulation The details of our approach are given elsewhere (Reed and Whiting, 1992). It is similar to the approach of Rubin (Diwekar and Rubin, 1991; Frey and Rubin, 1992) in studying the uncertainty of processes to design variable uncertainty. We use 100 simulation runs, which gives a 90% confidence region of 0.85-0.95 on the 0.90 fractile.
O S S S - ~ ~ ~ ~ I ~ ~ I ~ S 0~ 1993 ~ - ~American ~ ~ ~ ~Chemical O ~ . OSociety ~ I O
1368 Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993
This means that the probability that the true 0.90 fractile is either smaller than our calculated 0.85 fractile or greater than our calculated 0.95 fractile is only 10% (Morgan and Henrion, 1990). It is important to note that the number of uncertain input variables does not affect the number of simulation runs required for a given level of confidence. In our examples below, we use 18and 25 uncertain inputs, but we could look at a larger number of input variables without increasing the computational effort. We use the Latin hypercube sampling technique for choosing the values of the model parameters for each run (McKay et al., 1979; Iman and Shortencarier, 1984). In this technique, the probability distribution for each input variable is divided into 100 regions of 0.01 probability each. One value is chosen from each region at random, weighted by the distribution. Then, 100 sets of input variables are constructed by choosing one value for each variable at random from the values obtained above,without replacement. The advantage of this type of so-called stratified sampling technique is that it produces a better approximation of the output probability distribution than would normal Monte Carlo simulation for the same number of runs. This is especially true if many input parameters are being studied, but a few are much more important than the others (Morgan and Henrion, 1990;Ripley, 1988).
Statistical Measures of Results The results of the Monte Carlo simulation are presented in several ways. First, the cumulative frequency distribution for the output variable is shown (Figure 1, for example). This distribution quantifies the uncertainty of the output variable (in this case reflux ratio) and allows one to determine a safety factor for a chosen level of risk. For example, if we want to be 90% certain that we can meet the product specification,we need to be able to handle a reflux ratio of 13, even though the nominal design reflux ratio is about 3.5. If the uncertainty in the input model parameters is less, this curve will be steeper, and our expected range of reflux ratios will be smaller for a chosen level of risk. Other statistical measures are used to study the importance of the input variable uncertainties to the output variable uncertainties. Both the input variables (Y) and are standardized by subtracting the output variables (Xi) the mean (Yor Xi) and dividing by the standard deviation (Sy or SX,).A normal multilinear regression is then performed on a standardized output variable and the standardized input variables to obtain the best-fit regression equation:
The regression coefficients are known as standardized regression coefficients (SRCi), and they are measures of the first-order sensitivity of the output variable to the input variables. Because they are standardized by their uncertainties, the SRC’s are unitless and can be compared directly with one another. The SRC’s with the largest absolute values, thus, correspond to the most important variables in a statistical sense. In other words, if X i is changed by one standard deviation, Ywill change by SRCi times its standard deviation. Next, one input variable is ignored and another multilinear regression is performed. The relative improvement in the regression model obtained with all input variables over that obtained with all but one input variable is another measure of the importance of the excluded input variable. The square root of this value is called the partial correlation
Table I. Suwrfractionator S w i f i c a t i o n s number of theoretical stages 55 feed stage 36O column pressure 8.816 bar distillate purity 0.0008 mole fraction ethane bottoms purity 0.0053 mole fraction ethylene feed pressure 8.816 bar mole fraction vapor 0.27 methane 0.16 kmol/h ethylene 894.53 kmol/h ethane 717.58 kmol/h propylene 5.18 kmol/h a
Partial condenser is stage 1; reboiler is stage 55.
coefficient (PCC). The PCC is related to the SRC by (PCCi)2= (SRCJ2(1- R,;)/(l-
R?)
where Rx? is the coefficient of determination of the multilinear regression of Xi as a function of Y and all the other X s and Ry2 is the coefficient of determination of the multilinear regression of Y as a function of all of the XS.
Put another way, the PCC is a measure of the effect of the excluded variable, if all of the other variables were held constant. The input variables with the largest magnitude PCC are the most important variables in the sense that their contributions to the regression model cannot adequately be accounted for by the other input variables. That is, if we ignore the input variable with the largest PCC, we have a greater degradation of the model than we would have if we ignored any other single input variable. We calculate the SRC and the PCC for each input variable to decide which are the most important input parameters to determine more accurately, i.e., to which input variables is the output variable most sensitive. The greater the absolute value of SRC, the more important is the input variable, and similarly for PCC. Because of the subtle differences between these two parameters, slight differences in the ranking of the input variables can result. (For details, see Edwards (19851.) These results can be the basis for a traditional Pareto analysis for model improvements. The SRC’s and PCC’s are based on linear regressions. As such, they tend toward less reliability for nonlinear problems. Therefore, we take the raw data for each input or output variable and substitute its rank (1-100) and then perform the regression calculationsagain. The results are called standardized rank regression coefficients (SRRC) and partial rank correlation coefficients (PRCC),and they are used in the same way as are the SRC and PCC above. The SRRC and PRCC, however, have been shown to be more reliable measures for systems that either are highly nonlinear or may contain rogue data (Iman and Conover, 1979). The details of the calculations of the SRC, PCC, SRRC, and PRCC are given by Iman et al. (1985). Superfractionator Example This example is originally from Hernandez et al. (1983). Macchietto et al. (1986) used this example to demonstrate their efficient technique for determining first-order sensitivities. The detailed specifications are given in Table I. There are four components to be separated in a distillation column. The purities of the top and bottom products are specified. We perform tray-by-tray material and energy balances to determine the required reflux ratio and duties for the condenser and reboiler. The results of the Monte Carlo simulation are shown in Figures 1-3. The
Ind. Eng. Chem. Res., Vol. 32,No. 7,1993 1369 1W
1
Table 11. Superfractionator Regression Results' (Output Variable: Reflux Ratio)
input variable Tc(CH4) Pc(CH4) w(CH4) Tc(C2Hd Pc(C2Hd w(CnH4) Tc(C2Hd Pc(C2Hd w(CzHd Tc(Cd%) Pc(C3Hd w(CsHs) k(CHJC2H4) k(CHJC2He) k(CHJCsH6) k(C&/C&) k(C&/C&) k(C2HdC&)
Reflux Rauo
Figure 1. Uncertainty of reflux ratio for superfractionator.
R2
PCC -0.174 (3) -0.025 (16) -0.090 (9) 0.382 (2) -0.083 (10) -0.108 (7) -0.388 (1) -0.069 (11) 0.051 (13) 0.124 (4) 0.018 (17) 0.036 (15) -0.096 (8) 0.122 (5) 0.053 (12) -0.051 (14) -0.016 (18) 0.177 (6)
SRC -0.147 (3) -0.021 (16) -0.075 (9) 0.347 (2) -0.070 (10) -0.091 (7) -0.351 (1) -0.058 (11) 0.042 (14) 0.105 (4) 0.015 (17) 0.030(15) -0.080 (8) 0.103 (5) 0.044 (12) -0.043 (13) -0.014 (18) 0,099(6) 0.315
PRCC -0.017 (18) -0.136 (7) 0.130 (9) 0.786 (2) -0.210 (4) 0.131 (8) -0.808 (1) 0.188 (6) -0.209 (5) -0.044 (15) -0.102 (11) 0.033(16) -0.019 (17) -0.069 (13) 0.099 (12) -0.340 (3) -0.049 (14) -0.126 (10)
SRRC 0.008 (18) -0.063 (7) 0.061 (9) 0.588 (2) -0.099 (4) 0.061 (8) -0.634 (1) 0.088 (6) -0.099 (5) -0,020(15) -0,048(11) 0.015 (16) -0.009 (17) -0.032 (13) 0.046 (12) -0.167 (3) -0.023(14) -0.059(10) 0.789
a The values in parentheses are the ranks of the absolute values of the coefficients.
20
40 MI Q(Rebo1ler) (Glih)
80
1W
Figure 2. Uncertainty of reboiler heat duty for superfractionator.
0 0
20
40
60
Q(Candcarer) (Glh)
Bo
lW
Figure 3. Uncertainty of condenser heat duty for superfractionator.
cumulative frequency distributions are quite broad, indicating that we must be prepared for a wide range of reflux ratios and heat-exchanger duties if we are to be reasonably certain that we can meet the specifications. We have included in the simulation only uncertainties in the critical temperatures, critical pressures, acentric factors, and binary-interaction parameters. To investigate the meaning of these results and the path to a solution, we study the SRC, PCC, SRRC, and PRCC. Table I1 shows the results for the reflux ratio. The other two output variables (heat duties for the condenser and reboiler) show very similar results. We have examined three output variables while varying 18 uncertain input variables. However,three of the input variables seem most responsible for the uncertainties in the outputs: critical temperature of ethane, critical temperature of ethylene, and (ethane/ethylene)binary-interactionparameter. Thus, we can get the most improvement if we can get better data for these three model inputa. The binary-interaction
parameter requires binary phase-equilibrium data, and better data or more data of the same quality will improve this uncertainty. However, the effect of the uncertainty in the critical values is surprising. If we look at the quality codes used in the DIPPR compilation, we find that they are in broad ranges (e.g., 0.2%, 1% , 3% , 5 % 1. We chose for our illustrative examples the high end of the range. If we had chosen the low end of the ranges, these correlation and regression Coefficients would have been less significant in this example,but these critical temperatures would still be among the most significant input variables. This leads to a common problem in safety factor and uncertainty analysis: safety factors are added to safety factors. In this case, even the low range of the quality codes seems too high in some cases (that is, a safety factor has been applied). Then, we chose the high end of the range (because we had no basis for preferring any specific value within the range). That is, we added another safety factor. Notice that we did not add the safety factors in the case of the binary-interaction coefficients. We merely applied our logic to the data available. The number of data were not great, so we cannot be sure that our value of 0.015 for the average absolute deviation is especially accurate. However, we did not add in the bias of a safety factor. Three of the 100 simulations resulted in infeasible situations. The minimum number of stages was greater than the number of stages specified for the column. For the PCC and SRC calculations, these three simulations were excluded, as there was no reasonable value for the output variables. However, for the calculations that only involve the rank of the variables (PRCC and SRRC), these three simulations were included with ranks indicating how many additional stages would be necessary to meet the product specifications. This is a good example of the advantage of the rank coefficients. Also, the coefficient for the standardized rank regression of determination (R2) was 0.8, whereas it was only 0.3 for the standardized regression. This latter figure shows that only 30% of the variation in the output variables can be explained by a first-order sensitivity analysis. Although we treated all of the input variables the same in the correlation and regression coefficients calculations, there are two distinct classes of input parameters in this problem. Whereas the critical temperatures, critical pressures, and acentric factors are essentially directly measurable quantities, independent from the thermody-
Ind. Eng. Chem. Res., Vol. 32,No. 7, 1993 Table IV. Flash Regression Results. (Output Variable: VaDor Methane Flow Rate) ~
5
40 45 50 Lle‘hane Vapor F:ourare (kmo’lhl
is
60
Figure 4. Uncertainty of flash. Table 111. Flash Smcifications temperature pressure feed carbon dioxide hydrogen sulfide methane ethane propane
225 K 60.78 bar 6 kmol/h 24 kmol/h 66 kmol/h 3 kmol/h 1kmol/h
namic model, the binary-interaction parameters are determined using the chosen equation of state. The uncertainty of these parameters, then, is some (indirect) measure of the uncertainty of the chosen model. While we can improve the uncertainty of the first set of variables by better experiments,we need to improve the model (and, possibly, the vapor-liquid-equilibrium experiments) to improve the uncertainty of the binary-interaction parameters.
Flash Example This example is originally from Michelsen (19801, and a sensitivity analysis was reported by Macchietto et al. (1986). The details are given in Table 111. A fivecomponent system is flashed at a specified temperature and pressure. The output variable that we track is the methane vapor flow rate. Again, the Soave-RedlichKwong equation of state is used. The results of the Monte Carlo simulation are given in Figure 4. The cumulative frequency distribution indicates the broad range of possible outcomes of this flash unit. Either the downstream units need to be able to handle this uncertainty, or the conditions of the flash need to be studied for greater flexibility. Table IV shows the results of the regression study. For this case, 1 output variable and 25 input variables were considered. The most important variables are the (H& CHI) binary-interaction parameter, critical temperature of methane, critical pressure of hydrogen sulfide, and critical temperature of hydrogen sulfide. We note again that the uncertainties for the critical constants are probably overestimated by the broad quality code ranges. The coefficient of determination again was higher for the rank regression (0.949)than for the normal regression (0.843).
Conclusions The Monte Carlo simulation approach can provide insight into the uncertainty of the operation of a specified unit caused by the inaccuracies of the thermodynamic model, as a quantitative alternative to safety factors. This approach requires an unbiased (even if imprecise) assessment of the uncertainties of the input parameters to the
PCC -0.165 (8) 0.085 (17) -0.107 (12) 0.452 (4) 0.489(3) 0.372(6) -0.845 (1) 0.125 (9) -0.036 (24) 0.079(18) 0.108(11) 0.093 (14) -0.007 (25) -0.086 (16) -0.059 (19) -0.101 (13) 0.419 (5) 0.124 (10) 0.054 (22) 0.792 (2) 0.058(20) 0.057 (21) 0.286 (7) 0.091(15) 0.041(23)
SRC -0.066 (8) 0.034 (17) -0.043 (12) 0.201 (4) 0.222(3) 0.159(6) -0.626 (1) 0.050(9) -0.014 (24) 0.031(18) 0.043(11) 0.037(14) -0.003 (25) -0.034 (16) -0.023 (19) -0.040 (13) 0.183(5) 0.050 (10) 0.021 (22) 0.515 (2) 0.023 (20) 0.023 (21) 0.119 (7) 0.036(15) 0.016(23) 0.843
PRCC 4.010 (24) 0.448 (6) 0.158 (14) 0.771 (4) 0.828(3) 0.676(5) -0.931(2) 0.001(25) 0.036(21) -0.221 (9) -0,160(13) 0.066 (19) 0.020 (22) -0.103 (17) -0.176 (12) -0.275 (8) 0.383 (7) 0.080 (18) 0.011(23) 0.937 (1) -0.061 (20) -0.212 (10) 0.210 (11) 0.125 (16) -0.140(15)
SRRC -0.002(24) 0.114(6) 0.036 (14) 0.274 (4) 0.335 (3) 0.207 (5) -0.580 (2) O.Oo0 (25) 0.008 (21) -0.052 (9) -0.037 (13) 0.015 (19) 0.004(22) -0.023 (17) -0,041 (12) -0.065(8) 0.094(7) 0.018 (18) 0.002 (23) 0.612(1) -0.014 (20) -0.049 (10) 0.049 (11) 0.029(16) -0.032 (15) 0.949
a The values in parentheses are the ranks of the absolute values of the coefficients.
thermodynamic model. The results of the Monte Carlo simulation can be used also to direct further study to reduce those parameter uncertainties that are of greatest importance. This approach includes the first-order sensitivity analysis, but as is shown in the superfractionator example, reasonable uncertainty estimation can be made even when a first-order sensitivity analysis would only be able to explain 30% of the variation in the output variable. Many input variables can be considered, without increasing the number of simulations or compromising the quality of the result.
Literature Cited Danner, R. P. Personal Communication, 1992. Daubert, T. E.; Danner, R. P. Physical and Thermodynamic Properties of Pure Chemicals: Data Compilation (DIPPR); Hemisphere: New York, 1989. Daubert, T. E.; Sibul, H. M.; Steggins, C.; Kendall, R. Policies and Procedures DocumentingCompilation,Prediction, and Correlation for DIPPR Data Compilation Project. AIChE Symp. Ser. 1990, 86 (279),62. Diwekar, U. M.; Rubin, E. S. Stochastic Modeling of Chemical Processes. Comput. Chem. Eng. 1991,15, 105. Edwards, A. L. Multiple Regression and the Analysis of Variance and Covariance, 2nd ed.; W. H. Freeman: New York, 1985. Frey, H. C.; Rubin, E. S.Evaluate Uncertainties in Advanced Process Technologies. Chem. Eng. Prog. 1992,88(5),63. Hernandez, M. R.; Gani, R.; Romagnoli, J. A.; Brignole, E. A. The Prediction of Properties and Ita Influence on the Deaign and Modeling of Superfractionators. Proceedings, Second International Conference on Foundations of Computer-Aided Process Design, Snowmass, CO, June 1983;CACHE Publications: Ann Arbor, MI, 1984;pp 709-740. Iman, R. L.;Conover, W. J. The Use of the Rank Transform in Regression. Technometrics 1979,21,499. Iman, R. L.;Shortencarier, M. J. “A FORTRAN 77 Program and User’s Guide for the Generation of Latin Hypercube and Random Samples for Use with Computer Models”; Report NUREG/CR3624;National Technical Information Service: Springfield, VA, 1984. Iman, R. L.;Shortencarier, M. J.; Johnson, J. D. ‘A FORTRAN 77 Program and User’sGuide for the Calculationof PartialCorrelation
Ind. Eng. Chem. Res., Vol. 32, No. 7, 1993 1371 and Standardized Regression Coefficients";Report NUREG/CR4122; National Technical Information Service: SpringField, VA, 1985. Larsen, A. H. Data Quality for Process Design. Fluid Phase Equilib. 1986, 29, 47. Macchietto, S.; Maduabeuke, G.; Szczepanski, R. Exact Determination of Process Sensitivity to Physical Properties. Fluid Phase Equilib. 1986, 29, 59. McKav. M. D.: Beckman. R. J.: Conover. W. J. A ComDarison of Thiee Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code. Technometrics 1979, 21, 239. Michelsen, M. L. Calculation of Phase Envelopes and Critical Points for Multicomponent Mixtures. Fluid Phase Equilib. 1980,4, 1. Morgan, M. G.; Henrion, M. Uncertainty: A Guide t o Dealing with Uncertainty in Quuntitative Risk and Policy Analysis; Cambridge University Press: Cambridge, UK, 1990. Nelson, A. R.; Olson, J. H.; Sandler, S. I. Sensitivity of Distillation Process Design and Operation to VLE Data. Znd. Eng. Chem. Process Des. Dev. 1983,22, 547.
Reed, M. E.; Whiting, W. B. Sensitivity and Uncertainty of Process Designs to Thermodynamic Model Parameters: A Monte Carlo Approach. Chem. Eng. Commun. 1993, in press. Ripley, B. D. Uses and Abuses of Statistical Simulation. Math. Program. 1988,42, 53. Whiting, W. B.; Hwang, M.-J.; Anderson, T. F. Confidence Regions, MaximumLikelihood,and Phase-Equilibrium ModelUncertainty. Submitted for publication in Can. J. Chem. Eng. 1993. Zeck, S. Thermodynamics in Process Development in the Chemical Industry-Importance, Benefits, Current State and Future Development. Fluid Phase Equilib. 1991, 70, 125. Zudkevitch, D. Forensic Thermodynamics: Erroneous Decisions on Thermodynamic Data Can Cause Plant Failures. EFCE Publ. Ser. 1980, No. 11, 885. Received for review November 2, 1992 Revised manuscript received March 22, 1993 Accepted April 5, 1993
