Root Cause Analysis of Linear Closed-Loop Oscillatory Chemical

Article pubs.acs.org/IECR

Root Cause Analysis of Linear Closed-Loop Oscillatory Chemical Process Systems S. Babji,† U. Nallasivam,‡ and R. Rengaswamy*,† †

Department of Chemical Engineering, Texas Tech University, Lubbock, Texas, United States Department of Chemical Engineering, Clarkson University, Potsdam, New York, United States

‡

ABSTRACT: In a single-input single-output (SISO) closed-loop system, under constant or nonoscillatory set-point, oscillations in the output can occur mainly because of one or a combination of the following reasons: (i) presence of stiction in control valve, (ii) marginally stable control loops (because of aggressively tuned controller/changes in process time constant/gain/time delays or a combination of them), and (iii) disturbances external to the loop. The presence of these oscillations can propagate plantwide and force plants to deviate from optimal operating conditions. Therefore, it is essential to develop techniques that can diagnose the source of oscillations in control loops. Several data-driven methods have been developed to address the diagnosis problem by focusing on only one of the causes for oscillations. In the current study, an off-line data driven approach is developed to identify the root cause for oscillations in control loops using the routine plant operating data. Unlike the existing techniques, this approach identifies and distinguishes between the three major causes for oscillations in linear closed loop systems. The proposed methodology combines both parametric (Hammerstein-based approach) and nonparametric (Hilbert−Huang spectrum) schemes for performing oscillation diagnosis. Simulation and industrial case studies that demonstrate the utility and limitations of the proposed method for root cause diagnosis in closed loop systems are discussed.

1. INTRODUCTION AND LITERATURE REVIEW A process plant may have a few hundred control loops to several thousands depending on the complexity of the plant. Oscillations are a common type of plant-wide disturbance whose detection and diagnosis have generated considerable interest in recent years.1−4 Oscillations in control loops increase variability in product quality, accelerate equipment wear and may cause other issues that could potentially disrupt the operation.5 Therefore, it is important to diagnose the root cause for oscillations in process plants. Further, increasing emphasis on plant safety and profitability strongly motivates the search for techniques to detect and diagnose plant-wide oscillations. 1.1. Approaches for Plant-Wide Oscillation Detection and Diagnosis. The most important tasks in plant-wide oscillation detection and diagnosis are: (i) detection of presence of one or more oscillations (along with the periods of oscillations) indicated by a regular pattern in the data, (ii) detection of number of sources causing the oscillations and isolation of these oscillatory loops and, and (iii) diagnosis of causes for oscillations in the isolated control loops. Several methods have been developed for detection and determination of the periods of oscillations in control loops. Oscillation detection techniques can be broadly classified as (i) time domain techniques,6 (ii) auto-covariance function based techniques (ACF),7,8 and (iii) spectral peak detection techniques.9 Multivariate techniques like spectral principal component analysis (PCA),2 spectral nonnegative matrix factorization (NMF),10 spectral independent component analysis (ICA),11 spectral envelope method,9 and graph theoretical approaches12 have been proposed for identification of number of sources and isolation of oscillatory control loops. An excellent review of these multivariate techniques for © 2012 American Chemical Society

detection and isolation of oscillatory control loops can be found in Thornhill and Horch.3 Identification of root causes for oscillations is the final step in plant-wide oscillation diagnosis. Several surveys performed in different industrial plants (paper mills, refineries, cement, steel, and mining) at various times13,14 indicate that the major causes for oscillations in control loops are (i) stiction in control valves, (ii) marginally stable control loops (because of aggressively tuned controller/changes in time constant/gain/time delays or a combination of them), (iii) oscillatory external disturbances, (iv) valve saturation, and (v) nonlinearities in the process. In this work, an algorithm is developed to distinguish between the above first three major causes for oscillations. 1.2. Approaches for Stiction Detection. Stiction in control valves is a major cause for oscillations in control loop. Stiction in control valves introduces nonlinearities in the process output.15 Numerous techniques utilize this information for stiction identification in linear control loops. Stiction detection techniques can be broadly classified as (i) shapebased approaches,16−19 (ii) frequency domain based approaches,15,20 and (iii) model-based approaches21−25 In the proposed algorithm, Hammerstein model based approach21 is used for valve stiction detection. Though a variety of methods are available for stiction detection, only a few techniques identify the oscillations caused due to marginally stable control loop and external oscillatory disturbances.26 The difficulty involved in performing this task along with the various existing approaches for distinguishing Received: Revised: Accepted: Published: 13712

October 22, 2011 February 21, 2012 August 27, 2012 August 28, 2012 dx.doi.org/10.1021/ie2024323 | Ind. Eng. Chem. Res. 2012, 51, 13712−13731

Industrial & Engineering Chemistry Research

Article

process changes that result in marginally stable control loop and external oscillatory disturbances. In industrial plants, it is desired to have plant-wide oscillation diagnosis algorithms that meets the following conditions: (i) algorithm should perform well just with the available routine plant operating data (OP and PV) without exciting the plant using external signals and (ii) the algorithm must work without any prior model information. The proposed root cause analysis (RCA) algorithm for linear processes is an off-line diagnosis algorithm designed to satisfy these realistic requirements that make this problem hard to solve. Further, minimal and realistic assumptions are made regarding the disturbance and noise corrupting the process. As mentioned before, the proposed solution approach uses Hammerstein model identification and Hilbert−Huang Transform (HHT) for oscillation diagnosis. The use of Hammerstein model identification for stiction detection21−25 has been well established in the literature and will be briefly reviewed later in the solution approach section. However, the use of HHT in the area of root cause analysis has not been pursued before. In view of this, before presenting the use of HHT for root cause analysis, HHT is introduced in the context of oscillating signals.

between marginally stable control loops (due to presence of aggressively tuned controller/changes in process parameters) and disturbance caused oscillations are discussed in section 4.2. 1.3. Contribution of This Work. In this article, a method for diagnosis of three major causes for oscillations in closed loop systems using regular plant operating data is presented. Though a mini version of this work was presented by the authors at IFAC 2011 conference,27 theoretical frequency domain analysis arguments, detailed discussion of the algorithm with its possibility of extension to handle nonstationary disturbances and drifts along with testing on various industrial case studies are presented here. The key aspects of the proposed data-driven technique are (i) combinations of both parametric Hammerstein algorithm (stiction identification algorithm) and nonparametric (Hilbert−Huang Transform) technique (used to distinguish oscillations due to marginally stable control loop and external disturbance) for root cause analysis and (ii) identification of unique frequency domain signature for diagnosis of disturbance and controller related oscillations. The problem definition is presented first.

2. PROBLEM DEFINITION A typical closed loop system with stiction in control valve is shown in Figure 1. In Figure 1, u(t) (OP) is the controller

3. PRELIMINARIES: HILBERT HUANG TRANSFORM OF PROCESS SIGNALS Recently, Huang et al.28 developed HHT to analyze the timefrequency characteristics of time-dependent signals. HHT comprises of two distinct parts: (i) empirical model decomposition (EMD) of the signal to identify the so-called intrinsic mode functions (IMFs) and (ii) Hilbert transform of the IMFs to obtain instantaneous frequencies. EMD is an adaptive technique which is derived from the assumption that any signal consists of characteristic oscillations that are separated on a time-scale. The first step in HHT is the computation of IMFs from the data. A detailed explanation on the computation of IMFs using the EMD technique is provided in the Appendix A. In the second step, Hilbert transform of all the IMFs are computed, which yields the amplitude spectra along with the instantaneous frequencies (definition of Hilbert transform and a detailed discussion on computation of instantaneous frequencies are presented in Appendix A). The

Figure 1. General process control loop.

output, Gp and Gc denote the linear process and controller dynamics, yp(t) represents the uncorrupted process output, and y(t) (PV) is the process output corrupted with white noise e(t) and oscillatory disturbance d(t). The fundamental problem that is being addressed in this article is to distinguish between oscillations due to stiction, aggressively tuned controller/

Figure 2. (a) Sum of two sinusoidal signals (frequency −0.03 and 0.08). (b) Power spectrum obtained from Fourier transform. (c) One of the two individual IMFs obtained from EMD. (d) HH spectrum of the signal. 13713

dx.doi.org/10.1021/ie2024323 | Ind. Eng. Chem. Res. 2012, 51, 13712−13731


Article

Figure 3. (a) Nonlinear signal (sum of sinusoidal frequency of 0.03 and 0.08 is passed through a cubic function). (b) Power spectrum obtained from fourier transform. (c) One individual IMF obtained from EMD. (d) HH spectrum of the signal.

Figure 4. Amplitude-based discrimination analysis using HH spectrum for mixture of sinusoidal signals with different magnitude ratios I.

at all times and will be equal to the frequency of the sine wave. While this is obvious for a single sine, consider a signal which is a sum of two sine waves. Now at every time, there are two frequencies and hence a single instantaneous frequency cannot be defined. HHT solves this difficulty by first separating the signal into two components in the time domain and then calculating the instantaneous frequency for each individual time-resolved signal. This separation is achieved through the EMD process. This idea is explained in Figure 2, where Figures 2c and 2d depict the IMFs and the corresponding instantaneous frequencies of a signal that is a summation of two sine signals. The notion of instantaneous frequency results in nonlinear distortions to frequency components being handled in a completely different manner in HHT when compared to FT. In general, FT of the nonlinear transformed time signal will show several frequency components that are integer multiples of the original frequency component. This is usually referred to as

amplitude of the instantaneous frequencies from all IMFs are plotted as a function of time, either as a 3-D plot (z axis surface is amplitude) or as 2-D plot (color scale represents amplitude). This plot of superposition of all instantaneous frequencies obtained from all IMFs in the time-frequency plane along with their amplitudes is called as the Hilbert−Huang spectrum. A detailed discussion on the differences between HHT and all the other competing representations such as Fourier transform (FT), short time Fourier transform (STFT), and wavelet transform (WT) can be found in the seminal paper by Huang et al.28 A layman description of why HHT is important for oscillation characterization and how it is different from FT, and other techniques derived from FT such as STFT is presented here. The key difference between HHT and FT is the notion of instantaneous frequency at every time instant that HHT uses to describe signals. For example, for a single sine wave, the instantaneous frequency of the signal will be constant 13714



Article

Figure 5. Amplitude-based discrimination analysis using HH spectrum for mixture of sinusoidal signals with different magnitude ratios II.

formed signal compared to HH spectrum of the sum of sinusoidal signal. In fact, this observation in HHT is similar to the fact that Fourier Transform (FT) of nonlinearly transformed signal will show lower amplitudes at fundamental frequency, provided both the time domain sinusoid and nonlinearly transformed sinusoid have same power. According to Parseval’s theorem the total power in time domain should match that of frequency domain. Since FT uses harmonics to represent nonlinearities, the power is spread across both the fundamental and harmonics resulting in decreased power at fundamental frequency in the nonlinear transformed signal. Similarly, HHT uses intrawave modulations around fundamental frequency to represent nonlinearities resulting in spread of power around the fundamental frequency. Thus, the power of the sum of sinusoidal signal at fundamental frequencies is higher and clearly visible compared to that of the power spread (modulations) around the fundamental frequency. Since the amplitude at one of the frequencies in the sum of sinusoids is lowered, the modulation/power spread around the corresponding fundamental frequency is less clear in HHT of the nonlinearly transformed signal. Now, the magnitude ratio of the frequencies is changed to 10:1. The signals obtained are shown in Figure 5a and c. The corresponding HH spectrum plot is shown in Figure 5b and d. It can be clearly seen from the HH spectrum of nonlinear transformed signal that modulation is present around only one frequency. However, the sum of sinusoidal signal shows amplitude at two frequencies. Further the magnitude ratio of the frequencies is changed to 50:1 and the corresponding signals are shown in Figure 5e and g. The corresponding HH spectra are shown in Figure 5f and h. It can be clearly seen from the Figures that the HH spectrum of both the signals indicate that the low frequency content is absent compared to high frequency component in the signal. From this example, the key points to be observed are as follows: (1) HH spectrum of the nonlinearly transformed sum of signals (with magnitude of one signal higher compared to the other) shows modulations around the fundamental frequency, which has higher amplitude. (2) A suitable nonlinear transformation

harmonic distortion. However, when the same nonlinearly transformed time signal is analyzed using HHT, the frequency spectrum displays a drastically different behavior with oscillations/modulations around the fundamental frequency. The importance of this representation can be clearly seen when we look at a nonlinear transformation of an input signal that is a summation of two sines as shown in Figure 3. FT spreads the two frequencies in the frequency scale due to harmonic distortions and hence from a FT plot it becomes difficult to infer that the input signal had two frequency components when nonlinearities are involved. This inference is further exacerbated when there is noise present in the signals. However, the Hilbert−Huang (HH) spectrum shown in Figure 3 shows two waves around the fundamental frequencies that are clearly separated and easy to identify. 3.1. Amplitude-Based Analysis Using HHT. Now, an example to illustrate the use of HHT in the identification of nonlinearly transformed signals with varying amplitudes is presented. This is an important example which is related to the identification of disturbance and marginally stable control loop caused oscillations which is discussed in section 4.2. Let us consider: (i) sum of two sinusoidal signals of frequency (f1 = 2π0.05 and f 2 = 2π0.24) and (ii) nonlinear transformation (cubic nonlinearity) of the sum of two sinusoidal signals (frequency same as (i) which is f1 = 2π0.05 and f 2 = 2π0.24). The low and high frequency components in both the signals initially have a magnitude of one as shown in Figure 4a and c. The corresponding HH spectra of the two signals are shown in Figure 4b and d. Notice that the HHT spectrum of both the signals show considerable magnitude at both frequencies. Now, the amplitude of the high frequency component (0.24) is increased so that the magnitude ratio (ratio of magnitude of high and low frequency) is 5:1. The signals obtained are shown in Figure 4e and g. Now again, HHT is computed for the sum of sinusoids (with gain ratio 5:1) and the cubic nonlinear transformed signal. The corresponding HH spectrum plots are shown in Figure 4f and h. It can be clearly seen from the HH spectrum plot that the amplitude of the low frequency component is significantly reduced in the nonlinearly trans13715



Article

al.21 that the nonlinear stiction element can be decoupled from the linear process dynamics. This algorithm uses the following one parameter valve model given by

of the sum of the signals can be used to obtain a HH spectrum as described in the above observation. The above two key observations can be used to distinguish between oscillations in marginally stable control loop (due to aggressively tuned controller/changes in process parameters) and disturbance caused oscillations based on their amplitudes. In other words, when oscillations are generated at different parts of a control loop (controller, disturbance) as shown in Figure 1 and when they pass through a nonlinear element such as stiction, separating the different components based on their amplitudes becomes possible which is discussed in section 4.2.

⎧ if |u(t ) − x(t − 1)| ≤ d ⎪ x(t − 1) x(t ) = ⎨ ⎪ otherwise ⎩ u(t )

(1)

Here x(t) and x(t − 1) are past and present stem movements, u(t) is the present controller output, and d is the valve stiction band. d is expressed in terms of the percentage or fraction of valve movement corresponding to the amount of stiction present in the valve. The estimation of this parameter is achieved by decoupling the stiction parameter estimation from the estimation of the linear process dynamics. This method uses a linear data driven model along with nonlinear stiction parameter d to fit the data between OP and PV. The squared errors between the model predicted and process outputs are summed over a period of time to obtain the Root Mean Squared Error (RMSE). The value of stiction parameter ’d’ corresponding to the model with minimum RMSE is used for stiction detection. A nonzero value of d indicates stiction and quantifies its level while a zero value implies the absence of stiction in the control valve. Further, an Hammerstein based approach with two parameter model as discussed in Choudhury et al.29 is also used for stiction detection. From the industrial case studies, it is shown that stiction detection results obtained from use of one/two parameter model in the Hammerstein based approach remains the same. In fact, there are several data based models for stiction and a discussion of these models is provided in Kano et al.30 Once stiction detection is performed, it is necessary to distinguish between oscillations due to marginally stable control loop and external disturbances. This is discussed next. 4.2. Distinction between Oscillations in Marginally Stable Control Loop and Disturbance Caused Oscillations. Identification of external oscillatory disturbances and oscillations due to marginally stable control loop (caused by aggressively tuned controller, changes in process parameters, or both) is a challenging task because of the following constraints arising in industrial environment:26,31 (i) constant set-points, process and controller outputs are the only routine operating signals available for diagnosis, (ii) knowledge about the process is limited, and (iii) generally, no open/closed-loop experiments can be performed. There is hardly any solution in literature to differentiate between oscillations due to marginally stable control loop and external disturbances that will work with these constraints (Chapter 7 of ref 31). However, a few frameworks are listed in the literature that can be used (under certain assumptions) to distinguish between these two sources of oscillations. 4.2.1. Existing Frameworks. The various ways to identify the oscillations due to external oscillatory disturbances and marginally stable control loop can be classified as follows: (1) Process and data-driven model based methods.26 These methods use the process and controller transfer functions to identify the ultimate frequency of oscillation. However, in reality, process transfer function is seldom available. Data driven models are used to obtain disturbance models from which distinction can be made.25 However, information about the disturbance corrupting the process is rarely available, and therefore, disturbance models can seldom be obtained. Further, when there is no stiction in the control loop, there will not be enough persistence of excitation for identification of a reliable

4. SOLUTION APPROACH Our solution approach for root cause diagnosis of oscillating loops is depicted in Figure 6. In this work, it is assumed that

Figure 6. Flowchart of RCA algorithm.

only one root cause is present at a time. While multiple problems at the same time are definitely possible, these scenarios are probabilistically rarer than the single fault scenario. Further, with the available routine operating data, comprehensive off-line diagnosis solutions do not exist even for a single failure scenario at this time. This is because of the realistic constraint that is imposed in the problem definition; process models are not available a priori and the algorithm should work with the available PV-OP data with no external interruptions to the plant. This makes a solution to this problem harder. One of the key ideas in the solution approach is the initial diagnosis that separates stiction from the other two root causes. It has been shown21,23,24 that stiction detection can be decoupled effectively from marginally stable control loop and disturbance related problems. The algorithm terminates if stiction is detected due to the assumption of single root cause. If stiction is not detected then further analysis is needed. This is accomplished through the steps shown in Figure 6 and discussed in section 4.2. Before these steps are discussed, a brief description of the well established stiction detection idea is summarized. 4.1. Stiction Detection in Linear Closed Loop Systems. In this work, a Hammerstein model-based joint identification algorithm proposed by Srinivasan et al.21 is used for detection of stiction in control valves. It has been shown by Srinivasan et 13716



Article

control loop is high.26,31 If the disturbance oscillation occurs at a frequency which is moderately/highly attenuated by the process dynamics then the magnitude of oscillations because of disturbance will be lower than that of the magnitude of oscillations in marginally stable control loop. In general, chemical processes are usually low pass systems and attenuate high frequency components.34 Therefore, disturbances external to the control loop will pass through low pass systems before it manifests as oscillations in the other control loops. However, if the disturbances corrupting the process output have high amplitude because of (i) amplification of low frequency disturbances by low pass filter nature of process dynamics or (ii) high magnitude of disturbance then the proposed algorithm will not be able to correctly identify the cause for oscillations. This is discussed in section 4.3.1. Now, all of these ideas are brought together through a simple example. 4.2.2. Simulation Example. Let us consider a second-order plus dead time (SOPTD) system with the PID controller. The process (Gp(s)) and controller (Gc(s)) transfer functions are given by

model. In other words, when constant set points are present, driven models obtained from closed loop data are not reliable. (2) Frequency domain minimum variance index approaches.26,31 The dominant frequency of oscillation present in the process output is filtered and the minimum variance index is computed for the filtered data. The basic premise here is that, aggressively tuned controller performs poorly at all frequencies. Therefore, if controller is the cause for oscillation, the minimum variance index of the filtered data should be high. These methods are not popular, seldom used and require the knowledge of closed-loop time delay to estimate the minimum variance index. (3) Information and graph theoretic approaches:3 In these approaches, cause and effect analysis of the process signals are identified to obtain a qualitative process model which is then utilized to identify the source of disturbances.12 Transfer entropy is mainly used to identify these qualitative process models. (4) Open loop experiments.31 The control loop under analysis is operated in manual mode for a brief period of time. If the oscillations starts dying down, then the loop is operating in marginally stable conditions. The drawbacks with this approach are (i) manual intervention is required to operate the control loop in manual mode and (ii) this method is cumbersome if the number of isolated control loops are large. (5) Amplitude-based discrimination analysis.31 These methods utilize the fact that control loop operating at marginally stable conditions exhibit higher oscillation amplitude compared to external disturbance.26 The proposed HH spectrum based approach falls under this category.332632 According to linear systems theory, it has been shown that under noise-free conditions, sinusoidal disturbances and marginally stable control loop lead to identical PV and OP signals in control loops.26,31 Therefore, reasonable assumptions have to be made regarding the disturbance and noise corrupting the process output. In this work, the following assumptions are made to develop the RCA algorithm for distinguishing between marginally stable control loop and disturbance caused oscillations: (1) Sensor noise corrupting the process is small and has a broad band spectrum similar to that of white noise. It is assumed that the signal to noise ratio (SNR) in the process output is large. In other words, noise does not drown out the system response. (2) The amplitude of the oscillatory disturbance corrupting the process output is lower compared to that of the oscillations due to marginally stable control loop. The proposed methodology uses this assumption on the disturbance along with the fact that in linear systems, oscillations in a marginally stable control loop occur at a single frequency. Further, it is assumed that chemical processes are low pass systems. However, we do not make any assumptions on the frequency of the disturbances corrupting the process or the critical frequency except that they do not occur at the same frequency (or very close). Remark 1. The first assumption is generally reasonable because if the noise corrupting the system is very large then the oscillations caused by a sticky valve or marginally stable control loop will be buried in the noise. Therefore, a practical first task in this case will be to reduce the sensor noise corrupting the system. Remark 2. The second assumption implies: magnitude of the external disturbance at various frequencies and associated process gains should give rise to a disturbance magnitude (corrupting the PV data) which is lower than that of the oscillations due to marginally stable control loop. Further, it is to be noted that magnitude of oscillations in marginally stable

Gp(s) =

2e−3s 30s + 13s + 1 2

⎛ ⎞ 1 + τds⎟ Gc(s) = Kc⎜1 + τis ⎝ ⎠

(2)

PID controller settings that stabilize the process output are Kc = 1.1, τi = 11, and τd = 0.182. The variance of the white noise corrupting the process is 0.3 (which is 10 times less than that of the oscillatory disturbance added to the process output). PV and OP data for (i) a marginally stable control loop and (ii) a control loop with oscillatory disturbances are generated using this SOPTD model. Data for the latter case is obtained by introducing a sinusoidal disturbance (3 sin (ωt), where ω = 2π( f = 0.08)) in the control loop. The data for the former case is obtained by changing the PID controller settings to Kc = 2.5, τi = 11, and τd = 0.182 which resulted in a marginally stable control loop. All data sets are collected using a sampling frequency (Fs) value of 1. The data generated for these two scenarios are shown in Figure 7. The task now is to distinguish between the oscillations due to marginally stable control loop and external disturbances just from this data with no further information. The first step is to see if any discriminatory information is available in Fourier or

Figure 7. SOPTD system. (a and b) Oscillations due to marginally stable control loop. (c and d) Disturbance caused oscillations. 13717



Article

Figure 8. Spectrum of OP data. (a and c) Fourier spectrum of marginally stable loop and disturbance caused oscillations. (b and d) HH spectrum of marginally stable loop and disturbance caused oscillations.

(noise component as discussed in section 4.2.3) gets concentrated (intra wave modulation) around a particular frequency and (ii) high amplitude present at higher frequency (due to disturbance as discussed in section 4.2.3) is reduced in the HH spectrum of the transformed data. This is inline with the example discussed in section section 3.1. These two features helps in distinguishing between oscillations due to marginally stable control loops and external oscillatory disturbances. The OPind data for both oscillations due to a marginally stable control loop and external disturbance are shown in Figure 10. The FT spectrum and HH spectrum of the OPind

HH spectra of the OP and PV signals. Fourier and HH spectra of the OP signals for the two scenarios are shown in Figure 8. From this figure, it can be seen that the magnitude of oscillation in marginally stable loop is higher compared to that of the disturbance caused oscillation. Further, in case of disturbance caused oscillations, there is a small magnitude at low frequencies (compared to amplitude at disturbance frequency) which is not present in case of the controller caused oscillations. This is shown in both the HHT and FT plots in Figure 8. This information in OP signal is utilized for identifying the source of oscillations. Though, OP signals contain the required information for distinguishing between marginally stable control loops and disturbance caused oscillations, in general it is hard to come up with a threshold in the HHT/FFT of the OP signal because of the presence of noise magnitude at all frequencies. To overcome this situation and extract the required information, a nonlinear transformation and analysis using HHT is used to uncover diagnostic information from just the PV-OP data. The OP data is transformed into what is termed as the OPind (ind for induced) data by passing the OP data through a stiction nonlinearity. This transformation is shown in Figure 9. While

Figure 10. OPind: (a) oscillations due to marginally stable control loop and (b) disturbance caused oscillations.

signals in both cases is shown in Figure 11. Remarkably, notice that the HH spectrum of OPind signal for disturbance caused oscillation shows two distinct frequency regimes. However, the HH spectrum of OPind signal for oscillations due to marginally stable loop does not contain such distinct frequency regimes. This is similar to the example discussed in section 3.1. As discussed earlier, the amplitude of the oscillations due to marginally stable control loop is higher compared to the disturbance caused oscillation resulting in single significant frequency regime compared to two regimes for the disturbance caused oscillations (same scenario as the example in 3.1). Notice that this diagnostic feature is not clearly observable in the Fourier spectrum of the transformed variable. This is because nonlinearities in FT give rise to power at fundamental frequency and multiples of the fundamental frequency

Figure 9. Inducing stiction in the controller output.

many transformations could be used, the authors find that passing the OP data through a stiction nonlinearity is one of the means of identifying the required diagnostic information. Remark 3. One could view this nonlinear transformation as a device for feature extraction.35 As described in section 3.1, nonlinear transformation can be used to vary the amplitude of the signal and introduce modulations/power spread around the fundamental frequency in HH spectrum. The idea behind doing this transformation is (i) the amplitude (present scattered) in the HH spectrum of OP data at low frequencies 13718



Article

Figure 11. Spectrum of OPind data: (a and c) Fourier spectrum of of marginally stable and disturbance caused oscillations and (b and d) HH spectrum of marginally stable and disturbance caused oscillations.

(harmonics). In case of external oscillatory disturbances, it is hard to find whether the power at low frequency is due to fundamental frequency of disturbances or due to noise component. This is further complicated by the fact that (i) noise is present at all frequencies and (ii) disturbances could have multiple frequencies. In the latter case, it is almost impossible to obtain the low frequency diagnostic information from the FT spectrum. However, as discussed in section 3, even when multiple frequency signals are present only oscillation around their fundamental frequencies exist which makes it easy to identify the presence of low frequency components (because of noise as discussed in the section 4.2.3). Thus, a nonlinear transformation of the OP signal is performed to obtain the information in HH spectrum, which is utilized to distinguish between oscillations due to marginally stable control loop and external disturbances. Now, from the example it can be observed that it is possible to distinguish between oscillations due to marginally stable loops and disturbances. However, the following important questions are to be answered: (1) In case of oscillations due to external disturbances, why is there a significant amplitude at lower frequencies compared to the frequency at which maximum amplitude is present? (2) How is this information obtained clearly after passing through a nonlinear valve stiction model? These questions are addressed in the forthcoming paragraphs and the discussions are related to the results obtained from the SOPTD example discussed above. 4.2.3. Frequency Domain Analysis of Closed-Loop System. Here, “qualitative” frequency domain arguments are provided to support the observation that significant amplitude is present at lower frequencies compared to the frequency at which maximum amplitude is present. Consider a linear closed-loop system as shown in Figure 1. The error signal (er(t)) and the controller output u(t) for a linear process are given by

U (s ) =

∫0

U (ω) =

Gc(ω) |R(ω) − (E(ω) + D(ω))| 1 + Gp(ω)Gc(ω) (5)

The set point r(t) is usually constant (even if it is a step change, the steady state data after disappearance of transients is of interest) and therefore R(ω) = 0 for all ω except at ω = 0. Now, eq 5 can be rewritten as |U (ω)| =

Gc(ω) |E(ω) + D(ω)| 1 + Gp(ω)Gc(ω)

(6)

Notice that the negative sign is omitted because it just represents the 180° phase shift due to the feedback. Scenario I: Oscillations in Marginally Stable Control Loop. Consider scenario I (as discussed in the SOPTD example in section 4.2.2), where oscillations are the result of a marginally stable control loop. In this case, external oscillatory disturbances are not present and therefore, the Fourier transform of the controller output from eq 8 is given by |U (ω)| =

Gc′(ω) |E(ω)| 1 + Gp(ω)Gc′(ω)

(7)

Here, Gc′(ω) represents the frequency response of the aggressively tuned controller settings that leads to marginally stable control loop. When the control loop is marginally stable, small perturbations in the system leads to oscillations. The frequency of oscillation (ωcr) in the controller output is obtained by solving the closed loop characteristic equation Gp(ω)Gc′(ω) = −1 (condition for marginal stability). Further, it is important to note that the gain of the controller at marginally

t

gc(τ )er(t − τ ) dτ

(4)

Here R(s), U(s), E(s), and D(s) are the Laplace transform of set point, controller output u(t), broad band noise e(t), and disturbance signal d(t). In the above equation, Gc(s)/(1 + Gp(s) Gc(s)) is the closed loop transfer function between the controller output and other signals. Now, computing the Fourier transform of the controller output by replacing s = jω in eq 4 and considering only the magnitudes

er(t ) = r(t ) − (yp (t ) + e(t ) + d(t )) u(t ) =

Gc(s) [R(s) − E(s) − D(s)] 1 + Gp(s)Gc(s)

(3)

Here, gc(t) is the impulse response of the controller. Now, Laplace transform of the controller output obtained from the closed loop system is given by 13719



Article

Figure 12. (a) Bode magnitude plot SOPTD process. (b) Bode magnitude plot closed loop transfer function in eq 8. (c) Bode magnitude plot of the closed loop transfer function in eq 9.

frequencies in the controller output are (i) amplitude of the disturbance signal and (ii) the weighing function |Gc(ω)/(1 + Gp(ω)Gc(ω)|. It is assumed that amplitude of external oscillatory disturbance amplitude is lesser than the oscillations due to marginally stable control loop. This assumption addresses the first factor on the amplitude of external disturbance. Now, it is required to show that OP data contains significant amplitude at low frequencies compared to the maximum amplitude at disturbance frequency (as observed in the example SOPTD system). This is addressed by analyzing Bode magnitude plots of the process transfer function and closed loop transfer functions (shown in eqs 8 and 9) in case of marginally stable control loop and disturbance caused oscillations. The Bode magnitude plot (second order Pade approximation for the delay term) of the process transfer function is shown in Figure 12a. It can be clearly seen that the SOPTD process is a low pass system (in accordance with the assumption that chemical processes are low pass systems). The Bode magnitude plot of weighting transfer functions given in eqs 8 and 9 are shown in Figures 12b and 12c. It can be clearly seen that the magnitude of disturbance is attenuated by the overall closed loop transfer function while in the case of marginally stable loop, the magnitude of oscillations is higher (not attenuated by overall weighting function). Further, the noise component at low frequency is amplified by the closed loop transfer function which is depicted in Figure 12c. The large magnitude of oscillations due to marginally stable control loop masks the presence of low frequency noise component considerably compared to that of disturbance caused oscillations. This is the primary reason for presence of significant amplitude at frequency lower than the disturbance frequency which has the maximum amplitude in case of disturbance caused oscillations. In summary, (i) in case of disturbance caused oscillations, we have significant amplitude at lower frequencies compared to that of the disturbance frequency, while (Figure 12b) (ii) in case of marginally stable control loop, amplitude is present only at characteristic frequency of closed loop system ωcr. This is

stable condition represents the maximum gain that can be set in the controller without making the process output unstable.36 Therefore, the amplitude of oscillation at frequency ωcr in a marginally stable control loop is high.26 The magnitude of the controller output (at ωcr) in case of oscillations because of a marginally stable control loop is given U (ωcr) =

Gc′(ωcr) |E(ωcr)| 1 + Gp(ωcr)Gc′(ωcr)

(8)

Magnitude of the controller output at other frequencies due to broad band noise component can be neglected as the closed loop gain of the system is higher at the critical frequency. This is the reason why there is maximum power at one frequency in both Fourier and HH spectra of controller output for oscillations caused by marginally stable control loop (Figure 8). Scenario II: Oscillations due to external oscillatory disturbances. Consider scenario II (as discussed in the SOPTD example in section 4.2.2), where oscillations are caused due to external oscillatory disturbance. External oscillatory disturbance could have multiple frequency components. However, for the sake of clarity in explanation, it is currently assumed that the external disturbance has a single frequency of oscillation ωd. Now, the Fourier transform of the controller output is given by |U (ω)| =

Gc(ω) |D(ωd) + E(ω)| 1 + Gp(ω)Gc(ω)

(9)

Here Gc(ω) is the frequency response of the controller that stabilizes the process output. From the above equation, we make a qualitative observation that the controller output has two components: (i) amplitude at disturbance frequency weighed by a transfer function and (ii) amplitude at frequencies due to noise component weighed by the same transfer function as that of the disturbance. The critical factors which decides the significant amplitude (amplitudes above certain threshold value) at various 13720



Article

Figure 13. SOPTD: (a) Sinusoidal signal of frequency ω = 2π0.03 and nonlinearly transformed sinusoidal signal using one parameter valve stiction model (b) HH spectrum of sinusoidal signal (c) HH spectrum of nonlinearly transformed data showing oscillation around fundamental frequency.

clearly illustrated in Figure 12. Thus, the first question raised in end of section 4.2.2 is addressed. Now, the nonlinear transformation on OP data which helps in identifying the diagnostic information is explained. 4.2.4. Nonlinear Transformation of the Controller Output. Though HH spectrum of OP data contains the required diagnostic information it is hard to identify a threshold which can be used in practice to distinguish between marginally stable control loop and external oscillatory disturbances. As explained in section 3.1, nonlinear transformation of a signal results in intrawave modulations/power spread around the fundamental frequency in HH spectrum. For a pure sinusoidal signal, the output of the nonlinear one parameter valve stiction model (used in this work) for one cycle of a sinusoidal input A sin(ωt) with ϕ = ωt (all variables defined in the reference) is given by Srinivasan et al.:37 ⎧−d /2 ⎪ ⎪ d /2 ⎪ ⎪ 3d /2 ⎪ c(t ) = ⎨ d /2 ⎪ ⎪−d /2 ⎪ ⎪−3d /2 ⎪ ⎩−d /2

at the fundamental frequency is reduced after nonlinear transformation (refer Figure 4). This is due to (i) reduction in total power of the transformed signal after nonlinear transformation using one parameter valve model and (ii) spread of power around fundamental power in HH spectrum compared to large power at fundamental frequency for the signal before nonlinear transformation. The HH spectrum of the sinusoidal signal and its nonlinear transformation are shown in Figure 13b and c. Intrawave modulations are present in the HH spectrum (Figure 13c) of the nonlinearly transformed sinusoid using one parameter valve stiction model. Now, this is also true for the nonlinear transformation of oscillatory data caused due to external oscillatory disturbance and marginally stable control loops. For the SOPTD system, the data obtained from this nonlinear transformation of the OP data for the two scenarios are shown in Figures 10a and 10b. The HH spectrum of the corresponding signals shown in Figure 11b and d exhibit intrawave modulations/oscillations around the fundamental frequency. Further, the power at the fundamental frequencies in HH spectrum after nonlinear transformation is reduced compared to that before nonlinear transformation. (refer HH spectrum in Figures 11 and 8). In case of disturbance caused oscillations, significant amplitudes is present at low frequencies (after nonlinear transformation) resulting in power in a single regime at low frequency apart from intrawave modulations at disturbance frequency. This is mainly due to the reasons discussed in sections 4.2.3 and 4.2.4. Thus, HH spectrum of nonlinear transformed controller output OPind reveals the required diagnostic information and addresses the second question raised in the end of section 4.2.2. Remark 4. As mentioned earlier, several kinds of nonlinearities could be used for feature extraction. In this work, control valve stiction type of nonlinearity is used since control valve is the element which succeeds controller. Regarding the amount of stiction, RCA algorithm uses a fixed value of 30% of maximum value of oscillation in the OP data. However, simulation studies (discussed in section 5) are performed by varying the amount of stiction introduced in the OP data. It is

if 0 < ϕ < ϕ1 if ϕ1 ≤ ϕ < ϕ2 if ϕ2 ≤ ϕ < π − ϕ1 if π − ϕ1 ≤ ϕ < π + ϕ1 if π + ϕ1 ≤ ϕ < π + ϕ2 if π + ϕ2 ≤ ϕ < 2π − ϕ1 if 2π − ϕ1 ≤ ϕ < 2π

(10)

The output of the one parameter valve model (d = 30% of maximum amplitude of the signal) for a sinusoidal signal is shown in Figure 13a. From the output of the valve model, the following important observations can be made: (i) output of the one parameter valve model never exceeds the sinusoidal input signal and exhibits stick jump behavior as discussed in ref 21, and therefore, (ii) the power of the nonlinearly transformed signal is less than the power of the sinusoidal signal in both time and frequency domain (relation between power in time and frequency domain established by Parseval’s theorem). As discussed in section 3.1, the magnitude of the oscillatory signal 13721



Article

Figure 14. Power of OPind data: (a) Phht, Disturbance caused oscillation; (b) Phht, Oscillations due to marginally stable control loop.

diagnosis is that marginally stable control loop is the cause for oscillations in the control loop. (5) In practice, for comparison of Phht of OPind and PV signals, a threshold value is required. This is to neglect small power values at various frequencies because of the presence of noise. In this work, the threshold value for Phht of OPind signal is chosen to be 5% of the maximum value of Phht computed from process output (PV). This normalizes the threshold based on the data directly. Only the power values above this threshold are considered to be significant. This threshold can be raised or lowered, if a priori knowledge on noise corrupting the process is known. Remark 5. The normalized frequency range for a system is 0−0.5. In RCA algorithm, a value of 0.46 is chosen as the upper limit to avoid the high frequency white noise from corrupting the HH transform. However, this is not critical and entire frequency region can be considered. Moreover, the sampling frequency of chemical systems is usually 5−10 times greater than that of the fast time constant in the process and therefore, it is unlikely that the propagated disturbance frequency/ marginally stable control loop critical frequency lies within the neglected normalized frequency region (0.46−0.5). The above feature extraction algorithm is implemented on the SOPTD system discussed previously. In case of disturbance caused oscillation, Phht computed for OPind data contain significant values at lower frequencies compared to that of the frequency at which maximum amplitude is present. This can be observed from Figure 14a, where the chosen threshold is plotted as a bold horizontal line. Thus, from steps 3 and 4 of the above algorithm, the cause for oscillation is determined as external disturbance. In the case of marginally stable control loop oscillations, the feature extraction algorithm is implemented and the Phht of OPind is shown in Figure 14b. From the Figure, it can be clearly seen that there is no significant amplitude at lower frequencies compared to the frequency at which maximum amplitude is present. Therefore, using steps 3 and 4 of the algorithm the cause for oscillation is determined to be marginally stable control loop. A flowchart of the proposed RCA algorithm used for diagnosis of oscillations in closed loop systems is shown in Figure 6.

observed that the results are fairly insensitive to the amount of stiction introduced. In general, the percentage of stiction introduced does not paly a role as long as the OP data is nonlinearly transformed leading to an output signal, which has a square wave nature. While until now appeal is made to visual inspection of the HH spectrum plots, a quantification of the diagnostic information is needed for root cause diagnosis. The quantification approach is discussed next. 4.3. Feature Extraction Algorithm for Distinguishing Controller and Disturbance Caused Oscillations. (1) Compute the HH spectrum of the OPind signal. This involves computation of power as a function of frequency and time (please refer to appendix for computation of HH spectrum). (2) Divide the total frequency range (normalized frequency value f = (0−0.5) with ω = 2πf) into various frequency regimes with an interval of 0.02 and starting frequency as 0.02 and ending frequency as 0.46. Now the normalized frequency regimes range from f(i) = 0.02:0.02:0.46 (ω(i) = 2πf(i)). The interval of 0.02 (normalized frequency value F/Fs, where Fs is the sampling frequency) is just 6% in the total frequency scale of 0.5 and therefore does not pose any problems in the identification of root cause of oscillations. (3) Compute the normalized power Phht( f(i)) ( f(i) defined in previous step) from HH power spectrum at various frequency regimes for OPind and PV (used for threshold described in step 4) signals using the following equation: f (i)

Phht(f (i)) =

∑ f = f (i − 1) Phht 2(f ) N

(11)

(4) To distinguish between marginally stable loop and disturbance caused oscillations, Phht of OPind signal is compared using a threshold. If there is a significant amplitude at lower frequencies compared to that of the frequency at which maximum amplitude is present, then the diagnosis is external disturbance induced oscillation. However, if there is no significant amplitude at lower frequencies compared to that of the frequency at which maximum amplitude is present, the 13722



Article

clearly seen from this Figure that significant amplitude is present only at one frequency region. This can also be seen from the plot of Phht( f(i)) shown in Figure 15c. From Figure 15c, it can be observed that there is no significant amplitude at lower frequency compared to that of the frequency at which maximum amplitude is present. This indicates that the cause for oscillations is due to marginally stable control loop. However, the true cause for oscillations in this loop is due to large external disturbance. Now, consider the same SOPTD system (discussed in section 4.2.2) with external oscillating disturbance frequency changed to 0.04 (from 0.08). Since the frequency of oscillation is reduced, the disturbances at low frequencies will be amplified as per the Bode magnitude plot of closed loop transfer function shown in Figure 12c. The PV and OP data obtained from this case is shown in Figure 16a. The proposed RCA algorithm is applied on this process data. Th plot of Phht(f(i)) is shown in Figure 16b. Notice, that the power at the disturbance frequency is much higher compared to the case when frequency is 0.08 (Figures 14a and 16b). This is because, in this case, disturbance (with frequency f = 0.04) is considerably amplified by the closed loop transfer function. From the figure, it can be seen that there is no significant amplitude (above threshold) at lower frequencies compared to that of the frequency (disturbance) at which maximum amplitude is present. RCA predicted the root cause of oscillation to be marginally stable control loop while the true cause is disturbance at frequencies which are amplified by the process transfer function. To highlight the fact that amplitude of disturbance is the important factor, let the amplitude of disturbance be reduced from 3 to 1 while the frequency of oscillation is the same (0.04). The process data obtained from this scenario is shown in Figure 16c. The proposed RCA algorithm is applied to this case and the result (plot of phht( f(i))) is shown in Figure 16d. The amplitude at lower frequency (lower to the frequency (disturbance) at which maximum amplitude is present) is very close to the threshold. Now the amplitude of external disturbance is still decreased (amplitude is 0.3). The process data obtained from this case is

A summary of the results of the above RCA algorithm implemented on the SOPTD system discussed in 4.2 is presented in Table 1. It can be seen from the Table that the Table 1. Results for SOPTD System I G(s) (2e−3s)/(30s2 + 13s + 1)

C(s)

actual case

d

predicted case

Kc = 1.1

stiction

0.4

stiction

τi = 10 τd = 0.182 Kc = 2.5

marginally stable

0

marginally stable

disturbance frequency ( f = 0.08)

0

disturbance

τi = 10 τd = 0.182 Kc = 1.1 τi = 0.1 τd = 0.2

proposed RCA algorithm is able to diagnose the source of oscillations in all the cases. Before discussing the results obtained from simulation and several industrial case studies we provide a discussion on the assumptions (stated earlier at the end of 4.2.1) and their implications. 4.3.1. Discussion on the Assumptions Made in the Proposed Approach. In the proposed approach, it is assumed that the amplitude of the disturbance corrupting the process output is less than the amplitude of the oscillations caused due to controller. However, if this assumption is violated the proposed approach may fail to identify the root cause for oscillations. Consider the SOPTD system discussed in section 4.2.2 with controller settings that stabilizes the process output. The amplitude of the sinusoidal disturbance corrupting the process output is increased to 8 (from magnitude of 3). The PV data obtained from the process is shown in Figure 15a. The HH spectrum of the OPind data is shown in Figure 15b. It can be

Figure 15. Failure case: (a) Process output, (b) HH spectrum of OPind data, and (c) Phht(f(i)). 13723



Article

Figure 16. Failure case: (a, c,and e) Process output and (b, d, and e) Phht( f(i)) for disturbance frequency 0.04 and various amplitudes (3, 1, and 0.3).

shown in Figure 16e. RCA algorithm is applied to this scenario and the result (plot of phht( f(i))) is shown in Figure 16f. Notice that the magnitude at the lower frequency (compared to frequency at which maximum amplitude is present) is clearly visible (compare with Figure 16d). In this case, RCA algorithm identified the root cause as external disturbance. Thus, with reduction in amplitude of external disturbance, RCA algorithm is able to correctly identify the root cause as external disturbance (even if it occurs at various frequencies). This simulation study shows that RCA algorithm will be able to identify the root cause correctly if: amplitude of external disturbance is smaller (than the marginally stable controller caused oscillations) so that the low frequency (compared to that of the external disturbance) feature due to noise component can be identified. The large amplitude of the disturbance signal could be the result of (i) external high amplitude oscillatory disturbances that are not attenuated considerably by controller or (ii) the weighting function (given in eq 9) that amplifies the disturbance signal, which could happen if the disturbance signal occurs at considerably low frequencies (assuming that the process acts as low pass system). In general, controllers are usually tuned to reject disturbances31,36 and therefore is unlikely to amplify the disturbances though it could not attenuate the effect of disturbance considerably. In summary, if amplitude of external disturbance is high, the proposed approach cannot identify the root cause for oscillations. For the SOPTD system, Table 2 provides the results obtained from the proposed approach for external sinusoids of varying amplitudes. Therefore, if the RCA algorithm indicates marginally stable control loop, caution should be exercised as it could be due to high amplitude external disturbances. Results obtained using RCA algorithm on simulation and several industrial case studies are discussed in the following section. In this approach it is assumed that process output is corrupted with broad band noise, which is a true industrial scenario. In case noise is completely absent, the distinguishing features will not be present.

Table 2. Results Obtained from SOPTD System I with Various Disturbance Amplitudes sinusoidal disturbance amplitude 1 2 4 6 7 8

HH spectrum of OPind amplitude regions amplitude regions amplitude regions amplitude regions amplitude regions amplitude regions

RCA algorithm diagnosis result

at 2 frequency

disturbance

at 2 frequency

disturbance

at 2 frequency

disturbance

at 2 frequency

disturbance

at 2 frequency

disturbance

at 1 frequency

marginally stable

5. SIMULATION AND INDUSTRIAL CASE STUDIES While the RCA algorithm was extensively tested on several simulation examples, in this section, results are reported for a representative simulation study and on industrial data sets. As a first example, the proposed RCA algorithm is benchmarked on a SOPTD system that is different from the motivating example discussed in section 4.2 to test the generality of the proposed solution approach. While the algorithm has been tested on higher order systems, chemical process systems are usually well represented by second order time delay models. Using this case study, it is shown that the proposed RCA algorithm is insensitive to the amount of stiction introduced in OP signal to obtain OPind signal. It is also shown that the proposed algorithm can (i) identify different values of stiction, (ii) identify disturbances of varying frequencies corrupting the system, and (iii) identify aggressively tuned controller with various controller settings. Industrial case studies presented in this section demonstrate the utility of the proposed method in practical situations. 5.1. Simulation Studies. To show the versatility of the proposed approach, consider an SOPTD system with PI controller whose transfer functions Gc(s) and Gp(s) are given by 13724


Industrial & Engineering Chemistry Research 1e−5s 50s + 20s + 1

Gc(s) = Gp(s)

Article

the robustness of the proposed RCA approach. Case studies with two different controller tuning related oscillations shows that the proposed algorithm is insensitive to the actual controller and identifies marginally stable control loops. The sum of sine disturbance case shows that the proposed approach can handle a broader range of disturbance characteristics. Remark 6. In simulation studies, marginally stable SOPTD systems are simulated by changing the controller settings. In reality, changes in controller setting and/or process dynamics (time constants gain and delays) leads to marginally stable control loop. However, since the systems considered are linear, change of controller settings (integral and controller gain) will affect the overall closed loop transfer function in a manner similar to that of the change in process gain and time constants that results in marginally stable control loop. Remark 7. On the basis of the assumptions and frequency domain analysis provided in section 4.2, RCA algorithm can identify sources of oscillation as long as the broad band oscillatory disturbance signal at various frequencies does not have amplitudes large enough compared to controller caused oscillation. This ensures that the disturbance signal does not mask the presence of significant amplitudes (above threshold) at frequencies lower to that of the frequency at which disturbance amplitudes are present. 5.2. Industrial Case Studies. In this section, the proposed approach is demonstrated on seven industrial data sets. The first two loops are the data sets (provided by ref 26) from control loops with valve stiction. The last five data sets are obtained from a chemical industry and pulp and paper industry (name not revealed here). In these, two loops contains a sticky valve and oscillations in one loop is due to a marginally stable control loop. Further, cause for oscillations in two control loops are not identified by the industry. The proposed root cause analysis algorithm is applied to all these industrial loops and the results obtained are discussed here. 5.2.1. Flow Control Loops. Results obtained using the RCA algorithm on two flow control loops provided by26 (FC525 and FC145) are discussed. The PV and OP data obtained from the flow loop FC525 is shown in Figure 17a. Hammerstein based approach using one parameter model is applied to the particular data set and the result obtained is shown in Figure 17b. As explained earlier in section 4.1, the minimum value of RMSE

2

⎛ 1⎞ = K c ⎜1 + ⎟ τis ⎠ ⎝

(12)

Data are generated for different scenarios that include sticky valve with d = 0, two controller tunings leading to marginally stable control loop and two disturbance cases, one with a single frequency of oscillation and an another disturbance with two frequencies. The results of the RCA algorithm on these data sets are presented in Table 3. From the Table it can be seen Table 3. Results for SOPTD System II G(s)

C(s)

(1e−5s)/(50s2 + 20s + 1)

Kc = 0.8 τi = 10 Kc = 1.1 τi = 4.16 Kc = 0.8 τi = 3.33 Kc = 0.8 τi = 3.33 Kc = 0.8 τi = 10 Kc = 0.8 τi = 0.1

actual case

d

predicted case

stiction

0.6

stiction

marginally stable

0

marginally stable

marginally stable

0

marginally stable

marginally stable*

0

marginally stable

disturbance frequency ( f = 0.20)

0

disturbance

disturbance frequency ( f = 0.20 + 0.25)

0

disturbance

that all the cases were correctly identified by the RCA algorithm. In the marginally stable case denoted with an asterisk (*) in Table 3, the induced stiction amount was changed to 50% from the default value of 30% that is used in the RCA algorithm. This has no effect on the final result demonstrating

Figure 17. Flow loop−FC525: (a) process data and (band c) results from one and two parameter model based Hammerstein approach. 13725



Article

Figure 18. Flow loop−FC145: (a and b) Process data and (c and d) results from one and two parameter model-based Hammerstein approach.

Figure 19. Pressure loops: (a) Process data, loop I, (b) results from one parameter model-based Hammerstein approach, Loop I, (c) process data, loop II, and (d) results from one parameter model-based Hammerstein approach, loop II.

The PV and OP data obtained from the flow loop FC145 is shown in Figure 18a and b. A Hammerstein-based approach using one parameter model is applied to the particular data set, and the result obtained is shown in Figure 18b. Results indicate that stiction is present in this particular loop (d = 0.32). Hammerstein approach with two parameter valve model also indicates the presence of stiction (S = 0.5 and J = 0.5). This result is shown in Figure 18c. 5.3. Pressure Control Loops. In this section, the results of RCA algorithm on two pressure control loops are presented. The loops are found to contain intermediate and high sticky control valves. The PV and OP data obtained from the first loop is shown in Figure 19a with the corresponding Hammerstein-based approach result in Figure 19b. Results

obtained after model fitting occurs at a nonzero value of d indicating the presence of stiction (stiction band value d = 2.16). Results obtained from Hammerstein approach using two parameter valve model29 is shown in Figure 17c. In the approach provided in Choudhury et al.,29 instead of one parameter valve model, a two parameters valve model is used. Here both one and two parameter models are used for stiction detection. The two parameter model based Hammerstein approach also provides a nonzero value for the S and J (2 and 0.5) parameters indicating the presence of stiction. The stiction detection results did not depend on the use of one or two parameter valve models. According to the algorithm depicted in Figure 6, the root cause for oscillations in this loop is successfully identified as control valve stiction. 13726



Article

Figure 20. Pulp and paper industry loop 1: (a) Process data, (b) results from one parameter model-based Hammerstein approach, (c) HH spectrum of OPind data, and (d) Phht(f(i)) for OPind data.

Figure 21. Pulp and paper industry loop 2: (a) Process data, (b) results from one parameter model-based Hammerstein approach, (c) HH spectrum of OPind data, and (d) Phht(f(i)) for OPind data.

industry control loop is shown in Figure 20a. Hammerstein based approach results for this particular loop is shown in Figure 20b. It can be seen that the minimum value of RMSE (best model) has a corresponding d value of zero indicating that there is no stiction in this control loop. According to the root cause analysis (RCA) algorithm approach provided in Figure 6, OPind data is obtained. HH spectrum of OPind is shown in Figure 20c. From this figure, it can be clearly seen that there is no significant amplitude at lower frequencies compared to the frequency at which maximum amplitude is present. This can be clearly noticed from the plot of Phht(f(i)) shown in Figure 20d. Thus RCA algorithm concluded that the root cause

indicate that minimum RMSE value is achieved at a stiction band value greater than zero indicating the presence of sticky valve. PV and OP data set obtained from the second pressure loop is shown in Figure 19c. From Figure 19d, it can be seen that d value is not zero for the model corresponding to the minimum RMSE. This clearly indicates the presence of stiction in this control loop. 5.4. Data Sets from Pulp and Paper and Chemical Plants. In this section, results obtained from application of RCA algorithm on three industrial data sets are discussed. The plot of PV and OP data obtained from a pulp and paper 13727



Article

Figure 22. Chemical plant control loop: (a) Process data, (b) results from one parameter model-based Hammerstein approach, (c) FT spectrum of OP data, and (d) HH spectrum of OPind data.

Figure 23. Chemical plant control loop: Phht( f(i)).

concluded that the root cause for oscillations in this loop is due to marginally stable control loop. Data set from a chemical industry is shown in Figure 22a. After it was tested by the industry, valve stiction was ruled out. This is corroborated by our analysis as clearly seen from Figure 22b, as minimum RMSE occurs at d = 0. The Fourier transform of the OP data shown in Figure 22c indicates that there is significant amplitude at low frequencies compared to that of the frequency at which maximum amplitude is present. This information can be also seen from the HH of OPind data shown in Figure 22d. The results obtained from applying the proposed feature extraction algorithm is shown in Figure 23. From the figure, it can be clearly seen that there is a presence of significant amplitude at low frequency compared to that of the frequency at which maximum amplitude is present. RCA

for oscillations in this loop is due to marginally stable control loop. In fact as per the industry, the cause for oscillations in this loop is improperly tuned controller (cause for marginally stable control loop). Data set from another pulp and paper industry control loop is shown in Figure 21a. In this case, the true cause for oscillations is informed to be other than valve stiction. However, the root caused is not revealed by the industry. Hammerstein based approach identified that there is no stiction in the control loop (d = 0). This can be clearly seen from Figure 21b, as minimum RMSE occurs at d = 0. Thus, the proposed algorithm correctly identified that the root cause for oscillations is not due to stiction. HH spectrum of OPind data shown in Figure 21c indicates that there is significant amplitude at only one frequency region. This can also be seen clearly from the plot of Phht(f(i)) shown in Figure 21d. RCA algorithm 13728



Article

Figure 24. SOPTD system: (a) Process data with drift, (b)Phht( f(i)) of OPind data, (c) process data with nonstationarities, and (d) HH spectrum of OP data with nonstationarities.

0.25 corrupting the output of the closed loop SOPTD process discussed in section 4.2.2. The PV and OP data for the system corrupted with the chirp disturbance is shown in Figure 24c. The HH spectrum of the OP data is shown in Figure 24d. From this figure, it can be seen clearly that the disturbance contains frequencies which varies with respect to time. Thus, one can identify the presence of nonstationary disturbances from the HH spectrum of the data. This can be incorporated in the current RCA algorithm with little modifications. In the presence of control loop interactions, the proposed root cause analysis algorithm would be able to identify the causes for oscillations under the following conditions: (i) presence of stiction in one control loop does not exhibit nonlinearity in the other interacting control loop and (ii) oscillations due to a marginally stable loop in one of the interacting loops does not cause oscillations with large amplitude in the other loops. The former condition is because, if nonlinearity is present in both the control loops, stiction identification algorithm will indicate stiction in both control loops. This condition is likely to be met in industrial scenarios because most chemical processes act as low pass systems and attenuate the high frequency harmonic components in the signal. The latter condition helps the root cause analysis algorithm to identify the presence of marginally stable control loop in one of the interacting loops. Future work is directed toward a detailed study of the proposed root cause analysis algorithm for interacting control loops.

algorithm concluded that the root cause for oscillations in this control loop is due to external disturbances. While until now, we have discussed the use of RCA algorithm for linear oscillatory disturbances, a discussion of the proposed approach to handle drifts, nonstationary oscillatory disturbances, and interaction among control loops is discussed next. 5.5. Drifts, Nonstationary Disturbances and Interacting Control Loops. The EMD step in HHT is an adaptive data processing technique, which extracts (i) the oscillating components at various scales as IMFs and (ii) the trend component of the data as the last IMF (discarded during computation of HHT as discussed in appendix). This feature of EMD (component of HHT) enables the proposed RCA algorithm to identify the source for oscillations even in the presence of drifts (could be due to sensors, slow changes in process dynamics, etc.) in the control loop. The ability of the EMD to handle drifts is discussed in detail with experimental results in Subramanian et al.38 Consider the closed loop SOPTD system whose output is oscillatory due to oscillatory external disturbances (section 4.2.2). Now the disturbance signal is modified to have a time varying drift component. The disturbance signal is d(t) = 3 sin(2π0.08t) + n2 where n = 0.001:0.1. The process data obtained in this case is shown in Figure 24a. RCA algorithm is applied to this data set. The plot of Phht( f(i)) obtained from RCA algorithm is shown in Figure 24b. It can be clearly seen that there is a significant amplitude at lower frequencies compared to that of the frequency at which maximum amplitude is present. As per the algorithm, the cause for oscillations in this control loop is due to external disturbance. This example demonstrates the ability of RCA algorithm to identify source of oscillations in presence of control loops with drifts/time varying trends. HHT is a time-frequency analysis technique which is primarily used to analyze nonstationary and nonlinear signals. Therefore HHT can also be used to identify nonstationary disturbances corrupting the process output. For instance, consider a chirp signal28,39 with frequency range from 0.15 to

6. CONCLUSIONS A novel and robust method for diagnosis of cause of oscillation in closed-loop systems is developed. The proposed method combines both the parametric and nonparametric techniques for root cause analysis of oscillatory systems. The advantages of the proposed method are (i) a unique signature for distinguishing between controller and external disturbance is developed, (ii) nonstationary nature of the disturbance can be naturally handled since HHT is a time-frequency analysis tool, 13729



Article

where 9 is the real part of a complex number. This representation of the IMF is obtained by computing analytic signal zk(t) from the kth IMF.20 The analytic signal zk(t) is given by

and (iii) no assumptions on noise structure (we assume only the presence of white or broad band noise) is necessary. There are only two tuning parameters, namely, (i) introduction of 30% stiction in OP data and (ii) threshold value (5% of maximum power value in PV). These two values are maintained as the same for all the validation results that include simulation and industrial data. The results obtained from simulation and industrial data demonstrate the utility of the proposed method for root cause analysis. Future work will focus on the enhancement of the current algorithm to identify multiple causes of oscillation in closed-loop systems.

zk(t ) = dk(t ) + iH(dk(t )) = Ak (t )exp(iθ(t ))

with amplitude Ak(t) = (dk(t)2 + H(dk(t)2))1/2 and phase θ(t) = arctan H(dk(t))/dk(t) = ωt. In eq 15, H(dk(t)) is the Hilbert transform of the kth IMF, given by

■

H(dk(t )) =

APPENDIX A: A HILBERT−HUANG TRANSFORM Hilbert−Huang Transform (HHT) is a relatively new method of analyzing nonlinear and nonstationary data.28,39,40 HHT is a data-based approach constituting of the following two steps: (i) adaptively decomposing the signal x(t) into bases called intrinsic mode functions (IMF) through a procedure called empirical mode decomposition (EMD) and (ii) computation of analytic signal from IMFs to obtain the HH spectrum.28,41,42

ωk(t ) =

(16)

d θ (t ) dt

(17)

N k=1

(18)

where ωk is the instantaneous frequency given by dθ/dt with θ being θ(t) = arctan H(dk(t))/dk(t). Equation 1828 enables us to represent the amplitude and the instantaneous frequency as functions of time in a 3-D plot, in which the amplitude can be contoured over the frequency-time plane as the z-axis surface or by using a color scale in the 2-D frequency-time plane to represent amplitude. This plot of a superposition of all ωk(t) in the frequency-time plane is called as Hilbert spectrum H(ω,t) obtained from EMD of the signal (Hilbert−Huang spectrum). Hilbert−Huang transform (HHT) because of its ability to derive the basis functions from the data can handle nonlinear as well as non-stationary data. The Fourier domain representation for the signal x(t) in terms of time-independent amplitude AK and frequency ωK is given by M

x(t ) =

∑ Ak eiω t k

k=1

(19)

In contrast, the time-dependent amplitude A k (t) and instantaneous frequency ωk(t) representation in eq 18 allows to analyze data with amplitude and time-varying frequency (e.g. non-stationary data). The ability of HHT to handle nonstationarities and uniquely represent nonlinearities is well established in the work in ref 28.

■

AUTHOR INFORMATION

Corresponding Author

(13)

*E-mail: [email protected].

where dk(t) is the kth intrinsic mode function and mN(t) is the residual trend obtained after extraction of N IMFs. Subsequently, Hilbert Transform is applied to the N IMFs. The residual mN(t) is ignored since it represents the monotonic or constant data trend.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge partial support from NSF grant CBET-0934348 titled "GOALI: An Integrated Framework for Stiction Detection and Compensation in Control Loops" and a subsequent NSF-TRAC supplement for the project.

A.0.2. Computation of Analytic Signal

The kth IMF obtained from the signal can be represented as

dk(t ) = 9(Ak (t )eiωt )

d (t ′)

x(t ) = 9 ∑ Ak (t )e(i ∫ ωk(t )dt )

N k=1

∞

∫−∞ t k− t′ d(t′)

Now, the signal x(t) is given by

A simplified description of the EMD is provided here. EMD algorithm obtains the IMFs by decomposing the given signal x(t) into local high-frequency d(t) (detail) and local lowfrequency components m(t).20,28 The individual IMFs which are obtained by the EMD algorithm should satisfy the following two criteria: (i) the number of extrema and zero crossings must be either equal or differ at most by one and, (ii) at any instant in time, the mean value of the envelopes defined by the local maxima and the local minima is zero. Effectively the EMD algorithm can be outlined as follows:28 (1) Identify all extrema of signal x(t). (2) Create “envelope” of extrema by performing a cubic spline fit of the extrema points. (3) Compute the average of the “envelopes” m(t) = (emin(t) + emax(t))/2. (4) Compute the signal h(t) which is given by h(t) = x(t) − m(t). If h(t) satisfies the conditions of an IMF, proceed to next step. Otherwise, assign x(t) = h(t) and iterate steps 1−4. (5) Assign h(t) to dk(t), the kth IMF. (6) Obtain residual r(t) = x(t) − dk(t). (7) If r(t) contains at least two extrema (one maximum and one minimum) set x(t) = r(t) and iterate (steps 1−7); otherwise assign r(t) to mN(t), the data trend and terminate the decomposition process. The sifting process, iterating steps 1−4 upon detail d(t), is performed until the detail can be considered as zero mean according to stopping criterion discussed in refs 28 and 39. At the end of the sifting process, detail d(t) is the IMF and the residual is computed as in step 7. EMD operation on a signal x(t) is given by

∑ dk(t )

1 P π

with P as the Cauchy principal value of the integral. The instantaneous frequency of the kth IMF is defined as

A.0.1. Empirical Mode Decomposition (EMD)

x(t ) = mN (t ) +

(15)

(14) 13730



■

Article

(26) Horch, A. Condition Monitoring of Control Loops; Department of Signals, Sensors and Systems, Royal Institute of Technology: Stockholm, Sweden, 2000. (27) Babji, S.; Nallasivam, U.; Rengaswamy, R. Diagnosis of root cause for oscillations in closedloop chemical process systems. Proceeding of the IFAC Conference, Italy 2011; International Federation of Automatic Control: Laxenburg, Austria, 2011; pp 13145−13150. (28) Huang, N.; Shen, Z.; Long, S.; Wu, M.; Shih, H.; Zheng, Q.; Yen, N.; Tung, C.; Liu, H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary timeseries analysis. Proc. R. Soc. London, Ser. A 1996, 454, 903−995. (29) Choudhury, M. A. A. S.; Jain, M.; Shah, S. L. Stiction Definition, modelling, detection and quantification. J. Process Control 2008, 18, 232−243. (30) Kano, M.; Maruta, H.; Kugemoto, H.; Shimizu, K. Practical model and detection algorithm for valve stiction. IFAC Symp. Dyn. Control Process Syst. 2004, 18. (31) Ordys, A.; Uduehi, D.; Johnson, M. Process Control Performance Assessment: From Theory to Implementation; Springer Verlag: London, 2007. (32) Ljung, L. System Identification: Theory for the User, 2nd ed.; Prentice Hall: Englewood Cliffs, NJ, 1999. (33) Bauer, M.; Cox, J.; Caveness, M.; Downs, J.; Thornhill, N. Finding the direction of disturbance propagation in a chemical process using transfer entropy. IEEE Trans. Control Syst. Technol. 2007, 15, 12−21. (34) Stephanopoulos, G. Chemical Process Control: An Introduction to Theory and Practice, 3rd ed.; Prentice Hall: Englewood Cliffs, NJ, 1984. (35) Venkatasubramanian, V.; Rengaswamy, R.; Yin, K.; Kavuri, S. A review of process fault detection and diagnosisPart I: Quantitative model-based methods. Comput. Chem. Eng. 2003, 27, 293−311. (36) Nise, N. Control Systems Engineering, 3rd ed.; John Wiley and Sons: USA, 2000. (37) Srinivasan, R. Peformance monitoring of control loops: stiction detection, estimation and compensation. Phd Dissertation. Department of Chemical Engineering, Clarkson University, Potsdam, NY, USA, 2005. (38) Subramanian, S.; Babji, S.; Zhao, J.; Rengaswamy, R.; McKenna, G. Application of empirical mode decomposition in the field of polymer physics. J. Polym. Sci., Part B: Polym. Phys. 2011, 49, 277−290. (39) Huang, N.; S.S.P. Shen, E. Hilbert−Huang Transforms: Introduction and Applications; World Scientific: Singapore, 2004. (40) Wu, Z.; Huang, N. A study of the characteristcs of white noise using Empirical Mode Decomposition method. Proc. R. Soc. London, Ser. A 2004, 460, 1597−1611. (41) Rilling, G.; Flandrin, P.; Goncalves, P. Empirical Mode Decomposition, Fractional Gaussian noise and Hurst exponent estimation. Int. Conf. Acoustics, Speech, Signal Process. 2005, 489−492. (42) Flandrin, P.; Rilling, G.; Goncalves, P. Empirical Mode Decomposition as a filter bank. IEEE Signal Process. Lett. 2004, 1, 112−114.

REFERENCES

(1) Choudhury, M. Detection and diagnosis of control loop nonlinearities, valve stiction and data comprression management. Ph.D. thesis, University of Alberta, Edmonton, Alberta, Canada, 2004. (2) Thornhill, N.; Shah, S.; Huang, B.; Vishnubhotla, A. Spectral principal component analysis of dynamic process data. Control Eng. Pract. 2002, 10, 833−846. (3) Thornhill, N.; Horch, A. Advances and new directions in plantwide disturbance detection and diagnosis. Control Eng. Pract. 2007, 15, 1196−1206. (4) Thornhill, N.; Shah, S.; Huang, B. Detection and diagnosis of unit wide oscillations. Process Control Instrum. 2000, 26 (July), 26−28. (5) Tangirala, A.; Shah, S.; Thornhill, N. PSCMAP: A new tool for plant-wide oscillation detection. J. Process Control 2005, 15, 931−941. (6) Thornhill, N.; Hagglund, T. Detection and diagnosis of oscillation in control loops. Control Eng. Pract. 1997, 5, 1343−1354. (7) Miao, T.; Seborg, D. Automatic detection of excessive feedback oscillatory loops. IEEE Int. Conf. Control Appl. 1999, 1, 359−364. (8) Thornhill, N.; Huang, B.; Zhang, H. Detection of multiple oscillations in control loops. J. Process Control 2003, 13, 91−100. (9) Jiang, H.; Choudhury, M.; Shah, S.; Cox, J.; Paulonis, M. Detection and diagnosis of plantwide oscillations from industrial data using the spectral envelope method. J. Process Control 2007, 17, 143− 155. (10) Tangirala, A.; Kanodia, J.; Shah, S. Nonnegative matrix factorization for detection and diagnosis of plantwide oscillations. Ind. Eng. Chem. Res. 2007, 46, 801−817. (11) Xia, C.; Howell, J.; Thornhill, N. Detecting and isolating multiple plant-wide oscillatiion via spectral independent component analysis. Automatica 2005, 41, 2067−2075. (12) Jiang, H.; Patwardhan, R.; Shah, S. Root cause diagnosis of plant-wide oscillations using the concept of adjacency matrix. J. Process Control 2009, 19, 1347−1354. (13) Yang, J.; Clarke, D. The self-validating actuator. Control Eng. Pract. 1999, 7, 249−260. (14) Ruel, M. Stiction the hidden menace. Control Mag. 2000, 13, 69−76. (15) Choudhury, M. A. A. S.; Shah, S. L.; Thornhill, N. F. Diagnosis of poor control loop performance using higher order statistics. Automatica 2004, 40 (10), 1719−1728. (16) Horch, A. a. A simple method for detection of stiction in control valves. Control Eng. Pract. 1998, 7, 1221−1231. (17) Rengaswamy, R.; Hägglund, T.; Venkatasubramanian, V. A qualitative shape analysis formalism for monitoring control loop performance. Eng. Appl. Artif. Intell. 2001, 14, 23−33. (18) Singhal, A.; Salsbury, T. I. A simple method for detecting valve stiction in oscillating control loops. J. Process Control 2005, 15, 371− 382. (19) Jelali, M.; Huang, B. Detection and Diagnosis of Stiction in Control Loops; Springer: London, 2009. (20) Babji, S.; Gorai, P.; Tangirala, A. Detection and quantification of control valve non-linearities using Hilbert−Huang transform. J. Adv. Adapt. Data Anal. 2009, 1, 425−446. (21) Srinivasan, R.; Rengaswamy, R.; Narasimhan, S.; Miller, R. M. Control loop performance assessment 2: Hammerstein model approach for stiction diagnosis. Ind. Eng. Chem. Res. 2005, 44, 6719−6728. (22) Lee, K.; Ren, Z.; Huang, B. Novel Closed-Loop Stiction Detection and Quantification Method via System Identification; ADCONIP: Jasper, Canada, 2008, (23) Jelali, M. Estimation of valve stiction in control loops using separable least-squares and global search algorithms. J. Process Control 2008, 18, 632−642. (24) Ivan, L.; Lakshminarayanan, S. A new unified approach to valve stiction quantification and compensation. Ind. Eng. Chem. Res. 2009, 7, 3474−3483. (25) Karra, S.; Karim, M. N. Comprehensive methodology for detection and diagnosis of oscillatory control loops. Control Eng. Pract. 2009, 7, 939−956. 13731


Root Cause Analysis of Linear Closed-Loop Oscillatory Chemical

Recommend Documents