Patent Publication Number: US-7899652-B2

Title: Linear programming support vector regression with wavelet kernel

Description:
FIELD OF THE INVENTION 
     The present invention relates to nonlinear dynamical systems identification and modeling, in particular using linear programming support vector regression (LP-SVR) with wavelet kernel. 
     BACKGROUND OF THE INVENTION 
     Models of dynamical systems are of great importance in almost all fields of science and engineering and specifically in control, signal processing and information science. Most systems encountered in the real world are nonlinear and in many practical applications nonlinear models are required to achieve an adequate modeling accuracy. A model is always only an approximation of a real phenomenon so that having an approximation theory which allows for the analysis of model quality is a substantial concern. 
     Wavelet theory has been extensively studied in recent years, and has found many applications in various areas throughout science and engineering, such as numerical analysis and signal processing. Wavelets can capture the local behavior of signals both in frequency and time. 
     SUMMARY OF THE INVENTION 
     A support vector machine with wavelet kernel was developed for accurate modeling of nonlinear systems. A method of providing an optimized model of a nonlinear system includes using a support vector machine (SVM) having a wavelet kernel, where support vectors include a family of multidimensional wavelets. Training the SVM allows optimization of the number of support vectors, the weights of the support vectors, and the translation factors of the support vectors. Use of a novel linear programming approach reduces computational demands required for training, allowing optimized support vectors to give an optimized model of the nonlinear system. Further, in many applications on-line retraining is possible, so that the model can be adapted to changing conditions, new system parameters, and the like. 
     Applications include an electronic control unit (ECU) of a vehicle engine, the SVM providing a model of engine performance. The model can be changed as the engine ages, or other systems variations detected. Comparison of the model predictions with actual performance, which can be performed on-line using a sliding window of data, allow system (engine) changes to be detected. The ECU may be reprogrammed through retraining of the SVM using collected engine data. Diagnostic outputs may be obtained from a comparison of the model and actual engine performance. 
     Optimized support vectors for an engine type may be determined during an engine test process, and used as the initial configuration of the SVM in the ECU for engines of the same type in new vehicles. However, for each vehicle engine, the SVM may be optimized further by retraining on actual data from the individual engine. Hence, the performance of individual vehicle engines can be separately optimized. 
     By collecting system data from a nonlinear system, a predicted system performance using the support vector machine can be compared with actual performance so as to determine modeling error. If the error is greater than a predetermined threshold, the support vector machine may be retrained, for example using a computationally simpler linear programming method, to re-optimize feature vector parameters using the collected system data. Hence, an optimized model of the nonlinear system is maintained even if the system changes over time, as the system model can be re-optimized. Feature vector parameters may be optimized using a simple algorithm, such as a least squares algorithm, if the modeling error is less than a threshold value. 
     An improved electronic control unit for a vehicle engine comprises engine data inputs configured to receive engine data, the engine data including engine operation data and operator input data; engine control outputs, operable to control engine parameters (such as fuel injection quantity and valve timing), and an electronic control circuit. The electronic control circuit comprises at least a memory and a processor, the electronic circuit operable to use a support vector machine with a wavelet kernel as a model of the vehicle engine, the model being used to determine required control outputs from the engine data. The electronic control circuit may further be operable to determine a modeling error by comparing a modeled engine performance with an actual system performance, and to re-train the support vector machine if the modeling error is greater than a threshold value. Retraining the support vector machine may use a linear programming method to determine feature vector parameters from collected engine data, such as coolant temperature data, throttle position data, and the like. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  shows identification by LP-SVR with RBF kernel on a training set, where the solid line is observation, and the dotted line a model from LP-SVR; 
         FIG. 2  shows identification by LP-SVR with RBF kernel on a validation set, where the solid line is observation, and the dotted line a model from LP-SVR; 
         FIG. 3  shows identification by QP-SVR with RBF kernel on the training set, where the solid line is observation, and the dotted line a model from QP-SVR; 
         FIG. 4  shows identification by QP-SVR with RBF kernel on the validation set, where the solid line is observation, and the dotted line a model from QP-SVR; 
         FIG. 5  shows identification by LP-SVR with wavelet kernel on the training set, where the solid line is observation, and the dotted line a model from LP-SVR; 
         FIG. 6  shows identification by LP-SVR with wavelet kernel on the validation set, where the solid line is observation, and the dotted line a model from LP-SVR; 
         FIG. 7  shows identification by QP-SVR with wavelet kernel on the training set, where the solid line is observation, and the dotted line a model from QP-SVR; 
         FIG. 8  shows identification by QP-SVR with wavelet kernel on the validation set, where the solid line is observation, and the dotted line a model from QP-SVR; and 
         FIG. 9  is a simplified schematic of an engine control system. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     A support vector machine allows modeling of complex dynamic systems. The kernel of the support vector machine may be of various types. In examples of the present invention, the support vector machine has a wavelet kernel, which is well suited to modeling nonlinear systems, particularly to modeling dynamic nonlinear systems with transient behaviors. 
     Wavelet kernels improve accuracy over conventional SVM kernels, as a wavelet kernel is especially suitable for local signal analysis, signal-noise separation and detection of transients and jumping signals. The wavelet kernel also improves model sparsity, readily providing approximate orthogonality. 
     Examples of the present invention use a wavelet kernel in linear programming support vector regression for nonlinear dynamical systems identification. A novel approach was used to construct the wavelet networks by linear programming, allowing excellent generalization capability inherent in support vector learning and compact model expression. A linear programming support vector regression (SVR) was formulated, increasing computational efficiency and model sparsity over conventional approaches. 
     Applications include all types of nonlinear system modeling, and the general approach described has widespread applications including electronic engine control in automobiles, image compression, speech signal processing, process control generally, and other complex systems. 
     Support vector kernel expansion allows computational efficient representation of nonlinear dynamical systems and facilitates analysis. Using the support vector regression with the conventional quadratic programming technique, implementation is computationally expensive and sufficient model sparsity cannot be guaranteed. 
     In examples of the present invention, linear programming support vector regression (LP-SVR) with wavelet kernel is used for system modeling. Relationships between LP-SVR with wavelet kernel and wavelet networks are analyzed, in relation to nonlinear dynamical systems identification. Wavelet networks were constructed by linear programming, with excellent generalization capability inherent in support vector learning and compact model expression. 
     Wavelet networks can be used in various classification and identification problems Wavelet networks facilitate detecting essential features in “frequency-rich” signals, including transient data. In wavelet networks, both the position (translation parameter) and the dilation of the wavelets may be optimized, in addition to the weights accorded each wavelet in the network. A family of wavelets can be constructed by translations and dilations performed on a single fixed function called the mother wavelet. Each wavelet has a translation factor m and dilation factor d in relation to the mother wavelet. Families of wavelet functions are known to be universal approximators, allowing them to be used in function approximation and process modeling. Any function of the type L 2 (R) can be approximated to any prescribed accuracy with a finite sum of wavelets, each wavelet having a translation factor, dilation factor, and a weight. Nonlinear regression by wavelet networks can be used for both static modeling and for dynamic input-output modeling. 
     A data set can be modeled by a finite sum of wavelets, in particular a weighted sum of wavelets functions whose translation and dilation parameters (m j  and d j  respectively, see Equation 1) are adjustable real numbers, and these parameters are trained together with the weights. In this way, a wavelet family of parameterized nonlinear functions can be used for nonlinear regression, and their parameters are estimated by a training step. 
     Further Discussion of Wavelet Networks 
     Wavelet families are constructed by simple operations of translation and dilation of a single fixed function called the mother wavelet (sometimes called the analyzing wavelet). A wavelet φ j (x) is derived from its mother wavelet φ(z) by the relation 
                       ϕ   j     ⁡     (   x   )       =       ϕ   (       x   -     m   j         d   j       )     =     ϕ   ⁡     (     z   j     )                 (   1   )               
where the translation factor m j  and the dilation factor d j  are real numbers in   and    + *, respectively. The family of functions generated by φ can be defined as
 
     
       
         
           
             
               
                 
                   
                     Ω 
                     c 
                   
                   = 
                   
                     { 
                     
                       
                         
                           1 
                           
                             
                               d 
                               j 
                             
                           
                         
                         ⁢ 
                         
                           ϕ 
                           ( 
                           
                             
                               x 
                               - 
                               
                                 m 
                                 j 
                               
                             
                             
                               d 
                               j 
                             
                           
                           ) 
                         
                       
                       , 
                       
                         
                           m 
                           j 
                         
                         ∈ 
                         
                           ℜ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           and 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             d 
                             j 
                           
                         
                         ∈ 
                         
                           ℜ 
                           + 
                           * 
                         
                       
                     
                     } 
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
     
     A family Ω c  is said to be a frame of L 2 (R) if there exists two constants A and B, 0&lt;A≦B&lt;∞, such that for any square integrable function ƒ, the following inequalities hold: 
                     A   ⁢          f        2       ≤       ∑         j               ϕ   j     ∈     Ω   c               ⁢            〈       ϕ   j     ,   f     〉          2       ≤     B   ⁢          f        2               (   3   )               
where ∥ƒ∥ denotes the norm of function ƒ and  ƒ,g  the inner product of functions ƒ and g. Families of wavelet frames of L 2 (R) are universal approximators.
 
     For the modeling of multivariable processes, a multi-dimensional wavelet is defined. In examples described further below, a dictionary composed of tensor product wavelet functions was used, and the elements of the wavelet dictionaries were of the form 
                       Φ   j     ⁡     (     [           x   1             ⋮             x   n           ]     )       =       ∏     k   =   1     n     ⁢     ϕ   (         x   k     -     m   jk         d   jk       )               (   4   )               
where m j =[m j1 , m j1 , . . . , m jn ] T  and d j =[d j1 , d j1 , . . . , d jn ] T  are the translation and dilation vectors, respectively.
 
     Multi-dimensional wavelet networks can be written as follows: 
                         f   p     ⁡     (   x   )       =       ∑     j   =   1     p     ⁢       w   j     ⁢       Φ   j     ⁡     (   x   )             ,           (   5   )               
where x=[x 1 , x 2 , . . . , x n ] T  and w j  are the weights of the network.
 
Linear Programming SVM with Wavelet Kernel
 
     An example linear programming SVM with wavelet kernel was developed, using soft-constrained linear programming SVM. Both LP-SVR and QP-SVR algorithms may adopt the ε-insensitive loss function, and use kernel functions in feature space. 
     Consider regression in the following set of functions
 
ƒ( x )= w   T  φ( x )+ b    (6)
 
with given training data, {(x 1 ,y 1 ), . . . , (x l ,y l )} where l denotes the total number of exemplars, x i  ε R n  are the input and y i  ε R are the target output data. The nonlinear mapping φ: R n →R m  (m&gt;n) maps the input data into a so-called high dimensional feature space (which can be infinite dimensional) and w ε R m , b ε R. In ε-SV regression, the goal is to find a function ƒ(x) that has at most ε deviation from the actually obtained targets y i  for all the training data, and at the same time, is as smooth as possible.
 
     In a support vector method, the empirical risk is minimized subject to elements of the structure: 
                       minimize   ⁢           ⁢     1   2     ⁢          w        2       +     C   ⁢       ∑     l   =   1     ℓ     ⁢     (       ξ   i     +     ξ   i   *       )           ⁢     
     ⁢     subject   ⁢           ⁢   to   ⁢           ⁢     {               y   i     -     〈     w   ,     φ   ⁡     (     x   i     )         〉     -   b     ≤     ɛ   +     ξ   i                       〈     W   ,     φ   ⁡     (     x   i     )         〉     +   b   -     y   i       ≤     ɛ   +     ξ   i   *                     ξ   i     ,       ξ   i   *     ≥   0                         (   7   )               
where the ξ i  and ξ i * are the slack variables, corresponding to the size of the excess deviation for positive and negative deviation respectively. This is a classic quadratic optimization problem with inequality constraints, and the optimization criterion penalizes data points whose y-values differ from ƒ(x) by more than ε. The constant C&gt;0 determines the trade-off between the flatness of ƒ and the amount up to which deviations larger than ε are tolerated. By defining the loss function,
 
                     L   ⁡     (       y   i     -     f   ⁡     (     x   i     )         )       =     {           0   ,             if   ⁢           ⁢            y   i     -     f   ⁡     (     x   i     )                ≤   ɛ                          y   i     -     f   ⁡     (     x   i     )              -   ɛ     ,         otherwise                   (   8   )               
the optimization problem (7) is equivalent to the following regularization problem,
 
                     minimize   ⁢           ⁢       R   reg     ⁡     [   f   ]         =         ∑     l   =   1     ℓ     ⁢     L   ⁡     (       y   i     -     f   ⁡     (     x   i     )         )         +     λ   ⁢          w        2                 (   9   )               
where ƒ(x) is in the form of (6) and λ∥w∥ 2  is the regularization term. According to the well-known Representer Theorem, the solution to the regularization problem (9) can be written as the following SV kernel expansion provided k(x i ,x i )=1
 
                     f   ⁡     (   x   )       =       ∑     l   =   1     ℓ     ⁢       β   i     ⁢     k   ⁡     (       x   i     ,   x     )                   (   10   )               
where k(x i ,x) is the kernel function. Three commonly-used kernel functions in literature are:
         Gaussian radial basis function (GRBF) kernel:       

     
       
         
           
             
               
                 
                   
                     k 
                     ⁡ 
                     
                       ( 
                       
                         x 
                         , 
                         
                           x 
                           ′ 
                         
                       
                       ) 
                     
                   
                   = 
                   
                     exp 
                     ( 
                     
                       
                         - 
                         
                           
                              
                             
                               x 
                               - 
                               
                                 x 
                                 ′ 
                               
                             
                              
                           
                           2 
                         
                       
                       
                         σ 
                         2 
                       
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
           
         
       
         
         
           
             Polynomial kernel:
 
 k ( x,x ′)=(1+   x, x ′ ) 9    (12)
 
             Sigmoid kernel:
 
 k ( x,x ′)=tan  h (α   x, x ′ +γ)   (13)
 
where σ, q, α, γ are the adjustable parameters of the above kernel functions. The kernel function provides an elegant way of working in the feature space avoiding all the troubles and difficulties inherent in high dimensions, and this method is applicable whenever an algorithm can be cast in terms of dot products. Defining
 
β=[β 1  β 2 . . . β l ] T ,   (14)
 
LP-SVR replaces (9) by
 
           
         
       
    
                     minimize   ⁢           ⁢       R   reg     ⁡     [   f   ]         =         ∑     l   =   1     ℓ     ⁢     L   ⁡     (       y   i     -     f   ⁡     (     x   i     )         )         +     λ   ⁢          β        1                 (   15   )               
where ƒ(x) is in the form of (10) and ∥β∥ 1  denotes the l 1  norm in coefficient space. This regularization problem is equivalent to the following constrained optimization problem
 
     
       
         
           
             
               
                 
                   
                     
                       minimize 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         1 
                         2 
                       
                       ⁢ 
                       
                         
                            
                           β 
                            
                         
                         1 
                       
                     
                     + 
                     
                       C 
                       ⁢ 
                       
                         
                           ∑ 
                           
                             i 
                             = 
                             1 
                           
                           ℓ 
                         
                         ⁢ 
                         
                           ( 
                           
                             
                               ξ 
                               i 
                             
                             + 
                             
                               ξ 
                               i 
                               * 
                             
                           
                           ) 
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     subject 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     to 
                     ⁢ 
                     
                       
                           
                       
                       ⁢ 
                       
                           
                       
                     
                     ⁢ 
                     
                       { 
                       
                         
                           
                             
                               
                                 
                                   y 
                                   i 
                                 
                                 - 
                                 
                                   
                                     ∑ 
                                     
                                       j 
                                       = 
                                       1 
                                     
                                     ℓ 
                                   
                                   ⁢ 
                                   
                                     
                                       β 
                                       j 
                                     
                                     ⁢ 
                                     
                                       k 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             x 
                                             j 
                                           
                                           , 
                                           
                                             x 
                                             i 
                                           
                                         
                                         ) 
                                       
                                     
                                   
                                 
                               
                               ≤ 
                               
                                 ɛ 
                                 + 
                                 
                                   ξ 
                                   i 
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 
                                   
                                     ∑ 
                                     
                                       j 
                                       = 
                                       1 
                                     
                                     ℓ 
                                   
                                   ⁢ 
                                   
                                     
                                       β 
                                       j 
                                     
                                     ⁢ 
                                     
                                       k 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             x 
                                             j 
                                           
                                           , 
                                           
                                             x 
                                             i 
                                           
                                         
                                         ) 
                                       
                                     
                                   
                                 
                                 - 
                                 
                                   y 
                                   i 
                                 
                               
                               ≤ 
                               
                                 ɛ 
                                 + 
                                 
                                   ξ 
                                   i 
                                   * 
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 ξ 
                                 i 
                               
                               , 
                               
                                 
                                   ξ 
                                   i 
                                   * 
                                 
                                 ≥ 
                                 0 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   16 
                   ) 
                 
               
             
           
         
       
     
     From the geometric perspective, it can be followed that ξ iξ   i *=0 in SV regression. Therefore, it is sufficient to just introduce slack ξ i  in the constrained optimization problem (16). Thus, the following formulation of SV regression is obtained, with fewer slack variables 
     
       
         
           
             
               
                 
                   
                     
                       minimize 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         1 
                         2 
                       
                       ⁢ 
                       
                         
                            
                           β 
                            
                         
                         1 
                       
                     
                     + 
                     
                       C 
                       ⁢ 
                       
                         
                           ∑ 
                           
                             i 
                             = 
                             1 
                           
                           ℓ 
                         
                         ⁢ 
                         
                           ( 
                           
                             
                               ξ 
                               i 
                             
                             + 
                             
                               ξ 
                               i 
                               * 
                             
                           
                           ) 
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     subject 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     to 
                     ⁢ 
                     
                       
                           
                       
                       ⁢ 
                       
                           
                       
                     
                     ⁢ 
                     
                       { 
                       
                         
                           
                             
                               
                                 
                                   y 
                                   i 
                                 
                                 - 
                                 
                                   
                                     ∑ 
                                     
                                       j 
                                       = 
                                       1 
                                     
                                     ℓ 
                                   
                                   ⁢ 
                                   
                                     
                                       β 
                                       j 
                                     
                                     ⁢ 
                                     
                                       k 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             x 
                                             j 
                                           
                                           , 
                                           
                                             x 
                                             i 
                                           
                                         
                                         ) 
                                       
                                     
                                   
                                 
                               
                               ≤ 
                               
                                 ɛ 
                                 + 
                                 
                                   ξ 
                                   i 
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 
                                   
                                     ∑ 
                                     
                                       j 
                                       = 
                                       1 
                                     
                                     ℓ 
                                   
                                   ⁢ 
                                   
                                     
                                       β 
                                       j 
                                     
                                     ⁢ 
                                     
                                       k 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             x 
                                             j 
                                           
                                           , 
                                           
                                             x 
                                             i 
                                           
                                         
                                         ) 
                                       
                                     
                                   
                                 
                                 - 
                                 
                                   y 
                                   i 
                                 
                               
                               ≤ 
                               
                                 ɛ 
                                 + 
                                 
                                   ξ 
                                   i 
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 ξ 
                                 i 
                               
                               ≥ 
                               0 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   17 
                   ) 
                 
               
             
           
         
       
     
     The optimization problem above may be converted into a linear programming problem. The terms β i  and |β i | are decomposed as follows:
 
β i =α i   + −α i   −  |β i |=α i   + +α i   −   (18)
 
where α i   + , α i   − ≧0. The decompositions in (18) are unique, i.e. for a given β i  there is only one pair (α i   + , α i   − ) which fulfils both equations. Furthermore, both variables cannot be larger than zero at the same time, i.e. α i   + ·α i   − =0. In this way, the l 1  norm of β can be written as
 
                          β        1     =       (         1   ,   1   ,   …   ⁢           ,   1       ︸   ℓ       ,       1   ,   1   ,   …   ⁢           ,   1       ︸   ℓ         )     ⁢     (           α   +               α   -           )               (   19   )               
where α + =(α 1   30  , α 2   + , . . . , α l   + ) T  and α − =(α 1   − , α 2   − , . . . , α l   − ) T . Furthermore, the constraints in the formulation (17) can also be written in the following vector form
 
                       (         K         -   K           -   I               -   K         K         -   I           )     ·     (           α   +               α   -             ξ         )       ≤     (           y   +   ɛ               ɛ   -   y           )             (   20   )               
where K ij =k(x i , x j ), ξ=(ξ 1 , ξ 2 , . . . , ξ l ) T  and I is l×l identity matrix. Thus, the constrained optimization problem (17) can be implemented by the following linear programming problem with the variables
 
     
       
         
           
             
               
                 
                   
                     minimize 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         c 
                         T 
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             
                               
                                 α 
                                 + 
                               
                             
                           
                           
                             
                               
                                 α 
                                 - 
                               
                             
                           
                           
                             
                               ξ 
                             
                           
                         
                         ) 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       subject 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       to 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           ( 
                           
                             
                               
                                 K 
                               
                               
                                 
                                   - 
                                   K 
                                 
                               
                               
                                 
                                   - 
                                   I 
                                 
                               
                             
                             
                               
                                 
                                   - 
                                   K 
                                 
                               
                               
                                 K 
                               
                               
                                 
                                   - 
                                   I 
                                 
                               
                             
                           
                           ) 
                         
                         · 
                         
                           ( 
                           
                             
                               
                                 
                                   α 
                                   + 
                                 
                               
                             
                             
                               
                                 
                                   α 
                                   - 
                                 
                               
                             
                             
                               
                                 ξ 
                               
                             
                           
                           ) 
                         
                       
                     
                     ≤ 
                     
                       ( 
                       
                         
                           
                             
                               y 
                               + 
                               ɛ 
                             
                           
                         
                         
                           
                             
                               ɛ 
                               - 
                               y 
                             
                           
                         
                       
                       ) 
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       where 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       c 
                     
                     = 
                     
                       ( 
                       
                         
                           
                             1 
                             , 
                             1 
                             , 
                             … 
                             ⁢ 
                             
                                 
                             
                             , 
                             1 
                           
                           
                             ︸ 
                             l 
                           
                         
                         , 
                         
                           
                             1 
                             , 
                             1 
                             , 
                             … 
                             ⁢ 
                             
                                 
                             
                             , 
                             1 
                           
                           
                             ︸ 
                             l 
                           
                         
                         , 
                         
                           
                             
                               2 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               C 
                             
                             , 
                             
                               2 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               C 
                             
                             , 
                             … 
                             ⁢ 
                             
                                 
                             
                             , 
                             
                               2 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               C 
                             
                           
                           
                             ︸ 
                             l 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
           
         
       
     
     In the QP-SVR case, the set of points not inside the tube coincides with the set of SVs. While, in the LP context, this is no longer true—although the solution is still sparse, any point could be an SV, even if it is inside the tube. The sparse solution still can be obtained in LP-SVR even though the size of the insensitive tube was set to zero due to the usage of soft constraints; however, usually sparser solution can be obtained by using non-zero ε. 
     Linear Programming SVM with Wavelet Kernel 
     In a nonlinear support vector algorithm, the kernel provides a general framework to represent data. In this example, a wavelet kernel was used in linear programming support vector regression, and this also allows determination of the topological structure and translation parameters of wavelet networks by an efficient optimization approach. From equation (4), the wavelet kernel can be defined as: 
                       k   ⁡     (     x   ,     x   ′       )       =       ∏     k   =   1     n     ⁢           ⁢     ϕ   ⁡     (         x   k     -     x   k   ′       d     )           ,           (   22   )               
and the corresponding SV wavelet kernel expansion can be written as:
 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           f 
                           ⁡ 
                           
                             ( 
                             x 
                             ) 
                           
                         
                         = 
                         
                           
                             ∑ 
                             
                               i 
                               = 
                               1 
                             
                             l 
                           
                           ⁢ 
                           
                             
                               β 
                               i 
                             
                             ⁢ 
                             
                               
                                 ∏ 
                                 
                                   k 
                                   = 
                                   1 
                                 
                                 n 
                               
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 ϕ 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       
                                         x 
                                         k 
                                       
                                       - 
                                       
                                         x 
                                         k 
                                         ′ 
                                       
                                     
                                     d 
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                         
                           
                             ∑ 
                             
                               i 
                               = 
                               1 
                             
                             l 
                           
                           ⁢ 
                           
                             
                               ( 
                               
                                 
                                   α 
                                   i 
                                   + 
                                 
                                 - 
                                 
                                   α 
                                   i 
                                   - 
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               
                                 ∏ 
                                 
                                   k 
                                   = 
                                   1 
                                 
                                 n 
                               
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 
                                   ϕ 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       
                                         
                                           x 
                                           k 
                                         
                                         - 
                                         
                                           x 
                                           k 
                                           ′ 
                                         
                                       
                                       
                                         d 
                                         i 
                                       
                                     
                                     ) 
                                   
                                 
                                 . 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   23 
                   ) 
                 
               
             
           
         
       
     
     The model sparsity in equation (23) is measured by the ratio of non-zero components in the vector β=[β 1  β 2 . . . β] T , i.e., the ratio of support vectors. In the realm of nonlinear systems identification, training data are usually finite and non-uniformly sampled, so the problem is consequently ill-posed. 
     Conversion to a well-posed problem is typically achieved with some form of capacity control, which aims to balance the fitting of the data with constraints on the model flexibility, producing a robust model that generalizes successfully. In practice, such an optimization is accomplished by searching for the minimum number of the basis functions under the Occam&#39;s razor arguing that the model should be no more complex than is required to capture the underlying systems dynamics. Hence, in an attempt to achieve the highest generalization capability and the lowest system complexity, the l 1  norm of the weights in the model (23) was employed for model capacity control. 
     An optimal compact representation, including the number of support vectors, the weights and the translation factors, can be determined by solving the following linear programming problem: 
                     minimize   ⁢           ⁢       c   T     ⁡     (           α   +               α   -             ξ         )         ⁢     
     ⁢       subject   ⁢           ⁢   to   ⁢           ⁢       (         K         -   K           -   I               -   K         K         -   I           )     ·     (           α   +               α   -             ξ         )         ≤     (           y   +   ɛ               ɛ   -   y           )               (   24   )               
where the K denotes the wavelet kernels matrix. Explicitly, the obtained SV wavelet kernel expansion (23) is essentially consistent with the expression of wavelet networks (5), and therefore the construction of the optimal wavelet network is fulfilled by solving the linear programming problem.
 
     EXAMPLES 
     The problem of nonlinear black-box systems identification of a hydraulic robot arm is a benchmark problem in the realm of nonlinear systems identification, and it has been used widely for testing various identification methods. 
     In this nonlinear black-box dynamical system, the input u(t) represents the size of the valve through which oil flows into the actuator, and the output y(t) is a measure of oil pressure which determines the robot arm position. For the purpose of comparison, the same regressor is used, x(t)=[y(t−1) y(t−2) y(t−3) u(t−1) u(t−2)] T , as that used previously in the literature. See A. Gretton et al., “Support vector regression for black-box system identification”, Proceedings of the 11th IEEE Workshop on Statistical Signal Processing (2001), and J. Sjöberg et al., “Nonlinear black-box modeling in system identification: a unified overview”,  Automatica,  31, 1691-1724 (1995). Half of the data set containing 511 training data pairs were used for SVM training, and half as validation data. The generalization capability and accuracy of regression algorithms could be evaluated by making use of the root mean square (RMS) error 
                     E   RMS     =         1   N     ⁢       ∑     i   =   0     N     ⁢       [         f   ^     ⁡     (     x   i     )       -     f   ⁡     (     x   i     )         ]     2                   (   25   )               
where {circumflex over (ƒ)}(x i ) is the estimated value at point x i  from the SVR model.
 
     Example 1 
     Identification of Hydraulic Robot Arm Dynamics by SVR with RBF Kernel 
     In this example, the soft-constrained LP-SVR and QP-SVR with RBF kernel are applied to model the dynamic behavior of a hydraulic robot arm respectively. The same learning parameters for LP-SVR and QP-SVR are used: kernel width parameter—1; C (bound)—100; ε(tolerance)—0.03. 
       FIGS. 1 and 2  show the identification results by LP-SVR, where the RMS error on validation data is 0.1959 and the ratio of SVs to the training vectors is 6.5%. 
       FIGS. 3 and 4  visualize the identification results by QP-SVR, where the RMS error on validation set is 0.1597 and the ratio of SVs to the training vectors is 39.7%. In this example, the prediction accuracy from LP-SVR is comparable with that from QP-SVR, and the LP-SVR is around 25 times faster than QP-SVR for training, which means that the computing resources required by QP-SVR may become prohibitively expensive when the size of the training set is increased. 
     Example 2 
     Identification of Hydraulic Robot Arm Dynamics by SVR with Wavelet Kernel 
     For the sake of validating the performance of wavelet kernel in nonlinear identification, the soft-constrained LP-SVR and QP-SVR with wavelet kernel are used to model the dynamic behavior of a hydraulic robot arm respectively. Without loss of generality, a translation-invariant wavelet kernel is constructed by the following wavelet function 
     
       
         
           
             
               
                 
                   
                     ϕ 
                     ⁡ 
                     
                       ( 
                       x 
                       ) 
                     
                   
                   = 
                   
                     
                       cos 
                       ⁡ 
                       
                         ( 
                         
                           1.75 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           x 
                         
                         ) 
                       
                     
                     ⁢ 
                     
                       exp 
                       ⁡ 
                       
                         ( 
                         
                           - 
                           
                             
                               x 
                               2 
                             
                             2 
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   26 
                   ) 
                 
               
             
           
         
       
     
     Analogous to Example 1, the same learning parameters for LP-SVR and QP-SVR are used: C (bound)—100; ε (tolerance)—0.03. Since SVMs cannot optimize the parameters of kernels, it is difficult to determine l parameters d i , i=1,2,. . . ,l in equation (23). For simplicity, let d i =d such that the number of parameters becomes 1. In our experiment, the parameter d is set to 3 by the widely-used technique of cross validation. 
     The identification results by LP-SVR are illustrated in  FIGS. 5-6 , where the RMS error on the validation data is 0.1667 and the ratio of SVs to the training vectors is 3.1%.  FIGS. 7-8  visualize the identification results by QP-SVR, where the RMS error on the validation set is 0.1016 and the ratio of SVs to the training vectors is 33.7%. 
     The experimental results obtained in Example 1 and 2 are summarized in Table 1. 
     
       
         
           
               
             
               
                 TABLE 1 
               
             
            
               
                   
               
               
                 Comparison On Model Sparsity And RMS Error Of Different SVRs 
               
            
           
           
               
               
               
               
            
               
                   
                   
                   
                 RMS error on 
               
               
                   
                 SVR algorithm 
                 SV ratio 
                 validation set 
               
               
                   
                   
               
            
           
           
               
               
               
               
            
               
                   
                 LP-SVR with RBF kernel 
                 6.5% 
                 0.1959 
               
               
                   
                 LP-SVR with wavelet kernel 
                 3.1% 
                 0.1667 
               
               
                   
                 QP-SVR with RBF kernel 
                 39.7% 
                 0.1597 
               
               
                   
                 QP-SVR with wavelet kernel 
                 35.0% 
                 0.1076 
               
               
                   
                   
               
            
           
         
       
     
     The data in Table 1 shows that a much sparser model can be generated from LP-SVR than that from QP-SVR, with comparable prediction accuracy, and a much better sparsity can be obtained using the wavelet kernel. Particularly, the LP-SVR is around 25 times faster than QP-SVR for training, which means that the computing resources required by QP-SVR may become prohibitively-expensive with the increase of the size of training set. 
     The RMS error on validation sets acquired by using other identification algorithms are reported in Gretton et al. and J. Sjöberg et al (see above for full references), including 0.579 from wavelet networks, 0.467 from one-hidden-layer sigmoid neural networks with 10 hidden nodes, and 0.280 from v-support vector regression. Higher prediction accuracy was obtained by using the soft-constrained LP-SVR and QP-SVR. 
     An on-line iterative algorithm can be developed for linear programming support vector regression with wavelet kernel, and an efficient linear programming algorithm developed. 
     Fuzzy Support Vector Regression 
     Fuzzy support vector regression with wavelet kernel can be used for identifying a time-varying nonlinear systems. In fuzzy LP-SVR, the fuzzy membership 0≦s i ≦1 is introduced to measure the importance of the corresponding point. The approach can be summarized as: 
                       minimize   ⁢           ⁢     1   2     ⁢          β        1       +     2   ⁢           ⁢   C   ⁢       ∑     i   =   1     l     ⁢       s   i     ⁢     ξ   i             ⁢     
     ⁢     subject   ⁢           ⁢   to   ⁢           ⁢     {             y   i     =         ∑     j   =   1     l     ⁢       β   j     ⁢     k   ⁡     (       x   j     ,     x   i       )           ≤     ɛ   +     ξ   i                           ∑     j   =   1     l     ⁢       β   j     ⁢     k   ⁡     (       x   j     ,     x   i       )           -     y   i       ≤     ɛ   +     ξ   i                     ξ   i     ≥   0                       (   27   )               
where β=[β 1 , β 2 , β 3 , . . . β 1 ] T , and C is a constant.
 
     For a sequential learning problem, aσ&gt;0 was chosen as the lower bound of fuzzy membership. Second, as time is the main property of this kind of problem, fuzzy membership s i  was made a function of time t i , i.e. s i =h(t i ), where t t ≦ . . . ≦t 1  is the time the data arrived in the system. The last point x 1  is made the most important point, i.e. s 1 =h(t 1 )=1, the first point x 1  is made the least important, where s 1 =h(t 1 )=σ. 
     If fuzzy membership is made a linear function of time, we can select:
 
 s   i   =h ( t   i )=at i   +b    (28)
 
     By applying the boundary conditions, we obtain: 
     
       
         
           
             
               
                 
                   
                     s 
                     i 
                   
                   = 
                   
                     
                       h 
                       ⁡ 
                       
                         ( 
                         
                           t 
                           i 
                         
                         ) 
                       
                     
                     ⁢ 
                     
                         
                     
                     = 
                     
                       
                         
                           
                             1 
                             - 
                             σ 
                           
                           
                             
                               t 
                               l 
                             
                             - 
                             
                               t 
                               1 
                             
                           
                         
                         ⁢ 
                         
                           t 
                           i 
                         
                       
                       + 
                       
                         
                           
                             
                               t 
                               l 
                             
                             ⁢ 
                             σ 
                           
                           - 
                           
                             t 
                             1 
                           
                         
                         
                           
                             t 
                             l 
                           
                           - 
                           
                             t 
                             1 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   29 
                   ) 
                 
               
             
           
         
       
     
     Alternatively, fuzzy membership is made a quadric function of the time, then we can select:
 
 s   i   =h ( t   i )= a ( t   i   −b ) 2   +c    (30)
 
     By applying the boundary conditions, we obtain: 
     
       
         
           
             
               
                 
                   
                     s 
                     i 
                   
                   = 
                   
                     
                       h 
                       ⁡ 
                       
                         ( 
                         
                           t 
                           i 
                         
                         ) 
                       
                     
                     = 
                     
                       
                         
                           ( 
                           
                             1 
                             - 
                             σ 
                           
                           ) 
                         
                         ⁢ 
                         
                           ( 
                           
                             
                               
                                 t 
                                 i 
                               
                               - 
                               
                                 t 
                                 1 
                               
                             
                             
                               
                                 t 
                                 l 
                               
                               - 
                               
                                 t 
                                 1 
                               
                             
                           
                           ) 
                         
                       
                       + 
                       σ 
                     
                   
                 
               
               
                 
                   ( 
                   31 
                   ) 
                 
               
             
           
         
       
     
     A sliding window strategy can be used to develop an on-line learning algorithm for identifying a time-varying nonlinear system. A variable sliding window for SV learning may be implemented by defining a threshold for the violations (violations then corresponding to the model predictions differing from actual data by more than the threshold value), and counting the number of violations of within the sliding window(s). A variable-length sliding window on-line identification scheme can be adapted to cope with abrupt changes in the system to be identified. If the number of violations is below the defined threshold defined, fast updating of the linear weights of wavelet networks is performed, for example using a least squares algorithm. However, if the number of violations exceeds the defined threshold, SV learning is implemented to refresh the structure of wavelet networks using a linear programming support vector regression algorithm. This approach allows on-line training, and allows the model used to track both sudden and long-term changes in the modeled system. 
     In an example approach, sliding window blockwise SV learning is combined with a linear least square (LS) algorithm (or other fast algorithm). According to the available computational capacity and the characteristic of the identified system, a fixed sliding window size for SV learning can be defined. SV learning with the data within the sliding window is triggered by a change detection scheme. The linear LS algorithm could be used for adjusting the linear weights when the SV learning is not triggered, i.e. the model deviations from the actual system are below a triggering threshold for learning. 
     Change detection may be based on the output prediction error: 
     
       
         
           
             
               
                 
                   
                     ⅇ 
                     i 
                   
                   = 
                   
                     
                       y 
                       i 
                     
                     - 
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         l 
                       
                       ⁢ 
                       
                         
                           β 
                           j 
                         
                         ⁢ 
                         
                           k 
                           ⁡ 
                           
                             ( 
                             
                               
                                 x 
                                 j 
                               
                               , 
                               
                                 x 
                                 i 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   32 
                   ) 
                 
               
             
           
         
       
     
     Occurrence of a system change can be declared if the following accumulated detection index: 
                       d   e     ⁡     (   k   )       =       1     M   e       ⁢       ∑     j   =     k   -     M   e     +   1       k     ⁢       e   l   T     ⁢     e   l                   (   33   )               
exceeds a pre-set threshold ρe: i.e. d e ≧ρ e .
 
     Hence, a fuzzy linear programming support vector regression (LP-SVR) with wavelet kernel is well suited to handle a time-varying nonlinear systems, and provides a fast and reliable system change the system detection and modeling scheme. 
     Using a multi-dimensional wavelet kernel may require computation of the tensor product: 
     
       
         
           
             
               
                 
                   
                     
                       Φ 
                       j 
                     
                     ⁡ 
                     
                       ( 
                       x 
                       ) 
                     
                   
                   = 
                   
                     
                       ∏ 
                       
                         k 
                         = 
                         1 
                       
                       n 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       φ 
                       ( 
                       
                         
                           
                             x 
                             k 
                           
                           - 
                           
                             m 
                             jk 
                           
                         
                         
                           d 
                           jk 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   34 
                   ) 
                 
               
             
           
         
       
     
     In other examples, a multi-dimensional wavelet kernel without tensor product can be used to improve computational efficiency. 
     Applications 
     Applications of examples of the present invention include modeling of dynamic systems, such as process control, including chemical, physical, and/or mechanical process control. Applications further include image analysis (such as image cataloging, object recognition such as face detection, and the like), voice recognition, and computationally demanding applications. A particular application is engine management and control applications, described further below. 
     Engine Control 
     Vehicles, such as automobiles, often include an electronic control unit (ECU) for engine control. The ECU receives engine data from various engine sensors, and uses the engine data and an engine model to modify engine operational parameters. Engine sensors providing engine data may include a throttle position sensor, oxygen sensor, manifold absolute pressure (MAP) sensor, air temperature sensor, coolant temperature sensor, and the like. Engine parameters controlled by an ECU may include fuel injection quantities, ignition timing, variable valve timing (VVT), turbocharger boost (if applicable), and the like. In a typical example, the throttle position sensor provides a throttle position signal as part of the engine data received by the ECU, which in response increases the fuel injection quantities. 
     Engine operation is enhanced by using an ECU that uses an accurate model of engine operation. However, the number of model parameters may become large, so that the computing power necessary to set up (train) and run the model becomes substantial, particularly using conventional approaches. Engine operation is a dynamic process and response to transient data may need to be substantially immediate, such as less than a second. Operator inputs (such as throttle position) are received by the ECU, and result in control of engine operation parameters (such as fuel injection quantity). There may also be adverse safety consequences if an engine model used does not readily allow the ECU to detect or respond quickly to transient data. Hence, dynamic modeling methods using an SVM with wavelet kernel described herein are particularly useful for engine control applications. 
       FIG. 9  is a simplified schematic of an engine control system, comprising sources of engine data such as operator input sensors  10  (such as a throttle position sensor), engine operation sensors  12  (such as a coolant temperature sensor, fuel flow sensor, rpm sensor, and the like), and optional other data sources  14 , such as signals from an automatic cruise control system. Engine data is provided to the ECU  16 , and the ECU uses an engine model to provide control outputs to engine parameter controls  18  (such as fuel flow or quantity controls). The engine model is a SVM ( 20 ) within the ECU (or otherwise in communication with the ECU), having a wavelet kernel as described more fully above. An error calculation module ( 24 ) monitors the accuracy of the model, for example by comparing predicted engine behavior and actual response to engine parameter controls. 
     Data from a sliding time window may be stored in a memory within the ECU and used for comparison and retraining of the SVM. A linear programming module ( 22 ) is used to reprogram the SVM if errors are above a threshold value, for example optimizing vector parameters using stored data. Possible ECU functional modules are represented by dashed boxes in the figure, and may correspond to algorithms executable by a processor within the ECU. The ECU may further include a clock and allow rolling time gated storage of data in memory. Support vector parameters may be stored in non-volatile memory. 
     For a new vehicle engine, SVM model parameters may be determined from a previous engine testing and model training process using another of the same type of engine. However, these training parameter values may become more inaccurate as the engine wears. Hence, the dynamic model can be adjusted based on a comparison of model predictions and actual observed engine performance. This allows individual optimization of engine performance. The optimized model parameters may further be diagnostic of engine condition. 
     A support vector machine can be used to model any dynamic process, and is particularly well suited to engine management applications. In examples of the present invention, dynamic system modeling using a support vector machine with a wavelet kernel uses a computationally simpler linear programming approach. 
     Hence, embodiments of the present invention allow improved engine behavior modeling as a function of adjustable engine parameters. An electronic control unit for an engine can use an SVM-based engine model to improve engine performance. In some examples, improved engine modeling allows engine performance to be optimized in relation to one or more goals, such as exhaust emission reduction, fuel economy enhancement, acceleration, or other engine performance metric. 
     Conventional approaches become computationally challenging as the number of engine control inputs increase. Further, conventional modeling approaches, such as the use of look-up tables, are not adaptive, often use steady state operation to generate the model, and hence often do not properly describe transient behavior (non-steady state behavior) of the engine. In contrast, engine modeling according to embodiments of the present invention provides accurate representation of transient engine performance data. Further, training of an engine model according to an embodiment of the present invention may be much simpler, as there is no need to achieve numerous steady state operating conditions of the engine. Data-rich transient operation may be used for training purposes. Further, on-line training is possible, during normal operation of the vehicle. 
     Embodiments of the present invention may be used for other types of modeling, in particular modeling of any type of complex multi-input, multi-output nonlinear systems. Sparsity of the generated model facilitates real-time computation and on-line model adaptability. Embodiments of the present invention include modeling of transient responses in complex systems. 
     The invention is not restricted to the illustrative examples described above. Examples are not intended as limitations on the scope of the invention. Methods, apparatus, and the like described herein are exemplary and not intended as limitations on the scope of the invention. Changes therein and other uses will occur to those skilled in the art. The scope of the invention is defined by the scope of the claims.