Adaptive detection method for industrial loop oscillation under variable working conditions

Abstract

The invention discloses a self-adaptive detection method for industrial loop oscillation under variable working conditions, and belongs to the field of industrial loop performance evaluation. The method comprises the following steps: carrying out adaptive variational mode decomposition on PV of an industrial loop; calculating the relative growth rate of the correlation coefficient, and selecting an important mode; calculating a sparse index of an important mode, and removing a noise mode; estimating the confidence upper and lower limits of the regularity index of the residual mode, and judging the possibility of oscillation; calculating the relative error between the oscillation period of the residual mode and the oscillation period of the SP, and eliminating the interference of the periodic demand; and judging by using the lower confidence limit of the regularity index of the residual mode, and reporting the specific oscillation condition. According to the method provided by the invention, from the perspective of loop oscillation under variable working conditions, interference of non-stationary factors and periodic requirements can be effectively eliminated, the automation degree and accuracy of industrial loop oscillation detection are improved, and a method basis is provided for online monitoring and performance evaluation of an industrial loop.

Description

Technical Field

[0001] The invention belongs to the field of industrial circuit performance evaluation, and in particular relates to an adaptive detection method for industrial circuit oscillation under variable working conditions. Background Art

[0002] Industrial loop oscillation is a common fault phenomenon that affects the performance of industrial control systems. The performance degradation of the control system will lead to fluctuations in product quality, increased production costs and reduced economic efficiency. In severe cases, it may even threaten the stability and safety of the entire system. During operation, the system often has the need for variable operating conditions, such as system setting operating condition switching or load changes. At this time, non-stationary factors will be introduced. At the same time, the system itself will also be affected by random noise and process disturbances, resulting in extremely strong non-stationarity of the system signal under variable operating conditions. There are often trend terms, noise and multiple frequency components in the signal. In this case, the results of oscillation detection are seriously affected, and even misdiagnosis may occur. However, most oscillation detection methods do not consider the variable operating conditions in the industrial process when facing non-stationary oscillations, and most of the existing methods cannot adaptively detect oscillations for input signals. For example, the number of modes and the number of zero crossing points of variational mode decomposition need to be set in advance, which hinders the actual deployment of detection and cannot meet the requirements of high automation. Summary of the invention

[0003] The purpose of the present invention is to overcome the defects in the prior art and provide an adaptive detection method for industrial circuit oscillation under variable working conditions.

[0004] The specific technical solutions adopted by the present invention are as follows:

[0005] The present invention provides an adaptive detection method for industrial circuit oscillation under variable working conditions, which is as follows:

[0006] S1: Obtain the setting signal SP, process variable signal PV and sampling frequency f in the industrial circuit to be detected s ; After adaptive variational mode decomposition of PV by comprehensive decomposition performance index, K modal IMFs are obtained k , k=1,2,…,K;

[0007] S2: Calculate the relative growth rate δ of the correlation coefficient k ; If δ k If it is greater than the preset value a1, the mode is judged to be important and S3 is performed; otherwise, the corresponding mode is eliminated;

[0008] S3: Calculate the sparse index SI of important modes k If SI k If it is greater than the preset value a2, the mode will proceed to S4, otherwise the corresponding mode will be eliminated;

[0009] S4: Estimate the confidence limits of the regularity index of the remaining modes after S3 screening [R 1_IMFk ,R 2_IMFk ]; if the upper confidence limit of the regularity index R 2_IMFk If it is greater than the preset value a3, the mode is judged to be oscillating and S5 is performed; otherwise, the corresponding mode is eliminated;

[0010] S5: Calculate the relative error e between the oscillation period of the remaining modes after S4 screening and the oscillation period of SP k If e k If it is greater than the preset value a4, the mode proceeds to S6, otherwise the corresponding mode is eliminated;

[0011] S6: Determine the lower confidence limit R of the regularity index of the remaining modes after S5 screening 1_IMFk Is it greater than the preset value a3? If so, report that the mode is oscillating and output its oscillation period and oscillation amplitude; otherwise, report that the corresponding mode may oscillate and output the oscillation possibility, oscillation period, and oscillation amplitude.

[0012] Preferably, in S1, the objective function of the adaptive variational mode decomposition is:

[0013]

[0014] Among them, {u k} is the set of modal IMFs, {ω k} is the corresponding set of center frequencies, K is the number of decomposed modes, x(t) is the process variable signal PV, u k (t) is the modal IMF k , δ(t) is the unit impulse function, i is the imaginary unit, * represents convolution, It means to find partial derivative, t is time, ω k is the modal IMF k The center frequency of

[0015] Furthermore, the steps of the adaptive variational mode decomposition are as follows:

[0016] S11: Initialize the number of decomposed modes K=2;

[0017] S12: Decompose PV using adaptive variational mode decomposition to obtain K modal IMFs k and the corresponding center frequency ω k , where k = 1, 2, ..., K;

[0018] S13: Judgment k The maximum value max(ω k ) is greater than f s / 2.56; if so, the maximum number of decomposable modes is K max =K-1, otherwise K=K+1 and repeat S12;

[0019] S14: Determine the maximum number of decomposable modes K max Is it greater than 1? If so, S15 is performed, otherwise adaptive variational mode decomposition is not performed;

[0020] S15: Set the value range of the decomposed mode number K to [2, K max ];

[0021] S16: Calculate the comprehensive decomposition performance index CPI under different modal numbers K K , where K = 2, 3, ..., K max ; S17: Get the optimal number of decomposed modes K best , meet CPI Kbest =max(CPI K );

[0022] S18: Setting K best For the number of decomposed modes, adaptive variational mode decomposition is performed on PV to obtain K best Modal IMF k , where k = 1, 2, ..., K best .

[0023] Furthermore, the comprehensive decomposition performance index CPI K The calculation method is as follows:

[0024]

[0025] in, is the average mutual information entropy between K modes and PV, CV K is the coefficient of variation of the center frequencies of the K modes, ΔE K is the residual energy between the energy sum of K modes and the PV energy, SE K It is the smallest sample entropy among the sample entropies of K modes.

[0026] Furthermore, the The calculation formula is as follows:

[0027]

[0028] Among them, x belongs to the elements of the process variable signal x(t), and u belongs to the modal IMF k elements of , p(x,u) is the process variable signal and the modal IMF kThe joint probability density function of p(x) is the marginal probability density function of the process variable signal, and p(u) is the modal IMF. k The marginal probability density function of ;

[0029] CV K The calculation formula is as follows:

[0030]

[0031] ΔE K The calculation formula is as follows:

[0032]

[0033] Among them, E x(t) is the energy of the process variable signal x(t), E x(t) =∫x(t) 2 dt, is the process variable signal u k (t) energy,

[0034] SE K The calculation formula is as follows:

[0035] SE K =min(SampEn(u k (t))),k=1,2,…,K

[0036] Among them, SampEn() is the sample entropy function.

[0037] Preferably, in S2, the relative growth rate δ of the correlation coefficient k The specific calculation steps are as follows:

[0038] S21: Modal IMF obtained by adaptive variational mode decomposition k , according to the descending order of the correlation coefficient with the process variable signal, the IMF is obtained k ;

[0039] S22: According to the new order, the first k modes are superimposed to obtain the superimposed signal sc k ;in, k = 1, 2, 3, ... K;

[0040] S23: Calculate the superposition signal sc k Correlation coefficient with process variable signal:

[0041]

[0042] Among them, Cov is the covariance, σ x(t) is the standard deviation of the process variable signal, is the standard deviation of the superimposed signal;

[0043] S24: Calculate the relative growth rate δ of the correlation coefficient k ; If δ k >a1, the mode is judged to be important and S3 is performed, otherwise the corresponding mode is eliminated; a1 is selected as 0.05; relative growth rate δ k The expression is as follows:

[0044]

[0045] Preferably, in S3, the sparse index SI k The calculation steps are as follows:

[0046] S31: Standardization using the Z-score method for important modes;

[0047] S32: Based on the sampling frequency f s , the amplitude spectrum of the corresponding mode and the amplitude vector composed of amplitudes are obtained through fast Fourier transform;

[0048] S33: Calculate the sparse index SI of the amplitude vector k ; If there is a modal SI k If it is greater than the preset value a2, the mode is processed in S4, otherwise the corresponding mode is eliminated; a2 is selected as 0.58; the sparse index SI k The expression is as follows:

[0049]

[0050] in, The amplitude vector is obtained from the amplitude spectrum of the kth mode, N is the vector length, and f represents the fth sample point.

[0051] Preferably, in S4, the confidence limits of the regularity index are The calculation steps are as follows:

[0052] S41: Calculate the autocorrelation function ACF of the modal IMF;

[0053] S42: Obtain the position of the ACF zero crossing point;

[0054] S43: Eliminate zero crossing points within the oscillation range of ±1.96 / sqrt(Id); where I is the data length of the IMF and d is the data interval; retain the positions z(j) of the remaining zero crossing points;

[0055] S44: Obtain the periodic sequence T(j)=2(z(j+1)-z(j)) and the average period through z(j) Where j = 1, 2, ..., L; the average period Denoted as modal IMF k The oscillation period

[0056] S45: Calculate the confidence limits of regularity indicators Determine the upper confidence limit Is it greater than the preset value a3? If so, the mode may oscillate and S5 is performed, otherwise the corresponding mode is eliminated; a3 is selected as 3; the upper and lower confidence limits of the regularity index The expression is as follows:

[0057]

[0058] Where, σ is the standard deviation of the periodic sequence T(j); is the average value of the periodic sequence T(j); α is 0.05; is the 100α / 2 percentile of the chi-square distribution with L-1 degrees of freedom.

[0059] Furthermore, in S5, the relative error e k The calculation steps are as follows:

[0060] S51: Calculate the sparse index of SP and make a judgment; if the sparse index is greater than 0.58, proceed to S52; if the sparse index is less than 0.58, set its oscillation period to 0;

[0061] S52: Estimate the confidence upper limit R of the regularity index of SP 2_SP If R 2_SP >3, calculate the oscillation period T of the SP signal according to S41~S44 SP , otherwise set its oscillation period to 0;

[0062] S53: Calculate the oscillation period T of the residual mode IMFk The oscillation period T of SP SP The relative error of k If e k If it is greater than the preset value a4, S6 is performed, otherwise the corresponding mode is eliminated; a4 is selected as 0.1; the relative error e k The expression is as follows:

[0063]

[0064] Furthermore, in S6, the oscillation period Obtained by S44; the oscillation amplitude is obtained by Hilbert transform to obtain the average value b1 of the upper envelope and the average value b2 of the lower envelope, and the oscillation amplitude = (b1-b2) / 2.

[0065]

[0066] Compared with the prior art, the present invention has the following beneficial effects:

[0067] Starting from the perspective of loop oscillation under variable working conditions, the present invention can effectively eliminate the interference of non-stationary factors and periodic demands, and can adaptively select appropriate algorithm parameters for oscillation detection based on the input signal, thereby improving the automation and accuracy of industrial loop oscillation detection, and providing a methodological basis for online monitoring and performance evaluation of industrial loops. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 It is a schematic diagram of the process of the present invention;

[0069] Figure 2 A diagram of a process variable signal PV and a setting signal SP collected by an embodiment of the present invention;

[0070] Figure 3 It is a decomposition diagram of the process variable signal after variational mode decomposition in an embodiment of the present invention. DETAILED DESCRIPTION

[0071] In order to make the above-mentioned purpose, features and advantages of the present invention more obvious and easy to understand, the specific implementation mode of the present invention is described in detail below in conjunction with the accompanying drawings. In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. The technical features in each embodiment of the present invention can be combined accordingly without conflicting with each other.

[0072] like Figure 1 As shown, an adaptive detection method for industrial circuit oscillation under variable working conditions provided by the present invention comprises: (1) performing adaptive variational modal decomposition on the process variable signal PV of the industrial circuit to obtain several modes; (2) calculating the relative growth rate of the correlation coefficient and selecting the important modes; (3) calculating the sparsity index of the important modes and eliminating the noise modes; (4) estimating the upper and lower confidence limits of the regularity index of the remaining modes, and using the upper confidence limit to judge whether the mode has the possibility of oscillation; (5) calculating the relative error between the oscillation period of the remaining modes and the oscillation period of SP, and eliminating the interference of periodic demand; (6) using the lower confidence limit of the regularity index of the remaining modes to judge, and calculate the possibility of oscillation, oscillation period, oscillation amplitude and report the specific oscillation situation. The method proposed by the present invention can effectively eliminate the interference of non-stationary factors and periodic demand from the perspective of circuit oscillation under variable working conditions, improve the automation and accuracy of industrial circuit oscillation detection, and provide a methodological basis for online monitoring and performance evaluation of industrial circuits.

[0073] In a preferred implementation of the present invention, an adaptive detection method for industrial circuit oscillation under variable working conditions is specifically as follows:

[0074] S1: Obtain the setting signal SP, process variable signal PV and sampling frequency f in the industrial circuit to be detected s ; Through the comprehensive decomposition performance index, the process variable signal PV of the industrial loop is adaptively decomposed to obtain K modal IMFs k (k=1,2,…,K).

[0075] As a preferred embodiment of the present invention, the objective function of adaptive variational mode decomposition is:

[0076]

[0077] Among them, {u k} is the set of modal IMFs, {ω k} is the corresponding set of center frequencies, K is the number of decomposed modes, x(t) is the process variable signal PV, u k (t) is the modal IMF k , δ(t) is the unit impulse function, i is the imaginary unit, * represents convolution, It means to find partial derivative, t is time, ω k is the modal IMF k The center frequency of

[0078] As a preferred embodiment of the present invention, the steps of adaptive variational mode decomposition are as follows:

[0079] S11: Initialize the number of decomposed modes K=2;

[0080] S12: Decompose PV using adaptive variational mode decomposition (VMD) to obtain K modal IMFs k and the corresponding center frequency ω k , where k = 1, 2, ..., K;

[0081] S13: Judgment k The maximum value max(ω k ) is greater than f s / 2.56; if so, the maximum number of decomposable modes is K max =K-1, otherwise K=K+1 and repeat S12;

[0082] S14: Determine the maximum number of decomposable modes K max Is it greater than 1? If so, S15 is performed, otherwise adaptive variational mode decomposition is not performed;

[0083] S15: Set the value range of the decomposed mode number K to [2, K max ];

[0084] S16: Calculate the comprehensive decomposition performance index CPI under different modal numbers K K , where K = 2, 3, ..., K max ;

[0085] S17: Get the optimal number of decomposed modes K best , meet CPI Kbest =max(CPI K );

[0086] S18: Setting K best For the number of decomposed modes, adaptive variational mode decomposition is performed on PV to obtain K best Modal IMF k , where k = 1, 2, ..., K best .

[0087] In addition, as a preferred embodiment of the present invention, in step S1 of the present invention, the comprehensive decomposition performance index CPI K The calculation method is as follows:

[0088]

[0089] Among them, MI K is the average mutual information entropy between K modes and PV, CV K is the coefficient of variation of the center frequencies of the K modes, ΔE K is the residual energy between the energy sum of K modes and the PV energy, SE K It is the smallest sample entropy among the sample entropies of K modes.

[0090] The calculation formula is as follows:

[0091]

[0092] Among them, x belongs to the elements of the process variable signal x(t), and u belongs to the modal IMF k elements of , p(x,u) is the process variable signal and the modal IMF k The joint probability density function of p(x) is the marginal probability density function of the process variable signal, and p(u) is the modal IMF. k The marginal probability density function of ;

[0093] CV K The calculation formula is as follows:

[0094]

[0095] ΔE K The calculation formula is as follows:

[0096]

[0097] Among them, E x(t) is the energy of the process variable signal x(t), E x(t) =∫x(t) 2 dt, is the process variable signal u k (t) energy,

[0098] SE K The calculation formula is as follows:

[0099] SE K =min(SampEn(u k (t))),k=1,2,...,K

[0100] Among them, SampEn() is the sample entropy function. The specific solution process is based on the literature: Richman JS, Moorman J R. Physiological time-series analysis using approximate entropy and sampleentropy [J]. American journal of physiology-heart and circulatory physiology, 2000, 278 (6): H2039-H2049.

[0101] S2: Calculate the relative growth rate δ of the correlation coefficient k ; If δ k If it is greater than the preset value a1, the mode is judged to be important and S3 is performed; otherwise, the corresponding mode is eliminated.

[0102] As a preferred embodiment of the present invention, in step S2 of the present invention, the relative growth rate δ of the correlation coefficient is k The specific calculation steps are as follows:

[0103] S21: Modal IMF obtained by adaptive variational mode decomposition k , according to the descending order of the correlation coefficient with the process variable signal, the IMF is obtained k ;

[0104] S22: According to the new order, the first k modes are superimposed to obtain the superimposed signal sc k ;in, k = 1, 2, 3, ... K;

[0105] S23: Calculate the superposition signal sc k Correlation coefficient with process variable signal:

[0106]

[0107] Among them, Cov is the covariance, σ x(t) is the standard deviation of the process variable signal, is the standard deviation of the superimposed signal;

[0108] S24: Calculate the relative growth rate δ of the correlation coefficient k ; If δ k >a1, the mode is judged to be important and S3 is performed, otherwise the corresponding mode is eliminated; a1 is preferably selected as 0.05; relative growth rate δ k The expression is as follows:

[0109]

[0110] S3: Calculate the sparse index SI of important modes k If SI k If it is greater than the preset value a2, the mode proceeds to S4, otherwise the corresponding mode is eliminated.

[0111] As a preferred embodiment of the present invention, in step S3 of the present invention, the sparse index SI k The calculation steps are as follows:

[0112] S31: Standardization using the Z-score method for important modes;

[0113] S32: Based on the sampling frequency f s , the amplitude spectrum of the corresponding mode and the amplitude vector composed of amplitudes are obtained through fast Fourier transform;

[0114] S33: Calculate the sparse index SI of the amplitude vector k ; If there is a modal SI k If it is greater than the preset value a2, the mode is processed in S4, otherwise the corresponding mode is eliminated; a2 is preferably selected as 0.58; the sparse index SI k The expression is as follows:

[0115]

[0116] in, The amplitude vector is obtained from the amplitude spectrum of the kth mode, N is the vector length, and f represents the fth sample point.

[0117] S4: Estimate the confidence limits of the regularity index of the remaining modes after S3 screening If the upper confidence limit of the regularity index is If it is greater than the preset value a3, it is determined that the mode may oscillate and S5 is performed; otherwise, the corresponding mode is eliminated.

[0118] As a preferred embodiment of the present invention, in step S4 of the present invention, the confidence limits of the regularity index are The calculation steps are as follows:

[0119] S41: Calculate the autocorrelation function ACF of the modal IMF;

[0120] S42: Obtain the position of the ACF zero crossing point;

[0121] S43: Eliminate zero crossing points within the oscillation range of ±1.96 / sqrt(Id); where I is the data length of the IMF and d is the data interval; retain the positions z(j) of the remaining zero crossing points;

[0122] S44: Obtain the periodic sequence T(j)=2(z(j+1)-z(j)) and the average period through z(j) Where j = 1, 2, ..., L; the average period Denoted as modal IMF k The oscillation period T IMFk ;

[0123] S45: Calculate the confidence limits of regularity indicators Determine the upper confidence limit Is it greater than the preset value a3? If so, the mode may oscillate and S5 is performed, otherwise the corresponding mode is eliminated; a3 is preferably selected as 3; the upper and lower confidence limits of the regularity index The expression is as follows:

[0124]

[0125] Where, σ is the standard deviation of the periodic sequence T(j); is the average value of the periodic sequence T(j); α is generally a small positive real number, usually 0.05; is the 100α / 2 percentile of the chi-square distribution with L-1 degrees of freedom.

[0126] S5: Calculate the relative error e between the oscillation period of the remaining modes after S4 screening and the oscillation period of SP k If e k If it is greater than the preset value a4, the mode proceeds to S6, otherwise the corresponding mode is eliminated.

[0127] As a preferred embodiment of the present invention, in step S5 of the present invention, the relative error e k The calculation steps are as follows:

[0128] S51: Calculate the sparse index of SP and make a judgment; if the sparse index is greater than 0.58, proceed to S52; if the sparse index is less than 0.58, set its oscillation period to 0;

[0129] S52: Estimate the confidence upper limit R of the regularity index of SP 2_SP If R 2_SP >3, calculate the oscillation period T of the SP signal according to steps S41 to S44 SP , otherwise set its oscillation period to 0;

[0130] S53: Calculate the oscillation period of the residual mode The oscillation period T of SP SP The relative error of k If e k If it is greater than the preset value a4, S6 is performed, otherwise the corresponding mode is eliminated; a4 is preferably selected as 0.1; the relative error e k The expression is as follows:

[0131]

[0132] S6: Determine the lower confidence limit of the regularity index of the remaining modes after S5 screening Is it greater than the preset value a3? If so, report that the mode is oscillating and output its oscillation period and oscillation amplitude; otherwise, report that the corresponding mode may oscillate and output the oscillation possibility, oscillation period, and oscillation amplitude.

[0133] As a preferred embodiment of the present invention, in step S6 of the present invention, the oscillation period T IMFk Obtained through step S44; the oscillation amplitude is obtained by Hilbert transform to obtain the average value b1 of the upper envelope and the average value b2 of the lower envelope, and the oscillation amplitude = (b1-b2) / 2.

[0134] In order to better demonstrate the specific implementation and technical effect of the present invention, an adaptive detection method for industrial circuit oscillation under variable working conditions shown in steps S1 to S6 in the above preferred implementation is applied to a specific example.

[0135] Example

[0136] The specific implementation process of the adaptive detection method for industrial circuit oscillation under variable working conditions adopted in this embodiment is as described above and will not be repeated here.

[0137] Taking a flow control loop in a chemical industry as an example, the following describes in detail the oscillation detection method for a chemical process with valve sticking. For this loop, it is known a priori that the system has oscillations caused by valve sticking and the sampling frequency f sThe frequency is 0.1 Hz. The process variable signal PV and the setting signal SP are from the 35th chemical loop in the book Jelali M, Huang B. Detection and diagnosis of stiction in control loops: state of the art and advanced methods [M]. Springer London, 2009. The process variable signal data chemical.loop35.PV and the setting signal chemical.loop35.SP data can be obtained from the website https: / / sites.ualberta.ca / ~bhuang / Stiction-Book.htm.

[0138] The process variable signal PV and the setting signal SP collected in this embodiment are as follows: Figure 2 As shown, Figure 2 The horizontal axis is the sampling point, the sampling frequency is 0.1Hz, and the vertical axis is the flow rate. The modes obtained by adaptive decomposition are as follows: Figure 3 As shown, two modes IMF1 and IMF2 are extracted, and the relevant parameters for oscillation judgment are calculated, as shown in Table 1. The diagnosis result shows that mode IMF1 is an oscillation mode, and the oscillation amplitude and oscillation period of the oscillation mode are output. This embodiment demonstrates the advancement and effectiveness of the method of the present invention.

[0139] Table 1

[0140]

[0141] Starting from the perspective of loop oscillation under variable working conditions, the present invention can effectively eliminate the interference of non-stationary factors and periodic demands, and can adaptively select appropriate algorithm parameters for oscillation detection based on the input signal, thereby improving the automation and accuracy of industrial loop oscillation detection, and providing a methodological basis for online monitoring and performance evaluation of industrial loops.

[0142] The above-described embodiment is only a preferred solution of the present invention, but it is not intended to limit the present invention. A person skilled in the relevant technical field may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, any technical solution obtained by equivalent replacement or equivalent transformation falls within the protection scope of the present invention.

Claims

1. An adaptive detection method for industrial circuit oscillation under variable operating conditions, characterized in that: The details are as follows: S1: Obtain the setting signal SP, process variable signal PV and sampling frequency f in the industrial circuit to be detected s ; After adaptive variational mode decomposition of PV by comprehensive decomposition performance index, K modal IMFs are obtained k , k=1,2,…,K; S2: Calculate the relative growth rate δ of the correlation coefficient k ; If δ k If it is greater than the preset value a1, the mode is judged to be important and S3 is performed; otherwise, the corresponding mode is eliminated; S3: Calculate the sparse index SI of important modes k If SI k If it is greater than the preset value a2, the mode will proceed to S4, otherwise the corresponding mode will be eliminated; S4: Estimate the confidence limits of the regularity index of the remaining modes after S3 screening If the upper confidence limit of the regularity index is If it is greater than the preset value a3, the mode is judged to be oscillating and S5 is performed; otherwise, the corresponding mode is eliminated; S5: Calculate the relative error e between the oscillation period of the remaining modes after S4 screening and the oscillation period of SP k If e k If it is greater than the preset value a4, the mode proceeds to S6, otherwise the corresponding mode is eliminated; S6: Determine the lower confidence limit of the regularity index of the remaining modes after S5 screening Is it greater than the preset value a3? If so, report that the mode is oscillating and output its oscillation period and oscillation amplitude; otherwise, report that the corresponding mode may oscillate and output the oscillation possibility, oscillation period, and oscillation amplitude.

2. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 1, characterized in that: In S1, the objective function of adaptive variational mode decomposition is: Among them, {u k } is the set of modal IMFs, {ω k } is the corresponding set of center frequencies, K is the number of decomposed modes, x(t) is the process variable signal PV, u k (t) is the modal IMF k , δ(t) is the unit impulse function, i is the imaginary unit, * represents convolution, It means to find partial derivative, t is time, ω k is the modal IMF k The center frequency of 3. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 2, characterized in that: The steps of the adaptive variational mode decomposition are as follows: S11: Initialize the number of decomposed modes K=2; S12: Decompose PV using adaptive variational mode decomposition to obtain K modal IMFs k and the corresponding center frequency ω k , where k = 1, 2, ..., K; S13: Judgment k The maximum value max(ω k ) is greater than f s / 2.56; if so, the maximum number of decomposable modes is K max =K-1, otherwise K=K+1 and repeat S12; S14: Determine the maximum number of decomposable modes K max Is it greater than 1? If so, S15 is performed, otherwise adaptive variational mode decomposition is not performed; S15: Set the value range of the decomposed mode number K to [2, K max ]; S16: Calculate the comprehensive decomposition performance index CPI under different modal numbers K K , where K = 2, 3, ..., K max ; S17: Get the optimal number of decomposed modes K best , meet CPI Kbest =max(CPI K ); S18: Setting K best For the number of decomposed modes, adaptive variational mode decomposition is performed on PV to obtain K best Modal IMF k , where k = 1, 2, ..., K best .

4. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 3 is characterized in that: The comprehensive decomposition performance index CPI K The calculation method is as follows: in, is the average mutual information entropy between K modes and PV, CV K is the coefficient of variation of the center frequencies of the K modes, ΔE K is the residual energy between the energy sum of K modes and the PV energy, SE K It is the smallest sample entropy among the sample entropies of K modes.

5. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 4, characterized in that: Said The calculation formula is as follows: Among them, x belongs to the elements of the process variable signal x(t), and u belongs to the modal IMF k elements of , p(x,u) is the process variable signal and the modal IMF k The joint probability density function of p(x) is the marginal probability density function of the process variable signal, and p(u) is the modal IMF. k The marginal probability density function of ; CV K The calculation formula is as follows: ΔE K The calculation formula is as follows: Among them, E x(t) is the energy of the process variable signal x(t), E x(t) =∫x(t) 2 dt, is the process variable signal u k (t) energy, SE K The calculation formula is as follows: LAUGH K =min(CompEn(u k (t))),k=1,2,...,K Among them, SampEn() is the sample entropy function.

6. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 1, characterized in that: In S2, the relative growth rate of the correlation coefficient δ k The specific calculation steps are as follows: S21: Modal IMF obtained by adaptive variational mode decomposition k , according to the descending order of the correlation coefficient with the process variable signal S22: According to the new order, the first k modes are superimposed to obtain the superimposed signal sc k ;in, S23: Calculate the superposition signal sc k Correlation coefficient with process variable signal: Among them, Cov is the covariance, σ x(t) is the standard deviation of the process variable signal, is the standard deviation of the superimposed signal; S24: Calculate the relative growth rate δ of the correlation coefficient k ; If δ k >a1, the mode is judged to be important and S3 is performed, otherwise the corresponding mode is eliminated; a1 is selected as 0.05; relative growth rate δ k The expression is as follows:

7. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 1, characterized in that: In S3, the sparse index SI k The calculation steps are as follows: S31: Standardization using the Z-score method for important modes; S32: Based on the sampling frequency f s , the amplitude spectrum of the corresponding mode and the amplitude vector composed of amplitudes are obtained through fast Fourier transform; S33: Calculate the sparse index SI of the amplitude vector k ; If there is a modal SI k If it is greater than the preset value a2, the mode is processed in S4, otherwise the corresponding mode is eliminated; a2 is selected as 0.58; the sparse index SI k The expression is as follows: in, The amplitude vector is obtained from the amplitude spectrum of the kth mode, N is the vector length, and f represents the fth sample point.

8. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 1, characterized in that: In S4, the confidence limits of the regularity index The calculation steps are as follows: S41: Calculate the autocorrelation function ACF of the modal IMF; S42: Obtain the position of the ACF zero crossing point; S43: Eliminate zero crossing points within the oscillation range of ±1.96 / sqrt(Id); where I is the data length of the IMF and d is the data interval; retain the positions z(j) of the remaining zero crossing points; S44: Obtain the periodic sequence T(j)=2(z(j+1)-z(j)) and the average period through z(j) Where j = 1, 2, ..., L; the average period Denoted as modal IMF k The oscillation period S45: Calculate the confidence limits of regularity indicators Determine the upper confidence limit Is it greater than the preset value a3? If so, the mode may oscillate and S5 is performed, otherwise the corresponding mode is eliminated; a3 is selected as 3; the upper and lower confidence limits of the regularity index The expression is as follows: Where, σ is the standard deviation of the periodic sequence T(j); is the average value of the periodic sequence T(j); α is 0.05; is the 100α / 2 percentile of the chi-square distribution with L-1 degrees of freedom.

9. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 8, characterized in that: In S5, the relative error e k The calculation steps are as follows: S51: Calculate the sparse index of SP and make a judgment; if the sparse index is greater than 0.58, proceed to S52; if the sparse index is less than 0.58, set its oscillation period to 0; S52: Estimate the confidence upper limit R of the regularity index of SP 2_SP If R 2_SP >3, calculate the oscillation period T of the SP signal according to S41~S44 SP , otherwise set its oscillation period to 0; S53: Calculate the oscillation period of the residual mode The oscillation period T of SP SP The relative error of k If e k If it is greater than the preset value a4, S6 is performed, otherwise the corresponding mode is eliminated; a4 is selected as 0.1; the relative error e k The expression is as follows:

10. The method for adaptively detecting industrial circuit oscillation under variable working conditions according to claim 8, characterized in that: In S6, the oscillation period Obtained through S44; The oscillation amplitude is obtained by Hilbert transform to obtain the average value b1 of the upper envelope and the average value b2 of the lower envelope, and the oscillation amplitude = (b1-b2) / 2;

Images