[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 | |
---|
[131] | 16 | #include <stat/loggers.h> |
---|
| 17 | |
---|
[105] | 18 | using namespace itpp; |
---|
| 19 | |
---|
[131] | 20 | vec getPsi_a ( int t, mat &D, mat &Du , mat &X ); |
---|
| 21 | vec getPsi_b ( int t, mat &D, mat &Du , mat &X ); |
---|
[105] | 22 | |
---|
| 23 | int main() { |
---|
| 24 | // Kalman filter |
---|
| 25 | int Ndat = 90000; |
---|
| 26 | mat D; |
---|
[131] | 27 | mat Du; |
---|
| 28 | mat X; |
---|
[105] | 29 | |
---|
[131] | 30 | dirfilelog L ( "exp/sim_var_arx",1000 ); |
---|
| 31 | |
---|
[105] | 32 | it_file itf ( "sim_var.it" ); |
---|
| 33 | itf >> Name ( "D" ) >> D; |
---|
[131] | 34 | itf >> Name ( "Du" ) >> Du; |
---|
| 35 | itf >> Name ( "X" ) >> X; |
---|
[105] | 36 | |
---|
[131] | 37 | Array<std::string> Names = "{ia ib om dom r }"; |
---|
[105] | 38 | int rglen = Names.length(); |
---|
| 39 | //Regressor |
---|
[162] | 40 | RV rgr ( Names ); |
---|
[131] | 41 | mat V0 = 0.0001*eye ( rglen ); V0 ( 0,0 ) =200; |
---|
[105] | 42 | double nu0 = rglen+1; |
---|
| 43 | |
---|
| 44 | //Autoregressive model |
---|
[135] | 45 | ARX Ar_a ( rgr,V0,nu0 ,0.95 ); |
---|
| 46 | ARX Ar_b ( rgr,V0,nu0 ,0.95 ); |
---|
[170] | 47 | const epdf& pA= Ar_a._epdf(); |
---|
| 48 | const epdf& pB= Ar_b._epdf(); |
---|
[105] | 49 | |
---|
[162] | 50 | RV rta ( "{th_a }",vec_1 ( rglen ) ); |
---|
| 51 | RV rtb ( "{th_b }",vec_1 ( rglen ) ); |
---|
[131] | 52 | int tha_log = L.add ( rta,"" ); |
---|
| 53 | int thb_log = L.add ( rtb,"" ); |
---|
| 54 | L.init(); |
---|
| 55 | |
---|
[105] | 56 | vec Psi ( rglen ); |
---|
[131] | 57 | vec Save ( 13 ); |
---|
[105] | 58 | for ( int t=2; t<Ndat; t++ ) { |
---|
[131] | 59 | Psi =getPsi_a ( t, D,Du ,X ); |
---|
| 60 | Ar_a.bayes ( Psi ); |
---|
| 61 | Psi =getPsi_b ( t, D,Du ,X ); |
---|
| 62 | Ar_b.bayes ( Psi ); |
---|
| 63 | |
---|
| 64 | Save = pA.mean(); |
---|
| 65 | L.logit ( tha_log, Save ); |
---|
| 66 | Save = pB.mean(); |
---|
| 67 | L.logit ( thb_log, Save ); |
---|
| 68 | L.step(); |
---|
[105] | 69 | } |
---|
[162] | 70 | L.finalize(); |
---|
[105] | 71 | |
---|
[131] | 72 | ivec bestind = Ar_a.structure_est ( egiw ( rgr,V0,nu0 ) ); |
---|
| 73 | ivec bestind2 = Ar_b.structure_est ( egiw ( rgr,V0,nu0 ) ); |
---|
[105] | 74 | cout << bestind <<endl; |
---|
[131] | 75 | cout << bestind2 <<endl; |
---|
[105] | 76 | |
---|
| 77 | return 0; |
---|
| 78 | } |
---|
| 79 | |
---|
[131] | 80 | vec getPsi_a ( int t, mat &D, mat &Du, mat &X ) { |
---|
| 81 | vec Psi ( 5 ); |
---|
| 82 | Psi ( 0 ) = D ( t,0 )-Du ( t,0 ); // a=0 |
---|
[105] | 83 | |
---|
[131] | 84 | Psi ( 1 ) = D ( t,2 ); |
---|
| 85 | Psi ( 2 ) = D(t,2)-D(t-1,2 ); |
---|
| 86 | Psi ( 3 ) = D(t,0)-D(t-1,0 ); |
---|
| 87 | Psi ( 4 ) = X ( t,2 ) - X ( t-1,2 ); |
---|
[105] | 88 | |
---|
[131] | 89 | return Psi; |
---|
| 90 | } |
---|
[105] | 91 | |
---|
[131] | 92 | vec getPsi_b ( int t, mat &D, mat &Du, mat &X ) { |
---|
| 93 | vec Psi ( 5 ); |
---|
| 94 | Psi ( 0 ) = D ( t,1 )-Du ( t,1 ); //b=1 |
---|
| 95 | |
---|
| 96 | Psi ( 1 ) = D(t,3); |
---|
| 97 | Psi ( 2 ) = D(t,3)-D(t-1,3 ); |
---|
| 98 | Psi ( 3 ) = D ( t,1 )-D(t-1,1); |
---|
| 99 | Psi ( 4 ) = X ( t,2 ) - X ( t-1,2 ); |
---|
| 100 | |
---|
[105] | 101 | return Psi; |
---|
| 102 | } |
---|