Article pubs.acs.org/IECR
Discussion on Time Difference Models and Intervals of Time Difference for Application of Soft Sensors Hiromasa Kaneko and Kimito Funatsu* Department of Chemical System Engineering, The University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo 113-8656, Japan ABSTRACT: In chemical plants, soft sensors are widely used to estimate process variables that are difficult to measure online. The predictive accuracy of soft sensors decreases over time because of changes in the state of chemical plants, and soft sensor models based on time difference (TD) have been constructed. However, many details of models based on TD remain to be clarified. In this study, TD models are discussed in terms of noise in data, autocorrelation in process variables, predictive accuracy, and so on. We theoretically clarify and formulate the differences of predictive accuracy between normal models and TD models and the effects of noise, autocorrelation, TD intervals, and so on on the predictive accuracy. The relationships and the formulas were verified by analyzing simulation data. Furthermore, we analyzed dynamic simulation data and real industrial data and confirmed that the predictive accuracy of TD models increased when TD intervals were optimized.
1. INTRODUCTION Soft sensors have been widely used to estimate process variables that are difficult to measure online.1,2 In soft sensors, an inferential model is constructed between those variables that are easy to measure online and those that are not, and an objective variable or a dependent variable, y, is then estimated using that model. Through the use of soft sensors, the values of y can be estimated with a high degree of accuracy. Their practical use, however, involves some difficulties. One crucial difficulty is that their predictive accuracy gradually decreases as a result of changes in the state of chemical plants (e.g., reduction in catalyst performance) and sensor and process drift, among others. To reduce the degradation of the soft sensor model, the updating of regression models,3,4 just-in-time (JIT) modeling,5,6 and ensemble modeling7,8 have been proposed.9 While many excellent results based on the use of these methods have been reported, problems remain for the introduction of soft sensors in practice.10−12 First, if soft sensor models are reconstructed with the inclusion of abnormal data, their predictive ability can deteriorate.4 Although such abnormal data must be detected with high accuracy in real time, under present circumstances, it is difficult to accurately detect all abnormal data. Second, reconstructed models are often effective for making predictions over a narrow data range only.13 Subsequently, when rapid variations in the process variables occur, these models cannot predict the resulting variations in data with a high degree of accuracy. Third, if a soft sensor model is reconstructed, the parameters of the model, for example, the regression coefficients in linear regression modeling, can change dramatically in some cases. Without the operators’ understanding of a soft sensor model, the model cannot be practically applied. When soft sensor models are reconstructed, operators check the parameters of the models to ensure that they will be safe during operation. This is time-consuming and laborintensive because tens of soft sensors can be used in a plant.14 Fourth, the data used to reconstruct soft sensor models are also affected by drift. In the construction of the model, data must be © 2013 American Chemical Society
selected from a database that includes both data affected by the drift and data after correction of the drift. The use of time difference (TD) models has been proposed to solve these problems. In this approach, soft sensor models are constructed on the basis of the TD of explanatory variables or independent variables, X, and that of y (without reconstruction of the model) to reduce the effects of deterioration of model accuracy with time as a result of drift and gradual changes in the state of plants.10 In soft sensor modeling, models that are not affected by these temporal changes can be constructed not by using the values of process variables, but by using the TD. The problems inherent in model reconstruction as described above can be avoided because TD models do not have to be reconstructed; the data are represented as the TD, and the TD is not affected by drift. By analyzing industrial data, we showed previously that the TD model can maintain its predictive accuracy for a period of 3 years, without the model being reconstructed.10 TD models have been applied in cases with nonlinear process variables11 and process variations.12 However, the predictive accuracy of TD models is sometimes lower than that of updated models and JIT models,10,15 and many details relating to TD models remain to be clarified. For example, the TD interval applied is usually the smallest practicable TD interval, i.e., a measurement interval of y, because TD models with short intervals tend to have higher prediction accuracy. However, if the TD interval is too small, noise in the process variables will be amplified and the effect of noise on the TD model can be larger than that of significant variations in the process. The shortest TD interval is not necessarily the optimal TD interval. However, if the TD interval is too long, the TD of process variables cannot Received: Revised: Accepted: Published: 1322
September 22, 2012 December 23, 2012 January 3, 2013 January 3, 2013 dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
2.2. Elimination of Terms That Are Linear Functions of Time by Using TD. We will show that terms that are linear functions of time can be eliminated by TD using an example of a process variable, x ∈ Rm×1. First, x related to time t and t−i is separated into a term that is a linear function of time, xL ∈ Rm×1, and a term that is a nonlinear function of time, xNL ∈ Rm×1, as follows:
represent short-term variations. There are no accepted methods to determine the optimal TD interval. In this study, TD models are discussed in terms of the amount of noise in data, autocorrelation in process variables, and degree of model accuracy, among others. We theoretically clarify and formulate the difference of predictive accuracy between normal models and TD models. This approach enables the comparison of TD models and other models and facilitates consideration of an optimal TD interval and the construction of predictive TD models. The relationships and the formulas of the TD model were verified by analyzing numerical simulation data. We analyzed data obtained by dynamic simulation of an existing full-scale depropanizer distillation column and examined various TD intervals and the predictive ability of TD models. Furthermore, by analyzing two sets of data from an operational industrial plant, the increase of predictive accuracy of TD models was confirmed by changing the TD intervals of training data and test data.
(1)
Δy(t ) = y(t ) − y(t − i)
(2)
Δypred (t ′) = f (Δx(t ′))
(5)
= {xL(t ) + xNL(t )} − {xL(t − i) + xNL(t − i)} = ΔxNL(t ) + xL(t ) − xL(t − i) (10)
By using eq 9, Δx(t) is finally expressed as follows: Δx(t ) =ΔxNL(t ) + {p t(t ) + q − (p t(t − i) + q)} = ΔxNL(t ) + pi (11)
because the interval between t(t) and t(t−i) is a constant, i. pi is a constant and can be eliminated by a preprocessing method such as centering. Even when the interval is not constant, any terms depending on t do not remain, and thus, the xL term can be eliminated by using TD. The formulas above are applicable not only to x but also to y. Numerical experiments investigating TD models and terms that are linear functions of time are given in the work of Kaneko and Funatsu.10 By constructing TD models, the effects of deterioration of model accuracy with time, such as those caused by the drift and gradual changes in the state of plants, can be accounted for because data are represented as TDs that are not affected by these factors. 2.3. SN Ratios of a Process Variable and TD of a Process Variable. In time series analysis, the SN ratio of a time-differenced variable is affected by the autocorrelation of the variable.16 In this section, the SN ratios of a process variable and the TD of the process variable are investigated. When x is assumed to be a stationary process, we consider the SN ratios of the x-variable and the TD of the x-variable. First, x(t) is separated into significant variation, xV(t) ∈ Rm×1, and white noise, xE(t) ∈ Rm×1, as follows: x(t ) = xV (t ) + xE(t )
(12)
The SN ratio of x(t), ηx, is calculated as
ypred(t′) can be calculated as ypred (t ′) = Δypred (t ′) + y(t ′ − i)
(9)
Δx(t ) =x(t ) − x(t − i)
where f is a regression model and e ∈ Rm×1 is a vector of calculation errors. In terms of prediction, the constructed model, f, predicts the TD of y(t′), Δy(t′), using the TD of the new data, Δx(t′) ∈ R1 ×n, calculated as follows: (4)
(8)
where p and q are constants and t(t) ∈ Rm×1 is a time vector. From eqs 1, 7, and 8, TD of x(t), Δx(t), can be transformed as follows:
(3)
Δx(t ′) = x(t ′) − x(t ′ − i)
x(t − i) = xL(t − i) + xNL(t − i)
xL(t ) = p t(t ) + q
Then, the relationship between ΔX(t) and Δy(t) is modeled by a regression method. Δy(t ) = f (ΔX(t )) + e
(7)
xL is a linear function of t given by
2. METHODS After a brief introduction of the TD modeling method, the approach is used to eliminate a term that is a linear function of time. The effect of the signal-to-noise (SN) ratios of a process variable and the TD of the process variable are then explained. Finally, we formulate the predictive accuracy of TD models, taking into consideration the determination coefficient (r2), autocorrelation, and the SN ratios. 2.1. Time Difference Modeling Method. In a traditional procedure, the modeling relationship between explanatory variables or independent variables, X(t) ∈ Rm×n, with m rows of data samples and n columns of variables, and an objective variable or a dependent variable, y(t) ∈ Rm×1, is done by regression methods after preparing data, X(t) and y(t), related to time t. In terms of prediction, the constructed model predicts the value of y(t′) with the new data x(t′) ∈ R1×n. In TD modeling, the TD of X and that of y, ΔX(t) and Δy(t), are first calculated between the present values, X(t) and y(t), and those at some time i before the target time, X(t−i) ∈ Rm×n and y(t−i) ∈ Rm×1. ΔX(t ) = X(t ) − X(t − i)
x(t ) = xL(t ) + xNL(t )
ηx =
(6)
because y(t′−i) is given previously. TD models can be constructed and values of TD can be predicted even when interval i is not constant.
var(xV (t )) var(xE(t ))
(13)
where var(xV(t)) and var(xE(t)) are the variances of xV(t) and xE(t), respectively. Then, we consider TD of x(t), Δx(t) ∈ Rm×1, as 1323
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
Δx(t ) =x(t ) − x(t − i)
RMSEdiff =
= {xV (t ) + xE(t )} − {xV (t − i) + xE(t − i)} = ΔxV (t ) + ΔxE(t )
where (15)
ΔxE(t ) = xE(t ) − xE(t − i)
(16)
= =
+ y(t − i)} = Δy(t ) − Δycalc (t )
var(ΔxV (t )) var(ΔxE(t )) var(xV (t ) − xV (t − i)) var(xE(t ) − xE(t − i))
r =1−
∑ {y(t ) − y ̅ (t )}
RMSE =
(18)
(26)
rdiff2
∑ {y(t − i) − y ̅ (t − i)}2 ≈ ∑ {y(t ) − y ̅ (t )}2
(20)
(27)
(28)
applies if the TD interval is sufficiently small and there is little difference between y(t) and y(t−i). By using eq 28, eq 27 can be expressed as follows:
∑ {Δy(t ) − Δy(t )}2 = 2 ∑ {y(t ) − y ̅ (t )}2 − 2 ∑ {(y(t ) − y ̅ (t ))(y(t − i) − y ̅ (t − i))}2
∑ {Δy(t ) − Δycalc (t )}2 ∑ {Δy(t ) − Δy(t )}
∑ {y(t ) − y ̅ (t )}2
The following relationship between the first and third terms in eq 27
∑ {y(t ) − ycalc (t )}2
2
∑ {Δy(t ) − Δy(t )}2
∑ {Δy(t ) − Δy(t )}2 = ∑ {(y(t ) − y ̅ (t )) − (y(t − i) − y ̅ (t − i))}2 = ∑ {y(t ) − y ̅ (t )}2 − 2 ∑ {(y(t ) − y ̅ (t ))(y(t − i) − y ̅ (t − i))}2 + ∑ {y(t − i) − y ̅ (t − i)}2
where ycalc is the calculated y value. r2 represents the fitting accuracy of constructed models, and values close to unity for r2 are favorable. The lower the RMSE value, the more accurate the prediction obtained with the constructed model. The r2 values calculated between Δy and Δycalc (rdiff2) and the RMSE values calculated between Δy and Δycalc (RMSEdiff) are expressed as follows: rdiff 2 = 1 −
(25)
The relationships between r and and between RMSE and RMSEdiff are thus obtained. When the numerator and the denominator in eq 26 are divided by m − 1, it is evident that A is the ratio of the variances of Δy(t) and y(t). If the values of rdiff2 are the same, small values of A, i.e., a large variance of y(t) and a small variance of Δy(t), reflect large values of r2. In other words, even if rdiff2 values are small, r2 values become large when the variance of Δy(t) is small (in the extreme case, r2 = 1 if A = 0). In this case, RMSEdiff values are small since the variance of Δy(t) is small, but there is a possibility that values of the TD of y might not be calculated accurately. Not only r2 but also rdiff2 should be checked when constructing TD models. The approach used in this section is applied to validation data and test data of y. 2.5. Autocorrelation Coefficient for y and the Performance of TD Models. The numerator of A in eq 26 can be transformed as follows:
(19)
m
RMSE = RMSEdiff
2
∑ {y(t ) − ycalc (t )}2 2
(24)
A=
From eqs 13 and 18, although the large SN ratios are favorable, ηΔx is smaller than ηx if xV(t) is autocorrelated. This means that the SN ratio decreases upon using TD because a process variable and the i delayed variable do not usually have negative correlation in process data if i is relatively small. On the one hand, a term that is a linear function of time can be eliminated by using TD, as confirmed in section 2.2; on the other hand, TD decreases the SN ratio. Therefore, when TD models are used, we should pay attention to the effect of the autocorrelation on the SN ratio, although linear trends with time are removed by TD. Of course, the above formulas are applicable not only to x but also to y. 2.4. Performance of TD Models Relevant to Estimated y and Δy. The r2 and root mean squared error (RMSE) values are used as measures of accuracy of models and are defined as 2
1 − r 2 = A(1 − rdiff 2)
where A is given as follows:
var(xV (t )) + var(xV (t − i)) − 2cov(xV (t ), xV (t − i)) var(xE(t )) + var(xE(t − i)) (17)
var(xV (t )) − cov(xV (t ), xV (t − i)) var(xE(t ))
(23)
Therefore, eqs 19 and 21 and eqs 20 and 22 can be combined as
where cov(xV(t),xV(t−i)) is the covariance between xV(t) and xV(t−i), which means the autocorrelation16 of xV. Here, var(xV(t−i)) = var(xV(t)) and var(xE(t−i)) = var(xE(t)) are assumed, which is reasonable when i is relatively small. Thus, ηΔx is rewritten as follows: ηΔx =
(22)
y(t ) − ycalc (t ) ={Δy(t ) + y(t − i)} − {Δycalc (t )
Because xE(t) is not autocorrelated, the SN ratio of Δx(t), ηΔx, is calculated as ηΔx =
m
From eqs 2 and 6, in which ypred becomes ycalc, a formula of the difference between y(t) and ycalc(t) is derived as follows:
(14)
ΔxV (t ) = xV (t ) − xV (t − i)
∑ {Δy(t ) − Δycalc (t )}2
(21)
(29)
Therefore, A becomes 1324
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
⎛ ∑ {(y(t ) − y ̅ (t ))(y(t − i) − y ̅ (t − i))}2 ⎞ ⎟ A =2⎜1 − ∑ {y(t ) − y ̅ (t )}2 ⎠ ⎝
plant) we examine TD intervals and the predictive ability of TD models. Furthermore, through the analysis of two sets of real industrial data, TD intervals are adjusted to increase the predictive accuracy of TD models. 3.1. Analysis of numerical simulation data. In this study, the number of X variables was set as 1 and x(t) was prepared as follows:
= 2[1 − corrcoef(y(t ), y(t − i))] (30)
where corrcoef(y(t),y(t−i)) is the correlation coefficient between y(t) and y(t−i), i.e., the autocorrelation coefficient for y. From eqs 24 and 30, we confirmed that the larger the autocorrelation coefficient for y is, the smaller the values of A are and the larger the r2 values are. 2.6. SN Ratios and the Performance of TD Models. Equation 19 can be transformed as
x(t ) = sin(πt /250)
y(t ) = x(t ) + ns × N (0, 1)
∑ {y(t ) − y ̅ (t )}2
=1−
∑ {y(t ) − ycalc (t )}2 /(m − 1) ∑ {y(t ) − y ̅ (t )}2 /(m − 1) var(yF(t )) var(y(t ))
(31)
where yF(t) ∈ R is the model error of y(t). In addition, ηy can be transformed as m×1
ηy = =
var(y(t ) − yE(t )) var(yE(t )) var(y(t )) − 2 cov(y(t ), yE(t )) + var(yE(t ))
=1+
(36)
where N(0,1) is a random number from the normal distribution with a standard deviation of 1 and a mean of 0 and ns means the amount of noise. The number of data was set as 2000. The first 1000 data were used for training and the next 1000 data were the test data. We used the ordinary linear least-squares regression method as the regression approach. First, ns was changed from 0 to 1 in steps of 0.05, and regression modeling and prediction were carried out by using the data obtained. r2 for the test data in this case was represented by rpred,norm2. Then, the TD interval was changed from 1 to 120 in steps of 1, ns was varied as above, and regression modeling and prediction were performed. The r2 and RMSE values for the test data of Δy are represented by rpred,diff2 and RMSEpred,diff, respectively. The r2 and RMSE values for the test data after conversion back from Δy to y by using eq 6 are represented by rpred2 and RMSEpred, respectively. Figure 1 shows the relationship between rpred,diff2 and rpred2 and that between RMSEpred,diff and RMSEpred. As we confirmed
∑ {y(t ) − ycalc (t )}
=1−
(35)
After x(t) was autoscaled, y(t) was set as
2
r 2 =1 −
(t = 1, 2, ..., 2000)
var(yE(t )) var(y(t )) var(yE(t ))
(32)
because there is no correlation between y(t) and yE(t) ∈ Rm×1, which is the white noise of y(t). When the relationship between X and y or ΔX(t) and Δy(t) is modeled appropriately, yF(t) equals yE(t). In this case, from eqs 31 and 32, the relationship between r2 and ηy is given as follows: 1 1 =1+ 2 η r y (33) This can also be applied to rdiff2 and ηΔy as follows: 1 1 =1+ 2 ηΔy rdiff
Figure 1. Relationship between the prediction ability of Δy and that of y using numerical simulation data: (a) the plot of rpred,diff2 and rpred2 and (b) the plot of RMSEpred,diff and RMSEpred. (34)
These relationships between the SN ratios and the performance of TD models can be obtained only if yF(t) = yE(t), i.e., the relationships between X and y or ΔX(t) and Δy(t) are modeled appropriately. The formulas developed above are also applicable to validation data and test data of y. yF(t) = yE(t) is essential for the equality of eqs 33 and 34 to be true; however, the data for X include noise, and it will be difficult for yF(t) = yE(t) to be realized in practice. In addition, yE(t) will be unknown in real cases.
in eqs 24 and 25, RMSEpred,diff and RMSEpred were equal (Figure 1b), although rpred,diff2 did not relate to rpred2 (Figure 1a). Then A in eq 26 was calculated and the relationship between A(1− rpred,diff2) and 1 − rpred2 was plotted in Figure 2. As shown in eq 24, A(1−rpred,diff2) and 1 − rpred2 were equal. Figure 3 shows the relationship between the autocorrelation coefficient for y and A. The data were directly aligned in the straight line of eq 30. We could therefore confirm the accuracy of eqs 24, 25, and 30 with respect to the performance of TD models and the autocorrelation coefficient for y. The relationships between 1/rpred,diff2 and 1 + 1/ηΔy (blue dots) and between 1/rpred,norm2 and 1 + 1/ηy (pink asterisks) are shown in Figure 4. Figure 4a is the whole plot and Figure 4b is an expanded plot for 1 + 1/η in the range 0.5−5. Equations 33 and 34 were therefore verified as valid relationships between
3. RESULTS AND DISCUSSION First, we analyze numerical simulation data to verify the relationships and formulas of the TD model, and then, using data obtained by dynamic simulation of a depropanizer distillation column modeling (based on a real industrial 1325
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
Figure 2. Relationship between A(1−rpred,diff2) and 1 − rpred2 using numerical simulation data. For A see eq 26.
Figure 5. Relationship between rpred2 and var(ΔyE)/var(Δy) (blue dots) and that between rpred,norm2 and var(yE)/var(y) (pink asterisks) using numerical simulation data.
right of the line have A values that are lower than 1, and the points for rpred2 and var(ΔyE)/var(Δy) to the left of the line have A values that are higher than 1. When A is less than 1, i.e., the variation of Δy is small, the previous value of y has a significant impact on the value of y in eq 6, and thus ypred matches y after adding the previous y to Δypred. In contrast, when A is greater 1, i.e., the variation of y is small, the previous value of y has a small impact on the value of y in eq 6, and thus ypred does not match y even after adding the previous y to Δypred. In this numerical simulation, noise was added to y only. Almost the same results were obtained when noise was added to x only. However, when noise was added to both x and y, the relationships between 1/rpred,diff2 and 1 + 1/ηΔy, and rpred2 and var(ΔyE)/var(Δy) shown in Figures 4 and 5 could not be confirmed and the results were complex. This is because appropriate models could not be constructed with data where both x and y have noise and yF is different from yE. In addition, we do not always have quantitative information about noise, and therefore, it is quite difficult to calculate rpred2, rpred,diff2, and so on from eqs 33, 34, and others, although rpred2 and rpred,diff2 can be theoretically estimated from ηy and ηΔy with eqs 33 and 34. Thus, constructed models must be carefully validated in terms of their predictive accuracy in real cases. As shown in sections 2.2 and 2.3, noise in process variables did not affect the SN ratios, and linear trends with time are removed by TD, although the SN ratios increase by TD
Figure 3. Relationship between autocorrelation coefficient for y and A using numerical simulation data. For A see eq 26.
the performance of TD models and the SN ratio. Normal models and TD models can be integrated in terms of those relationships. When ns is large, the variance from the diagonal is also large in Figure 4. This is because the assumption that yF and yE are equal is not realized with large amounts of noise. The relationships between rpred2 and var(ΔyE)/var(Δy) (blue dots) and between rpred,norm2 and var(yE)/var(y) (pink asterisks) are shown in Figure 5. The straight line relationship of eq 31 between rpred,norm2 and var(yE)/var(y), where yF equals yE, was confirmed, indicating that appropriate models could be constructed. The data for rpred2 and var(ΔyE)/var(Δy) to the
Figure 4. Relationships between 1/rpred,diff2 and 1 + 1/ηΔy (blue dots), and 1/rpred,norm2 and 1 + 1/ηy (pink asterisks) using numerical simulation data. (a) the whole plot. (a) the plot whose ranges are both from 0.5 to 5. 1326
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
because of the autocorrelation in process variables. However, noise in process variables makes it difficult to calculate rpred2, rpred,diff2 and so on from eqs 33, 34, and others. This is a limitation not only for TD models but also for normal models. To reduce noise in process variables, dimensionality reduction methods such as principal component analysis 17 and independent component analysis method4,18 and a smoothing method such as the Savitzly−Goley19 method are useful for raw measured data. 3.2. Analysis of Dynamic Simulation Data. In this section, we discuss TD intervals and the predictive ability of TD models using dynamic simulation data. Visual Modeler (VM) developed by Omega Simulation Co., Ltd.20 was used as the dynamic simulator. Figure 6 shows the target plant (a
Table 2. Process Variables Measured in the Depropanizer Distillation Column Selected for Dynamic Simulation symbol A no. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
F1 F2 F3 F4 F5 F6 P1 P2 T1 T2 T3 T4 T5 T6 T7 L1 L2
feed flow bottom flow top flow steam flow flow to reboiler reflux flow pressure 1 pressure 2 feed temperature top temperature bottom temperature steam temperature temperature 1 temperature 2 temperature 3 liquid level 1 liquid level 2
molar concentrations of propane and n-butane and the steam temperature randomly walked within ±1%, ± 1%, and ±3 °C as noise. Process variables other than those in Table 2 were removed because the correlation coefficients of pairs of the variables were greater than 0.99. The setting values of each instrument and control system were given as default values. Two types of disturbances were generated in this study. One was changes in the steam temperature, i.e., 200 °C from 1 to 2880 min (48 h), 185 °C from 2881 to 5760 min (96 h), 200 °C from 5761 to 8640 min (144 h), and 185 °C from 8641 to 11 520 min (192 h). The other was changes in the feed concentrations of propane and n-butane; the concentration of n-butane was 35% from 1 to 2880 min (48 h), 30% from 2881 to 5760 min (96 h), 35% from 5761 to 8640 min (144 h), and 30% from 8641 to 11 520 min (192 h). The feed concentration of propane changed accordingly. In reality, concentrations are measured at intervals from minutes to hours with a time delay, but in this case study, all process variables were simulated at 1min intervals to allow the detailed investigation of TD intervals. Random numbers from the normal distribution with a mean of 0 and a standard deviation 0.03 times that of the standard deviation of each process variable were added to the process variable. This added noise may be relatively small compared with that of real industrial data, but the amount of the noise does not matter because relationships between TD intervals and the predictive ability of TD models were discussed in this case study. The first 5760 min (96 h) of data were used for the training and the next 5760 min (96 h) of data were the test data. To incorporate the dynamics of process variables into soft sensor models,21 X included each process variable delayed for durations ranging from 0 to 60 min in steps of 10 min, i.e., X consists of 7 time points (0, 10, 20, 30, 40, 50, 60 min) for each of 17 variables (7 × 17 = 119). Partial least-squares (PLS) regression22 was used for statistical analysis. When the steam temperature changed, the TD intervals of the training data were set as 1 and 30 min, the TD intervals of the test data were set as 1 and 30 min, and PLS models were constructed with training data and used for the prediction of
Figure 6. Schematic representation of the depropanizer distillation column selected for dynamic simulation.
depropanizer distillation column) that was used in the VM training program. In VM, its own calculation routine is used as solver, and we selected the Soave−Redlich−Kwong equation as an equation of state. Table 1 shows the process conditions and Table 1. Process Conditions and Parameters of the Column in Dynamic Simulation Process Conditions feed rate 7.5 × 103 kg/h feed temperature 87.8 °C external temperature 25.0 °C atmospheric pressure 101.325 kPa Parameters of the Column number of stages feed stage cross-sectional area of tray effective weight
objective variable mol concentration of bottom n-butane symbol explanatory variables
28 16 1.54 m2 1.00 × 104 kg
parameters of the column. The y variable is the molar concentration of the bottom n-butane and the X variables are the 17 variables shown in Table 2. The basic feed molar concentrations of ethane, propane, isobutene, n-butane, and isopentane were set as 0.2%, 65%, 4.7%, 30%, and 0.1%, respectively. The steam temperature was set as 185 °C. The 1327
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
Figure 7. Relationships between observed and predicted Δy of test data when steam temperature changed in the depropanizer distillation column.
Figure 8. Relationships between observed and predicted y of test data when the steam temperature changed in the depropanizer distillation column. The TD interval of test data is 30 min.
test data. Figure 7 shows the relationships between the observed and predicted Δy when the TD intervals of the training data and test data were 1 and 30 min. When the TD interval of the test data was 1 min, both rpred,diff2 values were low and the models were not effective in predicting the variations of Δy in test data. TD intervals of 1 and 30 min for the training data did not affect the predictive performance of the TD models; Δy with a TD interval of 1 min was more affected by amplified noise than by significant variation of y, and the SN ratio was quite small. Parts c and d of Figure 7 show the results when the TD interval of the training data was 1 and 30 min, respectively. In Figure 7c, there were outlying data away from
the diagonal, and so with a training data TD interval of 1 min, an appropriate model could not be constructed because the effect of noise was stronger than that of variations in y, i.e., the SN ratio was small. In Figure 7d, conversely, the plot showed a relatively tight cluster of predicted values along the diagonal and the rpred,diff2 value was higher than that of Figure 7c; thus, the TD model with a large TD interval of 30 min resulted in high predictive accuracy. It was thus confirmed that a small TD interval does not always give good results, and appropriate TD intervals need to be set for the accurate prediction of Δy. Figure 8 shows the relationships between the observed and predicted y values when the TD interval of the test data was 30 1328
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
Figure 9. Relationship between observed and predicted Δy of test data when the feed concentrations changed in the depropanizer distillation column.
Figure 10. Relationship between observed and predicted y of test data when the inlet concentrations changed in the depropanizer distillation column. The TD interval of test data is 30 min.
min. The rpred2 values were almost same, but when the TD interval of the training data was 1 min, the data for y values in the range from approximately 86 to 88 had large errors, as shown in Figure 8a. This was not unexpected, because the data in Figure 8a came from the results shown in Figure 7c, which had poor prediction accuracy of Δy. In contrast, the model trained with data whose TD interval was 30 min could predict values of y with high accuracy. For changes in the feed concentrations, TD intervals of training data and test data were set at 1 and 30 min, as was the
case for changes in steam temperature. PLS models were constructed with training data and were used to predict the test data. The details of PLS are shown in the Appendix. The relationships between the observed and predicted Δy values are shown in Figure 9. The rpred,diff2 values were low, data were outlying far from the diagonals, and the predictive accuracy was low when the test data TD interval was 1 min, just as happened when the steam temperature changed. Even when the TD interval of the test data was 30 min, the model constructed with a training data TD interval of 1 min could not predict Δy 1329
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
Figure 11. The rpred2 values for different training data TD intervals when test data TD intervals were 10, 30, and 120 min.
Chemical Corp. Figure 12 shows a schematic of the first distillation column, and Table 3 shows the process variables.
accurately. An appropriate model could not be obtained with such a small training data TD interval when the ratio of noise to the significant variation of Δy was large. In contrast, when both the training data and test data TD intervals were 30 min, the plot in Figure 9d shows a tight cluster of predicted values along the diagonal, meaning that the predictive accuracy was high. The rpred,diff2 value was higher than that of Figure 9c. By setting relatively large TD intervals, the ratio of variance of y and Δy to noise increased and the SN ratio became large. As a result, the predictive accuracy of the TD model increased on the basis of eqs 33 and 34. Figure 10 shows the relationships between observed and predicted y values when the TD interval of test data was 30 min. As was the case for changes in the steam temperature, the rpred2 values were almost same, but when the training data TD interval was 1 min, some predicted data were outlying from the diagonal in Figure 10a. However, the model with a training data TD interval of 30 min could predict values of y with high accuracy. From Figures 7−10, we can say that TD models with small TD intervals have poor prediction accuracy of Δy, and those with suitably large TD intervals have good prediction accuracy of Δy and can predict values of y with high reliability. By analyzing two types of disturbances in a depropanizer distillation column, it was confirmed that it is important to determine appropriate TD intervals for model construction and for predictions. In fact, if the measurement interval of y is constant and each measured value of y can be obtained, this interval can be used as the TD interval and a TD model predicts values of y. However, if the measurement interval of y is irregular or a y-analyzer is broken, a TD interval has to change for each prediction of Δy and y. Therefore, in the case of steam temperature changes and feed concentration changes, the training data TD interval changed continuously for each interval of test data and we investigated the predictive ability of the TD models. The TD interval of the test data was changed to 10, 30, and 120 min, and the rpred2 values were calculated for training data TD intervals of every minute from 1 to 60 min (see Figure 11). As the test data TD interval increased, larger training data TD intervals were needed to increase the rpred2 values. The SN ratio increases as a result of the increase of the test data TD interval, but long-term variation of Δy must be modeled, and therefore, also the training data TD intervals have to be set large. Appropriate training data TD intervals must be determined to accurately predict values of y with each TD interval of new data. 3.3. Analysis of Real Industrial Data. We analyzed data obtained from the operation of distillation columns at the Mizushima Works and the Kashima Works, Mitsubishi
Figure 12. Schematic representation of the first distillation column.
The y variable is the concentration of the bottom product with the lowest boiling point, and the X variables are the 19 variables given in Table 3. The measurement interval of y was 30 min. To incorporate the dynamics of process variables into soft sensor models, X was composed of each explanatory variable delayed for durations from 0 to 60 min in steps of 10 min, i.e., X consists of 7 time points (0, 10, 20, 30, 40, 50, 60 min) for each of 19 variables (7 × 19 = 133). The PLS method was used as the regression method. The details of PLS are shown in the Appendix. We collected data from January 2003 to December 2006 and divided them into periods of 2 days. Subsequently, repeated training and testing were carried out in series for periods of 2 days each. The ith period is training data and the (i + 1)th period is test data. Modeling and prediction were conducted for various TD intervals for training data and test data. Then, on incrementing i, the procedure was repeated. Data reflecting variations that could have been caused by troubles with the concentration analyzer and plant inspections were eliminated in advance.10 The TD intervals for training data and test data were independently changed from 30 to 300 min in steps of 30 min. Figure 13 shows the relationships between the rpred2 values where the TD interval for each training data was fixed at 30 min 1330
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
having lower variance (as a result of a shorter TD) tend to have larger rpred2 values. Table 4 shows the number of data for which rpred2 is the largest value for a range of training data TD intervals for test data TD intervals of 30 and 90 min. With a test data TD interval of 30 min, there were large numbers of optimum data points for training data TD intervals of 30 and 60 min. With some large TD intervals, the models had high predictive accuracy due to the noise in the data. With a test data TD interval of 90 min, the largest number of optimum data points was when the training data TD interval was 90 min; however, the numbers for training data TD intervals of 30 and 60 min were also large. Thus, an appropriate training data TD interval exists for each 2-day test data and each test data TD interval. Next, the predictive accuracy of the TD models and the TD intervals are discussed with respect to the data measured in the second distillation column; a schematic representation of this column is given in Figure 14, and Table 5 lists 23 process variables. A y-variable is the concentration of the bottom product with the lower boiling point, and X-variables are the 23 process variables. The measurement interval of y is 8 min. To incorporate the dynamics of process variables into soft sensor models, X included each explanatory variable that was delayed for durations ranging from 0 to 60 min in steps of 10 min. The PLS method was used as a regression method. The details of PLS are shown in the Appendix. The data were collected from July 2008 to March 2011 and divided into sets over 2-day periods. Then, sequential data sets, spanning two 2-day periods, were used in training and testing: the ith span being training data and the (i + 1)th span being test data, with modeling and prediction conducted. Next, i is incremented and the procedure repeated. Data reflecting variations that could have been caused by troubles with the concentration analyzer and plant inspections were eliminated in advance.10 The TD intervals of the training and test data were independently changed from 8 to 80 min in steps of 8 min. Figure 15 shows the relationship between the rpred2 values for an 8 min TD interval and those where the training data TD interval was optimized when the test data TD interval was 8 and 24 min. The ranges of each plot are 0 < rpred2 < 1. With data located above the line, the predictive accuracy of the TD models is seen to have increased by optimizing the TD
Table 3. Process Variables Measured in the First Distillation Column symbol
objective variable
A no.
symbol
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
F1 F2 F3 F4 F5 F6 L1 P1 P2 T1 T2 T3 T4 T5 T6 T7 T8 F4/F3 = R F1/F6 = F
bottom product concentration explanatory variables reflux flow reboiler flow feed 1 flow feed 2 flow bottom flow top flow liquid level pressure 1 pressure 2 temperature 1 temperature 2 temperature 3 temperature 4 bottom temperature feed 1 temperature feed 2 temperature top temperature reflux ratio feed flow ratio
and those where the TD interval for each training data was maximized for each rpred2 value for test data TD intervals of 30 and 90 min. The ranges of each plot are constrained between 0 and 1 inclusive. Many data are located in the upper left part of the graph, meaning that the predictive accuracy of the TD models was increased by optimizing the training data TD intervals. Robust TD models less affected by noise can be constructed with relatively large TD intervals. In this case, 77.7% and 83.6% of rpred2 values increased after training data TD interval optimization when test data TD intervals were 30 and 90 min, respectively. The minimum TD interval was relatively large at 30 min, and therefore, the difference between the two percentages is small. The values of rpred2 and the number of data points in the plot seemed to be larger in Figure 13a than in Figure 13b. This can be understood from eqs 24 and 26, which indicate that test data
Figure 13. Relationships between the rpred2 values when the TD interval was 30 min and the rpred2 values when the TD interval was optimized in the first distillation column. The TD intervals of test data were 30 and 90 min. The ranges were constrained to the interval 0 through 1. 1331
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
Table 4. Number of Data Intervals Where the rpred2 is the Largest Value with a Target TD Interval of Training Data for Each TD Interval of Test Data number of data intervals where the rpred2 is the largest value with a target TD interval of training data TD interval of test data
30 min 90 min
30 min
60 min
90 min
120 min
150 min
180 min
210 min
240 min
270 min
300 min
169 124
163 128
90 157
77 89
56 56
49 49
28 43
44 33
40 37
52 42
intervals. Robust TD models less affected by noise will be constructed with some large TD intervals. Comparing parts a and b of Figure 15, the number of data where rpred2 values were increased by optimizing the training data TD intervals was more than that for which the training data TD interval is fixed at 8 min for both test data TD intervals, namely, 8 and 24 min. In fact, only 65.0% compared to 83.6% of rpred2 had increased when test data TD intervals were 8 and 24 min, respectively. This is because, as verified with dynamic simulation data, the predictive accuracy of the TD models (rpred,diff2) is low and is not affected by training data TD intervals when test data TD intervals are small. Meanwhile, a small test data TD interval implies that rpred2 values tend to be large because of the small variance of Δy and eqs 24 and 26. The numbers of data where rpred2 has the largest value for a target training data TD interval for test data TD intervals of 8 and 24 min are shown in Table 6. With a test data TD interval of 8 min, the numbers for training data TD intervals of 8 and 16 min were large. With some large TD intervals, the models had high predictive accuracy due to the noise in the data. The numbers for the training data TD interval of 24 min was largest for a test data TD interval of 24 min, though those of 16 and 32 min were also large. The significance is that large training data TD intervals were needed for certain TD models because of the choice of some large TD intervals. We can say that appropriate training data TD intervals exist for each test data and corresponding TD interval. In fact, measured values of y sometimes do not exist because of, for example, a y-analyzer fault. The TD intervals of new data must be changed and large under such situations because measured values cannot be obtained. Through the analyses of dynamic simulation data and two real industrial data sets, we confirmed that appropriate training data TD intervals need to be chosen for performing high predictive accuracy with large TD intervals of new data. In this paper, the training data TD intervals could be optimized by checking the predictive accuracy of the test data. However, optimization of training data TD intervals for a test data TD interval should be performed using only training data in practice.
Figure 14. Schematic representation of the second distillation column.
Table 5. Process Variables Measured in the Second Distillation Column symbol
objective variable
A no.
symbol
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
F1 F2 F3 F4 F5 F6 F7 L1 P1 P2 P3 P4 P5 T1 T2 T3 T4 T5 T6 T7 T8 T9 T10
bottom product concentration explanatory variables feed flow reboiler 1 flow top flow 1 flow top flow 2 flow reflux flow reboiler 2 flow bottom flow liquid level pressure 1 pressure 2 pressure 3 pressure 4 pressure 5 feed temperature top temperature temperature 1 temperature 2 temperature 3 temperature 4 temperature 5 temperature 6 reboiler temperature bottom temperature
4. CONCLUSION In this study, we discussed TD models and optimal TD intervals to improve the predictive accuracy of soft sensor models. Various kinds of information on TD models and TD intervals could be obtained in terms of noise and variance of data, autocorrelation of process variables, accuracy of TD models, etc. For example, the performance in predicting Δy and y can be analyzed from eqs 24−26, whereas the SN ratios and the performance of models can be obtained from eqs 33 and 34. The equations derived in this study were confirmed by using numerical simulation data. In addition, we could integrate the performance of a normal model and that of a TD model. In the end, the autocorrelation coefficient for y, A, ηy, and ηΔy are indexes affecting r2, which represents the performance of TD models. 1332
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research
Article
Figure 15. Relationship between the rpred2 values when the TD interval is 8 min and the rpred2 when the TD interval is optimized in the second distillation column. The TD intervals of test data were 8 and 24 min.
Table 6. Number of Data Intervals Where the rpred2 is the Largest Value with a Target Training Data TD Interval for Each Test Data TD Interval number of data intervals where the rpred2 is the largest value with the TD interval of training data TD interval of test data
8 min 24 min
8 min
16 min
24 min
32 min
40 min
48 min
56 min
64 min
72 min
80 min
99 33
91 68
51 69
16 40
10 25
7 13
4 9
3 10
2 6
0 10
b = W(P′W)−1q
TD intervals and the performance of TD models were analyzed through dynamic simulation data based on a depropanizer distillation column and two sets of real industrial data measured in these columns. The performance of the TD models improved by setting appropriate TD intervals. In fact, when values of y cannot be measured for a period of time, because of y-analyzer faults, for example, a TD interval can be set large and a predictive TD model can be constructed even under such situations by selecting an appropriate TD interval. In this study, we did not consider unobserved disturbances. The comparison of TD models and normal models should be done with the unobserved disturbances, which are represented as nonlinear functions of time. Additionally, a task for the future is the determination of the optimal TD intervals with only training data. By using the information we obtained in this study, the predictive ability of soft sensors will improve with the goal of maintaining chemical plant stablity.
where W is an X-weight matrix and b is a vector of regression coefficients.
■
* Tel: +81-3-5841-7751. Fax: +813-5841-7771. E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
y = Tq + f
(A.2)
■
NOTATION A ratio of the variances of Δy and y corrcoef correlation coefficient cov covariance e error f regression model i interval of time difference m number of data n number of variables ns amount of noise p, q constants r2 determination coefficient rdiff2 determination coefficient of time difference rpred2 determination coefficient for test data
where T is a score matrix, P is an X-loading matrix, q is a yloading vector, E is a matrix of X residuals, and f is the vector of y residuals. The PLS regression model is y = Xb + const
ACKNOWLEDGMENTS
H.K. is grateful for financial support of the Japan Society for the Promotion of Science (JSPS) through a Grant-in-Aid for Young Scientists (B) (No. 24760629). The authors acknowledge the support of Mizushima Works and Kashima Works, Mitsubishi Chemical Co., and the financial support of Mizuho Foundation for the Promotion of Sciences and the Toyota Physical and Chemical Research Institute.
APPENDIX A. PLS PLS is a method for relating explanatory variables or independent variables, X, and an objective variable or an objective variable, y. Using a linear multivariate model, it goes beyond traditional regression methods in that it also models the structures of X and y. In PLS modeling, the covariance between y and the score vector ti is maximized. A PLS model consists of the following two equations (A.1)
AUTHOR INFORMATION
Corresponding Author
■
X = TP′ + E
(A.4)
(A.3) 1333
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334
Industrial & Engineering Chemistry Research rpred,diff2 rpred,norm2 RMSE RMSEdiff RMSEpred RMSEpred,diff SN t t TD x X xE xNL xL xV y yE yF(t) ycalc ypred var ΔX Δy η
■
Article
(12) Kaneko, H.; Funatsu, K. A Soft Sensor Method Based on Values Predicted from Multiple Intervals of Time Difference for Improvement and Estimation of Prediction Accuracy. Chemom. Intell. Lab. Syst. 2011, 109, 197. (13) Kaneko, H.; Arakawa, M.; Funatsu, K. Applicability Domains and Accuracy of Prediction of Soft Sensor Models. AIChE J. 2011, 57, 1506. (14) Ookita, K. Operation and Quality Control for Chemical Plants by Soft Sensors. CICSJ Bull. 2006, 24, 31 (in Japanese).. (15) Okada, T.; Kaneko, H.; Funatsu, K. Development of a Model Selection Method Based on Reliability of a Soft Sensor Model. Songklanakarin J. Sci. Technol. 2012, 34, 217. (16) Alkemade, C. T. J.; Snelleman, W.; Boutilier, G. D.; Pollard, B. D.; Winefordner, J. D.; Chester, T. L.; Omenetto, N. A Review and Tutorial Discussion of Noise and Signal-to-Noise Ratios in Analytical SpectrometryI. Fundamental Principles of Signal-to-Noise Ratios. Spectrochim. Acta, Part B. 1978, 33, 383. (17) Wold, S.; Esbensen, K.; Geladi, P. Principal Component Analysis. Chemom. Intell. Lab. Syst. 1987, 2, 37. (18) Comon, P. Independent Component Analysis, a New Concept? Signal Process. 1994, 36, 287. (19) Savitzky, A.; Golay, M. J. E. Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Anal. Chem. 1964, 36, 1627. (20) http://www.omegasim.co.jp/product/vm/index.htm. (21) Kaneko, H.; Funatsu, K. A New Process Variable and Dynamics Selection Method Based on Genetic Algorithm-Based Wavelength Selection Method. AIChE J. 2012, 58, 1829−1840. (22) Wold, S.; Sjöström, M.; Eriksson, L. PLS-Regression: A Basic Tool of Chemometrics. Chemom. Intell. Lab. Syst. 2001, 58, 109.
determination coefficient of time difference for test data determination coefficient for test data when a normal regression model is used root mean squared error root mean squared error of time difference root mean squared error for test data root mean squared error of time difference for test data signal-to-noise time time vector time difference process variable explanatory variable or independent variable white noise term that is a nonlinear function of time term that is a linear function of time significant variation objective variable or dependent variable white noise of an objective variable or a dependent variable model error of an objective variable or a dependent variable calculated value of an objective variable or a dependent variable predicted value of an objective variable or a dependent variable variance time difference of explanatory variables or independent variables time difference of an objective variable or a dependent variable signal-to-noise ratio
REFERENCES
(1) Kano, M.; Nakagawa, Y. Data-Based Process Monitoring, Process Control, and Quality Improvement: Recent Developments and Applications in Steel Industry. Comput. Chem. Eng. 2008, 32, 12. (2) Kadlec, P.; Gabrys, B.; Strandt, S. Data-Driven Soft Sensors in the Process Industry. Comput. Chem. Eng. 2009, 33, 795. (3) Qin, S. J. Recursive PLS Algorithms for Adaptive Data Modelling. Comput. Chem. Eng. 1998, 22, 503. (4) Kaneko, H.; Arakawa, M.; Funatsu, K. Development of a New Soft Sensor Method Using Independent Component Analysis and Partial Least Squares. AIChE J. 2009, 55, 87. (5) Cheng, C.; Chiu, M. S. A New Data-Based Methodology for Nonlinear Process Modeling. Chem. Eng. Sci. 2004, 59, 2801. (6) Fujiwara, K.; Kano, M.; Hasebe, S.; Takinami, A. Soft-Sensor Development Using Correlation-Based Just-in-Time Modeling. AIChE J. 2009, 55, 1754. (7) Kadlec, P.; Gabrys, B. Local Learning-Based Adaptive Soft Sensor for Catalyst Activation Prediction. AIChE J. 2010, 57, 1288. (8) Liu, Y.; Huang, D.; Li, Y. Development of Interval Soft Sensors Using Enhanced Just-in-Time Learning and Inductive Confidence Predictor. Ind. Eng. Chem. Res. 2012, 51, 3356. (9) Kadlec, P.; Grbic, R.; Gabrys, B. Review of Adaptation Mechanisms for Data-Driven Soft Sensors. Comput. Chem. Eng. 2011, 35, 1. (10) Kaneko, H.; Funatsu, K. Maintenance-Free Soft Sensor Models with Time Difference of Process Variables. Chemom. Intell. Lab. Syst 2011, 107, 312. (11) Kaneko, H.; Funatsu, K. Development of Soft Sensor Models Based on Time Difference of Process Variables with Accounting for Nonlinear Relationship. Ind. Eng. Chem. Res. 2011, 50, 10643. 1334
dx.doi.org/10.1021/ie302582v | Ind. Eng. Chem. Res. 2013, 52, 1322−1334