Changeset 131 for pmsm

Show
Ignore:
Timestamp:
07/07/08 15:48:31 (16 years ago)
Author:
smidl
Message:

Odhad Q s opravenym modelem synchronizace

Location:
pmsm
Files:
3 modified

Legend:

Unmodified
Added
Removed

  • pmsm/pmsm_mix.cpp

    r125 r131  
    3030        double h = 1e-6; 
    3131        int Nsimstep = 125; 
    32         int Npar = 2000; 
     32        int Npar = 10; 
    3333 
    3434        dirfilelog L("exp/pmsm_mix",1000); 
    … …  
    7979//      cout << Eevol.sample() <<endl; 
    8080 
    81         mepdf evolQR(rQR,rQR,&Ucom); 
     81        mepdf evolQR(rQR,rQR,&Eevol); 
    8282        MPF<EKFful_unQR> M ( rx,rQR, evolQR, evolQR, Npar, EKU ); 
    8383        M.set_est ( evolQR._epdf() ); 

  • pmsm/sim_var.cpp

    r117 r131  
    4040        vec dt ( 2 ); 
    4141        vec ut ( 2 ); 
     42        vec dut ( 4 ); 
     43        vec dit (2); 
    4244        vec xtm=zeros ( 4 ); 
    4345        vec xdif=zeros ( 4 ); 
    … …  
    7779        RV rR("11", "{R }", "4","0"); 
    7880        RV rUD("12 13 14 15", "{u_isa u_isb i_isa i_isb }", ones_i(4),zeros_i(4)); 
     81        RV rDu("16 17 18 19","{dux duy duxf duyf }",ones_i(4),zeros_i(4)); 
     82        RV rDi("20 21","{disa disb }",ones_i(2),zeros_i(2)); 
    7983        int X_log = L.add(rx,"X"); 
    8084        int Efix_log = L.add(rx,"XF"); 
    … …  
    8488        int R_log = L.add(rR,"R"); 
    8589        int D_log = L.add(rUD,"D"); 
     90        int Du_log = L.add(rDu,"Du"); 
     91        int Di_log = L.add(rDi,"Di"); 
    8692        L.init(); 
    8793 
    … …  
    97103                dt ( 0 ) = KalmanObs[2]; 
    98104                dt ( 1 ) = KalmanObs[3]; 
     105                dut ( 0 ) = KalmanObs[4]; 
     106                dut ( 1 ) = KalmanObs[5]; 
     107                dut ( 2 ) = KalmanObs[6]; 
     108                dut ( 3 ) = KalmanObs[7]; 
     109                dit ( 0 ) = KalmanObs[8]; 
     110                dit ( 1 ) = KalmanObs[9]; 
    99111 
    100112                xt = fxu.eval ( xtm,ut ); 
    … …  
    104116                xtm ( 0 ) =x[0];xtm ( 1 ) =x[1];xtm ( 2 ) =x[2];xtm ( 3 ) =x[3]; 
    105117                xdif = xtm-xt; 
     118                if (xdif(0)>3.0){ 
     119                        cout << "here" <<endl; 
     120                        } 
    106121                if ( xdif ( 3 ) >pi ) xdif ( 3 )-=2*pi; 
    107122                if ( xdif ( 3 ) <-pi ) xdif ( 3 ) +=2*pi; 
    … …  
    132147                L.logit(R_log,          vec(Rt._data(),4) );  
    133148                L.logit(D_log,  vec(KalmanObs,4) );  
     149                L.logit(Du_log, dut );  
     150                L.logit(Di_log, dit );  
    134151                 
    135152                L.step(false); 
    … …  
    137154 
    138155        L.step(true); 
    139         //L.itsave("sim_var.it");        
     156//      L.itsave("sim_var.it");  
    140157         
    141158 

  • pmsm/sim_var_arx.cpp

    r105 r131  
    1414#include <estim/arx.h> 
    1515 
     16#include <stat/loggers.h> 
     17 
    1618using namespace itpp; 
    1719 
    18 vec getPsi ( int t, mat &D, mat &Q ); 
     20vec getPsi_a ( int t, mat &D, mat &Du , mat &X ); 
     21vec getPsi_b ( int t, mat &D, mat &Du , mat &X ); 
    1922 
    2023int main() { 
    2124        // Kalman filter 
    2225        int Ndat = 90000; 
    23         mat Q; 
    24         mat R; 
    2526        mat D; 
     27        mat Du; 
     28        mat X; 
     29 
     30        dirfilelog L ( "exp/sim_var_arx",1000 ); 
    2631 
    2732        it_file itf ( "sim_var.it" ); 
    28         itf >> Name ( "Q" ) >> Q; 
    29         itf >> Name ( "R" ) >> R; 
    3033        itf >> Name ( "D" ) >> D; 
     34        itf >> Name ( "Du" ) >> Du; 
     35        itf >> Name ( "X" ) >> X; 
    3136 
    32         Array<std::string> Names = "{ia ib ua ub iam ibm uam ubm iamm ibmm uamm ubmm r }"; 
     37        Array<std::string> Names = "{ia ib om dom r }"; 
    3338        int rglen = Names.length(); 
    3439        //Regressor 
    3540        RV rgr ( linspace ( 1,rglen ),Names,ones_i ( rglen ),zeros_i ( rglen ) ); 
    36         mat V0 = 0.0001*eye ( rglen ); V0 ( 0,0 ) *=10; 
     41        mat V0 = 0.0001*eye ( rglen ); V0 ( 0,0 ) =200; 
    3742        double nu0 = rglen+1; 
    3843 
    3944        //Autoregressive model 
    40         ARX Ar ( rgr,V0,nu0 ); 
    41         epdf& ARep = Ar._epdf(); 
     45        ARX Ar_a ( rgr,V0,nu0 ,0.99 ); 
     46        ARX Ar_b ( rgr,V0,nu0 ,0.99 ); 
     47        epdf& pA= Ar_a._epdf(); 
     48        epdf& pB= Ar_b._epdf(); 
     49 
     50        RV rta ( "22","{th_a }",vec_1 ( rglen ),"0" ); 
     51        RV rtb ( "23","{th_b }",vec_1 ( rglen ),"0" ); 
     52        int tha_log = L.add ( rta,"" ); 
     53        int thb_log = L.add ( rtb,"" ); 
     54        L.init(); 
    4255 
    4356        vec Psi ( rglen ); 
     57        vec Save ( 13 ); 
    4458        for ( int t=2; t<Ndat; t++ ) { 
    45                 Psi =getPsi ( t, D,Q ); 
    46                 Ar.bayes ( Psi ); 
     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(); 
    4769        } 
     70        L.step ( true ); 
    4871 
    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 ) ); 
     72        ivec bestind = Ar_a.structure_est ( egiw ( rgr,V0,nu0 ) ); 
     73        ivec bestind2 = Ar_b.structure_est ( egiw ( rgr,V0,nu0 ) ); 
    5274        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; 
     75        cout << bestind2 <<endl; 
    6376 
    6477        return 0; 
    6578} 
    6679 
    67 vec getPsi ( int t, mat &D, mat &Q ) { 
    68         vec Psi ( 13 ); 
    69         Psi ( 0 ) = log(exp(1)+Q ( t,0 )); 
     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 
    7083 
    71         Psi ( 1 ) = D ( t,0 ); 
    72         Psi ( 2 ) = D ( t,1 ); 
    73         Psi ( 3 ) = D ( t,2 ); 
    74         Psi ( 4 ) = D ( t,3 ); 
     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 ); 
    7588 
    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 ); 
    8589        return Psi; 
    8690} 
     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}