252
Ind. Eng. Chem. Res. 2001, 40, 252-260
Extraction of Event Times in Batch Profiles for Time Synchronization and Quality Predictions Nitin Kaistha and Charles F. Moore* Department of Chemical Engineering, University of Tennessee, Knoxville, Tennessee 37996
A technique for the extraction of the time of occurrence of consistent features in batch profiles is presented. The event times are needed for the synchronization of the batch trajectories for meaningful analysis of the historical data using multivariate statistical techniques such as principal component analysis and partial least squares and also for studying the correlation of the event times with the final batch product quality. The event detection technique is based on the matched filter approach. A normalization strategy that compensates for the measurement axis variability in the profile snapshots is enunciated for accurate extraction of event times. The methodology is applied to profiles generated from a simulated methyl methacrylate polymerization batch reactor and an industrial batch polymerization reactor. The use of the event times for time synchronization is illustrated. The importance of studying the correlation of the event times with the final batch quality is also emphasized. Introduction Batch processes are used in the chemical industry for the production of various chemicals and are especially preferred when the product volume is small, market demand is uncertain, and flexibility is needed to produce various grades of the product as dictated by the market.1 Examples of batch manufacturing include specialty chemicals such as pharmaceuticals and biochemicals and commodity chemicals such as polymers. Almost all of the unit operations typically encountered in the continuous process industry have their batch counterparts (e.g., batch distillation columns, batch reactors, batch mixers, batch crystallization units, etc.). The batch cycle typically consists of a sequence of steps, also known as phases in a batch.1 For example, the cycle for a batch reactor may consist of (a) a charging phase, where specified amounts of ingredients are added to the batch vessel; (b) a preheat phase, where the vessel ingredients are heated to the reaction temperature; (c) a reaction phase, where the catalyst is added, signaling the beginning of the reaction; and (d) a product discharge phase, where the product is discharged from the vessel after completion of the reaction and prepared for the next batch. Once the batch is complete, a sample of the product is assessed for quality in the quality control laboratory as online product quality measurements are usually not available. The unavailability of product quality until after the quality tests are complete can lead to a sequence of offspecification batches because disturbances such as impurities in the raw material are not compensated for. Indirect information on the disturbances can however be obtained from the various online measurements on easily measurable variables such as temperatures, flows, etc. There is therefore considerable interest in the empirical modeling of these measurement profiles for batch monitoring and quality correlation purposes. Vast amounts of historical data on the measurement profiles, for both successful and unsuccessful batches, * To whom correspondence should be sent. Phone: (865) 974 2421. Fax: (865) 974 7076. E-mail:
[email protected].
are usually available in a historical database. Multiway methods such as multiway principal component analysis (MPCA),2,3 multiway partial least squares (MPLS),4 and, more recently, Tucker models5 have been used to empirically model the past successful batch operation for statistical process control (SPC) monitoring and quality correlation purposes. An implicit assumption in the application of all of these methods is that the profiles are already time-scaled so that all of the significant events (usually phases) in a batch are synchronized, resulting in equal duration batches. The time synchronization itself is a very crucial step in the meaningful application of these multivariate methods. To accomplish the synchronization, the use of monotonic measurements such as reaction conversion, also called indicator variables, has been suggested.2-4 It has been applied to an industrial emulsion batch polymerization process.6 As an alternative, dynamic time warping (DTW),7-9 a technique from speech recognition, has recently been applied to batch profiles.10 DTW was earlier applied for phase detection to aid in the supervision of bioprocesses.11 DTW takes into consideration the multivariate nature of all of the available measurements and extracts a nonlinear warp function that maps the time indices of a batch to those of a reference. However, it is computationally intensive and transforms the time axis variability to the warp function, which in turn must be properly characterized for a meaningful analysis. The warp function is similar to the indicator variable, so that in the reported batch application for time synchronization,10 more than 80% weight is given to an indicator variable. The techniques thus rely heavily on the availability of such measurements. Indicator variables usually require online analytical measurements, which are as yet quite uncommon in industrial batch processes. Even when available, their applicability is limited to only a single phase in the batch (usually the reaction phase) and not for the various other phases (charge, preheat, cooldown, etc.). The beginning and/or end of the various phases is usually characterized by features that occur in one or more of the profiles. Examples include a step in the inlet
10.1021/ie990937c CCC: $20.00 © 2001 American Chemical Society Published on Web 12/02/2000
Ind. Eng. Chem. Res., Vol. 40, No. 1, 2001 253
flow profile indicating the beginning of the loading phase or a ramp down in the level indicating the beginning of the product discharge phase. Events such as steps, peaks, valleys, etc., that consistently occur from batch to batch indicate similar underlying dynamic phenomena, and a tool is needed for the extraction of the event times. The times are critical for trajectory synchronization for a meaningful subsequent multivariate analysis such as MPCA or MPLS. More importantly, the time of occurrence of the consistent events characterizes the variability along the time axis in the profiles and provides useful information on the health of the process, making them good candidates for the application of simple SPC monitoring charts. Finally, the event times are vital from the point of view of the correlation that they may have with the final product quality. To the best of our knowledge, all previous literature on batch data analysis treats time synchronization only as a necessary pretreatment step and focuses solely on the subsequent application of the multivariate methods such as MPCA or MPLS to the synchronized profiles. The authors strongly feel that much useful information is embedded in the time axis variability and proper characterization of such variation, as demonstrated here, is probably the most important step in the overall analysis. Given the obvious importance of the time of occurrence of consistent events in batch profiles, this paper develops a methodology based on template matching for extracting the event times. The approach is based on the matched filter12 that uses convolution of the profile with a feature template and extracts the event times from the maximum in the output. An enhancement to the matched filter approach involving normalization of the profiles, for accurate event detection, is presented. The proposed methodology is applied for event detection in profiles from a simulated polymerization batch reactor and an industrial polymerization reactor. The paper demonstrates the use of the extracted event times for time synchronization and also for modeling the final product quality. A brief discussion of the various engineering issues and the conclusions from the work is presented at the end. Matched Filter The detection of characteristic features in a signal can be accomplished using digital filters. A filter convolves the signal s with the filter coefficient a to compute the correlation function r as follows: n
r(k) )
a(k - i)s(i) ∑ i)1
(1)
Convolution flips a and translates it over s, calculating the overlap at each translation point k. A matched filter typically uses the reversed (flipped) feature f* as its filter coefficients. If p denotes the I × 1 profile vector and f the J × 1 feature template vector, the matched filter convolves f*, the flipped f, with p so that
r(k) ) fkTpk
(2)
where the subscript k is used to highlight that the vectors in eq 2 are of varying length depending on k, I, and J. Assuming that I > J, i.e., the feature is shorter than the profile, we have
Figure 1. Zones in translation of a feature over the profile. Solid line: profile time axis. Box: feature time axis. Relevant time indices are indicated.
Figure 2. Matched filter output showing no maximum corresponding to the occurrence of the feature in the profile. The feature template is mean-centered for offset suppression.
pk ) p(1:k)
if k < J
p(k - J + 1:k)
if k g J and k e I
p(k - J + 1:I)
if k > I
fk ) f(J - k + 1:J)
if k < J
f(1:J)
if k g J and k e I
f(1:J + I - k)
if k > I
(3a)
(3b)
The three parts of the equation are pictorially depicted in Figure 1 and are seen to correspond to the moving in, the complete overlap, and the moving out of the feature upon translation. Peaks in the filter output r correspond to the occurrence of the feature in the profile. The matched filter has been used for phase detection and classification in biotechnological processes.13 For simple situations such as isolated peaks in a flat profile, the matched filter approach is adequate. In more complex situations, convolution of the feature of interest with the profile may not result in the filter output showing a maximum corresponding to the time of occurrence of the feature. This is illustrated in Figure 2, which plots the profile, the feature template, and the matched filter output. No maximum is seen around the time of occurrence of the feature with the filter coefficients mean-centered and normalized for offset suppression.13
254
Ind. Eng. Chem. Res., Vol. 40, No. 1, 2001
Furthermore, even if the filter output shows a maximum corresponding to the occurrence of the feature, the extracted event times may not be accurate because of the inherent variability in the rate of occurrence of the feature. As an example, the feature in Figure 2 can be described in terms of the profile decreasing and saturating at a value. The rate of settling however varies from batch to batch, so that the accuracy of the extracted event times is adversely affected. The traditional matched filter approach must therefore be suitably modified. An approach based on normalization at each time point k, especially suitable for event detection in batch profiles, is discussed in the next section. Normalized Matching The traditional matched filter output is the convolution of a pretreated signal with the flipped feature template. To guarantee the occurrence of a maximum in the output corresponding to the feature and also to compensate for the inherent variability in the shape of the feature, normalization at each translation step k can be performed. More specifically, consider the normalized feature fn where
fn ) (f - mf)/sf
(4a)
mf ) mean(f)
(4b)
sf ) ||f - mf||
(4c)
Figure 3. Schematic of event detection methodology.
with
At each translation point k, the windowed profile pk is normalized to obtain pnk as
pnk ) [pk - mp(k)]/sp(k)
(5a)
mp(k) ) mean(pk)
(5b)
sp(k) ) ||pk - mp(k)||
(5c)
where
The output at the kth translation point is thus
r(k) ) fnkTpnk
(6)
where fnk is the windowed part of the feature with indices as in eq 3b. Using such a normalization scheme, overlap values in r that are closest to 1 correspond to the time of occurrence of the feature. More significantly, the normalization compensates for the inherent variability in the rate of occurrence of the feature in the various batches. This leads to accurate extraction of the event times. The approach uses the matched filter philosophy but differs in that each translation step involves normalization of the windowed profile pk. The traditional matched filter generally uses a single signal pretreatment for the whole profile p. Note that multiple maxima in the output are possible and the maximum corresponding to the feature of interest must be distinguished from the other maxima. The mean mp and the normalization factor sp that are used to pretreat the windowed profile can be used effectively for this purpose. For all of the batch profiles, the feature of interest occurs in the same way as that
Figure 4. Normalized matching output for profile and feature in Figure 1.
of the template with nominal variability around the template. Thus, elements in the normalization factors mp and/or sp in the vicinity of the time of occurrence of the feature should be quite similar to those of the normalization factors mf and sf used for normalizing the feature template to zero mean and unit magnitude. In other words, the maximum of interest is the one for which the windowed profile normalization factors are similar to the feature normalization factors. Appropriate lower and upper bounds on mp and/or sp can be used for the purpose. Also, a feature generally occurs within a time snapshot for all of the batches. For example, the catalyst may be added between 30 and 90 min after the start of a batch for all of the batches in the historical database. The maximum of interest in r that corresponds to the event would therefore be within the 30-90 min snapshot. This can be exploited in conjunction with the bounds on the normalization factors for effectively distinguishing the maximum of interest from the rest. The complete feature time extraction methodology is depicted in Figure 3. It consists of two modules: the normalized matching that performs the translation and the classification of the maxima based on the normalization factors and the snapshot of interest. As an illustration, Figure 4 plots the output for normalized matching for the example profile used earlier. The output is set to 0 whenever division by 0 occurs in eq 5a. Note that the maximum around time 40 is the greatest in magnitude and also corresponds to the occurrence of the feature of interest in the profile. Other maxima also occur but are of lower intensity. The event time can then be easily extracted from the time index of the relevant maximum in r.
Ind. Eng. Chem. Res., Vol. 40, No. 1, 2001 255
Figure 5. Profiles for the simulated batch polymerization reactor (set 1).
Data Sets Batch profile data from a simulated polymerization reactor and an industrial polymerization reactor are used to demonstrate the event detection methodology. These data sets are briefly described below. A. Simulated Methyl Methacrylate Polymerization Reactor. The system modeled is a reactor with a cooling jacket for polymerizing methyl methacrylate in toluene solvent using azobisisobutyronitrile as the initiator. A proportional-integral-derivative controller is used to control the reactor temperature by manipulating the heat input to the jacket. The complete data set consists of 50 batches with four measurements per batch. The measurements are (1) reactor temperature, (2) jacket temperature, (3) coolant flow, and (4) electric heater power. The process is modeled using dynamic material and energy balance equations on the reactor and the cooling system. The reactor size, reactor and jacket holdups, and other physical properties have been taken from Souroush and Kravaris.14 The reaction mechanism is modeled as a free-radical polymerization as described in Appendix 1 where the kinetic parameters have been taken from Baillagou and Soong.15,16 For industrial relevance, a hypothetical impurity (not measurable) capable of significantly affecting the polymer quality by increasing the dead-chain formation is also included in the reaction mechanism. The polymerization kinetics are characterized by a strong gel effect at high monomer conversion. Other simulation details can be found in Kaistha.17 Autocorrelated noise is added to the various measurements in the simulation. A total of N ) 50 batches are simulated with nominal variability in the process parameters such as the jacket heat-transfer coefficient, ambient temperature, impurities in the monomer, initial charge to the reactor, etc. The profiles generated are plotted in Figure 5. This data set is subsequently referred to as set 1. B. Industrial Polymerization Batch Reactor. This data set, subsequently referred to as set 2, consists
Table 1. Event Descriptions for Sets 1 and 2a feature no.
set 1
set 2
1 2 3 4 5 6 7 8
reach minimum (2) settle to 0 (3) sudden rise up from 0 (3) slow rise up from 0 (4) reach maximum (4) sudden drop to 0 (4)
begin ramp down (1) settle to 0 (1) rise up from 0 (1) sudden drop (2) settle at value (3) sudden drop (3) sudden drop to 0 (3) sudden drop to 0 (4)
a Values in parentheses indicate the variable no. associated with the feature.
of profiles from an industrial polymerization batch reactor. Profiles from a total of 44 batches are available with four measurements per batch. These are scaled to be approximately between 0 and 1 and are plotted in Figure 6. Both of the batch data sets show batch-to-batch time misalignment in the various features in the profiles. For set 1, mild misalignment occurs toward the end of the batch in the gel effect period, while for set 2, severe misalignment in the features is seen. The next section demonstrates the application of the normalized matching scheme to the data sets for event time extraction. Event Detection The time of occurrence of features that consistently occur in the two data sets is extracted using the normalized matching approach described earlier. The feature templates (not normalized) superimposed on a profile are plotted in Figures 7 and 8 for sets 1 and 2, respectively. Numeric identifiers for the features are also shown in the figures. To quantify the accuracy of event detection, the event times are obtained manually using the description in Table 1 for the two data sets. If t denotes the calculated event times using the normalized matching approach and tr the actual feature times, where both t and tr are N × 1 column vectors, the root-mean-square error (RMSE)
RMSE ) [(t - tr)T‚(t - tr)/N]0.5
(7)
256
Ind. Eng. Chem. Res., Vol. 40, No. 1, 2001
Figure 6. Profiles from an industrial batch polymerization reactor (set 2).
Figure 7. Feature templates for set 1 used for event detection. The raw templates are mean-centered and normalized to unit norm before translation. Table 2. RMSE Values for Extracted Event Times for Sets 1 and 2 feature no.
set 1
set 2
feature no.
set 1
set 2
1 2 3 4
0.6851 0.8330 0.4738 0.5533
0 0.3492 0 0
5 6 7 8
0.6389 0.5714
0.4939 0.2209 0 0.7325
can be used to quantify the average error in the event detection. In the above, N is the total number of batches. The RMSEs for the different features in sets 1 and 2 are reported in Table 2. The low RMSE values (
