Detecting and Analyzing Nonlinearity-Caused Oscillations in Process Control Systems Using an Improved VNCMD

Nonlinearity-caused oscillations are a frequent issue in process control systems. Its incidence degrades the product quality, stability and safety of the plant. Therefore, it is important to detect and analyze the nonlinearity-caused oscillations to maintain the control performance. In this study, we propose a novel oscillation detection and analyze method based on an improved variational nonlinear chirp mode decomposition (VNCMD) algorithm. Specifically, the original VNCMD needs to manually set the mode number in advance, which is a challenging task in practice. To tackle this problem, an improved VNCMD (IVNCMD) is proposed by utilizing the approximate entropy of instantaneous frequency. Then a novel IVNCMD-based detector is developed to detect and analyze the nonlinearity-induced oscillations by revealing the harmonic content of process variable. Besides detecting the nonlinearity problem, the IVNCMD-based method can contribute in locating the root cause for nonlinearity-caused unit-wide oscillations. The proposed method is model-free and data-driven thus requiring no prior knowledge about the process dynamics. Compared with the latest related methods, the proposed method is able to process nonstationary oscillations and providing corresponding time-frequency information. The effectiveness and advantages are demonstrated through simulations as well as industrial applications.


I. INTRODUCTION
Oscillations are one of the most common problems in process control systems [1]. Oscillation behavior will seriously affect the control system performance, such as reducing the product quality and uniformity, increasing the energy consumption, wasting the raw materials and even causing the plant down [2]. It is reported that over 30% control loops are oscillating because of the presence of nonlinearity problems [3], [4]. The nonlinearity-induced oscillations can be confused with other causes of malfunction, such as poor controller tuning, multi loop interactions. In addition, this kind of oscillation cannot be completely eliminated by controller tuning or by the action of digital valve The associate editor coordinating the review of this manuscript and approving it for publication was Juntao Fei . positioners [5]. Therefore, it is essential to detect and analyze the nonlinearity-caused oscillations [6].
With the development of research in control performance assessment, many methods have been proposed to detect and analyze the nonlinearity-induced oscillations. Herein, a brief overview is provided as follows.
Horch [7] first presented an automatic approach for detecting nonlinearity-caused oscillations based on the cross-correlation function, which utilized the relationship between the process variable (PV) and controller output (OP). However, this method is limited to deal with non-integrating process. To address this limitation, Horch [8] utilized the probability distribution of the second-order derivative of the controlled variable to detect the nonlinearities in integrating processes. Because the valve position (MV) does not change even though OP changes if nonlinearity problems occur in control valves, it is effective to detect the  ApEn value versus mode number. The black dotted line represents the threshold. It is observed that when mode number is larger than 3, the ApEn value will be higher than the threshold. According to Algorithm 1, the mode number should be set as 3. nonlinearity through investigating the behaviour of the MV against the OP. For example, the sticking-valve positions-based [9], parallelogram-based [10], and qualitative analysis-based [11] methods are successively developed by utilizing the relationship between the OP and the Choudhury et al. [12] used an ellipse to fit the filtered PV and OP signals to detect the nonlinearity. Thornhill [13] formulated a statistic based on a three-sigma test, which considered the prediction errors for the surrogates as a reference probability distribution. Choudhury et al. [14] developed a bicoherence-based metric of nonlinearity measure to distinguish the linear and nonlinear oscillations. Besides, the pertaining approaches include the area peak [15], curve fitting [16], bihocerence [17], and so forth. However, these traditional methods are only applicable to stationary process, thus restricting their utility. Recently, signal decomposition methods have been widely introduced into detecting and analyzing oscillations [18]. These signal decomposition-based techniques are able to process signals contaminated by nonstationarity and noise. Although empirical mode decomposition (EMD) [19], [20], intrinsic time-scale decomposition [21], [22], and local mean decomposition [23] have been used for detecting oscillations, these methods do not diagnose the oscillation type. Babji et al. [24] applied the EMD to detect the nonlinearity-induced oscillations. Nevertheless, this method only provided the qualitative description. Aftab et al. [25] further presented an EMD-based measure to quantify the severity of nonlinearity. But it is subjected to mode-mixing and end-effect. To remedy these issues, Aftab et al. [26] utilized the dyadic filter bank property of multivariate EMD (MEMD) to automatically detect the nonlinearity-induced plant-wide oscillations [27]. The above signal decomposition methods are empirically established thus lacking theoretical basis.  Variational mode decomposition (VMD) [28] pioneered the optimization-based signal decomposition method, which has been a research hotspot. Chen et al. [2] proposed an improved VMD to detect the nonlinearity-induced oscillations. The VMD-based method displays satisfactory performance in noise robustness and anti-mode-mixing ability [28]- [30]. However, it is only suitable to analyzing the time-invariant oscillations and cannot provide the time-frequency information. Recently, Chen et al. [31] presented a novel time-frequency analysis approach, termed as variational nonlinear chirp mode decomposition (VNCMD). VNCMD is inspired by the principle of VMD but is able to process nonstationary, time-varying, and wide-band signals [32]. It first demodulates the wide-band signal into the narrow-band signal through demodulation techniques. Then, a variational objective function is established by estimating the bandwidth of demodulated signals. In the end, a joint-optimization scheme based on alternating direction method of multipliers (ADMM) is developed to accurately capture signal modes and the corresponding instantaneous frequencies. VNCMD has shown powerful effect in early fault-detection and multi-feature extraction [31]. However, its requires users to provide the mode number in advance, which is a challenging task in practice [33]. At present, there are few reports on this issue.
To tackle issue, this paper proposes an improved VNCMD (IVNCMD) algorithm by utilizing the approximate entropy of instantaneous frequency. The proposed IVNCMD is able to adaptively determine the mode number and providing the decomposition modes with instantaneous frequencies.
Based on the decomposition results of IVNCMD, a novel oscillation detector is presented by combining the normalized cross-correlation coefficient and regularity index. The cross-correlation coefficient is used to discard the spurious VOLUME 9, 2021 FIGURE 5. The decomposition performance of the original NCMD with K = 2. It is clear that the decomposition results suffer from mode-mixing issue, which is due to the fact that the total mode number is set too small.  modes and regularity index can measure the oscillation degree. Inspired by the fact that the nonlinearity-caused oscillations contain higher order harmonics [26], [34], thus the presence of harmonics can be used as an indicator of nonlinearity problems. Then, a nonlinearity-caused oscillation analysis algorithm is proposed through calculating the frequency relationship between modes. The contributions of this work are as follows,  (i) An improved VNCMD (IVNCMD) is proposed to adaptively determine the mode number parameter; (ii) An IVNCMD-based oscillation detector is developed by combing the normalized correlation coefficient and regularity index to capture the oscillating components of process variables; (iii) A novel nonlinearity-caused oscillation analysis approach is developed through studying the frequency relationship between modes obtained from IVNCMD; (iv) The proposed method is data-driven and can automatically detect and analyze the oscillation type for process control systems. Compared with the VMD-based [2] and EMD-based [26] methods, the proposed method can process both time-invariant and time-varying data, and improve the detection reliability and accuracy.
The rest of this paper is organized as follows. The original VNCMD is briefly described in Section II.
Section III elaborately introduces the proposed improved VNCMD and its application in detecting nonlinearity-caused oscillations. Simulations and industrial cases are provided in Section IV and V, respectively, followed by conclusions.

II. OVERVIEW OF VNCMD
VNCMD [31] assumes the nonstationary signal g (t) consists of several nonlinear chirp modes, expressed as where u k (t) = a k (t) cos 2π t 0 f k (τ ) dτ + ϕ k represents the nonlinear chirp mode. a (t) and f (t) stands for the VOLUME 9, 2021 TABLE 2. The detection and analysis results of the four-input-four-output control system based on IVNCMD.

FIGURE 9.
The structure diagram of a four-input-four-output control system under consideration. The nonlinearity unit is a valve stiction model [42].
instantaneous amplitude and instantaneous frequency, respectively. Both a (t) and f (t) are assumed to be smooth function and vary much slower than the phase function. ϕ is the initial phase. e typifies the decomposition error. Through demodulation techniques, (1) can be reformulated as wheref k stands for the estimated instantaneous frequency.
x k (t) and y k (t) are two demodulated parts, expressed as, The demodulated signals are assumed to be smooth function, whose bandwidth are relative narrow. There a variational objective function can be established by minimizing the bandwidth of x k (t) and y k (t) as where · 2 2 denotes the square of l 2 -norm, which is used to estimate the bandwidth. ε is an upper bound determined by the noise level. (5) can be effectively solved by the alternating direction method of multipliers, thus extracting the desired modes and instantaneous frequencies. The detailed solution procedures are provided in [31].
4: For a given F k (i), count the number for measures the frequency of patterns similar to the one given by the window d within a tolerance r, expressed as 6: Calculate the natural logarithm of each C d r (i) and average them over i, shown as 7: Let d = d + 1 and repeat step 2-6 to get the d+1 (r), 8: The ApEn value is calculated as ApEn = d+1 (r) − d (r).

III. PROPOSED METHOD
VNCMD is able to process both time-invariant and time-varying signals, and it has been successfully applied to mechanical fault diagnosis [31]. However, the original VNCMD requires users to know the mode number K in advance, which is a challenging task in practice. In order to tackle this issue, an improved VNCMD (IVNCMD) algorithm is developed in this section using the approximate entropy (AE) of instantaneous frequency. IVNCMD can adaptively determine the mode number K and decompose the process variable (PV) into a series of modes. It is effective to detect and analyze the process oscillations through investigating the oscillation information contained in these modes. Thus a novel nonlinearity-caused oscillation framework is developed based on the IVNCMD in the end.

A. APPROXIMATE ENTROPY
Approximate entropy (ApEn) can quantitatively describe the complexity of time series by calculating the marginal probability distribution [35]. It has been widely used to measure the randomness and irregularity of sequences. The general procedures for estimating the ApEn are listed in Algorithm 1, where d = 2 and r = 0.2 is recommended [36]. It is well-established that a smaller ApEn value corresponds to a more regular time series; while a larger ApEn value indicates the sequence is random and irregular. With respect to the instantaneous frequency, if it is extracted from the significant mode, its ApEn value should be small, because VNCMD assume the instantaneous frequency is smooth function; on the contrary, if the instantaneous frequency sequence is obtained from the noise, the corresponding curve must be irregular thus shows a relative larger ApEn value.

B. IMPROVED VNCMD
Because the ApEn value of instantaneous frequency can indicate whether this instantaneous frequency is extracted from the significant mode, it is reasonable to use this index to determine the mode number of VNCMD. More specifically, we first set the mode number K as 1 and apply the VNCMD to the input signal. Then, calculate the ApEn value of instantaneous frequency of the extracted mode. If the ApEn value is smaller than the predefined threshold, it is claimed that this mode is significant and repeat the VNCMD with K = k + 1; On the contrary, a larger ApEn value indicates this mode is likely to be obtained from noise, which means the mode number setting exceeds the actual value. The detailed procedures of the proposed IVNCMD is provided in Algorithm 2.

C. DETECTING AND ANALYZING NONLINEARITY-CAUSED OSCILLATIONS
The proposed IVNCMD is able to extract K modes from the input signal. The next procedure is to confirm whether these modes are oscillating components. It is reasonable to VOLUME 9, 2021 identify the significant modes before investigating oscillation behaviors, because IVNCMD may produce pseudo modes. According to the definition in (1), the generated modes are orthogonal, thus the significant mode will be strongly correlated with the input signal; while the correlation relationship between pseudo components with input signal should be weak. Herein, the cross-correlation coefficient [37], [38] is adopted to quantify the relationship, shown as for k = 1, 2, . . . , K , where Cov (·) and σ · stand for the covariance and standard deviation, respectively. For convenience, the cross-correlation coefficient is normalized using for k = 1, 2, . . . , K . Generally, if α k > α µ , the corresponding mode u k can be distinguished as significant component and retained for the following analysis, where α µ =0.15 is suggested by [26].
According to the definition in [3], oscillation should show a periodic variation that is not completely hidden in noise. Therefore, the regularity statistic is adopted to measure the periodic variation [39]. Because the auto-covariance function (ACF) displays the same oscillating period as the original oscillation signal, but is not influenced by noise, the significant modes are converted to the ACF. Denoting t as the time interval between two successive zero crossings of ACF, the average time period is calculated as where Z typifies the interval number and is set as 10. The regularity statistic is determined using the relationship, where σ P k represents the standard deviation of the time intervals between zero crossings. The modes with β k > β µ = 1 are regarded to be oscillating [26].
To sum up, the oscillations are reported if the decomposition modes satisfy the condition α k > α µ and β k > β µ .
When the oscillation is detected, we can analyze whether it is caused by nonlinearity problems by investigating the frequency relationship between modes. More specifically, the presence of nonlinearity-caused oscillations can be indicated by identifying the higher order harmonics [26]. Algorithm 3 is developed to automatically capture the harmonics to identify the oscillation type. The complete workflow for detecting and analyzing the nonlinearity-caused oscillations are provided in Figure 1.

IV. SIMULATIONS
In this section, a numerical example and a four-input-fouroutput control system are used to validate the effectiveness and advantages of the proposed method.

A. NUMERICAL EXAMPLE
Firstly, a numerical example is tested to show the effectiveness and advantages of the proposed IVNCMD. Signal (13) consists of three sub-signals, including two sinusoidal components and a time-varying mode. g (t) = cos (2πt)+2 cos 2t 2 + 10t +1.5 cos (2π35t)+e (13) where e ∼ N (0, 0.2) is noise. The trend of ApEn value versus mode number K is plotted in Figure 2. It is observed that when mode number is larger than 3, the ApEn value will be higher than the threshold. According to Algorithm 1, the mode number should be set as 3. The instantaneous frequencies and decomposition modes are provided in Figure 3 and 4, respectively. It can be seen from Figure 3 that the estimated instantaneous frequencies are completely consistent with the true frequency curves, which indicates the proposed IVNCMD correctly captures the time-frequency information. The decomposition results shown in Figure 4 also demonstrates the IVNMCD is able to extract the sub-signals contained in the original signal (13). Therefore, it is claimed that the proposed IVNCMD can process the complex nonstationary signal and providing accurate time-frequency information. Besides, we apply the original NCMD to signal (13) with inappropriate mode number. Taking k = 2 and 4 as examples, the decomposition performances of NCMD is shown in Figure 5 and 6, respectively. It can be seen that NCMD suffers from severe consequences, such as mode-mixing and endeffect, due to the inappropriate mode number setting. With  respect to the traditional methods, such as EMD [40] and VMD [28], the corresponding decomposition outputs are displayed in Figure 7 and 8, respectively. It is observed that both EMD and VMD meet fateful difficulties. Specifically, EMD is subjected to mode-splitting problem [41] and VMD is influenced by the time-varying factor. Therefore, it can be concluded that the proposed method outperforms NCMD, EMD and VMD in adaptively decomposing nonstationary signals, which would guarantee the detection performance of process oscillations in the following.

B. CONTROL SYSTEM
A four-input-four-output control system is taken from [42], [43] to validate the effectiveness of the IVNCMD-based detector. The structure diagram is shown in Figure 9. The corresponding transfer functions are listed in Table 1. In this system, four single-loop PI controllers are utilized with where K c is the proportional gain, T is the controller sampling time and T r is the integral time. The controller parameters are K c = [0.816, 0.625, 0.184, 0.37] and T r = [20, 16, 2.86, 5], respectively. In order to generate the nonlinearity-induced oscillations, a valve stiction model [42], [44] is embedded into the third loop to simulate the nonlinearity problem. Because the loops are connected with each other, the oscillations generated in Loop 3 will propagate to other loops. Herein, we apply the proposed IVNCMD to these four process variables in succession.
For the first loop, its process variable is plotted in the first row of Figure 11. According to the ApEn value in Figure 10, IVNCMD-Based Detector for Nonlinearity-Caused Oscillations Input: Process variable g (t); Output: Whether this oscillation is caused by nonlinearity problems; 1: Decompose the process variable g (t) into K modes using Algorithm 2; 2: Detect the oscillating modes according to α k > α µ and β k > β µ , for k = 1, 2, . . . , K ; Assume there are J significant oscillating modes retained. 3: Calculate the average oscillation period P j using (11) and the corresponding upper and lower limits are obtained as P U j = 1 P j − σ P j and P L j = 1 P j + σ P j , respectively; 4: Regard the most correlated mode (according to the normalized cross-correlation coefficient (10)) as the fundamental oscillation component and the corresponding average oscillation period is the base frequency; 5: Let n represent an integer. i corresponds to the index of the most correlated mode and j = i 6: if nP i ∈ P L j , P U j ,j = 1, 2, . . . , J then 7: The n-th order harmonic is detected, which demonstrates the presence of nonlinearity-caused oscillations 8: else 9: The oscillation is not caused by nonlinearity problems. 10: end if the mode number should be set as 1. Figure 11 provides the decomposition result of IVNCMD. The corresponding detection results are listed in the first row of Table 2. It is claimed that only one oscillating mode with frequency about 0.02 Hz is detected, i.e., there is no higher order harmonic. Therefore, Algorithm 3 reports this oscillation is not caused by nonlinearity problems. In light of the presetting fault, the oscillation in this loop is transmitted from other loops rather than caused by nonlinearity. Thus, the analysis result is consistent with the facts.
Similar to Loop 1, Figure 12 and 13 display the ApEn curve (according to Algorithm 2) and decomposition outputs of the process variable of Loop 2, respectively. The detailed information detected by Algorithm 3 is listed in the second case of Table 2. In terms of these information, the detected oscillation is identified as linear type by Algorithm 3 because no higher order harmonics are captured. This result is in accordance with the actual situation.
With respect to Loop 3, the mode number is determined as 2 because the ApEn value (shown in Figure 14) exceeds the threshold when the mode number is larger than 2. According VOLUME 9, 2021      to Figure 15 and Table 2, the process variable is distinguished as nonlinearity-caused oscillation, because a third order harmonic is detected. This analysis results conform to the fact that a valve nonlinearity is embedded into this loop, which leads to the loop oscillating.
The detection and analysis results for the last Loop are given in Figure 17 and Table 2. It is observed that the detected oscillation is not classified as nonlinearity-caused oscillation because it is an external disturbance from Loop 3. This case is similar to Loop 1 and 2.  In addition, because different parts of a plant tends to act as low pass mechanical filters, which can filter out the higher harmonics as we move away from the source of nonlinearity. Therefore, it is reasonable to regard Loop 3 as the oscillation source since it shows the highest order harmonic according to Table 2.
To conclude, the proposed method is able to accurately detect the process oscillation and correctly analyze whether the oscillation is caused by nonlinearity problems. The obtained results also can contribute to root cause analysis for unit-wide oscillations.

V. INDUSTRIAL CASE
This section considers a more complex case study, where time trends from an Australian refinery separation unit are tested [45]. The simplified structure diagram is shown in Figure 18. This data set includes the steam flow, analyser and temperature measurements (PV). Measurements from the upstream and downstream pressure controllers PC1 and PC2 are also contained. The sampling time is 20s. It is well-established that this plant displays a unit-wide oscillation in a distillation column. The time trend of the analyzer indicates the composition of the product leaving the top of column was varying in an undesirable behavior. It is known that there was a nonlinearity problem in the steam flow loop FC1 [46]. It was an orifice plate flow meter but there was no weep-hole in the plate, which had the effect that condensate collected on the upstream side until it reached a critical level, and the accumulated liquid would then periodically clear itself by siphoning through the orifice. Another disturbance is presented in PC1 and PC2, whose root cause is likely to be due to a controller interaction between the two pressure loops. The oscillations contained in PC1 and PC2 are linear [37]. The process variables and decomposition outputs of IVNCMD for these loops are provided in Figure 19, 20, 21, 22, and 23, respectively. It is observed that except for the loop FC1, which has two oscillating modes, the rest loops only have one mode. Based on the information provided in Table 3, it is reported that the second harmonic is detected in Loop FC1 and there is no higher order harmonic in other loops. According to Algorithm 3, the oscillations in loop FC1 should be caused by nonlinearity issues, which is consistent with the actual situation [46]. With respect to the oscillations in loop AC1 and TC1, their oscillation frequency is about 0.0024 Hz, which is the same as that of loop FC1. Due to the low-pass filter effect, the higher order harmonics contained in loop FC1 are filtered out during the transmission process. Therefore, it is referred that both oscillations contained in loop AC1 and TC1 is transmitted from loop FC1. In addition, measurements from the upstream and downstream pressure controllers PC1 and PC2 show different oscillation frequency, thus these two oscillations are likely to be caused by other factors rather than nonlinearity problem. The analysis results tally with the previous studies [46].
For comparisons, we apply the latest VMD-based method [2] to this dataset. The corresponding results are provided in Table 4. It is observed that the VMD-based method does not correctly decompose the signals collected from loop FC1 and TC1. Specifically, the second order harmonic contained in loop FC1 is falsely extracted and the oscillation frequency of loop TC1 obtained from VMD is not consistent with the prior. Figure 24 and 25 displays the incorrect decomposition results of VMD for loop FC1 and TC1, respectively. Therefore, it is claimed that the proposed method outperforms the VMD-based approach [2] in detecting and analyzing the nonlinearity-caused oscillations for process control systems.
In addition, the EMD-based decomposition results of this industrial case are also provided in figure 26, 27, 28, 29, and 30. It is evident that the EMD decomposition results suffer from mode-mixing and end-effect issues. These issues destroy the effective information extracted by EMD. Therefore, it can be concluded that the proposed method is effective and shows better performance than the EMD-based method in practice.

VI. CONCLUSION
In this paper, we propose an improved variational nonlinear chirp mode decomposition (IVNCMD) algorithm, which utilizes the approximate entropy of instantaneous frequency to adaptively determine the mode number of NCMD. Compared with EMD and VMD, the proposed IVNCMD shows better performance in decomposing nonstationary signals and providing accurate time-frequency information. Then, a novel IVNCMD-based detector is developed to extract the oscillations from process variables. Through investigating the instantaneous frequency relationship between modes, a novel nonlinearity-caused oscillation analysis algorithm is presented. In the end, the effectiveness and advantages of the presented method are demonstrated by simulations as well as industrial cases. More specifically, the industrial case study reports that the proposed method is able to correctly identify the nonlinearity-induced oscillations, while the VMD-based method fails. In addition, the results obtained from EMD are subject to severe mode-mixing and end-effect issues, which make the decomposition results meaningless. Therefore, it can be concluded that, in contrast to the existing methods, such as VMD and EMD, the results obtained from the proposed method is more reliable and accurate in practice. In the future, we will further study the plant-wide oscillation detection.