Automatic Detection of Nonstationary Multiple Oscillations by an

Article pubs.acs.org/IECR

Automatic Detection of Nonstationary Multiple Oscillations by an Improved Wavelet Packet Transform Zixu Guo,† Lei Xie,*,† Alexander Horch,‡ Yixiang Wang,¶ Hongye Su,† and Xu Wang§ †

State Key Laboratory of Industrial Control Technology and ¶State Key Laboratory of Fluid Power Transmission and Control, Zhejiang University, 310027 Hangzhou, China ‡ ABB Corporate Research Center Germany, 68526 Ladenburg, Germany § Shanghai Nuclear Engineering Research and Design Institute, 200231 Shanghai, China ABSTRACT: An automatic detection technique of nonstationary multiple oscillations in industrial control loops is proposed. The technique is based on an improved wavelet packet transform (WPT) which integrates Shannon Entropy, with a proposed non-Gaussianity test and quasi-intrinsic mode function index to automatically generate the wavelet packet decomposition structure. An optimal frequency resolution of the process data for oscillation detection purposes is obtained in a form of a binary tree. Comparative studies on simulated examples as well as industrial cases indicate that the proposed approach is robust to the nonstationary features of oscillations.

1. INTRODUCTION Motivated by the fact that many control loops in the process industry perform poorly, the last two decades have witnessed a rapid development of research in control system performance monitoring and assessment.1 Performance degradation in control loops manifests due to the following aspects: (i) poor setpoint tracking, (ii) oscillations, (iii) poor disturbance rejection, and (iv) excessive final control element variation. Multiple and plantwide oscillations have detrimental effects on control loop performance and force the plant to run at suboptimal conditions.2 Therefore, it is important to detect and characterize multiple oscillations in process plants. Existing oscillation detection techniques can be roughly classified into the following categories:1 (i) time-domain approaches, e.g. integral absolute error (IAE) and autocorrelation function (ACF) methods, (ii) frequency domain approaches, e.g. fast Fourier transform (FFT) method, and (iii) hybrid approaches including wavelet transform (WT) method. Hägglund3 first proposed a technique to detect oscillations in control loops based on the IAE between the zero-crossings of the control error. Forsman and Stattin4 improved the IAE method by regularizing the upper and lower IAEs. Miao and Seborg5 developed a method with decay ratio of the autocorrelation coefficients. Thornhill et al.6 introduced an oscillation index as regularity of the zero-crossings in the ACF, but its accuracy of oscillation detection is limited by the choice of bandpass filters. Matsuo and Sasaoka7 presented an oscillation detection approach with wavelet transform. More recent techniques include discrete cosine transform (DCT) by Li8 and empirical mode decomposition (EMD) by Srinivasan,9,10 which also provide solutions for nonstationary data. Zakharov11 introduced an oscillation detection index by identifying the peak positions of the dominant frequency component to obtain a robust measurement for oscillation. Despite the recent developments in oscillation detection techniques, there are a number of practically important issues that have not been adequately addressed. For instance, many © 2014 American Chemical Society

data sets recorded from industrial control loops contain slowly varying trends, high-frequency noises and nonpersistent oscillations. A proper oscillation detection technique should be robust to such scenarios1 and provide accurate characterizations of the oscillations for further diagnosis. To this end, this article proposes an improved WPT for automatic oscillation detection when the recorded data exhibits time-variant features including nonconstant mean and oscillation amplitude. In contrast to the discrete wavelet transform (DWT), which only focuses on the approximation subspaces, WPT decomposes both the detailed and approximation subspaces to create the decomposition structure in a full binary tree. With computation complexity considered, the proposed improved WPT approach for oscillation detection employs the traditional Shannon Entropy and provides further solutions to reduce the scale of the decomposition structure with nonGaussianity tests and quasi-IMF (intrinsic mode function) index to automatically generate a parsimonious wavelet packet tree. The proposed technique has the following promising features: (i) applicability to multiple oscillations with nonconstant mean and amplitude; (ii) improved frequency resolution with reduced computation complexity; and (iii) automatic implementation without process-specific knowledge. The main contribution of this paper is presenting a best decomposition structure oriented in the area of oscillation detection, which is within the frame of wavelet packets analysis. The reliability of verifying the oscillatory components is greatly improved under this specific structure. The rest of the article is organized as follows. Section 2 presents a general review of existing work. The proposed improved WPT based oscillation detection approach is presented in section 3. Simulation examples are presented for Received: Revised: Accepted: Published: 15686

May 26, 2014 August 12, 2014 September 9, 2014 September 9, 2014 dx.doi.org/10.1021/ie502129c | Ind. Eng. Chem. Res. 2014, 53, 15686−15697

Industrial & Engineering Chemistry Research

Article

would be better to reduce the computational efforts and to improve the quality of oscillation characterization. The motivation for the proposed automatic multiple oscillation detection method centers on the shortcomings of the existing techniques. An improved version of wavelet packets analysis provides a promising tool to describe the oscillations. In regard to the requirements for the oscillation detection, wavelet packets (WPs) analysis with entropy-based best basis selection provides more compact frequency resolution with data features reserved. The quasi-IMF index and nonGaussianity test significantly reduce the scale of the decomposition structure and guarantees the accuracy of oscillation characterization. As basic requirement for industrial applications, it is important that it is applicable to noisy and time-variant data without process knowledge required.

further explanation of proposed technique in section 4. An industrial case study is discussed in section 5. A comparative study is presented in section 6 with some recently developed methods. This paper ends with a concluding summary.

2. LITERATURE REVIEW AND MOTIVATION This section briefly reviews the pros and cons of the main existing work on the process oscillation detection. The motivation for the proposed automatic approach is presented at the end of this section. Hägglund3 first raised a industrial solution to oscillation detection technique based on the calculation the integral of absolute error (IAE) between zero-crossings of control error e. t IAE is defined as IAE = ∫ tii−1|e(t)| dt, where ti and ti−1 are two consecutive instants of zero-crossings. The IAE exceeding a threshold IAElim indicates the presence of process oscillation. Thornhill12 presented an improved version, which eliminated the necessity of process knowledge. Forsman and Stattin4 introduced an alternative IAE approach by calculating the IAEs for positive and negative control errors separately, and the threshold is calculated in a recursive manner. Another category of oscillation detection technique is in the ACF domain. The ACF of an oscillating signal is also oscillatory with the same period, and it is robust to process noise.6 The approach by Miao and Seborg5 utilizes the decay ratio of the ACF as an index to evaluate significant level of oscillation. Thornhill et al.6 suggested to use the zero-crossings of the ACF to detect the oscillations. Oscillations can be detected and characterized by calculating the regularity of the interval periods. Regularity r is defined as a dimensionless factor, r = T̅ p/3σTp. T̅ p and σTp denote the average and standard deviation of the interval periods, respectively. In spite of the successful applications of the above approaches, their performance may deteriorate due to the following characteristics of the industrial data: (i) nonstationary trends induced by setpoint change or external disturbances, (ii) non-Gaussian distributed noise, and (iii) multiple oscillations due to the control loops couplings.9 The quantified description of the above features is essential for further diagnosis procedure in a comprehensive oscillation detection method. IAE cannot determine whether the data is multioscillatory. More precisely, only the dominant frequency contributes to IAE value. The decay ratio approach primarily aims at detecting the presence of oscillation and is not able to distinguish multi oscillation frequencies either. Regularity of zero-crossings of ACF provides a possibility to indicate each oscillation pattern in a statistical way. But the presence of nonstationary trends and multiple oscillations may mislead the result unless proper data filters are carefully selected.1,6 Besides, there is a restriction on the Fourier transform that the inspected process must be piecewise stationary. Unfortunately, it is not always fulfilled in the real industrial applications. Both DCT and EMD, along with their improved versions share a common shortcoming when extracting the oscillation components. Their decompositions lack of specific improvements for the efficiency of oscillation detection purpose. Their decomposition structures remain the same with the original versions which may extract many unnecessary components. Although EMD and DCT are alternatives of oscillation detection techniques, a better decomposition structure under properly designed basis functions for oscillation extraction

3. DESCRIPTION OF THE METHODOLOGY This section introduces the proposed technique for multiple oscillation detection. A brief summary of the wavelet packet

Figure 1. Topology of the full rooted binary tree.

transform with entropy-based basis selection, non-Gaussianity test, and quasi-IMF index (QIMFI) is introduced as preliminaries before the description of the proposed technique. 3.1. Preliminaries. Wavelet Packet Transform with Entropy-Based Basis Selection. Wavelet packets (WPs) and general multirate filter banks have emerged as important signal representation schemes. The basis family of WPs is particularly appealing for the analysis of pseudostationary time series processes. With this promising feature, WPs are able to deal with multioscillatory data via extraction and decomposition. The top level of the WPT is the time representation of the signal, whereas the bottom level has a better frequency resolution with decomposition structure.13 WPs allows decomposing the observation space into subspaces associated with different frequency bands. Its basis family is characterized by a tree structure induced by its filter bank implementation, that recursively iterates through a two channel orthonormal filter bank. In the process of cascading a basic block of analysis, it is possible to generate a collection of orthonormal bases {ψij} for the space of finite energy sequences associated with different time-scale representation properties. The basic block of generating {ψij} is 1 2

ψ2ij+ 1(t ) = ψ2ij++11(t ) 15687

∑ h(k)ψji⎛⎝ t ⎜

k

1 = 2

⎞ − k⎟, ⎠ 2

∑ g(k)ψji⎛⎝ t ⎜

k

2

⎞ − k⎟ ⎠

(1)

dx.doi.org/10.1021/ie502129c | Ind. Eng. Chem. Res. 2014, 53, 15686−15697


Article

xji(t ) = c jiψ ji(t )

(3)

Let f(t) ∈ ? = 1 be the raw observation space. The application of the basic block of analysis decomposes ? into two subspaces ?10 and ?11, and orthonormal bases ensure ? = ?10 ⊕ ?11 (⊕ is direct sum of spaces). In any of the resulting sub-band spaces, ?10 and ?11, reapply the basic block of analysis to generate a new indexed basis. By iterating this process, it is possible to construct a binary tree-structured collection of indexed bases for ? . For instance, by iterating the decomposition recursively for l-times in every sub-band space, we can generate the indexed basis ∪j∈{0,...,2i−1}{ψij}, where ? = ⊕j∈{0,...,2i−1}? ij and ? ij = span{ψij} satisfying ? ij = ? i2+j 1 ⊕ ? i2+j +11

(4)

The generative process of producing a particular indexed basis in the WP family can be presented by a rooted binary tree. Let lmax be the maximum number of iteration times of this subband decomposition process, which denotes the maximum level of the binary tree. Let G = (E, V) be a graph with indexes E = {(0, 0), (1, 0), (1, 1), ..., (l, 0), ..., (l, 2l − 1), ..., (lmax, 0), ..., (lmax, l, 2lmax − 1)} and V the collection of arcs (degrees) on E × E that characterizes a full rooted binary tree (denoted by Tfull) with ROOT(0, 0) as illustrated in Figure 1. A rooted binary tree contains one root node with two degrees, internal nodes with three degrees, and leaf nodes with only one degree. Orthogonal WPs have a well-controlled time-frequency localization property, which is exactly determined by defining the wavelet packet tree structure. The intension of this study is to improve WP decomposition and make it optimal for further oscillation detection purpose, under the orthogonal basis of WPs with the minimum information required. One of the most popular and traditional way of achieving this goal is adopting cost functions to measure the information quantity in nodes with Shannon Entropy as15 Ei , j = −∑ c ji(k)2 log(c ji(k)2 ) k

In this way, a decomposition ? ij = ? i2+j 1 ⊕ ? i2+j +11 is considered to be necessary only if Ei,j > Ei+1,2j + Ei+1,2j+1, which indicates the information cost of the new collection of subspaces is reduced. Non-Gaussianity Test. A non-Gaussianity test is proposed which measures negentropy values of data in our previous work.16 It is proper alternative for the test in this article because entropy values are adopted during the WPs basis selection. Since Gaussian variable v has the largest entropy among all the random variables of unit variance, this allows the definition of the negentropy J(y) = H(v) − H(y). H(y) denotes to entropy function H(y) = −∫ f(y) log(f(y)) dy, where y and f(·) is a random variable and its density function, respectively. J(y) can be approximated by

Figure 2. Flowchart of identifying the terminal nodes.

where i is the modulation parameter or the decomposition level of WPT (i = 0, 1, 2, ..., lmax), j is the dilation parameter (j = 0, 1, 2, ..., 2i − 1) and k is the translation parameter. ψ00(t) is also known as mother wavelet. The two-channel discrete filters h(k) and g(k) are quadrature mirror filters (QMF) banks. A QMF pair splits an input signal into two bands with quadrature mirror property and power complementary property, which ensures the possible reconstruction. The selection of h(k) and g(k) is wavelet-dependent.14 The wavelet packet coefficients sequence cij corresponding to the signal f(t) can be calculated by projections onto spaces spanned by orthogonal wavelet bases {ψij} c ji =

∫ f (t )ψji(t ) dt

(5)

J(y) ≈ [E{G(y)} − E{G(v)}]2

(6)

where G(·) is the nonquadratic function and is preferred to be −exp(−y2/2).17 The negentropy function of y follows a Normal distribution y ∈ 5 (0,1), and the following holds as the sample number K → ∞,

(2)

xij(t)

The subportion of f(t) can be obtained through a reconstruction process

K ·J(y) ∼ var{G(v)}χ 2 (1) 15688

(7)


Industrial & Engineering Chemistry Research ⎡1 J (y ) = ⎢ ⎢⎣ K

⎤2 ∑ G(yi ) − E(G(v))⎥ ⎥⎦ i=1

Article

oscillating pattern. A symmetrical wavelet is preferred as a mother wavelet for industrial data to guarantee symmetrical and regular shaped decomposed components.18 Orthogonal wavelets, Daubechies and Coiflets wavelet of high orders (“db5” or “coif5”) as an example, show good symmetry and are suitable for the oscillating time series to be analyzed. In the case of the original WPT with a full rooted binary tree Tfull, the decomposition at level n contains 2n subspaces. Shannon Entropy enables a reduction of the decomposition tree complexity and yields a pruned subtree T of Tfull. If a split at a particular node results in an increase to the total information cost of T, i.e. E(i,j) < E(i+1,2j) + E(i+1,2j+1), it suggests that the split is costly and node (i, j) should be preserved as a leaf node. The difference in this article from traditional mechanism is that in our algorithm the entropy is calculated and compared along with the generation of T. Shannon Entropy is working as a necessary condition of decomposing a subspace rather than pruning Tfull when it is obtained. Note that leaf nodes in T, optimized by entropy do not necessarily exhibit oscillatory behaviors. For process data containing noise in the high frequency band, it is unnecessary to decompose the noise subspaces further for oscillation detection purpose. However, entropy test itself cannot tell the presence of a pure noise segment and the decomposition will continue until the compound noises are divided into narrower frequency bands with no oscillation exists. In order to reduce the scale of T further, it is suggested to evaluate the reconstructed signal with the accompanying criteria. The nonGaussianity test can identify whether the node contains information about oscillation or it is simply a noise. In contrast, quasi-IMF index can cease the decomposition where the reconstructed signal from the node shape like an IMF, which does not contain multiple oscillation frequencies. To summarize, a leaf node (i, j) in the wavelet packet decomposition tree splits further into (i + 1, 2j) and (i + 1, 2j + 1) during the generation of the WPT structure, if it satisfies the following criteria:

K

(8)

where v is a Gaussian variable statistically independent of y and var{·} denotes variance. In order to determine whether y follows a Gaussian distribution, a confidence limit for the statistic TG = (K/var{G(v)})·J(y) can be computed on the basis of a ?2 distribution with one degree of freedom, and a preferred significance level is 0.05. Quasi-IMF Index. The empirical mode decomposition (EMD) algorithm obtains the intrinsic mode functions (IMFs) of the time series data to identify the presence of single/multiple oscillations. The individual IMFs should satisfy the following properties: (i) the number of maxima and minima points must be equal or differ by one from the number of zerocrossings and (ii) at any instant in time, the mean value of the envelopes defined by the local maxima and local minima is zero or, in other words, the amplitude between each consecutive extrema points must be symmetric. A sine wave is a simple and intuitive example of an IMF which can be concluded as an oscillation. A pseudo random binary sequence (PRBS) is also an IMF in a binary form, but it is not a regular oscillation. Considering a subspace ? ij in eq 4 which exhibits timevarying amplitudes, it may not fulfill the second property to be a symmetric component. However, the same extrema points and zero-crossings in numbers can still imply a single oscillation pattern. It would be unnecessary to decompose such subspace because there is no possibility of multiple oscillation’s existence. In this study, a quasi-IMF index (QI) is proposed according to the first property of IMF

QI(i , j) = |Nex − N0|

(9)

where Nex and N0 denotes the number of extrema and zerocrossings of ? ij in time series. If QI(i, j) ≤ 1, the reconstructed subsignal from node (i, j) shows potential to exhibit an single oscillation pattern. It is unnecessary to split this node, and it should be marked as a “terminal” leaf node. In addition, it is the task of ACF evaluation, which will be illustrated in section 3.2, to identify whether the leaf node components are regular enough to be identified as oscillatory. Note that QI ≤ 1 is the necessary but not sufficient condition for the existence of oscillation. A counterexample is PRBS signal, whose QI(i, j) ≤ 1, but it is not regular enough to determine the existence of oscillation. 3.2. Main Idea of the WPT-Based Method. The main idea of the WPT-based method is to isolate different frequency components of a time series via an improved WPT algorithm and to detect oscillations by checking the regularity of ACF zero-crossings of these isolated components. The improved WPT ensures a computationally efficient decomposition of the original data for oscillation detection. Three sifting criteria, namely Shannon Entropy, Non-Gaussianity test, and quasi-IMF index, are adopted along with the iteration of WP analysis which results in the optimal decomposition of the time series for oscillation detection. Improved Wavelet Packet Decomposition. The raw observation space ? is decomposed into subspaces by improved WPT which is oscillation detection oriented. In order to obtain a compact decomposition, wavelet basis should provide a sparser representation of the source signal, which requires the waveform of mother wavelet to be similar to the

(i) Shannon Entropy values E(i,j) < E(i+1,2j) + E(i+1,2j+1) (ii) Quasi-IMF index QI(i, j) ≥ 1 (iii) Non-Gaussian distributed Otherwise, the node (i, j) is marked as a terminal leaf node. A flowchart of identifying the leaf nodes as terminal and representing the generation procedure of the decomposition tree is shown in Figure 2. No further evaluation process is needed on the “Discarded nodes”, as starred in the flowchart, since those leaf nodes have been confirmed to be Gaussian distributed noise. The maximum decomposition level should be determined by the predefined detection precision requirement. The decomposition level n corresponds to a frequency resolution of 2−n Hz (sampling time is 1 s by default). It is recommended that the maximum level is 5 or 6 for industrial oscillation detection. Oscillation Characterization with Cross-Checking Rule. Oscillation period estimation is essential to the root-cause diagnosis of the oscillation. Generally, an oscillation is considered to be regular if the standard deviation of the period is less than one-third of its mean value. The regularity statistic is defined as6 r= 15689

Tp̅ 3σTp

(10) dx.doi.org/10.1021/ie502129c | Ind. Eng. Chem. Res. 2014, 53, 15686−15697


Article

S1.4 Calculate the Shannon Entropy values of all the split nodes and their child nodes in step 2.3 to determine whether the split results in a decrease of the information cost. For a node at (i, j), if E(i,j) < E(i+1,2j) + E(i+1,2j+1), then label their parent node (i, j) as a terminal leaf and prune the child nodes (i, 2j) and (i, 2j + 1). Otherwise, preserve the child nodes. S1.5 Repeat steps 1.2−1.4 until all the leaf nodes are labeled terminal or the predefined decomposition level lmax has been reached. The best decomposition of WPT for oscillation detection is obtained as a pruned subtree T of Tfull. Step 2. Oscillation Characterization with Cross-Checking S2.1 Reconstruct the signals at all the leaf nodes of the optimal sub-band tree structure, which are non-Gaussian distributed and denoted as xij. S2.2 Evaluate the zero-crossing regularity of ACF of xij, denoted as rij. If r0 < 1, then ? ij does not contains an oscillation (rule 1). S2.3 Repeat S2.2 with a threshold denoising parameter THR at the same node (i, j). Obtain the reconstructed signal x̃ij from the sliced wavelet coefficients. S2.4 Evaluate the ACF regularity factor r̃ij of x̃ij and check if xij is an oscillation component with the cross-checking rule (rule 2). S2.5 Cluster the oscillations with similar period via a classification procedure by Chatfield.19 3.4. Threshold Selection for Cross-Checking. Note that a certain WP tree structure is a prerequisite to the calculation and application of the threshold. Birgé and Massart20 provide an algorithm to obtain THR by a wavelet packet coefficient selection rule using a penalization method. THR minimizes the penalized criterion and lets t* be the minimizer of M(t),

After a signal is decomposed via WPT, oscillations with different frequencies are distributed in separated bands wavelet domain. The challenge is that the industrial loops are usually contaminated by noise. Although a Gaussianity test on nodes will exclude some of the Gaussian-distributed subspaces and greatly reduced the scale of T, some subsignals will still contain noises. This is often the case when a noise segment is decomposed from an oscillation-contained subspace which usually locates in a narrow band and behaves non-Gaussian but irregular. The wavelet coefficients of a colored noise are random sequences, which also act like noise. In order to reduce the noise effect, the adaptive global threshold (THR), which will be elaborated in section 3.4, is introduced to reveal the most significant components in wavelet coefficients of the leaf nodes. The sliced wavelet coefficients c̃ij are obtained from cij by softthresholding as ⎧ c ji(k) ⎪ c ji(k) − THR i , |c ji(k)| ≥ THR ⎪ |c j(k)| c jĩ (k) = ⎨ ⎪ ⎪ 0, |c ji(k)| < THR ⎩

(11)

In other words, c̃ij preserves the most significant components of cij. The benefits of an adaptive global threshold is the accurate estimation of the noise variance can be obtained from the Gaussian distributed subspaces. It removes the noise segment at the same level from all leaf nodes. Apply the reconstruction process to cij and c̃ij according to (3), then xij and its counterpart x̃ij are obtained, respectively. Evaluate the corresponding ACF regularity factors via condition 10. And, the ACF evaluation rules are set up: Rule 1−Basic Rule: If the ACF of xij does not satisfy regularity condition (10), then xij is not an oscillation. But, it may be the nonstationary trend or simply the noise component. Rule 2−Cross-Checking Rule: If the ACFs of both xij and x̃ij satisfy the regularity condition 10, then conclude that xij is an oscillation component; If xij satisfies (10) but the counterpart x̃ij does not, then xij is not an oscillation but probably a bandpassed segment of the noise. The oscillation component is characterized as period Tp and regularity r in this step. Evaluate all the leaf nodes of the optimized wavelet packets tree, and all the oscillatory components of the input data will be identified with ACF evaluation rules. An illustrative example of cross-checking rule is presented in section 4 (example 2). 3.3. Steps of the WPT-Based Method. This section presents the proposed automated multiple oscillation detection method based on WPT. The method contains the following steps: Step 1. Improved Wavelet Packets Decomposition S1.1 Split the root node (0, 0) into binary child nodes (1, 0) and (1, 1) by WP analysis. S1.2 Apply a non-Gaussianity test to the current leaf nodes (without terminal label) and evaluate their quasi-IMF indexes. S1.3 Split a leaf node (i, j) into binary child nodes (i + 1, 2j) and (i + 1, 2j + 1), if (i) QI(i, j) ≥ 1 (ii) TG(i, j) statistic shows non-Gaussianity Otherwise, mark the leaf node (i, j) as terminal. During the iteration of the improved WPs, a leaf node in the WPT structure will not be split further if it is labeled as a terminal.

t * = argmin{M(t )} ⎛ ⎛ n ⎞⎞ M(t ) = −∑ c*j i(k) + ασ 2t ⎜2 + log⎜ ⎟⎟ ⎝ t ⎠⎠ ⎝ k≤t

(12)

where c*j i(k) are the wavelet packets coefficients sorted in decreasing order of their absolute value, n is the number of coefficients, and σ is the standard deviation (STD) of the noise segment, then THR = |c*j i(t*)|. α is a trade-off factor. The sparsity of the wavelet packets representation of the sliced signal grows with α, which means with larger α, more significant oscillation which acts more significantly will be preserved. Larger α applies to industrial scenarios with more significant signal-to-noise ratio. Let σij be the STD of the wavelet coefficients of subspace ? ij , and σ in the penalized criterion 12 is estimated as σ = mini,j{σij}. The improved wavelet packet decomposition structure, which identifies the Gaussian distributed subspaces during the decomposition, ensures the optimal estimation of the noise variance. It is preferred to choose α = 5 as default.

4. SIMULATION EXAMPLES In this section, four simulation examples are presented to illustrate the improved WPT-based oscillation detection method. The sampling periods are all selected as 1 s in the following examples. Example 1. This example illustrates the necessity of inducing the quasi-IMF index and the non-Gaussianity test to reduce the scale of the wavelet packet tree structure. A simple oscillating signal with noise corruption is generated as 15690



Article

Figure 3. Oscillating signals in time series for simulation study.

Figure 4. Reconstructed subsignals from original coefficients in example 2.

y1(k) = 2.5sin(0.1k) + v(k)

v(k) is a zero-mean white noise sequence with σv2 = 1.0. The original signal is shown in Figure 3-1. The improved WPT decomposed y1 into x30, x31, x21, and x11. Dominant oscillation is detected in y1 with a period around 62.9 s, as presented in Table 2. The quasi-IMF index successfully ceased the iteration of WP decomposition at the subspace ?30 and ?12 . ?11 did not split further because of its Gaussianity. It should be noticed that the maximum decomposition level is set at 5, but the level of improved WPT is only 3. Without the QIMFI and nonGaussianity test, a Shannon Entropy criterion will lead to a further decomposition, which deteriorates the regularity of reconstructed signal and increases the scale of wavelet packets tree structure. The decomposition tree with the Shannon Entropy test alone has 13 leaf nodes and 11 internal nodes. However, QIMFI and non-Gaussianity tests significantly reduce the tree structure into four leaf nodes(3, 0), (3, 1), (2, 1), and (1, 1)and two internal nodes. The results given in Table 1 also imply that the noise segment may be distorted due to the wavelet filters. The original signals show some regularity in zero-crossing distribution of their ACFs. But the existence of oscillations is denied in those narrow frequency bands according to crosschecking rules. Example 2. This example illustrates the effectiveness of the cross-checking rule in detecting multiple oscillations. Oscillating signal with two dominant frequencies is constructed as

Figure 5. Reconstructed subsignals from sliced coefficients in example 2.

oscillation components. The cross-checking rule with an adaptive threshold selection (α = 5) rejected the possibility of some oscillation existence, since the wavelet coefficients at these nodes are not significant compared with the noise estimation. r̃ = 0 implies the counterpart signal x̃ is almost zero in time series, because most of the wavelet coefficients are sliced by the soft-threshold. The sliced coefficients affect the noise-dominated subspace only, and four subsignals are concluded with no oscillation, namely x̃11, x̃21, x̃31, and x̃52. Their rotating behaviors are significantly reduced while the other three components, x̃50, x̃51, and x̃54 remain almost unchanged with

y2 (k) = 2.0sin(0.07k) + 1.5sin(0.2k) + v(k)

The signal is shown in Figure 3-2. There are two dominant oscillatory components in the original signal, and the improved WPT indicates at first that six subspaces probably contained Table 1. Results of the Proposed Method for Example 1 original

sliced

subspace

period

regularity

period

regularity

?30 ?13 ?12 ?11

62.9 ± 1.01

20.7

62.9 ± 1.01

20.7

conclusion

9.8 ± 2.9

1.1

0

not oscillation

4.2 ± 0.9

1.6

0

not oscillation

oscillation

Gaussian 15691



Article

Table 2. Results of the Proposed Method for Example 2 original

sliced

subspace

period

regularity

period

regularity

conclusion

?50 ?15 ?52 ?53 ?13 ?12 ?11

92.5 ± 10.6

2.9

92.5 ± 10.6

2.9

oscillation 1

37.3 ± 9.6

1.3

37.0 ± 10.3

1.2

oscillation 2

18.9 ± 3.7

1.7

31.3 ± 1.47

7.1

10.8 ± 0.97 5.6 ± 1.4

0.4

not oscillation

6.8

oscillation 2

3.7

0.6

not oscillation

1.3

0

not oscillation

31.3 ± 1.53

Gaussian

Table 3. Results of the Proposed Method for Example 3 original

sliced

subspace

period

regularity

period

regularity

?30 ?13 ?12 ?11

62.9 ± 1.01

20.7

62.9 ± 1.01

20.7

conclusion

10.2 ± 3.1

1.1

0.5

not oscillation

5.6 ± 1.6

1.2

0

not oscillation

oscillation

Gaussian

Figure 6. Wavelet packet tree structure of example 4.

Example 3. This example verified the proposed detection method with damping oscillation−time-varying amplitude. It illustrated the applicability of the improved WPT for one of the common time-varying features of the oscillation. A damping oscillation is constructed as

some minor effects on oscillation period estimation. The reconstructed signal from the original tree structure and the sliced one are presented in Figures 4 and 5, respectively. The result is shown in Table 2. Two dominant oscillations are extracted from the original signal and confirmed via ACF evaluation with the cross-checking rule. Oscillations x51 and x53 should be clustered together due to a classification procedure. Their scaled distance (|T5pl − T5p3|/max(σ51,σ53)) < 1, which indicate a dominant oscillation with a period around 31.3 s (more regular oscillation is recorded). Another dominant oscillation is detected in ?50 with period around 92.5 s.

y3 (t ) = 5e−0.0012t sin(0.1t ) + v(t )

It is shown in Figure 3-2. A similar result with example 1 was obtained and listed in Table 3. Only one dominant oscillation exists. Conventional techniques based on Fourier transform ignore the time-localization of the oscillation. Therefore, an inspection on the PSD of the original plot may not result in a 15692



Article

Table 4. Results of the Proposed Method for Example 4 original subspace

?50 ?15 ?52 ?53 ? 24 ? 34 ?12 ?11

sliced

period

regularity

period

regularity

50.5 ± 7.7

2.2

47.8 ± 13.3

1.2

16.3 ± 0.97

5.6

16.3 ± 1.01

5.4

oscillation 2

26.4 ± 3.52

2.5

0.4

not oscillation

0.0

not oscillation

2.8

oscillation 2

0.7

9.3 ± 1.48

2.1

15.7 ± 0.87

6.0

conclusion nonstationary trend

15.4 ± 1.83

oscillation 1

Gaussian Gaussian

significant peak value, which renders the oscillation presence to be ignored. Since the WPT analysis is scaled in both the time and frequency domain, the damped oscillations will be extracted through wavelet packet decomposition without information loss. The damping feature of the oscillation is preserved in one of the subspaces decomposed by improved WPT. The dominant oscillation has a period of 62.9 s. Example 4. The last example focuses on the most common time-variant feature of oscillations in industrial control loops nonconstant mean. A time series with two dominant oscillations and nonstationary trend w(k) is generated as y4 (k) = 2.0sin(0.07k) + 1.3sin(0.4k) + w(k) w(k) =

1 v (k ) 1 − 0.99z −1

Figure 7. Pressure control loop data PC3 in time series.

The nonstationary trend simulates the response of a loop which is excited by strong noise with a closed-loop pole on the stability boundary. The time series is shown in Figure 3-4. The structure obtained by the improved WPT method is given in Figure 6. Every node is marked with its entropy value and every terminal leaf node is annotated with the corresponding criteria for stopping further decomposition. The cross-checking rule excluded the possibility of oscillation existence at ? 24 and ?53 and confirmed that the other three subspaces contained oscillations. The result of ACF evaluation is presented in Table 4. There are two dominant oscillations in y4. One locates around 50.5 s. Oscillations x51 and x43 should be clustered as one, whose period is around 15.7 s. The nonstationary trend is usually localized at the lowest frequency band.

In the process of wavelet packet decomposition, the depth of wavelet decomposition is set to 5 and the db4 wavelet is selected. The wavelet packet decomposition of the data PC3 is shown in Figure 8. A non-Gaussainity test confirms that subspace ?11 is Gaussian distributed. Quasi-IMF indexes (QIs) indicate that ?13, ? 64 and ? 74 have the pattern of single oscillation. ?32 is marked terminal due to the criterion of Shannon Entropy. The result of the proposed method via ACF evaluation is presented in Table 5. It is noticed that x̃41 and x̃31 are adjacent in oscillation period, and calculation of their scaled distance, (|T̃ 4p1 − T̃ 3p1|/max(σ̃41, σ̃31)) < 1, indicates that they should be clustered as one oscillation. The reconstructed dominant oscillation along with its ACF is presented in Figure 9. Due to a set-point change of the pressure in PC3, the data has a nonconstant mean over the time domain. Wavelet packet transforms extracted the nonstationary trend from the raw data, as ?50 , which is shown in Figure 10. Case 2. This industrial case is presented to illustrate the applicability of the improved WPT to the multioscillatory data with a slow-varying trend. The data is provided by Horch,21 as presented in Figure 11, which is recorded from a flow loop FIC172 with proportional− integral−derivative (PID) control. The dominant oscillation is probably induced by couplings. There is another oscillation frequency with smaller amplitude. This case illustrates the benefit of a cross-checking rule for elimination of noise, while preserving the weak oscillation component. The wavelet packet decomposition is carried out with a depth of 6 and coif5 wavelet. The improved WPT structure is

5. APPLICATION TO INDUSTRIAL DATA SET Some benchmark industrial cases are presented to illustrate the improved WPT method. Analysis results indicate its applicability to control loops with nonstationary trend and multioscillations. Case 1. This benchmark industrial case illustrates the benefit of wavelet transform to deal with nonstationary data with a strong and sharp increase. The most common type of nonconstant mean of control loop data is the set-point change, which may demonstrate a sharp increase/decrease in the process variable. The benchmark data provided by Thornhill21 is a pressure measurement of control loop “PC3” in a refinery. The data record includes a possible set-point change introduced at t = 800 samples, whose oscillating frequency is high as shown in Figure 7. 15693



Article

Figure 8. Improved wavelet packet tree structure of PC3.

Table 5. Cross-checking Results of Subspaces in PC3 original subspace

period

sliced regularity

period

regularity

conclusion

?50

0

nonstationary trend

?15

1.0

not oscillation

?14 ?13 ?32 ? 64 ? 74 ?11

17.2 ± 3.0

1.9

18.2 ± 5.5

1.1

14.5 ± 2.8

1.7

14.1 ± 3.6

1.3

oscillation

oscillation

4.4 ± 0.86

1.7

0

not oscillation

5.0 ± 0.98

1.7

0

not oscillation

6.7 ± 1.2

1.9

0.8

not oscillation Gaussian

Figure 10. PC3 data and its nonconstant mean in time series. Figure 9. Dominant oscillation in PC3 data and its ACF.

presented in Figure 12. Since the process noise usually locates in a high frequency band, subspaces ?11, ?12 , and ?13 are confirmed to be Gaussian distributed, and it significantly reduced the computation effort of the evaluation process. Results of ACF cross-checking are presented in Table 6, with 6 ?1 and ?53 identified to be oscillation components of FIC172. The reconstructed oscillations along with their ACFs are

presented in Figure 13. The slow-varying trend in FIC217 is extracted by the improved WPT automatically in ? 60. The nonstationary trend x60 is presented in Figure 14. Case 3. This case study illustrates the availability of the proposed method in classifying nonoscillating behaviors. A nonoscillating signal recorded from a level control loop LC101 is presented in Figure 15. 15694



Article

the simulation examples detailed in section 4. All methods are carried out under default parameters without prior process knowledge. The IAE-based technique only successfully reports the single oscillation without description of the oscillation pattern. The results of ACF decay ratio and direct ACF zero-crossing analysis are highly dependent on the selection of band-pass filters to separate different oscillation frequencies. To implement the earlier modified version of EMD9 (EMD1) in oscillation detection, there is a basic step to remove the nonconstant mean from the signal. This method is effective to capture the most dominant oscillation only when a threshold is properly selected. When it deals with multiple oscillations, it requires different thresholds to recognize different oscillations by repeating this method on residuals. The more advanced EMD10 (EMD 2) for oscillation detection relies on the identification of key oscillatory IMFs of signal, which is applicable to mean-nonstationary data with multiple oscillation. However, because the decomposition structure is the same, there are many unnecessary IMFs included via corresponding indexes in the second step. EMD makes more efforts to decompose the nonconstant mean into several IMFs rather than giving a more detailed frequency resolution near the oscillating frequency. On the contrary, it is a compromise to use improved WPT to balance between number of subsignals and frequency resolution of the decomposition. In addition, the interpolations of EMD is more computationally costly than improved WPT, which enables a fast wavelet transform for industrial applications with limited processing capabilities. A typical execution time of the Matlab script for a time series with 1000 samples on a Intel Core i5 Processor is 0.43 s. DCT8 is able to detect oscillations in the simulation examples correctly as well. But the decomposition under the basis of discrete cosines is more redundant than needed. It is necessary to select the most significant DCT coefficients with

Figure 11. Flow control loop data FIC172 in time series.

The parameters for wavelet packet decomposition are chosen as 5 levels and db5 wavelet. The improved WPT structure is presented in Figure 16. It is exactly the same with DWT structure because the additional criteria proposed in this article successfully ceased the iteration on all of the detailed subspaces. It greatly reduced the computational cost to implement the proposed technique in real industrial applications. The cross-checking results are presented in Table 7. Most of the detailed subspaces are confirmed to be Gaussian and the only exception is ?14 . The decomposition at ?14 is stopped because of another criterionQIMFI, which indicates that x41 shapes like an potential oscillation. Cross-checking rule is applied to both x41 and x̃41 to confirm whether both ACF zerocrossings are regular enough to be oscillatory. Conclusion can be made that there is no oscillation in this signal.

6. COMPARATIVE STUDY A brief comparative study of the proposed oscillation detection method with some other methods is illustrated in Table 8 with

Figure 12. Improved wavelet packet tree structure of FIC172. 15695



Article

Table 6. Cross-Checking Results of Subspaces in FIC172 original subspace

? 60 ?16 ?15 ?52 ?53 ?13 ?12 ?11

sliced

period

regularity

54.0 ± 6.7

2.7

50.0 ± 13.9

1.2

period

regularity

53.0 ± 7.4

2.4

oscillation 1

0.9

not oscillation

0

17.1 ± 1.8

3.2

28.3 ± 3.5

2.7

conclusion nonstationary trend

21.1 ± 5.4

0.4

not oscillation

1.3

oscillation 2 Gaussian Gaussian Gaussian

Figure 15. Level control loop data LC101 in time series. Figure 13. Oscillations in FIC172 data and corresponding ACFs.

Figure 14. FIC172 data and its slow-varying trend in time series. Figure 16. Improved wavelet packet tree structure of LC101.

two noise level parameters, which are separately obtained from each set of coefficients. There is a similar idea in this work as the cross-checking rule when dealing with the noisy segments. It is highlighted that in this article, only one threshold is required under the optimal wavelet packet structure with a global estimation of the noise rather than denoising subsignals individually. By adjusting the trade-off factor α, it is easy for the user to handle how sparse the representation of the denoised

signal should be. In this way, improved WPT is able to provide either compact or detailed representation with tuning α. In DCT-based methods, there is no better solution for obtaining a more compact decomposition with the original algorithm of DCT, while extracting all the interested frequency bands with potential oscillation existence. This dilemma is solved in the 15696



■

ACKNOWLEDGMENTS The authors wish to thank sponsor and financial support from Natural Science Foundation of China under grant No. 61134007, No. 60904039, and the 111 Project under grant No. B07031.

Table 7. Cross-Checking Results of Subspaces in LC101 original subspace

period

sliced

regularity

? 50

period

regularity

0.7

conclusion

?15

56.38 ± 8.0

2.3

0.9

nonstationary trend not oscillation

?14 ?13 ?12 ?11

21.2 ± 6.9

1.0

0.3

not oscillation

■

Gaussian Gaussian

Table 8. Result Comparison of the Proposed Oscillation Detection Method with Existing Work method

signal 1

signal 2

signal 3

signal 4

yes yes no no yes yes yes yes

yes yes yes yes no yes yes yes

yes no yes yes yes yes yes yes

no no no no no yes yes yes

improved WPT with a mechanism of generating representation by three criteria. The theoretical maximum number of subsignals decomposed from improved WPT is only 16 at the fifth level.

7. CONCLUDING SUMMARY Based on an improved wavelet packet transform, a new automatic approach for oscillation detection is proposed. In contrast to existing oscillation detection approaches, (i) improved WPT represents an optimal decomposition tree structure for oscillation detection, which provides the possibility to be robust to the time-varying data features. (ii) The improved WPT is combined with entropy-based base selection and another two oscillation-detection-oriented sifting criteria, namely the non-Gaussianity test and quasi-IMF index. It reduces the computational cost of the decomposition, while the frequency resolution is improved where oscillations may occur. The isolated subspaces is the most sufficient and effective decomposition for oscillation detection. In this way, all the subspaces where oscillations potentially exist will be sifted out and evaluated via basic and cross-checking rules. (iii) The improved WPT method for oscillation detection can be automatically implemented without human intervention, since the sifting criteria determined the scale of the tree and the threshold selection is adaptive. In addition, there is no information loss during the decomposition procedure. Simulation examples and the implementation on industrial data verified the applicability of the proposed technique. These applications under common oscillation scenarios illustrate the advantages of the proposed approach over traditional techniques.

■

REFERENCES

(1) Jelali, M., Huang, B. Detection and diagnosis of stiction in control loops: state of the art and advanced methods; Springer: London, 2010. (2) Tangirala, A. K.; Kanodia, J.; Shah, S. L. Ind. Eng. Chem. Res. 2007, 46, 801−817. (3) Hägglund, T. Control Engineering Practice 1995, 3, 1543−1551. (4) Forsman, K.; Stattin, A. A new criterion for detecting oscillations in control loops. European Control Conference, Karlsruhe, Germany, Aug 31−Sep 3, 1999. (5) Miao, T.; Seborg, D. E. Automatic detection of excessively oscillatory feedback control loops. IEEE International Conference on Control Applications, Kohala Coast, HI, Aug 22−27, 1999. (6) Thornhill, N.; Huang, B.; Zhang, H. Journal of Process Control 2003, 13, 91−100. (7) Matsuo, T.; Sasaoka, H.; Yamashita, Y. Detection and diagnosis of oscillations in process plants. Knowledge-Based Intelligent Information and Engineering Systems; Lecture Notes in Computer Science; Springer, 2003; pp 1258−1264. (8) Li, X.; Wang, J.; Huang, B.; Lu, S. Journal of Process Control 2010, 20, 609−617. (9) Srinivasan, R.; Rengaswamy, R.; Miller, R. Control Engineering Practice 2007, 15, 1135−1148. (10) Srinivasan, B.; Rengaswamy, R. Control Engineering Practice 2012, 20, 733−746. (11) Zakharov, A.; Jämsä-Jounela, S.-L. Ind. Eng. Chem. Res. 2014, 53, 5973−5981. (12) Thornhill, N.; Hägglund, T. Control Engineering Practice 1997, 5, 1343−1354. (13) Learned, R. E.; Willsky, A. S. Applied and Computational Harmonic Analysis 1995, 2, 265−278. (14) Daubechies, I. Ten lectures on wavelets; SIAM, 1992; Vol. 61. (15) Shannon, C. E. ACM SIGMOBILE Mobile Computing and Communications Review 2001, 5, 3−55. (16) Liu, X.; Xie, L.; Kruger, U.; Littler, T.; Wang, S. AIChE J. 2008, 54, 2379−2391. (17) Hyvarinen, A. Neural Networks, IEEE Transactions on 1999, 10, 626−634. (18) Shinde, A.; Hou, Z. Structural Health Monitoring 2005, 4, 153− 170. (19) Chatfield, C.; Collins, A. Introduction to Multivariate Statistics; Chapman and Hall/CRC, 1980. (20) Birgé, L.; Massart, P. From model selection to adaptive estimation; Springer, 1997. (21) Thornhill, N.; Horch, A., et al. Appendix A. Detection and diagnosis of stiction in control loops: state of the art and advanced methods; Jelali, M., Huang, B., Eds.; Springer: London, 2010.

Gaussian

IAE by Hagglund IAE by Forsman decay ratio of ACF direct ACF EMD 1 EMD 2 DCT improved WPT

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. 15697


Automatic Detection of Nonstationary Multiple Oscillations by an

Recommend Documents