A Fault Diagnosis Method for Chemical Processes Based on Multivariate Correlation Analysis

Abstract

The invention discloses a chemical process fault diagnosis method based on multivariate correlation analysis, including constructing a sample data set, screening out variables with strong correlation with faults from the sample data set by using a ReliefF algorithm to form a key feature set; using an HVD decomposition algorithm to decompose the key feature set to obtain a number of IMF components, using the Spearman correlation coefficient to extract the IMF components with strong correlation with the key feature set, constructing an energy change rate feature vector after signal decomposition as a sample set; constructing a fault diagnosis model fused with TimesNet and OSRELM, using an improved NRBO algorithm to optimize model parameters, using a sample set to train the model; and using the trained model to perform fault diagnosis. The chemical process fault diagnosis method of the invention has high accuracy and reliability.

Description

Technical Field

[0001] The invention relates to a fault diagnosis method, and in particular to a chemical process fault diagnosis method based on multivariate correlation analysis. Background Art

[0002] Chemical processes are complex, and the equipment involved is diverse and the operating conditions are changeable. These factors make it difficult to accurately predict and diagnose chemical process faults. At present, chemical process fault diagnosis mainly includes model-based methods and data-based methods. Model-based methods rely on accurate mathematical models of chemical processes, but due to the nonlinearity and strong coupling of chemical processes, it is very difficult to establish accurate models. Data-based methods, such as principal component analysis (PCA) and support vector machines (SVM), can process complex chemical data to a certain extent, but they usually ignore the time correlation in the data and are difficult to accurately capture and predict process faults. In addition, when dealing with large-scale, high-dimensional, and strongly coupled chemical process data, existing fault diagnosis methods often face problems such as limited feature extraction capabilities (difficult to extract effective fault features from complex chemical data), insufficient model generalization capabilities (poor adaptability to different working conditions, difficult to be widely used in different chemical processes), and poor real-time performance (the speed of fault diagnosis cannot keep up with the needs of real-time monitoring). These problems make the accuracy and reliability of fault diagnosis poor, limiting the application effect of fault diagnosis in actual chemical production. Summary of the invention

[0003] Purpose of the invention: The purpose of the present invention is to provide a chemical process fault diagnosis method based on multivariate correlation analysis with high accuracy and reliability.

[0004] Technical solution: The chemical process fault diagnosis method based on multivariate correlation analysis of the present invention comprises:

[0005] (1) Collect data S D =[s1,s2,...,s n ], where s is the variable and n is the number of variables; a label l is added to each set of data to indicate a normal state or a specific fault state. z , is the number of label types; according to data S D and label l z Construct a sample dataset;

[0006] (2) Using the ReliefF feature selection algorithm, feature input variables that are strongly correlated with faults are selected from the sample data set to form a key feature set;

[0007] (3) The HVD decomposition algorithm is used to decompose the key feature set to obtain several IMF components. The Spearman correlation coefficient is used to extract the IMF components that have a strong correlation with the key feature set, and the energy change rate feature vector after signal decomposition is constructed as the sample set;

[0008] (4) Construct a fault diagnosis model TimesNet-OSRELM that integrates TimesNet and OSRELM, and use the sample set to train the TimesNet-OSRELM model;

[0009] (5) Use the trained TimesNet-OSRELM model to perform fault diagnosis and determine whether a fault has occurred and the fault type based on the fault diagnosis results and the set alarm threshold.

[0010] Furthermore, in step (1), the collected data includes process variable data and operating parameter data.

[0011] Furthermore, in step (2), the ReliefF feature selection algorithm screens feature input variables including:

[0012] (201) Initialization: Assign an initial weight to each feature in the sample data set;

[0013] (202) Iterative sampling: Perform multiple iterations and randomly select a sample point s in each iteration. g , in s g Find k nearest neighbor samples NHT in the same class, and find k nearest neighbor samples NMS in different classes;

[0014] (203) Update weights: For each feature, calculate s g The difference with NHT and NMS, and update the feature weights, the update formula is as follows:

[0015]

[0016] Among them, ω z+1 is the updated feature weight, ω z is the feature weight before updating, diff(z,s g ,NHT t ) indicates that feature z is at sample point s g and NHT t The difference between class(s g ) represents the sample point s g The fault type or status, P(class(s g )) represents the sample point s gThe prior probability of the category, P(C) is the prior probability of class C, NMS(C) is the tth nearest neighbor of the Cth different class, diff(z,s g , NMS(C) represents feature z at sample point s g and the difference on NMS(C); steps (201) to (203) are repeated m times;

[0017] (204) Calculate the average weight: After m iterations, calculate the average weight of each feature;

[0018] (205) According to the average weight of the features, the first H feature input variables with the largest weight are selected to form the key feature set.

[0019] Furthermore, in step (3), the HVD decomposition algorithm decomposes the key feature set to include:

[0020] (301) Calculate the instantaneous frequency of the maximum amplitude component:

[0021]

[0022] Wherein, F(t) represents a two-component non-stationary signal; a1(t) represents the instantaneous amplitude of the first component of the signal; a2(t) represents the instantaneous amplitude of the second component of the signal; i represents an imaginary unit; f1(t) represents the instantaneous frequency of the first component of the signal; f2(t) represents the instantaneous frequency of the second component of the signal; θ1 represents the initial phase of the first component of the signal; θ2 represents the initial phase of the second component of the signal;

[0023] When a1(t)>a2(t), the instantaneous frequency obtained by Hilbert transform is:

[0024]

[0025] Where, f(t) represents the instantaneous frequency; f1 represents the instantaneous frequency with the highest component amplitude, and f2 represents the instantaneous frequency of the component;

[0026] (302) Calculate the instantaneous amplitude and phase of the maximum amplitude component, multiply the original signal F(t) by two orthogonal signals respectively, and use a filter to filter out the oscillation part, so that the instantaneous amplitude a can be calculated. r (t) and instantaneous phase θ r :

[0027]

[0028] in, It represents the in-phase output signal obtained by filtering out the oscillating part with a low-pass filter; represents the quadrature phase output signal obtained by filtering out the oscillating part using a low-pass filter;

[0029] (303) Iterative calculation, by repeating the above steps (301) and (302), extracting the maximum amplitude component F1(t), and subtracting the two-component non-stationary signal F(t) from the maximum amplitude component F1(t), thereby obtaining a new initial signal F n-1 (t):

[0030] F n-1 (t) = F(t) - F1(t) (6)

[0031] Repeat steps (301) and (302) until an amplitude component that meets the conditions is obtained.

[0032] Furthermore, in step (3), extracting the IMF component using the Spearman correlation coefficient includes:

[0033] The key feature set is decomposed into q independent IMF components by HVD, and the undamaged signal E is calculated by equations (7) and (8): wp and the damaged signal E yp The optimal IMF component energy of is:

[0034]

[0035] Where T is the sampling time; F wp (t) is the IMF component of the signal decomposition in the undestructed state; F yp (t) is the IMF component of the signal decomposition under lossy state; p is the order of the optimal IMF component;

[0036] The optimal IMF energy before and after damage can be obtained from equations (7) and (8), and the energy change rate can be obtained from this:

[0037] The energy change rate feature vector after signal decomposition is constructed as the sample set, and divided into training set and test set.

[0038] Furthermore, in step (3), the number of iterations n or F of formula (6) n-1 The normalized standard deviation of (t) σ<0.001 is used as the termination condition of the iteration.

[0039] Furthermore, in step (4), the regularization parameter C and the ridge parameter L of the OSRELM model are optimized using the NRBO algorithm.

[0040] Furthermore, in step (4), the energy change rate feature vector is input into TimesNet, and TimesNet is used for feature selection to transform the input feature variable from a one-dimensional sequence X 1DGet deep features through the embedding layer That is, the input of the L-th layer TimesBlock is And the Lth layer extracts the two-dimensional temporal changes through 2D convolution, as shown in formula (10);

[0041]

[0042] in, is the output of TimesBlock of layer L;

[0043] The TimesBlock module in TimesNet is used to perform implicit feature extraction on the input feature vector, and the extracted data is sent to the input layer of OSRELM;

[0044] The Newton-Raphson optimization algorithm process is as follows:

[0045] ① Select the initial point x0: Randomly select an initial point x0 from the parameter definition domain as the starting point of the algorithm;

[0046] ② Calculate the fitness function f(x) and its derivative f'(x): Calculate the value of the function f(x) and the value of the first-order derivative f'(x) at the initial point x0; these values ​​will be used to update the iteration point;

[0047] ③ Update the iteration point: Use the Newton-Raphson formula to iterate and calculate the next iteration point x n+1 :

[0048]

[0049] where x n is the approximate root of the nth iteration; f(x n+1 ) is the function in x n The value at f'(x n+1 ) is the function in x n The derivative value at ;

[0050] ④ Check convergence conditions: When the maximum number of iterations of the algorithm is reached or the classification accuracy is met, the algorithm is considered to be converged and the iteration is stopped;

[0051] ⑤ Repeat iteration: If the stopping condition is not met, repeat steps ② and ③ until the stopping condition is met;

[0052] ⑥Output results: When the stopping condition is met, the optimal regularization parameter C and ridge parameter L are output.

[0053] Furthermore, in step (4), the elite reverse dynamic boundary strategy and adaptive sine-cosine search strategy are introduced to improve the NRBO algorithm.

[0054] Furthermore, the elite reverse dynamic boundary strategy includes:

[0055] 1) The elite reverse dynamic boundary continuously narrows the search space, improves the search space introduced in the position update, and improves the optimization efficiency of NRBO. The expression is as follows:

[0056]

[0057] In the formula, x n+1 is the position of the individual after iteration; r1∈(0,1) is a random number; ub is the upper limit of the optimized parameter; lb is the lower limit of the optimized parameter; is the optimal position of the individual;

[0058] 2) Adaptive sine and cosine search strategies include:

[0059] In the later stage of the Newton-Raphson optimization algorithm, the sine and cosine search strategy is added to disturb the optimal solution to obtain a new solution, reducing the possibility of the original algorithm falling into the local optimal solution. The formula is as follows:

[0060]

[0061] w(n)=w max *(w max -w min )*(Mit-n) / Mit(14)

[0062] in, is the position after the individual disturbance, sin(n) represents the sine function, cos(n) represents the cosine function, n is the nth iteration number, w(n) is the adaptive inertia weight; P is the random probability; w max is the maximum inertia weight; w min is the minimum inertia weight; Mit is the maximum number of iterations.

[0063] The present invention adopts the improved Newton-Raphson optimization algorithm (NRBO), which is an improvement on the traditional Newton method. By introducing more efficient search strategies and optimization mechanisms, the parameters of the TimesNet-OSRELM model are precisely tuned. This step ensures that the model can maintain the best diagnostic performance under various working conditions, and also improves the robustness of the model to abnormal data.

[0064] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages:

[0065] (1) Deep feature extraction and optimized variable selection: The present invention uses the ReliefF algorithm to conduct in-depth analysis of massive historical data of chemical processes, screen out variables with strong fault correlation, effectively remove redundant information between variables, and enhance the accuracy and generalization ability of the fault diagnosis model. This step not only improves the sensitivity of fault detection, but also enhances the generalization ability of the model by removing redundant information between variables. Furthermore, the Hilbert vibration decomposition (HVD) is used to decompose the selected variables in the time-frequency domain, and combined with the Spearman correlation coefficient, the feature components with the strongest correlation with the key feature set are accurately extracted, providing strong feature support for subsequent fault diagnosis.

[0066] (2) Fault diagnosis model construction: Based on feature extraction, the present invention constructs a composite model TimesNet-OSRELM that integrates time series analysis and machine learning technology. The TimesNet module is specially designed for time series data and can effectively capture the dynamic changes and long-term dependencies in the chemical process, while the OSRELM module uses its online learning ability to quickly adapt to process changes and perform real-time fault diagnosis. The construction of this model makes fault diagnosis no longer rely on empirical judgment, but based on scientific data-driven methods, greatly improving the accuracy and reliability of diagnosis.

[0067] (3) Real-time fault diagnosis and early warning: The present invention can respond to data changes in the chemical process in real time, quickly and accurately diagnose potential faults, and issue early warning signals in a timely manner. This real-time performance is due to the efficient computing power of the TimesNet-OSRELM model and the fast convergence characteristics of the NRBO algorithm. Through real-time monitoring and early warning, the present invention can minimize production losses caused by faults, and also provide a scientific basis for equipment maintenance. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 It is a flowchart of a chemical process fault diagnosis method based on multivariate correlation analysis provided by an embodiment of the present invention;

[0069] Figure 2 Schematic diagram of the TimesNet-OSRELM model structure in an embodiment of the present invention. DETAILED DESCRIPTION

[0070] The present invention will be further described below in conjunction with the accompanying drawings.

[0071] like Figure 1 As shown, an embodiment of the present invention provides a chemical process fault diagnosis method based on multivariate correlation analysis, comprising the following steps:

[0072] (1) Collect data S D=[s1,s2,...,s n ], where s is the variable and n is the number of variables; for each set of data S D Add labels indicating normal status or specific fault status z , is the number of label types; according to data S D and label l z Construct a sample dataset;

[0073] In step (1), historical data of the chemical process is collected by sensors and detection equipment, and the collected data include process variable data and operating parameter data. Taking the penicillin fed-batch fermentation process as an example, the process variable data include fermentation liquid volume, biomass concentration, penicillin concentration, substrate concentration, pH value, temperature, stirring rate, ventilation rate, stirring power, etc., and the operating parameter data include feeding rate, intermediate discharge interval time, dilution rate, metabolic water generation rate and water evaporation rate. Fault types include substrate concentration fault, pH and temperature control fault, agitator fault, accelerated flow concentration fault, etc.

[0074] (2) Using the ReliefF feature selection algorithm, we can select the characteristic input variables that are strongly correlated with the fault from the sample data set, remove the redundant information between the variables, and form a key feature set to improve the accuracy and efficiency of the fault diagnosis model;

[0075] The implementation process of the ReliefF feature selection algorithm is as follows:

[0076] (201) Initialization: Assign an initial weight to each feature in the sample data set;

[0077] (202) Iterative sampling: Perform multiple iterations and randomly select a sample point s in each iteration. g , in s g Find k nearest neighbor samples NHT in the same class, and find k nearest neighbor samples NMS in different classes;

[0078] (203) Update weights: For each feature, calculate s g The difference with NHT and NMS, and update the feature weights, the update formula is as follows:

[0079]

[0080] Among them, ω z+1 is the updated feature weight, ω z is the feature weight before updating, diff(z,s g ,NHT t ) indicates that feature z is at sample point s gand NHT t The difference between class(s g ) represents the sample point s g The fault type or status, P(class(s g )) represents the sample point s g The prior probability of the category, P(C) is the prior probability of class C, NMS(C) is the tth nearest neighbor of the Cth different class, diff(z,s g , NMS(C) represents feature z at sample point s g and the difference on NMS(C); steps (201) to (203) are repeated m times;

[0081] (204) Calculate the average weight: After m iterations, calculate the average weight of each feature;

[0082] (205) According to the average weight of the features, the first H feature input variables with the largest weight are selected to form the key feature set.

[0083] (3) Using the HVD (Hilbert vibration) decomposition algorithm, the key feature set of step (2) is decomposed to obtain several IMF components to reduce the computational difficulty. The Spearman correlation coefficient is used to extract the IMF components that have a strong correlation with the key feature set, and the energy change rate feature vector after signal decomposition is constructed as the sample set;

[0084] HVD decomposition generally uses an iterative method to decompose a non-stationary signal into several components of different amplitudes, and the time-frequency information of each component in the initial signal is obtained by iterative calculation. Taking a two-component non-stationary signal as an example, the decomposition process is as follows.

[0085] (301) Calculate the instantaneous frequency of the maximum amplitude component:

[0086]

[0087] Wherein, F(t) represents a two-component non-stationary signal; a1(t) represents the instantaneous amplitude of the first component of the signal; a2(t) represents the instantaneous amplitude of the second component of the signal; i represents an imaginary unit; f1(t) represents the instantaneous frequency of the first component of the signal; f2(t) represents the instantaneous frequency of the second component of the signal; θ1 represents the initial phase of the first component of the signal; θ2 represents the initial phase of the second component of the signal;

[0088] When a1(t)>a2(t), the instantaneous frequency obtained by Hilbert transform is:

[0089]

[0090] Where, f(t) represents the instantaneous frequency; f1 represents the instantaneous frequency with the highest component amplitude, and f2 represents the instantaneous frequency of the component;

[0091] It can be seen from formula (3) that the instantaneous frequency f(t) is obtained by adding the instantaneous frequency f1 with the highest component amplitude to the transformed high-frequency part corresponding to f1. The instantaneous frequency f1 with the highest amplitude is usually obtained by filtering with a low-pass filter.

[0092] (302) Calculate the instantaneous amplitude and phase of the maximum amplitude component, multiply the original signal F(t) by two orthogonal signals respectively, and use a low-pass filter to filter out the oscillation part, thereby calculating the instantaneous amplitude a r (t) and instantaneous phase θ r :

[0093]

[0094]

[0095] in, It represents the in-phase output signal obtained by filtering out the oscillating part with a low-pass filter; represents the quadrature phase output signal obtained by filtering out the oscillating part using a low-pass filter;

[0096] (303) Iterative calculation, by repeating steps (301) and (302), extracting the maximum amplitude component F1(t), and subtracting the two-component non-stationary signal F(t) from the maximum amplitude component F1(t), thereby obtaining a new initial signal F n-1 (t):

[0097] F n-1 (t) = F(t) - F1(t) (6)

[0098] Repeat steps (301) and (302) until the amplitude component that meets the conditions is obtained. n-1 The normalized standard deviation of (t) σ<0.001 is used as the termination condition of the iteration.

[0099] The key feature set is decomposed into several IMF components through HVD. In order to retain more features of the original signal, the Spearman correlation coefficient is used to perform correlation analysis on the decomposed IMF components to select the IMF component with the strongest correlation with the key feature set and extract the optimal IMF component energy eigenvector.

[0100] The key feature set is decomposed into q independent IMF components by HVD, and the undamaged signal E is calculated by equations (7) and (8): wp and the damaged signal E ypThe optimal IMF component energy of is:

[0101]

[0102] Where T is the sampling time; F wp (t) is the IMF component of the signal decomposition in the undestructed state; F yp (t) is the IMF component of the signal decomposition under lossy state; p is the order of the optimal IMF component.

[0103] The optimal IMF energy before and after damage can be obtained from equations (7) and (8), and the energy change rate can be obtained from this:

[0104]

[0105] Therefore, the energy change rate feature vector after signal decomposition is constructed as the sample set, and the training set and test set are divided and input into the subsequent TimesNet-OSRELM model for deep learning and fault diagnosis.

[0106] (4) Construct a fault diagnosis model TimesNet-OSRELM that integrates TimesNet and OSRELM (online sequential regularized learning machine), introduce the elite reverse dynamic boundary strategy and adaptive sine-cosine search strategy to improve the NRBO (Newton-Raphson optimization) algorithm, use the improved NRBO algorithm to optimize the regularization parameter C and ridge parameter L of the OSRELM model, and use the sample set generated in step (3) to train the TimesNet-OSRELM model to obtain the optimal model network parameters;

[0107] The chemical process fault diagnosis model TimesNet-OSRELM is constructed as follows Figure 2 As shown in the figure, the energy change rate feature vector is input into TimesNet, and the feature selection provided by TimesNet is used to convert the input feature variable from a one-dimensional sequence X 1D Get deep features through the embedding layer That is, the input of the L-th layer TimesBlock is And the Lth layer extracts two-dimensional temporal changes through 2D convolution, as shown in formula (10).

[0108]

[0109] in, It is the output of TimesBlock of the Lth layer.

[0110] In the TimesBlock module, firstly, the periodicity, local patterns and trends are mined based on the time dimension through the Fast Fourier Transform (FFT), which automatically learns and extracts the features in the data. Therefore, the TimesBlock module in TimesNet is used to perform implicit feature extraction on the input feature vector. The extracted data is fed into the input layer of OSRELM.

[0111] OSRELM is an online sequential regularized extreme learning machine, which mainly includes an input layer, a hidden layer and an output layer. The main idea is to combine the regularization method to further improve its generalization performance, and adopt a new update formula to eliminate the initialization stage. In order to enable OSRELM to adapt to new data effectively and reliably, an efficient cross-validation method is implemented. Finally, the matrix reconstruction method is used to solve the unstable update problem. Since the OSRELM model results are highly dependent on the regularization parameter C and the ridge parameter L. Therefore, the improved Newton-Raphson optimization algorithm is used to optimize these two parameters to improve the accuracy of the model.

[0112] The Newton-Raphson optimization algorithm process is as follows:

[0113] ① Select the initial point x0: Randomly select an initial point x0 from the parameter definition domain as the starting point of the algorithm;

[0114] ② Calculate the fitness function f(x) and its derivative f'(x): Calculate the value of the function f(x) and the value of the first-order derivative f'(x) at the initial point x0. These values ​​will be used to update the iteration point.

[0115] ③ Update the iteration point: Use the Newton-Raphson formula to iterate and calculate the next iteration point x n+1 :

[0116]

[0117] Among them, x n is the approximate root of the nth iteration; f(x n+1 ) is the function in x n+1 The value at f'(x n+1 ) is the function in x n+1 The derivative value at .

[0118] ④ Check convergence conditions: When the maximum number of iterations of the algorithm is reached or the classification accuracy is met, the algorithm is considered to be converged and the iteration is stopped.

[0119] ⑤ Repeat iteration: If the stop iteration condition is not met, repeat steps ② and ③ until the stop iteration condition is met.

[0120] ⑥Output results: When the conditions for stopping iteration are met, the optimal regularization parameter C and ridge parameter L are output.

[0121] Although the Newton-Raphson optimization algorithm performs well in many problems, it still has problems such as being easily trapped in local optimal solutions. Therefore, the elite reverse dynamic boundary strategy and adaptive sine-cosine search strategy are introduced to improve the original NRBO algorithm.

[0122] 1) Elite Reverse Dynamic Border Strategy

[0123] The elite reverse dynamic boundary continuously narrows the search space, improves the search space introduced in the position update, and improves the optimization efficiency of NRBO. The expression is as follows:

[0124]

[0125] In the formula, x n+1 is the position of the individual after iteration; r1∈(0,1) is a random number; ub is the upper limit of the optimized parameter; lb is the lower limit of the optimized parameter; is the optimal position of the individual;

[0126] 2) Adaptive sine and cosine search strategy

[0127] In the later stage of the Newton-Raphson optimization algorithm, the sine and cosine search strategy is added to disturb the optimal solution to obtain a new solution, reducing the possibility of the original algorithm falling into the local optimal solution. The formula is as follows:

[0128]

[0129] w(n)=w max *(w max -w min )*(Mit-n) / Mit(14)

[0130] in, is the position after the individual disturbance, sin(n) represents the sine function, cos(n) represents the cosine function, n is the nth iteration number, w(n) is the adaptive inertia weight; P is the random probability; w max is the maximum inertia weight; w min is the minimum inertia weight; Mit is the maximum number of iterations.

[0131] (5) Use the trained TimesNet-OSRELM model to perform fault diagnosis and determine whether a fault has occurred and the fault type based on the fault diagnosis results and the set alarm threshold.

[0132] Alarm threshold setting: Set the fault alarm threshold for each fault type according to the safety standards and operating requirements of the chemical plant.

[0133] Result evaluation: Compare the model’s diagnostic results with the preset alarm thresholds.

[0134] Fault type definition: If the diagnostic results show that the fault characteristics exceed the alarm threshold, the fault is classified according to the severity and type as follows:

[0135] Minor fault: The fault characteristics slightly exceed the threshold and may not require immediate shutdown, but close monitoring is required.

[0136] Moderate fault: The fault characteristics clearly exceed the threshold and may require adjustment of the production process or partial maintenance.

[0137] Serious fault: The fault characteristics far exceed the threshold and require immediate shutdown and comprehensive inspection.

[0138] Define specific fault types based on the specific indicators of fault characteristics and situations where the threshold is exceeded:

[0139] Sensor Fault: Abnormal sensor readings, possibly due to sensor damage or signal interference.

[0140] Equipment failure: The operating parameters of key equipment such as pumps and compressors are outside the normal range.

[0141] Process control failure: The process variable (such as temperature, pressure, flow) exceeds the predetermined control range.

[0142] Material imbalance: The quantity or quality of input or output material does not meet expectations.

[0143] The TimesNet-OSRELM model learns fault characteristics through training data and is able to predict the type of fault based on the input data. However, even after the model is trained, its prediction results are not always 100% accurate and there may be some uncertainty and error. The alarm threshold provides a mechanism to determine whether the model output is credible enough to trigger a fault alarm.

[0144] Through testing and verification, the present invention can accurately diagnose faults in chemical processes and timely determine the type of faults when they occur, thus providing strong support for fault prevention and maintenance.

Claims

1. A chemical process fault diagnosis method based on multivariate correlation analysis, characterized in that: include: (1) Collect data S D =[s1,s2,...,s n ], where s is a variable and n is the number of variables; a label l is added to each set of data to indicate a normal state or a specific fault state. z , is the number of label types; according to data S D and label l z Construct a sample dataset; (2) Using the ReliefF feature selection algorithm, feature input variables that are strongly correlated with faults are selected from the sample data set to form a key feature set; (3) The HVD decomposition algorithm is used to decompose the key feature set to obtain several IMF components. The Spearman correlation coefficient is used to extract the IMF components that have a strong correlation with the key feature set, and the energy change rate feature vector after signal decomposition is constructed as the sample set; (4) Construct a fault diagnosis model TimesNet-OSRELM that integrates TimesNet and OSRELM, and use the sample set to train the TimesNet-OSRELM model; (5) Use the trained TimesNet-OSRELM model to perform fault diagnosis and determine whether a fault has occurred and the fault type based on the fault diagnosis results and the set alarm threshold; In step (4), the regularization parameter C and the ridge parameter L of the OSRELM model are optimized using the NRBO algorithm; In step (4), the elite reverse dynamic boundary strategy and adaptive sine-cosine search strategy are introduced to improve the NRBO algorithm; the elite reverse dynamic boundary strategy includes: 1) The elite reverse dynamic boundary continuously narrows the search space, improves the search space introduced in the position update, and improves the optimization efficiency of NRBO. The expression is as follows: In the formula, x n+1 is the position of the individual after iteration; r1∈(0,1) is a random number; ub is the upper limit of the optimized parameter; lb is the lower limit of the optimized parameter; is the optimal position of the individual; 2) Adaptive sine and cosine search strategies include: In the later stage of the Newton-Raphson optimization algorithm, the sine and cosine search strategy is added to disturb the optimal solution to obtain a new solution, reducing the possibility of the original algorithm falling into the local optimal solution. The formula is as follows: in(n)=in max *(In max -In min )*(Myth-n) / Myth(14) in, is the position after the individual disturbance, sin(n) represents the sine function, cos(n) represents the cosine function, n is the nth iteration number, w(n) is the adaptive inertia weight; P is the random probability; w max is the maximum inertia weight; w min is the minimum inertia weight; Mit is the maximum number of iterations.

2. The chemical process fault diagnosis method based on multivariate correlation analysis according to claim 1 is characterized in that: In step (1), the collected data includes process variable data and operating parameter data.

3. The chemical process fault diagnosis method based on multivariate correlation analysis according to claim 1 is characterized in that: In step (2), the ReliefF feature selection algorithm selects feature input variables including: (201) Initialization: Assign an initial weight to each feature in the sample data set; (202) Iterative sampling: Perform multiple iterations and randomly select a sample point s in each iteration. g , in s g Find k nearest neighbor samples NHT in the same class, and find k nearest neighbor samples NMS in different classes; (203) Update weights: For each feature, calculate s g The difference with NHT and NMS, and update the feature weights, the update formula is as follows: Among them, ω z+1 is the updated feature weight, ω z is the feature weight before updating, diff(z,s g ,NHT t ) indicates that feature z is at sample point s g and NHT t The difference between class(s g ) represents the sample point s g The fault type or status, P(class(s g )) represents the sample point s g The prior probability of the category, P(C) is the prior probability of class C, NMS(C) is the tth nearest neighbor of the Cth different class, diff(z,s g , NMS(C) represents feature z at sample point s g and the difference on NMS(C); steps (201) to (203) are repeated m times; (204) Calculate the average weight: After m iterations, calculate the average weight of each feature; (205) According to the average weight of the features, the first H feature input variables with the largest weight are selected to form the key feature set.

4. The chemical process fault diagnosis method based on multivariate correlation analysis according to claim 1 is characterized in that: In step (3), the HVD decomposition algorithm decomposes the key feature set to include: (301) Calculate the instantaneous frequency of the maximum amplitude component: Wherein, F(t) represents a two-component non-stationary signal; a1(t) represents the instantaneous amplitude of the first component of the signal; a2(t) represents the instantaneous amplitude of the second component of the signal; i represents an imaginary unit; f1(t) represents the instantaneous frequency of the first component of the signal; f2(t) represents the instantaneous frequency of the second component of the signal; θ1 represents the initial phase of the first component of the signal; θ2 represents the initial phase of the second component of the signal; When a1(t)>a2(t), the instantaneous frequency obtained by Hilbert transform is: Among them, f(t) represents the instantaneous frequency; f1 represents the instantaneous frequency with the highest amplitude of the component, and f2 represents the instantaneous frequency of the component; (302) Calculate the instantaneous amplitude and phase of the maximum amplitude component, multiply the original signal F(t) by two orthogonal signals respectively, and use a filter to filter out the oscillation part, so that the instantaneous amplitude a can be calculated. r (t) and instantaneous phase θ r : in, It represents the in-phase output signal obtained by filtering out the oscillating part with a low-pass filter; represents the quadrature phase output signal obtained by filtering out the oscillating part using a low-pass filter; (303) Iterative calculation, by repeating the above steps (301) and (302), extracting the maximum amplitude component F1(t), and subtracting the two-component non-stationary signal F(t) from the maximum amplitude component F1(t), thereby obtaining a new initial signal F n-1 (t): F n-1 (t)=F(t)-F1(t)(6) Repeat steps (301) and (302) until an amplitude component that meets the conditions is obtained.

5. The chemical process fault diagnosis method based on multivariate correlation analysis according to claim 4 is characterized in that: In step (3), the Spearman correlation coefficient extracts the IMF components including: The key feature set is decomposed into q independent IMF components by HVD, and the undamaged signal E is calculated by equations (7) and (8): wp and the damaged signal E yp The optimal IMF component energy of: Where T is the sampling time; F wp (t) is the IMF component of the signal decomposition in the undestructed state; F yp (t) is the IMF component of the signal decomposition under lossy state; p is the order of the optimal IMF component; The optimal IMF energy before and after damage can be obtained from equations (7) and (8), and the energy change rate can be obtained from this: The energy change rate feature vector after signal decomposition is constructed as the sample set, and divided into training set and test set.

6. The chemical process fault diagnosis method based on multivariate correlation analysis according to claim 4 is characterized in that: In step (3), the number of iterations n or F n-1 The normalized standard deviation of (t) σ<0.001 is used as the termination condition of the iteration.

7. The chemical process fault diagnosis method based on multivariate correlation analysis according to claim 1 is characterized in that: In step (4), the energy change rate feature vector is input into TimesNet, and TimesNet is used for feature selection to transform the input feature variable from a one-dimensional sequence X 1D Get deep features through the embedding layer That is, the input of the L-th layer TimesBlock is And the Lth layer extracts the two-dimensional temporal changes through 2D convolution, as shown in formula (10); in, is the output of TimesBlock of layer L; The TimesBlock module in TimesNet is used to perform implicit feature extraction on the input feature vector, and the extracted data is sent to the input layer of OSRELM; The Newton-Raphson optimization algorithm process is as follows: ① Select the initial point x0: Randomly select an initial point x0 from the parameter definition domain as the starting point of the algorithm; ② Calculate the fitness function f(x) and its derivative f'(x): Calculate the value of the function f(x) and the value of the first-order derivative f'(x) at the initial point x0; these values ​​will be used to update the iteration point; ③ Update the iteration point: Use the Newton-Raphson formula to iterate and calculate the next iteration point x n+1 : where x n is the approximate root of the nth iteration; f(x n+1 ) is the function in x n The value at f'(x n+1 ) is the function in x n The derivative value at ; ④ Check convergence conditions: When the maximum number of iterations of the algorithm is reached or the classification accuracy is met, the algorithm is considered to be converged and the iteration is stopped; ⑤ Repeat iteration: If the stopping condition is not met, repeat steps ② and ③ until the stopping condition is met; ⑥Output results: When the stopping condition is met, the optimal regularization parameter C and ridge parameter L are output.

Images