root/tests/pmsm_sim.cpp @ 54

/*
  \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"

#include <netcdfcpp.h>
void write_to_nc ( NcFile &nc, mat &X, std::string Xn, Array<std::string> A ) {
        char tmpstr[200];
        int Len = X.rows();

        sprintf ( tmpstr,"%s.length",Xn.c_str() );
        NcDim* lengt = nc.add_dim ( tmpstr, ( long ) Len );
        for ( int j=0; j<X.cols(); j++ ) {
                if ( j<A.length() )
                        sprintf ( tmpstr,"%s_%s",Xn.c_str(), ( A ( j ) ).c_str() );
                else
                        sprintf ( tmpstr,"%s_%d",Xn.c_str(),j );
                // Create variables and their attributes
                NcVar* P = nc.add_var ( tmpstr, ncDouble, lengt );
                const double* Dp = X._data();
                P->put ( &Dp[j*Len],Len );
        }
}


using namespace itpp;
//!Extended Kalman filter with unknown \c Q
class EKF_unQ : public EKFCh , public BMcond {
public:
        //! Default constructor
        EKF_unQ ( RV rx, RV ry,RV ru,RV rQ ) :EKFCh ( rx,ry,ru ),BMcond ( rQ ) {};
        void condition ( const vec &Q0 ) {
                Q.setD ( Q0,0 );
                //from EKF
                preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() );
        };
};


int main() {
        // Kalman filter
        int Ndat = 30000;
        double h = 1e-6;
        int Nsimstep = 125;
        int Npart = 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.00001 0.00001" ); //zdenek: 0.01 0.01 0.0001 0.0001
        vec Rdiag ( "0.0005 0.0005" ); //var(diff(xth)) = "0.034 0.034"
        chmat Q ( Qdiag );
        chmat R ( Rdiag );
        EKFCh KFE ( rx,ry,ru );
        KFE.set_parameters ( &fxu,&hxu,Q,R );
        KFE.set_est ( mu0, chmat ( 1*ones ( 4 ) ) );

        RV rQ ( "100","{Q}","2","0" );
        EKF_unQ KFEp ( rx,ry,ru,rQ );
        KFEp.set_parameters ( &fxu,&hxu,Q,R );
        KFEp.set_est ( mu0, chmat ( 1*ones ( 4 ) ) );

        mgamma evolQ ( rQ,rQ );
        MPF<EKF_unQ> M ( rx,rQ,evolQ,evolQ,Npart,KFEp );
        // initialize
        evolQ.set_parameters ( 100.0 ); //sigma = 1/10 mu
        evolQ.condition ( "0.01 0.01" ); //Zdenek default
        epdf& pfinit=evolQ._epdf();
        M.set_est ( pfinit );
        evolQ.set_parameters ( 1000.0 );

        //

        epdf& KFEep = KFE._epdf();
        epdf& Mep = M._epdf();

        mat Xt=zeros ( Ndat ,9 ); //true state from simulator
        mat Dt=zeros ( Ndat,4+2 ); //observation
        mat XtE=zeros ( Ndat, 4 );
        mat XtM=zeros ( Ndat,6 ); //Q + x

        // 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 ) );
                M.bayes ( concat ( dt,ut ) );

                Xt.set_row ( tK,vec ( x,9 ) ); //vec from C-array
                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))) ) );
                XtE.set_row ( tK,KFEep.mean() );
                XtM.set_row ( tK,Mep.mean() );
        }

        it_file fou ( "pmsm_sim.it" );

        fou << Name ( "xth" ) << Xt;
        fou << Name ( "Dt" ) << Dt;
        fou << Name ( "xthE" ) << XtE;
        fou << Name ( "xthM" ) << XtM;
        //Exit program:

        ////////////////
        // Just Testing
        NcFile nc ( "pmsm_sim.nc", NcFile::Replace ); // Create and leave in define mode
        if ( ! nc.is_valid() ) {        std::cerr << "creation of NCFile failed."<<endl;}

        write_to_nc ( nc,Xt,"X","{isa isb om th }" );
        write_to_nc ( nc,XtM,"XtM","{q1 q2 isa isb om th }" );
        write_to_nc ( nc,XtE,"XE","{isa isb om th }" );
        write_to_nc ( nc,Dt,"Dt","{isa isb ua ub }" );
        return 0;
}