CN105629958A - Intermittence process fault diagnosis method based on sub-period MPCA-SVM - Google Patents

Abstract

The invention discloses an intermittence process fault diagnosis method based on a sub-period MPCA-SVM, and relates to the field of fault diagnosis based on pattern recognition. The method comprises the steps: firstly carrying out the unfolding of three dimensional data of a fermentation process, and carrying out slicing in a time direction; secondly carrying out the rough dividing and fine dividing of a time period for the intermittence process through employing an MPCA; and finally building an MPCA monitoring model and an SVM diagnosis model in each sub-period. The online fault diagnosis comprises the steps: carrying out the processing of the collected data according to a model, calculating a statistic quantity, and comparing the statistic quantity with a control limit. The production is carried out normally if the limit is exceeded, or else, the data is substituted into the SVM diagnosis model of the corresponding time period for fault diagnosis. The method just carries out the filling of data in a time period when a fault happens, and reduces the impact on the accuracy of the SVM fault diagnosis from manual excessive filling of known data. Meanwhile, the method also reduces the number of models, and solves a problem that a diagnosis process is complex because of the frequent updating of the model.

Description

A kind of batch process method for diagnosing faults based on sub-period MPCA-SVM

Technical field

The present invention relates to the fault diagnosis technology field based on pattern recognition, particularly to one for batch process on-line fault diagnosis technology. Namely the method based on pattern recognition of the present invention is the concrete application in typical batch process penicillin fermentation process malfunction monitoring.

Background technology

Batch process has the features such as specific function, high added value, small lot, multi items because of its product so that its proportion shared aborning is increasing, therefore that the fault diagnosis of batch process is also more and more important. But batch process has the features such as dynamic, strong nonlinearity and period characteristic and the high impact being vulnerable to the factors such as environment with product quality of operation complexity so that the research of its fault diagnosis is faced more challenges.

Fault diagnosis for batch process, method the more commonly used at present has the method for contribution plot and pattern recognition, some scholar utilizes multidirectional principal component analysis (Multi-wayprincipalcomponentanalysis, MPCA) batch process is carried out on-line monitoring at times, and utilize the method for contribution plot to review fault variable, but the method adopts normal data to carry out fault diagnosis, can not really react the information of fault, and have ignored the dependency between variable, single argument fault can only be diagnosed, and the method for pattern recognition can overcome the deficiency of MPCA contribution plot. support vector machine (SupportVectorMachine, SVM) has stronger nonlinear system study ability because of it under small sample so that SVM becomes a kind of mode identification method being widely used. some scholar utilizes MPCA and SVM that whole batch process is set up Off-line faults diagnosis model, the unknown data of artificial filling process is inevitably needed in diagnostic phases, but the data filled often there are differences with real process data, influences whether the accuracy rate of fault diagnosis to a certain extent, in order to solve this problem, sliding time window setting technique and SVM are conjointly employed in the fault diagnosis of batch process by other scholar, though can solve because a large amount of data of filling bring the problem of accuracy rate of diagnosis, but there is model modification shortcoming frequently so that diagnosis process is loaded down with trivial details.

Summary of the invention

In order to overcome the deficiencies in the prior art, the invention provides a kind of batch process on-line fault diagnosis method based on sub-period MPCA-SVM. Utilize MPCA to divide batch process, the search volume of fault diagnosis is confined to one by one in the specific period, and each period is set up fault detection and diagnosis model. The method is without filling the unknown data of whole production process, and the data of the period that is filled only with breaking down, therefore can reduce owing to artificial fills the impact that the accuracy rate of SVM fault diagnosis is brought by too much unknown data. It is also possible to reduce the quantity of modeling, thus solving the problem that the diagnosis process caused because of Renewal model continually is complicated.

Present invention employs following technical scheme and realize step:

A. the off-line modeling stage

1) gathering the historical data under sweat nominal situation, the I lot data under the same sweat same process that described historical data X is obtained by off-line test is constituted, X=(X₁,X₂,...,X_I)^T, wherein X_i(i=1,2 ..., I) represent the i-th lot data. Each batch comprises K sampling instant, i.e. X_i=(X_i,1,X_i,2,...,X_i,K), wherein X_i,kRepresent the data that i-th batch of kth sampling instant gathers. Each sampling instant gathers J process variable, i.e. X_i,k=(x_i,k,1,x_i,k,2,...,x_i,k,J), wherein x_i,k,jThe measured value of the jth process variable of kth sampling instant in representing i-th batch;

2) being standardized historical data X processing, processing mode is as follows:

Average and the standard variance of all process variables, the wherein average of the jth process variable of kth sampling instant is engraved when first calculating historical data X allComputing formula beWherein x_i,k,jThe measured value of the jth process variable of kth sampling instant in representing i-th batch, k=1 ..., K, j=1 ..., J; The standard variance s of the jth process variable of kth sampling instant_k,jComputing formula be, $s_{k,j}=\sqrt{\frac{1}{I-1}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{(x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j})}^2},$ K=1 ..., K, j=1 ..., J;

Then historical data X being standardized, wherein the standardized calculation formula of the jth process variable of kth sampling instant is as follows in i-th batch:

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}---\left(1\right)$

Wherein, i=1 ..., I, j=1 ..., J, k=1 ..., K;

3) X is by step 2) standardization obtain new two-dimensional matrix X', this matrix has (K �� J) individual column vector, i.e. X'=(X₁',X'₂,...,X'_K��J), wherein jth column vector X'_j=(X'_j,1,...,X'_j,K)^T, X'_j,k=(X'_j,k,1,...,X'_j,k,I)^T, wherein X'_j,k,iRepresent through step 2) value corresponding in i-th batch of jth process variable kth sampling instant after standardization, wherein i=1 ..., I, j=1 ..., J, k=1 ..., K;

4) multidirectional pivot analysis MPCA method is utilized to extract the main constituent in each moment in X'. If extracting kth sampling instant data X'_kMain constituent, concrete step is as follows:

4.1) two-dimensional matrix X' is obtained_kCovariance matrix COV;

$COV=\frac{1}{I-1}{X}_k^{\′T}{X}_k^{\′}---\left(2\right)$

4.2) Eigenvalues Decomposition to Matrix C OV;

COV=V �� V^T(3)

Formula (3) discloses the incidence relation of covariance matrix, wherein ��=(��₁,��₂,��,��_v) for diagonal matrix, v is two-dimensional matrix X'_kEigenvalue number, �� comprises the non-negative factual investigation (�� that m amplitude is successively decreased₁�ݦ�₂�ݡ��ݦ�_m>=0). V is orthogonal matrix (V^��V=E, E are unit battle array).

In definition, by MPCA method, matrix X'_kIt is equivalent to X'_k=QP^��+ R, and be updated to formula (3) and obtain:

$S=\frac{1}{I-1}{\mathrm{PQ}}^{T}QPQ---\left(4\right)$

Wherein, Q is score matrix, and P is load matrix, and R is residual matrix.

Corresponding (3) and (4) are every, can obtain:

$\begin{array}{c}P=V\\ \mathrm{\Λ}=\frac{1}{I-1}{Q}^{T}Q\end{array}---\left(5\right)$

If the sum equation more than one threshold value 0.85 of front A main constituent, then before just extracting, A pivot is used as aggregative indicator, then original J dimension space just becomes A dimension and A��J.

4.3) score matrix Q is obtained;

In order to optimally obtain the variable quantity of data, minimize random noise to the PCA impact produced simultaneously, retain the load vector corresponding with A eigenvalue of maximum, then X'_kProjection information at lower dimensional space is included in score matrix:

Q=X'_kP(6)

5) utilize MPCA that batch process is carried out Time segments division:

5.1) period slightly divides:

When neighbouring sample point has identical main constituent number (Q matrix column number), just these sampled points are divided into the same period, if have the Period Length L of identical main constituent number less than whole batch process 1/10 time, in this period data is grouped into that the adjacent main constituent number difference in left and right is relatively small one section. If this Period Length is less than the 1/10 of whole sweat and when the main constituent number of this period is equal with the main constituent number difference of adjacent two periods, then take main constituent number that front and back period main constituent number is this section respectively and calculate its contribution rate, this section was grouped in the period that contribution rate change is relatively small.

Final thick division can obtain F period, is expressed as S₁,S₂,��,S_F��

5.2) period carefully divides: at the arbitrarily thick period S divided_f(f=1,2 ..., F) in, utilize between load matrix angle information and range information to define similarity measurement formula (7).

$d_{b,c}^{S_f}=\underset{l=1}{\overset{a_f}{\mathrm{\Σ}}}{\mathrm{\γ}}_l\frac{\left|\right|P_{bl}^T{P}_{cl}\left|\right|}{\left|\right|P_{bl}\left|\right|\·\left|\right|P_{cl}\left|\right|}+(1-{\mathrm{\γ}}_l)e^{-\left|\right|P_{bl}-P_{cl}\left|\right|}---\left(7\right)$

${\mathrm{\γ}}_l=\frac{1}{l}/\underset{h=1}{\overset{a_f}{\mathrm{\Σ}}}\frac{1}{h},l=h=1,2,...,a_f---\left(8\right)$

Wherein b and c represents at period S_fTwo sampling times of interior arbitrary neighborhood, a_fFor period S_fInterior sampled point number, l and h represents that sampled point 1 arrives a_fIn any one value, | | P_bl| | with | | P_cl| | it is the modulus value of the load matrix of b and c sampling instant respectively,Represent the distance of b and c sampling instant load matrix, r_lFor weight coefficient, in order to emphasize the different importances of different projecting direction,Represent a in two timeslice load matrix_fThe weighted sum of the included angle cosine value of individual projecting direction,Represent the similarity of b and c sampling instant load matrix.

G period is got in final refinement, is expressed as S₁',S'₂,...,S'_G, and calculate any one sub-period S'_g(g=1,2 ..., G) obtain average load matrix

6) at S₁',S'₂,...,S'_GIn set up each sub-period MPCA monitoring model, such as formula (9).

Being respectively provided with close load matrix for each sub-period, therefore the MPCA model of each sub-period all can adopt average load matrixDescribe. For anyon period S'_gIn the MPCA model of any t can be expressed as:

$\left\{\begin{array}{c}{q}_t=S_g^{\′}{\stackrel{\‾}{P}}_{S_g^{\′}}\\ \widehat{S_g^{\′}}=q_t{\left({\stackrel{\‾}{P}}_{S_g^{\′}}\right)}^T\\ r_t=S_g^{\′}-\widehat{S_g^{\′}}\end{array}\right.---\left(9\right)$

Wherein q_tIt is sub-period S'_gIn the load vector of t,Sub-period S'_gCharacteristic vector, r_tIt it is residual vector.

7) when utilizing MPCA to set up monitoring model, it is necessary to first determining that two control limit, control to limit the use of to judge whether current data is in normal operating condition, the two controls limit and is called T²Statistic controls limit and SPE statistic controls limit, calculates T²The control limit of statistic meets F-distribution formula, and such as formula (10), wherein �� is main constituent number, and I is batch number of modeling, and �� is significance level.

$T_{t,\mathrm{\α}}^2:\frac{\mathrm{\β}(I-1)}{I-\mathrm{\β}}{F}_{\mathrm{\β},I-\mathrm{\β},a}---\left(10\right)$

The control limit of SPE statistic is calculated by formula (11), wherein mean_t, var_tRepresent that each sub-period batch measurement data in modeling data is in the average of t SPE and variance respectively.

$\left\{\begin{array}{c}{\mathrm{SPE}}_{t\mathrm{\α}}~\frac{{\mathrm{var}}_t}{2{\mathrm{mean}}_t}{\mathrm{\χ}}_{t,h_t,\mathrm{\α}}^2\\ h_t=\frac{2{\mathrm{mean}}_t^2}{{\mathrm{var}}_t}\end{array}\right.---\left(11\right)$

8) in each sub-period, adopt 1) and 2) in mode, take each I of data under dissimilar failure operation state respectively_faultGroup, utilizes 3) methodological standardization processes, and utilizes 4) in method extract the main constituent Q of fault data of each type₁,Q₂,��,Q_fault;

9) the main constituent Q that will extract₁,Q₂,��,Q_faultAs the input of supporting vector machine model, in each sub-period, set up a fault diagnosis model, finally can obtain G model, be SVM respectively₁, SVM₂..., SVM_G, wherein each SVM model adopts man-to-man multi-categorizer make.

9.1) any one SVM_g(g=1,2 ..., G) model is all a decision hyperplane.

9.2) man-to-man multi-categorizer building method, inhomogeneity sample combination of two, for u class training sample, is used for training by the methodIndividual two graders. Then test sample is substituted into all graders, adopt ballot method to carry out decision-making, if grader SVM_1,2Thinking that sample belongs to the 1st class, then add 1 in the ballot of the 1st class, otherwise add 1 in the ballot of the 2nd class, finally add up all kinds of polls, who gets the most votes's class is test sample generic.

B. the on-line fault diagnosis stage

10) current time data X is gathered_t(1 �� J), and judge S' during the operator residing for current time_g(g=1,2 ..., G) section;

11) X is calculated according to formula (9)_tScore vector Q corresponding to (1 �� J)_t(1 �� ��) and forecast error vector R_t(1 �� J), in formula, E is the unit matrix of (J �� J);

$\left\{\begin{array}{c}{Q}_t=X_t{\stackrel{\‾}{P}}_{S_g^{\′}}\\ R_t=X(E-{\stackrel{\‾}{P}}_{S_g^{\′}}{\left({\stackrel{\‾}{P}}_{S_g^{\′}}\right)}^T)\end{array}\right.---\left(12\right)$

12) X is calculated_t(1 �� J) correspondenceAnd SPE_tStatistic, computing formula such as formula (10);

$\left\{\begin{array}{c}{T}^2=Q_t{\left(S_g^{\′}\right)}^{-1}{Q_t}^T\\ SPE=R_t{R_t}^T\end{array}\right.---\left(13\right)$

Wherein, R_tIt is matrix X_tThe forecast error vector of (1 �� J).

13) judge nowOr SPE_tStatistic is no beyond controlling limit T²Or SPE, if both of which without departing from, then active procedure measurement data is normal, returns 1); Otherwise current operation process has fault, carries out 14);

14) at sub-period S'_gIn, the data that new lot is unknown are carried out data filling, filling data is current time X_t(1 �� J), finally can obtain the data matrix of this period is

15) to data matrixStandardization also extracts the main constituent of this data matrix, is updated to the SVM of this sub-period afterwards_g(g=1,2 ..., G) in model, carry out fault diagnosis.

Beneficial effect

Compared with prior art, the present invention utilizes MPCA to divide batch process, sets up fault detection and diagnosis model at each sub-period, reduces the search volume of fault diagnosis. When using the method on-line fault diagnosis, it is not necessary to fill the unknown data of whole process, and be the unknown data filling the period of breaking down, so can reduce the impact that too much filling unknown data brings, thus the accuracy rate of fault diagnosis can be improved. Simultaneously without updating fault diagnosis model frequently, and only need to select the diagnostic cast of period of breaking down, reduce the complexity of fault diagnosis.

Accompanying drawing explanation

Fig. 1 is the off-line modeling flow chart of the present invention;

Fig. 2 is on-line fault diagnosis flow chart;

Fig. 3 is normal data X (30 �� 10 �� 400) composition form schematic diagram;

Fig. 4 sub-period division result;

The MPCA T to fault batch in Fig. 5 sub-period 5²With SPE monitoring figure;

Fig. 6 design sketch that test data are not classified by SVM at times;

Fig. 7 present invention design sketch to test data classification;

Detailed description of the invention

The penicillin emulation platform PenSim2.0 of American I llinois state Institute of Technology process monitoring and technology group development, monitoring and fault diagnosis for batch process provide a standard platform, and this platform has become more influential penicillin emulation platform in the world.

The present invention is with this platform for simulation object, and the response time arranging penicillin fermentation each batch is 400h, and the sampling interval is 1 hour, chooses 10 process variables and carries out simulation study, as shown in table 1. Meanwhile, this platform can set that three kinds of faults: 1. air mass flow, 2. power of agitator, 3. bottoms stream rate of acceleration. The type of fault disturbance has step disturbance and slope disturbance two kinds, can set the amplitude of two kinds of disturbances, introducing time of disturbance and termination time further.

Model variable set up by table 1

30 batches of normal data X (30 �� 10 �� 400) of this experiment simulation are as dividing the period and setting up the data base of MPCA malfunction monitoring model, and wherein the initial condition of every batch is slightly changed, and all measurands is all added random measurement noise simultaneously.

Additionally, also simulate fault 1 (variable 1 both phase step fault data), fault 2 (variable 2 both phase step fault data) and each 30 batches of fault 3 (variable 3 both phase step fault data) to be used for training SVM fault diagnosis model, wherein each fault selects 15 batches to be used for setting up fault diagnosis model, and remaining 15 batches are used for testing this diagnostic cast. In order to make each fault data comprise more fully fault message, 15 batches of training datas of same fault are introduce disturbance in different response time and be all extended to reaction and terminate respectively.

The inventive method is applied to above-mentioned penicillin fermentation process simulation object and includes off-line modeling and the big step of inline diagnosis two, be specifically presented below:

A. the off-line modeling stage:

Step 1: above-mentioned 30 batches of normal data X (30 �� 10 �� 400) being launched along batch direction, concrete form is shown in Fig. 2. It figure is the matrix of 30 row 4000 row.

Step 2: be standardized X (30 �� 10 �� 400) processing.

First formula is pressedThe jth process variable of calculating kth sampling instant average on all batches, wherein x_i,k,jFor X_{30��10��400}The measured value of the jth process variable of kth sampling instant in i-th batch, k=1 ..., 400, j=1 ..., 10; The standard variance s of the jth process variable of kth sampling instant_k,jComputing formula be, $s_{k,j}=\sqrt{\frac{1}{30-1}\underset{i=1}{\overset{30}{\mathrm{\Σ}}}{(x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j})}^2},$ K=1 ..., 400, j=1 ..., 10.

Afterwards X (30 �� 10 �� 400) being standardized, wherein the standardized calculation formula of the jth process variable of kth sampling instant is as follows in i-th batch:

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}$

Wherein, i=1 ..., 30, j=1 ..., 10, k=1 ..., 400;

Step 3:X (30 �� 10 �� 400) is by step 2) to obtain two-dimensional matrix after standardization be X'(30 �� (10 �� 400));

Step 4: by X'(30 �� (10 �� 400)) 400 timeslice matrix X it are divided into along moment direction_k' (30 �� 10), k=1 ..., 400.

Step 5: utilize principal component analysis function (princomp) in MATLAB to extract two-dimensional matrix X_k' the main constituent Q of (30 �� 10)_k��

Step 6: work as Q_��=Q_��(��=��+1, �� < �� and ��=2 ..., 400) time, �� and �� can be divided into the same period, if Period Length L does not reachDuring sampled point, it is possible to this segment data is grouped in relatively small one section of the main constituent number difference that left and right is adjacent. If Period Length L less thanSampled point and equal with adjacent main constituent number difference time, then take the main constituent number that main constituent number is this period of former and later two periods of this period respectively and calculate its contribution rate, in this period be grouped into that contribution rate change is relatively small one period. Finally mark off 4 thick period S_f(f=1 ..., 4).

Step 7: at period S_fIn, utilize formula (7), when calculatingDuring more than 0.8, just b and the c sampling time was put in the period that same refinement divides. Batch process is divided into 7 sub-period S the most at last_g(g=1, ..., 7), the respectively period 1 (1-37h), the period 2 (38-85h), the period 3 (86-174h), period 4 (175-215h), period 5 (216-296h), the period 6 (297-377h), the period 7 (378-400h).

Step 8: off-line modeling stage by stage

This part is for the period 5, and MPCA monitoring model and SVM fault diagnosis model to sub-period are set up and be described.

Step 8.1: take out X'(30 �� (10 �� 400)) in the 5th period corresponding data set up the MPCA model of period 5, and calculate the average load matrix of period 5

Step 8.2: utilize formula 10) and 11) determine the T of period 5²Control limit with SPE statistic, control limit figure such as Fig. 5.

Fault 1,2 and 3 data in period 5 are carried out pretreatment, obtain matrix X after standardization by step 8.3: the method utilizing step 1, step 2, step 3 and step 4_5f1,X_5f2And X_5f3��

Step 8.4: utilize the principal component analysis function (princomp) in Matlab to X_5f1,X_5f2And X_5f3Decompose, main constituent Q can be obtained_5f1,Q_5f2And Q_5f3��

Step 8.5: utilize the svmtrain function in libsvm workbox in Matlab to Q_5f1,Q_5f2And Q_5f3It is trained, decision hyperplane F may finally be obtained₅(), this decision hyperplane is exactly fault diagnosis model SVM₅��

In like manner, according to the method for step 8.2 to 8.5, set up the fault diagnosis model of other periods, finally give model SVM₁, SVM₂, SVM₃, SVM₄, SVM₆, SVM₇��

B. the inline diagnosis stage:

This part introduces the both phase step fault of 5% for test data with variable 1 220 moment.

Step 9: gather the data x of 10 process variables of current sweat t sampling instant_t, and it is standardized obtaining by the average and standard variance according to the t obtained in step 2The wherein jth process variable of t sampling instantStandardization formula as follows:

${\tilde{x}}_{t,j}=\frac{x_{t,j}-{\stackrel{\&OverBar;}{x}}_{t,j}}{s_{t,j}}$

Wherein, x_t,jFor the jth process variable in current the gathered Fermentation Data of t sampling instant,It is the meansigma methods of the jth process variable of t sampling instant, s_t,jIt is the standard variance of the jth process variable of t sampling instant, j=1 ..., 10;

Step 10: according to formula 15), the t after normalized gathers dataCorresponding principal component vector Q_5tWith forecast error vector R_5t��

Step 11: calculate t according to the method in step 8.3And SPE_tStatistic.

Step 12: by above-mentioned calculated monitoring statisticAnd SPE_tThe control limit determined with modeling procedure 8.3 compares, if without transfiniting, is normal, then gathers the data in next moment and return step 9, and 10,11 and 12, until sweat is complete; Otherwise it is assumed that break down.

Step 13: if there being fault, then first determine whether that the stage that fault occurs was the 5th period, and transfer out the fault diagnosis model SVM of this period₅;

Step 14: within the 5th period, utilizes current time sampled data x_tUnknown data filled full, then the data matrix obtaining this period is X₅;

Step 15: propose principal component method to X according to PCA in the data normalization method and steps 5 in step 2₅Carry out data prediction, main constituent Q may finally be obtained₅��

Step 16: in Matlab, svmpredict function utilizes the model model parameter of svmtrain training in step 8.5, calculates Q₅Fault data be judged to other label be 1 and differentiate accuracy rate accurate be 86.67%.

Above-mentioned steps is the inventive method concrete application in penicillin fermentation emulation platform fault detection and diagnosis field. In order to verify effectiveness of the invention, three kinds of each 15 batch datas of fault are tested. Fig. 5 is that the failure monitoring figure in the 5th period occurs fault, and failure monitoring figure and Fig. 5 of other periods is similar no longer to be illustrated here. Line parallel with abscissa in figure represents that curve is monitor in real time value by the control limit that Density Estimator method is determined. If curve has fault to occur higher than controlling limit explanation, it it is otherwise normal operating condition.

Fig. 6 and Fig. 7 represents that fault occurs when sub-period 5 respectively, the fault diagnosis design sketch under not grading method and the inventive method, and in figure, circle represents the classification of actual test set, and asterisk represents prediction test set classification. Under same test set sample, when two kinds of marks coincide, then illustrate that test result is identical with actual result, otherwise test result mistake.

In order to significantly more efficient proof the inventive method is applied to the superiority in sweat fault diagnosis, other periods have been carried out again fault diagnosis research by the present invention respectively, and diagnostic result list contrast is as follows:

Each period fault diagnosis result of table 2

From upper table it is seen that, use the inventive method can be greatly improved the accuracy rate of fault when fault occurred in front several period, but along with the continuation produced, though the method is better than traditional method diagnosis effect, but the amplitude that accuracy rate improves is really limited. As a whole, the inventive method can significantly improve fault accuracy rate.

Claims (1)

1. based on a batch process on-line fault diagnosis method of sub-period MPCA-SVM, including " off-line modeling " and " inline diagnosis " two stages, it is characterised in that specifically comprise the following steps that

A. the off-line modeling stage

1) gathering the historical data under sweat nominal situation, the I lot data under the same sweat same process that described historical data X is obtained by off-line test is constituted, X=(X₁,X₂,...,X_I)^T, wherein X_i(i=1,2 ..., I) represent the i-th lot data; Each batch comprises K sampling instant, i.e. X_i=(X_i,1,X_i,2,...,X_i,K), wherein X_i,kRepresent the data that i-th batch of kth sampling instant gathers; Each sampling instant gathers J process variable, i.e. X_i,k=(x_i,k,1,x_i,k,2,...,x_i,k,J), wherein x_i,k,jThe measured value of the jth process variable of kth sampling instant in representing i-th batch;

2) being standardized historical data X processing, processing mode is as follows:

Average and the standard variance of all process variables, the wherein average of the jth process variable of kth sampling instant is engraved when first calculating historical data X allComputing formula beWherein x_i,k,jThe measured value of the jth process variable of kth sampling instant in representing i-th batch, k=1 ..., K, j=1 ..., J; The standard variance s of the jth process variable of kth sampling instant_k,jComputing formula be, $s_{k,j}=\sqrt{\frac{1}{I-1}\underset{i=1}{\overset{I}{\&Sigma;}}{(x_{i,k,j}-{\stackrel{\&OverBar;}{x}}_{k,j})}^2},k=1,...,K,j=1,...,J;$

Then historical data X being standardized, wherein the standardized calculation formula of the jth process variable of kth sampling instant is as follows in i-th batch:

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\stackrel{\&OverBar;}{x}}_{k,j}}{s_{k,j}}---\left(1\right)$

Wherein, i=1 ..., I, j=1 ..., J, k=1 ..., K;

3) X is by step 2) standardization obtain new two-dimensional matrix X', this matrix has (K �� J) individual column vector, i.e. X'=(X '₁,X��₂,...,X'_K��J), wherein jth column vector X'_j=(X'_j,1,...,X'_j,K)^T, X'_j,k=(X'_j,k,1,...,X'_j,k,I)^T, wherein X'_j,k,iRepresent through step 2) value corresponding in i-th batch of jth process variable kth sampling instant after standardization, wherein i=1 ..., I, j=1 ..., J, k=1 ..., K;

4) multidirectional pivot analysis MPCA method is utilized to extract the main constituent in each moment in X'; If extracting kth sampling instant data X'_kMain constituent, concrete step is as follows:

4.1) two-dimensional matrix X' is obtained_kCovariance matrix COV;

$COV=\frac{1}{I-1}{X}_k^{\&prime;T}{X}_k^{\&prime;}---\left(2\right)$

4.2) Eigenvalues Decomposition to Matrix C OV;

COV=V �� V^T(3)

Formula (3) discloses the incidence relation of covariance matrix COV, wherein ��=(��₁,��₂,��,��_v) for diagonal matrix, v is two-dimensional matrix X'_kEigenvalue number, �� comprises the non-negative factual investigation (�� that m amplitude is successively decreased₁�ݦ�₂�ݡ��ݦ�_m>=0); V is orthogonal matrix (V^��V=E, E are unit battle array);

In definition, by MPCA method, matrix X'_kIt is equivalent to X'_k=QP^��+ R, and be updated to formula (3) and obtain:

$COV=\frac{1}{I-1}{\mathrm{PQ}}^{T}QPQ---\left(4\right)$

Wherein, Q is score matrix, and P is load matrix, and R is residual matrix;

Corresponding (3) and (4) are every:

P=V

$\mathrm{\&Lambda;}=\frac{1}{I-1}{Q}^{T}Q---\left(5\right)$

If the sum equation more than one threshold value 0.85 of front A main constituent, then before just extracting, A pivot is used as aggregative indicator, then original J dimension space just becomes A dimension and A��J;

4.3) score matrix Q is obtained;

In order to optimally obtain the variable quantity of data, minimize random noise to the PCA impact produced simultaneously, retain the load vector corresponding with A eigenvalue of maximum, then X'_kProjection information at lower dimensional space is included in score matrix:

Q=X'_kP(6)

5) utilize MPCA that batch process is carried out Time segments division:

5.1) period slightly divides:

When neighbouring sample point has identical main constituent number, just these sampled points are divided into the same period, if have the Period Length L of identical main constituent number less than whole sweat 1/10 time, in this period data is grouped into that the adjacent main constituent number difference in left and right is relatively small one section; If this Period Length is less than the 1/10 of whole sweat and when the main constituent number of this period is equal with the main constituent number difference of adjacent two periods, then take main constituent number that front and back period main constituent number is this period respectively and calculate its contribution rate, this section is grouped in the period that contribution rate change is relatively small;

Final thick division obtains F period, is expressed as S₁,S₂,��,S_F;

5.2) period carefully divides: at the arbitrarily thick period S divided_f(f=1,2 ..., F) in, utilize between load matrix angle information and range information to define similarity measurement formula (7);

$d_{b,c}^{S_f}=\underset{l=1}{\overset{a_f}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_l\frac{\left|\right|P_{bl}^T{P}_{cl}\left|\right|}{\left|\right|P_{bl}\left|\right|\&CenterDot;\left|\right|P_{cl}\left|\right|}+(1-{\mathrm{\&gamma;}}_l)e^{-\left|\right|P_{bl}-P_{cl}\left|\right|}---\left(7\right)$

${\mathrm{\&gamma;}}_l=\frac{1}{l}/\underset{h=1}{\overset{a_f}{\mathrm{\&Sigma;}}}\frac{1}{h},l=h=1,2,...,a_f---\left(8\right)$

Wherein b and c represents at period S_fTwo sampling times of interior arbitrary neighborhood, a_fFor period S_fInterior sampled point number, l and h represents that sampled point 1 arrives a_fIn any one value, | | P_bl| | with | | P_cl| | it is the modulus value of the load matrix of b and c sampling instant respectively,Represent the inner product of b and c sampling instant load matrix, r_lFor weight coefficient in order to emphasize the different importances of different projecting direction,Represent a in two timeslice load matrix_fThe weighted sum of the included angle cosine value of individual projecting direction,Represent at period S_fThe similarity of interior b and c sampling instant load matrix;

G period is got in final refinement, is expressed as S '₁,S��₂,...,S'_G, and calculate any one sub-period S'_g(g=1,2 ..., G) average load matrix

6) at ready-portioned S '₁,S'₂,...,S'_GIn set up each sub-period MPCA monitoring model, such as formula (9);

Being respectively provided with close load matrix for each sub-period, therefore the MPCA model in each sub-period all adopts average load matrixDescribe; For anyon period S'_gIn the MPCA model of any t be expressed as:

$\left\{\begin{array}{c}{q}_t=S_g^{\&prime;}{\stackrel{\&OverBar;}{P}}_{S_g^{\&prime;}}\\ \widehat{S_g^{\&prime;}}=q_t{\left({\stackrel{\&OverBar;}{P}}_{S_g^{\&prime;}}\right)}^T\\ r_t=S_g^{\&prime;}-\widehat{S_g^{\&prime;}}\end{array}\right.---\left(9\right)$

Wherein q_tIt is sub-period S'_gIn the load vector of t,It is sub-period S'_gCharacteristic vector, r_tIt it is residual vector;

7) when utilizing MPCA to set up monitoring model, it is necessary to first determining that two control limit, control to limit the use of to judge whether current data is in normal operating condition, the two controls limit and is called T²Statistic controls limit and SPE statistic controls limit, calculates T²The control limit of statistic meets F-distribution formula, and such as formula (10), wherein �� is main constituent number, and I is batch number of modeling, and �� is significance level;

$T_{t,\mathrm{\&alpha;}}^2:\frac{\mathrm{\&beta;}(I-1)}{I-\mathrm{\&beta;}}{F}_{\mathrm{\&beta;},I-\mathrm{\&beta;},a}---\left(10\right)$

The control limit of SPE statistic is calculated by formula (11), wherein mean_t, var_tRepresent that each sub-period batch measurement data in modeling data is in the average of t SPE and variance respectively;

$\left\{\begin{array}{c}{\mathrm{SPE}}_{t\mathrm{\&alpha;}}~\frac{{\mathrm{var}}_t}{2{\mathrm{mean}}_t}{\mathrm{\&chi;}}_{t,h_t,\mathrm{\&alpha;}}^2\\ h_t=\frac{2{\mathrm{mean}}_t^2}{{\mathrm{var}}_t}\end{array}\right.---\left(11\right)$

8) in each sub-period, adopt 1) and 2) in mode, take each I of data under dissimilar failure operation state respectively_faultGroup, utilizes 3) methodological standardization processes, and utilizes 4) in method extract the main constituent Q of fault data of each type₁,Q₂,��,Q_fault;

9) the main constituent Q that will extract₁,Q₂,��,Q_faultAs the input of supporting vector machine model, in each sub-period, set up a fault diagnosis model, finally obtain G model, be SVM respectively₁, SVM₂..., SVM_G, wherein each SVM model adopts man-to-man multi-categorizer make;

9.1) any one SVM_g(g=1,2 ..., G) model is all a decision hyperplane;

9.2) man-to-man multi-categorizer building method, inhomogeneity sample combination of two, for u class training sample, is used for training by the methodIndividual two graders; Then test sample is substituted into all graders, adopt ballot method to carry out decision-making, if grader SVM_1,2Thinking that sample belongs to the 1st class, then add 1 in the ballot of the 1st class, otherwise add 1 in the ballot of the 2nd class, finally add up all kinds of polls, who gets the most votes's class is test sample generic;

B. the on-line fault diagnosis stage

10) current time data X is gathered_t(1 �� J), and judge S' during the operator residing for current time_gSection;

11) X is calculated according to formula (9)_tScore vector Q corresponding to (1 �� J)_t(1 �� ��) and forecast error vector R_t(1 �� J), in formula, E is the unit matrix of (J �� J);

$\left\{\begin{array}{c}{Q}_t=X_t{\stackrel{\&OverBar;}{P}}_{S_g^{\&prime;}}\\ R_t=X\left(E-{\stackrel{\&OverBar;}{P}}_{S_g^{\&prime;}}{\left({\stackrel{\&OverBar;}{P}}_{S_g^{\&prime;}}\right)}^T\right)\end{array}\right.---\left(12\right)$

12) X is calculated_t(1 �� J) correspondenceAnd SPE_tStatistic, computing formula such as formula (13);

$\left\{\begin{array}{c}{T_t}^2=Q_t{\left(S_g^{\&prime;}\right)}^{-1}{Q_t}^T\\ {\mathrm{SPE}}_t=R_t{R_t}^T\end{array}\right.---\left(13\right)$

Wherein, R_tIt is matrix X_tThe forecast error vector of (1 �� J);

13) judge nowOr SPE_tWhether statistic is beyond controlling limit T²Or SPE, if both of which without departing from, then active procedure measurement data is normal, returns 1); Otherwise current operation process has fault, carries out 14);

14) at sub-period S'_gIn, the data that new lot is unknown are carried out data filling, filling data is current time sampled data X_t(1 �� J), finally obtaining the data matrix of this period is

15) to data matrixStandardization also extracts the main constituent of this data matrix, is updated to the SVM of this sub-period afterwards_g(g=1,2 ..., G) in model, carry out fault diagnosis.