root/pmsm/TR2245/unitsteps.cpp @ 281

/*!
  \file
  \brief TR 2525 file for testing Toy Problem of mpf for Covariance Estimation
  \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 "../pmsm.h"
#include "simulator.h"
#include "../sim_profiles.h"

using namespace bdm;

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

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

        vec Qdiag;
        vec Rdiag;
        try {
                // Kalman filter
                F.lookupValue ( "ndat", Ndat );
                F.lookupValue ( "Npart",Npart );

                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"
        }
        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 ( zeros ( 4 ) ) );

        RV rQ ( "{Q }","4" );
        EKFCh_unQ KFEp ;
        KFEp.set_parameters ( &fxu,&hxu,Q,R );
        KFEp.set_est ( mu0, chmat ( zeros ( 4 ) ) );

        //mgamma_fix evolQ ( rQ,rQ );
        migamma_fix evolQ ;
        MPF<EKFCh_unQ> M ( &evolQ,&evolQ,Npart,KFEp );
        // initialize
        evolQ.set_parameters ( 0.1, 10*Qdiag, 1.0 ); //sigma = 1/10 mu
        evolQ.condition ( 10*Qdiag ); //Zdenek default
        M.set_est ( *evolQ._e() );
        evolQ.set_parameters ( 0.10, 10*Qdiag,0.999 ); //sigma = 1/10 mu
        //

        const epdf& KFEep = KFE.posterior();
        const epdf& Mep = M.posterior();

        dirfilelog *L; UIbuild(F.lookup("logger"), L);// ( "exp/mpf_test",100 );
        int l_X = L->add ( rx, "xt" );
        int l_D = L->add ( concat ( ry,ru ), "" );
        int l_XE= L->add ( rx, "xtE" );
        int l_XM= L->add ( concat ( rQ,rx ), "xtM" );
        int l_VE= L->add ( rx, "VE" );
        int l_VM= L->add ( concat ( rQ,rx ), "VM" );
        int l_Q= L->add ( rQ, "" );
        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 vecW=getvec(F.lookup("profile"));

        for ( int tK=1;tK<Ndat;tK++ ) {
                //Number of steps of a simulator for one step of Kalman
                for ( int ii=0; ii<Nsimstep;ii++ ) {
                        //simulator
                        sim_profile_vec01t ( Ww,vecW );
                        pmsmsim_step ( Ww );
                };
                ut ( 0 ) = KalmanObs[4];
                ut ( 1 ) = KalmanObs[5];
                xt = fxu.eval ( xtm,ut ) + diag ( sqrt ( Qdiag ) ) *randn ( 4 );
                dt = hxu.eval ( xt,ut );
                xtm = xt;

                //Variances
                if ( tK==1000 )  Qdiag ( 0 ) *=10;
                if ( tK==2000 ) Qdiag ( 0 ) /=10;
                if ( tK==3000 )  Qdiag ( 1 ) *=10;
                if ( tK==4000 ) Qdiag ( 1 ) /=10;
                if ( tK==5000 )  Qdiag ( 2 ) *=10;
                if ( tK==6000 ) Qdiag ( 2 ) /=10;
                if ( tK==7000 )  Qdiag ( 3 ) *=10;
                if ( tK==8000 ) Qdiag ( 3 ) /=10;

                //estimator
                KFE.bayes ( concat ( dt,ut ) );
                M.bayes ( concat ( dt,ut ) );

                vec mea=Mep.mean();
                if (max(mea)>1e3){
                        cout << "here"<<endl;
                }
                L->logit ( l_X,xt );
                L->logit ( l_D,concat ( dt,ut ) );
                L->logit ( l_XE,KFEep.mean() );
                L->logit ( l_XM, mea);
                L->logit ( l_VE,KFEep.variance() );
                L->logit ( l_VM,Mep.variance() );
                L->logit ( l_Q,Qdiag );
                L->step();
        }
        L->finalize();
        //Exit program:

        delete L;
        return 0;
}