Solving Kinetic Inversion Problems via a Physically Bounded Gauss

Apr 5, 2005 - An iterative physically bounded Gauss-Newton (PGN) method has been ... solving the Gauss-Newton update equations within the physical ...
0 downloads 0 Views 260KB Size
Download PDF
3626

Ind. Eng. Chem. Res. 2005, 44, 3626-3637

Solving Kinetic Inversion Problems via a Physically Bounded Gauss-Newton (PGN) Method Weiyong Tang, Libin Zhang, Andreas A. Linninger,* Robert S. Tranter,† and Kenneth Brezinsky Laboratory for Product and Process Design, Department of Chemical Engineering, University of Illinois at Chicago, Chicago, Illinois 60607

An iterative physically bounded Gauss-Newton (PGN) method has been formulated to estimate unknown kinetic parameters from experimental measurements. A physically bounded approach is adopted to reduce the size of the search space and ensure search within physically meaningful ranges of kinetic rates. First-order sensitivity information of state variables, with respect to unknown rate parameters, is computed simultaneously with the integration of the governing ordinary differential equations (ODEs). Optimal kinetic parameters are obtained by iteratively solving the Gauss-Newton update equations within the physical parameter range. Four different reaction systems, including both simple and complex networks, have been utilized for validating the performance of the proposed method. Comparisons with other state-of-the-art algorithms showed its efficiency, robustness, and accuracy. The methodology was also applied to analyze ethane oxidation based on the widely cited GRI-Mech 3.0 mechanism to demonstrate the algorithm’s performance in large-scale practical applications. The model predictions have been greatly improved by determining the optimal pre-exponential factors for five unimolecular reactions from experimental data obtained at elevated pressures. 1. Introduction Rate parameter information in kinetic expressions is essential for the design, optimization, and control of chemical reaction processes. The reaction rates are generally obtained by fitting model equations against time series data measured through chemical kinetic studies in a variety of techniques, such as absorption spectroscopy and laser photolysis.1 The problem of determining kinetic parameters from experiments is referred to as a kinetic inversion problem (KIP). The review of existing techniques for KIP shows that the most common methods are gradient-based leastsquares methods.2-5 Conventionally, the least-squares methods iterate toward the solution by means of successive evaluations of the error function and its gradients obtained analytically or numerically. As an example, Tadi and co-workers4,5 proposed a novel gradientbased algorithm to ensure convergence by treating the unknown parameters as time-dependent variables. An alternative to forward integration of the kinetic equations lies in global finite-element polynomial collocation for the discretization of the concentration profiles.6 Very recently, Arora and Biegler developed a trust region successive quadratic programming (SQP) method for equality-constrained parameter estimation with successful application to a polymerization reactor.7,8 In addition to the deterministic methods mentioned previously, researchers also have deployed various stochastic search techniques, including simulated annealing (SA),9,10 genetic algorithms (GA),11 and Monte Carlo (MC) simulations.12 However, these techniques with a random element require frequent evaluations of the governing kinetic equations and, hence, are com* To whom correspondence should be addressed. Phone: (312) 996-2581. Fax: (312) 996-0808. E-mail: [email protected]. † Current address: Chemistry Division, Argonne National Laboratory, Argonne, IL.

putationally extremely demanding. Furthermore, the results of stochastic methods are heavily dependent on tuning parameters such as annealing schedule (SA), population size and probability of mutation (GA), and sampling strategy (MC). The “un-robust” character of these “robust” methods can cause large errors in rate recovery, as will be demonstrated in this article. An alternative approach avoiding prohibitively large computational cost is to parametrize the time integration of chemical kinetic equations with a set of algebraic polynomial formulas called the response surface or map. Various names have been proposed for this technique, such as solution mapping13,14 and high-dimensional model representation,15 although their basic mathematical ideas are similar. The key to these methods is to find an algebraic map capable of representing the composition trajectory for all possible kinetic rates. Unfortunately, it is generally difficult and computationally expensive to extract such a map from symbolic information.15 The open literature reveals that there are several numerical methods, both local and global, devoted to studying reaction networks, but usually with 99%). 3. Applications The proposed PGN method was successfully applied to determine rate constants for both simple and complex reaction networks. Its efficiency and robustness are validated by four different reaction mechanisms. Its capability of handling complex reaction networks is demonstrated by the application to a practical largescale problem (e.g., ethane oxidation). Five optimal parameters were determined, leading to an improved GRI-Mech 3.0 kinetic model that is in better agreement with measured data. We found that our new approach was successful in rendering kinetic reaction information from experimental data, even in this large-scale problem. Note that, in all application studies in this work, the physical bounds were not chosen as a “tuning” parameter. It is known that classical transition-state theory can provide an order-of-magnitude accuracy in estimating the pre-exponential factor, of which estimates are derived from vibrational frequencies that correspond to the activated complexes and the reactants.1 Consequently, the physical bounds of a rate constant, [k0/10, k0 × 10], was used in most cases in the present work, including the practical problem of ethane oxidation, and is recommended for future application of the PGN. To test the robustness of the algorithm, wider bounds also were used. 3.1. Simple Chemical Reaction Systems. Firstorder reversible chemical reaction mechanisms are simple but occur in real processes, such as hydrocarbon cracking and the formation of the first aromatic ring in sooting flames.32 Two simple reaction systems were utilized to compare the performance of the present algorithm with other methods.

3630

Ind. Eng. Chem. Res., Vol. 44, No. 10, 2005

3.1.1. Three-Species Reversible Isomerization. The KIP of three-species isomerization has been the subject of previous studies.16,17 The first test of the PGN method with a linear kinetic system was performed on a three-species reversible reaction network that is given in Figure 1. True values of components of the rate constant vector are k0 ) [4.623 10.344 3.371 5.616 3.724 1.0]T. The boundary conditions were set as ki ∈ [0.01, 100]. Error-free experimental data were synthesized by numerical simulation. The concentrations of species A1 and A2 at a time interval of 0.05 s, computed from k0, are shown in Figure 2. The KIP that is intended to identify all six rate constants from 14 perfect “experimental” data points. For comparison with a previous study by Tadi and Cai,5 the same conditions were used in the current work with initial species concentrations as y0 ) [0.04 0.96 0.0]T and the initial guess vector as k0 ) [3.8979 4.0796 3.9987 4.0028 4.0810 3.9222]T. Mathematically, the KIP in this case study can be formulated as system (16) (the system defined by eqs 16). 7

min F(k) ) k

∑ i)1

[(

)

cal yexp A (ti) - yA (ti)

yexp A (ti)

(

Figure 1. Simple reaction network involving three-species isomerization.

2

+

)]

cal yexp B (ti) - yB (ti)

yexp B (ti)

2

(16)

s. t.

dycal A cal cal cal ) -k12ycal A + k21yB - k13yA + k31yC dt dycal B cal cal cal ) -k21ycal B + k12yA - k23yB + k32yC dt dycal C cal cal cal ) -k31ycal C + k13yA - k32yC + k23yB dt ycal A (t0 ) 0) ) 0.04 ycal B (t0 ) 0) ) 0.96 ycal C (t0 ) 0) ) 0.00 0.01 e k12 e 100 0.01 e k21 e 100 0.01 e k23 e 100 0.01 e k32 e 100 0.01 e k31 e 100 0.01 e k13 e 100 Using the PGN approach, the rate constants were identified in just eight iterations, as detailed and shown in Table 1. As mentioned earlier in the introduction, the KIP can also be approached using a traditional NLP code hooked to an ODE integrator that produces sensitivity information. For comparison, we also solved system (2) directly with a secant method (L-BFGS-B). In that test, the cost function value and its gradients, F(k) and ∂F/∂kj, with

Figure 2. Trajectories of state variables (solid lines) and errorfree experimental data points (symbols) for the simple reaction network in Figure 1.

respect to unknown rate constants, ki, must be encoded. However, the derivatives of the cost function, with respect to the unknown kinetic rates, cannot be computed explicitly. Therefore, the necessary numerical differentiation requires complete reintegration for each partial derivative in each iteration. When solving system (2) directly for the KIP expressed in system (16), convergence was achieved after 60 iterations, including 360 reintegration steps. Table 2 presents the slow convergence of the computational results of the direct solution approach. The proposed PGN method, by contrast, required only eight iterations, with one reintegration per iteration. The same KIP was also studied by Tadi and Cai,5 using their time-dependent treatment of unknown parameters. The number of iterations needed for convergence was extremely large. Even after 800 000 iterations, there still was an average deviation of ∼1% from true values. We also tested the sensitivity of the PGN approach to experimental errors. More realistic “experimental data” were generated by artificially introducing a 5% random error into the true concentrations. Despite errors, we found that the rate parameters were successfully identified with a very high accuracy (>99%). Errors in the experimental data may also complicate the KIP by introducing multiple local minima. Therefore, we wanted to validate the results of the PGN with global solvers. Note that the comparison between a global optimizer such as BARON and our local method is presented for practical purposes, only to report on the expected computational performance from the kinetic practitioners’ point of view. In this particular case, the global optimization problem consisted of a dynamic optimization problem with ODE constraints discretized by orthogonal collocation on finite elements (OCFE).33

Ind. Eng. Chem. Res., Vol. 44, No. 10, 2005 3631 Table 1. Convergence History To Determine Six Unknown Rate Constants in a Three-Species Reversible Isomerization Network from 12 Error-Free Data, Using the PGN Approach iteration

F(k)

k12

k21

k23

k32

k31

k13

0 1 2 3 4 5 6 7 8

2.29 × 100 1.24 × 100 3.21 × 10-2 2.18 × 10-3 6.67 × 10-4 6.80 × 10-6 2.20 × 10-9 1.13 × 10-12 2.05 × 10-16

3.8979 0.0100 0.7690 2.4716 6.3075 4.9192 4.6273 4.6229 4.6230

4.0796 9.5291 9.9747 10.2869 10.3711 10.3460 10.3438 10.3440 10.3440

3.9987 2.7791 3.6717 3.2521 3.4051 3.3710 3.3706 3.3710 3.3710

4.0028 10.3884 19.3415 12.1046 0.0100 4.5898 5.5999 5.6163 5.6160

4.0810 0.0100 0.0100 1.0959 6.5051 4.1033 3.7142 3.7234 3.7240

3.9222 1.3108 3.0166 2.3127 0.1456 0.8156 0.9929 0.9999 1.0000

4.623

10.344

3.371

5.616

3.724

1.0

actual value

Table 2. Convergence History To Determine Six Unknown Rate Constants in a Three-Species Reversible Isomerization Network from 12 Error-Free Data, by Solving System (2) Directly with L-BFGS-B iteration 0 1 2 3 l 17 18 19 20 l 37 38 39 40 l 57 58 59 60

F(k)

k12

k21

k23

k32

k31

k13

2.29 × 100 1.08 × 100 4.57 × 10-1 2.86 × 10-1 l 2.75 × 10-6 2.75 × 10-6 2.75 × 10-6 2.75 × 10-6 l 1.62 × 10-9 1.62 × 10-9 1.62 × 10-9 1.61 × 10-9 l 3.00 × 10-12 7.43 × 10-13 3.53 × 10-14 1.06 × 10-15

3.8979 3.3842 2.6485 2.3669 l 4.7062 4.7062 4.7063 4.7063 l 4.6104 4.6104 4.6104 4.6104 l 4.6229 4.6230 4.6230 4.6230

4.0796 5.0530 6.1123 6.7937 l 10.3344 10.3344 10.3344 10.3343 l 10.3441 10.3441 10.3441 10.3441 l 10.3440 10.3440 10.3440 10.3440

3.9987 4.5424 4.7591 4.6135 l 3.3777 3.3778 3.3779 3.3780 l 3.3715 3.3715 3.3715 3.3715 l 3.3710 3.3710 3.3710 3.3710

4.0028 3.7175 3.5909 3.8104 l 5.3118 5.3118 5.3118 5.3119 l 5.6628 5.6628 5.6628 5.6627 l 5.6163 5.6159 5.6162 5.6160

4.0810 4.2713 4.7265 5.1452 l 6.2128 6.2128 6.2126 6.2122 l 3.7329 3.7329 3.7328 3.7328 l 3.7224 3.7235 3.7241 3.7240

3.9222 3.7150 3.1533 2.5322 l 1.6055 1.6055 1.6054 1.6053 l 1.0152 1.0152 1.0152 1.0151 l 0.9996 0.9998 1.0001 1.0000

4.623

10.344

3.371

5.616

3.724

1.0

actual value

The fully OCFE-discretized differential equations converted the dynamic optimization problem to a set of nonlinear equations, in terms of the unknown concentration and the unknown kinetic rate constants (parameters). The nonlinear optimization problem was then solved by the global optimal solver “BARON” that is available in GAMS.34 With the same physical constraints and initial guesses, the same solution with the PGN approach was attained after 24 min of CPU time. 3.1.2. Four-Species Reversible Isomerization. Figure 3 describes another first-order kinetic network that involves four species in a reversible isomerization system. The same conditions as those in a previous paper5 were used: initial concentrations, y0 ) [0.4 0.2 0.1 0.3]T; true rate constant values, k0 ) [1.0 3.0 5.0 10.0 2.0 4.0 10.0 8.0]T; and initial guesses, k0 ) [19.72 0.781 0.516 29.38 11.88 1.328 2.003 7.934]T. Physical boundary conditions were set as ki ∈ [ki,0/20,ki,0 × 20]. Error-free data were also synthesized for time intervals

Figure 3. Simple reaction network involving four-species isomerization.

of 0.1 s. Concentration trajectories of all four species are shown in Figure 4. In this case, the KIP involved 36 observed data and 8 unknown rate constants. Table 3 displays the convergence history of the PGN approach to identify all 8 unknown rate constants in the aforementioned KIP. The true solution was attained after 30 iterations. Very good accuracy was already achieved after 20 iterations. In contrast, a local minimum was found with L-BFGS-B applied directly to system (2) after more than 70 iterations, as shown in Table 4. On the other hand, the time-dependent treatment of unknown parameters5 required 40 000 iterations to achieve convergence. Even at that point, the result was less accurate than the solution reported here.

Figure 4. Trajectories of state variables (solid lines) and errorfree experimental data points (symbols) for the simple reaction network in Figure 3.

3632

Ind. Eng. Chem. Res., Vol. 44, No. 10, 2005

Table 3. Convergence History To Determine Eight Unknown Rate Constants in a Four-Species Reversible Isomerization Network from 36 Error-Free Data, Using the PGN Approach iteration 0 1 2 3 l 17 18 19 20 l 27 28 29 30

F(k)

k12

k21

k23

k32

k34

k43

k41

k14

4.48 × 101 5.89 × 101 2.46 × 101 7.53 × 100 l 1.79 × 10-9 7.41 × 10-10 2.93 × 10-9 1.18 × 10-10 l 1.96 × 10-13 7.86 × 10-14 3.15 × 10-14 1.26 × 10-14

19.7200 0.0500 15.6290 1.4426 l 1.0000 1.0000 1.0000 1.0000 l 1.0000 1.0000 1.0000 1.0000

0.7810 6.1267 22.0580 0.1500 l 3.0016 2.9990 3.0006 2.9996 l 3.0000 3.0000 3.0000 3.0000

0.5160 6.9844 7.4023 5.8912 l 5.0005 4.9997 5.0002 4.9999 l 5.0000 5.0000 5.0000 5.0000

29.3800 0.5000 3.3821 5.6891 l 10.0031 9.9980 10.0013 9.9992 l 10.0000 10.0000 10.0000 10.0000

11.8800 0.1000 0.1000 0.1231 l 1.9980 2.0013 1.9992 2.0005 l 2.0000 2.0000 2.0000 2.0000

1.3280 2.7041 0.5164 0.2000 l 4.0003 3.9998 4.0001 3.9999 l 4.0000 4.0000 4.0000 4.0000

2.0030 20.8447 36.7336 104.3499 l 10.0133 9.9914 10.0054 9.9966 l 10.0001 9.9999 10.0001 10.0000

7.9340 46.1909 0.4000 24.9351 l 8.0100 7.9936 8.0040 7.9974 l 8.0001 7.9999 8.0000 8.0000

5.0

10.0

actual value

1.0

3.0

2.0

4.0

10.0

8.0

Table 4. Convergence History To Determine Eight Unknown Rate Constants in a Four-Species Reversible Isomerization Network from 36 Error-Free Data by Solving System (2) Directly with L-BFGS-Ba iteration 0 1 2 3 l 17 18 19 20 l 37 38 39 40 l 67 68 69 70

F(k)

k12

k21

k23

k32

k34

k43

k41

k14

4.48 × 101 7.60 × 101 2.46 × 101 1.83 × 101 l 8.98 × 10-2 9.75 × 10-2 7.57 × 10-2 4.05 × 10-2 l 6.63 × 10-4 6.57 × 10-4 6.53 × 10-4 6.51 × 10-4 l 3.98 × 10-2 3.98 × 10-2 3.98 × 10-2 3.98 × 10-2

19.7200 18.9219 19.4644 19.4379 l 6.5984 4.0597 5.3369 4.1005 l 0.5523 0.5032 0.4888 0.4925 l 0.2498 0.2498 0.2498 0.2498

0.7810 11.2726 4.1408 4.6044 l 9.4436 5.2722 7.3707 5.1532 l 2.6853 2.6603 2.6538 2.6575 l 3.2506 3.2506 3.2506 3.2506

0.5160 10.6492 3.7610 4.2024 l 16.0800 16.0062 16.0433 16.0314 l 14.8913 14.8648 14.8581 14.8599 l 14.1916 14.1916 14.1916 14.1916

29.3800 29.1040 29.2916 29.2857 l 24.2217 23.6957 23.9603 23.6752 l 24.6513 24.6734 24.6783 24.6760 l 25.0566 25.0566 25.0566 25.0566

11.8800 12.2876 12.0105 11.9393 l 4.4027 0.5178 2.4722 0.3902 l 2.5212 2.5544 2.5528 2.5403 l 1.4081 1.4081 1.4081 1.4081

1.3280 0.2000 0.9668 1.7501 l 3.6518 0.2000 1.9365 0.2000 l 4.7300 4.8046 4.8204 4.8102 l 4.9593 4.9593 4.9593 4.9593

2.0030 0.5000 1.5217 2.3001 l 11.4772 12.0885 11.7810 12.0861 l 12.9999 13.0098 13.0113 13.0083 l 12.3237 12.3237 12.3237 12.3237

7.9340 9.0847 8.3025 8.1838 l 7.6258 7.6423 7.6340 7.6211 l 10.2017 10.2363 10.2420 10.2368 l 10.4637 10.4637 10.4637 10.4637

4.0

10.0

actual value a

1.0

3.0

5.0

10.0

2.0

8.0

Only a local minimum was reached.

We also solved the problem with simulated annealing (SA) in a manner similar to the work by Terry and Messina.11 The SA algorithm was performed with an annealing schedule given by eqs 17:

T(n+1) ) 0.98 × T(n) T

(0)

) 55.0

(17a) (17b)

where T(n+1) is the temperature after N (N ) n × 10) function evaluations. The initial temperature and temperature reduction factor were chosen as 55.0 and 0.98, respectively. No convergence was achieved. Even after 50 000 evaluations of the governing ODEs, the result was still far from the true solution, with the resulting rate vector as kSA) [88.904 0.0441 46.437 0.2975 58.847 0.7083 37.128 0.0837]T and the cost function value of F(k) ) 0.487676. Other trials that included changing the cooling schedule and the probability of acceptance led to no improvements. These experiments suggested the shortcomings of stochastic methods for the KIP. Moreover, the simulations showed that SA was highly dependent on the choice of tuning parameters. 3.2. Complex Chemical Reaction Systems. Two different reaction networks were used to test the PGN method for the KIP of complex reaction systems: (i) the recombination of trifluoromethyl and (ii) H2/O2 combustion.

3.2.1. Recombination of Trifluoromethyl. This reaction system has been experimentally studied by Vakhtin.35 The reaction mechanism and rate constants in appropriate units are given in system (18) (the system defined by eqs 18).

CF3 + NO f CF3NO CF3 + CF3 f C2F6

(R1) (R2)

CF3 + CF3NO f (CF3)2NO NO + (CF3)2NO f (CF3)2(NO)2

(R3) (R4)

(18)

In this case study, the rate expressions for the trifluoromethyl recombination included second-order reversible reaction models, as exemplified in system (19) (the system defined by eqs 19):

d[CF3NO] ) k1[CF3][NO] - k3[CF3][CF3NO] dt d[CF3] ) - k1[CF3][NO] - k2[CF3]2 dt k3[CF3][CF3NO] (19) where k1, k3, and k4 are known values (7.83 × 1012, 8.43 × 1010, and 1.02 × 1013 cm3 mol-1 s-1, respectively). The

Ind. Eng. Chem. Res., Vol. 44, No. 10, 2005 3633 Table 5. Kinetic Mechanism for H2/O2 Combustiona reaction number

Figure 5. Iteration path of k2 and cost function F(k2): (- -4- -) the path starting from k02/k2,0 ) 206.4305, and (s×s) the path with k02/k2,0 ) 0.0036. The correct location of F(k2) is k2/k2,0 ) 1.

KIP consisted of finding k2 from measured [CF3NO] concentrations. In this study, the reported k2 value was used to generate error-free observed data, the same as that observed in the previous study by Tadi and Yetter.4 A total of six data points were taken with the initial mole fractions [CF3] ) 0.5 and [NO] ) 0.5. The boundary conditions were set as ki ∈ [ki,0/1000, ki,0 × 1000]. To study whether the solution by the PGN method is dependent on the initial guess, two arbitrary guess values were tested. One was taken close to the upper bound at k02/k2,0 ) 206.4305 and the other was close to the lower bound at k02/k2,0 ) 0.0036. Figure 5 shows the iteration paths of the convergence of k2 starting from these two arbitrary initial values. In both cases, the optimal solution of k2 was found after 5 iterations. In comparison, an accurate result could be achieved only after 100 iterations by others,4 even when better initial guess were used (k02/k2,0 ) 0.2547). 3.2.2. H2/O2 Combustion. In some KIP, the kinetic parameters of interest may be very sensitive to certain species concentrations, i.e., their values significantly affect the model predictions of species concentration, while being insensitive to others. The following sensitivity discussion reports an observation on a simple but realistic KIP, which exhibits features such as funneling that have been observed by other authors.36 The impacts of varying sensitivities on numerical treatments must be carefully examined. In this work, H2/O2 combustion was studied to test the performance of the PGN method, which involves both situations from measured [OH] data. The kinetic model was taken from the literature37 and is shown in Table 5, including 9 species in 20 reactions. Note that the reverse rate constant of every reaction had to be calculated via an equilibrium constant computed from thermodynamic properties of the involved species at the reaction temperature. The reactants were H2/O2/N2 (in a ratio of 2/1/3.76) and reacted under atmospheric pressure with an initial reaction temperature of 1000 K. Figure 6 demonstrates that reactions R11 and R1 are the most two sensitive reactions to the [OH] concentration. On the other hand, reactions R9 and R17 have little influence on [OH]. 3.2.2.1. [OH] Sensitivity Case. The PGN method was first used to identify the rate constants of reactions R11 and R1 from nine given points of [OH] concentration, which were computed at time intervals of 2.0 × 10-5 s, using the kinetic rates from the literature shown

reaction

A (mol cm s K)

b

E (cal/mol)

R1

H2O Enhanced by 1.86 × 101 H + O2 + M ) HO2 + M 3.61 × 1017

-0.7

0

R2 R3 R4 R5

H2 Enhanced by 2.86 × 100 H + H + M ) H2 + M 1.00 × 1018 9.20 × 1016 H + H + H2 ) H2 + H2 H + H + H2O ) H2 + H2O 6.00 × 1019 H + OH + M ) H2O + M 1.60 × 1022

-1 -0.6 -1.3 -2

0 0 0 0

R6

H2O Enhanced by 5.00 × 100 H + O + M ) OH + M 6.20 × 1016

-0.6

0

R7 R8 R9 R10 R11 R12 R13 R14 R15 R16 R17 R18 R19 R20

H2O Enhanced by 5.00 × 100 O + O + M ) O2 + M 1.89 × 1013 H2O2 + M ) OH + OH + M 1.30 × 1017 H2 + O2 ) 2OH 1.70 × 1013 OH+H2 ) H2O+H 1.17 × 109 O + OH ) O2 + H 3.61 × 1014 O + H2 ) OH + H 5.06 × 104 OH + HO2 ) H2O + O2 7.50 × 1012 H + HO2 ) 2OH 1.40 × 1014 O + HO2 ) O2 + OH 1.40 × 1013 2OH ) O + H2O 6.00 × 108 1.25 × 1013 H + HO2 ) H2 + O2 2.00 × 1012 HO2 + HO2 ) H2O2 + O2 H2O2 + H ) HO2 + H2 1.60 × 1012 1.00 × 1013 H2O2 + OH ) H2O + HO2

0 0 0 1.3 -0.5 2.7 0 0 0 1.3 0 0 0 0

-1788 45500 47780 3626 0 6290 0 1073 1073 0 0 0 3800 1800

a From ref 37. Parameters are given for the equation k ) ATb i exp[-E/(RT)].

Figure 6. Sensitivity of the [OH] concentration, with respect to various reaction constants. (See Table 5 for reaction number identification.) The calculation was performed with the literature values shown in Table 5.

Figure 7. Residual map of the cost function F(k11, k1); the minimum location is (k11/k11,0, k1/k1,0) ) (1, 1).

in Table 5. The boundary conditions were set as ki ∈ [ki,0/100, ki,0 × 100]. Figure 7 displays the residual map of F(k11, k1) that has multiple minima, which is more

3634

Ind. Eng. Chem. Res., Vol. 44, No. 10, 2005

Figure 8. Contour diagram showing multiple local minima in the cost function F(k11, k1). The numbers shown in the figure are the contour levels of F(k11, k1).

apparent in the contour diagram of Figure 8. Nevertheless, our approach successfully located the true 0 and solution from two arbitrary initial guesses of k11 0 0 0 0 k1: (k11/k11,0 ) 12.77, k1/k1,0 ) 0.0446) and (k11/k11,0 ) 0.0723, k01/k1,0 ) 15.54). The computational convergence history is shown in Table 6. In both cases, true values of k11 and k1 were determined within 20 iterations. However, although the PGN method was successful in locating the global minimum in this particular case study, it might fail on certain KIPs with saddle points. Lucia and Feng18,19 have developed global terrain methods that can overcome saddle point limitations. 3.2.2.2. [OH] Insensitivity Case. The PGN method was used to recover the insensitive rate-constant pair k9 and k17 from the same [OH] data set and the identical boundary conditions. The residual map of F(k9, k17) is much smoother than that in the previous problem, as shown in Figure 9. The contour diagram of F(k9,k17) in Figure 10 only has one minimum. Again, two arbitrary 0 /k17,0 ) initial conditions were tested: (k09/k9,0 ) 0.1, k17 0 0 10.0) and (k9/k9,0 ) 10.0, k17/k17,0 ) 0.1). In both cases, the correct values of k9 and k17 were determined within 10 iterations. The computational details are shown in Table 7.

Figure 9. Residual map of the cost function F(k9, k17); the minimum location is (k9/k9,0, k17/k17,0) ) (1, 1).

Figure 10. Contour diagram of the cost function F(k9, k17); the numbers shown in the figure are the contour levels of F(k9, k17).

Thus, the PGN method was successful in identifying both sensitive and insensitive rate constants in this demonstration. The funnel-shaped residual error maps (see Figures 7 and 9) of the KIP have also been observed previously.36 Thus, it seems to be a property inherent to the KIP. 3.3. Ethane Oxidation. The proposed PGN was applied to a large-scale complex reaction network. We intended to determine five optimal parameters in the

Table 6. Identification of Two Most-Sensitive Rate Constants k11 and k1 from Error-Free Experimental OH Data for H2/O2 Combustiona trial 1 iteration

F(k11,k1)

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

2.45 × 9.96 × 1011 9.69 × 1010 4.73 × 109 2.58 × 108 6.56 × 106 1.91 × 106 1.04 × 105 2.66 × 104 1.29 × 103 1.23 × 102 9.11 × 100 6.44 × 10-1 4.59 × 10-2 2.45 × 10-4 6.79 × 10-8 1.68 × 10-16

a

1013

trial 2 k11/k11,0 12.770083 3.538227 2.258532 1.532881 1.224017 1.013324 0.931136 0.844820 0.794321 0.713019 0.690609 0.705125 0.815346 1.029972 0.997452 1.000028 1.000000

k1/k1,0 0.044598 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.050186 0.010000 0.129706 0.296454 0.615485 1.067839 0.994848 1.000083 1.000000

iteration

F(k11,k1)

k11/k11,0

k1/k1,0

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

8.60 × 7.67 × 1013 8.58 × 100 1.27 × 1011 8.48 × 109 6.73 × 108 2.97 × 108 1.55 × 107 4.97 × 106 6.77 × 105 7.67 × 104 2.16 × 104 1.33 × 103 1.22 × 102 8.96 × 100 5.98 × 10-1 3.80 × 10-2 1.70 × 10-5 2.27 × 10-10

0.072299 52.523546 0.010000 2.408781 1.682659 1.352050 1.246343 1.058172 0.993823 0.898089 0.854183 0.784875 0.718393 0.702576 0.719612 0.825097 1.018643 0.999058 1.000000

15.540166 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.111723 0.037687 0.026906 0.161183 0.328366 0.633518 1.043906 0.997922 1.000000

The correct values were recovered starting from two arbitrary initial guesses.

100

Ind. Eng. Chem. Res., Vol. 44, No. 10, 2005 3635 Table 7. Identification of Two Less-Sensitive Rate Constants k9 and k17 from Error-Free Experimental OH Data for H2/O2 Combustion trial 1

trial 2

iteration

F(k9,k17)

k9/k9,0

k17/k17,0

iteration

F(k9,k17)

k9/k9,0

k17/k17,0

0 1 2 3 4 5 6

5.79 × 100 9.92 × 100 1.59 × 10-1 9.82 × 10-4 2.61 × 10-7 1.02 × 10-8 5.45 × 10-10

0.100000 1.133059 1.048882 0.999294 1.000118 1.000000 1.000000

10.000000 0.010000 0.970160 0.980640 1.000160 0.999920 1.000000

0 1 2 3 4 5 6 7 8 9 10

2.88 × 105 7.95 × 104 4.44 × 104 1.43 × 103 1.76 × 101 3.21 × 100 1.33 × 10-3 9.17 × 10-6 1.83 × 10-7 1.06 × 10-8 6.08 × 10-10

10.000000 7.711176 5.549765 4.109176 1.303471 1.040412 1.011353 0.999118 1.000118 1.000000 1.000000

0.100000 0.010000 0.010000 1.444960 0.010000 0.419136 1.036480 0.997200 1.000320 0.999920 1.000000

a

The correct values were recovered starting from two arbitrary initial guesses.

GRI-Mech 3.0 reaction mechanism,38 to interpret experimental measurements at elevated pressures in the presence of experimental errors. GRI-Mech 3.0 is a detailed chemical reaction mechanism that contains 53 species and 325 reactions capable of describing the chemistry of natural gas flames and ignitions. It was optimized for the combustion of light hydrocarbon fuels in the temperature range of 1000-2500 K with pressures in the range of 0.014-10 atm and an equivalence ratio in the range of 0.1-5 for premixed systems. Recent experimental work conducted with a high-pressure shock tube showed that GRI-Mech 3.0 failed to give good predictions for ethane oxidation at elevated pressures.39 Therefore, it is of great interest to determine whether better rate parameters can be extracted from our experimental measurements. The experimental data were obtained from stoichiometric reaction mixtures (φ ) 1) with slightly different reaction times of ∼0.0018 s at a nominal reaction pressure of 340 bar. For each shock, four stable species (C2H6, C2H4, C2H2, and CO) were quantitatively measured using GC/MS analytic techniques. A total of 21 experiments at different temperatures were performed, and 84 obtained data points were used for the KIP. The five reactions considered for optimization are listed in system (20) (the system defined by eqs 20), of which the series number in the original GRI-Mech 3.0 reaction model is indicated on the right-hand side and the letter “(M)” in the reaction expressions represents third-body gases.

H + CH3 (+ M) T CH4 (+ M) (R52 in GRI-Mech 3.0) H + C2H4 (+ M) T C2H5 (+ M) (R74 in GRI-Mech 3.0) 2CH3 (+ M) T C2H6 (+ M) (R158 in GRI-Mech 3.0) 2CH3 T H + C2H5

(R159 in GRI-Mech 3.0)

C2H5 + O2 T HO2 + C2H4 (R175 in GRI-Mech 3.0) (20) All of the rate constants have a modified Arrhenius expression of ki ) AiTbi exp[-Ei/(RT)] with Ai as the preexponential factor, bi as the temperature-dependent parameter, and Ei as the reaction activation energy. Ongoing work indicates that the pressure-dependent features of these reactions have not yet been fully

understood.40,41 Therefore, our purpose was to determine the optimal rate parameters to fit high-pressure shock tube data. The pre-exponential factors of the aforementioned five reactions were chosen for optimization, and the boundary condition was set to Ai ∈ [Ai,0/10, Ai,0 × 10], where Ai,0 is the literature value in GRI-Mech 3.0. Using the PGN method, the optimal pre-exponential factors were determined. Figure 11 compares the agreement between the model predictions of the improved and the original GRI-Mech 3.0 and the actual experimental data. The figure clearly indicates that the deviations between model predictions and experimental measurements are largely reduced for both C2H6 and C2H4. Considerable improvement also can be observed, with respect to C2H2 and CO over the entire temperature window. With the optimal rate parameters obtained by the proposed PGN method, the overall performance of GRI-Mech 3.0 was greatly improved. The computation converged in just 10 iterations, and the CPU time was ∼2.0 h. The time-controlling step was the integration of eq 1 and evaluation of F(k) in system (2), with ∼10 min of CPU time consumed in each iteration. This should not be a surprise, because the governing ODEs are constructed by a large reaction system, containing 325 reactions and 53 species. If other gradient-based methods2,3 were applied for this problem, months would be needed to solve the problem, because they generally demand an order of 104 function evaluations. The computational effort is even more prohibitive for a stochastic search (e.g., the SA, GA, and MC algorithms), because this method generally requires 105-106 iterations for F(k) evaluation. This case study provides a benchmark work for the computational efficiency of the PGN for large-scale kinetic reaction problems. 4. Multiple Solutions and Global Optimality It is of great significance to address the issue of multiple solutions and global optimality when applying numerical methods for KIP. It is well-known that, in KIP, multiple solutions may exist with equal error margins.42,43 Neither the parameter set nor the extracted reaction scheme is necessarily unique to a given experimental data sample. Different reaction schemes can generate the same species trajectories, even when different parameter sets are used. For example, Vajda and Rabitz42 demonstrated the unidentifiability and indistingushability in first-order reaction systems; Zsely and co-workers also demonstrated that very different combustion kinetic models could produce similar re-

3636

Ind. Eng. Chem. Res., Vol. 44, No. 10, 2005

Figure 11. Comparison of shock tube data and model predictions for ethane oxidation at 340 bar: (a) C2H6, (b) C2H4, (c) C2H2, and (d) CO. Legend is as follows: (9) experimental data, (- - -) GRI-Mech 3.0, and (s) GRI-Mech 3.0 with optimal parameters.

sults.43 The inherent errors in the experimental data caused by instruments and experimental procedures contribute to the multiplicity of solutions. Among those multiple solutions to the KIP, the “global” solution might be mathematically attractive by offering the smallest residual error. However, this global residual error minimum is not guaranteed to be physically better than other local minima. The proposed PGN approach, on the other hand, offers the advantage that it maintains rate constants within the specific physically meaningful bounds. However, it is still a local method, which means that the PGN method may be trapped in a local minimum. The main advantage of the PGN method is its fast convergence. This is essential for the solvability of the KIP of complex reaction systems. 5. Conclusions In this paper, the physically bounded Gauss-Newton (PGN) method has been formulated and tested on various nonlinear dynamically constrained kinetic inversion problems (KIPs). It distinguishes itself from other gradient-based methods by determining the update step length through solution of a quadratic problem with a physical boundary. This approach provides large step lengths while ensuring the parameters’ physical boundedness. The computational effort is reduced by solving repeatedly a constrained KIP directly, using accurate first-order response surfaces. The true rate constants were successfully recovered for the case studies presented in this article. The proposed algorithm was also compared against a global optimization algorithm. Using the PGN approach, the optimal preexponential factors of five target reactions were determined in a small number of iterations. The overall performance of GRI-Mech 3.0 was greatly improved in

fitting ethane oxidation data obtained under extreme reaction conditions. Acknowledgment K.B. gratefully acknowledges the support of the National Science Foundation, Combustion and Plasma System Program, under Contract No. CTS 0109053. Literature Cited (1) Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics, 2nd Edition; Prentice Hall: Upper Saddle River, NJ, 1999. (2) Lisy, J. M.; Simon, P. Evaluation of Parameters in Nonlinear Models by the Least Square Method. Comput. Chem. 1998, 22, 509. (3) Farinha, J. P. S.; Martinho, J. M. G.; Pogliani, L. Nonlinear Least Squares and Chemical Kinetics. J. Math. Chem. 1997, 21, 131. (4) Tadi, M.; Yetter, R. A. Evaluation of the Rate Constants in Chemical Reactions. Int. J. Chem. Kinet. 1998, 30, 151. (5) Tadi, M.; Cai, W. Inverse Matrix Evaluation for Linear Systems. Inverse Probl. Eng. 1998, 6, 15. (6) Cervantes, A. M.; Waechter, A.; Tutuncu, R.; Biegler, L. T. A Reduced Space Interior Point Strategy for Optimization of Differential Algebraic Systems. Comput. Chem. Eng. 2000, 24, 39. (7) Arora, N.; Biegler, L. T. A Trust Region SQP Algorithm for Equality Constrained Parameter Estimation with Simple Parameter Bounds. Comput. Optim. Appl. 2004, 28, 51. (8) Arora, N.; Biegler, L. T. Parameter Estimation for a Polymerization Reactor Model with a Composite-Step Trust-Region NLP Algorithm. Ind. Eng. Chem. Res. 2004, 43, 3616. (9) Locatelli, M. Simulated Annealing Algorithm for Continuous Global Optimization: Convergence Conditions. J. Optim. Theory Appl. 2000, 104, 121. (10) Eftaxias, A.; Font, J.; Fortuny, A.; Fabregat, A.; Stuber, F. Nonlinear Kinetic Parameter Estimation Using Simulated Annealing. Comput. Chem. Eng. 2002, 26, 1725.

Ind. Eng. Chem. Res., Vol. 44, No. 10, 2005 3637 (11) Terry, D. B.; Messina, M. Heuristic Search Algorithm for Determination of Rate Constants and Reaction Mechanisms from Limited Concentration Data. J. Chem. Inf. Comput. Sci. 1998, 38, 1232. (12) Alper, J. S.; Gelb, R. I. Monte Carlo Method for the Determination of Confidence Intervals: Analysis of Nonnormally Distributed Errors in Sequential Experiments. J. Phys. Chem. 1991, 95, 104. (13) Frenklach, M.; Wang, H.; Rabinowitz, M. J. Optimization and Analysis of Large Chemical Kinetic Mechanisms Using the Solution Mapping Method: Combustion of Methane. Prog. Energy Combust. Sci. 1992, 18, 47. (14) Eiteneer, B.; Frenklach, M. Experimental and Modeling Study of Shock Tube Oxidation of Acetylene. Int. J. Chem. Kinet. 2003, 35, 391. (15) Shenvi, N.; Geremia, J. M.; Rabitz, H. Nonlinear Kinetic Parameter Identification through Map Inversion. J. Phys. Chem. A 2002, 106, 12315. (16) Esposito, W. R.; Floudas, C. A. Global Optimization for the Parameter Estimation of Differential-Algebraic Systems. Ind. Eng. Chem. Res. 2000, 39, 1291. (17) Papamichail, I.; Adjiman, C. S. Global Optimization of Dynamic Systems. Comput. Chem. Eng. 2004, 28, 403. (18) Lucia, A.; Feng, Y. Global Terrain Methods. Comput. Chem. Eng. 2002, 26, 529. (19) Lucia, A.; Feng, Y. Multivariable Terrain Methods. AIChE J. 2003, 49, 2553. (20) Kearfott, R. B. Rigorous Global Search: Continuous; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996. (21) Esposito, W. R.; Floudas, C. A. Parameter Estimation in Nonlinear Algebraic Models via Global Optimization. Comput. Chem. Eng. 2002, 22, 213. (22) Caracotsios, M.; Stewart, W. E. Sensitivity Analysis of Initial Value Problems with Mixed ODEs and Algebraic Equations. Comput. Chem. Eng. 1985, 9, 359. (23) Li, S.; Petzold, L. R. Software and Algorithms for Sensitivity Analysis of Large-Scale Differential Algebraic Systems. J. Comput. Appl. Math. 2000, 125, 131. (24) Leis, J. R.; Kramer, M. A. Sensitivity Analysis of Systems of Differential and Algebraic Equations. Comput. Chem. Eng. 1985, 9, 93. (25) Saltelli, A.; Chan, K.; Scott, M. Sensitivity Analysis; Wiley: New York, 2000. (26) Dennis, J. E.; Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; SIAM: Philadelphia, PA, 1996. (27) Petzold, L. R. A Description of DASSL, Sandia National Laboratories Report SAND82-8637, Sandia National Laboratories, Albuquerque, NM, 1982. (28) Kee, R. J.; Rupley, F. M.; Miller, J. A. Chemkin-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase

Chemical Kinetics, Sandia National Laboratories Report SAND898009, Sandia National Laboratories, Albuquerque, NM, 1990. (29) Lin, C. J.; More, J. Newton’s Method for Large BoundConstrained Optimization Problems. SIAM J. Optim. 1999, 9, 1100. (30) Zhu, C.; Byrd, R. H.; Nocedal, J. L-BFGS-B: Fortran Routines for Large Scale Bound Constrained Optimization. ACM Trans. Math. Soft. 1997, 23, 550. (31) Byrd, R. H.; Lu, P.; Nocedal, J. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM J. Sci. Comput. 1995, 16, 1190. (32) Tranter, R. S.; Tang, W.; Anderson, K. B.; Brezinsky, K. Shock Tube Study Thermal Rearrangement of 1,5-Hexadiyne over Wide Pressure and Temperature Range. J. Phys. Chem. A 2004, 108, 3406. (33) Villadsen, J. V.; Michelsen, M. L. Solution of Differential Equations by Polynomial Approximation; Prentice Hall: Upper Saddle River, NJ, 1978. (34) Sahinidis,.V. BARON: A General Purpose Global Optimization Software Package. J. Glob. Optim. 1996, 8, 201. (35) Vakhtin, A. B. The Rate Constant for the Recombination of Trifluoromethyl at T ) 296 K. Int. J. Chem. Kinet. 1996, 28, 443. (36) Lucia, A.; DiMaggio, P. A.; Depa, P. Funneling Algorithms for Multiscale Optimization on Rugged Terrains. Ind. Eng. Chem. Res. 2004, 43, 3770. (37) Kim, T. J.; Yetter, R. A.; Dryer, F. L. New Results on Moist CO Oxidation. Proc. Combust. Inst. 1994, 25, 759. (38) Smith, G. P.; Golden, D. M.; Frenklach, M.; Moriarty, N. W.; Eiteneer, B.; Goldenberg, M.; Bowman, C. T.; Hanson, R. K.; Song, S.; Gardiner, W. C. Jr.; Lissianski, V. V.; Qin, Z. http:// www.me.berkeley.edu/gri_mech/, accessed in 2004. (39) Tranter, R. S.; Sivaramakrishnan, R.; Brezinsky, K.; Allendorf, M. D. High Pressure, High-Temperature Shock Tube Studies of Ethane Pyrolysis and Oxidation. Phys. Chem. Chem. Phys. 2002, 4, 2001. (40) Wagner, A. F. The Challenges of Combustion for Chemical Theory. Proc. Combust. Inst. 2002, 29, 1173. (41) Piling, M. J.; Robertson, S. H. Master Equation Models for Chemical Reaction of Importance in Combustion. Annu. Rev. Phys. Chem. 2003, 54, 247. (42) Vajda, S.; Rabitz, H. Identifiably and Distinguishability of First-Order Reaction System. J. Phys. Chem. 1988, 92, 701. (43) Zsely, I. G.; Zador, J.; Turanyi, T. How Can Very Different Combustion Models Produce Similar Results? Proc. Combust. Inst. 2004, 30, poster 1F2-11.

Received for review November 22, 2004 Revised manuscript received February 23, 2005 Accepted February 28, 2005 IE048872N