CN104679844A - Intermittent process batch data synchronizing method based on improved DTW (Dynamic Time Wrapping) algorithm - Google Patents

Abstract

The invention discloses an intermittent process batch data synchronizing method based on an improved DTW (Dynamic Time Wrapping) algorithm. The method comprises the following steps: selecting a plurality of batches of data of the same sampling point number from the batch data of a normal intermittent process, clustering the batch data by adopting a clustering algorithm, and selecting batch data which is closest to a clustering center from each cluster for serving as reference batch data, totaling D reference batch data; performing equal-length synchronization processing on testing batch data needing to be synchronized by taking the D reference batch data as a reference track respectively with the improved DTW algorithm to obtain D synchronization data, and averaging the D synchronization data to obtain a synchronization result of the testing batch data. By adopting the improved DTW algorithm, the search area is limited, and the algorithm efficiency is increased; moreover, the plurality of reference batch data is taken as the reference track, so that partial synchronization caused by a single reference sample is avoided, and the stability of the synchronization result is enhanced.

Description

Based on the batch process lot data synchronous method improving DTW algorithm

Technical field

The invention belongs to batch process field of fault detection, more specifically say, relating to a kind of batch process lot data synchronous method based on improving DTW algorithm.

Background technology

Batch process is the topmost mode of production of modern process industry, due to the dirigibility of itself, is widely used in the production of multi items, high value added product.But, due to the complicacy of its process itself and the impact of other interference, make often there is larger difference between actual motion track and desired trajectory, and finally cause product quality significantly to decline.In addition, this type systematic, once break down, not only can cause personnel and property loss, and environmental pollution is also much serious than other accidents.Therefore, the malfunction monitoring of batch process and prediction of quality are one of focuses of process control field research always.

Traditional batch process fault detection method is as MPCA ((Multilinear Principal Component Analysis, polyteny principal component analysis (PCA)), an important prerequisite of carrying out batch process fault detect is that the data of different batches have identical length.And actual batch process has very strong complicacy, so just cause can not reaching desirable duplication of production between the different batches of same batch process, therefore the length of process data also can not be identical.In multistage batch process, data asynchronous problem performance quite complicated, data asynchrony phenomenon all likely occurs in some or certain several specific stage.For above problem, Lakshminarayanan etc. propose all batches of tracks to expand to the longest, and in namely choosing batch, the longest track is standard, and expand other lot data to standard trajectory length, people is for adding measured value.Certainly also having by the shortest track is standard, and more artificial data points of deleting random length lot data, realize different batches data isometric.These methods all have ignored the feature of local mode, namely just data are become isometric simply by force, probably amplify or lose the local characteristics of raw data.Dynamic time warping (Dynamic Time Wraping, DTW) be thought based on dynamic programming (DP), it is the method for the seasonal effect in time series similarity that a kind of measurement two length are different, be mainly used in template matches, such as speech recognition, gesture identification, data mining and information retrieval etc., thus the method by regular for each for batch process lot data to standard batch sequence, can realize each batch synchronously.

DTW is a kind of method for mode matching proposed based on the thoery of dynamic programming (DP), Time alignment and distance measure calculate and combine by this algorithm, undertaken compressing, expand or change some vector, to obtain the minor increment between two tracks by the similar features searched between two tracks.

For the present invention for batch process, if T (t × N) and R (r × N) is 2 polynary tracks, represent respectively with reference to batch and with reference to batch, wherein t and r is respectively the sampling number of two batches, namely sampling number, the variable number of N for gathering at every turn.DTW uses dynamic programming principle, and non-linearly misplace 2 tracks, homotaxy event, makes wherein each vector of a track and each vector of another track corresponding, to obtain the bee-line between two tracks.If i and j is respectively the coordinate about the time on R (r × N) and T (t × N) track, namely which time sampling is represented, span is respectively 1≤i≤r, 1≤j≤t, R (i × N) represents R (r × N) track i-th sampled data, and T (j × N) represents T (t × N) track jth sampled data.DTW sets up the F* sequence of K point in r × t grid:

F*＝{c(1),c(2),...,c(k),...,c(K)} (1)

Wherein max (t, r)≤K≤t+r, c (k)=[i (k), j (k)], for representing the every bit that i and j mates in grid.

In DTW algorithm, F* sequence can regard the optimal path be in r × t grid making gauged distance between two tracks the shortest as.If d (i, j) represents the Euclidean distance value of R (i × N) and T (j × N):

$d(i,j)=\underset{c=1}{\overset{N}{\mathrm{\Σ}}}(w_c\·{(R(i,c)-T(j,c))}^2)---\left(2\right)$

Wherein, w _cfor the weights of each variable, reflect the relative importance of each measured variable.R (i, c) represents the value of c variable in i-th sampling in R (r × N), and T (j, c) represents the value of c variable in jth time sampling in T (t × N).

Construct a similarity matrix D _a(i.e. Cumulative Distance matrix), the diversity factor be used between description two tracks.D can be obtained by d (i, j) with the algorithm of dynamic programming _a(i, j), the local restriction that application Itakura proposes can obtain following stepping type:

$D_A(i,j)=\min \begin{array}{}\end{array}\left\{\begin{array}{c}{D}_A(i-1,j)+d(i,j)\\ D_A(i-1,j-1)+d(i,j)\\ D_A(i-1,j-2)+d(i,j)\end{array}\right.---\left(3\right)$

D in formula _a(1,1)=d (1,1).

DTW searches for optimal path in t × r rectangular node, and when rectangular area is larger, calculated amount is quite large, is difficult to reach calculate in real time and application on site.

Summary of the invention

The object of the invention is to overcome the deficiencies in the prior art, a kind of batch process lot data synchronous method based on improving DTW algorithm is provided, for multiple with reference to lot data, adopt improvement DTW algorithm to obtain synchrodata respectively, more on average obtain the synchronized result of test lot data.Improve efficiency of algorithm on the one hand, improve the stability of synchronized result on the other hand.

For achieving the above object, the present invention is based on the batch process lot data synchronous method improving DTW algorithm, comprise the following steps:

S1: the sampling number r of preset reference batch, some lot data that sampling number is r are chosen from the lot data of normal batch process, cluster is carried out to the lot data chosen, the number of clusters obtained is designated as D, from each cluster, select the lot data nearest with cluster centre as with reference to lot data, amount to D with reference to lot data;

S2: carry out synchronous test lot data for needing, D that obtains with step S1 respectively with reference to lot data for reference locus, adopt improvement DTW algorithm to carry out isometric synchronous process, obtaining D sampling number is the synchrodata of r, and the concrete grammar of isometric synchronous process comprises the following steps:

S2.1: calculate the distance matrix d between each sampled data in test lot data T (t × N) and reference lot data R (r × N), d (i, j) the Euclidean distance value of R (i × N) and T (j × N) is represented, R (i × N) represents with reference to lot data R (r × N) i-th sampled data, T (j × N) represents test lot data T (t × N) jth sampled data, the span of i and j is respectively 1≤i≤r, 1≤j≤t;

S2.2: initialization cumulant matrix D _amiddle D _a(1,1)=d (1,1), all the other elements are just infinite;

S2.3: make loop variable i=2;

S2.4: calculate this region of search [y _min, y _max], be divided into following three kinds of situations:

As i ∈ [1, X _l], y _min, y _maxcomputing formula be:

y _min＝ai-a+1

y _max＝bi-b+1

As i ∈ [X _l+ 1, X _r], y _min, y _maxcomputing formula be:

y _min＝ai-a+1

y _max＝ai+(b-a)x _l-b+1

As i ∈ [X _r+ 1, r], y _min, y _maxcomputing formula be:

y _min＝bi-(b-a)x _r-a+1

y _max＝ai+(b-a)x _l-b+1

Wherein:

$\left\{\begin{array}{c}{X}_l=\frac{t-\mathrm{ar}+b-1}{b-a}\\ X_r=\frac{\mathrm{br}-t-a+1}{b-a}\end{array}\right.$

X _land X _rall get the integer the most close with result of calculation, a, b are the slope of default path restriction parallelogram adjacent two edges, wherein 0 < a < k, b > k, k=t/r;

S2.5: successively to the region of search [y _min, y _max] in each value j, calculate Cumulative Distance D _a(i, j), computing formula is:

$D_A(i,j)=\min \begin{array}{}\end{array}\left\{\begin{array}{c}{D}_A(i-1,j)+d(i,j)\\ D_A(i-1,j-1)+d(i,j)\\ D_A(i-1,j-2)+d(i,j)\end{array}\right.$

S2.6: judge whether i < r, if so, make i=i+1, return step S2.4, otherwise enter step S2.7;

S2.7: backtracking obtains synchrodata S (r × N), and concrete steps comprise:

S2.7.1: make last data S (r × N)=T (t × N) in synchrodata S (r × N);

S2.7.2: initialization p=r-1, q=t;

S2.7.3: according to Cumulative Distance matrix D _astraight line D _a(p, q), D _a(p, q-1) and D _a(p, q-2), selects wherein minimum value, and the sampled point sequence number of the test lot data of its correspondence is designated as q ', makes S (p × N)=T (q ' × N);

S2.7.4: if p > 1, make p=p-1, q=q ', return step S2.7.3, otherwise terminate search, obtain synchrodata S (r × N);

S3: D the synchrodata obtained by step S2 is averaged, and obtains the synchronized result of these test lot data.

The present invention is based on the batch process lot data synchronous method improving DTW algorithm, first from the lot data of normal batch process, choose some lot data of identical sampling number, the lot data that employing clustering algorithm is selected after carrying out cluster to these lot data from each cluster and cluster centre is nearest, as with reference to lot data, amounts to D with reference to lot data; Synchronous test lot data are carried out for needing, respectively with D with reference to lot data for reference locus, adopt and improve DTW algorithm and carry out isometric synchronous process, obtain D synchrodata, again D synchrodata is averaged, obtains the synchronized result of these test lot data.

The present invention has following beneficial effect:

(1) adopt improvement DTW algorithm, restriction search area, thus reduce calculated amount, improve efficiency of algorithm;

(2) adopt multiple reference lot data as reference track, the mistake avoiding single reference sample to cause is partially synchronous, improves the stability of synchronized result.

Accompanying drawing explanation

Fig. 1 is DTW algorithm search path profile;

Fig. 2 is the process flow diagram that the present invention is based on the batch process lot data synchronous method improving DTW algorithm;

Fig. 3 improves the process flow diagram that DTW algorithm carries out isometric synchronous process;

Fig. 4 is synchrodata search routine figure;

Fig. 5 is the DTW path profiles of test lot data to 3 reference locus.

Embodiment

Below in conjunction with accompanying drawing, the specific embodiment of the present invention is described, so that those skilled in the art understands the present invention better.Requiring particular attention is that, in the following description, when perhaps the detailed description of known function and design can desalinate main contents of the present invention, these are described in and will be left in the basket here.

Embodiment

In order to better the present invention is described, first the DTW algorithm of improvement of the present invention is briefly described.

Improving DTW algorithm is on DTW algorithm basis, in order to reduce calculated amount, prevent from carrying out path overbend when DTW searches for, limits, thus improve the search speed in path to the direction of search.Fig. 1 is DTW algorithm search path profile.As shown in Figure 1, DTW algorithm carries out search optimal path in (r × t) grid rectangle, can search for along 3 directions according to any point in formula (3) known grid, be two net points directly over the neighbor mesh points of right, right neighbor mesh points respectively, in figure, namely the signal of black heavy line is optimal path along above-mentioned rule search.DTW algorithm carries out the search in above-mentioned three directions to each point in (r × t) grid rectangle and calculates Cumulative Distance matrix, and when r and t is larger, calculated amount can become very large, is unfavorable for practical application.Path is limited in the parallelogram shown in Fig. 1 dotted line by the improvement DTW algorithm that the present invention adopts, search area is significantly less than original rectangular search area, namely to the net point outside parallelogram, thinks that path is unreachable, not within computer capacity, thus greatly can reduce calculated amount.

In route matching process, the slope of bending is defined owing to improving DTW algorithm, many lattice points is actually and does not reach, therefore the frame matching distance that the lattice point outside parallelogram is corresponding is not calculative, also there is no need to preserve all frame matching distance matrixes and Cumulative Distance matrix, and the matching primitives on each row lattice point has only used three grids of previous column, make full use of these two features and can reduce calculated amount and save storage space.

As shown in Figure 1, the dynamic bending of reality is divided into 3 sections: [1, X _l], [X _l+ 1, X _r], [X _r+ 1, r].X _l, X _rvalue can obtain according to slope a, the b of parallelogram adjacent two edges preset, wherein 0 < a < k, b > k, k=t/r.X _l, X _rcomputing formula be:

$\left\{\begin{array}{c}{X}_l=\frac{t-\mathrm{ar}+b-1}{b-a}\\ X_r=\frac{\mathrm{br}-t-a+1}{b-a}\end{array}\right.---\left(4\right)$

X _land X _rall get the integer the most close with result of calculation, round up by result, when Searching point is when parallelogram is not inner, thinks that difference of them is too large, cannot dynamic time warping be carried out.Above-mentioned condition makes searching route become by original rectangular area the parallelogram region that an area is less than rectangle.Because matching process follows the principle of shortest path, can think, shortest path too tilts scarcely, therefore, in order to reduce calculated amount and improve matching speed, can Constrain Searching path slope within the specific limits.The span arranging slope in the present embodiment is, if k < 1,0.2k≤a < is k, and k < b≤5k; If k > 1,0.5k≤a < is k, k < b≤2k.

Fig. 2 is the process flow diagram that the present invention is based on the batch process lot data synchronous method improving DTW algorithm.As shown in Figure 2, the present invention is based on the batch process lot data synchronous method improving DTW algorithm to comprise the following steps:

S201: obtain D with reference to lot data:

The sampling number r of preset reference batch, some lot data that sampling number is r are chosen from the lot data of normal batch process, cluster is carried out to the lot data chosen, the number of clusters obtained is designated as D, from each cluster, select the lot data nearest with cluster centre as with reference to lot data, amount to D with reference to lot data.

S202: order is with reference to lot data sequence number λ=1;

S203: improve DTW algorithm synchronously isometric:

Using λ reference lot data as with reference to track, adopting improvement DTW algorithm to carry out isometric synchronous process to needing to carry out synchronous test lot data, obtaining λ synchrodata.Fig. 3 improves the process flow diagram that DTW algorithm carries out isometric synchronous process.As shown in Figure 3, improving the concrete steps that DTW algorithm carries out isometric synchronous process is:

S301: calculate distance matrix d:

Distance matrix d is tried to achieve according to formula (2), distance matrix d is the matrix of r × t, t represents the sampling number of test lot data, d (i, j) the Euclidean distance value of R (i × N) and T (j × N) is represented, i and j span is respectively 1≤i≤r, 1≤j≤t.

S302: initialization Cumulative Distance matrix D _a:

Initialization cumulant matrix D _amiddle D _a(1,1)=d (1,1), all the other elements are just infinite.

S303: make loop variable i=2;

S304: calculate this region of search [y _min, y _max]:

According to the explanation improving DTW algorithm, knownly can calculate the region of search corresponding to each i according to parallelogram four edges, its computing formula is divided into following three kinds of situations:

As i ∈ [1, X _l], y _min, y _maxcomputing formula be:

y _min＝ai-a+1

(5)

y _max＝bi-b+1

As i ∈ [X _l+ 1, X _r], y _min, y _maxcomputing formula be:

y _min＝ai-a+1

(6)

y _max＝ai+(b-a)x _l-b+1

As i ∈ [X _r+ 1, r], y _min, y _maxcomputing formula be:

y _min＝b _i-(b-a)x _r-a+1

(7)

y _max＝ai+(b-a)xl-b+1

S305: make loop variable j=y _min;

S306: calculate Cumulative Distance D _a(i, j):

$D_A(i,j)=\min \begin{array}{}\end{array}\left\{\begin{array}{c}{D}_A(i-1,j)+d(i,j)\\ D_A(i-1,j-1)+d(i,j)\\ D_A(i-1,j-2)+d(i,j)\end{array}\right.---\left(8\right)$

S307: judge whether j < y _max, if so, enter step S308, otherwise enter step S309.

S308: make j=j+1, returns step S306.

S309: judge whether i < r, if so, enter step S310, otherwise enter step S311.

S310: make i=i+1, returns step S304.

S311: backtracking obtains synchrodata S (r × N):

After Cumulative Distance matrix computations completes, now need back to search for from last element of matrix, obtain isometric synchrodata.Fig. 4 is synchrodata search routine figure.As shown in Figure 4, the concrete steps of synchrodata search comprise:

S401: make last data S (r × N)=T (t × N) in synchrodata S (r × N).

S402: initialization p=r-1, q=t.

S403: search obtains S (p × N):

Search procedure is, according to Cumulative Distance matrix D _aalong three direction straight line D of formula (3) _a(p, q), D _a(p, q-1) and D _a(p, q-2), selects wherein minimum value, and the sampled point sequence number (i.e. ordinate) of the test lot data of its correspondence is designated as q ', makes S (p × N)=T (q ' × N).

S404: judge whether p > 1, if so, enters step S405, otherwise terminates search.

S405: make p=p-1, q=q ', return step S403.

The synchrodata sequence S (r × N) of one and reference locus equal length can be generated from aforesaid operations, achieve the isometric synchronous process of lot data.

S204: judge whether λ < D, if so, enter step S205, otherwise enter step S206.

S205: make λ=λ+1, returns step S203.

S206: the D an obtained synchrodata is averaged, obtains the synchronized result of these test lot data.

In order to beneficial effect of the present invention is described, apply the present invention in the fault detect of typical TE chemical process.This TE process has 41 process variable, about each batch of duration be 1060 chronomeres.Select from the lot data of normal batch process sampling number be 1060 several lot data, 3 classes are divided into by clustering algorithm, from each cluster, select the lot data nearest with cluster centre as with reference to lot data, amount to 3 with reference to lot data.With the test lot data instance that length is 1080, carry out synchronously as with reference to track with reference to lot data with 3 respectively.Fig. 5 is the DTW path profiles of test lot data to 3 reference locus.As shown in Figure 5, a DTW path is respectively had for each reference locus, wherein two paths (two solid line) are close to and overlap, and Article 3 path (dotted line) differs larger with two other path, may make a world of difference in the DTW path of known same sample to different reference sample, thus add in the present invention to different ask the method for DTW synchronized result weighted fitting again to remove to a great extent respectively with reference to lot data mistake that single reference sample causes is partially synchronous, have good stability.

Although be described the illustrative embodiment of the present invention above; so that those skilled in the art understand the present invention; but should be clear; the invention is not restricted to the scope of embodiment; to those skilled in the art; as long as various change to limit and in the spirit and scope of the present invention determined, these changes are apparent, and all innovation and creation utilizing the present invention to conceive are all at the row of protection in appended claim.

Claims (3)

1., based on the batch process lot data synchronous method improving DTW algorithm, it is characterized in that, comprise the following steps:

S1: the sampling number r of preset reference batch, some lot data that sampling number is r are chosen from the lot data of normal batch process, cluster is carried out to the lot data chosen, the number of clusters obtained is designated as D, from each cluster, select the lot data nearest with cluster centre as with reference to lot data, amount to D with reference to lot data;

S2: carry out synchronous test lot batch secondary data for needing, D that obtains with step S1 respectively with reference to lot data for reference locus, improvement DTW algorithm is adopted to carry out isometric synchronous process, obtaining D sampling number is the synchrodata of r, and the concrete grammar of isometric synchronous process comprises the following steps:

S2.1: calculate the distance matrix d between each sampled data in test lot data T (t × N) and reference lot data R (r × N), d (i, j) the Euclidean distance value of R (i × N) and T (j × N) is represented, R (i × N) represents with reference to lot data R (r × N) i-th sampled data, T (j × N) represents test lot data T (t × N) jth sampled data, the span of i and j is respectively 1≤i≤r, 1≤j≤t;

S2.2: initialization cumulant matrix D _amiddle D _a(1,1)=d (1,1), all the other elements are just infinite;

S2.3: make loop variable i=2;

S2.4: calculate this region of search [y _min, y _max], be divided into following three kinds of situations:

As i ∈ [1, X _l], y _min, y _maxcomputing formula be:

y _min＝ai-a+1

y _max＝bi-b+1

As i ∈ [X _l+ 1, X _r], y _min, y _maxcomputing formula be:

y _min＝ai-a+1

y _max＝ai+(b-a)x _l-b+1

As i ∈ [X _r+ 1, r], y _min, y _maxcomputing formula be:

y _min＝bi-(b-a)x _r-a+1

y _max＝ai+(b-a)x _l-b+1

Wherein:

$\left\{\begin{array}{c}{X}_l=\frac{t-\mathrm{ar}+b-1}{b-a}\\ X_r=\frac{\mathrm{br}-t-a+1}{b-a}\end{array}\right.$

X _land X _rall get the integer the most close with result of calculation, a, b are the slope of default path restriction parallelogram adjacent two edges, wherein 0 < a < k, b > k, k=t/r;

S2.5: successively to the region of search [y _min, y _max] in each value j, calculate Cumulative Distance D _a(i, j), computing formula is:

$D_A(i,j)=\min \begin{array}{}\left\{\begin{array}{c}{D}_A(i-1,j)+d(i,j)\\ D_A(i-1,j-1)+d(i,j)\\ D_A(i-1,j-2)+d(i,j)\end{array}\right.\end{array}$

S2.6: judge whether i < r, if so, make i=i+1, return step S2.4, otherwise enter step S2.7;

S2.7: backtracking obtains synchrodata S (r × N), and concrete steps comprise:

S2.7.1: make last data S (r × N)=T (t × N) in synchrodata S (r × N);

S2.7.2: initialization p=r-1, q=t;

S2.7.3: according to Cumulative Distance matrix D _astraight line D _a(p, q), D _a(p, q-1) and D _a(p, q-2), selects wherein minimum value, and the sampled point sequence number of the test lot data of its correspondence is designated as q ', makes S (p × N)=T (q ' × N);

S2.7.4: if p > 1, make p=p-1, q=q ', return step S2.7.3, otherwise terminate search, obtain synchrodata S (r × N);

S3: D the synchrodata obtained by step S2 is averaged, and obtains the synchronized result of these test lot data.

2. batch process lot data synchronous method according to claim 1, is characterized in that, the clustering method in described step S1 is K means clustering method.

3. batch process lot data synchronous method according to claim 1, it is characterized in that, in described step S2.4, if k < 1, the span of slope a, b is respectively 0.2k≤a < k, k < b≤5k; If k > 1, the span of slope a, b is respectively 0.5k≤a < k, k < b≤2k.