A modified empirical mode decomposition (EMD) process for oscillation characterization in control loops

ARTICLE IN PRESS

Control Engineering Practice 15 (2007) 1135–1148 www.elsevier.com/locate/conengprac

A modified empirical mode decomposition (EMD) process for oscillation characterization in control loops Ranganathan Srinivasana,1, Raghunathan Rengaswamya,, Randy Millerb a

Department of Chemical Engineering, Clarkson University, P.O. 5705, NY 13699, USA b Honeywell Process Solutions, Thousand Oaks, CA, USA Received 11 May 2005; accepted 29 January 2007 Available online 30 March 2007

Abstract This paper focuses on a new method for characterizing oscillations in measurements from control loops in chemical processes. The proposed method uses a modified empirical mode decomposition (EMD) process to uncover all the zero-crossings in the data. It is shown that this procedure characterizes the following attributes of the oscillations: (i) non-constant mean, (ii) time instances when the oscillations are present, (iii) strength of the oscillations at different times, (iv) time period of each sweep of oscillation. The proposed approach is tested on more than 150 industrial control loop data and results obtained for 10 loops are discussed. The utility of the characterization approach in the area of oscillation analysis is also discussed. r 2007 Published by Elsevier Ltd. Keywords: Oscillation detection; Oscillation diagnosis; Performance monitoring; Signal analysis; Empirical mode decomposition

1. Introduction It is well known that performance degradation in control loops manifest as one or more of the following: (i) poor set point (SP) tracking, (ii) oscillations, (iii) poor disturbance rejection and (iv) excessive final control element variation. Industrial surveys over the last decade indicate that only about one-third of industrial controllers provide acceptable performance and about 30% (Bialkowski, 1993; Desborough & Miller, 2001; Ender, 1993) of all control loops oscillate. Since oscillations can lead to loss of energy and off-spec products, isolating the root cause for oscillations is important for improving the performance of the oscillating control loops. Automatic detection of oscillation can be used to focus the operators’ attention on control loops that might have performance problems and has become a standard activity. Analysis of oscillation in control loops is approached as a two-step process, namely, oscillation Corresponding author. Tel.: +1 315 268 4423; fax: +1 315 268 6654.

E-mail addresses: [email protected] (R. Srinivasan), [email protected] (R. Rengaswamy). 1 Current address: Automation and Control Labs, Honeywell Inc, Phoenix, AZ, USA. 0967-0661/$ - see front matter r 2007 Published by Elsevier Ltd. doi:10.1016/j.conengprac.2007.01.014

detection followed by oscillation diagnosis. The methods proposed for oscillation diagnosis assume that oscillation has, a priori, been detected. Oscillation diagnosis can suffer when the oscillation detection procedure fails or reports inaccurate oscillation periods. In this work, a novel oscillation characterization scheme that can directly provide metrics for oscillation diagnosis is proposed. First, a brief description of existing oscillation detection methods is provided.

1.1. Background Ha¨gglund (1995) proposed a technique to detect oscillating loops ‘‘on-line’’ using the IAE criterion and an ‘‘off-line’’ procedure to find the root cause of oscillation. This method does not assume any particular shape for oscillation and only requires that the measurement deviate significantly beyond a threshold. Ha¨gglund (1995) also proposed a diagnostic procedure for finding the source of oscillation and eliminating it. The diagnostic procedure is carried out by disconnecting the feedback (i.e. switching the controller to manual mode). This approach is simple and efficient, and probably the most comprehensive

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

1136

procedure available for diagnosing the root cause of oscillations. However, switching the controller to manual mode may not be always allowed. Further, it may not be possible for a control engineer to apply this approach on thousands of loops in a routine fashion. Thornhill and Ha¨gglund (1997) presented an off-line technique for detecting oscillations using a regularity factor. This method requires the user to specify the root mean square value of the noise and a threshold. Miao and Seborg (1999) used autocorrelation function and proposed an oscillation index (0—no oscillation, 1—highly oscillating). Industrial experience appears to be favorable and oscillations were detected with reasonable reliability. Thornhill and Ha¨gglund (1997) and Thornhill, Huang, and Zhang (2003) proposed a set of procedures to detect and diagnose oscillating loops using off-line data. Desborough and Harris (1992) combine the techniques of controller performance assessment along with operational signatures (OP-PV plots) and spectral analysis of the controller errors for diagnosis and have tested their approach on industrial loops. Recently, Tangirala, Shah, and Thornhill (2005) presented a new visualization tool termed as power spectral correlation map (PSCMAP) to detect plant wide oscillation. Matsuo and Sasaoka (2005) presented an interesting oscillation detection procedure using wavelet transforms. Wavelet analysis, currently a popular tool for analyzing non-stationary data, is essentially an adjustable window Fourier spectral analysis. It gives a time-frequency characterization of the signal and is useful in analyzing data with gradual frequency changes. Although wavelet analysis is promising, its success depends on several factors such as the number of levels, choice of mother wavelet and scaling filter, and

importantly, the level considered for analysis. Hence, automating wavelet analysis for thousands of industrial control loops to detect/diagnose oscillations will be a difficult task (Srinivasan, 2005). Figs. 1 and 2 show industrial loops that were diagnosed (by plant engineers) to limit cycle due to high static friction in control valves. It is seen that the oscillations have time varying frequency and magnitude, and is also intermittent. The use of FFT/ Spectral, correlation methods (Miao & Seborg, 1999; Pryor, 1982; Thornhill et al., 2003) or computing the area under the signal between successive zero crossings (Forsman & Stattin, 1999; Ha¨gglund, 1995; Thornhill & Ha¨gglund, 1997) can detect oscillations in the signal only under restrictive assumptions such as the data is stationary and the oscillations are periodic. However, these methods are prone to error in detecting and estimating the frequency or period when the process data has non-constant mean. Careful observation of industrial data often indicates more than one oscillation mode due to a combination of: slow changes in SP, constant effect of external disturbances and hardware faults.

1.2. Focus of this work The techniques presented so far were intended to detect the presence of oscillation in control loops. Some techniques also report the corresponding oscillation period. However, if by some means, from the data analysis one can extract the start and end time of oscillations along with their respective zero-crossings, then using these attributes of oscillation, the following tasks can be performed:

Process Output (PV)

200 PV SP

198 196 194 192

0

500

1000

1500

2000 2500 3000 Sampling instant

3500

4000

4500

5000

0

500

1000

1500

2000 2500 3000 Sampling instant

3500

4000

4500

5000

Controller Output (OP)

59.5

59

58.5

Fig. 1. Flow loop with slowly varying set-point (SP) and oscillating due to high static friction in the control valve.

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

1137

Process Output (PV)

2170

2150 2140 2130 2120 2110

Controller Output (OP)

PV SP

2160

0

500

1000

1500

2000 2500 3000 Sampling instant

3500

4000

4500

5000

0

500

1000

1500

2000 2500 3000 Sampling instant

3500

4000

4500

5000

32.2 32 31.8 31.6

Fig. 2. Flow loop exhibiting intermittent oscillation due to static friction.

1. Remove the non-constant mean from the oscillating data. This mean shifted data can improve the success rate of the available methods for oscillation diagnosis. 2. Isolate the portions of data where oscillations are present. For each extracted sweep of oscillation (obtained using the zero-crossing information), a shape analysis can be carried out. Shape analysis can be used to diagnose the presence of static friction in control valves. The zero-crossing information can also be used in other fault diagnosis techniques that are based on qualitative pattern matching techniques. Also, each extracted sweep of oscillation allows the diagnosis technique to handle varying amplitude and period. 3. Isolate the portions of data where oscillations are not present. This isolated portion of data may be used for more realistic controller performance calculations. 4. An iterative procedure on the extracted non-constant mean can be used to uncover and isolate other oscillation modes, if present. In this work, a procedure that will uncover the zerocrossing information in each oscillating loop is proposed. The uncovered zero-crossings in each oscillation mode give complete information regarding: (i) time period of each sweep of oscillation, (ii) detection of more than one oscillation mode by iteratively applying it on the extracted non-constant mean, (iii) amplitude and strength of each oscillation mode and (iv) time instances when oscillations are present. The proposed method uses the empirical mode decomposition (EMD) procedure introduced by Huang et al. (1998), the details of which are given later. Since the proposed method has to be routinely applied across

thousands of loops, the number of tuning parameters in the characterization algorithm has to be minimum (preferably zero). Special attention will be paid to this aspect of automation in this paper and the tuning parameters in this approach and how these can be automatically fixed will be discussed. This paper is organized as follows: the proposed methodology, which is a modified EMD technique, is discussed is Section 2. Illustrative examples and results for 10 industrial loops are presented in Section 3. The utility of the proposed oscillation characterization technique towards oscillation diagnosis is also discussed in Section 3. 2. Proposed method The proposed method has two main features: (a) the analysis is done in the time domain, and (b) the shape of the signal is preserved as only the non-constant mean is removed from the signal. There are two basic steps in the proposed oscillation characterization method: Step 1: The first step removes the non-constant mean (i.e. the low-oscillation mode) from the signal. A modified EMD procedure is employed in this step. Step 2: Cumulative area of the signal separated out using Step 1 is computed. The cumulative area averages the noise and reduces the number of spurious zero-crossings that may be reported. Extrema points of the cumulative area capture the zero-crossing points. These extrema points identified are reported as the zero-crossing points of the dominant oscillation mode. The data pre-processing step, i.e., replacing outliers and missing data is not discussed here. Quiet periods where

ARTICLE IN PRESS 1138

R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

oscillations are not present in an intermittently oscillating signal are identified using a qualitative trend approach (Rengaswamy & Venkatasubramanian, 1995). The quiet periods are removed before the proposed method is applied. 2.1. Step 1: removing the non-constant mean using modified EMD procedure In this step, the non-constant mean from the data is removed by an adaptive mean-shifting procedure using the modified EMD procedure. First, the EMD procedure is explained (Huang et al., 1998). Then the proposed approach is presented. 2.1.1. Empirical mode decomposition (EMD) Consider a sub-portion X 1 ðtÞ of a signal X ðtÞ, extracted between two consecutive extrema, e.g. between two minima occurring at t and tþ . This sub-portion corresponds to an oscillation starting and ending at a minimum and passes through a maximum. It represents the high frequency variations in the sub-portion considered and can be denoted as detail cðtÞ. For the same sub-portion considered, a low frequency or local trend rðtÞ such that we have X 1 ðtÞ ¼ rðtÞ þ cðtÞ between t ototþ is defined. Note that rðtÞ, termed as residual, can be a first-order trend or a slowly varying trend. The signal cðtÞ satisfies the properties of an intrinsic mode function (IMF). An IMF (Huang et al., 1998) must satisfy the following two conditions: (1) the number of maxima and minima points must be equal or differ by one from the number of zero-crossings and (2) the signal must have a zero-mean or in other words, the amplitude between each consecutive maxima and minima point must be symmetric. A simple example of an IMF is a sine wave with zero bias. The EMD method decomposes the data into several IMFs and serves two purposes: (a) separates the finest scale (or high frequency) from the signal and (b) makes the separated oscillation modes symmetric. From the original signal, say X ðtÞ, the first IMF ‘c1 ’ is identified. This represents the finest scale. This step involves finding the local maxima and minima points, respectively. After all extrema points are identified, an upper envelope that connects all the local maxima points and a lower envelope that connects all the local minima are obtained using a cubic spline interpolation between those respective points. Averaging the upper and the lower envelopes gives the nonconstant mean of the signal. The shifted signal c1 is then obtained by subtracting the computed non-constant mean from the actual signal. c1 is separated from the signal by

3. Similar to the upper envelope obtain the lower envelope by connecting all the minima points. The lower envelope is denoted as X l ðtÞ. 4. Compute the average r1 ðtÞ ¼ ðX u ðtÞ þ X l ðtÞÞ/2. 5. Compute the fastest oscillation mode (IMF) c1 ðtÞ ¼ X ðtÞ  r1 ðtÞ. Iterate steps 1–4 until c1 ðtÞ meets the two conditions of an IMF according to some stopping criterion. Huang et al. (1998) compare a threshold value (0.2–0.3) against the standard deviation computed between two consecutive shifting results. Once c1 is extracted, the residual r1 ðtÞ is decomposed as r1 ðtÞ ¼ c2 ðtÞ þ r2 ðtÞ, where c2 represents the second IMF. The separation into IMFs terminate when no further IMF can be extracted. This condition occurs when the longest component (or residual) becomes an increasing or a decreasing trend containing no extrema points. Effectively, EMD can be viewed as a process for separating the original signal into several components in descending order of frequency using a set of filter banks. The example trend shown in Fig. 3 is considered to demonstrate the EMD process. This signal was constructed from summation of two sine waves of 10 and 1 Hz of amplitude 2 and 3, respectively. A 2nd order bias was also added to the signal. The signal and their corresponding IMFs are given in Fig. 3. The IMFs effectively decompose the two frequency modes. From Fig. 3, it can be observed that the IMF extraction stopped once an increasing trend with no extrema in the residue r2 was found. However, this signal was noise free and in the presence of noise, first few details will represent the signal noise. For a detailed discussion on the EMD process, refer to Huang et al. (1998).

The residue ‘r1 ’ contains information related to longer period components. The steps for obtaining c1 are summarized below:

2.1.2. Modified EMD process Fig. 4 shows a loop that was diagnosed to oscillate due to high static friction (or in short stiction) in control valves. It can be observed that the effect of static friction was dominant over the unmeasured disturbance. Also, the unmeasured disturbance existed only for a short period (see Fig. 4 between sampling instants 2000 and 3110). From a truly comprehensive fault diagnosis perspective, it is desirable to isolate both the causes, i.e., unmeasured disturbance and stiction and their respective time occurrences. However, in general, a control engineer is often more interested in knowing the dominant cause(s) that affect the profitability of the loop rather than all the causes. This becomes a necessity when a control engineer becomes responsible for about 400 loops (Desborough & Miller, 2001). The PV data shown in Fig. 4 had three dominant oscillation modes:

1. Identify all the extrema points in X ðtÞ. 2. Use cubic spline interpolation to connect all the maxima points. This upper envelope is denoted as X u ðtÞ.

 Due to the unmeasured disturbance.  Stiction induced limit cycle oscillations.  Low frequency SP changes.

X ðtÞ  c1 ¼ r1 .

(1)

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

1139

X(t)

10 5 0 -5

0

0.5

1

1.5 Time (secs)

2

2.5

3

0

0.5

1

1.5 Time (secs)

2

2.5

3

0

0.5

1

1.5 Time (secs)

2

2.5

3

0

0.5

1

1.5 Time (secs)

2

2.5

3

C1

5 0 -5

C2

5 0 -5

r2

2 1 0 -1

Fig. 3. Signal X ðtÞ ¼ 2 sinð2p10tÞ þ 3 sinð2ptÞ þ 0:2t2 . Decomposition of X(t) into IMFs (C 1 ; C 2 ; r2 ). The residue r2 shows an increasing trend and hence no more IMFs can be extracted.

Process Output (PV)

34

33.5

33

32.5

Process output (PV) Set point (SP) 0

500

1000

1500

2000

2500 3000 Time (secs)

3500

4000

4500

5000

0

500

1000

1500

2000

2500 3000 Time (secs)

3500

4000

4500

5000

Controller Output (OP)

46.5 46 45.5 45 44.5

Fig. 4. Effect of static friction and unmeasured disturbance.

Since the oscillation was caused around the SP and stiction was diagnosed as the major cause, it will be of interest to see if the EMD process can extract the above dominant oscillation mode. However, it was assessed that extracting

the major cause with the current IMF process was not practically feasible.This can be readily explained using Fig. 5, which gives the IMFs of the PV data obtained using the standard EMD procedure (Huang et al., 1998). The

ARTICLE IN PRESS 1140

R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

Fig. 5. IMF decomposition for the PV data shown in Fig. 4.

IMFs 1–3 give the noise levels present in the PV data. IMFs 4 and 5 separate higher frequency contents that are due to external disturbances. IMF 6 separates the oscillation that originated due to stiction, IMF 7 gives the variation in PV due to the SP changes (see SP in Fig. 4) and the R gives the residual. Although the oscillations were separated out, ambiguity exists in choosing the level (similar to the problem faced when using Wavelets) for diagnosing the major cause. Even though an algorithm can be put together to assess all levels (or IMFs) for fault diagnosis, computationally this will be an expensive step and is very sensitive to noise. Since the focus of this work is on automation and computational simplicity, the option to analyze all the IMFs is not considered. Instead, a modification to the IMF extraction process that is geared towards identifying only the dominant oscillation modes is proposed. This is achieved by an adaptive thresholding to extract the dominant oscillation modes. The work of Huang et al. (1998) used the EMD procedure as a pre-processing step to obtain the instantaneous frequency (i.e. time-frequency) by applying Hilbert transform on the IMFs. The main objective in Huang et al. (1998) work was to find a tool that can extract the time-frequency information of a non-stationary, nonlinear signal.

However, in the proposed approach the idea of EMD procedure but with a meaningful threshold to extract all the dominant oscillation modes (different from IMFs) that are useful for oscillation diagnosis is used. Zero-crossings are then detected for each one of the dominant oscillation mode. A key difference from the original EMD procedure (Huang et al., 1998) is that the resulting output from the modified EMD procedure does not satisfy the properties of an IMF (due to the use of a threshold while identifying the extrema points, the IMF conditions get violated) and hence may have some extrema between consecutive zero-crossings. As will be seen from the results of this approach, removing this (IMF) restriction is an important shift that makes the proposed approach routinely applicable to a wide variety of control loops. In summary, to extract the dominant oscillation mode, it will be beneficial to remove the non-constant mean alone from the data in one single step rather than finding all IMFs present in the data. This can be done by modifying the first step in the standard EMD procedure (given in Section 2.1.1) with the following: instead of including all extrema points as suggested in the EMD process for extracting the IMFs, include only those extrema points that differ beyond a specified threshold value between successive extrema (minima followed by maxima). Using such an

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

approach, minor fluctuations in the data will be ignored and only the dominant cause will be extracted successfully. The proposed method was implemented on the data set shown in Fig. 4. The resulting mean shifted PV data is shown in Fig. 6. The non-constant mean is computed using a pre-specified threshold value (16% of the range of the signal). Since the value of threshold is critical in extracting the correct oscillation mode, an automated and meaningful approach for fixing the threshold value is discussed in the next section. 2.1.3. Automating the threshold value Current distributed control systems (DCS) allow access to the following parameters of a control loop through OPC standards (ProcessDoc, 2005): (1) PV/OP/SP values, (2) Operation range of the process output (PV), (3) Controller tuning constants and (4) the respective engineering units. For example, a temperature loop can have an operating range 32 to 100 F. The threshold for the modified EMD procedure indirectly depends on the operation range and is relative to the operating region dictated by the current SP. Hence, the threshold can be set in a meaningful way based on these two parameters and it will be easy to automate the threshold as a percentage of the operation range and the mean. In cases where the operational range information is unavailable, the threshold can be heuristically fixed based on the current operating range and the mean. 2.2. Step 2: identifying the zero-crossing points In step 1, the non-constant mean was removed to get the most dominant oscillation mode present in the signal. The

1141

next logical step is to identify the zero-crossings of this dominant oscillation. An auto-correlation test (Miao & Seborg, 1999; Pryor, 1982) is applied to ensure the presence of oscillation on the mean-shifted data obtained from Step 1. Computing the area between successive zero-crossings of the error signal between the SP and PV has been reported by Forsman and Stattin (1999). The ratio of successive areas of same sign and their time span is computed. However, this approach becomes ineffective for noisy data and reports spurious zero-crossings. To overcome this, a method based on the cumulative area of the data, similar to the CUSUM statistical charts that are used to detect shift (or change) in process data is proposed. The cumulative area for data X of length N is computed as follows: cusumi ¼

k¼i1 X

X ðkÞ;

i ¼ 1 . . . N.

(2)

k¼1

The cumulative area for the example data shown in Fig. 4 (that was mean-shifted in the Step 1, see Fig. 6) was computed using Eq. (2); this is shown in Fig. 7. The cumulative area changes slope at points where the signal crosses the zero. A local maxima or minima in the cumulative area corresponds to a zero-crossing point in the original signal. The cumulative area decreases (or increases) after the zero-crossing point till the next zerocrossing point. Hence, the extrema points in cumulative area plot capture the zero-crossing points of the original signal. The extrema points are obtained and a simple clustering technique that groups nearby zero-crossing points which do not differ by more than two sample instants is employed to further minimize the effect of fluctuations due to noise. The clustered zero-crossing

Process Output (PV)

34

33.5

33

32.5

PV data Computed non-constant mean 0

500

1000

1500

2000 2500 3000 Sampling instants

3500

4000

4500

5000

0

500

1000

1500

2000 2500 3000 Sampling instants

3500

4000

4500

5000

Mean shifted PV

0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8

Fig. 6. Modified mean-shifting process with a threshold value of 16% of the operating range of the signal.

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

1142

Cumulative area of the mean shifted PV data

300 250 200 150 100 50 0 -50

0

500

1000

1500

2000 2500 3000 Sampling instants

3500

4000

4500

5000

Fig. 7. Cumulative area computation for the mean shifted PV data (see Fig. 6b). The slope changing points represents the signal zero-crossing instants.

vector which characterizes the oscillating signal is finally reported. The two steps discussed above were implemented in MATLAB 6.5. Overall, the method is completely automated and the threshold value, which is the only tuning parameter, was fixed using the heuristics discussed in Section 2.1.3. A flowchart of the proposed approach is shown in Fig. 8. Table 1 gives the zero-crossings and the corresponding attributes obtained for the PV data shown in Fig. 4. In Table 1, the first column lists the identified zerocrossings, the second column gives the amplitude reporting the maximum and minimum values, the third column reports the period of oscillation (in sampling instants) and the last column gives the consistency of the identified oscillation period. The consistency attribute is defined as the ratio of total duration for which the oscillation occurred to the total duration of the signal. Here it is 4436 5000, i.e. 89%. The execution time in Intel P4 1.5 Ghz processor with 512 MB RAM on windows 2000 was about 320 CPU milliseconds (one-third of a second) per loop (5000 data points). It can be seen that the cumulative area approach averages the noise effects and reduces spurious zerocrossings. In this section, a modified EMD process to extract the most dominant oscillation mode was proposed. In the following section, a methodology to handle signals that have simultaneous presence of more than one dominant oscillating mode is given. 2.3. Addressing simultaneous occurrence of multiple timescale oscillations In practice, a control loop can exhibit more than one dominant oscillation mode. For ensuring correct diagnosis, it is important that all dominant oscillation modes are extracted. Fig. 9a shows data from a loop which has a fast oscillation mode imposed over a slow oscillation mode.

When simultaneous multiple oscillations are present in the data, the modified EMD process (described in Section 2.1) will extract the fastest oscillation mode lumping the other slower oscillation mode with the non-constant mean. Since the two oscillations are well separated in period (or frequency), an iterative procedure will allow these oscillation modes to be extracted. The iteration is performed by applying the modified EMD process on the residual (i.e. the non-constant mean separated from the previous step) to obtain the slower oscillation mode and its corresponding non-constant mean. For the example considered, both the oscillation modes are obtained in an iterative fashion (see Fig. 9b–e). It is important to note that the successive oscillation modes are obtained using an adaptive threshold, i.e., the threshold heuristic remains the same but the absolute values for the threshold changes since the signal being analyzed (i.e. the non-constant mean of the earlier step) is different at each iteration. Although this iterative approach appears very similar to the ‘sifting process’ in the Hilbert Huang Transform (Huang et al., 1998), there exists two important differences. The ‘sifting process’ in Hilbert Huang transform is performed without any threshold setting (equivalent to setting the threshold ¼ 0 in the modified EMD). The other important difference is that the extracted signal obtained using the modified procedure (in each step) does not possess the properties of an IMF. Also, it was demonstrated that ambiguities persists in choosing the correct level when the signal is decomposed into several IMFs using the EMD procedure outlined by Huang et al. (1998). Thus by specifying a threshold, it is seen that the dominant oscillation modes are generally obtained. 2.4. Comparison of EMD process with a high pass filter It is observed that the EMD process is equivalent to a set of filter banks separating frequency modes in a descending

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

1143

Fig. 8. Flowchart of the proposed algorithm.

Table 1 Oscillation characterization for the PV data shown in Fig. 4 Zero-crossings

[210 860 1193 1520 1960 2820 3210 3482 4017 4353 4646]

Amplitude

Period

Consistency (%)

Min

Max

Min

Max

0.28

0.52

629

1250

order. It may occur that a high pass filter with a cleverly chosen corner frequency may also yield the highest oscillation mode present in the signal. To study this a simple example signal ðX ðtÞ ¼ 3 sinð2ptÞ þ 0:5t2 Þ was considered. A corner frequency of 1 Hz (the actual frequency contained in the signal) was used to design a Butterworth filter of order 2. Fig. 10 gives the non-constant obtained using the EMD approach and the filtered signal. It is evident from Fig. 10 that the non-constant mean extracted by the EMD method is very close to the actual nonconstant mean. The non-constant mean obtained using the filter approach may yield comparable results with the EMD method only if the non-constant mean changes are very slow. However, for the filter approach, an exact corner frequency present in the signal has to be specified to yield the desired non-constant mean. This requirement cannot be met easily and a wrong choice of corner frequency may lead to incorrect oscillation attributes.

89

3. Results 3.1. Oscillation characterization 3.1.1. Non-constant mean Figs. 11a and b show the PV, SP and OP data of a selfregulating level loop. The level in the tank was controlled by regulating the inlet flow through a PI controller. The oscillation in the loop was diagnosed as being caused by pulsating inlet flow. Since the oscillations are external to the loop, the cross-correlation test suggested by Horch (1999) was applied. The cross-correlation test on the raw data failed due to the non-constant nature of the controller output (see Fig. 11c). However, after meanshifting both the OP and PV data (see Figs. 11d and e), the cross-correlation test (see Fig. 11f) correctly identified the source of oscillations as external to the loop. It can be seen that simple mean shifting process improved an already published

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

1144

410 400

0

500

Mean shifted data after 1st iteration

390

1000

1500

2000 2500 3000 Sampling instants

10

420

5

415

Non-constant mean after 1st iteration

Process Output (PV)

420

0 -5 -10

0

Non-constant mean after 2nd iteration

Mean shifted data after 2nd iteration

4500

5000

405 0

2000 4000 Sampling instants

6000

0

2000 4000 Sampling instants

6000

415

10 5 0 -5

4000

410

400

1000 2000 3000 4000 5000 Sampling instants

3500

0

2000 4000 Sampling instants

410 405 400

6000

Fig. 9. Extraction of multiple time-scale oscillations using two levels. (a) Simultaneous presence of multiple time-scale oscillations. (b) Mean shifted signal obtained after first iteration. (c) Non-constant mean (or residual) obtained after first iteration. (d) Mean shifted signal obtained after second iteration. (e) Non-constant mean (or residual) obtained after second iteration.

6 Signal 2 x(t) = 3*sin(2 π t) + 0.5t Actual non-constant 2 mean (0.5t ) Mean obtained using High-pass filter Mean obtained using EMD

5

Signal Amplitude

4 3 2 1 0 -1 -2 -3

0

0.5

1

1.5 Time (secs)

2

2.5

3

Fig. 10. Signal X ðtÞ ¼ 3 sinð2ptÞ þ 0:5t2 . Comparison of non-constant means obtained using a high pass filter and EMD methods.

method (Horch, 1999). Interestingly, as the loop was not generating oscillations due to improper controller settings, the variance caused due to this oscillation alone (i.e. A20 /2, where A0 is the oscillation amplitude) can be subtracted

from the output variance before calculating the actual controller performance index (Desborough & Harris, 1992), thus computing realistic performance indices for oscillating loops as suggested by Ettaleb (1999).

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

5

60

50 45

Level (PV) Set point (SP) 0

200

400 Time (secs)

600

80

10

60

5

40 20 0

0

200

400 Time (secs)

600

200

400 Time (secs)

600

800

0

200

400 Time (secs)

600

800

-5

Cross-correlation coefficients

1

0

-0.5 -500

0

0

-10

800

0.5 Cross-correlation coefficients

0

-5

800

Mean shifted OP data

Controller output (%)

Mean shifted PV data

Level (%)

55

40

1145

0 Lags

500

0.5 0 -0.5 -1 -500

0 Lags

500

(c) Fig. 11. Level loop oscillation diagnosis. (a) Raw PV and SP data. (b) Raw OP data. (c) OP/PV Cross-correlation coefficients computed using Raw data. (d) Mean shifted PV data. (e) Mean shifted OP data. (f) OP/PV Cross-correlation coefficients computed using mean shifted data.

3.1.2. Shape analysis Since stiction in control valves leave distinct qualitative shapes such as square, triangular and sinusoidal in the controller output (OP) and process variable (PV) data (Rengaswamy, Ha¨gglund, & Venkatasubramanian, 2001), each extracted sweep of oscillation (three consecutive zero crossings constitute one sweep) obtained from the proposed technique can be used to find the shapes in PV and OP to confirm the presence of stiction. Such an approach was applied in this work (Srinivasan, Rengaswamy, & Miller, 2005) for diagnosing stiction induced limit cycles. Also, the accuracy of other published methods for stiction diagnosis (Choudhury, Shah, & Thornhill, 2004; Singhal & Salsbury, 2005) can also improve by using the mean shifted data. 3.2. Application to a large data set To investigate the effectiveness of the proposed oscillation characterization approach, a collection of 150 industrial loops were considered. Fig. 12 shows the OP/

PV/SP data for the 10 loops. It can be seen from Fig. 12 that several loops oscillate and the SPs vary, making the data fairly non-constant. Results for the proposed oscillation characterization technique are summarized in Table 2. Table 2 lists the number of identified zero-crossings found manually versus the number of zero-crossings identified by the oscillation characterization algorithm. The average error between the identified position of the actual and identified zero-crossings is given in sampling units. The accuracy is defined as the ratio of the number of identified zero-crossings and the total number of zerocrossings in the data in percentage terms. It was seen that the average error in the zero-crossings were around 1–2 sampling units. Also more than 95% of the time the zerocrossings were identified correctly. The last column in Table 2 gives the threshold used in Step 1 (i.e. removing non-constant mean) of the oscillating technique. The threshold value was automated in a heuristic manner based on a linear combination of the mean and the range ð0:25  mean þ 0:3  rangeÞ of the signal under analysis. All

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

1146

62

200

60

L1

200.5

199.5

100

200

300

400

500

58 640

2

62

L2

2.50

L3

1.5

100

200

300

400

500

60

180

750

16

70

14

100

200

300

400

500

0

65

L4

100

200

300

400

500

62.6

103

74

L5

74.50

L6

102.5

100

200

300

400

500

73.5

2.50

650

2

60

1.5

L7

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

62.8

103.50

100

200

300

400

500

55

2.40

660

2.2

64

2

L8

200

630

0.8

100

200

300

400

500

62

2.50

700

2

65

1.5

L9

100

100

200

300

400

500

60

3.50

800

3

70

2.5

100

200

300

400

500

2.50

60

L 10

700 65

2

0

100

200

300

400

500

Process output Set-point

60

0

Controller output (%)

Fig. 12. (a) Set-point (SP), Process output (PV). (b) Controller Output (OP).

results reported in Table 2 were obtained with the heuristic unchanged. The new approach failed to characterize the oscillations for some of the loops that had very long oscillation periods with one or two extrema. For the oscillation characterization to succeed, at least one maxima and one minima point

is required for constructing the upper and lower envelopes for finding the non-constant mean. This requirement can be met by analyzing longer periods of data. Some loops in the 150 loop data set reported a small number of spurious zero-crossings. On further analysis, it was confirmed that a better threshold value would remove the spurious

ARTICLE IN PRESS R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

1147

Table 2 Tabulation of results using oscillation characterization algorithm Loop No

Actual no. of zero crossings

Identified no. of zero crossings

Average error with each zero crossing point

Accuracy (%)

Threshold as % of signal range

1 2 3 4 5 6 7 8 9 10

42 21 23 28 20 61 55 77 33 61

42 21 23 27 20 59 52 77 32 61

1 1 0 3 0 2 2 0 2 1

100 100 100 96 100 97 94 100 97 100

15.6 23.9 28.6 26.4 24.4 18.9 19.1 20 18 19.4

zero-crossings in these loops. However, the proposed characterization algorithm found both the number and position of zero-crossings accurately given that the tuning parameter heuristics was unchanged.

4. Conclusion A new oscillation characterization algorithm was proposed and demonstrated successfully on a large industrial data set. The proposed approach uses a modified empirical mode decomposition (EMD) procedure. All the dominant oscillation modes that are useful for oscillation diagnosis are extracted using this approach. The idea of zero-crossing detection on each dominant oscillation mode was proposed. A main advantage of this approach is that the analysis can be done completely in time domain and is automated with a single intuitive tuning parameter. The new approach opens up avenues for the plant engineers to look beyond detecting oscillations. As seen in Section 3.1, removing non-constant mean by itself improves already proposed methods for oscillation diagnosis. It was pointed out that the shape analysis of each sweep of oscillation can be useful for stiction diagnosis. These applications indicate the usefulness of the proposed approach. As an extension, an oscillation diagnosis framework that utilizes the oscillation attributes obtained using the proposed approach along with the existing oscillation diagnosis techniques can be pursued.

Acknowledgments The research was financially supported by Department of Chemical Engineering, Clarkson University. We also thank Mark Cook, Senior technician for setting the level loop experiment. The authors also thank Dr. Sridhar Narasimhan for his valuable feedback on the manuscript. The authors thank the National Science Foundation for partial support through the Grant CTS-0553992.

References Bialkowski, W. L. (1993). Dreams versus reality: A view from both sides of the gap. Pulp and Paper Canada, 94(11), 19–27. Choudhury, M. A. A. S., Shah, S. L., & Thornhill, N. F. (2004). Detection and quantification of control valve stiction. Boston, USA. DYCOPS. Desborough, L. D., & Harris, T. J. (1992). Performance assessment measures for univariate feedback control. Canadian Journal of Chemical Engineering, 70, 1186–1197. Desborough, L. D., & Miller, R. M. (2001). Increasing customer value of industrial control performance monitoring—Honeywell’s experience. Arizona, USA. CPC-VI. Ender, D. B. (1993). Process control performance: Not as good as you think. Control Engineering, 9, 180–190. Ettaleb, L. (1999). Control loop performance assessment, oscillation detection. Ph.D. Dissertation, Department of Electrical and Computer Engineering, University of British Columbia, Canada. Forsman, K., & Stattin, A. (1999). A new criterion for detecting oscillations in control loops, Karlsruhe, Germany. Proceedings of the European control conference. Ha¨gglund, T. (1995). A control-loop performance monitor. Control Engineering Practice, 5(11), 1543–1551. Horch, A. (1999). A simple method for the detection of stiction in control valves. Control Engineering Practice, 7, 1221–1231. Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, E. H., Zheng, Q., et al. (1998). The empirical mode decomposition and the Hilbert Spectrum for non-linear and non-stationary time series. Proceedings of Royal Society of London A, 454, 903–995. Matsuo, T., & Sasaoka, H. (2005). Application of wavelet transform to process behavior monitoring. Journal of Society of Instrument and Control Engineers, 44, 139–142. Miao, T., & Seborg, D.E. (1999). Automatic detection of excessively oscillating feedback control loops, Hawaii. IEEE Conference on Control Applications (pp. 359–364). ProcessDoc, Matrikon Inc. (2005). hhttp://www.matrikon.comi. Pryor, C. (1982). Autocovariance and power spectrum analysis—derive new information from process data. Control Engineering, 2, 103–106. Rengaswamy, R., Ha¨gglund, T., & Venkatasubramanian, V. (2001). A qualitative shape analysis formalism for monitoring control loop performance. Engineering Applications of Artificial Intelligence, 14, 23–33. Rengaswamy, R., & Venkatasubramanian, V. (1995). A syntactic patternrecognition approach for process monitoring and fault diagnosis. Engineering Applications of Artificial Intelligence, 8(1), 35–51. Singhal, A., & Salsbury, T. I. (2005). A simple method for detecting valve stiction in oscillating control loops. Journal of Process Control, 15, 371–382.

ARTICLE IN PRESS 1148

R. Srinivasan et al. / Control Engineering Practice 15 (2007) 1135–1148

Srinivasan, R. (2005). Control loop performance monitoring: Modeling, diagnosing and compensating stiction phenomenon in process control valves. Ph.D. Dissertation Department of Chemical Engineering, Clarkson University, USA. Srinivasan, R., Rengaswamy, R., & Miller, R. M. (2005). Control loop performance assessment 1: A qualitative pattern matching approach for stiction diagnosis. Industrial and Engineering Chemistry Research, 44, 6708–6718.

Tangirala, A. K., Shah, S. L., & Thornhill, N. F. (2005). PSCMAP: A new tool for plant-wide oscillation detection. Journal of Process Control, 15, 931–941. Thornhill, N. F., & Ha¨gglund, T. (1997). Detection and diagnosis of oscillation in control loops. Control Engineering Practice, 5, 1343–1354. Thornhill, N. F., Huang, B., & Zhang, H. (2003). Detection of multiple oscillations in control loops. Journal of Process Control, 13(1), 91–100.