Approximating Kalman filtering in large dynamical systems Harri - - PowerPoint PPT Presentation

Approximating Kalman filtering in large dynamical systems

Harri Auvinen, Heikki Haario and Tuomo Kauranne Department of Applied Mathematics Lappeenranta University of Technology tuomo.kauranne@lut.fi

4 January 2006

– Typeset by FoilT EX –

Outline

Introduction to data assimilation Kalman filter Extended Kalman filter 3D-Var 4D-Var Variational Kalman filter Simulated assimilation results Conclusions

– Typeset by FoilT EX – 1

Introduction to data assimilation

The formulation of the general data assimilation (state estimation) problem for discrete time steps t = 1, 2, 3, ..., n contains an evolution or prediction equation and an observation equation: x(t) = Mt(x(t − 1)) + Et (1) y(t) = Kt(x(t)) + et, (2) where Mt is a model of the evolution operator, xt is the state of the process at time t and Et is a vector valued stochastic process. The second equation connects the state estimate x(t) to the measurements y(t), with their associated observation errors et, by way

• f an observation operator Kt.

– Typeset by FoilT EX – 2

Popular assimilation methods:

Kalman filter – Extended Kalman filter (EKF) – Fast Kalman filter (FKF) – Reduced Rank Kalman filter (RRKF) – Ensemble Kalman filter (EnKF) 3D-Var 4D-Var (used in operational weather forecasting)

– Typeset by FoilT EX – 3

Basic linear Kalman filter

The basic Kalman filter operates on a linear version of the general assimilation problem, equations (1) and (2): x(t) = Mtx(t − 1) + Et (3) y(t) = Ktx(t) + et, (4) where Mt is the linear evolution operator and similarly Kt is the linear

• bservation operator. The Kalman filter operates sequentially in time.

– Typeset by FoilT EX – 4

Kalman filter algorithm (1/2)

Let xest(t − 1) be an estimate of state x(t − 1) and Sest(t − 1) be the corresponding error covariance matrix of the estimate. At time t the evolution operator is used to produce an a priori estimate xa(t) and its covariance Sa(t): xa(t) = Mtxest(t − 1) (5) Sa(t) = MtSest(t − 1)MT

t + SEt,

(6) where SEt is the covariance of the prediction error Et.

– Typeset by FoilT EX – 5

Kalman filter algorithm (2/2)

The next step is to combine xa(t) with the observations y(t) made at time t to construct an updated estimate of the state and its covariance: Gt = Sa(t)KT

t (KtSa(t)KT t + Set)−1

(7) xest(t) = xa(t) + Gt(y(t) − Ktxa(t)) (8) Sest(t) = Sa(t) − GtKtSa(t), (9) where Gt is the Kalman gain matrix, which is functionally identical to the maximum a posteriori estimator. In a more general case, when the evolution model and/or the observation model is non-linear, the Extended Kalman Filter (EKF) is required.

– Typeset by FoilT EX – 6

Extended Kalman filter algorithm (1/3)

The filter uses the full non-linear evolution model equation (1) to produce an a priori estimate: xa(t) = Mt(xest(t − 1)) . In order to obtain the corresponding covariance Sa(t) of the a priori information, the prediction model is linearized about xest(t − 1): Mt = ∂Mt(xest(t − 1)) ∂x (10) Sa(t) = MtSest(t − 1)MT

t + SEt.

(11) The linearization in equation (10) is produced by taking finite differences between evolved small state vector perturbations, which is a computationally very expensive operation for large models, unless an

– Typeset by FoilT EX – 7

adjoint code to the model M has been produced. In both cases, integrating all the columns in M forward in time in equation (11) is expensive.

– Typeset by FoilT EX – 8

Extended Kalman filter algorithm (2/3)

The observation operator is linearized at the time of the observation about the a priori estimate xa(t) in order to obtain Kt, which is then used for calculating the gain matrix: Kt = ∂Kt(xa(t)) ∂x (12) Gt = Sa(t)KT

t (KtSa(t)KT t + Set)−1.

(13)

– Typeset by FoilT EX – 9

Extended Kalman filter algorithm (3/3)

After this, the full non-linear observation operator is used to update xa(t) and this is then used to produce a current state estimate and the corresponding error estimate: xest(t) = xa(t) + Gt(y(t) − Kt(xa(t))). (14) Sest(t) = Sa(t) − GtKtSa(t). (15) If the linearization of the observation operator at xa(t) is not good enough to construct xest(t), it is necessary to carry out some iterations of the last four equations.

– Typeset by FoilT EX – 10

Three-dimensional variational method (3D-Var)

The idea of 3D-Var method is to avoid the computation of the gain completely by looking for the analysis as an approximate solution to the equivalent minimization problem defined by the cost function: J(xt) = (xa(t) − xt)T(Sa(t))−1(xa(t) − xt) (16) + (yt − K(xt))TSe−1

t (yt − K(xt)),

(17) where K is the nonlinear observation operator. Opposite to the EKF method, the solution is retrieved iteratively by performing several evaluations of the cost function and of its gradient using a suitable descent algorithm. The approximation lies in the fact that

• nly a small number of iterations are performed.

– Typeset by FoilT EX – 11

Illustration of the variational cost-function minimization

– Typeset by FoilT EX – 12

The four-dimensional variational method (4D-Var)

4D-Var is a generalization of the less expensive 3D-Var method and it also supersedes earlier statistical interpolation methods, which do not take into account model dynamics. In its general form, it is defined as the minimization of the following cost function: J(x0) = Jb + Jo (18) = (xb − x0)TSa(t0)−1(xb − x0) (19) +

T

• t=0

(y(t) − Kt(Mt(x0)))TSe−1

t (y(t) − Kt(Mt(x0))).(20)

4D-Var method assumes that the model is perfect over the assimilation period T. This leads to the so-called strong constraint formalism.

– Typeset by FoilT EX – 13

Illustration of the 4D-Var assimilation

– Typeset by FoilT EX – 14

Approximating EKF

(Strong constraint) 4DVAR: Sest(t) is a static background error matrix Sa(t0) and there is no model error Sest(t) ≡ Sa(t0)

– Typeset by FoilT EX – 15

Approximating EKF

Reduced Rank Kalman Filtering RRKF: RRKF behaves like the Extended Kalman Filter in a low dimensional subspace determined at the beginning of the assimilation window Sest(t) = MtLT ˆ S F T F I

• LM T

t

L is the orthogonal matrix that transforms model variables into control

• variables. ˆ

S is the model error covariance matrix in the chosen k-dimensional subspace of the control space and F is its cross-correlation matrix with the complement subspace

– Typeset by FoilT EX – 16

Approximating EKF

Fast Kalman Filtering FKF: FKF produces results numerically identical to the full Kalman filter EKF but has one order of magnitude lower complexity FKF complexity is still superlinear, but it suits medium sized systems well

– Typeset by FoilT EX – 17

Approximating EKF

Ensemble Kalman Filtering EnKF: A low dimensional approximation En − Sest to Sa(t0) is propagated by the full nonlinear model and there is no model error En − Sest(t) =

• si∈V

(M(si)(t) − M(si)(t))(M(si)(t) − M(si)(t))T V is a sample of vectors, such as a multinormal sample modulated by a low rank approximation to the analysis error covariance matrix Sa(t0) at initial time

– Typeset by FoilT EX – 18

Variational Kalman filter

The idea of the VKF method is to replace the most time consuming part of the Extended Kalman Filter (EKF) method, the updating of the prediction error covariance matrix Sest(t − 1) equation (11) in order to

• btain Sa(t), by a low rank update that only requires a few dozen

model integrations. The Variational Kalman Filter can be viewed as a predictor corrector scheme in the Cartesian sum of the state and covariance spaces. In the scheme, a 3D-Var method is used to find both the state estimate and the error covariance matrix. Dominant Hessian singular vectors at time t + 1 are approximated by evolving the search directions from time t. The rank of the error covariance approximation can be regulated.

– Typeset by FoilT EX – 19

Variational Kalman filter algorithm (1/3)

The VKF method uses the full non-linear prediction equation (1) to construct an a priori estimate from the previous state estimate: xa(t) = Mt(xest(t − 1)). (21) The corresponding approximated covariance Sa(t) of the a priori information is available from previous time step of VKF method. In order to avoid the computation of the Kalman gain we perform a 3D-Var optimization with equivalent cost function. As a result of the

• ptimization we get the state estimate xest(t) for present time t.

– Typeset by FoilT EX – 20

Variational Kalman filter algorithm (2/3)

The error estimate Sest(t) of the state estimate follows from the formula: Sest(t) = 2(Hess(t))−1, (22) where the matrix Hess(t) can be approximated by using search directions

• f the optimization process. During the optimization another task is done:

the full non-linear evolution model equation (1) is used to transfer search directions to next time step t + 1. These evolved search directions are then used to update the approximation of the present covariance Sa(t) in order to approximate Sa(t + 1).

– Typeset by FoilT EX – 21

Variational Kalman filter algorithm (3/3)

In practice the updates are carried out by using the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) formula in order to maintain an upper limit on the rank of the Hessian approximation used in the L-BFGS Quasi-Newton method. After a 3D-Var optimization we have an updated covariance ˆ Sa(t) which is then used to produce: Sa(t + 1) = ˆ Sa(t) + SEt. (23) At each iteration of the optimization method the evolved search directions and the evolved gradients of J are stored and the approximation of the ˆ Sa(t) is updated as the optimization method updates it’s own approximation of the Hessian ∇∇J.

– Typeset by FoilT EX – 22

Simulated assimilation results

Simple linear evolution model 100 observation times 100 grid points data is obtained at the last 2 grid points in each set of 5 (i.e. the

• bserved grid point indexes were 4, 5, 9, 10, 14, 15, 19, 20, 24, 25, ....)

– Typeset by FoilT EX – 23

Simulated assimilation results (1/4)

10 20 30 40 50 60 70 80 90 100 −10 −5 5 10 15

Relative error: 100*(EKF−True)./True Relative error [%] time t

– Typeset by FoilT EX – 24

Simulated assimilation results (2/4)

10 20 30 40 50 60 70 80 90 100 −10 −5 5 10 15

Relative error: EKF vs. VKF Relative error [%] time t

Red lines: EKF Blue lines: VKF – Typeset by FoilT EX – 25

Simulated assimilation results (3/4)

10 20 30 40 50 60 70 80 90 100 −20 −15 −10 −5 5 10 15 20

Relative error: 100*(EKF−True)./True Relative error [%] time t

– Typeset by FoilT EX – 26

Simulated assimilation results (4/4)

10 20 30 40 50 60 70 80 90 100 −20 −15 −10 −5 5 10 15 20

Relative error: EKF vs. VKF Relative error [%] time t

Red lines: EKF Blue lines: VKF – Typeset by FoilT EX – 27

Conclusions

EKF is computationally unfeasible for large dynamical systems FKF is a good, exact supplement to EKF for medium sized systems VKF is a good approximation to EKF for large systems It allows model error to be incorporated It has a dynamic error covariance matrix It retains the linear complexity of 4DVAR

– Typeset by FoilT EX – 28