CN111626099A - Industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM - Google Patents

Abstract

The invention discloses an industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM, which comprises the following steps: calculating mutual information between process output signals of two oscillation loops to be analyzed to perform characteristic selection, reserving loop signal pairs with relatively high correlation degree, removing noise and period items by using EMD and DFA methods, calculating cross mapping indexes of reconstructed sub-signal pairs under different time lags by using CCM methods, preliminarily judging the truth of causal relationship according to the positive and negative of the optimal time lag, calculating the cross mapping indexes of each target sub-signal pair under different sample lengths corresponding to the optimal time lag, judging whether convergence occurs according to a convergence threshold value, and finally obtaining a final causal relationship network of the oscillation loop so as to determine an oscillation propagation path and position an oscillation source. The invention can realize the rapid and accurate diagnosis of the oscillation propagation path in the industrial control system and the positioning of the oscillation source, and provides a new idea for carrying out fault diagnosis on the oscillation circuit of the industrial control system by using a causal analysis method.

Description

Industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM

Technical Field

The invention belongs to the field of fault diagnosis of industrial control systems, and particularly relates to an industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM.

Background

Due to controller over-tuning, valve sticking, external disturbances, and process non-linearities, approximately 30% of industrial control loops have oscillatory faults. Hunting is one of the main causes of the degradation of the control performance of the industrial control system. Continuous oscillation increases the abrasion of industrial equipment, reduces the product quality, seriously influences the production benefit of a factory and even generates potential safety hazard. There may be thousands of loops in a large industrial control system, and due to the interaction and high coupling between loops, the oscillation of any one loop may be transferred to other loops of the system, eventually resulting in global oscillation.

The method comprises the following steps of carrying out causal analysis on an oscillation circuit in the industrial process, determining an oscillation propagation path and positioning an oscillation source so as to take corresponding measures more pertinently, maintaining the damaged machine in time, shortening the shutdown time and reducing related losses. Therefore, an effective causal analysis method is designed for the oscillation circuit, the oscillation propagation path is diagnosed efficiently and accurately, the fault source is located accurately, and the method has important significance and reference value for fault diagnosis of the industrial control circuit.

In recent years, the diagnosis of industrial process faults based on causal analysis has been rapidly developed. Wakefield B J et al use Granger and the transfer entropy method to diagnose faults in a simulated milling circuit of a mineral processing plant. Bauer M et al estimate the time delay between process variables using cross-correlation functions, give a qualitative model of the oscillation propagation path in the form of a causal graph and use this method for Eastman chemical processes. Duan P et al propose the concept of Direct Transfer Entropy (DTE) for detecting the presence of a Direct flow of information between variables and a normalization method of differential Direct transfer entropy to determine the strength of a Direct causal relationship. However, the Granger method has a large model fitting error for nonlinear variables, unstable variables and variables in an inseparable system with complex relationships, and may present a causal relationship detection result of missing judgment and erroneous judgment. The transfer entropy method relates to the estimation of probability density function, and is sensitive to parameters and low in calculation efficiency.

The Convergent Cross Map (CCM) method provided by Sugihara can be used for an inseparable dynamic system with weak to medium coupling strength, the defect that the Granger method cannot be used for the inseparable system is overcome, and the method is simple in principle, easy to implement and high in calculation efficiency. However, the CCM method can only obtain its causal relationship by pairwise analysis of variables. If a large industrial process has global oscillation, in order to diagnose the oscillation, one method is to carry out pairwise causal analysis on all oscillation loops to obtain a related causal network. However, this approach has significant limitations: 1) the related variables are excessive, the calculated amount is extremely large, the process is complicated, and the efficiency is extremely low; 2) unnecessary cause and effect detection may be involved, which is time consuming; 3) because of a plurality of variables and large span, the analysis result is not necessarily accurate. In addition, the method is provided based on an ecosystem and is mainly used for ecology and neuroscience, and most industrial processes have the characteristics of high dimension, strong coupling and the like, which are obviously different from the ecosystem which is suitable for the traditional CCM method. Particularly, part of the oscillation signals are quasi-periodic sinusoidal signals, phase differences exist between the sinusoidal signals, and specific lag information between the sinusoidal signals and the sinusoidal signals cannot be determined directly according to the phase differences, so that the accuracy of causal analysis is reduced. The periodic signals may also cause the signals to become similar and coordinated, interfering with causal analysis. Besides periodic signals, unwanted noise, which is almost unavoidable in industrial processes, may also have an impact on the causal analysis.

Based on the background, aiming at multi-loop oscillation signals in an industrial control system, the CCM method is improved to greatly reduce unnecessary detection workload and improve causal analysis efficiency, meanwhile, interference items in the signals are eliminated, so that the accuracy and reliability of causal analysis results are greatly improved, the oscillation propagation path is quickly and accurately diagnosed, an oscillation source is positioned, and the method has important application value for fault diagnosis in the industrial control system.

Disclosure of Invention

The invention provides an industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM, which is suitable for diagnosis of an oscillation propagation path in a multi-loop oscillation signal and positioning of an oscillation source, has high analysis efficiency and high accuracy, and only needs to obtain relevant oscillation data of an industrial process without process mechanism knowledge.

An industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM comprises the following steps:

(1) collecting process output signals of oscillation loops of the industrial control systems to be analyzed;

(2) calculating mutual information values between output signals of two oscillating circuits in the process, and obtaining a target signal pair through feature selection;

(3) denoising and cycle item removing processing are carried out on signals in each target signal pair, and the rest part is reconstructed to obtain corresponding target sub-signals;

(4) determining the optimal embedding dimension of each target sub-signal, calculating cross mapping indexes of each target sub-signal pair under different time lags by using a CCM (continuous current measurement) method, and preliminarily judging the truth of a causal relationship according to the positive and negative of the optimal time lags;

(5) calculating cross mapping indexes of each target sub-signal pair corresponding to the optimal time lag under different sample lengths, and judging whether convergence occurs according to a convergence threshold value to obtain a causal relationship of each signal pair;

(6) and determining an oscillation propagation path and positioning an oscillation source according to the final causal relationship network.

The method comprises the steps of performing characteristic selection by calculating mutual information between process output signals of every two oscillation circuits to be analyzed, reserving a circuit signal pair with relatively high correlation degree, removing noise (optional) and period items by using EMD and DFA methods, calculating cross mapping indexes of a reconstructed sub signal pair under different time lags by using a CCM method, preliminarily judging the truth of a causal relationship according to the positive and negative of the optimal time lag, calculating the cross mapping indexes of each target sub signal pair under different sample lengths corresponding to the optimal time lag, judging whether convergence occurs according to a convergence threshold value, and finally obtaining a final causal relationship network of the oscillation circuit, thereby determining an oscillation propagation path and positioning an oscillation source.

Aiming at multi-loop oscillation signals in an industrial control system, the invention greatly reduces unnecessary causal analysis workload, improves analysis efficiency, eliminates interference items in the signals, greatly improves the accuracy and reliability of causal analysis results, can quickly and accurately diagnose oscillation propagation paths and position oscillation sources, provides a new idea for carrying out fault diagnosis on the oscillation loops of the industrial control system by using a causal analysis method, and has certain reference significance and application value.

The invention directly adopts measurable variables of the chemical process as process output signals, and all oscillation loop process output signals to be analyzed are acquired in real time on site.

The specific process of the step (2) is as follows:

(2-1) calculating a mutual information value between process output signals of every two oscillating circuits by using a K-order nearest neighbor method;

(2-2) carrying out 95% confidence test on the mutual information value by using a normalized substitution data method;

(2-3) retaining the signal pair having the mutual information value greater than 0 as the target signal pair.

The specific process of the step (3) is as follows:

(3-1) decomposing the target signal by using an EMD method to obtain a series of modal components Imf and residual components;

(3-2) calculating the DFA indexes of each layer of components by using a DFA method, if all the DFA indexes are more than 0.5, judging that the DFA indexes do not contain noise, and directly executing the step (3-3); if the DFA index of a certain layer is less than or equal to 0.5, judging the layer to be useless noise, and if multiple layers exist, accumulating the layers to be final noise;

(3-3) judging the period items in the modal components, and accumulating all the period items to be used as final period items;

and (3-4) if the signal contains noise, performing denoising and cycle item removing processing on the target signal, if the signal does not contain noise, directly performing cycle item removing processing on the target signal, and reconstructing the rest part after the processing to obtain a target sub-signal.

The specific process of the step (4) is as follows:

(4-1) determining the optimal embedding dimension E of each target sub-signal by using a G-curve method, wherein the determination standard is as follows: and when the G function is less than 0.5, the corresponding E is the optimal embedding dimension.

(4-2) determining two variables X and Y for each target sub-signal pair, setting a time lag λ ∈ [ -c: Δ c: c) of X with respect to Y as a preset cause variable and a result variable]. For each time lag, ρ between X (t + λ) and Y (t) can be calculated as (4-3) to (4-8)_X→Y(X (t) is replaced by X (t + lambda); similarly, rho between Y (t + lambda) and X (t)_Y→X(X (t) is replaced by Y (t + lambda) and Y (t) is replaced by X (t));

(4-3) for the target sub-signal pair X (t) and Y (t) with the sample length P, respectively establishing corresponding time lag coordinate vectors, wherein the calculation formula is as follows:

in the formula, E is the optimal embedding dimension determined in the step (4-1), and if the E of the target sub-signal pair is different, the maximum value of the E and the maximum value is selected as the final E; τ is the time lag, and is 1 by default;x(t)、y(t) sets of X (t), Y (t) corresponding shadow manifolds M_X、M_Y；

(4-4) finding the shadow manifold M_YPoint ofy(t) obtaining a distance sequence { d) in descending order from the E +1 nearest neighbors at the time t₁,d₂,…,d_E+1And a sequence of time indices t corresponding to the sequence of distances₁,t₂,…,t_E+1And (c) the step of (c) in which,xthe distance between (t) is calculated as follows:

d_i＝D(y(t_i),y(t_j))

wherein D (a, b) represents the euclidean distance between vectors a and b;

(4-5) calculating the related weights by using the E +1 nearest neighbors of each point y (t), wherein the weight calculation formula is as follows:

w_i＝u_i/U

in the formula (I), the compound is shown in the specification,

normalization factor

(4-6) calculating an estimated value of X by using the weight

The calculation formula is as follows:

in the formula, w_iIs a weight, t_iIs a time index corresponding to the weight;

(4-7) calculating X and

standard pearson correlation coefficient p between_X→YAs a cross-map correlation coefficient, the causal strength of X → Y is represented by the following calculation formula:

(4-8) pairing ρ with a substitution data method_X→YCarrying out 95% confidence test;

(4-9) obtaining a series of cross mapping indexes of the X → Y direction changing along with time lag by adopting the steps, and obtaining a series of cross mapping indexes of the Y → X direction changing along with time lag by using the same method as the steps (4-2) to (4-8); wherein the time lag corresponding to the maximum rho value is the optimal time lag lambda_X→Y、λ*_Y→XAnd (judging by a computer), and preliminarily judging the authenticity of the causal relationship according to the positive and negative of the optimal time lag.

The concrete process of preliminarily judging the truth of the causal relationship according to the positive and negative of the optimal time lag is as follows: the optimal time lag is less than or equal to 0, which indicates that the result variable can better estimate the past value but not the future value of the cause variable, and the cause and effect in the direction is a real cause and effect relationship; an optimal skew greater than 0 indicates that the effect precedes the cause, and the cause and effect of that direction is pseudo-cause and effect.

The specific process of the step (5) is as follows:

(5-1) for each target signal pair, if the optimal time lags of the target sub-signals are all larger than 0, judging that no causal relationship exists between the target signal pairs, and judging the next target signal pair; if both sides in the optimal time lag are less than or equal to 0, (5-2) is executed, and if only one side is less than or equal to 0, (5-3) is executed;

(5-2) dividing the sample length N of the target signal pair into a plurality of intervals of Δ N N ∈ [ N_min:ΔN:N_max]Calculating the length n according to steps (4-3) to (4-8) respectively

And Y (t) and

and X (t) are cross-mapped to an index ρ_X→Y、ρ_Y→XUntil N is equal to N_maxObtaining rho_X→YCurve sum p_Y→XA curve, and then executing the step (5-4);

(5-3) setting the time lag with the optimal time lag less than or equal to 0 as lambda^*The sample length N of a signal pair is divided into a number of intervals of Δ N N ∈ [ N_min:ΔN:N_max]Calculating X (t + lambda) with the length of n according to the steps (4-3) to (4-8) respectively^*) And Y (t) and Y (t + lambda)^*) And X (t) are cross-mapped to an index ρ_X→Y、ρ_Y→XUntil N is equal to N_maxObtaining rho_X→YCurve sum p_Y→XSelecting a rho curve in a direction corresponding to the optimal time lag of less than or equal to 0;

(5-4) calculation of ρ by the Z-transform method_X→Y、ρ_Y→XIf | Z | > 1.96 and the ρ curve increases with increasing length n until it no longer fluctuates, the ρ curve for that direction is considered to converge and the causal relationship for that direction is true, otherwise it does not.

Compared with the prior art, the invention has the following beneficial effects:

1. the invention only utilizes the oscillation circuit signal in the industrial control system, does not need external additional signal excitation when acquiring the signal, does not introduce additional disturbance to the control system, and can realize non-invasive detection and diagnosis.

2. The invention completely adopts a data driving type method, does not need prior process knowledge and does not need manual intervention.

3. The invention has no stationarity requirement on the acquired oscillation loop signal.

4. The algorithm principle of the invention is simple, easy to realize and convenient to operate.

5. The method performs feature selection through mutual information calculation based on the K-order nearest neighborhood, greatly reduces the workload of CCM causal analysis, and improves the analysis efficiency.

6. The invention carries out denoising (optional) and periodization item removing treatment on the original oscillation signal subjected to characteristic selection, improves the accuracy and reliability of CCM causal analysis results, and widens the application range of the CCM causal analysis results.

7. According to the CCM causal analysis method, the preliminary causal relationship is obtained according to the optimal time lag corresponding to the cross mapping index of the signal, and whether the cross mapping curve is converged is judged based on the optimal time lag, so that the accuracy and the reliability of the CCM causal analysis result are improved.

8. The method is suitable for causal analysis of the multi-loop oscillation signals of the industrial control system, determines an oscillation propagation path and accurately positions oscillation sources including linear and nonlinear oscillation sources according to a causal relationship network, and provides a new thought and a reliable basis for diagnosis of linear and nonlinear oscillation of the multi-loop oscillation signals of the industrial control system by using a CCM causal analysis method.

9. The method is also suitable for causal analysis of the global oscillation signal of the industrial control system, the oscillation propagation path is determined according to the causal relationship network, the oscillation source is accurately positioned and comprises a linear oscillation source and a nonlinear oscillation source, and a new thought and a reliable basis are provided for linear and nonlinear oscillation diagnosis of the global oscillation signal of the industrial control system by using a CCM causal analysis method.

Drawings

Fig. 1 is a schematic diagram of a TE process flow and a control strategy based on a distributed control strategy according to an embodiment of the present invention;

FIG. 2 is a schematic flow chart of a multi-loop oscillation causal relationship analysis method of an industrial control system based on improved CCM according to the present invention;

FIG. 3 is a diagram of process output signals of 17 to-be-analyzed industrial control system oscillation loops collected in an embodiment of the present invention;

FIG. 4 is a diagram of the results of 17 loop mutual information calculations in an embodiment of the present invention;

FIG. 5 is a diagram illustrating the EMD decomposition results of the 16 th loop in the embodiment of the present invention;

FIG. 6 is a graph showing the results of the optimal embedding dimension E of the 11 th, 16 th and 17 th loops in the embodiment of the present invention;

FIG. 7 is a diagram illustrating skew CCM detection results of the 16 th loop and the 17 th loop in the embodiment of the present invention;

FIG. 8 is a diagram illustrating the convergence of the optimal time-skewed CCM for the 16 th loop and the 17 th loop in an embodiment of the present invention;

FIG. 9 is a final causal relationship network in an embodiment of the present invention.

Detailed Description

The invention will be described in further detail below with reference to the drawings and examples, which are intended to facilitate the understanding of the invention without limiting it in any way.

The following describes in detail an improved CCM causal relationship analysis method for a chemical process with multiple oscillation loops, taking the oscillation loops in the TE process based on a distributed control strategy as an example.

The TE Process in this embodiment adopts a distributed Control strategy, and the flow and Control strategy are shown in FIG. 1 according to "Ricker N L. decentralized Control of the Tennessee Eastman cell Process [ J ]. Journal of Process Control,1996,6(4): 205-.

The TE process is a chemical process model designed by Eastman chemical company based on actual industrial processes. The whole process mainly comprises five operation units: the system comprises a reactor, a product condenser, a gas-liquid separator, a circulating compressor and a product desorption tower. There are four feed streams, one product stream and one waste gas stream. Two products G, H were produced from four reactants A, C, D, E, one inert component B and byproduct F, as follows:

a (gas) + C (gas) + D (gas) → G (liquid)

A (gas) + C (gas) + E (gas) → H (liquid)

A (gas) + E (gas) → F (liquid)

3D (gas) → 2F (liquid)

All reactions are exothermic and irreversible. The reaction rate is mainly influenced by the reaction temperature and the concentration of the reactants.

In fig. 1, RCL1, RCL2, RCL3, RCL4, RCL5, RCL6, RCL7, FCL8, LCL9, LCL10, LCL11, PCL12, XCL13, XCL14, XCL15, TCL16, and TCL17 are controlled variables of the 17 main circuits to be analyzed in the present embodiment. In fig. 1, V1, V2, V3, V4, V6, V7, V8, V10, and V11 are operation variables (actuator valves) in the respective circuits. Table 1 is the information about 17 major loops of the TE process based on the distributed control strategy.

TABLE 1

Loop numbering	Controlled variable	Operating variables
			1	A flow rate	V3
2	Flow rate of D	V1
			3	E flow rate	V2
4	Flow rate of C	V4
			5	Discharge rate of circulating compressor	V6
6	Liquid flow of gas-liquid separator	V7
			7	Liquid flow of product stripper	V8
8	Flow rate of product	Fp
			9	Liquid level of product resolution tower	Ratio of Loop 7 (V8)
10	Liquid level of gas-liquid separator	Loop				6 ratio (V7)
			11	Reactor level	Arrangement of circuit 17Fixed value (V11)
12	Reactor pressure	Ratio of Loop 5 (V6)
			13	Mod% G in product stream 11	Eadj(V1,V12)
14	Ratio of A to A + C	Loop				1 ratio, r1(V3)
			15	A+C	r1+r4(V3,V4)
16	Reactor temperature	V10
			17	Temperature of gas-liquid separator	V11

TE process equipment and reaction process parameter settings are implemented according to "downloads J, Vogel E F.A plant-wireless independent process control project [ J ]. Computers & chemical engineering,1993,17(3): 245-.

The TE process has a sample length of 721 for 17 loops (referred to as the loop's controlled variable) and a sampling interval of 0.1 hour (6 minutes). After the system operates stably for 20 hours, in a 16 th loop (L for short)₁₆) Adding oscillation disturbance (oscillation period is 105 sampling points), and global oscillation occurs in 17 loops.

As shown in fig. 2, an industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM includes:

step 1, collecting output signals of the oscillation loop process of each industrial control system to be analyzed.

The method for acquiring the process output signal comprises the following steps: the process data in the control loop to be detected is recorded in each preset sampling interval, and the process data collected in each sampling interval is added to the tail end of the process data collected previously.

The sampling interval refers to the sampling interval of the performance evaluation system. The process data is continuously updated over time, with new process data added to the end of the previously acquired process data for each length of time that a sampling interval has elapsed. The sampling interval of the performance evaluation system is generally the same as the control period in the industrial control system, and can also be selected as an integral multiple of the control period, and is specifically determined according to the real-time requirements and data storage capacity limitations of performance monitoring and industrial sites.

The raw data of the 17 process output signals of the oscillation Loop collected in this embodiment is shown in fig. 3, the abscissa is a sampling point ordinal number, and the unit is Sample (1 Sample corresponds to one sampling time, and the time corresponding to adjacent samples is a sampling interval), Loop is a Loop number, and the ordinate is a measurement value of a controlled variable of the Loop, and the unit is a unit of the controlled variable.

In this embodiment, in order to simulate the actual industrial process more ideally by adding the oscillation disturbance manually, the 250 th to 720 th sampling points of the output signal of the oscillation loop process are used as the loop data to be analyzed, that is, the data does not include the signal of the oscillation starting time.

And 2, calculating mutual information between output signals of the two oscillating circuits in the process, and obtaining a target signal pair through feature selection.

Mutual information is the reduction of uncertainty of another variable given the information of a variable, and can be used as a measure of the degree of correlation between the two variables. Theoretically, mutual information is positive, which indicates that correlation exists among variables, and the larger the value is, the higher the correlation degree among the variables is; if the mutual information is 0, no correlation exists between the variables. And (3) calculating mutual information values among the 17 oscillating circuits by using a K-order nearest neighbor method to select characteristics, removing circuits which are irrelevant or have low relevance degree, and screening out circuits with relatively high relevance degree (MI is larger than 0). A 95% confidence test was performed using the normalized substitution data method.

In this example, the K-nearest neighbor method is performed according to "Schreiber T.measuring information transfer [ J ]. Physical review letters,2000,85(2):461 ]. In this embodiment, the K-order neighborhood of the main parameter in the K-order nearest neighborhood method is set to 8.

In this embodiment, the normalized substitution data method is implemented according to "Sprite light string.

The mutual information calculation results of the 17 loops are shown in fig. 4, where the abscissa and ordinate of fig. 4 represent each loop, and the color patch represents the mutual information value. After mutual information characteristic selection, 40 pairs of loops (target signal pairs) enter the next step, and the number of the loops is 29% of the original analysis amount.

And 3, performing denoising and cycle item removing treatment on the signals in each target signal pair, and reconstructing the rest part to obtain corresponding target sub-signals.

The oscillation signals in the embodiment are periodic sinusoidal signals, phase differences often exist between the sinusoidal signals, and specific lag information between the sinusoidal signals and the sinusoidal signals cannot be determined directly according to the phase differences, so that the accuracy of causal analysis is reduced. The periodic signals may also cause the signals to become similar and coordinated, interfering with causal analysis. Besides periodic signals, unwanted noise, which is almost unavoidable in industrial processes, may also have an impact on the causal analysis.

The target signal is decomposed using the EMD method to obtain a series of modal components Imf and residual components. And calculating the DFA index of each layer of component by using a DFA method. If the DFA index of a certain layer is less than or equal to 0.5, the layer is judged to be useless noise, and if multiple layers exist, the layers are accumulated to be final noise. The period terms in Imf are determined by observation and all period terms are accumulated as the final period term. And removing noise (optional) and periodic terms in the target signal, and reconstructing the rest part to obtain a target sub-signal.

In this example, the EMD method is performed according to "Huang N E, Shen Z, Long S R, et al, the empirical mode decomposition and the Hilbert spectrum for non-linear and non-stationary time series analysis [ J ]. Proceedings of the Royal Society of London series A. physical, physical and engineering series, 1998,454(1971):903 and 995 ].

In this example, the DFA method was performed according to "Peng C K, Havlin S, Stanley H E, et.Quantification of scaling exponents and cross phenyl genomic extracted time series [ J ]. Chaos: an inter scientific family of nonlinear science,1995,5(1):82-87 ].

With L in this embodiment₁₆For example, each component obtained by performing EMD decomposition is shown in fig. 5, and the DFA index corresponding to each component is shown in table 2. In fig. 5, the abscissa is the ordinal number of a sampling point, the unit is Sample (1 Sample corresponds to a sampling time, and the corresponding time between adjacent samples is a sampling interval), and the first layer Original Data is L₁₆The second layer Imf1 through the seventh layer Imf6 are the 6 modal components in the decomposition, the last layer Residual is the Residual component, and the ordinate units are all in degrees c. As is clear from table 2, the DFA index of Imf1 is less than 0.5, and therefore, it is judged to be noise. As can be seen from FIG. 5, Imf4, Imf5 and Imf6 are periodic terms, and the period is around 105, which corresponds to the period of the applied oscillation disturbance in the present embodiment. Thus for L₁₆Removing noise Imf1, removing period term Imf4+ Imf5+ Imf6, and reconstructing the rest part to obtain the target sub-signal. The target sub-signals of the other 16 loops can be obtained in the same way.

TABLE 2

Modal component	DFA index
		Imf1	0.2985
Imf2	0.6427
		Imf3	1.5554
Imf4	3.4493
		Imf5	3.8271
Imf6	3.9703
		Imf7	4.0525

And 4, determining the optimal embedding dimension of each target sub-signal, calculating the cross mapping indexes of each target sub-signal pair under different time lags by using a CCM (continuous current cycle) method, and preliminarily judging the truth of the causal relationship according to the positive and negative of the optimal time lags.

And 4-1, determining the optimal embedding dimension E of each target sub-signal by using a G-cure method. With L₁₁、L₁₆、L₁₇For example, the result of determining the optimal embedding dimension E by the G-cure method is shown in FIG. 6. In FIG. 6, the abscissa Embedding dimension is the Embedding dimension E, the ordinate G (k) is the value of the G function, E_L11、E_L16、E_L17Are respectively L₁₁、L₁₆、L₁₇The embedding dimension of (a). The determination criteria for the optimal embedding dimension E are as follows: and when the G function is less than 0.5, the corresponding E is the optimal embedding dimension. In this embodiment, the optimizationThe embedding dimension E is 5.

Step 4-2, determining two variables X and Y of each target sub-signal pair, and setting the time lag lambda ∈ [ -c: deltac: c) of X relative to Y as a preset cause variable and a result variable]. For each time lag, ρ between X (t + λ) and Y (t) can be calculated in steps 4-3 to 4-8_X→Y(X (t) is replaced by X (t + lambda); similarly, rho between Y (t + lambda) and X (t)_Y→X(X (t) is replaced by Y (t + lambda) and Y (t) is replaced by X (t)). In this embodiment, c is 50.

Step 4-3, for a pair of target sub-signal pairs X (t) and Y (t) with the sample length of P, respectively establishing corresponding time lag coordinate vectors, wherein the calculation formula is as follows:

wherein E is the optimal embedding dimension determined in the step 4-1, if the E of the target sub-signal pair is different, the maximum value of the two is selected as the final E, tau is the time lag and is generally regarded as 1 by default,x(t)、ythe set of (t) is the corresponding shadow manifold M of X (t), Y (t)_X、M_Y. In this embodiment, E is 5, τ is 1, and P is 370.

Step 4-4, finding the shadow manifold M_YPoint ofy(t) obtaining a distance sequence { d) in descending order from small to large from (E +1) nearest neighbors at time t₁,d₂,...,d_E+1And a sequence of time indices t corresponding to the sequence of distances₁,t₂,...,t_E+1And (c) the step of (c) in which,xthe distance between (t) is calculated as follows:

d_i＝D(y(t_i),y(t_j))

where D (a, b) refers to the Euclidean distance between vectors a and b.

Step 4-5, with each pointy(E +1) nearest neighbors of (t) calculating relevant weights, and calculating the weightsThe formula is as follows:

w_i＝u_i/U

in the formula

Normalization factor

Step 4-6, calculating the estimated value of X by weight

The calculation formula is as follows:

in the formula w_iIs a weight, t_iIs a time index corresponding to the weight.

Step 4-7, calculating X and

standard Pearson's correlation coefficient between

(ρ_X→Y) As a cross-map correlation coefficient, the causal strength of X → Y is represented by the following calculation formula:

4-8, utilizing a substitution data method to correct rho_X→YCarrying out 95% confidence test;

step 4-9, respectively obtaining a series of cross mapping indexes changing along with time lag in the X → Y direction and the Y → X direction from the steps 4-2 to 4-8, wherein the time lag corresponding to the maximum rho value is the optimal time lag lambda_X→Y、λ*_Y→XThe truth of the causal relationship is preliminarily judged according to the positive and negative of the optimal time lag (judged by a computer), and the concrete mode is as follows: optimal skew is less than or equal to 0, tableThe cause variable can better estimate the past value of the cause variable but not the future value, and the cause and effect of the direction is a real cause and effect relationship; an optimal skew greater than 0 indicates that the effect precedes the cause, which is clearly contrary to the definition of causality, so causality in this direction is pseudo-causality.

In this embodiment, the G-cure method is implemented according to "Wang Y, Hu F, Cao Y, et al, improved CCM facilitated availability detection in complex systems [ J ]. Control engineering practice,2019,83:67-82.

In this embodiment, the substitution data method is implemented according to "Sprite, substitution data and application [ J ]. university of east China, 2011:1-2 ]. But the normalization step in which the mean and euclidean distances are calculated is left out. Since the CCM method involves convergence, the causal strength related indicator theoretically takes whether 1 is reached as a basis for convergence, and the indicator ρ itself is a normalization coefficient, the detection result is expressed by an absolute value rather than a relative value.

With L₁₆As X, L₁₇As an example of Y, fig. 7 shows the detection result of the skew CCM. In fig. 7, the abscissa time lag λ is the time lag between variables, and the ordinate ρ is the cross-mapping index between corresponding variables. The red curve represents the X → Y direction, and the blue curve represents the Y → X direction. As can be seen from fig. 7, in the present embodiment,

then L is₁₆→L₁₇For real cause-and-effect, L₁₇→L₁₆Is a pseudo-causal relationship. Thus preliminarily judging L₁₆Is L₁₇Reason for (L)₁₆→L₁₇) I.e. L₁₆Influence of one way L₁₇. The same results can be obtained for the preliminary causal relationship between the other 39 pairs of target sub-signals.

And 5, calculating cross mapping indexes of each target sub-signal pair corresponding to the optimal time lag under different sample lengths, and judging whether convergence occurs according to a convergence threshold value to obtain the causal relationship of each signal pair.

Step 5-1, for each target signal pair, if the optimal time lags of the target sub-signals are all larger than 0, judging that no causal relationship exists between the target signal pairs, and judging the next target signal pair; and if both sides in the optimal time lag are less than or equal to 0, directly executing the step 5-2, and if only one side is less than or equal to 0, directly executing the step 5-3.

Step 5-2, dividing the sample length N of the signal pair into a plurality of intervals with the interval delta N, N ∈ [ N_min:ΔN:N_max]Calculating the length n according to steps 4-3 to 4-8 respectively

And Y (t) and

and X (t) are cross-mapped to an index ρ_X→Y、ρ_Y→XUntil N is equal to N_maxObtaining rho_X→YCurve sum p_Y→XCurve and then perform step 5-4. In this embodiment, N is 471, N_min＝10，ΔN＝10，Nmax＝450。

Step 5-3, setting the time lag with the optimal time lag less than or equal to 0 as lambda^*The sample length N of a signal pair is divided into a number of intervals of Δ N N ∈ [ N_min:ΔN:N_max]Calculating X (t + lambda) of length n according to steps 4-3 to 4-8, respectively^*) And Y (t) and Y (t + lambda)^*) And X (t) are cross-mapped to an index ρ_X→Y、ρ_Y→XUntil N is equal to N_maxObtaining rho_X→YCurve sum p_Y→XAnd (4) a curve, namely a rho curve with the optimal time lag being less than or equal to 0 corresponding direction is focused. In this embodiment, N is 471, N_min＝10，ΔN＝10，Nmax＝450。

Step 5-4, calculating rho by using a Z conversion method_X→Y、ρ_Y→XIf Z > 1.96 and the ρ curve increases with increasing length n until it no longer fluctuates, the ρ curve for that direction is considered to converge and the causal relationship for that direction holds, otherwise it does not hold.

In this embodiment, the Z transformation method is implemented according to "Aftab M F, Hovd M, Sivalingam S.Convergentcross mapping (CCM)" based approach for approximating the source of plant-wide disorders [ C ]//2017 IEEE Conference on Control Technology and Applications (CCTA). IEEE,2017: 1492-.

With L₁₆As X, L₁₇As an example of Y, obtained by step 4

Then the cross mapping indexes between X (t-2) and Y (t-2) and X (t) at different sample lengths are calculated. Fig. 8 is a diagram of optimal skew CCM convergence. In fig. 8, the abscissa L is the sample length and the ordinate ρ is the cross-mapping index between the corresponding variables. In the present embodiment, the first and second electrodes are,

the curve converges and the corresponding Z index is 9.814, greater than 1.96, so L₁₆→L₁₇Is established in a causal relationship of L₁₆Is L₁₇The reason for (1). The final causal relationship judgment result between other 39 pairs of target sub-signals can be obtained in the same way.

And 6, determining an oscillation propagation path and positioning an oscillation source according to the final causal relationship network.

When the variables are causal to each other, the variable that first gives an influence can be considered a causal variable. In this embodiment, the final causal relationship network is shown in fig. 9. The rows in FIG. 9 are result variables and the columns are cause variables. For one-way causality, a black square (value 1) represents that this direction causality exists; for bi-directional causality, the black box (value 1) represents the direction causality before the reverse direction causality, and the vermilion box (value 0.5) represents the direction causality after the reverse direction causality. If both directions are black squares, it is indicated that both affect each other almost simultaneously. If the coordinates of the 17 th row and the 16 th column are (17,16), the corresponding color blocks are black squares, which represents L₁₆For the reason, L₁₇As a result, the direction cause and effect relationship exists.

Loop oscillation cause and effect in this exampleThe results of the relationship analysis show that L₁₆By direct or indirect influence of L₂、L₃And L₁₃Other variables than L, and only L₁₂Influence (known from the theoretical knowledge, L)₁₆First acts on L₁₂Then L is₁₂Immediately counteracts L₁₆) Thus determining L₁₆The source of the TE process oscillation coincides well with the location in this embodiment where the oscillation disturbance is manually added. And the whole oscillation propagation path conforms to the process flow information flow and the control strategy of the TE process.

By utilizing the method, the multi-loop oscillation signals in the industrial control system can be quickly and accurately diagnosed, the oscillation propagation path can be quickly and accurately diagnosed, the oscillation source can be positioned, a new thought is provided for carrying out fault diagnosis on the oscillation loop of the industrial control system by using a causal analysis method, and the method has certain reference significance and application value.

The embodiments described above are intended to illustrate the technical solutions and advantages of the present invention, and it should be understood that the above-mentioned embodiments are only specific embodiments of the present invention, and are not intended to limit the present invention, and any modifications, additions and equivalents made within the scope of the principles of the present invention should be included in the scope of the present invention.

Claims (6)

1. An industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM is characterized by comprising the following steps:

(1) collecting process output signals of oscillation loops of the industrial control systems to be analyzed;

(2) calculating mutual information between output signals of two oscillating circuits in the process, and obtaining a target signal pair after feature selection;

(3) denoising and cycle item removing processing are carried out on signals in each target signal pair, and the rest part is reconstructed to obtain corresponding target sub-signals;

(4) determining the optimal embedding dimension of each target sub-signal, calculating cross mapping indexes of each target sub-signal pair under different time lags by using a CCM (continuous current measurement) method, and preliminarily judging the truth of a causal relationship according to the positive and negative of the optimal time lags;

(5) calculating cross mapping indexes of each target sub-signal pair corresponding to the optimal time lag under different sample lengths, and judging whether convergence occurs according to a convergence threshold value to obtain a causal relationship of each signal pair;

(6) and determining an oscillation propagation path and positioning an oscillation source according to the final causal relationship network.

2. The improved CCM-based industrial control system multi-loop oscillation causal relationship analysis method of claim 1, wherein the specific process of step (2) is as follows:

(2-1) calculating a mutual information value between process output signals of every two oscillating circuits by using a K-order nearest neighbor method;

(2-2) carrying out 95% confidence test on the mutual information value by using a normalized substitution data method;

(2-3) retaining the signal pair having the mutual information value greater than 0 as the target signal pair.

3. The improved CCM-based industrial control system multi-loop oscillation causal relationship analysis method of claim 1, wherein the specific process of step (3) is as follows:

(3-1) decomposing the target signal by using an EMD method to obtain a series of modal components and residual components;

(3-2) calculating the DFA indexes of each layer of components by using a DFA method, if all the DFA indexes are more than 0.5, judging that the DFA indexes do not contain noise, and directly executing the step (3-3); if the DFA index of a certain layer is less than or equal to 0.5, judging the layer to be useless noise, and if multiple layers exist, accumulating the layers to be final noise;

(3-3) judging the period items in the modal components, and accumulating all the period items to be used as final period items;

and (3-4) if the signal contains noise, performing denoising and cycle item removing processing on the target signal, if the signal does not contain noise, directly performing cycle item removing processing on the target signal, and reconstructing the rest part after the processing to obtain a target sub-signal.

4. The improved CCM-based industrial control system multi-loop oscillation causal relationship analysis method of claim 1, wherein the specific process of step (4) is as follows:

(4-1) determining the optimal embedding dimension E of each target sub-signal by using a G-curve method, wherein the determination standard is as follows: when the G function is less than 0.5, the corresponding E is the optimal embedding dimension;

(4-2) determining two variables X and Y of each target sub-signal pair, and setting the time lag lambda epsilon [ -c: deltac: c ] of X relative to Y as a preset cause variable and a preset result variable;

(4-3) for the target sub-signal pair X (t) and Y (t) with the sample length P, respectively establishing corresponding time lag coordinate vectors, wherein the calculation formula is as follows:

in the formula, E is the optimal embedding dimension determined in the step (4-1), and if the E of the target sub-signal pair is different, the maximum value of the E and the maximum value is selected as the final E; τ is the time lag, and is 1 by default; the set of x (t), y (t) is X (t), Y (t) corresponding to the shadow manifold M_X、M_Y；

(4-4) finding the shadow manifold M_YPoint ofy(t) obtaining a distance sequence { d) in descending order from the E +1 nearest neighbors at the time t₁,d₂,...,d_E+1And a sequence of time indices t corresponding to the sequence of distances₁,t₂,…,t_E+1And (c) the step of (c) in which,xthe distance between (t) is calculated as follows:

d_i＝D(y(t_i),y(t_j))

wherein D (a, b) represents the euclidean distance between vectors a and b;

(4-5) use of each doty(t) calculating the relevant weights of the E +1 nearest neighbors, wherein the calculation formula of the weights is as follows:

w_i＝u_i/U

in the formula (I), the compound is shown in the specification,

normalization factor

(4-6) calculating an estimated value of X by using the weight

The calculation formula is as follows:

in the formula, w_iIs a weight, t_iIs a time index corresponding to the weight;

(4-7) calculating X and

standard pearson correlation coefficient p between_X→YAs a cross-map correlation coefficient, the causal strength of X → Y is represented by the following calculation formula:

(4-8) pairing ρ with a substitution data method_X→YCarrying out 95% confidence test;

(4-9) obtaining a series of cross mapping indexes of the X → Y direction changing along with time lag by adopting the steps, and obtaining a series of cross mapping indexes of the Y → X direction changing along with time lag by using the same method as the steps (4-2) to (4-8); wherein the time lag corresponding to the maximum rho value is the optimal time lag lambda_X→Y、λ*_Y→XAnd preliminarily judging the truth of the causal relationship according to the positive and negative of the optimal time lag.

5. The improved CCM-based industrial control system multi-loop oscillation causal relationship analysis method of claim 4, wherein in the step (4-9), the detailed process of preliminarily judging the authenticity of the causal relationship according to the positive and negative of the optimal time lag is as follows:

the optimal time lag is less than or equal to 0, which indicates that the result variable can better estimate the past value but not the future value of the cause variable, and the cause and effect in the direction is a real cause and effect relationship; an optimal skew greater than 0 indicates that the effect precedes the cause, and the cause and effect of that direction is pseudo-cause and effect.

6. The improved CCM-based industrial control system multi-loop oscillation causal relationship analysis method of claim 5, wherein the specific process of step (5) is as follows:

(5-1) for each target signal pair, if the optimal time lags of the target sub-signals are all larger than 0, judging that no causal relationship exists between the target signal pairs, and judging the next target signal pair; if both sides in the optimal time lag are less than or equal to 0, (5-2) is executed, and if only one side is less than or equal to 0, (5-3) is executed;

(5-2) dividing the sample length N of the target signal pair into a plurality of intervals of Δ N N ∈ [ N_min:ΔN:N_max]Calculating the length n according to steps (4-3) to (4-8) respectively

And Y (t) and

and X (t) are cross-mapped to an index ρ_X→Y、ρ_Y→XUntil N is equal to N_maxObtaining rho_X→YCurve sum p_Y→XA curve, and then executing the step (5-4);

(5-3) setting the time lag with the optimal time lag less than or equal to 0 as lambda^*The sample length N of a signal pair is divided into a number of intervals of Δ N N ∈ [ N_min:ΔN:N_max]Calculating X (t + lambda) with the length of n according to the steps (4-3) to (4-8) respectively^*) And Y (t) andand Y (t + lambda)^*) And X (t) are cross-mapped to an index ρ_X→Y、ρ_Y→XUntil N is equal to N_maxObtaining rho_X→YCurve sum p_Y→XSelecting a rho curve in a direction corresponding to the optimal time lag of less than or equal to 0;

(5-4) calculation of ρ by the Z-transform method_X→Y、ρ_Y→XIf | Z | > 1.96 and the ρ curve increases with increasing length n until it no longer fluctuates, the ρ curve for that direction is considered to converge and the causal relationship for that direction is true, otherwise it does not.

Images