Numerical Results of Nonlinear Filtering Problem ... - Semantic Scholar
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
971
Numerical Results of Nonlinear Filtering Problem from Yau-Yau Method Zhen Liu Department of Mathematics, Zhejiang University of Technology， Hangzhou 310023，China Email: [email protected]
Fangfang Dong School of Statistics and Mathematics ,Zhejiang Gongshang University， Hangzhou 310018， China Email: [email protected]
Luwei Ding School of Science Huazhong Agricultural University， Wuhan 430070， China Email: [email protected] Abstract— In the paper, we introduce a kind of method to solve the nonlinear filtering problem. Firstly, we review the basic filtering problem and the reduction from robust Duncan-Mortensen-Zakai equation to Kolmogorov equation. Then we use the difference discrete method to solve the Kolmogorov equation. The result is given to prove that the solution of the difference scheme convergences pointwise to the solution of the initial-value problem of the Kolmogorov equation. At last, the numerical results show that the numerical method can give the exact result. Index Terms—Nonlinear filtering equation, DMZ equation, Kolmogorov equation
I. INTRODUCTION One of the key goals of data assimilation is to obtain, recursively in time, good estimates of the state of a stochastic dynamical system based on noisy partial observations of the same. The major breakthrough in this classical problem of signal analysis was the landmark work of Kalman and Bucy [1]. The power of KalmanBucy filtering, which can be gauged from the broad spectrum of disciplines where it is used, lies in its elegant simplicity. Since then, the Kalman-Bucy filtering has proved useful in many areas such as navigational and guidance systems, radar tracking, solar mapping, and satellite orbit determination. Ever since the technique of the Kalman-Bucy filtering was popularized, there has been an intense interest in developing nonlinear filtering theory. In the 1960s and early 1970s, the basic approach to nonlinear filtering theory was via innovation methods originally proposed by Kailath [2] and Frost and Kailath[3] and subsequently rigorous developed by Fujisaki, Kallianpur, and Kunita [4] in 1972. But the difficulty with this approach is that innovation process is not, in general, explicitly computable. In the late 1970s, Brockett and Clark [5], Brockett [6] and Mitter [7]
© 2012 ACADEMY PUBLISHER doi:10.4304/jcp.7.4.971-976
proposed the idea of using estimation algebras to construct a finite-dimensional nonlinear filtering. Yau and his coworkers have finally classified all finite dimensional estimation algebras of Maximal rank ([8][19]). Although it is an interesting and challenging problem to classify all finite dimensional filtering, it appears that from Yau’s previous works, finite dimensional filtering simply does not exist for many practical situations. In ([20-21]), Yau and Yau proved the existence and decay estimates of the solution to the DMZ(Duncan-MortensenZakai) equation under the assumption that the drift terms of the signal dynamic and observation dynamic respectively have linear growths. Later they [22] showed that the real-time solution of DMZ equation can be reduced to off-time solution of Kolmogorov equation if the drift terms have linear growths. Recently Yau and Yau ([23-24]) finally proved that the real-time solution of DMZ equation can be reduced to off-time solution of Kolmororov equation if the growth of the drift term of observation dynamic at infinity is faster than the growth of the drift term of signal dynamic at infinity. In section 2, we recall the basic filtering problem. In section 3, we review the reduction from robust DMZ equation to Kolmogorov equation of Yau-Yau method. In section 4, we shall use finite difference method to derive the solution of Kolmororov equation and the numerical simulation of Yau-Yau method is shown. II. INTRODUCTION OF NONLINEAR FILTERING PROBLEM The filtering problem considered here is based on the signal observation model dx(t ) = f ( x(t ))dt + g ( x(t ))dv(t ) , x(0) = x0 , (1) dy (t ) = h( x(t ))dt + dw(t ) ，
y (0) = 0 , n
p
m
(2) m
in which x, v, y and w are R , R , R and R valued processes and v, w have components that are independent,
972
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
standard Brownian processes. We further assume that n = p . f , g and h are vector-valued, matrix-valued and vector-valued C ∞ smooth functions. We shall refer to x(t ) as the state of the system and y (t ) as the observation at time t . Let ρ ( x, t ) denote the conditional probability density of the state given the observation { y ( s ) : 0 ≤ s ≤ t} . It is well known that ρ ( x, t ) is given by normalizing a function σ (t , x) that satisfies the following DuncanMortensen-Zakai equation m
dσ (t , x) = L0σ (t , x)dt + ∑ Liσ (t , x) dyi (t ), i =1
σ (0, x) = σ 0 ( x) , Where n n ∂f 1 m 1 n ∂2 ∂ L0 = ∑ 2 − ∑ f i − ∑ i − ∑ hi2 . 2 i =1 ∂xi i =1 ∂xi i =1 ∂xi 2 i =1
(3)
(4)
And for i = 1, 2, ⋅⋅⋅, m, Li is the zero degree differential operator of multiplication by hi , σ 0 is the probability density of the initial point. Equation (3) is a stochastic partial differential equation. In real application, we are interested in constructing robust state estimators from observed sample paths with some property of robustness. From the robust algorithms proposed by Davis [25], it define a new unormalized density m
u (t , x) = exp(−∑ hi ( x) yi (t ))σ (t , x).
(5)
i =1
Then we can reduce the equation (3) to the following equation, which is called robust DMZ equation, 1 ut = Δu + F (t , x) ⋅∇u + V (t , x)u , u (0, x) = σ 0 ( x) .(6) 2 where drift function F (t , x) = − f ( x) + ∇k (t , x) , source function 1 1 2 V (t , x) = −divf ( x) − h( x) + Δk (t , x) 2 2 1 2 − f ( x ) ⋅ ∇k (t , x) + ∇k (t , x) . 2
uk (τ k , x ) uniformly in x
(10)
Where C is the constant and Tk = sup1≤i ≤ n { τ i − τ i −1 } .
So the solution of (7) approximate the solution of robust D-M-Z equation very well in the pointwise sense. Therefore it remains to describe an algorithm to compute u k (τ k , x), which is based on the following property. Proposition 1. [23-24] u (t , x) satisfies the following Kolmogorov equation n ∂u 1 ∂u (t , x) = Δu (t , x) − ∑ fi ( x) (t , x) ∂t 2 ∂xi i =1 (11) n ∂f 1 m − (∑ i ( x) + ∑ hi2 ( x))u (t , x), 2 i =1 i =1 ∂xi for τ l −1 ≤ t ≤ τ l , if and only if m
u (t , x) = exp(−∑ hi ( x) yi (τ l ))u (t , x) i =1
satisfies the robust DMZ equation with observation being fixed at y (τ l ) ∂u 1 (t , x) = Δu (t , x) + F (τ l , x) ⋅∇u + V (τ l , x)u. (12) ∂t 2 From the property, u1 (τ 1 , x) can be computed by u1 (τ 1 , x) where u1 (t , x) satisfies the following Kolmogorov equation n ∂u1 ∂u 1 (t , x) = Δu1 (t , x) − ∑ f i ( x) 1 (t , x) 2 ∂t ∂xi i =1 , n m ∂fi 1 2 ( x) + ∑ hi ( x))u1 (t , x), − (∑ 2 i =1 i =1 ∂xi m
u1 (0, x) = σ 0 ( x) exp(∑ hi ( x) y i (τ 1 )) .
(13)
m
The fundamental problem of nonlinear filtering theory is how to solving the robust DMZ equation(6) in real time and in memoryless manner. In the section, we introduce the Yau-Yau algorithm which achieves this goal for a large class of filtering system with arbitrary initial distribution by reducing it to solve Kolmogorov equation. Suppose that u (t , x) is the solution of robust D-M-Z equation and we want to compute u (τ , x) .Let a
partition
of
[0,τ ] .And let ui (t , x) be a solution of the following
partial differential equation for τ i −1 ≤ t ≤ τ i
© 2012 ACADEMY PUBLISHER
k →0
In fact,
EQUATION
be
and u (τ , x) = lim T
i =1
III. YAU-YAU METHODS FOR THE ROBUST DMZ
Tk = {0 = τ 0 < τ 1 < ⋅⋅⋅ < τ k = τ }
∂ui 1 (t , x) = Δu (t , x) + F (τ i , x) ⋅ ∇u + V (τ i , x)u , (7) 2 ∂t ui (τ i −1 , x) = ui −1 (τ i −1 , x) . (8) Theorem 1. [23-24] Let u (t , x) and uk (t , x) be the solutions of (6) and (7) respectively. For any given ε > 0 and sufficiently large n , u (t , x) − uk (t , x) ≤ ε TeCT , (9)
u1 (τ 1 , x) = exp(−∑ hi ( x) y i (τ 1 ))u1 (τ 1 , x).
(14)
i =1
For i ≥ 2 , ui (τ i , x) can be computed by ui (τ i , x) where ui (t , x ) for τ i −1 ≤ t ≤ τ i satisfies the following Kolmogorov equation n ∂ui ∂u 1 (t , x) = Δui (t , x) − ∑ f i ( x) i (t , x) ∂t ∂xi 2 i =1 n
− (∑ i =1
m
∂fi 1 m ( x) + ∑ hi2 ( x))ui (t , x), ∂xi 2 i =1
ui (τ i −1 , x) = exp(∑ hi ( x)( y i (τ i )) − yi (τ i −1 ))ui −1 (τ i −1 , x). (15) i =1
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
973
Observe that in the algorithm at step i ,we only need the observation at time τ i and τ i −1 . We do not need any other previous observation data. Observe also that the Kolmogorov equation (15) is uniform for all time steps and it depends on observation y (t ) only initial condition.
convergences pointwise to the solution of the initialvalue problem 1 ∂ 2u ∂u ∂u (24) (t , x) = (t , x) + p ( x) (t , x) + q ( x)u (t , x). 2 2 ∂x ∂t ∂x Proof. We denote the exact solution to the initial-value problem by u = u (t , x) and set
IV. NUMERICAL METHOD FOR THE KOLMOGOROV
zkn = u (k Δx, nΔt ) − vkn . (25) If we insert u into equation (20), and multiply through by Δt , we see that ukn = u (k Δx, nΔt ) satisfies
EQUATION
Consider the Kolmogorov equation (11), we use finite difference method to approximate the solution of the equation. We take the Kolmogorov equation ∂u 1 ∂ 2u ∂u (t , x) = (t , x) + p( x) (t , x) + q ( x)u (t , x) , (16) ∂t 2 ∂x 2 ∂x as an example to explain the method, where ∂f 1 m p( x) = − f ( x), q( x) = −( ( x) + ∑ hi2 ( x)) . ∂x 2 i =1 To approximate the equation (16) by finite difference, we discretize the spatial and time domain by placing a grid over the domain. For convenience, we use a uniform grid, with grid spacing Δx and Δt .The space-time domain of our problem then is approximated by the lattice of points. We will attempt to approximate the solution to our problem at the points on the lattice. Notationally, we will define vkn to be a function at the point (k Δx, nΔt ) .To arrive at the finite difference equation used to approximation the Kolmogorov equation, we use the approximation v n +1 − vkn , (17) ut (nΔt , k Δx ) ≈ k Δt v n − 2vkn + vkn−1 u xx (nΔt , k Δx) ≈ k +1 , (18) Δx 2 v n − vkn−1 . (19) u x ( n Δt , k Δ x ) ≈ k Δx We see that 1 u t − u xx − p( x)u x − q( x)u 2 n +1 u − ukn 1 ukn+1 − 2ukn + ukn−1 = k − 2 Δt Δx 2 (20) u n − ukn−1 − p ( xk ) k − q( xk )ukn , Δx where the above expression assumes that the higher order derivatives of u at (k Δx, nΔt ) are bounded. Then we have the following theorem. Theorem 2. Suppose q( x k ) < 0 ,and 1 2Δt Δt + p ( xk ) + Δtq( xk ) ≥ 0, 2 Δx 2 Δx 1 Δt Δt − p( xk ) ≥ 0, 2 Δx 2 Δx then the solution of the difference scheme 1 Δt n vkn +1 = vkn + (vk +1 − 2vkn + vkn−1 ) 2 Δx 2 Δt + p ( xk )(vkn − vkn−1 ) + Δtq( xk )v kn , Δx 1−
© 2012 ACADEMY PUBLISHER
(21) (22)
(23)
ukn +1 = ukn +
1 Δt n (uk +1 − 2ukn + ukn−1 ) 2 Δx 2
Δt (26) p ( xk )(ukn − ukn−1 ) + Δtq ( xk )ukn Δx + O(Δt 2 + ΔxΔt ). Then by subtracting equation (26) from equation (23), we see that zkn satisfies +
1 2Δt Δt + p( xk ) + Δtq ( xk )) zkn 2 Δx 2 Δx 1 Δt Δt (27) +( − p( xk )) zkn−1 2 Δx 2 Δx 1 Δt n + zk +1 + O(Δt 2 + ΔxΔt ). 2 Δx 2 In order to obtain bounds for zkn , we need to ensure that the coefficients of the three terms on the right of this equation are all nonnegative and have a sum no greater than unity. We are therefore led to the condition 1 2Δt Δt + (28) 1− p ( xk ) + Δtq( xk ) ≥ 0, 2 Δx 2 Δx 1 Δt Δ t − (29) p ( xk ) ≥ 0, 2 Δx 2 Δx as well as q( xk ) < 0 . Then we have 1 2Δt Δt zkn +1 ≤ (1 − + p( xk ) + Δtq ( xk )) zkn 2 Δx 2 Δx 1 Δt Δt +( − p ( xk )) zkn−1 (30) 2 Δx 2 Δx 1 Δt n 2 + zk +1 + A(Δt + ΔxΔt ) 2 Δx 2 ≤ Z n + A(Δt 2 + ΔxΔt ), where A is the constant associated with the " O " terms and depends on the assumed bound of the higher order zkn +1 = (1 −
{ }
derivative of u , and Z n = sup k zkn . Then taking the sup over k on the left side of the above equation yields Z n +1 ≤ Z n + A(Δt 2 + ΔxΔt ) . (31) Applying the above equation repeatedly yields Z n +1 ≤ Z n + A(Δt 2 + ΔxΔt ) ≤ " (32) ≤ Z 0 + (n + 1) A(Δt 2 + ΔxΔt )
974
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
f 2 ( x, y ) = x 2 + y 2 , h1 ( x, y ) = x − y, h2 ( x, y ) = x + y,
Since Z 0 = 0 , u (k Δx, (n + 1)Δt ) − vkn +1 ≤ z n +1 , and (n + 1)Δt → t , then we have u (k Δx, (n + 1)Δt ) − vkn +1 ≤ (n + 1) A(Δt 2 + ΔxΔt ) → 0 , (33)
as Δt , Δx → 0 . ,
A. Numerical simulation Let f ( x) = cos( x), h( x) = sin( x), and f ( x) = x 2 , h( x) = sin( x) respectively, we can have figure 1 and 2,
and f1 ( x, y ) = − cos( x) − sin( y ), f 2 ( x, y ) = cos( x ) + sin( y ), h1 ( x, y ) = x − y, h2 ( x, y ) = x + y ,
σ 0 = exp(−5( x 2 + y 2 )),
where σ 0 = exp(−10 ∗ x ). The dashed line denotes the mean of the conditional probability density ρ and the solid line is generated from discrete extended Kalman filter equations. Figure 3 and 4 show the numerical behavior of the conditional probability density ρ in different time. Let f1 ( x, y ) = − x 2 − y 2 ,
we have figure (5,6), in which figure(5,6) are the means about the x, y respectively. The dashed line denotes the mean of the conditional probability density ρ and the solid line is generated from discrete extended Kalman filter equations. Figure (7,8) show the numerical behavior of the conditional probability density ρ in different time.
Fig. 1
Fig. 2
2
Fig. 3
© 2012 ACADEMY PUBLISHER
Fig. 4
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
975
Fig.5
Fig.6
Fig.7
Fig.8
National Natural Science Foundation of China (Grant No.10801045). V. CONCLUSION In the paper, we use the Yau-Yau method to solve the nonlinear filtering problem, which reduce the robust DMZ equation to Kolmogorov equation. We give a finite difference schemes to derive the numerical solution of the Kolmogorov equation. A proof is given to prove that the solution of the difference scheme convergences pointwise to the solution of the initial-value problem of the Kolmogorov equation. At last, we give some numerical experiments to test the method. ACKNOWLEDGMENT The authors are grateful to the anonymous referees for their valuable remarks and helpful suggestions, which have significantly improved the paper. This work is supported by Zhejiang Provincial Natural Science Foundation of China (Grant No. Y6100526 ) and by
© 2012 ACADEMY PUBLISHER
REFERENCES [1] R.E. Kalman and R.S. Bucy, New results in linear filtering and prediction theory,Trans. ASME Series D, J. Basic Engineering 83(1961), 95-108. [2] T. Kailath, An innovations approach to least-squares estimation, Part I: Linear filtering in additive white noise, IEEE Transactions on Automatic Control,13(1986), pp.646655. [3] P.A. Frost and T. Kailath, An innovations approach to leastsquares estimation II, IEEE Transactions on Automatic Control, 16(1971), pp.217-226. [4] M.Fujisaki, G. Kallianpur, and H. Kunita, stochastic differential equations for the nonlinear filtering problem, Osaka Journal of Mathematics, 2(1972), PP.19-40. [5] R.W. Brockett and J.M.C. Clark, The geometry of the conditional density functions, in: Analysis and Optimization of Stochastic Systems, O.L.R. Jacobs eral (Eds.), Academic Press, New York, (1980) 299-309. [6] R.W. Brockett, Nonlinear systems and nonlinear estimation theory, in: The mathematics of Filtering and Identification
976
and Application, M.Hazewinkel and J.C. Williams (Eds.), Reidel, Dordrecht, (1981), 479-504. [7] S.K. Mittar, On the analogy between mathematical problems of nonlinear filtering and quantum physics, Ricerche Automat., 10(1979), pp. 163-216. [8] S. S.-T. Yau, Recent results on nonlinear filtering: New class of finite dimensional filters, in Proceedings of the 29th Conf. Decision and Control, Honolulu, 1990, pp. 231-233. [9] S. S.-T. Yau, Finite dimensional filters with nonlinear drift I: A class of filters including both Kalman-Bucy filters and Benes filters, J. Math. Syst., Estimation and Control, 4(2) (1994), pp. 181-203. [10]Yau, S.-T., and Hu, G.-Q. New direct method for Yau filtering system with arbitrary conditions. In Proceedings of 35th IEEE Conference on Design and Control, Kobe, Japan, (1996) ,pp. 2539-254 [11] Chen, J., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift VI: Linear structure of - matrix. Mathematics of Control, Signals and Systems, 9 (1996), pp. 370-385. [12] S.-T. Yau and G.-Q. Hu, Direct method without Riccati equation for Kalman-Bucy filtering system with arbitrary initial conditions, in Proceedings of the 13th World Congress IFAC, vol. II, San Francisco, CA,(1996), pp. 469474. [13] Chen, J., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift VII: Mitter conjecture and structure of η. SIAM J. Control and Optimization, 35, 4 (1997), pp. 11161131. [14] Chen, J., Yau, S. S.-T., and Leung, C.-W. Finite dimensional filters with nonlinear drift VIII: Classification of finite dimensional estimation algebra of maximal rank with state space dimension 4. SIAM J. Control and Optimization, 35, 4 (1997), pp. 1132-1141. [15] Rasoulian, A., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift IX: Construction of finite dimensional estimation algebra of non-maximal rank. Systems Control Letters, 30 (1997), pp. 109-118. [16] S. S.-T. Yau, X. Wu and W. S. Wong, Hessian Matrix nondecomposition theorem, Math. Res. Lett., 6 (1999), pp. 663-673. [17] S. S.-T. Yau, Brockett’s problem on nonlinear filtering theory, in Lectures on Systems, Control and Information, AMS/IP, Stud. Adv. Math., 17 (2000), pp. 177-212. [18] Hu, G.-Q., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift X: Explicit solution of DMZ equation. IEEE Transactions on Automatic Control, 46, 1 (2001), pp.142-148. [19]S.S.-T. Yau and G.-Q. Hu, Classification of finite dimensional estimation algebras of maximal rank with arbitrary state-space dimension and Mitter Conjecture, International Journal of Control, 78:10(2005), pp.689-705. [20]Yau, S. S.-T., and Yau, S.-T. Finite dimensional filters with nonlinear drift III: Duncan-Mortensen-Zakai equation with arbitrary initial condition for linear filtering system and the
© 2012 ACADEMY PUBLISHER
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
Benes filtering system. IEEE Transactions on Aerospace and Electronic Systems, 33, 4 (1997), pp.1277-1294. [21] Yau, S.-T., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift XI: Explicit solution of the generalized Kolmogorov equation in the Brockett-Mitter program. Advances in Mathematics, 140 (1998), pp.156-189. [22] S.-T. Yau and S.S.-T Yau, Existence and uniqueness and decay estimates for the time dependent parabolic equation with application to Duncan-Mortensen-Zakai equation, Asian Journal of Mathematics, 2(1998), pp.1079-1149. [23] S.-T. Yau and S.S.-T Yau, Real time solution of nonlinear filtering problem without memory I, mathematical Research Letter, 7(2000), pp.671-693. [24] S.-T. Yau and S.S.-T Yau, Real time solution of nonlinear filtering problem without memory II, SIAM. J. Control. Optim. 47(1)(2008). pp. 163-195. [25]Davis, M.H.A., On a multiplicative functional transformation arising in nonlinear filtering theory, Z.Wahrsch Verw. Gebiete, 54(1980), pp.125-139.
Zhen Liu received his Ph.D. degree in the Institute of Theoretical Physics from Chinese Academy of Sciences in 2005, and then worked at Zhejiang University and Harvard University as post-doctoral research fellow. Liu is currently an associate professor in the College of Science, Zhejiang University of Technology. His research interests include image processing, structurepreserving numerical algorithm, and discrete dynamic system
Fangfang Dong received the B.C. degree in the College of Mathematics and Information Science from Henan Normal University and the Ph.D. degree in the Center of Mathematical Sciences at Zhejiang University. Dong is currently an instructor in School of Statistics and Mathematics , Zhejiang Gongshang University. Her research interests include image processing, partial differential equation.
Luwei Ding, lecturer of College of Science, Huazhong Agricultural University, graduated from Henan University with the bachelor's degree in 1999, acquired master’s degree of economics from Wuhan University of Technology in 2009, is currently working at the doctorate degree of economics management, Huazhong Agricultural university. His research interests include quantitative economics.
971
Numerical Results of Nonlinear Filtering Problem from Yau-Yau Method Zhen Liu Department of Mathematics, Zhejiang University of Technology， Hangzhou 310023，China Email: [email protected]
Fangfang Dong School of Statistics and Mathematics ,Zhejiang Gongshang University， Hangzhou 310018， China Email: [email protected]
Luwei Ding School of Science Huazhong Agricultural University， Wuhan 430070， China Email: [email protected] Abstract— In the paper, we introduce a kind of method to solve the nonlinear filtering problem. Firstly, we review the basic filtering problem and the reduction from robust Duncan-Mortensen-Zakai equation to Kolmogorov equation. Then we use the difference discrete method to solve the Kolmogorov equation. The result is given to prove that the solution of the difference scheme convergences pointwise to the solution of the initial-value problem of the Kolmogorov equation. At last, the numerical results show that the numerical method can give the exact result. Index Terms—Nonlinear filtering equation, DMZ equation, Kolmogorov equation
I. INTRODUCTION One of the key goals of data assimilation is to obtain, recursively in time, good estimates of the state of a stochastic dynamical system based on noisy partial observations of the same. The major breakthrough in this classical problem of signal analysis was the landmark work of Kalman and Bucy [1]. The power of KalmanBucy filtering, which can be gauged from the broad spectrum of disciplines where it is used, lies in its elegant simplicity. Since then, the Kalman-Bucy filtering has proved useful in many areas such as navigational and guidance systems, radar tracking, solar mapping, and satellite orbit determination. Ever since the technique of the Kalman-Bucy filtering was popularized, there has been an intense interest in developing nonlinear filtering theory. In the 1960s and early 1970s, the basic approach to nonlinear filtering theory was via innovation methods originally proposed by Kailath [2] and Frost and Kailath[3] and subsequently rigorous developed by Fujisaki, Kallianpur, and Kunita [4] in 1972. But the difficulty with this approach is that innovation process is not, in general, explicitly computable. In the late 1970s, Brockett and Clark [5], Brockett [6] and Mitter [7]
© 2012 ACADEMY PUBLISHER doi:10.4304/jcp.7.4.971-976
proposed the idea of using estimation algebras to construct a finite-dimensional nonlinear filtering. Yau and his coworkers have finally classified all finite dimensional estimation algebras of Maximal rank ([8][19]). Although it is an interesting and challenging problem to classify all finite dimensional filtering, it appears that from Yau’s previous works, finite dimensional filtering simply does not exist for many practical situations. In ([20-21]), Yau and Yau proved the existence and decay estimates of the solution to the DMZ(Duncan-MortensenZakai) equation under the assumption that the drift terms of the signal dynamic and observation dynamic respectively have linear growths. Later they [22] showed that the real-time solution of DMZ equation can be reduced to off-time solution of Kolmogorov equation if the drift terms have linear growths. Recently Yau and Yau ([23-24]) finally proved that the real-time solution of DMZ equation can be reduced to off-time solution of Kolmororov equation if the growth of the drift term of observation dynamic at infinity is faster than the growth of the drift term of signal dynamic at infinity. In section 2, we recall the basic filtering problem. In section 3, we review the reduction from robust DMZ equation to Kolmogorov equation of Yau-Yau method. In section 4, we shall use finite difference method to derive the solution of Kolmororov equation and the numerical simulation of Yau-Yau method is shown. II. INTRODUCTION OF NONLINEAR FILTERING PROBLEM The filtering problem considered here is based on the signal observation model dx(t ) = f ( x(t ))dt + g ( x(t ))dv(t ) , x(0) = x0 , (1) dy (t ) = h( x(t ))dt + dw(t ) ，
y (0) = 0 , n
p
m
(2) m
in which x, v, y and w are R , R , R and R valued processes and v, w have components that are independent,
972
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
standard Brownian processes. We further assume that n = p . f , g and h are vector-valued, matrix-valued and vector-valued C ∞ smooth functions. We shall refer to x(t ) as the state of the system and y (t ) as the observation at time t . Let ρ ( x, t ) denote the conditional probability density of the state given the observation { y ( s ) : 0 ≤ s ≤ t} . It is well known that ρ ( x, t ) is given by normalizing a function σ (t , x) that satisfies the following DuncanMortensen-Zakai equation m
dσ (t , x) = L0σ (t , x)dt + ∑ Liσ (t , x) dyi (t ), i =1
σ (0, x) = σ 0 ( x) , Where n n ∂f 1 m 1 n ∂2 ∂ L0 = ∑ 2 − ∑ f i − ∑ i − ∑ hi2 . 2 i =1 ∂xi i =1 ∂xi i =1 ∂xi 2 i =1
(3)
(4)
And for i = 1, 2, ⋅⋅⋅, m, Li is the zero degree differential operator of multiplication by hi , σ 0 is the probability density of the initial point. Equation (3) is a stochastic partial differential equation. In real application, we are interested in constructing robust state estimators from observed sample paths with some property of robustness. From the robust algorithms proposed by Davis [25], it define a new unormalized density m
u (t , x) = exp(−∑ hi ( x) yi (t ))σ (t , x).
(5)
i =1
Then we can reduce the equation (3) to the following equation, which is called robust DMZ equation, 1 ut = Δu + F (t , x) ⋅∇u + V (t , x)u , u (0, x) = σ 0 ( x) .(6) 2 where drift function F (t , x) = − f ( x) + ∇k (t , x) , source function 1 1 2 V (t , x) = −divf ( x) − h( x) + Δk (t , x) 2 2 1 2 − f ( x ) ⋅ ∇k (t , x) + ∇k (t , x) . 2
uk (τ k , x ) uniformly in x
(10)
Where C is the constant and Tk = sup1≤i ≤ n { τ i − τ i −1 } .
So the solution of (7) approximate the solution of robust D-M-Z equation very well in the pointwise sense. Therefore it remains to describe an algorithm to compute u k (τ k , x), which is based on the following property. Proposition 1. [23-24] u (t , x) satisfies the following Kolmogorov equation n ∂u 1 ∂u (t , x) = Δu (t , x) − ∑ fi ( x) (t , x) ∂t 2 ∂xi i =1 (11) n ∂f 1 m − (∑ i ( x) + ∑ hi2 ( x))u (t , x), 2 i =1 i =1 ∂xi for τ l −1 ≤ t ≤ τ l , if and only if m
u (t , x) = exp(−∑ hi ( x) yi (τ l ))u (t , x) i =1
satisfies the robust DMZ equation with observation being fixed at y (τ l ) ∂u 1 (t , x) = Δu (t , x) + F (τ l , x) ⋅∇u + V (τ l , x)u. (12) ∂t 2 From the property, u1 (τ 1 , x) can be computed by u1 (τ 1 , x) where u1 (t , x) satisfies the following Kolmogorov equation n ∂u1 ∂u 1 (t , x) = Δu1 (t , x) − ∑ f i ( x) 1 (t , x) 2 ∂t ∂xi i =1 , n m ∂fi 1 2 ( x) + ∑ hi ( x))u1 (t , x), − (∑ 2 i =1 i =1 ∂xi m
u1 (0, x) = σ 0 ( x) exp(∑ hi ( x) y i (τ 1 )) .
(13)
m
The fundamental problem of nonlinear filtering theory is how to solving the robust DMZ equation(6) in real time and in memoryless manner. In the section, we introduce the Yau-Yau algorithm which achieves this goal for a large class of filtering system with arbitrary initial distribution by reducing it to solve Kolmogorov equation. Suppose that u (t , x) is the solution of robust D-M-Z equation and we want to compute u (τ , x) .Let a
partition
of
[0,τ ] .And let ui (t , x) be a solution of the following
partial differential equation for τ i −1 ≤ t ≤ τ i
© 2012 ACADEMY PUBLISHER
k →0
In fact,
EQUATION
be
and u (τ , x) = lim T
i =1
III. YAU-YAU METHODS FOR THE ROBUST DMZ
Tk = {0 = τ 0 < τ 1 < ⋅⋅⋅ < τ k = τ }
∂ui 1 (t , x) = Δu (t , x) + F (τ i , x) ⋅ ∇u + V (τ i , x)u , (7) 2 ∂t ui (τ i −1 , x) = ui −1 (τ i −1 , x) . (8) Theorem 1. [23-24] Let u (t , x) and uk (t , x) be the solutions of (6) and (7) respectively. For any given ε > 0 and sufficiently large n , u (t , x) − uk (t , x) ≤ ε TeCT , (9)
u1 (τ 1 , x) = exp(−∑ hi ( x) y i (τ 1 ))u1 (τ 1 , x).
(14)
i =1
For i ≥ 2 , ui (τ i , x) can be computed by ui (τ i , x) where ui (t , x ) for τ i −1 ≤ t ≤ τ i satisfies the following Kolmogorov equation n ∂ui ∂u 1 (t , x) = Δui (t , x) − ∑ f i ( x) i (t , x) ∂t ∂xi 2 i =1 n
− (∑ i =1
m
∂fi 1 m ( x) + ∑ hi2 ( x))ui (t , x), ∂xi 2 i =1
ui (τ i −1 , x) = exp(∑ hi ( x)( y i (τ i )) − yi (τ i −1 ))ui −1 (τ i −1 , x). (15) i =1
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
973
Observe that in the algorithm at step i ,we only need the observation at time τ i and τ i −1 . We do not need any other previous observation data. Observe also that the Kolmogorov equation (15) is uniform for all time steps and it depends on observation y (t ) only initial condition.
convergences pointwise to the solution of the initialvalue problem 1 ∂ 2u ∂u ∂u (24) (t , x) = (t , x) + p ( x) (t , x) + q ( x)u (t , x). 2 2 ∂x ∂t ∂x Proof. We denote the exact solution to the initial-value problem by u = u (t , x) and set
IV. NUMERICAL METHOD FOR THE KOLMOGOROV
zkn = u (k Δx, nΔt ) − vkn . (25) If we insert u into equation (20), and multiply through by Δt , we see that ukn = u (k Δx, nΔt ) satisfies
EQUATION
Consider the Kolmogorov equation (11), we use finite difference method to approximate the solution of the equation. We take the Kolmogorov equation ∂u 1 ∂ 2u ∂u (t , x) = (t , x) + p( x) (t , x) + q ( x)u (t , x) , (16) ∂t 2 ∂x 2 ∂x as an example to explain the method, where ∂f 1 m p( x) = − f ( x), q( x) = −( ( x) + ∑ hi2 ( x)) . ∂x 2 i =1 To approximate the equation (16) by finite difference, we discretize the spatial and time domain by placing a grid over the domain. For convenience, we use a uniform grid, with grid spacing Δx and Δt .The space-time domain of our problem then is approximated by the lattice of points. We will attempt to approximate the solution to our problem at the points on the lattice. Notationally, we will define vkn to be a function at the point (k Δx, nΔt ) .To arrive at the finite difference equation used to approximation the Kolmogorov equation, we use the approximation v n +1 − vkn , (17) ut (nΔt , k Δx ) ≈ k Δt v n − 2vkn + vkn−1 u xx (nΔt , k Δx) ≈ k +1 , (18) Δx 2 v n − vkn−1 . (19) u x ( n Δt , k Δ x ) ≈ k Δx We see that 1 u t − u xx − p( x)u x − q( x)u 2 n +1 u − ukn 1 ukn+1 − 2ukn + ukn−1 = k − 2 Δt Δx 2 (20) u n − ukn−1 − p ( xk ) k − q( xk )ukn , Δx where the above expression assumes that the higher order derivatives of u at (k Δx, nΔt ) are bounded. Then we have the following theorem. Theorem 2. Suppose q( x k ) < 0 ,and 1 2Δt Δt + p ( xk ) + Δtq( xk ) ≥ 0, 2 Δx 2 Δx 1 Δt Δt − p( xk ) ≥ 0, 2 Δx 2 Δx then the solution of the difference scheme 1 Δt n vkn +1 = vkn + (vk +1 − 2vkn + vkn−1 ) 2 Δx 2 Δt + p ( xk )(vkn − vkn−1 ) + Δtq( xk )v kn , Δx 1−
© 2012 ACADEMY PUBLISHER
(21) (22)
(23)
ukn +1 = ukn +
1 Δt n (uk +1 − 2ukn + ukn−1 ) 2 Δx 2
Δt (26) p ( xk )(ukn − ukn−1 ) + Δtq ( xk )ukn Δx + O(Δt 2 + ΔxΔt ). Then by subtracting equation (26) from equation (23), we see that zkn satisfies +
1 2Δt Δt + p( xk ) + Δtq ( xk )) zkn 2 Δx 2 Δx 1 Δt Δt (27) +( − p( xk )) zkn−1 2 Δx 2 Δx 1 Δt n + zk +1 + O(Δt 2 + ΔxΔt ). 2 Δx 2 In order to obtain bounds for zkn , we need to ensure that the coefficients of the three terms on the right of this equation are all nonnegative and have a sum no greater than unity. We are therefore led to the condition 1 2Δt Δt + (28) 1− p ( xk ) + Δtq( xk ) ≥ 0, 2 Δx 2 Δx 1 Δt Δ t − (29) p ( xk ) ≥ 0, 2 Δx 2 Δx as well as q( xk ) < 0 . Then we have 1 2Δt Δt zkn +1 ≤ (1 − + p( xk ) + Δtq ( xk )) zkn 2 Δx 2 Δx 1 Δt Δt +( − p ( xk )) zkn−1 (30) 2 Δx 2 Δx 1 Δt n 2 + zk +1 + A(Δt + ΔxΔt ) 2 Δx 2 ≤ Z n + A(Δt 2 + ΔxΔt ), where A is the constant associated with the " O " terms and depends on the assumed bound of the higher order zkn +1 = (1 −
{ }
derivative of u , and Z n = sup k zkn . Then taking the sup over k on the left side of the above equation yields Z n +1 ≤ Z n + A(Δt 2 + ΔxΔt ) . (31) Applying the above equation repeatedly yields Z n +1 ≤ Z n + A(Δt 2 + ΔxΔt ) ≤ " (32) ≤ Z 0 + (n + 1) A(Δt 2 + ΔxΔt )
974
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
f 2 ( x, y ) = x 2 + y 2 , h1 ( x, y ) = x − y, h2 ( x, y ) = x + y,
Since Z 0 = 0 , u (k Δx, (n + 1)Δt ) − vkn +1 ≤ z n +1 , and (n + 1)Δt → t , then we have u (k Δx, (n + 1)Δt ) − vkn +1 ≤ (n + 1) A(Δt 2 + ΔxΔt ) → 0 , (33)
as Δt , Δx → 0 . ,
A. Numerical simulation Let f ( x) = cos( x), h( x) = sin( x), and f ( x) = x 2 , h( x) = sin( x) respectively, we can have figure 1 and 2,
and f1 ( x, y ) = − cos( x) − sin( y ), f 2 ( x, y ) = cos( x ) + sin( y ), h1 ( x, y ) = x − y, h2 ( x, y ) = x + y ,
σ 0 = exp(−5( x 2 + y 2 )),
where σ 0 = exp(−10 ∗ x ). The dashed line denotes the mean of the conditional probability density ρ and the solid line is generated from discrete extended Kalman filter equations. Figure 3 and 4 show the numerical behavior of the conditional probability density ρ in different time. Let f1 ( x, y ) = − x 2 − y 2 ,
we have figure (5,6), in which figure(5,6) are the means about the x, y respectively. The dashed line denotes the mean of the conditional probability density ρ and the solid line is generated from discrete extended Kalman filter equations. Figure (7,8) show the numerical behavior of the conditional probability density ρ in different time.
Fig. 1
Fig. 2
2
Fig. 3
© 2012 ACADEMY PUBLISHER
Fig. 4
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
975
Fig.5
Fig.6
Fig.7
Fig.8
National Natural Science Foundation of China (Grant No.10801045). V. CONCLUSION In the paper, we use the Yau-Yau method to solve the nonlinear filtering problem, which reduce the robust DMZ equation to Kolmogorov equation. We give a finite difference schemes to derive the numerical solution of the Kolmogorov equation. A proof is given to prove that the solution of the difference scheme convergences pointwise to the solution of the initial-value problem of the Kolmogorov equation. At last, we give some numerical experiments to test the method. ACKNOWLEDGMENT The authors are grateful to the anonymous referees for their valuable remarks and helpful suggestions, which have significantly improved the paper. This work is supported by Zhejiang Provincial Natural Science Foundation of China (Grant No. Y6100526 ) and by
© 2012 ACADEMY PUBLISHER
REFERENCES [1] R.E. Kalman and R.S. Bucy, New results in linear filtering and prediction theory,Trans. ASME Series D, J. Basic Engineering 83(1961), 95-108. [2] T. Kailath, An innovations approach to least-squares estimation, Part I: Linear filtering in additive white noise, IEEE Transactions on Automatic Control,13(1986), pp.646655. [3] P.A. Frost and T. Kailath, An innovations approach to leastsquares estimation II, IEEE Transactions on Automatic Control, 16(1971), pp.217-226. [4] M.Fujisaki, G. Kallianpur, and H. Kunita, stochastic differential equations for the nonlinear filtering problem, Osaka Journal of Mathematics, 2(1972), PP.19-40. [5] R.W. Brockett and J.M.C. Clark, The geometry of the conditional density functions, in: Analysis and Optimization of Stochastic Systems, O.L.R. Jacobs eral (Eds.), Academic Press, New York, (1980) 299-309. [6] R.W. Brockett, Nonlinear systems and nonlinear estimation theory, in: The mathematics of Filtering and Identification
976
and Application, M.Hazewinkel and J.C. Williams (Eds.), Reidel, Dordrecht, (1981), 479-504. [7] S.K. Mittar, On the analogy between mathematical problems of nonlinear filtering and quantum physics, Ricerche Automat., 10(1979), pp. 163-216. [8] S. S.-T. Yau, Recent results on nonlinear filtering: New class of finite dimensional filters, in Proceedings of the 29th Conf. Decision and Control, Honolulu, 1990, pp. 231-233. [9] S. S.-T. Yau, Finite dimensional filters with nonlinear drift I: A class of filters including both Kalman-Bucy filters and Benes filters, J. Math. Syst., Estimation and Control, 4(2) (1994), pp. 181-203. [10]Yau, S.-T., and Hu, G.-Q. New direct method for Yau filtering system with arbitrary conditions. In Proceedings of 35th IEEE Conference on Design and Control, Kobe, Japan, (1996) ,pp. 2539-254 [11] Chen, J., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift VI: Linear structure of - matrix. Mathematics of Control, Signals and Systems, 9 (1996), pp. 370-385. [12] S.-T. Yau and G.-Q. Hu, Direct method without Riccati equation for Kalman-Bucy filtering system with arbitrary initial conditions, in Proceedings of the 13th World Congress IFAC, vol. II, San Francisco, CA,(1996), pp. 469474. [13] Chen, J., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift VII: Mitter conjecture and structure of η. SIAM J. Control and Optimization, 35, 4 (1997), pp. 11161131. [14] Chen, J., Yau, S. S.-T., and Leung, C.-W. Finite dimensional filters with nonlinear drift VIII: Classification of finite dimensional estimation algebra of maximal rank with state space dimension 4. SIAM J. Control and Optimization, 35, 4 (1997), pp. 1132-1141. [15] Rasoulian, A., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift IX: Construction of finite dimensional estimation algebra of non-maximal rank. Systems Control Letters, 30 (1997), pp. 109-118. [16] S. S.-T. Yau, X. Wu and W. S. Wong, Hessian Matrix nondecomposition theorem, Math. Res. Lett., 6 (1999), pp. 663-673. [17] S. S.-T. Yau, Brockett’s problem on nonlinear filtering theory, in Lectures on Systems, Control and Information, AMS/IP, Stud. Adv. Math., 17 (2000), pp. 177-212. [18] Hu, G.-Q., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift X: Explicit solution of DMZ equation. IEEE Transactions on Automatic Control, 46, 1 (2001), pp.142-148. [19]S.S.-T. Yau and G.-Q. Hu, Classification of finite dimensional estimation algebras of maximal rank with arbitrary state-space dimension and Mitter Conjecture, International Journal of Control, 78:10(2005), pp.689-705. [20]Yau, S. S.-T., and Yau, S.-T. Finite dimensional filters with nonlinear drift III: Duncan-Mortensen-Zakai equation with arbitrary initial condition for linear filtering system and the
© 2012 ACADEMY PUBLISHER
JOURNAL OF COMPUTERS, VOL. 7, NO. 4, APRIL 2012
Benes filtering system. IEEE Transactions on Aerospace and Electronic Systems, 33, 4 (1997), pp.1277-1294. [21] Yau, S.-T., and Yau, S. S.-T. Finite dimensional filters with nonlinear drift XI: Explicit solution of the generalized Kolmogorov equation in the Brockett-Mitter program. Advances in Mathematics, 140 (1998), pp.156-189. [22] S.-T. Yau and S.S.-T Yau, Existence and uniqueness and decay estimates for the time dependent parabolic equation with application to Duncan-Mortensen-Zakai equation, Asian Journal of Mathematics, 2(1998), pp.1079-1149. [23] S.-T. Yau and S.S.-T Yau, Real time solution of nonlinear filtering problem without memory I, mathematical Research Letter, 7(2000), pp.671-693. [24] S.-T. Yau and S.S.-T Yau, Real time solution of nonlinear filtering problem without memory II, SIAM. J. Control. Optim. 47(1)(2008). pp. 163-195. [25]Davis, M.H.A., On a multiplicative functional transformation arising in nonlinear filtering theory, Z.Wahrsch Verw. Gebiete, 54(1980), pp.125-139.
Zhen Liu received his Ph.D. degree in the Institute of Theoretical Physics from Chinese Academy of Sciences in 2005, and then worked at Zhejiang University and Harvard University as post-doctoral research fellow. Liu is currently an associate professor in the College of Science, Zhejiang University of Technology. His research interests include image processing, structurepreserving numerical algorithm, and discrete dynamic system
Fangfang Dong received the B.C. degree in the College of Mathematics and Information Science from Henan Normal University and the Ph.D. degree in the Center of Mathematical Sciences at Zhejiang University. Dong is currently an instructor in School of Statistics and Mathematics , Zhejiang Gongshang University. Her research interests include image processing, partial differential equation.
Luwei Ding, lecturer of College of Science, Huazhong Agricultural University, graduated from Henan University with the bachelor's degree in 1999, acquired master’s degree of economics from Wuhan University of Technology in 2009, is currently working at the doctorate degree of economics management, Huazhong Agricultural university. His research interests include quantitative economics.
