Changeset 54 for tests

Show
Ignore:
Timestamp:
03/25/08 13:40:27 (17 years ago)
Author:
smidl
Message:

funkcni verze testu + pmsm

Location:
tests
Files:
1 added
3 modified

Legend:

Unmodified
Added
Removed
  • tests/CMakeLists.txt

    r42 r54  
    4040 
    4141# PMSM stuff 
     42add_executable (testbidiff testbidiff.cpp) 
     43target_link_libraries (testbidiff ${BdmLibs}) 
     44 
    4245include_directories (${BDM_SOURCE_DIR}/simulator_zdenek) 
    4346link_directories (${BDM_BINARY_DIR}/simulator_zdenek) 
    4447add_executable (pmsm_sim pmsm_sim.cpp) 
    45 target_link_libraries (pmsm_sim ${BdmLibs} pmsmsim) 
     48target_link_libraries (pmsm_sim ${BdmLibs} pmsmsim netcdf_c++) 
  • tests/pmsm.h

    r48 r54  
    11#ifndef PMSM_H 
    22#define PMSM_H 
     3 
     4#include <stat/libFN.h> 
     5 
    36 
    47//TODO hardcoded RVs!!! 
  • tests/pmsm_sim.cpp

    r48 r54  
    1919#include "simulator.h" 
    2020 
     21#include <netcdfcpp.h> 
     22void write_to_nc ( NcFile &nc, mat &X, std::string Xn, Array<std::string> A ) { 
     23        char tmpstr[200]; 
     24        int Len = X.rows(); 
     25 
     26        sprintf ( tmpstr,"%s.length",Xn.c_str() ); 
     27        NcDim* lengt = nc.add_dim ( tmpstr, ( long ) Len ); 
     28        for ( int j=0; j<X.cols(); j++ ) { 
     29                if ( j<A.length() ) 
     30                        sprintf ( tmpstr,"%s_%s",Xn.c_str(), ( A ( j ) ).c_str() ); 
     31                else 
     32                        sprintf ( tmpstr,"%s_%d",Xn.c_str(),j ); 
     33                // Create variables and their attributes 
     34                NcVar* P = nc.add_var ( tmpstr, ncDouble, lengt ); 
     35                const double* Dp = X._data(); 
     36                P->put ( &Dp[j*Len],Len ); 
     37        } 
     38} 
     39 
     40 
    2141using namespace itpp; 
    2242//!Extended Kalman filter with unknown \c Q 
     
    2949                //from EKF 
    3050                preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() ); 
    31         };  
     51        }; 
    3252}; 
    3353 
     
    3555int main() { 
    3656        // Kalman filter 
    37         int Ndat = 10000; 
     57        int Ndat = 30000; 
    3858        double h = 1e-6; 
    3959        int Nsimstep = 125; 
    40         int Npart = 100; 
    41          
     60        int Npart = 1000; 
     61 
    4262        // internal model 
    4363        IMpmsm fxu; 
     
    4868 
    4969        vec mu0= "0.0 0.0 0.0 0.0"; 
    50         vec Qdiag ( "0.05 0.05 0.001 0.001" ); //zdenek: 0.05 0.05 0.001 0.001 
    51         vec Rdiag ( "0.03 0.03" ); //var(diff(xth)) = "0.034 0.034" 
     70        vec Qdiag ( "0.01 0.01 0.00001 0.00001" ); //zdenek: 0.01 0.01 0.0001 0.0001 
     71        vec Rdiag ( "0.0005 0.0005" ); //var(diff(xth)) = "0.034 0.034" 
    5272        chmat Q ( Qdiag ); 
    5373        chmat R ( Rdiag ); 
    5474        EKFCh KFE ( rx,ry,ru ); 
    5575        KFE.set_parameters ( &fxu,&hxu,Q,R ); 
    56         KFE.set_est ( mu0, chmat ( 1000*ones ( 4 ) ) ); 
     76        KFE.set_est ( mu0, chmat ( 1*ones ( 4 ) ) ); 
    5777 
    58         RV rQ ( "100","{Q}","4","0" ); 
     78        RV rQ ( "100","{Q}","2","0" ); 
    5979        EKF_unQ KFEp ( rx,ry,ru,rQ ); 
    6080        KFEp.set_parameters ( &fxu,&hxu,Q,R ); 
    61         KFEp.set_est ( mu0, chmat ( 1000*ones ( 4 ) ) ); 
     81        KFEp.set_est ( mu0, chmat ( 1*ones ( 4 ) ) ); 
    6282 
    6383        mgamma evolQ ( rQ,rQ ); 
     
    6585        // initialize 
    6686        evolQ.set_parameters ( 100.0 ); //sigma = 1/10 mu 
    67         evolQ.condition ( "0.05 0.05 0.001 0.001" ); //Zdenek default 
     87        evolQ.condition ( "0.01 0.01" ); //Zdenek default 
    6888        epdf& pfinit=evolQ._epdf(); 
    6989        M.set_est ( pfinit ); 
    70         evolQ.set_parameters ( 1000.0 );  
     90        evolQ.set_parameters ( 1000.0 ); 
    7191 
    7292        // 
     
    7595        epdf& Mep = M._epdf(); 
    7696 
    77         mat Xt=zeros ( 9,Ndat ); //true state from simulator 
    78         mat Dt=zeros ( 4,Ndat ); //observation 
    79         mat XtE=zeros ( 4,Ndat ); 
    80         mat XtM=zeros ( 8,Ndat ); //Q + x 
     97        mat Xt=zeros ( Ndat ,9 ); //true state from simulator 
     98        mat Dt=zeros ( Ndat,4+2 ); //observation 
     99        mat XtE=zeros ( Ndat, 4 ); 
     100        mat XtM=zeros ( Ndat,6 ); //Q + x 
    81101 
    82102        // SET SIMULATOR 
     
    86106        static long k_rampa_tmp=0; 
    87107        vec dt ( 2 ); 
    88         vec dtVS =zeros( 2 ); 
    89         vec xtVS =zeros(4); 
    90         vec et ( 4 ); 
    91         vec wt ( 2 ); 
    92108        vec ut ( 2 ); 
    93         mat XtV=zeros ( 4,Ndat ); 
    94109 
    95110        for ( int tK=1;tK<Ndat;tK++ ) { 
     
    112127                dt ( 0 ) = KalmanObs[2]; 
    113128                dt ( 1 ) = KalmanObs[3]; 
    114                  
    115                 // My own simulator for testing : Asuming ut is OK 
    116                 NorRNG.sample_vector ( 4,et ); 
    117                 NorRNG.sample_vector ( 2,wt ); 
    118                 xtVS = fxu.eval ( xtVS,ut ) + Q.sqrt_mult ( et ); 
    119                 dtVS = hxu.eval ( xtVS,ut ) + R.sqrt_mult ( wt ); 
    120129 
    121130                //estimator 
     
    123132                M.bayes ( concat ( dt,ut ) ); 
    124133 
    125                 Xt.set_col ( tK,vec ( x,9 ) ); //vec from C-array 
    126                 Dt.set_col ( tK, concat ( dt,ut ) ); 
    127                 XtE.set_col ( tK,KFEep.mean() ); 
    128                 XtM.set_col ( tK,Mep.mean() ); 
    129                 XtV.set_col ( tK,xtVS ); 
     134                Xt.set_row ( tK,vec ( x,9 ) ); //vec from C-array 
     135                Dt.set_row ( tK, concat ( dt,ut,vec_1(sqrt(pow(ut(0),2)+pow(ut(1),2))), vec_1(sqrt(pow(dt(0),2)+pow(dt(1),2))) ) ); 
     136                XtE.set_row ( tK,KFEep.mean() ); 
     137                XtM.set_row ( tK,Mep.mean() ); 
    130138        } 
    131139 
     
    136144        fou << Name ( "xthE" ) << XtE; 
    137145        fou << Name ( "xthM" ) << XtM; 
    138         fou << Name ( "xthV" ) << XtV; 
    139146        //Exit program: 
     147 
     148        //////////////// 
     149        // Just Testing 
     150        NcFile nc ( "pmsm_sim.nc", NcFile::Replace ); // Create and leave in define mode 
     151        if ( ! nc.is_valid() ) {        std::cerr << "creation of NCFile failed."<<endl;} 
     152 
     153        write_to_nc ( nc,Xt,"X","{isa isb om th }" ); 
     154        write_to_nc ( nc,XtM,"XtM","{q1 q2 isa isb om th }" ); 
     155        write_to_nc ( nc,XtE,"XE","{isa isb om th }" ); 
     156        write_to_nc ( nc,Dt,"Dt","{isa isb ua ub }" ); 
    140157        return 0; 
     158} 
    141159 
    142 }