[105] | 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/arx.h> |
---|
| 15 | |
---|
| 16 | using namespace itpp; |
---|
| 17 | |
---|
| 18 | vec getPsi ( int t, mat &D, mat &Q ); |
---|
| 19 | |
---|
| 20 | int main() { |
---|
| 21 | // Kalman filter |
---|
| 22 | int Ndat = 90000; |
---|
| 23 | mat Q; |
---|
| 24 | mat R; |
---|
| 25 | mat D; |
---|
| 26 | |
---|
| 27 | it_file itf ( "sim_var.it" ); |
---|
| 28 | itf >> Name ( "Q" ) >> Q; |
---|
| 29 | itf >> Name ( "R" ) >> R; |
---|
| 30 | itf >> Name ( "D" ) >> D; |
---|
| 31 | |
---|
| 32 | Array<std::string> Names = "{ia ib ua ub iam ibm uam ubm iamm ibmm uamm ubmm r }"; |
---|
| 33 | int rglen = Names.length(); |
---|
| 34 | //Regressor |
---|
| 35 | RV rgr ( linspace ( 1,rglen ),Names,ones_i ( rglen ),zeros_i ( rglen ) ); |
---|
| 36 | mat V0 = 0.0001*eye ( rglen ); V0 ( 0,0 ) *=10; |
---|
| 37 | double nu0 = rglen+1; |
---|
| 38 | |
---|
| 39 | //Autoregressive model |
---|
| 40 | ARX Ar ( rgr,V0,nu0 ); |
---|
| 41 | epdf& ARep = Ar._epdf(); |
---|
| 42 | |
---|
| 43 | vec Psi ( rglen ); |
---|
| 44 | for ( int t=2; t<Ndat; t++ ) { |
---|
| 45 | Psi =getPsi ( t, D,Q ); |
---|
| 46 | Ar.bayes ( Psi ); |
---|
| 47 | } |
---|
| 48 | |
---|
| 49 | vec th = ARep.mean(); |
---|
| 50 | th.del ( th.length()-1 ); //remove est of r |
---|
| 51 | ivec bestind = Ar.structure_est ( egiw ( rgr,V0,nu0 ) ); |
---|
| 52 | cout << bestind <<endl; |
---|
| 53 | cout << th <<endl; |
---|
| 54 | |
---|
| 55 | //Reconstruction |
---|
| 56 | vec Rrec ( Ndat ); |
---|
| 57 | for ( int t=2; t<Ndat; t++ ) { |
---|
| 58 | Psi =getPsi ( t, D,R ); |
---|
| 59 | Rrec ( t ) = th(bestind-1)*Psi ( bestind); |
---|
| 60 | } |
---|
| 61 | it_file itfr ( "Qrec.it" ); |
---|
| 62 | itfr <<Name ( "Rrec" ) <<Rrec; |
---|
| 63 | |
---|
| 64 | return 0; |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | vec getPsi ( int t, mat &D, mat &Q ) { |
---|
| 68 | vec Psi ( 13 ); |
---|
| 69 | Psi ( 0 ) = log(exp(1)+Q ( t,0 )); |
---|
| 70 | |
---|
| 71 | Psi ( 1 ) = D ( t,0 ); |
---|
| 72 | Psi ( 2 ) = D ( t,1 ); |
---|
| 73 | Psi ( 3 ) = D ( t,2 ); |
---|
| 74 | Psi ( 4 ) = D ( t,3 ); |
---|
| 75 | |
---|
| 76 | Psi ( 5 ) = D ( t-1,0 ); |
---|
| 77 | Psi ( 6 ) = D ( t-1,1 ); |
---|
| 78 | Psi ( 7 ) = D ( t-1,2 ); |
---|
| 79 | Psi ( 8 ) = D ( t-1,3 ); |
---|
| 80 | |
---|
| 81 | Psi ( 9 ) = D ( t-2,0 ); |
---|
| 82 | Psi ( 10 ) = D ( t-2,1 ); |
---|
| 83 | Psi ( 11 ) = D ( t-2,2 ); |
---|
| 84 | Psi ( 12 ) = D ( t-2,3 ); |
---|
| 85 | return Psi; |
---|
| 86 | } |
---|