17 Modified Continuous-Time Particle Filter Algorithm …
251
ˆ
X k =
M
i=1
ω
i
k X
i
k
and the estimate of the posterior covariance matrix
ˆ
R k =
M
i=1
ω
i
k (X
i
k − ˆ
X k )(X
i
k − ˆ
X k )
T
.
Verify the condition T − t k = 0. If it is met, then stop. Obtain a realization of the
estimated state vector and the corresponding measurement vector at t k+1 :
X k+1 = F(t k , X k , h), Y k+1 = C(t k , X k , Y k , h), Z k =
Y k+1 − Y k
h
.
3. See Step 3 of Algorithm 1.
4. See Step 4 of Algorithm 1.
The given standard normalization procedure is involved in both the discretetime and continuous-time particle filters [3]. However, it might be insufficient when
μ(t k , X
i
k , Z k )h is greater or lower than some threshold value. In this situation, the
resampling procedure is also ineffective due to underflow or overflow errors. The
threshold value is determined by the data type used for storage of the floating-point
number. For example, when using the double-precision floating-point format (for
storage of the floating-point number it is required 64 bit) this threshold value is
about 706.893. Certainly, the value μ(t k , X
i
k , Z k )h can be lowered by the choice of
the integration step but decreasing h proportionally increases the calculation time.
And if it is required to solve the optimal filtering problem in real-time then such a
way may be unrealizable. If the single-precision floating-point format is used (for
storage of the floating-point number it is required 32 bit), then the threshold value is
just 88.722. For the extended precision floating point format (for storage of floatingpoint numbers it is required 80 bit) the threshold value is 11,356.523, but this data
type is used much more rarely in practice.
When overflow errors appear it is proposed to switch from the exponential function e
μ(t k ,X
i
k ,Z k )h to the expression involving the exponent μ(t k , X
i
k , Z k )h and the
quantity that provides the correctness of this procedure. In other words, it provides
the exponential function calculation without overflow errors. For that, instead of the
expression
ω
i
k+1 = ω
i
k e
μ(t k ,X
i
k ,Z k )h
it is used the relation
ω
i
k+1 = exp{ln ω
i
k + μ(t k , X
i
k , Z k )h − γ k },
251
ˆ
X k =
M
i=1
ω
i
k X
i
k
and the estimate of the posterior covariance matrix
ˆ
R k =
M
i=1
ω
i
k (X
i
k − ˆ
X k )(X
i
k − ˆ
X k )
T
.
Verify the condition T − t k = 0. If it is met, then stop. Obtain a realization of the
estimated state vector and the corresponding measurement vector at t k+1 :
X k+1 = F(t k , X k , h), Y k+1 = C(t k , X k , Y k , h), Z k =
Y k+1 − Y k
h
.
3. See Step 3 of Algorithm 1.
4. See Step 4 of Algorithm 1.
The given standard normalization procedure is involved in both the discretetime and continuous-time particle filters [3]. However, it might be insufficient when
μ(t k , X
i
k , Z k )h is greater or lower than some threshold value. In this situation, the
resampling procedure is also ineffective due to underflow or overflow errors. The
threshold value is determined by the data type used for storage of the floating-point
number. For example, when using the double-precision floating-point format (for
storage of the floating-point number it is required 64 bit) this threshold value is
about 706.893. Certainly, the value μ(t k , X
i
k , Z k )h can be lowered by the choice of
the integration step but decreasing h proportionally increases the calculation time.
And if it is required to solve the optimal filtering problem in real-time then such a
way may be unrealizable. If the single-precision floating-point format is used (for
storage of the floating-point number it is required 32 bit), then the threshold value is
just 88.722. For the extended precision floating point format (for storage of floatingpoint numbers it is required 80 bit) the threshold value is 11,356.523, but this data
type is used much more rarely in practice.
When overflow errors appear it is proposed to switch from the exponential function e
μ(t k ,X
i
k ,Z k )h to the expression involving the exponent μ(t k , X
i
k , Z k )h and the
quantity that provides the correctness of this procedure. In other words, it provides
the exponential function calculation without overflow errors. For that, instead of the
expression
ω
i
k+1 = ω
i
k e
μ(t k ,X
i
k ,Z k )h
it is used the relation
ω
i
k+1 = exp{ln ω
i
k + μ(t k , X
i
k , Z k )h − γ k },