CN109270907B - Process monitoring and fault diagnosis method based on hierarchical probability density decomposition - Google Patents

Process monitoring and fault diagnosis method based on hierarchical probability density decomposition Download PDF

Info

Publication number
CN109270907B
CN109270907B CN201811243457.XA CN201811243457A CN109270907B CN 109270907 B CN109270907 B CN 109270907B CN 201811243457 A CN201811243457 A CN 201811243457A CN 109270907 B CN109270907 B CN 109270907B
Authority
CN
China
Prior art keywords
fault
probability density
data
value
output
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811243457.XA
Other languages
Chinese (zh)
Other versions
CN109270907A (en
Inventor
应安妮
曾九孙
蔡晋辉
姚燕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Jiliang University
Original Assignee
China Jiliang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Jiliang University filed Critical China Jiliang University
Priority to CN201811243457.XA priority Critical patent/CN109270907B/en
Publication of CN109270907A publication Critical patent/CN109270907A/en
Application granted granted Critical
Publication of CN109270907B publication Critical patent/CN109270907B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B19/00Programme-control systems
    • G05B19/02Programme-control systems electric
    • G05B19/418Total factory control, i.e. centrally controlling a plurality of machines, e.g. direct or distributed numerical control [DNC], flexible manufacturing systems [FMS], integrated manufacturing systems [IMS] or computer integrated manufacturing [CIM]
    • G05B19/4184Total factory control, i.e. centrally controlling a plurality of machines, e.g. direct or distributed numerical control [DNC], flexible manufacturing systems [FMS], integrated manufacturing systems [IMS] or computer integrated manufacturing [CIM] characterised by fault tolerance, reliability of production system
    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05BCONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B2219/00Program-control systems
    • G05B2219/30Nc systems
    • G05B2219/31From computer integrated manufacturing till monitoring
    • G05B2219/31457Factory remote control, monitoring through internet
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02PCLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/02Total factory control, e.g. smart factories, flexible manufacturing systems [FMS] or integrated manufacturing systems [IMS]

Landscapes

  • Engineering & Computer Science (AREA)
  • General Engineering & Computer Science (AREA)
  • Manufacturing & Machinery (AREA)
  • Quality & Reliability (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Testing And Monitoring For Control Systems (AREA)

Abstract

The invention discloses a process monitoring and fault diagnosis method based on hierarchical probability density decomposition. Constructing a layered probability map model, decomposing the joint probability density of all process variables into the product of a conditional probability density and a low-dimensional probability density, and estimating and statistically testing the decomposed probability densities by a kernel conditional density estimation method and a kernel density estimation method respectively; and (3) obtaining an optimal fault value for the input variable and the output variable in each possible fault direction in an iterative mode by using an iterative diagnosis method based on fault reconstruction, and performing corresponding correction to minimize the fault influence on the variables so as to realize the separation of the fault variables. The invention provides reliable and effective technical support for the evaluation of industrial production control behaviors and the diagnosis of fault variables.

Description

Process monitoring and fault diagnosis method based on hierarchical probability density decomposition
Technical Field
The invention belongs to the field of performance evaluation, process monitoring and fault diagnosis in an industrial control system, and particularly relates to a process monitoring and fault diagnosis method based on hierarchical probability density decomposition.
Background
With the rapid popularization of artificial intelligence and the rapid development of computer science, the industrial process changes with the earth, the complexity of the industrial process is increasing day by day, and the importance of process monitoring and fault diagnosis research is higher and higher. The complexity of complex industrial processes is mainly reflected in: the characteristics of the multi-variable system are interfered by various uncertain factors such as raw material components, operation conditions, equipment states and the like, and are difficult to describe by a mechanism model. For example, in complex production processes such as metallurgy and chemical engineering, certain faults can occur in the production process due to the complexity of the composition and structure of equipment, the change of external environment, the aging of the equipment, human misoperation and other factors, the probability of the faults is increased along with the lapse of time, and the destructiveness of major faults is increased. When a fault occurs, if the fault cannot be diagnosed accurately in time, the serious property loss and casualties may be brought.
Therefore, in order to ensure the reliability and safety of the production process, improve the product quality, reduce the production cost, and realize the timely monitoring and accurate positioning of process faults and the separation of fault variables, the method is a problem which needs to be solved urgently.
Typical fault diagnosis methods can be roughly classified into two categories: data-based methods and knowledge-based methods. The data-based method utilizes a large amount of data acquired in production, does not depend or only depends on a small amount of process knowledge, and can realize effective monitoring on faults through analysis of the production data; the main disadvantage is that it is difficult to locate and diagnose a fault after it has been monitored. The knowledge-based method utilizes qualitative knowledge in the production process to diagnose, and can obtain better effect under the condition of complete knowledge; however, due to the high complexity of industrial processes, acquiring complete expert knowledge requires a great deal of effort, which is also a bottleneck limiting the importance of expert knowledge-based methods. The graph model is very suitable for the process of processing uncertainty and incomplete knowledge, can convert the process qualitative knowledge into the prior distribution of a model structure or parameters, and comprehensively calculates probability likelihood or posterior probability with a data sample, thereby realizing the inference of unknown parameters and the diagnosis of faults. Compared with the traditional method, the method of the graph model can more fully utilize precious and available process knowledge and production data, and realize more accurate fault location and diagnosis.
In practical application, the process monitoring and fault diagnosis method based on the hierarchical probability density decomposition can effectively monitor whether faults occur in the industrial production process in time, estimate the fault amplitude and accurately separate fault variables, and has very important practical value for accurately diagnosing the faults in the industrial production process.
Disclosure of Invention
In order to solve the technical problems in the background art, the invention provides a process monitoring and fault diagnosis method based on hierarchical probability density decomposition, which can be used for carrying out fault monitoring on an uncertain industrial process with incomplete process knowledge, estimating a fault amplitude value and further realizing the separation of fault variables.
The technical scheme adopted by the invention is as follows:
step 1, acquiring data of two process variables of an input measurement value and an output measurement value of a process industrial production process under a normal working condition by using a sensor as training data, respectively using the input measurement value and the output measurement value as an input variable and an output variable, carrying out normalization processing on the training data, and calculating the mean value and the variance of the training data;
the data of the process industrial production process specifically comprises output and input measurement values acquired by a sensor arranged on assembly line process operation equipment of the process industrial production process, and the input measurement value is uk,uk-1,uk-2… time series of measurements, the output measurement being yk,yk-1,yk-2… time series of measured values, wherein k represents the ordinal number of the sampling time, the sensors include temperature sensors installed at the outlet of the reboiler, the top and middle of the distillation column, velocity sensors at the feed inlet, the overhead gas and the reflux drum, pressure sensors inside the reboiler, and concentration sensors at the overhead gas
Step 2, identifying the data subjected to normalization processing in the step 1 by using an N4SID subspace identification method to obtain a state space matrix and an initial process state vector of an input and output state space model, and establishing the input and output state space model by taking the process state vector as a state variable to describe the process of the process industrial production;
step 3, constructing a layered probability graph model according to the causal dependence among the three process variables of the input variable, the output variable and the state variable, and decomposing a joint probability density function of all the process variables for monitoring the process state of the process industrial production into a product of a conditional probability density function and a low-dimensional probability density function;
step 4, estimating the data normalized in the step 1 by a non-parameter density estimation method to obtain a decomposed conditional probability density function and a decomposed low-dimensional probability density function, calculating the probability density value of each sampling moment of the training data, and respectively setting the control limit of each probability density;
step 5, acquiring data of the process industrial production process under the working condition to be detected as data to be detected through a sensor, carrying out normalization processing on the data to be detected by using the mean value and the variance of the training data obtained in the step 1, calculating the probability density value of each sampling moment of the data to be detected by using the conditional probability density function and the low-dimensional probability density function which are obtained in the step 4, and monitoring the data to be detected by using the control limit obtained in the step 4;
step 6, if the fault is not monitored, no processing is carried out; if the fault is monitored, calculating by using an iterative diagnosis method based on fault reconstruction to obtain an optimal fault amplitude value;
and 7, positioning and separating the fault variable by using the optimal fault amplitude to complete the monitoring and fault diagnosis of the process industrial production process.
The step 2 specifically comprises:
step 2.1: establishing the following input and output state space model, which is expressed as:
xk=Axk-1+Buk-1+Kek
yk=Cxk+Duk+ek
wherein, ykRepresenting the output measurement value, u, monitored by the sensor at the kth sampling instantk,uk-1Respectively representing the input measured values, x, monitored by the sensors at the k, k-1 th momentk,xk-1Respectively representing the process state vectors at the k, k-1 th instant, ekA random noise vector representing the process industrial production process at the kth moment; in specific implementation, random noise vectors e at all sampling moments are randomly generated in normal distribution;
a, B, C, D and K respectively represent a first state space matrix, a second state space matrix, a third state space matrix, a fourth state space matrix and a fifth state space matrix, wherein the state space matrices are represented by A, B, C, D and K, and the process industrial production process is set to be gradually stable, namely all characteristic values of the first state space matrix A are in a unit circle; k represents the ordinal number of the sampling moment, and k is more than or equal to 1;
step 2.2: and then based on the output measurement yk,yk-1,yk-2…, inputting measured value uk,uk-1,uk-2… identifying and obtaining five state space matrixes and initial process state vector x of the I/O state space model by using the N4SID subspace identification method0(ii) a Defining initial input measurements u0Initial output measurement value y of 00Calculating to obtain a process state vector x as 0k,xk-1,xk-2…. Output measured value yk,yk-1,yk-2… (hereinafter all using the output variable y)kRepresentation), input measured values uk,uk-1,uk-2… (all hereafter with input variable u)kRepresentation) and a process state vector xk,xk-1,xk-2… (hereinafter all using the state variable x)kExpressed) as a process variable of a process industrial process.
In the step 3, the specific process is as follows:
step 3-1, firstly, obtaining the causal dependency relationship among the process variables based on the input and output state space model in the step 2, wherein the process variables comprise the output measured value ykInputting the measured value ukProcess state vector xkThe causal dependency between process variables is specifically: the input measurement value influences the process state vector, and further influences the output measurement value;
then, constructing a layered probability map model, wherein the layered probability map model specifically comprises the following steps: the output measurement value, the input measurement value and the process state vector are represented by nodes, the causal dependency relationship among the process variables is represented by a directed arrow, the input measurement value is used as a father node, points to the process state vector used as an intermediate node, and points to the output measurement value used as a leaf node;
the specific implementation is shown in fig. 1, wherein nodes represent process variables, and directional arrows between the nodes represent causal dependencies between the process variables. The hierarchical probabilistic graph model contains 5 process variables, represented by nodes in the graph, of which the output variable ykAnd inputting variables: u. ofk,uk-1The state variable is as follows: x is the number ofk,xk-1(ii) a The causal relationship between process variables is: input variable u at time k-1k-1Influencing the state variable x at the time k-1k-1And further influences the input variable u at the next time (k-th time)kAnd a state variable xkAnd further influences the output variable y at the k-th timekIn the figure, causal relationships are represented by directional arrows
Step 3-2, establishing the joint probability density of the error items in the input and output state space model, and expressing the joint probability density as p (e)k)=p(yk,xk,uk,xk-1,uk-1),p(ek) Representing the random noise vector e at the kth sampling instantkProbability value of p (y)k,xk,uk,xk-1,uk-1) Representing probability values of all relevant process variables at the kth sampling moment;
the joint probability density is decomposed into the product of three conditional probability densities and one low-dimensional probability density, namely p (e) by using a chain rule and combining a hierarchical probability graph modelk)=p1(yk|xk,uk)p2(xk,uk|xk-1,uk-1)p3(xk-1|uk-1)p4(uk-1),p1(yk|xk,uk) Denotes xk,ukOutput measurement y under known conditionskFirst conditional probability density of (1), p2(xk,uk|xk-1,uk-1) Denotes xk-1,uk-1Process state vector x under known conditionskAnd inputting the measured value ukSecond conditional probability density, p3(xk-1|uk-1) Represents uk-1Process state vector x under known conditionsk-1Third conditional probability density, p4(uk-1) Representing the low dimensional probability density of the input measurement at the (k-1) th sampling instant.
In the step 4, each decomposed probability density is estimated by a nonparametric density estimation method, wherein the conditional probability density is estimated by a kernel conditional density estimation method, the low-dimensional probability density is calculated by a kernel density estimation method, the significance level α is 0.05 when the probability density value is calculated, and 100 x α% quantiles are selected as control limits.
In the step 5, each probability density value of each sampling time of the data to be measured is specifically compared with the corresponding control limit:
a. if the four probability density values are higher than the respective corresponding control limits, the monitored process industrial process is considered to be normal;
b. if probability density p4(uk-1) Above the control limit and with a probability density p1(yk|xk,uk)、p2(xk,uk|xk-1,uk-1)、p3(xk-1|uk-1) If the measured values are all lower than the control limit, the system of the monitored process industrial process is considered to have a fault, and the output measured value is possibly faulty;
c. if probability density p2(xk,uk|xk-1,uk-1)、p3(xk-1|uk-1)、p4(uk-1) Above the control limit, probability density p1(yk|xk,uk) If the output measurement value is lower than the control limit, the output measurement value of the monitored process industrial process is considered to have a fault;
d. if the four probability densities are all lower than the control limit, the input measurement value of the monitored process industrial process is considered to have a fault, and the system and the output measurement value are considered to have faults.
The step 6 comprises the following steps:
step 6-1, setting fault directions of process variables, wherein one fault direction corresponds to one variable type in one process variable, and J different fault directions are totally expressed as { theta [ ]j,j=1,2,…,J},ΘjRepresenting the jth fault direction, the total number of fault directions J being the same as the total number of all species in all process variables; the fault amplitude corresponding to the fault direction is f, and the data containing the fault in the jth fault direction is expressed as sf=s+Θjf, sf being at fault amplitude fData collected under fault conditions, s represents nominal normal data, sf refers to output measured value y under fault conditionsfOr inputting measured values uf(ii) a s is the nominal output measured value y or the nominal input measured value u obtained if the measured value is under the normal working condition;
step 6-2, for input variables, in each possible fault direction implementation, e.g. u1~u6Based on the kernel condition density estimation method, the optimal fault amplitude is iteratively solved for data by adopting the following fixed point iterative optimization algorithm
Figure BDA0001839952010000051
Figure BDA0001839952010000052
In the formula, f (t) represents the fault amplitude obtained by the t-th iteration, f (t +1) represents the fault amplitude obtained by the t + 1-th iteration, the initial value f (0) of the fault amplitude is taken as a zero vector, t is an iteration index, y is the output measured value of the data to be measured, and y is the output measured value of the data to be measurediIs the output measurement in the ith training data, [ theta ]jIndicating the jth fault direction, βiThe coefficient of the ith kernel conditional density estimation function in the kernel conditional density estimation method is sigma, and the sigma is the kernel bandwidth; for convenient representation, the state variable x of the data to be measured and the input measured value u collected under the fault working condition containing the fault amplitude f are usedfThe composed data matrix is wfIs represented by, i.e. wf=[x uf]By wiRepresenting a process state vector x from the ith training dataiAnd inputting the measured value uiConstituent data matrices, i.e. wi=[xiui]I represents the sampling time of the training data, and N represents the total number of the training data;
the fault amplitude obtained by iterating the optimization iterative algorithm to convergence is the optimal fault amplitude
Figure BDA0001839952010000061
Obtaining the optimal fault amplitude value according to the solution
Figure BDA0001839952010000062
First conditional probability density p for data collected under fault conditions including fault amplitude f1(y|x,uf) Performing fault reconstruction by calculating the difference between the input measurement and the optimum fault amplitude
Figure BDA0001839952010000067
Then obtaining nominal normal data, and obtaining a first conditional probability density p after fault reconstruction for the nominal normal data by using a kernel conditional density estimation method1(y | x, u), i.e.
Figure BDA0001839952010000063
Step 6-3, for the output variable, in each possible fault direction implementation, e.g. y1~y6Based on the kernel condition density estimation method, the optimal fault amplitude is iteratively solved for data by adopting the following fixed point iterative optimization algorithm
Figure BDA0001839952010000064
Figure BDA0001839952010000065
Wherein f (t) represents the fault amplitude obtained by the t iteration, f (t +1) represents the fault amplitude obtained by the t +1 iteration, the initial value f (0) of the fault amplitude is taken as a zero vector, t is an iteration index, y is the sum of the error coefficients of the t iteration and the t iteration indexfOutput measured value, y, collected under fault condition containing fault amplitude fiIs the output measurement in the ith training data, [ theta ]jIndicating the jth fault direction, βiThe coefficient of the ith kernel conditional density estimation function in the kernel conditional density estimation method is sigma, and the sigma is the kernel bandwidth; for convenience of representation, a data matrix formed by a state variable x and an input variable u of the data to be measured is represented by w, namely w is equal to (x, u), and a process state vector x containing ith training data is represented by wiiAnd inputting the measured value uiOf data matrices, i.e.wi=[xiui]I represents the sampling time of the training data, and N represents the total number of the training data;
the fault amplitude obtained from the last iteration to convergence is the optimal fault amplitude
Figure BDA0001839952010000066
Obtaining the optimal fault amplitude value according to the solution
Figure BDA0001839952010000071
First conditional probability density p for data collected under fault conditions including fault amplitude f1(yf| x, u) fault reconstruction, i.e. subtracting the optimal fault amplitude from the output measurement value
Figure BDA0001839952010000074
Then obtaining nominal normal data, and obtaining a first conditional probability density p after fault reconstruction for the nominal normal data by using a kernel conditional density estimation method1(y | x, u) is
Figure BDA0001839952010000072
And 6-4, recording each probability density value obtained after the input measurement value and the output measurement value are reconstructed in each possible fault direction.
The step 7 comprises the following steps:
7-1, calculating the average value of the probability density of the data to be detected before the fault reconstruction as M for the probability density function of the monitored fault;
step 7-2, calculating the average value of the probability density of the data to be tested after the fault reconstruction in each possible fault direction based on the probability density value after the fault reconstruction obtained in the step 6-4, and recording the average value as Mj
7-3, selecting the maximum value of the average values of the probability densities of the data to be detected after the fault reconstruction in all possible fault directions, and recording the maximum value as the average valuemax(Mj);
Step 7-4, calculating a fault reconstruction index zeta by adopting a formulaj
Figure BDA0001839952010000073
If fault reconstruction index ζjIf the current direction is smaller than the preset reconstruction threshold, the fault is not reconstructed according to the correct direction, namely the jth direction is not influenced by the fault;
conversely, if the fault reconstruction index ζjThe reconstruction threshold value is larger than or equal to the preset reconstruction threshold value, so that the fault is reconstructed along the correct direction, namely the jth direction is influenced by the fault, and the positioning and the separation of fault variables are realized.
The method comprises the steps of extracting causal dependence among process variables, constructing a layered probability graph model according to the causal dependence, decomposing joint probability density of all the process variables into products of a plurality of conditional probability densities and a low-dimensional probability density, estimating the decomposed probability densities by a kernel conditional density estimation method and a kernel density estimation method respectively, designing statistical tests by using estimated probability density functions, and realizing process monitoring and pre-analysis of fault variables; in order to provide more detailed information for experienced process operators, based on the obtained fault information, an iterative diagnosis method based on fault reconstruction is utilized to obtain an optimal fault value for each possible fault variable in an iterative mode, and the fault value is correspondingly corrected to minimize the fault influence on the variable, so that the fault variables are separated.
For the monitoring process, a hierarchical probability graph model is constructed by utilizing the causal dependence among process variables obtained by process knowledge, so that a combined probability density function is hierarchically decomposed, the decomposed conditional probability density and low-dimensional probability density are estimated by a nonparametric density estimation method, and the fault monitoring is realized by analyzing the distribution of the probability density function; after the fault is monitored, an iterative diagnosis method based on fault reconstruction is applied to the input variable and the output variable in each possible fault direction, and the fault variables are further positioned and separated.
The process monitoring and fault diagnosis method based on the hierarchical probability density decomposition is mainly based on the following considerations: the control limit of the probability distribution is set through the training data collected under the normal working condition, and if the data to be detected is collected under the normal working condition, the probability distribution is in the acceptance domain. Whether the data to be detected has faults or not can be judged by drawing a probability distribution map of the data to be detected, so that fault monitoring is realized; and the fault data is changed on the basis of the normal data, when a fault occurs, the probability of the fault data is inevitably distributed in a rejection region, if the fault amplitude can be obtained, nominal normal data corresponding to the fault data is obtained after correction, and then the probability distribution of the nominal normal data is obtained, the closer the obtained fault amplitude is to the real fault amplitude, the closer the obtained probability distribution of the nominal normal data is to the distribution in a receiving region, so that the positioning and the separation of the fault variables are realized.
Compared with the prior art, the invention has the following beneficial effects:
1. and the causal dependency relationship among the variables is fully mined by utilizing process knowledge and the like, and a hierarchical probability graph model is constructed, so that the process is more visual.
2. The joint probability density function is converted into the product of a plurality of conditional probability density functions and a low-dimensional probability density function, and a nonparametric density estimation method is introduced for solving, so that the operation amount is greatly reduced, and the problem of dimension disaster is solved.
3. The smearing effect (smearing effect) of the fault variable on the normal variable is greatly reduced, and the method has good noise and interference resisting capability and is beneficial to engineering technicians to accurately and effectively confirm the fault variable.
4. The method has the advantages that precious and available process knowledge and production data are utilized more fully, more accurate fault positioning and diagnosis is realized, and reliable and effective technical support is provided for evaluation of industrial production control behaviors and diagnosis of fault variables.
Drawings
FIG. 1 is a hierarchical probability map model corresponding to a state space in accordance with the present invention;
FIG. 2 is a schematic view of a distillation column according to the present invention;
FIG. 3 is a flow chart of the method of the present invention;
FIG. 4 is a process monitoring graph obtained based on the hierarchical probability density decomposition method of the present invention;
FIG. 5 is a diagram of process monitoring derived based on the DICA method for comparison with the method of the present invention;
FIG. 6 is a diagram of a variable fault isolation obtained based on the iterative diagnostic method of fault reconstruction of the present invention.
FIG. 7 is a graph of ICA contribution for comparison with the method of the present invention.
Detailed Description
The invention is further illustrated by the following figures and examples.
The embodiment of the invention and the implementation process thereof are as follows:
the following description will be made in detail of process monitoring and fault diagnosis and separation methods based on data recorded in an operation process, taking a process of purifying butane in an industrial distillation apparatus as an example.
As shown in fig. 2, the commercial distillation process, which employs a tray column containing 33 trays, achieves the purification of butane from a hydrocarbon mixture containing butane, hexane, and propane impurities as the major components. In the figure, circles indicate sensors, T indicates a temperature sensor, F indicates a flow sensor, P indicates a pressure sensor, a indicates a concentration sensor, and arrows indicate the liquid or gas flow direction.
The industrial distillation process is specifically as follows: the reboiler heats the liquid at the bottom of the tower to partially gasify the liquid, the liquid enters the distillation tower from the tower inlet of the reboiler at the bottom of the tower, the vapor and the descending liquid are in countercurrent two-phase contact, the volatile (low boiling point) component in the descending liquid is continuously transferred into the vapor, the difficult volatile (high boiling point) component in the vapor is continuously transferred into the descending liquid, the concentration of the volatile component is higher as the vapor is closer to the tower top, and the difficult volatile component is enriched as the descending liquid is closer to the tower bottom, so that the aim of separating and purifying the butane by the components is fulfilled. The vapor ascending from the tower top enters a condenser, one part of the condensed liquid is returned to the tower top as reflux liquid, and the rest part is sent out as distillate. And feeding a part of liquid flowing out of the tower bottom into a reboiler to heat and return the liquid to the tower, and taking the other part of liquid as kettle residual liquid. The product quality refers to the concentration of hexane at the top of the tower, i.e. the hexane concentration needs to be lower than a preset value.
In order to reduce waste and ensure economic benefit, the concentration of butane in the still residual liquid extracted from the bottom of the tower is lower than a certain specific value. Because there is no corresponding controller at the top and bottom of the tower to adjust the concentration of each component, the whole process can only be controlled by adjusting the feeding amount and the temperature of the feeding material, and the quality of the product can be affected by the change of any variable. To study the process, 8 variables relating to temperature, pressure, flow rate and concentration were selected for analysis, and the descriptions and numbers of these variables for process monitoring and fault diagnosis and isolation are shown in table 1.
TABLE 1
Figure BDA0001839952010000091
Figure BDA0001839952010000101
As shown in fig. 3, in order to apply the process monitoring and fault diagnosis method based on the layered probability density decomposition to monitor, diagnose and separate the butane purification process of the industrial distillation apparatus, the following steps are established:
step 1, collecting 1000 data of input and output two process variables of an industrial distillation process under normal working conditions through a sensor to serve as training data, carrying out normalization processing on the training data, and calculating the mean value and variance of the training data;
the sensors comprise temperature sensors arranged at the outlet of the reboiler, the top and the middle part of the distillation tower, speed sensors at the feed inlet, the overhead gas and the reflux tank, a pressure sensor inside the reboiler and a concentration sensor for the overhead gas; the input variables are the propane impurity concentration and the feed flow; the output variables are the temperature at the top and middle of the distillation column, reflux flow, distillate flow, reboiler pressure and reboiler outlet temperature.
Step 2, identifying the data subjected to the normalization processing in the step 1 by using an N4SID subspace identification method to obtain a state space matrix and an initial process state vector of an input and output state space model, thus obtaining a process variable which is a state variable, and establishing an input and output state space model to describe the industrial distillation process;
the input and output state model is as follows:
xk=Axk-1+Buk-1+Kek
yk=Cxk+Duk+ek
the embodied industrial distillation process comprises 5 process variables, of which the output variable ykAnd inputting variables: u. ofk,uk-1The state variable is as follows: x is the number ofk,xk-1
Step 3, constructing a layered probability graph model according to the causal dependence among the three process variables of the input variable, the output variable and the state variable, and monitoring the joint probability density function p (y) of all the process variables of the industrial distillation process statek,xk,uk,xk-1,uk-1) Decomposition into a conditional probability density function p1(yk|xk,uk)、p2(xk,uk|xk-1,uk-1)、p3(xk-1|uk-1) And a low dimensional probability density function p4(uk-1) The product of (a).
And 4, estimating the data subjected to the normalization processing in the step 1 by using a non-parameter density estimation method to obtain a decomposed conditional probability density function and a decomposed low-dimensional probability density function, calculating the probability density value of each sampling moment of the training data, selecting the significance level α to be 0.05, and selecting 100 x α% quantile as a control limit.
Step 5, collecting 2000 data of the industrial distillation process under the fault working condition as the disturbance of unit material balance and thermodynamic balance caused by the sudden change of the feeding amount of the data to be detected by a sensor, introducing the fault from the 1001 st sampling point, carrying out normalization processing on the data to be detected by using the mean value and the variance of the training data obtained in the step 1, calculating the probability density value of each sampling moment of the data to be detected by using the conditional probability density function and the low-dimensional probability density function obtained in the step 4, and drawing a probability density distribution diagram to monitor the data to be detected by combining the probability density value and the control limit of the training data obtained in the step 4, wherein the fault type of the data to be detected is the disturbance of unit material balance and thermodynamic balance caused;
step 6, as shown in fig. 4, near the 1001 st sampling point, the four probability density values are all lower than the control limit, so that it is considered that at least the input variable of the distillation process is certain to be failed; the black straight line segment in the figure is a control limit, and the part exceeding the control limit in the black dashed line box indicates that a fault is detected.
Then for the input and output variables, in each possible fault direction, i.e. u1,u2,y1,y2,y3,y4,y5,y6Calculating the eight fault directions based on an iterative diagnosis method of fault reconstruction to obtain an optimal fault amplitude value;
step 7, performing fault reconstruction on the first conditional probability density of the first 10, 20, 30 and 50 fault data of the acquired data under the fault working condition by using the optimal fault amplitude to obtain the reconstructed first conditional probability density, calculating a Fault Reconstruction Index (FRI) by using a formula, and drawing, wherein as shown in fig. 6, the fault variable which can be processed by the method shown in fig. 6 is u2,y1,y2,y3And the actual fault condition is met.
For comparison, the same data set is monitored by a traditional dynamic independent principal component analysis (DICA) method, the monitoring result is shown in FIG. 5, a black straight line segment in the graph is a control limit, and a part exceeding the control limit in a black dotted line box indicates that a fault is monitored.
By comparison, DICA can also monitor faults, but the sensitivity is not high in the method provided by the invention, and more fault information cannot be provided; meanwhile, the fault is separated by using the ICA contribution diagram method, and as a result, as shown in FIG. 7, it can be found that the ICA contribution diagram does not realize the correct positioning and separation of fault variables.
In conclusion, the method provided by the invention can complete fault monitoring on an industrial process based on the hierarchical probability density decomposition and the iterative diagnosis method based on fault reconstruction, realize the positioning and on-line separation of fault variables, and effectively improve the sensitivity of fault monitoring and the accuracy of fault information positioning.

Claims (6)

1. A process monitoring and fault diagnosis method based on hierarchical probability density decomposition is characterized by comprising the following steps:
step 1, acquiring data of two process variables of an input measurement value and an output measurement value of a process industrial production process under a normal working condition by using a sensor as training data, respectively using the input measurement value and the output measurement value as an input variable and an output variable, carrying out normalization processing on the training data, and calculating the mean value and the variance of the training data;
step 2, identifying the data subjected to normalization processing in the step 1 by using an N4SID subspace identification method to obtain a state space matrix and an initial process state vector of an input and output state space model, and establishing the input and output state space model by taking the process state vector as a state variable to describe the process of the process industrial production;
step 3, constructing a layered probability graph model according to the causal dependence among the three process variables of the input variable, the output variable and the state variable, and decomposing a joint probability density function of all the process variables for monitoring the process state of the process industrial production into a product of a conditional probability density function and a low-dimensional probability density function;
step 4, estimating the data normalized in the step 1 by a non-parameter density estimation method to obtain a decomposed conditional probability density function and a decomposed low-dimensional probability density function, calculating the probability density value of each sampling moment of the training data, and respectively setting the control limit of each probability density;
step 5, acquiring data of the process industrial production process under the working condition to be detected as data to be detected through a sensor, carrying out normalization processing on the data to be detected by using the mean value and the variance of the training data obtained in the step 1, calculating the probability density value of each sampling moment of the data to be detected by using the conditional probability density function and the low-dimensional probability density function which are obtained in the step 4, and monitoring the data to be detected by using the control limit obtained in the step 4;
step 6, if the fault is monitored, calculating by using an iterative diagnosis method based on fault reconstruction to obtain an optimal fault amplitude value;
the step 6 comprises the following steps:
step 6-1, setting fault directions of process variables, wherein one fault direction corresponds to one variable type in one process variable, and J different fault directions are totally expressed as { theta [ ]j,j=1,2,…,J},ΘjRepresenting the jth fault direction, the total number of fault directions J being the same as the total number of all species in all process variables; the fault amplitude corresponding to the fault direction is f, and the data containing the fault in the jth fault direction is expressed as sf=s+Θjf, sf is data collected under the fault working condition with the fault amplitude value f, and s represents nominal normal data;
6-2, for the input variable, in each possible fault direction, based on the kernel condition density estimation method, iterating the data by the following fixed point iteration optimization algorithm to obtain the optimal fault amplitude
Figure FDA0002478895330000021
Figure FDA0002478895330000022
In the formula, f (t) represents the fault amplitude obtained by the t-th iteration, f (t +1) represents the fault amplitude obtained by the t + 1-th iteration, the initial value f (0) of the fault amplitude is taken as a zero vector, t is an iteration index, y is the output measured value of the data to be measured, and y is the output measured value of the data to be measurediIs the output measurement in the ith training data, [ theta ]jIndicating the jth fault direction, βiThe coefficient of the ith kernel conditional density estimation function in the kernel conditional density estimation method is sigma, and the sigma is the kernel bandwidth; for convenient representation, of the data to be measuredState variable x and input measured value u collected under fault working condition containing fault amplitude ffThe composed data matrix is wfIs represented by, i.e. wf=[x uf]By wiRepresenting a process state vector x from the ith training dataiAnd inputting the measured value uiConstituent data matrices, i.e. wi=[xiui]I represents the sampling time of the training data, and N represents the total number of the training data;
the fault amplitude obtained by iterating the optimization iterative algorithm to convergence is the optimal fault amplitude
Figure FDA0002478895330000023
Obtaining the optimal fault amplitude value according to the solution
Figure FDA0002478895330000024
First conditional probability density p for data collected under fault conditions including fault amplitude f1(y|x,uf) Performing fault reconstruction by calculating the difference between the input measurement and the optimum fault amplitude
Figure FDA0002478895330000025
Then obtaining nominal normal data, and obtaining a first conditional probability density p after fault reconstruction for the nominal normal data by using a kernel conditional density estimation method1(y|x,u);
6-3, for the output variable, in each possible fault direction, based on the kernel condition density estimation method, iterating the data by the following fixed point iteration optimization algorithm to obtain the optimal fault amplitude
Figure FDA0002478895330000026
Figure FDA0002478895330000031
Wherein f (t) represents the fault amplitude obtained by the t-th iteration, and f (t +1) represents the fault amplitude obtained by the t + 1-th iteration, soThe initial value f (0) of the barrier amplitude is taken as a zero vector, t is an iteration index, yfOutput measured value, y, collected under fault condition containing fault amplitude fiIs the output measurement in the ith training data, [ theta ]jIndicating the jth fault direction, βiThe coefficient of the ith kernel conditional density estimation function in the kernel conditional density estimation method is sigma, and the sigma is the kernel bandwidth; for convenience of representation, a data matrix formed by a state variable x and an input variable u of the data to be measured is represented by w, namely w is equal to (x, u), and a process state vector x containing ith training data is represented by wiiAnd inputting the measured value uiOf a data matrix, i.e. wi=[xiui]I represents the sampling time of the training data, and N represents the total number of the training data;
the fault amplitude obtained from the last iteration to convergence is the optimal fault amplitude
Figure FDA0002478895330000032
Obtaining the optimal fault amplitude value according to the solution
Figure FDA0002478895330000033
First conditional probability density p for data collected under fault conditions including fault amplitude f1(yf| x, u) fault reconstruction, i.e. subtracting the optimal fault amplitude from the output measurement value
Figure FDA0002478895330000034
Then obtaining nominal normal data, and obtaining a first conditional probability density p after fault reconstruction for the nominal normal data by using a kernel conditional density estimation method1(y|x,u);
6-4, recording each probability density value obtained after input measurement values and output measurement values are reconstructed in each possible fault direction;
and 7, positioning and separating the fault variable by using the optimal fault amplitude to complete the monitoring and fault diagnosis of the process industrial production process.
2. The method of claim 1 for process monitoring and fault diagnosis based on hierarchical probability density decomposition, wherein: the step 2 specifically comprises:
step 2.1: establishing the following input and output state space model, which is expressed as:
xk=Axk-1+Buk-1+Kek
yk=Cxk+Duk+ek
wherein, ykRepresenting the output measurement value, u, monitored by the sensor at the kth sampling instantk,uk-1Respectively representing the input measured values, x, monitored by the sensors at the k, k-1 th momentk,xk-1Respectively representing the process state vectors at the k, k-1 th instant, ekA random noise vector representing the process industrial production process at the kth moment;
a, B, C, D and K respectively represent a first state space matrix, a second state space matrix, a third state space matrix, a fourth state space matrix and a fifth state space matrix, and the process industrial production process is set to be gradually stable, namely all characteristic values of the first state space matrix A are in a unit circle; k represents the ordinal number of the sampling moment, and k is more than or equal to 1;
step 2.2: and then based on the output measurement yk,yk-1,yk-2…, inputting measured value uk,uk-1,uk-2… identifying and obtaining five state space matrixes and initial process state vector x of the I/O state space model by using the N4SID subspace identification method0(ii) a Defining initial input measurements u0Initial output measurement value y of 00Calculating to obtain a process state vector x as 0k,xk-1,xk-2…, outputting the measured value yk,yk-1,yk-2…, inputting measured value uk,uk-1,uk-2… and a process state vector xk,xk-1,xk-2… as a process variable in a process industrial process.
3. The method of claim 1 for process monitoring and fault diagnosis based on hierarchical probability density decomposition, wherein: in the step 3, the specific process is as follows:
step 3-1, firstly, obtaining the causal dependency relationship among the process variables based on the input and output state space model in the step 2, wherein the process variables comprise the output measured value ykInputting the measured value ukProcess state vector xkThe causal dependency between process variables is specifically: the input measurement value influences the process state vector, and further influences the output measurement value;
then, constructing a layered probability map model, wherein the layered probability map model specifically comprises the following steps: the output measurement value, the input measurement value and the process state vector are represented by nodes, the causal dependency relationship among the process variables is represented by a directed arrow, the input measurement value is used as a father node, points to the process state vector used as an intermediate node, and points to the output measurement value used as a leaf node;
step 3-2, establishing the joint probability density of the error items in the input and output state space model, and expressing the joint probability density as p (e)k)=p(yk,xk,uk,xk-1,uk-1),p(ek) Representing the random noise vector e at the kth sampling instantkProbability value of p (y)k,xk,uk,xk-1,uk-1) Representing probability values of all relevant process variables at the kth sampling moment;
the joint probability density is decomposed into the product of three conditional probability densities and one low-dimensional probability density, namely p (e) by using a chain rule and combining a hierarchical probability graph modelk)=p1(yk|xk,uk)p2(xk,uk|xk-1,uk-1)p3(xk-1|uk-1)p4(uk-1),p1(yk|xk,uk) Denotes xk,ukOutput measurement y under known conditionskFirst conditional probability density of (1), p2(xk,uk|xk-1,uk-1) Denotes xk-1,uk-1Process state vector x under known conditionskAnd inputting the measured value ukSecond conditional probability density, p3(xk-1|uk-1) Represents uk-1Process state vector x under known conditionsk-1Third conditional probability density, p4(uk-1) Representing the low dimensional probability density of the input measurement at the (k-1) th sampling instant.
4. The method as claimed in claim 1, wherein in step 4, each decomposed probability density is estimated by a non-parametric density estimation method, wherein the conditional probability density is estimated by a kernel conditional density estimation method, the lower dimensional probability density is calculated by a kernel density estimation method, the significance level α is 0.05 when a control limit is set during the calculation of the probability density values, and 100- α% quantiles are selected as the control limit.
5. The method of claim 3 for process monitoring and fault diagnosis based on hierarchical probability density decomposition, wherein: in the step 5, each probability density value of the data to be measured is specifically compared with the corresponding control limit:
a. if the four probability density values are higher than the respective corresponding control limits, the monitored process industrial process is considered to be normal;
b. if probability density p4(uk-1) Above the control limit and with a probability density p1(yk|xk,uk)、p2(xk,uk|xk-1,uk-1)、p3(xk-1|uk-1) If the measured values are all lower than the control limit, the system of the monitored process industrial process is considered to have a fault, and the output measured value is possibly faulty;
c. if probability density p2(xk,uk|xk-1,uk-1)、p3(xk-1|uk-1)、p4(uk-1) Above the control limit, probability density p1(yk|xk,uk) If the output measurement value is lower than the control limit, the output measurement value of the monitored process industrial process is considered to have a fault;
d. if the four probability densities are all lower than the control limit, the input measurement value of the monitored process industrial process is considered to have a fault, and the system and the output measurement value are considered to have faults.
6. The method for process monitoring and fault diagnosis based on hierarchical probability density decomposition of claim 1, wherein: the step 7 comprises the following steps:
7-1, calculating the average value of the probability density of the data to be detected before the fault reconstruction as M for the probability density function of the monitored fault;
step 7-2, calculating the average value of the probability density of the data to be tested after the fault reconstruction in each possible fault direction based on the probability density value after the fault reconstruction obtained in the step 6-4, and recording the average value as Mj
7-3, selecting the maximum value of the average values of the probability densities of the data to be detected after the fault reconstruction in all possible fault directions, and recording the maximum value as max (M)j);
Step 7-4, calculating a fault reconstruction index zeta by adopting a formulaj
Figure FDA0002478895330000051
If fault reconstruction index ζjIf the current direction is smaller than the preset reconstruction threshold, the fault is not reconstructed according to the correct direction, namely the jth direction is not influenced by the fault;
conversely, if the fault reconstruction index ζjThe reconstruction threshold value is larger than or equal to the preset reconstruction threshold value, so that the fault is reconstructed along the correct direction, namely the jth direction is influenced by the fault, and the positioning and the separation of fault variables are realized.
CN201811243457.XA 2018-10-24 2018-10-24 Process monitoring and fault diagnosis method based on hierarchical probability density decomposition Active CN109270907B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811243457.XA CN109270907B (en) 2018-10-24 2018-10-24 Process monitoring and fault diagnosis method based on hierarchical probability density decomposition

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811243457.XA CN109270907B (en) 2018-10-24 2018-10-24 Process monitoring and fault diagnosis method based on hierarchical probability density decomposition

Publications (2)

Publication Number Publication Date
CN109270907A CN109270907A (en) 2019-01-25
CN109270907B true CN109270907B (en) 2020-07-28

Family

ID=65194290

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811243457.XA Active CN109270907B (en) 2018-10-24 2018-10-24 Process monitoring and fault diagnosis method based on hierarchical probability density decomposition

Country Status (1)

Country Link
CN (1) CN109270907B (en)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109933953B (en) * 2019-04-11 2023-04-07 东南大学 Power distribution network switch state identification method based on probability map model
DE102019112099B3 (en) * 2019-05-09 2020-06-18 Dürr Systems Ag Monitoring method for an application system and corresponding application system
DE112020002294A5 (en) 2019-05-09 2022-03-03 Dürr Systems Ag Analysis methods and devices for this
JP2022531886A (en) * 2019-05-09 2022-07-12 デュール システムズ アーゲー Analytical method and equipment for it
US11928628B2 (en) 2019-05-09 2024-03-12 Dürr Systems Ag Method for checking workpieces, checking facility and treatment facility
CN112651534A (en) * 2019-10-10 2021-04-13 顺丰科技有限公司 Method, device and storage medium for predicting resource supply chain demand
CN112818564B (en) * 2021-02-26 2022-07-15 北京航空航天大学 Multi-process manufacturing system functional health state modeling method based on coupling operation factors

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012050262A1 (en) * 2010-10-15 2012-04-19 한국전력공사 Method and system for monitoring the performance of plant instruments using ffvr and glrt
CN105700518A (en) * 2016-03-10 2016-06-22 华中科技大学 Fault diagnosis method during industrial process
CN107272625A (en) * 2017-07-12 2017-10-20 温州大学 A kind of industrial process method for diagnosing faults based on bayesian theory
CN108388232A (en) * 2018-03-20 2018-08-10 江南大学 A kind of operational mode fault monitoring method of crude oil desalting process

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012050262A1 (en) * 2010-10-15 2012-04-19 한국전력공사 Method and system for monitoring the performance of plant instruments using ffvr and glrt
CN105700518A (en) * 2016-03-10 2016-06-22 华中科技大学 Fault diagnosis method during industrial process
CN107272625A (en) * 2017-07-12 2017-10-20 温州大学 A kind of industrial process method for diagnosing faults based on bayesian theory
CN108388232A (en) * 2018-03-20 2018-08-10 江南大学 A kind of operational mode fault monitoring method of crude oil desalting process

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Monitoring of dynamic process using hierarchical probability density decomposition;Anni Ying等;《2017 Chinese Automation Congress (CAC)》;20180101;第4793-4797页 *

Also Published As

Publication number Publication date
CN109270907A (en) 2019-01-25

Similar Documents

Publication Publication Date Title
CN109270907B (en) Process monitoring and fault diagnosis method based on hierarchical probability density decomposition
EP0906593B1 (en) Industrial process surveillance system
TWI648609B (en) Program monitoring system and method
US5586066A (en) Surveillance of industrial processes with correlated parameters
Lau et al. Fault diagnosis of Tennessee Eastman process with multi-scale PCA and ANFIS
Wise et al. The process chemometrics approach to process monitoring and fault detection
US7640145B2 (en) Automated model configuration and deployment system for equipment health monitoring
D’Angelo et al. A new fault classification approach applied to Tennessee Eastman benchmark process
Lahdhiri et al. Supervised process monitoring and fault diagnosis based on machine learning methods
Hajihosseini et al. Fault detection and isolation in the challenging Tennessee Eastman process by using image processing techniques
CN112904810B (en) Process industry nonlinear process monitoring method based on effective feature selection
CN112000081B (en) Fault monitoring method and system based on multi-block information extraction and Mahalanobis distance
EP3712728A1 (en) Apparatus for predicting equipment damage
KR20210017651A (en) Method for Fault Detection and Fault Diagnosis in Semiconductor Manufacturing Process
CN111639304B (en) CSTR fault positioning method based on Xgboost regression model
CN112199409A (en) Method and device for monitoring real-time working condition of catalytic reforming device
Madakyaru et al. Monitoring distillation column systems using improved nonlinear partial least squares-based strategies
CN116380445A (en) Equipment state diagnosis method and related device based on vibration waveform
CN109683594B (en) Method for accurately identifying and positioning abnormal variable
Mao et al. Comparative study on prediction of fuel cell performance using machine learning approaches
Liu et al. Structured sequential Gaussian graphical models for monitoring time-varying process
Dong et al. Dynamic-inner canonical correlation analysis based process monitoring
Huang et al. Distributed SFA-CA monitoring approach for nonstationary plant-wide process and its application on a vinyl acetate monomer process
Zhou et al. Overview of fault detection and identification for non-linear dynamic systems
KR102486463B1 (en) Method and Apparatus for Real Time Fault Detection Using Time series data According to Degradation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant