CN111626099B - 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: and calculating mutual information between process output signals of every two oscillation loops to be analyzed, carrying out feature selection, reserving loop signal pairs with relatively high correlation degree, removing noise and period items by using an EMD (empirical mode decomposition) method and a DFA (complementary mode decomposition) method, calculating cross mapping indexes of reconstructed sub-signals under different time lags by using a CCM (continuous mode correlation) method, primarily judging authenticity of causal relations according to positive and negative of optimal time lags, calculating cross mapping indexes of each target sub-signal pair corresponding to the optimal time lags under different sample lengths, judging whether convergence is carried out according to a convergence threshold value, and finally obtaining a final causal relation network of the oscillation loops, thereby determining an oscillation propagation path and positioning an oscillation source. The invention can realize rapid and accurate diagnosis of the oscillation propagation path in the industrial control system and position the oscillation source, and provides a new thought for fault diagnosis of the industrial control system oscillation loop 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

About 30% of industrial control loops have oscillating faults due to controller overscaling, valve sticking, external disturbances, and process nonlinearities. Oscillation is one of the main causes of degradation in control performance of industrial control systems. The continuous oscillation increases the abrasion of industrial equipment, reduces the quality of products, seriously influences the production benefit of factories and even generates potential safety hazards. There may be thousands of loops in a large industrial control system, and due to the interaction and high coupling between loops, the oscillations of any one loop may be transferred to the other loops of the system, ultimately resulting in global oscillations.

And carrying out causal analysis on an oscillating circuit in the industrial process, determining an oscillating propagation path and positioning an oscillating source so as to take corresponding measures more specifically, maintain a damaged machine in time, shorten the downtime and reduce the related loss. Therefore, the method is designed for the effective causal analysis method of the oscillating circuit, efficiently and accurately diagnoses the oscillating propagation path and accurately locates the fault source, and has important significance and reference value for fault diagnosis of the industrial control circuit.

In recent years, industrial process fault diagnosis based on causal analysis has progressed rapidly. Wakefield B J et al used Granger and transfer entropy methods to perform fault diagnosis on simulated mineral processing plant milling circuits. Bauer M et al estimate the time delay between process variables using a cross-correlation function, 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 a concept of direct transfer entropy (Direct Transfer Entropy, DTE) for detecting whether a direct information flow path exists between variables and a normalization method of differential direct transfer entropy to determine direct causal relationship strength. However, the Granger method has larger model fitting error for nonlinear variables, unstable variables and variables in an inseparable system with complex relation, and can present causal relation detection results of missed judgment and misjudgment. The transfer entropy method involves estimation of probability density functions and is sensitive to parameters and low in calculation efficiency.

The Convergent Cross Map (CCM) method proposed by Sugihara can be used for an inseparable dynamic system with weak to medium coupling strength, overcomes the defect that the Granger method cannot be used for the inseparable system, and has the advantages of simple principle, easy realization and higher calculation efficiency. However, CCM methods can only obtain their causal relationship by pairwise analysis of variables. If a large industrial process generates global oscillation, for the purpose of carrying out oscillation diagnosis on the large industrial process, one method is to carry out causal analysis on all oscillation loops two by two to obtain a related causal network. However, this method has significant limitations: 1) The method has the advantages of excessive related variables, extremely large calculated amount, complex process and extremely low efficiency; 2) May involve unnecessary causal relationship detection, which is time consuming and labor intensive; 3) Because of the large number of variables and large span, the analysis result is not necessarily accurate. In addition, the method is proposed based on an ecosystem and is mostly used for ecology and neuroscience, and industrial processes mostly have high-dimensional, strong-coupling and other characteristics, which are remarkably different from the ecosystem applicable to the traditional CCM method. Especially, part of the oscillating signals are periodic-like sinusoidal signals, phase differences exist between the sinusoidal signals, and specific hysteresis information between the sinusoidal signals cannot be determined directly according to the phase differences, so that the accuracy of causal analysis is reduced. Periodic signals may also cause the signals to become similar, synergistic, and interfere with causal analysis. In addition to periodic signals, almost unavoidable unwanted noise in industrial processes can also have an impact on causal analysis.

Based on the background, aiming at the multi-loop oscillation signals in the industrial control system, the CCM method is improved, so that unnecessary detection workload is greatly reduced, causal analysis efficiency is improved, meanwhile, interference items in the signals are eliminated, and therefore accuracy and reliability of causal analysis results are greatly improved, the rapid and accurate diagnosis of oscillation propagation paths and the positioning of oscillation sources are realized, 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 an improved CCM, which is suitable for diagnosing an oscillation propagation path and positioning an oscillation source in a multi-loop oscillation signal, has high analysis efficiency and high accuracy, and only needs to acquire industrial process related oscillation data without process mechanism knowledge.

An industrial control system multi-loop oscillation causal relationship analysis method based on improved CCM, comprising:

(1) Collecting process output signals of each industrial control system oscillation circuit to be analyzed;

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

(3) Denoising and removing period item 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 the cross mapping index of each target sub-signal pair under different time lags by using a CCM method, and preliminarily judging the authenticity of the 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, judging whether convergence is carried out according to a convergence threshold value, and obtaining causal relation of each signal pair;

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

The invention performs characteristic selection by calculating mutual information between output signals of the process of the two-by-two oscillation loops to be analyzed, reserves loop signal pairs with relatively higher correlation degree, removes noise (optional) and period items by using an EMD (empirical mode decomposition) method and a DFA (differential mode decomposition) method, calculates cross mapping indexes of reconstructed sub-signals under different time lags by using a CCM (continuous mode decomposition) method, preliminarily judges the authenticity of causality according to the positive and negative of the optimal time lag, calculates cross mapping indexes of each target sub-signal pair corresponding to the optimal time lag under different sample lengths, judges whether convergence is carried out according to a convergence threshold, and finally obtains a final causality network of the oscillation loops, 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, simultaneously eliminates interference items in the signals, greatly improves the accuracy and reliability of causal analysis results, can rapidly and accurately diagnose an oscillation propagation path and position an oscillation source, and provides a new thought for fault diagnosis of the industrial control system oscillation loop by using a causal analysis method, thereby having certain reference significance and application value.

The invention directly adopts the measurable quantity of the chemical process as the process output signal, and all the process output signals of the oscillating circuit 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 the process output signals of every two oscillating circuits by using a K-order nearest neighbor method;

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

(2-3) reserving a signal pair having a mutual information value greater than 0 as a 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 index of each layer component by using the DFA method, and if all the DFA indexes are greater than 0.5, judging that the component does not contain noise, and directly performing the step (3-3); if the DFA index of a certain layer is less than or equal to 0.5, judging the layer as useless noise, and if a plurality of layers exist, accumulating the layers as final noise;

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

and (3-4) if the signal contains noise, carrying out denoising and cycle term removal on the target signal, if the signal does not contain noise, directly carrying out cycle term removal on the target signal, and reconstructing the rest part after the processing is finished 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-cure method, wherein the determination standard is as follows: and when the G function is smaller 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, setting a time delay λε [ -c: Δc: c of X relative to Y as preset cause and result variables]. For each time lag, ρ between X (t+λ) and Y (t) can be calculated as (4-3) to (4-8) _X→Y (X (t+λ) is used instead of X (t)); similarly, ρ between Y (t+λ) and X (t) can be obtained _Y→X (Y (t+λ) is used instead of X (t), and X (t) is used instead of Y (t));

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

wherein E is the optimal embedding dimension determined in the step (4-1), and if E of the target sub-signal pair is different, the maximum value of the two is selected as the final E; τ is the time lag, defaulting to 1;x(t)、ythe set of (t) is shadow manifold M corresponding to X (t) and Y (t) _X 、M _Y ；

(4-4) finding shadow manifold M _Y Points ony(t) E+1 nearest neighbors at time t, obtaining a distance sequence { d ] arranged in order of small to large ₁ ,d ₂ ,...,d _E+1 Time index sequence { t } corresponding to a distance sequence ₁ ,t ₂ ,...,t _E+1 And } wherein,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) with each dotyE+1 nearest neighbors of (t) calculate the correlation weights, the weight calculation formula is as follows:

w _i ＝u _i /U

in the method, in the process of the invention,

normalization factor->

(4-6) calculating the estimated value of X with weights

The calculation formula is as follows:

wherein w is _i Is weight, t _i Is a time index corresponding to the weight;

(4-7) calculating X and

standard pearson correlation coefficient ρ therebetween _X→Y As a cross-map correlation coefficient, representing the causal intensity of x→y, the calculation formula is as follows:

(4-8) Using alternative data method vs ρ _X→Y Performing 95% confidence test;

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

The concrete process for 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, and the result variable can better estimate the past value of the reason variable instead of the future value, and the cause and effect of the direction is a real cause and effect relationship; an optimal time lag greater than 0 indicates that the result occurs before the cause, and the cause and effect of the direction is a pseudo-cause and effect relationship.

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 greater than 0, judging that there is no causal relationship between the target signal pairs, and judging the next target signal pair; if both sides in the optimal time lag are smaller than or equal to 0, executing (5-2), and if only one side is smaller than or equal to 0, executing (5-3);

(5-2) dividing the sample length N of the target signal pair into a plurality of sections with an interval Δn: n E [ N ] _min :ΔN:N _max ]Calculating the length n according to the steps (4-3) to (4-8)

And Y (t)/(t)>

Cross-map index ρ with X (t) _X→Y 、ρ _Y→X Up to n=n _max Obtain ρ _X→Y Curve sum ρ _Y→X Curve, then performing step (5-4);

(5-3) setting the time lag of 0 or less to lambda ^* The sample length N of a signal pair is divided into several intervals of interval Δn: n E [ 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) ^* ) Cross-map index ρ with X (t) _X→Y 、ρ _Y→X Up to n=n _max Obtain ρ _X→Y Curve sum ρ _Y→X A curve, wherein a rho curve with the optimal time lag smaller than or equal to 0 in the corresponding direction is selected;

(5-4) calculating ρ using Z-transform method _X→Y 、ρ _Y→X If |Z| > 1.96, and the ρ curve increases with increasing length n until no longer fluctuates, then the ρ curve in that direction is considered to converge, and the causal relationship in that direction is true, otherwise it is not.

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

1. the invention only uses the oscillation circuit signal in the industrial control system, does not need external additional signal excitation when collecting signals, 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 method, does not need prior knowledge of the process and does not need manual intervention.

3. The invention has no stability requirement on the collected oscillation loop signals.

4. The algorithm of the invention has simple principle, easy realization and convenient operation.

5. According to the CCM causal analysis method, the characteristic selection is performed through mutual information calculation based on the K-order nearest neighbor, so that the CCM causal analysis workload is greatly reduced, and the analysis efficiency is improved.

6. The invention carries out denoising (optional) cycle term removal treatment on the original oscillation signal subjected to feature 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 primary causal relation 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 invention is suitable for causal analysis of multi-loop oscillation signals of an industrial control system, determines an oscillation propagation path according to a causal relation network and accurately positions an oscillation source, comprises linear and nonlinear oscillation sources, and provides a new thought and a reliable basis for linear and nonlinear oscillation diagnosis of the multi-loop oscillation signals of the industrial control system by utilizing a CCM causal analysis method.

9. The invention is also suitable for causal analysis of global oscillation signals of an industrial control system, determines an oscillation propagation path according to a causal relation network and accurately positions oscillation sources, including linear and nonlinear oscillation sources, and provides a new thought and a reliable basis for diagnosing linear and nonlinear oscillations of the global oscillation signals of the industrial control system by using a CCM causal analysis method.

Drawings

FIG. 1 is a schematic diagram of TE process flow and control strategy based on distributed control strategy in 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 an improved CCM according to the present invention;

FIG. 3 is a graph of process output signals of 17 industrial control system tank circuits to be analyzed collected in an embodiment of the present invention;

FIG. 4 is a diagram showing the calculation results of 17 mutual information of loops in the embodiment of the invention;

FIG. 5 is a graph showing the EMD decomposition result of the 16 th loop in the embodiment of the present invention;

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

FIG. 7 is a graph showing the results of the time-lapse CCM detection of the 16 th loop and the 17 th loop in the embodiment of the present invention;

FIG. 8 is a graph showing the convergence of the optimal time-lapse CCM for the 16 th loop and the 17 th loop in the 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 with reference to the drawings and examples, it being noted that the examples described below are intended to facilitate the understanding of the invention and are not intended to limit the invention in any way.

Taking an oscillation circuit in a TE process based on a distributed control strategy as an example, an improved CCM causal relationship analysis method for a chemical process with a plurality of oscillation circuits is described in detail below.

The TE process in this embodiment employs a distributed control strategy, the flow and control strategy is shown in FIG. 1, according to "Ricker N L. Decentized control of the Tennessee Eastman challenge process [ J ]. Journal of Process Control,1996,6 (4): 205-221 ].

TE process is a chemical process model designed by Eastman chemical company based on actual industrial process. The whole flow mainly comprises five operation units: a reactor, a product condenser, a gas-liquid separator, a recycle compressor and a product analysis tower. There are four feed streams, one product stream and one waste stream. Two products G, H are produced from four reactants A, C, D, E, together with an inert component B and by-product 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 affected by the reaction temperature and the concentration of the reactants.

In fig. 1, 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 of the present embodiment. In fig. 1, V2, V3, V4, V6, V7, V8, V10, V11 are operation variables (execution valves) in each circuit, respectively. Table 1 is information about 17 main loops of TE process based on distributed control strategy.

TABLE 1

Loop numbering	Controlled variable	Operating variables
			1	A flow rate	V3
2	D flow rate	V1
			3	E flow rate	V2
4	C flow rate	V4
			5	Compressor discharge rate	V6
6	Liquid flow rate of gas-liquid separator	V7
			7	Liquid flow rate of product analysis tower	V8
8	Product flow rate	Fp
			9	Liquid level of product analysis tower	Ratio of loop 7 (V8)
10	Liquid level of gas-liquid separator	Ratio of Loop 6 (V7)
			11	Reactor liquid level	Setting value of loop 17 (V11)
12	Reactor pressure	Ratio of loop 5 (V6)
			13	Mod% G in product stream 11	E adj (V1,V12)
14	Ratio of A to A+C	Ratio of loop 1, r1 (V3)
			15	A+C	r1+r4(V3,V4)
16	Reactor temperature	V10
			17	Temperature of gas-liquid separator	V11

The setting of the parameters of the equipment and the reaction process of the TE process is carried out according to the methods of 'Down J J, vogel E F.A plant-wide industrial process control problem [ J ]. Computers & chemical engineering,1993,17 (3): 245-255'.

The TE process has a sample length 721 of 17 loops (referring to the controlled variables of the loops) with a sampling interval of 0.1 hours (6 minutes). After the system is stably operated for 20 hours, the system is started in a 16 th loop (L for short ₁₆ ) Adding oscillation disturbance (oscillation period is 105 sampling points), and generating global oscillation for 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 process output signals of each industrial control system oscillation circuit to be analyzed.

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

The sampling interval refers to the sampling interval of the performance evaluation system. The process data is updated over time, with new process data being added to the end of the previously acquired process data every time a sampling interval elapses. The sampling interval of the performance evaluation system is generally the same as the control period in the industrial control system, and can be selected to be an integral multiple of the control period, and is specifically determined according to the performance monitoring and the real-time requirements of the industrial field and the limitation of data storage.

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

In this embodiment, the oscillation disturbance is added manually, and in order to more ideally simulate the actual industrial process, 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, i.e. the data does not include the signal at the oscillation starting time.

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

The mutual information is the reduction of uncertainty of another variable given variable information, and can be used as a measure of the degree of correlation between two variables. Theoretically, the mutual information is positive, which means that there is a correlation between variables, and the greater the value, the higher the degree of correlation between variables; and if the mutual information is 0, no correlation exists between the variables. And calculating mutual information values among 17 oscillation loops by using a K-order nearest neighbor method, performing feature selection, removing loops with irrelevant or low correlation degree, and screening loops with relatively high correlation degree (MI is more than 0). A 95% confidence test was performed using the normalized surrogate data method.

In this example, the K-th order nearest neighbor method is implemented 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 neighbor method is set to 8.

In this example, the normalized substitution data method was performed according to "data on substitution data and application [ J ]. University of eastern China, 2011:1-2.

The calculation result of mutual information of 17 loops is shown in fig. 4, the abscissa and the ordinate of fig. 4 are the loops, and the color block is the mutual information value. After mutual information feature selection, 40 pairs of loops (target signal pairs) are added to the next step, wherein the total number of loops is 29% of the original analysis amount.

And step 3, carrying out denoising and cycle term removal processing on the signals in each target signal pair, and reconstructing the rest part to obtain corresponding target sub-signals.

In this embodiment, the oscillation signals are periodic sinusoidal signals, and there is often a phase difference between the sinusoidal signals, and specific hysteresis information between the sinusoidal signals cannot be determined directly according to the phase difference, so that accuracy of causal analysis is reduced. Periodic signals may also cause the signals to become similar, synergistic, and interfere with causal analysis. In addition to periodic signals, almost unavoidable unwanted noise in industrial processes can also have an impact on causal analysis.

The target signal is decomposed using an EMD method to obtain a series of modal components Imf and residual components. And then 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 a plurality of layers exist, the layers are accumulated to be the final noise. The period entries in Imf are determined by observation and all period entries are accumulated as the final period entry. Noise (optional) and period terms in the target signal are removed, and the remaining part is reconstructed to obtain a target sub-signal.

In this example, the EMD method is performed in accordance with "Huang N E, shen Z, long S R, et al empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis [ J ]. Proceedings of the Royal Society of London. Series A: matrical, physical and engineering sciences,1998,454 (1971):903-995.

In this example, the DFA method is carried out in accordance with "Peng C K, havlin S, stanley H E, et al Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series [ J ]. Chaos: an interdisciplinary journal of nonlinear science,1995,5 (1): 82-87 ].

In the present embodimentL of (3) ₁₆ For example, the components obtained by EMD-decomposing the components are shown in fig. 5, and the DFA index corresponding to each component is shown in table 2. In FIG. 5, the abscissa indicates the sampling point number in samples (1 Sample corresponds to one sampling time, and the corresponding time between adjacent samples is the sampling interval), and the first layer Original Data is L ₁₆ The second layer Imf to the seventh layer Imf6 are the 6 modal components of the decomposition result, the last layer Residual is the Residual component, and the ordinate units are all in degrees celsius. As is clear from table 2, since the DFA index of Imf1 is smaller than 0.5, it is determined as noise. As can be seen from fig. 5, imf, imf and Imf6 are periodic terms, and the period is about 105, which corresponds to the period of the oscillation disturbance applied in the present embodiment. Thus for L ₁₆ Noise Imf is removed, the period term Imf4+ Imf5+ Imf6 is removed, and the remaining part is reconstructed to obtain a 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 index of each target sub-signal pair under different time lags by using a CCM method, and preliminarily judging the authenticity 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 results of determining its optimal embedding dimension E using the G-cure method are shown in FIG. 6. In FIG. 6, the abscissa Embedding dimension is the embedding dimension E, the ordinate G (k) is the G function value, E _L11 、E _L16 、E _L17 Respectively is L ₁₁ 、L ₁₆ 、L ₁₇ Is embedded in the embedded dimension of (a). The determination criteria for the optimal embedding dimension E are as follows: and when the G function is smaller than 0.5, the corresponding E is the optimal embedding dimension. In this embodiment, the optimal embedding dimension e=5.

Step 4-2, determining two variables X and Y of each target sub-signal pair, and setting time delay lambda E [ -c: delta c: c of X relative to Y as preset reason variable and 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 used instead of X (t)); similarly, ρ between Y (t+λ) and X (t) can be obtained _Y→X (Y (t+λ) is used instead of X (t), and X (t) is used instead of Y (t)). In this embodiment, c=50.

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

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

Step 4-4, finding the shadow manifold M _Y Points ony(t) at time t, (E+1) nearest neighbors, obtaining a distance sequence { d ] arranged in order of small to large ₁ ,d ₂ ,...,d _E+1 Time index sequence { t } corresponding to a distance sequence ₁ ,t ₂ ,...,t _E+1 And } wherein,xthe distance between (t) is calculated as follows:

d _i ＝D(y(t _i ),y(t _j ))

where D (a, b) points to the Euclidean distance between the quantities a and b.

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

w _i ＝u _i /U

in the middle of

Normalization factor->

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

The calculation formula is as follows:

w in _i Is weight, t _i Is a time index corresponding to the weight.

Step 4-7, calculating X and

standard pearson correlation coefficient between +.>

As a cross-map correlation coefficient, representing the causal intensity of x→y, the calculation formula is as follows:

step 4-8, p is compared with the alternative data method _X→Y Performing 95% confidence test;

step 4-9, respectively obtaining a series of cross mapping indexes of the X-Y direction and the Y-X direction along with time delay change by the steps 4-2 to 4-8, wherein the time delay corresponding to the maximum rho value is the optimal time delay lambda _X→Y 、λ* _Y→X The judgment is carried out by a computer, and the authenticity of the causal relationship is preliminarily judged according to the positive and negative of the optimal time lag, and the specific mode is as follows: the optimal time lag is less than or equal to 0, and the result variable can better estimate the past value of the reason variable instead of the future value, and the cause and effect of the direction is a real cause and effect relationship; an optimal time lag greater than 0 indicates that the result occurs before the cause, which obviously is contrary to the definition of causal relationships, so the causal relationship of this direction is a pseudo-causal relationship.

In this example, the G-cure method is performed in accordance with "Wang Y, hu F, cao Y, et al improved CCM for variable causality detection in complex systems [ J ]. Control Engineering Practice,2019,83:67-82.

In this example, the substitution data method is carried out according to "headland data and application thereof [ J ]. University of east China, 2011:1-2. But the normalization step in which the mean and euclidean distance are calculated is omitted. Because the CCM method involves convergence, the causal intensity related index theoretically takes whether 1 is reached as the basis of convergence, and the index ρ itself is a normalized coefficient, the detection result is represented by an absolute value instead of a relative value.

With L ₁₆ As X, L ₁₇ As an example of Y, fig. 7 shows the detection result of time-lapse CCM. In fig. 7, the abscissa time lag λ is the time lag between variables, and the ordinate ρ is the cross-map 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 this embodiment,

then L is ₁₆ →L ₁₇ L is true causal relationship ₁₇ →L ₁₆ Is a pseudo-causal relationship. Thus preliminarily judging L ₁₆ Is L ₁₇ Reasons (L) ₁₆ →L ₁₇ ) I.e. L ₁₆ Unidirectional influence L ₁₇ . And the same can obtain the primary causal relationship judgment result between other 39 pairs of target sub-signal pairs.

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

Step 5-1, for each target signal pair, if the optimal time lag of the target sub-signals is greater than 0, judging that no causal relation exists between the target signal pairs, and judging the next target signal pair; if both of the optimal time lags are smaller than or equal to 0, directly executing the step 5-2, and if only one of the optimal time lags is smaller 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 an interval Δn: n E [ N ] _min :ΔN:N _max ]Calculating the length n according to steps 4-3 to 4-8

And Y (t)) And +.>

Cross-map index ρ with X (t) _X→Y 、ρ _Y→X Up to n=n _max Obtain ρ _X→Y Curve sum ρ _Y→X The curve is then followed by step 5-4. In this embodiment, n=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 several intervals of interval Δn: n E [ 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) ^* ) Cross-map index ρ with X (t) _X→Y 、ρ _Y→X Up to n=n _max Obtain ρ _X→Y Curve sum ρ _Y→X The curve is focused on the ρ curve with the optimal time lag smaller than or equal to 0 corresponding direction. In this embodiment, n=471, N _min ＝10，ΔN＝10，Nmax＝450。

Step 5-4, calculating ρ by using Z-transformation method _X→Y 、ρ _Y→X If |Z| > 1.96, and the ρ curve increases with increasing length n until no longer fluctuates, then the ρ curve in that direction is considered to converge, and the causal relationship in that direction is true, otherwise it is not.

In this embodiment, the Z-transform method is implemented according to "Aftab M F, hovd M, sivalingam S.convergent cross-mapping (CCM) based approach for isolating the source of plant-wide disturbances [ C ]//2017IEEE Conference on Control Technology and Applications (CCTA) & IEEE,2017:1492-1498.

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

Then X (t-2) and Y (t) are calculated for different sample lengthsAnd (3) a cross mapping index between Y (t-2) and X (t). Fig. 8 is an optimal time-lapse CCM convergence graph. The abscissa L of fig. 8 is the sample length, and the ordinate ρ is the cross-map index between the corresponding variables. In the present embodiment ρ _L16→L17 The curve converges and the corresponding Z index is 9.814, which is greater than 1.96, therefore L ₁₆ →L ₁₇ Is established, L ₁₆ Is L ₁₇ For reasons of (2). The same applies to the final causal relationship determination between the other 39 pairs of target sub-signal pairs.

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 affects can be considered as the causative variable. In this embodiment, the final causal relationship network is shown in fig. 9. In fig. 9, the behavior is the result variable, and the column is the cause variable. For unidirectional causality, a black square (value 1) represents that this directional causality exists; for bi-directional cause and effect, the black box (value 1) represents the direction cause and effect before the anti-direction cause and effect, and the gray box (value 0.5) represents the direction cause and effect after the anti-direction cause and effect. If both directions are black squares, it means 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 block is a black square, which represents L ₁₆ For the reason L ₁₇ As a result, and this directional causality exists.

The analysis result of the causal relationship of the loop oscillation in the embodiment shows that L ₁₆ Direct or indirect influence on the division of L ₂ 、L ₃ And L ₁₃ Other variables than L ₁₂ Influence (combined with theoretical knowledge, L ₁₆ First act on L ₁₂ And then L ₁₂ Immediate reaction to L ₁₆ ) Thus judge L ₁₆ The source of the TE process oscillation is exactly in line with the location of the artificially added oscillation disturbance in this embodiment. And the overall oscillation propagation path is in accordance with the flow of process flow information and control strategy of the TE process.

The method of the invention realizes the rapid and accurate diagnosis of the oscillation propagation path and the positioning of the oscillation source for the multi-loop oscillation signal in the industrial control system, provides a new thought for the fault diagnosis of the industrial control system oscillation loop by using a causal analysis method, and has certain reference significance and application value.

The foregoing embodiments have described in detail the technical solution and the advantages of the present invention, it should be understood that the foregoing embodiments are merely illustrative of the present invention and are not intended to limit the 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 invention.

Claims (3)

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

(1) Collecting process output signals of each industrial control system oscillation circuit to be analyzed;

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

(3) Denoising and removing period item 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 the cross mapping index of each target sub-signal pair under different time lags by using a CCM method, and preliminarily judging the authenticity of the causal relationship according to the positive and negative of the optimal time lags; the specific process is as follows:

(4-1) determining the optimal embedding dimension E of each target sub-signal by using a G-cure method, wherein the determination standard is as follows: when the G function is smaller 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 time delay lambda epsilon [ -c: delta c: c ] of X relative to Y as preset reason variables and result variables;

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

wherein E is the optimal embedding dimension determined in the step (4-1), and if E of the target sub-signal pair is different, the maximum value of the two is selected as the final E; τ is the time lag, defaulting to 1; the set of X (t) and Y (t) is shadow manifold M corresponding to X (t) and Y (t) _X 、M _Y ；

(4-4) finding shadow manifold M _Y E+1 nearest neighbors of the point y (t) on t at the moment t are obtained to obtain a distance sequence { d ] arranged in order from small to large ₁ ,d ₂ ,...,d _E+1 Time index sequence { t } corresponding to a distance sequence ₁ ,t ₂ ,…,t _E+1 -wherein the distance between x (t) is calculated as:

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 a correlation weight using E+1 nearest neighbors of each point y (t), the weight calculation formula is as follows:

w _i ＝u _i /U

in the method, in the process of the invention,

normalization factor->

(4-6) calculating the estimated value of X with weights

The calculation formula is as follows:

wherein w is _i Is weight, t _i Is a time index corresponding to the weight;

(4-7) calculating X and

standard pearson correlation coefficient ρ therebetween _X→Y As a cross-map correlation coefficient, representing the causal intensity of x→y, the calculation formula is as follows: />

(4-8) Using alternative data method vs ρ _X→Y Performing 95% confidence test;

(4-9) obtaining a series of cross mapping indexes of the X-Y direction changing with time delay by adopting the steps, and obtaining a series of cross mapping indexes of the Y-X direction changing with time delay by adopting the same method from the step (4-2) to the step (4-8); wherein the time lag corresponding to the maximum rho value is the optimal time lag lambda _X→Y 、λ* _Y→X The authenticity of the causal relationship is primarily judged according to the positive and negative of the optimal time lag, and the specific process is as follows:

the optimal time lag is less than or equal to 0, and the result variable can better estimate the past value of the reason variable instead of the future value, and the cause and effect of the direction is a real cause and effect relationship; if the optimal time lag is greater than 0, the result is shown to occur before the reason, and the cause and effect of the direction is a pseudo cause and effect relationship;

(5) Calculating cross mapping indexes of each target sub-signal pair corresponding to the optimal time lag under different sample lengths, judging whether convergence is carried out according to a convergence threshold value, and obtaining causal relation of each signal pair; the specific process is as follows:

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

(5-2) dividing the sample length N of the target signal pair into a plurality of sections with an interval Δn: n E [ N ] _min :ΔN:N _max ]Calculating the length n according to the steps (4-3) to (4-8)

And Y (t)/(t)>

Cross-map index ρ with X (t) _X→Y 、ρ _Y→X Up to n=n _max Obtain ρ _X→Y Curve sum ρ _Y→X Curve, then performing step (5-4);

(5-3) setting the time lag of 0 or less to lambda ^* The sample length N of a signal pair is divided into several intervals of interval Δn: n E [ 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) ^* ) Cross-map index ρ with X (t) _X→Y 、ρ _Y→X Up to n=n _max Obtain ρ _X→Y Curve sum ρ _Y→X A curve, wherein a rho curve with the optimal time lag smaller than or equal to 0 in the corresponding direction is selected;

(5-4) calculating ρ using Z-transform method _X→Y 、ρ _Y→X If |Z| > 1.96, and the ρ curve increases with increasing length n until no fluctuation occurs, then the ρ curve in that direction is considered to converge, and the causal relationship in that direction is established, otherwise it is not established;

(6) And determining an oscillation propagation path and positioning an oscillation source according to the final causality 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 the process output signals of every two oscillating circuits by using a K-order nearest neighbor method;

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

(2-3) reserving a signal pair having a mutual information value greater than 0 as a 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 index of each layer component by using the DFA method, and if all the DFA indexes are greater than 0.5, judging that the component does not contain noise, and directly performing the step (3-3); if the DFA index of a certain layer is less than or equal to 0.5, judging the layer as useless noise, and if a plurality of layers exist, accumulating the layers as final noise;

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

and (3-4) if the signal contains noise, carrying out denoising and cycle term removal on the target signal, if the signal does not contain noise, directly carrying out cycle term removal on the target signal, and reconstructing the rest part after the processing is finished to obtain a target sub-signal.

Images