Abstract:
The invention provides a fuzzy curve analysis based soft sensor modeling method using time difference Gaussian process regression, it is suitable for application in chemical process with time delay characteristics. This method can extract stable delay information from the historical database of process and introduce more relevant modeling data sequence to the dominant variable sequence. First of all, the method of fuzzy curve analysis (FCA) can intuitively judge the importance of the input sequence to the output sequence, estimate the time-delay parameters of process, and such offline time-delay parameter set can be utilized to restructure the modeling data. For the new input data, based on the historical variable value before a certain time, the current dominant value can be predicted by time difference Gaussian Process Regression (TDGPR) model. This method does not encounter the problem of model updating and can effectively track the drift between input and output data. Compared with steady-state modeling methods, this invention can achieve more accurate predictions of the key variable, thus improving product quality and reducing production costs.

Description:
CROSS-REFERENCES AND RELATED APPLICATIONS 
       [0001]    This application claims the benefit of priority to Chinese Application No. 201510541727.5, entitled “A fuzzy curve analysis based soft sensor modeling method using time difference Gaussian process regression”, filed on Aug. 28, 2015, which is herein incorporated by reference in its entirety. 
       BACKGROUND OF THE INVENTION 
       [0002]    Field of the Invention 
         [0003]    The invention relates to a fuzzy curve analysis based time difference Gaussian process regression soft sensor modeling method (FCA-TDGPR), and belongs to the field of complex industrial process modeling and soft sensing. 
         [0004]    Description of the Related Art 
         [0005]    The traditional soft sensor modeling methods mostly consider the characteristics of zero delay, that is, considering the input and output with the same sampling interval and input variable collected at time t corresponds to the t-th dominant variable sample in the database. However, there is a significant time lag between the input data collected by each sensor and the output data obtained through the laboratory analysis or online instrumentation. If we continue using steady-state modeling approaches, the established model will be unable to fully explain the characteristics of the process, and it does not meet the causality of actual process. In order to ensure that the soft sensor model can achieve accurate predictions of the key variables in a long time, it is necessary to take measures to introduce the time delay and dynamic information. 
         [0006]    Essentially speaking, most of the existing time delay estimation methods are trying to find the auxiliary variables which are most closely related to the dominant variables for modeling. When applying such methods to practical applications, it is needed to get a tradeoff between the complexity and the accuracy of the algorithm. In view of the process time delay estimation problem, at present, the existing methods include mutual information (MI) method, correlation coefficient method, etc. The invention adopts fuzzy curve analysis method to introduce variable time delay information into the soft sensor model, and the characteristic of this method is low computational complexity while being easy to understand. It is possible to visually and effectively determine the importance degree of the input variables. 
         [0007]    The performance of soft sensor model needs to be maintained by periodic reconstruction in order to dynamically track and effectively take control of the process. The main methods include Moving Window (MW), iterative approach and Just-in-time Learning (JITL), but these methods often need to update the model frequently. Time difference model not only can deal with the problem of model performance degradation due to time lapse and achieve the effect of tracking process dynamics exactly as model updating it also is able to minimize the likelihood of model reconstruction. 
         [0008]    Up to now, data-driven based soft sensor modeling methods have received more and more attention. Some commonly used methods like partial least squares (PLS) and principal component analysis (PCA). They can describe well the linear relationship between input variables and output variables. And artificial neural networks (ANN), support vector machine (SVM) and least squares support vector machine (LS-SVM) can effectively deal with the nonlinear characteristics of the process. In recent years, as a non-parametric probabilistic model, GPR not only can give the predictive value, but also can get the uncertain degree of predictive value. Therefore, the GPR model is selected as the basic soft sensor modeling algorithm in the present invention, and combined with the TD modeling strategy to effectively deal with the drift between input and output data in the process. 
         [0009]    In summary, considering the establishment of a time delay online soft sensor model for the strict control of key variables in the process is of great significance. 
       DETAILED DESCRIPTION 
       [0010]    For time delay process, the soft sensor model without delay can&#39;t be modeled by the sequence of auxiliary variables which are the most relevant to the dominant variable. When such a model is employed to do estimation, the estimation accuracy of the dominant variable will be greatly affected. In order to effectively extract the process delay information and set up an online soft sensor model without the case of frequent model updating, it is necessary to provide a more effective online strategy for the prediction of key variable. Therefore, this invention provides a soft sensor modeling method based on FCA-TDGPR. 
         [0011]    The invention is realized by the following technical scheme: 
         [0012]    Soft sensor modeling method based on FCA-TDGPR which comprises the following steps: aiming at the time delay process, first of all, collect enough samples of data and constitute the historical database with the historical values of auxiliary variables and dominant variable. 
         [0013]    The FCA method is used to determine the time delay parameters of each auxiliary variable with respect to the dominant variable sequence, which is used to reconstruct the soft sensor modeling data; 
         [0014]    Then, When new input samples are available, a time difference Gaussian process regression (TDGPR) model is employed for current time online prediction based on historical variable values collected certain moments ago, thus making it possible to realize real-time estimation and control of the key variable, obtain more accurate results, increase the yield and reduce production costs. 
     
    
     
       BRIEF DESCRIPTION OF DRAWINGS 
         [0015]      FIG. 1  is a flow chart of online soft sensor method based on FCA-TDGPR; 
           [0016]      FIG. 2  is a schematic diagram of debutanizer process; 
           [0017]      FIG. 3  is a schematic diagram of TDGPR modeling approach; 
           [0018]      FIG. 4  contains fuzzy curve distribution diagrams of the original variables and optimal time delay variables; 
           [0019]      FIG. 5  contains scatterplots of butane concentration predictions with different j values. 
       
    
    
     EXAMPLES 
       [0020]    The modeling flow chart, which is shown in  FIG. 1  below, is further detailed in the present invention: 
         [0021]    Take the actual chemical process as an example, debutanizer is an important part of naphtha desulfurization and separation device of oil refining production process, and one of the dominant variables needed to be controlled for this process is the concentration of the bottom butane (C4). The schematic diagram of the process is shown in  FIG. 2 , due to the value of C4 cannot be directly measured, therefore, there is a delay issue in analyzing and obtaining C4 concentration values. At the same time, different auxiliary variables show different degrees of time delay. Experimental data is derived from the actual industrial process which contains 2394 samples, a total of 7 auxiliary variables. As shown in  FIG. 2 , x 1  is the top temperature; x 2  is the top pressure; x 3  is the reflux flow; x 4  is the top product outflow; x 5  is the 6th tray temperature; x 6  is the bottom temperature 1; x 7  is the bottom temperature 2. The 1 dominant variable is the bottom butane concentration, which in this invention is predicted as the key variable of the process. 
         [0022]    Step1: Collect historical input and output data to form a training database which contains N continuous samples. Assuming that the data is expressed as {X(t),y(t)},t 1, 2, . . . , N and is preprocessed, and the 2 bottom temperature variables are averaged as 1 auxiliary variable, then, X(t)=[x 1 (t), x 2 (t), x 3 (t), x 4 (t), x 5 (t), x 6 (t)] T . The maximum time delay T max  of 6 variables is set to 19. 
         [0023]    Step2: For each of the original variables x i , i∈{1, 2, . . . , 6}, they are extended to the input variables with time delay {x i (t−λ), λ=0, 1, . . . , T max } by formula (1), and a set of 120 dimensional delay variables will be obtained for subsequent analysis 
         [0000]    
       
                 
         
             
             
         
       
     
         [0024]    Step3: Determine the importance of each variable in the time delay input variable set by FCA, for (x i (t),y(t)), fuzzy membership function of variable x i  is defined as: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       Φ 
                       it 
                     
                      
                     
                       ( 
                       
                         x 
                         i 
                       
                       ) 
                     
                   
                   = 
                   
                     exp 
                      
                     
                       [ 
                       
                         - 
                         
                           
                             ( 
                             
                               
                                 
                                   
                                     x 
                                     i 
                                   
                                    
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   x 
                                   i 
                                 
                               
                               b 
                             
                             ) 
                           
                           2 
                         
                       
                       ] 
                     
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
     
         [0025]    For each x i , {Φ it , y(t)} provides a fuzzy rule which is described as {if x i  is Φ it (x i ), then y is y(t)}, and Φ it  is a fuzzy membership function of input variable x i  at t-th data point; In formula (2) a Gaussian fuzzy membership function is selected; b is determined as 20% the range of variable x i . As a result, For N training samples, each sample corresponding to each variable has N fuzzy rules. In the fuzzy membership function, Φ it =1 holds true at each point {x i (t),y(t)}. 
         [0026]    For time delay process, by introducing time delay information the original variable x i  becomes (T max +1)-dimensional, which can be expressed as x i (t−λ), λ=0, 1, . . . , T max , λ is a variable delay value to be introduced; fuzzy curve C i,λ  with the condition that λ is the i-th variable delay value can be obtained by making centroid defuzzification of each new expanded variable using formula (3); as shown in the formula (4), d i  is the λ which can make the maximum coverage of fuzzy curve C i,λ ; C i,λ (λ) max  is the maximum value of the fuzzy curve point range, while C i,λ (λ) min  is the minimum value of the fuzzy curve point range; 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       C 
                       
                         i 
                         , 
                         λ 
                       
                     
                      
                     
                       ( 
                       λ 
                       ) 
                     
                   
                   = 
                   
                     
                       
                         ∑ 
                         
                           t 
                           = 
                           1 
                         
                         N 
                       
                        
                       
                         
                           
                             Φ 
                             it 
                           
                            
                           
                             [ 
                             
                               
                                 x 
                                 i 
                               
                                
                               
                                 ( 
                                 
                                   t 
                                   - 
                                   λ 
                                 
                                 ) 
                               
                             
                             ] 
                           
                         
                         · 
                         
                           y 
                            
                           
                             ( 
                             t 
                             ) 
                           
                         
                       
                     
                     
                       
                         ∑ 
                         
                           t 
                           = 
                           1 
                         
                         N 
                       
                        
                       
                         
                           Φ 
                           it 
                         
                          
                         
                           [ 
                           
                             
                               x 
                               i 
                             
                              
                             
                               ( 
                               
                                 t 
                                 - 
                                 λ 
                               
                               ) 
                             
                           
                           ] 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
             
               
                 
                   
                     d 
                     i 
                   
                   = 
                   
                     
                       argmax 
                       λ 
                     
                      
                     
                       [ 
                       
                         
                           
                             
                               C 
                               
                                 i 
                                 , 
                                 λ 
                               
                             
                              
                             
                               ( 
                               λ 
                               ) 
                             
                           
                           max 
                         
                         - 
                         
                           
                             
                               C 
                               
                                 i 
                                 , 
                                 λ 
                               
                             
                              
                             
                               ( 
                               λ 
                               ) 
                             
                           
                           min 
                         
                       
                       ] 
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
         [0027]    If the scope of the C i,λ (λ) range is closer to that of y, then the input variable x i (t−λ) is more important. In view of this point, the importance degree of each variable can be determined by sorting the coverage of C i,λ (λ). Finally, the optimal delay parameter d i  as well as time delay variable x i (t−d i ) can thus be obtained by FCA method, which later on can be used for soft sensor modeling data reconstruction. 
         [0028]    Step4: Based on the previous step, the time delay parameters d 1 , d 2 , . . . , d m  are used to reconstruct the training input sample set for on-line modeling, the reconstructed input dataset is denoted as X d (t), X d (t)[x 1 (t−d 1 ), x 2 (t−d 2 ), x 3 (t−d 3 ), x 4 (t−d 4 ), x 5 (t−d 5 ), x 6 (t−d 6 )]. If there is a new input sample X(t+1), then the delay input set could be restructured based on historical database samples with the same parameters, then go to step 5, otherwise, wait for the arrival of new data. 
         [0029]    Step 5: After the reorganization procedure, the training set and the new data are processed by j order time difference treatment (the value of j can be determined according to the sampling period and property of dominant variable): 
         [0000]      Δ X   d,j ( t ) X   d ( t )− X   d ( t−j )
 
         [0000]      Δ y   j ( t )= y ( t )− y ( t−j )  (5)
 
         [0030]    Next, make a regression of the relationship between ΔX d,j (t) and Δy j (t) by GPR, which satisfies Δ(t)=f(ΔX j (t))+e(t). The GPR method can obtain the mapping relationship through the given training input and output samples. In this way, the corresponding predictive value and the uncertainty degree can be obtained given the new input data, which means the result will be probabilistic. The GPR algorithm is shown as below. 
         [0031]    In general, the relationship between the observed output value y and noise e satisfies: 
         [0000]        y   i   =f ( x   i )+ε
 
         [0000]      ε˜ N (0,σ n   2 )  (6)
 
         [0032]    If the mean function and covariance function are determined, then the distribution of the Gaussian process is well-determined. For simplicity, the mean function is usually preprocessed into 0. Covariance function can transform the correlation of output data into the function of input data. As similar inputs produce similar outputs, the covariance function can be selected according to the characteristics of the sample distribution. One condition which must be satisfied is that the closer the distance of samples is, the more correlated the two samples are, and vise versa. The covariance function form of this invention is shown in formula (7): 
         [0000]    
       
         
           
             
               
                 
                   
                     k 
                      
                     
                       ( 
                       
                         
                           x 
                           p 
                         
                         , 
                         
                           x 
                           q 
                         
                       
                       ) 
                     
                   
                   = 
                   
                     v 
                      
                     
                         
                     
                      
                     
                       exp 
                        
                       
                         [ 
                         
                           
                             - 
                             
                               1 
                               2 
                             
                           
                            
                           
                             
                               ∑ 
                               
                                 d 
                                 = 
                                 1 
                               
                               D 
                             
                              
                             
                                 
                             
                              
                             
                               
                                 
                                   π 
                                   d 
                                 
                                  
                                 
                                   ( 
                                   
                                     
                                       x 
                                       p 
                                       d 
                                     
                                     - 
                                     
                                       x 
                                       q 
                                       d 
                                     
                                   
                                   ) 
                                 
                               
                               2 
                             
                           
                         
                         ] 
                       
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
         [0033]    In the formula, x p , x q ∈R D , v controls the magnitude of the covariance function. π d  describes the relative importance of each input attribute x d . The determination of the hyper-parameter Θ gp =(v, π 1 , . . . , π D , σ R   2 ) in the Gaussian process is generally estimated by the MLE method. The optimization of the parameters can be realized by using the conjugate gradient method. Based on test sample and training data, the posterior distribution of test data x *  can be calculated, and its predictive value obey the joint Gaussian distribution described in formula (9), where K(X,X) is n-dimensional covariance matrix of training samples; k(x * ,X) is the covariance vector of test sample and training samples; k(x * ,x * ) is the autocovariance of test sample, and f gp  is a predictive value of GPR. 
         [0000]    
       
         
           
             
               
                 
                   
                     L 
                      
                     
                       ( 
                       
                         Θ 
                         gp 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         - 
                         
                           1 
                           2 
                         
                       
                        
                       
                         
                           
                             y 
                             T 
                           
                            
                           
                             [ 
                             
                               
                                 K 
                                  
                                 
                                   ( 
                                   
                                     X 
                                     , 
                                     X 
                                   
                                   ) 
                                 
                               
                               + 
                               
                                 
                                   σ 
                                   n 
                                   2 
                                 
                                  
                                 I 
                               
                             
                             ] 
                           
                         
                         
                           - 
                           1 
                         
                       
                        
                       y 
                     
                     - 
                     
                       
                         1 
                         2 
                       
                        
                       log 
                        
                       
                           
                       
                        
                       
                         det 
                          
                         
                           [ 
                           
                             
                               K 
                                
                               
                                 ( 
                                 
                                   X 
                                   , 
                                   X 
                                 
                                 ) 
                               
                             
                             + 
                             
                               
                                 σ 
                                 n 
                                 2 
                               
                                
                               I 
                             
                           
                           ] 
                         
                       
                     
                     - 
                     
                       
                         n 
                         2 
                       
                        
                       log 
                        
                       
                           
                       
                        
                       2 
                        
                       π 
                     
                   
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
             
               
                 
                   
                       
                   
                    
                   
                     
                       
                         
                           f 
                           gp 
                         
                          
                         X 
                       
                       , 
                       y 
                       , 
                       
                         
                           x 
                           * 
                         
                         ~ 
                         
                           N 
                            
                           
                             ( 
                             
                               
                                 
                                   f 
                                   gp 
                                 
                                 _ 
                               
                               , 
                               
                                 cov 
                                  
                                 
                                   ( 
                                   
                                     f 
                                     gp 
                                   
                                   ) 
                                 
                               
                             
                             ) 
                           
                         
                       
                     
                      
                     
                         
                     
                      
                     
                       
 
                     
                      
                     
                         
                     
                      
                     
                       
                         s 
                         . 
                         t 
                         . 
                         
                             
                         
                          
                         
                           
                             f 
                             gp 
                           
                           _ 
                         
                       
                       = 
                       
                         
                           
                             
                               k 
                                
                               
                                 ( 
                                 
                                   
                                     x 
                                     * 
                                   
                                   , 
                                   X 
                                 
                                 ) 
                               
                             
                              
                             
                               [ 
                               
                                 
                                   K 
                                    
                                   
                                     ( 
                                     
                                       X 
                                       , 
                                       X 
                                     
                                     ) 
                                   
                                 
                                 + 
                                 
                                   
                                     σ 
                                     n 
                                     2 
                                   
                                    
                                   
                                     I 
                                     n 
                                   
                                 
                               
                               ] 
                             
                           
                           
                             - 
                             1 
                           
                         
                          
                         y 
                       
                     
                      
                     
                       
 
                     
                      
                     
                         
                     
                      
                     
                       
                         cov 
                          
                         
                           ( 
                           
                             f 
                             gp 
                           
                           ) 
                         
                       
                       = 
                       
                         
                           k 
                            
                           
                             ( 
                             
                               
                                 x 
                                 * 
                               
                               , 
                               
                                 x 
                                 * 
                               
                             
                             ) 
                           
                         
                         - 
                         
                           
                             
                               
                                 
                                   k 
                                    
                                   
                                     ( 
                                     
                                       X 
                                       , 
                                       
                                         x 
                                         * 
                                       
                                     
                                     ) 
                                   
                                 
                                 T 
                               
                                
                               
                                 [ 
                                 
                                   
                                     K 
                                      
                                     
                                       ( 
                                       
                                         X 
                                         , 
                                         X 
                                       
                                       ) 
                                     
                                   
                                   + 
                                   
                                     
                                       σ 
                                       n 
                                       2 
                                     
                                      
                                     
                                       I 
                                       n 
                                     
                                   
                                 
                                 ] 
                               
                             
                             
                               - 
                               1 
                             
                           
                           · 
                           
                             k 
                              
                             
                               ( 
                               
                                 X 
                                 , 
                                 
                                   x 
                                   * 
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
         [0034]    when the new input data arrives at time t+1, the calculating formula of predictive value y j,pred (t+1) with TDGPR method is: 
         [0000]      Δ X ( t+ 1)= X ( t+ 1)− X ( t+ 1− j )
 
         [0000]      Δ y   j,pred ( t+ 1)= f   GPR (Δ X   j ( t+ 1))
 
         [0000]        y   j,pred ( t+ 1)= y   j ( t+ 1− j )+Δ y   j,pred ( t+ 1)  (10)
 
         [0035]    In the actual industrial process, there will be the case of instrument damage or laboratory analysis with delay, and the circumstance that the time interval of obtaining dominant variable is large and the quantity is small or there is a lack of dominant analysis value in the database. Thus, shown in  FIG. 3 , for new incoming test data X(t+1), based on y i (t+1−j) stored in the database j moment ago, the predictive value of the dominant variable at time t+1 can be obtained. The predicted output y j,pred (t+1) of the online model is calculated by formula (10), and the predicted result of bottom butane concentration can be obtained. 
         [0036]    As shown in  FIG. 4 , compared the original variables without time delay, 6 reconstructed variables contribute more to the dominant variable, which introduce more relevant modeling data for online modeling. At the same time, in order to verify the effectiveness of this invention for on-line estimation, the first 1519 samples are selected in 2394 samples to reconstruct 1500 training samples. The final 875 samples are then used as test samples, and a soft sensor is established for on-line prediction of butane concentration. 
         [0037]      FIG. 5  contains scatterplots of butane concentration prediction results with two methods respectively denoted as the FCA-TDGPR method (present method) which involves delay estimation, and t-TDGPR method which is without time delay estimation. 
         [0038]    From  FIG. 5 , when the time difference order j increases from 1 to 10, the time interval of prediction based on historical database is gradually increasing, and the prediction accuracy is declining. This is because the more recent the analysis value is, the better tracking ability the model has for current process dynamics. 
         [0039]    Although the accuracy of the two methods are in decline, compared with the TDGPR method without considering the time delay, the predicted results of the present invention can be better close to the true value of the butane concentration when the time difference increases. This suggests that extracted delay information is in line with the actual causal relationship of the process, and the soft sensor model with variable time delay estimation is more accurate. 
         [0040]    After fuzzy curve analysis method is taken to determine the optimal parameters, reconstructed data is proved to be capable of enhancing the accuracy of online model significantly by introducing more contributing auxiliary variables to dominant variable sequence. At the same time, it reflects that the GPR method can explain the dynamic change of the process well. The online soft sensor model based on TDGPR method can adaptively estimate real-time butane concentration with historical variables collected j time ago. 
         [0041]      FIGS. 4 and 5  have jointly validated that fuzzy curve analysis based time difference Gaussian process regression soft sensor modeling method has good accuracy for on-line prediction of bottom butane concentration. 
         [0042]    While the present invention has been described in some detail for purposes of clarity and understanding, one skilled in the art will appreciate that various changes in form and detail can be made without departing from the true scope of the invention. All figures, tables, appendices, patents, patent applications and publications, referred to above, are hereby incorporated by reference.