Abstract:
The method for Kalman filter state estimation in bilinear systems provides for state estimation in dynamic systems, and is a bilinear extension of the Kalman filter and the Kalman smoother. The method for Kalman filter state estimation in bilinear systems introduces a nonlinear state equation coupled with a linear measurements equation. The specific nonlinearity is of the bilinear form, depending upon the system dynamics.

Description:
BACKGROUND OF THE INVENTION 
       [0001]    1. Field of the Invention 
         [0002]    The present invention relates to state estimation in dynamic systems, and particularly to a bilinear extension of the Kalman filter and the Kalman smoother. 
         [0003]    2. Description of the Related Art 
         [0004]    The Kalman filter is an estimator for the linear-quadratic problem, which is the problem of estimating the instantaneous state of a linear dynamic system perturbed by white noise using measurements that are linearly related to the state, but are also corrupted by white noise. The Kalman filter produces values that tend to be closer to the true values of the measurements and their associated calculated values by predicting an estimate of uncertainty of the predicted value via a weighted average of the predicted and measured values. The Kalman filter is also an algorithm for efficient performance of the exact inference in a linear state space model that has some statistical properties. The resulting estimator is statistically optimal with respect to some quadratic function of the estimation error. 
         [0005]    Viewed mathematically, the Kalman filter is a set of equations that provides an efficient recursive solution of the least square method. It provides estimates of the past, present, and future states, and it can do so when the precise nature of the system model is unknown. 
         [0006]    The Kalman estimator may also be adapted as a smoother. A smoother estimates the state of a system at time k, using measurements made before and after time k. The accuracy of a smoother is generally superior to that of a filter, since it uses more measurements for its estimate. The Kalman smoother can be derived from the Kalman filter model. 
         [0007]    The general derivation methodology for the Kalman smoother uses the Kalman filter for measurements up to (each) time k that the state is to be estimated, combined with another algorithm also derived from the Kalman filter for the measurements beyond that time. This second algorithm for the smoother can be derived by running the Kalman filter backward from the last measurement to just past k, and then optimally combining the two independent estimates (forward and backward) of the state at time k based on the two independent sets of measurements. 
         [0008]    In the following, the Kalman filter for a linear-discrete Gaussian state space model is considered. Initially, the nonlinear Kalman filter is examined. A discrete-time linear state space model within a dynamical system may be represented as: 
         [0000]        x   k+1   =Ax   k   +w   k ,  (1)
 
         [0000]      with measurements 
         [0000]        y   k   =Cx   k   +v   k ,  (2)
 
         [0000]    where x k ε           n  is the system state vector at time k (i.e., x k  represents an n-dimensional real vector), A ε           n×n  is the transition matrix (an n×n matrix), y k ε           p  is the corresponding measurement vector at time k, C ε           p×n  is the observation (or “measurement”) matrix (a p×1 matrix), w k ε           n  is the dynamical (or system) noise at time k, and v k ε           p  is the observation (or measurement) noise at time k. 
         [0009]    w k  and v k  are both uncorrelated, white and Gaussian, with zero mean and covariance Q and R, respectively. In other words, w k ˜N(0, Q); v k ˜N(0, R); E(w k w l   T )=Q for k=1 and E(w k  w l   T )=0 for k≠1; E(v k v l   T )=R for k=1 and E(v k v l   T )=0 for k≠1; and E(w i v j   T )=0 for ∀i, j, where E represents the expectation. 
         [0010]    From the above, for the Kalman filter, the predicted state estimate (a priori) is given by equation (1) as x k+1 =Ax k +w k , and the predicted (a priori) estimate covariance would be given by P k   f =AP f   f A T +Q, wherein the superscript “f” represents the forecast (e.g., x k   f  represents the forecast state of x k ). The updated (a posteriori) state estimate is given by x k   a =x k   f +K k (y k −Cx k   f ), where the superscript “a” represents the estimation (e.g., x k   a  is the estimation of x k ) and K k  is the optimal Kalman gain, given by K k =P k   f C T (CP k   f C T +R) −1 . The updated (a posteriori) estimate covariance is given by P k   a =P k   f −K k CP k   f . 
         [0011]    With regard to the Kalman smoother, x k  is obtained based on data samples y 1 , . . . , y t , where k≦t (denoted as x k   t ). Such estimators are generally referred to as “smoothers” since a time plot of the sequence {x k   t : k=1, . . . , t} is typically smoother than the forecasts {x k   k-1 : k=1, . . . , t} or the filter {x k   k : k=1, . . . , t}. Thus, smoothing implies that each estimated value is a function of the past, present and future, whereas the filter estimator depends on the past and the present. The forecast depends on only the past. For the Kalman smoother, the updated estimate covariance is given by P k   t =P k   a +J k └P k+1   t −P k+1   f ┘J k   T , where J k =P k   a A T (P k+1   f ) −1 . The lag-one covariance smoother is a type of smoother problem, which is a set of recursions for obtaining P k+1,k   t . For the lag-one covariance smoother, P k,k−1   t =P k   a J k−1   T +J k └P k+1,k   t −AP k   a ┘J k−1   T . 
         [0012]    The Kalman filter is an algorithm, commonly used since the 1960s, for improving vehicle navigation (among other applications, although aerospace is typical), that yields an optimized estimate of the system&#39;s state (e.g., position and velocity). The algorithm works recursively in real time on streams of noisy input observation data (typically, sensor measurements) and filters out errors using a least-squares curve-fit optimized with a mathematical prediction of the future state generated through a modeling of the system&#39;s physical characteristics. The model estimate is compared to the observation and this difference is scaled by the Kalman gain, which is then fed back as an input into the model for the purpose of improving subsequent predictions. The gain can be “tuned” for improved performance. With a high gain, the filter follows the observations more closely. With a low gain, the filter follows the model predictions more closely. This method produces estimates that tend to be closer to the true unknown values than those that would be based on a single measurement alone or the model predictions alone. 
         [0013]    The Kalman filter has many applications in technology, and is an essential part of space and military technology development. A very simple example (and perhaps the most commonly used type of Kalman filter) is the phase-locked loop, which is now ubiquitous in FM radios and most electronic communications equipment. Another common application is with guidance, navigation and control of vehicles, like aircraft or spacecraft. Sensors are used to make measurements of the vehicle&#39;s state (its position and velocity at the time of the measurement), but such measurements are intermittent and have significant stretches of time between measurements. Also, the measurements are corrupted with a certain amount of error, including noise. 
         [0014]    The Kalman Filter algorithm is an optimized method for determining the best estimation of the vehicle&#39;s state. The basic concept is similar to simple mathematical curve fitting of data points using a least-squares approximation (where the deviation is squared so that negative errors will not cancel out positive ones) and enables predictions of the state into future time steps. The most basic concepts of the filter are related to interpolation and extrapolation, where data estimates are filled in between given points and the latter involves data estimates being extended beyond the given set (as with future estimates). 
         [0015]    In each time step, the Kalman filter produces estimates of the true unknown values, along with their uncertainties. Once the outcome of the next measurement is observed, these estimates are updated using a weighted average, with more weight being given to estimates with lower uncertainty. From a theoretical standpoint, the main assumption of the Kalman filter is that the underlying system is a linear dynamical system and that all error terms and measurements have a Gaussian distribution (often a multivariate Gaussian distribution). Extensions and generalizations to the method have also been developed. The underlying model is a Bayesian model similar to a hidden Markov model, but where the state space of the latent variables is continuous and where all latent and observed variables have Gaussian distributions. 
         [0016]    The Kalman filter may be viewed as a type of “observer,” that helps to estimate the state variables of a dynamical system. For example, in general, individual state variables of a dynamic system cannot be determined exactly by direct measurements. Rather, measurements associated with the state variables, e.g., sensor data in a control system, may be a function of the state variables. In addition, such measurements are generally subject to measurement error, which, in one embodiment, may be modeled as random noise. Further, the system itself may also be subjected to random or other disturbances. The Kalman filter estimates the state variables based on the noisy measurements. 
         [0017]    In  FIG. 8 , a simplified system  200  is illustrated. The system  200  includes one or more sensors  202  that provide data to an electronic state observer  100 . The state observer  100  outputs a state (e.g., a real valued array of state variables) to a control system  206 , which outputs signals to one or more controls  208  based on the state output by the state observer  100 . As an example, the system  200  may be a control system of an aircraft or spacecraft. In such a system, the sensors  202  may include, for example, one or more of attitude, velocity, position, altitude, etc., while the controls  208  may include control surfaces, engine controls, etc. The state observer  100  is the Kalman filter, having a corresponding covariance matrix for which positive definiteness is to be maintained in order to maintain filter stability. In this example, the control system  206  receives aircraft state data from the Kalman filter  100 , which uses data from the sensors  202  to update its estimate of one or more state variables associated with the aircraft. The control system  206  generates control signals to the controls  208  of the aircraft to achieve a selected control operation. 
         [0018]      FIG. 9  illustrates the Kalman filter algorithm, using the “observer” view noted above. At step  10 , the method begins with the system  200  generating an observer, e.g., a Kalman filter, for estimation of at least one state variable (contained within the state vector) associated with the system. The Kalman filter  100  provides periodically updated state variables given a sequence of observations (i.e., sensor data) about those state variables. The data from the sensors  202  generally includes some error or noise, such as measurement errors that can depend on the sensor, the operational environment, the nature of the measured value, etc. 
         [0019]    The Kalman filter helps to remove the effects of the errors to estimate the state variables of the system. The Kalman filter is a recursive estimator or state observer. Hence, the estimated state from the previous step and a new sensor measurement are used to compute the estimate for the current state. The Kalman filter is used in two steps, i.e., a prediction step and an update step. The prediction step uses the state estimate from the previous time step to produce an estimate of the state at the current time step. In the update step, sensor measurement information at the current time step is used to refine this prediction to arrive at a new state estimate, again for the current time step. The state observer  100  is represented by the Kalman filter defined above. 
         [0020]    In step  12 , the state observer  100  receives the data from the sensors  202 . Next, at step  14 , the observer calculates the initial estimates of the state vector x k−1  and the covariance P k−1 . At step  16 , the state vector and covariance are projected ahead by one time increment, as x k+1 =Ax k +w k  and P k   f =AP f   f A T +Q, respectively. At step  18 , the Kalman gain is computed as K k =P k   f C T (CP k   f C T +R) −1 . At step  20 , a new set of measurements y k  are taken from sensors  202  to further update the state vector and covariance as x k   a =x k   f +K k (y k −Cx k   f ) and P k   a =P k   f −K k CP k   f , respectively. These updated estimates are then fed back to step  16  for incremental projection by the next time increment. 
         [0021]    With regard to identification of the linear Gaussian state-space model, “system identification” is a general term to describe the mathematical tools and algorithms that build dynamical models from measurement data. System identification plays an important role in uncertain dynamical systems. The Expectation-Maximization (EM) algorithm is a broadly applicable approach to the iterative computation of the maximum likelihood (ML) estimates, which are useful in a variety of incomplete data problems. On each iteration of the EM algorithm, there are two steps, the Expectation step (E-step) and the Maximization step (M-step). The EM algorithm is a technique that can be use to estimate the parameters after filling in the initial values for the missing data. The latter are then updated by their predicted values using these initial parameter estimates. The parameters are then re-estimated, and so on, proceeding iteratively until convergence. This technique is generally considered to be intuitive and natural. 
         [0022]    The Kalman filter and the Kalman smoother are the basic tools for calculating the expectation in the E-step. The discrete-time linear state space model within a dynamical system is given above in equations (1) and (2). In this model, the initial state vector x 0  is assumed to be a Gaussian random vector with a mean μ and a covariance V, where x 0 ˜N(μ, V). 
         [0023]    The maximum likelihood estimator for the parameters in this model is given by application of the EM algorithm. Using the Kalman filter and the Kalman smoother, x k   t =E t (x k )=E(x k |{y} 1   t ) (where E represents the expectation); P k   t =Cov(x k |{y} 1   t ), (where Cov represents the covariance); and P k,k−1   t =Cov(x k ,x k−1 |{y} 1   t ). For the M-step, the parameters of the system (A, C, Q, R and μ) are collected in a vector α={A, C, Q, R, μ}. To estimate these parameters, the corresponding partial derivatives of the expected log-likelihood function are taken, and then set to zero. Thus, α i  is updated as α i+1 . In other words, α i+1 ={A(i+1), C(i+1), Q(i+1), R(i+1), μ(i+1)} is obtained by maximizing the log likelihood function with respect to each parameter. 
         [0024]    The value of μ is simply given by μ=x 0   t =E t (x 0 ). A is expressed as A(i+1)=γβ −1 , where 
         [0000]    
       
         
           
             γ 
             = 
             
               
                 
                   ∑ 
                   
                     k 
                     = 
                     0 
                   
                   
                     t 
                     - 
                     1 
                   
                 
                  
                 
                   
                     ( 
                     
                       
                         P 
                         
                           
                             k 
                             + 
                             1 
                           
                           , 
                           k 
                         
                         t 
                       
                       + 
                       
                         
                           
                             x 
                             
                               k 
                               + 
                               1 
                             
                             t 
                           
                            
                           
                             ( 
                             
                               x 
                               k 
                               t 
                             
                             ) 
                           
                         
                         T 
                       
                     
                     ) 
                   
                    
                   
                       
                   
                    
                   and 
                    
                   
                       
                   
                    
                   β 
                 
               
               = 
               
                 
                   ∑ 
                   
                     k 
                     = 
                     0 
                   
                   
                     t 
                     - 
                     1 
                   
                 
                  
                 
                   
                     ( 
                     
                       
                         P 
                         k 
                         t 
                       
                       + 
                       
                         
                           
                             x 
                             k 
                             t 
                           
                            
                           
                             ( 
                             
                               x 
                               k 
                               t 
                             
                             ) 
                           
                         
                         T 
                       
                     
                     ) 
                   
                   . 
                 
               
             
           
         
       
     
         [0000]    C is given by C(i+1)=εβ −1 , where 
         [0000]    
       
         
           
             ɛ 
             = 
             
               
                 ∑ 
                 
                   k 
                   = 
                   0 
                 
                 
                   t 
                   - 
                   1 
                 
               
                
               
                 
                   
                     E 
                     t 
                   
                    
                   
                     ( 
                     
                       
                         y 
                         k 
                       
                        
                       
                         x 
                         k 
                         T 
                       
                     
                     ) 
                   
                 
                 . 
               
             
           
         
       
     
         [0000]    Q is given by Q(i+1)=t −1 (σ−γβ −1 γ T ), where 
         [0000]    
       
         
           
             σ 
             = 
             
               
                 ∑ 
                 
                   k 
                   = 
                   0 
                 
                 
                   t 
                   - 
                   1 
                 
               
                
               
                 
                   ( 
                   
                     
                       P 
                       k 
                       t 
                     
                     + 
                     
                       
                         x 
                         k 
                         t 
                       
                        
                       
                         x 
                         k 
                         
                           t 
                           T 
                         
                       
                     
                   
                   ) 
                 
                 . 
               
             
           
         
       
     
         [0000]    Lastly, R is given by R(i+1)=t −1 (ξ−εβ− 1 ε T ), where 
         [0000]    
       
         
           
             ξ 
             = 
             
               
                 ∑ 
                 
                   k 
                   = 
                   0 
                 
                 
                   t 
                   - 
                   1 
                 
               
                
               
                 
                   
                     E 
                     t 
                   
                    
                   
                     ( 
                     
                       
                         y 
                         k 
                       
                        
                       
                         y 
                         k 
                         T 
                       
                     
                     ) 
                   
                 
                 . 
               
             
           
         
       
     
         [0025]    In order to identify the parameters of the state-space model using the EM algorithm to obtain the maximum likelihood estimation vector {circumflex over (α)}, first the initial values of A(0), C(0), Q(0), R(0) and μ(0) are selected, and then x k   t , P k   t  and P k,k−1   t  are computed using the Kalman smoother. The conditional expectation of the log likelihood function is then calculated, and then {A(i+1),C(i+1),Q(i+1),R(i+1), μ(i+1)} are computed to find the next iterative estimated parameters that maximize the conditional expectation of the log likelihood. These new parameters are then inserted back into the state-space model, and computed using the Kalman smoother. These steps are repeated until the log likelihood converges. 
         [0026]    When the dynamical system and measurements system of the state-space model are linear, the Kalman filter is the optimal estimator. However, in most “real world” processes of interest, the dynamical and/or the measurements systems are nonlinear. Thus, a suitable extension of the Kalman filter is needed. Both the Lorenz-96 model and the Lotka-Volterra model, which are used for atmospheric dynamics, biology, chemistry, and control systems, are examples of bilinear models. Generally, bilinear systems are nonlinear systems that are of great interest, as they represent a wide variety of important physical processes. Bilinear systems can also be used in approximation or alternate representations for a range of other nonlinear systems. Further, bilinear models may be used to model nonlinear processes in signal, image and communication systems. Bilinear systems may be represented as state-space models. However, bilinear models do not work with the traditional Kalman filter. It would, therefore, be desirable to extend the Kalman filter and Kalman smoother to be able to handle such bilinear systems. 
         [0027]    Thus, a method for Kalman filter state estimation in bilinear systems solving the aforementioned problems is desired. 
       SUMMARY OF THE INVENTION 
       [0028]    The present invention relates to state estimation in dynamic systems, and particularly to a bilinear extension of the Kalman filter and the Kalman smoother. The method for Kalman filter state estimation in bilinear systems introduces a nonlinear state equation coupled with a linear measurements equation. The specific nonlinearity is of the bilinear form, depending upon the system dynamics. 
         [0029]    The method for Kalman filter state estimation in bilinear systems includes the steps of: (a) generating an observer in an electronic device for estimating a state vector x k , associated with a bilinear system, where the observer defines a covariance matrix P k  statistically relating the state vector x k , where k represents an integer time increment, and the bilinear system is defined by x k+1 =Ax k +B(x k           x k )+w k , where A is a transition matrix defined by the bilinear system, B is a measurement matrix defined by parameters of the bilinear system, and w k  is system noise at time k; (b) receiving sensor data y k  at the electronic device from a sensor, the sensor data being associated with the state vector x k  as y k =Cx k +v k , where C is a measurement matrix defined by the parameters of the bilinear system and v k  is measurement noise at time k; (e) generating initial estimates of the state vector and the covariance matrix from the sensor data and storing the initial estimates of the state vector and the covariance matrix in non-transitory computer readable memory of the electronic device; (d) calculating a projected state vector x k+1   k  representing the state vector projected ahead by one time increment as x k+1   l =Ax k   k +Bz k   k  and storing the projected state vector x k+1   k  in the non-transitory computer readable memory of the electronic device, where z k =x k           x k  and z k   k =E k  (z k ); (e) calculating a projected covariance matrix P k+   1  representing the covariance matrix projected ahead by one time increment as P k+1   k =AP k   k A T +A{umlaut over (P)} k   k B T +B({umlaut over (P)} k   k ) T A T +B{dot over (P)} k   k B T +Q, where Q represents a covariance of the system noise w k , and storing the projected covariance matrix P k+1   k  in the non-transitory computer readable memory of the electronic device; (f) calculating a Kalman gain K k+1  as K k+1 =P k+1   k C T (CP k+1   k C T +R) −1 , where R is a covariance of the measurement noise v k , and storing the Kalman gain K k+1  in the non-transitory computer readable memory of the electronic device; (g) receiving a new set of measurement data y k  from the sensor and updating the state vector and the covariance matrix based upon the new set of measurement data as x k+1   k+1 =K k+1   k +K k+1 (y k −Cx +1   k ) and P +1   k+ =(I−K k+1 C)P k+1   k , respectively, wherein x k+1   k+1  represents the updated state vector for time (k+1) and P k+1   k+1  represents the updated covariance matrix for time (k+1), where I represents the identity matrix, and storing the updated state vector and the updated covariance matrix in the non-transitory computer readable memory of the electronic device; and (h) returning to steps (d) and (e) with the updated state vector and the updated covariance matrix, respectively. 
         [0030]    These and other features of the present invention will become readily apparent upon further review of the following specification and drawings. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0031]      FIG. 1  is a block diagram illustrating an overview of a system for performing a method for Kalman filter state estimation in bilinear systems according to the present invention. 
           [0032]      FIG. 2  is a graph illustrating calculated error between actual states and states estimated by the method for Kalman filter state estimation in bilinear systems according to the present invention. 
           [0033]      FIG. 3  is a graph illustrating calculated error between actual states and states estimated by a Kalman smoother extension of the method for Kalman filter state estimation in bilinear systems according to the present invention. 
           [0034]      FIG. 4  is a graph illustrating calculated error using parameters estimated by the method for Kalman filter state estimation in bilinear systems according to the present invention produced via the EM algorithm. 
           [0035]      FIG. 5  is a graph illustrating calculated error between simulated states and states estimated by the method for Kalman filter state estimation in bilinear system according to the present inventions. 
           [0036]      FIG. 6  is a graph illustrating calculated error between simulated states and states estimated by the ensemble Kalman filter algorithm for purposes of comparison with  FIG. 5 . 
           [0037]      FIG. 8  is block diagram illustrating a system overview for implementation of a method for Kalman filter state estimation according to the prior art. 
           [0038]      FIG. 9  is a flow chart illustrating method steps of the Kalman filter algorithm according to the prior art. 
       
    
    
       [0039]    Similar reference characters denote corresponding features consistently throughout the attached drawings. 
       DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
       [0040]    The discrete-time linear state space model within a dynamical system is represented by equations (1) and (2) above. The bilinear Gaussian discrete state space model is a variant on equations (1) and (2) and is represented by: 
         [0000]        x   k+1   =Ax   k   +B ( x   k             x   k )+ w   k ,  (3)
 
         [0000]      with measurements 
         [0000]        y   k   =Cx   k   +v   k ,  (4)
 
         [0000]    where x k ε           n  is the system state vector at time k (i.e., x k  represents an n-dimensional real vector), A ⊖           n×n  is the transition matrix (an n×n matrix), y k ε           p  is the corresponding measurement vector at time k, B ε           n×[n(n+1)/2]  and C ε           p×n  are the observation (or “measurement”) matrices (i.e., the parameters of the model), w k ε           n  is the dynamical (or system) noise at time k, and v k ε           p  is the observation (or measurement) noise at time k. 
         [0041]    w k  and v k  are both uncorrelated, white and Gaussian, with zero mean and covariance Q and R, respectively. In other words, w k ˜N(0, Q); v k ˜N(0, R); E(w k w l   T )=Q for k=1 and E(w k w l   T )=0 for k≠1, E(v k v l   T )=R for k=1 and E(v k v l   T )=0 for k≠1; and E(w k v j   t )=0, where E represents the expectation and (x k           x k ) is the Kronecker product of the state x k  with itself. 
         [0042]    Due to the nonlinearity in the bilinear model, the estimation is performed as a set of differential equations. For the bilinear state space model of equations (3) and (4), the state vector and estimator are respectively given as: 
         [0000]        x   k+1   k   =Ax   k   k   +Bz   k   k   (5)
 
         [0000]        P   k+1   k   =AP   k   k   A   T   A{umlaut over (P)}   k   k   B   T   +B ( {umlaut over (P)}   k   k ) T   A   T   +B{dot over (P)}   k   k   B   T   +Q,   (6)
 
         [0000]    where z k =x k           x k  and z k   k =E k (z k ), where E represents the expectation, and 
         [0000]        x   k+1   k+1   =x   k+1   k   +K   k+1 ( y   k   −Cx   k+1   k )  (7)
 
         [0000]        P   k+1   k+1 =( I−K   k+1   C ) P   k+1   k ,  (8)
 
         [0000]    where I represents the identity matrix and K k  is the optimal Kalman gain. In the above, x k   t =E{x k |{y} l   t }=E t x t ), P k   t =E{(x k −x k   t )(x k −x k   t ) T }, z k   t =E{z k |{y} l   t }=E t (z t ), {dot over (P)} k   t =E{(z k −z k   t )(z k −z k   t ) T }, and {umlaut over (P)} k   t =E t {(x k −x k   t )(z k −z k   t ) T }, where 1≦k≦t, 1≦t≦n, and {y} l   t  is the measurement sequence; i.e., {y} l   t ={y j , . . . , y t }. 
         [0043]    Applying a Taylor expansion to z k  yields 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         z 
                         k 
                       
                       = 
                       
                         
                           z 
                           k 
                           j 
                         
                         + 
                         
                           
                             
                               f 
                               ′ 
                             
                              
                             
                               ( 
                               
                                 x 
                                 k 
                                 j 
                               
                               ) 
                             
                           
                            
                           
                             ( 
                             
                               
                                 x 
                                 k 
                               
                               - 
                               
                                 x 
                                 k 
                                 j 
                               
                             
                             ) 
                           
                         
                         + 
                         
                           
                             1 
                             2 
                           
                            
                           
                             H 
                              
                             
                               ( 
                               
                                 
                                   x 
                                   k 
                                 
                                 , 
                                 
                                   x 
                                   k 
                                   j 
                                 
                               
                               ) 
                             
                           
                            
                           
                             ( 
                             
                               
                                 x 
                                 k 
                               
                               - 
                               
                                 x 
                                 k 
                                 j 
                               
                             
                             ) 
                           
                         
                       
                     
                     , 
                     
                       
 
                     
                      
                     where 
                   
                    
                   
                     
 
                   
                    
                   
                     
                       
                         z 
                         k 
                         j 
                       
                       = 
                       
                         
                           E 
                           j 
                         
                          
                         
                           ( 
                           
                             z 
                             k 
                           
                           ) 
                         
                       
                     
                     , 
                     
                       
                         
                           f 
                           ′ 
                         
                          
                         
                           ( 
                           x 
                           ) 
                         
                       
                       = 
                       
                         
                           [ 
                           
                             
                               
                                 ∂ 
                                 
                                   x 
                                   i 
                                 
                               
                                
                               
                                 x 
                                 j 
                               
                             
                             
                               ∂ 
                               
                                 x 
                                 l 
                               
                             
                           
                           ] 
                         
                         
                           i 
                           , 
                           j 
                           , 
                           
                             l 
                             = 
                             1 
                           
                           , 
                           2 
                           , 
                           … 
                            
                           
                               
                           
                           , 
                           m 
                         
                       
                     
                     , 
                     and 
                   
                    
                   
                     
 
                   
                    
                   
                     
                       
                         H 
                          
                         
                           ( 
                           
                             
                               x 
                               
                                 k 
                                 + 
                                 1 
                               
                             
                             , 
                             
                               x 
                               
                                 k 
                                 + 
                                 1 
                               
                               j 
                             
                           
                           ) 
                         
                       
                       = 
                       
                         ( 
                         
                           
                             
                               
                                 
                                   
                                     ( 
                                     
                                       
                                         x 
                                         
                                           k 
                                           + 
                                           1 
                                         
                                       
                                       - 
                                       
                                         x 
                                         
                                           k 
                                           + 
                                           1 
                                         
                                         j 
                                       
                                     
                                     ) 
                                   
                                   T 
                                 
                                  
                                 
                                   D 
                                   1 
                                 
                               
                             
                           
                           
                             
                               
                                 
                                   
                                     ( 
                                     
                                       
                                         x 
                                         
                                           k 
                                           + 
                                           1 
                                         
                                       
                                       - 
                                       
                                         x 
                                         
                                           k 
                                           + 
                                           1 
                                         
                                         j 
                                       
                                     
                                     ) 
                                   
                                   T 
                                 
                                  
                                 
                                   D 
                                   2 
                                 
                               
                             
                           
                           
                             
                               ⋮ 
                             
                           
                           
                             
                               
                                 
                                   
                                     ( 
                                     
                                       
                                         x 
                                         
                                           k 
                                           + 
                                           1 
                                         
                                       
                                       - 
                                       
                                         x 
                                         
                                           k 
                                           + 
                                           1 
                                         
                                         j 
                                       
                                     
                                     ) 
                                   
                                   T 
                                 
                                  
                                 
                                   D 
                                   m 
                                 
                               
                             
                           
                         
                         ) 
                       
                     
                     , 
                     where 
                   
                    
                   
                     
 
                   
                    
                   
                     D 
                     = 
                     
                       ( 
                       
                         
                           
                             
                               D 
                               1 
                             
                           
                         
                         
                           
                             ⋮ 
                           
                         
                         
                           
                             
                               D 
                               m 
                             
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
         [0000]    is the m 2 ×m matrix of second derivatives. Having defined z k   j , f′(x) and H(x k+1 , x k+1   j ) the bilinear state space model Kalman filter can be completed with the expressions: 
         [0000]        {umlaut over (P)}   k+1   k+1   =P   k+1   k+1 ( V   k+1   k+1 ) T   (10)
 
         [0000]        P   k+1   k+1 =( V   k+1   k+1 ) {umlaut over (P)}   k+1   k+1 ,  (11)
 
         [0000]    where the Kalman gain is given by: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       K 
                       
                         k 
                         + 
                         1 
                       
                     
                     = 
                     
                       
                         P 
                         
                           k 
                           + 
                           1 
                         
                         k 
                       
                        
                       
                         
                           
                             C 
                             T 
                           
                            
                           
                             ( 
                             
                               
                                 
                                   CP 
                                   
                                     k 
                                     + 
                                     1 
                                   
                                   k 
                                 
                                  
                                 
                                   C 
                                   T 
                                 
                               
                               + 
                               R 
                             
                             ) 
                           
                         
                         
                           - 
                           1 
                         
                       
                     
                   
                   , 
                   and 
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       V 
                       
                         k 
                         + 
                         1 
                       
                       
                         k 
                         + 
                         1 
                       
                     
                     = 
                     
                       
                         
                           f 
                           ′ 
                         
                          
                         
                           ( 
                           
                             x 
                             
                               k 
                               + 
                               1 
                             
                             
                               k 
                               + 
                               1 
                             
                           
                           ) 
                         
                       
                       + 
                       
                         
                           1 
                           2 
                         
                          
                         
                           H 
                            
                           
                             ( 
                             
                               
                                 x 
                                 
                                   k 
                                   + 
                                   1 
                                 
                                 k 
                               
                               , 
                               
                                 x 
                                 
                                   k 
                                   + 
                                   1 
                                 
                                 
                                   k 
                                   + 
                                   1 
                                 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                    
                   
                     
 
                   
                    
                   
                     
                       
                         for 
                          
                         
                             
                         
                          
                         k 
                       
                       = 
                       0 
                     
                     , 
                     … 
                      
                     
                         
                     
                     , 
                     
                       t 
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
         [0044]    The Kalman smoother is an estimator which operates in the time domain where t&gt;k+1. For the interval of k=t−1, . . . , 1, the Kalman smoother is given by: 
         [0000]        x   k   t   =x   k   k   +J   k ( x   k+1   t   −x   k+1   k )  (14)
 
         [0000]        P   k   t   =P   k   k   +J   k ( P   k+1   t   −P   k+1   k ) J   k   T ,  (15)
 
         [0000]    where J k =(P k   k A T +{umlaut over (P)} k   k B T )(P k+1   k ) −1 . The bilinear lag-one covariance smoother is simply given by: 
         [0000]        P   k+1,k   t   =AP   k   k   +B ( {umlaut over (P)}   k   k ) T , and  (16)
 
         [0000]        {umlaut over (P)}   k+1,k   t   =P   k+1,k   t ( V   k   t ) T .  (17)
 
         [0045]    As noted above, the term system identification describes the mathematical tools and algorithms that build dynamical models from measured data. In the following, the bilinear state space model defined by equations (3) and (4) is identified using the expectation maximization algorithm. The expectation maximization (EM) algorithm is an iterative technique used for obtaining the maximum likelihood estimation. In the EM algorithm, first q(θ, θ(i))=E{log p(θ, X, Y)|Y} is computed, where p(θ, X, Y) is the maximum likelihood function. This is the E-step. Next, in the M-step, q(θ, θ(i)) is maximized with respect to θ, which is the parameter vector, given by θ={A, B, C, Q, R, V, μ}; X={x} 1   t ={x 1 , . . . , x t } and Y={y} 1   t ={y 1 , . . . , y t }. 
         [0046]    For the bilinear state space model of equations (3) and (4), the following assumptions are applied: x 0 ˜N(μ, V); w k ˜N(0, Q); and v k ˜N(0, R). Then, q(θ, θ(i)) is given by: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         
                           q 
                            
                           
                             ( 
                             
                               θ 
                               , 
                               
                                 θ 
                                  
                                 
                                   ( 
                                   i 
                                   ) 
                                 
                               
                             
                             ) 
                           
                         
                         = 
                           
                          
                         
                           E 
                            
                           
                             { 
                             
                               
                                 log 
                                  
                                 
                                     
                                 
                                  
                                 
                                   p 
                                    
                                   
                                     ( 
                                     
                                       θ 
                                       , 
                                       X 
                                       , 
                                       Y 
                                     
                                     ) 
                                   
                                 
                               
                               | 
                               Y 
                             
                             } 
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                          
                         
                           
                             
                               - 
                               
                                 1 
                                 2 
                               
                             
                              
                             log 
                              
                             
                                
                               V 
                                
                             
                           
                           - 
                           
                             
                               1 
                               2 
                             
                              
                             Tr 
                              
                             
                               { 
                               
                                 
                                   V 
                                   
                                     - 
                                     1 
                                   
                                 
                                  
                                 
                                   ( 
                                   
                                     Δ 
                                     - 
                                     
                                       
                                         
                                           x 
                                           ^ 
                                         
                                         0 
                                       
                                        
                                       
                                         μ 
                                         T 
                                       
                                     
                                     - 
                                     
                                       μ 
                                        
                                       
                                           
                                       
                                        
                                       
                                         
                                           x 
                                           ^ 
                                         
                                         0 
                                         T 
                                       
                                     
                                     + 
                                     
                                       μμ 
                                       T 
                                     
                                   
                                   ) 
                                 
                               
                               } 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   18 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         
                           - 
                           
                             t 
                             2 
                           
                         
                          
                         log 
                          
                         
                            
                           Q 
                            
                         
                       
                        
                       
                         
 
                       
                       - 
                       
                         
                           1 
                           2 
                         
                          
                         Tr 
                          
                         
                           { 
                           
                             
                               Q 
                               
                                 - 
                                 1 
                               
                             
                              
                             
                               ( 
                               
                                 
                                   
                                     
                                       Θ 
                                       - 
                                       
                                         Ψ 
                                          
                                         
                                             
                                         
                                          
                                         
                                           A 
                                           T 
                                         
                                       
                                       - 
                                       
                                         Π 
                                          
                                         
                                             
                                         
                                          
                                         
                                           B 
                                           T 
                                         
                                       
                                       - 
                                       
                                         A 
                                          
                                         
                                             
                                         
                                          
                                         
                                           Ψ 
                                           T 
                                         
                                       
                                       + 
                                     
                                   
                                 
                                 
                                   
                                     
                                       
                                         A 
                                          
                                         
                                             
                                         
                                          
                                         Φ 
                                          
                                         
                                             
                                         
                                          
                                         
                                           A 
                                           T 
                                         
                                       
                                       - 
                                       
                                         B 
                                          
                                         
                                             
                                         
                                          
                                         
                                           Π 
                                           T 
                                         
                                       
                                       + 
                                       
                                         B 
                                          
                                         
                                             
                                         
                                          
                                         Λ 
                                          
                                         
                                             
                                         
                                          
                                         
                                           B 
                                           T 
                                         
                                       
                                     
                                   
                                 
                               
                               ) 
                             
                           
                           } 
                         
                       
                        
                       
                         
 
                       
                       - 
                       
                         
                           t 
                           2 
                         
                          
                         log 
                          
                         
                            
                           R 
                            
                         
                       
                       - 
                       
                         
                           1 
                           2 
                         
                          
                         Tr 
                          
                         
                           { 
                           
                             
                               R 
                               
                                 - 
                                 1 
                               
                             
                              
                             
                               ( 
                               
                                 δ 
                                 - 
                                 
                                   Ω 
                                    
                                   
                                       
                                   
                                    
                                   
                                     C 
                                     T 
                                   
                                 
                                 - 
                                 
                                   C 
                                    
                                   
                                       
                                   
                                    
                                   
                                     Ω 
                                     T 
                                   
                                 
                                 + 
                                 
                                   C 
                                    
                                   
                                       
                                   
                                    
                                   Φ 
                                    
                                   
                                       
                                   
                                    
                                   
                                     C 
                                     T 
                                   
                                 
                               
                               ) 
                             
                           
                           } 
                         
                       
                       + 
                       
                         const 
                         . 
                       
                     
                     , 
                     
                       
 
                     
                      
                     
                       
                         
                           where 
                            
                           
                               
                           
                            
                           Δ 
                         
                         = 
                         
                           
                             E 
                             t 
                           
                            
                           
                             ( 
                             
                               
                                 x 
                                 0 
                               
                                
                               
                                 x 
                                 0 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                       ; 
                       
                         
                           
                             x 
                             ^ 
                           
                           0 
                         
                         = 
                         
                           
                             E 
                             t 
                           
                            
                           
                             ( 
                             
                               x 
                               0 
                             
                             ) 
                           
                         
                       
                       ; 
                       
                         Θ 
                         = 
                         
                           
                             ∑ 
                             
                               k 
                               = 
                               0 
                             
                             
                               t 
                               - 
                               1 
                             
                           
                            
                           
                             
                               E 
                               t 
                             
                              
                             
                               ( 
                               
                                 
                                   x 
                                   
                                     k 
                                     + 
                                     1 
                                   
                                 
                                  
                                 
                                   x 
                                   
                                     k 
                                     + 
                                     1 
                                   
                                   T 
                                 
                               
                               ) 
                             
                           
                         
                       
                       ; 
                     
                   
                    
                   
                     
 
                   
                    
                   
                     
                       Φ 
                       = 
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             0 
                           
                           
                             t 
                             - 
                             1 
                           
                         
                          
                         
                           
                             E 
                             t 
                           
                            
                           
                             ( 
                             
                               
                                 x 
                                 k 
                               
                                
                               
                                 x 
                                 k 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                     
                     ; 
                     
                       Ψ 
                       = 
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             0 
                           
                           
                             t 
                             - 
                             1 
                           
                         
                          
                         
                           
                             E 
                             t 
                           
                            
                           
                             ( 
                             
                               
                                 x 
                                 
                                   k 
                                   + 
                                   1 
                                 
                               
                                
                               
                                 x 
                                 k 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                     
                     ; 
                     
                       Π 
                       = 
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             0 
                           
                           
                             t 
                             - 
                             1 
                           
                         
                          
                         
                           
                             E 
                             t 
                           
                            
                           
                             ( 
                             
                               
                                 x 
                                 
                                   k 
                                   + 
                                   1 
                                 
                               
                                
                               
                                 z 
                                 k 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                     
                     ; 
                   
                    
                   
                     
 
                   
                    
                   
                     
                       Γ 
                       = 
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             0 
                           
                           
                             t 
                             - 
                             1 
                           
                         
                          
                         
                           
                             E 
                             t 
                           
                            
                           
                             ( 
                             
                               
                                 x 
                                 k 
                               
                                
                               
                                 z 
                                 k 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                     
                     ; 
                     
                       Λ 
                       = 
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             0 
                           
                           
                             t 
                             - 
                             1 
                           
                         
                          
                         
                           
                             E 
                             t 
                           
                            
                           
                             ( 
                             
                               
                                 z 
                                 k 
                               
                                
                               
                                 z 
                                 k 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                     
                     ; 
                     
                       Ω 
                       = 
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             0 
                           
                           
                             t 
                             - 
                             1 
                           
                         
                          
                         
                           
                             E 
                             t 
                           
                            
                           
                             ( 
                             
                               
                                 y 
                                 k 
                               
                                
                               
                                 x 
                                 k 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                     
                     ; 
                   
                    
                   
                     
 
                   
                    
                   
                     
                       δ 
                       = 
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             0 
                           
                           
                             t 
                             - 
                             1 
                           
                         
                          
                         
                           
                             E 
                             t 
                           
                            
                           
                               
                           
                            
                           
                             ( 
                             
                               
                                 y 
                                 k 
                               
                                
                               
                                 y 
                                 k 
                                 T 
                               
                             
                             ) 
                           
                         
                       
                     
                     ; 
                     
                       
                         and 
                          
                         
                             
                         
                          
                         
                           z 
                           k 
                         
                       
                       = 
                       
                         
                           x 
                           k 
                         
                         ⊗ 
                         
                           
                             x 
                             k 
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                     
                 
               
             
           
         
       
     
         [0047]    The function q(θ, θ(i)), given above in equation (18), is maximized over θ when: 
         [0000]    
       
         
           
             
               μ 
               = 
               
                 
                   x 
                   ^ 
                 
                 0 
               
             
             ; 
             
               V 
               = 
               Δ 
             
             ; 
             
               A 
               = 
               
                 Ψ 
                  
                 
                     
                 
                  
                 
                   Φ 
                   
                     - 
                     1 
                   
                 
               
             
             ; 
             
               B 
               = 
               
                 Π 
                  
                 
                     
                 
                  
                 
                   Λ 
                   
                     - 
                     1 
                   
                 
               
             
             ; 
           
         
       
       
         
           
             
               Q 
               = 
               
                 
                   1 
                   t 
                 
                  
                 
                   ( 
                   
                     Θ 
                     - 
                     
                       
                         Ψ 
                         T 
                       
                        
                       
                         Φ 
                         
                           - 
                           1 
                         
                       
                        
                       Ψ 
                     
                     - 
                     
                       
                         ΠΛ 
                         
                           - 
                           1 
                         
                       
                        
                       
                         Π 
                         T 
                       
                     
                   
                   ) 
                 
               
             
             ; 
           
         
       
       
         
           
             
               C 
               = 
               
                 Ω 
                  
                 
                     
                 
                  
                 
                   Φ 
                   
                     - 
                     1 
                   
                 
               
             
             ; 
             
               
                 and 
                  
                 
                     
                 
                  
                 R 
               
               = 
               
                 
                   t 
                   
                     - 
                     1 
                   
                 
                  
                 
                   
                     { 
                     
                       δ 
                       - 
                       
                         
                           ΩΦ 
                           
                             - 
                             1 
                           
                         
                          
                         
                           Ω 
                           T 
                         
                       
                     
                     } 
                   
                   . 
                 
               
             
           
         
       
     
         [0048]    In the following, the above bilinear Kalman model is applied to the nonlinear estimation and identification of the Lotka-Volterra nonlinear model. The Lotka-Volterra nonlinear model has wide ranging applications in various domains of life science, biology, chemistry, economics and neural networks. In order to show the efficiency and accuracy of the present bilinear Kalman filter, the bilinear Kalman filter and smoother are applied to simultaneously estimate states and parameters from noise data of a Lotka-Volterra system. 
         [0049]    The Lotka-Volterra predator-prey equations are a pair of first-order, nonlinear differential equations frequently used to describe the dynamics of biological systems in which two species interact; i.e., predator and prey. The Lotka-Volterra model has the following form: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                        
                       x 
                     
                     
                        
                       t 
                     
                   
                   = 
                   
                     
                       λ 
                        
                       
                           
                       
                        
                       x 
                     
                     - 
                     
                       β 
                        
                       
                           
                       
                        
                       xy 
                     
                   
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                          
                         y 
                       
                       
                          
                         t 
                       
                     
                     = 
                     
                       
                         β 
                          
                         
                             
                         
                          
                         xy 
                       
                       - 
                       
                         γ 
                          
                         
                             
                         
                          
                         y 
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where y represents the number of predators and x represents the number of the corresponding prey, 
         [0000]    
       
         
           
             
               
                  
                 x 
               
               
                  
                 t 
               
             
              
             
                 
             
              
             and 
              
             
                 
             
              
             
               
                  
                 y 
               
               
                  
                 t 
               
             
           
         
       
     
         [0000]    represent the growth of the two populations against time t, and λ, β and γ are parameters representing the interaction of the two species. 
         [0050]    The Lotka-Volterra competing species model is a system of ordinary differential equations of the following form: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                        
                       x 
                     
                     
                        
                       t 
                     
                   
                   = 
                   
                     
                       λ 
                        
                       
                           
                       
                        
                       
                         x 
                          
                         
                           ( 
                           
                             1 
                             - 
                             x 
                           
                           ) 
                         
                       
                     
                     - 
                     
                       β 
                        
                       
                           
                       
                        
                       xy 
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
             
               
                 
                   
                     
                        
                       y 
                     
                     
                        
                       t 
                     
                   
                   = 
                   
                     
                       γ 
                        
                       
                           
                       
                        
                       
                         y 
                          
                         
                           ( 
                           
                             1 
                             - 
                             y 
                           
                           ) 
                         
                       
                     
                     - 
                     
                       β 
                        
                       
                           
                       
                        
                       
                         xy 
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   22 
                   ) 
                 
               
             
           
         
       
     
         [0000]    Equation (21) indicates that the population of the species x grows according to a logistic law in absence of species y (i.e., when y=0). Additionally, the rate of growth of x is negatively proportional to y, representing competition between members of x and members of y. The larger the population y, the smaller the growth rate of x. Similarly, equation (22) describes the rate of growth for population y. 
         [0051]    The Lotka-Volterra competing species model can be written in the form of the present bilinear model (of equations (3) and (4)), with 
         [0000]    
       
         
           
             A 
             = 
             
               
                 
                   ( 
                   
                     
                       
                         
                           λ 
                            
                           
                               
                           
                         
                       
                       
                         0 
                       
                     
                     
                       
                         0 
                       
                       
                         γ 
                       
                     
                   
                   ) 
                 
                  
                 
                     
                 
                  
                 and 
                  
                 
                     
                 
                  
                 B 
               
               = 
               
                 
                   ( 
                   
                     
                       
                         
                           - 
                           λ 
                         
                       
                       
                         
                           - 
                           β 
                         
                       
                       
                         0 
                       
                     
                     
                       
                         0 
                       
                       
                         
                           - 
                           β 
                         
                       
                       
                         
                           - 
                           γ 
                         
                       
                     
                   
                   ) 
                 
                 . 
               
             
           
         
       
     
         [0000]    The same can be performed for the Lotka-Volterra predator-prey model of equations (19) and (20). 
         [0052]    The following simulation using the present bilinear Kalman filter and bilinear Kalman smoother was performed on the nonlinear Lotka-Volterra competition model in state-space form, which is given by: 
         [0000]    
       
         
           
             
               
                 
                   
                     x 
                     
                       k 
                       + 
                       1 
                     
                   
                   = 
                   
                     
                       Ax 
                       k 
                     
                     + 
                     
                       Bz 
                       k 
                     
                     + 
                     
                       w 
                       k 
                     
                   
                 
               
               
                 
                   ( 
                   23 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         y 
                         
                           k 
                           + 
                           1 
                         
                       
                       = 
                       
                         
                           Cx 
                           k 
                         
                         + 
                         
                           v 
                           k 
                         
                       
                     
                     , 
                     
                       
 
                     
                      
                     
                       
                         where 
                          
                         
                             
                         
                          
                         A 
                       
                       = 
                       
                         ( 
                         
                           
                             
                               1 
                             
                             
                               0.1 
                             
                           
                           
                             
                               
                                 - 
                                 0.1 
                               
                             
                             
                               1 
                             
                           
                         
                         ) 
                       
                     
                     , 
                     
                       B 
                       = 
                       
                         ( 
                         
                           
                             
                               0.01 
                             
                             
                               0.01 
                             
                             
                               0 
                             
                           
                           
                             
                               0 
                             
                             
                               
                                 - 
                                 0.01 
                               
                             
                             
                               0 
                             
                           
                         
                         ) 
                       
                     
                   
                    
                   
                     
 
                   
                    
                   
                     
                       
                         and 
                          
                         
                             
                         
                          
                         C 
                       
                       = 
                       
                         ( 
                         
                           
                             
                               1 
                             
                             
                               0 
                             
                           
                           
                             
                               0 
                             
                             
                               1 
                             
                           
                         
                         ) 
                       
                     
                     , 
                   
                 
               
               
                 
                   ( 
                   24 
                   ) 
                 
               
             
           
         
       
     
         [0000]    and the bilinear term z k =x k           x k , where x k  is the state vector. The random noise w k  and v k  are uncorrelated with w k ˜(0, W) and v k ˜(0, V), where W=0.0004I 2  and V=0.0004I 2 . The initial state is 
         [0000]    
       
         
           
             
               x 
               0 
             
             = 
             
               
                 ( 
                 
                   
                     
                       1 
                     
                   
                   
                     
                       1 
                     
                   
                 
                 ) 
               
               . 
             
           
         
       
     
         [0053]    The error for the estimated quantities is required in order to state the reliability of the results. The error is provided by the covariance matrix to compute the difference between the true states and the estimated states. In other words, if x k  is the true state, and x k   t  is the estimated state, then the estimation error can be computed as: 
         [0000]      ε=∥ x   k   −x   k   t ∥=( x   k   −x   k   t ) T ( x   k   −x   k   t ).  (25)
 
         [0000]    In the bilinear Kalman filter, the estimated state is x k+1 . The value of x k+1   k  can be obtained from equation (5). Having run the simulation with the present Kalman filter and Kalman smoother, it is found that the error between the true states and the estimated states via the bilinear Kalman filter is very small. The error in the case of the bilinear filter is shown in  FIG. 2 . 
         [0054]    In the bilinear Kalman smoother, the estimated state is x k+1   t  for t&gt;k+1. The estimated state for the Lotka-Volterra model can be computed by using the bilinear Kalman smoother of equation (14). By applying the equations of this estimator to estimate the state x j , the error between the true state and the estimated state is found to be very small, as shown in  FIG. 3 . These results show that the present bilinear Kalman filter and bilinear Kalman smoother work well and are applicable to a bilinear model. 
         [0055]    The utility of the bilinear Kalman filter and bilinear Kalman smoother have also been simulated for estimation of bilinear system parameters via the EM approach. In this simulation, the Lotka-Volterra state space model is, once again, considered, with the parameters A, B, C, W and V being unknown. The initial value for the state in this simulation is 
         [0000]    
       
         
           
             
               
                 x 
                 0 
               
               = 
               
                 ( 
                 
                   
                     
                       1 
                     
                   
                   
                     
                       1 
                     
                   
                 
                 ) 
               
             
             , 
           
         
       
     
         [0000]    such that x 0 ˜N(μ, P). The additive noise w k  and v k  are uncorrelated with w k ˜(0, W) and v k ˜(0, V), with initial estimates of W and V. The estimation for the system of the model given by equations (23) and (24) is determined from this information about the model via the EM algorithm. The simulation begins with initial guesses for such parameters, and these are updated recursively until convergence to the true system; i.e., if the estimated parameters are very close to the true parameters, then a small error is obtained between the estimated state and the true state. The E-step and M-step are both outlined above. 
         [0056]      FIG. 4  shows the error between the true state and the estimated state using the estimated values of the parameters that were obtained from the simulation via the bilinear EM algorithm. The resulting small errors indicate that using the bilinear EM approach to estimate the exact values of the parameters is very reliable. 
         [0057]    The extended Kalman filter (EKF) and the ensemble Kalman filter (EnKF) are conventional prior art filters used for nonlinear models. Application of these filters to the nonlinear Lotka-Volterra model produces relatively good results. However, when compared against the present bilinear Kalman filter, the error between the true states and the estimated states is found to be smaller with the present bilinear Kalman filter than either of the prior art filters. The error between the true states and the estimated states for the present bilinear Kalman filter is shown in  FIG. 5 . The same error for the extended Kalman filter is shown in  FIG. 6 , and the error for the ensemble Kalman filter is shown in  FIG. 7 . Comparing the errors, the error is found to be smallest with the present bilinear Kalman filter, thus showing that the bilinear Kalman filter produces more accurate results than either conventional prior art variant on the Kalman filter. 
         [0058]    It should be understood that the calculations of state observer  100  may be performed by any suitable computer system, such as that diagrammatically shown in  FIG. 1 . Data is fed into system  100  by sensors  202 , and may further be input via any suitable type of user interface  116 , and may be stored in memory  112 , which may be any suitable type of computer readable and programmable memory. Calculations are performed by processor  114 , which may be any suitable type of computer processor and are fed, as control signals, to controls  208 . Controls  208  may also include a display for the user, which may be any suitable type of computer display. 
         [0059]    Processor  114  may be associated with, or incorporated into, any suitable type of computing device, for example, a personal computer or a programmable logic controller. The controls  208 , the processor  114 , the memory  112  and any associated computer readable recording media are in communication with one another by any suitable type of data bus, as is well known in the art. 
         [0060]    Examples of computer-readable recording media include a magnetic recording apparatus, an optical disk, a magneto-optical disk, and/or a semiconductor memory (for example, RAM, ROM, etc.). Examples of magnetic recording apparatus that may be used in addition to memory  112 , or in place of memory  112 , include a hard disk device (HDD), a flexible disk (FD), and a magnetic tape (MT). Examples of the optical disk include a DVD (Digital Versatile Disc), a DVD-RAM, a CD-ROM (Compact Disc-Read Only Memory), and a CD-R (Recordable)/RW. 
         [0061]    Alternatively, the present Kalman filter may be implemented by a digital signal processor, a microcontroller, an application specific integrated circuit (ASIC), or any other suitable circuit or device programmed or configured to carry out the steps of the method. 
         [0062]    It is to be understood that the present invention is not limited to the embodiments described above, but encompasses any and all embodiments within the scope of the following claims.