Determination of Operating Regimes for Chemical Reactors Using

3390

Ind. Eng. Chem. Res. 2007, 46, 3390-3399

PROCESS DESIGN AND CONTROL Determination of Operating Regimes for Chemical Reactors Using Incomplete Mechanistic Knowledge and a Mutual Information-Based Dependency Test Aidong Yang,*,† Julian Morris,‡ Gary Montage,† and Elaine Martin† Centre for Process Analytics and Control Technology and School of Chemical Engineering and AdVanced Materials, Newcastle UniVersity, Merz Court, Newcastle upon Tyne NE1 7RU, United Kingdom

A chemical process that involves the interplay between chemical reaction(s) and transport phenomena may be operated in one of several possible operating regimes. In the development of such a process, the operating regime which yields desirable performance at a smaller scale should be identified to help realize successful scale-up of the process. In the past, experimental procedures for identifying operating regimes have been proposed, which utilize the qualitative trends of the response of a process to changes in operating conditions. In this study, an alternative approach which utilizes both experimental data and mathematical modeling is explored. The rationale behind this approach is that an incorrect assumption with regard to an operating regime may result in a false dependency between particular process variables. On the basis of this hypothesis, operating regime determination can be formulated around a dependency test. In this paper, the computation of a dependency measure is based on the mutual information approach as defined in information theory. The methodology proposed is able to handle the situation where mechanistic knowledge such as chemical reaction kinetics is missing. In practice, this information is often not available during the early stages of process development. The proposed methodology is tested and compared with alternative approaches on a simulated toluene nitration process. Its potential in terms of achieving a reliable discrimination between different assumptions with regard to the operating regime with fewer constraints on the experimental requirements is demonstrated. 1. Introduction In the development of a new chemical process it is important to gain an understanding of how the process performs with respect to a wide range of equipment design parameters and operating conditions. When a process involves chemical reactions as well as transport phenomena, such as intra-phase mixing and/or inter-phase mass transfer, the performance of the process is typically determined by a controlling step which places limits on the yield, selectivity, and so forth.1 When the design and operating conditions change, the controlling step may vary, and consequently a number of operating regimes may materialize.2 For example, a heterogeneous reaction system can be operated under either a slow or a fast reaction regime, depending on how rapid the inter-phase mass transfer is in comparison to the rate of the chemical reaction. In the situation where a desirable level of performance has been achieved in smaller scale experiments for a particular operating regime, this regime should be identified and the information utilized when making scale up decisions. This has been discussed in the literature, particularly with respect to mixing.3-5 For example, Atherton3 showed that the effect of agitator speed on overall conversion rates is qualitatively different between fast-reaction and slow-reaction regimes. Consequently the identification of operating regimes is an important issue that needs to be addressed in the early stages of process development. * Corresponding author. Tel.: +44 (0)191-222-5331. Fax: +44(0)191-222-5748. E-mail: [email protected]. † School of Chemical Engineering and Advanced Materials. ‡ Centre for Process Analytics and Control Technology.

Because operating regimes usually correspond to the relative rate of different phenomena, a particular operating regime can be characterized by a specific range of a dimensionless number which denotes the ratio of two characteristic times. For example, the Damkoehler number (Da) of a reactive-mixing process is the ratio of the engulfment (mixing) time constant and the halflife of the chemical reaction. Likewise, the Hatta number (Ha) of a heterogeneous reaction process denotes the comparison of the diffusion time and the reaction time.6,5 When sufficient information about a process is available, such a dimensionless number can be computed, and the identification of the operating regime becomes straightforward. However, this is often not the case when a new process is being developed because there is usually only limited knowledge about the scaled-up process at this stage, especially with respect to the chemical reaction kinetics. In such a situation, an experimental procedure, specific to the type of process being investigated, can be implemented to determine qualitatively the operating regime, although particular types of equipment might be needed to support the determination process.3,5 For example, the determination of the operating regime of a heterogeneous reaction system can be achieved by applying both a baffled turbine agitated reactor and a constant interfacial area cell and by changing operating conditions such as agitation speed and phase volume.3 Johansen and Foss7 proposed a methodology whereby they decomposed a global model of the operation of a chemical process into a number of local models, each of which corresponds to a particular operating regime. Solutions for determining the number of local models (or operating regimes) and for

10.1021/ie0609721 CCC: $37.00 © 2007 American Chemical Society Published on Web 04/11/2007

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3391

identifying each local model were developed using a data-based approach. The research presented in this current paper differs in that the goal is to characterize a chemical reactor utilizing one of a number of possible models that represent different operating regimes. In this study, quantitative alternatives to the existing qualitative approaches for determining the operating regime for a reactor are explored. Of particular focus is the case where the transport phenomena in a process can be characterized by generic laws but the process-specific knowledge of the chemical reaction kinetics is incomplete. This is one possible scenario in the early stages of process development. The aim is thus to establish a generic numerical approach to address the problem of operating regime determination. In contrast to relying on the use of different types of equipment, the proposed approach utilizes the well-known overall structure of the mathematical model of the chemical process (usually derived from universal conservation laws) with a particular formulation of the model, corresponding to each operating regime. In the process of developing such quantitative alternatives, an approach based on model discrimination has been previously proposed by the authors,8 which incorporates hybrid modeling and optimal experimental design. In this approach, a black-box sub-model was used to represent the missing chemical kinetics knowledge with the other part of the model being formulated from the mechanistic laws of conservation and transport phenomena. Adopting this approach enabled the issue of operating regime determination to be presented in terms of model discrimination. However, a black-box sub-model has the potential to approximate arbitrary relationships, thereby artificially reducing the detectable difference between two operating regimes. In this paper, an alternative approach is proposed with the aim of better utilizing the “intrinsic” distance between the models for different operating regimes. The distance between two models, as described here, refers to the difference between the inferential results for the two models. The basic idea is to only consider the known mechanistic parts of these models to avoid an artificial reduction in the distances between the different models which may occur if the hybrid models (representing different regimes) are fitted to the same set of measurement data. By not using a black-box model to capture the chemical reaction kinetics, a complete model of the process will not be available, due to the missing mechanistic knowledge about the kinetics. Consequently the model discrimination approach as reported previously by the authors8 is not be applicable. The approach proposed in this paper is based on a dependency test as opposed to model discrimination. The basic hypothesis is as follows. In principle, the intrinsic chemical reaction kinetic relationship (as opposed to the value of the reaction rate) is not influenced by agitation or other conditions that are used to control the transport phenomena in the reactor. That is, the former is independent of the latter. Thus, ideally the mechanistic (though possibly incomplete) model corresponding to the true operating regime (termed the “true model” hereafter) should reflect this independency unlike other models (termed the “false model” hereafter). More specifically, given a set of noisy experimental data, a specific dependency measure calculated according to the true model should be the lowest of those calculated, that is, when compared with “false” models. Thus, the true operating regime may be identified by ranking the computed dependency measures.

In section 2, a detailed description of the proposed approach and the supporting techniques are introduced, with an application of the methodology to a simulation of toluene nitration being presented in sections 3 and 4. Section 5 draws together the main conclusions of this study. 2. Dependency Test Based Approach In this section, the proposed approach is first described by explaining how false models, as defined in the introduction, will lead to the false dependency of the intrinsic chemical reaction kinetics on the transport-influencing conditions. The key underpinning techniques for computing a dependency measure then follow. 2.1. Origin of Dependency. In general, a chemical process system can be represented by

y(t) ) h(x(t), u(t))

(1)

where y, x, and u are vectors of the output, state, and input variables, respectively. For the purpose of illustration, the following explanation will focus on a stirred chemical reactor, where the overall rate of conversion is influenced by the agitation. More specifically, two variables, the intrinsic reaction rate (r0) and the stirrer speed (N), are considered. With reference to eq 1, r0 ∈ x and N ∈ u. As part of the chemical kinetics, r0 is determined from a function of the state variables through a relationship that is assumed in this study to be unknown:

r0 ) g(xr)

(2)

where xr are the state variables not including the intrinsic reaction rate, r0, and from which the value of r0 can be determined. Depending on the nature of the underlying chemical kinetics, r0 may rely on quantities such as the reaction temperature and the concentration of the chemical components. In this study, it is assumed that the elements in xr can either be measured or be computed from measurements. A “true” model, M0, that is, the model that represents the true operating regime of the chemical process, is now defined:

M0: y ) h0(r0, xr, N, u0)

(3)

where u0 represents the input variables other than the stirrer speed, N, and the other variables are as described previously. With the appropriate measurement information and when h0 satisfies the continuity requirements as stated in the implicit function theorem,9 r0 may be computed as follows:

r00 ) h0-1(y, xr, N, u0)

(4)

where h0-1 denotes a relationship (i.e., an explicit expression) or an algorithm (i.e., an iterative numerical procedure) derived from M0 and r00 is the intrinsic reaction rate (i.e., r0) computed through h0-1 for given outputs, inputs, and state variables. Practically, the computation of r00 can be complicated by the presence of noise in the measurements. This issue is considered further in section 2.2. Although eq 4 does not use the fundamental relationship between r0 and xr (eq 2), the resulting values of r0 should potentially satisfy this relationship, because M0 reflects the true operating regime and therefore should describe the chemical kinetics. This therefore ensures that

r00 ) g(xr) whatever form g(‚) takes.

(5)

3392

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007

Now, an incorrect model which represents an operating regime that differs from the true one is considered. It is denoted by

M1: y ) h1(r0, xr, N, u0)

(6)

between a variable (i.e., N) and a relationship (i.e., (ri0, xr)). For this purpose, it is proposed to adopt an information-theoretic approach and use the mutual information19 between N and (ri0, xr) as the measure of dependency. The general definition of mutual information between the m and n dimensional random variables X and Y is given by20

where h1 differs from h0. As for M0, the true operating regime, the values of r0 can be calculated according to this model through

r10 ) h1-1(y, xr, N, u0)

(7)

where h1-1 and r10 are the counterparts of h0-1 and r00, respectively, in eq 4. When the same measurement information is used, the difference between the two models will materialize in a difference between r00 and r10. By combining eqs 4, 5, and 7, the following relationship is obtained: -1

-1

r1 ) g(x ) + [h1 (y, x , N, u ) - h0 (y, x , N, u )] 0

r

r

0

r

0

Di: N f (ri , x ), i )1, 2, ..., n r

(9)

Because it is not known a priori which model is the right one, the value of the dependency, Di, should be determined for n candidate models to identify the true model and, hence, the true operating regime. 2.2. Realization of the Dependency Test. Given the measurement data, the dependency, Di, for a particular model Mi can be determined as follows. The first step is to obtain the pair (ri0, xr). It is assumed that xr is directly or indirectly observable and is therefore model independent. Hence, the key task is to compute ri0, which is a standard state estimation problem. State estimation is an area that has been widely studied in engineering.10,11 Of particular interest is the estimation of the chemical reaction rate from concentration measurements, which is an important task in chemical engineering. It involves the calculation of the derivatives of the concentrations and is well-known as an ill-posed problem because small errors in the concentration measurements can be significantly amplified in the computed values of the reaction rate. Existing methods for dealing with this problem are the use of smoothing functions12 and, in particular, regularized smoothing splines.13-15 Once ri0 is computed by using one of the state estimation techniques, the dependency test can be carried out, as described in the next section (section 2.2.1). Section 2.2.2 summarizes the stages of the entire procedure for realizing the proposed dependency test based approach. 2.2.1. Dependency Test. Detecting the dependencies between a set of variables is a common task in many areas, including engineering. When linear relationships are assumed, statistical methods such as correlation analysis and analysis of variance (ANOVA) are common techniques.16 For more general situations that involve nonlinear relationships, information theoretic approaches have been proposed (e.g., Ihler et al.17 and Dionisio et al.18). In this study, the dependency to be investigated is

r(x, y)

∫Rm+n r(x, y) log p(x) q(y) dx dy

(10a)

where r(x, y), (x, y) ∈ Rm+n, p(x), x ∈ Rm, and q(y), y ∈ Rn, are the probability density functions of (X, Y), X, and Y, respectively. On the basis of this definition, the measure of dependency considered in this paper can be computed as follows:

Di )

∫N∫r ∫x ...∫x p(N, ri0, xr) 0

i

r

m

r

1

p(N, ri0, xr)

log

(8)

From the second term on the right-hand side of eq 8, calculating r0 according to the false model, unlike when the true model is utilized, will not satisfy the underlying relationship between r0 and xr. A further consequence is that this relationship, as opposed to satisfying the fundamental chemical kinetics and thus being independent of the stirrer speed (N), will now generally be dependent on N. The value of r0 computed through a model Mi is generally denoted by ri0, and the dependency is given by 0

I(X, Y) )

0

r

p(N) p(ri , x )

dx1r ... dxmr dri0 dN (10b)

where p(‚) denotes a probability density function and x1r, ..., xmr are the elements of xr. In this formulation, it is assumed that all the variables are from continuous distributions. However, if the stirrer speed, N, is treated as being from a discrete distribution, eq 10b can be redefined as

Di )

∑l ∫r ∫x ...∫x pl(ri0, xr) 0

i

r

m

r

1

pl(ri0, xr) log

0

r

dx1r ... dxmr dri0 (11)

P(Nl) pl(ri , x ) where P(Nl) is the probability that N takes one of its possible values, Nl, and pl(ri0, xr) is the conditional distribution function of (ri0, xr) when N takes the value Nl:

pl(ri0, xr) ) pl(ri0, xr|N ) Nl)

(12)

Applying eq 10 or 11 requires the estimation of three probability density functions. A number of methods are available in the literature for this task (e.g., ref 21). In this work the Gaussian mixture model and the EM (expectation maximization) algorithm are adopted.22,23 With regard to the computed dependency measure, Di, the mutual information between N and (ri0, xr) should theoretically be zero when they are completely independent. However, the computational inaccuracies in the state estimation and density estimation steps will not allow such a result to be realized. Therefore, a reasonable expectation, which is adopted in the case study, would be to identify the true model as being that which yields the smallest value of the dependency measure among all the candidate models. If, among the computed dependency measures, there is no single dependency measure that is statistically smaller than the rest, it would have to be concluded that the true model does not exhibit a significant distance from the other candidate models (within the range of operating conditions considered), and therefore the true operating regime cannot be uniquely determined. 2.2.2. Stages of the Dependency Test for the Determination of the Operating Regime. A procedure for performing the proposed approach is summarized in Figure 1. With regard to how the experimental data should be generated, step 2, a number of different settings that affect the transport phenomenon require to be identified for the operating condition. Their influence on

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3393

according to the intrinsic kinetic relationship can either be measured or be derived from the measurements so that the influence of the process condition affecting the transport phenomena on the intrinsic chemical kinetics can be assessed through the analysis of the results from the dependency test. 3. Case Study: The Toluene Nitration Process

Figure 1. Procedure for dependency test based determination of operating regimes.

the computed relationship between the selected dependent kinetics state and the independent variables that affect the state can then be effectively determined. In doing so, different and representative values for this operating condition should be applied to cover its possible range. Furthermore, the number of experiments should be aligned with the number of values identified for the particular operating condition. One exemplary experimental setting following these guidelines is given in the case study (section 3). In step 6, shown in Figure 1, it is implied that the significance of the comparison between the dependency measures attained from different models should be examined. This is discussed further in section 4.1. Note that although stirrer speed (N) and intrinsic reaction rate (r0) have been used as examples for describing the approach, the applicability of this procedure is not limited to these specific quantities. The process conditions to vary in an experiment (cf. step 2) can be N or any other variables that significantly affects the transport phenomena. Similarly, the chemical reaction kinetics states involved in step 3 can be r0 or other quantities, for example, apparent reaction rate constant as in the following case study. In general, the principle for applying the proposed method is that (1) existing (mechanistic) knowledge should allow for a “partial” model to be developed which in turn allows for a chemical reaction kinetics state to be computed from the measurement data and (2) the variables that determine this state

Toluene nitration is a member of the family of aromatic nitrations in mixed acid and is an important unit process in the chemicals industry.24 The behavior of the batch and semi-batch realizations has resulted in significant research being undertaken in chemical reaction engineering, as a consequence of the interplay between the chemical reaction and the mass transfer as well as the changing (apparent) reaction rate constant during the nitration process.5,25-27 The toluene nitration process which forms the basis of the study is assumed to be run in a 2-L batch stirred-tank laboratory reactor. The first stage in the process is the charging of a particular amount of mixed acid, comprising HNO3, H2SO4, and H2O into the empty reactor prior to rapidly feeding a certain amount of toluene into the reactor. No further additions are made during the process. The reactor contains an organic and an aqueous phase, with a constant temperature being maintained. The nitration product includes three mononitrotoleuene isomers, namely o-, m-, and p-MNT. Unlike in the application of the qualitative experimental approach for the determination of the operating regimes (as reported for example in Atherton3 and in Bourne5), a constant interfacial area is not required in this study. In the remainder of this section, the selected operating regimes and the underlying mathematical models are described (section 3.1), followed by how the dependency test procedure is implemented in this case study (section 3.2). 3.1. Definition and Modeling of the Operating Regimes. Toluene nitration is a typical liquid-liquid heterogeneous reaction process. The possible operating regimes for this type of process originate from the comparative speed of the chemical reaction and the inter-phase mass transfer and have been wellstudied in the literature with respect to chemical reactor modeling and design (e.g., Doraiswamy and Sharma;2 Westerterp et al.;6 Bourne5). In this study, three possible operating regimes of the toluene nitration process were considered, namely a Very slow reaction, a slow reaction, and a fast reaction. These three regimes were also those addressed by Atherton4 in his general consideration of liquid-liquid heterogeneous reactors. The amount of H2SO4 charged into the reactor is used as the manipulating condition to control which regime the process operates within, because the amount of H2SO4 can effectively change the chemical reaction rate. From a modeling perspective, each of the operating regimes corresponds to a particular mathematical model of the overall conversion rate, r, defined as the rate of consumption of toluene within a unit volume in the total reaction medium. This leads to the conservation of toluene in the reactor being denoted by

dntoluene ) -rVreaction dt

(13)

where ntoluene is the mole number of toluene in the reactor and Vreaction is the volume of the entire reaction medium. Equation 13 applies to all operating regimes, but they can be differentiated between with respect to how the overall conservation rate is calculated. The overall conversation rate models are derived from film theory (Figure 2). In all cases, nitration is assumed to occur in the aqueous phase only, and the mass transfer

3394

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007

Figure 2. Film theory applied to toluene nitration. Figure 3. Simulation results of toluene nitration process.

resistance in the (nonreacting) organic phase is ignored. Furthermore, for any experimental run the organic phase is always dispersed into the continuous aqueous phase. Regarding the intrinsic reaction kinetics, the models assume a second-order reaction but with a reaction rate constant which, as explained in more detail in section 3.2, is dependent on the reaction conditions according to an unknown relationship. (1) Very Slow Reaction Regime.5 The nitration process is completely controlled by the reaction with the resistance for toluene to transfer from the phase interface to the bulk of the aqueous phase being neglected. This essentially yields

where kL is the mass transfer coefficient of toluene and a is the specific surface area available for inter-phase mass transfer. Both kL and a are agitation dependent and can be calculated through the models as described by Zaldivar et al.25 (3) Fast Reaction Regime.5 The nitration process is controlled by both the reaction and the mass transfer, but the reaction is so fast that it is completed within the film on the aqueous phase side, leaving the bulk of the aqueous phase with no reaction. In this case, the overall conversion rate can be expressed as

Figure 3 presents illustratively the results of the simulation for each operating regime according to the above models. 3.2. Implementation of the Dependency Test Procedure. The procedure presented in section 2 is customized for the application under consideration. First, it is not appropriate in this case study to select the intrinsic reaction rate as the dependent reaction kinetics variable. This is because the independent process variables that influence the reaction rate include the concentration of toluene in the aqueous phase. This concentration value is only available for the case where the process is operated in the very slow reaction regime (eq 14). Consequently, it is not always possible to properly investigate the relationship between the intrinsic reaction rate and all its influential factors, and hence it is not possible to conduct the dependency tests proposed earlier. For this reason, the apparent reaction rate constant, k, instead of the intrinsic reaction rate is selected as the dependent reaction kinetics variable. Other modifications made to the procedure are summarized below. (1) Factors that Influence k. When the reaction temperature is kept constant, the state variables that determine the value of the apparent reaction rate constant k are the composition of the mixed acid, that is, the concentrations of HNO3, H2SO4, and H2O.25,26 These concentrations, after being normalized, show almost identical absolute profiles. This is a consequence of the stoichiometric relationship between HNO3 and H2O and the linear relationship between the component molar numbers and the volume of the mixed acid under a constant temperature that is used in the simulation model to generate the measurement data. Consequently, only the concentration of HNO3 (arbitrarily chosen from the three components in the mixed acid) is taken to form xr, that is, xr is a scalar in this case study. (2) Estimation of k. For a particular operating regime, the value of k can be calculated according to the kinetic expression corresponding to this regime (eqs 15-17). All the quantities required for the calculation are obtained by numerical simulation using gPROMS.28 The correlations and parameters required for the simulation, that is, those for computing a, the specific surface area, and for modeling the mass transfer are taken from the literature.5,25,26 Of specific note is that r (the overall conversion rate) has to be derived from the concentration of toluene, which itself is obtained from the simulation. 5% random (Gaussian) noise was added to the simulation results to reflect the measurements that would be recorded from actual physical experiments (further discussions regarding noise levels are given in section 4.2). A smoothing function-based approach was then applied to deal with the noisy data. As the shape of the trajectory for the amount of toluene suggests that a function of exponential form would be appropriate (cf. Figure 3), the following smoothing function was adopted and fitted to the noisy data:

r ) axDtoluenekCHNO3Ctoluene,interface

y(t) ) n0ef(t)

Ctoluene,aqueous ≈ Ctoluene,interface ) mtolueneCtoluene,organic

(14)

where Ctoluene,aqueous, Ctoluene,interface, and Ctoluene,organic are the molar concentrations of toluene in the bulk of the aqueous phase, at the interface, and in the organic phase, respectively, and mtoluene is the distribution coefficient of toluene. The value of mtoluene can be found in Zaldivar et al.25 The overall conversion reaction rate, r, can then be expressed as

r ) kφCHNO3Ctoluene,interface

(15)

where k is the apparent reaction rate constant, φ is the volume ratio of the aqueous phase (i.e., the reacting phase), and CHNO3 is the concentration of HNO3. (2) Slow Reaction Regime.5 The nitration process is controlled by both the reaction and the mass transfer. The resistance for toluene to transfer from the phase interface to the bulk of the aqueous phase is considered, and it is further assumed that nitration takes place completely in the bulk of the aqueous phase with no reaction occurring within the film on the aqueous phase side. In this case, the derivation of the overall conversion rate gives

r)

Ctoluene,interface 1 1 + kLa kφCHNO3

(16)

(17)

(18)

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3395 Table 1. (a) Typical Experimental Information and (b) Experimental Information Specific to the Different Operating Regimes

Table 2. Dependency Measures Computed for Three Different Models

manipulated variables

measured variables

reaction temperature initial amount of toluene stirrer speed

remaining amount of toluene in the reactor molar concentration of HNO3

unchanged

changed in different experiments

318 K 2.0 mol 5, 10, 15 s-1 each applied once for each operating regime

true regime

D1

D2

D3

hypothesis

P value

1

5.98 ( 0.55

8.11( 0.38

11.29 ( 0.34

H 0 : D1 ) D 2

0.0007

H1: D1 < D2 2

9.12 ( 0.30

7.70 ( 0.21

9.33 ( 0.25

H 0 : D2 ) D 1

0.0001

H1: D2 < D1 3

12.25 ( 0.45

infeasible state estimation

b. Experimental Information Specific to the Different Operating Regimes mass of sulfuric acid (kg) measured length of period of reaction (s) number of sampling points (evenly distributed)

hypothesis test (R ) 0.01)

dependency measure

a. Typical Experimental Information

regime 1

regime 2

regime 3

0.6 20 000

0.7 20 000

1.9 2400

100

100

40

where n0 is the (known) initial amount of toluene and f(t) is a polynomial of t (time) with order up to four. (3) Calculation of the Dependency Measure. In this study, the stirrer speed is not expected to change continuously in any of the experimental runs. Thus, eq 11 is used for computing the dependency measure. In particular, it is assumed that for every operating regime three speeds (5/10/15 s-1) are applied, each for 1/3 of the experiments conducted within this regime. The estimation of the density functions is performed using subroutines from BNT, a MATLAB tool developed by Murphy.29 The rest of the computations, primarily the integration of the probability density functions over a two-dimensional domain, was again undertaken using MATLAB. The assumed experimental scheme is summarized in Tables 1a,b. 4. Results and Discussion Three sets of tests were performed. In each set, one of the three operating regimes was taken as the true regime, and three experimental runs were conducted, each with a distinct stirrer speed as specified in Table 1a. To present the results, the operating regimes are sequentially numbered as 1, very slow reaction; 2, slow reaction; and 3, fast reaction. The dependency measures (D) computed according to these models are numbered in the same way, which gives D1, D2, and D3, respectively. 4.1. Treatments for Reducing Inherent Uncertainties and the Results. From the results of some preliminary runs undertaken prior to the major computational work, two sources of uncertainty in the computational procedure (Figure 1) were identified. First when executing step 3, the reaction rate has to be recovered from the noisy concentration measurements. The underlying issue is that different (simulated) runs of an experiment for a particular set of operating conditions will yield different measurements, due to the presence of noise. Consequently, different sets of smoothed concentration data will be associated with the same set of experimental conditions. Consequently, the resulting mutual information and thus the dependency estimates corresponding to a particular set of operating regime assumptions may differ between different runs. To address this uncertainty, it was decided to generate multiple sets of (simulated) noise-affected measurements for a particular set of experimental specification. In the case study, five sets of measurements were actually generated. The dependency measure was then calculated for each series of measurements using the

8.47 ( 0.23

H 0 : D3 ) D 1

0.0000

H1: D3 < D1

proposed computational procedure (Figure 1). Both the sample mean (i.e., averaged value) of the computed dependency measures and the estimated standard error of the sample mean were computed and were used to draw conclusions on the outcome of the study. More specifically, a hypothesis-test framework was utilized to investigate whether the averaged dependency measure computed from one model is different from that from another model, to gain statistical evidence that the averaged dependency measure computed from the true model is smaller than that from the false model that appears to be closest to the true model in terms of the computed dependency measure. Adopting this formal approach, it is assumed that the averaged dependency measure follows a normal distribution with the computed sample mean as its mean and the (squared) estimated standard error as its variance. The difference between the averaged dependency measures from the two models thus follows a normal distribution with the difference between the two corresponding sample means as its mean and the sum of the variances of the two averaged dependency measures as its variance. With this setup, the null hypothesis (H0) that the two averaged dependency measures do not differ can be readily tested. A level of significance of 0.01 was utilized for the test. Another source of uncertainty relates to step 4 (Figure 1) where the key computational task is the estimation of the probability density distributions. The method adopted essentially estimates the parameters and weights of individual Gaussian mixture models by maximizing the likelihood. The EM algorithm implemented for this purpose is sensitive to the initial guess of the parameters and weights. Thus, different initial guesses presented to this algorithm can result in different estimates of the probability density distribution and, hence, different estimates of the dependency measures. To address this uncertainty, 100 sets of randomly generated initial guesses were considered, and it was observed that a portion of these initial guesses always gave rise to the highest likelihood value when benchmarked against other local maxima. The set of parameters and weights corresponding to this highest likelihood value were subsequently used for the computation of the dependency measure. The results of the dependency test generated following the above computational treatments are summarized in Table 2. Applying the hypothesis test framework described previously to the averaged values of the computed dependency measures, it can be observed that the proposed dependency test based approach has resulted in the correct identification of the true operating regimes with appropriate statistical confidence. This confidence was gained from the fact that all the computed P values are less than the preset level of significance (i.e., R ) 0.01). With regard to the estimation of ki, the computed apparent reaction rate constant corresponding to the ith assumed operating

3396

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007

Figure 4. Measurements with 5% noise and data derived from the measurements. (a) Noisy measurements: original and smoothed. (b) Apparent rate constant: true and derived values.

regime, no feasible result was obtained when the slow reaction model was used and the true regime was regime 3, that is, the fast reaction regime. This means that, given the measurements from a batch operated within the fast reaction regime, the solution of the slow reaction model was unable to produce a physically reasonable estimate of the reaction rate constant, ki. In this case, such an outcome allows the immediate exclusion of the possibility that the slow reaction regime is the true regime. In practice, the identification of a very slow reaction regime is fairly easy because the behavior of a process operating in this regime does not vary with stirrer speed. This piece of qualitative knowledge was not used in this case study, because the goal here was to demonstrate the potential of the proposed quantitative approach. 4.2. Influence of Noise on the Measurements. The results presented in Table 2 were obtained from simulated measurements to which 5% (in terms of the standard deviation of the

relative errors) random (Gaussian) noise was added. It can be observed from Figure 4 that although the smoothed measurement curve is similar to the noise-free curve, a significant deviation is evident between the apparent reaction rate constants derived from the smoothed measurements and their true values. While the distinction between the operating regimes has been made successfully even with the presence of such a deviation, one should bear in mind that this deviation typically increases as the level of noise increases, thereby possibly materializing in unusable derivations of the apparent reaction rates. To illustrate this, Figure 5 shows as an example the case where 10% random noise has been added to the simulated measurements. In this case the difference between the derived apparent reaction rate constants and their true values is significantly larger than for the case of 5% measurement noise. Therefore where the magnitude of noise is high, a method more advanced than that adopted in this paper, for example, that of smoothing splines,6

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3397

Figure 5. Measurements with 10% noise and data derived from the measurements. (a) Noisy measurements: original and smoothed. (b) Apparent rate constant: true and derived values.

will be required, to ensure that the proposed computational procedure generates useful results. 4.3. Comparison with Other Methodologies. The proposed approach can be compared with other methodologies that have been used to investigate the same type of reaction systems and regimes. According to the literature,4,8 the distinction between the slow-reaction (2) and the fast-reaction (3) regimes can be made through a qualitative approach, facilitated by a constant interfacial area cell to maintain the interface area between the two phases. In this work, the same task of operating regime identification has been successfully carried out in a quantitative manner through the combined use of experimental data, mechanistic modeling, and statistical analysis. This latter quantitative approach does not require the use of the constant interfacial area cell to distinguish between the slow and the fast reaction regimes. For the operating regime identification problem addressed in this paper, the results attained through the implementation of the dependency test based approach are comparable to those achieved when an alternative quantitative approach previously

reported by the authors is utilized. This method was based on model discrimination.24 The key difference between the two approaches is that the latter is based on a hybrid model for the characterization of a particular operating regime and the unknown chemical kinetics are modeled by means of a blackbox sub-model, while the approach described in this paper uses only a mechanistic (though incomplete) model. The adoption of a black-box sub-model may potentially impair the possibility of discriminating between two hybrid models, depending on the role the black-box sub-model plays in the entire model as well as how flexible the black-box sub-model is in terms of approximating arbitrary relations. Therefore, the dependency test based approach should be regarded as more reliable, although its advantage in comparison with the model discrimination based approach is still to be demonstrated through actual cases. In the paper of Yang et al.8 which describes the model discrimination based approach, a fourth operating regime, slowfast transition regime, was considered in addition to the three regimes investigated in this paper. For this regime, the nitration process is controlled by both the reaction and the mass transfer,

3398

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007

Figure 6. Simulation results for two complete, mechanistic models at conditions N ) 5, mass of H2SO4 ) 1.9 kg.

but the reaction is sufficiently fast for some toluene to be consumed within the film on the aqueous phase side, while the rest of reaction still takes place in the bulk of the aqueous phase. It was shown in Yang et al.8 that difficulties in terms of differentiation tend to occur when this regime was distinguished from its neighbors (i.e., the slow- and fast-reaction regimes). Although not reported in detail in this paper, similar difficulties were encountered when the dependency test based approach was applied to differentiate between this transition regime and its neighbors. In fact, it was observed from numerical experiments that such difficulties arose primarily as a result of the similarity of this transition regime and its neighbors under the operating conditions applied. For example, Figure 6 illustrates how similar the simulation results are for the transition and fastreaction regimes as given by the rigorous mechanistic models. This outcome implies that if a formal test was performed to distinguish between these two models (cf. Asprey and Macchietto),30 a negative result would potentially be obtained. From a physical perspective, this is understandable. When the operating conditions give rise to a very fast reaction, operating the reactor in the slow-to-fast transient regime will result in most of the toluene being consumed within the film and therefore only a very limited amount of the reaction will take place in the bulk of the mixed acid. This behavior is therefore reflective of the fast reaction regime, where the reaction is completed within the film and no reaction occurs in the bulk liquid. The implications resulting from the difficulties faced by trying to distinguish between this transition regime and its neighbors when examined in the context of scale-up studies is not serious. In practice in operating regime identification studies, if this transition regime is the true regime of the operation, then it may be identified as one of its neighbors, that is, the slow reaction or the fast reaction regime, depending on which of these two regimes it is most closely aligned with. Such an outcome will still constitute a useful input with respect to making the right decision in scale-up. 5. Conclusions In this paper the problem of operating regime determination has been formulated within the context of dependency testing.

The rationale for adopting this approach is based on the hypothesis that an incorrect assumption about the operating regime can lead to a false dependency between (1) the intrinsic chemical reaction kinetics relationship that links, for example, the apparent rate constant to its influencing factors, and (2) a variable that affects the transport phenomena such as the stirrer speed. The detection of such a dependency requires both state estimation, which in this case utilizes a smoothing function to handle measurement noise, and the calculation of the dependency measure, the basis of which was an information-theoretic dependency measure, that is, mutual information. A significant benefit of the proposed approach is that it does not rely on any special experimental requirements such as the maintenance of a constant interfacial area as is the case for qualitative approaches. The proposed methodology and techniques were investigated on a simulation of a toluene nitration process, and it was observed that a true operating regime was always correctly identified. In comparison with another quantitative approach which is based on model discrimination, the dependency test based approach potentially is able to avoid the reduction of the intrinsic difference between two compared operating regimes which may be caused by the introduction of black-box models as required by the model discrimination based approach. The effectiveness of the proposed dependency test is affected by the distance that exists between the candidate models. For example, when the models to be compared are too simplistic (due to a lack of understanding about the process) and/or are not reliable, such dependency tests may be unreliable or unfeasible. Computationally, the performance of the state and distribution density estimation will have a significant impact on the accuracy of the numerical value of the dependency measure and hence on the conclusions drawn. Therefore, care should be taken when dealing with noisy measurements (in connection with state estimation) and when applying existing statistical software tools for estimating the probability distribution densities, as has been discussed in the paper. Finally, it has been assumed in the current study that only part of the chemical kinetics knowledge is missing, while in reality the parameters that describe the transport phenomena may also not be available. An initial exploration with respect

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3399

to model identification has been undertaken to tackle this situation31 with extensions being proposed to handle operating regime determination. Acknowledgment The authors acknowledge the financial support of the EPSRC Grant GR/R64407/01 “Vertical Integration of Product Development and Manufacturing”. A.Y. is grateful to T. Chen and X.L. Ou for useful discussions on probability density estimation. Nomenclature a ) interfacial area, m2 m-3 C ) molar concentration, mol L-1 Di ) dependency measure computed according to model i Dtoluene ) Diffusivity of toluene, m2 s-1 k ) apparent reaction rate constant, l mol-1 s-1 kL ) mass-transfer coefficient, m s-1 mtoluene ) distribution coefficient of toluene n0 ) initial amount of toluene, mol ntoluene ) amount of toluene, mol N ) stirrer speed, s-1 p ) probability density P ) probability r ) overall conversion rate, mol L-1 s-1 r0 ) intrinsic reaction rate t ) time, s u ) input variables Vreaction ) volume of reaction medium, L x ) state variables xr ) state variables not including r0 y ) output variables Greek Letters φ ) phase ratio Literature Cited (1) Froment G. F.; Bischoff, K. B. Chemical reactor analysis and design; John Wiley & Sons: New York, 1990. (2) Doraiswamy, L. K.; Sharma, M. M. Heterogeneous Reactions: Analysis, Examples, and reactor Design. Volume 2: Fluid-Fluid-Solid Reactions; John Wiley & Sons: New York, 1984. (3) Atherton, J. H. Methods for study of reaction mechanisms in liquidliquid and liquid-solid reaction systems and their relevance to the development of fine chemical processes. Trans. Inst. Chem. Eng., Part A 1993, 71, 111. (4) Atherton, J. H. Chemical Aspects of Scale-up. In Pilot Plants and Scale-up of Chemical Process II; Hoyle, W., Ed.; The Royal Society of Chemistry: Cambridge, U.K., 1999. (5) Bourne, J. R. Mixing and selectivity of chemical reactions. Org. Process Res. DeV. 2003, 7 (4), 471. (6) Westerterp K. R.; van Swaaij, W. P. M.; Beenackpers, A. A. C. M. Chemical Reactor Design and Operation; John Wiley & Sons: New York, 1990. (7) Johansen, T. A.; Foss, B. A. ORBIT - operating-regime-based modeling and identification toolkit. Control Engineering Practice 1998, 6, 1277-1286. (8) Yang, A.; Martin, E.; Montague, G.; Morris, J. Numerical approaches to the determination of operating regimes with incomplete mechanistic knowledge. In Proceedings of 16th European Symposium on Computer Aided Process Engineering and 9th International Symposium on Process

Systems Engineering, July 9-13, 2006, Garmisch-Partenkirchen, Germany; Marquardt, W., Pantelides, C., Eds.; Elsevier: Amsterdam, 2006; pp 273278. (9) Goursat, E. (translated by Hedrick, E. R.). A Course in mathematical analysis; Ginn and Company: Boston, 1904; Vol. I, p 35. (10) Alvarez, J.; Lopez, T. Robust dynamic state estimation of nonlinear plants. AIChE J. 1999, 45 (1), 107. (11) Anderson, B. D. O.; Moore J. B. Optimal Filtering; Prentice Hall: New York, 1979. (12) Hosten, L. H. A comparative study of short cut procedures for parameter estimation in differential equations. Comput. Chem. Eng. 1979, 3, 117. (13) Dunfield, L. G.; Read, J. F. Determination of reaction rates by use of cubic spline interpolation. J. Chem. Phys. 1972, 57 (5), 2178. (14) Bardow, A.; Marquardt, W. Incremental and simultaneous identification of reaction kinetics: methods and comparison. Chem. Eng. Sci. 2004, 59, 2673. (15) Yeow, Y. L.; Wickramasinghe, S. R.; Han, B.; Leong, Y.-K. A new method of processing the time-concentration data of reaction kinetics. Chem. Eng. Sci. 2003, 58, 3601. (16) Haerdle W.; Simar, L. Applied MultiVariate Statistical Analysis; Springer: Berlin, 2003. (17) Ihler, A. T.; Fisher, J. W.; Willsky A. S. Nonparametric hypothesis tests for statistical dependency. IEEE Transactions on Signal Processing 2004, 52 (8), 2234. (18) Dionisio, A.; Menezes, R.; Mendes, D. A. Mutual information: a measure of dependency for nonlinear time series. Physica A 2004, 344, 326. (19) Cover, T. M.; Thomas, J. A. Elements of Information Theory; Wiley: New York, 1991. (20) Ihara, S. Information theory for continuous systems; World Scientific: London, 1993, p38. (21) Scott, D.W. MultiVariate Density Estimation; John Wiley & Sons, Inc.: New York, 1992. (22) Dempster, A. P.; Laird, N. M.; Rubin, D. B. Maximum likelihood from incomplete data via the EM algorithm (with discussion). J. Royal Stat. Soc.: Series B 1977, 39, 1. (23) Bilmes, J. A. A gentle tutorial of the EM Algorithm and its application to parameter estimation for Gaussian mixture and hidden MarkoV models; Technical report (TR-97-021); International Computer Science Institute: Berkeley, CA, 1998. (24) Schofield, K. Aromatic Nitration; Cambridge University Press: Cambridge, 1980. (25) Zaldivar, J. M.; Alo´s, M. A.; Molga, E.; Herna´ndez, H.; Westerterp, K. R. Aromatic nitrations by mixed acid: slow liquid-liquid reaction regime. Chem. Eng. Process 1995, 34, 543. (26) Zaldivar, J. M.; Molga, E.; Alo´s, M. A.; Herna´ndez, H.; Westerterp, K. R. Aromatic nitrations by mixed acid: fast liquid-liquid reactions. Chem. Eng. Process 1996, 35, 91. (27) D’Angelo, F. A.; Brunet, L.; Cognet, P.; Cabassud, M. Modelling and constraint optimisation of an aromatic nitration in liquid-liquid medium. Chem. Eng. J 2003, 91, 75. (28) Process Systems Enterprise. gPROMS AdVanced User Guide; Process Systems Enterprise, Ltd.: 2004. (29) Murphy, K. Bayes Net Toolbox for MATLAB. http://www.ai.mit.edu/∼murphyk/Software/BNT/bnt.html (accessed online on February 18, 2005). (30) Asprey, S. P.; Macchietto, S. Statistical tools for optimal dynamic model building. Comput. Chem. Eng. 2000, 24, 1261. (31) Yang, A.; Martin, E.; Morris, J. A two-level approach to the identification of hybrid models involving unknown physical parameters. Presented at 7th Chemical Process Control Conference, Dec 1-8, 2006, Lake Louise, Alberta, Canada.

ReceiVed for reView July 26, 2006 ReVised manuscript receiVed January 22, 2007 Accepted February 16, 2007 IE0609721