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