root/pmsm/mpf_u_delta.cpp @ 232

/*!
  \file
  \brief TR 2525 file for testing Toy Problem of mpf for Covariance Estimation
  \author Vaclav Smidl.

  \ingroup PMSM

  -----------------------------------
  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 "sim_profiles.h"

using namespace itpp;

//!Extended Kalman filter with unknown \c Q and delta u
class EKFCh_du_kQ : public EKFCh , public BMcond {
        chmat Qref;
public:
        //! Default constructor
        EKFCh_du_kQ ( RV rx, RV ry,RV ru,RV rC ) :EKFCh ( rx,ry,ru ),BMcond ( rC ),Qref(rx.count()) {};
        void set_ref(const chmat &Qref0){Qref=Qref0;}
        void condition ( const vec &val ) {
                pfxu->condition ( val.left(2) );
                
                Q = Qref*std::abs(val(2));
                preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() );
        };
};

class IMpmsm_delta :  public IMpmsm {
protected:
        vec ud;
public:
        IMpmsm_delta() :IMpmsm(),ud ( 2 ) {ud=zeros ( 2 );};
        //! Set mechanical and electrical variables

        void condition ( const vec &val ) {ud = val;}
        vec eval ( const vec &x0, const vec &u0 ) {
                // last state
                double iam = x0 ( 0 );
                double ibm = x0 ( 1 );
                double omm = x0 ( 2 );
                double thm = x0 ( 3 );
                double uam = u0 ( 0 );
                double ubm = u0 ( 1 );

                vec xk=zeros ( 4 );
                //ia
                xk ( 0 ) = ( 1.0- Rs/Ls*dt ) * iam + Ypm/Ls*dt*omm * sin ( thm ) + ( uam+ud ( 0 ) ) *dt/Ls;
                //ib
                xk ( 1 ) = ( 1.0- Rs/Ls*dt ) * ibm - Ypm/Ls*dt*omm * cos ( thm ) + ( ubm+ud ( 1 ) ) *dt/Ls;
                //om
                xk ( 2 ) = omm + kp*p*p * Ypm/J*dt* ( ibm * cos ( thm )-iam * sin ( thm ) ) - p/J*dt*Mz;
                //th
                xk ( 3 ) = thm + omm*dt; // <0..2pi>
                if ( xk ( 3 ) >pi ) xk ( 3 )-=2*pi;
                if ( xk ( 3 ) <-pi ) xk ( 3 ) +=2*pi;
                return xk;
        }

};


int main() {
        // Kalman filter
        int Ndat = 9000;
        double h = 1e-6;
        int Nsimstep = 125;
        int Npart = 20;

        mat Rnoise = randn ( 2,Ndat );

        // internal model
        IMpmsm2o fxu0;
        IMpmsm_delta 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 );
        fxu0.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.4 0.4 0.001 0.000001" ); //zdenek: 0.01 0.01 0.0001 0.0001
        vec Rdiag ( "0.1 0.1" ); //var(diff(xth)) = "0.034 0.034"
        chmat Q ( Qdiag );
        mat Q2o=diag(Qdiag);
        chmat R ( Rdiag );
        EKFCh KFE ( rx,ry,rU );
        KFE.set_parameters ( &fxu0,&hxu,Q,R );
        KFE.set_est ( mu0, chmat ( ones ( 4 ) ) );

        RV rUd ( "{ud k}", "2 1" );
        EKFCh_du_kQ KFEp ( rx,ry,ru,rUd );
        KFEp.set_parameters ( &fxu,&hxu,Q,R );
        KFEp.set_ref(Q);
        KFEp.set_est ( mu0, chmat ( ones ( 4 ) ) );

        mlnorm<ldmat> evolUd ( rUd,rUd );
        MPF<EKFCh_du_kQ> M ( rx,rUd,evolUd,evolUd,Npart,KFEp );
        // initialize
        vec Ud0="0 0 100.0";
        evolUd.set_parameters ( eye ( 3 ), zeros(3), ldmat ( vec( "1 1 10" )));
        evolUd.condition ( Ud0 );
        epdf& pfinit=evolUd._epdf();
        M.set_est ( pfinit );
        evolUd.set_parameters ( eye ( 3 ), zeros(3), ldmat ( vec(" 9e-8 9e-8 9e-7" )));

        mat Xt=zeros ( Ndat ,4 ); //true state from simulator
        mat Dt=zeros ( Ndat,2+2+2 ); //observation
        mat XtE=zeros ( Ndat, 4 );
        mat Qtr=zeros ( Ndat, 4 );
        mat XtM=zeros ( Ndat, 3+4 ); //W + x
        mat XtMTh=zeros ( Ndat,1 );

        // 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 ( 4 );
        vec utm ( 4 );
        vec xt ( 4 );
        vec xtm=zeros ( 4 );
        double Ww=0.0;
        vec vecW="1 2 4 8 4 2 0 -4 -9 -16 -4 0 0 0";
        vecW*=10.0;

        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[0];
                ut ( 1 ) = KalmanObs[1];
                ut ( 2 ) = utm(0)-ut(0);
                ut ( 3 ) = utm(1)-ut(1);
                dt ( 0 ) = KalmanObs[2]+0.3*Rnoise ( 0,t );
                dt ( 1 ) = KalmanObs[3]+0.3*Rnoise ( 1,t );
                xt = vec ( x,4 );

                utm =ut;
                
                Q2o*=0.9;
                Q2o+=outer_product(x2o_dbg,x2o_dbg);
                //estimator
                KFE.bayes ( concat ( dt,ut ) );
                for (int ii=0; ii<Npart; ii++){dynamic_cast<EKFCh_du_kQ*>(M._BM(ii))->set_ref(chmat(Q2o));}
                M.bayes ( concat ( dt,ut ) );

                Xt.set_row ( tK, xt ); //vec from C-array
                Dt.set_row ( tK, concat ( dt,ut.left(2), vec_2 ( KalmanObs[4],KalmanObs[5] ) ) );
                Qtr.set_row ( tK, Qdiag );
                XtE.set_row ( tK,KFE._e()->mean() );
                XtM.set_row ( tK,M._e()->mean() );
                // correction for theta

                {
                        double sumSin=0.0;
                        double sumCos=0.0;
                        vec mea ( 4 );
                        vec* _w;

                        for ( int p=0; p<Npart;p++ ) {
                                mea = M._BM ( p )->_e()->mean();
                                _w = M.__w();
                                sumSin += ( *_w ) ( p ) *sin ( mea ( 3 ) );
                                sumCos += ( *_w ) ( p ) *cos ( mea ( 3 ) );
                        }
                        double Th = atan2 ( sumSin, sumCos );

                        XtMTh.set_row ( tK,vec_1 ( Th ) );
                }

        }

        it_file fou ( "mpf_u_delta.it" );

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

        return 0;
}