Distributed System Design Under Uncertainty - Industrial

However, many technically important reaction processes involve chemical reactions with strong heat and mass transfer coupling. This article will study...
0 downloads 5 Views 344KB Size
Download PDF
8352

Ind. Eng. Chem. Res. 2006, 45, 8352-8360

Distributed System Design Under Uncertainty Libin Zhang, Celia Xue, Andres Malcolm, Kedar Kulkarni, and Andreas A. Linninger* Laboratory for Product and Process Design, Department of Chemical Engineering and BioEngineering, UniVersity of Illinois at Chicago, Chicago, Illinois 60607

Many chemical reactions are intrinsically coupled with transport phenomena. In realistic systems, both the chemical reactions as well as the transport phenomena are functions of space. Despite spatial distribution, it is customary to fit transport and kinetic properties in experiments to simplified lumped models. This paper outlines a method to acquire unknown system properties in reacting systems with transport coupling. Rigorous metrics to quantify confidence levels of these distributed models and their parameters to represent the experimental data will also be computed. We will advocate the advantage of deploying confidence information together with apparent physical and chemical properties for flexible process design. Hence, we will find the most cost-effective design subject to the accurate statistical error information of the model parameters. The paper will introduce mathematical foundations and nonlinear programming solutions based on trust region methods with inexact function and gradient evaluations specifically designed for inversion problem of distributed systems. A comprehensive example for the optimal design of catalytic pellet reactor will illustrate the proposed approach. 1. Introduction Chemical process development traverses through two stages. In the first stage, experiments provide data about chemical and physical state changes so that proper mathematical models and their parameters can be established. This stage is referred to as model identification. The second stage uses the mathematical model to design the process for optimal performance. Process optimization maximizes the expected performance to determine the best design variables such as reaction temperature, operating pressure, and reactor dimensions. Current chemical reactors models often focus on lumped models without considering the spatial variation of the velocity, temperature, and concentration fields. Model parameters obtained from such experiments do not account for mass and heat transport. However, many technically important reaction processes involve chemical reactions with strong heat and mass transfer coupling. This article will study process design based on experimental data using two-dimensional distributed reaction and transport models. The advantages of estimating reaction and transport properties in distributed systems simultaneously will be demonstrated. Problem formulation and numerical solutions techniques will be discussed. This activity is termed problem inVersion. Many researchers have worked in inversion problems, but most approaches only focus on the kinetic mechanisms without considering their spatial dependence.1-4 Most parameter estimation work does not address simultaneously reaction and transport phenomena. Prior work in distributed systems concerned inverse heat conduction, air pollution, and data assimilation for meteorological and/or oceanic models.5 Inverse heat transfer problems have been studied because of the importance of aerodynamic heating of space vehicles, whose thermal shell temperature during reentry into the atmosphere is higher than can be measure with thermocouples directly.6 Borggaard and Burns7 presented a trust region approach to determine the shape parameters for matching * To whom correspondence concerning this article should be addressed. Phone: (312) 413-7743. Fax: (312) 996-5921. E-mail: [email protected].

the flow to aerodynamic design. Truncated Newton’s optimization with finite elements was implemented to reconstruct the tissue optical property maps from exterior measurements of frequency-domain photon migration (FDPM).8 Another shortcoming in existing distributed process development is the lack in mathematical methods to quantify uncertainty. Most designs optimized for nominal conditions are subsequently relaxed to accommodate uncertainty. If uncertainty sources are considered, the prevailing assumption is to treat them as statistically independent of each other. This parametric independence, however, cannot be justified for coupled transport and kinetic properties. The coupling of the process parameters can be quantified with the help of the experimental data.9,10 In this article, we propose to extract the optimal parameters along with their correlations for distributed systems from experimental data. Once the parameters and the information regarding their uncertainty are obtained rigorously, we propose to conduct flexible design using probabilistic optimization to maximize expected performance. The article is organized as follows. Section 2 will introduce the mathematical background to formulate inversion problems for distributed systems. Transport problem discretization by the finite volume approach will be shown. We will propose inexact trust region methods to solve the parameter estimation problem as well as to obtain the desired uncertainty information. In Section 3, we will illustrate the application of a transport and kinetic inversion problem using a tubular flow reactor for catalyzed desulfurization. After identifying unknown reaction rates and transport properties, rigorous design optimization under uncertainty will be the subject of Section 4. The paper closes with conclusions. 2. Design for the Distributed Systems Consider the proposed process development methodology in two stages shown in Figure 1. The first stage seeks to identify unknown reaction and transport parameters from experimental observations. Ideally, we would like to discover reaction and transport phenomena with limited concentration or temperature

10.1021/ie060082l CCC: $33.50 © 2006 American Chemical Society Published on Web 09/09/2006

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8353

Figure 2. Typical balance envelope is considered by the finite volume method in two dimensions.

permeability as well as reaction kinetic rates appearing in source terms S(x,θ,φ). The unknown parameters, θ, in system 1 will be recovered by least-squares optimization against the experimental data set. F is the covariance matrix of experimental data accounting for measurement accuracy.

minψ(P) ) (φ(x,θ) - φˆ (x))TF-1(φ(x,θ) - φˆ (x))

(1)

F∇(u bφ) ) ∇(Γ∇φ) + S(x,θ,φ)

(2)

θ

s.t.

Figure 1. Overall information flowsheet of the proposed flexible design methodology.

measurements collected at convenient positions of the reaction experiment. The optimal parameter set is obtained by solving a least-squares-fitting problem between the observations and the proposed distributed model. Statistical information characterizing the confidence in the models’ parameters with respect to experimental data is computed simultaneously. In the second stage, we propose a flexible design method based on accurate probability distributions for each parameter. The design optimization maximizes performance within the uncertainty space delineated in the first stage. Both subproblems entail mathematical programs with partial differential equation constraints. Three main ideas help solve the distributed inversion problem: (i) discretization of transport equations using a finite volume approach; (ii) inexact trust region methods to solve large distributed optimization problems; and (iii) beneficial use of response maps to assemble efficiently and robustly parameter sensitivities and the Hessian matrix for large-scale transport and kinetic inversion problems. 2.1. Parameter Estimation in Distributed System. Systems 1 and 2 define the steady-state parameter estimation problem in the distributed system. Experiments yield known measurements of state variables φˆ (x) at locations x. They can include species concentrations, temperatures, or flow velocity measurements. First principles mathematical models quantify the distributed states, φ(x,θ), in terms of unknown reaction and transport properties, θ. The system of constraints (eq 2) represents, generically, conservation laws for mass (φ ) 1), momentum (φ ) u), species (φ ) C), or energy transfer (φ ) H). The set of unknown system parameters, θ, may include unknown transport properties such as diffusivity, porosity, and

The transport and kinetic inversion problem presented above constitutes a nonlinear optimization problem with partial differential equation (PDE) constraints.5 Before solving the problem numerically, the PDE constraints need to be converted into algebraic equations, as discussed in the next section. 2.2. Discretization of the PDE Constraints: The Finite Volume Method. For a large-scale transport and kinetic inversion problem (TKIP), we propose the finite volume method for the spatial discretization of the transport equations. This process will produce algebraic equations instead of eq 2. In the finite volume method, the spatial domain is divided into a discrete number of nonoverlapping control volumes. The conservation balances of eq 2 can be formulated over each volume shown in Figure 2 according to eq 3.

∂φ ∂ ∂φ ∂ ∂φ + FV - Γ - Γ - S) dV ) 0 ∫V (Fu ∂φ ∂x ∂y ∂x( ∂x ) ∂y( ∂y ) i

∀ Vi ∈ V (3)

Integration of eq 3 along the main coordinate directions with proper profile assumptions yields algebraic expressions for each transport of quantity, φ, across each control volume. Similar balances for all other control volumes convert the continuous PDE into a set of algebraic constraints for all states, φ, as given in eq 4.

f(φ) ) aPφP - aEφE - aWφW - aNφN - aSφS S‚∆x‚∆y ) 0 ∀ Vi θ ) {Γ,k} (4) An attractive feature of the finite volume formulation is the integral conservation of quantities such as mass, momentum, and energy over each control volume individually as well as over the entire problem domain. This property is harder to enforce in finite difference and finite element methods.11 After the constraints’ discretization, the problem becomes a large-scale nonlinear optimization problem (NLP). To solve this NLP, we propose a reduced space optimization approach based on the inexact trust region method. 2.3. Mathematical Programming Solution. Two challenges pertaining to the transport and kinetic inversion problem require special attention: (i) its size due to the discretization of the problem domain and (ii) difficulty in convergence, especially

8354

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006

Figure 3. Approximate quadratic error surface in each iteration for trust region method.

when large variations in the parameter set cause physically meaningless departures. To attack these problems, we propose an inexact trust region method with the following advantages: (i) a drastic design variable set reduction by separating the search for unknown parameter updates from the solution of the distributed transport equations; (ii) a rigorous enforcement of physically meaningful bounds on parameter updates in each iteration; and (iii) efficient computation of parameter sensitivity information to provide necessary gradient information. As a result, our approach converges globally by sequentially solving linear optimization subproblems in which the true error surface is replaced by first-order response maps leading to quadratic subproblems, as shown in Figure 3. Section 2.4 will discuss the efficient computation of sensitivity information. Section 2.5 describes the implementation of the proposed inexact trust region method. 2.4. Response Maps of Transport Equations. Solving objective 1 subject to discretized transport eq 4 entails finding the derivatives of each transport quantity, φi, with respect to each unknown parameter, θ. The computation of the sensitivity information requires analytical differentiation of the transport equations with the help of the chain rule. Implicit differentiation produces additional M × N sensitivity equations as in eq 5, where M is the total number of unknown state variables of the transport equation and N is the number of optimization parameters.

∑ | i)1 ∂φ M

∂h

i θ,φ

∂φi ∂θj

+

∂h

|

∂θj

)0

j ) 1, ..., N

(5)

θ,φ

The unknown sensitivity variables, ∂φ/∂θj ) ∇θjφ, can be expressed in matrix form as a linear algebraic problem given in eq 6.

∇φh‚∇θjφ ) -∇θjh

j ) 1, ..., N

(6)

The known M × M matrix,∇φh, contains the derivatives of the transport equations with respect to the state variables φ. The right-hand side vector,∇θjh, represents the partial derivatives of transport equations with respect to the unknown parameters, θ. The same derivative matrix,∇φh, also appears in the transport problem when solved with Newton’s methods. Thus, the state Jacobian matrix,∇φh, has to be factorized only once; the computation of the sensitivity equations can then be done inexpensively.12-15

2.5. Implementation of Inversion Problems with Inexact Newtons Methods. The trust region method is a gradient-based nonlinear optimization routine with step-size control.16-18 It chooses the search direction for parameter updates as a combination of the Newton and the steepest descent direction as depicted in Figure 4a. Accordingly, it can also be deployed when the Jacobian or Hessian matrix is singular or near singular. This advantage is crucial in ill-conditioned transport and kinetic inversion problems. Even when using accurate sensitivity information, the numerical solution of TKIP is time-consuming and often fails to converge. As a remedy, the proposed trust region method limits the maximum permissible parameter update from one iteration to the next. The name trust region derives from the concept of the scalar radius, γ, sizing the maximum permissible parameter change (|∆θ| e γ). Moreover, it is not necessary to solve for state variables and derivative information exactly when the algorithm is still far from the final solution. Therefore, the algebraic state eq 2 is not solved exactly, and only an approximation φ˜ is computed. Similarly, gradient computations requiring the solution of linear subproblems are converged only partially to yield approximate sensitivity information, J˜. The inexact trust region method can then be formulated as in eqs 7 and 8. Its iterative search directions in comparison to the classical trust region methods are shown in Figure 4b.

min ψ ˜ (θk + ∆θ) ) ∆θ

(φ˜ (θk) + J˜‚∆θ - φˆ (x))T(φ˜ (θk) + J˜‚∆θ - φˆ (x)) (7) s.t. |∆θ| e γ

(8)

The proposed inexact trust region method uses approximation for the state, φˆ , and the sensitivity, J˜ ) ∂φ˜ /∂θ, but is otherwise identical to the traditional trust region method. What are the accuracy requirements on evaluations of functions and derivative information? To keep the strong global convergence for the inexact trust region method, the relative error of the Jacobian, J˜, has to be bounded in each iteration. Some researchers have shown global convergence of the inexact trust region algorithm provided that the relative gradient error satisfies eq 9.16,17

|ek| eσ |gk|

(9)

The expression |ek| is the norm of residual errors of the discretized transport functions in each iteration. The scalar, |gk|, represents the first-order derivative of the least-squares objective function, which is (φ(x,θk) - φˆ (x))TJ˜. The user-specified constant σ has a typical value, σ ) 0.8.16 The inexact algorithm summarized in Table 1 uses the same steps as the classical trust region method. The detailed discussion and implementation of the algorithm can be found elsewhere.18-20 A successful application of the proposed methodology is demonstrated by the design of a catalytic pellet reactor under uncertainty. 3. Application In this section, we present the design of a catalytic pellet reactor under uncertainty. First, we recover the optimal parameters and their variability from experimental data using transport and kinetic inversion. We propose to sample the uncertain space using the accurate probability information. Finally, using

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8355

Figure 4. Trust region step, pTR, for transport optimization problem. (a) Trust region: The trust region step, ∆θTR, is a linear combination of the Newton, ∆θN, and steepest descent, ∆θsd, direction inside the trust region sphere of radius γ. (b) Inexact trust region: The inexact trust region step, ∆θˆ TR, is a linear combination of the inexact Newton, ∆θˆ N, and an inexact steepest descent, ∆θˆ sd, direction.

Table 1. Algorithm of Inexact Trust Region Method Give initial guesses for the unknown parameters, θ0, and trust region radius γ. Do: Solve the distributed model eq 4 and its sensitivity equation eq 6 inexactly. Calculate ∆θˆ sd, β, and ∆θˆ N. Repeat: if: γ e β|∆θˆ sd|, then: ∆θˆ ) γ∆θˆ sd/|∆θˆ sd|. “Cauchy point outside trust region” Else if: γ g |∆θˆ N||, then: ∆θˆ ) ∆θˆ N. “Newton inside trust region radius” Else: |∆θˆ N| < γ < β|∆θˆ sd|, then: ∆θˆ ) η∆θˆ N + (1 - η)β∆θˆ sd, where η ) (γ - β|∆θˆ sd|)/(|∆θˆ N| - β|∆θˆ sd|). “Trust region radius between Cauchy point and Newton” End Update θk+1 ) θk + ∆θˆ . Determine whether θk+1 is acceptable by checking whether |ψ(θk+1)| < |ψ(θk)|, and calculate a new value of γ. Until θk+1 is an acceptable next point.

stochastic optimization techniques, the best design will be obtained with minimum “expected” cost. The experimental setting for this example generates concentration measurements for all species at the inlet and the outlet of the reactor with an accuracy of 10%. In addition, temperature measurements are conveniently available along the reactor axis with a variance of 5%, as depicted in Figure 5. The goal of the experiment is to determine the diffusivity of the reactor, DA; the reaction constant, Kref; the activation energy, E; and the mass transfer coefficient, hs, based only on these limited experimental data. 3.1. Transport and Kinetic Inversion Problem. The reaction rate in the catalytic fixed-bed reactor is distributed with spatial and temporal variation. In this section, four parameters of a distributed reaction model in the fixed-bed reactor, such as reactor diffusivity, reaction rate, activation energy, and mass transfer coefficient between bulk flow and pellet, are recovered by problem inversion. 3.2. Catalytic Pellet Reactor. Catalytic pellet reactors have many commercial applications for reducing automotive emission

8356

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 Table 2. Physical Properties of the Catalytic Pellet Reactor Physical Properties of the Pellet  ) 0.4 (porosity) De ) 5.2 × 10-6 m2/s R0 ) 0.003 m Fc ) 2.8 × 103 kg/m3 Physical Properties of The Bed b ) 0.4 Fb ) Fc(1 - b) Physical Properties of The Inlet Gas To ) 773 K Reaction Conditions ∆H ) 97 200 J/mol

The reactions in the bulk are coupled with the kinetics and mass transfer of the porous pellets via the effectiveness factor, η. The effectiveness factor accounts for mass transfer and reaction inside the pellet and expresses the overall reaction rate in terms of the bulk concentration and temperature. The rate of reaction is expressed by an Arrhenius-type equation, as shown in eq 14.

[ (

rA ) -ηKref exp Figure 5. Multiple length scale of a catalytic pellet reactor (a) and a schematic of measurements (b) to determine reaction, mass transfer, and heat transfer properties.

gas, complete oxidation of hydrocarbons, paraffin dehydrogenation, hydrocracking, and dehydrocyclization. A typical tubular reactor uses a packed bed of spherical catalyst pellets, shown in Figure 5. The pellets are porous media supporting diffusion and heterogeneous chemical reactions on their outer and inner surfaces.21 If the chemical reaction is accompanied with a heat effect, significant temperature gradients can develop within the pellet. Diffusive and convective transport of species also occur in the bulk. The whole reactor model traverses multiple length scales from the reactor scale to the pellets down to the micropores. In this case study, a sulfur dioxide (SO2) oxidation catalyst for production of sulfuric acid is chosen. The overall reaction given in eq 10 is first order with respect to SO2. Catalyst 1 SO2(g) + O2(g) 798 SO3(g) 2

(10)

3.3. Two-Dimensional Reactor Model with Porous Catalyst Packing. The species and heat balances in axial and radial directions for a cylindrical control volume of the catalytic reactor filled with pellet packing, as shown in Figure 5, are given in eqs 11 and 12, respectively.

∇(u bCA) ) ∇(DA‚∇CA) + rA

(11)

FgCp∇(u bT) ) ∇(ke‚∇T) + rA∆H

(12)

where rA is the overall reaction rate within the catalyst per unit volume. Since the reactor is packed with small catalytic pellets, its content can be treated as a porous material. As a first approximation, the momentum transfer in the pellet reactor is described by Darcy’s law as in eq 13.

-∇P B)

µ b u km

(13)

)]

E Tref - 1 CA RTref T

(14)

The external effectiveness factor η relates the actual conversion in the pellet to a homogeneous reference.22 All constants are listed in Table 2. The effectiveness factor links the bulk conversion of reactant to the transport and catalytic reactions inside the pellet. 3.4. Catalytic Pellet Model. The catalytic pellets can be treated as a porous medium with coupled reactions, mass and heat transfer. The reactions in the pellet are assumed to be fast enough to justify the steady-state assumption for the pellets. The fully coupled, nonisothermal concentration profile inside the pellet can be expressed in dimensionless form as in eq 15, with boundary conditions in eq 16:

˜A d2C dR2 De

+

(

)

1- C ˜A 2 -Φ2 dC ˜ A exp γβ ) 0 (15) R 1 + β(1 - C ˜ A)

|

dC ˜A dR

R)1

) R0hs(1 - C ˜ A) and

|

dC ˜A dR

R)0

)0

(16)

The dimensionless variables represent the Thiele modulus, Φ2 ) R02Kref/De; the species concentration, C ˜ A ) CA/C0; the pellet radius, R ) r/R0; the activation energy, γ ) E/RgT0; and the heat generation, β ) C0∆HDe/keT0.23,24 Here, C0 and T0 are taken to be the concentration of the reactant and the temperature on bulk flow, respectively. The boundary condition in eq 16 shows that reactant concentrations at the pellet surface are not identical to bulk concentration. The mass flux can be determined with the help of the mass transfer coefficient, hs. At the center of the pellet, fluxes are zero, dC ˜ A/dR ) 0. Orthogonal collocation on finite elements is employed to solve for C ˜ A in each volume.18 Figure 6 displays typical dimensionless concentration and temperature profiles in catalytic pellets. 3.5. Results. The concentrations of reactants are measured at the inlet and outlet of the reactor. Temperatures are obtained at different locations along the reactor. With the measured inlet and outlet concentration and the temperatures inside the reactor, the diffusivity of the reactor, DA; the reaction constant, Kref; the activation energy, E; and the mass transfer coefficient, hs, were determined according to the inversion problem (eq 17).

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8357

Figure 6. (a) Concentration and (b) temperature profiles in a catalytic pellet.

Measurement errors in the covariance matrix F were 5% for temperature and 10% for concentration measurements.

min ψ(θ) ) θ,C,T

([ ] [ ]) ([ ] [ ]) C C ˆ T Tˆ

T

F-1

C C ˆ T Tˆ θ ) {DA,Kref,E,hs} (17)

s.t. eqs 11-16 The optimal parameters, such as the bulk diffusivity of species A in the carrier gas, DA; the reaction constant, Kref; the activation energy, E, for the heterogeneous reaction; and the mass transfer coefficient from the bulk to the catalytic pellet, hs, can be computed from experimental data with help from the solution of eqs 1-5 with the inexact trust region method described in Section 2. After 30 iterations, the following results were obtained: DA ) 7.32 × 10-5 (m2/s); Kref ) 0.82 × 10-4 (1/s); E ) 5268.62 (J/mol); and hs ) 1.52 (m/s). Figure 7 depicts the concentration and temperature fields predicted for the optimal parameters. 3.6. Assessment of Parameter Variability. Rigorous computation of the parameter confidence aims at determining the range of the uncertain parameter variations as a function of the quality of the underlying experimental data. The Jacobian, J, accounting for the sensitivity information permits the construction of probabilistic parameter ranges. Probability density functions can be formulated with the help of confidence regions. Confidence regions, in turn, are related to the parameter covariance matrix, Fθˆ , which can be expressed with the help of the Jacobian as in eq 18, where n is the number of experimental data, m is the number of estimation parameters, and e is the error vector between the experimental and the model values.25

Fθ ) s2(J(θ)TJ(θ))-1

Figure 7. Concentration and temperature distribution in the catalytic pellet reactor: (a) steady-state axial concentration profile and (b) steady-state axial temperature distribution.

8a. Similarly, joint confidence region for the parameters can be constructed according to eq 20. The joint confidence region has the advantage of accounting for cross-correlation among uncertain variables. The joint confidence region for the parameters of the pellet reactor with a confidence of R ) 0.99 is Figure 8b. Section 4 will illustrate how to draw relevant samples from the hypercube or hyperellipsoid model of the uncertain space.

|Θj - Θ ˆ j| e Fjj1/2tn-m,1-R/2

(19)

ˆ ) e mFm,n-m,1-R (Θ - Θ ˆ )TF-1(Θ - Θ

(20)

T

e e with s2 ) n-m

(18)

The covariance matrix, Fθˆ , helps construct probability density functions for uncertain parameter realizations. For flexible design, it is advantageous to limit the range of infinite uncertain space according to a tolerable confidence level. The scalar confidence or R-leVel defines a cutoff point outside which the risk associated with highly unlikely parameter realization is deemed tolerable. After specifying appropriate confidence levels (like R ) 0.99), eq 19 delineates an indiVidual confidence region for each parameter, Θj, according to its variance, Fjj, and the outcome of a student t-distribution. When applying eq 19 to all uncertain parameters, a hypercube centered at the nominal point and with specific parameter ranges emerges, as given in Figure

3.7. Result for the Pellet Reactor. The variance-covariance matrix of three estimated parameters, the diffusivity of the reactor, DA; the reaction constant, Kref; and the activation energy, E, and the mass transfer coefficient, hs, is summarized in eq 21.

[

0.37 1.81 Fθˆ ) 0.48 1.73

1.81 0.014 9.27 0.95

0.48 9.27 12.89 7.35

1.73 0.95 7.35 0.054

]

(21)

Accordingly, the bulk diffusivity of species A, DA; the reaction constant, Kref, and the activation energy, E, for the heterogeneous reaction and the mass transfer coefficient from

8358

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006

Figure 8. Random sampling in the confidence interval and joint confident region: (a) sampling within individual confidence region for each parameter (θ1, θ2); and (b) random sampling in joint confidence region of two parameters θ1, θ2. Table 3. Uncertain Parameters for the Catalytic Pellet Reactor optimal (mean) value DA Kref E bulk gas flowrate

10-5

7.32 × 0.82 × 10-4 (1/s) 5268.62 (J/mol) 0.6 m3/s

(m2/s)

variance (0.71 × 10-5 (m2/s) (0.12 × 10-4 (1/s) (3.11 (J/mol) (10%

the bulk to the catalytic pellet vary between DA ) 7.32 ( 0.71 × 10-5 (m2/s); Kref ) 0.82 ( 0.12 × 10-4 (1/s); E ) 5268.62 ( 3.11 (J/mol); and hs ) 1.52 ( 0.062 (m/s). Table 3 lists the mean value and the standard deviations of these three uncertainty parameters of the distributed model and uncertain inlet flowrate. 3.8. Solution Multiplicity. Figure 9 depicts the convergence trajectory of parameters diffusivity, DA, and reaction constant, Kref, on the error surface. In this case, there are two different solutions; each one was successfully identified when starting from different initial guesses (Kref ) 2.2, DA ) 5.0) and (Kref ) 2.2, DA ) 1.0). In other inversion problems, entire local solution clusters were observed.15 Solution multiplicity is frequently the case in kinetic inversion problems when based on experimental data with unavoidable measurement errors. A naı¨ve approach would unconditionally declare the global minimum the true solution of any inversion problem. However, the global residual error minimum is not guaranteed to represent the true physical phenomena. The selection of the true solution from multiple candidates based on physical principles is a challenge that requires further investigation. For completeness, there exists a family of computational algorithms capable of

Figure 9. Convergence trajectory of inversion dynamic reactor. With experimental errors, solution multiplicity of TKIP is common: (a) shows the convergent path from a starting point (Keff ) 2.2, DA ) 5.0); (b) shows the convergent trajectory from a starting point (Keff ) 2.2, DA ) 1.0).

locating the global minimum as well as all other local minima. These methods include global terrain method,26-29 interval analysis,30 BARON,31 and RBB.1,2 Figure 7 shows the concentration and temperature profiles along the main axis of the catalytic pellet reactor. This case study demonstrates that detailed kinetic and transport properties can be extracted from concentration and temperature data. 4. Design of Distributed Systems under Uncertainty In classical deterministic design optimization, the system parameters are considered at their nominal values only; probabilistic information associated with each parameters is usually ignored. Design decisions such as residence times and cooling rates oVerdesign to account for parameter uncertainty. Unfortunately, overdesign bounds set arbitrarily do not guarantee robust performance. We propose to deploy a more mathematically rigorous design procedure which uses both nominal values as well as the confidence information. As a consequence, flexible design optimization is a probabilistic approach that maximizes expected performance with respect to the uncertainty space described by the confidence regions. Prior research introduced mathematical programming for rigorous design under uncertainty.32-37 The stochastic problem in eqs 22-24 formulates the robust design problem for the catalytic pellet reactor by searching for optimal reactor diameter

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8359 Table 4. Optimal Design under Different Sampling Strategies normal design

stochastic design with individual confidence region

stochastic design with joint confidence region

1.39 31.50 84.53

2.63 45.80 147. 70

2.35 42.72 134.90

length of reactor (m) flowrate of cooling solvent (mol/h) minimum cost ($ 103)

and length, D and L, as well as cooling policy, Fcool, that minimizes expected cost over the entire uncertainty space, Ω. Equality constraints 11-16 represent the coupled reaction and the transport equations. Inequality constraints 23 and 24 ensure minimal required conversion and hot spot temperatures inside the reactor below a critical threshold.

s.t. eqs 11-16 χ g 0.95 T e Tmax ) 815 K

conversion threshold

(23)

maximum admissible temperature (24)

4.1. Probabilistic Description of the Uncertain Space. Problem inversion produces estimates of nominal parameter values, Θ ˆ k, as well as their covariance matrix, FΘˆ . As discussed in Section 3, the covariance matrix and appropriate confidence levels define confidence regions. This subsection will construct models of the uncertain space based on probability density functions and samples9,10,38 To compute the expected performance in objective 22, probabilities of parameter realizations need to be computed. In the first method, Monte Carlo samples39 are drawn for each parameter’s individual confidence regions (ICRs). The probability of the outcome for each parameter, pk(Θk) is evaluated according to its mean and its variance. The scenario probability for the entire uncertain parameter vector ΘS follows from eq 25.

Ps(Θ) )

pk(Θk) ∏ k∈M

(25)

The second sampling strategy used the joint confidence region (JCR) to address parameter cross-correlations, as shown in Figure 8. The crude probability for the entire parameter vector ΘkS follows from the joint probability density function of eq 26; the final probability is normalized according to the likelihood of the total ensemble.

Ps )

{

}

1 1 ˆ - ΘS)T‚FΘˆ -1‚(Θ exp - (Θ ˆ - ΘS) 1/2 2 (2π) |FΘˆ | (26) N/2

4.2. Results. The flexible design of the pellet reactor involves choices for the reactor diameter, length, and cooling policy of the reactor for minimum expected cost, while not compromising product quality or violating safety constraints such as a maximum permissible hot-spot temperature. The cheapest design results from considering only nominal values. However, slight parameter variations would violate product and safety constraints, inequalities 23 and 24.

The most robust reactor with flexibility against different uncertainty sources is given in Table 4. The expected annualized cost is $134.9 × 103 for the diameter 0.5 m, where the optimized length is 2.35 m and the necessary cooling flow is 42.7 mol/h; the average conversion of the optimal design is found to be 0.985. Note that the design constraint enforces the conversion to be at least 95%, even in worst-case scenarios. The maximum expected temperature is 814.5 K. In this case study, a more conservative design is obtained with the individual confidence region model, perhaps because of the larger area covered by the uncertain space using ICR. The optimal design with JCR is 9.4% cheaper than the flexible design using ICR. However, the difference between two models for uncertain space is the same of the order variance of the uncertainty sources. General guidelines for determining which of the two metrics to use require additional studies. All computational results were obtained using 30 samples inside the three-dimensional confidence regions, requiring no more than 17 optimization iterations in less than 2 CPU h. The solution of the transport equations for all scenarios in each iteration consumes the largest amount of CPU time. 5. Conclusions and Significance This article introduces a mathematical programming approach for the inversion of coupled reaction and transport problems in distributed problem domains. In addition, robust design of distributed reactors was addressed based on uncertainty information extracted from the problem inversion. The presented methodology for transport and kinetic inVersion problems (TKIPs) simultaneously determines distributed reaction and transport coefficients. We have proposed a numerical optimization approach to quantify the unknown transport and reaction parameters from limited experimental data. The methodology uses advanced optimization techniques such as inexact trust region optimization, sparse linear algebra, and response maps for the solution of the PDE constrained nonlinear mathematical program. The example of a catalytic tubular reactor demonstrates a realistic application for determining unknown parameters as well as their probability to perform optimal design under uncertainty. We found that the individual confidence region gives a more conservative design than the joint confidence region in the catalytic example. Acknowledgment Financial support for the undergraduate research experience at the University of Illinois at Chicago for C.X. during the summer of 2005 was provided by an REU supplement of Grant NSF-DMI 0328134-REU and the NSF-EEC-REU 0453432. Partial support for this work was provided by NSF-DMI 0328134. Nomenclature aP, aE, aW, aS, aN ) The coefficients of the discretized transport model C1, C2, C3 ) price CA ) concentration of reactant A, mol/m3

8360

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006

C ˜ A ) scaled pellet surface conditions Cp ) heat capacity of gas reactant, J/(kg‚K) DA ) diffusion coefficient, m2/s E ) Arrhenius activation energy, J/mol ke ) heat conductivity coefficient of the bed, J/(m‚s‚K) Keff ) reaction rate coefficient, 1/s n ) total number of experimental observations m ) number of parameters M ) total number of scenarios r ) radial position along the pellet, m R ) scaled pellet radius R0 ) pellet radius, m rA ) overall rate of reaction, mol/(s‚m3) S ) source term in transport equations T ) pellet temperature, K V ) total domain Vi ) control volume i in domain V u,V ) velocity, m/s ∆H ) heat of reaction, J/mol Greek Symbols β ) dimensionless variable in pellet equation γ ) dimensionless variable in pellet equation ζ ) scaled pellet temperature η ) overall effectiveness factor θ ) unknown parameters Fb ) density of the bed, kg/m3 Fg ) density of reactant gas, kg/m3 Φ ) Thiele modulus Literature Cited (1) Esposito, W. R.; Floudas, C. A. Global Optimization for the Parameter Estimation of Differential-Algebraic Systems. Ind. Eng. Chem. Res. 2000, 39, 1291. (2) Esposito, W. R.; Floudas, C. A. Parameter Estimation in Nonlinear Algebraic Models via Global Optimization. Comput. Chem. Eng. 2002, 22, 213. (3) 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. (4) 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. (5) Biegler, L.; Ghattas, O.; Heinkenschloss, M.; van Bloemen Waanders, B. Large-Scale PDE-Constrained Optimization. Lecture Notes in Computational Science and Engineering; Springer-Verlag: New York, 2003; Vol. 30. (6) Ozisik, M. N.; Orlande, H. R. B. InVerse Heat Transfer; Taylor and Francis: New York, 2000. (7) Borggaard, J.; Burns, J. A PDE Sensitivity Equation Method for Optimal Aerodynamic Design. J. Comput. Phys. 1997, 136, 366. (8) Roy, R.; Sevick-Muraca, E. M. Truncated Newton’s Optimization Scheme for Absorption and Fluorescence Optical Tomography. Part I. Theory and Formulation. Opt. Express 1999, 4, 353. (9) Rooney, W. C.; Biegler, L. T. Incorporating Joint Confidence Regions into Design under Uncertainty. Comput. Chem. Eng. 1999, 23, 1563. (10) Rooney, W. C.; Biegler, L. T. Design for Model Parameter Uncertainty using Nonlinear Confidence Regions. AIChE J. 2001, 47, 1794. (11) Patankar S. V. Numerical Heat Transfer and Fluid Flow; Hemisphere Publishing Corporation: Washington, DC, 1980. (12) Leis, J. R.; Kramer, M. A. Sensitivity Analysis of Systems of Differential and Algebraic Equations. Comput. Chem. Eng. 1985, 9, 93. (13) Li, S.; Petzold, L. R. Software and Algorithms for Sensitivity Analysis of Large-Scale Differential Algebraic Systems. J. Comput. Appl. Math. 2000, 125, 131.

(14) Ozyurt, D. B.; Barton, P. I. Large-Scale Dynamic Optimization Using the Directional Second-Order Adjoint Method. Ind. Eng. Chem. Res. 2005, 44, 1804. (15) Tang, W.; Zhang, L.; Linninger, A. A.; Tranter, R.; Brezinsky, K. Solving Kinetic Inversion Problems via a Physically Bounded GaussNewton (PGN) Method. Ind. Eng. Chem. Res. 2005, 44, 3626. (16) Carter, R. G. On the Global Convergence of Trust Region Algorithms Using Inexact Gradient Information. SIAM J. Numer. Anal. 1991, 28, 251. (17) Heinkenschloss, M.; Vicente, L. N. Analysis of Inexact TrustRegion SQP Algorithms. SIAM J. Optim. 2001, 12, 283. (18) Zhang, L.; Linninger, A. A. Discovery of Transport and Reaction Properties in Distributed Systems. AIChE J. 2006, submitted for publication. (19) Conn, A. R.; Gould, N. I. M.; Toint, P. L. Trust Region Methods. SIAM: Philadelphia, PA, 2000. (20) Biegler, L. T.; Grossmann, I. E.; Westerberg, A. W. Systematic Methods of Chemical Process Design; Prentice Hall PTR: Upper Saddle River, NJ, 1997. (21) Weisz, P. B.; Hicks, J. S. The Behavior of Porous Catalyst Particles in View of Internal Mass and Heat Diffusion Effects. Chem. Eng. Sci. 1962, 17, 265. (22) Fogler, H. S. Elements of Chemical Reaction Engineering, 3rd ed.; Prentice Hall Inc.: Upper Saddle River, NJ, 1999. (23) Pita, J.; Balakotaiah, V.; Luss, D. Thermoflow Multiplicity in a Packed-Bed Reactor: Condition and Cooling Effects. AIChE J. 1989, 35, 373. (24) Villadsen, J. V.; Michelsen, M. L. Solution of Differential Equations by Polynomial Approximation; Prentice Hall: New York, 1978. (25) Bard, Y. Non-Linear Parameter Estimation; Wiley: New York, 1979. (26) Lucia, A.; Fang, Y. Global Terrain Methods. Comput. Chem. Eng. 2002, 26, 529. (27) Lucia, A.; Fang, Y. Multivariable Terrain Methods. AIChE J. 2003, 49, 2553. (28) Lucia, A.; DiMaggio, P. A.; Depa, P. Funneling Algorithms for Multiscale Optimization on Rugged Terrains. Ind. Eng. Chem. Res. 2004, 43, 3770. (29) Lucia, A.; Xu, J.; Layn, K. M. Nonconvex Process Optimization. Comput. Chem. Eng. 1996, 20, 1375. (30) Kearfott, R. B. Rigorous Global Search: Continuous; Kluwer Academic Publishers: Dordrecht, Netherlands, 1996. (31) Sahinidis, V. BARON: A General Purpose Global Optimization Software Package. J. Global Optim. 1996, 8, 201. (32) Halemaane, K. P.; Grossmann, I. E. Optimal Process Design under Uncertainty. AIChE J. 1983, 29, 425. (33) Swaney, R. E.; Grossmann, I. E. An Index for Operational Flexibility in Chemical Process Design. Part I: Formulation and Theory. AIChE J. 1985, 31, 621. (34) Swaney, R. E.; Grossmann, I. E. An Index for Operational Flexibility in Chemical Process Design. Part II: Computational Algorithms. AIChE J. 1985, 31, 631. (35) Pistikopoulos, E. N.; Mazzuchi, T. A. A Novel Flexibility Analysis Approach for Processes with Stochastic Parameters. Comput. Chem. Eng. 1990, 991. (36) Pistikopoulos, E. N.; Grossmann, I. E. Optimal Retrofit Design for Improving Process Flexibility in Linear System. Comput. Chem. Eng. 1988, 12, 719. (37) Pistikopoulos, E. N.; Grossmann, I. E. Stochastic Optimization of Flexibility in Retrofit Design of Linear System. Comput. Chem. Eng. 1988, 12, 1215. (38) Chakraborty, A.; Linninger, A. A. Plant-Wide Waste Management. 2. Decision Making under Uncertainty. Ind. Eng. Chem. Res. 2003, 42, 357. (39) Niederreiter, H. Random Number Generation and Quasi-Monte Carlo Methods; SIAM: Philadelphia, PA, 1992.

ReceiVed for reView January 18, 2006 ReVised manuscript receiVed June 6, 2006 Accepted June 26, 2006 IE060082L