Ind. Eng. Chem. Res. 1997, 36, 23-41
23
Mechanism and Temperature-Dependent Kinetics of the Dehydration of tert-Butyl Alcohol in Hot Compressed Liquid Water Xiaodong Xu and Michael Jerry Antal, Jr.* Hawaii Natural Energy Institute and the Department of Mechanical Engineering, University of Hawaii at Manoa, Honolulu, Hawaii 96822
Donald G. M. Anderson Aiken Computation Laboratory, Harvard University, Cambridge, Massachusetts 02138
In the presence of acid or base, or under neutral conditions, isobutylene is the only observed product of the reactions of tert-butyl alcohol in compressed liquid water at 225, 250, and 320 °C. On the basis of a kinetic analysis of limited results at 250 °C, an earlier paper concluded that tert-butyl alcohol dissociates as an Arrhenius acid at 250 °C and thereby catalyzes its own dehydration to isobutylene. Kinetic analyses of the data sets contained in this paper do not corroborate the alleged ability of tert-butyl alcohol to protonate water. Instead, we find that hydronium ions formed by the ordinary dissociation of water are the primary catalytic agents for the dehydration reaction. In agreement with the earlier work, all three data sets are consistent with a heterolytic reaction mechanism involving protonated alcohol, carbocation, ditert-butyl ether, and protonated ether as intermediates. The kinetics still suggest that tertbutyl alcohol weakly dissociates as a Bronsted acid at these conditions. Values of Kw determined by the kinetic model enjoy good agreement with electrochemical values available in the literature. Values of the other parameters (rate constants) which compose the kinetic model are evaluated with less precision, and a novel method is proposed to estimate the uncertainty associated with each parameter. Introduction During a visit to Honolulu 7 years ago, Professor Maitland Jones, Jr., of Princeton University inquired of us if a kinetic study of the acid-catalyzed dehydration of tert-butyl alcohol in near-critical water could discriminate between various possible mechanistic models for the reaction chemistry and yield the accepted mechanism. His question was motivated by a desire to understand the depth of mechanistic insight that can be gleaned from studies of reaction rates. At that time, the reaction chemistry of tert-butyl alcohol in water was thought to be well understood (see Xu and Antal, 1994 for a review); consequently, a study of its kinetics was certain to offer insight into the utility of kinetic studies. We were unable to answer Jones’ question with confidence. Students of chemical engineering are taught that almost any reasonable model with enough free parameters should be able to fit available experimental data. For example, the authoritative textbook by Hill (1977) teaches “...with the flexibility given by these (rate) constants it is not hard to get a good fit of data to a mathematical expression. In the words of one of my colleagues, ‘with six constants you can fit a charging rhinoceros.’ ” On the other hand, chemists learn that “the most important single tool for mechanistic investigation is probably kinetics. Measurements of reaction rates are made for the purpose of verifying the consistency of the proposed pattern of reaction steps through dependence of rates on concentration, and also in order to derive activation parameters, ...” (Lowry and Richardson, 1987). Thus, this study was initiated with the goal of shedding light on the ability of kinetics to detail mechanism. A closely related secondary goal was to * Author to whom all correspondence should be addressed at the Hawaii Natural Energy Institute, Holmes Hall 305, University of Hawaii at Manoa, Honolulu, HI 96822. Phone: 808/956-7267. Fax: 808/956-2336. E-mail: antal@ wiliki.eng.hawaii.edu. S0888-5885(96)00349-1 CCC: $14.00
assess the reliability of the values of the parameters associated with the best-fitting model. As described in an earlier paper (Xu and Antal, 1994), the experimental work associated with this study was successful, but the results were unexpected. We found that isobutylene forms from tert-butyl alcohol in the absence of acid in hot (250 °C), compressed (34.5 MPa) liquid water (see also Ogi et al., 1990). The reaction is unusually selective: isobutylene is the only product detected by gas chromatography and mass spectrometry. It is also surprisingly fast: equilibrium is attained after about 30 s. Addition of as little as 10-5 mol/dm3 sulfuric acid to the reactant mixture noticeably enhanced the reaction rate without introducing any detectable byproducts. With the intention of explicating these unanticipated results, we developed and elaborated three chemically motivated, root models for the observed reaction chemistry (Xu and Antal, 1994). Each model presupposed acid-catalyzed, heterolytic reaction schemes wherein the needed protons were supplied by the dissociation of either water or tert-butyl alcohol. Two of the models assumed the dehydration chemistry to be governed by an E1 type mechanism, and one of the models posited an E2 mechanism. Although many combinations of many elementary reaction steps (see Table 2 of Xu and Antal, 1994) were examined as subsets of these three root models, only one model with one relatively simple combination of steps (see Figure 1) was able to fit the available data. This model involved protonated tert-butyl alcohol, carbocation, ditert-butyl ether, and protonated ether as intermediates in the reaction network. Protons were supplied by the dissociation of tert-butyl alcohol as both an Arrhenius (eq 2 in Figure 1) and a Bronsted (eq 4) acid. Isobutylene was formed by the unimolecular (eq 11) decomposition of the ether. Although the literature contains information concerning synthesis of the ether (Petrov, 1889; Evans and Edlund, 1936) and its thermodynamic properties (Fenwick et al., 1975), it is not commercially © 1997 American Chemical Society
24
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997
Figure 1. Kw/Auto/E1/AdE2/Uni/Bi model for the dehydration of tert-BuOH in near-critical water.
available. Consequently, we did not undertake studies of its reaction chemistry. To obtain agreement between the model and the experimental data, the pKa of tertbutyl alcohol at 250 °C was predicted to be about 9, whereas its value at normal temperature and pressure (NTP) is 18. In summarizing this earlier work, we would be remiss if we did not call attention to the startling ability of this model to fit unexpected and unusual trends in the experimental data (see Figure 4 of Xu and Antal, 1994). Dramatic variations of the pKa of various acids with increasing temperature are well-known. For example, the pKa of HF increases from 3.2 to about 6.7 as the temperature increases from 20 to 300 °C, whereas the pKa of NH4+ decreases from 9.5 to about 4.2 over the same temperature range (Ellis and Mahon, 1977). Nevertheless, we were skeptical of the 9 order of magnitude decrease in the Ka of tert-butyl alcohol that was required to enable the kinetic model to fit the dehydration data at 250 °C. To check this result, we undertook studies of the reaction kinetics at 320 and 225 °C. Both the reaction rate and the extent of conversion were large at the higher temperature; consequently, some experiments were executed in the presence of base NaOH to slow the rate of conversion. These
experiments were successful, but the model was unable to fit the new data. Moreover, we discovered that the model could not fit data at 250 °C that incorporated new NaOH results with the earlier measurements. This paper describes the steps we took to comprehend the kinetic and mechanistic implications of these new findings. The primary objective of this study was to develop a kinetic model based on elementary reactions that could describe the effects of temperature, residence time, reactant concentration, and acid and base concentration on the yields of isobutylene from tert-butyl alcohol in hot liquid water. In light of the successful outcome of earlier studies of the reactions of alcohols (Antal et al., 1987; Ramayya et al., 1987; Narayan and Antal, 1989, 1990; Xu et al., 1990, 1991), acids (Mok et al., 1989; Carlsson et al., 1994), and sugars (Antal et al., 1990a,b, 1991) in hot liquid water and supercritical water, this objective did not seem to be unrealistic. Apparatus and Experimental Procedures Four plug flow reactors fabricated from Hastelloy C276 tubing were employed in this work. Experiments conducted before Sept 9, 1993 involved two reactors that were described in detail in an earlier paper (Xu and Antal, 1994). Experiments conducted after Sept 9, 1993
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997 25
used two newly fabricated reactors. The new, annular flow reactor employed a 61 cm (24 in.) long piece of 0.635 cm (0.25 in.) o.d. × 0.457 cm (0.18 in.) i.d. tube for its outer shell and a 0.397 cm (5/32 in.) o.d. sintered alumina tube for its annulus. A 27 cm (10.63 in.) long section of 0.159 mm (0.0625 in.) o.d. × 0.108 mm (0.0425 in.) i.d. capillary tubing constituted the fourth reactor. These four reactors enabled studies of a broad range of residence times at reaction temperatures of 225, 250, and 320 °C. Following design practices developed early in our work (Antal et al., 1987), good temperature profiles (i.e., measured maximum excursions of no more than 3 °C from the desired reaction temperature over the length of the isothermal zone), with abrupt transitions from room temperature to the reaction temperature at the entrance and exit of the reactor, were always obtained. Criteria which ensure the legitimate use of the plug flow idealization for these reactors were discussed in Appendix I of the paper by Antal and Xu (1994), which was based on the earlier work of Cutler et al. (1988) and Ramayya and Antal (1989). Occasionally, colleagues have expressed concern that an isobutylene-rich phase could form in the reactor (Schneider, 1973), thereby invalidating the plug flow assumption. Such an occurrence would presumably give rise to nonreproducible yields, particularly among different reactors, but such nonreproducibilities have never been observed. NIST researchers (Ochel, 1994, Ochel et al., 1994) have observed one-phase equilibrium mixtures in a similar system at temperatures between 200 and 300 °C with an initial tert-butyl alcohol concentration of 1 M/dm3. Since all the work presented in this paper involved concentrations below this level, there is no reason to believe that phase separation has compromised our results. Liquid products from experiments conducted before Sept 9, 1993 were analyzed by high-performance liquid chromatography (HPLC) and gaseous products by gas chromatography (GC) according to procedures described earlier (Xu and Antal, 1994). Product analysis for the more recent experiments was performed using a HewlettPackard Model 5890 series II GC equipped with a thermal conductivity detector and a flame ionization detector. The liquid products were analyzed using a J & W Scientific DB-1 fused-silica capillary column with a 8% hydrogen in helium carrier gas mixture. These analyses used the following temperature program: 3 min at 35 °C followed by a 30 °C/min ramp to 150 °C and a 1 min soak at 150 °C. Gaseous products were separated by an Alltech Carbosphere column that employed the same carrier gas mixture. The temperature program for gas analysis involved a 4.2 min hold at 35 °C, followed by a 15 °C/min ramp to 227 °C and a second 70 °C/min ramp to 350 °C with a subsequent 11.2 min soak at 350 °C. Some qualitative analyses of the gas and liquid samples were accomplished using a Hewlett-Packard Model 5790 GC equipped with a Hewlett-Packard Model 5970 MSD and a J&W FSOT capillary column. In most cases six liquid samples (three from each of two different sampling loops) taken at each reaction condition were analyzed by GC. Two injections were made into the GC from each sample, and the area of each detected peak was recorded. The typical ratio of the sample standard deviation of a peak area measurement to its mean value for a selected reaction condition was less than 3%. Laboratory-grade tert-butyl alcohol from Hawaii Chemical & Scientific was used as the reactant. No
impurities were detected in this reagent by either HPLC or GC analysis. Fisher certified-grade 10 N sulfuric acid was employed as catalyst. Isobutylene supplied by Aldrich Chemical Co. and air standards were used to calibrate the GC. Model Formulation Because many different models were examined during the course of this work and because of their complexity, it is not possible to fully detail each one here. The Kw/ Auto/E1/AdE2/Uni/Bi Model displayed as eqs 1-14 in Figure 1 is representative of the models we considered, and it realizes an effectively perfect fit to all of our experimental data. Consequently, the mathematics of its formulation are presented in full detail here to serve as an illustrative example of model development. Equations 1-14 define a set of six ordinary differential equations (ODE’s) which govern the time-dependent behavior of the concentrations of tert-butyl alcohol [(CH3)3COH or BuOH], di tert-butyl ether [(CH3)3COC(CH3)3 or BuOBu], protonated tert-butyl alcohol [(CH3)3COH2+ or BuOH2+], tert-butyl carbocation [(CH3)3C+ or Bu+], tert-butyloxy anion [(CH3)3CO- or BuO-], and protonated ether [(CH3)3COH+C(CH3)3 or BuOH+Bu]. These six ODE’s are d[BuOH]/dt ) -k1[BuOH][H2O] + k-1[BuO-][H3O+] k2A[BuOH][H3O+] + k-2A[BuOH2+][H2O] 2k2B[BuOH]2 + 2k-2B[BuOH2+][BuO-] + k4B[Bu+] × [BuO-] - k-4B[C4H8][BuOH] - k5[Bu+][BuOH] + k-5[BuOH+Bu] + k6B[BuOH+Bu][BuO-] k-6B[BuOBu][BuOH] + k7[BuOBu] + 2k8[BuOBu] × [H2O] - 2k-8[BuOH]2 + k10[Bu+][OH-] - k-10[BuOH] (1) d[BuOBu]/dt ) k6A[BuOH+Bu][H2O] - k-6A[BuOBu] × [H3O+] + k6B[BuOH+Bu][BuO-] - k-6B[BuOBu] × [BuOH] - k7[BuOBu] - k8[BuOBu][H2O] + k-8[BuOH]2 (2) d[BuOH2+]/dt ) k2A[BuOH][H3O+] - k-2A[BuOH2+] × [H2O] + k2B[BuOH]2 - k-2B[BuOH2+][BuO-] k3[BuOH2+] + k-3[Bu+][H2O] (3) d[Bu+]/dt ) k3[BuOH2+] - k-3[Bu+][H2O] k4A[Bu+][H2O] + k-4A[C4H8][H3O+] - k4B[Bu+] × [BuO-] + k-4B[C4H8][BuOH] - k5[Bu+][BuOH] + k-5[BuOH+Bu] - k9[Bu+][OH-] + k-9[C4H8][H2O] k10[Bu+][OH-] + k-10[BuOH] (4) d[BuO-]/dt ) k1[BuOH][H2O] - k-1[BuO-][H3O+] + k2B[BuOH]2 - k-2B[BuOH2+][BuO-] - k4B[Bu+] × [BuO-] + k-4B[C4H8][BuOH] - k6B[BuOH+Bu] × [BuO-] + k-6B[BuOBu][BuOH] (5) d[BuOH+Bu]/dt ) k5[Bu+][BuOH] - k-5[BuOH+Bu] k6A[BuOH+Bu][H2O] + k-6A[BuOBu][H3O+] k6B[BuOH+Bu][BuO-] + k-6B[BuOBu][BuOH] (6)
26
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997
In these equations the concentration of hydronium cations is [H3O+], the concentration of hydroxyl anions is [OH-], and the concentration of water is [H2O]. To avoid unnecessary integrations, three algebraic equations are employed to represent the concentrations of isobutylene, hydronium ion, and water. Carbon conservation results in the following expression for isobutylene:
[C4H8] ) [BuOH]0 - [BuOH] - [BuO-] [BuOH2+] - [Bu+] - 2[BuOBu] - 2[BuOH+Bu] (7) where [BuOH]0 is the initial concentration of tert-butyl alcohol that enters the reactor. When acid is used as a catalyst, a charge balance offers the following equation:
[H3O+] + [BuOH2+] + [Bu+] + [BuOH+Bu] ) [OH-] + [HSO4-] + [BuO-] (8) Since H2SO4 fully dissociates into H3O+ and HSO4- at the conditions discussed in this paper and HSO4- does not further dissociate (Xu, 1992), we have [HSO4-] ) [H2SO4]0. Employing the relationship Kw ) [H3O+][OH-], we find
1 [H3O+] ) (-S + xS2 + 4Kw) 2
(9)
where
S ) [BuOH2+] + [Bu+] + [BuOH+Bu] - [BuO-] [HSO4-] (10) When NaOH is used (in the absence of H2SO4) to lower the background concentration of [H3O+] and thereby retard the reaction rate, we have [Na+] ) [NaOH]0 and [H3O+] is given by eq 9 with
S ) [Na+] + [BuOH2+] + [Bu+] + [BuOH+Bu] [BuO-] Finally, the concentration of water is assumed to remain constant throughout the reaction with a value given by
[H2O] ) (1 - instantaneous mole fraction FH2O,RTP (11) of tert-butyl alcohol) × 55.55 FH2O,NTP where FH2O,NTP is the density of water at 0.1 MPa and 25 °C, and FH2O,RTP is the density of water at reaction temperature and pressure (RTP). The initial conditions (t ) 0 at the entrance of the reactor) are as follows:
[BuOH]0 ) (NTP concentration of tert-butyl alcohol) ×
FH2O,RTP
(12)
FH2O,NTP [C4H8]0 ) [BuOBu]0 ) [BuOH2+]0 ) [Bu+]0 ) [BuOH+Bu]0 ) [BuO-]0 ) 0 (13)
[H2O]0 ) (1 - initial mole fraction FH2O,RTP (14) of tert-butyl alcohol) × 55.55 FH2O,NTP and
1 [H3O+]0 ) ([H2SO4]0 + x4Kw + [H2SO4]02) (15) 2 when H2SO4 is present or
1 [H3O+]0 ) (-[NaOH]0 + x4Kw + [NaOH]02) (16) 2 when NaOH is present. Since the intermediate species included in the model are undetected by our analytic instrumentation, some assumptions must be made concerning their disposition at the exit of the reactor. A best fit is achieved when both ether and protonated ether are converted to two molecules of tert-butyl alcohol, the tert-butyloxy anion reverts to tert-butyl alcohol, and the protonated alcohol and its carbocation form isobutylene at the exit of the reactor. Since model values of the concentrations of these undetected intermediates are always low, these assumptions do not strongly influence the fit of the model to the data. Finally, to calculate the residence time of the reactant at RTP, the density of the tert-butyl alcohol solution is assumed to be that of pure water at the same conditions. Modeling studies indicate that this approximation has no significant effect on the results (Xu and Antal, 1994). Parameter Estimation The rate constants kj, j ) 1, 2, ..., n, are parameters whose values are chosen to give a best fit (in the leastsquares sense) of the model to the experimental data. For any given condition, the quality of this fit is measured by the residual ei ) (Ciexp - Cimod), where Ci is the concentration of tert-butyl alcohol exiting the reactor at reaction condition i ) 1, 2, ..., m as measured experimentally (exp) or calculated by the model (mod). In this paper we optimize the fit of the model to the data by minimizing the value of the least-squares objective function χν2: m
χν2 )
(ei/σi)2/ν ∑ i)1
(17)
where σi is the sample standard deviation associated with the measurement of Ciexp and ν is the number of degrees of freedom. As is evident from eq 17, the value of σi plays a pivotal role in the evaluation of χν2. During the past decade a large number of statistical studies (Xu et al., 1991) of the accuracy of both HPLC and GC instruments with automatic injectors have consistently shown σi ) µCiexp , where µ ) 0.04 for experiments with inlet reactant concentration [BuOH]0 > 0.05 mol/dm3 and µ ) 0.08 when [BuOH]0 e 0.05 mol/dm3. The larger value of µ reflects our recent experience, which indicates that there is a significant loss of precision when low concentrations ( 0, contains the actual parameter vector, with a probability increasing with . Assuming a
Gaussian distribution for the measurement errors, the traditional choice is
)
n F(n,m-n;R) m-n
associated with probability 1 - R, where F is a known, tabulated distribution function (see, for example, Myers, 1990). For R ) 0.01 and typical values of m and n, this yields an of order unity, which may exceed the range of validity of the simplifying assumptions. Bates and Watts (1988) discuss maximum likelihood, sampling theory, and Bayesian interpretations of the inference region, and also several alternative approaches which have been deemed too cumbersome and uneconomical for use in the present circumstances. For moderately large n, it is difficult to come to grips with such influence regions. We therefore consider probing the boundary of the region in selected directions. For given > 0, and for any direction d, with ||d|| ) 1, define δˆ (d|) as the smallest positive δ such that φ(xˆ +δd) ) (1 + )φ(xˆ ). In principle, it is not difficult to calculate δˆ (d|) by solving a simple univariate nonlinear root-finding problem. We have
1 φ(xˆ +δd) ) φ(xˆ ) + δ2d*H(xˆ ) d + O{δ3} 2 and we than obtain, for sufficiently small δˆ (d|)
δˆ (d|) ≈ {2φ(xˆ )/d*H(xˆ ) d}1/2 or
δˆ (d|) ≈ x||f(xˆ )||/{||J(xˆ ) d||2 + d* S(xˆ ) d}1/2 If we further have
|d*S(xˆ ) d| , ||J(xˆ ) d||2 we see that
δˆ (d|) ≈ x||f(xˆ )||/||J(xˆ ) d|| For a given d, we define the insensitivity indicator
34
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997
η(d) :)
1 ||f(xˆ )||/{||J(xˆ ) d||2 + d* S(xˆ ) d}1/2 10
which approximates the distance one must move from xˆ in the direction d to increase φ by 1%, provided η(d) is sufficiently small. If η(d) is not small enough to invoke this interpretation, we may still regard η(d) as a measure of insensitivity to change in direction d. Note that η(d) is well-defined for all d if H(xˆ ) is positive definite. Because of the difficulty of evaluating S(xˆ ), we are motivated to consider a simplified insensitivity indicator
η˜ (d) :)
1 ||f(xˆ )||/||J(xˆ ) d|| 10
which is well-defined for all d if J(xˆ ) has maximal rank. We shall focus primarily on η˜ (d) hereafter. Note that while δˆ (-d|) * δˆ (d|), in general, we have η˜ (-d) ) η˜ (d) and also η(-d) ) η(d). Two sets of orthogonal directions are of particular interest. The first set consists of the coordinate directions, ej, j ) 1, 2, ..., n, whereby we measure insensitivity to changes in a single parameter, with the others held fixed. The second consists of the right singular vectors of J(xˆ ) vj, j ) 1, 2, ..., n. (For η(d), the eigenvectors of H(xˆ ) wj, j ) 1, 2, ..., n, play a similar role.) Observe that v1 is the direction in which η˜ (d) is maximal and vn that in which it is minimal. We have
σn e ||J(xˆ ) d|| e σ1 Hence, we find
||f(xˆ )||/10σ1 e η˜ (d) e ||f(xˆ )||/10σn If J(xˆ ) is ill-conditioned, η˜ (d) varies markedly with d. Moreover, the assumption that |d*S(xˆ )d| , ||J(xˆ ) d||2 will usually not be valid when ||J(xˆ ) d|| is small. In the problems at hand, σ1 is usually significantly larger than 1, and σn is usually much smaller than 1 (see below). Since we have n
||J(xˆ ) d||2 )
2 σi2|v* ∑ i d| i)1
we see that, while the bounds above are sharp, for most d, ||J(xˆ ) d|| will lie closer to its upper rather than its lower bound, so η˜ (d) will lie closer to its lower rather than its upper bound. Only for special d such that σi|v*i d| is of order σn for all i will η˜ (d) lie close to its upper bound. Before discussing the implications of the foregoing, we must pause to consider the effect of parameter transformations. Recall that we defined
xj ) ln(1 + kj/sj),
j ) 1, 2, ..., n
The remarks above have been phrased in terms of the transformed variables xj and the quantities
and
e* j g(x) )
∂φ (x) ∂xj
J(x) ej )
∂f (x) ∂xj
e* j H(x) el )
∂2φ (x) ∂xj∂xl
These quantities and the related remarks can be rephrased in terms of the original variables kj by observing that
∂ ∂ ) (kj + sj)-1 ∂kj ∂xj Because of the wide range of values of the kj, expressing things in terms of relative variables ln(kj) is more meaningful, and we have
∂ ∂ ) (kj / (kj + sj)) ∂xj ∂ ln(kj) Note that when sj ≈ kj, the derivatives with respect to the ln(kj) and xj variables are comparable, differing by roughly a factor of 2. Since η˜ (ej) is inversely proportional to ||J(xˆ ) ej||, we can derive the corresponding insensitivity indicators for the original kj or relative ln(kj) variables from those for the transformed variables xj by multiplying by (kj + sj) and (kj + sj)/kj, respectively. This is not true for η˜ (d) in directions other than ej, j ) 1, 2, ..., n, since the effects of all variables are coupled. However, it is a simple matter to calculate the gradient, Jacobian and Hessian in the original or relative variables from those in the transformed variables. In particular, the singular values and singular vectors of the Jacobian, and the eigenvalues and eigenvectors of the Hessian will depend on the variables used, though those for the relative variables will be closely related to those for the transformed variables when sj ≈ kj for j ) 1, 2, ..., n. We have found in the problems at hand that when x ≈ xˆ , choosing sj ≈ kˆ j for j ) 1, 2, ..., n, which necessarily leads to comparable values for the xj, leads also to comparable values for ||J(x) ej|| and hence for e*j H(xˆ ) ej. The condition number of J(xˆ ) and by extension H(xˆ ), is necessarily large if the columns of the Jacobian are of disparate sizes, as is the case for the original variables because of the disparate sizes of the variables. However, the condition numbers may still be large if the Jacobian columns are of comparable size. (See examples below.) When x is not close to xˆ , the choice sj ≈ kj may not yield ||J(x) ej|| of comparable size. To illustrate the strongly nonlinear and ill-conditioned character of the problems at hand, we quote the following results for two particular cases. Case one corresponds to a temperature of 250 °C with m ) 40 and n ) 10. The original kj parameters range from 9.3 × 105 to 1.5 × 10-12, the norms of the columns of the corresponding Jacobian range from 3.6 × 106 to 4.3 × 10-6, and the singular values range from 3.6 × 106 to 3.8 × 10-8, for a condition number of order 1014. The transformed xj variables are of order unity, the norms of the columns of the corresponding Jacobian range from 88.7 to 4.9, and the singular values range from 1.4 × 102 to 2.5 × 10-2, for a condition number of order 104. Note the mollification wrought by the transformation reducing the condition number from order 1014 to order 104, which is still substantial. We achieve χν2 ) 0.21 corresponding to ||f(xˆ )|| ) 2.5 and φ(xˆ ) ) 6.2, with variations in φ of less than 0.01% in the search directions and a scaled gradient whose elements range in magnitude from 2.2 × 10-2 to 5.5 × 10-3. The elements of the approximate Hessian 2J(xˆ )*J(xˆ ) and the
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997 35
Hessian H(xˆ ) are of order 104, while the elements of 2S(xˆ ) are of order 102. The eigenvalues of 2J(xˆ )*J(xˆ ) range from 4.0 × 104 to 6.4 × 10-4 for a condition number of order 108. The eigenvalues of H(xˆ ) range from 4.1 × 104 to -2.1 × 102. Evidently, 2S(xˆ ) is not positive definite and suffices to render H(xˆ ) not positive definite by perturbing the smaller eigenvalues of 2J(xˆ )*J(xˆ ) by amounts comparable to the elements of 2S(xˆ ), as is to be expected. The η˜ (ej) range from 5.1 × 10-2 to 2.8 × 10-3 for the transformed xj parameters, from 5.9 × 104 to 5.9 × 10-14 for the original kj parameters, and from 8.4 × 10-2 to 6.1 × 10-3 for the relative ln(kj) parameters; the η˜ (vj) range from 9.8 to 1.7 × 10-3 for the xj variables. Since H(xˆ ) is not positive definite, it is not meaningful to calculate η(d). Case two corresponds to a temperature of 320 °C with m ) 28 and n ) 11. The original kj parameters range from 1.6 × 1012 to 3.8 × 10-12, the norms of the columns of the corresponding Jacobian range from 5.9 × 106 to 3.9 × 10-12, and the singular values range from 5.9 × 106 to 4.4 × 10-15, for a condition number of order 1021. The transformed xj variables are of order unity, the norms of the columns of the corresponding Jacobian range from 1.1 × 102 to 4.3, and the singular values range from 1.9 × 102 to 5.6 × 10-3, for a condition number of order 104. Note the mollification wrought by the transformation reducing the condition number from order 1021 to order 104, which is still substantial. We achieve χν2 ) 0.58 corresponding to ||f(xˆ )|| ) 3.1 and φ(xˆ ) ) 9.6, with variations in φ of less than 0.01% in the search directions and a scaled gradient whose elements range in magnitude from 4.6 × 10-3 to 1.3 × 10-4. The elements of the approximate Hessian 2J(xˆ )*J(xˆ ) are of order 104, while the elements of 2S(xˆ ) and H(xˆ ) are of order 105. The eigenvalues of 2J(xˆ )*J(xˆ ) range from 7.7 × 104 to 1.2 × 10-4 for a condition number of order 109. The eigenvalues of H(xˆ ) range from 1.9 × 105 to 9.8 for a condition number of order 104. Evidently, 2S(xˆ ) is positive definite and dominates sufficiently to render H(xˆ ) positive definite and better conditioned than 2J(xˆ )*J(xˆ ). The η˜ (ej) range from 7.3 × 10-2 to 2.8 × 10-3 for the transformed xj parameters, from 7.9 × 1010 to 9.0 × 10-14 for the original kj parameters, and from 1.5 × 10-1 to 5.6 × 10-3 for the relative ln(kj) parameters; the η˜ (vj) range from 5.6 × 101 to 1.6 × 10-3 for the xj variables. Since H(xˆ ) is positive definite, η(ej) for the xj variables can be calculated and range from 6.5 × 10-2 to 1.1 × 10-3. In both cases we observe that the problem is highly nonlinear and severely ill-conditioned, increasingly so for the larger temperature. Moving to the transformed xj variables mollifies both aspects. It is sometimes argued that 2S(xˆ ) can be neglected in comparison with 2J(xˆ )*J(xˆ ). In case one, 2S(xˆ ) is 2 orders of magnitude smaller than 2J(xˆ )*J(xˆ ) but still affects d*H(xˆ )d in directions d associated with the small singular values of J(xˆ ). In case two, 2S(xˆ ) is comparable in magnitude to 2J(xˆ )*J(xˆ ) and dominates d*H(xˆ )d. We have identified η˜ (d), or η(d), as a measure of the insensitivity of φ to changes in the parameter vector in direction d, provided the effects of nonlinearity and illconditioning do not confound the interpretation as a measure of uncertainty in the parameter vector. In particular, η˜ (ej)sor η(ej) when applicablesfor j ) 1, 2, ..., n, are analogous to standard deviations derived from the diagonal elements of the covariance matrix in the sense that they reflect the rate of change of φ with respect to the jth parameter, holding the others fixed.
In this sense, we can assign a measure of uncertainty to relative changes in individual parameters by computing η˜ (ej), or η(ej), with respect to the ln(kj) variables. If η˜ (ej) is small, a quantitative interpretation is possible; if η˜ (ej) is large, only qualitative information is provided. Identification of possible outliers in the data requires indices of the sensitivity of xˆ and φ(xˆ ) to changes in data values and to removal of data sets. We shall defer consideration of such indices to another occasion and offer only brief remarks here. We shall term the removal of a putative outlier inconsequential if the resulting decrease in φ(xˆ ) is accompanied by no significant change in xˆ and consequential if xˆ is substantially altered. For example, in Table 2, the removal is inconsequential in the 320 °C case, but substantial changes in xˆ result in the 225 °C case and significant (though less dramatic) changes in xˆ occur in the 250 °C case, so both removals may be regarded as consequential. If a data set has high leverage (for example, because it corresponds to extreme experimental conditions), then the xˆ chosen to fit the data will necessarily attempt to make the residual for that data set relatively small. Consequently, upon removal of that data set, xˆ may change substantially but need not do so if the data set is consistent with the xˆ which fits the rest of the data. Observe in Table 2 that there is a systematic decrease in the insensitivity indicators η˜ (ej) upon removal of one or more candidate outliers. The factor ||f(xˆ )|| necessarily decreases upon such removal, but the factor ||J(xˆ ) d|| could conceivably decrease so much more for some directions d that the ratio, and hence η˜ (d), increases. The sensitivity indicator η˜ (ej) is inherently incapable of providing quantitative information about changes in xˆ j upon removal of data, since φ changes in unpredictable ways, or even of providing a reliable indication whether such a removal might be consequential: compare the results for the 250 and 320 °C cases. The striking behavior in the 225 °C case may be an accidental artifact attributable to the fact that the coordinate direction associated with Kw happens to nearly coincide with the right singular vector of J(xˆ ) associated with the smallest singular value, hence with the maximal value of η˜ (d). However, an anomalously large value of η˜ (ej) for some j, or even of η˜ (vn), may suggest exploration of the possibility of removing outliers on the grounds that such removals may be consequential, though the removal would have to be justified on other grounds. There are, of course, other reasons for large insensitivity indicators which have nothing to do with outliers. In summary, we have not been able to resolve the issue of parameter quality entirely to our satisfaction but have been able to put forward readily computable insensitivity indicators which are a first step. A second step might be to actually compute δˆ ((ej|0.01) and possibly δˆ ((vj|0.01) or to pursue other lines of attack discussed in Bates and Watts (1988), but we have not undertaken these tasks. Conclusions 1. When tert-butyl alcohol undergoes dehydration to isobutylene in hot (T > 220 °C), compressed (P > Psat) liquid water, the water itself is the primary catalytic agent in the dehydration chemistry. At these temperatures the concentration of hydronium ions (protons), that result from the dissociation of water, is sufficient to catalyze the dehydration reaction. One outcome of this finding is that tert-butyl alcohol does not dissociate as an Arrhenius acid, as reported earlier (Xu and Antal,
36
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997
1994). Hydronium ions needed to catalyze the reaction derive from the dissociation of water, not reactant tertbutyl alcohol. 2. The use of kinetic models to discriminate between differing mechanisms and active catalytic agents continues to support the earlier finding (Xu and Antal, 1994) that tert-butyl alcohol acts as a Bronsted acid and thereby catalyzes its own dehydration to isobutylene in hot, liquid water. 3. Over widely ranging temperatures and reactant and catalyst concentrations, the modeling effort also continues to confirm the roles of E1 and AdE2 heterolytic reaction mechanisms in the dehydration chemistry, involving the protonated reactant, its carbocation, protonated ether, and ether as reactive intermediates. Isobutylene formation occurs via elementary unimolecular and bimolecular (with water) decompositions of the unprotonated ether. In the absence of base OH-, isobutylene is not formed directly from the carbocation. 4. As indicated in 3 above, a cascade of reaction pathways leading from a strong acid (H3O+) and base (BuOH) to a weaker acid (BuOH+) and base (H2O) ultimately results in the formation of carbocation (CH3)3C+. The carbocation is so weak an acid that it is unable to protonate water; consequently (in the absence of base OH-) it is not the immediate source of product isobutylene. Thus, the elementary rules of acid-base chemistry rationalize the surprising conclusion that isobutylene forms via the decomposition of ether and not the deprotonation of carbocation. 5. Many other mechanistic models were examined during the course of this work. No others were found to be consistent with the experimental data. 6. In addition to their ability to discriminate between realistic, candidate chemical mechanisms, kinetic studies of the kind illustrated in this paper can provide precise measures of some of the parameters employed in the model. For example, values of Kw determined by the kinetic model used in this work enjoy good agreement with electrochemical measurements. 7. Insensitivity indicators, derived and employed in this paper, offer a measure of the quality of best-fitting parameters that result from the nonlinear least-squares kinetic analysis. These insensitivity indicators also can signal problems within the experimental data set that might otherwise go unnoticed. 8. Without the parameter transformation employed in this paper, the simulation effort would have been numerically intractable. After this transformation was made, all simulations were easily accomplished on an IBM compatible 486 personal computer. 9. Nonlinear least-squares methods customarily neglect the contribution of the S matrix to the Hessian. As indicated in one of the examples herein, it is our experience that the S matrix can make a significant contribution to the Hessian, even near the minimum of the objective function. 10. Future experimental work should focus on the role of di-tert-butyl ether in the reaction mechanism. Kinetic studies of the decomposition of di-tert-butyl ether in hot, liquid water are needed to corroborate our findings. Special attention should be given to the role of elementary unimolecular and bimolecular decompositions of the ether in water. Information concerning the synthesis (Petrov, 1889; Evans and Edlund, 1936) and thermodynamic properties (Fenwick et al., 1975) of this ether is available in the literature.
11. Industry currently employs heterogeneous catalysis to manufacture some isobutylene from tert-butyl alcohol (Abraham and Prescott, 1992). Results reported in this paper suggest that water alone is adequate to catalyze the dehydration chemistry in a commercial setting. Appendix A We review here salient topics from matrix theory (see, for example, Stewart, 1973). If H ∈ Rn×n is symmetric (H* ) H), then H has real eigenvalues λi and associated eigenvectors wi ∈Rn such that Hwi ) λiwi, for i ) 1, 2, ..., n. The eigenvalues can be labeled so that λ1 g λ2 g ... g λn, and the eigenvectors can be chosen to be orthonormal: w*i wj ) δij. An eigenvalue is repeated according to its geometric multiplicity: the number of orthonormal eigenvectors which can be associated therewith. H has a spectral or eigenvalue representation n
H)
λiwiw*i ∑ i)1
The rank and nullity of H are the number of λi which are nonzero and zero, respectively. H is nonsingular if the rank is n, in which case H-1 has eigenvalues λi-1 and associated eigenvectors wi, for i ) 1, 2, ..., n, and the spectral representation n
H-1 )
λi-1wiw* ∑ i i)1
Assume henceforth that we have λ1 > 0. Then, we can classify H as positive definite, positive semidefinite, or indefinite according to λn > 0, λn ) 0, or λn < 0. If H is positive definite, we define its condition number as κ ) λ1/λn. H is ill-conditioned if κ . 1. The significance of ill-conditioning is that H is nearly singular in the sense that a relatively small change in H (for example, omission of the term λnwnw*n in its spectral representation) will produce a singular matrix. Moreover, a relatively small change (for example, replacing the term λnwnw* n by -λnwnw* n) will produce an indefinite matrix. For H ) I, we have λi ) 1, for i ) 1, 2, ..., n; hence, κ ) 1, the minimum possible value. In this case, the eigenvalues can be chosen as any orthonormal set. Define |z| ) xz*z. We have n
x*Hx )
2 λi|w* ∑ i x| i)1 n
x*H-1x )
λi-1|w*i x|2 ∑ i)1
and n
x*x )
2 2 |w* ∑ i x| ) ||x|| i)1
so we obtain, for x * 0,
λn e x*Hx / |x|2 e λ1 and
λ1-1 e x*H-1x / |x|2 e λn-1
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997 37
and these bounds are sharp. We see that
κ ) max (x*Hx/|x|2)/min (x*Hx/|x|2) x*0
x*0
Under our assumption that J has maximal rank, J*J is positive definite and κ(J*J) ) κ(J)2. If J is rank deficient, then J*J is positive semidefinite; if m > n, then JJ* is positive semidefinite. We have
and
n
2
x*0
σi(v* ∑ i x)ui i)1
Jx )
2
κ ) max (x*Hx/|x| ) max (x*H x/|x| ) -1
x*0
n
In particular, we have
J*x )
n
n
j)1
j)1
κ g max (e* j Hej)/min (e* j Hej)
σi(u* ∑ i x)vi i)1
and
and
n
n
n
j)1
j)1
J+x )
κ g max (e*j Hej) max (e*j H-1ej) If Hp ) g, so p ) H-1g, we obtain
so we obtain, for x * 0,
σn e |Jx|/|x| e σ1
n
p)
λi ∑ i)1
-1
(w*i g)wi
and, for g * 0, g*p/|g| |p| ) (g*H-1g/|g|2)(|g|/|p|) > 0. J ∈ Rm×n, m g n, has real singular values, σi, i ) 1, 2, ..., n, which are customarily labeled so that σ1 g σ2 g ... g σn g 0. Associated with these singular values are left singular vectors ui ∈ Rm and right singular vectors vi ∈ Rn, for i ) 1, 2, ..., n, such that Jvi ) σiui and J*ui ) σvi, where u* i uj ) δij and v* i vj ) δij. We see that J*Jvi ) σi2vi and JJ*ui ) σi2ui, so the right singular vectors vi are orthonormal eigenvectors of J*J associated with eigenvalues σi2, and the left singular vectors ui are orthonormal eigenvectors of JJ* associated with eigenvalues σi2. (The other eigenvectors of JJ* correspond to a zero eigenvalue.) J has a singular value representation
σn e |J*x|/|x| e σ1 and
σ1-1 e |J+x|/|x| e σn-1 and these bounds are sharp. We see that
κ ) max (|Jx|/|x|)/min (|Jx|/|x|) x*0
x*0
and that
κ ) |J| |J+| where
|J| ) max (|Jx|/|x|) x*0
n
σiuiv*i ∑ i)1
J)
and
|J+| ) max (|J+x|/|x|)
J* has a singular value representation n
J* )
σiviu* ∑ i i)1
The rank and nullity of J are the number of σi which are nonzero and zero, respectively. Assume hereafter that J has maximal rank, so σn > 0. The pseudoinverse J+ ∈ Rn×m has singular values σi-1 and left and right singular vectors vi and ui; hence, the singular value representation n
J+ )
σi-1(u* ∑ i x)vi i)1
σi-1viu*i ∑ i)1
We define the condition number of J as κ ) σ1/σn, and say that J is ill-conditioned if κ . 1. The significance of ill-conditioning is that J is nearly rank deficient in the sense that a relatively small change in J (for example, omission of the term σnunv*n in its singular value representation) will produce a rank deficient matrix. If m ) n and J is positive definite, the eigenvalues and singular values coincide, and the left and right singular vectors can be chosen to coincide with the eigenvectors, so the condition numbers coincide.
x*0
In particular, we have n
n
j)1
j)1
κ g max |Jej|/min |Jej| and n
n
j)1
j)1
κ g max |Jej| max |J+ej| If g ) J*f and p ) J+f, for f * 0, then p ) 0 iff g ) 0; and if g * 0,
g*p/|g||p| ) (|f|/|g|)(|f|/|p|) > 0 Observe that n
J+J )
viv*i ) I ∑ i)1
but that, for m > n, n
JJ+ )
uiu* ∑ i * I i)1
38
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997
Writing
We have n
Jp ) JJ+f )
x(l+1) ) x(l) - p(l)
(u* ∑ i f)ui i)1
so Jp is the projection of f onto the range of J and is the corresponding projector. We see that
we obtain JJ+
n
J*Jp )
H(x(l)) p(l) ) g(x(l)) or
σi(u* ∑ i f)vi ) J*f i)1
and p is the minimizer of |f - Jp|, or equivalently of |f - Jp|2 ) (f - Jp)*(f - Jp). The residual f - Jp ) (I JJ+)f is orthogonal to the range of J.
p(l) ) H(x(l))-1g(x(l)) Assuming that J(x(l)) has maximal rank, the GaussNewton method generates x(l+1) by minimizing the local quadratic approximation to φ(x) given by
ψ(x|x(l)) ) |f(x|x(l))|2
Appendix B For smooth f : Rn f R, we wish to minimize
where
φ(x) ) f(x)*f(x) ) |f(x)|2
f(x|x(l)) ) f(x(l)) + J(x(l)) (x - x(l))
subject to the nonnegativity restrictions x g 0: that is, xi g 0, for i ) 1, 2, ..., n. A strong local minimizer xˆ g 0 is characterized by the fact that φ(xˆ + d) > φ(xˆ ) for all sufficiently small > 0 and for all directions d such that xˆ + d g 0, with |d| ) 1 for the sake of definiteness. Necessary conditions that xˆ be such a minimizer are that g(xˆ )*d g 0 for all admissible d and d*H(xˆ ) d g 0 for all admissible d such that g(xˆ )*d ) 0. Sufficient conditions are obtained by strengthening the last inequality to d*H(xˆ ) d > 0. These sufficient conditions are equivalent to the assertions (i) that gi(xˆ ) ) 0 for all i such that xˆ i > 0, (ii) that gi(xˆ ) g 0 for all i such that xˆ i ) 0, and (iii) that the submatrix of H(xˆ ) obtained by deleting all rows and columns corresponding to i such that xˆ i ) 0 and gi(xˆ ) > 0 is positive definite. For xˆ > 0, these sufficient conditions are simply that g(xˆ ) ) 0 and that H(xˆ ) is positive definite, the familiar conditions for an unrestricted minimizer. See Gill et al. (1981) for more details. We shall focus initially on unconstrained minimization, tacitly assuming that xˆ > 0 and that we can confine attention to a neighborhood of xˆ in which x > 0. Later we will consider the case where xˆ i ) 0 for some i, or in which we stray from such a neighborhood. The case where xˆ i ) 0 for some i is significant in the problems of interest; it corresponds to a rate constant which cannot be distinguished from zero with the data at hand, permitting the model to be simplified by eliminating the corresponding reaction. It is important to be able to recognize this situation, even though we would ordinarily repeat the calculation with the simplified model, returning to the case xˆ > 0. It is likewise important to avoid difficulties occasioned by accidentally encountering the boundary of the domain, since a negative rate constant is meaningless. The choice of the transformation ln(1 + k/s), rather than the more conventional ln(k/s), is motivated by the desire to admit k ) 0 gracefully. For given x(0), we seek to generate a sequence of approximations x(l), l ) 1, 2, ..., converging to xˆ . Assuming that H(x(l)) is positive definite, the Newton method generates x(l+1) by minimizing the local quadratic approximation to φ(x) given by
is a local affine approximation to f(x). It is easily shown that
φ(x|x(l)) ) φ(x(l)) + g(x(l))*(x - x(l)) + 1 (x - x(l))*H(x(l)) (x - x(l)) 2
φ(x|x(l)) ) ψ(x|x(l)) + (x - x(l))*S(x(l)) (x - x(l)) We see that the Gauss-Newton method can be regarded as a modified version of the Newton method in which the Hessian H(x(l)) ) 2 [J(x(l))*J(x(l)) + S(x(l))] is replaced by the approximate Hessian 2J(x(l))*J(x(l)). Writing
x(l+1) ) x(l) - p(l) again, we obtain
J(x(l))*J(x(l))p(l) ) J(x(l))*f(x(l)) or
p(l) ) (J(x(l))*J(x(l)))-1J(x(l))*f(x(l)) This is equivalent to solving the overdetermined linear system
J(x(l)) p(l) ) f(x(l)) in the least-squares sense, to obtain
p(l) ) J(x(l))+f(x(l)) (See Appendix A.) Given the difficulty and expense involved in estimating S(x(l)) by differencing, we have used the Gauss-Newton method to generate p(l), except in special circumstances. The matter will be examined further below. If the assumption that H(x(l)) is positive definite is invalid or if the condition number of H(x(l)) is too large, we replace the Newton method by the regularized Newton method obtained by minimizing
1 φ˜ (x|x(l)) ) φ(x|x(l)) + τl2|x - x(l)|2 2 which results in solving
H ˜ (x(l)) p˜ (l) ) g(x(l)) where
H ˜ (x(l)) ) H(x(l)) + τl2I
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997 39
for
p˜ (l) ) H ˜ (x(l))-1g(x(l)) The parameter τl > 0 is chosen large enough so that H ˜ (x(l)) is positive definite and has a condition number equal to a specified upper bound. If the assumption that J(x(l)) has maximal rank is invalid or if the condition number of J(x(l)) is too large, we replace the GaussNewton method by the regularized Gauss-Newton method obtained by minimizing
ψ ˜ (x|x(l)) ) ψ(x|x(l)) + τl2|x - x(l)|2 which results in
p˜ (l) ) J ˜ (x(l))+˜f(x(l)) where
J ˜ (x(l)) )
[ ]
˜f(x(l)) )
[ ]
J(x(l)) τ lI
and
f(x(l)) 0
This is equivalent to solving
J ˜ (x(l))*J ˜ (x(l)) p˜ (l) ) J ˜ (x(l))*f˜(x(l)) or
(J(x(l))*J(x(l)) + τl2I)p˜ (l) ) J(x(l))*f(x(l)) The parameter τl > 0 is chosen large enough so that J˜ (x(l)) has maximal rank and has a condition number equal to a specified upper bound. Assuming that g(x(l)) * 0, the rate of decrease of φ as we move away from x(l) in direction d(l) is maximal for d(l) ) -g(x(l))/|g(x(l))|, the direction of steepest descent. Assuming that J(x(l)) has maximal rank, the Cauchy method chooses x(l+1) ) x(l) - q(l) as the minimizer of ψ(x|x(l)) as we move away from x(l) along the ray with direction d(l). We obtain
q(l) ) [ |g(x(l))|2 / |J(x(l))g(x(l))|2 ]g(x(l)) A regularized Cauchy method could be obtained by replacing ψ(x|x(l)) by ψ ˜ (x|x(l)). In all of the foregoing methods, we have x(l+1) ) x(l) r(l), where r(l) satisfies a linear equation Ar(l) ) g(x(l)) for a positive definite matrix A characteristic of the method. We have (l)
(l)
(l)
(l)
r *g(x ) ) r *Ar > 0 for r(l) * 0. Therefore, -r(l)/|r(l)| is a descent direction in which φ decreases locally as we move away from x(l) a short distance. This does not guarantee, however, that φ(x(l+1)) < φ(x(l)). In the damped version of each method, one takes instead x(l+1) ) x(l) - Rlr(l), where Rl approximates the smallest positive minimizer of φ(x(l) - Rr(l)). We know that φ(x(l) - Rr(l)) < φ(x(l)) for r(l) * 0 and R > 0 sufficiently small, so we can choose Rl so that we satisfy the descent condition φ(x(l) - Rlr(l)) < φ(x(l)). The process of choosing Rl is called a line search. We choose to scan along the ray x(l) - Rr(l), seeking three
values of R which bracket the minimizer; and then take Rl as the minimizer of the quadratic interpolating polynomial approximating φ(x(l) - Rr(l)) defined by these three values of R. The motivation for introducing the damping parameter Rl suggests the restriction 0 < Rl e 1, with R ) 0 and R ) 1 being two candidate sample values for the interpolation; but when φ(x(l) - r(l)) < φ(x(l)), we can profit by permitting Rl > 1. Though r(l) has been constructed so that the choice R ) 1 is natural, the line search rarely chooses Rl ) 1; indeed, Rl is commonly significantly smaller or larger than 1. Let p(l) be the step suggested by the Newton or GaussNewton method, regularized as necessary, and q(l) that for the Cauchy method. We choose to conduct line searches for r(l) ) p(l) and r(l) ) q(l) and choose the method which results in the greatest decrease in φ at that stage. An insufficient decrease in both methods signals termination of the iteration. The line searches involve only a fraction of the cost of evaluating J(x(l)) and even less when evaluating H(x(l)). More elaborate line search processes, more sophisticated combinations of Cauchy with Newton or GaussNewton methods, and more elegant regularization procedures may be found in Dennis and Schnabel (1983). We have found the simpler method outlined above both adequate and more convenient for our purposes than the IMSL codes embodying some of these more compllicated techniques. The ability to monitor and control the course of the calculation is valuable for nonroutine minimization problems. The nonnegativity restrictions can now be incorporated as follows: If x(l) > 0 and r(l) e 0, the restrictions are not operative, and x(l+1) may be determined as discussed above. If x(l) > 0 and there is at least one i such that ri(l) > 0, let βl be the smallest of the ratios xi(l)/ri(l) for all i such that ri(l) > 0. We then confine the search for Rl to the interval (0, βl]. If no minimizer of φ(x(l) - Rr(l)) is detected in (0, βl), we take Rl ) βl and we will have xi(l+1) ) 0 for at least one i. If xi(l) ) 0 for at least one i, let S(l) be the set of all i such that xi(l) ) 0 and gi(x(l)) > 0. If S(l) is empty, x(l+1) may be determined as discussed above. If S(l) is nonempty, we consider the reduced problem of minimizing φ subject to the restriction xi(l+1) ) 0 for i ∈ S(l). The procedure given previously can be applied to the reduced problem to determine xi(l+1) for i ∉ S(l). Whether S(l) is empty or nonempty, the fact that r(l) is a descent direction guarantees that x(l+1) g 0 and φ(x(l+1)) < φ(x(l)). Practical termination criteria must be based on quantities like the relative rate of decrease of the objective function (φ(x(l)) - φ(x(l+1)))/φ(x(l)) and the scaled gradient elements
g˜ i(x(l)) ) (J(x(l))ei)*f(x(l)) / |J(x(l))ei||f(x(l))| for i ) 1, 2, ..., n. We see that g˜ i(x(l)) is the cosine of the angle between f(x(l)) and J(x(l))ei and should be close to zero for all i such that xi(l) > 0 and not significantly less than zero for all i such that xi(l) ) 0. To provide intuition as to the potential effects of nonlinearity and ill-conditioning on the iteration and the parameter sensitivity, envision the shape of the objective function surface φ(x) in the vicinity of xˆ as approximated by
ψ(x|xˆ ) ) |f(x|xˆ )|2 ) |f(xˆ ) + J(xˆ ) (x - xˆ )|2 for n ) 2. The description involves the singular values σ1 and σ2 and corresponding right singular vectors v1
40
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997
and v2 of J(xˆ ), with κ ) σ1/σ2 . 1. We have an elongated bowl which is virtually a steep-sided trough. The alongaxis direction of the “trough” is aligned with v2, and the cross-axis direction, with v1. The rates of change of ψ(x|xˆ ), and by extension φ, in the along-axis and crossaxis directions are proportional to σ2 and σ1, respectively. By virtue of nonlinearity, the trough may be expected to “curve” as one moves away from xˆ in the along-axis direction. For n . 2 and κ ) σ1/σn . 1, one can proceed by analogy in what follows by identifying the along-axis direction with the subspace spanned by the right singular vectors corresponding to singular values less than or equal to xσ1σn and the cross-axis direction with the subspace spanned by the right singular vectors corresponding to singular values greater than xσ1σn. Any vector can be expressed as the sum of orthogonal cross-axis and along-axis components. Separation into components based on xσ1σn is arbitrary, of course, but a plausible heuristic. Because of nonlinearity, any such picture is a caricature except in the immediate vicinity of xˆ , possibly marginally improved by replacing ψ(x|xˆ ) with φ(x|xˆ ) and replacing the singular values σi of J(xˆ ) by the square roots of the eigenvalues λi of H(xˆ ) and the right singular vectors vi with the eigenvectors wi. During the course of the iteration, ψ(x|x(l)) and J(x(l)), or φ(x|x(l)) and H(x(l)), play analogous roles. Let us examine salient features of the iteration and parameter sensitivity in light of this geometric picture. For the sake of brevity, we shall refer to the Newton and Gauss-Newton methods, and their regularized counterparts, by invoking the p(l) search direction and to the Cauchy method by invoking the q(l) search direction. Unless x(l) lies near the bottom of the trough, the gradient and hence both p(l) and q(l) will be dominated by their cross-axis components. When x(l) lies closer to the bottom, the gradient will have a significant along-axis component which will be reflected in q(l) and in strongly magnified form in p(l), to the extent that p(l) and q(l) may be nearly orthogonal. When x(l) lies very near the bottom, the gradient may be small, with comparable cross-axis and along-axis components or even a dominant (though small) along-axis component. Without regularization or damping an overly aggressive excursion along the trough in the p(l) search direction may result. Regularization restricts the along-axis component of p(l), while damping restricts the extent of motion in the search direction, though in favorable circumstances one finds Rl > 1 so the extent of motion is actually accentuated. We can anticipate rapid reduction of φ in the early stages of the iteration, slowing to a crawl in later stages, especially if significant motion along the trough is required to reach xˆ . Since we must terminate the iteration based on the relative rate of reduction of φ and the size of the scaled gradient, we cannot exclude the possibility of apparent convergence even when x(l) is some distance from xˆ , but hopefully φ(x(l)) is close to φ(xˆ ). Moreover, such minimization methods are inherently capable only of finding local minimizers; we cannot exclude the possibility that the global minimizer lies in some unprobed region of parameter space. We are confident, based on our monitoring and tuning of the iterative process, that the final value of φ is close to φ(xˆ ), and the excellent values of χν2 suggest that φ(xˆ ) is indeed a global minimum. However, we have encountered situations wherein misleading results could emerge in the absence of careful monitoring and tuning.
Typically the directions ej, j ) 1, 2, ..., n, all have substantial cross-axis components, so the η˜ (ej) are small and comparable in size, for the transformed or relative variables. For the original variables, scaling considerations intrude. There are compensation effects when coordinated changes of several parameters are considered. For directions in parameter space whose alongaxis component is dominate, φ is highly insensitive to changes in the parameter vector. Nonlinearity dictates that such directions are simply tangential to a curved surface in parameter space along which φ varies slowly. Thus, we cannot pin down the parameters yielding a good fit entirely, because dramatic compensation effects are possible in the presence of ill-conditioning and nonlinearity. However, it should be appreciated that such directions in parameter space are very special for large n. The parameters are well-determined with respect to changes in most directionssthe case n ) 2 is potentially misleading in this regard. Acknowledgment This work was supported by the National Science Foundation under Grants CBT88-12954 and BCS9111743 and the Coral Industries Endowment. We thank Professor Maitland Jones, Jr. (Princeton University), and Dr. Wilbur Hurst (NIST) for helpful discussions and advice and Dr. Fred Heineken and Dr. Maria Burka (NSF) for their continuing interest in this work. We also thank three anonymous reviewers for their helpful comments. Literature Cited Abraham, O. C.; Prescott, G. F. Make isobutene from TBA. Hydrocarbon Proc. 1992, Feb, 51-54. Antal, M. J.; Anderson, D. G. M. Inverse Problem Solver, Department of Mechanical Engineering, University of Hawaii at Manoa, Honolulu, HI, 1994. Antal, M. J.; Brittain, A.; DeAlmeida, C. Ramayya, S.; Roy, J. C. Heterolysis and Homolysis in Supercritical Water. In Supercritical Fluids; Squires, T. G., Paulaitis, M., Eds.; ACS Symposium Series 329; American Chemical Society: Washington, DC, 1987; pp 77-87. Antal, M. J.; Mok, W. S. L.; Richards, G. N. Mechanism of Formation of 5-(Hydroxymethyl)-2-Furaldehyde from D-Fructose and Sucrose. Carbohydr. Res. 1990a, 199, 91-109. Antal, M. J.; Mok, W. S. L.; Richards, G. N. Four-Carbon Model Compounds for the Reactions of Sugars in Water at High Temperature. Carbohydr. Res. 1990b, 199, 111-115. Antal, M. J.; Leesomboon, T. C.; Mok, W. S. L.; Richards, G. N. Mechanism and Kinetics of Formation of 2-Furaldehyde from D-Xylose. Carbohydr. Res. 1991, 217, 71-85. Bard, Y. Nonlinear Parameter Estimation, Academic Press: New York, 1974. Bates, D. M.; Watts, D. G. Nonlinear Regression Analysis and Its Applications; John Wiley: New York, 1988. Carlsson, M.; Habenicht, C.; Kam, L. C.; Antal, M. J.; Bian, N.; Cunningham, R. J.; Jones, M. Study of the Sequential Conversion of Citric to Itaconic to Methacrylic Acid in Near-Critical and Supercritical Water. Ind. Eng. Chem. Res. 1994, 33, 19891996. Cutler, A. H.; Antal, M. J.; Jones, M. A Critical Evaluation of the Plug Flow Idealization of Tubular-Flow Reactor Data. Ind. Eng. Chem. Res. 1988, 27, 691-697. Dennis, J. E. Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; PrenticeHall: Englewood Cliffs, NJ, 1983. Ellis, A. J.; Mahon, W. A. J. Chemistry and Geothermal Systems; Academic Press: New York, 1977; p 128. Evans, T. W.; Edlund, K. R. Tertiary Alkyl Ethers Preparation and Properties. Ind. Eng. Chem. 1936, 28, 1186-1188.
Ind. Eng. Chem. Res., Vol. 36, No. 1, 1997 41 Fenwick, J. O.; Harrop, D.; Head, A. J. Thermodynamic Properties of Organic Oxygen Compounds. J. Chem. Thermodyn. 1975, 7, 943-954. Gill, P. E.; Murray, W.; Wright, M. H. Practical Optimization; Academic Press: New York, 1981. Habenicht, C.; Kam, L. C.; Wilschut, M. J.; Antal, M. J. Homogeneous Catalysis of Ethyl tert-Butyl Ether Formation from tertButyl Alcohol in Hot, Compressed Liquid Ethanol. Ind. Eng. Chem. Res. 1995, 34, 3784-3792. Hill, C. G., Jr. An Introduction to Chemical Engineering Kinetics and Reactor Design; John Wiley & Sons: New York, 1977. Lowry, T. H.; Richardson, K. S. Mechanism and Theory in Organic Chemistry, 3rd ed.; Harper & Row: New York, 1987. Marshall, W. L.; Franck, E. U. Ion Product of Water Substance, 0-1000 °C, 1-1000 bars; New International Formulation and Its Background. J. Phys. Chem. Ref. Data 1981, 10, 295-304. Mok, W. S. L.; Antal, M. J.; Jones, M. The Formation of Acrylic Acid from Lactic Acid in Supercritical Water. J. Org. Chem. 1989, 54, 4596-4602. Myers, R. H. Classical and Modern Regression with Applications, 2nd ed.; PWS-Kent: Boston, 1990. Narayan, R.; Antal, M. J. A Kinetic Elucidation of the AcidCatalyzed Mechanism of 1-Propanol Dehydration in Supercritical Water. In Supercritical Science and Technology; Johnston, K. P., Penninger, J. M. L., Eds.; ACS Symposium Series 406; American Chemical Society: Washington, DC, 1989; pp 226241. Narayan, R.; Antal, M. J. Influence of Pressure on the AcidCatalyzed Rate Constant for 1-Propanol Dehydration in Supercritical Water. J. Am. Chem. Soc. 1990, 112, 1927-1931. Ochel, H. Private communication, 1994. Ochel, H.; Bean, V. E.; Bowers, W. J.; Davis, R. W.; Hurst, W. S. Optical Probing of Hydrothermal Reactions and Flow Patterns in a High-Pressure Flow Reactor. Proceedings of the 12th International Conference on the Properties of Water and Steam, Orlando, FL, Sept 11-16, 1994. Ogi, T.; Yokoyama, S.; Minowa, T.; Dote, Y. Role of Butanol Solvent in Direct Liquefaction of Wood. J. Jpn. Petrol. Inst. 1990, 33, 383-389.
Petrov, J. J. Russ. Phys.-Chem. Ges. 1889, 21, 348. Ramayya, S.; Antal, M. J. An Evaluation of Systematic Error Incured in the Plug Flow Idealization of Laminar Flow Reactor Data. Energy Fuels 1989, 3, 105-108. Ramayya, S.; Brittain, A.; DeAlmeida, C.; Mok, W. S. L.; Antal, M. J. Acid Catalyzed Dehydration of Alcohols in Supercritical Water. Fuel 1987, 66, 1364-1371. Schneider, G. M. Phase Behavior of Aqueous Solutions at High Pressures. In WatersA Comprehensive Treatise; Franks, F., Ed.; Plenum Press: New York, 1973; Vol. 2. Stewart, G. W. Introduction to Matrix Computations; Academic Press: New York, 1973. Xu, X. Autocatalytic Dehydration of tert-Butanol in Near-Critical Water. Ph.D. Thesis, University of Hawaii, Honolulu, HI, 1992. Xu, X.; Antal, M. J. Kinetics and Mechanism of Isobutene Formation from t-Butanol in Hot Liquid Water. AIChE J. 1994, 40, 1524-1534. Xu, X.; DeAlmeida, C.; Antal, M. J. Mechanism and Kinetics of the Acid-Catalyzed Dehydration of Ethanol in Supercritical Water. J. Supercrit. Fluids 1990, 3, 228-232. Xu, X.; DeAlmeida, C.; Antal, M. J. Mechanism and Kinetics of the Acid-Catalyzed Formation of Ethene and Diethyl Ether from Ethanol in Supercritical Water. Ind. Eng. Chem. Res. 1991, 30, 1478-1485.
Received for review June 17, 1996 Revised manuscript received September 6, 1996 Accepted September 6, 1996X IE960349O
Abstract published in Advance ACS Abstracts, November 1, 1996. X
