Nonlinear Characterization of Pressure Fluctuations in Fluidized Beds

Ind. Eng. Chem. Res. 2008, 47, 9497–9507

9497

Nonlinear Characterization of Pressure Fluctuations in Fluidized Beds Reza Zarghami, Navid Mostoufi,* and Rahmat Sotudeh-Gharebagh Multiphase Systems Research Laboratory, Department of Chemical Engineering, UniVersity of Tehran, P.O. Box 11155/4563, Tehran, Iran

Nonlinear time series analysis techniques were applied to predict pressure fluctuation data in fluidized beds in two different hydrodynamic states. The method of delays was used to reconstruct the state space attractor to carry out analysis in the reconstructed state space. The state space reconstruction parameters, i.e., time delay and embedding dimension, were determined and the results shown that their values were different for various types of methods introduced in the literature. Chaotic behavior and predictability of fluidized system were determined by introducing two nonlinear dynamic invariants, correlation dimension and entropy, in different ways. The traditional linear autoregression method and state space based prediction methods (SSBPMs), i.e., nearest neighbors and locally linear, and global linear methods, were applied to predict the pressure fluctuation signals. The quality of prediction was assessed by comparison of the predicted data with its original benchmark. In addition, the dynamic invariants of measured and predicted attractor of the pressure signals were compared. The results showed that SSBPMs are preferred to the traditional linear methods. Finally, a continuous uncertainty band of pressure signals of single and multiple bubble regimes for the prediction methods was presented. 1. Introduction Several researchers have found characteristics of lowdimensional deterministic chaos from time series of pressure fluctuations in fluidized beds.1,2 Developed nonlinearity tests applied on bubbling fluidized bed data have confirmed the data to be nonlinear.3 Although some researchers have attributed the chaotic behavior of the pressure fluctuations in the fluidized beds to the noise in the measurements,4,5 in the present study, it was assumed that they are best approximated as chaotic, rather than purely stochastic. The state of a fluidized bed, which is a dissipative system, at a certain time, could be determined by projecting all variables governing the system into a multidimensional space, i.e., the state space. However, it is practically impossible to know all governing variables of a fluidized bed. Takens6 proved that the dynamic state of a system can be reconstructed from the time series of only one characteristic variable, such as the local pressure in a fluidized bed. Autocorrelation function8,9 and mutual information function10 have been used to determine the optimum value of the time delay for the state space reconstruction. The false nearest neighbors addressed by Kennel et al.11 and use of correlation dimension12,13 are the methods of choosing the minimum embedding dimension of a time series. An alternative method has been proposed to select the optimum value of the time window (multiple values of embedding dimension and time delay).7,9,14,15 Prediction of future observations is an important problem in the analysis of the time series. Linear autoregressive moving average (ARMA) method and some existing state space based prediction methods (SSBPMs), such as nearest neighbors,11,16 locally linear prediction,8,17,18 and global method,8 have been applied for this purpose. Usually, the first part of the data, about 50% for the single bubble regime and 80% for the multiple bubble regime, is used as the training section, and the remaining segment would be used for prediction validation. The prediction quality could be determined by normalized error. In the present * To whom correspondence should be addressed. Tel.: (98-21)66967797. Fax: (98-21)6646-1024. E-mail: [email protected].

study, maximum likelihood estimation of correlation dimension and Kolmogorov entropy12-14,19-22 of predicted time series were compared with actual pressure fluctuations in the fluidized bed in order to characterize the predicted attractor of this system. 2. Experimental Data Time series of pressure measurements in a fluidized bed extracted from the work of Johnsson et al.7 for single and multiple bubble regimes were employed in the present work. The experiments were carried out in a fluidized bed, operated under ambient conditions. Two noncirculating conditions were considered in this work: multiple bubble and single bubble regimes. The sampling frequency for pressure fluctuation signals was 400 Hz in both cases, satisfying the Nyquist criterion. The total number of data points in each sample was 786 432, corresponding to about 33 min of sampling time. They indicated that no significant distortion of the signal can be expected in the range of frequencies studied, and their measuring technique is suitable for getting dynamic information from pressure signals. A great advantage of the pressure signals is that they include the effect of many different dynamic phenomena taking place in the bed, such as bubble formation, bubble coalescence, and bubble passage. 3. State Space and Embedding Theory The state of a dynamical system can be defined in a state, such that a point in this state space specifies the state of the system. For a dissipative dynamical system such as fluidized bed,2 a set of initial conditions, after some transient time, is attracted to a subset of state spaces which is called attractor of the system.8,9 To characterize the structure of the underlying attractor from an observed time series, it is necessary to reconstruct a state space from the time series. The process of reconstructing the attractor of the system is commonly referred to as embedding.6 There are several methods of reconstructing the state space from the observed quantity which may differ in the quality of the resulting coordinates.6,23,24 A number of these methods are

10.1021/ie800460f CCC: $40.75  2008 American Chemical Society Published on Web 11/05/2008

9498 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 N-τ

ACF )

∑ [x(i) - jx][x(i + τ) - jx] i)1

N-τ

∑ [x(i) - jx]

(2)

2

i)1

where N

jx )

Figure 1. Attractor reconstruction using the method of time delays.

used to reconstruct the attractor of the fluidized bed system from pressure and void fluctuations.1,2,7,15 It is still not clear; however, which method is the most efficient. The lack of a unique solution is a source of uncertainty in nonlinear study of chaotic behavior of fluidized bed systems. The simplest method to embed scalar data, such as pressure fluctuations, is the method of delays6 which was used in this work. In this method, the state-space is reconstructed from a scalar time series, by using delayed copies of the original time series as components of the reconstructed state space (RSS). The RSS is sometimes referred to as the reconstructed state vector. Letting x(i), with i ) 1, 2, ..., N, represent the measured pressure signal at equidistant time intervals, ∆t (i.e., with a sampling frequency of fs ) 1/∆t) and N be the total number of samples, the RSS could be represented as s(i) ) (x(i), x(i + τ), ... , x(i + (d - 1)τ)) (1) where d is the number of elements of the reconstructed state vector which is equal to the number of coordinates in the reconstructed state space and is called the embedding dimension and τ is the embedding delay. It should be noted that the total number of state vectors (points) is equal to N - (d - 1)τ. A sample of the reconstructed state space method is illustrated in Figure 1. In this figure, d ) 3 and the time delay is equal to the time between the samples, i.e., τ ) 1 (corresponding to 1/fs). 4. Embedding Parameters Determinations the embedding parameters, time delay, and embedding dimension are the most important steps in the nonlinear time series modeling. Several methods have been developed for determining the time delay and the embedding dimension. However, no rigorous method exists for determining the optimal value of the time delay and embedding dimension. A number of these methods for selection of these parameters are described below and applied to pressure fluctuations time series of fluidized bed. 4.1. Time Delay. The first step in the state space reconstruction is to choose the optimum delay parameter, τ. Generally, time delay is not the subject of the embedding theorem. Therefore, from a mathematical point of view, its determination is arbitrary.2,8,9 In application, however, the proper choice of the delay time is quite important and based on geometrical arguments the attractor should be unfolded in the optimum delay time. Different prescriptions have appeared in literature to choose τ. Nevertheless, since τ has no relevance in the mathematical framework, they are all empirical and do not necessarily provide appropriate unique estimate. Some of these methods are described below: 4.1.1. Autocorrelation Function. From the basic mathematics, the autocorrelation function, ACF, compares linear dependence of two time series separated by delay and is defined as

∑

1 x(i) N i)1

(3)

The delay for the attractor reconstruction, τ, is then taken at a specific threshold value of ACF where a linear independence is found. The recommendation for this threshold is the first value of τ for which ACF equals one-half or zero or the first inflection point of ACF.8,9 4.1.2. Mutual Information. While the autocorrelation function measures the linear dependence of two variables, Fraser and Swinney10 suggested using the mutual information (MI) function as a kind of nonlinear correlation function to determine when the values of x(i) and x(i + τ) are independent enough of each other to be useful as coordinates in a time delay vector but not so independent as to have no connection which each other at all. The MI of the attractor reconstruction coordinates is defined as N-(d-1)τ

MI )

∑

P(x(i), x(i + τ), ... , x(i + (d - 1)τ)) ×

i)1

P(x(i), x(i + τ), ... , x(i + (d - 1)τ)) [ P(x(i))P(x(i + τ)), ... , P(x(i + (d - 1)τ)) ]

log

(4)

where P(x(i)) refers to the individual probability of the time series variable x(i) and P(x(i), x(i+τ), ..., x(i + (d - 1)τ)) is the joint probability density of attractor coordinate s(i) ) (x(i), x(i + τ), ..., x(i + (d - 1)τ))T. A suitable choice of time delay requires the mutual information to be minimum.1,2,8,9 In this case, the attractor is unfolded and stretched out. This condition for the choice of the delay time is known as the minimum mutual information criterion. In general, the time delay provided by the MI criteria is normally smaller than what calculated by the ACF and provides appropriate characteristic time scales for the motion. As mentioned above, the MI presents a kind of nonlinear correlation concept, while the ACF provides an optimum linear correlation criterion. 4.2. Embedding Dimension. As mentioned earlier, the dimension of a reconstructed state space is called the embedding dimension, d. However, different methods lead to different values for the embedding dimension. Working in a dimension larger than the minimum required by the data will lead to excessive requirements in terms of the number of data points and computation times necessary when investigating different issues such as nonlinear dynamic invariants calculation, prediction, etc. According to Whitney’s embedding theorem,25 every m dimensional manifold can be embedded within an embedding space of dimension d ) 2m + 1. Therefore, an attractor may, in general, be reconstructed within an embedding space of dimension d ) 2m + 1 where m is the dimension of the state space of the underlying real attractor. Alternatively, Sauer et al.26 generalized Whitney’s embedding theorem with connection to the correlation dimension. If the correlation dimension, D, of the attractor is known, then a closer estimate of the minimum embedding would be d > 2D. Several practical methods have been addressed in literature for determining the embedding dimension from an observed time series. The usual method for

Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9499

choosing the minimum embedding dimension is to calculate some invariants (correlation dimension, entropy, etc.) of the reconstructed attractor. Since these invariants are geometric properties of the attractor, they become independent of the dimension after the geometry is unfolded.8,9 Embedding dimension can be then determined when the value of the invariant stops changing and reached to its steady state status. A number of these methods are explained and applied in this work as follows. 4.2.1. Correlation Dimension. Practical method to determine a minimum dimension of the embedding space is to calculate the correlation dimension for the reconstructed attractor. The correlation dimension can be calculated from the power-law relation between the correlation integral, C, of an attractor and a neighborhood radius, ε. Generally, correlation dimension which is the slope of the log plot of C versus ε, initially increases with the embedding dimension and reaches a limiting value in which further increase in the embedding dimension would not increase the correlation dimension.12,13 The embedding dimension can be specified when the correlation dimension reach to its saturation state. Detail description of the correlation dimension is introduced in the next section. 4.2.2. False Nearest Neighbors. The method of false nearest neighbors (FNN) was developed by Kennel et al.11 This method states that if the dynamics of a system is to be reconstructed successfully in Rd, then all the neighbor points in Rd would be also neighbors in Rd+1. This method checks the neighbors in successively higher embedding dimensions until only a negligible number of false neighbors are found when increasing the dimension from d to d + 1. This dimension is chosen as the embedding dimension. 4.2.3. Time Window. Another method to select embedding parameters is to find an optimum value for the time window Tw ) dτ · ∆t. In this case, dτ can be determined for a specified sampling frequency. For the attractors systems with a dominant periodic characteristic, the time window is a segment [t, t + Tw] of the time series, and its optimum value can be chosen as one or one-quarter of the dominant period or average cycle time.7,9,14,15 Average cycle time is defined as the length of the time series (in time units) divided by the number of cycles. Number of cycles is the number of time that the time series crosses its average value, divided by two. The dominant period of a dynamic system is identified from the spectral analysis of the system in the frequency domain. However, average cycle time is calculated directly in the time domain and can be used in this work. 5. Nonlinear Dynamical Invariants After reconstructing the attractor with proper selection of the embedding parameters, the most common methods to characterize the reconstructed attractors and compare their geometric structure are the evaluation of their nonlinear dynamical invariant such as the correlation dimension and entropy. The correlation dimension expresses the spatial complexity of the attractor in state space, whereas the entropy is a measure of the predictability of the system. It should be noted that as it is mentioned before, the correlation dimension is also used for selection of embedding dimension. 5.1. Correlation Dimention. Strange attractors of chaotic systems are objects with complex geometrical structure which are self-similar at various resolutions. Self-similarity in a geometrical structure of a system is a strong signature of the chaotic behavior of underlying system. There are several methods to quantify the self-similarity of a geometrical object

by its dimension. Grassberger and Procaccia12,13 introduced the correlation dimension, D, for practical applications where the geometrical object has to be reconstructed from a finite sample of data points. Noninteger value of the correlation dimension shows self-similarity in geometrical structure of a state space attractor. In fact, when the correlation dimension is noninteger, the corresponding system has chaotic (nonlinear) behavior. The correlation dimension for simple systems (linear) is an integer. The power-law relation between the correlation integral of a reconstructed attractor and the neighborhood radius, ε, of the assumed hypersphere, C(ε) ∝ εD, can be used to provide an estimate of the correlation dimension:8,9 q(ε, M) )

∂ ln C(ε, M) ∂ ln ε

D ) lim lim q(ε, M) Mf∞ εf0

(5) (6)

where C(ε) is the correlation integral and is defined as M

C(ε) )

M

∑∑

2 Θ(ε - |s(i) - s(j)|) M(M - 1) i)1 j)i+1

(7)

in which s is a point on the attractor which has M such points; M ) N - (d - 1)τ. Θ is the Heaviside function which is equal to unity/zero if the value inside the parentheses is positive/ negative. s(i) are the points on the reference trajectory, and s(j) are other points on the attractor in the vicinity of s(i). The correlation integral is essentially a measure of the number of points within a neighborhood of radius ε, averaged over the entire attractor. 5.2. Entropy. Entropy is a well-known measure used to quantify the amount of disorder in a system. It has also been associated with the loss of information along the attractor. In fact two very close points on separate trajectories on the attractor (i.e., different initial conditions) are evolving into two different trajectories, thus, two distinguished states due to the divergence of the nearby trajectories. Therefore, the initial information would be lost after a certain time. When expressed in bits per second, entropy indicates the amount of information lost in the time unit. On the basis of the information theory, Grassberger and Procaccia19-21 stated that the required information to predict a system during the time interval [t1, t2] is I[t1,t2] ) I[t1] + K(t2 -t1)

(8)

where, I[t1] is information in bits at time t1 and K is the entropy. The value of K is zero for a predictable system (e.g., fully periodic), infinity for a stochastic system, and finite and positive for a chaotic system. The entropy can be estimated as the second order Renyi entropy (K2) and can be related to the correlation integral of the reconstructed attractor:8 K2 ∼

Cd(ε) 1 lim ln τ εf0 Cd+1(ε)

(9)

df∞

where d is the embedding dimension and τ is the time delay used for the attractor reconstruction. In a practical situation, the values of ε and d are restricted by the resolution of the attractor and the length of the time series. 6. Prediction Methods When the chaotic behavior and predictability of the system were confirmed by introducing correlation dimension and entropy then any nonlinear analysis, such as prediction, model-

9500 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008

ing, or control, would have been possible. Nonlinear prediction of future state of a system from its observable variable (e.g., pressure fluctuation in fluidized bed) is an interesting and challenging problem especially in applications where the system exhibits chaotic behavior. A variety of approaches to modeling nonlinear dynamic systems and predicting nonlinear signals have been addressed in the literature. In the following sections linear and several state space based prediction methods (SSBPMs) are explained and applied for prediction of pressure fluctuations of two different regime of fluidization system. 6.1. Autoregressive Moving Average Model. In the traditional linear models, such as the autoregressive moving average (ARMA) model, a future observation point x(i) is taken to be a linear combination of a certain number of previous observations and random, mostly Gaussian, disturbances:

Table 1. Embedding Parameters for Pressure Fluctuations Time Series of Single and Multiple Bubble Regimes regime parameters time delay embedding dimension correlation dimension entropy a

MAR

x(i) ) a0 +

MMA

∑ a x(i - k) + ∑ b η(i - j) k

j

k)1

multiple bubble

ACF mutual information correlation dimension

50 (145)a 61 not obtained

12 (24)a 22 not obtained

false nearest neighbors time window slope of log C vs log ε

16 195 1.6 ( 0.1

7 72 1.3 ( 0.2

DML Renyi (K2) Kolmogorov (KML)

2.7 not obtained 9.2 (bits/s)

7.1 not obtained 40.8 (bits/s)

The time delay at the first pass of the ACF from one-half (zero).

(10)

j)1

where η(i) are independent Gaussian random numbers with zero mean and unit variance. The constant a0 is the offset and is often omitted; ai and bj are adjusted during the training period. The first term in eq 10 is called the autoregressive (AR) term and describes the dynamics of the process. The second term is referred to the moving average (MA) term over the random input values. The numbers MAR and MMA of adjustable parameters are called the order of the process. 6.2. Nearest Neighbor Method. In the SSBPMs, the last known state of the system, which is represented by the vector s(M), is determined to predict the point x(N + τ) where M ) N - (d -1)τ is the number of points on the attractor, N is the total number of samples in the time series, and τ is the time delay. Many nearest neighbor algorithms have been proposed for predicting the state space. A simple prediction scheme, to predict the future state s(M + 1) of a given state s(M) at time tM, is to search for the closest point to s(M) with respect to a nearness norm, s(i), at time ti in the past. If the state of the system at time ti was similar to the present point, thus close in the state space, then s(i + 1) would be also close to s(M + 1). This algorithm does not introduce any modeling, and if used for more prediction steps, it would start regenerating a past part of the time series in a periodic orbit form. Another proposed nearest neighbor method is to use an averaging over all close neighbors within a specified distance.8,11,16 In this case, the time series is searched to find k similar states that have occurred in the past, where similarity is determined by evaluating the distance between vector s(M) and its neighbor vector s(i) in the d-dimensional state space. For all k points s(i), that is all points closer than ε to s(M), the future observations points in the past, x(i + τ) are looked at and used to predict x(N + τ): x(N + τ) )

single bubble

method

∑

1 x(i + τ) k s(i)∈U (s(M))

Figure 2. Autocorrelation function against delay time.

(11)

ε

where Uε(s(M)) denotes the neighborhood of radius ε around the point s(M). 6.3. Local Linear Model. In the local linear method,8,17,18 a local map can be reconstructed from the data by looking at its behavior in the neighborhood of s(M). s(M + 1) )

F

s(i)∈Uε(s(M))

(s(i))

(12)

This map of F can be approximated by fitting a polynomial which maps k nearest neighbors of s(M) onto its next immediate values, s(M + 1). The square of the Euclidean distance is chosen

Figure 3. Mutual information against delay time.

as the appropriate nearness measure. The pair (k, ε) would be selected by minimizing the error. 6.4. Global Model. Although local models generate predictions by finding local neighborhood relation and map it forward in time, global models act on the whole attractor. In the case of

Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9501

Figure 5. Percent of false neighbor against embedding dimension for single and multiple bubble regimes.

Figure 4. Correlation integral for pressure fluctuations versus ε with discarding points which are less than 500 time steps apart: (a) single bubble regime, (b) multiple bubble regime.

the global models,8 usually an appropriate function F is taken to be a linear superposition of k basic functions Φi, k

F)

∑RΦ i

i

(13)

i)1

These basic functions are fixed during fitting and only the coefficients, Ri, are fitted by minimizing the prediction error. The dynamics of a time series in global modeling can be described by a polynomial. Since F acts on a d-dimensional state space attractor, it has to be a multivariate n-order polynomial with k ) (d + n)!/d!n! independent coefficients. 7. Results and Discussion As mentioned earlier, there are different aspects that should be carefully studied before using any nonlinear analysis of an underlying system (fluidized bed in this work). These aspects are reconstructing the attractor by proper selection of embedding

Figure 6. Local slopes of correlation integral for pressure fluctuations of single bubble regime versus ε: (a) single bubble regime, (b) multiple bubble regime.

parameters, then evaluation of chaotic behavior and predictability of the system from geometrical properties of its reconstructed attractor. These aspects were calculated and discussed for two different regimes of fluidization.

9502 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008

Figure 7. Original and predicted pressure fluctuations for single bubble regime. (a) Original with global method. (b) Original with locally linear method. (c) Original with nearest neighbors method. (d) Original with linear method.

7.1. Embedding Parameters. While the selection of the time delay and the embedding dimension is a basic step for any type of nonlinear dynamic study, different methods in literature lead to different results. This conclusion has been derived from the results described below whose summary is given in Table 1. This shows that application of nonlinear time series (state space) analysis on the experimental time series (fluidized bed pressure fluctuations) is not straightforward and is still subject to research. 7.1.1. Time Delay. Two different methods of autocorrelation function and mutual information profile were used to identify the proper value of the time delay. The results are summarized in Table 1 for comparison. (a) Autocorrelation Function. The trends of the autocorrelation function (ACF) for the pressure fluctuations in single and multiple bubble regimes are shown in Figure 2. As could be seen in this figure, there is a fast decay in the autocorrelation with time lag in the multiple bubble regime while in the single bubble regime the behavior is more periodic. The first pass of the ACF from one-half and the time delay at which the ACF becomes zero occur at 50 and 145, respectively, for the single bubble regime. These values are 12 and 24 for the multiple bubble regime. These values of autocorrelation function which provide a measure of linear dependence may be chosen as the time delay in the both regimes of fluidization.

(b) Mutual Information. The mutual information profiles for both flow regimes are illustrated in Figure 3. The first minima of the mutual information for single and multiple bubble regimes occur at a delay time of 61 and 22, respectively. While the autocorrelation function compares linear dependence, these values of mutual information, provide a measure of nonlinear dependence and provide a better estimate of time delay for single and bubble fluidization regimes. 7.1.2. Embedding Dimension. The methods of correlation dimension, false nearest neighbor, and time window were used to estimate the embedding dimension of the fluidized bed attractors for the two regimes. The results are presented in Table 1 for comparison. (a) Correlation Dimension. As set before, the correlation dimension can be obtained from the slope of lines in log-log plot of C versus ε when lines are approaching to a clear saturation state. Correlation integral of pressure fluctuations in single and multiple bubble regimes were plotted vs hypersphere radius for a range of embedding dimensions. The log-log plot of the correlation sum, C(ε), versus ε (not shown here) was not quite straight line but more or less S-shaped, especially in large embedding dimensions. This trend is due to the temporal correlation problem in the definition of correlation sum, eq 7. According to Theiler et al.,27 the pairs of points in the time

Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9503

Figure 8. Original and predicted pressure fluctuations for multiple bubble regime. (a) Original with global method. (b) Original with locally linear method. (c) Original with nearest neighbors method. (d) Original with linear method.

Figure 9. Graphical representation of variable and time uncertainty.

series which are measured within a short ∆t lead to close points in the state space, thus, introduce a bias in estimation of the

correlation sum. In fact, these points are close not because of the attractor geometry, which should be considered in definition

9504 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008

Figure 10. Continuous uncertainty bands of pressure signals: (a) single bubble regime, (b) multiple bubble regime.

of correlation sum, but because they are correlated. Temporal correlations are seen in almost every data set. However, their effect on the correlation sum is more severe for high resolution time series. To solve this problem, the pairs of points which are close, not due to the attractor geometry, should be excluded. Therefore, the second summation in eq 7 should be started after a typical minimum correlation time, ∆t (including mct points in its discrete form), has passed as follows: M

C(ε) )

M

∑ ∑

2 Θ(ε - |s(i) - s(j)|) (M - mct)(M - mct - 1) i)1 j)i+m ct

(14) This correlation time, sometimes referred as Theiler window, can be estimated by introducing time separation plot as addressed by Provenzale et al.28 and Kantz and Schreiber.8 They explained that in the presence of temporal correlation, the probability that a given pair of points on the reconstructed attractor has a distance smaller than ε (geometry distance on the attractor) does depend not only on the position of the points but also on the time that has elapsed between them, ∆t. This dependence can be detected by plotting the number of neighbor points as a function of two variables, the time separation and the spatial distance. The pairs which are close in time should be excluded only for values of the time separation at saturation condition. The space time separation plot for the single and multiple bubble regimes were obtained (not shown at this work) and following results were observed. For the pressure fluctuations of the multiple bubble regime, this saturation occurs at ∆t > 100 points, therefore, a minimum correlation time of 100 points would be enough for calculating the correlation sum. The space time separation plot for the single bubble regime was more complicated. Initially (∆t < 100 points), the curves were increased fast. However, later on, a stable oscillation was found with the dominant period of the system in this regime. In the present study, to be in the safe margin, a value of 500 points was chosen for the correlation time for both regimes. Figure 4a and b illustrates the correlation integral obtained for pressure fluctuations of single and multiple bubble regimes when all pair points closer than 500 time step are discarded. These figures show straight lines as an indicator of self-similar geometry. However, these curves still do not show any rigid

saturation even for higher embedding dimensions. Some temporal saturation regions can be found and lines are separated again and do not specify any clear embedding dimension for either of the systems. This phenomenon may be connected to the noise effect on the experimental time series and the problem might be solved after the data set would be cleaned by the noise reduction scheme which is not addressed here. (b) False Nearest Neighbors. The percent of false nearest neighbors is shown versus embedding dimension in Figure 5 for both regimes of fluidization. As it can be seen, the minimum values of false neighbors percent occur in the embedding dimensions of 16 and 7 for single and multiple bubble regimes, respectively. In addition it is worth mentioning that there is a fast decay in the percent of false nearest neighbors with embedding dimension in the multiple bubble regime while in the single bubble regime the behavior is more periodic. (c) Time Window. The average cycle time was shown to be a good choice for the size of the time window for pressure fluctuation of fluidized beds.14 This is due to the fact that one cycle approximately corresponds on the main physical characteristics of the bed like bubble passage or bubble coalescence.15 In the present work, the time windows of 487.5 and 180 ms were selected for single and multiple bubble regimes, respectively, corresponding to their average cycle times. These time windows contain 195 and 72 values per cycle, respectively, for the assumed sampling frequency of 400 Hz. Therefore, the embedding dimensions of the reconstructed attractor for these regimes would be 195 and 72, respectively. As it can be seen from Table 1, the values of embedding dimension are different for various types of methods. Alternatively, Sauer et al.26 suggested that a closer estimate of the minimum embedding would be close to 2D, it was found in this work that this value could be considerably higher in some cases. 7.2. Nonlinear Dynamical Invariants. Nonlinear dynamic invariants are used to characterize the geometrical behavior of the reconstructed attractors of the system. Two nonlinear dynamic invariants of the correlation dimension and entropy were calculated in this work and are discussed below. 7.2.1. Correlation Dimension. As described above, the correlation dimension can be used for selection of a proper value for the embedding dimension. Although no specific and clear embedding dimension could be obtained from the correlation

Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9505

dimension method, local slopes, q(ε, M), of the correlation integral versus ε for the pressure fluctuations of single were used to study the chaotic behavior of corresponding time series. These slopes are shown in Figure 6a and b. Different curves in these figures represent different embedding dimensions, d. The scaling behavior, which shows self-similarity, is clearly visible in these figures. For large scales (ε > 1000), the correlation integral does not obey the power low. On the smaller length scales (90 < ε < 500), the curves are relatively flat and suggest a correlation dimension of 1.6 ( 0.1 for single bubble regime and 1.3 ( 0.2 for multiple bubble regime which confirms a low and chaotic dimensional behavior for both regimes. At the smaller scales, curves are separated again and self-similarity is obviously destroyed. There is another method of estimating the correlation dimension from the correlation integral7,14 which is simple in use for the experimental time series analysis of data such as pressure fluctuations of fluidized beds. In this method, a maximum likelihood estimation of the correlation dimension, DML, is used to characterize the geometrical structure of the attractor in the state space and the correlation dimension would be estimated from the distribution of distances between points on the reconstructed attractor.14 The maximum likelihood estimation of the correlation dimension, DML, shows that the fluctuations corresponding to passage of single bubbles through the bed are nonlinear chaos with a low noninteger dimension, 2.7, whereas for multiple bubble this value is 7.1. These values are in agreement with the values reported in the literature.1,2 It should be noted that these dimensions are different from those obtained by the direct method of using local slopes of the correlation integral versus ε. This phenomenon may be related to the noise effect on the experimental time series. The correlation dimension calculated by the maximum likelihood method gives an overestimation of the correct value for noisy data.14 These results are also summarized in Table 1 for comparison. 7.2.2. Entropy. Predictability of the pressure fluctuations can be studied from estimated entropy of chaotic fluidized bed system. Generally chaotic systems are between periodic systems, which are predictable for long time, and the stochastic systems, which are completely unpredictable. The numerical values of the Renyi entropy (K2) for both single and multiple bubble regimes were plotted (not shown here). However, a clear plateau was not visible at all and curves were not approaching to a specific limit. This might be attributed to the noise of the measurements. Another method for estimation the entropy of the fluidized beds is proposed by Schouten et al.22 This method has a simple algorithm and easy to apply for pressure fluctuations signals. Schouten et al.22 explained that entropy, as Kolmogorov entropy, can be calculated by the maximum likelihood method. The results of this method are also summarized in Table 1. The maximum likelihood estimation of the Kolmogorov entropy, KML, equal to 9.2 and 40.8 bits/s for single and multiple bubble regimes, respectively, shows an overall deterministic behavior or predictability rather than the stochastic behavior. However, it is worth mentioning that due to the large approximate entropy in the multiple bubble regime, its predictability is poor. As discussed earlier, typical characteristics of pressure fluctuations in fluidized beds are low dimension chaotic behavior for both regimes, good predictability for single bubble regime and low predictability for multiple bubble regime. Due to inherent predictability condition, it is impossible to predict the long time behavior of nonlinear dynamic systems.

Table 2. Properties of Original and Predicted Time Series multiple bubble regime

single bubble regime

method

RMSE

DML

KML

RMSE

DML

KML

original global locally linear nearest neighbors linear

25.35 655.81 684.22 950.85

2.7 2.9 3.0 2.9 1.0

9.2 10.9 9.5 7.7 9.1

24.74 107.35 409.45 850.24

7.1 7.2 6.8 6.1 5.3

40.8 40.9 40.9 35.9 32.2

7.3. Prediction. In order to build a useful phenomenological model based on the known information to estimate future observation points, the first part of the time series was considered as the training section. The last segment of series was then used for prediction. While 786 432 samples (about 33 min) were used for the attractor reconstruction, 4096 samples were taken for the prediction, corresponding to 10 s of the total sampling time. Since multiple bubble pressure fluctuations are less periodic, more training duration (about 3500 points) was considered compared to the pressure fluctuations of the single bubble for which 2000 points was enough. Since the mutual information provides a measure of nonlinear dependence, it is chosen for the time delay. While the correlation dimension do not specify any clear embedding dimension and estimated dimensions of time window method seems to give an overestimation values in both regimes, the false nearest neighbors method is used for selection of the minimum embedding dimension in the state space based prediction methods (SSBPMs). Prediction profiles of pressure signals are presented in Figures 7a-d and 8a-d. In general, the quality of prediction was found to be satisfactory. Figure 7a-d demonstrates the predicted signal versus measured time series in single bubble regime for different methods of prediction. As shown in these figures, semiperiodicity is the same as the measured values with equal frequency for single bubble regime pressure fluctuations. Figure 8a-d shows the same trends for the multiple-bubble regime which has more stochastic than the single bubble regime. The root mean squared error (RMSE) was calculated for the predicted time series compared with the original signal for each prediction method. In addition, a maximum likelihood estimation of the correlation dimension, DML, and Kolmogorov entropy, KML, of measured and predicted time series of pressure signals were compared. Although the maximum likelihood method gives an overestimation of the correct value for noisy data, the maximum likelihood estimation of the correlation dimension and entropy were used due to simplicity of these methods. The prediction errors, correlation dimensions, and Kolmogorov entropies corresponding to each method and time series are given in the Table 2. The prediction errors shown in this table suggest that the global model has a clear advantage over other methods for the fluidized bed time series. This is mostly because the fluidized bed time series is a good example of a chaotic time series and the overall character of the global model could best capture its dynamics. The locally linear predictor and the nearest neighbors show comparable performance while the linear method is the worst. In general, the nonlinear computation of predictions is relatively fast; however, it is rather slow in comparison with the linear method. Calculations were performed with a typical PC (Pentium 4, CPU 3.6 GHz, 1.0 GB RAM) for which the calculation time was less than 5 s for linear method, less than 10 s for locally linear and nearest neighbors methods, and less than 50 s for global method in the single bubble regime (2000 points for training and 2096 points for prediction); less than 5 s for linear method, less than 10 s for locally linear and nearest neighbors methods, and less than 25 s for global method in the

9506 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008

multiple bubble regime (3500 points for training and 596 points for prediction). Figure 9a-d shows that there is a continuous variable and time uncertainty for the predicted time series. Definition of the uncertainty in time and variable can be derived from Figure 9a-d. It is shown in Figure 9a and b that each point on the curve is affected by a variable error (Uv) and a time error (Ut). Owing to the uncertainty, each point may take any value within the rectangle identified by the variable and the time errors shown in Figure 9c. The amount of error, bounded by the edges of the rectangle, and the way to combine these rectangles are shown in Figure 9d. Continuous uncertainty bands of pressure signals of single and multiple bubble regimes for the above-mentioned prediction methods are presented in Figure 10a and b. Three major sources of uncertainty are reconstructed attractor (embedding dimension and time delay), measurement noise, and the prediction method. 8. Conclusions The method of reconstruction in the state space has been applied theoretically based on the fact that the dynamic state of a system can be reconstructed from the time series of only one characteristic variable such as the local pressure in a fluidized bed. For practical applications it also depends on the measurement techniques, which is considered properly. First, the state space reconstruction parameters, i.e., time delay and embedding dimension, were determined for two different hydrodynamic states. Autocorrelation function and mutual information function have been used to determine the optimum value of the time delay. The false nearest neighbors and use of correlation dimension and time window have been considered for choosing the minimum embedding dimension. It was shown that the values of the time delay and embedding dimensions were different for various types of methods. Next, the nonlinear dynamical invariants, such as the correlation dimension and entropy, were evaluated in different proposed methods. The results were different for different methods. However, in case of the single bubble regime, a nonlinear behavior with a low dimension chaotic system was observed from the correlation dimension statistics, whereas for the multiple regime, it was shown that the correlation dimensions are slightly higher. The entropy values evaluated in this work showed a deterministic (short-term predictable) behavior rather than stochastic, especially for the single bubble regime. Autoregressive moving average model and several state space based prediction methods, i.e., nearest neighbors, locally linear, and global methods were applied to predict the pressure fluctuation signals. The prediction error corresponding to each method and values of the dimension and entropy were evaluated for each method. Although the corresponding prediction errors and values of the dimension and entropy showed a good prediction, still a continuous variable and time uncertainty was observed for the predicted time series. The results showed a better prediction profile for the global model and a reasonable prediction for the nearest neighbors and locally linear models. However, for the autoregressive moving average model the prediction results was found to be poor. These results confirm the suggestion that for the chaotic time series, the SSBPMs are preferred to the traditional linear methods. Finally, it was shown that there is always a continuous variable and time uncertainty bands for the predicted time series. The major sources of this uncertainty are connected to improper selection of the embedding parameters, effect of noise in the pressure measurement, and prediction methods.

Notation ai ) constants in eq 10 AAD ) average absolute deviation ACF ) autocorrelation function AR ) autoregressive bj ) constants in eq 10 C ) correlation integral d ) embedding dimension D ) correlation dimension DML ) maximum likelihood estimation of the correlation dimension F ) map (function) I[t1] ) information in bits at time t1 m ) manifold dimension MA ) moving average k ) number of points in neighborhood of a point on attractor k ) number of independent coefficient for n order polynomial in global model K ) Kolmogorov entropy, bits/s K2 ) Renyi entropy KML ) maximum likelihood estimation of the entropy M ) number of points on attractor MI ) mutual information N ) total number samples q ) local slope of correlation integral RMSE ) root mean squared error fs ) sampling frequency, Hz s ) state vector, point on state space attractor t ) time, s tct ) minimum correlation time, s Tw ) time window Ut ) time error Uv ) variable error x ) (pressure) signal (Pa) Greek Letters Ri ) fitted coefficients in global model ∆t ) time resolution, s Θ ) Heaviside function η ) independent Gaussian random number Φi ) basic functions in global model ε ) neighborhood radius around a point on attractor ε ) geometry distance on attractor τ ) embedding time delay, s

Literature Cited (1) Van den Bleek, C. M.; Schouten, J. C. Deterministic chaos: a new tool in fluidized bed design and operation. Chem. Eng. J. 1993, 53, 75. (2) Van der Stappen, M. L. M.; Schouten, J. C.; van den Bleek, C. M. Application of deterministic chaos theory in understanding the fluid dynamic behaviour of gas-solids fluidization. AIChE Symp. Ser. 1993, 89 (296), 91– 02. (3) Van der Stappen, M. L. M. Chaotic hydrodynamics of fluidized beds. Thesis, Delft University Press, Delft, 1996. (4) Tam, S. W., Devine, M. K. Is there a strange attractor in a fluidized bed? In Measures of Complexity and Chaos; Abraham, N. B., Albano, A. M., Passamante, A., Rapp, P. E., Eds.; Plenum: New York, 1989; pp 193197. (5) Letaief, N.; Roze, C.; Gouesbet, G. Noise/chaos distinction applied to the study of a fluidized bed. J. Phys. II Fr. 1995, 5, 1883–1899. (6) Takens, F. Detecting strange attractors in turbulence. In Proceedings of Dynamical Systems and Turbulence, Lecture Notes in Mathematics 898; Rand, D. A., Yong, L. S., Eds.; Springer Verlag: Berlin, 1981; pp 366381. (7) Johnsson, F.; Zijerveld, R. C.; Schouten, J. C.; van den Bleek, C. M.; Leckner, B. Characterization of fluidization regimes by time-series analysis of pressure fluctuations. Int. J. Multiphase Flow 2000, 26, 663. (8) Kantz, H.; Schreiber, T. Nonlinear time series analysis, 2nd ed.; Cambridge University Press: New York, 2003.

Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9507 (9) Addison, P. S. Fractals and Chaos: An Illustrated Course; IOP Publishing Ltd.: London, UK, 2005. (10) Fraser, A.; Swinney, H. Independent coordinates for strange attractors from mutual information. Phys. ReV. A 1986, 33, 1134. (11) Kennel, M. B.; Isabelle, S. Method to distinguish possible chaos from colored noise and to determine embedding parameters. Phys. ReV. A 1992, 46, 3111. (12) Grassberger, P.; Procaccia, I. Characterization of strange attractors. Phys. ReV. Lett. 1983, 50, 346. (13) Grassberger, P.; Procaccia, I. Measuring the strangeness of strange attractors. Physica D 1983, 9, 189. (14) Schouten, J. C.; Takens, F.; van den Bleek, C. M. Estimation of the dimension of a noisy attractor. Phys. ReV. E 1994, 50, 3. (15) Van Ommen, J. R.; Coppens, M.; van den Bleek, C. M. Early warning of agglomeration in fluidized beds by attractor comparison. AIChE J. 2000, 46, 2183. (16) McNames, J. A nearest trajectory strategy for time series prediction. Proceedings of the International Workshop on AdVanced Black-box Techniques for Nonlinear Modeling, Katholieke Universiteit, Leuven, Belgium, July 1998; pp 112-128. (17) Casdagli, M.; Des Jardins, D.; Eubank, S.; Doyne Farmer, J.; Gibson, J.; Theiler, J. Nonlinear modeling of chaotic time series: Theory and applications. In Applied Chaos; Kim, J. H., Stringer, J., Eds.; John Wiley & Sons Inc.: New York, 1992; pp 335-380. (18) Liu, Z.; Ren, X.; Zhu, Z. Equivalence between different local prediction methods of chaotic time series. Phys. Let. A 1997, 227, 37. (19) Grassberger, P.; Procaccia, I. Estimation of the Kolmogorov entropy from a chaotic signal. Phys. ReV. A 1983, 28, 2591.

(20) Grassberger, P.; Procaccia, I. Dimensions and entropies of strange attractors from a fluctuating dynamics approach. Physica D 1984, 13, 34. (21) Grassberger, P. Estimating the fractal dimensions and entropies of strange attractors. In Chaos-Nonlinear Science: Theory and Applications; Holden, A. V., Ed.; Manchester University Press: Manchester, 1986; pp 291-311. (22) Schouten, J. C.; Takens, F.; van den Bleek, C. M. Maximum likelihood estimation of the entropy of an attractor. Phys. ReV. E 1994, 49, 126. (23) Breeden, J. L.; Packard, N. H. A learning algorithm for optimal representation of experimental data. Int. J. Bif. Chaos 1994, 4, 311. (24) Broomhead, D. S.; King, G. P. Extracting qualitative dynamics from experimental data. Physica D 1986, 20, 217. (25) Whitney, H. Differentiable manifolds. Ann. Math. 1936, 37, 645. (26) Saure, T.; Yorke, J. A.; Casdagli, M. Embedology. J. Stat. Phys. 1991, 65, 578. (27) Theiler, J.; Eubank, S.; Longtin, A.; Galdrikian, B.; Farmer, J. D. Testing for nonlinearity in time series: the method of surrogate data. Physica D 1992, 58, 77. (28) Provenzale, A.; Smith, L. A.; Vio, R.; Murante, G. Distinguishing between low dimensional dynamics and randomness in measured time series. Physica D 1992, 58, 31.

ReceiVed for reView March 21, 2008 ReVised manuscript receiVed September 15, 2008 Accepted September 23, 2008 IE800460F