| 1 | /* |
|---|
| 2 | \file |
|---|
| 3 | \brief Models for synchronous electric drive using IT++ and BDM |
|---|
| 4 | \author Vaclav Smidl. |
|---|
| 5 | |
|---|
| 6 | ----------------------------------- |
|---|
| 7 | BDM++ - C++ library for Bayesian Decision Making under Uncertainty |
|---|
| 8 | |
|---|
| 9 | Using IT++ for numerical operations |
|---|
| 10 | ----------------------------------- |
|---|
| 11 | */ |
|---|
| 12 | |
|---|
| 13 | #include <itpp/itbase.h> |
|---|
| 14 | #include <estim/libKF.h> |
|---|
| 15 | #include <estim/libPF.h> |
|---|
| 16 | #include <stat/libFN.h> |
|---|
| 17 | |
|---|
| 18 | #include "pmsm.h" |
|---|
| 19 | |
|---|
| 20 | using namespace itpp; |
|---|
| 21 | /* |
|---|
| 22 | // PMSM with Q on Ia and Ib given externally |
|---|
| 23 | class EKF_unQ : public EKF<chmat> , public BMcond { |
|---|
| 24 | public: |
|---|
| 25 | EKF_unQ( rx,ry,ru,rQ):EKF<chmat>(rx,ry,ru),BMcond(rQ){}; |
|---|
| 26 | void condition(const vec &Q0){}; |
|---|
| 27 | };*/ |
|---|
| 28 | |
|---|
| 29 | |
|---|
| 30 | int main() { |
|---|
| 31 | // Kalman filter |
|---|
| 32 | int Ndat = 1000; |
|---|
| 33 | |
|---|
| 34 | // cout << KF; |
|---|
| 35 | // internal model |
|---|
| 36 | IMpmsm fxu; |
|---|
| 37 | // Rs Ls dt Fmag(Ypm) kp p J Bf(Mz) |
|---|
| 38 | fxu.set_parameters ( 0.28, 0.003465, 20*1e-6, 0.1989, 1.5 ,4.0, 0.04, 0.0 ); |
|---|
| 39 | // observation model |
|---|
| 40 | OMpmsm hxu; |
|---|
| 41 | |
|---|
| 42 | vec mu0= "100 100 100 1"; |
|---|
| 43 | vec Qdiag ( "0.1 0.1 0.01 0.01" ); |
|---|
| 44 | vec Rdiag ( "0.02 0.02" ); |
|---|
| 45 | vec vQ = "0.01:0.01:100"; |
|---|
| 46 | chmat Q ( Qdiag ); |
|---|
| 47 | chmat R ( Rdiag ); |
|---|
| 48 | EKFCh KFE ( rx,ry,ru ); |
|---|
| 49 | KFE.set_est ( mu0, chmat ( 1000*ones ( 4 ) ) ); |
|---|
| 50 | KFE.set_parameters ( &fxu,&hxu,Q,R ); |
|---|
| 51 | |
|---|
| 52 | mat ll(100,Ndat); |
|---|
| 53 | |
|---|
| 54 | EKFCh* kfArray[100]; |
|---|
| 55 | |
|---|
| 56 | for ( int i=0;i<100;i++ ) { |
|---|
| 57 | vec Qid ( Qdiag ); |
|---|
| 58 | Qid ( 0 ) = vQ ( i ); Qid ( 1 ) = vQ ( i ); |
|---|
| 59 | kfArray[i]= new EKFCh ( rx,ry,ru ); |
|---|
| 60 | kfArray[i]->set_est ( mu0, chmat ( 1000*ones ( 4 ) ) ); |
|---|
| 61 | kfArray[i]->set_parameters ( &fxu,&hxu,chmat ( Qid ),R ); |
|---|
| 62 | } |
|---|
| 63 | |
|---|
| 64 | epdf& KFEep = KFE._epdf(); |
|---|
| 65 | |
|---|
| 66 | //simulator values |
|---|
| 67 | vec dt ( 2 ); |
|---|
| 68 | vec wt ( 2 ); |
|---|
| 69 | vec ut ( 2 ); |
|---|
| 70 | vec xt=mu0; |
|---|
| 71 | vec et ( 4 ); |
|---|
| 72 | |
|---|
| 73 | mat Xt=zeros ( 4,Ndat ); |
|---|
| 74 | mat XtE=zeros ( 4,Ndat ); |
|---|
| 75 | Xt.set_col ( 0,KFEep.mean() ); |
|---|
| 76 | |
|---|
| 77 | for ( int t=1;t<Ndat;t++ ) { |
|---|
| 78 | //simulator |
|---|
| 79 | UniRNG.sample_vector ( 2,wt ); |
|---|
| 80 | |
|---|
| 81 | if ( rem ( t,500 ) <200 ) ut = rem ( t,500 ) *ones ( 2 ); |
|---|
| 82 | else |
|---|
| 83 | ut=zeros ( 2 ); |
|---|
| 84 | |
|---|
| 85 | NorRNG.sample_vector ( 4,et ); |
|---|
| 86 | NorRNG.sample_vector ( 2,wt ); |
|---|
| 87 | xt = fxu.eval ( xt,ut ) + Q.sqrt_mult ( et ); |
|---|
| 88 | dt = hxu.eval ( xt,ut ) + R.sqrt_mult ( wt ); |
|---|
| 89 | |
|---|
| 90 | //estimator |
|---|
| 91 | KFE.bayes ( concat ( dt,ut ) ); |
|---|
| 92 | for ( int i=0;i<100;i++ ) {kfArray[i]->bayes( concat ( dt,ut ) );ll(i,t)=ll(i,t-1) + kfArray[i]->_ll(); |
|---|
| 93 | } |
|---|
| 94 | |
|---|
| 95 | Xt.set_col ( t,xt ); |
|---|
| 96 | XtE.set_col ( t,KFEep.mean() ); |
|---|
| 97 | } |
|---|
| 98 | |
|---|
| 99 | it_file fou ( "pmsm.it" ); |
|---|
| 100 | |
|---|
| 101 | fou << Name ( "xth" ) << Xt; |
|---|
| 102 | fou << Name ( "xthE" ) << XtE; |
|---|
| 103 | fou << Name ( "ll" ) << ll; |
|---|
| 104 | //Exit program: |
|---|
| 105 | return 0; |
|---|
| 106 | |
|---|
| 107 | } |
|---|