Wavelet Transform Based Methodology for Detection and

Article pubs.acs.org/IECR

Wavelet Transform Based Methodology for Detection and Characterization of Multiple Oscillations in Nonstationary Variables Elham Naghoosi and Biao Huang* Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta T6G 2G6, Canada ABSTRACT: Diagnosing the root cause of a propagated oscillation in the operation requires detection of all process variables that are oscillating with similar frequencies followed by application of an appropriate root cause diagnosis procedure. Oscillations in chemical processes are usually caused by controller tuning, valve problems, or external oscillatory disturbances. There are several methods proposed in literature for root cause diagnosis of oscillations within the system. However, most of the methodologies can only work for a specific type of oscillation. For example, the methodologies based on quantifying the nonlinearity of variables can help with root cause diagnosis of a valve-induced oscillation but cannot help if the oscillation actually has occurred due to aggressive controller tuning or due to an external oscillatory disturbance. Therefore, before trying to find out which loop within the system has caused the oscillation, it is important to categorize the oscillation meaning to learn if the oscillation is caused by a nonlinear valve within the system, controller tuning or an external disturbance to choose an appropriate diagnostic procedure. The different characteristics of these three oscillation types are studied in the literature with methodologies to distinguish them from each other. However, the proposed methodologies can work reliably when there is only one oscillatory component present in the variables and cannot help in cases of multiple oscillations. Also, nonstationary trends and noise within variables are yet a challenging issue in detection and diagnosis of oscillations. This paper presents a comprehensive oscillation detection and characterization procedure based on wavelet transform. The methodology is capable of both detection and independent characterization of multiple oscillation frequencies in variables as well as implementing automatic noise and nonstationary trend removal algorithms. Advantages of the proposed method are illustrated through analysis of data sampled from an industrial process as well as simulations.

1. INTRODUCTION Oscillations are common faults in chemical processes mostly due to nonlinearity in the valves, aggressive controller tuning, or external disturbances. Oscillations are detrimental to operation of the process, and thus their detection and diagnosis is an important task. Root cause diagnosis of an oscillatory fault propagated through several loops within the system requires a few steps. First, all the variables affected by the oscillation should be detected because the root cause could be any one of the variables oscillating with similar frequencies. After detection, an appropriate diagnostic procedure is required to find the loop and the element within the system that is causing the oscillation among all the possibilities. There are several algorithms proposed in the literature for oscillation detection and diagnosis.1,2 However, there is one step yet missing. Even though there are several methods proposed to diagnose the root cause of an oscillation, most of them can work reliably only for a specific type of oscillation. For example, methodologies based on nonlinearity ranking of variables can work only if the oscillation is caused due to nonlinearity within the process and cannot help if the oscillation is actually due to controller tuning or an external oscillatory disturbance. In the same way, if the oscillation is caused due to an aggressively tuned controller, it is possible to use causality analysis procedures such © XXXX American Chemical Society

as Granger causality to diagnose its root cause because the nature of the oscillation allows us to model it as the response of a linear model. However, application of Granger causality may not lead to reliable results if the oscillation is due to valve stiction, because nonlinear behavior cannot be captured through a linear model. External oscillatory disturbance if it has a deterministic nature requires a different diagnostic algorithm, which will not be discussed in this paper. Thus, an important step before root cause diagnosis is to learn the category of the oscillatory fault to be able to select a proper diagnostic procedure. Learning the oscillation category, which is called oscillation characterization in this paper, means to learn if the oscillation is likely to be due to controller tuning, valve nonlinearity, or an external disturbance. This task is possible due to the different properties of oscillations caused by the mentioned three different sources. Oscillatory faults more or less maintain their properties when they propagate through the process. Therefore, this paper proposes to check the properties of a specific oscillation Received: Revised: Accepted: Published: A

August 23, 2016 November 7, 2016 February 1, 2017 February 1, 2017 DOI: 10.1021/acs.iecr.6b03075 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

of using wavelet bicoherence over classical bicoherence is the possibility to study the changes in data properties with time. Therefore, the proposed methodology can potentially be applied to study the structural changes in variables with time, which is valuable when dealing with nonstationary data. The paper is organized as follows. First, a preliminary introduction to wavelet transform is given. The wavelet power spectrum (wps) and its properties are reviewed in section 3 followed by development of the methodology for diagnosis of oscillations between controller tuning and external disturbances based on its properties. Section 4 reviews wavelet bicoherence and proposes a hypothesis test for detection of nonlinearity in oscillations. Section 5 presents a case study, and section 6 draws the conclusions.

frequency in all the affected variables and, after learning its category, choose the proper diagnostic algorithm. Choudhury et al.3 have discussed extensively the properties of valve stiction-induced oscillations and proposed to use bicoherence as a nonlinearity measure for their diagnosis. Naghoosi et al.4 discussed the unique properties of controllerinduced oscillations and the methodology to distinguish them from sinusoidal disturbances and valve-induced oscillations. The mentioned methodologies can work reliably for stationary variables containing a single oscillatory component. This paper transfers these concepts to wavelet domain to achieve a comprehensive algorithm capable of both detection and independent characterization of multiple oscillations hidden in the usual nonstationary trends and noises of industrial variables. Wavelet transform provides a time dependent, multiresolution decomposition of signals that facilitates handling nonstationary trends, analyzing intermittent oscillations, and dealing with the presence of multiple oscillations. It also facilitates individual characterization of various oscillatory components within one variable due to its intrinsic nature in decomposing the variables. Wavelet transform is extensively applied in signal and image processing especially for the purpose of noise filtering and data compression.5 It is also applied for fault detection and diagnosis in engineering applications. Peng et al.6 reviewed the application of wavelet transform in monitoring machine conditions. Yang et al.7 proposed a fault detection algorithm based on the use of wavelet transform for noise filtering and preprocessing of data followed by kernel entropy component analysis. Alabi et al.8 also integrated wavelet transform with a generic dissimilarity measure to achieve a better performance monitoring tool. Zhu et al.9 presents a review on the literature regarding monitoring tool conditions based on wavelet transform. Bakshi10 integrated wavelet transform with PCA to monitor the slow and fast features of the data individually along with better noise removal. Lin et al.11 proposed a wavelet denoising method to extract the weak features in variables for fault diagnosis purposes. Wavelet transform has also been utilized for the purpose of oscillation detection. Matsu et al.12 proposed visual inspection of wavelet transform of the variables as a method to detect oscillations in noisy environment. Guo et al.13 used wavelet packet transform for detection of multiple nonstationary oscillations in process variables. Plett14 proposed a hypothesis test based on using cross wavelet transform and coherence to detect transient oscillations in variables. Several other methodologies have also been proposed in the literature for the purpose of oscillation detection in the presence of multiple oscillations and nonstationary trends such as refs 15−17. A simple, easily implementable and yet effective oscillation detection methodology based on wavelet transform is proposed here that does not require the implementation of a separate oscillation detection algorithm. However, the main novelty is to facilitate the independent characterization of different oscillation frequencies within one variable. It should be mentioned that the idea of using wavelet transform for diagnosis of oscillations between controller tuning and external disturbances was very briefly introduced in ref 18, which is elaborated on in this work. There are two hypothesis tests proposed in this paper. One is based on examining the wavelet power spectrum to distinguish between controller-induced oscillations from sinusoidal external disturbances. The other is based on examining wavelet bicoherence to detect nonlinearity in oscillations. An advantage

2. REVIEW ON WAVELET TRANSFORM The continuous wavelet transform (CWT) of a signal x(t) at time b and scale a (wx(b,a)) is defined as +∞

wx(b ,a) =

∫−∞

x(t ) Ψ*a , b(t ) dt

(1)

The wavelet function Ψa,b is a scaled and translated version of mother wavelet function Ψ as 1 ⎛t − b⎞ ⎟ Ψ⎜ a ⎝ a ⎠

Ψa , b =

(2)

The mother wavelet function Ψ is usually chosen as a decaying oscillatory function. Changes in the scale a and time b provide three-dimensional representation of the signal in the time− frequency domain. The values of wavelet coefficients are a function of position in the time domain, scale, and the choice of mother wavelet function in addition to the signal itself. The relation between scale and frequency is approximately fa =

fc (3)

aΔ

where fc is the central frequency of the wavelet function, Δ is the sampling period, and a is the scale. The choice of the mother wavelet function can be done subjectively depending on which feature of the data is of interest to be analyzed. Because this work is mainly concerned with using the wavelet transform for oscillation diagnosis, the Morlet wavelet function is selected as defined in the following equation. Ψ(t ) = π −1/4 exp(iωct ) exp( −t 2/2)

(4)

The continuous wavelet transform (CWT) has difficulties in estimation due to the continuous change over both time and frequency domains and thus it is usually estimated at a finite number of time samples and scales. The discretized CWT of a signal is obtained as shown in eq 5. N−1

wx(b ,a) =

⎛ Δ ⎞1/2 ⎡ (b′ − b)δt ⎤ ⎟ Ψ*⎢ ⎥ ⎣ ⎦ a⎠ a

∑ xn′⎜⎝

n ′= 0

(5)

where a is the scaling factor and b is the translation factor. Changing a between 0 and N − 1 changes the width of the wavelet function and therefore provides different frequency resolution. b varies between 1 and N (the number of data samples), moving the wavelet function in the time horizon and providing a time dependent representation of the signal. B

DOI: 10.1021/acs.iecr.6b03075 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research As an example to illustrate wavelet transform, consider a variable constructed as ⎛ 2π ⎞ ⎛ 2π ⎞ x(t ) = sin⎜ t ⎟ + u(t − 250) sin⎜ t ⎟ ⎝ 20 ⎠ ⎝ 35 ⎠

where u(t) is the unit step function. The variable is plotted in the top panel of Figure 1, which shows that it consists of a persistent

Figure 2. Wavelet coefficients of x(t) for scales from 1 to 50.

wavelet function, which is a damped sinusoidal in the case of the Morlet wavelet function, in the same way that Fourier transform measures the similarity between the signal and sinusoidal waves of different frequencies. If the similarity is high, the absolute value of the coefficients will be high, which usually happens if the scale of the wavelet function is close to the period of the signal’s oscillation. Therefore, the frequencies corresponding to the scales in which the coefficients have larger absolute values can be taken as oscillation frequencies within the variable. To detect those scales, it is possible to take an average of absolute wavelet coefficients (or squared values) through time to obtain the map of wavelet power versus scale as in the following equation: 1 awps x(a) = ΣbN= 1 |wx(b ,a)|2 (6) N The curve of awps of the signal plotted in Figure 1 versus scale is shown in Figure 3. Local peaks in this curve correspond to the

Figure 1. Top: x(t). Middle: wavelet coefficietns at a scale approximately coresponding to oscillation period of 20. Bottom: wavelet coefficietns at a scale approximately coresponding to oscillation period of 35.

oscillation with a period of 20 and an additional oscillation with a period of 35 from sample 250 onward. The wavelet coefficients at the scale of 16, which is the closest scale to an oscillation period of 20, calculated using Morlet wavelet function, are plotted in the middle panel. The coefficients show a persistent oscillation during the whole time interval, as is expected. Wavelet coefficients at the scale of 28, approximately corresponding to an oscillation period of 35, are also plotted at the bottom panel. The coefficients are very close to 0 approximately up to the time sample of 250 when the oscillation actually appears in the data. The wavelet coefficients correctly indicate the presence of oscillation with a period of 35 for the second half of the whole period. The fluctuations at the beginning and the end of the time interval are due to a phenomenon explained as a cone of influence in the literature. Usually the wavelet coefficients are estimated for a range of scales and the absolute value of coefficients is color coded in a figure, as shown in Figure 2, which shows the wavelet transform of x(t) for scales from 1 to 50. The scales close to 16 and 28 where the signal actually has its highest energy have brighter colors, which indicates larger wavelet coefficients. The map of wavelet coefficients shows which scales or frequencies have the highest contribution to the signal and at which time samples. A simple method to detect oscillations within a variable based on wavelet transform is to detect the scales where the coefficients have larger values than the neighboring scales. CWT coefficients are actually the inner product of the signal with shifted and scaled versions of the mother wavelet function. Thus, wavelet coefficients quantify the similarity between the signal and the

Figure 3. Averaged wavelet power spectrum versus scale for signal plotted in Figure 1.

socillation periods within the variable. In practice, this curve should be plotted after removing the nonstationary trend and noise of the data, which is possible through implementation of a wavelet based filtering and trend removal algorithm.

3. DIAGNOSIS OF OSCILLATIONS BETWEEN CONTROLLER TUNING AND EXTERNAL DISTURBANCE BASED ON WAVELET POWER SPECTRUM Wavelet power spectrum (wps) is utilized in this work for the purpose of oscillation diagnosis. First, there is a quick review on the properties of estimated wps followed by the methodology to use it for the diagnosis purpose. C


Article

Industrial & Engineering Chemistry Research 3.1. Properties of Wavelet Power Spectrum. wps of a variable x(t) at time b and scale a is estimated as 2

wps x (b ,a) = |wx(b ,a)|

ρτ = pr

The confidence interval for wps of a normally distributed random variable is available in the literature. Zhang et al.19 mathematically proved that the distribution of the wavelet power spectrum of white noise obtained on the basis of Morlet wavelet function can be obtained as

pm

2

(8)

2

be 0 and the distribution simplifies to |wx(b,a)|2 ⇒ σ (X12 + X22) or 2

χ22

σ2 2 χ 2 2

(9)

where has χ distribution with two degrees of freedom. It should be noted that if the wavelet function is real with no imaginary part, the χ2 distribution will have one degree of freedom. Torrence et al.20 heuristically proposed that the probability distribution of the wavelet spectrum of the response of a general linear process to Gaussian innovations is

|wx(b ,a)|2 ∼

2

1 λk χ2 2 2

N−t cos(ωt ) N − tm

(13)

where pm is the peak value at time tm, which is in the middle of the wavelet coefficient’s time domain. The peak at the middle time is considered to avoid the uncertainty in the estimated CWT coefficients at both the beginning and end due to the cone of influence phenomena. In the case of oscillation due to aggressively tuned controller, the output of the loop can be modeled with an ARMA structure that has complex poles. A pair of complex poles generate the oscillation in the response of the system in a noisy environment. Considering a second-order AR process with imaginary poles, the response to noise can be shown as in eq 14 where amplitude (βt) and phase (ϕt) of the oscillation change as functions of time:

where X1 and X2 are independent standard Gaussian random variables and σ2 is the noise variance. If the wavelet central frequency is sufficiently large (ωc ≥ 6), the exponential terms will

|wx(b ,a)|2 ⇒

(12)

where N is the sample size and τr is the time lag corresponding to the peak chosen as the reference. When CWT of this ACF is estimated, because CWT at the scale corresponding to ω is just the convolution of the deterministic wavelet function with this deterministic ACF, it will also have a deterministic behavior, which can be written as

(7)

2 2 σ2 σ2 |wx(b ,a)| ⇒ (1 + e−ωc )X12 + (1 − e−ωc )X 2 2 2 2

N−τ cos(ωτ ) N − τr

yt = βt r t cos(tω + ϕt ) + εt

(14)

The ACF of this response will also be oscillatory. Again, in the case of a second-order AR process with complex poles, ACF can be written as in eq 15.

(10)

where λk is the power spectrum of the data at frequency k corresponding to scale s. A confidence interval for the true value of |Wx(b,a)|2 is derived on the basis of the estimated value and the distribution in eq 10, as shown in eq 11.20 2 2 |wx(b ,a)|2 ≤ |Wx(b ,a)|2 ≤ 2 |wx(b ,a)|2 2 χ2 (p /2) χ2 (1 − p /2)

ρj =

a 2 j /2(sin(j + 1)ω − a 2 sin(j − 1)ω) (1 + a 2) sin ω

(15)

where a2 is a model parameter. This equation shows that the ACF is oscillatory with the same frequency as the model response but with a decay ratio of a2j/2. Therefore, the unique characteristic of the ACF of the response of a system oscillating due to controller tuning is that the ACF converges to 0 after a few time lags due to the decay ratio. It should be mentioned that when this ACF is estimated from data, the oscillation amplitude will be random, similar to the model response. To illustrate the difference between the two mentioned oscillation types, Figure 4 plots the sampled output of a feedback

(11)

where p is the significance level and λk in eq 10 is replaced by the true wavelet spectrum value. Even though this confidence interval is utilized a lot in literature (the corresponding paper was highly cited), it is questionable. Because the expected value of wps of colored noise is not 0, the distribution cannot be central χ2. There is no proof that the distribution is χ2 in general, but even if assuming it is χ2, it has to be noncentral χ2. However, eq 11 can be used for Gaussian noise and its use in oscillation diagnosis is described in the following section. 3.2. Diagnosis Methodology. This section proposes a hypothesis test based on using wps for automatic diagnosis between oscillations due to controller tuning and external disturbances. The diagnosis is done on the basis of the different properties of the ACFs of oscillations caused by controller tuning compared with sinusoidal disturbances as described in ref 4. As a quick review, ACF of a harmonic with added white noise (xt = Σki=1Ai cos(ωit+ϕi) + εt) has a deterministic behavior theoretically described as ρτ = 12 Σki=10.5Ai2 cos(ωiτ) where σx is

2.9s 2 , which is controlled 34s 2 + 17s + 2 0.2 by a PID controller implemented as 3.3 + + 0.8 100100 . There is s 1+

loop (yt) with a process model as

s

(

also a sinusoidal disturbance as 2 sin

2π 0.0419

)

entering the loop

through the output. The ACF of the response along with CWT of the ACF is also plotted. From the ACF in the second panel, it can be seen that there is a high frequency oscillation present in the variable in addition to the sinusoidal wave. The awps curve of yt shown in Figure 5 indicates that two scales have higher energy within the spectrum, which indicates the presence of oscillations in the variable with periods close to the two scales. From the ACF itself it is not possible to correctly characterize both of these oscillations. The cwt of the ACF of y(t) at two scales corresponding to the two oscillations is plotted in Figure 6. As can be seen at the bottom panel of Figure 6, the CWT coefficients corresponding to the harmonic part of the signal show a deterministic oscillation amplitude that linearly decays with time as is expected from the

σx

the standard deviation of xt and τ is time lag. When ACF is estimated from data, it will have a linear decay with increasing time lag due to the numerical issues. However, the estimated ACF of a single sinusoidal function as xt = A cos(ωt+ϕ) + εt also has a deterministic behavior that can be determined by the knowledge of its frequency and the value of one of the peaks in the ACF (Pr) using the following equation. D


Article

Industrial & Engineering Chemistry Research N−τ

using ρτ = pr N − τ cos(ωτ ). This prediction will result in large r

errors if the oscillation is actually caused by controller tuning. The hypothesis test can be developed by verifying the ratio between the calculated wavelet power spectrum of the ACF and the predicted values. From eq 10 we can derive |w(b ,a)|2 1 ∼ χ12 2 |W (b ,a)|2

(16)

In eq 16, w(b,a) is taken as the calculated value of wavelet coefficients and W(b,a) is taken as the predicted value by eq 13 assuming that the oscillation is harmonic. This hypothesis test is equivalent to assuming a confidence interval for the ratio between the calculated wavelet spectrum and the predicted ones. The ratio corresponding to truly harmonic signals should be equal to 1 with a minor error because ⎡ ⎤ ⎢ ⎥ |wρ(b ,a)|2 ⎥=1 E⎢ 2 N−t ⎢ ⎥ ⎢⎣ pm N − tm cos(ωt ) ⎥⎦

Figure 4. Simulated process output (yt) with its ACF and CWT of its ACF.

(17)

If the ratio is significantly different from 1 at some samples, then it is inferred that the oscillation is caused by the system due to controller tuning. The significance level of the distribution determines the confidence interval. The top panel of Figure 7 plots this ratio corresponding to the controller-induced oscillation in the data plotted in Figure 6 and

Figure 5. awps of process output plotted in Figure 4 versus scale.

Figure 7. Ratio between calculated peak values of the CWT of the ACF of yt at two scales.

the bottom panel plots the ratio corresponding to the sinusoidal disturbance. Here only the peak values are analyzed. The ratio is estimated by taking all the peak values as the reference one by one and predicting the rest of the peaks. As can be seen in Figure 7, the ratio corresponding to the harmonic oscillation is close to 1 and within 95% confidence interval and the ratio corresponding to the controller tuning-induced oscillation is far from one for many samples. Thus, the proposed test can characterize each oscillatory component in the variable independently, which is not possible by other methodologies in literature so far. One point to mention here is regarding the possibility of a sinusoidal disturbance occurring in the system. There are different ways a sinusoidal disturbance can occur within the process. It could be due to high resolution of a sensor or electricity

Figure 6. CWT coefficients of ACF of y(t) at two scales.

ACF of a harmonic itself (the beginning and end should be disregarded due to the cone of influence phenomena). Also, the coefficients at the top panel corresponding to the oscillation due to the controller tuning show a decaying oscillation with varying amplitude as is expected from the respective ACF. The predictable behavior of CWT of the ACF at the scale corresponding to harmonic oscillation makes it possible to develop a hypothesis test to distinguish these two from each other. The null hypothesis is that the oscillation is a harmonic process and not due to controller tuning. On the basis of this assumption, CWT coefficients of the ACF can be predicted with negligible error by taking a nominal peak as the reference and E


Article


It should be noted that bicoherence has 12 symmetric regions in the f1, f 2 plane.21 The principal domain is defined as the region where f1 ≤ f 2 and f1 + f 2 ≤ fs. The principal domain contains all the information and the other regions are redundant. It is proved in the literature that even though a linear time series with non-Gaussian innovations will have non-zero bicoherence, the value of bicoherence will be constant and independent of bifrequencies.21 However, if there is nonlinearity in the structure of data generation process, bicoherence will have non-zero as well as nonconstant values. These facts have been used in the literature to propose methods for detection of nonGaussianity and nonlinearity in signals.21,24 4.1. Wavelet Bicoherence. Wavelet bicoherence25 is the generalization of the classical bicoherence. Similar to other quantities based on wavelet transform, wavelet bicoherence is specifically suitable for analysis of nonstationary signals containing temporary behaviors, pulses, and abrupt changes.26 Wavelet bicoherence has been applied to detect phase coupling and nonilinearites in variables. The squared wavelet bicoherence between two scales a1 and a2 is defined as

disturbance or the other case is when an oscillation due to valve stiction propagates far from its source. Nonlinear oscillations (oscillations caused due to the nonlinearity within the process) contain several high frequency harmonics that produce the nonlinear shape (ex. square or triangular oscillations). These harmonics get damped while going through the process, which usually acts as a low pass filter. Thus, valve-induced oscillations become more linear and similar to sinusoidal disturbances when getting farther from the source. However, a controller-induced oscillation will also keep its characteristics when traveling to other loops. As mentioned in the Introduction, when it is determined that a specific oscillation frequency shows the characters of controller-induced oscillation in all the affected variables, then proper diagnostic procedure can be selected to find which variable is actually the source of that oscillation.

4. BICOHERENCE AND ITS APPLICATION IN OSCILLATION DIAGNOSIS Bicoherence measures the interaction between two frequencies in one variable or between two variables. The interaction between two frequencies, which is also called phase coupling, is usually a sign of nonlinearity or existence of some kind of structure in the variable. Phase coupling occurs when oscillations with two frequencies, f1 and f 2, are present in the signal as well as the frequency f = f1 + f 2, while the phases of these frequency components also have the relation ϕf1 + ϕf 2 = ϕf + const. The value of the bicoherence estimated for two frequencies is close to 1 if these conditions are satisfied and is close to 0 otherwise. Bicoherence has been used to detect the presence of nonlinearity in the signal generation process.21 Because the oscillations due to valve stiction are nonlinear, it can be used to diagnose oscillations in control systems. To understand bicoherence, it should be mentioned that bicoherence is the normalized version of bispectrum, which is the frequency domain counterpart of the third-order cumulant C3. The third-order cumulant of variable x(t) with 0 mean value is defined as C3x(m ,n) = E[x(t ) x(t +m) x(t +n)]

wbicx(a1 ,a 2)2 =

|Σmk = 1wx(k ,a1) wx(k ,a 2) wx*(k ,a)|2 Σmk = 1 |wx(k ,a1) wx(k ,a 2)|2 Σmk = 1 |wx(k ,a)|2 (20)

where the summation is over m samples around the desired time sample. wbicx(a1,a2)2 has a value between 0 and 1. Squared bicoherence will be non-zero if the variable does not follow normal distribution. There are few papers in the literature proposing methods to estimate the significance level of squared wavelet bicoherence. All of the proposed methodologies are based on comparing the estimated wavelet bicoherence with the one that could be estimated from Gaussian noise due to numerical issues. Ge27 has mathematically derived the probability distribution function (PDF) of the wavelet bicoherence estimated from white Gaussian noise as

(18)

A linear variable x(t) (linear in the sense that it can be expanded as x(t) = Σ∞ s=−∞h(t−s) ε(s) where ε(t) are independent and identically distributed (i.i.d.) innovations), is normally distributed if the distribution of innovations ε(t) is normal. C3 of a linear time series can be non-zero only if the innovations ε(t) are not normal distributed and have non-zero third-order cumulant. When the innovations are normal distributed, x(t) may have non-zero C3 only if it is a nonlinear function of the innovations. Therefore, non-zero third-order cumulant indicates that the variable does not have a symmetric distribution, which could be due to the nonlinearity in signal generation process or the nonGaussian distribution of the innovations. Bicoherence is defined on the basis of Fourier transform as

fwbic2 (z) =

Ne(f1 ,f2 ,f )

2

2πz (1 − z)

e−(arctanh

z)

/2[Ne(f1 ,f2 ,f )] (21)

where Ne(f1 ,f2 ,f ) =

m ⎡ fc fs max⎢ ⎣

f12

+ f2 2

f1 f2

,

⎤

2 ⎥ f1 + f2 ⎦

ln d (22)

where d is the allowable correlation between wavelet coefficients. The significance at the level of α (Dα) can be obtained as ∞

bic(f1 ,f2 ) =

Dα =

E[X(f1 ) X(f2 ) X *(f1 + f2 )] 2

2

E[|X(f1 ) X(f2 )| ]E[|X(f1 + f2 )| ]

∫D

α

(19)

fwbic2 (z) dz = α

(23)

Squared bicoherence significantly different from 0 is a sign of nonGaussian distributed variable. Again, the nonconstant value of squared bicoherence can be taken as a sign of nonlinear structure in the variable. To observe the advantage of wavelet bicoherence over the bicoherence estimated on the basis of Fourier method, consider

where X(f) is the discrete Fourier transform of time series x(t) and * denotes the complex conjugate. Bicoherence itself is a complex quantity whereas squared bicoherence has a real value between 0 and 1. For more details on the methods on estimation of bicoherence, readers can refer to refs 22 and 23. F


Article

Industrial & Engineering Chemistry Research a variable as

Pr(|x−μ|≥kσ ) ≤

⎛ ⎛ π⎞ π ⎞⎟ y(t ) = sin⎜2π × 0.04 + ⎟ + sin⎜2π × 0.025 + ⎝ ⎝ 3⎠ 12 ⎠ ⎛ π⎞ + u(t − 1000) sin⎜2π × 0.065 − ⎟ ⎝ 4⎠

1 k2

(26)

In this case μ = wbic2(a≠aosc) and σ can be estimated as the standard deviation of wbic2(a1≠aosc,a2≠aosc). k is a tuning parameter that defines the significance level. If wbicmax2(aosc,a2) is significantly larger than wbic2(a≠aosc), it implies the oscillation has a nonlinear structure. In terms of oscillation diagnosis, if a feedback loop is oscillating and no nonlinearity is detected in the oscillation, the cause could be due to a tightly tuned controller, external disturbance, or sensor. The oscillations caused by valve problems such as stiction or backlash are naturally nonlinear. Therefore, bicoherence can diagnose if the fault is caused due to the nonlinear behavior of the valve or the process itself. The flowchart in Figure 9 shows the steps required for oscillation diagnosis. It should be noted that this work diagnoses

where u(t − 1000) is a step function with initial value of 1 and final value of 0. The step function eliminates the oscillation with the frequency of 0.065 from the signal from time sample 1000 onward. The bicoherence between scales 32 and 20 of this signal is plotted in Figure 8. From Figure 8 the change in the

Figure 8. Bicoherence of signal y(t) between scales 32 and 20 in time.

bicoherence at time sample equal to 1000 is easily observable. Therefore, the wavelet method is a much more suitable tool to detect the changes in signal’s behavior. Also, to obtain a general representation of the signal, total bicoherence based on averaging the bicoherence values over the time range can be used. 4.2. Diagnosis of Oscillations Due to Nonlinearity. The purpose of using bicoherence is to detect if the oscillation is nonlinear, which indicates the source of the oscillation to be valve problem. After detection of the oscillation period, it is required to verify if wbicx2(a1,a2) estimated for the scale corresponding to the oscillation has a larger value compared to wbicx2(a1,a2) estimated for the rest of the scales. A hypothesis test is required to detect nonlinear oscillations. The null hypothesis is that the oscillation is not nonlinear, and therefore, the bicoherence has a constant value at all the scales. In this test, the maximum of wbicx2(a1, a2) when a1 equals to the oscillation scale aosc and a2 changes over the feasible range of scales is taken as wbic max 2(aosc ,a 2) = max a2(wbicx 2(aosc ,a 2))

Figure 9. Required steps for oscillaton diagnosis.

the oscillation category for each oscillation frequency individually thanks to the use of wavelet transform. After analysis of all the variables of the process, the ones carrying similar oscillations should be grouped together and further analysis is required to find the loop that propagates the oscillation to the other ones. It is also noteworthy that due to the low pass filtering nature of chemical processes, nonlinear oscillations become more linear further away from the oscillation source. This has been used as a method to diagnose which loop has caused the oscillation when several loops carry the same kind of oscillation. Because bicoherence quantifies the nonlinearity, it can be said that the loop with the highest bicoherence value is propagating the oscillation to other loops. Diagnosing the source of the oscillation among several loops with similar oscillations requires application of other methodologies if the oscillation is due to controller tuning or due to a sinusoidal disturbance. The appropriate methodologies are being considered in the ongoing work.

(24)

The average of the wbicx2(a1,a2) when neither a1 nor a2 is equal to aosc is also estimated as wbic 2(a≠aosc) =

1 Σa1≠ aoscΣa2 ≠ aoscwbicx 2(a1,a 2) N

(25)

where N is the number of summands. If wbicmax2(aosc,a2) is significantly different from wbic2(a≠aosc), then it is concluded that the oscillation is due to nonlinearity. Because the probability distribution of wbicx2(a1,a2) values is unknown, the hypothesis test is performed on the basis of the Chebyshev’s inequality. Chebyshev’s inequality states that for a random variable x with mean μ and standard deviation σ the following inequality holds. G


Article

Industrial & Engineering Chemistry Research 4.3. Note Regarding the Selection of the Scales for Which the Bicoherence Should Be Estimated. A correct inference about nonlinearity in the process based on estimated bicoherence values requires careful selection of the scales for which the bicoherence is estimated. Specially estimation of bicoherence for scales larger than what is required may lead to wrong inference about the presence of nonlinearity. To illustrate this fact, consider ⎛ ⎛ π⎞ π ⎞⎟ x(t ) = 2sin⎜2π × 0.04 + ⎟ + 2sin⎜2π × 0.025 + ⎝ ⎝ 3⎠ 12 ⎠ ⎛ π⎞ + 2sin⎜2π × 0.065 − ⎟ + ε(t ) ⎝ 4⎠

and its bicoherence plotted in Figure 10 with consideration of scales up to 80.

Figure 11. Bicoherence of signal ⎛ ⎛ π⎞ π ⎞⎟ x(t ) = sin⎜2π × 0.04 + ⎟ + sin⎜2π × 0.025 + ⎝ ⎝ 3⎠ 12 ⎠ for scales up to 80.

Figure 10. Bicoherence of signal x(t) for scales up to 80.

From Figure 10 it can be observed that the bicoherence value again increases considerably for scales around 60−70. These scales approximately correspond to periods between 75 and 85. Periods 75−85 are approximate multiples of the two periods present in the signal. As a result, the wavelet coefficients at these scales have larger values compared to neighboring scales leading to a larger bicoherence value between periods 40 (frequency 0.025) and 75 (frequency 0.013). In this case the sinusoidal wave with frequency 0.04 is causing the coupling between these two frequencies. These large values of bicoherence that can wrongly indicate the presence of nonlinearity in the signals can always exist in frequencies that are multiples of the two frequencies present in the signal. The same bicoherence is again plotted in Figure 11 when the sinusoidal wave with frequency of 0.065 is removed from the signal. Even though there is no real phase coupling in signal

Figure 12. Wavelet coeficients of signal ⎛ ⎛ π⎞ π ⎞⎟ x(t ) = 2 sin⎜2π × 0.04 + ⎟ + 2 sin⎜2π × 0.025 + ⎝ ⎝ ⎠ 3 12 ⎠ ⎛ π⎞ + 2 sin⎜2π × 0.065 − ⎟ + ε(t ) ⎝ 4⎠

examination of the map of wavelet coefficients. For example, from Figure 12 it can be observed that the energy of the signal is mostly concentrated in scales up to 41. Bicoherence is better to be estimated for these scales. If no high value bicoherency is detected at these scales, it means the oscillations are linear.

⎛ ⎛ π⎞ π ⎞⎟ x(t ) = sin⎜2π × 0.04 + ⎟ + sin⎜2π × 0.025 + ⎝ ⎝ ⎠ 3 12 ⎠

5. CASE STUDY The introduced tools are used for diagnosis of a unit-wide oscillation through an oil sands industrial data set. The analysis consists of detecting oscillatory variables, estimating multiple oscillation frequencies and finally diagnosing the source of the oscillation. The data set contains the controller outputs as well as the process values (PV) of 24 loops. The detailed analysis of 5 of the loops carrying similar oscillations is presented here. The variables are denoted as Loop 1 to Loop 4 and Loop 30. Figure 13 plots 4 PVs along with their awps.

high value bicoherency is observed in periods from 75 to 85 similar to Figure 10. To avoid wrong detection of phase coupling in a variable, it is important to avoid estimating bicoherence values for periods that are approximately multiples of the oscillation periods present in the signal. Bicoherence can be estimated for scales up to the scale corresponding to the largest oscillation period within the variable. The largest oscillation period can be detected by running the oscillation detection algorithm on the data or visual H


Article


Figure 15. Wavelet bicoherence of down-sampled Loop 1 by a factor of 2.

Figure 13. Four industrial variables along with their awps.



From the right column of Figure 13, we can observe that all four variables have similar patterns of the power spectrum. Similar patterns indicate that all four loops are carrying similar oscillations and suffer from the same fault. The source of the oscillation could be one of these loops or other variables within the process that have similar patterns.

Figure 14. ACFs of the four industrial variables. I


Article



To categorize the oscillation, we first examine the ACF of the variables that are plotted in Figure 14. The ACFs of the variables show the presence of an oscillation. Even though the amplitude of the oscillation is very small in the ACFs, it does not decay to 0. The oscillation is steady with almost constant amplitude at all the time lags, which is not similar to the ACF of an oscillation induced by controller tuning. Therefore, the oscillation could be due to an external disturbance to these loops or could be due to nonlinearity in one of them. The next step is to examine the presence of nonlinearity in these 4 loops. Figures 15−18 plot the bicoherence of the four loops. Because scale 300 is too large for estimation of bicoherence, the variables are down-sampled by a factor of 2. Again, similar patterns in the plot of averaged bicoherence values of the four loops are observed. The average is with respect to time. All the bicoherence plots of the five variables show a peak at the scale around 120, which approximately corresponds to an oscillation period of 300 samples. Figure 19 plots the controller output of Loop 30 and its wavelet coefficients. Oscillation detection algorithm finds an


compared to those of the other four loops. This shows that Loop 30 is closer to the source of the oscillation whereas the other four loops receive the oscillation from the external source (propagated by Loop 30). Also, similar bicoherence values for Loop 1 to Loop 4 indicates that they are influenced by Loop 30 in a similar way.

6. CONCLUSION In this paper, a comprehensive algorithm is proposed that is capable of both detection and characterization of individual oscillatory components of variables in the presence of multiple oscillations and nonstationary signals. The independent characterization of multiple oscillations is viable due to the capability of wavelet transform in decomposing the variables to its components of different frequencies. Two hypothesis tests are proposed on the basis of the properties of wavelet bicoherence and wavelet power spectrum for the purpose of oscillation diagnosis.

■

ACKNOWLEDGMENTS This work is supported in part by Natural Sciences and Engineering Research Council of Canada and Alberta Innovates Technology Futures.

■

AUTHOR INFORMATION

Corresponding Author

*B. Huang. E-mail: [email protected]. ORCID

Biao Huang: 0000-0001-9082-2216 Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Jelali, M., Huang, B. Detection and diagnosis of stiction in control loops: state of the art and advanced methods, 1st ed.; Springer: London, 2010. (2) Horch, A. Benchmarking control loops with oscillation and stiction. In Process Control Performance Assessment, 1st ed.; Springer: London, 2006. (3) Choudhury, M. A. A.; Shah, S. L; Thornhill, N. F. Diagnosis of process nonlinearities and valve stiction: data driven approaches, 1st ed.; Springer: London, 2008.

Figure 19. Top: controller output of Loop 30 in time. Bottom: wavelet coefficients.

oscillation with period of 300 samples with 0 standard deviation from the ACF of the variable. The bicoherence estimated for the controller output of Loop 30, plotted in Figure 20, is much larger J


Article

Industrial & Engineering Chemistry Research (4) Naghoosi, E.; Huang, B. Diagnosis of Oscillations Between Controller Tuning and Harmonic External Disturbances. IEEE Transactions on Control Systems Technology 2015, 23 (4), 1283−1293. (5) Pan, Q.; Zhang, L.; Dai, G.; Zhang, H. Two denoising methods by wavelet transform. IEEE Transactions on Signal Processing 1999, 47 (12), 3401−3406. (6) Peng, Z. K.; Chu, F. L. Application of the wavelet transform in machine condition monitoring and fault diagnostics: a review with bibliography. Mechanical systems and signal processing 2004, 18 (2), 199− 221. (7) Yang, Y.; Li, X.; Liu, X.; Chen, X. Wavelet kernel entropy component analysis with application to industrial process monitoring. Neurocomputing 2015, 147 (5), 395−402. (8) Alabi, S. I.; Morris, A. J.; Martin, E. B. On-Line Dynamic Process Monitoring Using Wavelet-Based Generic Dissimilarity Measure. Chem. Eng. Res. Des. 2005, 83 (6), 698−705. (9) Zhu, K.; Wong, Y. S.; Hong, G. S. Wavelet analysis of sensor signals for tool condition monitoring: A review and some new results. International Journal of Machine Tools and Manufacture 2009, 49 (7), 537−553. (10) Bakshi, B. Multiscale PCA with application to multivariate statistical process monitoring. AIChE J. 1998, 44 (7), 1596−1610. (11) Lin, J.; Qu, L. Feature extraction based on Morlet wavelet and its application for mechanical fault diagnosis. Journal of sound and vibration 2000, 234 (1), 135−148. (12) Matsuo, T.; Sasaoka, H.; Yamashita, Y. Detection and Diagnosis of Oscillations in Process Plants. Lecture Notes in Computer Science 2003, 2773, 1258−1264. (13) Guo, Z.; Xie, L.; Horch, A.; Wang, Y.; Su, H.; Wang, X. Automatic detection of nonstationary multiple oscillations by an improved wavelet packet transform. Ind. Eng. Chem. Res. 2014, 53 (40), 15686−15697. (14) Plett, M. I. Transient detection with cross wavelet transforms and wavelet coherence. IEEE Transactions on Signal Processing 2007, 55 (5), 1605−1611. (15) Thornhill, N. F.; Huang, B.; Zhang, H. Detection of multiple oscillations in control loops. J. Process Control 2003, 13 (1), 91−100. (16) Naghoosi, E.; Huang, B. Automatic detection and frequency estimation of oscillatory variables in the presence of multiple oscillations. Ind. Eng. Chem. Res. 2014, 53 (22), 9427−9438. (17) Srinivasan, R.; Rengaswamy, R.; Miller, R. A modified empirical mode decomposition (EMD) process for oscillation characterization in control loops. Control Engineering Practice 2007, 15, 1135−1148. (18) Naghoosi, E., Huang, B. Diagnosis between sinusoidal oscillations and oscillatory colored noise. In American Control Conference (ACC); IEEE, 2014; pp 3518−3523. (19) Zhang, Z.; Moore, J. C. A Comment on ”Significance tests for the wavelet power and the wavelet power spectrum by Ge (2007). Annales Geophysicae-Atmospheres Hydrospheres and Space Sciences 2012, 30 (12), 1743. (20) Torrence, C.; Compo, G. P. A practical guide to wavelet analysis. Bull. Am. Meteorol. Soc. 1998, 79 (1), 61−78. (21) Hinich, M. J. Testing for Gaussianity and linearity of a stationary time series. Journal of Time Series Analysis 1982, 3, 169−176. (22) Rosenblatt, M.; van Ness, J. W. Estimation of the bispectrum. Ann. Math. Stat. 1965, 36, 1120−1136. (23) Nikias, C. L., Petropulu, A. Higher-order spectra: A nonlinear signal processing frame work; Prentice-Hall: Englewood Cliffs, NJ, 1993. (24) Choudhury, M. A. A.; Shah, S. L.; Thornhill, N. F. Diagnosis of poor control-loop performance using higher-order statistics. Automatica 2004, 40 (10), 1719−1728. (25) Van Milligen, B. P.; Hidalgo, C.; Sanchez, E. Nonlinear phenomena and intermittency in plasma turbulence. Phys. Rev. Lett. 1995, 74 (3), 395. (26) Van Milligen, B. P.; Sanchez, E.; Estrada, T.; Hidalgo, C.; Branas, B.; Carreras, B.; Garcia, L. Wavelet bicoherence: A new turbulence analysis tool. Phys. Plasmas 1995, 2 (8), 3017−3032. (27) Ge, Z. Significance testing for wavelet bicoherence and its application in analyzing nonlinearity in turbulent shear flows. Phys. Rev. 2010, 81 (5), 056311. K


Wavelet Transform Based Methodology for Detection and

Recommend Documents