Detecting and Characterizing Multiple Oscillations in Process Control Loop Based on PSO-NCMD

Oscillations are a common abnormal phenomenon in the process control system. They can degrade control performance or even cause plant shutdown. It is crucial to accurately detect and characterize the process oscillations. In this article, a novel oscillation detector is proposed by combining the particle swarm optimization (PSO) and nonlinear chirp mode decomposition (NCMD). Because the performance of NCMD relies on the selection of mode number <inline-formula> <tex-math notation="LaTeX">$Q$ </tex-math></inline-formula> and bandwidth parameter <inline-formula> <tex-math notation="LaTeX">$\alpha $ </tex-math></inline-formula>, PSO is utilized to search the optimal parameter pairs. Then, the multiple oscillations contained in the process variables (PV) can be extracted by NCMD with the optimal parameters. The normalized correlation index and sparseness index are used to discarding the spurious modes and quantifying the degree of oscillations, respectively. After detecting, by utilizing the time-frequency information provided by the oscillating modes, multiple oscillation types can be accurately characterized. Comparisons are provided to show the advantages of the proposed method. The effectiveness and utility are validated by simulations as well as various industrial cases.


I. INTRODUCTION
Oscillations are a common adverse factor in process control loops, which often result in the waste of raw materials and energy, product quality fluctuation, or even threatening the stability and safety of the plant. External disturbances, poor controller tuning, and valve stiction are the main causes of oscillations. The oscillations caused by the first two factors are generally regarded as linear oscillations, while the stiction-induced oscillations are classified as nonlinear oscillations [1], [2]. In order to accurately eliminate the fault, it is crucial to detect the oscillations and characterize their types as early as possible.
In the past two decades, many oscillation detection techniques have been developed [3]. In 1995, Hägglund [4] proposed the first workflow of detecting and diagnosing oscillations. It utilized the integral absolute errors (IAEs) to measure the regularity of control errors. If IAE values exceed the predefined limit, it is conclude the oscillations have generated. This method is simple and intuitive; while it only could detect single oscillations. Its corresponding diagnosis The associate editor coordinating the review of this manuscript and approving it for publication was Salman Ahmed . framework is not data-driven, which needs to manually stop the control loop. Then Forsman and Stattin [5] modified the IAE-based methods. Miao and Seborg [6] used the decay ratio property of auto-covariance function (ACF) to detect oscillation. Regularity of zero crossings [7] is also the common index to indicate the presence of oscillations. Nevertheless, the above methods are limited to detecting single oscillations [8].
Due to the presence of nonlinearity, nonstationarity and noise, etc., the detection and characterization of multiple oscillations is still a challenging problem at present. Naghoosi and Huang [9] combined clustering algorithm and ACF to detect multiple oscillations and estimated the corresponding oscillating frequencies. Because the multiple oscillations can be regarded as the frequency components contained in the data, researchers and engineers tended to use signal processing tools to solve this problem. Fourier power spectrum was adopted in [10] to detect plant-wide oscillations. But Fourier transform is only suitable for processing linear and stationary signals, which is difficult to meet in practice. Matsuo and Sasaoka [11] introduced the wavelet transform into detecting multiple oscillations. However, this method required to define the decomposition level and mother wavelet in advance. With the development of signal decomposition techniques, they have been a powerful tool for oscillation detection. The basic assumption of signal decomposition is that the complex signals are composed of a series of modes, which are relatively simple signals. Through extracting these modes, it would be easier to mine the information in the data. Huang et al. [12] first proposed the empirical mode decomposition (EMD), which extracted the intrinsic mode functions (IMFs) through a recursive sifting process. Srinivasan et al. [13] applied EMD to characterizing oscillations. Nevertheless, EMD often suffers from mode-mixing and end-effect, which damage the decomposition performance and reduce the accuracy and reliability of the detection results. Although some modified version of EMD [14], such as EEMD [15] and CEEMD [16], are developed successively to remedy these issues, the selection of parameters of these method is subjective and lacks criterion. Guo et al. [17] proposed an improved intrinsic time-scale decomposition (ITD) to detect the oscillations online. ITD has the advantage of low computational complexity, but its time resolution is poor. Then, Xie et al. [18] developed a detector based on local mean decomposition (LMD) and Lempel-Ziv complexity to monitor time-varying oscillations. However, LMD is subjected to mode-mixing, which degrades the detection performance [19].
The above signal decomposition methods adopted in the above detection schemes are empirical and lack theoretical basis. Starting from variational mode decomposition (VMD) [20], those methods based on optimization theory have become a hot topic. More and more related practical applications have emerged [21]. Although VMD shows good performance in anti-mode-mixing and robustness, it is limited to decomposing narrow-band signals. It only could provide the center frequencies (scalar), which indicates a lack of time-frequency representation [22]. Recently, Chen et al. [23] proposed a nonlinear chirp mode decomposition (NCMD) algorithm, which is the latest development in the field of signal decomposition. NCMD is established on the fact that a wide-band signal can be transformed into a set of narrow-band signals [24]. It formulates the signal decomposition problem as a demodulation issue and extracts multiple modes concurrently. Because NCMD shows excellent decomposition performance even when the modes are close or crossed, it has been widely used in various applications, such as bearing fault diagnosis [25] and wind speed forecasting [26], etc. However, the performance of NCMD relies on the selection of mode number Q and bandwidth parameter α. How to select the parameters is an important issue, but it has not been studied yet.
In this article, in order to achieve the best performance of NCMD, particle swarm optimization (PSO) [27] is used to search these two parameters, thus a PSO-NCMD algorithm is successfully developed. Then, process variables are decomposed by PSO-NCMD with the optimal parameters, and the corresponding modes are obtained. Following, a novel oscillation detector is proposed based on correlation index and sparseness index. Specifically, the spurious modes can be discarded by the normalized correlation index. And oscillation properties of the retained modes are quantified by sparseness index. After detecting, a power-weighted frequency is developed to estimate the oscillation periods, which also can be used to characterizing the type of multiple oscillations.
The contributions of this work are threefold, (i) A PSO-NCMD algorithm is proposed to improve the performance of the original NCMD; (ii) A novel multiple oscillation detection and characterization framework is developed based on PSO-NCMD and power-weighted frequency; (iii) Compared with the the existing detection methods based on EMD or VMD, the proposed method shows better performance in anti-mode-mixing, detection accuracy and reliability.
This article hereafter is organized as follows. The NCMD and PSO are briefly introduced in Section II. The construction of PSO-NCMD and oscillation detection and characterization framework are described in Section III. Simulations and industrial cases are provided in Section IV and V, respectively, followed by conclusions.

II. THEORETICAL FOUNDATION
A brief introduction of NCMD and PSO algorithm is presented in this section.

A. NONLINEAR CHIRP MODE DECOMPOSITION
In NCMD [23], the complex signals are assumed to be composed of several sub-signals, called as nonlinear chirp modes. More specifically, this relationship is shown as where Q is the mode number; a i , f i , and φ i are the instantaneous amplitudes (IAs), instantaneous frequencies (IFs), and initial phases, respectively; η (t) stands for the noise. Generally, the fluctuation of IAs and IFs is much smoother than that of phases.
The goal of NCMD is to extract these modes m i from the given signal s. Firstly, NCMD transforms the signal s into its de-chirped form through demodulation techniques, where u i (t) and v i (t) are two demodulated signals, wheref i denotes the frequency functions of cos(2π t 0f i (z)dz) and sin(2π t 0f i (z)dz), which are the demodulation operators. The corresponding IA can be recovered as a i = (u i (t)) 2 + (v i (t)) 2 . Iff i is matched to that of the true IF f i , the demodulated signals will have narrowest bandwidth. Therefore, an optimal demodulation problem is established on the basis of minimizing the bandwidth of u i (t) and v i (t), shown as where ε is the reconstruction error tolerance, which can be estimated by noise level; The square of the l 2 -norm of the second-order derivative is adopted to estimate the signal bandwidth. (4) is expressed in a continuous form; while the signals sampled in practice are discrete. Thus, it is necessary to discretize (4). Assuming time vector is t = t 0 , . . . , t N −1 where N is the number of samples. The corresponding discrete version of (4) can be written as where and is the second-order difference operator By introducing an auxiliary variable ω, the inequality constraints in (5) can be transformed into equality constraints, where I C ε (ω) is the indicator function of an Euclidean ball [23]. The corresponding augmented Lagrangian of (14) is where α and λ are the quadratic penalty coefficient and Lagrangian multiplier, respectively. (15) can be effectively solved by the alternating direction method of multipliers (ADMM). Refer to [23] for more details.
As can be seen from (15), there are two parameters (mode number Q and quadratic penalty coefficient α) that need to be specified in advance. Q determines the number of modes obtained from decomposition. A large or small Q value will lead to over-decomposition or under-decomposition issues, which degrade the performance of NCMD a lot. NCMD can be viewed as a time-frequency filter, which has similar filter bank property with Vold-Kalman filter [28]. α is used to regulate the filter bandwidth. However, There is no criterion for the selection of α value. Therefore, it is crucial to develop a method to search both parameters.

B. PARTICLE SWARM OPTIMIZATION
PSO [29] is a kind of swarm intelligence optimization algorithm, which imitates the birds' foraging. It has been widely used in various applications, such as controller design [30], [31], optimal scheduling [32], [33], to name a few. It first initializes a group of particles in the solution space, each of which represents a potential optimal solution. The particle characteristics are described by three indexes: position, velocity and fitness value. The fitness value is calculated by fitness function, and its value indicates the quality of particles. These particles move in solution space and update their position by tracking individual and group extremum. More specifically, PSO supposes that there are J particles in D-dimension solution space. In each process of iteration, the particles update their own velocity and position by the information of the individual optimal value and the global optimal value in accordance with (16) and (17) and where V d j stands for the amount in which will change in magnitude and direction in the next iteration; P d j represents the current position for each particle. The subscripts · b and · g denotes the best position that has been discovered by a particular particle and the best global solution of the population in the current iteration, respectively. The lower subscript · j and the upper subscript · d are the index of particle and dimension, respectively. r 1 and r 2 are random numbers distributed in [0, 1]. The weighting factor , c 1 and c 2 are empirically set to be 0.7, 1.4, and 1.4, respectively.

III. PROPOSED METHOD A. PSO-NCMD
Defining an appropriate fitness function is the key to find the parameters of NCMD with PSO. There are many various common measurement indexes, such as sparsity measurement [34], smoothness index [35], correlation coefficient [23], and permutation entropy [36], etc. In this work, NCMD is expected to be applied to decompose the oscillations in process control loop, which often show almost regular periodic variation property. Both correlation coefficient (CC) and permutation entropy (PE) are widely adopted in detecting process oscillations [37]. Following, these two indexes are briefly described.
Correlation coefficient can be used to quantify the correlation between the original signal and the decomposed modes.
where cov (·) and σ represent the covariance and standard deviation, respectively. If the correlation coefficient is too small, it means that the corresponding mode may be spurious. Permutation entropy is developed from information theory [38] and can be used to measure the complexity for time series. Compared with other nonlinear monotonous transformations, the computation of permutation entropy is efficient and robust. Given a one-dimensional time series {x t , t = 1, 2, . . . , N }, it can be reconstructed in phase space with embedding dimension d and time delay τ , Each row of matrix Y i = x i x i+τ · · · x i+(d−1)τ is a reconstruction component. The reconstruction components are arranged in ascending order according to the index i, where i 1 , i 2 , . . . , i d represent the sequence number of the column in which each element in the reconstruction component is located. If two elements in the reconstruction component have the same value, namely, x i+(i 1 −1)τ = x i+(i 2 −1)τ , they will be sorted by the size of i 1 and i 2 , i.e., if Following, the reconstruction components in Y can be mapped into a group of Symbol sequence, where l = 1, 2, . . . L and L ≤ d!. Sym (l) is one of the d! permutations of d distinct symbols. Denoting the probability of occurrence of each permutation as P 1 , P 2 , . . . P L , the permutation entropy H (P) can be defined based on Shannon entropy, For the convenience of application, H (P) is usually normalized into the range from 0 to 1 .
The more regular the time series, the smaller the corresponding PE. If the parameters of NCMD are appropriate, the obtained modes will be regular time series and the PE value will be small. Therefore, PE can be used to evaluate the decomposition performance of NCMD. Permutation entropy is only dependent on the regularity of the decomposed mode itself. If a mode is regular but does not contain valuable information, for example, this mode is extracted from noise, the corresponding PE will be small, too. Thus, it is not proper to only use the permutation entropy as the fitness function of NCMD. Correlation coefficient is able to measure the similarity between signals, which can be utilized to characterize the information between the original signal and modes. Therefore, a synthetic index, termed PECC, is developed by combining PE and CC to serve as the fitness function for the NCMD, where NH and ρ are the normalized permutation entropy and the correlation coefficient, respectively. In PECC index, the modes' PE values are weighted by their corresponding CC.
To sum up, the proposed method can be described as where i is the mode index; Q and α are the mode number and bandwidth parameter of NCMD. PECC (24) is mode's VOLUME 8, 2020 fitness function, which is used to quantify the decomposition performance of NCMD. The optimal parameter pairs of Q and α are obtained by PSO. The procedures of the proposed method is summarized in Algorithm 1.

B. DETECTING AND CHARACTERIZING MULTIPLE OSCILLATIONS
After obtaining the optimal parameters of NCMD, the process variables can be well decomposed into a series of modes, which contain valuable information about multiple oscillations. Herein, the normalized correlation coefficient ξ [23] and sparseness index [39] are combined to detect the presence of oscillating modes. More specifically, the normalized correlation coefficient ξ i can be calculated as where ρ i is the correlation coefficient obtained from (18). This indicator is used to identify and discard the spurious modes, which are a common issue in signal decomposition techniques. Only modes with ξ i > T ξ is regarded as significant modes and retained for the following analysis. The sparseness index ς i is where F stands for the Fourier transform; m i is decomposition modes and N is the data length. If ς i value approaches 1, it means the frequency response of m i shows evident peaks, which correspond to the oscillation components. On the contrary, ς i value close to 0 indicates the corresponding modes are DC components or noise. Only modes with ς i > T ς are considered as an oscillating component. The thresholds of ξ i and ς i are empirically set to be 0.15 and 0.58, respectively [40].
Having detected the oscillations, the next task is to characterize the oscillations' properties, such as the frequency information and oscillation type. In regard of oscillation frequency, previous studies [41], [42] only consider a few zero-crossing internals, which do not fully mining the information. Although [40] uses the full frequency information to reveal the oscillation properties, its oscillation frequency is calculated as the mean value of mode's instantaneous frequency. It is not proper to make a simple arithmetic mean of instantaneous frequency, because this way loses the amplitude information. When the oscillations are timevarying, the arithmetic mean will lead to big deviation. Herein, inspired by [43], we propose a power-weighted average to estimate the oscillating frequency f osci , as shown in where f i (t) and a i (t) are the instantaneous frequency and instantaneous amplitude obtained from NCMD, respectively. This power-weighted framework is able to correctly calculate the frequency information of both time-invariant oscillation and time-variant oscillation. Based on the oscillation frequency information, it is feasible to further identify the type of oscillations. For example, the nonlinear oscillations can be characterized if there are higher order harmonics among various modes [42]. The detection and characterization procedures are provided in Algorithm 2. The complete flowchart of the proposed framework is shown in Figure 1, which integrates parameter optimization, oscillation detection and characterization. Apart from processing single oscillations, it can effectively detect and characterize the multiple oscillations.

Algorithm 2 Detecting and Characterizing Multiple Oscillations
Input: Collect process variable as input signal Output: The oscillation frequency (f osci ) and type (nonlinear or linear) 1: Decompose the input signal into a set of modes using NCMD with the optimal parameters 2: Calculate the oscillation indices for each mode 3: Regard modes whose normalized correlation coefficient and sparseness index satisfy ξ > T ξ and ς > T ς as oscillating modes 4: Calculate the power-weighted frequency of oscillating modes according to (28)  The higher order harmonics are detected, which indicates the presence of nonlinear oscillations 8: else 9: The oscillation type is linear 10: end if

IV. SIMULATIONS
In order to verify the improvement of PSO on the optimization of NCMD parameters, a synthesized signal (29) is tested  as the subject. s = cos (2π t)+cos (10πt)+2 cos(16πt +40πt 2 )+η (29) where η ∈ N (0, 0.1) is noise. The PSO settings are as follows: the maximum iteration number and particle number are 20 and 30, respectively. The fitness value of PSO is plotted in Figure 2. It is observed that the fitness value converges after about 12 iterations. Figure 3 and 5 display the decomposition performance of NCMD in the initial and final stages of optimization, respectively. It can be seen that the decomposition performance in the case of optimal parameters is much better than that in the initial iteration. Their corresponding instantaneous frequencies are plotted in Figure 4 and 6, respectively. It can be seen that the true time-frequency curve  is in good agreement with the extracted time-frequency curve under the appropriate parameters; while if the parameters are not proper, the extracted IF curve would deviate greatly. Therefore, the proposed PSO-NCMD can effectively improve the performance of the original NCMD.
For the sake of contrast, the decomposition performance of VMD and EMD are provided in Figure 7 and 8, respectively. It is observed that both VMD and EMD cannot correctly decompose signal (29), because they are not good at dealing with time-varying signals. Therefore, compared with VMD and EMD, the proposed PSO-NCMD has better performance in processing complex signals.
Following, a simple single-input and single-output closed loop is tested to show the effectiveness of the proposed detection and characterization framework. The feedback system is shown in Figure 9. Its transfer function is  This loop is regulated by a PI controller with P = 0.1 and I = 0.5, which can enable the system run stably. External disturbances, poor controller tuning, and valve stiction are the main causes of oscillations in process control system. Totally four cases, namely, external disturbances (linear), poorly tuned controller-caused oscillations (linear), stiction-induced oscillations (nonlinear), and both external disturbances and stiction-induced oscillations (nonlinear and linear), are investigated. Their data and the corresponding decomposition performance are displayed in Figure 10, 11, 12, and 13, respectively. It is observed that all oscillating signals are well decomposed. The detailed descriptions on this experiment are provided as follows.
Firstly, this single-input and single-output loop is contaminated by an external oscillatory disturbance with amplitude 1 and frequency 1.15 rad/s, which is approximately equal to 0.1830 Hz. The parameters of PI controller are P = 0.1 and   I = 0.5. According to Algorithm 1, the optimal parameters of NCMD are Q = 1 and α = 4e − 8. It is observed that the predefined oscillation can be accurately extracted by NCMD with optimal parameters. The obtained power-weighted frequency is 0.1830 Hz (in Table 1), which is consistent with the preset disturbance. Because only one mode is extracted, no higher order harmonics are characterized. It is inferred that this is a linear oscillation.
Secondly, P remains the same while the integral parameter of controller is tuned to I = 1.5, which leads to an excessive integral action. The loop outputs apparent oscillations, as shown in the first row of Figure 11. Note that the oscillations induced by external disturbances or poorly tuned controller are significant. The optimal mode number Q and bandwidth parameter α obtained from PSO-NCMD are 1 and 3.5e − 6, respectively. The decomposition performance is  displayed in the second row of Figure 11. The second case in Table 1 clearly indicates that the poorly tuned controller leads to a linear oscillation with 0.2389 Hz.
Thirdly, a simulated valve stiction model [44] is embedded into the control loop, which causes nonlinearity-induced oscillations. Using the PSO-NCMD, the process variable is decomposed into two modes, which are shown in Figure 12. The third case of Table 1 provided the oscillating frequencies VOLUME 8, 2020  of these two modes, which indicates the presence of higher order harmonics in this loop. Thus the nonlinear oscillations are confirmed.
The last simulation is a mixture of linear and nonlinear oscillations, which is more complex than single oscillation type. Although processing such multiple oscillations is challenging, the proposed method still correctly extracted three modes by PSO-NCMD. According to the oscillation detection and characterization algorithm, it can be judged that the first mode and the second mode are induced by nonlinearity, while the last one is caused by linear oscillation source. Therefore, it is concluded that besides single oscillations, the proposed method is able to hand with complex multiple oscillations.
The corresponding detection and characterization results are listed in Table 1. The oscillation frequency and characterized oscillation type of various faults are consistent with the predefined settings, which demonstrates the proposed framework can be used to detect and characterize various complex oscillations in the process control system.

V. INDUSTRIAL CASES
Since the above simulations on the proposed methodology acts well, more widespread industrial applications can be expected. In this section, the proposed method is applied to data sets sampled from various industries [1], such as commercial buildings, chemicals, pulp and paper mills, power plants, mining, and metal processing.
More specifically, the first industrial data are sampled from a temperature control loop in the commercial building. Its sampling time is 1s. Figure 14 displays the oscillation signal and the corresponding decomposition modes. The characterization results are listed in the first case of Table 2. The higher order harmonics are detected in this loop, thus it can be concluded that this oscillation is introduced by nonlinearity problems. The second case is collected from a flow control loop in chemicals and the sampling time is 12s. The decomposition performance of the proposed method is provided in Figure 15. The second case of Table 2 reports the presence of higher order harmonics, thus the nonlinear oscillations are confirmed. With respect to the signals obtained from a liquid  flow control loop in the pulp and paper mills, the corresponding decomposition outputs and characterization results are shown in Figure 16 and Table 2, respectively. The sampling time is 1s. It can be seen that the third order harmonic is captured by the proposed method, thus the oscillation is regarded as a nonlinear type. The fourth control loop is used to regulate the level in the power plant. Its sampling time is 5s. It is observed from Figure 17 and Table 2 that the proposed method captures the second order harmonic, thus identifying the nonlinearity-induced oscillations. The fifth industrial data set is collected from a temperature control loop in the mining process, whose sampling time is 60s. Its corresponding decomposition and detection results are depicted in Figure 18 and Table 2, respectively. It is apparent that the presence of third order harmonic indicates the oscillation is caused by nonlinearity. The last example is utilized to control thickness for metals and its sampling time is 0.05s. Figure 19 and Table 2 display the decomposition modes and characterization results, which indicate this signal is a linear oscillation,  because no higher order harmonic is detected. According to the above experiments, it is concluded that the proposed   method could elegantly extract the oscillation modes from complex signals in various practice environments. The detection and characterization results listed in Table 2 are well consistent with the prior knowledge, which demonstrates the effectiveness of the proposed method.
In order to highlight the advantages of the proposed method, the EMD-based [13] and VMD-based [22] methods  are used to detect the oscillations contained in these data. The decomposition plots for the three misclassified cases using EMD (for pulp and papers, metals), or VMD (for mining) are provided in Figure 20, 21, and 22, respectively. The corresponding detection and characterization results are listed in Table 3. It is observed that both EMD-based and VMD-based methods make some mistakes. For example, EMD-based method regards the oscillations contained in the pulp and paper mills are induced by a linear source, because there are no higher order harmonics are detected. This result is inconsistent with the facts. EMD-based method also meets difficulties in the case of metal process, because the decomposition performance of EMD is unsatisfied. In addition, VMD-based method detects two oscillation modes for the temperature control loop in the mining plant, but the corresponding oscillation frequencies obtained by VMD are not accurate, thus resulting in a false report. Therefore, it is concluded that the proposed framework can be applied to detect and characterize the oscillations in various actual process control systems, and it shows better performance than other methods.

VI. CONCLUSION
In this article, the PSO algorithm is utilized to search the optimal parameters of NCMD. The correlation coefficient and permutation entropy are combined to establish a reasonable fitness function, which can quantify the decomposition performance of NCMD. Following, based on the decomposition results of PSO-NCMD, a novel oscillation detector is proposed by using the normalized correlation coefficient and the sparseness index. Then a power-weighted frequency is developed to indicate the multiple oscillating frequencies.
The oscillation type can be characterized through investigating the frequency relationship among various modes. Simulations shows that the proposed PSO-NCMD is able to effectively improve the decomposition performance of the original NCMD. In addition, PSO-NCMD outperforms EMD and VMD in decomposing complex signals. In the end, the SISO experiments and various industrial cases demonstrate the utility and advantages of the proposed oscillation detection and characterization framework. Compared with the the existing detection techniques based on EMD [13] and VMD [22], the proposed method is more accurate and reliable. Therefore, it can be concluded that this study has great potential values for control performance monitoring.
In the future, we will focus on detecting plant-wide oscillations and locating the corresponding root causes. For example, some causality analysis methods, such as Granger causality, transfer entropy, and Bayesian network, can be used to find the root cause of oscillations.