# Discrete Kalman filter

**URL:** <https://discourse.julialang.org/t/discrete-kalman-filter/109014>\
**Category:** Modelling & Simulations\
**Created:** [January 19, 2024, 3:51pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014 "2024-01-19T15:51:12Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 19, 2024, 3:51pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/1 "2024-01-19T15:51:12Z")

</div>

## Kalman filter

A wind turbine is modeled under the assumption that the rotor is directly attached to the generator and that the generator torque can be controlled directly. Then we can write (without noise):  
\dot\omega\_\mathrm{r} = \frac{T\_\mathrm{a}-T\_\mathrm{g}}{J\_\mathrm{r}},  
where \omega\_\mathrm{r} is the rotor speed, T\_\mathrm{a} the aerodynamic torque, T\_\mathrm{g} the generator torque and J\_\mathrm{r} the combined inertia of the rotor and the generator.

#### Model equations

\left[\begin{array}{ccc} \dot\omega\_\mathrm{r} \\ \dot{T}\_\mathrm{a} \\ \end{array} \right] = \left[\begin{array}{ccc} 0 & 1/J\_\mathrm{r} \\ 0 & 0 \\ \end{array} \right] \left[\begin{array}{ccc} \omega\_\mathrm{r} \\ {T}\_\mathrm{a} \\ \end{array} \right] + \left[\begin{array}{ccc} -1/J\_\mathrm{r} \\ 0 \\ \end{array} \right] T\_\mathrm{g}+ \left[\begin{array}{ccc} w\_{\omega\_\mathrm{r}} \\ w\_\mathrm{T\_\mathrm{a}} \\ \end{array} \right],  
where w\_{\omega\_\mathrm{r}} and w\_\mathrm{T\_\mathrm{a}} are the process noise of \omega\_\mathrm{r} and \mathrm{T\_\mathrm{a}} respectively.

The output equation is defined by  
 y = \omega\_\mathrm{r} + v\_\mathrm{\omega\_\mathrm{r}},  
where v\_\mathrm{\omega\_\mathrm{r}} is the measurement noise.

**Question:**  
How can I transform this system into a discrete time system and determine the matrices A, B, C and D that I need for implementing it with [LowLevelParticleFilters.jl](https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/#Kalman-filter) ?

Input will be T\_\mathrm{g} and the noisy measurement of \omega\_\mathrm{r}, output filtered estimates for T\_\mathrm{a} and \omega\_\mathrm{r}.

---

<div class="post-metadata">

**Author:** ![tom-plaa](https://avatars.discourse-cdn.com/v4/letter/t/d26b3c/32.png) [@tom-plaa](https://discourse.julialang.org/u/tom-plaa)\
**Post date:** [January 19, 2024, 5:00pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/2 "2024-01-19T17:00:06Z")

</div>

Without having actually used the package, I’m just going to try to rewrite the equations here in the hope that it can help you see the problem in a different way and make it easier for you to reach the solution.

I would write the continuous version of the differential equation for \omega\_r as  
d\omega\_r(t) = \frac{T\_a - T\_g}{J\_r}dt + dw\_{\omega\_r}(t)

For this stochastic differential equation we can use different approximating numerical schemes, but if we stick with the basic Euler method we can write it as  
\omega\_r(t+\delta t) \approx \omega(t) + \frac{T\_a - T\_g}{J\_r}\delta t + w\_{\omega\_r}(t + \delta t) - w\_{\omega\_r}(t)

In this case you can define a random variable Z(\delta t) that is supposed to represent the increment distribution of the noise: w\_{\omega\_r}(t + \delta t) - w\_{\omega\_r}(t), for which the distribution depends on the elapsed time \delta t (the time step).

Would this help?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 19, 2024, 11:14pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/3 "2024-01-19T23:14:02Z")

</div>

Have you seen [Discretization · LowLevelParticleFilters Documentation](https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/discretization/)?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 20, 2024, 2:26am UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/4 "2024-01-20T02:26:42Z")

</div>

No, I overlooked that. Thanks for pointing out!

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 21, 2024, 4:33pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/5 "2024-01-21T16:33:49Z")

</div>

OK, I am on step further. I wrote the following code:

```julia
using ControlSystemsBase, LowLevelParticleFilters
function kalman()
    J = 4.047e7

    A = [0 1/J; 0 0]
    B = [-1/J; 0]
    C = [1 0]
    D = 0

    Ts = 0.1

    sys = ss(A,B,C,D)
    sysd = c2d(sys, Ts) # Discretize the dynamics

    σ1 = 0.5
    N = σ1*[1; 0]
    sys_w = ss(A,N,C,D)
    sys_wd = c2d(sys_w, Ts) # Discretize the noise system
    Nd = sys_wd.B # The discretized noise input matrix
    R1 = Nd*Nd' # The final discrete-time covariance matrix
    R2 = [1;;]

    kf = KalmanFilter(sysd.A, sysd.B, sysd.C, sysd.D, R1, R2)
 end

```

It fails with the error message:

```julia
julia> kalman()
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
  [1] checkpositivedefinite
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:67 [inlined]
  [2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::NoPivot; check::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:269
  [3] cholesky!
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:267 [inlined]
  [4] cholesky!(A::Matrix{Float64}, ::NoPivot; check::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:301
  [5] cholesky! (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:295 [inlined]
  [6] cholesky(A::Matrix{Int64}, ::NoPivot; check::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:401 [inlined]
  [7] cholesky (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:401 [inlined]
  [8] PDMats.PDMat(mat::Matrix{Float64})
    @ PDMats ~/.julia/packages/PDMats/cAM9h/src/pdmat.jl:19
  [9] MvNormal
    @ ~/.julia/packages/Distributions/UaWBm/src/multivariate/mvnormal.jl:201 [inlined]
 [10] MvNormal
    @ ~/.julia/packages/Distributions/UaWBm/src/multivariate/mvnormal.jl:218 [inlined]
 [11] KalmanFilter(A::Matrix{Float64}, B::Vector{Float64}, C::Matrix{Int64}, D::Int64, R1::Matrix{Float64}, R2::Matrix{Int64})
    @ LowLevelParticleFilters ~/.julia/packages/LowLevelParticleFilters/rO2as/src/kalman.jl:42
 [12] kalman()
    @ Main ~/repos/WindSpeedEstimators/src/kalman.jl:24
 [13] top-level scope
    @ REPL[1]:1

```

Any idea?

In addition, which matrix is not positive definite but has to be?

UPDATE:  
The following lines cause (inside KalmanFilter) the problem:

```julia
julia> using Distributions
julia> MvNormal(R1)
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.

```

```julia
julia> R1
2×2 Matrix{Float64}:
 0.0025 0.0
 0.0 0.0

```

But why?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 22, 2024, 4:10am UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/6 "2024-01-22T04:10:48Z")

</div>

You can solve this problem by either

- Adding a tiny scalar times the identity matrix to R\_1
- Use the Square-root Kalman filter [`SqKalmanFilter`](https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/api/#LowLevelParticleFilters.SqKalmanFilter) type. This type does not have a user-friendly constructor that accepts pre-constructed Cholesky factors, so the option above is simpler to get you started.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 22, 2024, 4:23am UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/7 "2024-01-22T04:23:58Z")

</div>

Looking a bit closer at your system definition, the system you have posted is non-minimal (not controllable), with both a pole and a zero in the origin. Is this expected?

Further, you mention in your first post that you have two process-noise sources, w\_{\omega\_r} and w\_{T\_a}, this is not reflected in your modeled continuous-time covariance matrix, which would for such a setup be diagonal and non-singular

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 22, 2024, 10:17am UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/8 "2024-01-22T10:17:18Z")

</div>

> [@baggepinnen](#):
>
> Looking a bit closer at your system definition, the system you have posted is non-minimal (not controllable), with both a pole and a zero in the origin. Is this expected?

The system definition is a simplified version of Eq. (4) in [1]. They write:  
“The … model given by (2) is further augmented by appending the Ta as an additional state.”

I assume that augmenting the model made it non-minimal, and I think it was done to take the noise  
of Ta into account.

The estimator is using \omega and T\_r as noisy inputs and estimates the real \omega.

[1] [Kalman Filter-based Wind Speed Estimation for Wind Turbine Control](https://link.springer.com/article/10.1007/s12555-016-0537-1)

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 22, 2024, 10:20am UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/9 "2024-01-22T10:20:11Z")

</div>

> [@baggepinnen](#):
>
> Further, you mention in your first post that you have two process-noise sources, w\_{\omega\_r}wωrw\_{\omega\_r} and w\_{T\_a}wTaw\_{T\_a}, this is not reflected in your modeled continuous-time covariance matrix, which would for such a setup be diagonal and non-singular

What is the continues-time covariance matrix in the following piece of code:

```julia
σ1 = 0.5
N = σ1*[1; 0]
sys_w = ss(A,N,C,D)
sys_wd = c2d(sys_w, Ts) # Discretize the noise system
Nd = sys_wd.B # The discretized noise input matrix
R1 = Nd*Nd' # The final discrete-time covariance matrix
R2 = [1;;]

```

It cannot be N, because N is a column vector, I cannot replace it with a diagonal matrix…

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 22, 2024, 11:12am UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/10 "2024-01-22T11:12:03Z")

</div>

> [@ufechner7](#):
>
> What is the continues-time covariance matrix in the following piece of code:

`N*N'`

In the model you have at the moment, you are assuming that the dynamics of the second state variable is completely deterministic, this is what is causing the singular exception. Does this make sense? If so, what is the reason for adding this state variable rather than treating it as a known input? If you look through your equaitons, there is nothing at all affecting the second state variable, the seond row of A is all zeros, and both `B[2]` and `N[2]` are zero.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 22, 2024, 11:44am UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/11 "2024-01-22T11:44:38Z")

</div>

> [@baggepinnen](#):
>
> If so, what is the reason for adding this state variable rather than treating it as a known input?

Well, Ta is not known, the only thing we know about Ta is the mean and the variance. We  
want to estimate it based on omega and Tg, which are the (noisy) inputs…

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 22, 2024, 4:25pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/12 "2024-01-22T16:25:50Z")

</div>

So why does the covariance matrix not reflect that T\_a is affected by noise?

> [@ufechner7](#):
>
> `N = σ1*[1; 0]`

you only model noise acting on \omega\_r.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 23, 2024, 12:18pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/13 "2024-01-23T12:18:06Z")

</div>

I use now this code:

```julia
using ControlSystemsBase
import LowLevelParticleFilters as ll

J = 4.047e7
sigma_p = 0.01 # process noise
sigma_m = 1000 # measurement noise
Ts = 0.1 # sampling time [s]

A = [0 1/J; 0 0]
B = [-1/J; 0;;]
C = [1 0]
D = [0;;]

sys = ss(A, B, C, D) # Define the linear system
sysd = c2d(sys, Ts) # Discretize the dynamics

N = sigma_p * [1; 1;;]
sys_w = ss(A, N, C, D) # Define the noise system
sys_wd = c2d(sys_w, Ts) # Discretize the noise system
Nd = sys_wd.B # The discretized noise input matrix
R1 = Nd*Nd' + 1e-6*I # The final discrete-time covariance matrix
R2 = [sigma_m;;]

kf = ll.KalmanFilter(sysd.A, sysd.B, sysd.C, sysd.D, R1, R2)

```

Did I derive R1 and R2 now correctly, assuming “the process and measurement noises are set to 10e-2 and 1000, respectively”?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 23, 2024, 12:39pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/14 "2024-01-23T12:39:48Z")

</div>

> [@ufechner7](#):
>
> Did I derive R1 and R2 now correctly, assuming “the process and measurement noises are set to 10e-2 and 1000, respectively”?

Not quite, this

```julia
N = sigma_p * [1; 1;;]

```

indicates that _the same_ source of noise acts on _both_ state variables. What you probably intend to write is

```julia
N = sigma_p * I(2)

```

indicating that there are two independent sources of dynamics noise, each affecting one of the state variables. If this is what you mean, you have a full-rank covariance matrix which can be discretized directly without ZoH assumption using

```julia
R1c = sigma_p^2 * I(2) # Notice the square here since this is a covariance matrix
R1d = c2d(sysd, R1c, Ts) # Discrete version of R1

```

A second point is that the two state variables are describing two completely different things, so it would probably make sense to give them unique variances, rather than a single variance parameter `sigma_p^2` for both of them, the state variables do not even have the same unit so neither does the noise acting on them 🙂

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 23, 2024, 12:59pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/15 "2024-01-23T12:59:53Z")

</div>

Next try:

```julia
using ControlSystemsBase
import LowLevelParticleFilters as ll

J = 4.047e7 # total inertia on the side of the rotor
sigma_p = 0.01 # process noise
sigma_m = 1000 # measurement noise
Ts = 0.1 # sampling time [s]

A = [0 1/J; 0 0]
B = [-1/J; 0;;]
C = [1 0]
D = [0;;]

sys = ss(A,B,C,D) # linear system
sysd = c2d(sys, Ts) # discretize the dynamics

R1c = sigma_p^2 * I(2) # linear covariance matrix
R1d = c2d(sys, R1c, Ts) # discretized covariance matrix
R2d = [sigma_m;;]

kf = ll.KalmanFilter(sysd.A, sysd.B, sysd.C, sysd.D, R1d, R2d)

```

Is this correct now?

I think you had a typo in this line: `R1d = c2d(sysd, R1c, Ts) `. The c2d function does not accept a discrete system as input.

You further wrote: “A second point is that the two state variables are describing two completely different things, so it would probably make sense to give them unique variances, rather than a single variance parameter…”

I completely agree, but in the moment I just try to reproduce the results of the paper mentioned above, and they used a single variance parameter.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 23, 2024, 1:09pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/16 "2024-01-23T13:09:47Z")

</div>

Looks alright

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 23, 2024, 1:17pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/17 "2024-01-23T13:17:51Z")

</div>

Extra question:  
I use:

```julia
    kf.x .= [omega_0, ta_0]

```

to define the initial state of the Kalman filter. Is that a good way to do it, or am I using an internal property of the filter object that can change in future versions and there is a better way of doing it?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 23, 2024, 3:24pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/18 "2024-01-23T15:24:48Z")

</div>

> [@ufechner7](#):
>
> is a better way of doing it?

The intended way is to provide the distribution of the initial state, you do this by providing the `d0` argument [here](https://baggepinnen.github.io/LowLevelParticleFilters.jl/stable/api/#LowLevelParticleFilters.KalmanFilter), for example

```julia
d0 = MvNormal(x0, R0)

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [January 23, 2024, 3:49pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/19 "2024-01-23T15:49:16Z")

</div>

And what would be R0?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 23, 2024, 4:32pm UTC](https://discourse.julialang.org/t/discrete-kalman-filter/109014/20 "2024-01-23T16:32:02Z")

</div>

The covariance matrix of the initial state estimate