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