Abstract
A new approach which uses the joint probability matrix computation of noisy time series is proposed to construct a phase space portrait which reflects the orbit visitation frequency of the different regions of the phase space. The resulting representation provides a clear cut of the dynamical reconstructed attractor giving both quantitative information and qualitative information about the attractor structure. The orbital distribution recovered in the map is studied by an angular first-return map where the orbital time used for the reconstruction is obtained from the magnitude information of the complex representation of the data belonging to the probability phase portrait. The resulting phase delay coordinates serve to identify phase intermittency. The Lorenz-like Shimizu-Morioka model and the Rossler model are used to present the methodology. Finally, some experimental pressure time series measured on gas-solid fluidized beds operated at different dynamical regimes are presented to analyze the reliability of the proposed method to deal with experimental noise time series.
1. Introduction
Experimental noisy time series are of paramount importance in industrial instrumentation and process control, which has a crucial role to keep the operation of the plant under control according to some specified conditions. Since a control system is essentially an information collection and transmission network, the amazing development occurring in communication and microprocessor technologies leads to several conceptual revolutions of the control system characteristics. Besides the necessity of safety, sustainable development demands tighter control. Thus, the modern day control strategies include multiple independent “control” layers to facilitate both collection and transmission of process information, and most of that information travels in the form of noisy time series.
The process control system is the first layer of safety, and especially for some complex multiphase flow operation units such as distillation column, trickled beds, bubble columns, and fluidized bed systems, the control becomes a complex issue. In particular, the complexity of the hydrodynamics gas-solid fluidized bed makes the scale-up of these systems difficult, in which many research efforts have been conducted. Most gas-solid fluidized beds operate in the bubbling regime, which is characterized by bubbles that continuously pass through the bed interacting by means of coalescing and splitting processes, to finally burst at the surface. Accordingly, the mass, heat transfer, and gas-solid mixing characteristics of this regime are a consequence of the bubble flow pattern inside the fluidized system. The low-dimensional nonlinear behavior of gas-solid fluidized beds has been extensively reported in literature by Daw et al. [1]. However, fluidized bed dynamics currently exhibit correlation dimension larger than 2, being difficult to visualize in two and even in three coordinate axes. The visualization of their dynamics reconstructing the corresponding attractors by conventional embedding techniques (i.e., differential, singular value of decomposition, PCA, or the Hilbert transform method) apparently does not provide enough information for template identification (UPOs and phase space regions). The presence of nonwandering sets can be helpful for controlling purposes as in Boccaletti et al. [2]. Moreover, the presence of experimental noise due to the system, as well as to the measurement process, makes that little useful information be extracted from these dynamical representations as well as from the well-known conventional first-return maps. For these complex systems, advanced control strategies are often designed and used to reject process disturbances in order to keep the unit under control. However, the underlying nonlinearities of multiphase flow dynamics and uncontrolled process disturbances occurring during plant performance can lead the process to become inherently unstable over a range of operational conditions, changing the process state to an “unstable” operational region. In order to deal with these unwanted situations, advanced time series analysis is commonly applied. Dynamic Matrix Control (DMC), Artificial Neural Networks (ANNs), and Fuzzy Logic (FL) models are currently used to provide a model to help process operation. Moreover, statistical tools derived from the information and from the deterministic chaos theories have shown their capability to extract dynamical information of the process dynamics. Nevertheless, there are at least two major restrictions for time series that apparently exhibit chaotic features, as follows: (i) the limitation of the black box model provided by the DMC, ANNs, and FL approaches and (ii) the unclear relation between the physical phenomena that manages the process and the nonlinear measurements used for the dynamic diagnosis. Therefore, it is necessary to get insight into the dynamical mechanisms behind the data set.
Topological analysis can be used to study the stretching and squeezing mechanism responsible for generating the measured chaotic data as in Gilmore [3]. However, it is well known that experimental time series can be contaminated by process and measurement noise, which can destroy or obscure the dynamical “skeleton” on which the true strange attractor is built, hiding the information of the unstable periodic orbits (UPOs) needed to compute the topological invariants for control purposes [2].
In this paper, a novel method based on the computation of the joint probability matrix of noisy time series is used to construct a two-dimensional picture of the orbit visitation frequency of the different regions of the phase space. To analyze the structure of the dynamical attractor, the information provided by the angular first-return map is used [4]. Thus a complex representation of the probability phase portrait serves to define a phase angle, , corresponding to the angle, in polar coordinates, between the original times series and the delay time series. Instead of constructing the first-return map from the local maxima of the phase vector, the orbital time used for the reconstruction is obtained from the magnitude information of the complex representation of the data belonging to the probability phase portrait. The angular first-return map obtained serves to clearly identify phase intermittency. Pressure time series collected from different gas-solid fluidized beds are analyzed to show the reliability of the proposed approach.
1.1. Cutting the Attractor
The dynamic of a system is characterized by its attractor. One of the most used mapping techniques is the first-return plot, which is constructed from a surrogated time series. This approach consists of obtaining the local maxima of the time series and their location in the data set, , and, once the maxima have been extracted; the first-return plot is obtained by representing versus . However, most of the real systems under study present correlation dimension larger than two, and, for the experimental time series collected in these systems, the information extracted from the first-return map is poor. In contrast, the proposed mapping methodology, based on the joint probability matrix, , reduces the -dimensional initial problem concentrating the dynamical information of the -dimensional attractor into a three-dimensional portrait, providing an optimal representation of the attractor section (Figure 1), and preserving quantitative and qualitative information about orbit visitation frequency.
1.2. The Hypothesis
The joint probability, , is defined as the probability that measurement falls into bin and its delayed version falls into bin . Then, while computing , the phase space defined by and its surrogate is filled by points . The resulting projection provides a two-dimensional picture of the dynamics defined by ; the visitation frequency over is later accounted by . Thus, if the bin size, , is small enough, such as 1% of the nominal attractor diameter [3], then the probability matrix and the plane defined by can be seen as a hypersurface, , that intersects the trajectories of the -dimensional attractor, which contains the information about the visitation frequencies of orbits crossing . The choice of the delay is critical to obtain a clean cut of the phase space. According to that, it is well known that the first minimum of the mutual information function, , gives a measurement of the decorrelation time and, consequently, its use as time lag, , guarantees the maximum independence between the coordinates and . Consequently, the resulting hypersurface provides an “optimal” cut of the trajectories of the -dimensional attractor, providing a clear cut and giving quantitative and qualitative information about the attractor structure. As shown later through the paper, any other time lag different from the will provide a “transversal” cut and therefore the resulting surface would be more contaminated by wandering or transient orbits.
1.3. The Cutting Methodology
The joint probability matrix at time lag is computed by setting a bin size, , which is estimated using the expression commonly applied for computing the noise threshold level of time series during wavelet shrinkage noise removal methods [5]. The resulting bin size has a value around 1% to 3% of the nominal attractor diameter, satisfying the condition for to be seen as a hypersurface. Moreover, as a consequence of the bin size chosen, the noise is removed from the constructed projection, and thus artificial orbits caused by the measurement noise would present a probability value close to zero. This effect influences the visualization of the mapping projection, fading out the noise to the background. Furthermore, in order to reduce the effect of the noise, the method of the relative frequency is used; such method requires the additional computation of the averaged probability value of the neighborhood of each bin (Figure 1).
1.4. The Complex Surrogates Approach
The methodology was inspired by the method reported in Daw and Halow [6], where symbolic time series were constructed accounting the crossing of the phase space trajectories. Thus, as stated above, first the -dimensional problem is reduced into a three-dimensional portrait; later the challenge is to separate the “temporal” variation (magnitude) from the spatial variation (phase). In order to concentrate, respectively, both the orbit visitation frequency and the spatial orbital information during the computation of the probability projection, the “phase space” is divided into probability regions, and the dynamical information is then retrieved by performing a component-plane crossing analysis over the surface. Thus, the dynamical information encoded within the portrait is unfolded into a spatial and a temporal orbital evolution terms by constructing a surrogate time series of “complex numbers” (Figures 2(a) and 2(b)). The approach is similar to the approach followed in [4] to investigate free liquid jets from phase portrait reconstructed with delay coordinates. Later, the new data are stored in rectangular notation preserving the filling bin order (as being located on a complex plane). During the mapping conversion, one of the bin size axes acts as the real number line while the other bin size axis acts as the imaginary part (Figure 2(c)), reducing the -dimensional problem into a two-dimensional problem.
Once the complex surrogate time series is obtained, the polar notation is preferred (in order to separate both “temporal” and the “spatial” dynamic information). On the one hand, the temporal variation of the time series is carried out in the magnitude vector; . According to that, the corresponding variations through the different spatial directions are disregarded and the information concerning visitation frequencies is enhanced by concentrating the magnitude fluctuation into just one plane (magnitude versus time). On the other hand, the spatiotemporal evolution followed by the attractor orbits through the different regions in the phase is accounted by the phase vector quantity; .
Finally, as shown below, the computation of the close return histogram over the magnitude time series is used to identify the characteristic orbital cycle time, , needed to generate a surrogated time series, , sampled from the phase vector points at , which is finally used to construct a phase-return plot representing versus ().
2. Experimental
In order to show the possibilities and reliability of the proposed close return plots construction methodology, pressure time series are analyzed. Moreover, the Rossler and Lorenz model are used for comparison, and the effect of uniform random noise on the constructed angular first-return map is investigated.
2.1. Experimental Time Series
The pressure fluctuations time series were measured by means of two AP32 KEYENCE pressure gauges. The signal was digitized through a PCI 6023E I/O board (National Instrument Co.). The experimental conditions as well as the fluidized bed used during the experimental runs are shown in [7]. Moreover, pressure time series collected by means of Kistler type 7261 transducers are shown for comparison [8].
2.2. Simulated Time Series
The Lorenz-like Shimizu-Morioka equations reported in Gilmore [3] have been integrated by using an explicit Runge-Kuta method of 4th order (Figure 3(a)). At initial conditions, . Consider The Rossler model (Figure 3(b)) has been solved by the same method at the same initial conditions compared to the Lorenz model; the equations are expressed as follows:
(a)
(b)
3. Results
The reliability of the probability projection methodology is illustrated over the Rossler and Lorenz models (Figures 4 and 5); thus several cuts at different “optimal” delay times are computed. The visual comparison of Figures 4 and 5 with the reconstructed attractors presented in Figure 3 shows that the best cut is the one provided at . The next step of the analysis is to reduce the dynamical information stored within the three-dimensional portrait into a two-dimensional problem by constructing the surrogate time series of complex numbers. The new data, which as stated above is initially stored in rectangular notation, is converted into polar form to provide the magnitude and phase information separately. Figure 6 compares the time delay histograms constructed, respectively, from the -coordinate of the simulated Lorenz and Rossler models.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
According to the literature [3], the time delay histograms are constructed by counting close returns as a function of delay; the peaks in the histograms indicate the cycle time and the approximate period of unstable periodic orbits (UPOs). The peak at the smallest time delay provides a good estimate of the dynamical cycle time. It is worth to point out that, due to the measurement and process noise contained in experimental time series, these tests cannot be used to prove that a data set is chaotic, but in contrast they can be used to reject the alternative hypothesis that it is stochastic.
As a result of concentrating the magnitude fluctuation just into one plane, the peaks detected by the time-delay histogram computed from the magnitude vectors, and , are larger than the equivalent peaks obtained from the time-delay histogram computed from the corresponding simulated -coordinate. Moreover, additional peaks appearing at the smaller time delay which are not detected by conventional time delay histogram appear in the histograms computed from the magnitude vectors (Figures 6(c) and 6(d)). When using the lag value at those local maxima as orbital time to construct a phase portrait from the simulated time series, a clear picture of the underlying attractor can be obtained (Figures 7(a) and 7(b)). Therefore, such delay time provides a good estimate of the orbital cycle time and is subsequently used to construct the corresponding angular first-return map for the Rossler and Lorenz models (Figures 7(c) and 7(d)); when no local maxima are not detected in the histogram, the first minimum of the mutual information function is used. It can be observed that, for the Lorenz model (Figure 7(d)), the angular map is around four times tangent to the main diagonal on the interval , which is relevant since as reported previously in literature the period-one orbit can be identified in the first-return map as the points closest to the diagonal [3]. In addition, it can be observed that the points are primarily confined in the curve near the main diagonal and close to a curve in the upper left corner of the angular map. According to the literature, those results could indicate the presence of an intermittent mechanism on the phase angle involving the presence of fast (regions far from the main diagonal) and slow dynamics (regions near the main diagonal) on the phase angle interval [4]. For the Rossler map, the angular points are also confined around two curves above and below the main diagonal of the first-return map and the dynamics exhibit a similar intermittent behavior. In contrast, a traditional first-return map construction obtained from the local maxima over the phase information shows no meaningful results (Figures 7(e) and 7(f)). As a result, the usefulness of the proposed method to recover information about the underlying attractor structure is clear.
(a)
(b)
(c)
(d)
(e)
(f)
To illustrate the usefulness of the proposed angle first-return map to characterize the deterministic features exhibited by the dynamics of interest, the effect of uniform random noise is explored in Figure 8. Thus, the effects of additive uniform random noise (i.e., due to measurement or instrumentation noise) are considered. Accordingly, four time series are produced by a superposition of uniform random noise with the simulated data corresponding with the -coordinate of the Lorenz model. Prior to the superposition, the noise is scaled to attain a certain signal-to-noise ratio, SNR. The SNR is defined as the ratio of the power in the noise-free signal and that of the pure noise signal. Figure 8 shows the angle first-return plot recovered from noise signal having, respectively, a SNR of 1000 (Figure 8(a)), 100 (Figure 8(b)), 10 (Figure 8(c)), and 1 (Figure 8(d)). According to the angular return maps, a signal-to-noise ratio lower than or about 10 can be considered as high noise. Thus it can be observed that the points of the first-return map are not confined on any structure for Figures 8(c) and 8(d). Moreover, for the SNR = 1 case, it can be observed that the points visit all regions of the angular first-return map (Figure 8(d)).
(a)
(b)
(c)
(d)
Finally, the potential application of the proposed methodology is illustrated, studying experimental noisy time series. Figure 9 shows the corresponding phase portrait and the angular first-return map constructed from measured pressure time series. The signals were collected from different gas-solid fluidized beds operating at multiple bubble regimes, exploding regime, single bubble regime, and slugging regime. The characteristics of those regimes have been frequently addressed in literature [7–9]. The results shown in Figure 9 are in agreement with the dynamical features exhibited for gas-solid fluidized beds when operating at those conditions. Thus, for the multiple bubble regime, the dynamics is characterized by a continuous passage of bubbles through the bed, and several bubbles can erupt simultaneously at the surface of the bed, which gives a high turbulence to the bed dynamics. Consequently, the angular first-return plot reconstructed at those conditions (Figure 9(a)) does not exhibit a clear intermittent structure similar to the structure previously observed for the Lorenz or Rossler models. Also, according to previous result shown in literature [4], the nonvisited regions on the angular return map shown in Figure 9(a) are due to a property of the phase angle and are not related to the bed dynamics. In contrast, when moving towards single bubble dominated dynamics, that is, for exploding, single bubble, and slugging regime (Figures 9(b), 9(c), and 9(d)), phase intermittency dominates. Thus, even in the most turbulent “single” bubble situation, such as the exploding regime in Figure 9(b), a characteristic curve, around which the points of the angular return map are spread, is recognizable. The best situation where the deterministic features are very well identified is for the slugging regime, which is a dynamical regime characterized by a low-dimensional bulk behavior [9]. Therefore, from the results shown in Figure 9, it is clear that although the noise time series do not exhibit a clear annular structure in the reconstructed phase portrait, the proposed angular first-return plot can reveal the underlying dynamical structure of the measured dynamics.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Moreover, how the angular first-return map shown in Figures 9(d) and 9(h) presents strong similarities between the structures around which most of the points are concentrated can be observed, such fact is remarkable since the compared dynamics were collected at very different bubbling regimes. Thereby, the exploding regime is a regime with very large gas velocities in the bottom bed, which is usually found in circulating gas-solid fluidized bed boilers. In contrast, slugging regime is a dynamical situation strongly dependent on geometrical constraints such as bed aspect ratio (bed height/bed diameter) and therefore can be found at considerable lower velocities than the exploding regime. That fact apparently reveals the existence of similar stretching and squeezing mechanism responsible for the underlying dynamics behind the measured experimental time series; such similarities might be addressed by some further symbolic sequence statistics, which is not yet fully developed for gas-solid fluidized bed, and apparently the proposed angular first-return plot might be useful to help the symbolic analysis of fluidized bed dynamics.
4. Conclusions
The results shown are very promising and can be very useful to help the topological analysis of experimental noisy time series. The possibility to construct first-return map without linear artifact correlations or measurement noise can be critical for the identification of the “true” attractor structure that characterizes the fluidized bed dynamics. Such information can be later used for symbolic analysis, for instance, to perform a template identification for modeling purposes.
The proposed methodology to construct angular first-return map uses the first minimum of the mutual information function as orbital time. Due to that, the approach is able to identify phase intermittency even from nonannular phase portraits which, as shown through the paper, is critical to identify the attractor structure hidden within the experimental measurement noise.
Notation
| : | Magnitude vector time series |
| : | Surrogate time series of local maxima |
| : | Surrogate phase vector time series |
| : | Joint probability matrix |
| : | Hypersurface |
| : | Resample time series |
| : | -coordinate. |
| : | Bin size (V) |
| : | Time lag (s). |
| : | Lorenz model |
| : | Rossler model |
| : | Sampling. |
| FB: | Fluidized bed |
| MIF: | Mutual information function |
| PCA: | Principal component analysis |
| UPOs: | Unstable periodic orbits. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors would like to specially thank Professor Filip Johnsson (Chalmers University) who kindly provided the Kistler pressure time series.