CN103116306B - Automatic stepping type ordered time interval dividing method - Google Patents

Automatic stepping type ordered time interval dividing method Download PDF

Info

Publication number
CN103116306B
CN103116306B CN201310046432.1A CN201310046432A CN103116306B CN 103116306 B CN103116306 B CN 103116306B CN 201310046432 A CN201310046432 A CN 201310046432A CN 103116306 B CN103116306 B CN 103116306B
Authority
CN
China
Prior art keywords
time
matrix
mrow
monitoring
new
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201310046432.1A
Other languages
Chinese (zh)
Other versions
CN103116306A (en
Inventor
赵春晖
李文卿
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
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 Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN201310046432.1A priority Critical patent/CN103116306B/en
Publication of CN103116306A publication Critical patent/CN103116306A/en
Application granted granted Critical
Publication of CN103116306B publication Critical patent/CN103116306B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Testing And Monitoring For Control Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses an automatic stepping type ordered time interval dividing method which can be used for dividing an interval of a multi-operation procedure into different sub time intervals according to change of process characteristics. Each time interval has similar process characteristic and can be represented by the same model, while different time intervals have different characteristics so that different models are needed to be built, and accordingly, the number of the models and complexity are reduced greatly. Online monitoring indexes are used as the judgment basis for time interval dividing, and thus, precision of the models for the time interval dividing is improved, and subsequent process online monitoring performance is improved greatly. In addition, the automatic stepping type ordered time interval dividing method is easy to implement and is applied successfully in the mold injection process, assists people to learn the specific process characteristics, enhances reliability and confidence of the actual online process monitoring, assists engineers to make correct judgment for the process running states and find out fault in time, and accordingly guarantees safe and reliable running of actual production and high quality of products.

Description

Automatic stepping type sequential time interval dividing method
Technical Field
The invention belongs to the field of intermittent process statistical monitoring, and particularly relates to a method for automatically dividing time intervals in an intermittent production process in multiple operation time intervals.
Background
As an important production mode in industrial production, the intermittent process is closely related to the life of people and is widely applied to the fields of fine chemical engineering, biological pharmacy, food, polymer reaction, metal processing and the like. In recent years, with the more urgent market demand of modern society for various products, various specifications and high-quality products, industrial production is more focused on the intermittent process for producing small-lot and high-value-added products. The safe and reliable operation of batch production and the high quality pursuit of products have become the focus of attention. Multivariate statistical analysis techniques based on data are increasingly favored by researchers and field engineers because they only require process data under normal operating conditions to build models, and they have unique advantages in processing high-dimensional, highly coupled data. Statistical modeling, online monitoring, fault diagnosis, and quality prediction of intermittent processes have become a subject of extensive research. However, the statistical analysis methods of multi-way principal component analysis (MPCA) and multi-way minimum-deviation-two-times regression (MPLS) treat all data of one intermittent operation as a whole and cannot reflect the change of the correlation of variables in the time direction. In addition, unknown measurement data must be estimated in online application, so that online implementation is difficult to realize, and wide application of the method in actual production is prevented.
It has been found through research that intermittent operation processes have a time-of-day characteristic. The process variable correlation in the intermittent operation does not change along with time and moment, but changes regularly along with the change of the process operation progress or the process mechanism characteristic, and shows sectionalization. In different time periods, each time period has different process variable tracks, operation modes and correlation characteristics, and the correlation of the variables is obviously different. The correlation of the process variables at different sampling times is approximately the same during the same time period. In consideration of the multi-period characteristics of the intermittent process, people expect that the whole intermittent process can be divided into different periods, so that a plurality of models based on the periods can be established and used for on-line process monitoring, so that the occurrence of faults can be timely and effectively detected, and the safety of industrial production is ensured. Take a typical multi-time batch process-injection molding-as an example. The injection molding process is a typical batch industrial process and is also a specific application background process of the time division method proposed in the patent. The construction of the injection molding machine and the basic operating principle of the injection molding process will be briefly described here. Injection molding is one of the major molding methods for thermoplastic or thermoset plastic articles. About 1/3, the plastic products in our lives are manufactured by injection molding. A universal injection molding machine mainly comprises an injection device, a mold closing device, a hydraulic system and an electric control system. As a typical batch process with multiple operation stages, a complete injection molding process consists of procedures such as mold closing, advancing of an injection seat, injection, pressure maintaining, plasticizing, cooling, mold opening, ejection of a part, and the injection stage, the pressure maintaining stage and the cooling stage are the most important three operation stages for determining the quality of the part. In the injection section, the hydraulic system pushes the screw to inject the viscous fluid of plastic into the mold cavity until the mold cavity is filled with fluid. When the process is in the pressure maintaining section, a small amount of viscous fluid is still extruded into the mold cavity under high pressure to compensate the volume shrinkage of the viscous fluid of the plastic during cooling and plasticizing. The pressure maintaining stage is continued until the gate of the mold cavity is frozen, and the process enters a cooling stage. While the fluid in the die cavity is solidified in the cooling stage, the plastic particles in the machine barrel realize the change of the physical state under the action of a heating device outside the machine barrel and the shearing heat generated by the rotation of the screw rod, become a plastic viscous state and are conveyed to the head part of the screw rod. When the melt at the head of the screw is gradually increased and the pressure of the melt is greater than the back pressure of the injection oil cylinder, the screw retreats and simultaneously starts the volume calculation. After the head melt reaches a certain injection quantity, the screw stops moving back and rotating, and the process state of the time is also called a plasticizing stage. As the melt in the mold cavity continues to cool, the plastic returns from a viscous state to a glassy state and sets. When the plastic part is completely solidified, the mold is opened, and the plastic part is ejected, thereby completing a working cycle.
How to reasonably divide an intermittent process into different sub-periods is the basis and key of subsequent statistical modeling and real-time fault detection based on the periods, and is directly related to the accuracy and reliability of process monitoring. The predecessors have made corresponding research and discussion on the above, and have proposed corresponding time interval identification methods based on different angles. The following three types can be summarized: (a) a method based on process mechanism knowledge and expert experience (b) a characteristic signal variable analysis method (c) an automatic identification method. The time interval identification methods have various application occasions and advantages and disadvantages. The clustering method is a relatively common time interval automatic division method, and focuses on analyzing and tracking the change of relevant characteristics of the process. However, the clustering-based time interval division method does not take the time sequence of the time intervals into consideration. If only the similar characteristics of the process are considered, the sampling moments in different time regions may be wrongly divided into the same sub-period, and the time points in the same time region may also be divided into different time periods, which may cause the division result to be complicated and difficult to understand. In addition, these time interval division methods do not take into account the transition characteristics between time intervals, but strictly divide time slices into a certain time interval. In the transition region where the process dependency changes, if the process mode belonging to the time is divided into two subclasses, it may cause "misclassification". On one hand, the inclusion of the transition mode in the subinterval can affect the accuracy of the subinterval model; on the other hand, if the transition pattern is annihilated in a subinterval, the probability of the first type of detection error is increased when the subinterval transition region is monitored online. In other words, the previous time interval automatic division method does not consider the time sequence and the time interval transitivity of the process operation, so that the subsequent process modeling precision and the monitoring performance are directly or intermittently influenced.
The invention further considers the time-varying property of the potential characteristics of the intermittent process, the time sequence of the actual process operation and the influence of the time interval division result on the subsequent monitoring performance, and provides a novel automatic sub-time interval division method. So far, no research report related to the invention is seen.
Disclosure of Invention
The invention aims to provide an automatic stepping type sequential time interval dividing method aiming at the defects of the existing time interval dividing technology aiming at the intermittent production process. The method can automatically capture the development change of potential process characteristics according to the running sequence of the intermittent production process, determine local time blocks in the time direction, namely a sub-period and a transition mode, and model based on the period division result, can improve the model precision and improve the on-line process monitoring performance, and can be finally applied to an actual industrial production field to ensure the safe and reliable running of intermittent production and the high-quality pursuit of products.
The purpose of the invention is realized by the following technical scheme: an automatic stepwise sequential time interval division method, comprising the steps of:
step 1: acquiring process analysis data: if an intermittent operation is provided with J measurement variables and K sampling points, each measurement batch can obtain a K multiplied by J matrix, and after the measurement steps of the I batch are repeated, the obtained data can be expressed as a three-dimensional matrixWherein the measurement variables are temperature, speed, pressure, displacement and other state parameters which can be measured in the batch operation process;
step 2: data preprocessing: combining three-dimensional matricesExpanding according to the collection batch direction, namely arranging the variables on each sampling point in an operation batch according to the time sequence to obtain a two-dimensional matrix X (I multiplied by KJ), and obtaining a matrix X by K time slicesk(I × J), where subscript k is a time indicator;
let two-dimensional matrix XkThe variable at any point in the table is xk,i,jThe variable is normalized by subtracting the mean value and dividing by the standard deviation, wherein the subscript i represents the batch and j represents the variable, and the calculation formula of the normalization process is as follows:
<math> <mrow> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>x</mi> <mo>&OverBar;</mo> </mover> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow> <msub> <mi>s</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mfrac> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein: k is a time slice indicator.Is XkMean of any column of the matrix, sk,jIs XkThe standard deviation of the moment for the corresponding column,
<math> <mrow> <msub> <mover> <mi>x</mi> <mo>&OverBar;</mo> </mover> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>I</mi> </mfrac> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>I</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>,</mo> </mrow> </math>
<math> <mrow> <msub> <mi>s</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msqrt> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>I</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>x</mi> <mo>&OverBar;</mo> </mover> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>/</mo> <mrow> <mo>(</mo> <mi>I</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msqrt> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow> </math>
and step 3: time slice PCA modeling, which is implemented by the following sub-steps:
(3.1) normalizing each time slice matrix X processed in the step 2k(I × J) performing PCA decomposition to establish a time slice PCA model, wherein the PCA decomposition formula is as follows:
<math> <mrow> <msub> <mi>X</mi> <mi>k</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>k</mi> </msub> <msup> <msub> <mi>P</mi> <mi>k</mi> </msub> <mi>T</mi> </msup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>t</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>r</mi> </mrow> </msub> <msub> <mi>p</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>r</mi> </mrow> </msub> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein: t is tk,nBeing orthogonal principal component vectors, pk,nThe load vector is orthogonal normalization, r represents different PCA decomposition directions, and superscript T represents the transposition of the matrix; t isk(I × J) represents a score matrix that retains all of the pivot elements, Pk(J × J) represents the corresponding load matrix.
(3.2) selecting the number of the principal elements, and restating the formula (3) into the following form:
<math> <mrow> <msub> <mi>X</mi> <mi>k</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>k</mi> </msub> <msup> <msub> <mi>P</mi> <mi>k</mi> </msub> <mi>T</mi> </msup> <mo>+</mo> <msub> <mi>E</mi> <mi>k</mi> </msub> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>R</mi> </munderover> <msub> <mi>t</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>r</mi> </mrow> </msub> <msub> <mi>p</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>r</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>E</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein: r represents different PCA decomposition directions; t isk(I×Rk) And Pk(J×Rk) Representing the load matrix as reserved RkScore matrix and load matrix after an individual principal element, EkIs a residual matrix. Through the transformation, the multi-directional principal component analysis model decomposes an original data space into a principal component space and a residual error space, wherein the principal component space represents main system process fluctuation information; number of principal elements R reserved herekCan reflect 90% of the process fluctuation information in the original process. Comprehensively considering the whole process, selecting the R with the most occurrence timeskThe values are unified into the final modeling pivot number R, namely, the same pivot number is reserved for all time slice models.
(3.3) calculating SPE indexes corresponding to each batch in each time slice k in the residual error space:
SPEk,i=ek,i Tek,i (5)
wherein the index i denotes the batch in the time slice, ek,iIs the residual column vector for the ith batch at time k. Obeying chi with weighting coefficient according to SPE values of different batches at the same time2Distribution to determine the control limit Ctr at each time pointkIt reflects the time slice PCA modelType of reconstitution capability.
And 4, step 4: determining SPE index control limit based on variable expansion model: from the initial point of the intermittent process, the next time slice is combined with the previous time slice in turn and expanded in a variable mode to obtain a time block Xv,k(Ik J), where the subscript v represents the way the variables expand. Carrying out PCA analysis on the new time block matrix to extract a load matrix Pv,k(J × R). Calculating SPE value and obeying chi with weighting coefficient according to SPE values of different batches at the same time2Distribution to determine the control limit Ctr at each time pointk
And 5: determining a first time period division point k*: comparison of Ctr at each time point within the same time regionkAnd Ctrv,kIf three consecutive samples are found to exhibit Ctrv,k>αCtrkThen the newly added time slice has a significant impact on the PCA monitoring model and the corresponding monitoring performance for that time block. The time before adding the new time slice is recorded as k*. Where α is dependent on CtrkIs called the mitigation factor, which reflects the extent to which the time block model allows monitoring of the loss of accuracy compared to the time slice model. Then k is*The time slice before the moment of time may be considered a sub-period.
Step 6: the process analyzes the data updates and determines all of the divided periods: according to the time k obtained in step 5*Removing the first sub-period, bringing the remaining intermittent process data as new input data into step 5 and repeating the above steps 5-6 for different time periods until no data remains.
And 7: and (3) statistical modeling based on time interval division results: according to the time interval division result in the step 6, time slices in each time interval are combined into a sub-time interval representative modeling data group in a variable expansion mode, Xc(IKcX J), where subscript c is a period indicator. Similar process potential characteristics in each time interval can be obtained by pairing Xc(IKcxJ) carrying out PCA divisionPerforming solution extraction:
Tc=XcPc
X ^ c = T c P c T = X c P c P c T - - - ( 6 )
E c = X c - X ^ c
wherein, Pc(J×Rc) Is the PCA load matrix of the sub-period, revealing the main fluctuation direction, R, within the sub-periodcRepresenting the number of PCA principal components retained by the subinterval model. T isc(IKc×Rc) A pivot score matrix representing the subinterval. Thus, the main fluctuation of the period is through TcPc TThe characterization represents the main fluctuation information in the sub-period; the remainder is the sub-period residual matrix, Ec(IKcX J) represents the noise information within the sub-period. Each time slice principal score, Tk(I×Rc) The matrix T can be easily scored from the sub-periodsc(IKc×Rc) According to the corresponding extraction of the process time,thereby obtaining the covariance relation S at each sampling moment by calculationk. Residual error matrix E of each time slicek(I × J) may also be derived from the sub-period residual matrix Ec(IKcXJ) in (B).
And 8: calculating a real-time monitoring statistical index: according to the principal component score time slice T obtained from the result of the calculation of the formula (6)k(I×Rc) And residual matrix time slice Ek(I × J) two monitoring statistics can be calculated at each time instant: Hotelling-T2Statistical indicators and SPE statistics.
Hotelling-T2The statistical indicator is used for measuring the distance of the process variable deviating from the average track under the normal working condition at each sampling moment, and the control limit under the significance level alpha is calculated as follows:
<math> <mrow> <msubsup> <mi>T</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> </mrow> <mn>2</mn> </msubsup> <mo>-</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>&OverBar;</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <msubsup> <mi>S</mi> <mi>k</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>&OverBar;</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>~</mo> <mfrac> <mrow> <msub> <mi>R</mi> <mi>c</mi> </msub> <mrow> <mo>(</mo> <mi>I</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>I</mi> <mo>-</mo> <msub> <mi>R</mi> <mi>c</mi> </msub> </mrow> </mfrac> <msub> <mi>F</mi> <mrow> <msub> <mi>R</mi> <mi>c</mi> </msub> <mo>,</mo> <msub> <mrow> <mi>I</mi> <mo>-</mo> <mi>R</mi> </mrow> <mi>c</mi> </msub> <mo>,</mo> <mi>&alpha;</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein R iscIs the number of principal elements reserved in the PCA model in the time period; t is ti,k(RcX 1) is the principal component score of the ith batch at time k, i.e. the corresponding time slice score matrix Tk(I×Rc) Row i of (1); whileIs Tk(I×Rc) For mean vectors of different batches, since the measured data of each time slice is already centered to zero mean in the data preprocessing, the mean vector of the time slice isIt is a zero vector.
For the residual subspace, the SPE statistics for different batches at each time are calculated as:
SPE i , k ( x i , k - x ^ i , k ) T ( x i , k - x ^ i , k ) - - - ( 8 )
wherein x isi,kThe process measurement vector representing the ith lot at time k,then the corresponding result is reconstructed from the PCA model. The calculation result in equation (8) may form an I × 1 vector [ SPE ] at each time1,k,SPE2,k,...,SPEI,k]TAnd the vector approximation is subject to a weighting χ2And distributing to obtain the monitoring control limit of the SPE at each moment.
And step 9: on-line process monitoring based on a time period model: and (4) on the basis of the time period divided in the step (6), the time period model monitoring system established in the step (7) and the two monitoring statistics obtained in the step (8), the state of a new operation intermittent process such as injection molding and the like can be monitored on line. This step is realized by the following substeps:
(9.1) acquiring new measurement data and preprocessing the new measurement data: during on-line monitoring, new process measurement data x are collectednewAfter (JX 1), whereinThe new label represents the new sample, J is the measured variable, the same as the measured variable in step 1; first, data preprocessing is required. And (3) calling the mean value and the standard deviation corresponding to the moment according to the mean value and the standard deviation obtained in the step (2) and carrying out standardization preprocessing on the collected process measurement data as shown in a formula (1) according to the indication of the process time.
(9.2) calculating new monitoring statistics: after data preprocessing, calling a model P corresponding to the time interval of the new sampling moment according to a PCA sub-time interval model calculated by the formula (6)c(J×Rc) Calculating to obtain principal component score, estimating residual error and its corresponding Hotelling-T2And SPE two monitoring statistical indexes:
tnew T=xnew TPc
x ^ new = P c t new ( 9 )
<math> <mrow> <msubsup> <mi>T</mi> <mi>new</mi> <mn>2</mn> </msubsup> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>new</mi> </msub> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>&OverBar;</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <msup> <msub> <mi>S</mi> <mi>k</mi> </msub> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>new</mi> </msub> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>&OverBar;</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> </mrow> </math>
SPE new = ( x new - x ^ new ) T ( x new - x ^ new )
wherein,is the new process measurement data that is,is T previously obtained from training datakMean vector of (S)kIs XkThe covariance matrix of (2).
(9.3) judging the process running state on line: comparing the two monitoring indexes with respective statistical control limits in real time, and if the two monitoring indexes are both within the statistical control limits, indicating that the process is normal in operation; if more than one monitoring index exceeds the normal control limit, the abnormal condition of the process is indicated.
Further, in step 1, the measured variables are the following 9: the temperature control device comprises a pressure valve opening, a flow valve opening, a screw stroke, a screw speed, injection pressure, a nozzle temperature, a barrel head temperature, a barrel middle temperature and a barrel tail temperature.
Compared with the prior art, the invention has the beneficial effects that: the invention provides a new research idea for multi-stage process time interval division, modeling and monitoring under the condition of no process prior knowledge. The automatic sub-period division method can be applied to a class of intermittent production processes with multiple operation periods, is divided into different sub-periods according to the change of process characteristics, and can improve the accuracy of the established sub-period model and the performance of on-line monitoring of the subsequent process. The method has the advantages that the method carries out detailed experimental research in the injection molding industrial process, and achieves successful application, improves the understanding of the specific process operation characteristics through automatic division of multiple operation time periods of the intermittent process, improves the monitoring efficiency of the process monitoring process and the accuracy of a fault detection result, can be finally applied to an actual industrial production field, and ensures safe and reliable operation of intermittent production and high quality pursuit of products.
Drawings
FIG. 1 is a flow chart of an automated step-wise sequential time interval division method of the present invention;
FIG. 2 is a schematic diagram of a three-dimensional data expansion model of the automated step-by-step sequential time interval division method of the present invention;
FIG. 3 is a graph comparing the time interval division result of the method of the present invention with the division result of the conventional method in the embodiment of the present invention.
Detailed Description
The invention is further described with reference to the following drawings and specific examples.
The injection molding process is a typical multi-stage batch process, and generally comprises three stages of injection, pressure maintaining and cooling. In addition, the plasticizing process is completed at the initial stage of cooling. Specifically, during the injection phase, the hydraulic system pushes the screw to inject the viscous fluid of plastic into the mold cavity until the mold cavity is filled with fluid. When the process is in the pressure maintaining stage, a small amount of viscous fluid is still extruded into the mold cavity under high pressure to compensate the volume shrinkage of the viscous fluid of the plastic during cooling and plasticizing. The pressure maintaining stage is continued until the gate of the mold cavity is frozen, and the process enters a cooling section. When the melting material at the head of the screw rod is gradually increased and reaches a certain injection amount, the screw rod stops retreating and rotating, and the process state of the period is also called a plasticizing section. As the molten material in the mold cavity continues to cool, the plastic part is fully solidified, the mold is opened, and the plastic part is ejected, thereby completing a work cycle.
The invention discloses an automatic step-by-step ordered time interval dividing method, which comprises the following steps:
1. obtaining process analysis data
If an intermittent operation is provided with J measurement variables and K sampling points, each measurement batch can obtain a K multiplied by J matrix, and after the measurement steps of the I batch are repeated, the obtained data can be expressed as a three-dimensional matrixWherein the measurement variables are temperature, speed, pressure, displacement and other state parameters which can be measured in the batch operation process; in this example, 526 samples were collected, with 9 measurement variables: the temperature control device comprises a pressure valve opening, a flow valve opening, a screw stroke, a screw speed, injection pressure, a nozzle temperature, a barrel head temperature, a barrel middle temperature and a barrel tail temperature. In this example, 30 normal batches are used to test the time interval division method proposed by the present invention and establish a corresponding on-line monitoring system. I.e., three-dimensional modeling data matrix ofIn addition, three batches of faults were experimentally obtained for verifying the online fault detection performance of the established monitoring system, wherein each fault included 49 batches. Three of theseThe seed faults are: the failure of 90% efficiency of the thermocouple, the failure of 80% efficiency of the heating ring under open loop, and the failure of compounding (i.e. adding blue polypropylene to the original high density polyethylene raw material).
Step 2: data preprocessing: combining three-dimensional matricesExpanding according to the collection batch direction, namely arranging the variables on each sampling point in an operation batch according to the time sequence to obtain a two-dimensional matrix X (I multiplied by KJ), and obtaining a matrix X by K time slicesk(I × J), where subscript k is a time indicator;
let two-dimensional matrix XkThe variable at any point in the table is xk,i,jThe variable is normalized by subtracting the mean value and dividing by the standard deviation, wherein the subscript i represents the batch and j represents the variable, and the calculation formula of the normalization process is as follows:
<math> <mrow> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>x</mi> <mo>&OverBar;</mo> </mover> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mrow> <msub> <mi>s</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mfrac> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein: k is a time slice indicator.Is XkMean of any column of the matrix, sk,jIs XkThe standard deviation of the moment for the corresponding column,
<math> <mrow> <msub> <mover> <mi>x</mi> <mo>&OverBar;</mo> </mover> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>I</mi> </mfrac> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>I</mi> </munderover> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>,</mo> </mrow> </math>
<math> <mrow> <msub> <mi>s</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <msqrt> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>I</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>x</mi> <mo>&OverBar;</mo> </mover> <mrow> <mi>k</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>/</mo> <mrow> <mo>(</mo> <mi>I</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msqrt> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow> </math>
and step 3: time slice PCA modeling, which is implemented by the following sub-steps:
(3.1) normalizing each time slice matrix X processed in the step 2k(I × J) performing PCA decomposition to establish a time slice PCA model, wherein the PCA decomposition formula is as follows:
<math> <mrow> <msub> <mi>X</mi> <mi>k</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>k</mi> </msub> <msup> <msub> <mi>P</mi> <mi>k</mi> </msub> <mi>T</mi> </msup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>J</mi> </munderover> <msub> <mi>t</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>r</mi> </mrow> </msub> <msub> <mi>p</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>r</mi> </mrow> </msub> <mo>;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein: t is tk,nBeing orthogonal principal component vectors, pk,nThe load vector is orthogonal normalization, r represents different PCA decomposition directions, and superscript T represents the transposition of the matrix; t isk(I × J) represents a score matrix that retains all of the pivot elements, Pk(J × J) represents the corresponding load matrix.
(3.2) selecting the number of the principal elements, and restating the formula (3) into the following form:
<math> <mrow> <msub> <mi>X</mi> <mi>k</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>k</mi> </msub> <msup> <msub> <mi>P</mi> <mi>k</mi> </msub> <mi>T</mi> </msup> <mo>+</mo> <msub> <mi>E</mi> <mi>k</mi> </msub> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>R</mi> </munderover> <msub> <mi>t</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>r</mi> </mrow> </msub> <msub> <mi>p</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>r</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>E</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein: r represents different PCA decomposition directions; t isk(I×Rk) And Pk(J×Rk) Representing the load matrix as reserved RkScore matrix and load matrix after an individual principal element, EkIs a residual matrix. Through the transformation, the multi-directional principal component analysis model decomposes an original data space into a principal component space and a residual error space, wherein the principal component space represents main system process fluctuation information; number of principal elements R reserved herekCan reflect 90% of the process fluctuation information in the original process. Comprehensively considering the whole process, selecting the R with the most occurrence timeskThe values are unified into the final modeling pivot number R, namely, the same pivot number is reserved for all time slice models.
(3.3) calculating SPE indexes corresponding to each batch in each time slice k in the residual error space:
SPEk,i=ek,i Tek,i (5)
wherein the index i denotes the batch in the time slice, ek,iIs the residual column vector for the ith batch at time k. Obeying chi with weighting coefficient according to SPE values of different batches at the same time2Distribution to determine the control limit Ctr at each time pointkIt reflects the reconstruction capability of the time slice PCA model.
And 4, step 4: determining SPE index control limit based on variable expansion model: from the initial point of the intermittent process, the next time slice is combined with the previous time slice in turn and expanded in a variable mode to obtain a time block Xv,k(Ik J), where the subscript v represents the way the variables expand. Carrying out PCA analysis on the new time block matrix to extract a load matrix Pv,k(J × R). Calculating SPE value and obeying chi with weighting coefficient according to SPE values of different batches at the same time2Distribution to determine the control limit Ctr at each time pointk
And 5: determining a first time period division point k*: comparison of Ctr at each time point within the same time regionkAnd Ctrv,kIf three consecutive samples are found to exhibit Ctrv,k>αCtrkThen the newly added time slice has a significant impact on the PCA monitoring model and the corresponding monitoring performance for that time block. The time before adding the new time slice is recorded as k*. Where α is dependent on CtrkIs called the mitigation factor, which reflects the extent to which the time block model allows monitoring of the loss of accuracy compared to the time slice model. Then k is*The time slice before the moment of time may be considered a sub-period.
Step 6: the process analyzes the data updates and determines all of the divided periods: according to the time k obtained in step 5*Removing the first sub-period, bringing the remaining intermittent process data as new input data into step 5 and repeating the above steps 5-6 for different time periods until no data remains.
And 7: and (3) statistical modeling based on time interval division results: according to the time interval division result in the step 6, time slices in each time interval are combined into a sub-time interval representative modeling data group in a variable expansion mode, Xc(IKcX J), where subscript c is a period indicator. Similar process potential characteristics in each time interval can be obtained by pairing Xc(IKcX J) PCA extraction:
Tc=XcPc
X ^ c = T c P c T = X c P c P c T - - - ( 6 )
E c = X c - X ^ c
wherein, Pc(J×Rc) Is the PCA load matrix of the sub-period, revealing the main fluctuation direction, R, within the sub-periodcRepresenting the number of PCA principal components retained by the subinterval model. T isc(IKc×Rc) A pivot score matrix representing the subinterval. Thus, the main fluctuation of the period is through TcPc TThe characterization represents the main fluctuation information in the sub-period; the remainder is the sub-period residual matrix, Ec(IKcX J) represents the noise information within the sub-period. Each time slice principal score, Tk(I×Rc) The matrix T can be easily scored from the sub-periodsc(IKc×Rc) According to the corresponding process time, the covariance relation S at each sampling moment can be calculated and obtainedk. Residual error matrix E of each time slicek(I × J) may also be derived from the sub-period residual matrix Ec(IKcXJ) in (B).
And 8: calculating real-time monitoring statistical index
According to the principal component score time slice T obtained from the result of the calculation of the formula (6)k(I×Rc) And residual matrix time slice Ek(I × J) two monitoring statistics can be calculated at each time instant: Hotelling-T2Statistical indicators and SPE statistics.
Hotelling-T2The statistical indicator is used for measuring the distance of the process variable deviating from the average track under the normal working condition at each sampling moment, and the control limit under the significance level alpha is calculated as follows:
<math> <mrow> <msubsup> <mi>T</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> </mrow> <mn>2</mn> </msubsup> <mo>-</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>&OverBar;</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <msubsup> <mi>S</mi> <mi>k</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>&OverBar;</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>~</mo> <mfrac> <mrow> <msub> <mi>R</mi> <mi>c</mi> </msub> <mrow> <mo>(</mo> <mi>I</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>I</mi> <mo>-</mo> <msub> <mi>R</mi> <mi>c</mi> </msub> </mrow> </mfrac> <msub> <mi>F</mi> <mrow> <msub> <mi>R</mi> <mi>c</mi> </msub> <mo>,</mo> <msub> <mrow> <mi>I</mi> <mo>-</mo> <mi>R</mi> </mrow> <mi>c</mi> </msub> <mo>,</mo> <mi>&alpha;</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein R iscIs the number of principal elements reserved in the PCA model in the time period; t is ti,k(RcX 1) is the principal component score of the ith batch at time k, i.e. the corresponding time slice score matrix Tk(I×Rc) Row i of (1); whileIs Tk(I×Rc) For mean vectors of different batches, since the measured data of each time slice is already centered to zero mean in the data preprocessing, the mean vector of the time slice isIt is a zero vector.
For the residual subspace, the SPE statistics for different batches at each time are calculated as:
SPE i , k ( x i , k - x ^ i , k ) T ( x i , k - x ^ i , k ) - - - ( 8 )
wherein x isi,kThe process measurement vector representing the ith lot at time k,then the corresponding result is reconstructed from the PCA model. The calculation result in equation (8) may form an I × 1 vector [ SPE ] at each time1,k,SPE2,k,...,SPEI,k]TAnd the vector approximation is subject to a weighting χ2And distributing to obtain the monitoring control line of SPE at each moment.
And step 9: online process monitoring based on time period model
And (4) on the basis of the time period divided in the step (6), the time period model monitoring system established in the step (7) and the two monitoring statistics obtained in the step (8), the state of a new operation intermittent process such as injection molding and the like can be monitored on line. This step is realized by the following substeps:
(9.1) acquiring and preprocessing new measurement data
During on-line monitoring, new process measurement data x are collectednew(J × 1) (where the subscript new represents a new sample), data preprocessing is first required. And (3) calling the mean value and the standard deviation corresponding to the moment according to the mean value and the standard deviation obtained in the step (2) and according to the indication of the process time to carry out standardization preprocessing on the existing data as shown in a formula (1).
(9.2) calculating New monitoring statistics
After data preprocessing, calling a model P corresponding to the time interval of the new sampling moment according to a PCA sub-time interval model calculated by the formula (6)c(J×Rc) Calculating to obtain principal component score, estimating residual error and its corresponding Hotelling-T2And SPE two monitoring statistical indexes:
tnew T=xnew TPc
x ^ new = P c t new ( 9 )
<math> <mrow> <msubsup> <mi>T</mi> <mi>new</mi> <mn>2</mn> </msubsup> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>new</mi> </msub> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>&OverBar;</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <msup> <msub> <mi>S</mi> <mi>k</mi> </msub> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>new</mi> </msub> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>&OverBar;</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> </mrow> </math>
SPE new = ( x new - x ^ new ) T ( x new - x ^ new )
whereinIs the new process measurement data that is,is T previously obtained from training datakMean vector of (S)kIs XkThe covariance matrix of (2).
(9.3) judging the running state of the process on line
Comparing the two monitoring indexes with respective statistical control limits in real time, and if the two monitoring indexes are both within the statistical control limits, indicating that the process is normal in operation; if more than one monitoring index exceeds the normal control limit, the abnormal condition of the process is indicated. Appropriate fault diagnosis methods, such as contribution graph methods, can then be used to analyze and isolate possible fault variables.
According to the monitoring model established by the sub-period division result, an engineer can obtain the online monitoring result of the sampling data of the new process in real time and compare the corresponding monitoring statistic control line, so that the state of the new running process can be automatically judged, and whether a fault occurs or not can be detected. When 5 continuous monitoring statistics exceed the normal control line, a certain fault occurs in the process operation, and the certain fault represents the First Alarm Time (FAT), and the index can be used for evaluating the sensitivity of fault detection. In the application example of injection molding, there are 49 batches of data for each fault type, and the Mean value (Mean) and the Mean Absolute Deviation (MAD) of the FAT can be calculated for comprehensively evaluating the fault detection performance of the built monitoring model for 49 batches. As shown in table 1, the performance of fault detection based on the method of the present invention was compared with that based on conventional cluster partitioning. Firstly, for a monitoring model based on clustering time interval division, the best monitoring performance is selected from different theta values (as shown in fig. 3) for display. The monitoring model established by the partitioning method displays the alpha value range superior to the clustering method in the obtained monitoring result. In addition, the best and worst monitoring results are selected from the results and compared with the monitoring results based on the clustering method. Generally, the monitoring model established based on the time interval division strategy provided by the invention has the monitoring performance far superior to that of the traditional clustering-based time interval division method, the reliability and the credibility of the actual online process monitoring are greatly improved, the accurate judgment of the process running state by an industrial engineer is facilitated, and the safe and reliable running of the actual production process is ensured.
Table 1: the time interval division method based on the invention is compared with the traditional division method based on the online SPE monitoring comparison result.
In the table, the number of the first and second,*: the online monitoring performance obtained based on the time interval division method is superior to the value range of alpha of the clustering method;
h: selecting the best SPE online monitoring result from different values of the second row alpha for displaying;
o: and selecting the worst SPE online monitoring result from different values of the second row alpha for displaying.
The invention relates to a stepping type sequential time interval automatic division method, which captures the change of process characteristics by analyzing the influence on the model reconstruction precision and the monitoring performance, and successfully proves that the stepping type sequential time interval automatic division method can automatically divide the intermittent production process represented by injection molding and the like into different time intervals by being applied to multi-time interval intermittent production processes such as injection molding and the like, avoids the complicated subsequent processing steps of the traditional method, and greatly improves the precision of the established time interval monitoring model and the accuracy of the subsequent on-line process monitoring. The method comprises the steps of firstly constructing a time slice model, then continuously conducting time slice fusion from the initial time of the process, establishing a sub-period model based on variable expansion in a period of time region to be compared with the time slice model, and analyzing whether the process characteristics of the time slices in the period of time region are similar or not; after the sub-periods comprising similar time slices are determined, iteration is repeated continuously to obtain subsequent sub-periods. The process monitoring system established based on the time interval division result can provide a high-precision online process monitoring result for a technical management department in an actual industrial production field, provides a reliable basis for judging the state of the production process in real time and identifying whether a fault occurs, and finally lays a foundation for safe and reliable operation of production and high-quality pursuit of products.
It should be understood that the present invention is not limited to the injection molding process of the above-described embodiments, and that equivalent modifications or substitutions can be made by those skilled in the art without departing from the spirit of the present invention, and are intended to be included within the scope of the appended claims.

Claims (2)

1. An automatic step-by-step sequential time interval division method, comprising the steps of: step 1: acquiring process analysis data: if an intermittent operation is provided with J measurement variables and K sampling points, each measurement batch can obtain a K multiplied by J matrix, and after the measurement steps of the I batch are repeated, the obtained data can be expressed as a three-dimensional matrixX(I × J × K), wherein the measured variables are state parameters that can be measured during the batch run, including temperature, velocity, pressure, and displacement;
step 2: data preprocessing: will be three-dimensionalMatrix arrayXExpanding according to the collection batch direction, namely arranging the variables on each sampling point in an operation batch according to the time sequence to obtain a two-dimensional matrix X (I multiplied by KJ), and obtaining a matrix X by K time slicesk(I × J), where subscript k is a time indicator;
let two-dimensional matrix XkThe variable at any point in the table is xk,i,jThe variable is normalized by subtracting the mean value and dividing by the standard deviation, wherein the subscript i represents the batch and j represents the variable, and the calculation formula of the normalization process is as follows:
wherein: k is an index of the time slice,is XkMean of any column of the matrix, sk,jIs XkThe standard deviation of the moment for the corresponding column,
and step 3: time slice PCA modeling, which is implemented by the following sub-steps:
(3.1) normalizing each time slice matrix X processed in the step 2k(I × J) performing PCA decomposition to establish a time slice PCA model, wherein the PCA decomposition formula is as follows:
wherein: t is tk,nBeing orthogonal principal component vectors, pk,nThe load vector is orthogonal normalization, r represents different PCA decomposition directions, and superscript T represents the transposition of the matrix; t isk(I X J) representational securityLeaving the score matrix, P, of all principal elementsk(J × J) represents a corresponding load matrix;
(3.2) selecting the number of the principal elements, and restating the formula (3) into the following form:
wherein: r represents different PCA decomposition directions; t isk(I×Rk) And Pk(J×Rk) Representing the load matrix as reserved RkScore matrix and load matrix after an individual principal element, EkIs a residual error matrix; decomposing an original data space into a principal component space and a residual error space by a multidirectional principal component analysis model through a formula (4), wherein the principal component space represents main system process fluctuation information; number of principal elements R reserved herekCan reflect 90% of process fluctuation information in the original process; comprehensively considering the whole process, selecting the R with the most occurrence timeskThe values are unified into the final modeling principal element number R, namely the same principal element number is reserved for all time slice models;
(3.3) calculating SPE indexes corresponding to each batch in each time slice k in the residual error space:
SPEk,i=ek,i Tek,i (5)
wherein the index i denotes the batch in the time slice, ek,iIs the residual column vector of the ith batch corresponding to the k time; obeying chi with weighting coefficient according to SPE values of different batches at the same time2Distribution to determine the control limit Ctr at each time pointkIt reflects the reconstruction capability of the time slice PCA model;
and 4, step 4: determining SPE index control limit based on variable expansion model: from the initial point of the intermittent process, the next time slice is combined with the previous time slice in turn and expanded in a variable mode to obtain a time block Xv,k(Ik × J), wherein the subscript v represents the way the variables expand; carrying out PCA analysis on the new time block matrix to extract a load matrix Pv,k(J × R); calculating SPE value and obeying the authority according to the SPE values of different batches at the same timeChi of the weight coefficient2Distribution to determine the control limit Ctr at each time pointk
And 5: determining a first time period division point k*: comparison of Ctr at each time point within the same time regionkAnd Ctrv,kIf three consecutive samples are found to exhibit Ctrv,k>αCtrkThen, the newly added time slice has a significant influence on the PCA monitoring model and the corresponding monitoring performance of the time block; the time before adding the new time slice is recorded as k*Where α is dependent on CtrkIs a constant, called the mitigation factor, which reflects the extent to which the time block model allows monitoring of the loss of accuracy compared to the time slice model; then k is*A time slice before the time of day may be considered a sub-period;
step 6: the process analyzes the data updates and determines all of the divided periods: according to the time k obtained in step 5*Removing the first sub-period, taking the rest intermittent process data as new input data into the step 5, repeating the step 5-6, and dividing different time periods until no data remains;
and 7: and (3) statistical modeling based on time interval division results: according to the time interval division result in the step 6, time slices in each time interval are combined into a sub-time interval representative modeling data group in a variable expansion mode, Xc(IKcXj), where subscript c is a period indicator; similar process potential characteristics in each time interval can be obtained by pairing Xc(IKcX J) PCA extraction:
Tc=XcPc
wherein, Pc(J×Rc) Is the PCA load matrix of the sub-period, revealing the sub-periodMain direction of fluctuation in time interval, RcRepresenting the number of PCA principal components reserved by the sub-period model; t isc(IKc×Rc) A principal score matrix representing the sub-period; thus, the main fluctuation of the period is through TcPc TThe characterization represents the main fluctuation information in the sub-period; the remainder is the sub-period residual matrix, Ec(IKcxJ) represents noise information within the sub-period; each time slice principal score, Tk(I×Rc) The matrix T can be easily scored from the sub-periodsc(IKc×Rc) According to the corresponding process time, the covariance relation S at each sampling moment can be calculated and obtainedk(ii) a Residual error matrix E of each time slicek(I × J) may also be derived from the sub-period residual matrix Ec(IKcXJ) is correspondingly obtained;
and 8: calculating a real-time monitoring statistical index: according to the principal component score time slice T obtained from the result of the calculation of the formula (6)k(I×Rc) And residual matrix time slice Ek(I × J) two monitoring statistics can be calculated at each time instant: Hotelling-T2Counting indexes and SPE statistics;
Hotelling-T2the statistical indicator is used for measuring the distance of the process variable deviating from the average track under the normal working condition at each sampling moment, and the control limit under the significance level alpha is calculated as follows:
wherein R iscIs the number of principal elements reserved in the PCA model in the time period; t is ti,k(RcX 1) is the principal component score of the ith batch at time k, i.e. the corresponding time slice score matrix Tk(I×Rc) Row i of (1); whileIs Tk(I×Rc) For the mean vector of different batches, the measurement data of each time slice is already measuredCentering to zero mean in data preprocessing, hereIt is a zero vector in nature;
for the residual subspace, the SPE statistics for different batches at each time are calculated as:
wherein x isi,kThe process measurement vector representing the ith lot at time k,then is the corresponding result reconstructed by the PCA model; the calculation result in equation (8) may form an I × 1 vector [ SPE ] at each time1,k,SPE2,k,...,SPEI,k]TAnd the vector approximation is subject to a weighting χ2Distributing to obtain the monitoring control limit of the SPE at each moment;
and step 9: on-line process monitoring based on a time period model: monitoring the state of the new operation intermittent process on line based on the time period divided in the step 6, the time period model monitoring system established in the step 7 and the two monitoring statistics obtained in the step 8; this step is realized by the following substeps:
(9.1) acquiring new measurement data and preprocessing the new measurement data: during on-line monitoring, new process measurement data x are collectednew(J × 1), where the subscript new represents the new sample and J is the same measured variable as in step 1; firstly, data preprocessing is needed; calling the mean value and the standard deviation corresponding to the moment according to the indication of the process time to carry out standardization preprocessing on the acquired process measurement data as shown in a formula (1) according to the mean value and the standard deviation obtained in the step 2;
(9.2) calculating new monitoring statistics: after data preprocessing, the corresponding new process measurement data x is called according to the PCA sub-period model calculated by the formula (6)new(JX 1) model P of the time periodc(J×Rc) Calculating to obtain principal component score, estimating residual error and its corresponding Hotelling-T2And SPE two monitoring statistical indexes:
tnew T=xnew TPc
(9)
wherein,is the new process measurement data that is,is T previously obtained from training datakMean vector of (S)kIs XkThe covariance matrix of (a);
(9.3) judging the process running state on line: comparing the two monitoring indexes with respective statistical control limits in real time, and if the two monitoring indexes are both within the statistical control limits, indicating that the process is normal in operation; if more than one monitoring index exceeds the normal control limit, the abnormal condition of the process is indicated.
2. The automated step-by-step sequential time interval division method according to claim 1, wherein in said step 1, said measurement variables are the following 9: the temperature control device comprises a pressure valve opening, a flow valve opening, a screw stroke, a screw speed, injection pressure, a nozzle temperature, a barrel head temperature, a barrel middle temperature and a barrel tail temperature.
CN201310046432.1A 2013-02-05 2013-02-05 Automatic stepping type ordered time interval dividing method Active CN103116306B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310046432.1A CN103116306B (en) 2013-02-05 2013-02-05 Automatic stepping type ordered time interval dividing method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310046432.1A CN103116306B (en) 2013-02-05 2013-02-05 Automatic stepping type ordered time interval dividing method

Publications (2)

Publication Number Publication Date
CN103116306A CN103116306A (en) 2013-05-22
CN103116306B true CN103116306B (en) 2015-06-17

Family

ID=48414719

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310046432.1A Active CN103116306B (en) 2013-02-05 2013-02-05 Automatic stepping type ordered time interval dividing method

Country Status (1)

Country Link
CN (1) CN103116306B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108664009A (en) * 2017-08-03 2018-10-16 湖州师范学院 Divided stages based on correlation analysis and fault detection method

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103777627B (en) * 2014-01-24 2016-03-30 浙江大学 A kind of batch process on-line monitoring method based on a small amount of batch
CN104298187B (en) * 2014-06-12 2017-03-29 东北大学 Golden hydrometallurgy whole process three-decker process monitoring method
CN104714537B (en) * 2015-01-10 2017-08-04 浙江大学 A kind of failure prediction method based on the relative mutation analysis of joint and autoregression model
CN104699075B (en) * 2015-02-12 2017-04-12 浙江大学 Unequal time period automatic ordered partition-based process monitoring method
CN105353607B (en) * 2015-11-26 2017-10-27 江南大学 A kind of batch process self study dynamic optimization method driven by data difference
EP3420424B1 (en) * 2016-02-22 2024-04-03 Kistler Holding AG Method for carrying out a cyclical production process
CN108803531B (en) * 2018-07-17 2019-10-15 浙江大学 Closed-loop system process monitoring method based on sound feature Cooperative Analysis and orderly Time segments division
CN109932908B (en) * 2019-03-20 2022-03-01 杭州电子科技大学 Multi-directional principal component analysis process monitoring method based on alarm reliability fusion
CN117032996B (en) * 2023-10-09 2023-12-22 湖南中青能科技有限公司 Power metadata management method and system

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102426439A (en) * 2011-04-02 2012-04-25 东北大学 Pierced billet quality forecasting and control method based on data driving
CN102431136A (en) * 2011-09-16 2012-05-02 广州市香港科大***研究院 Multi-phase batch process phase dividing method based on multiway principal component analysis method

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE102005008975A1 (en) * 2005-02-28 2006-08-31 Robert Bosch Gmbh Processor system monitoring method, involves newly starting cyclic permutation of two timers, if one of processes is started, and displaying error signal, if time interval detected by timer exceeds pre-determined maximum time interval

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102426439A (en) * 2011-04-02 2012-04-25 东北大学 Pierced billet quality forecasting and control method based on data driving
CN102431136A (en) * 2011-09-16 2012-05-02 广州市香港科大***研究院 Multi-phase batch process phase dividing method based on multiway principal component analysis method

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于多时段MPCA模型的间歇过程监测方法研究;常玉清等;《自动化学报》;20100930;第36卷(第9期);1312-1320 *
时段划分的多向主元分析间歇过程监测及故障变量追溯;玉姝等;《控制理论与应用》;20110228;第28卷(第2期);150-156 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108664009A (en) * 2017-08-03 2018-10-16 湖州师范学院 Divided stages based on correlation analysis and fault detection method

Also Published As

Publication number Publication date
CN103116306A (en) 2013-05-22

Similar Documents

Publication Publication Date Title
CN103116306B (en) Automatic stepping type ordered time interval dividing method
CN103777627B (en) A kind of batch process on-line monitoring method based on a small amount of batch
CN103336507B (en) Automatically the statistical modeling divided based on the multi-modal collaborative period and on-line monitoring method
EP3809220B1 (en) Method and system for semi-supervised deep anomaly detection for large-scale industrial monitoring systems based on time-series data utilizing digital twin simulation data
CN108803531B (en) Closed-loop system process monitoring method based on sound feature Cooperative Analysis and orderly Time segments division
CN104166787B (en) A kind of aero-engine method for predicting residual useful life based on multistage information fusion
CN103970092B (en) Multi-stage fermentation process fault monitoring method based on self-adaption FCM algorithm
CN105955241B (en) A kind of quality fault localization method based on joint data-driven production process
CN110245460B (en) Intermittent process fault monitoring method based on multi-stage OICA
CN104699075B (en) Unequal time period automatic ordered partition-based process monitoring method
CN104699077B (en) A kind of failure variable partition method based on nested iterations Fei Sheer discriminant analyses
CN108664009B (en) Stage division and fault detection method based on correlation analysis
CN105629958B (en) A kind of batch process method for diagnosing faults based on sub-period MPCA SVM
CN105574587B (en) A kind of online operating mode course monitoring method of plastic injection molding process
CN103838217B (en) A kind of sweat fault monitoring method based on MICA-OCSVM
CN103310095A (en) Intermittent process quality index soft measuring method
CN103473540A (en) Vehicle track incremental modeling and on-line abnormity detection method of intelligent traffic system
CN108334674A (en) A kind of steam turbine high-pressure cylinder method for monitoring operation states based on parameter association intellectual analysis
CN105334823A (en) Supervision-based industrial process fault detection method of linear dynamic system model
CN103926919A (en) Industrial process fault detection method based on wavelet transform and Lasso function
CN106897505A (en) A kind of structure monitoring data exception recognition methods for considering temporal correlation
Zhao et al. Inter-batch-evolution-traced process monitoring based on inter-batch mode division for multiphase batch processes
CN112749623B (en) Processing and feature extraction method and system for high-frequency sensor data of injection molding process
CN109932908B (en) Multi-directional principal component analysis process monitoring method based on alarm reliability fusion
CN104731955A (en) Methods and systems for diagnostic standard establishment and intelligent diagnosis of wind generation set oil monitoring

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant