root/pmsm/TR2245/pmsm_wishart.cpp @ 301

/*!
  \file
  \brief TR 2525 file for testing Wishart random walk on PMSM simulator
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/



#include <estim/libPF.h>
#include <estim/ekf_templ.h>
#include <stat/libFN.h>

#include <stat/loggers_ui.h>
#include <stat/libEF_ui.h>

#include "../pmsm_ui.h"

using namespace bdm;

int main ( int argc, char* argv[] ) {
        const char *fname;
        if ( argc>1 ) {fname = argv[1]; }
        else { fname = "pmsm_wishart.cfg"; }
        UIFile F ( fname );

        int Ndat;
        int Npart;
        double h = 1e-6;
        int Nsimstep = 125;

        vec Qdiag;
        vec Rdiag;

        pmsmDS* DS;

        double k;
        double l;
//      mpdf* evolQ ;
        try {
                // Kalman filter
                F.lookupValue ( "ndat", Ndat );
                F.lookupValue ( "Npart",Npart );
                
                F.lookupValue ( "k", k);
                F.lookupValue ( "l",l);

//              UIbuild ( F.lookup ( "Qrw" ),evolQ );
                Qdiag= getvec ( F.lookup ( "dQ" ) ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001
                Rdiag=getvec ( F.lookup ( "dR" ) );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034"
                
                UIbuild(F.lookup("system"),DS);
        }
        catch UICATCH;
// internal model

IMpmsm fxu;
//                  Rs    Ls        dt       Fmag(Ypm)    kp   p    J     Bf(Mz)
fxu.set_parameters ( 0.28, 0.003465, Nsimstep*h, 0.1989, 1.5 ,4.0, 0.04, 0.0 );
        // observation model
        OMpmsm hxu;

        vec mu0= "0.0 0.0 0.0 0.0";
        chmat Q ( Qdiag );
        chmat R ( Rdiag );
        EKFCh KFE ;
        KFE.set_parameters ( &fxu,&hxu,Q,R );
        KFE.set_est ( mu0, chmat ( 0.01*eye( 4 ) ) );
        KFE.set_rv ( rx );

        RV rQ ( "{Q }","16" );
        EKFCh_chQ KFEp ;
        KFEp.set_parameters ( &fxu,&hxu,Q,R );
        KFEp.set_est ( mu0, chmat (0.01* eye( 4 ) ) );

        rwiWishartCh* evolQw = new rwiWishartCh; 
        evolQw->set_parameters(4, k, sqrt(Qdiag),l);
        MPF<EKFCh_chQ> M;
        M.set_parameters ( evolQw,evolQw,Npart );
        // initialize
        chmat Ch0(diag(Qdiag));
        evolQw->condition ( vec(Ch0._Ch()._data(),16) ); //Zdenek default
        M.set_statistics ( evolQw->_e() , &KFEp );
        //

        M.set_rv ( concat ( rQ,rx ) );

        dirfilelog *L; UIbuild ( F.lookup ( "logger" ), L );// ( "exp/mpf_test",100 );
        
        KFE.set_options ( "logbounds" );
        KFE.log_add ( L,"KF" );
        M.set_options ( "logbounds" );
        M.log_add ( L,"M" );
        DS->log_add(*L);
        L->init();

        // SET SIMULATOR
        pmsmsim_set_parameters ( 0.28,0.003465,0.1989,0.0,4,1.5,0.04, 200., 3e-6, h );
        vec dt ( 2 );
        vec ut ( 2 );
        vec xt ( 4 );
        vec xtm=zeros ( 4 );
        double Ww=0.0;
        
        vec obs(6);
        for ( int tK=1;tK<Ndat;tK++ ) {
                //Number of steps of a simulator for one step of Kalman
                DS->step();
                DS->getdata(obs);
                
                dt=obs.left(2);
                ut=obs.right(2);
                //estimator
                KFE.bayes ( concat ( dt,ut ) );
                M.bayes ( concat ( dt,ut ) );

                DS->logit (*L);

                KFE.logit ( L );
                M.logit ( L );
                L->step();
        }
        L->finalize();
        //Exit program:

        delete L;
        return 0;
}