Application of Variable Structure Models to Oscillatory and Spatially

Jun 15, 1996 - line and off-line parameter estimation for a single CSTR as well as for a train of reactors with the Belousov-Zhabotinski or polynomial...
0 downloads 0 Views 142KB Size
+

+

2590

Ind. Eng. Chem. Res. 1996, 35, 2590-2595

PROCESS DESIGN AND CONTROL Application of Variable Structure Models to Oscillatory and Spatially-Dependent Processes B. E. Dainson, M. Sheintuch, and D. R. Lewin* Department of Chemical Engineering, Technion I.I.T., Haifa 32000, Israel

In this paper, the variable structure model (VSM) approach is applied to modeling and parameter estimation in processes exhibiting simple and complex relaxation oscillations. The approach reduces the full model to two simple submodels and permanently chooses the submodel which has a better performance in the description of the input-output behavior of the process. The degree to which the approach will adequately approximate the true process dynamics increases with the system stiffness. Excellent performance of the approach is demonstrated in both online and off-line parameter estimation for a single CSTR as well as for a train of reactors with the Belousov-Zhabotinski or polynomial kinetics. 1. Introduction Recently, a method has been proposed for modeling, state estimation, and control of processes characterized by sequential phases or sequential transformation steps (Dainson et al., 1995). Using the proposed method, such processes are described by a battery of alternative submodels, each of which qualitatively represents one of the process phases. The submodels cannot themselves track transitions from one phase to another. This is handled by switching to an alternative submodel when its predictive performance is significantly better than that of the submodel presently in use. This set of submodels together with a supervisor which performs switches is referred to as a variable structure model (VSM). The main advantage of this approach is that each submodel is significantly simpler than the full model of the system and it typically incorporates fewer parameters. Banerjee et al. (1995) presented a conceptually similar approach in which multiple linear dynamic models obtained either by identification or by linearization of the full nonlinear model are proposed to explain plant behavior at different operating points. The dynamics of the overall composed system is given by the weighted sum of the dynamics of the linear submodels. The weights are interpreted as probabilities of the submodels which are evaluated on the basis of the last measurements. Simplicity of the submodels is achieved since they are linear dynamic systems. This facilitates future use of well-developed linear control methods. The VSM method used in this paper is free from the requirement of linearity in submodels’ dynamics and assumes that at each time instant a single submodel describes the process. Simplicity of the submodels is understood here as lower order and fewer parameters of the submodels as compared with the full model. It has been shown that the VSM can serve as a basis for design of state observers and feedback controllers. The application proposed in Dainson et al. (1995) was for * To whom all correspondence should be addressed. email: [email protected]. Telephone: +972-4-8292006. Fax: +972-4-230476.

the description of batch processes that incorporate several phases like diauxic growth of bacteria on two limiting substrates or sequential polymerization reactions. The number of temporal phases in such systems is small (typically two or three), and the transitions are not always sharp. Periodic relaxation oscillations offer a simple example of time-dependent processes with periodic switching between two or more phases. Relaxation oscillations are caused by interaction of a fast activator, that shows several stable branches for a range of inhibitor values, and a slowly changing inhibitor. The ratio of the characteristic time scales of the activator and the inhibitor is generally small, and the smaller it is, the sharper is the switching between the phases. Relaxation oscillations appear in a wide variety of processes like catalytic oxidation reactions (Middya et al., 1993), highly exothermic reactions, the colorful BelousovZhabotinsky reaction (Field and Noyes, 1974), and various biological processes (Murray, 1989). The main objective of this paper is to apply the VSM approach for parameter estimation to systems that exhibit simple or complex relaxation oscillations. The method will also be extended to deal with distributed parameter systems, such as a train of CSTRs with or without backmixing between the units. The latter application is especially promising, since distributed systems with oscillatory kinetics exhibit an interesting array of spatiotemporal patterns. The problem of parameter estimation from such observations, or the control of patterned systems, has been scarcely addressed in the literature. 2. Methodology It would be useful to develop a criterion that will indicate when the VSM approach can be applied to an investigated process. The approach can be thought of as model reduction for a complex dynamical model. The novelty is that for different time phases or similarly for different areas of the state space, different reduced models are constructed. Therefore, the first check should be whether or not model reduction can be performed for the system under consideration. Stiff

+

+

Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 2591

systems of differential equations are good candidates for model reduction. In such systems, fast dynamics can be approximated by jumps between quasi-steady states. Time evolution along these quasi-steady-state branches occurs on a slower time scale and can be described by a reduced system of differential equations. Essentially, this will be a zero-order approximation of the original system obtained by singular perturbation methods. If the system equations can be normalized and a small parameter can be explicitly written out, its value can serve as a measure of stiffness and suggests how effective the VSM approach could be when applied to this system. For many stiff systems of differential equations the state variables can be divided into fast and slow motions and, after appropriate time rescaling, the system can be written in the form:

dxj 

dt

) hf (xj,zj)

dzj/dt ) gj (xj,zj)

(1)

with hf and gj of order unity, and  , 1 is a small parameter which is a measure of stiffness. In such systems, the fast xj variable can be assumed to follow in quasi-steady state the slow zj variable (inhibitor) according to hf(xj,zj) ) 0. The slow dynamics are then described by

dzj/dt ) gj (xj0,zj) ≡ G h (zj) where xj0 is a solution of hf(xj,zj) ) 0. If this equation is single-valued for all zj, then the reduction is simple and, except for the initial jump in xj, all fast time scales have been eliminated from the system. Problems arise when the fast subsystem is multivalued. One can solve hf(xj,zj) ) 0 along every stable branch separately to find xjj(zj). The dynamics is then governed by

h j(zj) dzj/dt ) gj (xjj,zj) ≡ G

(2)

In the examples below there are only two such stable branches, but the methodology can be extended to any number of branches. Relaxation oscillations emerge when there exists a periodic sequence of jumps and slow motions. This is always the case when the steady state hf ) gj ) 0 belongs to an unstable branch of hf ) 0. The submodels of the VSM are then described by eq 2 which typically incorporates a smaller number of parameters and only one time scale (which scales the time coordinate). Switching between submodels occurs at the boundaries of the stable branches, i.e., the limit points of hf(xj,zj) ) 0, which are characterized by

det

( ) ∂fh(xj,zj) ∂xj

)0

(3)

Condition (3) together with hf(xj,zj) ) 0 defines all variables (xj and zj) at the switching points. On the other hand, for a system of differential equations of a more general form of

dxj/dt ) h h (xj)

(4)

stiffness can be evaluated directly as the condition number of the matrix of the linearized system:

κ(xj) )

|

λmax(A C (xj)) λmin(A C (xj))

|

(5)

where A C (xj) ) ∂h h (xj)/∂xj is the Jacobian matrix and κ(xj) is the ratio of its largest absolute value of the eigenvalues to the smallest eigenvalue. As seen from definition (5), the condition number is a function of the system’s state. It can be computed in an investigated region in the state space or along typical time trajectories. In order to obtain a scalar measure of the system’s stiffness, the condition number can be averaged over the chosen region or trajectory or its minimum over these sets can be found. The stiffer the system is, the better candidate it is for model reduction and, consequently, for application of the VSM approach. Performance of a VSM can be defined in a standard way as integral square error (ISE) between its output and the actual output of the process. Moreover, it is of interest to keep track of how often switching between submodels occurs. Thus, the number of switches between submodels on a chosen time trajectory is used as an additional performance measure of the VSM. For a periodic process, it is convenient to divide the number of switches during a long enough interval by the number of cycles exhibited by the process. If this value is more than the number of submodels used to describe the process, then the VSM is “ringing” in the transition between time phases. Naturally, extensive ringing will make further utilization of the VSM, for purposes of parameter/state estimation and regulatory control, more difficult. Extensive ringing is undesirable, therefore, and indicates that application of the VSM approach to the investigated system is not valid. If a VSM results in the number of switches per cycle equal to the number of submodels used, then the VSM applies each submodel once a cycle and it can be evaluated as perfect from this point of view. 3. Application to the Belousov-Zhabotinsky (BZ) Reaction The proposed approach will now be applied to the characterization of the Oregonator model of the Belousov-Zhabotinsky (BZ) reaction (Field and Noyes, 1974):

px (dx/dt) ) f1 ≡ x(1 - x) - y(x - q)

(6a)

py (dy/dt) ) f2 ≡ -y(x + q) + bgz

(6b)

dz/dt ) f3 ≡ 2x - bz

(6c)

where x, y, and z denote scaled concentrations of HBrO2, Br-, and Ce(IV), respectively. The model parameters are

q ) 4 × 10-4, b ) 10, g ) 0.35, py ) px/24 , 1 It will be assumed that px is a variable and its influence on the performance of the VSM can be studied. The fast subset is the surface described by f1 ) f2 ) 0 or

F)

x(1 - x)(x + q) - bgz ) 0; x-q

y)

x(1 - x) x-q

(7)

The limit points xL of eq 7 are given by eq 3 or by ∂F/∂x ) 0, i.e.,

+

+

2592 Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996

(1 - 2xL)(xL + q)(xL - q) - 2qxL(1 - xL) ) 0 (8) The slow dynamics along f1 ) f2 ) 0 is described by eq 6c after substituting x(z) from eq 7. The latter is rather cumbersome, and hence it is more convenient to describe the system in terms of x or to approximate x(z). In terms of x the simplified model can be written as follows:

dz dz dx ) ) dt dx F)0 dt 1 (1 - 2x)(x + q)(x - q) - 2qx(1 - x) dx ) bg dt (x - q)2 1 x(1 - x)(x + q) 2x g x-q

( )

which gives

2gx(x - q) - x(1 - x)(x + q) dx ; ) b(x - q) dt (1 - 2x)(x + q)(x - q) - 2qx(1 - x) for x > xL1, x < xL2 (9) where xL1 and xL2 are the largest and smallest positive roots of eq 8, respectively. If q is sufficiently small, two simple approximations to the solutions of F ) 0 can be suggested (Crowley and Field, 1984). Plotting z(x) obtained from this equation, it can be verified that the smallest (and stable) root will be on the order of q. Then the first term of eq 6c can be assumed to be constant as compared with the second term and is approximated by 2q. This will lead to the equation which will constitute the first submodel in terms of the z variable only:

Submodel 1: dz/dt ) 2q - bz

(10)

The second stable root will be on the order of unity. Then q may be negleced as compared with x, and this root is approximated by an appropriate solution of a quadratic equation obtained from F ) 0 (if for a specific value of z this quadratic equation had no solution, x ) 1 was assumed). That leads to the following second submodel:

Submodel 2: dz/dt ) 1 + x1 - 4bgz - bz

(11)

If z is measured, the submodels (10) and (11) can be simulated in parallel and the leading submodel is chosen by comparison of the errors between the submodels and the measurement. Details of this procedure together with the description of correction of trailing submodels in order to prevent their divergence far from the process can be found in Dainson et al. (1995). The integral square error (ISE), the number of switches between submodels (10) and (11) (both computed over one cycle), and px are plotted against the stiffness of the system in Figure 1. The latter was computed as an average of the local stiffness, defined by eq 5, over one period of oscillations. Evidently, for very stiff systems, with very low px values, low ISE and no unnecessary switches are obtained. Note that, above a certain stiffness, the VSM exhibits exactly two switches per cycle which is equal to the number of submodels used. The performance of the VSM deteriorates with increasing px.

Figure 1. Performance of the VSM for the Oregonator model: 1000 × ISE (solid line), switches per cycle (dashed line), 1000 × px (dotted line).

Note that the stiffness of a dynamical model indicates only that some kinds of model reduction can be applied. If a model reduction yields different reduced models in different regions of the state space or in different time phases, a VSM is obtained. If the reduced models were designed properly, the stiffness of the obtained VSM can be significantly lower than that of the full model. On the example above it is easily seen, as the submodels do not include small parameters px and py. This will result in much faster simulations and, more significantly, may facilitate future utilization of the VSM for process control, optimization, and parameter estimation. 4. On-Line Parameter Estimation The simplicity of the submodels facilitates on-line parameter estimation with the help of VSM. Since each submodel consists of fewer equations than the full model of the process with hopefully fewer unknown parameters, the VSM provides an efficient approach for parameter estimation. It is worthwhile to perform parameter correction of a submodel only during the time phase where this submodel is a leading one. According to the accepted approach, trailing submodels do not represent dynamics of the system and there is no reason to correct their parameters on the basis of current measurements. This approach has been tested on the described Oregonator model of the Belousov-Zhabotinski reaction. It is assumed that the values of the parameters b and g are unknown and should be estimated based on the measurement z(t). Here, a general procedure is used which can be applied to systems which depend nonlinearly on unknown parameters. Namely, the state vector of the submodel is augmented by the parameters to be estimated, and a nonlinear observer is used to estimate the new state on the basis of available measurements. The state vectors of the submodels are

u j 1 ) [z bˆ ]T

and

u j 2 ) [z gˆ ]T

with bˆ and gˆ denoting the current estimates of the corresponding parameters. Note that the second submodel uses the last estimate of the parameter b generated by the first submodel. Thus the original parameter estimation problem has been decoupled into two problems, each involving a single dynamic equation with a state variable that is measured and a single parameter to be estimated. This task is significantly simpler than the estimation of two or even four parameters (if px and py are not known) of the original stiff model of three equations measuring a single-state variable. The non-

+

+

Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 2593

Figure 2. On-line parameter estimation for the Oregonator model. True values of the parameters are b ) 10 and g ) 0.35.

linear observer proposed by Gauthier et al. (1992) has been used. Numerous simulations have shown that good estimates for the unknown parameters are obtained after a few periods of oscillations (see Figure 2 for a typical run). Incorrect initial estimates for the parameters do not lead to frequent erroneous switches between the submodels. 5. Off-Line Parameter Estimation A similar procedure can be applied for off-line parameter estimation. Namely, the VSM is first run with the initial guesses for parameter values in order to determine the time intervals where each submodel is the leading one. Data from these time intervals are used to identify parameters of the leading submodels in the next step. Now the determination of the time phases corresponding to different leading submodels can be corrected with the help of the latest estimates of the parameters. This way several iterations can be performed if needed. If some of the parameters enter several submodels, the most general identification procedure will search for all of the parameters using all the data available. Nevertheless, it can be worthwhile to use a stepwise approach: at each step identify the submodel with the minimal number of unknown parameters and substitute the obtained estimates into the others. Once estimated on a step of the procedure, parameters are assumed to be known for the subsequent steps. This method is motivated by the widely used practice of estimation of different groups of model parameters from separate experiments and/or separate parts of experimental data. This can ease convergence and significantly decrease computational load. If required, final tuning of parameter estimates obtained by the stepwise procedure can be performed using data corresponding to different time phases. If some parameter appears in different submodels, then it can be used to test the quality of the VSM. Namely, estimates of such a parameter found by several submodels can be compared. Close values confirm that the VSM application is successful, while a significant difference between the estimates indicates problems. Similarly, estimates of a submodel can be compared with the results of an identification procedure which was applied to the whole model. The stepwise parameter identification can be demonstrated on the Oregonator model of BZ oscillations. The first submodel (eq 10) contains a single parameter b, while both of the parameters enter the second submodel (eq 11). Thus, in the first step, b can be identified from the first phase and its estimate can be substituted into the second submodel. In the second step, g is estimated from the data corresponding to the second

submodel. This procedure results in the values b ) 9.84 and g ) 0.364. For comparison, simultaneous identification of these parameters on the basis of all the data available was performed by the standard simplex algorithm and gave similar results: b ) 9.95 and g ) 0.358. The parameter identification procedure for the full model of the investigated process can utilize information obtained by identification of the VSM. Namely, parameter estimates of the VSM can serve as values for the full model, or they are expected to be a good initial point for any identification procedure. If a parameter of the full model does not appear in any of the submodels (as when  f 0), this indicates either that the sensitivity of the full model to the parameter is low or that it mainly controls the transition between phases. For the Oregonator model of the BZ reaction, parameters px and py have been considered as the small parameters when the singular perturbation technique has been applied. The submodels obtained are zero-order approximations and, therefore, do not contain px and py. These parameters mainly influence the accuracy of the zero-order approximation and control the speed of the transition from one time phase to another, which is assumed by the VSM to be instantaneous. Therefore, it can be concluded that while these parameters remain small, the full model has a low sensitivity to their variations. As a consequence, such parameters cannot be accurately determined from the measurements available. 6. A Train of CSTRs The VSM approach will now be applied to problems of modeling and parameter estimation of a train of n CSTRs in which the BZ reaction takes place and where certain communication, due to flow or backmixing, occurs between the reactors. Such a system can also be considered as an approximation of a distributed parameter system. Specifically, there has recently been interest in pattern formation due to convection, diffusion (or backmixing), and reaction in a fixed-bed or a plugflow reactor (Sheintuch and Shvartsman, 1994; Rovinsky and Menzinger, 1994). The simplest approximation for such systems is a cascade of mixed reactors with flow down the reactor and certain backmixing between reactors. In the example below it will be assumed that x and y are immobilized and are not flowing or backmixing, while z is flowing and backmixing. Phase plane analysis demonstrates that a flow term of the fast autocatalytic components x and/or y will cause the reactor to synchronize (xi ) xj, yi ) yj, zi ) zj) since a transition from one f ) 0 branch to another in one reactor will induce a similar transition downstream. When z is flowing, a phase difference will develop between adjacent reactors. Attention will be restricted to this case which will produce wavelike patterns (Figure 3). Equations for ith CSTR will take the following form:

px (dxi/dt) ) f1i ≡ xi(1 - xi) - yi(xi - q) + qxf(xi-1 - xi) + qxb(xi+1 - xi) (12a) py (dyi/dt) ) f2i ≡ -yi(xi + q) + bgzi + qyf(yi-1 - yi) + qyb(yi+1 - yi) (12b)

+

+

2594 Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996

Figure 3. Wavelike patterns of the train of four CSTRs with the BZ reaction.

Figure 4. On-line estimates of the bi parameters in a train of different CSTRs.

dzi/dt ) f3i ≡ 2xi - bzi + qzf(zi-1 - zi) + qzb(zi+1 - zi) (12c) Here qmf and qmb denote forward and backward flows of component m, respectively. Analysis of each unit of the train is similar to the analysis of a single CSTR. Two solutions corresponding to stable branches of fli ) f2i ) 0 obtained analytically or numerically for a current value of zi should be substituted into eq 12c. This will form two submodels for each section. Therefore, n pairs of submodels are formed and, at each time instant, a leading submodel from each pair is chosen. The solutions explicitly use inlet concentrations of the fast components (xi-1, xi+1, yi-1, and yi+1) which may be unmeasured. If so, their estimates by the leading submodel of the corresponding section can be utilized. If all the parameters of the full model are known, a typical run of the VSM of the train of CSTRs looks as shown in Figure 3. The parameters are qxf ) qxb ) qyf ) qyb ) 0, qzf ) 0.05, qzb ) 0.0125, and b ) 10 for all of the reactors (with all other parameters as before). The VSM based on measurements of zi perfectly follows the full model with exactly two switches between alternative submodels for one cycle of oscillations. Trains with more CSTRs were found to exhibit the same features, but they are a bit more difficult for representation. Distributed systems may be nonuniform and have different values of kinetic parameters in different units due to differences in catalytic activity, in heat transfer coefficient, etc. Similarly, parameters in distributed systems may have spatial variations. This significantly increases the number of unknown parameters and makes the identification procedure much more complicated. As in the case of a single reactor, the developed VSM can serve as a basis for parameter estimation. To demonstrate this, it will be assumed that units of the described train have different true values of the b parameter: b1 ) 11, b2 ) 10, b3 ) 9, b4 ) 8, where the subscript denotes the number of the unit in the train. Estimates of these parameters are augmented to the state vectors of corresponding submodels as it is described in section 4, and the nonlinear observer of Gauthier et al. (1992) has been applied. Simulation results (Figure 4) show that the estimates converge during one cycle and remain relatively close to the true values of the parameters. The initial estimates were taken as bˆ i(0) ) 12 for all of the reactors. Usually, submodels will have a lower accuracy in the regions near the critical points (which correspond to switches between submodels). Therefore, parameter estimation based on the dynamics of submodels in these regions may be incorrect. This effect is clearly seen on Figure 4 as a fast periodic deviation at the switch points with

Figure 5. Simulation of VSM for the cubic-nonlinearity model (13a-c): (a) x(t); (b) yˆ (t) coincides with y(t).

rapid convergence to the true values in between. If the switching areas can be specified, this feature can be eliminated simply by ceasing parameter correction in these areas. 7. Systems with More Complex Dynamics The proposed approach is not limited to simple relaxation oscillations, as demonstrated above, and can be applied to complex (multipeak, chaotic, or quasiperiodic) behavior. The method will now be applied to describe a system with a two-dimensional slow dynamics manifold on which the motion is oscillatory. The model equations are



dx ) f1 ≡ -x3 + px2 + qx - ry dt

(13a)

dy/dt ) λx + ay + bz

(13b)

dz/dt ) cy + dz

(13c)

with the parameters a ) 1, b ) -2, c ) 2, d ) -0.8, p ) 5, q ) -6.3, r ) 0.1, λ ) 10, and  ) 10-4. Models with cubic nonlinearity have long been employed as a simple representation showing complex dynamics in mechanical (Van-der-Pol oscillations) as well as in chemical systems (Sheintuch and Luss, 1989). Analysis of this model can be performed according to section 2. Equation f1 ) 0 has two stable branches of solutions: the low (∞ < x e xL1 ) 0.84) and high (2.49 ) xL2 e x < ∞) branches. For x values on the low branch, the slow dynamics described by eqs 13b and 13c has the character of an unstable focus. When the amplitude of the oscillations increases and y reaches the limit point yL1 ) -23.57, x jumps to the high branch. At this branch y increases to the second limit point at yL2 ) -1.25, which is followed by the jump of x back to the low branch, and the trajectory is injected to the unstable focus on the low branch.

+

+

Ind. Eng. Chem. Res., Vol. 35, No. 8, 1996 2595

This process can be described by the VSM comprised of two submodels. Each of the submodels has two state variables yˆ and zˆ which represent the slow dynamics of eqs 13b and 13c. The x variable in these equations was substituted by xˆ 1 or by xˆ 2 which are the stable roots of f1(xˆ ,yˆ ) ) 0 from the low and high branches, respectively. These roots were found numerically at each step of integration. When the equation f1(xˆ ,yˆ ) ) 0 had a single root which belonged to the low (high) branch, the second root was assigned the limit value of the opposite branch: xˆ 2 ) xL2 (xˆ 1 ) xL1). Simulation of the model (13) together with the VSM which was based on the measurement of y only is given in Figure 5. The VSM follows the full model perfectly without any redundant switches between the submodels. As in the examples above, this VSM can serve for the estimation of unknown model parameters. 8. Conclusion It has been shown that the VSM approach is appropriate for the modeling and parameter estimation of systems featuring relaxation oscillations. The degree to which the approach will adequately approximate the true process dynamics increases with the system stiffness. Excellent performance of the approach has been demonstrated in both on-line and off-line parameter estimation for the Oregonator model of the BelousovZhabotinsky reaction for a single CSTR as well as for a train of reactors. These results suggest that the VSM approach should be successful in state/parameter estimation and control of distributed parameter systems featuring relaxation oscillations.

Tartakovsky from the Biotechnology Research Institute, NRC, Montreal, for his contribution to this work. Literature Cited Banerjee, A.; Arkun, Y.; Pearson, R.; Ogunnaike, B. H∞ control of nonlinear processes using multiple linear models. Proc. 3rd Eur. Control Conf. 1995, 2671. Crowley, M. F.; Field, R. J. Asymptotic solutions of a reduced Oregonator model of the Belousov-Zhabotinsky reaction. J. Phys. Chem. 1984, 88, 762. Dainson, B. E.; Tartakovsky, B.; Sheintuch, M.; Lewin, D. R. Variable structure models in process observation and control. Ind. Eng. Chem. Res. 1995, 34, 3008. Field, R. J.; Noyes, R. M. J. Chem. Phys. 1974, 60, 1877. Gauthier, J. P.; Hammouri, H.; Othman, S. A simple observer for nonlinear systems applications to bioreactors. IEEE Trans. Automat. Control 1992, AC-37, 87. Middya, U.; Graham, M. D.; Luss, D.; Sheintuch, M. Pattern selection in controlled catalytic wires. J. Chem. Phys. 1993, 96, 2823. Murray, J. D. Mathematical Biology; Springer-Verlag: Berlin, 1989. Rovinsky, A. B.; Menzinger, M. Differential flow instability in dynamical system without an unstable (activator) subsystem. Phys. Rev. Lett. 1994, 72, 2017. Sheintuch, M.; Luss, D. Identification of bifurcation centers in systems with complex dynamic behaviour. J. Phys. Chem. 1989, 93, 5227. Sheintuch, M.; Shvartsman, S. Patterns due to convectiondiffusion-reaction interaction in a fixed-bed catalytic reactor. Chem. Eng. Sci. 1994, 49, 5315.

Received for review February 27, 1996 Revised manuscript received May 8, 1996 Accepted May 10, 1996X IE960118H

Acknowledgment This work was supported by the U.S.-Israel Binational-Science-Foundation. The authors thank Dr. B.

X Abstract published in Advance ACS Abstracts, June 15, 1996.