CN111160776A - Method for detecting abnormal working condition in sewage treatment process by utilizing block principal component analysis - Google Patents

Method for detecting abnormal working condition in sewage treatment process by utilizing block principal component analysis Download PDF

Info

Publication number
CN111160776A
CN111160776A CN201911394625.XA CN201911394625A CN111160776A CN 111160776 A CN111160776 A CN 111160776A CN 201911394625 A CN201911394625 A CN 201911394625A CN 111160776 A CN111160776 A CN 111160776A
Authority
CN
China
Prior art keywords
data
sub
matrix
variables
principal component
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.)
Pending
Application number
CN201911394625.XA
Other languages
Chinese (zh)
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.)
East China University of Science and Technology
Original Assignee
East China University of Science and Technology
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 East China University of Science and Technology filed Critical East China University of Science and Technology
Priority to CN201911394625.XA priority Critical patent/CN111160776A/en
Publication of CN111160776A publication Critical patent/CN111160776A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
    • G06F18/24155—Bayesian classification
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/10—Services
    • G06Q50/26—Government or public services

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Strategic Management (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Tourism & Hospitality (AREA)
  • Economics (AREA)
  • Educational Administration (AREA)
  • Development Economics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Marketing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Entrepreneurship & Innovation (AREA)
  • General Business, Economics & Management (AREA)
  • Evolutionary Biology (AREA)
  • Artificial Intelligence (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Probability & Statistics with Applications (AREA)
  • Game Theory and Decision Science (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)
  • Testing And Monitoring For Control Systems (AREA)

Abstract

The invention relates to a method for identifying abnormal working conditions in a sewage treatment process by utilizing block principal component analysis. The method takes data under the normal working condition of the sewage treatment process as a training data set, blocks variables by calculating mutual information of the variables in the training data set, establishes principal component models in each sub-block according to a principal component analysis method, and establishes T in the principal component model of each sub-block2And (4) statistics and threshold values thereof, establishing a principal component model of the real-time data by using the real-time data, and judging whether a fault occurs according to whether the detection statistics of the real-time data in the whole data space exceeds a control limit. The method can establish the principal component model according to the training data set in the sewage treatment process, identify the fault point and the fault variable, and has the characteristics of high efficiency, rapidness, accuracy and the like.

Description

Method for detecting abnormal working condition in sewage treatment process by utilizing block principal component analysis
Technical Field
The invention belongs to the field of sewage treatment, and particularly relates to a method for detecting abnormal working conditions in a sewage treatment process by utilizing block principal component analysis.
Background
Municipal sewage contains large particle solid suspensions, various pathogens, nitrogen compounds, phosphorus compounds, carbohydrates, etc., and is a very complex mixed liquor. FIG. 1 shows a flow diagram of a wastewater treatment process. The sewage treatment process is a typical multivariable, multi-coupling, long-flow and nonlinear complex process, and a plurality of interference factors and uncertain factors exist at the same time. At present, the method adopted by urban sewage treatment mainly comprises the steps of oxidizing degradable organic matter components in the sewage by virtue of microbial population adsorption and decomposition capacity, and degrading and separating organic matter from the sewage by using complex biological and chemical reactions and physical treatment to purify the sewage. Figure 2 shows a simplified flow diagram of activated sludge. The existing urban sewage treatment system widely uses an activated sludge method and a derived improved process thereof, but the activated sludge can be damaged to a certain degree due to the influence of working environment, water consumption, weather, poisonous water, inlet water quality and water quantity fluctuation, so that the fault of the sewage treatment process is caused, and the quality of outlet water of the whole sewage treatment system is finally influenced to be not up to the standard. Because of the characteristics of long flow, continuity and irreplaceability of the sewage treatment process, once a certain process of the sewage treatment system breaks down, the whole sewage treatment system breaks down, and huge economic loss and great environmental pollution are brought to sewage plants and society. Therefore, it is very important to perform online fault monitoring and fault diagnosis on abnormal conditions in the sewage treatment process, and further take necessary measures to reduce or suppress the occurrence of the abnormal conditions.
In the current field of process monitoring, fault detection methods are generally divided into three categories: analytical model-based methods, knowledge-driven-based methods, and data-driven-based methods. The method based on the analytic model is that the analytic model of the system is utilized, and the acquired input and output signals are added to construct a residual signal, so that the deviation degree between the actual behavior and the expected state of the system is reflected, and the fault detection can be carried out by setting the threshold value of the residual signal. However, the variables, scales, etc. involved in the actual industrial process are very large, and it becomes very difficult to establish a specific analytical model, even cannot be completed. Although the diagnosis and identification of abnormal working conditions of special conditions can be realized by a fault detection method of a mechanism model, the method excessively depends on the mechanism model of biological reaction, and the existing mechanism model involves too many parameters which are difficult to identify, so that the uncertainty of the abnormal working conditions is increased. Meanwhile, the parameters of the model are related to the running state of the actual plant, a unified biochemical mechanism model does not exist at present, and abnormal diagnosis results obtained based on different mechanism models are possibly contradictory, so that the accuracy of fault diagnosis based on the mechanism model is seriously influenced.
The knowledge-based driving method is mainly used for classifying the running state of the system by means of the prior knowledge and expert experience, and fault behaviors can be accurately and intuitively recognized without establishing an accurate mathematical model. But the collection of knowledge and the establishment of an expert system are very difficult, and the reliability and richness of the expert experience are decisive factors for the identification accuracy of the system. In the field of fault diagnosis of sewage treatment, two methods, namely a decision tree method and an expert system, are widely applied in a knowledge-based research method, but the types of the knowledge-based fault diagnosis method are few, sludge load, chemical oxygen demand, SVI (singular value index), dissolved oxygen, pH (potential of hydrogen) and the like are mainly analyzed, a decision tree or an expert database is established through indexes, and a risk assessment model of abnormal working conditions is established. The model can realize risk assessment on abnormal working conditions, but cannot perform quantitative analysis on the abnormal working conditions.
Based on the data driving method, a specific analytical model and rich expert knowledge are not required, variable data obtained in an industrial process can be directly utilized, and fault detection can be carried out on the running state of the system by extracting characteristic information in the data. With the development of artificial intelligence, more fault diagnosis methods are used at the present stage: a deep learning model, a staged analysis model, a BP (Back propagation) neural network model, a principal component analysis model and the like. The deep learning model utilizes the classifier to classify for multiple times, and the accuracy is high. After the fault hypothesis is introduced, the staged analysis model can quickly detect the suspicious elements and then screen the suspicious elements. The BP neural network model utilizes the fault characteristics and the diagnosis target to establish a fault model of the three-layer neural network, and the extracted fault signal is used as an input signal of the BP neural network, so that the diagnosis speed is low, and the accuracy is high. The principal component analysis model has the advantage that a large amount of data can be used for rapid distinction without building a model.
However, most of the actual sewage treatment processes are complex variable non-Gaussian dynamic processes, and the variables are various and are coupled with each other. A wide range of multivariable dynamically varying processes is much more complex than a single process variable, and the current measurements in a local cell depend not only on the past few measurements in the local cell, but also on the current or past measurements in neighboring cells. The conventional Principal Component Analysis (PCA) is a method of performing matrix operation on all variables, the dimensionality of a matrix is usually large, a large amount of time is spent in operation, and once a characteristic value is improperly selected, a large error is caused to the construction of a model and the generation of a principal Component space, so that the conventional principal Component Analysis has the defects of low recognition rate and large calculation amount.
Disclosure of Invention
In view of the above problems, the present invention provides a method for detecting abnormal conditions in a sewage treatment process by using a block principal component analysis. The method divides acquired preprocessed data under normal working conditions into training sub-blocks by utilizing mutual information, establishes principal component models for each training sub-block according to a traditional principal component analysis method, and calculates T for each training sub-model2Statistics and determining T from the kernel density estimate2A threshold value of the statistic; then collecting real-time data, carrying out the same modeling on the real-time data, and calculating the T of each detection sub-model2Statistics and determining T according to Bayesian inference2Statistical failure probability, for T2And weighting the fault probability of the statistics to obtain the detection statistics of the real-time data in the whole data space, and if the detection statistics exceed the control limit, determining that abnormal working conditions possibly occur.
Specifically, the invention provides a method for detecting abnormal working conditions in a sewage treatment process by utilizing block principal component analysis, which is characterized by comprising the following steps of:
the method comprises the following steps: selecting a monitoring variable in the sewage treatment process, and collecting data under normal working conditions to serve as a training data set;
step two: preprocessing training data;
step three: respectively calculating mutual information between every two variables in the preprocessed training data by using a mutual information method;
step four: appointing a threshold value of mutual information, and dividing variables in the training data into different sub-blocks according to a mutual information rule, namely training sub-blocks;
step five: establishing a principal component analysis model, namely a training sub-model, in each training sub-block;
step six: calculating T for each training sub-model2Statistics, determining T for each training sub-block based on kernel density estimates2A threshold value of the statistic;
step seven: acquiring real-time data in the sewage treatment process, preprocessing the real-time data, dividing variables in the preprocessed real-time data into different subblocks according to the partitioning mode determined in the step four, namely detection subblocks, and establishing a principal component analysis model, namely a detection submodel, in each detection subblock;
step eight: calculating T for each detection submodel2Statistics, determining T of each detection submodel according to Bayesian inference2A failure probability of the statistics;
step nine: according to T of each training sub-block2Threshold value of statistic, T of each detection submodel2And calculating the detection statistic of the real-time data in the whole data space and identifying fault points.
In one or more embodiments, in step one, the wastewater treatment process should substantially satisfy the dynamic process in the activated sludge model No.1, which essentially comprises: (1) aerobic growth of heterotrophs, (2) anoxic growth of heterotrophs, (3) aerobic growth of autotrophs, (4) attenuation of heterotrophs, (5) attenuation of autotrophs, (6) ammoniation of soluble organic nitrogen, (7) hydrolysis of adsorbed slow-degrading organic carbon, and (8) hydrolysis of adsorbed slow-degrading organic nitrogen.
In one or more embodiments, in step one, the wastewater treatment process conforms to the long-term model No.1 reference simulation model, comprising two anoxic tanks, three aerobic tanks, and a secondary sedimentation tank.
In one or more embodiments, in step one, the monitoring variable is selected to reflect the operating conditions of the wastewater treatment process.
In one or more embodiments, in step one, the selected monitoring variable is selected from the group consisting of dissolved oxygen concentration, water inflow, sludge return, water outflow, effluent ammonia nitrogen content, chemical oxygen demand, biological oxygen demand, pH, suspended solids concentration, water pressure, and water temperature.
In one or more embodiments, in step two or step seven, preprocessing the data comprises: removing samples with data missing, removing samples with overlarge deviation with other samples, and carrying out zero-averaging on the data to eliminate the influence of different dimensions; wherein zero-averaging the data comprises: suppose there are M sets of data { XmEach set of data is N-dimensional, thus constituting a matrix Xm×nThe formula for zero-averaging the data is:
Figure BDA0002345962270000041
wherein, i is 1,2 … M, j is 1,2 … N.
In one or more embodiments, in step three, two variables x1And x2Mutual information between I (I)1，x2) The following calculations were made:
Figure BDA0002345962270000051
wherein, p (x)1，x2) Is a joint probability density function, p (x)1) And p (x)2) Are each x1And x2The marginal probability density function of (a); wherein, p (x)1，x2)、p(x1) And p (x)2) The following calculations were made:
Figure BDA0002345962270000052
Figure BDA0002345962270000053
Figure BDA0002345962270000054
wherein, f (x)1，x2) Is x1And x2The joint distribution function of (1).
In one or more embodiments, in step four, the process of partitioning is as follows:
suppose there are n samples, each containing m variables, constituting X ∈ Rn×mThe data set of (a);
step a: calculating two variables xi，xjMutual information I (x) between (I ═ 1,2, …, m; j ═ 1,2, …, m)i,xj) Thus obtaining a mutual information matrix I e Rm×m；
Step b, a threshold η of mutual information is given according to experience if I (x)i,xj) Satisfy I (x)i，xj) Equal to η, the variable x is setiAnd xjDivided into one sub-block;
step c: and for the variables of which the mutual information values corresponding to other variables are all lower than the threshold value, collecting the variables into a new sub-block or distributing the variables into other existing sub-blocks.
In one or more embodiments, the threshold η for mutual information is 0.5.
In one or more embodiments, in step five or step seven, the establishing a principal component analysis model includes:
suppose a sub-block contains M sets of data { XmEach data is N-dimensional, thus constituting a matrix Xm×n；
Step 1: for matrix Xm×nPerforming zero equalization processingThe formula is as follows:
Figure BDA0002345962270000055
wherein, i is 1,2 … M, j is 1,2 … N; matrix Xm×nMarking the matrix obtained after zero averaging as X, and solving the covariance matrix S of the matrix X after zero averagingT；
Step 2: find STCharacteristic value λ ofiAnd corresponding unitized orthogonal eigenvectors ai；
And step 3: the feature vector aiAccording to its corresponding characteristic value lambdaiIs arranged in a matrix from large to small [ p ]1，p2，…，pN]；
And 4, step 4: determining dimension k of the submodel, comprising calculating CPV (k) according to:
Figure BDA0002345962270000061
wherein λ isiSetting a CPV threshold value for the ith characteristic value arranged from large to small according to experiments and experiences, and taking a k value when the CPV (k-1) < the CPV threshold value and the CPV (k) > or less than the CPV threshold value as the dimension k of the submodel;
and 5: according to the dimension k determined in the step 4, selecting the matrix [ p ] obtained in the step 31，p2，…，pN]Forming a matrix u by the k groups, namely, reducing the dimension of the matrix X to k dimensions; matrix [ p ] obtained in step 31，p2，…，pN]The first k columns of vectors form a principal component space, and the rest vectors form a residual error space;
step 6: principal component analysis model
Figure BDA0002345962270000062
The following calculations were made:
Figure BDA0002345962270000063
wherein, [ t ]1，t2，…，tk]＝[Xp1，Xp2，…Xpk]，[p1，p2，…，pk]Is the matrix u obtained in step 5, and X is the matrix X obtained in step 1.
In one or more embodiments, in step six or step eight, T2The statistical quantity is calculated by the formula:
T2＝XTPΛ-1PTX，
wherein, T2Is T2Statistic, Λ ═ diag (λ)1，λ2，...λk)，λ1,λ2,…λkFor the eigenvalues λ calculated in step 2 of establishing principal component analysis models described hereiniThe first k eigenvalues arranged from large to small, P is the matrix u obtained in the step 5 of establishing a principal component analysis model described herein, and X is the principal component analysis model calculated in the step 6 of establishing a principal component analysis model described herein
Figure BDA00023459622700000611
In one or more embodiments, the CPV threshold is 85%.
In one or more embodiments, in step six, T is determined2The process of thresholding the statistics is as follows:
t is calculated according to the following formula2Distribution function of statistics
Figure BDA0002345962270000064
Figure BDA0002345962270000065
Wherein the content of the first and second substances,
Figure BDA0002345962270000066
represents T2Row ith and column ith elements of statistic, k represents T2Dimension k, h of the statistic represents the bandwidth; satisfy the requirement of
Figure BDA0002345962270000067
Threshold value of
Figure BDA0002345962270000068
Q (q is 1,2 … k) of (a)
Figure BDA0002345962270000069
T as the training sub-block2Threshold value of statistic
Figure BDA00023459622700000610
The KDE threshold is determined experimentally and empirically.
In one or more embodiments, the bandwidth h is 0.1.
In one or more embodiments, the KDE threshold is 95%.
In one or more embodiments, in step seven, the monitored variables of the real-time data are consistent with the categories and numbers of the monitored variables in the training data set.
In one or more embodiments, in step eight, T of the b-th detector sub-block is calculated according to the following formula2Failure probability of statistics
Figure BDA0002345962270000071
B is 1,2, …, B is the number of sub-blocks:
Figure BDA0002345962270000072
wherein the content of the first and second substances,
Figure BDA0002345962270000073
n and F represent normal and abnormal conditions, respectively;
Figure BDA0002345962270000074
and
Figure BDA0002345962270000075
respectively normal and abnormal processesA priori probability;
Figure BDA0002345962270000076
and
Figure BDA0002345962270000077
the calculation formula of (a) is as follows:
Figure BDA0002345962270000078
Figure BDA0002345962270000079
wherein the content of the first and second substances,
Figure BDA00023459622700000710
is T of the training subblock corresponding to the detector subblock2Threshold value of the statistic, T2Is T of the detector block2Statistics are obtained.
In one or more embodiments, step nine is based on T of each training sub-block2Threshold value of statistic, T of each detection submodel2Statistics and fault probability thereof, detection statistics of real-time data in whole data space
Figure BDA00023459622700000711
The following calculations were made:
Figure BDA00023459622700000712
wherein B is the number of the sub-blocks,
Figure BDA00023459622700000713
and
Figure BDA00023459622700000714
as described in any embodiment herein; if it is
Figure BDA00023459622700000715
Greater than differentPrior probability of constant process
Figure BDA00023459622700000716
It is assumed that an abnormal condition may occur.
Drawings
FIG. 1 is a flow chart of a sewage treatment process;
FIG. 2 is a simplified flow diagram of activated sludge;
FIG. 3 is a general flowchart of the abnormal condition detection method of the sewage treatment process using the block principal component analysis according to the present invention;
FIG. 4 is a simplified diagram of a Long-Term Model No.1 reference Simulation Model (Long-Term Benchmark Simulation Model No.1, BSM 1);
FIG. 5 shows the results of the test of the sewage treatment process by the block principal component analysis (block PCA) method in example 1
Figure BDA0002345962270000081
A statistical quantity graph;
FIG. 6 shows T obtained by detecting abnormal conditions in the wastewater treatment process of example 1 by the conventional Principal Component Analysis (PCA)2A statistical quantity graph;
FIG. 7 shows T obtained by detecting abnormal conditions in the wastewater treatment process in example 1 by Slow Feature Analysis (SFA)2And (5) a statistical quantity graph.
Detailed Description
To make the features and effects of the present invention comprehensible to those skilled in the art, general description and definitions are made below with reference to terms and expressions mentioned in the specification and claims. Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs.
In this context, for the sake of brevity, not all possible combinations of features in the various embodiments or examples are described. Therefore, the respective features in the respective embodiments or examples may be arbitrarily combined as long as there is no contradiction between the combinations of the features, and all the possible combinations should be considered as the scope of the present specification.
The invention discovers that the traditional principal component analysis method is to perform matrix operation on all variables, neglects the correlation among the variables, ensures that the dimensionality of a matrix is usually large, takes a large amount of time during operation, and causes large errors on the construction of a model and the generation of a principal component space once characteristic values are improperly selected. The invention provides a block principal component analysis method on the basis of a principal component analysis method, considers the correlation of variables and overcomes the defects of low recognition rate and large calculated amount of the traditional principal component analysis method.
The invention comprises a method for detecting abnormal working conditions in a sewage treatment process by utilizing block principal component analysis (block PCA method for short). Fig. 3 shows a general flow chart of the method of the invention. In the present invention, the abnormal condition of the sewage treatment process generally refers to an abnormal condition recognized in the art, for example, the abnormal condition includes a fault of increased water intake. In some embodiments, the abnormal condition of the wastewater treatment process may also include an abnormal condition that is specifically identified for the wastewater treatment process to be tested.
The block PCA method of the invention comprises the following steps:
the method comprises the following steps: selecting a monitoring variable in the sewage treatment process, and collecting data under normal working conditions to serve as a training data set;
step two: preprocessing training data;
step three: respectively calculating mutual information between every two variables in the preprocessed training data by using a Mutual Information (MI) method;
step four: specifying a threshold value of MI, and dividing variables in the training data into different sub-blocks (also called subspaces herein) according to mutual information rules, namely training sub-blocks;
step five: establishing a principal component analysis model, namely a training sub-model, in each training sub-block;
step six: calculating T for each training sub-model2Statistics, determining T for each training sub-block based on kernel density estimates2A threshold value of the statistic;
step seven: acquiring real-time data in the sewage treatment process, preprocessing the real-time data, dividing variables in the preprocessed real-time data into different subblocks according to the partitioning mode determined in the step four, namely detection subblocks, and establishing a principal component analysis model, namely a detection submodel, in each detection subblock;
step eight: calculating T in each detection submodel2Statistics, determining corresponding failure probability according to Bayesian inference;
step nine: computing detection statistics of real-time data across data space
Figure BDA0002345962270000091
A point of failure is identified.
These steps will be described in detail below. It is understood that within the scope of the present invention, the above-described technical features of the present invention and the technical features described in detail below (e.g., the embodiments) can be combined with each other to constitute a preferred technical solution.
Firstly, selecting monitoring variables and constructing a training data set
In step one, the sewage treatment process preferably substantially satisfies the dynamic process in Activated sludge model 1 (ASM 1). The activated sludge model No.1 is well known in the art, and the dynamic process mainly comprises the following steps: (1) aerobic growth of heterotrophs, (2) anoxic growth of heterotrophs, (3) aerobic growth of autotrophs, (4) attenuation of heterotrophs, (5) attenuation of autotrophs, (6) ammoniation of soluble organic nitrogen, (7) hydrolysis of adsorbed slow-degrading organic carbon, and (8) hydrolysis of adsorbed slow-degrading organic nitrogen. It will be understood by those skilled in the art that the wastewater treatment process satisfying the activated sludge model No.1 means that the wastewater treatment process includes the above-described 8 dynamic processes. In certain embodiments, the wastewater treatment process follows the flow scheme shown in fig. 1.
In certain embodiments, the wastewater treatment process conforms to the Long-term Model No.1 benchmark Simulation Model (Long-term benchmark Simulation Model No.1, BSM 1). FIG. 4 shows a simplified diagram of a long-term model No.1 baseline simulation model. The long-term model No.1 reference simulation model is well known in the art and comprises two anoxic tanks, three aerobic tanks and a secondary sedimentation tank.
It is understood that, in step one, the selected monitoring variable reflects the operation condition of the sewage treatment process. In the present invention, the monitoring variable refers to the kind of parameter to be collected. Optional monitoring variables include, but are not limited to, dissolved Oxygen concentration, water inflow, sludge return, water outflow, effluent ammonia nitrogen content, Chemical Oxygen Demand (COD), Biological Oxygen Demand (BOD), pH, solids concentration, water pressure, water temperature, total nitrogen concentration in the wastewater, total COD in the wastewater, NH in the wastewater4 -N and NH3 -The concentration of N, the total amount of solids in the wastewater, the concentration of BOD in the wastewater, the concentration of easily degradable substrates in the influent water, the concentration of ammonia nitrogen in the influent water, the dissolved oxygen content of one or more reaction tanks, and the NO of one or more reaction tanks3 -N and NO2 -N concentration, NH of one or more reaction cells4 -N and NH3 -N concentration, concentration of easily degradable substrate in one or more reaction cells, etc.
In some embodiments, the wastewater treatment process conforms to the long-term model No.1 reference simulation model, and the selected monitoring variables include total nitrogen concentration in the wastewater, total COD amount in the wastewater, and NH in the wastewater4 -N and NH3 -The concentration of N, the total solid content in the sewage, the concentration of BOD in the sewage, the concentration of easily-degradable substrate in the inlet water, the concentration of ammonia nitrogen in the inlet water, the dissolved oxygen content in the 3 rd reaction tank, and NO in the 3 rd reaction tank3 -N and NO2 -N concentration, NH of 3 rd reaction cell4 -N and NH3 -N concentration, 4 th reaction tank dissolved oxygen content, 5 th reaction tank NO3 -N and NO2 -N concentration, NH of the 5 th reaction cell4 -N and NH3 -N concentration and concentration of easily degradable substrate in the 5 th reaction tank.
When data under normal working conditions are collected, it is generally required that the data are not less than 200 groups, and preferably, the data contain as many operating conditions as possible, such as dry weather, wet weather, rainy weather and the like. As will be appreciated by those skilled in the art, a set of data is generally referred to herein as data collected at a certain point in time for a selected monitored variable.
The invention takes the collected data under the normal working condition as training data, and the training data form a training data set. It will be understood by those skilled in the art that when performing matrix operations on data, a set of collected data is usually used as a row vector of a matrix, and data belonging to the same monitoring variable is used as a column vector of the matrix. In the present invention, unless otherwise specified, a variable refers to a column vector (corresponding to a monitoring variable) of a matrix, and the number of the column vectors is the dimension of data.
Secondly, preprocessing the data
In step two, the preprocessing of the data generally includes: and removing the samples with data missing, removing the samples with overlarge deviation with other samples, and carrying out zero-averaging on the data to eliminate the influence of different dimensions.
Methods of zero-averaging the data may be known in the art. In some embodiments, zero-averaging the data comprises: assume that there are M sets of sample data { Xm }, each data sample being N-dimensional, from which X is composedm×nThe dimensional matrix, for the data zero-mean formula, is:
Figure BDA0002345962270000111
wherein, i is 1,2 … M, j is 1,2 … N.
The data preprocessing mode in the second step is also suitable for preprocessing the real-time data in the seventh step.
Thirdly, mutual information between the calculated variables
The Mutual Information (MI) method considers both the cross relationship between variables and the high-order statistical Information between variables, and realizes automatic division of blocks.
The MI method is a relatively new statistical technique derived from information science that can quantitatively measure the degree of mutual independence of two variables from the perspective of entropy. Two variables x1And x2The MI in between can be calculated as:
Figure BDA0002345962270000112
wherein, I (x)1，x2) Is MI, p (x)1，x2) Is a joint probability density function, p (x)1) And p (x)2) Are each x1And x2The marginal probability density function of (a); wherein, p (x)1，x2)、p(x1) And p (x)2) The following calculations were made:
Figure BDA0002345962270000113
Figure BDA0002345962270000121
Figure BDA0002345962270000122
wherein, f (x)1，x2) Is x1And x2The joint distribution function of (1).
Fourthly, setting a mutual information threshold value, and blocking the variable
And in the fourth step, a threshold value of mutual information is appointed, and the variables are divided into different sub-blocks by using a mutual information rule according to the mutual information among the calculated variables.
In some embodiments, the specific algorithm steps for blocking with Mutual Information (MI) are as follows:
suppose there are n samples, each containing m variables, constituting X ∈ Rn×mThe data set of (a);
step 1: calculating two variables xi，xjMutual information I (x) between (I ═ 1,2, …, m; j ═ l, 2, …, m)i,xj) Thus, a mutual information matrix I ∈ R can be obtainedm×m；
Step 2, a threshold η of mutual information is given according to experience if I (x)i,xj) Satisfy I (x)i，xj) Equal to η, the variable x is setiAnd xjDivided into one sub-block;
and step 3: for the variable with the mutual information value smaller than the threshold value of the mutual information, the variable can be collected into a new sub-block, and can also be distributed into other existing sub-blocks.
One of the main principles of chunking is that variables in different blocks should be equally independent, while variables in the same block are not necessarily highly correlated.
In step four, the threshold for MI is determined in combination with empirical and actual requirements. In certain embodiments, the threshold value for MI is set to 0.5.
In certain embodiments, the variables are divided into sub-blocks according to the MI blocking method described herein as follows:
X＝[X1,X2,...,Xb]，
wherein b is the number of sub-blocks.
In the invention, the subblocks into which the variable in the training data is divided are called training subblocks, and the subblocks into which the variable in the real-time data is divided are called detection subblocks. In the seventh step, the partitioning of the detection sub-blocks is performed in the partitioning manner established in the fourth step (i.e. the training sub-blocks are partitioned), that is, as will be understood by those skilled in the art, if the variables in the training data are partitioned into n training sub-blocks, the variables in the real-time data are correspondingly partitioned into n detection sub-blocks, and if some variables in the training data are partitioned into the same training sub-block, the variables in the real-time data belonging to the same monitoring variables as those in the training data are correspondingly partitioned into one detection sub-block.
Fifthly, principal component models are respectively constructed in each sub-block
In step five, the method for constructing the principal component model (also called principal component analysis model) may be a method known in the art. In this context, a sub-model refers to a principal component model built in a sub-block.
In certain embodiments, the process of constructing a principal component analysis model is as follows:
assume that the sample data (e.g., a sub-block) contains M sets of data { X }mEach data is N-dimensional, thus constituting a matrix Xm×n；
Step 1: for matrix Xm×nCarrying out zero equalization processing, wherein the zero equalization processing formula is as follows:
Figure BDA0002345962270000131
wherein, i is 1,2 … M, j is 1,2 … N; matrix Xm×nMarking the matrix obtained after zero averaging as X, and solving the covariance matrix S of the matrix X after zero averagingT；
Covariance matrix STThe calculation of (c) is known in the art and can be calculated, for example, by the following formula:
ST＝XTX，
wherein X is the matrix X after zero equalization, and X is X1，x2，…，xn；
Step 2: find STCharacteristic value λ ofiAnd corresponding unitized orthogonal eigenvectors ai；
STCharacteristic value λ ofiAnd corresponding unitized orthogonal vector aiThe calculation of (c) is known in the art and can be calculated, for example, by the following formula:
|λE-ST|＝0，
where E is an identity matrix, from which S can be determinedTCharacteristic value λ of1,λ2,…,λn；
For each lambdaiA basic solution ξ of the homogeneous linear equation set is obtained by1，ξ2，…，ξn：
|λE-ST|X＝0；
By the following formula pairBasic cell line ξ1，ξ2，…，ξnOrthogonalizing and unitizing to obtain a unitized orthogonal vector ai：
a1＝ζ1，
Figure BDA0002345962270000141
……
Figure BDA0002345962270000142
And step 3: unitizing the orthogonal vector aiArranging the eigenvalues into a matrix [ p ] according to the magnitude of the corresponding eigenvalue and the order of the eigenvalue from large to small1,p2,…,pN]；
And 4, step 4: the value of the dimensionality k of the principal component analysis model is determined according to the cumulative variance Contribution (CPV) method using the following equation:
Figure BDA0002345962270000143
wherein λ isiThe ith characteristic value is arranged from big to small, N is the total number of the characteristic values, and k is a positive integer less than or equal to N;
in the CPV method, the dimension k is the pair matrix Xm×nPerforming a target dimension for dimension reduction; a CPV threshold is typically determined based on extensive experimentation and experience; sequentially adding k from 1 and bringing the k into CPV (k) one by one for calculation, and determining the k at the moment as the dimension k of the principal component model when the value of the CPV (k) is greater than or equal to the CPV threshold value for the first time;
in certain embodiments, the CPV threshold is set to 85%, i.e., the value of dimension k is determined in terms of CPV (k) ≧ 85%;
and 5: according to the dimension k determined in the step 4, selecting the matrix [ p ] obtained in the step 31,p2,…,pN]Form a matrix u (i.e., p)1，p2，…，pk]) Namely, the matrix X is data after dimensionality reduction to k dimension; matrix [ p ]1,p2,…,pN]The first k column vectors of (a) constitute the principal component space, the matrix p1,p2,…,pN]The rest vectors form a residual error space; the data is subjected to dimensionality reduction by the method;
step 6: principal component analysis model
Figure BDA0002345962270000144
The following calculations were made:
Figure BDA0002345962270000145
wherein, [ t ]1，t2，…，tk]＝[Xp1，Xp2，…Xpk]，[p1，p2，…，pk]Is the matrix u obtained in step 5, and X is the matrix X obtained in step 1.
The method for establishing the principal component model in the step five is also suitable for establishing the principal component model in the detection sub-block in the step seven.
Sixth, construct T2Statistic, determining T2Threshold value of statistic
In the sixth step, Hotelling T is utilized2Statistic (T for short)2Statistics) to characterize the variation of the sample vector in the principal component space. T is2The way statistics are calculated is known in the art and is calculated as:
T2＝XTPΛ-1PTX，
wherein, T2Is T2Statistic, a ═ diag (λ)1，λ2，…，λi)，λ1,λ2,…λkFor the characteristic value lambda calculated in the process of constructing the principal component modeliThe first k characteristic values are arranged from large to small, and P is a matrix [ P ] obtained in the process of constructing the principal component model1,p2,…,pN]The 'X' is a sub-model of principal component analysis calculated in the process of constructing the principal component model
Figure BDA0002345962270000151
By T2＝XTPΛ-1PTT calculated by X2Statistic T2Is as follows
Figure BDA0002345962270000152
A diagonal matrix of forms.
Herein, T of a submodel2The statistic is also called T of the corresponding sub-block2Statistics are obtained.
Considering that there is no prior knowledge of the data distribution, the residuals may not be gaussian distributed, so there is a Kernel Density Estimation (KDE). KDE is a method of estimating the density of unknown functions in probability theory, and can determine a threshold for residual statistics. The invention utilizes T of training sub-block2Statistics, determining T for each training sub-block based on kernel density estimates2A threshold value for the statistic. The invention trains T of the sub-block2Threshold value of statistic is used as T of corresponding detection subblock2A threshold value for the statistic.
In the sixth step, T is determined by using a nuclear density estimation method2A threshold value for the statistic. In certain embodiments, the present invention determines T using a Gaussian kernel density estimation method2The threshold of the statistic is specifically processed as follows:
t is calculated according to the following formula2Distribution function of statistics
Figure BDA0002345962270000153
Figure BDA0002345962270000154
Wherein the content of the first and second substances,
Figure BDA0002345962270000155
represents T2Row ith and column ith elements of statistic, k represents T2The dimension k, h of the statistic represents the bandwidth.
Figure BDA0002345962270000156
Ordinate of function
Figure BDA0002345962270000157
Is a value between 0 and 1, representing the probability; the abscissa x is an integer between 1 and k, k being T2The dimension k of the statistic.
Computing
Figure BDA0002345962270000161
q is a positive integer from 1 to k, k is T2The dimension k of the statistic; when q (q ═ 1,2 … k) satisfies
Figure BDA0002345962270000162
Threshold value of
Figure BDA0002345962270000163
Then, the value q corresponds to
Figure BDA0002345962270000164
(i.e. T)2Q row and q column elements of statistics) as T of the training sub-block2Threshold value of statistic
Figure BDA0002345962270000165
The bandwidth is an empirical value. In certain embodiments, the bandwidth h is 0.1.
In general, the threshold for a KDE may be determined based on a number of experiments and experience. In certain embodiments, the threshold for KDE is 95%.
Collecting real-time data, preprocessing the real-time data, dividing variables in the preprocessed real-time data into different sub-blocks, namely detection sub-blocks, according to the partitioning mode determined in the step four, and establishing a principal component analysis model in each detection sub-block
It is understood that, in step seven, the monitored variables of the real-time data should be consistent with the categories and the number of the monitored variables of the training data.
In the seventh step, the method for preprocessing the real-time data may be the same as the preprocessing method in the second step, and the method for establishing the principal component analysis model in the detection sub-block may be the same as the method for establishing the training sub-model in the fifth step.
Eighthly, calculating T for each detection submodel2The statistics, i.e. the detection statistics, determine their respective failure probabilities according to Bayesian inference
Step eight, calculating T of the detection submodel2The statistic can be compared with T of the calculation training submodel in step six2The method of statistics is the same.
The invention utilizes T of the detection submodel2Statistics and T of training submodels2Determining T of the detection submodel according to Bayesian inference2The probability of failure of the statistic.
Bayesian inference is typically used to combine the results with probabilities, the bayesian conditional formulation being:
Figure BDA0002345962270000166
where P (a | B) is the probability of event a occurring in the event of event B.
In step eight, T of the B-th detector sub-block (B is 1,2, …, B, and B is the number of sub-blocks)2Failure probability of statistics
Figure BDA0002345962270000171
Can be calculated as:
Figure BDA0002345962270000172
wherein the content of the first and second substances,
Figure BDA0002345962270000173
n and F are normal and abnormal conditions, respectively;
Figure BDA0002345962270000174
and
Figure BDA0002345962270000175
the prior probabilities of normal and abnormal processes, respectively;
Figure BDA0002345962270000176
and
Figure BDA0002345962270000177
the calculation formula of (a) is as follows:
Figure BDA0002345962270000178
Figure BDA0002345962270000179
wherein the content of the first and second substances,
Figure BDA00023459622700001710
is T of the training subblock corresponding to the detection subblock2Threshold value of the statistic, T2Is T of the detector block2Statistics are obtained.
Those skilled in the art will appreciate that the prior probabilities of normal and abnormal processes may be determined based on historical experience with actual conditions. The prior probabilities of the normal and abnormal processes can be determined by adopting a conventional method in the field, for example, the time that the sewage treatment process is in the abnormal working condition (fault) within a period of time (such as within one year) can be counted, and the time in the abnormal working condition is divided by the total time (such as one year), so that the prior probability of the abnormal process is obtained; counting the time of the sewage treatment process in normal working condition (normal operation) within a period of time (such as within one year), and dividing the time in the normal working condition by the total time (such as one year), wherein the time is the prior probability of the normal process.
In some embodiments of the present invention, the substrate is,
Figure BDA00023459622700001711
taking out 99 percent of the raw materials,
Figure BDA00023459622700001712
1 percent of the total weight is taken.
Nine, T for B sub-blocks2Weighting the fault probability of the statistic to obtain the detection statistic of the real-time data in the whole data space, and identifying the fault
The invention is based on the T of each training subblock2Threshold value of statistic, T of each detection submodel2Statistics and fault probability thereof, and calculating detection statistics (detection statistics for short) of real-time data in whole data space
Figure BDA00023459622700001713
Statistics); and if the detection statistic of the real-time data in the whole data space is larger than the control limit of the detection statistic, the abnormal working condition possibly occurs in the sewage treatment process when the real-time data is collected.
In the ninth step, according to T of different sub-blocks2Statistics, their thresholds and failure probabilities, detection statistics of real-time data across data space
Figure BDA00023459622700001714
The following calculations were made:
Figure BDA0002345962270000181
wherein B is the number of the sub-blocks,
Figure BDA0002345962270000182
t for the b-th detector block2The probability of failure of the statistics is,
Figure BDA0002345962270000183
wherein the content of the first and second substances,
Figure BDA0002345962270000184
is T of training sub-block corresponding to the b-th detection sub-block2Threshold value of the statistic, T2Is T of the b-th detector block2Statistics;
if it is
Figure BDA0002345962270000185
If the control limit is larger than the control limit of the detection statistic, the fault is possibly generated in the sewage treatment process when the real-time data is collected. The control limit of the detection statistic is preset according to historical experience of actual working conditions. In certain embodiments, the control limit of the detection statistic is a prior probability of an abnormal process as described herein
Figure BDA0002345962270000186
In some embodiments of the present invention, the substrate is,
Figure BDA0002345962270000187
1 percent of the total weight is taken.
Tenthly, adjusting the model to obtain better detection effect
In certain embodiments, the present invention further comprises the step of: and adjusting the parameters in the first to ninth steps according to the detection result of the abnormal working condition, the actual running condition of the sewage treatment process and/or the detection precision requirement of the abnormal working condition so as to obtain better detection effect. Parameters that may be adjusted include the type and number of monitoring variables, the number of sets of training data, the MI threshold, the CPV threshold, the bandwidth, the KDE threshold, the prior probabilities of normal and abnormal processes, the control limits of the detection statistics, and the like.
For example, the model established in step five can be modified according to the data and accuracy requirements of the actual industrial process, such as changing the MI threshold, the CPV threshold, and the like, so as to obtain better detection effect.
When the abnormal working condition is detected, firstly, the method in the step two of the invention is adopted to preprocess data, then the method in the step three is adopted to calculate the correlation between variables by utilizing mutual information, and the variables are divided into blocks according to the correlation, then the method in the step five is adopted to establish a principal component model in each sub-block, and the T of each model is calculated2Statistics and determining T according to the method described in step six2A threshold value of the statistic; after acquiring real-time monitoring data, establishing the same subblocks and submodels, and calculating T2Statistics and fault probability thereof are calculated and added according to the method of the step nineAnd obtaining the detection statistic of the real-time data in the whole data space, and identifying the fault point.
The invention has the following beneficial effects:
the abnormal working condition detection method is simple, quick and practical, and the data are partitioned by utilizing the correlation among mutual information calculation data. Compared with the traditional principal component analysis method, the method greatly shortens the detection time, considers the correlation among variables and improves the detection accuracy. Mechanism knowledge is not needed in the experimental process, and the industrial process is not needed to be suspended; meanwhile, the method can change the model according to actual needs and data changes so as to realize real-time tracking of the working conditions, obtain better prediction precision and reduce the cost of model maintenance.
The present invention will be specifically described below by way of examples. It should be noted that the following examples are only for illustrating the present invention and should not be construed as limiting the scope of the present invention, and any insubstantial modifications and adaptations by those skilled in the art based on the teachings of the present invention are still within the scope of the present invention. Algorithms and methods not specifically described in the examples are those well known in the art or those described herein. Algebra not explicitly described in the examples have the meaning known in the art or described herein.
Example 1
The method for detecting the abnormal working condition of the sewage treatment process by utilizing the block principal component analysis is described by taking an example of detecting the abnormal working condition of the sewage treatment process as follows, and comprises the following specific steps:
the method comprises the following steps: selecting 15 variables shown in the table 1 as monitoring variables, collecting data under normal working conditions in the sewage treatment process, collecting data of two weeks, and taking 1344 groups as training data sets;
step two: by using
Figure BDA0002345962270000191
Preprocessing training data;
step three: mutual information I (x) of calculated variables1，x2)：
Figure BDA0002345962270000192
Wherein, p (x)1，x2) Is a joint probability density function, p (x)1) And p (x)2) Are each x1And x2The marginal probability density function of (a); wherein, p (x)1，x2)、p(x1) And p (x)2) The following calculations were made:
Figure BDA0002345962270000193
Figure BDA0002345962270000194
Figure BDA0002345962270000195
wherein, f (x)1，x2) Is x1And x2A joint distribution function of (a);
step four: setting the threshold value of mutual information to be 0.5, and partitioning variables according to mutual information rules to obtain 15 training sub-blocks;
step five: establishing a principal component model (training submodel) for each training subblock according to a principal component analysis method, and setting a CPV threshold value to be 85%;
step six: calculating T in each training submodel2Statistics;
determining T using the kernel density estimate described herein2A threshold value of the statistic, wherein the bandwidth h is 0.1, and a threshold value of the kernel density estimation (KDE threshold value) is 0.95;
step seven: in the sewage treatment process, 1344 groups of real-time data are collected, wherein each group of data consists of 15 variables shown in the table 1; acquiring the fault that the water inflow is increased from the 672 th group of data according to the actual operation condition; respectively modeling each group of real-time data according to a modeling mode of training data, dividing each group of real-time data into 15 detection subblocks corresponding to 15 training subblocks one by one, establishing a principal component model (detection submodel) for each detection subblock according to a principal component analysis method, and setting a CPV threshold value to be 85%;
step eight: calculating T of each detection submodel2Statistics, determining T of each detection submodel according to Bayesian inference2Failure probability of statistics
Figure BDA0002345962270000201
Figure BDA0002345962270000202
Wherein the content of the first and second substances,
Figure BDA0002345962270000203
n and F are normal and abnormal conditions, respectively;
Figure BDA0002345962270000204
and
Figure BDA0002345962270000205
the prior probabilities of normal and abnormal processes, respectively, are taken
Figure BDA0002345962270000206
Figure BDA0002345962270000207
Figure BDA0002345962270000208
And
Figure BDA0002345962270000209
the calculation formula of (a) is as follows:
Figure BDA00023459622700002010
Figure BDA00023459622700002011
wherein the content of the first and second substances,
Figure BDA00023459622700002012
is T of the training subblock corresponding to the detector subblock2Threshold value of the statistic, T2Is T of the detector block2Statistics;
step nine: t for each detector sub-block of a set of real-time data2Weighting the failure probability of the statistic to obtain the detection statistic of the group of real-time data in the whole data space
Figure BDA00023459622700002013
According to
Figure BDA00023459622700002014
Whether the prior probability (1% in the embodiment) of the abnormal process is exceeded or not, whether the sewage treatment process fails or not when the set of real-time data is acquired is judged,
Figure BDA0002345962270000211
the following calculations were made:
Figure BDA0002345962270000212
wherein B is the number of sub-blocks (15 in this embodiment),
Figure BDA0002345962270000213
t of the b-th detection sub-block of the set of real-time data2The probability of failure of the statistics is,
Figure BDA0002345962270000214
wherein the content of the first and second substances,
Figure BDA0002345962270000215
t of training sub-block corresponding to the b-th detection sub-block of the set of real-time data2Threshold value of the statistic, T2Is T of the b-th detector sub-block of the set of real-time data2Statistics；
If a certain set of real-time data
Figure BDA0002345962270000216
If the probability is greater than the prior probability (1% in this embodiment) of the abnormal process, it is determined that the abnormal condition occurs in the sewage treatment process when the set of real-time data is collected.
The experimental results of example 1 for detecting abnormal conditions in the sewage treatment process are shown in fig. 5 and table 2.
The abnormal condition detection test was performed by Principal Component Analysis (PCA) and Slow Feature Analysis (SFA) which are conventional in the art using the same training data and detection data as in example 1, and the test results are shown in fig. 6, 7 and table 2.
Table 1: selection of variables
Figure BDA0002345962270000217
Figure BDA0002345962270000221
FIG. 5 shows the results of abnormal condition detection in example 1 by the block PCA method of the present invention
Figure BDA0002345962270000222
And (5) a statistical quantity graph. As for the detection result of the failure, it can be found from fig. 5 that,
Figure BDA0002345962270000223
the statistic value is greatly increased from the 672 th sampling point and far exceeds the specified control limit (example 1 is
Figure BDA0002345962270000224
) This shows that the method of the present invention can find out faults in time and ensure the safety of the production process.
Fig. 6 shows the result of detecting abnormal conditions by a conventional principal component analysis method. Fig. 7 shows the result of abnormal condition detection using the slow signature analysis method. PCA and SFA are widely used fault diagnosis methods in the industry. The best detection results of the block PCA method of the present invention can be found by comparing fig. 5-7.
Table 2 shows the detection rate, false alarm rate, and missing report rate of the experimental results obtained by the block PCA method, the conventional PCA method, and the SFA method of example 1, and the calculation methods of the detection rate, false alarm rate, and missing report rate are as follows:
Figure BDA0002345962270000225
Figure BDA0002345962270000226
Figure BDA0002345962270000227
where TP represents the number of results detected as normal and actually normal, FN represents the number of results detected as failed and actually normal, FP represents the number of results detected as normal and actually failed, and TN represents the number of results detected as failed and actually failed.
Table 2: block PCA, PCA and SFA detection effect comparison
Method of producing a composite material Detection rate False alarm rate Rate of missing reports
Block PCA 95.8333％ 4.1667％ 0
PCA 20.0297％ 0.3720％ 39.8065％
SFA 60.5655％ 0.5208％ 77.7488％
As can be seen from table 2, the false alarm rate of the block PCA method of example 1 is 4.1667%, and the false alarm rate is 0%; and comparing the block PCA, the PCA and the SFA, wherein the block PCA has the lowest missing report rate and the best detection effect.

Claims (10)

1. A method for detecting abnormal working conditions in a sewage treatment process by utilizing block principal component analysis is characterized by comprising the following steps:
the method comprises the following steps: selecting a monitoring variable in the sewage treatment process, and collecting data under normal working conditions to serve as a training data set;
step two: preprocessing training data;
step three: respectively calculating mutual information between every two variables in the preprocessed training data by using a mutual information method;
step four: appointing a threshold value of mutual information, and dividing variables in the training data into different sub-blocks according to a mutual information rule, namely training sub-blocks;
step five: establishing a principal component analysis model, namely a training sub-model, in each training sub-block;
step six: calculating T for each training sub-model2Statistics, rootDetermining T for each training sub-block based on kernel density estimation2A threshold value of the statistic;
step seven: acquiring real-time data in the sewage treatment process, preprocessing the real-time data, dividing variables in the preprocessed real-time data into different subblocks according to the partitioning mode determined in the step four, namely detection subblocks, and establishing a principal component analysis model, namely a detection submodel, in each detection subblock;
step eight: calculating T for each detection submodel2Statistics, determining T of each detection submodel according to Bayesian inference2A failure probability of the statistics;
step nine: according to T of each training sub-block2Threshold value of statistic, T of each detection submodel2And calculating the detection statistic of the real-time data in the whole data space and identifying abnormal working conditions.
2. The method of claim 1,
in the first step, the sewage treatment process meets the dynamic process in the activated sludge model No.1, and comprises the following steps: (1) aerobic growth of heterotrophs, (2) anoxic growth of heterotrophs, (3) aerobic growth of autotrophs, (4) attenuation of heterotrophs, (5) attenuation of autotrophs, (6) ammoniation of soluble organic nitrogen, (7) hydrolysis of adsorbed slow-degrading organic carbon, and (8) hydrolysis of adsorbed slow-degrading organic nitrogen; and/or
In the first step, the sewage treatment process conforms to a long-term type No.1 standard simulation model and comprises two anoxic tanks, three aerobic tanks and a secondary sedimentation tank; and/or
In the first step, the selected monitoring variable can reflect the running condition of the sewage treatment process; and/or
In the first step, the selected monitoring variables are selected from dissolved oxygen concentration, water inflow, sludge reflux, water yield, effluent ammonia nitrogen content, chemical oxygen demand, biological oxygen demand, pH value, solid suspended matter concentration, water pressure and water temperature.
3. The method of claim 1,
in the second or seventh step, the preprocessing of the data comprises: removing samples with data missing, removing samples with overlarge deviation with other samples, and carrying out zero-averaging on the data to eliminate the influence of different dimensions; wherein zero-averaging the data comprises: suppose there are M sets of data { XmEach set of data is N-dimensional, thus constituting a matrix Xm×nThe formula for zero-averaging the data is:
Figure FDA0002345962260000021
wherein, i is 1,2 … M, j is 1,2 … N.
4. The method of claim 1,
in step three, two variables x1And x2Mutual information I (x) between1，x2) The following calculations were made:
Figure FDA0002345962260000022
wherein, p (x)1，x2) Is a joint probability density function, p (x)1) And p (x)2) Are each x1And x2The marginal probability density function of (a); wherein, p (x)1，x2)、p(x1) And p (x)2) The following calculations were made:
Figure FDA0002345962260000023
Figure FDA0002345962260000024
Figure FDA0002345962260000025
wherein, f (x)1，x2) Is x1And x2The joint distribution function of (1).
5. The method of claim 1,
in step four, the process of partitioning is as follows:
suppose there are n samples, each containing m variables, constituting X ∈ Rn×mThe data set of (a);
step a: calculating two variables xi，xjMutual information I (x) between (I ═ 1,2, …, m; j ═ 1,2, …, m)i,xj) Thus obtaining a mutual information matrix I e Rm×m；
Step b, a threshold η of mutual information is given according to experience if I (x)i,xj) Satisfy I (x)i，xj) Equal to η, the variable x is setiAnd xjDivided into one sub-block;
step c: and for the variables of which the mutual information values corresponding to other variables are all lower than the threshold value, collecting the variables into a new sub-block or distributing the variables into other existing sub-blocks.
6. The method of claim 1,
in step five or step seven, the establishing of the principal component analysis model comprises:
suppose a sub-block contains M sets of data { XmEach data is N-dimensional, thus constituting a matrix Xm×n；
Step 1: for matrix Xm×nCarrying out zero equalization processing, wherein the zero equalization processing formula is as follows:
Figure FDA0002345962260000031
wherein, i is 1,2 … M, j is 1,2 … N; matrix Xm×nMarking the matrix obtained after zero averaging as X, and solving the covariance matrix S of the matrix X after zero averagingT；
Step 2: find outSTCharacteristic value λ ofiAnd corresponding unitized orthogonal eigenvectors ai；
And step 3: the feature vector aiAccording to its corresponding characteristic value lambdaiIs arranged in a matrix from large to small [ p ]1，p2，…，pN]；
And 4, step 4: determining dimension k of the submodel, calculating CPV (k) according to:
Figure FDA0002345962260000032
wherein λ isiSetting a CPV threshold value for the ith characteristic value arranged from large to small according to experiments and experiences, and taking a k value when the CPV (k-1) < the CPV threshold value and the CPV (k) > or less than the CPV threshold value as the dimension k of the submodel;
and 5: according to the dimension k determined in the step 4, selecting the matrix [ p ] obtained in the step 31，p2，…，pN]Forming a matrix u by the k groups, namely, reducing the dimension of the matrix X to k dimensions; matrix [ p ] obtained in step 31，p2，…，pN]The first k columns of vectors form a principal component space, and the rest vectors form a residual error space;
step 6: principal component analysis model
Figure FDA0002345962260000041
The following calculations were made:
Figure FDA0002345962260000042
wherein, [ t ]1，t2，…，tk]＝[Xp1，Xp2，…Xpk]，[p1，p2，…，pk]Is the matrix u obtained in step 5, and X is the matrix X obtained in step 1.
7. The method of claim 6,
step six or stepEight middle, T2The statistical quantity is calculated by the formula:
T2＝XTPΛ-1PTX，
wherein, T2Is T2Statistic, Λ ═ diag (λ)1，λ2，…λk)，λ1,λ2,…λkFor the eigenvalues λ calculated in step 2 in claim 6iThe first k eigenvalues in descending order, P being the matrix u obtained in step 5 in claim 6, and X being the principal component analysis model calculated in step 6 in claim 6
Figure FDA0002345962260000043
And/or
In the sixth step, T is determined2The process of thresholding the statistics is as follows:
t is calculated according to the following formula2Distribution function of statistics
Figure FDA0002345962260000044
Figure FDA0002345962260000045
Wherein the content of the first and second substances,
Figure FDA0002345962260000046
represents T2Row ith and column ith elements of statistic, k represents T2Dimension k, h of the statistic represents the bandwidth; satisfy the requirement of
Figure FDA0002345962260000047
Q (q is 1,2 … k) of (a)
Figure FDA0002345962260000048
T as the training sub-block2Threshold value of statistic
Figure FDA0002345962260000049
The KDE threshold is determined experimentally and empirically.
8. The method of claim 1, wherein in step seven, the monitored variables of the real-time data are consistent with the categories and numbers of the monitored variables in the training data set.
9. The method of claim 1,
in step eight, T of the b-th detection sub-block is calculated according to the following formula2Failure probability of statistics
Figure FDA0002345962260000051
B is the number of sub-blocks:
Figure FDA0002345962260000052
wherein the content of the first and second substances,
Figure FDA0002345962260000053
n and F represent normal and abnormal conditions, respectively;
Figure FDA0002345962260000054
and
Figure FDA0002345962260000055
the prior probabilities of normal and abnormal processes, respectively;
Figure FDA0002345962260000056
and
Figure FDA0002345962260000057
the calculation formula of (a) is as follows:
Figure FDA0002345962260000058
Figure FDA0002345962260000059
wherein the content of the first and second substances,
Figure FDA00023459622600000510
is T of the training subblock corresponding to the detector subblock2Threshold value of the statistic, T2Is T of the detector block2Statistics are obtained.
10. The method of claim 9,
in the ninth step, according to T of each training sub-block2Threshold value of statistic, T of each detection submodel2Statistics and fault probability thereof, detection statistics of real-time data in whole data space
Figure FDA00023459622600000511
The following calculations were made:
Figure FDA00023459622600000512
wherein B is the number of the sub-blocks,
Figure FDA00023459622600000513
and
Figure FDA00023459622600000514
as claimed in claim 9; if it is
Figure FDA00023459622600000515
Greater than a priori probability of an abnormal process
Figure FDA00023459622600000516
It is assumed that an abnormal condition may occur.
CN201911394625.XA 2019-12-30 2019-12-30 Method for detecting abnormal working condition in sewage treatment process by utilizing block principal component analysis Pending CN111160776A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911394625.XA CN111160776A (en) 2019-12-30 2019-12-30 Method for detecting abnormal working condition in sewage treatment process by utilizing block principal component analysis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911394625.XA CN111160776A (en) 2019-12-30 2019-12-30 Method for detecting abnormal working condition in sewage treatment process by utilizing block principal component analysis

Publications (1)

Publication Number Publication Date
CN111160776A true CN111160776A (en) 2020-05-15

Family

ID=70559125

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911394625.XA Pending CN111160776A (en) 2019-12-30 2019-12-30 Method for detecting abnormal working condition in sewage treatment process by utilizing block principal component analysis

Country Status (1)

Country Link
CN (1) CN111160776A (en)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112101731A (en) * 2020-08-18 2020-12-18 华南理工大学 Sewage plant online fault monitoring method based on mixed kernel canonical correlation analysis
CN112311662A (en) * 2020-10-22 2021-02-02 南京富岛信息工程有限公司 Intelligent oilfield edge computing gateway
CN112346338A (en) * 2020-10-10 2021-02-09 北京工业大学 Sewage treatment process layered model prediction control method based on fuzzy neural network
CN112417765A (en) * 2020-12-02 2021-02-26 华东理工大学 Sewage treatment process fault detection method based on improved teacher-student network model
CN112763678A (en) * 2020-12-30 2021-05-07 佛山科学技术学院 PCA-based sewage treatment process monitoring method and system
CN113092907A (en) * 2021-04-02 2021-07-09 长春工业大学 System fault detection method based on block slow characteristic analysis
CN113848307A (en) * 2021-11-19 2021-12-28 华南理工大学 Feature extraction principal component analysis online monitoring method for sludge bulking
CN114355846A (en) * 2021-12-07 2022-04-15 华南理工大学 Fault diagnosis method for papermaking sewage treatment process based on SBR simulation model
CN116384158A (en) * 2023-05-26 2023-07-04 广东合诚环境工程有限公司 Sewage treatment equipment operation monitoring method and system based on big data
CN113848307B (en) * 2021-11-19 2024-06-07 华南理工大学 On-line monitoring method for feature extraction principal component analysis aiming at sludge expansion

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105425779A (en) * 2015-12-24 2016-03-23 江南大学 ICA-PCA multi-working condition fault diagnosis method based on local neighborhood standardization and Bayesian inference
CN105676833A (en) * 2015-12-21 2016-06-15 海南电力技术研究院 Power generation process control system fault detection method
CN105955219A (en) * 2016-05-30 2016-09-21 宁波大学 Distributed dynamic process fault detection method based on mutual information
CN109447187A (en) * 2018-12-25 2019-03-08 中南大学 Method of Motor Fault Diagnosis and system

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105676833A (en) * 2015-12-21 2016-06-15 海南电力技术研究院 Power generation process control system fault detection method
CN105425779A (en) * 2015-12-24 2016-03-23 江南大学 ICA-PCA multi-working condition fault diagnosis method based on local neighborhood standardization and Bayesian inference
CN105955219A (en) * 2016-05-30 2016-09-21 宁波大学 Distributed dynamic process fault detection method based on mutual information
CN109447187A (en) * 2018-12-25 2019-03-08 中南大学 Method of Motor Fault Diagnosis and system

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
刘敏 等: "基于邻域保持嵌入-主成分分析的高压 电缆状态数据异常检测及分析", 《科学技术与工程》, pages 1 *
潘珍如: "基于主成分分析和草图的网络异常流量检测研究", 《中国优秀硕士学位论文全文数据库》, pages 4 *

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112101731B (en) * 2020-08-18 2024-03-12 华南理工大学 Sewage plant online fault monitoring method based on mixed kernel typical correlation analysis
CN112101731A (en) * 2020-08-18 2020-12-18 华南理工大学 Sewage plant online fault monitoring method based on mixed kernel canonical correlation analysis
CN112346338B (en) * 2020-10-10 2022-03-11 北京工业大学 Sewage treatment process layered model prediction control method based on fuzzy neural network
CN112346338A (en) * 2020-10-10 2021-02-09 北京工业大学 Sewage treatment process layered model prediction control method based on fuzzy neural network
CN112311662A (en) * 2020-10-22 2021-02-02 南京富岛信息工程有限公司 Intelligent oilfield edge computing gateway
CN112311662B (en) * 2020-10-22 2022-08-30 南京富岛信息工程有限公司 Intelligent oilfield edge computing gateway
CN112417765A (en) * 2020-12-02 2021-02-26 华东理工大学 Sewage treatment process fault detection method based on improved teacher-student network model
CN112417765B (en) * 2020-12-02 2024-02-02 华东理工大学 Sewage treatment process fault detection method based on improved teacher-student network model
CN112763678A (en) * 2020-12-30 2021-05-07 佛山科学技术学院 PCA-based sewage treatment process monitoring method and system
CN113092907A (en) * 2021-04-02 2021-07-09 长春工业大学 System fault detection method based on block slow characteristic analysis
CN113092907B (en) * 2021-04-02 2023-02-03 长春工业大学 System fault detection method based on block slow characteristic analysis
CN113848307A (en) * 2021-11-19 2021-12-28 华南理工大学 Feature extraction principal component analysis online monitoring method for sludge bulking
CN113848307B (en) * 2021-11-19 2024-06-07 华南理工大学 On-line monitoring method for feature extraction principal component analysis aiming at sludge expansion
CN114355846A (en) * 2021-12-07 2022-04-15 华南理工大学 Fault diagnosis method for papermaking sewage treatment process based on SBR simulation model
CN114355846B (en) * 2021-12-07 2023-10-31 华南理工大学 Fault diagnosis method for papermaking sewage treatment process based on SBR simulation model
CN116384158B (en) * 2023-05-26 2023-08-18 广东合诚环境工程有限公司 Sewage treatment equipment operation monitoring method and system based on big data
CN116384158A (en) * 2023-05-26 2023-07-04 广东合诚环境工程有限公司 Sewage treatment equipment operation monitoring method and system based on big data

Similar Documents

Publication Publication Date Title
CN111160776A (en) Method for detecting abnormal working condition in sewage treatment process by utilizing block principal component analysis
CN107025338B (en) Recursive RBF neural network-based sludge bulking fault identification method
US10570024B2 (en) Method for effluent total nitrogen-based on a recurrent self-organizing RBF neural network
CN111126870B (en) Sewage treatment process abnormal condition detection method by utilizing integrated principal component analysis
CN109408774B (en) Method for predicting sewage effluent index based on random forest and gradient lifting tree
CN111291937A (en) Method for predicting quality of treated sewage based on combination of support vector classification and GRU neural network
CN110232256B (en) KPLS (kernel principal component system) and RWFCM (wireless remote control unit) -based sewage treatment process monitoring method
CN112417765B (en) Sewage treatment process fault detection method based on improved teacher-student network model
CN109492265B (en) Wastewater effluent index prediction method based on dynamic nonlinear PLS soft measurement method
CN108665119B (en) Water supply pipe network abnormal working condition early warning method
CN110232062B (en) KPLS (kernel principal component plus minor component plus) and FCM (fiber channel model) -based sewage treatment process monitoring method
Li et al. Ensemble model of wastewater treatment plant based on rich diversity of principal component determining by genetic algorithm for status monitoring
CN112904810B (en) Process industry nonlinear process monitoring method based on effective feature selection
Sanchez-Fernández et al. Fault detection in wastewater treatment plants using distributed PCA methods
CN108830006B (en) Linear-nonlinear industrial process fault detection method based on linear evaluation factor
Xu et al. A complex-valued slow independent component analysis based incipient fault detection and diagnosis method with applications to wastewater treatment processes
CN109934334B (en) Disturbance-based chlorophyll a content related factor sensitivity analysis method
CN111191855A (en) Water quality abnormal event identification and early warning method based on pipe network multi-element water quality time sequence data
Yang et al. Teacher-student uncertainty autoencoder for the process-relevant and quality-relevant fault detection in the industrial process
CN112819087B (en) Method for detecting abnormality of BOD sensor of outlet water based on modularized neural network
CN201330211Y (en) Working parameter self-optimizing simulation system for sewage treatment plant
CN111204867B (en) Membrane bioreactor-MBR membrane pollution intelligent decision-making method
CN115983534A (en) Method and system for evaluating state of sewage treatment process
CN110542748B (en) Knowledge-based robust effluent ammonia nitrogen soft measurement method
CN117170221A (en) Artificial intelligence control system for sewage treatment

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