An adaptive detection method for industrial circuit oscillation under varying operating conditions

Abstract

The present invention discloses an adaptive detection method for industrial circuit oscillation under variable working conditions, which belongs to the field of industrial circuit performance evaluation. The method includes: performing adaptive variational modal decomposition on the PV of the industrial circuit; calculating the relative growth rate of the correlation coefficient and selecting important modes; calculating the sparsity index of the important modes and eliminating the noise modes; estimating the upper and lower confidence limits of the regularity index of the remaining modes and judging the possibility of oscillation; calculating the relative error between the oscillation period of the remaining modes and the oscillation period of the SP and eliminating the interference of periodic demand; using the lower confidence limit of the regularity index of the remaining modes to judge and report the specific oscillation situation. The method proposed by the present invention starts from the perspective of circuit oscillation under variable working conditions, can effectively eliminate the interference of non-stationary factors and periodic demand, improve the automation level and accuracy of industrial circuit oscillation detection, and provide a methodological basis for online monitoring and performance evaluation of industrial circuits.

Description

Technical Field

[0001] The present 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. Degraded control system performance can lead to product quality fluctuations, increased production costs, and reduced economic efficiency. In severe cases, it can even threaten the stability and security of the entire system. During operation, systems often experience variable operating conditions, such as switching between set operating conditions or changes in load. This introduces non-stationarity. Furthermore, the system itself is subject to random noise and process disturbances, resulting in extremely non-stationary system signals under variable operating conditions. Signals often contain trends, noise, and multiple frequency components. This severely impacts oscillation detection results and can even lead to misdiagnosis. However, most oscillation detection methods fail to account for the variable operating conditions of industrial processes when dealing with non-stationary oscillations. Furthermore, most existing methods are unable to adaptively detect oscillations based on the input signal. For example, the number of modes and zero-crossing points in variational mode decomposition must be pre-set, hindering practical detection deployment and preventing the achievement of highly automated requirements. 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 in 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 performing adaptive variational mode decomposition on PV by comprehensive decomposition performance indicators, 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 sparsity 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 have the possibility of oscillation 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 will proceed to S6, otherwise the corresponding mode will be 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 center frequency set, 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, Indicates 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 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, proceed to S15, otherwise do not perform adaptive variational mode decomposition;

[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 , meeting CPI Kbest =max(CPI K );

[0022] S18: Setting K best For the number of decomposed modes, perform adaptive variational mode decomposition on PV and 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 of K modes and PV, CV K is the coefficient of variation of the center frequencies of the K modes, ΔE K is the residual between the sum of the energies of the K modes and the PV energy, SE K It is the minimum 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, p(x,u) is the process variable signal and 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] As an advantage, in S2, the relative growth rate δ of the correlation coefficient is k The specific calculation steps are as follows:

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

[0039] S22: According to the new order, superimpose the first k modes 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 the 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 the value 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, where N is the vector length and f represents the fth sample point.

[0051] As a preference, in said S4, the upper and lower 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 position 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 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 taken as 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 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, the oscillation period is set to 0;

[0062] S53: Calculate the oscillation period T of the residual mode IMFk The oscillation period T of SP SP The relative error e k If e k If it is greater than the preset value a4, then 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] The present invention starts from the perspective of loop oscillation under variable working conditions, can effectively eliminate the interference of non-stationary factors and periodic demand, and can adaptively select appropriate algorithm parameters for oscillation detection based on the input signal, thereby improving the degree of 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 Schematic diagram of the process of the present invention;

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

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

[0071] In order to make the above-mentioned objects, features and advantages of the present invention more clearly understood, the specific embodiments of the present invention are described in detail below with reference to 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 than 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 the various embodiments of the present invention can be combined accordingly without conflicting with each other.

[0072] like Figure 1 As shown in the figure, 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 oscillation possibility, 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 level 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, a method for adaptively detecting industrial circuit oscillation under variable operating 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 center frequency set, 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, Indicates 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 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, proceed to S15, otherwise do not perform adaptive variational mode decomposition;

[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: Obtain the optimal number of decomposed modes K best , meeting CPI Kbest =max(CPI K );

[0086] S18: Setting K best For the number of decomposed modes, perform adaptive variational mode decomposition on PV and 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 of K modes and PV, CV K is the coefficient of variation of the center frequencies of the K modes, ΔE K is the residual between the sum of the energies of the K modes and the PV energy, SE K It is the minimum 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, p(x,u) is the process variable signal and 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] 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 , and the IMF is obtained by arranging them in descending order according to the correlation coefficient with the process variable signal. k ;

[0104] S22: According to the new order, superimpose the first k modes 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 sparsity 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 the 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 the value 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, where 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, the mode is judged to have the possibility of oscillation 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 position 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 bound 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, the oscillation period is set to 0;

[0130] S53: Calculate the oscillation period of the residual mode The oscillation period T of SP SP The relative error e k If e k If it is greater than the preset value a4, then 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 effects of the present invention, the 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] The following is a flow control loop in a chemical industry as an example to describe 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 setpoint signal (SP) are from Chemical Loop 35 in the book Jelali M and 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 setpoint signal data (chemical.loop35.SP) are available at https: / / sites.ualberta.ca / ~bhuang / Stiction-Book.htm.

[0138] The process variable signal PV and 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 mode obtained by adaptive decomposition is as follows Figure 3 As shown, two modal IMFs, IMF1 and IMF2, were extracted, and the relevant parameters for oscillation determination were calculated, as shown in Table 1. The diagnostic results indicate that modal IMF1 is an oscillating mode, and the oscillation amplitude and oscillation period of the oscillating mode are output. This example demonstrates the advanced nature and effectiveness of the method of the present invention.

[0139] Table 1

[0140]

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

[0142] The embodiment described above is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Persons skilled in the art may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, any technical solution obtained by equivalent substitution or equivalent transformation falls within the scope of protection 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 performing adaptive variational mode decomposition on PV by comprehensive decomposition performance indicators, 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 sparsity 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 have the possibility of oscillation 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 will proceed to S6, otherwise the corresponding mode will be 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 center frequency set, 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, Indicates 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 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, proceed to S15, otherwise do not perform adaptive variational mode decomposition; 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: Obtain the optimal number of decomposed modes K best , meeting CPI Kbest =max(CPI K ); S18: Setting K best For the number of decomposed modes, perform adaptive variational mode decomposition on PV and 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, characterized in that: The comprehensive decomposition performance index CPI K The calculation method is as follows: in, is the average mutual information entropy of K modes and PV, CV K is the coefficient of variation of the center frequencies of the K modes, ΔE K is the residual between the sum of the energies of the K modes and the PV energy, SE K It is the minimum 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: described 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, p(x,u) is the process variable signal and 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, superimpose the first k modes 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 the 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 the value 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, where 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 position 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 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 taken as 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 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 e k If e k If it is greater than the preset value a4, then 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 the 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. The oscillation amplitude = (b1-b2) / 2;

Images