root/tests/pmsm_sim.cpp @ 45

/*
  \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 <itpp/itbase.h>
#include <estim/libKF.h>
#include <estim/libPF.h>
#include <stat/libFN.h>

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

using namespace itpp;
/*
// PMSM with Q on Ia and Ib given externally
class EKF_unQ : public EKF<chmat> , public BMcond {
public:
        EKF_unQ( rx,ry,ru,rQ):EKF<chmat>(rx,ry,ru),BMcond(rQ){};
        void condition(const vec &Q0){};
};*/


int main() {
        // Kalman filter
        int Ndat = 10000;
        double h = 1e-6;
        int Nsimstep = 20;
        
        const it Nll = 10;

//      cout << KF;
        // 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.03 0.03 0.001 0.00001" );
        vec Rdiag ( "0.000001 0.000001" ); //var(diff(xth)) = "0.034 0.034"
        vec vQ = "0.01:0.01:0.1";
        chmat Q ( Qdiag );
        chmat R ( Rdiag );
        EKFCh KFE ( rx,ry,ru );
        KFE.set_est ( mu0, chmat ( 1000*ones ( 4 ) ) );
        KFE.set_parameters ( &fxu,&hxu,Q,R );

        mat ll ( Nll,Ndat );

        EKFCh* kfArray[Nll];


        for ( int i=0;i<Nll;i++ ) {
                vec Qid ( Qdiag );
                Qid ( 0 ) = vQ ( i ); Qid ( 1 ) = vQ ( i );
                kfArray[i]= new EKFCh ( rx,ry,ru );
                kfArray[i]->set_est ( mu0, chmat ( 100*ones ( 4 ) ) );
                kfArray[i]->set_parameters ( &fxu,&hxu,chmat ( Qid ),R );
        }

        epdf& KFEep = KFE._epdf();

        mat Xt=zeros ( 9,Ndat ); //true state from simulator
        mat Dt=zeros ( 4,Ndat ); //observation
        mat XtE=zeros ( 4,Ndat );

        // 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;
        static int k_rampa=1;
        static long k_rampa_tmp=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
                        Ww+=k_rampa*2.*M_PI*2e-4;    //1000Hz/s
                        if ( Ww>2.*M_PI*150. ) {
                                Ww=2.*M_PI*150.;
                                if ( k_rampa_tmp<500000 ) k_rampa_tmp++;
                                else {k_rampa=-1;k_rampa_tmp=0;}
                        };
                        if ( Ww<-2.*M_PI*150. ) Ww=-2.*M_PI*150.; /* */

                        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 ) );
                for ( int i=0;i<Nll;i++ ) {
                        kfArray[i]->bayes ( concat ( dt,ut ) );
                        ll ( i,tK ) =ll ( i,tK-1 ) + kfArray[i]->_ll();
                }

                Xt.set_col ( tK,vec ( x,9 ) );
                Dt.set_col ( tK, concat ( dt,ut ) );
                XtE.set_col ( tK,KFEep.mean() );
        }

        it_file fou ( "pmsm_sim.it" );

        fou << Name ( "xth" ) << Xt;
        fou << Name ( "Dt" ) << Dt;
        fou << Name ( "xthE" ) << XtE;
        fou << Name ( "ll" ) << ll;
        fou << Name ( "llgrid" ) << vQ;
        //Exit program:
        return 0;

}