704
Ind. Eng. Chem. Fundam. 1980, 2 5 , 704-712
GENERAL ARTICLES Isotopic Assessment of Fundamental Catalytic Mechanisms by Kinetic Modeling John Happei" Columbia University, New York, New York 10027
Eric Walter and Yves Lecourtler Laboratoire des Signaux et Systemes, CNRS-ESE, Plateau du Moulon, 9 1190, Gif-sur-Yvette,France
Information about reaction mechanisms in heterogeneous catalysis is often obtained by modeling of rate data. This paper discusses the advantages of using transient isotopic experimentation for obtaining fundamental information as compared with the widely used Langmuir-Hinshelwood-Hougen-Watson method. Methods are also discussed to study important structural properties of postulated models that have direct consequences on our a b i l i to estimate their parameters. These properties are known as identifiability and distinguishability. The procedure is illustrated by practical examples.
Introduction The object of this paper is to show the value of transient tracing for the study of heterogeneous catalytic reactions in order to determine the fundamental parameters that appear in models characterizing them. Complex chemical reactions involve mechansims consisting of sets of elementary steps. The velocities of these steps and the concentrations of intermediate species which enter into the mechanisms constitute fundamental parameters for understanding catalytic behavior. These parameters characterize a reaction apart from the construction of overall rate expressions, which require additional assumptions. Transient tracing furnishes direct access to these basic parameters as compared with the usual methods of experimentation employing overall reaction rate data. A step which is preliminary to the procedures to be discussed here is the selection of a plausible list of possible reaction mechanisms corresponding to a given choice of elementary steps. A systematic procedure for accomplishing this has been presented by Happel and Sellers (1982, 1983). Once a list of reaction mechanisms has been prepared, a logical next step is the construction of a set of models based on the mechanisms. From the mathematical structure of such models it is possible without experimental data to determine whether the desired parameters are obtainable for a given model. To accomplish this requires a study of two closely related problems, namely those of structural identifiability and structural distinguishability. If examination of the equations describing a given model shows that there exists only one set of values for these parameters that yields one given output behavior, the model is said to be structurally identifiable. This means that if indeed actual data are satisfactorily fitted by the
model chosen, then there is only one set of parameters that corresponds to the same fit. When the parameters prove unidentifiable, they can sometimes be combined or limits placed on their values to make them identifiable. Such restrictions should, of course, only be introduced on the basis of well-defined supplementary information. In order to establish the uniqueness of a given model once it has been shown to be identifiable, it is necessary to determine whether other possible models with different mathematical structures could also have the same behavior. Indeed, information from input-output behavior may not be sufficient to distinguish a proposed model from others. General methods for treating problems of structural identifiability and distinguishability have been presented by Walter et al. (1984a,b) and are applied to heterogeneous catalytic systems in this paper. These methods are applicable both to the study of models based on the Langmuir-Hinshelwood-Hougen-Watson methodology as well as to transient tracing. In the present treatment we will devote most attention to the study of models based on tracing after comparing both methods. In this paper we have employed a deterministic rather than a stochastic framework, so that the conclusions drawn apply to ideal experimental conditions. If under this hypothesis proposed models are not identifiable, it is clear that neither would they be under noisy conditions. The LHHW Method The development of rate equations based on early studies of Langmuir and Hinshelwood was first extended and formalized by Hougen and Watson (1947). The resulting rate expressions are often termed LHHW models, and the procedure for obtaining them has been well documented in standard texts on chemical reaction kinetics
0196-4313/86/1025-0704$01.50/0@ 1986 American Chemical Society
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986 705
(Carberry, 1976; Froment and Bischoff, 1979; Hill, 1977; Smith, 1981). These models form the basis for the majority of chemical engineering treatments of mechanisms in heterogeneous catalysis as well as for deriving rate equations for reactor design, optimization, and control. The difficulty in applying these models to the obtaining of the fundamental reaction velocities and concentrations of reaction intermediates is the need to postulate relationships governing the rates themselves. Usually, to obtain tractable rate expressions it is also necessary to assume that the overall reaction rate is governed by a single rate-controlling step, with other reaction steps being at equilibrium. The following simple example illustrates the LHHW approach for the case where a single rate-controlling step is not assumed, following essentially the notation of Froment and Bischoff (1979). Example 1. It is desired to develop a rate expression for the overall reaction
Equations 6 and 7 contain three separate groupings of these parameters that will correspond to three constants, namely
available in principle from rate data in which the partial pressures of species A and B are varied. The equilibrium constant K , provides one additional relationship among the reaction velocity constants
k+lk+2 = F4 k-1k-p Assuming knowledge of F, to F4, it is possible to compute uniquely k,, and k,, k+, =
(&+
1):
k-, = 1
F3F4
A=B (1) with the following assumed mechanism for heterogeneous catalysis "+I
A+l=Al "-1
"+2
Al=B+l u-2
where 1 represents catalyst sites for chemisorption and A1 is a chemisorbed intermediate. The reaction velocities for the two mechanistic steps a t steady state are given by u , ~ (i = 1, 2). The Langmuir hypothesis and law of mass action enable the velocities of individual mechanistic steps to be expressed as follows U-1 = k-lCA1 U + 1 = k+lPAC' u+p = k+2CA' U - 2 = k-.$BC' where the reaction velocity constants are denoted by and C' refers to surface concentrations of species. steady state = u+1 - u-1 = u+, - u-2
v
(3) k,i At (4)
To obtain C1,we use the equation ct = C' CA'
and determine-whether it is possible to solve uniquely for
+
(5) The total area C, must be obtained from independent information. It is taken as constant, and if it is only desired to obtain a rate equation, it is convenient to take it equal to unity since the units of the reaction velocity constants are arbitrary. These assumptions lead to the final expression for reaction velocity
Lip)
V = k+lk+.$AC'( l - k-1 + k+2 with
Ct (7) IZ+lPA+ k-Z.PB 1+ k-1 + k+2 K , = (k+,k+,)/(k-,k-,) is the equilibrium constant for the reaction. The factor in parentheses in eq 6 will approach zero as the reaction approaches equilibrium. For extracting mechanistic information from this LHHW model, we are interested in estimating u,,, v , ~ , C', and CA',and these are the parameters whose identifiability we wish to study. For that purpose, we shall first prove that the parameters k,, and k,, are identifiable.
C' =
so that these parameters are identifiable. In some cases such a procedure can be unsatisfactory, because it may happen that some relationships exist between the Fi's that are not taken into account. In order to avoid any risk of ambiguity, a recommended procedure is to assume that more than one set of k's is possible and examine the consequences. InJhe present case we could distinguish between k,, and k,, as follows. Instead of the four equations corresponding to the constants F,, we write k - , + k,, - k., + ff+, k+lk+2 L+lk+2 1 1 _ -k+z ff+2 k-2 = - 12-2 k+lk+2 k+lk+2 k+lk+2 k + l f f + 2 -7 k-,k-2 k-lk-z
k,, in terms of k,,. Clearly, from the second equation above k + z = k+2. Substitution_ of-this result into the third equation gives k-,/k+, = k - p / k + l . The last equation then yields k-, = k-,. Using these values for k,, and k-l in the first equation gives k + , = k+l. Finally, from either of the two last equations we find k-z = kw2. Since the parameters k,, and k,z are identifiable, it is now possible to prove that u * ~ u, * ~ C', , and CA1are identifiable too, assuming that C, can be independently determined. C1is calculated from eq 7. Equation 5 then gives CA'. Finally, the u,, are obtained from eq 3. The mechanistic parameters of the LHHW model are therefore identifiable. In the section about the identifiability of models obtained using transient tracing, we shall give a more precise terminology. It is not always possible even for simple systems to identify these parameters. Thus, Happel (1968) has discussed the following slightly more complicated mechanisms for the overall reaction A = B, given by eq 1:
706
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986
Following the same procedure for this model as in the previous example, we can also obtain just four relationships from overall rate data. However, since it is now necessary to evaluate six reaction velocity constants, k*t, (i = 1, 2, 31, the model expressed by eq 8 is not identifiable. With the more complicated models often encountered in practice the application of LHHW kinetics becomes correspondingly more involved. It is often not possible to estimate values for the reaction velocity constants, even with the additional simplifying assumption of a single rate-controlling step. In order to choose which model is most acceptable among a set of possible models, the usual procedure is to obtain parameter estimates by computer modeling (see e.g. Beck and Arnold, 1977), using statistical tests to govern the choice. If the structure of the models being tested is such that they cannot be distinguished in principle, such techniques will not enable the best model to be selected. As an example of this situation, consider the case discussed by Shub et al. (1979). Example 2. Consider the heterogeneous catalytic process which corresponds to the overall reaction A+B=C (9) Suppose that there are two models possible for the reaction, which for simplicity is considered to be unidirectional, as follows:
A + E ~ A ~ ~i
+ B 2c + 1
(10)
and
B + I ~ B I
BL+ A %
c+E
(11)
The overall reaction velocities corresponding to these mechanisms may be derived following the LHHW method. For eq 10 we obtain L’+~ =
k+gAC1
u += ~ k+&&’*‘
(12)
At steady state
v = u+1 = u+2
(13)
and we have
The constants k+i (i = 1, 2) are identifiable. However, if we go through the same derivation for eq 11, we find
V=
k+l’k+2’C$flB
k+l’pA + k+2’PB In this case the model is again identifiable, but the models corresponding to eq 14 and 15 are indistinguishable on the basis of steady-state kinetics. Shub et al. suggest that the choice between the two mechanisms could be made by a study of the reaction rate under unsteady-state conditions. In this case, of course, the remaining LHHW assumptions are still necessary.
Transient Isotope Substitution The fundamental parameters, surface concentrations of adsorbed species and velocities of individual mechanistic steps, can be quantitatively studied by the use of isotopic tracers. Substitution of a traced species will cause no change in reaction rates provided that isotopic kinetic effects are negligible. Advantage was first taken of this
fact many years ago by Enemoto and Horiuti (1953) in steady-state tracing studies. We have employed it in both steady-state and transient studies, most recently in catalytic methanation (Happel et al., 1980, 1982; Otarod et al., 1983). Step velocities appear in steady-state transfer in terms of sets of algebraic equations involving tracer marking levels. Both step velocities and surface concentrations of intermediates appear in transient tracer transfer. The marking levels of traced species change with time and are expressible in terms of linear differential equations. This is an important advantage of transient tracing because complicated nonlinear rate expressions are transformed into linear differential models. Use of transient tracing has also been widely discussed by Le Cardinal et al. (1977), Happel (1978), and Le Cardinal (1978) as well as by Walter (1975) in the general context of identifiability of the parameters. A useful procedure is to conduct tests in an open system using a so-called gradientless recirculating reactor. After the reaction system is brought to steady state, the feed stream is switched to one containing the same proportions of reactants except that one or more of them are labeled by isotopes. Thus, a step function of isotopes is introduced into the system by means of isotopic substitution. The isotope transients in feed and product species occur while the kinetics of the steady-state reaction itself remains unchanged unless isotopic kinetic effects are important. The parameters involving adsorption behavior and mass action kinetics are eliminated from consideration. Only velocities of mechanistic steps and surface concentrations of intermediates appear as parameters. Furthermore, it is possible to avoid the assumption of a rate-controlling step which leads to a proliferation of models with the LHHW procedure. Intermediates may sometimes be identified by appropriate surface spectroscopies. But the presence of a surface species does not guarantee that it is involved in the mechanism under consideration. Transient tracing models are similar to those frequently employed in the biosciences and described as compartmental. In the biosciences data are usually obtained by employing very small concentrations of radioactive isotopes. It is therefore possible to describe them with a linear model. On the other hand, in our method linearization is accomplished by tracer substitution, which allows us to use large concentrations of tracers. Although the reaction velocities of individual steps are often nonlinear, the redistribution of tracer is linear. Provided that the dependent variables involving traced species are expressed as fractions of tracer in a given total atomic species, the differential equations modeling the system are linear with constant coefficients. This is useful for establishing the conditions for identifiability and distinguishability though not needed for computer modeling. The general procedure for determining parameters begins by writing material balances to describe the transient concentrations of marked species. Such models can often be expressed in terms of linear differential equations with constant coefficients as follows: M,:
dz/dt = A(8)z + Bu
(16)
y = cz
If the differential equations (16) can be solved analytically, it is possible to determine the groupings containing the parameters of 8, where 8 is defined by 0 = [C,, U*,lT
(17)
However, except for the very simple heterogeneous reaction such as that considered in example 1, explicit so-
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986
lutions of the equations is not possible. It is thus necessary to consider the structure of the systems as discussed in the following sections. Study of Identifiability Using Transient Tracing The problem of structural identifiability has been discussed by many authors. For a detailed bibliography, see e.g. Walter (1982) or Godfrey and DiStefano (1985). Most of the papers originate from econometrics or the biosciences. The notion has been discussed in a chemical context by Park and Himmelblau (1982). These authors described criteria for models represented by algebraic equations such as the LHHW rate equations as well as models represented by differential equations. An application to models of the methanation of carbon monoxide is given by Walter et al. (1986). Among methods that are suitable for testing state-space models for identifiability including linear time-invariant models described by eq 16, the transfer function approach is especially applicable for models with a small number of inputs and outputs and a highly constrained structure, with many nonexistent exchanges between compartments. This is the situation encountered in the modeling of catalytic systems by transient tracing of terminal species. The method was originally proposed by Bellman and Astrom (1970) and takes advantage of the fact that the transfer function employed in control theory yields the same parameter groupings which would appear in the solution of eq 16 without actually requiring its explicit solution. This method is employed in the present paper. The transfer function is obtained by taking the Laplace transform of both sides of eq 16. In Laplace transform space, a t zero initial conditions, the response becomes y ( s ) = G(s,8) u(s) (18) where G(s,8)is the transfer matrix, which can be computed as G(s,8) = C(s1 - A)-'B (19) where I is the identity matrix. The elements of G(s,B)are ratios of polynomials in s. The unknown parameter vector 8 is structurally globally identifiable (sgi), if for almost any value of 8 belonging to 0
c FA
Figure 1. Simple isomerization model-example 3 (identifiability).
pure component A and that all the forward and reverse step velocities are constant and nonzero (at steady state but not at equilibrium). Note that the compartments need not be physically distinct phases. Thus, the compartments A and B coexist in the gas phase. Material balances following the notation of Happel (1974) are as follows. =
W
-($ + ,,,>, +
LIZA'
The parameter vector 8 is sli for almost any value of 8 belopging to 8; the fact that (20) and (21) hold implies that 8 = 8. In that case more than one set of parameter values may exist rather than a single unique set. The actual implementation of the transfer function approach can be best understood by following worked-out examples. First we will consider the simple system given in example 1. Example 3. For the reaction considered in example 1, with the mechanism shown by eq 2, assume that the species A and B contain a single atom that can be traced. If a step function input of traced A is introduced at time t = 0 and the marking level in A is observed, we would like to determine whether the system is sgi. The model may be sketched as shown in Figure 1. The symbol resembling an eye refers to the compartment whose specific activity can be obtained from the analysis of the effluent. It is assumed that for this example the feed is
+ FFA -u W
The assumption that the overall reaction is a t steady state allows us to write the following equation:
It is convenient to combine known constants separately, designating the parameters of interest as
o1 = u - ~ , o2 = u - ~ , 0 , = CA' with ko =
FB
-,
W
kl =
FA
-, W
k2 =
PCA PCB 7, k3 = 7
Equation 23 gives
v += ~ ko
+ 01,
u+2 = ko
FFA
- ko + k1 + 02, W
Equation 22 may be written as r
i
C(O)[sI - AB]-'B(O) = C(8)[sI - A8]-'B(8) for all s (20) implies that 8 = 8. That is, the vector 8 is sgi, if different values of 8 do not give the same input-output behavior. Structural local identifiability (sli) implies a restricted parameter space such that 110 - 811 e 6, 6 > 0 (21)
707
r
1
z+
y = [l 0 012 (24)
It is possible to determine ~ ( s )by taking the Laplace transforms of eq 24 and solving the three simultaneous linear algebraic equations for zA. Using eq 16, we obtain the solution Y(s) = C(s1 - A)-'BU(s)
(25)
In view of the particular form of the B and C matrices, eq 25 becomes
where all(-l)(s)is the entry in the first row and the first column of the matrix (SI - A)-1. Thus, the transfer function may be easily deduced from eq 24
with
708
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986
i F-4v
Figure 2. Deuterium tracing of methanation (identifiability). F = rate of atomic H or D feed; u = rate of atomic H or D transfer.
and
-0,+++ k
[
(6,
+ k o ) ( b + ko)
k3
(81
(0,
i
4v
+
83123
+ 0 2 + ko)(ko + k1 + 01) - O,(ko + 81) +
+ ko)(ki + ko + 81) k2k3
j.:
(01 + k0)(02 + ko)(ki + ko) 0 3 M 3
(29)
Identifying_G(s,O) with G(s,8) as in eq 20, we find that M ( e ) and M ( e ) are output indistinguishable, if and only if the following set of nonlinear algebraic equations is satisfied. These equations are derived by extracting the coefficients of the various polynomials in s from eq 28 and 29 assuming that the k , are known. 4 0, + k , 0, 8, + 8, + k o 8, +-= + - (30) 03 k3 $3 k3 (81 + ko)(O, + ko) - (8, + ko)(8, + k,) -(31) 0, 83 8, -01+ 81 + $2 + k , 0, 81 8, + 8, + ko + (32) k2 03 + -k ,= - k+, k3 (0, + ko)(k, + hi + 0,) + 0,h,
+
s,
+
_ I _
03
(8,+ ko)(k1 + k , + 8,) -
k,
(8, + k,)(ko + h , + 8,) + 8,kl
(8, + k o ) ( k l + k o + 8,)
4
k,
-+
(33) Note that the constant term of D(s,e) in eq 29 leads to the same equation as the constant term of N(s,O),Le., eq 31. Also, eq 33 follows because of the term of degree one in s in D(s,e). In this equation the term ( 8 , + ko)(O,+ k,)/O, does not appear because it is already present in eq 31. The test for structural identifiability is thus a problem of solving the four equations (30)-(33), which can easily be treated by hand here. In more complicated cases solution by hand can become extremely complicated and symbolic manipulation routines may be required (Walter and Lecourtier, 1982). Here, the difference between eq 32 and eq 30 leads to o1 = 8, (34) Then, eq 34 and 31 yield (35) If these two relationships are substituted into eq 33, we obtain, taking eq 30 into account 0, = 8, (36) and finally, by substitution of eq 36 into eq 35 0, =
8,
(37)
Thus this model is sgi. Models of this type are of general interest because they apply not only to isomerization reactions but to any system that involves a path of transfer of marked species without addition or removal of tracers from intermediate compartments in a chain. They are often described as catenary. It can be shown that catenary models with any number of compartments are sgi, provided that all steps are reversible, and that a few conditions, enumerated in the supplementary material, are satisfied. As noted by eq 8, even a system which would correspond to a three-compartment model is not identifiable by kinetic experiments, using LHHW modeling. More complicated situations in which one or more steps are unidirectional are also treated in the supplementary material. Different situations arise when the tracer species added moves stepwise to a reacting molecule in a series of steps. Such a system is shown in Figure 2. In this system deuterium substitution is employed to study the hydrogenation of carbon monoxide over a nickel catalyst to produce methane. Since the velocity refers to H or D transfer in this figure, the output 4u corresponds to the four atoms of H or D in a molecule of CHID,-,. Details of the experiments are given by Happel et al. (1982). This model can be shown to be sgi. Values have been reported for the parameters consisting of surface concentrations of adsorbed intermediates in terms of observed fractional markings of terminal deuteriomethanes. Modeling was actually conducted using concentrations of the individual species, rather than only fractional marking. S t u d y of Distinguishability Using Transient Tracing Given a set of identifiable models, discrimination often involves experimental design to select the model that best fits the data. Methods for accomplishing this have been discussed (see e.g. Beck and Arnold, 1977). It often happens that discrimination among a large number of proposed models is not clear. One is confronted by the problem that even if it is possible to identify the parameters precisely for each model, it may also be possible from the structure of the models themselves that more than one model may fit a given set of data equally well. Even careful experimentation would then not serve to decide which model is best. An approach to this problem is to examine the structure o f the proposed mechanisms by using principles similar to those discussed in the previous section for studying identifiability. Thus, one may seek to determine whether two given models can be distinguished from each other. The procedure is then repeated for all possible pairs of an assumed set to determine the extent that they can be distinguished. The basic procedure is set forth by Walter et al. (1984, 1985) as follows. Assume that there are two models M(0)and *(a) with different structures but still the same number of inputs and-outputs. We designate the outputs as y(t,O,u)^and Yit,O,u), respectively. The parameter vectors 0 and 0 are no longer supposed to belong to the samf: spa_ce,and we denote the admissible parameter space for 0 by 0. We wish
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986
4
V+ 2
709
with Model 1
and
4
V+ 1
Di(s,O) = D(s$)Is,=o
Model 2
* Figure 3. Simple isomerization model-example 4 (distinguishability).
io know whether it is possible to compute any parameters 8 for M , assuming that we already have available the parameters 8 for M , so that the input-output behavior of M and M would be identical. In that case there would be no possibility by computer modeling to discover the inadequacy of the models from discrepapcies between inputoutput behaviors of M and M . M will be structurally output distinguishable sd from M i f for almost any chosen value of 8 E 0 the equation 9(t,d,u) = y(t,8,u) for any positive t and piecewise continuous u (38) has-no solution for 8. It does not necessarily follow that if M is sd fr0m-M that the conyerse is true, namely that M is s i from y. For M and M to be sd, M must be sd from M and M must be sd from M . As in the case of eq 20, the use of the tranfer function approach allows elimination of time and inputs. As shown by Walter et al. (1984), structural identifiability of the models is neither necessary nor sufficient for them to be sd. Thus, it may be possible to determine the most likely model without being able to determine all the parameters comprising the model. We will illustrate the method by considering the distinguishability of two models similar to example 3. Example 4. For the reaction considered in examples 1 and 3 we wish to determine whether models are distinguishable in which the first and second steps, respectively, are unidirectional as shown in Figure 3. In both cases the tracer is introduced into compartment A while compartment B is observed. In order to test the distinguishability of the two models, it is necessary to calculate the transfer function for each of them in the same manner as in example 3. Model 1 is a special case of the model of example 3, where u - ~= 02 is taken equal to zero.
+ ko + 01) 3
ko + k l
where D(s,8) is defined by eq 29. Similarly, the transfer function associated with model 2 may be written as
with
and
D J s h = D(s,8)Ia,=O Note that G depends on and e3 and that G depends on and 8,. In order to test the distinguishability of model 1 and model 2, we must write_the equality of the two transfer functions G(s,8)= G(s,8)for all s, which is equivalent to the following set of equations:
*
(01 + ko)ko p
+
(4 + ko)(ko + k J
+ ko)
+
(82
Z=[?kZ
+ ko)(ki + ko)
Note that the equality of constant terms of D1 and D2 repeats eq 43. In order to test whether model 1 is distinguishable from model 2 , we must look for solutions of eq 43-45 with respect to and &, and in order to test whether model 2 is distinguishable from model 1,we mustJook fqr solutions of the same equations with respect to tI2 and 03. Here, it is easy to simplify the set of equations. By subtracting eq 43 from eq 44, we obtain
:"; ;J+[7]
y = [O 0 112
(82
(45)
81 k2
82
-=-
(39)
01) -
+
k2k3
k283
83123
+ ko + 01) k2k3
+ ko)(ko + k J
k2
u
ko(k1
k2°3
3
- ko(02
+
k3
(46)
By subtracting the product of
-ko+ 123
ko + kl k2
Note that because of observation of the B species y = [ 1 0 O]z in eq 24 is replaced by y = [0 0 l ] z in eq 39 and then eq 26 becomes
times eq 43 from eq 45 and simplifying, we obtain
The transfer function associated with model 1 may be written as
Except for special values of ko, k l , k2, and k3 eq 46 and 47 haye no solution either with respect to 0, or with respect to 02. Therefore, model 1 and model 2 are structurally distinguishable. The identifiability of these two models is studied in the supplementary material. For simple cases like this it is often not necessary to resort to development of the transfer function. Thus in
ko81 = (ko +
(47)
710
Ind. Eng. Chem. Fundam., Vol. 25,
No. 4, 1986
m1
0
n1‘
4)
m, m4
mi b
o
0 1
0
1
-1 0 2 2
7
0
1 1 3
0 -2
2 0 3 2 0 0
1 3
0 0 2
0
example 2, the problem of distinguishing whether the intermediate is AI or Bl could be solved by marking either A or B and observing the delay in appearance of marking in C. The procedures given in this paper for the study of identifiability and distinguishability can be applied in a systematic way to determine to what extent proposed isotopic tracer experiments could lead to the unambiguous choice of a reaction mechanism from those considered possible on the basis of an enumeration such as that advanced by Happel and Sellers (1982, 1983). T o show how this can be accomplished, the following example follows that detailed for ammonia synthesis by Happel and Sellers (their example 3 in Adu. Catal.). Example 5. For the catalytic synthesis of ammonia, all possible mechanisms were developed for the reaction N2
+ 3H,
N2 + 1 + N21
Sp:
N,1+ H,
S2:
N,H&
S4: Sg:
S,: S-:
Sa:
Sq:
+1
N2H2l 2NH1
N2 + 21 + 2NI N1 + H1 e NHI
+1 +1
+ HI “21 NHl + H2 + NH, + 1 H, + 21 + 2H1 N H & + HI + NH3 + 21 NHl
+1
N2
d
b
t
+
+
= 2NH3
consistent with the choice of elementary reaction steps given in eq 48 which are considered as plausible candidates to participate in any mechanism. As before, the symbol 1 in eq 48 refers to active surface sites on the catalyst. Every species with an 1 in it is an intermediate, and the rest are terminal species. It was shown that there are six possible mechanisms for the synthesis reaction, as listed in Table I. m2 and m4 are mechanisms proposed by Horiuti (1978) and Temkin 1973), who did extensive studies of ammonia synthesis. The steps chosen for eq 48 are a combination of all those proposed by these two scientists. SI:
4
\
Table I. Mechanisms for Ammonia Synthesis stew
(48)
We wish to determine how tracing with 15Ncould serve to discriminate among these six mechanisms. Figure 4 shows the path of atomic nitrogen transfer for the six mechanisms. It is assumed that traced N, species is injected into the system in all cases and that the rate of production of ammonia as well as concentrations of species N,, H,, and NH, is known. The transient of 15Nmarking in nitrogen, zN, is observed. Note that the models are all catenary or have modified catenary paths with bridging over some of the steps. A general treatment for such systems is given in the supplementary material available for this paper. Table I1 summarizes the results for the six models given in Figure 4.
Figure 4. Models for nitrogen tracing of ammonia synthesisexample 5. Table 11. Identifiability and Distinguishability of Ammonia Synthesis Models wrocess model M I M2 M3 M4 M5 sgi sd sd sd sd M, M2 sd sgi sd sd sd __ A, sd sd sgi sd sd M, sd sd sd sgi sd &, sd sd sd sd sgi M, sd sd sd sd sd
M6 sd sd sd
sd sd
7
Conclusions regarding identifiability are given in the principal diagonal in Table 11. For the ith case, sgi indicates that the model is structurally globally identifiable. The designation “?” for model 6 with bridged structure indicates that this model requires further study to determine whether it is identifiable. The remainder of the table gives the properties of distinguishability of the models, taken by pairs. Thus, sd in the position ij (ith row, j t h column) indicates that the “model” M, is structurally distinguishable from a “process” MI. The opposite is indicated by sd. Thus, for example, the entry sd corresponding to fi6, M , means that, if a set of parameters can be found for $f2, it will also be possible to find a set of parameters for M , . This case of indistinguishability in which the number of compartments in one model M 6 exceeds that in the other M , simply involves forcing parameters of the more general model to equal zero. From a practical viewpoint, if these parameters are numerically negligible, it would be reasonable to accept the simpler model. A different kind of indistinguishability is that in which the same number of compartments is involved such as distinguishing between Horiuti’s M , and Temkin’s M , models. If nitrogen tracing were employed, experiments in which nitrogen isotope scrambling was tested
Ind. Eng. Chem. Fundam., Vol. 25, No. 4, 1986 711
reaction component concentration of unoccupied active sites per unit mass of catalyst matrix of constants corresponding to observed terminal species concentration of total active sites per unit mass of catalyst concentration of intermediate A1 on the catalyst surface per unit mass of catalyst concentration of A and B in the gas phase, vol fraction rate of feed of A to the system rate of removal of A from the system rate of removal of B from the system constants in kinetic relationship transfer matrix identity matrix reaction velocity constants for mechanistic steps combinations of known constants catalyst surface site equilibrium constant refers to a model designation for mechanisms partial pressures of components A and B Laplace transform variable step in a mechanism structurally globally identifiable structurally locally identifiable structurally distinguishable time fractional marking of feed = zFA vector of input fractional tracer levels (u is usually a scalar since there is only a single input) overall reaction velocity per unit mass of catalyst reaction velocities of individual steps i or j per unit mass of catalyst total mass of catalyst in the system Laplace transform of y ( t ) vector of observed output tracer fractions dzfdt vector of state variables, corresponding to the atomic fraction of tracer in each species fractional marking of A, Al, and B
could be employed to distinguish between these mechanisms. Finally, of course, deuterium tracing could serve to model steps involving hydrogen adsorption, which are not modeled by the transients of nitrogen transfer. Discussion In this paper two different but complementary techniques have been advanced with the objective of obtaining more fundamental information on chemical kinetic parameters than is available by usual methods. First, it has been shown that obtaining rate data by transient isotope substitution enables one to avoid the postulation of parameters which are often difficult to justify. Only concentrations of intermediate species and velocities of mechanistic steps are sought. By use of recognized numerical methods estimates of parameters can often be obtained. Appropriate additional physical and chemical tests can serve to limit acceptable values of parameters. Second, the structure of transient tracing models can be examined mathematically, independent of numerical data, to consider the problems of global identifiability and distinguishability. Methods have been presented for determining a t the outset whether the parameters can be uniquely determined and whether it is possible to distinguish models from each other. For cases in which models are unidentifiable or indistinguishable a further step in model building would be, as shown in Walter e t al. (1986), to generate the set of local solutions consistent with given data and an appropriate model. From such local solutions it is then possible to generate sets of local solutions for corresponding indistinguishable models. Physical limitations posed by other experimental data and theoretical considerations can then be employed to narrow the choice of acceptable parameters and models. The present study constitutes a step in the elucidation of the parameters which are fundamental in models of reaction mechanisms. It should be useful to scientists and engineers concerned with the development of catalysts by enabling a more exact description to be made of the significant steps in a reaction mechanism. By establishing the individual velocities of reaction steps for a range of process variables, one should be able to assess the validity of various theories involving reaction rate and chemisorption constants. Since it will enable determination of rate-controlling steps to be made, transient tracing should furnish a logical basis for the construction of rate equations of overall reactions that employ this concept. The methodology discussed here should also find application in related fields such as enzyme and electrochemical catalysis. Acknowledgment This work was supported by the US.-France Joint Scientific Exchange Program, National Science Foundation (Grant INT-81 15613), and Centre National de la Recherche Scientifique (Grant 82 0412). John Happel also received support from the National Science Foundation (Grant CFE-81-16015). Notation A reaction component A1 adsorbed intermediate matrix consisting of constants appearing in the A material balances for each compartment B reaction component B1 adsorbed intermediate matrix determined from tracer input B
fractional marking of A in feed Greek Symbols dead volume of the system e parameter vector e admissible parameter space Other Symbols to distinguish the “model” from its “process” * implies that E belongs to Literature Cited
P
Beck, J. V.; Arnold, K. J. Parameter Estimation in Engineering and Science; Wiley: New York, 1977. Beilman, R.; Astrom. K. J. Math. Biosci. 1970, 7, 329. Carberry, J. J. Chemical and Catalytic Reaction Engineering; McGraw-Hill: New York, 1976. Enemoto, S.;Horiuti, J. J. Res. Inst. Catal., Hokkaido Univ. 1953, 2 , 87. Froment, G. F.; Bischoff, K. B. Chemical Reaction Analysis and Design; Wiley: New York, 1979. Godfrey, K. R.; DiStefano, J. J., 111 I n Zdentifiability and System Parameter Estimation; Barker, H. A,, Young, P. C., Eds.; Pergamon: Oxford, 1985; p 89. Happei, J. J. Res. Inst. Catal., Hokkaido Univ. 1988, 16. 305. Happel, J. J. Res. Inst. Catal., Hokkaido Univ. 1974, 22,206. Happel, J. Chem. Eng. Sci. 1978, 33, 1567. Happel, J.; Suzuki, I.; Kokayeff, P.; Fthenakis, V. J. Catal. 1980, 68. 59. Happel, J.; Cheh. H. Y.; Otarod, M.; Ozawa, S.; Severdia, A. V.; Yoshida, T.; Fthenakis, V. J. Catal. 1982, 72,314. Happei, J.; Sellers, P. H. Znd. Eng. Chem. Fundam. 1982, 27, 67. Happel, J.; Sellers, P. H. Adw. Catal. 1983, 32,273. Hill, C. G., Jr. An Introduction to Chemical Engineering Kinetics and Reactor Design; Wiley: New York, 1977. Horiuti, J. Ann. N . Y . Acad. Sci. 1973, 273,5. Hougen, 0.A,; Watson, K. M. Chemical Process Principles; Wiley: New York, 1947; Vol. 111. Le Cardinal, G.; Walter, E.; Bertrand, P.; Zoulaiian, A.; Geius, M. Chem. Eng. Sci. 1977, 32,733.
Ind. Eng. Chem. Fundam. 1986, 25, 712-718
712
Le Cardinal, G. Chem. f n g . Sci. 1978, 33, 1568. Lecourtier, Y.; Walter, E. If€€ Trans. Autom. Control 1981, AC-26. 800. Otarod, M.; Ozawa, S.; Yin, F.; Chew, M.; Cheh, H. Y.; Happel, J. J Cats/. 1983, 84, 156. Park. S. W.: Himmeiblau. D. M. Chem. €no. J. 1962, 2 5 , 163. Shub, F. S.;Zyskin, A. G.; Slinko, M. G.; hagovski, Yu. S . ; Temkin, M I , Kinet. Katal. 1979, 20, 334. Smith, J. M. Chemical Engineering Kinetics, 3rd ed.; McGraw-Hill. New York,
Walter, E.; Lecourtier, Y.; Raksanyi. A.; Happel. J. I n Mathematics and Computers in Biomedical Application; Eisenfeld, J., De Lisi, C., Eds.; Elsevier: Amsterdam, 1985; pp 145-180. Walter, E.; Lecourtier, Y.; Happel, J.; Kao, J.-Y. AIChE J. 1986, 32, 1360.
Received for review May 8, 1984 Accepted December 20, 1985
1QR1
Tem&:’M. I. Ann. N . Y . Acad. Sci. 1973, 213, 7 9 . Walter, E. These de 3’ Cycle, Universite de Paris Sud, 1975. Walter, E. Lect. Notes Biomath. 1982, 46. Walter, E.; Lecourtier, Y. Math. Comput. Simul. 1982. 24, 472. Walter, E.; Lecourtier, Y.; Happel. J. I€€€ Trans. Aufom. Control 1984, AC-
29,56.
Supplementary Material Available: Information concerning the structural identifiability and distinguishability Of catenary and quasi-catenary (19pages)* Ordering information is given on any current masthead page.
Dynamic Matrix Control: A Nonstochastic, Industrial Process Control Technique with Parallels in Applied Statistics Babatunde A. Ogunnaike’ Chemical Engineering Department, University of Lagos, Lagos, Nigeria
Concepts which are traditionally associated with certain aspects of statistical analysis (specifically, least-squares estimation, the Levenberg-Marquart nonlinear regression method, Ridge regression, and Bayesian sequential decisions) are shown to feature prominently in the structure of dynamic matrix control (DMC). This control technique, representative of a new class of predictive-cum-optimizing control techniques that have recently made a noticeable impact on the process control community, is however strictly nonstochastic in philosophy. I n drawing out both the obvious as well as the obscure statistical parallels, DMC is presented in a new light: as an innovative technique within whose structure resides a collection of the same principles upon which these time-proven statistical concepts are based, thus providing another vehicle for understanding its success where others have failed.
1. Introduction The application of probability theory and statistics to the treatment of stochastic control problems is not new. A notable example is the Kalman filter (Kalman, 1960; Kalman and Bucy, 1961), frequently used in state estimation problems and in refining system measurements corrupted with noise. The Box-Jenkins approach to practical feedfonvard/feedback control schemes (Box and Jenkins, 1976) stems from applications of time-series linear filter models. The theory for self-tuning regulators (cf., for example, Astrom and Wittenmark, 1973) involves a combination of on-line parameter estimation and minimum variance control. In each of these cases, the dynamical system is represented by a deterministic part in conjunction with a “noise model”. The individual problems are then solved by utilizing probability and statistics in dealing with the nondeterministic aspect of the prescribed two-part model. In certain industrial processes, the random fluctuations which endow the system with its stochastic nature are sometimes negligible compared with those variations which can be explained-to a greater or lesser extent-by a deterministic mathematical model. Ignoring such random variations leads to control laws which therefore do not ordinarily employ statistical ideas. The dynamic matrix control (DMC) technique (Cutler and Ramaker, 1979) is one of such nonstochastic techniques. Despite this fact, however, it has a structure which can be related to certain aspects of applied statistics. To be sure, some of the re-
* Formerly with Shell Development Company, P.O. Box 1380, Houston, TX 77082. 0196-43 1318611025-071 2$01.50/0
semblance did not occur by design. The purpose of this paper is to highlight those aspects of the DMC framework that have statistical parallels, in addition to demonstrating the interesting equivalence between statistical parameter estimation and control in dynamical systems. This linking of the DMC framework with equivalent principles found in time-proven statistical concepts is designed to provide another vehicle for understanding DMC’s success in solving industrial process control problems where others have failed. The general control problem is discussed first. The parameter estimation problem in linear models is next addressed in a general fashion and linked in particular with estimation of parameters for dynamical systems modeled like linear filters. We then show how a dual control problem can be posed from the same parameter estimation problem. A discussion of the dynamic matrix control (DMC) technique is then presented: first the basic structure and then the additional modifications necessary for satisfactory performance. Comparisons are made with relevant statistical concepts at the appropriate places. 2. Problem Formulation In order that we may focus on the essentials, we limit this discussion to the problem of controlling a single-input, single-output system. However, we note that the technique extends quite straightforwardly to the multivariable case, as has been discussed in the original publications (Cutler and Ramaker, 1979; Prett and Gillette, 1979; Cutler, 1983). The control problem is then summarized as follows: The system output, yk,at any arbitrary time instant h is desired to be kept as close as possible to a set target value, y*. (This could mean driving the output from a different state
0 1986 American Chemical Society
