root/pmsm/sim_var_arx.cpp @ 171

Revision 170, 2.2 kB (checked in by smidl, 16 years ago)

Mixtures of EF and related changes to libEF and BM

  • Property svn:eol-style set to native
Line 
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#include <stat/loggers.h>
17
18using namespace itpp;
19
20vec getPsi_a ( int t, mat &D, mat &Du , mat &X );
21vec getPsi_b ( int t, mat &D, mat &Du , mat &X );
22
23int main() {
24        // Kalman filter
25        int Ndat = 90000;
26        mat D;
27        mat Du;
28        mat X;
29
30        dirfilelog L ( "exp/sim_var_arx",1000 );
31
32        it_file itf ( "sim_var.it" );
33        itf >> Name ( "D" ) >> D;
34        itf >> Name ( "Du" ) >> Du;
35        itf >> Name ( "X" ) >> X;
36
37        Array<std::string> Names = "{ia ib om dom r }";
38        int rglen = Names.length();
39        //Regressor
40        RV rgr ( Names );
41        mat V0 = 0.0001*eye ( rglen ); V0 ( 0,0 ) =200;
42        double nu0 = rglen+1;
43
44        //Autoregressive model
45        ARX Ar_a ( rgr,V0,nu0 ,0.95 );
46        ARX Ar_b ( rgr,V0,nu0 ,0.95 );
47        const epdf& pA= Ar_a._epdf();
48        const epdf& pB= Ar_b._epdf();
49
50        RV rta ( "{th_a }",vec_1 ( rglen ) );
51        RV rtb ( "{th_b }",vec_1 ( rglen ) );
52        int tha_log = L.add ( rta,"" );
53        int thb_log = L.add ( rtb,"" );
54        L.init();
55
56        vec Psi ( rglen );
57        vec Save ( 13 );
58        for ( int t=2; t<Ndat; t++ ) {
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();
69        }
70        L.finalize();
71
72        ivec bestind = Ar_a.structure_est ( egiw ( rgr,V0,nu0 ) );
73        ivec bestind2 = Ar_b.structure_est ( egiw ( rgr,V0,nu0 ) );
74        cout << bestind <<endl;
75        cout << bestind2 <<endl;
76
77        return 0;
78}
79
80vec 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
83
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 );
88
89        return Psi;
90}
91
92vec 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
101        return Psi;
102}
Note: See TracBrowser for help on using the browser.