4992
J. Phys. Chem. 1995,99, 4992-5000
Control of Chemical Chaos and Noise: A Nonlinear Neural Net Based Algorithm D. Lebender, J. Muller, and F. W. Schneider* Institut j?ir Physikalische Chemie der Universitat Wurzburg, Marcusstrasse 911 1, 0-97070 Wurzburg, Germany Received: October 28, 1994; In Final Form: January 26, 1995@
A nonlinear chaos control method is presented which is based on a feedforward neural network (NN). It represents a nonlinear extension of the linear map based version of the control method by Ott, Grebogi, and Yorke (OGY). We compare the results of the nonlinear NN method with the linear map based version of the OGY method and another linear control method by Pyragas, which uses a continuous time delayed feedback algorithm. All three methods are shown by cross correlation to stabilize the same U P 0 in the same chaotic attractor. The nonlinear NN method turns out to be the fastest method to control chaos followed by the linear map based method and the linear Pyragas algorithm. Using the NN control, we also demonstrate the suppression of interactive Gaussian noise in a periodic attractor. All numerical simulations are performed with the three- and four-variable model (Montanator) of Gyorgyi and Field.
method by using the so-called Montanator m ~ d e l , ~ ~which -~O 1. Introduction is a simplified mechanism for the BZ reaction. We further The sensitive dynamical behavior of chaotic systems may be demonstrate the application of the nonlinear NN method to the controlled by imposing external perturbations or using feedback suppression of interactive noise which has been imposed on methods.'-I0 Ott,Grebogi, and Yorke proposed a new linear periodic states of the Montanator model. This emphasizes the control method (OGY method)' I which stabilizes unstable great variability of the NN method. periodic orbits (UPOS)'~contained in a chaotic attractor by In the laboratory chemical chaos is usually studied in a flow applying small time-dependent perturbations to a system reactor, for which the flow rate is a readily accessible bifurcation parameter and by monitoring any number of variables. Because parameter. Since chemical reactions may respond relatively of its simplicity, the OGY method has been implemented in s!owly to changes in flow rate, transient dynamics may be various experimental system^.'^-'^ An overview of literature observed. Long-lived transients occur in the BZ reaction as on control and synchronization of chaotic dynamical systems well as in the present Montanator model. Since transient is given in ref 17. In chemical reactions a map-based reduction behavior cannot be described by Poincark sections or oneof the OGY method was applied in theoretical and experimental dimensional maps, we have modified the map-based OGY studies of the Belousov-Zhabotinsky (BZ) reaction'8-2' and method by the use of a variation method. The variation method of electrochemical chaos,22where it is sufficient to monitor a detennines a proportionality factor between the observable and single variable. In the map-based reduction of the OGY method the controlling parameter in analogy to the map-based OGY a linear control mechanism is applied which forces a trajectory method. Chaos control by the map-based OGY method is also to move along a UPO. The relationship between the observed possible for transient behavior observed in "slow" chaotic variable and the controlling parameter is approximated by a reactions. The Pyragas method as well as the present nonlinear linear equation. Since the validity of the linearization is extension of the reduced OGY method may be applied to slow restricted to a narrow region around the UPO, it may take a chaotic systems as well. relatively long time for a trajectory to reach this region. After The important question arises whether the periodic oscillations the system's trajectory has entered this region it converges of a controlled chaotic system are stabilized UPOs or new toward the UP0 during control. Experimental chaos in the BZ states generated by the controlling algorithm itself. This and horseradish-peroxide reactions has also been ~ o n t r o l l e d ~ ~ . ~dynamic ~ question cannot be decided analytically. Since deterministic by another linear algorithm, a continuous control method of chaotic solutions can be obtained only by numerical methods, time-delayed feedback developed by P y r a g a ~ . ~The ~ . *method ~ it may be answered approximately by comparing the results of of Pyragas applies a weighted difference between the real time different control algorithms applied to the same chaotic system. and the time-delayed value of an observable, where the time If the controlled oscillations are identical for different control delay is equal to the period of the U P 0 and the weighted techniques, then they most likely represent a stabilized UP0 of difference approaches 0 during control. the chaotic system and not a newly generated dynamic state. In this work we propose a method that is a nonlinear extension The identity of the stabilized UPOs may be shown directly by of the reduced OGY method based on a small neutral net (NN). calculating their respective cross correlation functions. In independent work2' a neural net has been trained in the vicinity of the UP0 in order to control chaos in the Henon map. The latter method is based on the standard scheme of OGY, 2. Model including a variant that utilizes previous perturbations. In our The properties of the different control algorithms are tested nonlinear NN method the restriction of linearization in the in numerical simulations. As a common model, we use the vicinity of the UP0 is relaxed and the control operates over an Montanator, which was developed by Gyorgyi and Field28.29to arbitrary region around the UPO. In order to test this new simulate deterministic chaos observed in the B Z reaction. An nonlinear method, we compare its results with those of the linear account of the experimental discovery in the continuous-flow map based reduction of the OGY method and the Pyragas stirred tank reactor (CSTR) and numerical simulations based Abstract published in Advance ACS Ahstructs, March 15, 1995. on the BZ chemistry are given in ref 30. The model consists @
0022-3654/95/2099-4992$09.00/0 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 14, 1995 4993
Control of Chemical Chaos and Noise TABLE 1: Mechanism of the Four-Variable Model
+ + + + + + + + - + + + + -
Br03- Br- 2H+ HBr02 BrMA H B a 2 Br- H+ 2BrMA 1/2HBr02 Br03- H+ HBrOz f Ce4+ HBrO2 Ce4+ '/2HBr02 2HBr02 BrMA Br03- H+ Ce4+ MA inert products Ce4+ BrMA BrBrMA Br-
-
-
kl = 2.0 M-3 S-' kz = 2.0 x lo6 M-2 s-' k, = 6.2 x lo4M-* s-l k4 = 7.0 x lo3M-* s-I ks = 3.0 x lo3 M-I s-l k b = 0.3 M-' S-' ki = 30.0 M-' S-I ks = 2.4 x lo4M-I s-I
5.968-06 I
5.95e-06
+
r
-
5.948-06
0
o_ 5.939-06
5.92e-06 5.929-06 5.93e-06 5.94e-06 5.958-06 5.96e-06
[WVln
Figure 1. One-dimensional map of the chaotic time series at kf = 3.28 x
of two autocatalytic cycles, one describing the production of HBrO2 by the reduction of bromate with Ce(1V) and the other the formation of Br- from bromomalonic acid. The results of numerical simulations based on the Montanator model, which consists of 12 variables, are in good agreement with experiment. For the sake of simplicity we use reduced versions of the 12variable model, namely, a four- and a three-variable model. These simplified models show similar dynamical features shifted toward higher values of the flow rate in comparison with the full The chemical equations of the four-variable Montanator are given in Table 1. The variables of the system are bromous acid, bromide, bromomalonic acid, and cerium(IV). They are determined by a system of nonstoichiometric differential equations derived from the chemical mechanism29 referred to as model D Q ~ s The . numerical integrations of the differential equations were performed using the Gear meth0d.3'53~ The four-variable model can be simplified by a quasi-steadystate assumption for the concentration of Br-. One obtains the three-variable Montanator, as discussed in detail in ref 28. For example, the model shows deterministic chaos with a maximal positive Lyapunov exponent of 9.97 x s-' and a Hausdorff s-". This dimension of 2.3 at a flow rate kf = 3.28 x low-dimensional chaotic attractor is used in all our numerical experiments for controlling chemical chaos. In general, chaos may be reached via period d o ~ b l i n gor~ via ~ , ~the~ wrinkling of a t o r ~ s . Both ~ ~ . routes ~ ~ are based on the principle that stable orbits may bifurcate to lose their stability when the motion becomes chaotic. These unstable orbits are contained in the chaotic attractor. I * How can the UPOs be detected? A simple method is to look at the PoincarB section. If the system moves on a W O ,a point of the PoincarB plane will be mapped onto itself. For this point the corresponding one-dimensional map crosses the bisectrix, and the crossing point is an element of a UP0 (Figure 1). The aim of a control algorithm is to force the system to stay on a UP0 by methodically varying a parameter which pulls the system toward the crossing point. The control methods we present here differ in the way in which they minimize the difference between the actual trajectory and the trajectory to be stabilized. 3. Discontinuous Linear Control (Map-Based Method)
For applying the discontinuous control algorithm to a chaotic system one first determines the one-dimensional map after transients have disappeared. The one-dimensional map is obtained from the Ce(1V) values of the PoincarB section when the [HBrOz] concentration is constant (1.5 x M). Let [Ce4+], be a Ce(1V) concentration in the PoincarB section, then the kth crossing point of the system's trajectory with the
I
s-].
The crossing point with the bisectrix is an element
of a UPO.
Poincark plane is given by the one-dimensional map fk: [ce4+i,+, =fk([ce4+1,)
(1)
If [Ce4+], is a point on the UPO, then [Ce4+], will be mapped byfk onto itself [ce4+],+, = [ce4+1,
(2)
[Ce4+], is called a stable point [Ce4+],of the mapfk (Figure 1). The control algorithm involves the variation of a parameter in such a way that every point [Ce4+], is mapped onto the stable point [Ce4+],byfk. If this mapping is successful, the system is forced to remain on the UP0 and it will show a stabilized periodic motion. For controlling a chemical reaction in a CSTR, the flow rate kf is a convenient controlling parameter. In experiments it can be easily varied by adjusting the pump speed of the inflow feedstreams. In this case the map fk is a function of the flow rate kf. In a linear control algorithm the difference between an observed value of [Ce4+], and the desired value [Ce4+ls, A[Ce4+], = [Ce4'],
- [Ce4'],
(3)
has to be proportional to the variation Akf,, which is added to the flow rate kfo of the uncontrolled system: Akf, = uA[Ce4+],
(4)
Now the constant of proportionality a has to be determined. The experimental and the model system respond only slowly to a variation of the flow rate. In fact, the slow response of the system shows transient times which are much longer than the time between two consecutive crossings through the PoincarB section. Therefore, the constant of proportionality a cannot be derived from the PoincarB maps, but it has to be evaluated numerically. For controlling the chaotic system Akf should be chosen in such a way that the deviation A[Ce4+],+1 will be minimized. Therefore, the minimization problem, A[Ce4+],+, = minfk(A& A[Ce4+l,)
(5)
Akfn
has to be solved for every A[Ce4+],. The numerical solution of (5)leads to optimal changes of the flow rate A& as a function of A[Ce4+], (Figure 2). Figure 2 shows the necessary change of the flow rate when the system's trajectory crosses the PoincarB section at a certain value of [Ce4+],. In the region of A[Ce4+], FX 0 the graph in
Lebender et al.
4994 J. Phys. Chem., Vol. 99, No. 14, 1995 1.38-05 1.2e-05 Be-06
~
I1 \\
I
I
I
t
1.18-05
I -le08
0
le-OB
28-08
Start of Control
Ill I l l l,/lillIi 0
1000
A [Ce(lv)l Figure 2. Corrections of the flow rate A& versus A[Ce4+], as
2000 Time
3000
4000
2000 Time [s]
3000
4000
[SI
38-08 0.000329 I
calculated according to eq 5, which minimizes A[Ce4+],-,.
Figure 2 can be approximated by a straight line. The gradient of this line gives the proportionality constant a of eq 4. After the stable point [Ce4+],and the proportionality constant a have been determined, chaos control is turned on in the following manner: Every time the trajectory of the system crosses the Poincar6 plane, the values of A[Ce4+], and of are calculated using eqs 3 and 4, respectively. Now the flow rate is changed by adding Akf,, to the value of kf,~. The new flow rate is kept constant until the trajectory of the system crosses the Poincar6 plane for the next time. Then the same process is repeated. This control method is called discontinuous, because control takes place in pulses, i.e. when the trajectory crosses the Poincar6 plane. Due to its relatively "slow" dynamics the system is hardly affected by pulses which last for only a fraction of a whole cycle. For the three-variable Montanator a value of a = 840 M-l s-I can be derived from the graph in Figure 2. An application of eq 4 stabilizes a period 1 orbit. Figure 3 provides the result of a numerical controlling experiment. The motion of the uncontrolled system is deterministic chaotic. After 2000 s (Figure 3) the linear map based control algorithm is turned on. However, it is necessary to wait until the system crosses the Poincar6 plane with a value of Ce(IV) for which the linearization is valid. The approximate region of validity of linearity is marked with a rectangle in Figure 2. Then the flow rate can be regulated by eq 4. After a transient time the system shows a stable periodic motion. It tums out that the variations in the flow rate which are necessary for stabilizing the motion of the system are less than one-tenth of a percent. To quantify the effectiveness of convergence toward the UPO, we measured the average convergence time in 100 numerical experiments to be 381 s (Table 2).
wn
4. Continuous Linear Control (Pyragas Method)
0.0003278
0.0003274
I'
0
1000
Figure 3. Chaotic time series (top) and flow rate (bottom) of the discontinuous map based method for the three-variable model. In this example control starts at t = 2000 s and the system shows a periodic motion after a transience time of about 250 s. The period of the rectangular pulses (bottom) is equal to an orbit. For the stabilized orbit, T = 41.0 s.
TABLE 2: Effectiveness of the Convergence toward the U P 0 of the Different Control MethodsO method
mean transient time (s)
standard deviation (s)
OGY Pyragas neural net
38 1 359 131
241 98 101
a Average time to reach control (convergence time) and the standard deviation calculated in 100 numerical experiments are shown.
frequency of 2.44 x Hz, corresponding to T = 41 s. For a stabilized periodic orbit the mean quadratic deviation of Ce(lV) converges to 0: 1 h[Ce4+l2= -~"([Ce4+](t) t, - tz f l
- [Ce4+](t - ZJ)2 dt
-
0
(7)
The continuous control method developed by P y r a g a ~ ~is~ * * ~In order to find an appropriate value of /3, the mean quadratic deviation of Ce(1V) has been calculated as a function of ,6 based on the observation that if the trajectory of the system (Figure 5 ) . For values of j3 in the range between 6 and 8 M-' runs on a UPO, the difference between the observed value s-' the graph (Figure 5 ) shows a minimum. For performing a [Ce4+] at time t and the [Ce4+] value of a previous time t - T, control experiment we set ,6 = 6.0 M-' s-I and T = 41 s. Figure where Tis the period, is nearly 0. This small difference is used 6 displays the result of a numerical experiment using a chaotic for controlling the system in a linear way: state displayed by the three-variable Montanator (kf= 3.28 x s-l). Pyragas control is initiated after 2000 s. After a A@t) = j3([Ce4+l(t) - [Ce4+l(t- ZJ) (6) convergence time of 358 s (Table 2) the difference between [Ce4+](t)and [Ce4+](t - T ) declines to a low value (eq 7) and Equation 6 contains two unknown parameters: the period T of the U P 0 of period 41 s is stabilized according to eq 6. The a U P 0 and the constant of proportionality ,6. For example, a value for T can be obtained from the Fourier spectrum (Figure variations of the flow rate after stabilization are less than 5 per 4) of the uncontrolled chaotic reaction. It shows a main mille (Figure 6 ) . Since the adjustment of the flow rate takes
J. Phys. Chem., Vol. 99,No. 14, 1995 4995
Control of Chemical Chaos and Noise 88-07 t
1.38-05 '
~tartof'control
'
I
2000
3000
4000
2000
3000
4000
I
1.18-05
zi-
z
99-06
7e-06
59-06
A .
0 0
0.02
0.04 0.06 Frequency [Hz]
0.08
lo00
0
0.1
Figure 4. Fourier spectrum of the uncontrolled chaos at kf= 3.28 x s-' in the three-variable model. The main frequency is 2.44 x lo-* Hz (T = 41.0 s).
Time [SI
0.00035
0.00034
19-09 le10 0.00032
0.00031
'
0
I
0
"
I
10
20
30 40 a [l/s1/Ml
50
60
Figure 5. Pyragas method: mean quadratic deviation of Ce(IV) (eq 7) versus the feedback strength /3. Control is possible at /3 = 6-8
M-1
s-I
place continuously, the Pyragas control algorithm is termed to be continuous. 5. Neural Net Based Control
The main reason for the observed long times to reach control (Figures 3 and 6, Table 2) in both linear control methods is the linearization procedure in eqs 4 and 6 . One has to wait until the system reaches a region, for which the linearization is valid. Therefore, it is desirable to replace eqs 4 and 6 by a nonlinear equation which is valid over the total range of Ce(1V) concentrations. The dependency of the flow rate on the value of Ce(1V) is given in Figure 2. We will use a simple neural net to interpolate this dependency. Neural nets are dynamic structures made of a large number of highly interlinked processing units. They were developed according to structures found in the brain that consist of a great number of neurons interconnected by axons and synapses. The most frequently used type of neural net is the feedforward network,37historically based on Rosenblatt' s pre~ e p t r o n .The ~ ~ whole net consists of three layers of neurons (Figure 7). The input layer feeds the input pattems into the neural net. There is one intermediate layer of neurons that evaluate the weighted sum of all inputs by a certain sigmoid activation function which is independent of the specific task the net has to perform. The output signal of the previous layer
1000
Time [SI
Figure 6. Chaotic time series (top) and flow rate (bottom) for the Pyragas method (three-variable model). Control starts at r = 2000 s. The system shows a stable periodic motion (PIat /3 = 6) after a convergence time of about 700 s. For the stabilized UPO, T = 41.O s.
0 A0 0/ 0 \
I
Input Layer
HiddenLayer
/
\CY
Output Layer
Figure 7. Structure of the neural net consisting of one input, three neurons in the hidden layer, and one output neuron. Altogether six connectivities (arrows) are determined in the learning phase.
is sent forward to the next layer. This process continues until the output layer is reached. There are two phases of using a neural network: a learning phase and a working phase. The aim of the leaming phase is to adjust the weights of the connections between the neurons in such a way that a given input pattem is mapped on a corresponding output pattem. In order to train a neural net for controlling chaos, a number of points of the curve in Figure 2 are chosen arbitrarily as a training set. The A[Ce4+], values serve as input and the Awn values as output, respectively. The activation received in the input layer is processed forward through the hidden layer to the output layer. Let the neurons i be elements of one layer and the neurons j elements of the following layer. Then, the input u, to the neurons j is a linear function of the outputs vi of the neurons i and the weights wij of these connections:
49% J. Phys. Chem., Vol. 99, No. 14, 1995
Lebender et al. 1.39-05
Start of Control
The input u, to neuron j is processes by a sigmoid activation function to evaluate the output vj of this neuron: (9)
This output is propagated forward to the next layer of neurons. The adjustment of the weights is evaluated by back-propagating the calculation given by the difference between the calculated Awn and the desired Awn of the curve. We tested neural nets with one and two hidden layers, each of them consisting of up to eight neurons. The most efficient configuration for our problem is a small net consisting of an input and an output layer of one neuron each and one hidden layer of three neurons. On a common workstation the training of the net requires 2-5 min. The control process itself occurs in the microsecond range, which is faster than any other nonlinear interpolation procedure. The trained network is able to map a value of A[Ce4+], on a value of Awn according to the curve in Figure 2 as obtained by a variational method. When the trajectory of the system crosses the Poincart plane, the corresponding A[Ce4+], value is fed into the neural net. The net calculates the best value for Awn. For this case the system is forced to run on the UPO. In a numerical control experiment (Figure 8) the control of the chaotic motion is initiated after 2000 s. After a transient time of about 100 s a U P 0 of period 41 s is stabilized and the variations of the flow rate are less than 1 per mille (Figure 8). The advantage of the nonlinear neural net based control is a much shorter average convergence time of 131 s (Table 2) compared with the linear control methods. Furthermore, the method is more flexible than the linear control algorithms. It can easily be extended to a control using more than one parameter. Next we demonstrate that the method can be applied to noisy periodic time series to diminish interactive noise.
0
lo00
0
lo00
2000
3000
4000
2000
3000
4000
Time [s]
0.000329
F r
Y
8 0.0003286 4 ii 0.0003282
0.0003278
zeit [SI
Figure 8. Chaotic time series (top) and flow rate (bottom) of the neural net control (three-variable model). Control starts at t = 2000 s, and the system shows a periodic motion (Po after a convergence time of about 100 s. For the stabilized UPO, T = 41.0 s.
is the intensity of the noise. For charactxizing the effects of the added noise on the nonlinear system we draw our attention to the Poincart section. The single point in the Poincart plane which is expected for periodic motion is changed to a cloud of 6. Noise Control points by the imposed interactive noise. The size and the distribution of the points within the cloud contain information Fluctuations which exist in experimental systems can be roughly divided into interactive and noninteractive n ~ i s e . ~ O - ~ ~about the interaction between the noise and the nonlinear system. Figure 9 shows the Ce(IV) and the bromous acid values of a Noninteractive noise is produced by the detection system itself Poincark section (Br- = 2.8 x lo-' M) obtained from the and is thus not processed by the mechanism of the chemical attractor of a noisy periodic orbit ( D = 4.0 x M2 s-l). reaction. Interactive noise, however, represents fluctuations of The standard deviation 0 is a measure of the diameter of the system parameters (e.g. flow rate, temperature) and variables cloud of points. (e.g. c~ncentration).~~ The latter kind of noise may be amplified by the nonlinear system. Therefore, it is desirable to diminish noise by the method of noise control. As a model for our numerical investigations we use a periodic state of the four-variable M ~ n t a n a t o at r ~kf= ~ 5.2 x s-I. First, the noisy time series are simulated according to an The normalized third mean moment of the deviation of the algorithm developed by During the numerical integraCe(IV) values, the skewness, characterizes the degree of tion Gaussian noise is added at every single integration step to asymmetry of the cloud of points: the variables of the system: yi = y i ( l
+ ti)
i = 1, 2, ..., 4
I
1
n
\3
(10)
Since the random variables & are Gaussian distributed, their properties are sufficiently determined by
W E ; ) = ((5, where dt is the time interval (
= 2D& s) of the added noise and D
A skewness greater than 0 indicates that the points'scatter more into the region of high Ce(IV) values than into the region of low Ce(1V) values. A negative skewness indicates a distribution which is shifted toward lower values of Ce(1V). Since the skewness of the imposed Gaussian noise is equal to 0, any observed asymmetry in the distribution of the Ce(IV) values is
J. Phys. Chem., Vol. 99,No. 14,1995 4997
Control of Chemical Chaos and Noise
I
4.528-07
I
l.le-05
0
1.128-05
1.148-05
I
0
1.16~05
1.188-05
1.168-05
1.188-05
W V ) [MI
4.528-07
4.58-07 1.18-05
1.128-05
1.148-05
CdlV) [MI Figure 9. Noise control with a neural net: (a) Ce(1V) and bromous acid values of a Poincark section (Br- = 0.28 M) obtained from the attractor MZs-I) in the four-variable model; (b) Ce(IV) and bromous acid values of a Poincark section (Brof a noise-uncontrolled orbit (D = 4.0 x = 0.28 M) obtained from the attractor of the noise-controlled orbit (D = 4.0 x M2 s-') in the four-variable model.
caused by the interaction between the noise and the nonlinear system. Therefore, the skewness can be used to quantify the effect of interactive noise in a periodic time series. In our calculations the strength of noise D is varied between and 1.0 x lo-* M2 s-' in steps of 1.0 x M2 1.0 x s-I. For all simulated time series the standard deviations and skewness are calculated. The standard deviation (Figure loa, solid line) of the uncontrolled system increases with increasing strength of noise D. For D greater than 5.0 x M2 s-l the standard deviation reaches an upper limit. The skewness (Figure lob, solid line) of the distribution of the Ce(IV) concentrations is nearly 1 for all values of the noise strength. This indicates a strong asymmetric point cloud in the Poincar6 section, which is caused by the interaction between the noise and the chemical system. With increasing noise strength the value of the skewness decreases slightly. Equation 5 can be adapted to the suppression of noise. After every crossing of the trajectory with the Poincar6 section, the flow rate is varied in such a way that the next crossing of the trajectory is in the immediate neighborhood of the crossing point given by the noise-free system: (A [Ce4+],+ I ; A [HBrO,] + ,
I)
= minMA [Ce4+l,; kf
The map f describes the evolution of a point (A[Ce4+],;A-
[HBrOzl,) to its succeeding point (A[Ce4+1,+, ;A[HBrOzI,+l) of the Poincar6 section. The argument A[Ce4+], is the difference between the actual value and the value of the noise-free Ce(IV) concentration of the crossing point of the system trajectory and the Poincar6 plane. A[Ce4+], = [Ce4+], - [Ce4+l,
(15)
The concentration of bromous acid is given in an analogous way:
AtHBfl,l, = [HBQI,
- [HBfl,],
(16)
The numerical solution of eq 19 leads to a set of flow rates which are dependent on both A[Ce4+], and A[HBr02]n concentrations. The hyperplane in Figure 12 shows how the flow rate has to be chosen in order that the next point of the Poincar6 section is as near as possible to the desired point [Ce4+],. We use the hyperplane of Figure 11 to train a neural net. The Ce(1V) and bromous acid concentrations serve as input for the neural net, while the output is given by the corresponding flow rate. It turns out that a network consisting of two input neurons, a hidden layer of four neurons, and a hidden layer of two neurons followed by one ouptput neuron is best suited to this particular problem of two monitored variables. We trained the net using the back-propagation algorithm to calculate the optimal flow rate from the Ce(IV) and the bromous acid values of the Poincark section.
Lebender et al.
4998 J. Phys. Chem., Vol. 99, No. 14, I995 1.4e-07
z
I
1.2e-05
1.2e-07
E
8e-08
5
0
a
le-09 0
4e-09
7e-09
l.le-05
1e-06
1e-08
0
StMgUlofNoise
b)
2000
4000
Time [SI
6000
I
8000
I
0.6
0
0.4
I
-0.4' 1e49
"
"
"
4e-09 7e-09 strength of Noise
"
'
le-08
Figure 10. (a) Standard deviation of the uncontrolled (solid line) and neural net controlled (dotted line) system versus the strength of the interactive noise D. (b) Skewness of the uncontrolled (solid line) and neural net controlled (dotted line) system as a function of the strength of the interactive noise D. Fbwrate [l/s] 0.00058
a
0.00054 0.0005
0.00046
0
2000
4000
Time [SI
6000
8000
Figure 12. Uncontrolled time series (top) and NN controlled time series (bottom) showing the suppression of long-time correlated noise for D = 4.0 x 10-9 MZ s-1.
diameter of the cloud of points is half as large as for the uncontrolled system (Figure 9), which corresponds to a decreased standard deviation. The distribution of the points within the cloud in the controlled system is nearly symmetric, as indicated by the low skewness. This is the effect of the control of interactive noise. 7. Discussion
Figure 11. Hyperplane showing how the flow rate has to be chosen as a function of the Ce(IV) and bromous acid concentrations to optimally suppress noise of a noisy data set.
After each crossing of the Poincar6 section the noisy system is controlled by calculating a new flow rate using the neural net. This flow rate is kept constant until the system crosses the Poincar6 plane for another time. The result of a control experiment is shown in Figure 12. The upper graph contains the time series of a uncontrolled noisy periodic oscillator (D= 4.0 x M2 s-'). The dynamics of the noisy chemical reaction is dominated by long-time correlation effects. For example, between 5000 and 7000 s, the observed amplitude is lower than that expected without noise. The controlled time series (Figure 12, bottom) does not display these correlation effects. The amplitude in this time series is characterized by deviations which are statistical and independent of each other. The standard deviations of the controlled time series (Figure loa, dashed line) are smaller than those of the uncontrolled system (Figure loa, solid line). The skewness of the controlled time series (Figure lob, dashed line) decreases almost to 0. The
All chaos control methods represent approximations, since deterministic chaos cannot be described by any closed form algebraic equation. By cross correlating the results of all three control methods for the same UPO, it is possible to decide with certainty whether the state to be stabilized is a UP0 of a chaotic attractor or a newly generated state characteristic of the control method. As a result, all cross correlation functions are shown to provide maximal correlations of 1.O,which indicates that the stabilized UPOs are identical (Figure 13a,b). Since the stabilized states of all three control methods show the same period (=41 s) and the same amplitude, it is obvious that these states represent the same UP0 of the chaotic attractor (Figure 14). The amplitude of a stabilized periodic orbit is lower than the maximum amplitudes of the chaotic motion; thus the U P 0 is contained in the strange attractor. In the Pyragas method the control function (eq 7) approaches 0 when the delay time becomes equal to the period of the UPO. This is a manifestation of nonlinear resonance between the nonlinear system and its own time-delayed dynamics. The linear map based control and the Pyragas method are restricted to the region in which the linearization is valid. On the other hand, the nonlinear NN control interpolates the curve in Figure 2 also outside the linear region. Therefore, it can be
Control of Chemical Chaos and Noise
J. Phys. Chem., Vol. 99, No. 14, 1995 4999
1
a) 6
'P
1"
0
-11
0
'
'
'
'
'
300 Time [SI
-1
Time [s]
'
600
600
'
'
1
900
s
Figure 13. (a) Cross correlation between the stabilized UPOs of the map-based and the NN control method. (b) Cross correlation between the stabilized UPOs of the Pyragas and the NN control method. The period of the cross correlation function is 41 .O s. Maxima at 1.0 show that the stabilized UPOs are identical.
B 0.001175
0.001169 -0s
le-08
Figure 14. Stabilized period 1 orbit (solid) embedded in the chaotic attractor (dotted) of the Montanator model, where BM = bromomalonic acid and kf = 3.28 x s-'.
applied to all points in the Poincard section. It tums out that the NN control is relatively rapid (Table 2). However, any other efficient interpolation method would do as well. As an example, we tested a polynomial fit on the basis of Neville's algorithm and a cubic spline.46 For both methods the convergence times toward the U P 0 were nearly identical with that of the NN control. To further decrease the time for convergence, one may extend the NN algorithm to more than one controling parameter such as the inflow concentrations of selected species. By adding additional neurons, this extension is simpler to perform with a neural net than with the tested polynomial fit methods. Transient behavior is observed when a perturbation decays more slowly than the period of a cycle in the chaotic dynamics. In this case the time between two succeeding control events is less than the time for the transients to disappear. Thus, the control parameter has to be chosen in such a way as to force
the transient trajectory to coincide with the targeted UPO. This is done by the present variation method, which is used in the discontinuous linear map based method and the NN algorithm. Likewise, the Pyragas method is also suitable for transient trajectories. Noise control has been carried out in a noisy periodic orbit of the four-variable Montanator. Here the flow rate for all species is determined by two variables, the Ce(1V) and the HBr02 concentrations, which are monitored during noise control. In the noisy periodic system there is no UP0 to be stabilized. Our results show that the control algorithm reduces the extent as well as the asymmetry of the point cloud in the PoincarC section. Since the skewness of the imposed noise is 0 the skewness of the uncontrolled system is obviously generated by interactive noise. Noise control reduces only the long time correlation effects, leading to a reduction of the observed skewness. The remaining high-frequency fluctuations are a consequence of the nonamplified, noninteractive noise. Since control takes place only in the PoincarC section (about every 140 s in the chosen periodic oscillation) and the fluctuations are imposed every s, any noninteractive noise is too fast to be suppressed by the control algorithm. For an experimental implementation of chaos control both linear methods require two constants which have to be determined: a constant of proportionality and the stable point [Ce4+], in the map based method or the period of the UP0 in the Pyragas method. The stable point can be obtained from the onedimensional map after all transients have decayed, whereas the period of the U p 0 is obtained from the Fourier spectrum of the uncontrolled chaotic time series. The experimental choice of a given control method depends on the availability of the respective instrumentation. The Pyragas method tums out to be very robust toward noise, as observed in our calculations of chaos control including interactive noise, which is the subject of further work. It may be suspected that the use of more than one parameter in the OGY method would give an even better control. However, in experimental systems this is often not possible although desirable. Moreover, Figure 11 shows that monitoring an additional observable does not improve control in this case. In relatively slow chemical chaotic reactions the perturbation must be applied for a longer duration than a fraction of a cycle. In fact, in the present system the perturbation has to be applied for almost one cycle since shorter perturbations do not influence the slow chemical reaction at all. If an additional recursive feedback term is included in the OGY method, the time to reach control is likely to be increased. All discussed control methods are able to stabilize many UPOs contained in the chaotic attractor. For this purpose the parameters of the methods (stable point, delay time, constant of proportionality) have to be changed accordingly. The present methods are particularly suitable for low-dimensional chaos such as chemical chaos. In chaos control in higher dimensional systems a unique linear region may not exist. This may make it difficult to control high-dimensional chaos in a similarly simple fashion. Acknowledgment. We thank the Volkswagen Stiftung and the Deutschen Forschungsgemeinschaft for partial financial support. References and Notes (1) Hubler, A. W. Helv. Phys. Acta 1989, 62, 343. (2) Hiibler, A. W.; Georgii, R.; Kuchler, M.; Stebel, W.; Luscher, E. Helv. Phys. Acta 1989, 61, 203. (3) Singer, J.; Wang, Y.-Z.; Bau, H. H. Phys. Rev. Lett. 1991,66, 1123. (4) Minh, D. V . Phys. Rev. E 1993, 47, 714.
Lebender et al.
SO00 J. Phys. Chem., Vol. 99, No. 14, 1995 ( 5 ) Genesio, R.; Tesi, A. Auromarica 1992, 28, 531. (6) Braiman, Y.; Goldhirsch, I. Phys. Rev. Lett. 1991, 66, 2545. (7) Azevedo, A.; Rezende, S. M. Phys. Rev. Lett. 1991, 66, 1342. (8) Fahy, S.;Hamann, D. R. Phys. Rev. Lett. 1992, 69, 761. (9) Stone, E. F. Phvs. Lett. A 1992, 163, 367. (IO) Fowler, T. B. IEEE Trans. Autom. Control 1989, 34, 201. (11) Ott, E.; Grebogi, C.; Yorke, J. A. Phys. Rev. Lert. 1990, 64, 1196.
(12) Argoul, F.; Ameodo, A,; Richetti, P.; Roux, J. C. J . Chem. Phys. 1987, 86, 3325.
(13) Ditto, W. L.; Rauseo, S. N.; Spano, M. L. Phys. Rev. Lett. 1990, 65, 3211.
(14) Roy, R.; Murphy, T. W.; Maier, T. D.; Gills, Z.; Hunt, E. R. Phys. Rev. Lett. 1992, 68, 1259. (15) Hubinger, B.; Doemer, R.; Martienssen, W. Z. Phys. B 1993, 90, 103. (16) Gills, 2.; Iwata, C.; Roy, R.; Schwartz, I. B.; Triandaf, I. Phys. Rev. Lett. 1992, 69, 3169. (17) Chen, G. Control and synchronization of chaotic systems (a bibliography), EE Dept, University of Houston, TX. Available from ftp: “uhoop.egr.uh.edu/pub/reX/chaos.tex” (login and passward: “anonymous”). (18) Peng, B.; Petrov, V.; Showalter, K. J. Phys. Chem. 1991, 95,4957. (19) Peng, B.; Petrov, V.; Showalter, K. J. Chem. Phys. 1992, 96, 7506. (20) Peng, B.; Petrov. V.; Showalter, K. Physica A 1992, 188, 210. (21) Petrov, V.; GBspk V.; Masere, J.; Showalter, K. Nature 1993, 361, 240. (22) Parmananda. P.: Sherad. P.: Rollins. R. W.: Dewald. H. D. Phvs. Rev. E 1993, 47, R3003. (23) Schneider, F. W.; Blittersdorf, R.; Forster, A,; Hauck, T.; Lebender, D.; Muller, J. J . Phys. Chem. 1993, 97, 12244. (24) Forster, A.; Lekebusch, A.; Schneider, F. W. Submitted. (25) Pyragas, K. Phys. Lett. A 1992, 170, 421. (26) Pyragas, K. Naturforsch. 1993, 48a. 629.
(27) Alsing, P. M.: Gavrielides, A,: Kovanis, V. Phys. Rev. E . 1994, 49, 1225.
(28) Gyorgyi, L.; Field, R. J . Nature 1992, 355, 808. (29) Gyorgyi, L.; Field, R. J. J. Phys. Chem. 1991, 95, 6594. (30) Zhang, D.; Gyorgyi, L.: Peltier, W. R. Chaos 1993, 3, 723. (31) Gear, C. W. Numerical Initial Value Problems in Ordinary Diferential Equations: Prentice-Hall, Englewood Cliffs, NJ, 1971. (32) Shampine, L. F.; Gear, C. W. SlAM Rev. 1979, 21, 1. (33) Feigenbaum, M. J. J. Srat. Phys. 1978, 19, 25. (34) Feigenbaum, M. J. J. Stat. Phys. 1979, 21, 669. (35) Feigenbaum, M. J.; Kadanoff, L.; Shenker, S. Physica D 1982,5, 370. (36) Ostlund, S.;Rand, D.; Sethna, J. P.; Siggia, E. D. Physica D 1983, 8, 303. (37) Rumelhat, D. E.; McClelland, J. L. Parallel Distributed Processing; MIT Press: Cambridge, MA, 1986; Vols. 1, 2. (38) Rosenblatt, F. Psychol. Rev. 1955, 65, 386. (39) Rumelhart, D. E.; Hinton, G. E.; Williams, R. J. Nature 1986,323, 533. (40) Crutchfield, J. P.; Farmer, J. D.; Huberman, B. A. Phys. Rep. 1982, 92, 45. (41) Crutchfield, J. P.; Packard, N. H. Physica D 1983, 7, 201. (42) b e l , Th.-M.; Freund, A.; Schneider, F. W. J. Phys. Chem. 1990, 93, 416. (43) Fox, R. F.; Keizer, J. Phys. Rev. A 1991, 43, 1709. (44) Fox, R. F.; Gatland, I.: Roy, R.; Vemuri, G. Phys. Rev. A 1988, 38, 5938. (45) Fox, R. F. J. Stat. Phys. 1989, 54, 1353. (46) Stoer, J.; Bulirsch, R. Introduction to Numerical Analysis SpringerVerlag: New York, 1980.
JP9429 107
