CN105607631A - Batch process weak fault model control limit establishment method and weak fault monitoring method - Google Patents

Abstract

The invention discloses a batch process weak fault monitoring model establishment method based on tendency turning points, comprising the following steps: 1, acquiring I batches of normal operating data; 2, sectioning the batch process; 3, extracting turning points by adopting a method based on the tendency turning point, and ranking the turning points according to the orthogonal distance; 4, checking the turning points extracted with variables among the batches; 5, aggregating the tendency turning points of all variables in the batches and verifying the truth or falsity; 6, constructing a time sequence augmented matrix by using a time lag expanded analysis data array, and performing PCA (Principal Component Analysis) on the time sequence augmented matrix; 7, clustering all established models according to the similarity of the load direction at each phase; and 8, calculating T2 and SPE (Synergistic Processing Element) statistics and a control limit of each type. The invention further provides a batch process weak fault monitoring method based on tendency turning points.

Description

The weak fault model control limit method for building up of batch process and weak fault monitoring method

Technical field

The present invention relates to the weak fault monitoring method of a kind of batch process, especially a kind of based on trend reverse pointThe weak fault model of batch process set up and monitoring method.

Background technology

Along with developing rapidly of modern industry and science and technology, the structure of modern comfort becomes increasingly complex, fromThe scale of movingization system is also increasing, and this has increased the possibility that system breaks down greatly. Once itsBreak down, not only cause shutting down stopping production, also can cause the massive losses of personnel and property, cause disasterSexual behavior part. In order to improve security, reliability and the validity of industrial processes and control system, withIn time, improves the quality of products, reduces costs, in the urgent need to process of establishing control monitoring system. Process monitoring systemThe task of system is exactly the running status of work under supervision process, constantly variation and the fault message of testing process,To prevent the generation of catastrophic failure, reduce the fluctuation of product quality, for the economy that improves enterprise simultaneouslyThe competitiveness tool of benefit and raising enterprise is of great significance.

Multivariate statistical process monitoring technology has obtained many successful Application, Qi Zhongduo in industrial process at intermittenceExpand and improve with multidirectional PCA (MPCA) and multidirectional PLS (MPLS).But it is not the method based on MPCA is an overall modeling method, very responsive to the identification of Weak fault.After propose again that way moving window MPCA method solves that operating condition changes and slow time-varying characteristics with adaptiveAnswer the more model data between new lot of ground; Lee has proposed online updating MPCA method; Lu etc. have proposedA kind of two-dimentional dynamic principal component analysis 2D-DPCA strategy, by choosing suitable Data support territory for rightAnalytic unit carries out time and batch two-way expansion, can directly extract batch and time orientation on localDynamical correlation relation. But the analyst coverage of dependency relation and extract performance and be limited to that " data are propped up between batchHold territory " concrete setting. Although solve to a certain extent the malfunction monitoring problem of slow time-varying process,But the method is only for linear change process. Recurrence based on kernel function has been proposed again in recent yearsKPCA (RKPCA) algorithm becomes the adaptive process monitoring of non-linear process while being used for, but all main pins of these methodsTo the slow variation situation of systematic parameter.

In actual intermittently industrial process, in the time that the fluctuation of production primary condition is larger, poor between lot dataDifferent more obvious, then adopt said method to carry out modeling to process, None-identified is occurred in such casesWeak fault, has descended because weak fault has been submerged in the fluctuation of primary condition.

Therefore, a kind of can monitoring the serious weak fault of growth process that do not wait, and can be greatlyA little less than the batch process based on trend reverse point of the computation complexity of reduction algorithm, fault monitoring method becomes solutionThe certainly key of problem.

Summary of the invention

The present invention has designed and developed the weak fault model control limit of a kind of batch process method for building up, has overcome existingThere is algorithm complexity in technology, be difficult to carry out not waiting the defect of the weak fault detect of growth process, solved initialThe weak malfunction monitoring problem of conditional fluctuation under more greatly.

Another object of the present invention is to provide the weak fault monitoring method of a kind of batch process, has greatly improvedThe speed of on-line monitoring and efficiency.

Technical scheme provided by the invention is:

A little less than batch process based on trend reverse point, a malfunction monitoring method for establishing model, is characterized in that,Comprise the steps:

The data of step 1, I batch of normal operation of collectionWherein i=1 ..., I, K_iIBatch sampling number, J be gather variable number;

Step 2, adopt CPV1 discrete method to carry out segmentation each batch process lot data

Step 3, the method for employing based on trend reverse point are carried out turning point extraction, and are entered by orthogonal distanceLine ordering;

Step 4, by batch between the turning point that extracts of each variable carry out parallelism inspection, when the trend of extracting turnsWhen break both sides variation tendency has obvious difference, need to replace this trend reverse point, again extract;

Step 5, by batch in all variablees the polymerization of trend reverse point and carry out true and false verification;

Step 6, set up single-MDPCA model for each lot data based on trend reverse point,Use time lag data extending to analyze data matrix structure time series augmented matrix, and to time series augmentation squareBattle array is carried out PCA processing;

Step 7, at every one-phase, by set up all single-MDPCA models press loading direction phaseCarry out cluster like property;

Step 8, will belong to of a sort lot data and launch to re-establish model by variable, calculate T²WithThe control limit of SPE statistic and each class.

Preferably, between step 1 and step 2, also comprise: each lot data is adopted by variableThe threshold values denoising method denoising of Donoho.

Preferably, in step 3, turning point extracting method is: for time series, first by starting pointSave as turning point with terminal; Then between starting point and terminal, search by maximum orthogonality distance methodNext trend reverse point; In the new turning point both sides of finding, again find out respectively and there is maximum orthogonality distanceFrom turning point, and by distance compare, the point with maximum orthogonality distance is extracted as the next onePoint.

Preferably, in step 5, the process of true and false verification is: when the trend reverse point extracting is exceedingWhen having repetition in 2/3 total batch of number or approaching very much, these points are rejected, only leave one of them point;When duplicating or approach very much in the trend reverse point extracting being less than 1/3 total batch of number, repetitionPoint replace again and search.

Preferably, in step 6, time series augmented matrix is

$X_d=\left[\begin{array}{cccc}{x}^T\left(k\right)& x^T(k-1)& L& x^T(k-d)\\ x^T(k-1)& x^T(k-2)& L& x^T(k-d-1)\\ M& M& M& M\\ x^T(k+d-K)& x^T(k+d-K-1)& L& x^T(k-K)\end{array}\right]$

In formula, x (k)=[x_1,kx_2,kLx_J,k]^TThe observation vector of J variable at sampling instant k place, at kThe turning point vector of J variable of moment extraction, K is the sampled point of each batch of modeling data matrixNumber, d is time lag length.

Preferably, in step 7, all single-MDPCA models are undertaken by loading direction similitudeThe method of cluster comprises:

Adopt cosine angle between two model projection vectors to measure the similarity degree of two models:

$\begin{array}{c}\cos (P_i^m,P_k^m)=\frac{<P_i^m,P_k^m>}{\sqrt{\underset{m=1}{\overset{p_{\max }}{\&Sigma;}}{\left(P_i^m\right)}^2}\sqrt{\underset{m=1}{\overset{p_{\max }}{\&Sigma;}}{\left(P_k^m\right)}^2}}\\ i,k=1,2,L,I;m=1,2,L,p_{\max }.\end{array}$

In formula, P represents load matrix,Represent P_iM row,Represent P_kM row, p_maxForMaximum pivot number;

Introduce weight coefficient ω_m：

${\mathrm{\&omega;}}_m=\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_m}/\underset{m=1}{\overset{p_{\max }}{\&Sigma;}}\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_m}$

Wherein, λ is nonzero eigenvalue;

Can calculate two angles between model projection vector:

${\mathrm{\&theta;}}_{i,k}=\underset{m=1}{\overset{p_{\max }}{\&Sigma;}}{\mathrm{\&omega;}}_m{\mathrm{\&theta;}}_m;$

Finally adopt k-means algorithm cluster that the load of weighting is divided into different classes.

Preferably, in step 8,

T²Statistic is:

$T^2=[t_1,L,t_{p_{\max }}]{\mathrm{\&Lambda;}}^{-1}{[t_1,L,t_{p_{\max }}]}^T$

Here t is score vector, Λ^-1P before representing_maxThe diagonal matrix of individual characteristic vector institute character pair value compositionContrary;

T²The control limit of statistic can distribute and draw according to F:

$T^2~\frac{p_{max}(K-1)}{K-p_{max}}{F}_{p_{max},K-p_{max},\mathrm{\&alpha;}}$

In formula,Be to be α corresponding to insolation level, the free degree is p_max，K-p_maxF under conditionThe critical value distributing.

Preferably, in step 8,

SPE statistic is:

$SPE=\left|\right|\mathrm{\&Phi;}\left(x\right)-{\widehat{\mathrm{\&Phi;}}}_a\left(x\right)|{|}^2=\underset{k=1}{\overset{K}{\&Sigma;}}{t_k}^2-\underset{k=1}{\overset{p_{max}}{\&Sigma;}}{t_k}^2$

The control limit of SPE statistic can be according to x²Distribution draws:

${\mathrm{SPE}}_{\mathrm{\&alpha;}}~{\mathrm{g\&chi;}}_h^2;g=\frac{u}{2v},h=\frac{2v^2}{u}$

In formula, u and v are average statistical and the variances of modeling data SPE statistics value.

A fault monitoring method a little less than batch process based on trend reverse point, comprises the steps:

Step 1, Real-time Collection process data;

Step 2, in 10% stage length range, new lot data are monitored with all class models,And calculate overall target

SI＝SPE/CONs＋T²/CON_T，

Wherein, CON_sAnd CON_TBe respectively SPE and T²Control limit; By the model of the minimum SI of correspondenceAs best monitoring model;

The best monitoring model of choosing in step 3, use step 2 is monitored process data, works as numberAccording to SPE value and T²Exceed and controlled limit, in process, fault has occurred.

The invention has the beneficial effects as follows: weak fault monitoring method of the present invention is compared with existing weak fault methodAdvantage be: 1,, for the seriously not isometric characteristic of batch process, solved isometricization of lot data and askedTopic; 2, reduce computation complexity, improved speed and the efficiency of on-line monitoring; 3, successfully solvedWeak malfunction monitoring problem under primary condition fluctuation is large; 4, be successfully applied to penicillin fermentation processMalfunction monitoring; 5, the present invention is also applicable to the intermittently event of generative process such as pharmacy, semiconductor machining, chemical industryIn barrier monitoring.

Brief description of the drawings

Fig. 1 is normal modeling data schematic diagram.

Fig. 2 is orthogonal distance schematic diagram.

Fig. 3 is for extracting trend reverse point schematic diagram.

Fig. 4 is that power of agitator data are extracted trend reverse point schematic diagram.

Fig. 5 extracts trend reverse point schematic diagram for producing thermal data.

Fig. 6 is that dissolved oxygen concentration data are extracted trend reverse point schematic diagram.

Fig. 7 is that dense carbon dioxide degrees of data is extracted trend reverse point schematic diagram.

Fig. 8 is that temperature data extracts trend reverse point schematic diagram.

Fig. 9 is monitoring system structural representation of the present invention.

Figure 10 is historical data T²Control limit monitoring result schematic diagram.

Figure 11 is that historical data SPE controls limit monitoring result schematic diagram.

Figure 12 is weak On-line Fault monitoring T²Result schematic diagram.

Figure 13 is weak On-line Fault monitoring SPE result schematic diagram.

Detailed description of the invention

Below in conjunction with accompanying drawing, the present invention is described in further detail, to make those skilled in the art's referenceDescription word can be implemented according to this.

As shown in Figure 1, the invention provides the foundation side of the weak fault model occurring in a kind of observation processMethod, comprises the steps:

Data under step 1, collection nominal situation.

Under nominal situation condition, gather the data of I batch of normal operationi＝1,…,I，K_iRepresent the sampling number of i batch, J is the variable number gathering, as shown in Figure 1. I batchEach variable all gathers K_iIndividual sampled point. As front four batches of sampling numbers are respectively 400,420,420,460, K₁＝400，K₂＝420，K₃＝420，K₄＝460。

Step 2: adopt the threshold values denoising method of Donoho to carry out denoising by variable each lot data.

Step 3: by lot data segmentation.

Adopt CPV1 discrete method to carry out segmentation each batch process lot data. For each batchData, first form the data segment corresponding to this moment by certain sampling instant to the data of last sampling instant;Then each data segment is carried out to pivot analysis, obtain the first pivot contribution rate; Finally, adopt the first masterThe variation of unit's contribution rate characterizes the transformation in batch inner stage, and batch process is divided to the stage.

Step 4: adopt the method based on trend reverse point (TTP) to carry out turning point extraction.

As shown in Figure 2, the extracting method based on TTP can be passed through to calculate: a some x₂＝(v₂,t₂) to its twoEnd points x₁＝(v₁,t₁) and x₃＝(v₃,t₃) orthogonal distance D₀As shown in the formula:

$D_o=\sqrt{{(t_c-t_2)}^2+{(v_c-v_2)}^2}$

t_c＝(t₂+s×v₂＋s×v₃-s²×t₃)/(1＋s²)

v_c＝s×t_c-s×t₃＋v₃

$s=\frac{v_3-v_1}{t_3-t_1}.$

As shown in Figure 3, using the point with maximum orthogonality distance as trend reverse point (TTP). By TTPSort by orthogonal distance, adopt priority query's method to guarantee that each point extracting is all in remaining pointThere is the turning point of maximum orthogonality distance. Specific as follows:

For time series, first starting point and terminal are saved as turning point; Then in starting point andBetween terminal, search next trend reverse point by maximum orthogonality distance method; In the new turning point both sides of finding,Again find out respectively the turning point with maximum orthogonality distance, and compare by distance, will there is maximumThe point of orthogonal distance extracts a little as the next one; By that analogy, remaining point is all remained at every turn,With the distance-taxis of pressing of new extraction, to find out next trend reverse point.

For example a certain lot data gathers 8 temperature spot (T in chronological order₁,t₁)-(T₈,t₈), first by starting point(T₁,t₁) and terminal (T₈,t₈) as turning point, 6 point (T in the middle of then calculating successively₂,t₂)—(T₇,t₇) arriveThe orthogonal distance of two-end-point, and will be apart from sorting, as point (T₅,t₅) the distance maximum that obtains, will(T₅,t₅) as next turning point; At point (T₅,t₅) both sides, obtain two interval (T₁,t₁)—(T₅,t₅) and(T₅,t₅)—(T₈,t₈), and using these four points respectively as two interval starting points and terminal, in computation intervalEach point to the orthogonal distance of two-end-point separately, and sequence, finds out the point of maximum orthogonality distance as nextIndividual turning point, as point (T₃,t₃); Process has been divided into (T₁,t₁)—(T₃,t₃)、(T₃,t₃)—(T₅,t₅) and (T₅,t₅)—(T₈,t₈) three intervals, in the same way, then calculate respectively interval separately in each point to two-end-pointOrthogonal distance. Due to point (T₃,t₃) be the new characteristic point of extracting, the interval (T of these both sides₁,t₁)—(T₃,t₃)、 (T₃,t₃)—(T₅,t₅) need to recalculate orthogonal distance, and interval (T₅,t₅)—(T₈,t₈) orthogonal distance upperWhen an extract minutiae, calculated, only need to and just have the orthogonal distance of new calculatingHand over distance to sort relatively, can obtain new turning point, by that analogy.

Step 5: the turning point that between batch, each variable extracts carries out parallelism inspection;

Successively the trend reverse point of each batch of extraction is carried out to comparison test, when the trend reverse point two extractingWhen side variation tendency has obvious difference, need to replace this trend reverse point, again extract. Newly becomeThe extraction of gesture turning point will contrast the result of parallelism inspection and search, according to other batch of corresponding turnoverThe both sides variation tendency of point is extracted.

Step 6: by batch in all variablees the polymerization of trend reverse point and carry out true and false verification;

(1) polymerization: get union, after polymerization total trend reverse count into:

N_ReasonThe trend reverse that=variable number J × single variable the extracts f that counts;

In fact trend reverse count into

WhereinRepresent the number of times that the same time repeats, q represents the characteristic time number duplicating.

(2) true and false verification: after polymerization, the phenomenon that there will be indivedual trend reverse points to repeat. When becoming of extractingWhen gesture turning point has repetition or approaches very much in the total batch of number that exceedes 2/3, illustrate in these pointsThere is false point, these points need to be rejected, only leave one of them point; When the trend reverse point extracting is fewWhen duplicating or approach very much in 1/3 total batch of number, explanation has trend reverse point here really,Be true point, the point of repetition need to be replaced again and search, and again searches by variation tendency. For example, before4 batches each batch gathers two variablees, and data are respectively data (400 × 2), and (420 × 2),(420 × 2), (460 × 2), each variable extracts respectively 6 turning points, except starting point, and each variableBe equivalent to extract 4 turning points, because same lot data initial point position is identical, carry for each batchIt is 10 that the total characteristic of getting is counted,

There is position in 10 turning points of first secondary data:

1,20,57,96,178,235,235,279,357,400

There is position in 10 turning points of second batch secondary data:

1,34,65,123,196,240,240,279,357,420

There is position in 10 turning points of the 3rd lot data:

1,38,76,135,186,247,247,280,358,420

There is position in 10 turning points of the 4th lot data:

1,56,89,145,145,269,298,346,397,460

Wherein in first, second, and third batch, the time of occurrence of the 6th and the 7th turning point is identical, and repeat number of times exceeded 2/3 total batch of number, one in two turning points will be from thisIn three batches, remove. Become

1,20,57,96,178,235,279,357,400

1,34,65,123,196,240,279,357,420

1,38,76,135,186,247,280,358,420

1,56,89,145,145,269,346,397,460

In the 4th batch, the time of occurrence of the 4th and the 5th turning point is identical, but the number of times repeatingBe less than 2/3 total batch of number, need in the 4th batch, fill out a turning point, with fill out this turnBreak replaces the turning point of origin-location, can re-start and search by both sides variation tendency.

Step 7: at every one-phase, set up based on trend reverse point for each lot dataSingle-MDPCA model.

Dynamic characteristic in consideration process, first constructs ARMAEX time series models, uses time lagData extending is analyzed data matrix structure time series augmented matrix:

$X_d=\left[\begin{array}{cccc}{x}^T\left(k\right)& x^T(k-1)& L& x^T(k-d)\\ x^T(k-1)& x^T(k-2)& L& x^T(k-d-1)\\ M& M& M& M\\ x^T(k+d-K)& x^T(k+d-K-1)& L& x^T(k-K)\end{array}\right]$

In formula, x (k)=[x_1,kx_2,kLx_J,k]^TThe observation vector of J variable at sampling instant k place, at kThe turning point vector of J variable of moment extraction, K is the sampled point of each batch of modeling data matrixNumber, d is time lag length.

(2) augmented matrix is carried out to PCA processing, due to single batch of modeling stage by stage, therefore be calledSingle-MDPCA model.

Step 8: at every one-phase, all single-MDPCA models of setting up are pressed to loading direction phaseCarry out cluster like property, can adopt two cosine angles between model projection vector to measure the phase of two modelsSeemingly degree:

$\cos (P_i^m,P_k^m)=\frac{<P_i^m,P_k^m>}{\sqrt{\underset{m=1}{\overset{p_{\max }}{\&Sigma;}}{\left(P_i^m\right)}^2}\sqrt{\underset{m=1}{\overset{p_{\max }}{\&Sigma;}}{\left(P_k^m\right)}^2}},i,k=1,2,L,I;m=1,2,L,p_{\max }.$

Can obtain like this two angles between load matrix correspondence position projection vector is

${\mathrm{\&theta;}}_m=\left\{\begin{array}{cc}\frac{180\&CenterDot;\arccos (P_i^m,P_k^m)}{\mathrm{\&pi;}}& \left(cos\right(P_i^m,P_k^m)\&GreaterEqual;0)\\ 180-\frac{180\&CenterDot;\arccos (P_i^m,P_k^m)}{\mathrm{\&pi;}}& \left(cos\right(P_i^m,P_k^m)<0)\end{array}\right.$

In formula, P represents load matrix,Represent P_iM row,Represent P_kM row, p_maxForMaximum pivot number.

Consider the importance of each row in load, introduce weight coefficient ω_m：

${\mathrm{\&omega;}}_m=\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_m}/\underset{m=1}{\overset{p_{\max }}{\&Sigma;}}\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_m}$

Wherein, λ_mIt is nonzero eigenvalue.

Can calculate two angles between model projection vector:

Adopt k-means algorithm cluster that the load of weighting is divided into different classes;

Step 9: will belong to of a sort lot data and launch to re-establish model by variable, and form multimodeType, is TTP-basedMDPCA model, calculates T²Control limit with SPE statistic and each class:

T²Statistic is defined as:

$T^2=\&lsqb;t_1,L,t_{p_{max}}\&rsqb;{\mathrm{\&Lambda;}}^{-1}{\&lsqb;t_1,L,t_{p_{max}}\&rsqb;}^T$

Here t is score vector, Λ^-1P before representing_maxThe diagonal matrix of individual characteristic vector institute character pair value compositionContrary.

T²The control limit of statistic can distribute and draw according to F:

$T^2~\frac{p_{max}(K-1)}{K-p_{max}}{F}_{p_{max},K-p_{max},\mathrm{\&alpha;}}$

In formula,Be to be α corresponding to insolation level, the free degree is p_max，K-p_maxF under conditionThe critical value distributing.

SPE statistic is defined as:

$SPE=\left|\right|\mathrm{\&Phi;}\left(x\right)-{\widehat{\mathrm{\&Phi;}}}_a\left(x\right)|{|}^2=\underset{k=1}{\overset{K}{\&Sigma;}}{t_k}^2-\underset{k=1}{\overset{p_{max}}{\&Sigma;}}{t_k}^2$

The control limit of SPE statistic can be according to x²Distribution draws:

${\mathrm{SPE}}_{\mathrm{\&alpha;}}~{\mathrm{g\&chi;}}_h^2;g=\frac{u}{2v},h=\frac{2v^2}{u}$

In formula, u and v are average statistical and the variances of modeling data SPE statistics value.

After establishing weak fault model, can use this model to carry out on-line monitoring to process data, specifically stepRapid as follows:

Step 1, collection new lot data;

Step 2, in the starting stage, with all class models in 10% stage length range to new lotData are monitored, and are every model I data at the new lot data length of starting stage online acquisition10% of length, then brings these data in every model I into, calculates respectively SPE and T²And controlLimit, and adopt an overall target as monitoring standard:

SI＝SPE/CONs＋T²/CON_T，

CON_sAnd CON_TBe respectively SPE and T²Control limit. According to monitoring result, correspondence is minimumThe model of SI is as best monitoring model;

Step 3, best image data successively substitution monitoring model is calculated to T²With SPE statistic. Work as numberAccording to SPE value and T²Exceed and controlled limit, in process, fault has occurred.

Set forth process of the present invention below as an example of penicillin fermentation example:

Obtain the sweat data of penicillin from pharmaceutical factory, therefrom choose 30 batches of normal data formations and buildModulus certificate. Choose following process variables as primary variables: pH value, temperature, dissolved oxygen solubility,CO₂Concentration, power of agitator, culture medium capacity, ventilation rate, substrate concentration, cooling water flow, monitoringSystem as shown in Figure 9.

At the scene in practical application, due to instrumentation reliability, certainty of measurement and in-site measurement environmentEtc. the impact of factor, in measurement data, inevitably contain various noises. For improving signal to noise ratioWith the accuracy that ensures to extract trend reverse point, must first carry out denoising to signal. Small echo is at signalIn denoising, be widely used, and obtained good effect. Wavelet transformation is a kind of differentiates moreAnalyze, when yardstick hour, temporal resolution is high, is suitable for analysis of high frequency signal; In the time that yardstick is larger,Frequency resolution is high, is suitable for analysing low frequency signal, can see process overall picture, and this makes to become based on small echoThe detector changing has good robustness. This characteristic makes it aspect signal denoising, show obvious advantage.Therefore, first adopt the small echo threshold values denoising method that Donoho proposes to remove the noise in signal.

Using extraction algorithm based on trend reverse point to carry out turning point to variable the data acquisition after denoising carriesGet, as shown in accompanying drawing 4-8. Sweat is divided into two stages: thalli growth stage and penicillinSynthesis phase. The time lag length in two stages elects respectively 2 and 3 as. At every one-phase, by PCA andARMAEX time series models combine, and all lot data are carried out to single batch of modeling, set upSingle-MDPCA model. Then all single-MDPCA models of setting up are pressed to loading direction similarProperty is carried out cluster. To belong to of a sort lot data and launch to re-establish model by variable, formTTP-basedMDPCA model, calculates T²Control limit with SPE statistic and each class.

In order to verify that the present invention sets up the accuracy of model, we adopt and put menstruation with the model of above-mentioned foundationPut down is that 95% control limit is monitored historical production data. Figure 10, Figure 11 have provided wherein one groupThe monitoring result of historical production data. All SPE value and T²Value all, controlling under limit, shows this batchInferior be normal within the whole service time.

Figure 12, Figure 13 have provided the monitoring result of a fault batch. Can from monitoring result figureTo SPE value and T²Value has obviously exceeded during 100h to 300h controls limit, and fault has occurred in process.Find on inspection it is that small step saltus step has occurred power of agitator during 100h to 300h, make processData are offset, and cause process dependency relation that variation has occurred. Visible the present invention can accurately identifyThe generation of this weak fault. In sum, the present invention set up based on trend reverse point (TTP)The weak malfunction monitoring model of MDPCA, can meet the actual production process requirement of monitoring in real time.

Although embodiment of the present invention are open as above, it is not restricted to description and enforcement sideListed utilization in formula, it can be applied to various applicable the field of the invention completely, for being familiar with abilityThe personnel in territory, can easily realize other amendment, therefore do not deviate from claim and etc. homotypeEnclose under limited universal, the present invention is not limited to specific details and illustrates here and the figure describingExample.

Claims (9)

1. the weak fault model control limit of a batch process method for building up, is characterized in that, comprises following stepRapid:

The data of step 1, I batch of normal operation of collectionWherein i=1 ..., I, K_iForThe sampling number of i batch, J is the variable number gathering, R is real number set;

Step 2, adopt CPV1 discrete method to carry out segmentation each batch process lot data;

Step 3, the method for employing based on trend reverse point are carried out turning point extraction, and are entered by orthogonal distanceLine ordering;

Step 4, by batch between the turning point that extracts of each variable carry out parallelism inspection, when the trend of extracting turnsWhen break both sides variation tendency has obvious difference, need to replace this trend reverse point, again extract;

Step 5, by batch in all variablees the polymerization of trend reverse point and carry out true and false verification;

Step 6, set up single-MDPCA model for each lot data based on trend reverse point,Use time lag data extending to analyze data matrix structure time series augmented matrix, and to time series augmentation squareBattle array is carried out PCA processing;

Step 7, at every one-phase, by set up all single-MDPCA models press loading direction phaseCarry out cluster like property;

Step 8, calculating T²Control limit with SPE statistic and each class.

2. the weak fault model control limit of batch process according to claim 1 method for building up, its featureBe, between step 1 and step 2, also comprise: each lot data is adopted to Donoho's by variableThe denoising of threshold values denoising method.

3. the weak fault model control limit of batch process according to claim 1 method for building up, its featureBe, in step 3, turning point extracting method is: for time series, first starting point and terminal are doneFor turning point saves; Then between starting point and terminal, searching the next one by maximum orthogonality distance method becomesGesture turning point; In the new turning point both sides of finding, again find out respectively the turnover with maximum orthogonality distancePoint, and compare by distance, the point with maximum orthogonality distance is extracted a little as the next one.

4. the weak fault model control limit of batch process according to claim 1 method for building up, its featureBe, in step 5, the process of true and false verification is: when the trend reverse point extracting exceed 2/3 alwaysWhen having repetition or approach very much in batch number, these points are rejected, only leave one of them point; When carryingWhen the trend reverse point of getting duplicates or approaches very much in the total batch of number that is less than 1/3, the point of repetitionReplace again and search.

5. the weak fault model control limit of batch process according to claim 1 method for building up, its featureBe, in step 6, time series augmented matrix is

$X_d=\left[\begin{array}{cccc}{x}^T\left(k\right)& x^T(k-1)& L& x^T(k-d)\\ x^T(k-1)& x^T(k-2)& L& x^T(k-d-1)\\ M& M& M& M\\ x^T(k+d-K)& x^T(k+d-K-1)& L& x^T(k-K)\end{array}\right]$

In formula, x (k)=[x_1,kx_2,kLx_J,k]^TThe observation vector of J variable at sampling instant k place, at kThe turning point vector of J variable of moment extraction, K is the sampled point of each batch of modeling data matrixNumber, d is time lag length.

6. the weak fault model control limit of batch process according to claim 1 method for building up, its featureBe, in step 7, all single-MDPCA models carry out the side of cluster by loading direction similitudeMethod comprises:

Adopt cosine angle between two model projection vectors to measure the similarity degree of two models:

$\begin{array}{c}cos(P_i^m,P_k^m)=\frac{<P_i^m,P_k^m>}{\sqrt{\underset{m=1}{\overset{p_{max}}{\&Sigma;}}{\left(P_i^m\right)}^2}\sqrt{\underset{m=1}{\overset{p_{max}}{\&Sigma;}}{\left(P_k^m\right)}^2}}\\ i,k=1,2,L,I;m=1,2,L,p_{\max }.\end{array}$

In formula, P represents load matrix,Represent P_iM row,Represent P_kM row, p_maxForMaximum pivot number;

Introduce weight coefficient ω_m：

${\mathrm{\&omega;}}_m=\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_m}/\underset{m=1}{\overset{p_{max}}{\&Sigma;}}\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_m}$

Wherein, λ is nonzero eigenvalue;

Can calculate two angles between model projection vector:

${\mathrm{\&theta;}}_{i,k}=\underset{m=1}{\overset{p_{max}}{\&Sigma;}}{\mathrm{\&omega;}}_m{\mathrm{\&theta;}}_m;$

Finally adopt k-means algorithm cluster that the load of weighting is divided into different classes.

7. the weak fault model control limit of batch process according to claim 1 method for building up, its featureBe, in step 8,

T²Statistic is:

$T^2=\&lsqb;t_1,L,t_{p_{max}}\&rsqb;{\mathrm{\&Lambda;}}^{-1}{\&lsqb;t_1,L,t_{p_{max}}\&rsqb;}^T$

Here t is score vector, Λ^-1P before representing_maxThe diagonal matrix of individual characteristic vector institute character pair value compositionContrary;

T²The control limit of statistic can distribute and draw according to F:

$T^2~\frac{p_{max}(K-1)}{K-p_{max}}{F}_{p_{max},K-p_{max},\mathrm{\&alpha;}}$

In formula,Be to be α corresponding to insolation level, the free degree is p_max，K-p_maxF under conditionThe critical value distributing.

8. the weak fault model control limit of batch process according to claim 1 method for building up, its featureBe, in step 8,

SPE statistic is:

$SPE=\left|\right|\mathrm{\&Phi;}\left(x\right)-{\widehat{\mathrm{\&Phi;}}}_a\left(x\right)|{|}^2=\underset{k=1}{\overset{K}{\&Sigma;}}{t_k}^2-\underset{k=1}{\overset{p_{\max }}{\&Sigma;}}{t_k}^2$

The control limit of SPE statistic can be according to χ²Distribution draws:

${\mathrm{SPE}}_{\mathrm{\&alpha;}}~{\mathrm{g\&chi;}}_h^2;g=\frac{u}{2v},h=\frac{2v^2}{u}$

In formula, u and v are average statistical and the variances of modeling data SPE statistics value.

9. the weak fault monitoring method of batch process, is characterized in that, right to use requires to appoint in 1-8The weak fault model control limit of the batch process method for building up providing in one, and comprise the steps:

Step 1, Real-time Collection process data;

Step 2, in 10% length range, new lot data are monitored with all class models, and meterCalculate overall target

SI＝SPE/CONs+T²/CON_T

Wherein, CON_sAnd CON_TBe respectively SPE and T²Control limit; By the model of the minimum SI of correspondenceAs best monitoring model;

The best monitoring model of choosing in step 3, use step 2 is monitored process data, works as numberAccording to SPE value and T²Exceed and controlled limit, in process, fault has occurred.