1 | /*! |
---|
2 | \file |
---|
3 | \brief Bayesian Filtering for linear Gaussian models (Kalman Filter) and extensions |
---|
4 | \author Vaclav Smidl. |
---|
5 | |
---|
6 | ----------------------------------- |
---|
7 | BDM++ - C++ library for Bayesian Decision Making under Uncertaint16y |
---|
8 | |
---|
9 | Using IT++ for numerical operations |
---|
10 | ----------------------------------- |
---|
11 | */ |
---|
12 | |
---|
13 | #ifndef EKFfix_H |
---|
14 | #define EKFfix_H |
---|
15 | |
---|
16 | |
---|
17 | #include <estim/kalman.h> |
---|
18 | #include "fixed.h" |
---|
19 | #include "matrix.h" |
---|
20 | #include "matrix_vs.h" |
---|
21 | #include "reference_Q15.h" |
---|
22 | #include "parametry_motoru.h" |
---|
23 | #include "mpf_double.h" |
---|
24 | #include "fast_exp.h" |
---|
25 | |
---|
26 | using namespace bdm; |
---|
27 | |
---|
28 | double minQ(double Q); |
---|
29 | |
---|
30 | void mat_to_int16(const imat &M, int16 *I); |
---|
31 | void vec_to_int16(const ivec &v, int16 *I); |
---|
32 | void UDtof(const mat &U, const vec &D, imat &Uf, ivec &Df, const vec &xref); |
---|
33 | |
---|
34 | #ifdef XXX |
---|
35 | /*! |
---|
36 | \brief Extended Kalman Filter with full matrices in fixed point16 arithmetic |
---|
37 | |
---|
38 | An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean. |
---|
39 | */ |
---|
40 | class EKFfixed : public BM { |
---|
41 | public: |
---|
42 | void init_ekf(double Tv); |
---|
43 | void ekf(double ux, double uy, double isxd, double isyd); |
---|
44 | |
---|
45 | /* Declaration of local functions */ |
---|
46 | void prediction(int16 *ux); |
---|
47 | void correction(void); |
---|
48 | void update_psi(void); |
---|
49 | |
---|
50 | /* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/ |
---|
51 | int16 Q[16]; /* matrix [4,4] */ |
---|
52 | int16 R[4]; /* matrix [2,2] */ |
---|
53 | |
---|
54 | int16 x_est[4]; |
---|
55 | int16 x_pred[4]; |
---|
56 | int16 P_pred[16]; /* matrix [4,4] */ |
---|
57 | int16 P_est[16]; /* matrix [4,4] */ |
---|
58 | int16 Y_mes[2]; |
---|
59 | int16 ukalm[2]; |
---|
60 | int16 Kalm[8]; /* matrix [5,2] */ |
---|
61 | |
---|
62 | int16 PSI[16]; /* matrix [4,4] */ |
---|
63 | |
---|
64 | int16 temp15a[16]; |
---|
65 | |
---|
66 | int16 cA, cB, cC, cG, cH; // cD, cE, cF, cI ... nepouzivane |
---|
67 | |
---|
68 | int32 temp30a[4]; /* matrix [2,2] - temporary matrix for inversion */ |
---|
69 | enorm<fsqmat> E; |
---|
70 | mat Ry; |
---|
71 | |
---|
72 | public: |
---|
73 | //! Default constructor |
---|
74 | EKFfixed ():BM(),E(),Ry(2,2){ |
---|
75 | int16 i; |
---|
76 | for(i=0;i<16;i++){Q[i]=0;} |
---|
77 | for(i=0;i<4;i++){R[i]=0;} |
---|
78 | |
---|
79 | for(i=0;i<4;i++){x_est[i]=0;} |
---|
80 | for(i=0;i<4;i++){x_pred[i]=0;} |
---|
81 | for(i=0;i<16;i++){P_pred[i]=0;} |
---|
82 | for(i=0;i<16;i++){P_est[i]=0;} |
---|
83 | P_est[0]=0x7FFF; |
---|
84 | P_est[5]=0x7FFF; |
---|
85 | P_est[10]=0x7FFF; |
---|
86 | P_est[15]=0x7FFF; |
---|
87 | for(i=0;i<2;i++){Y_mes[i]=0;} |
---|
88 | for(i=0;i<2;i++){ukalm[i]=0;} |
---|
89 | for(i=0;i<8;i++){Kalm[i]=0;} |
---|
90 | |
---|
91 | for(i=0;i<16;i++){PSI[i]=0;} |
---|
92 | |
---|
93 | set_dim(4); |
---|
94 | E._mu()=zeros(4); |
---|
95 | E._R()=zeros(4,4); |
---|
96 | init_ekf(0.000125); |
---|
97 | }; |
---|
98 | //! Here dt = [yt;ut] of appropriate dimensions |
---|
99 | void bayes ( const vec &yt, const vec &ut ); |
---|
100 | //!dummy! |
---|
101 | const epdf& posterior() const {return E;}; |
---|
102 | |
---|
103 | }; |
---|
104 | |
---|
105 | UIREGISTER(EKFfixed); |
---|
106 | |
---|
107 | #endif |
---|
108 | |
---|
109 | //! EKF for testing q44 |
---|
110 | class EKFtest: public EKF_UD{ |
---|
111 | void bayes ( const vec &yt, const vec &cond ) { |
---|
112 | EKF_UD::bayes(yt,cond); |
---|
113 | vec D = prior()._R()._D(); |
---|
114 | |
---|
115 | if (D(3)>10) D(3) = 10; |
---|
116 | |
---|
117 | prior()._R().__D()=D; |
---|
118 | } |
---|
119 | }; |
---|
120 | UIREGISTER(EKFtest); |
---|
121 | |
---|
122 | /*! |
---|
123 | \brief Extended Kalman Filter with UD matrices in fixed point16 arithmetic |
---|
124 | |
---|
125 | An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean. |
---|
126 | */ |
---|
127 | class EKFfixedUD : public BM { |
---|
128 | public: |
---|
129 | LOG_LEVEL(EKFfixedUD,logU, logG, logD, logA, logP); |
---|
130 | |
---|
131 | void init_ekf(double Tv); |
---|
132 | void ekf(double ux, double uy, double isxd, double isyd); |
---|
133 | |
---|
134 | /* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/ |
---|
135 | int16 Q[16]; /* matrix [4,4] */ |
---|
136 | int16 R[4]; /* matrix [2,2] */ |
---|
137 | |
---|
138 | int16 x_est[4]; /* estimate and prediction */ |
---|
139 | |
---|
140 | int16 PSI[16]; /* matrix [4,4] */ |
---|
141 | int16 PSIU[16]; /* matrix PIS*U, [4,4] */ |
---|
142 | |
---|
143 | int16 Uf[16]; // upper triangular of covariance (inplace) |
---|
144 | int16 Df[4]; // diagonal covariance |
---|
145 | int16 Dfold[4]; // temp of D |
---|
146 | int16 G[16]; // temp for bierman |
---|
147 | |
---|
148 | int16 cA, cB, cC, cG, cH; // cD, cE, cF, cI ... nepouzivane |
---|
149 | |
---|
150 | enorm<fsqmat> E; |
---|
151 | mat Ry; |
---|
152 | |
---|
153 | public: |
---|
154 | //! Default constructor |
---|
155 | EKFfixedUD ():BM(),E(),Ry(2,2){ |
---|
156 | int16 i; |
---|
157 | for(i=0;i<16;i++){Q[i]=0;} |
---|
158 | for(i=0;i<4;i++){R[i]=0;} |
---|
159 | |
---|
160 | for(i=0;i<4;i++){x_est[i]=0;} |
---|
161 | for(i=0;i<16;i++){Uf[i]=0;} |
---|
162 | for(i=0;i<4;i++){Df[i]=0;} |
---|
163 | for(i=0;i<16;i++){G[i]=0;} |
---|
164 | for(i=0;i<4;i++){Dfold[i]=0;} |
---|
165 | |
---|
166 | for(i=0;i<16;i++){PSI[i]=0;} |
---|
167 | |
---|
168 | set_dim(4); |
---|
169 | dimy = 2; |
---|
170 | dimc = 2; |
---|
171 | E._mu()=zeros(4); |
---|
172 | E._R()=zeros(4,4); |
---|
173 | init_ekf(0.000125); |
---|
174 | }; |
---|
175 | //! Here dt = [yt;ut] of appropriate dimensions |
---|
176 | void bayes ( const vec &yt, const vec &ut ); |
---|
177 | //!dummy! |
---|
178 | const epdf& posterior() const {return E;}; |
---|
179 | void log_register(logger &L, const string &prefix){ |
---|
180 | BM::log_register ( L, prefix ); |
---|
181 | |
---|
182 | L.add_vector ( log_level, logG, RV("G",16), prefix ); |
---|
183 | L.add_vector ( log_level, logU, RV ("U", 16 ), prefix ); |
---|
184 | L.add_vector ( log_level, logD, RV ("D", 4 ), prefix ); |
---|
185 | L.add_vector ( log_level, logA, RV ("A", 16 ), prefix ); |
---|
186 | L.add_vector ( log_level, logP, RV ("P", 16 ), prefix ); |
---|
187 | |
---|
188 | }; |
---|
189 | //void from_setting(); |
---|
190 | }; |
---|
191 | |
---|
192 | UIREGISTER(EKFfixedUD); |
---|
193 | |
---|
194 | /*! |
---|
195 | * \brief Extended Kalman Filter with UD matrices in fixed point16 arithmetic |
---|
196 | * |
---|
197 | * An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean. |
---|
198 | */ |
---|
199 | class EKFfixedUD2 : public BM { |
---|
200 | public: |
---|
201 | LOG_LEVEL(EKFfixedUD2,logU, logG, logD, logA, logC, logP); |
---|
202 | |
---|
203 | void init_ekf2(double Tv); |
---|
204 | void ekf2(double ux, double uy, double isxd, double isyd); |
---|
205 | |
---|
206 | /* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/ |
---|
207 | int16 Q[4]; /* matrix [4,4] */ |
---|
208 | int16 R[4]; /* matrix [2,2] */ |
---|
209 | |
---|
210 | int16 x_est[2]; /* estimate and prediction */ |
---|
211 | int16 y_est[2]; /* estimate and prediction */ |
---|
212 | int16 y_old[2]; /* estimate and prediction */ |
---|
213 | |
---|
214 | int16 PSI[4]; /* matrix [4,4] */ |
---|
215 | int16 PSIU[4]; /* matrix PIS*U, [4,4] */ |
---|
216 | int16 C[4]; /* matrix [4,4] */ |
---|
217 | |
---|
218 | int16 Uf[4]; // upper triangular of covariance (inplace) |
---|
219 | int16 Df[2]; // diagonal covariance |
---|
220 | int16 Dfold[2]; // temp of D |
---|
221 | int16 G[4]; // temp for bierman |
---|
222 | |
---|
223 | int16 cA, cB, cC, cG, cH; // cD, cE, cF, cI ... nepouzivane |
---|
224 | |
---|
225 | enorm<fsqmat> E; |
---|
226 | mat Ry; |
---|
227 | |
---|
228 | public: |
---|
229 | //! Default constructor |
---|
230 | EKFfixedUD2 ():BM(),E(),Ry(2,2){ |
---|
231 | int16 i; |
---|
232 | for(i=0;i<4;i++){Q[i]=0;} |
---|
233 | for(i=0;i<4;i++){R[i]=0;} |
---|
234 | |
---|
235 | for(i=0;i<2;i++){x_est[i]=0;} |
---|
236 | for(i=0;i<2;i++){y_est[i]=0;} |
---|
237 | for(i=0;i<2;i++){y_old[i]=0;} |
---|
238 | for(i=0;i<4;i++){Uf[i]=0;} |
---|
239 | for(i=0;i<2;i++){Df[i]=0;} |
---|
240 | for(i=0;i<4;i++){G[i]=0;} |
---|
241 | for(i=0;i<2;i++){Dfold[i]=0;} |
---|
242 | |
---|
243 | for(i=0;i<4;i++){PSI[i]=0;} |
---|
244 | for(i=0;i<4;i++){C[i]=0;} |
---|
245 | |
---|
246 | set_dim(2); |
---|
247 | dimc = 2; |
---|
248 | dimy = 2; |
---|
249 | E._mu()=zeros(2); |
---|
250 | E._R()=zeros(2,2); |
---|
251 | init_ekf2(0.000125); |
---|
252 | }; |
---|
253 | //! Here dt = [yt;ut] of appropriate dimensions |
---|
254 | void bayes ( const vec &yt, const vec &ut ); |
---|
255 | //!dummy! |
---|
256 | const epdf& posterior() const {return E;}; |
---|
257 | void log_register(logger &L, const string &prefix){ |
---|
258 | BM::log_register ( L, prefix ); |
---|
259 | |
---|
260 | L.add_vector ( log_level, logG, RV("G2",4), prefix ); |
---|
261 | L.add_vector ( log_level, logU, RV ("U2", 4 ), prefix ); |
---|
262 | L.add_vector ( log_level, logD, RV ("D2", 2 ), prefix ); |
---|
263 | L.add_vector ( log_level, logA, RV ("A2", 4 ), prefix ); |
---|
264 | L.add_vector ( log_level, logC, RV ("C2", 4 ), prefix ); |
---|
265 | L.add_vector ( log_level, logP, RV ("P2", 4 ), prefix ); |
---|
266 | |
---|
267 | }; |
---|
268 | //void from_setting(); |
---|
269 | }; |
---|
270 | |
---|
271 | UIREGISTER(EKFfixedUD2); |
---|
272 | |
---|
273 | /*! |
---|
274 | * \brief Extended Kalman Filter with UD matrices in fixed point16 arithmetic |
---|
275 | * |
---|
276 | * An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean. |
---|
277 | */ |
---|
278 | class EKFfixedUD3 : public BM { |
---|
279 | public: |
---|
280 | LOG_LEVEL(EKFfixedUD3,logU, logG, logD, logA, logC, logP); |
---|
281 | |
---|
282 | void init_ekf3(double Tv); |
---|
283 | void ekf3(double ux, double uy, double isxd, double isyd); |
---|
284 | |
---|
285 | /* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/ |
---|
286 | int16 Q[9]; /* matrix [4,4] */ |
---|
287 | int16 R[4]; /* matrix [2,2] */ |
---|
288 | |
---|
289 | int16 x_est[3]; /* estimate and prediction */ |
---|
290 | int16 y_est[2]; /* estimate and prediction */ |
---|
291 | int16 y_old[2]; /* estimate and prediction */ |
---|
292 | |
---|
293 | int16 PSI[9]; /* matrix [4,4] */ |
---|
294 | int16 PSIU[9]; /* matrix PIS*U, [4,4] */ |
---|
295 | int16 C[6]; /* matrix [4,4] */ |
---|
296 | |
---|
297 | int16 Uf[9]; // upper triangular of covariance (inplace) |
---|
298 | int16 Df[3]; // diagonal covariance |
---|
299 | int16 Dfold[3]; // temp of D |
---|
300 | int16 G[9]; // temp for bierman |
---|
301 | |
---|
302 | int16 cA, cB, cC, cG, cF, cH; // cD, cE, cF, cI ... nepouzivane |
---|
303 | |
---|
304 | enorm<fsqmat> E; |
---|
305 | mat Ry; |
---|
306 | |
---|
307 | public: |
---|
308 | //! Default constructor |
---|
309 | EKFfixedUD3 ():BM(),E(),Ry(2,2){ |
---|
310 | int16 i; |
---|
311 | for(i=0;i<9;i++){Q[i]=0;} |
---|
312 | for(i=0;i<4;i++){R[i]=0;} |
---|
313 | |
---|
314 | for(i=0;i<3;i++){x_est[i]=0;} |
---|
315 | for(i=0;i<2;i++){y_est[i]=0;} |
---|
316 | for(i=0;i<2;i++){y_old[i]=0;} |
---|
317 | for(i=0;i<9;i++){Uf[i]=0;} |
---|
318 | for(i=0;i<3;i++){Df[i]=0;} |
---|
319 | for(i=0;i<4;i++){G[i]=0;} |
---|
320 | for(i=0;i<3;i++){Dfold[i]=0;} |
---|
321 | |
---|
322 | for(i=0;i<9;i++){PSI[i]=0;} |
---|
323 | for(i=0;i<6;i++){C[i]=0;} |
---|
324 | |
---|
325 | set_dim(3); |
---|
326 | dimc = 2; |
---|
327 | dimy = 2; |
---|
328 | E._mu()=zeros(3); |
---|
329 | E._R()=zeros(3,3); |
---|
330 | init_ekf3(0.000125); |
---|
331 | }; |
---|
332 | //! Here dt = [yt;ut] of appropriate dimensions |
---|
333 | void bayes ( const vec &yt, const vec &ut ); |
---|
334 | //!dummy! |
---|
335 | const epdf& posterior() const {return E;}; |
---|
336 | void log_register(logger &L, const string &prefix){ |
---|
337 | BM::log_register ( L, prefix ); |
---|
338 | }; |
---|
339 | //void from_setting(); |
---|
340 | }; |
---|
341 | |
---|
342 | UIREGISTER(EKFfixedUD3); |
---|
343 | |
---|
344 | /*! |
---|
345 | * \brief Extended Kalman Filter with Chol matrices in fixed point16 arithmetic |
---|
346 | * |
---|
347 | * An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean. |
---|
348 | */ |
---|
349 | class EKFfixedCh : public BM { |
---|
350 | public: |
---|
351 | LOG_LEVEL(EKFfixedCh,logCh, logA, logP); |
---|
352 | |
---|
353 | void init_ekf(double Tv); |
---|
354 | void ekf(double ux, double uy, double isxd, double isyd); |
---|
355 | |
---|
356 | /* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/ |
---|
357 | int16 Q[16]; /* matrix [4,4] */ |
---|
358 | int16 R[4]; /* matrix [2,2] */ |
---|
359 | |
---|
360 | int16 x_est[4]; /* estimate and prediction */ |
---|
361 | |
---|
362 | int16 PSI[16]; /* matrix [4,4] */ |
---|
363 | int16 PSICh[16]; /* matrix PIS*U, [4,4] */ |
---|
364 | |
---|
365 | int16 Chf[16]; // upper triangular of covariance (inplace) |
---|
366 | |
---|
367 | int16 cA, cB, cC, cG, cH; // cD, cE, cF, cI ... nepouzivane |
---|
368 | |
---|
369 | enorm<chmat> E; |
---|
370 | mat Ry; |
---|
371 | |
---|
372 | public: |
---|
373 | //! Default constructor |
---|
374 | EKFfixedCh ():BM(),E(),Ry(2,2){ |
---|
375 | int16 i; |
---|
376 | for(i=0;i<16;i++){Q[i]=0;} |
---|
377 | for(i=0;i<4;i++){R[i]=0;} |
---|
378 | |
---|
379 | for(i=0;i<4;i++){x_est[i]=0;} |
---|
380 | for(i=0;i<16;i++){Chf[i]=0;} |
---|
381 | |
---|
382 | for(i=0;i<16;i++){PSI[i]=0;} |
---|
383 | |
---|
384 | set_dim(4); |
---|
385 | dimc = 2; |
---|
386 | dimy =2; |
---|
387 | E._mu()=zeros(4); |
---|
388 | E._R()=zeros(4,4); |
---|
389 | init_ekf(0.000125); |
---|
390 | }; |
---|
391 | //! Here dt = [yt;ut] of appropriate dimensions |
---|
392 | void bayes ( const vec &yt, const vec &ut ); |
---|
393 | //!dummy! |
---|
394 | const epdf& posterior() const {return E;}; |
---|
395 | void log_register(logger &L, const string &prefix){ |
---|
396 | BM::log_register ( L, prefix ); |
---|
397 | |
---|
398 | L.add_vector ( log_level, logCh, RV ("Ch", 16 ), prefix ); |
---|
399 | L.add_vector ( log_level, logA, RV ("A", 16 ), prefix ); |
---|
400 | L.add_vector ( log_level, logP, RV ("P", 16 ), prefix ); |
---|
401 | |
---|
402 | }; |
---|
403 | //void from_setting(); |
---|
404 | }; |
---|
405 | |
---|
406 | UIREGISTER(EKFfixedCh); |
---|
407 | |
---|
408 | /*! |
---|
409 | * \brief Extended Kalman Filter with UD matrices in fixed point16 arithmetic |
---|
410 | * |
---|
411 | * An approximation of the exact Bayesian filter with Gaussian noices and non-linear evolutions of their mean. |
---|
412 | */ |
---|
413 | class EKFfixedCh2 : public BM { |
---|
414 | public: |
---|
415 | LOG_LEVEL(EKFfixedCh2,logCh, logA, logC, logP); |
---|
416 | |
---|
417 | void init_ekf2(double Tv); |
---|
418 | void ekf2(double ux, double uy, double isxd, double isyd); |
---|
419 | |
---|
420 | /* Constants - definovat jako konstanty ?? ?kde je vyhodnejsi aby v pameti byli?*/ |
---|
421 | int16 Q[4]; /* matrix [4,4] */ |
---|
422 | int16 R[4]; /* matrix [2,2] */ |
---|
423 | |
---|
424 | int16 x_est[2]; /* estimate and prediction */ |
---|
425 | int16 y_est[2]; /* estimate and prediction */ |
---|
426 | int16 y_old[2]; /* estimate and prediction */ |
---|
427 | |
---|
428 | int16 PSI[4]; /* matrix [4,4] */ |
---|
429 | int16 PSICh[4]; /* matrix PIS*U, [4,4] */ |
---|
430 | int16 C[4]; /* matrix [4,4] */ |
---|
431 | |
---|
432 | int16 Chf[4]; // upper triangular of covariance (inplace) |
---|
433 | |
---|
434 | int16 cA, cB, cC, cG, cH; // cD, cE, cF, cI ... nepouzivane |
---|
435 | |
---|
436 | enorm<fsqmat> E; |
---|
437 | mat Ry; |
---|
438 | |
---|
439 | public: |
---|
440 | //! Default constructor |
---|
441 | EKFfixedCh2 ():BM(),E(),Ry(2,2){ |
---|
442 | int16 i; |
---|
443 | for(i=0;i<4;i++){Q[i]=0;} |
---|
444 | for(i=0;i<4;i++){R[i]=0;} |
---|
445 | |
---|
446 | for(i=0;i<2;i++){x_est[i]=0;} |
---|
447 | for(i=0;i<2;i++){y_est[i]=0;} |
---|
448 | for(i=0;i<2;i++){y_old[i]=0;} |
---|
449 | for(i=0;i<4;i++){Chf[i]=0;} |
---|
450 | |
---|
451 | for(i=0;i<4;i++){PSI[i]=0;} |
---|
452 | for(i=0;i<4;i++){C[i]=0;} |
---|
453 | |
---|
454 | set_dim(2); |
---|
455 | dimc = 2; |
---|
456 | dimy = 2; |
---|
457 | E._mu()=zeros(2); |
---|
458 | E._R()=zeros(2,2); |
---|
459 | init_ekf2(0.000125); |
---|
460 | }; |
---|
461 | //! Here dt = [yt;ut] of appropriate dimensions |
---|
462 | void bayes ( const vec &yt, const vec &ut ); |
---|
463 | //!dummy! |
---|
464 | const epdf& posterior() const {return E;}; |
---|
465 | void log_register(logger &L, const string &prefix){ |
---|
466 | BM::log_register ( L, prefix ); |
---|
467 | |
---|
468 | L.add_vector ( log_level, logCh, RV ("Ch2", 4 ), prefix ); |
---|
469 | L.add_vector ( log_level, logA, RV ("A2", 4 ), prefix ); |
---|
470 | L.add_vector ( log_level, logC, RV ("C2", 4 ), prefix ); |
---|
471 | L.add_vector ( log_level, logP, RV ("P2", 4 ), prefix ); |
---|
472 | |
---|
473 | }; |
---|
474 | //void from_setting(); |
---|
475 | }; |
---|
476 | |
---|
477 | UIREGISTER(EKFfixedCh2); |
---|
478 | |
---|
479 | |
---|
480 | //! EKF for comparison of EKF_UD with its fixed-point16 implementation |
---|
481 | class EKF_UDfix : public BM { |
---|
482 | protected: |
---|
483 | //! logger |
---|
484 | LOG_LEVEL(EKF_UDfix,logU, logG); |
---|
485 | //! Internal Model f(x,u) |
---|
486 | shared_ptr<diffbifn> pfxu; |
---|
487 | |
---|
488 | //! Observation Model h(x,u) |
---|
489 | shared_ptr<diffbifn> phxu; |
---|
490 | |
---|
491 | //! U part |
---|
492 | mat U; |
---|
493 | //! D part |
---|
494 | vec D; |
---|
495 | |
---|
496 | mat A; |
---|
497 | mat C; |
---|
498 | mat Q; |
---|
499 | vec R; |
---|
500 | |
---|
501 | enorm<ldmat> est; |
---|
502 | |
---|
503 | |
---|
504 | public: |
---|
505 | |
---|
506 | //! copy constructor duplicated |
---|
507 | EKF_UDfix* _copy() const { |
---|
508 | return new EKF_UDfix(*this); |
---|
509 | } |
---|
510 | |
---|
511 | const enorm<ldmat>& posterior()const{return est;}; |
---|
512 | |
---|
513 | enorm<ldmat>& prior() { |
---|
514 | return const_cast<enorm<ldmat>&>(posterior()); |
---|
515 | } |
---|
516 | |
---|
517 | EKF_UDfix(){} |
---|
518 | |
---|
519 | |
---|
520 | EKF_UDfix(const EKF_UDfix &E0): pfxu(E0.pfxu),phxu(E0.phxu), U(E0.U), D(E0.D){} |
---|
521 | |
---|
522 | //! Set nonlinear functions for mean values and covariance matrices. |
---|
523 | void set_parameters ( const shared_ptr<diffbifn> &pfxu, const shared_ptr<diffbifn> &phxu, const mat Q0, const vec R0 ); |
---|
524 | |
---|
525 | //! Here dt = [yt;ut] of appropriate dimensions |
---|
526 | void bayes ( const vec &yt, const vec &cond = empty_vec ); |
---|
527 | |
---|
528 | void log_register ( bdm::logger& L, const string& prefix ){ |
---|
529 | BM::log_register ( L, prefix ); |
---|
530 | |
---|
531 | if ( log_level[logU] ) |
---|
532 | L.add_vector ( log_level, logU, RV ( dimension()*dimension() ), prefix ); |
---|
533 | if ( log_level[logG] ) |
---|
534 | L.add_vector ( log_level, logG, RV ( dimension()*dimension() ), prefix ); |
---|
535 | |
---|
536 | } |
---|
537 | /*! Create object from the following structure |
---|
538 | |
---|
539 | \code |
---|
540 | class = 'EKF_UD'; |
---|
541 | OM = configuration of bdm::diffbifn; % any offspring of diffbifn, bdm::diffbifn::from_setting |
---|
542 | IM = configuration of bdm::diffbifn; % any offspring of diffbifn, bdm::diffbifn::from_setting |
---|
543 | dQ = [...]; % vector containing diagonal of Q |
---|
544 | dR = [...]; % vector containing diagonal of R |
---|
545 | --- optional fields --- |
---|
546 | mu0 = [...]; % vector of statistics mu0 |
---|
547 | dP0 = [...]; % vector containing diagonal of P0 |
---|
548 | -- or -- |
---|
549 | P0 = [...]; % full matrix P0 |
---|
550 | --- inherited fields --- |
---|
551 | bdm::BM::from_setting |
---|
552 | \endcode |
---|
553 | If the optional fields are not given, they will be filled as follows: |
---|
554 | \code |
---|
555 | mu0 = [0,0,0,....]; % empty statistics |
---|
556 | P0 = eye( dim ); |
---|
557 | \endcode |
---|
558 | */ |
---|
559 | void from_setting ( const Setting &set ); |
---|
560 | |
---|
561 | void validate() {}; |
---|
562 | // TODO dodelat void to_setting( Setting &set ) const; |
---|
563 | |
---|
564 | }; |
---|
565 | UIREGISTER(EKF_UDfix); |
---|
566 | |
---|
567 | |
---|
568 | class MPF_pmsm_red:public BM{ |
---|
569 | double qom, qth, r; |
---|
570 | |
---|
571 | |
---|
572 | public: |
---|
573 | MPF_pmsm_red(){ |
---|
574 | dimy=2; |
---|
575 | dimc=2; |
---|
576 | qom=1e-1; |
---|
577 | qth=1e-6; |
---|
578 | r=1e-1; |
---|
579 | }; |
---|
580 | void bayes ( const vec &val, const vec &cond ) { |
---|
581 | /* const double &isa = val(0); |
---|
582 | const double &isb = val(1); |
---|
583 | const double &usa = cond(0); |
---|
584 | const double &usb = cond(1);*/ |
---|
585 | mpf_bayes((floatx)val(0),(floatx)val(1),(floatx)cond(0), (floatx)cond(1));//isa,isb,usa,usb); |
---|
586 | } |
---|
587 | |
---|
588 | class mp:public epdf{ |
---|
589 | LOG_LEVEL(mp,logth); |
---|
590 | public: |
---|
591 | mp():epdf(){set_dim(3); log_level[logth]=true; |
---|
592 | } |
---|
593 | vec sample() const {return zeros(3);} |
---|
594 | double evallog(const vec &v) const {return 0.0;} |
---|
595 | vec mean() const {vec tmp(3); |
---|
596 | floatx es,ec,eo; |
---|
597 | mpf_mean(&es, &ec, &eo); |
---|
598 | tmp(0)=es;tmp(1)=ec;tmp(2)=eo; |
---|
599 | return tmp; |
---|
600 | } |
---|
601 | vec variance() const {return zeros(3);} |
---|
602 | void log_register ( bdm::logger& L, const string& prefix ) { |
---|
603 | epdf::log_register ( L, prefix ); |
---|
604 | if ( log_level[logth] ) { |
---|
605 | int th_dim = N; // dimension - dimension of cov |
---|
606 | L.add_vector( log_level, logth, RV ( th_dim ), prefix ); |
---|
607 | } |
---|
608 | } |
---|
609 | void log_write() const { |
---|
610 | epdf::log_write(); |
---|
611 | if ( log_level[logth] ) { |
---|
612 | floatx Th[N]; |
---|
613 | mpf_th(Th); |
---|
614 | vec th(N); |
---|
615 | for(int i=0;i<N;i++){th(i)=Th[i];} |
---|
616 | log_level.store( logth, th ); |
---|
617 | } |
---|
618 | } |
---|
619 | }; |
---|
620 | |
---|
621 | mp mypdf; |
---|
622 | const mp& posterior() const {return mypdf;} |
---|
623 | |
---|
624 | void from_setting(const Setting &set){ |
---|
625 | BM::from_setting(set); |
---|
626 | UI::get ( log_level, set, "log_level", UI::optional ); |
---|
627 | |
---|
628 | UI::get(qom,set,"qom",UI::optional); |
---|
629 | UI::get(qth,set,"qth",UI::optional); |
---|
630 | UI::get(r,set,"r",UI::optional); |
---|
631 | } |
---|
632 | void validate(){ |
---|
633 | mpf_init((floatx)qom,(floatx)qth,(floatx)r); |
---|
634 | |
---|
635 | } |
---|
636 | }; |
---|
637 | UIREGISTER(MPF_pmsm_red); |
---|
638 | |
---|
639 | //! EKF with covariance R estimated by the VB method |
---|
640 | class EKFvbR: public EKFfull{ |
---|
641 | LOG_LEVEL(EKFvbR,logR,logQ); |
---|
642 | |
---|
643 | //! Statistics of the R estimator |
---|
644 | mat PsiR; |
---|
645 | //! Statistics of the Q estimator |
---|
646 | mat PsiQ; |
---|
647 | //! forgetting factor |
---|
648 | double phi; |
---|
649 | //! number of VB iterations |
---|
650 | int niter; |
---|
651 | //! degrees of freedom |
---|
652 | double nu; |
---|
653 | //! stabilizing element |
---|
654 | mat PsiR0; |
---|
655 | //! stabilizing element |
---|
656 | mat PsiQ0; |
---|
657 | |
---|
658 | void from_setting(const Setting &set){ |
---|
659 | EKFfull::from_setting(set); |
---|
660 | if (!UI::get(phi,set,"phi",UI::optional)){ |
---|
661 | phi = 0.99; |
---|
662 | } |
---|
663 | if (!UI::get(niter,set,"niter",UI::optional)){ |
---|
664 | niter = 3; |
---|
665 | } |
---|
666 | PsiQ = Q; |
---|
667 | PsiQ0 = Q; |
---|
668 | PsiR = R; |
---|
669 | PsiR0 = R; |
---|
670 | nu = 3; |
---|
671 | } |
---|
672 | void log_register ( logger &L, const string &prefix ) { |
---|
673 | EKFfull::log_register(L,prefix); |
---|
674 | L.add_vector(log_level, logR, RV("{R }",vec_1<int>(4)), prefix); |
---|
675 | L.add_vector(log_level, logQ, RV("{Q }",vec_1<int>(4)), prefix); |
---|
676 | } |
---|
677 | void bayes ( const vec &val, const vec &cond ) { |
---|
678 | vec diffx, diffy; |
---|
679 | mat Psi_vbQ; |
---|
680 | mat Psi_vbR; |
---|
681 | nu = phi*nu + (1-phi)*2 + 1.0; |
---|
682 | |
---|
683 | //save initial values of posterior |
---|
684 | vec mu0=est._mu(); |
---|
685 | fsqmat P0=est._R(); |
---|
686 | vec xpred = pfxu->eval(mu0,cond); |
---|
687 | |
---|
688 | for (int i=0; i<niter; i++){ |
---|
689 | est._mu() = mu0; |
---|
690 | est._R() = P0; |
---|
691 | |
---|
692 | EKFfull::bayes(val,cond); |
---|
693 | |
---|
694 | diffy = val - fy._mu(); |
---|
695 | Psi_vbR = phi*PsiR + (1-phi)*PsiR0+ outer_product(diffy,diffy)/*+C*mat(est._R())*C.T()*/; |
---|
696 | R = Psi_vbR/(nu-2); |
---|
697 | |
---|
698 | diffx = est._mu() - xpred; |
---|
699 | Psi_vbQ = phi*PsiQ + (1-phi)*PsiQ0+ outer_product(diffx,diffx) /*+mat(est._R()) /*+ A*mat(P0)*A.T()*/; |
---|
700 | Q = Psi_vbQ/(nu-2); |
---|
701 | } |
---|
702 | PsiQ = Psi_vbQ; |
---|
703 | PsiR = Psi_vbR; |
---|
704 | // cout <<"==" << endl << Psi << endl << diff << endl << P0 << endl << ":" << Q; |
---|
705 | log_level.store(logQ, vec(Q._M()._data(),4)); |
---|
706 | log_level.store(logR, vec(R._M()._data(),4)); |
---|
707 | |
---|
708 | } |
---|
709 | }; |
---|
710 | UIREGISTER(EKFvbR); |
---|
711 | |
---|
712 | |
---|
713 | #endif // KF_H |
---|
714 | |
---|