Model-Based Experimental Analysis of a Fixed ... - ACS Publications

Ind. Eng. Chem. Res. 2009, 48, 5165–5176

5165

KINETICS, CATALYSIS, AND REACTION ENGINEERING Model-Based Experimental Analysis of a Fixed-Bed Reactor for Catalytic SO2 Oxidation J. C. Scho¨neberger,* H. Arellano-Garcia, and G. Wozny Chair of Process Dynamics and Operation, Berlin Institute of Technology, Sekr. KWT-9, Strasse des 17. Juni 135, D-10623 Berlin, Germany

S. Ko¨rkel DFG Research Center MATHEON, Rudower Chaussee 25, D-12489 Berlin, Germany

H. Thielert ThyssenKrupp Uhde GmbH, Friedrich-Uhde-Strasse 15, D-44141 Dortmund, Germany

In this work, a systematic approach based on nonlinear experimental design as an efficient tool for the validation of kinetic models is proposed which has the objective of computing experimental layouts, setups, and controls in order to optimize the statistical reliability of parameter estimates from the resulting experimental data. Furthermore, because the optimum experimental design is dependent on the values of the parameters, a sequential design is proposed. This represents an interaction between proposed experiments (e.g., inlet temperature and concentration), their evaluation, parameter estimation, and new design information such as reactor layout, sampling points, and temperature measurement positions. This follows an iterative procedure until all parameters achieve a prescribed statistical quality. Because of the fact that the overlaying optimization problem contains a highly nonlinear objective function, a hybrid optimization framework is also proposed to overcome the problem of local minima. To show the efficiency of the developed framework, it has been applied to the contact reaction. This reaction is a basic step in the sulfuric acid production and is carried out in a catalytic fixed-bed reactor called the contact horde. 1. Introduction Mathematical models are widely used to represent individual unit operations, as well as entire processes and, consequently, their behavior. They can be applied for the simulation, process design, and/or improvement of existing understanding about specific phenomenon. On the other hand, accurate and reliable representation of the basic phenomena strongly relies upon model structure and the parameter values, which are to be adjusted to correspond with the real process. Therefore, experimental data are required to estimate the model parameters and, thus, validate the developed model. However, the estimation error is dependent on the design of experiments. In particular, for nonlinear systems, the decision on the experiments to be executed is decisive because of cost and time. This issue represents a compromise between experimental attempt and data quality. Methods for the determination of the data quality (i.e., the model parameter accuracy) were first developed for linear systems (see, e.g., ref 1) and then expanded to nonlinear problems.2,3 The descriptiveness of these properties allows the so-called “experimental design” to perform experiments that improve these characteristics.4,5 The application of these methods is especially useful when experiments are aligned with great financial and time effort. In this work, the methods of nonlinear optimal experimental design are applied to the process of sulfur dioxide oxidation. Because sulfuric acid is, worldwide, one of the most produced * To whom correspondence should be addressed. Fax: +49 30 31426915. E-mail address: [email protected].

chemicals,6 the reaction has been the subject of many experimental studies. This is also the reason why plenty of catalysts are commercially available and there are still new catalysts entering the market. On the other hand, the experimental effort for new catalyst materials represents an expensive and challenging task, because high temperatures are required for this particular reaction and corrosive and toxic gases are involved, which not only form sulfuric acid but also corrode almost every type of steel upon condensing. This makes an optimal experimental design necessary. In a first step, the experimental setup can be designed by applying the method of dynamic experimental design to a steady-state problem. Instead of time, space is considered to be an independent variable. Experiments are designed with the position of a measurement left as a free variable. With repetition of this procedure, the most desired measurement positions can be determined. The optimal control factors (e.g., inlet temperature, pressure, and concentrations) can then be calculated for each experiment, so as to reach maximal model parameter accuracy (i.e., maximal model quality). Finally, it is possible to determine the number of experiments necessary to achieve a desired or required parameter accuracy. The remainder of the paper is structured as follows. First, the sulfuric acid production process is introduced noting the importance of an accurate model for the sulfur dioxide conversion. The model used for the fixed-bed reactor is derived including its solution techniques and gradient calculation. Afterward, the approach employed for experimental design is presented, followed by the introduction of the developed hybrid optimization framework.

10.1021/ie801288d CCC: $40.75  2009 American Chemical Society Published on Web 05/11/2009

5166

Ind. Eng. Chem. Res., Vol. 48, No. 11, 2009

Figure 1. Flowsheet of a simplified sulfuric acid production process (ChemCAD).

the inlet condition changes in the contact tower and new available catalysts, the reaction rate models must be validated further. 3. Reactor Model

Figure 2. Control volume for a fixed-bed plug-flow reactor.

Finally, the results for some scenarios of the experimental procedure are given and discussed in detail. 2. Process Description All established processes for the production of sulfuric acid operate with the intermediate product sulfur dioxide. This can be produced by the combustion of brimstone, sulfur ore, or hydrogen sulfide. The entire process is depicted in Figure 1. The further oxidation after the combustion to sulfur trioxide occurs in fixed-bed reactors and is catalyzed by a vanadium pentoxide-containing catalyst. Because of the strong heat production, the reaction beds are separated allowing a quenching with air. They are allocated together in the contact tower. The sulfur dioxide is converted to sulfur trioxide, which can then be absorbed in sulfuric acid, increasing its concentration. The input of makeup-water allows the steady-state production of sulfuric acid with a desired concentration. The main task concerning the design of the process considered is the determination of the number of beds, their length, and the amount of quench air. With these steps, on the one hand, the conversion of the sulfur dioxide is determined, which is important to hold waste gas directives. On the other hand, the main stream amount is fixed, which leads to smaller or larger apparatus size. Thus, to determine an optimal design and fulfill the restrictions, a reliable reactor model is essential. The use of nonsuitable models in the design phase may result in critical situations, such as, for example, the destruction of the catalyst or the violation of toxic gas concentrations in the waste gas outlet. The present work is part of an ongoing industrial cooperation project, in which new operation conditions are to be validated for wet catalysis sulfuric acid processes. However, because of

To describe the behavior of the process unit, a model for a fixed-bed plug-flow reactor has been used. The control volume is shown in Figure 2. The steady-state model equations include the component balances in eq 1, where the concentrations are replaced by the molar flow rates (Fi). The volume of the gas phase is given as a function of the reactor diameter (D), the catalyst density (FK), and the void fraction (ε). The reaction rate for each component i is calculated as the product of its stoichiometry coefficient νi and the overall reaction rate r˙, as follows: dFi π ) νir˙(1 - ε)FK D2 dz 4

()

(1)

The energy balance in eq 2 is not only a strong function of the reaction rate, but also of the temperature dependency of the heat of reaction ∆HRx and the molar heat capacity cP,i of the components: -r∆HRx(1 - ε)FK(π/4)D2 dT ) dz (FicP,i)

(2)

∑

The pressure drop along the reactor is modeled using the Ergun equation, given in eq 3. This is necessary because the pressure drop in the fixed bed cannot be neglected a priori, and the overall reaction rate depends on the pressure:

∑

( )

(FiMi) 1 - ε dP )· dz FGas(π/4)D2DP ε3

[

∑

]

(FiMi) 150(1 - ε)η (3) + 1.75 DP (π/4)D2

where DP is the diameter of the catalyst particle, FGas the density of the gas phase, η the viscosity of the gas phase, and Mi the molar weight of component i. Using the laws of stoichiometry, the component balances in eq 1 can be reformulated to eq 4,


by considering the conversion of sulfur dioxide (as X ) 0 - FSO2)/FSO ) over the reactor length. 2

0 (FSO 2

r˙(1 - ε)FK(π/4)D2 dX ) 0 dz FSO 2

(4)

The reaction that occurs over the vanadium catalyst, 1 SO2+ O2 a SO3 2

(5)

can be described in eq 6, with the kinetic equation given by Eklund,7 as a function of the component partial pressures Pi, where KP is the equilibrium constant related to the reactant partial pressures. Equation 6 can be reformulated as a function

of the conversion with the initial component ratio Θi:

r˙ ) k

1-X ΘSO3 + X

[∑

ΘO2 - νO2X Θi - νO2X

P-

(

)]

ΘSO3 + X

2

(1 - X)KP

(7)

The parameters to be estimated during the model identification and validation procedure are the kinetic constants k1, k2, and k3, whereas the equilibrium constant KP can be calculated from property data and is considered to be known. However, note that the structure of eq 6 is an empirical one. It corresponds to the generic equation for power-law reaction rates: r˙ ) k+

∑

Pfi i

- k-

∑

Pbi i

)k

(∑

1 Pfi i KP

∑

)

Pbi i

(8)

The exponent values found by Eklund are given in Table 1. The same coefficients have been reported by Fogler.8 Carberry reported a slightly different reaction rate.9 Nevertheless, they all are empirical and should be treated with care. Here, the model published by Eklund (eq 6) is considered to be the true model. The complete set of equations is discretized with the orthogonal collocation method using three collocation points and five intervals, which shows a good convergence behavior.10 The model parameters, as well as the resulting algebraic equation system, which is composed of 60 equations and variables, are summarized in Appendix A. To obtain better convergence behavior, the state variables are normalized by scaling them on their inlet values, except for the conversion, which is, by nature, in the range of 0-1. Parameters are scaled in a logarithmic manner by including them in the exponential term in eq 6. The system is solved using a Newton-Ralphson step. The sensitivities of the state variables x ) [X,T,P]T concerning the parameters P ) [k1,k2,k3]T can then be calculated using the Jacobian of the equation system (Jx) and the internal numerical differentiation as follows: ∂x ) -Jx-1JP ∂P

(9)

where JP is the Jacobian of the equation system with respect to the parameters. Equation 9 gives the exact derivation of the numerical solution, with respect to the parameters.11 The sensitivities are used to calculate the covariance matrix, which will then be used in the objective function of the later introduced optimization problem in eq 12:

C)

T

5167

-1

[( ∂P∂x ) C ( ∂P∂x )] -1

(10)

M

If the state variables x are measured directly, the covariance matrixes C can be calculated from eq 10, where the sensitivities are weighted with the inverse variances of the measured variables CM.12 Therefore, independent-additive, normally distributed measurement errors are assumed. However, note that CM has a large influence on the result of the optimal experimental design. In this work, three state variables are measured: conversion (X), temperature (T), and pressure (P). The measurement accuracy of these variables differs considerably. The errors resulting from the measurement procedure can be taken from the devices, such as gas chromatograph, resistance thermometers, and pressure gauges. In addition to these measurement errors, a value for the model inaccuracy must be added, which, in many cases, is larger than the measurement errors. This can be, in the case of the temperature measurement, for instance, a penalty for the assumption of an “adiabatic reactor”. A real reactor is not adiabatic, and therefore, this assumption leads to model inaccuracy. Both the total standard deviations and the covariance matrix of the measured variables are given in eq 11. The high error for the temperature measurement is a result of the uncertainty of the related properties, i.e., the heat capacities of the components and the assumption of an adiabatic reactor.

(

σX ) 2.5%, CM )

σT ) 5 K, 2

σx

2

σT

σP2

)(

σP ) 200 Pa,

)

0 6.25 × 10-4 0 ) (11) 0 25 0 0 0 40000

The profiles for the state variables are depicted on the top of Figure 3. They show the typical behavior of an equilibrium-limited reaction. After a certain distance z in the reactor, the reaction stops because the equilibrium state is reached. An additional length represents a waste of catalyst and an unnecessary increase in the pressure drop. However, it can also be used as a safety margin to ensure the maximal possible conversion (i.e., equilibrium conversion). On the bottom of Figure 3, the sensitivities of the state variables, with regard to the parameters, are shown. A change in the kinetic constants ki does not lead to a change in the equilibrium conversion and temperature, because the term in brackets in eq 7 approaches zero when reaching chemical equilibrium, and the velocity constant k does not influence the reaction rate anymore. Thus, the sensitivity trajectories approach zero when the equilibrium state is reached. The pressure drop is a weak function of the reaction rate. The pressure drop is increased when the forward reaction is faster in the front part of the reactor. Nevertheless, a change of 0.1% of parameter k1 leads to a change of 0.04 Pa in the pressure drop. This is beyond the detectable range of the pressure measurement. Therefore, the pressure measurement is not used in the parameter identification procedure. 4. Nonlinear Optimal Experimental Design For the nonlinear optimal experimental design, a nonlinear constrained optimization problem must be solved. In this case, the optimization variables are the experimental design variables, namely, the inlet conditions, and the measurement positions for

5168


Table 1. Coefficients of the Power Law Approach Taken from ref 8 k1 [K] -97782

k2 [1/ln(K)] k3 [-] fN2 fSO2 fO2 -110.1

848.5

0

0.5

1

fSO3

bN2

bSO2

-0.5

0

-1.5

bO2 bSO3 0

1.5

Table 2. Bounds and Optimal Values of the Design Variables inlet condition

Reynolds number, Re

xSO2 [mol/mol]

xO2 [mol/mol]

Tin [°C]

lower bound upper bound reference value optimal value

300 600 300 300

0.002 0.080 0.080 0.080

0.002 0.100 0.100 0.100

400 500 500 447

the state variables. The resulting optimization problem is defined as follows: min φ(C) q,w

(12)

s.t. F(x˙, x, p, q, w, xs) ) 0 G(x, p, q, w, xs) ) 0 a e Ψ(q, w) e b with

{

trace(C) ∨ NP φ(C) ) D: det(C)1/NP∨ E: max(eig(C)) A:

where q represents the experimental conditions (i.e., inlet temperature, concentrations, etc.); w the measurement position weights, which indicate if a measurement is taken (w ) 1) or not (w ) 0); and xs the discrete state variables. F represents the model equations, and G is the nonlinear constraints with the state variables x and the parameters p, while Ψ defines the bounds of the optimization variables. NP denotes the number of parameters to be estimated. The objective function φ in eq 12 is a function of the covariance matrix C of the parameter estimation problem. There are different criteria that can be minimized (see, e.g., ref 2). In Figure 4, these criteria are interpreted graphically for the estimation of two parameters, making use of a linearization of the confidence region. The inner point gives the estimated value of the parameters p1 and p2. The ellipsoid surrounding this point gives the linearized confidence region for a certain probability. According to this, the real values of the parameters are inside the ellipsoid with a probability of, e.g., R ) 95%. The objective of the experimental design is to shrink this ellipsoid to an acceptable size. Therefore, the D-criterion minimizes the determinant of C. This can be interpreted as the surface area of the ellipsoids. This criterion can lead to very large stretched confidence regions with a high correlation of parameters, but to a small surface area. Because of the structure of the approach for the velocity constant k in eq 6, the parameters k1-k3 are strong correlated. Thus, working with the D-criterion will not lead to an improvement in model accuracy. The A-criterion minimizes the area of the circle formed with the average main axis length of the ellipsoid. It is calculated as the trace of C. The advantage of this criterion is a low computational effort, because the trace calculation is quite fast, compared to the other criteria. The E-criterion represents the length of the longest main axis of the ellipsoid. Its value is given from the largest eigenvalue of C. This includes an eigenvalue calculation and a discontinuity. In any case, the result is similar to the A-criterion for strong

correlated parameters. Because of these facts, in this work, the A-criterion is chosen, i.e., minimizing the trace of C. The optimization problem is solved with an approach based on SQP with a tailored derivative computation. For the derivative calculation procedure of the different objective functions, we refer to the work in ref 13. Here, the commercial solvers SNOPT and NPSOL14 from Tomlab are used in a Matlab environment. It has been found that NPSOL is more robust and fasterconverging. The model equation system is solved sequentially and not simultaneously, which leads to an optimization problem that contains few free variables but has a dense Jacobian. In the open literature, NPSOL is recommended for such problems (see references given in ref 15). 5. Hybrid Optimization Framework The objective function of the considered optimization problem in eq 12 is highly nonlinear. Furthermore, the functions of the covariance matrix exhibit several local minima within the considered region (see Figure 9). However, the surface of the objective function shows a trend toward a certain direction. This tendency can be exploited using the stochastic search algorithm PSO (particle swarm optimization).16 We refer to ref 17 for a detailed description of the PSO algorithm. It simulates particles that move over the surface of the objective function. The particles share their actual and historical information (i.e., the values of the objective function) and change the direction of their movement toward the best particle position. The main problem with this algorithm is that it is very difficult to find any minima, local or global, because it does not use gradient information. This leads to huge convergence times. To overcome this problem, the algorithm is combined with a gradient-based optimization algorithm. The key idea is that each particle should move to the next optimum, using gradient information, before it shares information with the other particles (see Figure 5). Because of the particle movement, the initial points of the gradient-based algorithm can be in the constrained region. This requires a robust algorithm such as NPSOL that can handle these types of problems.18 An illustrative example that highlights the advantage of the proposed algorithm can be found in Appendix B. 6. Results 6.1. Experimental Setup. To determine the best experimental setup, a specific scenario is given. As shown in Figure 6, the experimental reactor is divided into five sections. In all of these sections, the temperature is measurable. In addition, the concentration of the reactor outlet can be measured. Because concentration measurements are expensive and time-consuming (e.g., gas chromatography), only two additional concentration measurements are allowed. The positions at which positions these measurements are the most informative are yet to be determined. To answer this question, 10 optimal experiments are designed, leaving the measurement weights w of the positions 1-4 as free variables in eq 12. Furthermore, a constraint must be included that explicitly ensures that the sum of these variables must be equal to 2. The optimal positions are related to the experiment for which they have been determined. This is the reason why a set of experiments is designed. In nine of these, positions 3 and 4 are determined to be the best (see Figure 6). This can be interpreted as follows: the concentration profile can be “moved” through the reactor


5169

Figure 3. State variables conversion (X), temperature (T), and pressure drop (dP), and their sensitivities, for a sample test experiment.

Figure 4. Graphical interpretation of the optimization criteria (linearization).

with the other design variables (i.e., the inlet conditions). The latter positions have an advantage, because of the final concentration measurement, which is fixed at the reactor outlet. This leads to optimal experiments with a concentration change close to the reactor outlet. 6.2. Optimal Single Experiment. In this section, one optimal experiment is designed. Here, the design variables are the inlet conditions of the experimental reactor (i.e., the Reynolds number (Re), the inlet concentrations (xSO2), and xO2, as well as the inlet temperature Tin). They can be influenced with the controls shown in Figure 6. The optimal single experiment is compared to a reference experiment that has been derived intuitively. Based on the linear experimental design, it can be concluded that the optimum of a constrained optimization problem must lie on the bound of the constraints. The optimal bounds can be identified easily: the temperature and the educts concentrations are to be maximal to obtain a high conversion. Moreover, the total flow rate represented by the Re parameter must be minimal in order to maximize the temperature increase. Table 2 summarizes the derived

Figure 5. Hybrid optimization procedure.

inlet conditions, which concern the Reynolds number Re, the mole fractions xi of sulfur dioxide and oxygen, the inlet temperature Tin as well as their bounds. As shown in Figure 7, because of the nonlinearity of the system, the reference experiment is not optimal. The conversion reaches the equilibrium very fast, and no further parameter information can be extracted from the flat profiles. The optimal single experiment (i.e., the solution of the optimization problem given in eq 12), results in a conversion profile, in which the equilibrium conversion is reached at

5170


Figure 6. Experimental setup of the PFR system. Legend: FCR, flow control and registration; PCR, pressure control and registration; TCR, temperature control and registration; TIR, temperature indication and registration; QIR, quality indication and registration (conversion X).

Figure 7. Trajectories and confidence region of an intuitive experiment (top) versus an optimal single experiment (bottom).

the end of the reactor allowing an informative concentration measurement in the position 5 (see Figure 6). Moreover, Figure 7 shows the profiles of the two experimentssthe intuitive and the single optimalsas well as the confidence ellipsoid for their corresponding parameters k1 and k2. The confidence region is calculated based on the covariance matrix C with a multivariable normal distribution sampling method. For a better comparison, the parameters are scaled to unity and only the two parameters (k1 and k2) are considered. The two additional confidence regions (k1-k3 and k2-k3) appear similar. With the optimal single experiment, the spreading of the parameter values can, nevertheless, be reduced considerably. Note that both the reference and the optimal experiment may cause problems in the experimental implementation. In both cases, the temperature inside the reactor crosses the bound value

of T ) 600 °C. This is the highest temperature allowed for the vanadium pentoxide catalyst.19 Thus, in this work, the maximal temperature in the reactor is limited to Tmax ) 550 °C. This constraint is included in the optimization problem, which leads to a new optimal single experiment. The resulting trajectories and the confidence region are given in Figure 8. The new optimal decision variables are given in Table 3, and the resulting values for the objective function and the parameter spreading in Table 4. The last column in Table 4 represents the standard deviation including the correlation of the parameters. It can be concluded that, with only one single optimal experiment, no acceptable parameter accuracy can still be reached. Thus, additional experiments must be designed. In the sequel of this work, the designed experiments are restricted to the temperature bound. Results for an unrestricted experimental design have been reported in ref 20.


5171

Figure 8. Optimal single experiment with a restricted reactor temperature. Table 3. Optimal Values of the Design Variables in the Temperature-Constrained Single Experiment inlet condition

Reynolds number, Re

xSO2 [mol/mol]

xO2 [mol/mol]

Tin [°C]

optimal value

300

0.080

0.100

421

Table 4. Optimization Criteria and the Related Scattering of the Parameter Values experiment

A

D

E

σk1 [%]

σk2 [%]

σk3 [%]

σmax [%]

reference optimal constrained

0.9342 0.0021 0.0058

2.00 × 10-6 9.08 × 10-8 1.36 × 10-7

2.803 0.006 0.017

92.1 4.2 6.8

98.7 4.7 8.0

98.9 4.7 7.9

167.4 7.9 13.2

6.3. Multiple Experiments. Based on the presented solution strategies, it is now possible to design a set of multiple experiments. For this propose, a sequential procedure is recommended, because the estimated values of the parameters can change with each performed experiment, therefore leading to a different optimal experiment in the next step. In this study, the parameter values are assumed to be constant from experiment to experiment, because of the fact that no real experiments are performed so far. Anyhow, in the case of real performed experiments, there will be a change of parameter values and a sequential proceeding becomes necessary. In Figure 9, the surface of the objective function for the fifth experiment is plotted for the two decision variables: the inlet temperature Tin and the Reynolds number Re. The striped region in the image on the right in Figure 9 represents the infeasible combinations of these variables. The problem concerning the local minima can be seen on the left side of Figure 9, where the surface is plotted. The objective function is highly nonconvex. Thus, a tailored

Figure 10. Progress of the optimization criteria according to the number of experiments.

optimization algorithm is required. To overcome the issue of ending up in a local optimum, the proposed hybrid optimization algorithm is applied. In this case study, five particles and five iterations are used, which leads to 25 calls of the optimization problem. As a goal for the planning procedure, certain parameter accuracy can be assigned. Moreover, based on the sequential experimental design, it can be determined how many experiments are required to reach the requested accuracy. This gives an idea of the load of experimental work that must be performed to reach the given objective. However, this is only an approximation. The optimal experiments and parameter accuracies will change with each new parameter estimation, which must be done after each performed experiment, because the parameters converge to their “true” values. Besides, based on this procedure, the evolution of the confidence region can be

Figure 9. Objective function (A-criterion) for the fifth designed experiment with respect to the two decision variables by fixed inlet concentrations at xSO2 ) 0.05 and xO2 ) 0.1.

5172


Figure 11. Improvement of the confidence region with respect to the number of experiments.

Figure 12. First 15 optimal experiment settings and their trajectories.

Figure 13. Analytical calculated confidence regions for the parameter combinations (k1, k2, or k3).

monitored, and thus, it can be decided from which point on, further experiments will not increase parameter accuracy significantly. In Figure 10, the progress of the parameter value accuracies is shown and the horizontal line in that figure delineates the 2% standard deviation. Figure 11 visualizes the shrinking of the confidence ellipsoid for the parameters k1 and k2. After 10 experiments, curves that represent the optimization criteria are getting quite flat. This means that no significant improvement of the parameter value accuracies can be reached with more experiments. In Figure 12, one can notice that some high informative experiments are repeated with an increasing number of experiments (e.g., experiments 4, 5, 8, and 12).

These recurring experiments are optimal because they reduce the measurement error in the repeated experiment. In Figure 12, the first 15 optimal experiment settings and their corresponding trajectories are plotted. After 14 experiments, the standard deviation including the correlation of the parameters is lower than 2% (see Figure 10). This could also be an objective of the experimental design procedure. Nevertheless, lower deviation values imply higher experimental effort. On the other hand, this value is not satisfactory, in particular when it is compared to Figure 3, where a change of only 1% of k3 results in a change of up to 3% in the conversion. One way to overcome this problem and to reduce the model inaccuracy is to change the model structure. To


select the “responsible” parameter, the analytically derived confidence ellipsoids for the single experiment are calculated with a 95% confidence region, using a function stated in Appendix C (see also Figure 13). The parameters k2 and k3 are highly correlated, leading to the largest ellipsoid in Figure 13. The correlations of k1-k2 and k1-k3 are almost the same. Thus, it can be concluded that k2 or k3 could be eliminated. Mathematically, it makes no difference which variable is eliminated. Anyhow, from a physical point of view, it is recommended to remove the parameter k2, because this reduces the approach for the rate constant k in eq 7 to the well- known Arrhenius equation. 7. Conclusions The application of the methods of nonlinear optimal experimental design to steady-state space-variant reaction systems gives fairly good results and enables new options in the planning of experiments. The developed hybrid optimization framework improves the results of the experimental design, because it helps to overcome the problem of local minima. While a single experiment can be designed intuitively, in the case of multiple experiments, this is not possible. It is shown that the confidence region can be reduced considerably with the optimal designed experiments. Furthermore, control of the parameter accuracy is possible, and thus, an estimation of the experimental effort can be determined. Moreover, restrictions related to the state variables can be introduced. Based on this, it can be assured that the reactor temperature will not exceed a critical value. The results show an interesting change in the decision variables and trajectories of the optimal single experiment. This constraint represents also an important safety aspect, which must be considered in the laboratory work. During the model validation process highly correlated parameters can be identified such that the model structure can be improved in the sense of parameter identifiability, which also means model reliability. The impact of parameters, which are identified with a low accuracy on a process model, is demonstrated with the related sensitivities. From this point of view, desired model accuracy can be defined and the related experimental effort can be estimated. Based on this, the important procedure of model identification and validation can be accomplished in a more systematic way using the proposed approach based on the methods of nonlinear optimal experimental design.

Appendix A. Equation System and Jacobians The differential equations (eqs 2-4) are discretized with the orthogonal collocation on finite elements. In this work, NC ) 3 collocation points are used in NI ) 5 intervals, resulting in an

parameter D L FK ε DP

Table A2. Gas Properties and Further Parameters parameter

value DIPPR 107,22 ideal mixing ηN2, from Perry and Green22 from Perry and Green22 ideal gas heat of formations,22 cP see eq A823

cP η Mi Fgas ∆HRx Kp

algebraic equation system of 60 equations. The mass balance is given exemplarily in eqs A1-A4 for the first interval. F(1) ) X(1)

value 0.1 m 0.5 m 1082 kg/m3 56.7% 6.9 mm

(A1)

F(j) ) A(j, 1:NC + 1)X(1:NC + 1) 2 L r˙(1 - ε)FK(π/4)D NIzNC F0

j ) 2:NC (A2-A4)

SO2

where A and zNCrepresent the collocation matrix A for three collocation points and the coordinate of the last collocation point, respectively (see, e.g., ref 21). In eqs A5, A6, and A7, Jx, the resulting Jacobians Jx, and JP are defined.

[

Jx,I)1 )

[

∂F(1) ∂X(2) ∂F(2) ∂F(2) ∂X(1) ∂X(2) l l ∂F((NC + 1) · 3) ∂F((NC + 1) · 3) ∂X(1) ∂X(2) 1

··· ··· ·

··

···

∂F(1) ∂P(NC + 1) ∂F(2) ∂P(NC + 1)

]

l ∂F((NC + 1) · 3) ∂P(NC + 1) (A5)

]

Jx )

Jx,I)1 0 0 0 0 Jx,I)2 0 0 0 [-1 · · · -1 · · · -1] Jx,I)3 0 [-1 · · · -1 · · · -1] 0 0 Jx,I)4 0 0 [-1 · · · -1 · · · -1] 0 0 0 0 [-1 · · · -1 · · · -1] Jx,I)5 (A6)

The Jacobian versus the state variables Jx shows a block structure, which enables a faster inverse calculation. The elements below the blocks result from the consistency conditions which couple the intervals. They can be eliminated reducing the equation system to 45 equations with 45 variables.

[

]

JP )

∂F(1) ∂F(1) ∂F(1) ∂k1 ∂k2 ∂k3 l l l ∂F((Nc + 1) · NI · 3) ∂F((Nc + 1) · NI · 3) ∂F((Nc + 1) · NI · 3) ∂k1 ∂k2 ∂k3 (A7)

Tables A1 and A2 give all of the used parameter values and their corresponding units with KP (atm-1/2) ) exp

Table A1. Reactor and Catalyst Parameters

5173

( 11852T K - 1.3091)

(A8)

Appendix B. Illustrative Example and Validation of the Hybrid Optimization Algorithm The proposed hybrid algorithm is tested here based on an objective function, f, which features several local minima. The

5174


Figure B1. Test function for the hybrid optimization algorithm (left) and particle positions (right).

Figure B2. Comparison of the hybrid approach with a random seed of initial values.

function in eq B1 with a global optimum at X ) [6.8, 8]T is shown on the left side of Figure B1. f ) -(sin(X1) · X1 · X2 + cos(X2) · X22 + X1 + X2) (B1) On the right side of Figure B1, both the starting positions of the particles (squares) and the detected optima (points) are denoted. Here, five particles and five iterations are used. It can be shown that, although the start positions are distant from the global optimum (squares on the bottom of Figure B1), it can be found anyway by the new framework. In this example, 10 local minima were detected. It is quite evident that, with this procedure, one cannot be assured that the global optimum of the objective function will be found. However, it helps the gradient-based algorithm not to get stuck in a bad (i.e., high value) local minimum. For the sake of comparison, in the case of nonlinear experimental design problems, optimal experiments are calculated with a random seed of 25 initial points taking over always the best optimum. The development of the A-criterion for two approaches is shown in Figure B2. Occasionally, the same optimal experiments are found and the final difference in the A-criterion is only 3 × 10-5. Anyhow, this leads to a mean

increase in the parameter standard deviation of 0.2%, and another experiment would be necessary to reach the desired accuracy. However, the authors are aware of the fact that this is not a final proof of the advantages of the proposed approach. It has been used straightforward in this work to get improved results. A detailed examination of the approach on a set of selected optimization problems is underway.24

Appendix C. Analytic Calculation of the 2D Confidence Ellipsoid The calculation of the two-dimensional (2D) confidence ellipsoid is based on the probability density function (PDF) for a multivariable normal distribution (see eq C1).

R(p) )

√det(C-1) dim(p)/2

(2π)

(

exp -

pTC-1p 2

)

(C1)

The confidence region G can be calculated with the integral in eq C2, assuming an ellipsoidal form, which leads then to the inequality in eq C3.


∫ ∫ √det(C 2π

(

-1

)

exp -

G

(

1 - exp -

)

pTC-1p dp1 dp2 e R 2

)

pTC-1p eR 2

(C2)

(C3)

The ellipsoid now is positioned on the bound of the confidence region and can be calculated as the solution of eq C3. Because of the quadratic nature of the ellipsoid, eq C3 has two solutions for p1. They are given in eq C4. p2 ) -p1

c ( c2

( p21

)

c1 c2 ln(1 - R) 2 c2 c2 c2

(C4)

with C-1 )

[ ] c1 c c c2

To exclude complex solutions of eq C4, p1 must be bounded. The allowed interval is given in eq C5. Following this step, the ellipsoidal can be plotted (see Figure 13).

{

p1 ∈ -

ln(1 - R)c2 c - c1c2 2

;+

ln(1 - R)c2 c2 - c1c2

}

List of Notation Process Variables bi ) exponent in the reaction rate cP,i ) component heat capacity (J/(mol K)) D ) reactor diameter (m) DP ) particle diameter (m) dP ) pressure drop (Pa) dz ) differential length (m) fi ) exponent in the reaction rate F ) total flow rate (mol/s) Fi ) component flow rate(mol/s) Fi0 ) initial composition flow rate (mol/s) ∆HRx ) heat of reaction (J/mol) k, k+, k- ) velocity constants k1, k2, k3 ) constants in the extended Arrhenius equation KP ) equilibrium constant (1/atm1/2) L ) reactor length (m) Mi ) component molar mass (kg/mol) P ) pressure (atm) Pi ) partial pressure (atm) r˙ ) reaction rate (mol/(kg s)) Re ) Reynolds number T ) temperature (K or °C) xi ) component mole fraction X ) conversion of SO2 z ) length coordinate (m) ε ) void fraction η ) dynamic viscosity (Pa s) θi ) initial composition ratio νi ) composition stoichiometric coefficient FGas ) gas-phase density (kg/m3) FK ) catalyst density (kg/m3) Mathematical Symbols A ) collocation matrix a, b ) bounds C ) covariance matrix CM ) covariance matrix of the measurement errors

(C5)

5175

c ) element of the inverted covariance matrix f ) function value F ) (discretized) differential equations G ) algebraic equations JP, Jx ) Jacobians NC ) number of collocation points NI ) number of intervals NP ) number of parameters to be estimated P ) parameter vector p, p1, p2 ) model parameters q ) control w ) weight x ) state variable xS ) discrete state variable zNC ) position of the last collocation point σ ) standard deviation φ ) function of the covariance matrix Ψ ) constraint function

Literature Cited (1) Lazic, Z. R. Design of Experiments in Chemical Engineering; Wiley: New York, 2004. (2) Macchietto, S.; Asprey, S. P. Curvature-based methods for designing optimally informative experiments in multiresponse nonlinear dynamic situations. Ind. Eng. Chem. Res. 2005, 44 (18), 7120. (3) Franceschini, G.; Macchietto, S. Validation of a model for biodiesel production through model-based experiment design. Ind. Eng. Chem. Res. 2007, 46, 220. (4) Ko¨rkel, S.; Bauer, I.; Bock, H. G.; Schlo¨der, J. P. Scientific Computing in Chemical Engineering; Springer: Berlin, 1999. (5) Bardow, A.; Go¨ke, V.; Koss, H.-J.; Marquardt, W. Concentrationdependent diffusion coefficients from a single experiment using modelbased Raman spectroscopy. AIChE J. 2006, 52 (12), 4004. (6) McCoy, M. Facts & Figures. Chem. Eng. News 2003, 81, 27. (7) Eklund, R. B. The Rate of Oxidation of Sulfur Dioxide with a Commercial Vanadium Catalyst; Almqvist & Wiksells Boktr.: Uppsala, Sweden, 1956. (8) Fogler, H. S. Elements of Chemical Reaction Engineering; Pearson International: Upper Saddle River, NJ, 2006. (9) Carberry, J. J. Chemical and Catalytic Reaction Engineering; Dover Publications: Mineola, NY, 2001. (10) Scho¨neberger, J. C.; Arellano-Garcia, H.; Thielert, H.; Wozny, G. An efficient approach to robust simulation of claus processes in coking plants. Comput.-Aided Chem. Eng. 2007, 24, 521. (11) Kudryavtsev, L. D. Implicit Function. In Encyclopedia of Mathematics; Kluwer: Dordrecht, The Netherlands, 1990. (12) Faber, R.; Arellano-Garcia, H.; Li, P.; Wozny, G. An optimization framework for parameter. estimation of large-scale systems. Chem. Eng. Process. 2007, 46, 1085. (13) Bauer, I.; Bock, H. G.; Ko¨rkel, S.; Schlo¨der, J. P. Numerical Methods for Optimum Experimental Design in DAE Systems. J. Comput. Appl. Math. 2000, 120, 1. (14) Gill, P. E.; Murray, W.; Saunders, M. A.; Wright, M. H. User’s Guide for NPSOL 5.0: A Fortran Package for Nonlinear Programming, Technical Report SOL 86-2, Systems Optimization Laboratory, Department of Operations Research, Stanford University, Palo Alto, CA, 1998. (15) Stein, O.; Oldenburg, J.; Marquardt, W. Continuous reformulations of discrete-continuous optimization problems. Comput. Chem. Eng. 2004, 28, 1951. (16) Kennedy, J.; Eberhart, R. Particle swarm optimization. In Proceedings of the IEEE International Conference on Neuronal Networks, Volume 4; IEEE: Piscataway, NJ, 1995; p 1942. (17) Schwaab, M.; Biscaia, E. C.; Monteiro, J. L.; Pinot, J. C. Nonlinear parameter estimation through particle swarm optimization. Chem. Eng. Sci. 2008, 63 (6), 1542. (18) Gill, P. E.; Murray, W.; Saunders, M. A.; Wright, M. H. Some issues in implementing a sequential quadratic programming algorithm. ACM Signum Newsl. 1985, 20, 2. 13. (19) BASF Product Information: High Performance Catalysts for Chemicals; 2005.

5176


(20) Arellano-Garcia, H.; Scho¨neberger, J. C.; Ko¨rkel, S. Chem. Ing. Tech. 2007, 79 (10), 1625. (21) Finlayson, B. A. Nonlinear Analysis in Chemical Engineering; McGraw-Hill: New York, 1980. (22) Perry, R. H.; Green, D. W. Perry’s Chemical Engineers’ Handbook, 7th Ed.; McGraw-Hill: New York, 1998. (23) Duecker, W. W.; West, J. R. The Manufacture of Sulfuric Acid; Reinhold Publishing Corporation: New York, 1959.

(24) Schöneberger, J.; Arellano-Garcia, H.; Wozny, G. Local Optima in Model-Based Optimal Experimental Design. To be submitted to Ind. Eng. Chem. Res.

ReceiVed for reView August 25, 2008 ReVised manuscript receiVed January 6, 2009 Accepted February 11, 2009 IE801288D

Model-Based Experimental Analysis of a Fixed ... - ACS Publications

Recommend Documents