1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#![allow(non_snake_case)]
use nalgebra as na;
use na::{allocator::Allocator, DefaultAllocator, Dim, U1, MatrixN, RealField, VectorN};
use crate::linalg::rcond;
use crate::matrix::{check_non_negativ};
use crate::models::{KalmanEstimator, KalmanState, ExtendedLinearObserver, ExtendedLinearPredictor, Estimator};
use crate::noise::{CorrelatedNoise};
use nalgebra::MatrixMN;
impl<N: RealField, D: Dim> KalmanState<N, D>
where
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
{
pub fn new_zero(d: D) -> KalmanState<N, D> {
KalmanState {
x: VectorN::zeros_generic(d, U1),
X: MatrixN::zeros_generic(d, d),
}
}
}
impl<N: RealField, D: Dim> Estimator<N, D> for KalmanState<N, D>
where
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
{
fn state(&self) -> Result<VectorN<N, D>, &'static str> {
Ok(self.x.clone())
}
}
impl<N: RealField, D: Dim> KalmanEstimator<N, D> for KalmanState<N, D>
where
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>,
{
fn init(&mut self, state: &KalmanState<N, D>) -> Result<(), &'static str> {
self.x = state.x.clone();
self.X = state.X.clone();
let rcond = rcond::rcond_symetric(&self.X);
check_non_negativ(rcond, "X not PSD")?;
Ok(())
}
fn kalman_state(&self) -> Result<KalmanState<N, D>, &'static str> {
Ok(self.clone())
}
}
impl<N: RealField, D: Dim> ExtendedLinearPredictor<N, D> for KalmanState<N, D>
where
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
{
fn predict(
&mut self,
x_pred: &VectorN<N, D>,
fx: &MatrixN<N, D>,
noise: &CorrelatedNoise<N, D>,
) -> Result<(), &'static str> {
self.x = x_pred.clone();
self.X.quadform_tr(N::one(), &fx, &self.X.clone(), N::zero());
self.X += &noise.Q;
Ok(())
}
}
impl<N: RealField, D: Dim, ZD: Dim> ExtendedLinearObserver<N, D, ZD> for KalmanState<N, D>
where
DefaultAllocator: Allocator<N, D, D> + Allocator<N, ZD, ZD> + Allocator<N, ZD, D> + Allocator<N, D, ZD> + Allocator<N, D> + Allocator<N, ZD>
{
fn observe_innovation(
&mut self,
s: &VectorN<N, ZD>,
hx: &MatrixMN<N, ZD, D>,
noise: &CorrelatedNoise<N, ZD>,
) -> Result<(), &'static str> {
let XHt = &self.X * hx.transpose();
let S = hx * &XHt + &noise.Q;
let SI = S.clone().cholesky().ok_or("S not PD in observe")?.inverse();
let W = XHt * SI;
self.x += &W * s;
self.X.quadform_tr(N::one().neg(), &W, &S, N::one());
Ok(())
}
}