Reduced-Rank Unscented Kalman Filtering Using Cholesky-Based ...
2008 American Control Conference Westin Seattle Hotel, Seattle, Washington, USA June 11-13, 2008
WeB17.4
Reduced-Rank Unscented Kalman Filtering Using Cholesky-Based Decomposition J. Chandrasekar, I. S. Kim, D. S. Bernstein, and A. J. Ridley I. I NTRODUCTION Data assimilation for large-scale systems has gained increasing attention due to nonlinear and computationally intensive applications such as weather forecasting [1, 2]. These problems require algorithms that are computationally tractable despite the enormous dimension of the state. Reduced-order variants of the classical Kalman filter have been developed [3, 4]. A comparison of several techniques is given in [5]. An alternative technique for reducing the computational requirements of data assimilation for high-dimensional systems is the reduced-rank filter [6, 7]. In this method, the error-covariance matrix is factored to obtain a square root, whose rank is then reduced through truncation. The truncated square root is then propagated by the data assimilation algorithm. The primary technique for truncating the error-covariance matrix is the singular value decomposition (SVD), wherein the singular values are used to determine which entries of the error covariance matrix are most relevant to the accuracy of the state estimates [6, 7]. In related work [8], it is observed that the Kalman filter estimate update depends on the product Ck Pk , where Ck is the measurement map and Pk is the error covariance. In particular, it is shown in [8] that approximation of Ck Pk leads directly to truncation based on the Cholesky decomposition. Filter reduction based on the Cholesky decomposition provides state-estimation accuracy that is competitive with, and in many cases superior to, that of the SVD. An additional advantage of using the Cholesky decomposition in place of the SVD for reducedrank filtering is the fact that the Cholesky decomposition is computationally less expensive than the SVD. To assimilate data in nonlinear systems, particle filters are used to propagate a collection of state estimates from which statistics can be computed. These techniques include the ensemble Kalman filter (EnKF) [10, 11], which uses a stochastic construction, as well as the unscented Kalman filter (UKF) [12, 13], which deterministically constructs the collection of state estimates by perturbing the nominal state estimate. Specifically, UKF constructs the ensemble members by using the columns of the square root of the error covariance to perturb the nominal state estimate. For a model of order n, the n columns and their negatives result in 2n+1 ensemble members and thus 2n + 1 model updates. This research was supported by the National Science Foundation, under grants CNS-0539053 and ATM-0325332 to the University of Michigan, Ann Arbor, USA. J. Chandrasekar, I. S. Kim, and D. S. Bernstein are with The University of Michigan, Ann Arbor, MI, [email protected]
978-1-4244-2079-7/08/$25.00 ©2008 AACC.
A straightforward approach to reducing the UKF ensemble size is to use a factorization-and-truncation method to truncate n − q columns of the square root of the error covariance and construct the ensemble members using the remaining q columns. In [14], SVD-based decomposition-and-truncation is used to construct reduced-rank approximations of the square root of the error covariance, which are then used to construct the ensemble members resulting in an ensemble of size 2q + 1. In this paper, we use the Cholesky-based decomposition technique developed in [8] to construct the reduced-ensemble members. Specifically, we use the Cholesky decomposition to obtain a square root of the error covariance and select columns of the Cholesky factor to approximate Ck Pk . The retained columns of the Cholesky factor are used to construct the ensemble members. We compare the performance of the Cholesky-decomposition-based reduced-rank UKF and the SVD-based reduced-rank UKF on a linear advection model and a nonlinear system with chaotic dynamics. II. T HE U NSCENTED K ALMAN F ILTER We consider the discrete-time system with nonlinear dynamics xk+1 = f (xk , uk , k) + wk
(2.1)
and linearly dependent measurements yk = Ck xk + vk , n
m
(2.2) p
where xk , wk ∈ R , uk ∈ R , and yk , vk ∈ R . The input uk and output yk are assumed to be measured, and wk and vk are uncorrelated zero-mean white noise processes with covariances Qk and Rk , respectively. We assume that Rk is positive definite. The objective is to obtain estimates of the state xk using measurements yk . When the dynamics (2.1) are linear, the Kalman filter provides estimates that minimize the mean-square-error (MSE) in the state estimates [16]. Let x ∈ Rn , and let P ∈ Rn×n be positive semidefinite. The unscented transformation provides 2n + 1 ensembles Xi ∈ Rn and corresponding weights γi ∈ R, for 0 = 1, . . . , 2n, such that the weighted mean and weighted variance of the ensembles are x and P , respectively. Specifically, let S ∈ Rn×n satisfy SS T = P , and, for all i = 1, . . . , n, let coli (S) denote the ith column of S. For α > 0, the unscented transformation X = Ψ(x, S, α) ∈ Rn×(2n+1) of x with covariance P = SS T is defined by £ ¤ X , X0 · · · X2n , (2.3)
1274
where x, √ Xi , x + αcoli (S), √ x − αcoli−n (S),
i = 0, i = 1, . . . , n, i = n + 1, . . . , 2n.
III. SVD-BASED R EDUCED -R ANK U NSCENTED K ALMAN F ILTER (2.4)
The parameter α determines the spread of the ensembles around x. Next, define the weights γi ∈ R by 1 α−n , γi , , i = 1, . . . , 2n. (2.5) γ0 , α 2α Then, 2n 2n X X γi Xi = x, γi (Xi − x)(Xi − x)T = P. (2.6) i=0
i=0
UKF uses the unscented transformation to approximate the error covariance and estimate the state xk . Letting xf0 denote an initial estimate of x0 with error covariance P0f , UKF is given by the following steps. UKF data assimilation step: f f xda k = xk + Kk (yk − yk ),
(2.7)
= Ck xfk , (2.8) ¤ £ da da da X · · · X , (2.9) = Ψ(xda , S , α) = k k 0,k 2n,k f f = Sk Hk , (2.10) f T f T −1 Kk = Pk Ck (Ck Pk Ck + Rk ) , (2.11) f n×n where Hk ∈ R satisfies f f T Hk (Hk ) =In−(Ck Skf )T (Ck Skf (Ck Skf )T+Rk )−1Ck Skf(2.12) ykf Xkda Skda
To reduce the ensemble size, we use a reduced-rank apf f proximation Pˆs,k of Ps,k . The reduced-rank approximations f f are chosen such that kPˆs,k − Ps,k kF is minimized subject to f ˆ rank(Ps,k ) = q, where k·kF denotes the Frobenius norm. Let P ∈ Rn×n be positive semidefinite, let σ1 > · · · > σn > 0 be the singular values of P , and let u1 , . . . , un ∈ Rn be the corresponding mutually orthogonal singular vectors. Next, define Uq ∈ Rn×q and Σq ∈ Rq×q by £ ¤ Uq , u1 · · · uq , Σq , diag(σ1 , . . . , σq ). (3.1)
With this notation, the singular value decomposition of P is given by P = Un Σn UnT , where Un ∈ Rn×n is orthogonal. For q 6 n, let ΦSVD (P, q) ∈ Rn×q denote the SVD1/2 based rank-q approximation of the square root Un Σn of 1/2 P given by ΦSVD (P, q) , Uq Σq . As noted in [18], Pˆ = SS T , where S , ΦSVD (P, q), is the best rank-q approximation of P with respect to the Frobenius norm. The SVD-based reduced-rank square-root unscented Kalman filter (SVDRRUKF) is given by the following steps. SVDRRUKF data assimilation step: f f xda s,k = xs,k + Ks,k (yk − ys,k ),
f ys,k da Xs,k da Ss,k
and Skf ∈ Rn×n satisfies
Skf (Skf )T = Pkf .
(2.13)
Finally, define Pkda ∈ Rn×n by
Pkda , Skda (Skda )T .
(2.14)
Ks,k
=
xfk+1
=
f = Pk+1 f ∆Xi,k
da f (Xi,k , uk , k), i = 0, . . . , 2n, 2n X f γi Xi,k+1 , i=0 2n X f f (∆Xi,k+1 )T + Qk , γi ∆Xi,k+1 i=0 f Xi,k
xfk .
= =
£ da ¤ da = Xs,0,k · · · Xs,2q,k ,
(3.3) (3.4)
(3.5) ¢ ¡ −1 f f f f )T + Rk , (3.6) (Ck Ss,k )T Ck Ss,k (Ck Ss,k = Ss,k
SVDRRUKF forecast step:
f da Xs,i,k+1 = f (Xs,i,k , uk , k), i = 0, . . . , 2q,
(2.15) xfs,k+1 = (2.16) f Ps,k+1 =
(2.17)
2q X
i=0 2q X
(3.7)
f , γq,i Xs,i,k+1
(3.8)
f f γq,i ∆Xs,i,k+1 (∆Xs,i,k+1 ) T + Qk ,
(3.9)
i=0
f f Ss,k+1 = ΦSVD (Ps,k+1 , q).
where , − When the dynamics in (2.1) are linear, UKF is equivalent to the Kalman filter [12]. Furthermore, in the linear case, Pkda and Pkf are the covariances of f the error xk − xda k and xk − xk , respectively. However, in f da the nonlinear case, Pk and Pk are pseudo-error covariances. The case in which the measurements depend nonlinearly on the state is discussed in [13]. Note that Hkf and Skf satisfying (2.12) and (2.13) are not unique. Moreover, all square Hkf and Skf satisfying (2.12) and (2.13) are related by an orthogonal transformation. Specifically, the following result is given in [18, p. 188]. Lemma 2.1: Let S, Sˆ ∈ Rn×n . Then, SS T = SˆSˆT if and only if there exists an orthogonal matrix U ∈ Rn×n such that Sˆ = SU .
(3.2)
f f where Hs,k ∈ Rq×q satisfies (2.12) with Skf replaced by Ss,k .
UKF forecast step: f Xi,k+1
=
Ck xfs,k , da Ψq (xda s,k , Ss,k , α) f f Ss,k Hs,k ,
f ∆Xs,i,k
f Xs,i,k
(3.10)
xfs,k .
, − f da Next, define Pˆs,k , Pˆs,k ∈ Rn×n by
where
da T da da f f f ) . (Ss,k , Ss,k )T , Pˆs,k (Ss,k , Ss,k Pˆs,k
(3.11)
It then follows from (3.5) that da f f f f Pˆs,k = Pˆs,k − Pˆs,k CkT (Ck Pˆs,k CkT + Rk )−1 Ck Pˆs,k . (3.12) Furthermore, (3.6) and (3.12) imply that f f CkT + Rk )−1 . Ks,k = Pˆs,k CkT (Ck Pˆs,k
(3.13)
Since the SVD in (3.10) is computationally intensive [9], we introduce an alternative method to obtain a reducedrank approximation of a square root of the pseudo-error covariance.
1275
CDRRUKF forecast step:
IV. C HOLESKY-D ECOMPOSITION -BASED R EDUCED -R ANK U NSCENTED K ALMAN F ILTER
f da Xc,i,k+1 = f (Xc,i,k , uk , k), i = 0, . . . , 2q
The filter gain Kk of UKF depends on a particular subspace of the forecast error covariance Pkf . Specifically, Kk depends only on the correlation Ck Pkf between the error in the measured states and unmeasured states. Since rank(Ck ) = p, there exists a state space basis with respect to which Ck has the form £ ¤ Ck = Ip 0 . (4.1) The following result is given in [8]. Lemma 4.1: Partition Pkf as · f ¸ f Pp,k (Ppp,k )T f Pk = , f f Ppp,k Pp,k
(4.2)
f f ∈ Rp×p , and assume that Ck where Pp,k ∈ Rp×p and Pp,k has the form (4.1). Then, · f ¸ Pp,k f (Pp,k + Rk )−1 . (4.3) Kk = f Ppp,k To reduce the ensemble size, we construct a filter that f f uses a reduced-rank approximation Pˆc,k of Pc,k such that f f f ˆ ˆ rank(Pc,k ) < n and kCk (Pc,k − Pc,k )kF is minimized. To f obtain Pˆc,k , we perform a Cholesky decomposition of the f pseudo-error covariance Pc,k at each time step. Assuming n×n that P ∈ R is positive definite, the Cholesky decomposition of P yields a unique lower triangular Cholesky factor L ∈ Rn×n satisfying
LLT = P.
£
Truncating the last n − q columns of L = L1 · · · yields the rank-q Cholesky factor £ ¤ ΦCHOL (P, q) , L1 · · · Lq ∈ Rn×q .
(4.4) ¤ Ln (4.5)
The following result is given in [8]. Lemma 4.2: Let P ∈ Rn×n be positive definite, define S , ΦCHOL (P, q), where 0 < q 6 n, and partition P and Pˆ , SS T as ¸ · · ¸ Pq Pqq Pˆqq Pˆq ˆ , P = P = , (4.6) (Pqq )T Pq (Pˆqq )T Pˆq
where Pq , Pˆq ∈ Rq×q and Pq , Pˆq ∈ Rq×q . Then, £ ¤ £ ¤ (4.7) Pˆq Pˆqq = Pq Pqq . Lemma 4.2 implies that, if S = ΦCHOL (P, q), then the first q columns and rows of SS T and P are equal. The Cholesky-decomposition-based reduced-rank unscented Kalman filter (CDRRUKF) is summarized as follows. CDRRUKF data assimilation step: f f xda c,k = xc,k + Kc,k (yk − yc,k ),
f yc,k da Xc,k da Sc,k
=
= =
Kc,k =
Ck xfc,k , da Ψq (xda c,k , Sc,k , α), f f Sc,k Hc,k , f f (Ck Sc,k )T Sc,k
(4.8) (4.9) (4.10)
(4.11) ¡ ¢−1 f f T , (4.12) Ck Sc,k (Ck Sc,k ) +Rk
f f where Hc,k ∈ Rq×q satisfies (2.12) with Skf replaced by Sc,k .
xfk+1 = f Pc,k+1 =
2q X
i=0 2q X
f γq,i Xc,i,k+1 ,
(4.13) (4.14)
f f γq,i ∆Xc,i,k+1 (∆Xc,i,k+1 )T + Qk , (4.15)
i=0
f f Sc,k+1 = ΦCHOL (Pc,k+1 , q), f ∆Xc,i,k
(4.16)
f = Xc,i,k − xfc,k . da ˆ f Pˆc,k , Pc,k ∈ Rn×n
where f f f Next, define by Pˆc,k , Sc,k (Sc,k )T and da f f f f Pˆc,k , Pˆc,k − Pˆc,k CkT (Ck Pˆc,k CkT + Rk )−1 Ck Pˆc,k . (4.17)
da da T da It then follows from (4.11) that Sc,k (Sc,k ) = Pˆc,k . Furthermore, (4.12) and (4.17) imply that Kc,k = Pˆ f C T (Ck Pˆ f C T + Rk )−1 . (4.18) c,k
k
c,k
k
Hence, like the estimator gain Ks,k of SVDRRUKF, the estimator gain Kc,k of CDRRUKF depends on a reducedf rank approximation Pˆc,k of the pseudo-error covariance f Pc,k . Due to the rank-reduction step (4.16), CDRRUKF is generally not equivalent to UKF. However, we now discuss cases in which the performance of CDRRUKF is close to that of UKF. A. Basis Selection for CDRRUKF The following result given in [8] shows that CDRRUKF is equivalent to UKF for a single time step when Ck has the form (4.1). Proposition 4.1: Assume that p 6 n, q = p, and Ck has the structure in (4.1). Let Pkf ∈ Rn×n be positive semidefinite and Kk be given by (2.11). Furthermore, define f f Sc,k , ΦCHOL (Pkf , q) and let Pˆc,k be given by (4.17). Then, f f ˆ Ck Pc,k = Ck Pk and hence, Kc,k = Kk . If the dynamics (2.1) are linear and time-invariant, that is, for all k > 0, xk+1
=
Axk + Buk + wk ,
(4.19)
yk
=
Cxk + vk ,
(4.20)
then a basis for the state x can be chosen so that CDRRUKF is equivalent to UKF for r > 0 time steps. To construct such a basis, we define the observability matrix O(A, C) ∈ Rpn×n by ¤T △ £ O(A, C) = C T (CA)T · · · (CAn−1 )T . (4.21) The following result is given in [8]. Proposition 4.2: Assume that O(A, C) has the form ¸ · In . (4.22) O(A, C) = 0(p−1)n×n Let r > 0 be an integer such that pr < n, and let q = f pr. Furthermore, assume that Pc,0 = P0f . Then, for all k = 0, . . . , r, Kc,k = Kk . If, in addition, xfc,0 = xf0 , then for all k = 0, . . . , r, xfc,k = xfk . £ ¤T x1,k · · · xn,k Let xk have entries xk = . If O(A, C) has the form (4.22), and uk = wk = vk = 0,
1276
for all k > 0, then for all k > 0, and for every integer r > 0 such that pr 6 n, £ T ¤T £ T ¤T T yk · · · yk+r−1 = x1,k · · · xT . (4.23) pr,k
Next, we consider the nonlinear system (2.1) and assume that the dynamics in (2.1) can be expressed as xi,k+1 = fi (xφL (i,b) , . . . , xφR (i,b) , uk , k) + wi,k ,
(4.24)
where i = 1, . . . , n, b > 0 and φL (i, b) , max(1, i − b), φR (i, b) , min(n, i + b). (4.25) For example, in systems modeled by finite volume schemes, the next value of a physical variable in a given cell depends only on the present values of the physical variables in its neighboring cells. Next, let yk denote a measurement of a specific component of the state, so that yk = xj,k +vk , where j ∈ {1, . . . , n}, j− rb > 1, and j + rb 6 n. It follows from (4.24) that, if wk = vk = 0, for all k > 0, then yk , . . . , yk+r−1 depends on only first 2rb components of the state vector x ˜k at time step k. Hence, a state space basis can be chosen such that the outputs yk , . . . , yk+r−1 depend on only the first few components of the state vector. V. L INEAR A DVECTION M ODEL Consider a linear advection model [2] with n cells, and let xi,k be the energy in the ith cell at time k. The energy flow satisfies ( xi−1,k , if i = 2, . . . , n, (5.1) xi,k+1 = xn,k , if i = 1. Hence, energy in the ith cell flows to the (i + 1)th cell, while the periodic boundary condition ensures that the energy circulates continually. We choose n = 100 and assume that the disturbance wk enters selected cells, where wk ∈ Rn is a white noise process with covariance Qk = Q for all k > 0, and Q ∈ Rn×n is diagonal with nonzero entries Qi = 1 only for i ∈ {10, 20, . . . , 100}. Next, we assume that measurements of the energy in cells 50 and 51 are available so that yk = [x50,k x51,k ]T + vk , where vk is white noise process with covariance Rk = 0.1I2 . First, we use the measurements yk to estimate the energy in the remaining cells using UKF. In all three cases, the initial estimates xf0 , xfs,0 , and xfc,0 are not equal to the initial f f state x0 . Moreover, we choose P0f = Ps,0 = Pc,0 = 0.1In . Finally, we choose α = 0.6 for all three filters. As shown in Figure 1 and Figure 2, data assimilation is performed using SVDRRUKF and CDRRUKF for several values of q between 5 and 100. It can be seen that the performance of SVDRRUKF with 111 ensemble members (q = 55) is close to optimal, whereas the performance of CDRRUKF is close to optimal with 11 ensemble members (q = 5). The steady-state MSE of state estimates for various values of q is plotted in Figure 3 and Figure 4. The performance of SVDRRUKF is poor when q < 55, and close to optimal when q > 55. Thus the ensemble size can be reduced from 201 to 111 with negligible change in the
performance. However, the ensemble size can be reduced from 211 to 11 with negligible performance deterioration. Next, we repeat the same procedure except with a poor estimate of the process noise covariance for data assimilation. ˆ k , where Specifically, we replace Qk in (3.9) and (4.15) by Q ˆ Qk = I for all k > 0. The steady-state MSE of state estimates for different choices of q is plotted in Figure 3 and Figure 4. SVDRRUKF with a poor estimate of the error covariance is unstable for all q 6 95 (indicated by the X’s). However, Figure 4 shows that even with q = 5 and a poor estimate of the process noise covariance, the performance of CDRRUKF is close to optimal. ˆ k , where Q ˆ k = αI Finally, we replace Qk in (4.16) with Q for all k > 0, and perform state estimation using CDRRUKF. The steady-state MSE of the state estimates is shown in Figure 5 for several values of α. For all the cases, the performance of CDRRUKF is close to optimal when α > 1. This suggests that it is advantageous to overestimate the process noise covariance. SVDRRUKF is unstable for all choices of α = 0.005, . . . , 50. Hence, these simulations suggest that CDRRUKF is more robust than SVDRRUKF with respect to uncertainties in the process noise covariance. All of the results for CDRRUKF in figures 1-6 are obtained using a state space basis with respect to which the observability matrix has the form (4.22). VI. L96 M ODEL The L96 model mimics the propagation of an unspecified meteorological quantity along a latitude circle [17]. The dynamics are governed by d xi (t) = (xi+1 (t) − xi−2 (t))xi−1 (t) − xi (t) + ui (t), (6.1) dt where xi (t) ∈ R denotes the meteorological quantity at the ith grid point at time t, ui ∈ R denotes an external forcing term, and wi denotes unknown disturbances affecting the ith grid point. For all t > 0, the boundary conditions are defined by x0 (t) = xn (t), x−1 (t) = xn−1 (t), xn+1 (t) = x1 (t). We choose ui (t) = 8 for all i = 1, . . . , n and all t > 0. Using fourth-order Runge-Kutta discretization with a sampling time of 0.05 s, we obtain a discrete-time model of (6.1) that can be expressed as (2.1). Furthermore, we assume that the discretized model is corrupted by an unknown external disturbance that affects specified cells. We choose n = 40, and assume that wk is white noise process with covariance Qk = Q for all k > 0, where Q ∈ Rn×n is diagonal with nonzero entries Qi,i = 0.1 only for i ∈ {5, 15, 25, 35}. Next, we assume that measurements from cells with 20 and 21 are T available so that yk = [x20,k x21,k ] + vk , where vk is white noise process with covariance Rk = 0.01I2 . We use the measurements yk to estimate the state in the cells where measurements are not available. Next, as shown in Figure 6 and Figure 7, we reduce the ensemble size and use SVDRRUKF and CDRRUKF with q = 10, 20, 30. It can be seen that the performance of SVDRRUKF is poor compared to the performance of CDRRUKF for both q = 20 and q = 30. Moreover, the performance of CDRRUKF with
1277
3.6
log10(MSE of the state estimates)
3.2
3.6
SVDRRUKF q=5 SVDRRUKF q=15 SVDRRUKF q=25 SVDRRUKF q=35 SVDRRUKF q=45 SVDRRUKF q=55 UKF
3.2
3
2.8
2.6
2.4
2.2
3
2.8
2.6
2.4
2.2
2
2
1.8
1.8
1.6
50
100
CDRRUKF q=5 UKF
3.4
log10(MSE of the state estimates)
3.4
150 200 Time step k
250
300
1.6
350
Fig. 1. MSE of the state estimates obtained from SVDRRUKF for several values of q. SVDRRUKF with q = 5 is unstable, while the performance of SVDRRUKF with q = 55 is close to the performance of full-order UKF.
250
300
350
SVDRRUKF SVDRRUKF with a poor estimate of Qk
Steady−state MSE of the state estimates
2500
2000
1500
1000
UKF 500
0
X
X X 5
15
X
X
X
X
X
X
X
25 35 45 55 65 75 85 q : Rank of the approximation of the error covariance
X 95
100
Fig. 3. Steady-state performance of SVDRRUKF. For each value of q, we perform data assimilation with the true value of the process noise covariance and with a poor estimate of the process noise covariance. The X’s indicate cases in which the filter is unstable.
3000 CDRRUKF CDRRUKF with a poor estimate of Qk 2500
2000
1500
1000
UKF 500
0
R EFERENCES [1] J. Lewis, S. Lakshmivarahan, and S. Dhall, Dynamic Data Assimilation: A Least Squares Approach, Cambridge, 2006. [2] G. Evensen, Data Assimilation: The Ensemble Kalman Filter, Springer, 2006. [3] S. F. Farrell and P. J. Ioannou, “State Estimation Using a ReducedOrder Kalman Filter,” J. Atmos. Sci., vol. 58, pp. 3666–3680, 2001.
150 200 Time step k
3000
Steady−state MSE of the state estimates
We developed a reduced-rank square-root unscented Kalman filter based on the Cholesky decomposition of the pseudo-error covariance. We compare the performance of the Cholesky-based filter with an analogous filter that uses the singular value decomposition for a linear advection model and a nonlinear system that exhibits chaotic behavior. Although the computational requirement of the Choleskybased reduced-rank filter is less than that of the SVDbased reduced-rank filter, the results presented here suggest that the estimation accuracy of the Cholesky-based reducedrank filter is significantly better than that of SVD-based reduced-rank filter. Moreover, simulation results suggest that the Cholesky-based reduced-rank filter is more robust to uncertainties in the process noise covariance than the SVDbased reduced-rank filter.
100
Fig. 2. MSE of the state estimates obtained from CDRRUKF with q = 5. The performance of CDRRUKF with q = 5 is close to the full-order UKF performance.
61 (q = 30) ensemble members is close to the performance of UKF with 81 ensemble members. ˆ k , where Next, we replace Qk in (3.9) and (4.15) by Q ˆ k = αI for all k > 0. Figure 8 shows the time-averaged Q MSE of state estimates obtained using SVDRRUKF and CDRRUKF with q = 10 and q = 20 for several values of α between 0.001 and 100. It can be seen that, for all values of α, the performance of CDRRUKF is superior to the performance of SVDRRUKF. In fact, CDRRUKF with 21 ensemble members (q = 10) consistently outperforms SVDRRUKF with 41 ensemble members (q = 20). Finally, note that (6.1) can be expressed as (4.24) with b = 1. All of the results for CDRRUKF in figures 8-12 are obtained using the state space basis discussed in Section IV-A. VII. C ONCLUSIONS
50
5
15
25 35 45 55 65 75 85 q: Rank of the approximation of the error covariance
95
100
Fig. 4. Steady-state performance of CDRRUKF for values of q between 5 and 100. Note that for q = 5, the performance of CDRRUKF is close to optimal, irrespective of the value of the process noise covariance used for data assimilation.
1278
3.5 CDRRUKF q=5 CDRRUKF q=15 CDRRUKF q=25
4.5
No data assimilation CDRRUKF q=20 CDRRUKF q=30 UKF
3.3
log10(MSE of the state estimates)
log10(Steady−state MSE of the state estimates)
3.4
3.2
3.1
3
2.9
2.8
2.7
4
3.5
3
2.5
2
2.6 −2.5
−2
−1.5
−1
−0.5 0 log (α)
0.5
1
1.5
2
1.5
10
Fig. 5. Steady-state performance of CDRRUKF with q = 5, 15, 25. In all three cases, we use a poor estimate of the process noise covariance for data assimilation. For a fixed level of uncertainty in the process noise covariance, the performance of CDRRUKF improves when the ensemble size increases. The performance of SVDRRUKF is not shown since SVDRRUKF is unstable for all values of α and q = 5, 15, 25.
5
10
15
20
25 Time (sec)
30
35
40
45
50
Fig. 7. Performance of CDRRUKF with n = 40 and q = 20, 30. Note that the performance of CDRRUKF with q = 20 is better than the performance of SVDRRUKF with q = 30.
4 No data assimilation SVDRRUKF q=10 SVDRRUKF q=20 CDRRUKF q=10 CDRRUKF q=20 UKF
3.8 No data assimilation SVDRRUKF q=20 SVDRRUKF q=30 UKF
log10(Time−averaged MSE of the state estimates)
log10(MSE of the state estimates)
4.5
4
3.5
3
2.5
3.6
3.4
3.2
3
2.8
2.6
2.4
2.2 2
1.5
2 −3
5
10
15
20
25 Time (sec)
30
35
40
45
−2.5
−2
−1.5
−1
−0.5 log10(α)
0
0.5
1
1.5
2
50
Fig. 6. MSE of the state estimates obtained using SVDRRUKF with q = 20, 30. The error in state estimates based on UKF and data-free simulation is shown for comparison. SVDRRUKF with q = 20 and q = 30 sometimes yields estimates that are worse than estimates based on data-free simulation.
[4] L. Scherliess, R. W. Schunk, J. J. Sojka, and D. C. Thompson, “Development of a Physics-Based Reduced State Kalman Filter for the Ionosphere,” Radio Science, vol. 39, RS1S04, pp. 1-12, 2004. [5] I. S. Kim, J. Chandrasekar, H. J. Palanthandalam-Madapusi, A. Ridley, and D. S. Bernstein, “State Estimation for Large-Scale Systems Based on Reduced-Order Error-Covariance Propagation,” Proc. Amer. Contr. Conf., New York, pp. 5700-5705, June 2007. [6] M. Verlaan and A. W. Heemink, “Tidal Flow Forecasting Using Reduced-Rank Square-Root Filters,” Stochastic Hydrology and Hydraulics, vol. 11, pp. 349–368, 1997. [7] S. Gillijns, D. S. Bernstein, and B. D. Moor, “The Reduced Rank Transform Square Root Filter for Data Assimilation,” Proc. 14th IFAC Symposium on System Identification (SYSID2006), Newcastle, Australia, pp. 349–368, 2006. [8] J. Chandrasekar, I. S. Kim, and D. S. Bernstein, “Cholesky-Based Reduced-Rank Kalman Filtering,” submitted to Amer. Contr. Conf, 2008. [9] G. W. Stewart, Matrix Algorithms Volume 1: Basic Decompositions, SIAM, 1998. [10] G. Evensen, “The Ensemble Kalman Filter: Theoretical Formulation and Practical Implementation,” Ocean Dynamics, vol. 53, pp. 343-36,
Fig. 8. Time-averaged MSE of state estimates between 35 sec and 50 sec. The state estimates are obtained using SVDRRUKF and CDRRUKF with q = 10 and q = 20, and a poor estimate of the process noise covariance. For all values of α, the performance of CDRRUKF is better than the performance of SVDRRUKF.
2003. [11] J. S. Whitaker and T. M. Hamill, “Ensemble Data Assimilation without Perturbed Observations,” Mon. Wea. Rev., vol. 130, pp. 1913-1924, 2002. [12] S. Julier, J. Uhlmann, and H. F. Durrant-Whyte, “A New Method for the Nonlinear Transformation of Means and Covariances in Filters and Estimators,” IEEE Trans. Autom. Contr., vol. 45, pp. 477-482, 2000. [13] B. Ristic, S. Arulampalam and N. Gordon, Beyond the Kalman Filter, Artech House, 2004. [14] M. K. Tippett, J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker, “Ensemble Square Root Filters,” Mon. Wea. Rev., vol. 131, pp. 1485–1490, 2003. [15] D. S. Bernstein and D. C. Hyland, “Compartmental Modeling and Second-Moment Analysis of State Space Systems,” SIAM J. Matrix Anal. Appl., vol. 14, no. 3, pp. 880–901, 1993. [16] B. D. O. Anderson and J. B. Moore, Optimal Filtering, Prentice-Hall, 1979, reprinted by Dover, 2005. [17] E. N. Lorenz, “Predictability - A Problem Partly Solved,” Predictability of Weather and Climate, Cambridge University Press, 2006. [18] D. S. Bernstein, Matrix Mathematics, Princeton University Press, 2005.
1279
WeB17.4
Reduced-Rank Unscented Kalman Filtering Using Cholesky-Based Decomposition J. Chandrasekar, I. S. Kim, D. S. Bernstein, and A. J. Ridley I. I NTRODUCTION Data assimilation for large-scale systems has gained increasing attention due to nonlinear and computationally intensive applications such as weather forecasting [1, 2]. These problems require algorithms that are computationally tractable despite the enormous dimension of the state. Reduced-order variants of the classical Kalman filter have been developed [3, 4]. A comparison of several techniques is given in [5]. An alternative technique for reducing the computational requirements of data assimilation for high-dimensional systems is the reduced-rank filter [6, 7]. In this method, the error-covariance matrix is factored to obtain a square root, whose rank is then reduced through truncation. The truncated square root is then propagated by the data assimilation algorithm. The primary technique for truncating the error-covariance matrix is the singular value decomposition (SVD), wherein the singular values are used to determine which entries of the error covariance matrix are most relevant to the accuracy of the state estimates [6, 7]. In related work [8], it is observed that the Kalman filter estimate update depends on the product Ck Pk , where Ck is the measurement map and Pk is the error covariance. In particular, it is shown in [8] that approximation of Ck Pk leads directly to truncation based on the Cholesky decomposition. Filter reduction based on the Cholesky decomposition provides state-estimation accuracy that is competitive with, and in many cases superior to, that of the SVD. An additional advantage of using the Cholesky decomposition in place of the SVD for reducedrank filtering is the fact that the Cholesky decomposition is computationally less expensive than the SVD. To assimilate data in nonlinear systems, particle filters are used to propagate a collection of state estimates from which statistics can be computed. These techniques include the ensemble Kalman filter (EnKF) [10, 11], which uses a stochastic construction, as well as the unscented Kalman filter (UKF) [12, 13], which deterministically constructs the collection of state estimates by perturbing the nominal state estimate. Specifically, UKF constructs the ensemble members by using the columns of the square root of the error covariance to perturb the nominal state estimate. For a model of order n, the n columns and their negatives result in 2n+1 ensemble members and thus 2n + 1 model updates. This research was supported by the National Science Foundation, under grants CNS-0539053 and ATM-0325332 to the University of Michigan, Ann Arbor, USA. J. Chandrasekar, I. S. Kim, and D. S. Bernstein are with The University of Michigan, Ann Arbor, MI, [email protected]
978-1-4244-2079-7/08/$25.00 ©2008 AACC.
A straightforward approach to reducing the UKF ensemble size is to use a factorization-and-truncation method to truncate n − q columns of the square root of the error covariance and construct the ensemble members using the remaining q columns. In [14], SVD-based decomposition-and-truncation is used to construct reduced-rank approximations of the square root of the error covariance, which are then used to construct the ensemble members resulting in an ensemble of size 2q + 1. In this paper, we use the Cholesky-based decomposition technique developed in [8] to construct the reduced-ensemble members. Specifically, we use the Cholesky decomposition to obtain a square root of the error covariance and select columns of the Cholesky factor to approximate Ck Pk . The retained columns of the Cholesky factor are used to construct the ensemble members. We compare the performance of the Cholesky-decomposition-based reduced-rank UKF and the SVD-based reduced-rank UKF on a linear advection model and a nonlinear system with chaotic dynamics. II. T HE U NSCENTED K ALMAN F ILTER We consider the discrete-time system with nonlinear dynamics xk+1 = f (xk , uk , k) + wk
(2.1)
and linearly dependent measurements yk = Ck xk + vk , n
m
(2.2) p
where xk , wk ∈ R , uk ∈ R , and yk , vk ∈ R . The input uk and output yk are assumed to be measured, and wk and vk are uncorrelated zero-mean white noise processes with covariances Qk and Rk , respectively. We assume that Rk is positive definite. The objective is to obtain estimates of the state xk using measurements yk . When the dynamics (2.1) are linear, the Kalman filter provides estimates that minimize the mean-square-error (MSE) in the state estimates [16]. Let x ∈ Rn , and let P ∈ Rn×n be positive semidefinite. The unscented transformation provides 2n + 1 ensembles Xi ∈ Rn and corresponding weights γi ∈ R, for 0 = 1, . . . , 2n, such that the weighted mean and weighted variance of the ensembles are x and P , respectively. Specifically, let S ∈ Rn×n satisfy SS T = P , and, for all i = 1, . . . , n, let coli (S) denote the ith column of S. For α > 0, the unscented transformation X = Ψ(x, S, α) ∈ Rn×(2n+1) of x with covariance P = SS T is defined by £ ¤ X , X0 · · · X2n , (2.3)
1274
where x, √ Xi , x + αcoli (S), √ x − αcoli−n (S),
i = 0, i = 1, . . . , n, i = n + 1, . . . , 2n.
III. SVD-BASED R EDUCED -R ANK U NSCENTED K ALMAN F ILTER (2.4)
The parameter α determines the spread of the ensembles around x. Next, define the weights γi ∈ R by 1 α−n , γi , , i = 1, . . . , 2n. (2.5) γ0 , α 2α Then, 2n 2n X X γi Xi = x, γi (Xi − x)(Xi − x)T = P. (2.6) i=0
i=0
UKF uses the unscented transformation to approximate the error covariance and estimate the state xk . Letting xf0 denote an initial estimate of x0 with error covariance P0f , UKF is given by the following steps. UKF data assimilation step: f f xda k = xk + Kk (yk − yk ),
(2.7)
= Ck xfk , (2.8) ¤ £ da da da X · · · X , (2.9) = Ψ(xda , S , α) = k k 0,k 2n,k f f = Sk Hk , (2.10) f T f T −1 Kk = Pk Ck (Ck Pk Ck + Rk ) , (2.11) f n×n where Hk ∈ R satisfies f f T Hk (Hk ) =In−(Ck Skf )T (Ck Skf (Ck Skf )T+Rk )−1Ck Skf(2.12) ykf Xkda Skda
To reduce the ensemble size, we use a reduced-rank apf f proximation Pˆs,k of Ps,k . The reduced-rank approximations f f are chosen such that kPˆs,k − Ps,k kF is minimized subject to f ˆ rank(Ps,k ) = q, where k·kF denotes the Frobenius norm. Let P ∈ Rn×n be positive semidefinite, let σ1 > · · · > σn > 0 be the singular values of P , and let u1 , . . . , un ∈ Rn be the corresponding mutually orthogonal singular vectors. Next, define Uq ∈ Rn×q and Σq ∈ Rq×q by £ ¤ Uq , u1 · · · uq , Σq , diag(σ1 , . . . , σq ). (3.1)
With this notation, the singular value decomposition of P is given by P = Un Σn UnT , where Un ∈ Rn×n is orthogonal. For q 6 n, let ΦSVD (P, q) ∈ Rn×q denote the SVD1/2 based rank-q approximation of the square root Un Σn of 1/2 P given by ΦSVD (P, q) , Uq Σq . As noted in [18], Pˆ = SS T , where S , ΦSVD (P, q), is the best rank-q approximation of P with respect to the Frobenius norm. The SVD-based reduced-rank square-root unscented Kalman filter (SVDRRUKF) is given by the following steps. SVDRRUKF data assimilation step: f f xda s,k = xs,k + Ks,k (yk − ys,k ),
f ys,k da Xs,k da Ss,k
and Skf ∈ Rn×n satisfies
Skf (Skf )T = Pkf .
(2.13)
Finally, define Pkda ∈ Rn×n by
Pkda , Skda (Skda )T .
(2.14)
Ks,k
=
xfk+1
=
f = Pk+1 f ∆Xi,k
da f (Xi,k , uk , k), i = 0, . . . , 2n, 2n X f γi Xi,k+1 , i=0 2n X f f (∆Xi,k+1 )T + Qk , γi ∆Xi,k+1 i=0 f Xi,k
xfk .
= =
£ da ¤ da = Xs,0,k · · · Xs,2q,k ,
(3.3) (3.4)
(3.5) ¢ ¡ −1 f f f f )T + Rk , (3.6) (Ck Ss,k )T Ck Ss,k (Ck Ss,k = Ss,k
SVDRRUKF forecast step:
f da Xs,i,k+1 = f (Xs,i,k , uk , k), i = 0, . . . , 2q,
(2.15) xfs,k+1 = (2.16) f Ps,k+1 =
(2.17)
2q X
i=0 2q X
(3.7)
f , γq,i Xs,i,k+1
(3.8)
f f γq,i ∆Xs,i,k+1 (∆Xs,i,k+1 ) T + Qk ,
(3.9)
i=0
f f Ss,k+1 = ΦSVD (Ps,k+1 , q).
where , − When the dynamics in (2.1) are linear, UKF is equivalent to the Kalman filter [12]. Furthermore, in the linear case, Pkda and Pkf are the covariances of f the error xk − xda k and xk − xk , respectively. However, in f da the nonlinear case, Pk and Pk are pseudo-error covariances. The case in which the measurements depend nonlinearly on the state is discussed in [13]. Note that Hkf and Skf satisfying (2.12) and (2.13) are not unique. Moreover, all square Hkf and Skf satisfying (2.12) and (2.13) are related by an orthogonal transformation. Specifically, the following result is given in [18, p. 188]. Lemma 2.1: Let S, Sˆ ∈ Rn×n . Then, SS T = SˆSˆT if and only if there exists an orthogonal matrix U ∈ Rn×n such that Sˆ = SU .
(3.2)
f f where Hs,k ∈ Rq×q satisfies (2.12) with Skf replaced by Ss,k .
UKF forecast step: f Xi,k+1
=
Ck xfs,k , da Ψq (xda s,k , Ss,k , α) f f Ss,k Hs,k ,
f ∆Xs,i,k
f Xs,i,k
(3.10)
xfs,k .
, − f da Next, define Pˆs,k , Pˆs,k ∈ Rn×n by
where
da T da da f f f ) . (Ss,k , Ss,k )T , Pˆs,k (Ss,k , Ss,k Pˆs,k
(3.11)
It then follows from (3.5) that da f f f f Pˆs,k = Pˆs,k − Pˆs,k CkT (Ck Pˆs,k CkT + Rk )−1 Ck Pˆs,k . (3.12) Furthermore, (3.6) and (3.12) imply that f f CkT + Rk )−1 . Ks,k = Pˆs,k CkT (Ck Pˆs,k
(3.13)
Since the SVD in (3.10) is computationally intensive [9], we introduce an alternative method to obtain a reducedrank approximation of a square root of the pseudo-error covariance.
1275
CDRRUKF forecast step:
IV. C HOLESKY-D ECOMPOSITION -BASED R EDUCED -R ANK U NSCENTED K ALMAN F ILTER
f da Xc,i,k+1 = f (Xc,i,k , uk , k), i = 0, . . . , 2q
The filter gain Kk of UKF depends on a particular subspace of the forecast error covariance Pkf . Specifically, Kk depends only on the correlation Ck Pkf between the error in the measured states and unmeasured states. Since rank(Ck ) = p, there exists a state space basis with respect to which Ck has the form £ ¤ Ck = Ip 0 . (4.1) The following result is given in [8]. Lemma 4.1: Partition Pkf as · f ¸ f Pp,k (Ppp,k )T f Pk = , f f Ppp,k Pp,k
(4.2)
f f ∈ Rp×p , and assume that Ck where Pp,k ∈ Rp×p and Pp,k has the form (4.1). Then, · f ¸ Pp,k f (Pp,k + Rk )−1 . (4.3) Kk = f Ppp,k To reduce the ensemble size, we construct a filter that f f uses a reduced-rank approximation Pˆc,k of Pc,k such that f f f ˆ ˆ rank(Pc,k ) < n and kCk (Pc,k − Pc,k )kF is minimized. To f obtain Pˆc,k , we perform a Cholesky decomposition of the f pseudo-error covariance Pc,k at each time step. Assuming n×n that P ∈ R is positive definite, the Cholesky decomposition of P yields a unique lower triangular Cholesky factor L ∈ Rn×n satisfying
LLT = P.
£
Truncating the last n − q columns of L = L1 · · · yields the rank-q Cholesky factor £ ¤ ΦCHOL (P, q) , L1 · · · Lq ∈ Rn×q .
(4.4) ¤ Ln (4.5)
The following result is given in [8]. Lemma 4.2: Let P ∈ Rn×n be positive definite, define S , ΦCHOL (P, q), where 0 < q 6 n, and partition P and Pˆ , SS T as ¸ · · ¸ Pq Pqq Pˆqq Pˆq ˆ , P = P = , (4.6) (Pqq )T Pq (Pˆqq )T Pˆq
where Pq , Pˆq ∈ Rq×q and Pq , Pˆq ∈ Rq×q . Then, £ ¤ £ ¤ (4.7) Pˆq Pˆqq = Pq Pqq . Lemma 4.2 implies that, if S = ΦCHOL (P, q), then the first q columns and rows of SS T and P are equal. The Cholesky-decomposition-based reduced-rank unscented Kalman filter (CDRRUKF) is summarized as follows. CDRRUKF data assimilation step: f f xda c,k = xc,k + Kc,k (yk − yc,k ),
f yc,k da Xc,k da Sc,k
=
= =
Kc,k =
Ck xfc,k , da Ψq (xda c,k , Sc,k , α), f f Sc,k Hc,k , f f (Ck Sc,k )T Sc,k
(4.8) (4.9) (4.10)
(4.11) ¡ ¢−1 f f T , (4.12) Ck Sc,k (Ck Sc,k ) +Rk
f f where Hc,k ∈ Rq×q satisfies (2.12) with Skf replaced by Sc,k .
xfk+1 = f Pc,k+1 =
2q X
i=0 2q X
f γq,i Xc,i,k+1 ,
(4.13) (4.14)
f f γq,i ∆Xc,i,k+1 (∆Xc,i,k+1 )T + Qk , (4.15)
i=0
f f Sc,k+1 = ΦCHOL (Pc,k+1 , q), f ∆Xc,i,k
(4.16)
f = Xc,i,k − xfc,k . da ˆ f Pˆc,k , Pc,k ∈ Rn×n
where f f f Next, define by Pˆc,k , Sc,k (Sc,k )T and da f f f f Pˆc,k , Pˆc,k − Pˆc,k CkT (Ck Pˆc,k CkT + Rk )−1 Ck Pˆc,k . (4.17)
da da T da It then follows from (4.11) that Sc,k (Sc,k ) = Pˆc,k . Furthermore, (4.12) and (4.17) imply that Kc,k = Pˆ f C T (Ck Pˆ f C T + Rk )−1 . (4.18) c,k
k
c,k
k
Hence, like the estimator gain Ks,k of SVDRRUKF, the estimator gain Kc,k of CDRRUKF depends on a reducedf rank approximation Pˆc,k of the pseudo-error covariance f Pc,k . Due to the rank-reduction step (4.16), CDRRUKF is generally not equivalent to UKF. However, we now discuss cases in which the performance of CDRRUKF is close to that of UKF. A. Basis Selection for CDRRUKF The following result given in [8] shows that CDRRUKF is equivalent to UKF for a single time step when Ck has the form (4.1). Proposition 4.1: Assume that p 6 n, q = p, and Ck has the structure in (4.1). Let Pkf ∈ Rn×n be positive semidefinite and Kk be given by (2.11). Furthermore, define f f Sc,k , ΦCHOL (Pkf , q) and let Pˆc,k be given by (4.17). Then, f f ˆ Ck Pc,k = Ck Pk and hence, Kc,k = Kk . If the dynamics (2.1) are linear and time-invariant, that is, for all k > 0, xk+1
=
Axk + Buk + wk ,
(4.19)
yk
=
Cxk + vk ,
(4.20)
then a basis for the state x can be chosen so that CDRRUKF is equivalent to UKF for r > 0 time steps. To construct such a basis, we define the observability matrix O(A, C) ∈ Rpn×n by ¤T △ £ O(A, C) = C T (CA)T · · · (CAn−1 )T . (4.21) The following result is given in [8]. Proposition 4.2: Assume that O(A, C) has the form ¸ · In . (4.22) O(A, C) = 0(p−1)n×n Let r > 0 be an integer such that pr < n, and let q = f pr. Furthermore, assume that Pc,0 = P0f . Then, for all k = 0, . . . , r, Kc,k = Kk . If, in addition, xfc,0 = xf0 , then for all k = 0, . . . , r, xfc,k = xfk . £ ¤T x1,k · · · xn,k Let xk have entries xk = . If O(A, C) has the form (4.22), and uk = wk = vk = 0,
1276
for all k > 0, then for all k > 0, and for every integer r > 0 such that pr 6 n, £ T ¤T £ T ¤T T yk · · · yk+r−1 = x1,k · · · xT . (4.23) pr,k
Next, we consider the nonlinear system (2.1) and assume that the dynamics in (2.1) can be expressed as xi,k+1 = fi (xφL (i,b) , . . . , xφR (i,b) , uk , k) + wi,k ,
(4.24)
where i = 1, . . . , n, b > 0 and φL (i, b) , max(1, i − b), φR (i, b) , min(n, i + b). (4.25) For example, in systems modeled by finite volume schemes, the next value of a physical variable in a given cell depends only on the present values of the physical variables in its neighboring cells. Next, let yk denote a measurement of a specific component of the state, so that yk = xj,k +vk , where j ∈ {1, . . . , n}, j− rb > 1, and j + rb 6 n. It follows from (4.24) that, if wk = vk = 0, for all k > 0, then yk , . . . , yk+r−1 depends on only first 2rb components of the state vector x ˜k at time step k. Hence, a state space basis can be chosen such that the outputs yk , . . . , yk+r−1 depend on only the first few components of the state vector. V. L INEAR A DVECTION M ODEL Consider a linear advection model [2] with n cells, and let xi,k be the energy in the ith cell at time k. The energy flow satisfies ( xi−1,k , if i = 2, . . . , n, (5.1) xi,k+1 = xn,k , if i = 1. Hence, energy in the ith cell flows to the (i + 1)th cell, while the periodic boundary condition ensures that the energy circulates continually. We choose n = 100 and assume that the disturbance wk enters selected cells, where wk ∈ Rn is a white noise process with covariance Qk = Q for all k > 0, and Q ∈ Rn×n is diagonal with nonzero entries Qi = 1 only for i ∈ {10, 20, . . . , 100}. Next, we assume that measurements of the energy in cells 50 and 51 are available so that yk = [x50,k x51,k ]T + vk , where vk is white noise process with covariance Rk = 0.1I2 . First, we use the measurements yk to estimate the energy in the remaining cells using UKF. In all three cases, the initial estimates xf0 , xfs,0 , and xfc,0 are not equal to the initial f f state x0 . Moreover, we choose P0f = Ps,0 = Pc,0 = 0.1In . Finally, we choose α = 0.6 for all three filters. As shown in Figure 1 and Figure 2, data assimilation is performed using SVDRRUKF and CDRRUKF for several values of q between 5 and 100. It can be seen that the performance of SVDRRUKF with 111 ensemble members (q = 55) is close to optimal, whereas the performance of CDRRUKF is close to optimal with 11 ensemble members (q = 5). The steady-state MSE of state estimates for various values of q is plotted in Figure 3 and Figure 4. The performance of SVDRRUKF is poor when q < 55, and close to optimal when q > 55. Thus the ensemble size can be reduced from 201 to 111 with negligible change in the
performance. However, the ensemble size can be reduced from 211 to 11 with negligible performance deterioration. Next, we repeat the same procedure except with a poor estimate of the process noise covariance for data assimilation. ˆ k , where Specifically, we replace Qk in (3.9) and (4.15) by Q ˆ Qk = I for all k > 0. The steady-state MSE of state estimates for different choices of q is plotted in Figure 3 and Figure 4. SVDRRUKF with a poor estimate of the error covariance is unstable for all q 6 95 (indicated by the X’s). However, Figure 4 shows that even with q = 5 and a poor estimate of the process noise covariance, the performance of CDRRUKF is close to optimal. ˆ k , where Q ˆ k = αI Finally, we replace Qk in (4.16) with Q for all k > 0, and perform state estimation using CDRRUKF. The steady-state MSE of the state estimates is shown in Figure 5 for several values of α. For all the cases, the performance of CDRRUKF is close to optimal when α > 1. This suggests that it is advantageous to overestimate the process noise covariance. SVDRRUKF is unstable for all choices of α = 0.005, . . . , 50. Hence, these simulations suggest that CDRRUKF is more robust than SVDRRUKF with respect to uncertainties in the process noise covariance. All of the results for CDRRUKF in figures 1-6 are obtained using a state space basis with respect to which the observability matrix has the form (4.22). VI. L96 M ODEL The L96 model mimics the propagation of an unspecified meteorological quantity along a latitude circle [17]. The dynamics are governed by d xi (t) = (xi+1 (t) − xi−2 (t))xi−1 (t) − xi (t) + ui (t), (6.1) dt where xi (t) ∈ R denotes the meteorological quantity at the ith grid point at time t, ui ∈ R denotes an external forcing term, and wi denotes unknown disturbances affecting the ith grid point. For all t > 0, the boundary conditions are defined by x0 (t) = xn (t), x−1 (t) = xn−1 (t), xn+1 (t) = x1 (t). We choose ui (t) = 8 for all i = 1, . . . , n and all t > 0. Using fourth-order Runge-Kutta discretization with a sampling time of 0.05 s, we obtain a discrete-time model of (6.1) that can be expressed as (2.1). Furthermore, we assume that the discretized model is corrupted by an unknown external disturbance that affects specified cells. We choose n = 40, and assume that wk is white noise process with covariance Qk = Q for all k > 0, where Q ∈ Rn×n is diagonal with nonzero entries Qi,i = 0.1 only for i ∈ {5, 15, 25, 35}. Next, we assume that measurements from cells with 20 and 21 are T available so that yk = [x20,k x21,k ] + vk , where vk is white noise process with covariance Rk = 0.01I2 . We use the measurements yk to estimate the state in the cells where measurements are not available. Next, as shown in Figure 6 and Figure 7, we reduce the ensemble size and use SVDRRUKF and CDRRUKF with q = 10, 20, 30. It can be seen that the performance of SVDRRUKF is poor compared to the performance of CDRRUKF for both q = 20 and q = 30. Moreover, the performance of CDRRUKF with
1277
3.6
log10(MSE of the state estimates)
3.2
3.6
SVDRRUKF q=5 SVDRRUKF q=15 SVDRRUKF q=25 SVDRRUKF q=35 SVDRRUKF q=45 SVDRRUKF q=55 UKF
3.2
3
2.8
2.6
2.4
2.2
3
2.8
2.6
2.4
2.2
2
2
1.8
1.8
1.6
50
100
CDRRUKF q=5 UKF
3.4
log10(MSE of the state estimates)
3.4
150 200 Time step k
250
300
1.6
350
Fig. 1. MSE of the state estimates obtained from SVDRRUKF for several values of q. SVDRRUKF with q = 5 is unstable, while the performance of SVDRRUKF with q = 55 is close to the performance of full-order UKF.
250
300
350
SVDRRUKF SVDRRUKF with a poor estimate of Qk
Steady−state MSE of the state estimates
2500
2000
1500
1000
UKF 500
0
X
X X 5
15
X
X
X
X
X
X
X
25 35 45 55 65 75 85 q : Rank of the approximation of the error covariance
X 95
100
Fig. 3. Steady-state performance of SVDRRUKF. For each value of q, we perform data assimilation with the true value of the process noise covariance and with a poor estimate of the process noise covariance. The X’s indicate cases in which the filter is unstable.
3000 CDRRUKF CDRRUKF with a poor estimate of Qk 2500
2000
1500
1000
UKF 500
0
R EFERENCES [1] J. Lewis, S. Lakshmivarahan, and S. Dhall, Dynamic Data Assimilation: A Least Squares Approach, Cambridge, 2006. [2] G. Evensen, Data Assimilation: The Ensemble Kalman Filter, Springer, 2006. [3] S. F. Farrell and P. J. Ioannou, “State Estimation Using a ReducedOrder Kalman Filter,” J. Atmos. Sci., vol. 58, pp. 3666–3680, 2001.
150 200 Time step k
3000
Steady−state MSE of the state estimates
We developed a reduced-rank square-root unscented Kalman filter based on the Cholesky decomposition of the pseudo-error covariance. We compare the performance of the Cholesky-based filter with an analogous filter that uses the singular value decomposition for a linear advection model and a nonlinear system that exhibits chaotic behavior. Although the computational requirement of the Choleskybased reduced-rank filter is less than that of the SVDbased reduced-rank filter, the results presented here suggest that the estimation accuracy of the Cholesky-based reducedrank filter is significantly better than that of SVD-based reduced-rank filter. Moreover, simulation results suggest that the Cholesky-based reduced-rank filter is more robust to uncertainties in the process noise covariance than the SVDbased reduced-rank filter.
100
Fig. 2. MSE of the state estimates obtained from CDRRUKF with q = 5. The performance of CDRRUKF with q = 5 is close to the full-order UKF performance.
61 (q = 30) ensemble members is close to the performance of UKF with 81 ensemble members. ˆ k , where Next, we replace Qk in (3.9) and (4.15) by Q ˆ k = αI for all k > 0. Figure 8 shows the time-averaged Q MSE of state estimates obtained using SVDRRUKF and CDRRUKF with q = 10 and q = 20 for several values of α between 0.001 and 100. It can be seen that, for all values of α, the performance of CDRRUKF is superior to the performance of SVDRRUKF. In fact, CDRRUKF with 21 ensemble members (q = 10) consistently outperforms SVDRRUKF with 41 ensemble members (q = 20). Finally, note that (6.1) can be expressed as (4.24) with b = 1. All of the results for CDRRUKF in figures 8-12 are obtained using the state space basis discussed in Section IV-A. VII. C ONCLUSIONS
50
5
15
25 35 45 55 65 75 85 q: Rank of the approximation of the error covariance
95
100
Fig. 4. Steady-state performance of CDRRUKF for values of q between 5 and 100. Note that for q = 5, the performance of CDRRUKF is close to optimal, irrespective of the value of the process noise covariance used for data assimilation.
1278
3.5 CDRRUKF q=5 CDRRUKF q=15 CDRRUKF q=25
4.5
No data assimilation CDRRUKF q=20 CDRRUKF q=30 UKF
3.3
log10(MSE of the state estimates)
log10(Steady−state MSE of the state estimates)
3.4
3.2
3.1
3
2.9
2.8
2.7
4
3.5
3
2.5
2
2.6 −2.5
−2
−1.5
−1
−0.5 0 log (α)
0.5
1
1.5
2
1.5
10
Fig. 5. Steady-state performance of CDRRUKF with q = 5, 15, 25. In all three cases, we use a poor estimate of the process noise covariance for data assimilation. For a fixed level of uncertainty in the process noise covariance, the performance of CDRRUKF improves when the ensemble size increases. The performance of SVDRRUKF is not shown since SVDRRUKF is unstable for all values of α and q = 5, 15, 25.
5
10
15
20
25 Time (sec)
30
35
40
45
50
Fig. 7. Performance of CDRRUKF with n = 40 and q = 20, 30. Note that the performance of CDRRUKF with q = 20 is better than the performance of SVDRRUKF with q = 30.
4 No data assimilation SVDRRUKF q=10 SVDRRUKF q=20 CDRRUKF q=10 CDRRUKF q=20 UKF
3.8 No data assimilation SVDRRUKF q=20 SVDRRUKF q=30 UKF
log10(Time−averaged MSE of the state estimates)
log10(MSE of the state estimates)
4.5
4
3.5
3
2.5
3.6
3.4
3.2
3
2.8
2.6
2.4
2.2 2
1.5
2 −3
5
10
15
20
25 Time (sec)
30
35
40
45
−2.5
−2
−1.5
−1
−0.5 log10(α)
0
0.5
1
1.5
2
50
Fig. 6. MSE of the state estimates obtained using SVDRRUKF with q = 20, 30. The error in state estimates based on UKF and data-free simulation is shown for comparison. SVDRRUKF with q = 20 and q = 30 sometimes yields estimates that are worse than estimates based on data-free simulation.
[4] L. Scherliess, R. W. Schunk, J. J. Sojka, and D. C. Thompson, “Development of a Physics-Based Reduced State Kalman Filter for the Ionosphere,” Radio Science, vol. 39, RS1S04, pp. 1-12, 2004. [5] I. S. Kim, J. Chandrasekar, H. J. Palanthandalam-Madapusi, A. Ridley, and D. S. Bernstein, “State Estimation for Large-Scale Systems Based on Reduced-Order Error-Covariance Propagation,” Proc. Amer. Contr. Conf., New York, pp. 5700-5705, June 2007. [6] M. Verlaan and A. W. Heemink, “Tidal Flow Forecasting Using Reduced-Rank Square-Root Filters,” Stochastic Hydrology and Hydraulics, vol. 11, pp. 349–368, 1997. [7] S. Gillijns, D. S. Bernstein, and B. D. Moor, “The Reduced Rank Transform Square Root Filter for Data Assimilation,” Proc. 14th IFAC Symposium on System Identification (SYSID2006), Newcastle, Australia, pp. 349–368, 2006. [8] J. Chandrasekar, I. S. Kim, and D. S. Bernstein, “Cholesky-Based Reduced-Rank Kalman Filtering,” submitted to Amer. Contr. Conf, 2008. [9] G. W. Stewart, Matrix Algorithms Volume 1: Basic Decompositions, SIAM, 1998. [10] G. Evensen, “The Ensemble Kalman Filter: Theoretical Formulation and Practical Implementation,” Ocean Dynamics, vol. 53, pp. 343-36,
Fig. 8. Time-averaged MSE of state estimates between 35 sec and 50 sec. The state estimates are obtained using SVDRRUKF and CDRRUKF with q = 10 and q = 20, and a poor estimate of the process noise covariance. For all values of α, the performance of CDRRUKF is better than the performance of SVDRRUKF.
2003. [11] J. S. Whitaker and T. M. Hamill, “Ensemble Data Assimilation without Perturbed Observations,” Mon. Wea. Rev., vol. 130, pp. 1913-1924, 2002. [12] S. Julier, J. Uhlmann, and H. F. Durrant-Whyte, “A New Method for the Nonlinear Transformation of Means and Covariances in Filters and Estimators,” IEEE Trans. Autom. Contr., vol. 45, pp. 477-482, 2000. [13] B. Ristic, S. Arulampalam and N. Gordon, Beyond the Kalman Filter, Artech House, 2004. [14] M. K. Tippett, J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker, “Ensemble Square Root Filters,” Mon. Wea. Rev., vol. 131, pp. 1485–1490, 2003. [15] D. S. Bernstein and D. C. Hyland, “Compartmental Modeling and Second-Moment Analysis of State Space Systems,” SIAM J. Matrix Anal. Appl., vol. 14, no. 3, pp. 880–901, 1993. [16] B. D. O. Anderson and J. B. Moore, Optimal Filtering, Prentice-Hall, 1979, reprinted by Dover, 2005. [17] E. N. Lorenz, “Predictability - A Problem Partly Solved,” Predictability of Weather and Climate, Cambridge University Press, 2006. [18] D. S. Bernstein, Matrix Mathematics, Princeton University Press, 2005.
1279
