root/pmsm/old/pmsm_sim2.cpp @ 303

/*
  \file
  \brief Models for synchronous electric drive using IT++ and BDM
  \author Vaclav Smidl.

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

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


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

#include <stat/loggers.h>

#include "pmsm.h"
#include "simulator.h"
#include "sim_profiles.h"

using namespace bdm;
//!Extended Kalman filter with unknown \c Q
class EKF_unQ : public EKFCh , public BMcond {
public:
        //! Default constructor
        void condition ( const vec &Q0 ) {
                Q.setD ( Q0,0 );
                //from EKF
                preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() );
        };
        void bayes(const vec dt){
        EKFCh::bayes(dt);
                
                vec xtrue(4);
                //UGLY HACK!!! reliance on a predictor!!
                xtrue(0)=x[0];
                xtrue(1)=x[1];
                xtrue(2)=x[2];
                xtrue(3)=x[3];
                
                ll = -0.5* ( 4 * 1.83787706640935 +_P.logdet() +_P.invqform(xtrue));
        }
};
class EKF_unQful : public EKFfull , public BMcond {
public:
        void condition ( const vec &Q0 ) {
                Q=diag(Q0);
        };
        void bayes(const vec dt){
        EKFfull::bayes(dt);
                
                vec xtrue(4);
                //UGLY HACK!!! reliance on a predictor!!
                xtrue(0)=x[0];
                xtrue(1)=x[1];
                xtrue(2)=x[2];
                xtrue(3)=x[3];
                
                BM::ll = -0.5* ( 4 * 1.83787706640935 +log(det(P)) +xtrue* ( inv(P)*xtrue ) );
        }
};

int main() {
        // Kalman filter
        int Ndat = 90000;
        double h = 1e-6;
        int Nsimstep = 125;
        int Npart = 50;
        
        dirfilelog L("exp/pmsm_sim2",1000);
        
        // 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";
//      vec Qdiag ( "0.01 0.01 0.01 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001
        vec Qdiag ( "10 10 10 0.001" ); //zdenek: 0.01 0.01 0.0001 0.0001
        vec Rdiag ( "100 100" ); //var(diff(xth)) = "0.034 0.034"
        chmat Q ( Qdiag );
        chmat R ( Rdiag );
        EKFCh KFE ;
        KFE.set_est ( mu0, chmat( 1*eye ( 4 ) ) );
        KFE.set_parameters ( &fxu,&hxu,Q,R);

        RV rQ ( "{Q}","4" );
        EKF_unQful KFEp ;
        KFEp.set_est ( mu0,  1*ones ( 4 ) );
        KFEp.set_parameters ( &fxu,&hxu,diag(Qdiag),diag(Rdiag) );

        mgamma_fix evolQ ;
        MPF<EKF_unQful> M ( &evolQ,&evolQ,Npart,KFEp );
        // initialize
        evolQ.set_parameters ( 1000.0 ,Qdiag, 0.5); //sigma = 1/10 mu
        evolQ.condition ( Qdiag ); //Zdenek default
        epdf& pfinit=evolQ._epdf();
        M.set_est ( pfinit );
        evolQ.set_parameters ( 100000.0, Qdiag, 0.9999 );
        //
        const epdf& KFEep = KFE.posterior();
        const epdf& Mep = M.posterior();

        int X_log = L.add(rx,"X");
        int Efix_log = L.add(rx,"XF");
        RV tmp=concat(rQ,rx);
        int M_log = L.add(tmp,"M");
        L.init();

        // SET SIMULATOR
        pmsmsim_set_parameters ( 0.28,0.003465,0.1989,0.0,4,1.5,0.04, 200., 3e-6, h );
        double Ww=0.0;
        vec dt ( 2 );
        vec ut ( 2 );

        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_steps1(Ww);
                        pmsmsim_step ( Ww );
                };
                // collect data
                ut ( 0 ) = KalmanObs[0];
                ut ( 1 ) = KalmanObs[1];
                dt ( 0 ) = KalmanObs[2];
                dt ( 1 ) = KalmanObs[3];

                //estimator
                KFE.bayes ( concat ( dt,ut ) );
                M.bayes ( concat ( dt,ut ) );
                
                L.logit(X_log,  vec(x,4)); //vec from C-array
                L.logit(Efix_log, KFEep.mean() ); 
                L.logit(M_log,  Mep.mean() ); 
                
                L.step();
        }

        L.finalize(); 
        
        return 0;
}