- Email: [email protected]
Chaotic behaviour of a controlled CSTR
0098-1354190$3.00+0.00
Compuferschem.Engng,Vol.14,No. II,pp. 1237-1247, 1% Printed in Great Britain. All rights reserved
CHAOTIC
BEHAVIOUR L.
Department
of Industrial
(Received 23 Mar&
Chemistry
Copyright
OF A CONTROLLED
PELLEGRINI
Q
1990 Pergamon
Press plc
CSTR
and G. BIARDI
and Chemical Engineering I-20133, Milan, Italy
“G. Natta”,
Politecnico
di Milano,
1990; final revisionreceived4 July 1990; received for publicarion 18 July
1990)
Abstract-The influence of a PI controller on the dynamic behaviour of a continuous stirred tank reactor (CSTR) where an exothermic reaction takes placehas been investigated. The occurrence of chaotic phenomena is singled out by means of both topological features and of sensitivity analysis to initial conditions. The Feigenbaum route to chaos is recognized, since period doubling bifurcation cascades are observed: the universal characteristics of such a “scenario” are well pointed out. The possible destabilizing effect of a PI controller on such nonlinear systems is demonstrated. A parametric analysis in terms of proportional and integral control constants (kr vs k,) shows how no general guidance can be given in tuning control parameters, within the region where chaotic phenomena may occur.
INTRODUCTION shown in various Chemical reactors, typically, are among the systems which have been most studied for detecting unusual and/or pathological behaviour in chemical engineering. Since the earliest, classical, studies in the fifties by Aris and coworkers (1958) on bifurcation phenomena in chemical reacting systems, much has been written on steady state multiplicity, orbitally stable solutions, etc. both from the experimental and from the theoretical point of view. However, only in the last decade have researchers in chemical engineering paid attention to the aperiodic or chaotic behaviour of deterministic systems, that were already known in fluid-mechanics and atmospheric science (Saltzman, 1962; Lorenz, 1963; Ruelle and Takens, 1971) from studies on “turbulent motions”. An extended review regarding chaos in deterministic systems and its applications in the chemical engineering field has been made by Doherty and Ottino (1988) where a complete description of the chaotic phenomena and of the available methods to detect them is reported. In addition extensive literature relevant to the study of aperiodic phenomena from the experimental and/or theoretical point of view is cited. CSTR has certainly been the most widely studied reacting system. Kahlert ez al. (1981) and Jorgensen and Aris (1983) pointed out that consecutive reactions in a stirred tank can induce pathological behaviour leading to chaos. A cascade of perioddoubling bifurcations in an isothermal CSTR where quadratic kinetics take place has been analyzed by Smith and coworkers (1983). Recently great interest has fields in nonlinear dynamics.
been
Always in the field of stirred tank reactors Mankin and Hudson (1984) showed how chaotic phenomena can occur under forced operation. Chang and Chen (1984) have presented the possible destabilizing effects of feedback controllers, studying a substrate inhibition model and showing how exotic dynamic behaviour can be induced by adopting a simple PI controller. Prairie and Bailey (1987) have reported the occurrence of feedbackinduced Hopf bifurcation in a closed-loop CSTR, studied for the testing and development of a kinetic model for ethylene hydrogenation. In this paper where a nonisothermal, but kinetically simple CSTR is analyzed, we give further evidence that feedback controllers which are designed, according to linear theories, to stabilize a desired operating point, may, on the contrary, induce complicated dynamic phenomena which culminate in chaos.
THE SYSTEM Let’s consider at first the open-loop system. The dynamic model of a cooled CSTR, where a first-order exothermic reaction A + B takes place, is represented by the material balance for the species A and the enthalpy balance. By adopting the usual notation (see, however, the end of the paper), we can write the following autonomous system of ordinary differential equations:
Vz =Q (CA, - C,) - k,emEiRTCAV, Vpcp~=Q~cpVo--T) -AH,k,e-EtRTCAV-UA(T 1237
(1) - T,),
L.
1238 By means of the following
V 7=--,
ad
e
T-
=
ara.q
+
koTe-EiR(~~,r+Ar,dO)(]
T,,
1 _
_
form:
5 ),
(2)
0,
-k,
AT dt, s0
to the dimensionless
form:
f At?, = -kpAe
<),
ATa, ’
The dynamic features of these equations have been theoretically studied by Uppal et al. (1974). Different operating conditions may occur: one steady state, two steady states, periodic oscillation around a singular point [each being characterized by its region of asymptotic stability: see Pellegrini et al. (1988)]. However, as discussed in the basic paper by Ruelle and Takens (1971), such a system cannot exhibit a chaotic behaviour since it shows only two state variables or, in other words, two degrees of freedom. Let’s now consider the same system with a closedloop controlling outlet temperature. It is well known that a CSTR is usually “stabilized” at the desired set point (<,, T,) by means of a SISO PI controller whose integral action allows us to asymptotically remove the offset related to the proportional control. As for the feedback control: the reactor inlet temperature T, stands for the manipulated variable whilst the outlet temperature T is the controlled variable. The proportional-integral controller presents the following mode:
equivalent
_
7
k,Te -~/~(rrcr+
AT, = -kpAT
koTe-“,R(r~f+Ar~d”)(]
TT--4%
de dr = 8,, - NTUB, - (1 + NTU)B
1
+
de
the system can be put into its adimensional
+
_(
df
f=-,
Q
=(-*HR)c*o,
-<
_= I
t
PCP
TT=
and G. BIARDI
d5
assumptions:
UA cA, - cA , NTU=epe,; A0
5=.
*T
PELLEGRINI
- k,s
that now presents three degrees of freedom. As it has been proved by Packard et al. (1980), the two formulations are perfectly equivalent. It’s worthwhile stressing that the assumed set point is the only singular point of the closed-loop system. For the set of parameters listed in Table 1 the closed-loop system exhibits an aperiodic behaviour which seems, to every appearance, chaotic from a phase space analysis as shown in Fig. I. The aim of this paper is to clearly detect the onset of chaos and its occurrence in such a simple system and to investigate the possible road to chaos or “scenario”, by means both of the analysis of geometrical features and of the methods based on sensitivity to initial conditions. THEORY
Floquet theory has been applied to find out how periodic solutions degenerate towards chaos. The theory analyzes the stability of periodic orbits; a complete description is reported by Berg& et al. (1984). Since infinitesimal perturbations are considered, a linear stability analysis is sufficient. Let x(t) be the periodic solution of period T for a nonlinear autonomous system x = F(x, p) where p is the bifurcation parameter. In order to determine the orbital stability, what happens to an initial infinitesimal disp!acement Sx from the solution is analyzed. Linearizing the system along the limit cycle, we find that the initial condition x0 + 6x is mapped after one period into x + W(T) 6x where W(T) is the Floquet matrix.
A9 df, s0 Table I
where: Ae,=e,-e,,, *e=e-e,. So the original ODE system (2) is transformed into an integro-differential one or, after differentiation, into the following ODE system:
7=36OOr k,r=4.807908 E/R=SOOOK AT,,=ZOOK T,= r,<,=430 K Th= 298.42K T, = 373.16K NTU =O.S k, = 0.3 k.r=0.3
Chaotic behaviour of a controlledCSTR
380.00
400.00
420.00
440.00
460.00
1239
480.00
500.00
520.00
Temperature (K)
520.00
_
500.00 1 480.00 E r
460.00:
J P E
440.00 7
E I-
420.00400.00<
Fig. 1. Phase plane portrait and temperature vs time diagram of the chaotic behaviour for k,r = 0.3.
The problem of the linear stability for a periodic solution is reduced to the study of the eigenvalues of W(T) that are called the characteristic Floquet multipliers (CFM). For an autonomous system one of the CFM is always unity: this corresponds to a displacement 6x along the trajectory. The periodic orbit is stable if all the other Floquet multipliers are located inside the unit circle of the complex plane. When, by continuous variation of the system bifurcation parameter, one of the eigenvalues exits in the unit circle, loss of stability of the periodic solution occurs. The type of crossing ( - 1, + 1, complex conjugated eigenvalues) has different consequences on the system behaviour and is related to the “scenario” of the route to chaos. Swinney (I 983) has reviewed in his paper the most known among such routes to chaos.
A very useful tool for detecting chaotic phenomena is also provided by the Lyapunov characteristic exponents defined as: 1, = lim 1-00
log ‘ i(t) zP,@)’
where pa is the length of the principal axis of the n-ellipsoid into which an infinitesimal n-sphere of initial conditions has been transformed due to the locally-defined nature of the flow. In particular the largest Lyapunov characteristic exponent is a measure of chaos: a positive exponent shows the system sensitivity to initial conditions (two nearby trajectories diverge in the phase-space); moreover the value of the largest Lyapunov exponent is related to the rate of divergence: the greater LCE is, the faster the system leaves its initial conditions.
L.
1240
PELLEGRINI
and G.
IO'
lo2
10
1
BIARDI
Time
-4 0
20
40
60
60
100
, 120
I
I
t
I
1
I
20
40
60
60
100
120
-60 0
Time
Time
Fig. 2. Lyapunov characteristic exponents for
In our computations we followed the method suggested by Benettin et al. (1976) that has proved itself really valid and reliable though sometimes very time consuming. SIMULATlON
RESULTS
As a first step an accurate analysis was performed for the case of the operating conditions shown in Table 1 whose exotic behaviour is well illustrated by Fig. 1. The Runge-Kutta-Merson method with selfadjusting step was employed to follow the sharp peaks of 5 vs time and of T vs time. Because of the delicate nature of the dynamic problem under study, special attention was paid to the numerical computation; to this purpose all the routines used were in double precision. In our opinion the figure represents a good portrait of chaos; however, to dispel possible doubts, the three Lyapunov exponents were computed and reported in Fig. 2a-c. Figure 2a shows that, asymptotically, the largest Lyapunov characteristic exponent 1, exhibits a positive value (1, z 0.2); the other two Lyapunov exponents AZ and A, are, respectively, 0 and ~0 according to the theory. Moreover, the Poincari: section (Fig. 3), constructed following the suggestions by Henon (1982), is typical of a chaotic situation since no preferential points of accumulation can be noticed.
k,r =
0.3.
A parametric exploration of the adimensional integral control constant k,r was performed in the phase space in order to establish the chaotic range and the possible route to chaos: the simulation results are shown in Figs 4-7. Figure 4 represents a limit cycle corresponding to a period-l solution for k,r =O.l; at k,r =0.17 (Fig. 5) one still single, but of increased amplitude, limit cycle is found (on the other hand, moving towards lower k,s values, k,r +O, more and more shrinking limit cycles are encountered). As the bifurcation parameter is increased further, the solution becomes period-2, period-4, period-8 (Fig. 5) up to k,r z 0.23 where the phase space analysis seems to indicate the occurrence of chaos.
480 r
ConversIon
Fig. 3. PoincarC map for k,r =0.3.
Chaotic
behaviour
of a controlled
CSTR
1241
o.a2= 0.80;
430.00
435.00
440.00
Temperature
(K
445.00
450.00
I
425.00
420.00 Time
Fig. 4. Phase plane portrait
and temperature
Such behaviour can be noticed in the range ~0.23 + 0.43 (see Fig. 6 for k,r = 0.35) whilst computation at higher parameter values shows limit cycles of decreasing complexity and periodicity till one single limit cycle appears beyond k,r z 0.56 (Fig. 7). The picture looks like the classic one of the Feigenbaum period doubling cascade as a route to chaos. Floquet multipliers were computed solving the system (3) simultaneously with the linear sensitivity equations given by: k=Jlw
W(O)=U,
where JI stands for the Jacobian matrix evaluated along the curve x(t, p) and 0 for the identity matrix. As suggested also by Jorgensen and Aris (1983) a check was made on the unit Floquet multiplier and on the period computed between two successive
vs time diagram
for k,~ = 0.1.
temperature minima with an error, according to Nandapurkar et al. (!986), 0) we have always found three real characteristic Floquet multipliers, the first always equals unity while the third one is close to zero. The second one, related to the orbital stability, moves out from the unit cycle always along the real axis at - 1 (+ Feigenbaum “scenario”). This happens at each loss of orbital stability giving rise to another stable cycle whose period is twice the initial one. Table 2 shows the parameter values corresponding to a period doubling bifurcation and the relevant period (the value is the last one before bifurcation). Numerical problems did not allow period doublings beyond those listed in the table. About the same results have been achieved from the analysis of the k,r range from 0.4 to 0.6,
L.
I242
PELLEGRINI
and G. BIARDI
Fig. 5. Period doubling
showing the existence of a symmetrical period doubling cascade. The results of the Floquet multipliers computation for this last case are reported in Table 3. The use of AUTO package (Doedel, 1986) based on the continuation techniques as from Keller (1977) to check the validity of our computational analysis has given quite similar results. The period doubling bifurcations diagram from AUTO is reported in Fig. 8.
DISCUSSION
Paradoxical as it seems, it is now well known that some kind of order exists within chaos.
cascade.
Feigenbaum (1978) has shown that period doubling bifurcation occurs in a quite universal way and this may be tested for the case at hand. The first measurable property of the Feigenbaum route we tried to find out from our analysis is the f ratio:
* =LI-P”
” A-P”+,’ where p, is the parameter value corresponding to the doubling; its value, as n increases, approaches the universal constant 6 = 4.66920 1609 1. In Table 4 the k, parameter value and the calculated S, are listed for each period doubling: as it can be easily seen 6, approaches the Feigenbaum constant
Chaotic
0.40
I 380.00
I
I
51
I
I
400.00
I
behaviour
I [ 420.00
I
I
I
of a controlled
I
1 I
1 I
440.00
11
1 I
460.00
Temperature
CSTR
I
1243
I
1 I 9 1 I
480.00
1 v I
500.00
1 I 520.00
(K)
500.00 490.00 480.00 470.00 460.00 450.00 440.00 430.00 420.00 410.00 400.00 390.00 Time
Fig. 6. Phase plane portrait
and temperature
very well, also when compared with other authors results (see Smith et al., 1983). Another universal property relevant to the periodic regime preceding the chaotic region along the Feigenbaum route is that: K = pc, - const S -n, which allows, once two parameter p-values corresponding to period doubling are known, to evaluate the period doubling accumulation point pK. Table 5 reports the accumulation point values for the left to right cascade computed from each available pair of period doubling parameters; all the listed values don’t differ very much one from the other so that a mean value of pI = 0.22076 can be assumed without large errors. Similar analysis has been performed for the symmetrical cascade (see Table 6), so enabling us to
vs time diagram
for k,s = 0.35.
define the correct range of chaotic behaviour: k,r = 0.2207 + 0.437; such values are supported by punctual investigation by means of Lyapunov exponents as reported in the previous section. A remark that, in our opinion, may be of interest is that a window of periodicity is encountered for k,r z 0.40 + 0.4025. It is a narrow range of p within which the trajectories are strictly periodic. Just period-3 trajectory (Fig. 9) have been detected from our analysis due to the narrow range. Perhaps it would be possible to show how the universal behaviour can be extended (even if with different constants) to this window but it would be, in our opinion, just a piece of cleverness to go on along this way. Much more importance from the engineering viewpoint, should be accorded to the influence of the control action on the extent of the chaotic region.
L. F%LLEGRINI and G. BIARDI
Fig. 7. Symmetrical
period doubling cascade.
A parametric analysis has been performed in the k, vs k, plane and the relevant results are sketched in Fig. 10; point-by-point the determination of the phase space portrait together with the largest Lyapunov exponential computation has allowed the occurrence of chaos or periodic situations to be detected. It is clearly shown how the chaotic region (dotted one in the figure) can be represented by a spot with very irregular borders; the surrounding region
is a zone of periodic operation characterized, in some cases, by very complex behaviour that makes them unacceptable from the operational point of view. The region at the top of the figure is characterized by a point operation, i.e. the set point of system (3) is an asymptotically stable singular point. What strikes the mind is that the closed-loop system behaviour is completely unpredictable: no guidance can be given to safe operation of the system
Table 2 k-7
Period/r
0.1825 0.2126 0.2190 0.2203
4.218.42 8.42116.84 16.84133.68 33.681
k.r
0.4383 0.443 0.465 0.57
Table 3 Period/r 35.8/ 17.90/35.8 8.95/17.90 4.4318.95
Chaotic
425.0 - 0.30
behaviour
, -a20
of a controlled
,
1
-0.10
0.00 KI
Y P
I
I 0.10
a20
1245
I 0.30
(7)
465.0 -_
CSTR
r
465.0 4600 t
4406 0,100
I 0.125
, 0.150
I 0175
I 0.200
I 0.225
1 0.250
I Q275
445%:
0.225
0250
0.275
KI (~1
KI (~1 Fig. 8. Period doubling bifurcations diagrams (the three diagrams are relevant to different k,r ranges).
without a previous analysis; influence of k,, k, parameters cannot he unequivocally predetermined.
CONCLUSIONS
We have proved how the presence of the PI control action in a CSTR can induce a period doubling cascade (two symmetrical ones have been detected) Table
4
k,r
cl
0.1825 0.2126 0.2190 0.2203
4.7031 4.9230
0.57 0.465
4.7727
0.433 0.4383
Table
4.68085
5
n
P*(n)
l-2 2-3 l-3 IL4 24
0.22080 0.22074 0.22075 0.22068 0.22067
34
0.22065
that, according to the Feigenbaum route, leads to chaos. All the universal features of this “scenario” have been pointed out, together with the occurrence of periodicity windows. From an engineering viewpoint this kind of thing looks so exotic that it merits no special attention at first glance. However, on closer consideration the fact that the integral action of a PI controller can give rise to a destabilizing effect of such magnitude that chaos possibly occurs, is of paramount importance. Moreover, one realizes that no suggestion can be given to the operator about the ways of modifying the control parameters, as the results shown in Fig. 10 clearly indicate that, within the parameter ranges where the chaos is experienced, no direction seems to be helpful apart from going far away. Therefore two operational conclusions stem from this work: the former is that in the domain of determinist chaos much is still unknown and remains to be done, because we are dealing with a “new science” and we are able to forecast only a few of the possible consequences, the latter lies in the awareness we must have that such things do exist, and have to be taken into due account when designing chemical reactors and control systems.
Table6 n l-2 2-3 1-3 1-Q 24 3-4
a,(n) 0.43638 0.43700 0.43689 0.43699 0.43701 0.43702
NOMENCLATURE
A = Heat transfer area (m2) C,, = Reactant concentration (kmol m-‘) c, = Specific heat (J kg-’ K-l) E = Activation energy (J kmol-‘)
L.
1246
PELLEGRINI
and G. BIARDI
0.60
380.00
400.00
420.00
440.00
460.00
Temperature
460.00
500.00
520.00
1K 1
L
420.00
Time Fig. 9. Periodicity W(T) = Floquet matrix 0= Identity matrix 9 = Jacobian matrix k, = Frequency factor (h-l) k, = Integral control constant (h-’ ) k, = Proportional control constant NTU = Number of transport units p,(r) = Length of principal axis of an n-ellipsoid Q = Volumetric flowrate (m3 h-‘) R = Gas constant (J kmol-’ K-‘) I = Time (h) r= f adimensional r T= (/= V= x=
time
Temperature, period (K), (h) Overall heat transfer coefficient (J h-’ m-r K-‘) Reaction volume (m3) Periodic solution of the autonomous system x=F(x,p) 6 = Feigenbaum constant AH, = Heat of reaction (J kmol-‘)
=(-AK&o
AT
rd
PCP
window
for k,r ~0.4. = Adiabatic temperature rise of reaction mixture after complete conversion (K) dx = Array of infinitesimal displacements from the periodic solution de
I,, ,I>, i., = Lyapunov characteristic p = Bifurcation parameter (j_ T - T,, = adimensional AT,, p = Density (kg m-l) r = Contact time (h) c = Conversion
Subscripts 0 e ref s n
= = = = =
Inlet conditions Coolant Reference conditions Set point Doubling
exponents
temperature
Chaotic KP
of a controlled
1247
CSTR
,
CL34
0.32
behaviour
...... .......... ............ ............. ............... ................ ............... ............... .............. ............ ......... ........ ....... ....... ........ ....... ........
-
0.28-
..... ....... ....... ........ ....... ........ ........ ........ ........ ........ ........ ....... ....... ....... ...... ...... ...... ..... ____
0.26 -
0.24
0.22
1
! 0.19
0.24
Fig. IO. Chaotic
1 0.29
region (dotted
REFERENCES
Aris R. and N. R. Amundson, An analysis of chemical reactor stability and control. Chem. Engng Sci. 7, 121-131 (1958). Benettin G., L. Galgani and J. M. Strelcyn, Kolmogorov experiments. Phys. Rev. A14, entropy and numerical 2338-2345 (1976). Berg& P., Y. Pomeau and C. Vidal, Order Within Chaos. Hermann-Wiley, Paris (1984). Chang H. C. and L. H. Chen, Bifurcation characteristics of nonlinear systems under conventional PID control. Chem. Engng Sri. 39, 1127-1142 (1984). Doedel E., AUTO: software for continuation and bifurcation problems in ordinary differential equations. User Manual. California Institute of Technology Pasadena, Calif. (1986). Doherty M. F. and J. M. Ottino, Chaos in deterministic systems: strange attractors, turbulence, and applications in chemical engineering. Chem. Engng Sci. 43, 139-183 (1988). Feigenbaum M. J., Quantitative universality for a class of nonlinear transformations. J. Srarisr. Phys. 19, 25 (1978). Htnon M., On the numerical computation of Poincare maps. Physica 5D, 412-414 (1982). Jorgensen D. V. and R. Aris, On the dynamics of a stirred tank with consecutive reactions. Chem. Engng Sci. 38, 44-53 (1983). Kahlert C., 0. E. Rossler and A. Varma, Modelling Chemical Reacfion Sysfems. Springer-Verlag. Heidelberg (1981). Keller H. B.. Numerical solution of bifurcation and nonlinear eigenvalue problems. Applications of Bifurcation
I
0.34
zone) in the k, r
a39
- k,
t
3
plane.
Theory (P. H. Rabinowitz, Ed.), pp. 359-384. Academic Press, New York (1977). Lorenz E. N., Deterministic nonpcriodic flow. J. Afmos. Sci. 20, 130-141 (1963). Mankin J. C. and J. L. Hudson, Oscillatory and chaotic behaviour of a forced exothermic chemical reaction. Chem. Engng Sci. 39, 1807-1814 (1984). Nandapurkar P. J., V. Hlavacek and P. Van Rompay, Chaotic behaviour of a diffusion-reaction system. Chem. Engng Sci. 41, 2747-2760 (1986). Packard N. H., J. P. Crutchfield, J. D. Farmer and R. S. Show, Geometry from a time series. Phys. RPU. Lefrs 45, 712 (1980). Pellegrini L., G. Biardi and M. G. Grottoli, Determination of the region of asymptotic stability for a CSTR. Compur ers them. Engng 12, 237-241 (1988). Prairie M. R. and J. E. Bailey, Experimental and modeling investigations of a steady-state and dynamic characteristics of ethylene hydrogenation on Pt/A&O,. Chem. Engng Sci. 42, 2085-2102 (1987). Ruelle D. and F. Takens, On the nature of turbulence. Commun. Math. Phys. 20, 167-192 (1971). Saltzman B., Finite amplitude free convection as an initial value problem--I. J. Afmos. Sci. 19, 329 341 (1962). Smith C. B., B. Kuszta, G. Lyberatos and J. E. Bailey, Period doubling and complex dynamics in an isothermal chemical reaction system. Chem. Engng Sci. 38.425430 (1983). Swinney H. L., Observation of order and chaos in nonlinear systems. Physica ‘ID, 3-15 (1983). Uppal A., W. H. Ray and A. B. Poore, On the dynamic behaviour of continuous stirred tank reactors. Chem. Engng Sci. 29, 967-985 (1974).
Compuferschem.Engng,Vol.14,No. II,pp. 1237-1247, 1% Printed in Great Britain. All rights reserved
CHAOTIC
BEHAVIOUR L.
Department
of Industrial
(Received 23 Mar&
Chemistry
Copyright
OF A CONTROLLED
PELLEGRINI
Q
1990 Pergamon
Press plc
CSTR
and G. BIARDI
and Chemical Engineering I-20133, Milan, Italy
“G. Natta”,
Politecnico
di Milano,
1990; final revisionreceived4 July 1990; received for publicarion 18 July
1990)
Abstract-The influence of a PI controller on the dynamic behaviour of a continuous stirred tank reactor (CSTR) where an exothermic reaction takes placehas been investigated. The occurrence of chaotic phenomena is singled out by means of both topological features and of sensitivity analysis to initial conditions. The Feigenbaum route to chaos is recognized, since period doubling bifurcation cascades are observed: the universal characteristics of such a “scenario” are well pointed out. The possible destabilizing effect of a PI controller on such nonlinear systems is demonstrated. A parametric analysis in terms of proportional and integral control constants (kr vs k,) shows how no general guidance can be given in tuning control parameters, within the region where chaotic phenomena may occur.
INTRODUCTION shown in various Chemical reactors, typically, are among the systems which have been most studied for detecting unusual and/or pathological behaviour in chemical engineering. Since the earliest, classical, studies in the fifties by Aris and coworkers (1958) on bifurcation phenomena in chemical reacting systems, much has been written on steady state multiplicity, orbitally stable solutions, etc. both from the experimental and from the theoretical point of view. However, only in the last decade have researchers in chemical engineering paid attention to the aperiodic or chaotic behaviour of deterministic systems, that were already known in fluid-mechanics and atmospheric science (Saltzman, 1962; Lorenz, 1963; Ruelle and Takens, 1971) from studies on “turbulent motions”. An extended review regarding chaos in deterministic systems and its applications in the chemical engineering field has been made by Doherty and Ottino (1988) where a complete description of the chaotic phenomena and of the available methods to detect them is reported. In addition extensive literature relevant to the study of aperiodic phenomena from the experimental and/or theoretical point of view is cited. CSTR has certainly been the most widely studied reacting system. Kahlert ez al. (1981) and Jorgensen and Aris (1983) pointed out that consecutive reactions in a stirred tank can induce pathological behaviour leading to chaos. A cascade of perioddoubling bifurcations in an isothermal CSTR where quadratic kinetics take place has been analyzed by Smith and coworkers (1983). Recently great interest has fields in nonlinear dynamics.
been
Always in the field of stirred tank reactors Mankin and Hudson (1984) showed how chaotic phenomena can occur under forced operation. Chang and Chen (1984) have presented the possible destabilizing effects of feedback controllers, studying a substrate inhibition model and showing how exotic dynamic behaviour can be induced by adopting a simple PI controller. Prairie and Bailey (1987) have reported the occurrence of feedbackinduced Hopf bifurcation in a closed-loop CSTR, studied for the testing and development of a kinetic model for ethylene hydrogenation. In this paper where a nonisothermal, but kinetically simple CSTR is analyzed, we give further evidence that feedback controllers which are designed, according to linear theories, to stabilize a desired operating point, may, on the contrary, induce complicated dynamic phenomena which culminate in chaos.
THE SYSTEM Let’s consider at first the open-loop system. The dynamic model of a cooled CSTR, where a first-order exothermic reaction A + B takes place, is represented by the material balance for the species A and the enthalpy balance. By adopting the usual notation (see, however, the end of the paper), we can write the following autonomous system of ordinary differential equations:
Vz =Q (CA, - C,) - k,emEiRTCAV, Vpcp~=Q~cpVo--T) -AH,k,e-EtRTCAV-UA(T 1237
(1) - T,),
L.
1238 By means of the following
V 7=--,
ad
e
T-
=
ara.q
+
koTe-EiR(~~,r+Ar,dO)(]
T,,
1 _
_
form:
5 ),
(2)
0,
-k,
AT dt, s0
to the dimensionless
form:
f At?, = -kpAe
<),
ATa, ’
The dynamic features of these equations have been theoretically studied by Uppal et al. (1974). Different operating conditions may occur: one steady state, two steady states, periodic oscillation around a singular point [each being characterized by its region of asymptotic stability: see Pellegrini et al. (1988)]. However, as discussed in the basic paper by Ruelle and Takens (1971), such a system cannot exhibit a chaotic behaviour since it shows only two state variables or, in other words, two degrees of freedom. Let’s now consider the same system with a closedloop controlling outlet temperature. It is well known that a CSTR is usually “stabilized” at the desired set point (<,, T,) by means of a SISO PI controller whose integral action allows us to asymptotically remove the offset related to the proportional control. As for the feedback control: the reactor inlet temperature T, stands for the manipulated variable whilst the outlet temperature T is the controlled variable. The proportional-integral controller presents the following mode:
equivalent
_
7
k,Te -~/~(rrcr+
AT, = -kpAT
koTe-“,R(r~f+Ar~d”)(]
TT--4%
de dr = 8,, - NTUB, - (1 + NTU)B
1
+
de
the system can be put into its adimensional
+
_(
df
f=-,
Q
=(-*HR)c*o,
-<
_= I
t
PCP
TT=
and G. BIARDI
d5
assumptions:
UA cA, - cA , NTU=epe,; A0
5=.
*T
PELLEGRINI
- k,s
that now presents three degrees of freedom. As it has been proved by Packard et al. (1980), the two formulations are perfectly equivalent. It’s worthwhile stressing that the assumed set point is the only singular point of the closed-loop system. For the set of parameters listed in Table 1 the closed-loop system exhibits an aperiodic behaviour which seems, to every appearance, chaotic from a phase space analysis as shown in Fig. I. The aim of this paper is to clearly detect the onset of chaos and its occurrence in such a simple system and to investigate the possible road to chaos or “scenario”, by means both of the analysis of geometrical features and of the methods based on sensitivity to initial conditions. THEORY
Floquet theory has been applied to find out how periodic solutions degenerate towards chaos. The theory analyzes the stability of periodic orbits; a complete description is reported by Berg& et al. (1984). Since infinitesimal perturbations are considered, a linear stability analysis is sufficient. Let x(t) be the periodic solution of period T for a nonlinear autonomous system x = F(x, p) where p is the bifurcation parameter. In order to determine the orbital stability, what happens to an initial infinitesimal disp!acement Sx from the solution is analyzed. Linearizing the system along the limit cycle, we find that the initial condition x0 + 6x is mapped after one period into x + W(T) 6x where W(T) is the Floquet matrix.
A9 df, s0 Table I
where: Ae,=e,-e,,, *e=e-e,. So the original ODE system (2) is transformed into an integro-differential one or, after differentiation, into the following ODE system:
7=36OOr k,r=4.807908 E/R=SOOOK AT,,=ZOOK T,= r,<,=430 K Th= 298.42K T, = 373.16K NTU =O.S k, = 0.3 k.r=0.3
Chaotic behaviour of a controlledCSTR
380.00
400.00
420.00
440.00
460.00
1239
480.00
500.00
520.00
Temperature (K)
520.00
_
500.00 1 480.00 E r
460.00:
J P E
440.00 7
E I-
420.00400.00<
Fig. 1. Phase plane portrait and temperature vs time diagram of the chaotic behaviour for k,r = 0.3.
The problem of the linear stability for a periodic solution is reduced to the study of the eigenvalues of W(T) that are called the characteristic Floquet multipliers (CFM). For an autonomous system one of the CFM is always unity: this corresponds to a displacement 6x along the trajectory. The periodic orbit is stable if all the other Floquet multipliers are located inside the unit circle of the complex plane. When, by continuous variation of the system bifurcation parameter, one of the eigenvalues exits in the unit circle, loss of stability of the periodic solution occurs. The type of crossing ( - 1, + 1, complex conjugated eigenvalues) has different consequences on the system behaviour and is related to the “scenario” of the route to chaos. Swinney (I 983) has reviewed in his paper the most known among such routes to chaos.
A very useful tool for detecting chaotic phenomena is also provided by the Lyapunov characteristic exponents defined as: 1, = lim 1-00
log ‘ i(t) zP,@)’
where pa is the length of the principal axis of the n-ellipsoid into which an infinitesimal n-sphere of initial conditions has been transformed due to the locally-defined nature of the flow. In particular the largest Lyapunov characteristic exponent is a measure of chaos: a positive exponent shows the system sensitivity to initial conditions (two nearby trajectories diverge in the phase-space); moreover the value of the largest Lyapunov exponent is related to the rate of divergence: the greater LCE is, the faster the system leaves its initial conditions.
L.
1240
PELLEGRINI
and G.
IO'
lo2
10
1
BIARDI
Time
-4 0
20
40
60
60
100
, 120
I
I
t
I
1
I
20
40
60
60
100
120
-60 0
Time
Time
Fig. 2. Lyapunov characteristic exponents for
In our computations we followed the method suggested by Benettin et al. (1976) that has proved itself really valid and reliable though sometimes very time consuming. SIMULATlON
RESULTS
As a first step an accurate analysis was performed for the case of the operating conditions shown in Table 1 whose exotic behaviour is well illustrated by Fig. 1. The Runge-Kutta-Merson method with selfadjusting step was employed to follow the sharp peaks of 5 vs time and of T vs time. Because of the delicate nature of the dynamic problem under study, special attention was paid to the numerical computation; to this purpose all the routines used were in double precision. In our opinion the figure represents a good portrait of chaos; however, to dispel possible doubts, the three Lyapunov exponents were computed and reported in Fig. 2a-c. Figure 2a shows that, asymptotically, the largest Lyapunov characteristic exponent 1, exhibits a positive value (1, z 0.2); the other two Lyapunov exponents AZ and A, are, respectively, 0 and ~0 according to the theory. Moreover, the Poincari: section (Fig. 3), constructed following the suggestions by Henon (1982), is typical of a chaotic situation since no preferential points of accumulation can be noticed.
k,r =
0.3.
A parametric exploration of the adimensional integral control constant k,r was performed in the phase space in order to establish the chaotic range and the possible route to chaos: the simulation results are shown in Figs 4-7. Figure 4 represents a limit cycle corresponding to a period-l solution for k,r =O.l; at k,r =0.17 (Fig. 5) one still single, but of increased amplitude, limit cycle is found (on the other hand, moving towards lower k,s values, k,r +O, more and more shrinking limit cycles are encountered). As the bifurcation parameter is increased further, the solution becomes period-2, period-4, period-8 (Fig. 5) up to k,r z 0.23 where the phase space analysis seems to indicate the occurrence of chaos.
480 r
ConversIon
Fig. 3. PoincarC map for k,r =0.3.
Chaotic
behaviour
of a controlled
CSTR
1241
o.a2= 0.80;
430.00
435.00
440.00
Temperature
(K
445.00
450.00
I
425.00
420.00 Time
Fig. 4. Phase plane portrait
and temperature
Such behaviour can be noticed in the range ~0.23 + 0.43 (see Fig. 6 for k,r = 0.35) whilst computation at higher parameter values shows limit cycles of decreasing complexity and periodicity till one single limit cycle appears beyond k,r z 0.56 (Fig. 7). The picture looks like the classic one of the Feigenbaum period doubling cascade as a route to chaos. Floquet multipliers were computed solving the system (3) simultaneously with the linear sensitivity equations given by: k=Jlw
W(O)=U,
where JI stands for the Jacobian matrix evaluated along the curve x(t, p) and 0 for the identity matrix. As suggested also by Jorgensen and Aris (1983) a check was made on the unit Floquet multiplier and on the period computed between two successive
vs time diagram
for k,~ = 0.1.
temperature minima with an error, according to Nandapurkar et al. (!986),
L.
I242
PELLEGRINI
and G. BIARDI
Fig. 5. Period doubling
showing the existence of a symmetrical period doubling cascade. The results of the Floquet multipliers computation for this last case are reported in Table 3. The use of AUTO package (Doedel, 1986) based on the continuation techniques as from Keller (1977) to check the validity of our computational analysis has given quite similar results. The period doubling bifurcations diagram from AUTO is reported in Fig. 8.
DISCUSSION
Paradoxical as it seems, it is now well known that some kind of order exists within chaos.
cascade.
Feigenbaum (1978) has shown that period doubling bifurcation occurs in a quite universal way and this may be tested for the case at hand. The first measurable property of the Feigenbaum route we tried to find out from our analysis is the f ratio:
* =LI-P”
” A-P”+,’ where p, is the parameter value corresponding to the doubling; its value, as n increases, approaches the universal constant 6 = 4.66920 1609 1. In Table 4 the k, parameter value and the calculated S, are listed for each period doubling: as it can be easily seen 6, approaches the Feigenbaum constant
Chaotic
0.40
I 380.00
I
I
51
I
I
400.00
I
behaviour
I [ 420.00
I
I
I
of a controlled
I
1 I
1 I
440.00
11
1 I
460.00
Temperature
CSTR
I
1243
I
1 I 9 1 I
480.00
1 v I
500.00
1 I 520.00
(K)
500.00 490.00 480.00 470.00 460.00 450.00 440.00 430.00 420.00 410.00 400.00 390.00 Time
Fig. 6. Phase plane portrait
and temperature
very well, also when compared with other authors results (see Smith et al., 1983). Another universal property relevant to the periodic regime preceding the chaotic region along the Feigenbaum route is that: K = pc, - const S -n, which allows, once two parameter p-values corresponding to period doubling are known, to evaluate the period doubling accumulation point pK. Table 5 reports the accumulation point values for the left to right cascade computed from each available pair of period doubling parameters; all the listed values don’t differ very much one from the other so that a mean value of pI = 0.22076 can be assumed without large errors. Similar analysis has been performed for the symmetrical cascade (see Table 6), so enabling us to
vs time diagram
for k,s = 0.35.
define the correct range of chaotic behaviour: k,r = 0.2207 + 0.437; such values are supported by punctual investigation by means of Lyapunov exponents as reported in the previous section. A remark that, in our opinion, may be of interest is that a window of periodicity is encountered for k,r z 0.40 + 0.4025. It is a narrow range of p within which the trajectories are strictly periodic. Just period-3 trajectory (Fig. 9) have been detected from our analysis due to the narrow range. Perhaps it would be possible to show how the universal behaviour can be extended (even if with different constants) to this window but it would be, in our opinion, just a piece of cleverness to go on along this way. Much more importance from the engineering viewpoint, should be accorded to the influence of the control action on the extent of the chaotic region.
L. F%LLEGRINI and G. BIARDI
Fig. 7. Symmetrical
period doubling cascade.
A parametric analysis has been performed in the k, vs k, plane and the relevant results are sketched in Fig. 10; point-by-point the determination of the phase space portrait together with the largest Lyapunov exponential computation has allowed the occurrence of chaos or periodic situations to be detected. It is clearly shown how the chaotic region (dotted one in the figure) can be represented by a spot with very irregular borders; the surrounding region
is a zone of periodic operation characterized, in some cases, by very complex behaviour that makes them unacceptable from the operational point of view. The region at the top of the figure is characterized by a point operation, i.e. the set point of system (3) is an asymptotically stable singular point. What strikes the mind is that the closed-loop system behaviour is completely unpredictable: no guidance can be given to safe operation of the system
Table 2 k-7
Period/r
0.1825 0.2126 0.2190 0.2203
4.218.42 8.42116.84 16.84133.68 33.681
k.r
0.4383 0.443 0.465 0.57
Table 3 Period/r 35.8/ 17.90/35.8 8.95/17.90 4.4318.95
Chaotic
425.0 - 0.30
behaviour
, -a20
of a controlled
,
1
-0.10
0.00 KI
Y P
I
I 0.10
a20
1245
I 0.30
(7)
465.0 -_
CSTR
r
465.0 4600 t
4406 0,100
I 0.125
, 0.150
I 0175
I 0.200
I 0.225
1 0.250
I Q275
445%:
0.225
0250
0.275
KI (~1
KI (~1 Fig. 8. Period doubling bifurcations diagrams (the three diagrams are relevant to different k,r ranges).
without a previous analysis; influence of k,, k, parameters cannot he unequivocally predetermined.
CONCLUSIONS
We have proved how the presence of the PI control action in a CSTR can induce a period doubling cascade (two symmetrical ones have been detected) Table
4
k,r
cl
0.1825 0.2126 0.2190 0.2203
4.7031 4.9230
0.57 0.465
4.7727
0.433 0.4383
Table
4.68085
5
n
P*(n)
l-2 2-3 l-3 IL4 24
0.22080 0.22074 0.22075 0.22068 0.22067
34
0.22065
that, according to the Feigenbaum route, leads to chaos. All the universal features of this “scenario” have been pointed out, together with the occurrence of periodicity windows. From an engineering viewpoint this kind of thing looks so exotic that it merits no special attention at first glance. However, on closer consideration the fact that the integral action of a PI controller can give rise to a destabilizing effect of such magnitude that chaos possibly occurs, is of paramount importance. Moreover, one realizes that no suggestion can be given to the operator about the ways of modifying the control parameters, as the results shown in Fig. 10 clearly indicate that, within the parameter ranges where the chaos is experienced, no direction seems to be helpful apart from going far away. Therefore two operational conclusions stem from this work: the former is that in the domain of determinist chaos much is still unknown and remains to be done, because we are dealing with a “new science” and we are able to forecast only a few of the possible consequences, the latter lies in the awareness we must have that such things do exist, and have to be taken into due account when designing chemical reactors and control systems.
Table6 n l-2 2-3 1-3 1-Q 24 3-4
a,(n) 0.43638 0.43700 0.43689 0.43699 0.43701 0.43702
NOMENCLATURE
A = Heat transfer area (m2) C,, = Reactant concentration (kmol m-‘) c, = Specific heat (J kg-’ K-l) E = Activation energy (J kmol-‘)
L.
1246
PELLEGRINI
and G. BIARDI
0.60
380.00
400.00
420.00
440.00
460.00
Temperature
460.00
500.00
520.00
1K 1
L
420.00
Time Fig. 9. Periodicity W(T) = Floquet matrix 0= Identity matrix 9 = Jacobian matrix k, = Frequency factor (h-l) k, = Integral control constant (h-’ ) k, = Proportional control constant NTU = Number of transport units p,(r) = Length of principal axis of an n-ellipsoid Q = Volumetric flowrate (m3 h-‘) R = Gas constant (J kmol-’ K-‘) I = Time (h) r= f adimensional r T= (/= V= x=
time
Temperature, period (K), (h) Overall heat transfer coefficient (J h-’ m-r K-‘) Reaction volume (m3) Periodic solution of the autonomous system x=F(x,p) 6 = Feigenbaum constant AH, = Heat of reaction (J kmol-‘)
=(-AK&o
AT
rd
PCP
window
for k,r ~0.4. = Adiabatic temperature rise of reaction mixture after complete conversion (K) dx = Array of infinitesimal displacements from the periodic solution de
I,, ,I>, i., = Lyapunov characteristic p = Bifurcation parameter (j_ T - T,, = adimensional AT,, p = Density (kg m-l) r = Contact time (h) c = Conversion
Subscripts 0 e ref s n
= = = = =
Inlet conditions Coolant Reference conditions Set point Doubling
exponents
temperature
Chaotic KP
of a controlled
1247
CSTR
,
CL34
0.32
behaviour
...... .......... ............ ............. ............... ................ ............... ............... .............. ............ ......... ........ ....... ....... ........ ....... ........
-
0.28-
..... ....... ....... ........ ....... ........ ........ ........ ........ ........ ........ ....... ....... ....... ...... ...... ...... ..... ____
0.26 -
0.24
0.22
1
! 0.19
0.24
Fig. IO. Chaotic
1 0.29
region (dotted
REFERENCES
Aris R. and N. R. Amundson, An analysis of chemical reactor stability and control. Chem. Engng Sci. 7, 121-131 (1958). Benettin G., L. Galgani and J. M. Strelcyn, Kolmogorov experiments. Phys. Rev. A14, entropy and numerical 2338-2345 (1976). Berg& P., Y. Pomeau and C. Vidal, Order Within Chaos. Hermann-Wiley, Paris (1984). Chang H. C. and L. H. Chen, Bifurcation characteristics of nonlinear systems under conventional PID control. Chem. Engng Sri. 39, 1127-1142 (1984). Doedel E., AUTO: software for continuation and bifurcation problems in ordinary differential equations. User Manual. California Institute of Technology Pasadena, Calif. (1986). Doherty M. F. and J. M. Ottino, Chaos in deterministic systems: strange attractors, turbulence, and applications in chemical engineering. Chem. Engng Sci. 43, 139-183 (1988). Feigenbaum M. J., Quantitative universality for a class of nonlinear transformations. J. Srarisr. Phys. 19, 25 (1978). Htnon M., On the numerical computation of Poincare maps. Physica 5D, 412-414 (1982). Jorgensen D. V. and R. Aris, On the dynamics of a stirred tank with consecutive reactions. Chem. Engng Sci. 38, 44-53 (1983). Kahlert C., 0. E. Rossler and A. Varma, Modelling Chemical Reacfion Sysfems. Springer-Verlag. Heidelberg (1981). Keller H. B.. Numerical solution of bifurcation and nonlinear eigenvalue problems. Applications of Bifurcation
I
0.34
zone) in the k, r
a39
- k,
t
3
plane.
Theory (P. H. Rabinowitz, Ed.), pp. 359-384. Academic Press, New York (1977). Lorenz E. N., Deterministic nonpcriodic flow. J. Afmos. Sci. 20, 130-141 (1963). Mankin J. C. and J. L. Hudson, Oscillatory and chaotic behaviour of a forced exothermic chemical reaction. Chem. Engng Sci. 39, 1807-1814 (1984). Nandapurkar P. J., V. Hlavacek and P. Van Rompay, Chaotic behaviour of a diffusion-reaction system. Chem. Engng Sci. 41, 2747-2760 (1986). Packard N. H., J. P. Crutchfield, J. D. Farmer and R. S. Show, Geometry from a time series. Phys. RPU. Lefrs 45, 712 (1980). Pellegrini L., G. Biardi and M. G. Grottoli, Determination of the region of asymptotic stability for a CSTR. Compur ers them. Engng 12, 237-241 (1988). Prairie M. R. and J. E. Bailey, Experimental and modeling investigations of a steady-state and dynamic characteristics of ethylene hydrogenation on Pt/A&O,. Chem. Engng Sci. 42, 2085-2102 (1987). Ruelle D. and F. Takens, On the nature of turbulence. Commun. Math. Phys. 20, 167-192 (1971). Saltzman B., Finite amplitude free convection as an initial value problem--I. J. Afmos. Sci. 19, 329 341 (1962). Smith C. B., B. Kuszta, G. Lyberatos and J. E. Bailey, Period doubling and complex dynamics in an isothermal chemical reaction system. Chem. Engng Sci. 38.425430 (1983). Swinney H. L., Observation of order and chaos in nonlinear systems. Physica ‘ID, 3-15 (1983). Uppal A., W. H. Ray and A. B. Poore, On the dynamic behaviour of continuous stirred tank reactors. Chem. Engng Sci. 29, 967-985 (1974).