root/pmsm/pmsm_sim2.cpp @ 78

Revision 72, 6.2 kB (checked in by smidl, 18 years ago)

presun pmsm

  • Property svn:eol-style set to native
Line 
1/*
2  \file
3  \brief Models for synchronous electric drive using IT++ and BDM
4  \author Vaclav Smidl.
5
6  -----------------------------------
7  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
8
9  Using IT++ for numerical operations
10  -----------------------------------
11*/
12
13#include <itpp/itbase.h>
14#include <estim/libKF.h>
15#include <estim/libPF.h>
16#include <stat/libFN.h>
17
18#include "pmsm.h"
19#include "simulator.h"
20
21#include "iopom.h"
22
23using namespace itpp;
24//!Extended Kalman filter with unknown \c Q
25class EKF_unQ : public EKFCh , public BMcond {
26public:
27        //! Default constructor
28        EKF_unQ ( RV rx, RV ry,RV ru,RV rQ ) :EKFCh ( rx,ry,ru ),BMcond ( rQ ) {};
29        void condition ( const vec &Q0 ) {
30                Q.setD ( Q0,0 );
31                //from EKF
32                preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() );
33        };
34        void bayes(const vec dt){
35        EKFCh::bayes(dt);
36               
37                vec xtrue(4);
38                //UGLY HACK!!! reliance on a predictor!!
39                xtrue(0)=x[0];
40                xtrue(1)=x[1];
41                xtrue(2)=x[2];
42                xtrue(3)=x[3];
43               
44                ll = -0.5* ( 4 * 1.83787706640935 +_P->logdet() +xtrue* ( inv(_P->to_mat())*xtrue ) );
45        }
46};
47class EKF_unQful : public EKFfull , public BMcond {
48public:
49        //! Default constructor
50        EKF_unQful ( RV rx, RV ry,RV ru,RV rQ ) :EKFfull ( rx,ry,ru ),BMcond ( rQ ) {};
51        void condition ( const vec &Q0 ) {
52                Q=diag(Q0);
53        };
54        void bayes(const vec dt){
55        EKFfull::bayes(dt);
56               
57                vec xtrue(4);
58                //UGLY HACK!!! reliance on a predictor!!
59                xtrue(0)=x[0];
60                xtrue(1)=x[1];
61                xtrue(2)=x[2];
62                xtrue(3)=x[3];
63               
64                BM::ll = -0.5* ( 4 * 1.83787706640935 +log(det(P)) +xtrue* ( inv(P)*xtrue ) );
65        }
66};
67
68void set_simulator_t(double &Ww) {
69
70        if (t>0.2) x[8]=1.2;    // 1A //0.2ZP
71        if (t>0.4) x[8]=10.8;   // 9A
72        if (t>0.6) x[8]=25.2;  // 21A
73
74        if (t>0.7) Ww=2.*M_PI*10.;
75        if (t>1.0) x[8]=1.2;    // 1A
76        if (t>1.2) x[8]=10.8;   // 9A
77        if (t>1.4) x[8]=25.2;  // 21A
78
79        if (t>1.6) Ww=2.*M_PI*50.;
80        if (t>1.9) x[8]=1.2;    // 1A
81        if (t>2.1) x[8]=10.8;   // 9A
82        if (t>2.3) x[8]=25.2;  // 21A
83
84        if (t>2.5) Ww=2.*M_PI*100;
85        if (t>2.8) x[8]=1.2;    // 1A
86        if (t>3.0) x[8]=10.8;   // 9A
87        if (t>3.2) x[8]=25.2;  // 21A
88
89        if (t>3.4) Ww=2.*M_PI*150;
90        if (t>3.7) x[8]=1.2;    // 1A
91        if (t>3.9) x[8]=10.8;   // 9A
92        if (t>4.1) x[8]=25.2;  // 21A
93
94        if (t>4.3) Ww=2.*M_PI*0;
95        if (t>4.8) x[8]=-1.2;    // 1A
96        if (t>5.0) x[8]=-10.8;   // 9A
97        if (t>5.2) x[8]=-25.2;  // 21A
98
99        if (t>5.4) Ww=2.*M_PI*(-10.);
100        if (t>5.7) x[8]=-1.2;    // 1A
101        if (t>5.9) x[8]=-10.8;   // 9A
102        if (t>6.1) x[8]=-25.2;  // 21A
103
104        if (t>6.3) Ww=2.*M_PI*(-50.);
105        if (t>6.7) x[8]=-1.2;    // 1A
106        if (t>6.9) x[8]=-10.8;   // 9A
107        if (t>7.1) x[8]=-25.2;  // 21A
108
109        if (t>7.3) Ww=2.*M_PI*(-100.);
110        if (t>7.7) x[8]=-1.2;    // 1A
111        if (t>7.9) x[8]=-10.8;   // 9A
112        if (t>8.1) x[8]=-25.2;  // 21A
113        if (t>8.3) x[8]=10.8;   // 9A
114        if (t>8.5) x[8]=25.2;  // 21A
115       
116        if (t>9) Ww=2.*M_PI*0;
117        x[8]=0.0;
118}
119
120int main() {
121        // Kalman filter
122        int Ndat = 90000;
123        double h = 1e-6;
124        int Nsimstep = 125;
125        int Npart = 50;
126       
127        //StrSim:06:
128        vec SSAT(Ndat);
129
130        // internal model
131        IMpmsm fxu;
132        //                  Rs    Ls        dt           Fmag(Ypm) kp   p    J     Bf(Mz)
133        fxu.set_parameters ( 0.28, 0.003465, Nsimstep*h, 0.1989,   1.5 ,4.0, 0.04, 0.0 );
134        // observation model
135        OMpmsm hxu;
136
137        vec mu0= "0.0 0.0 0.0 0.0";
138//      vec Qdiag ( "0.01 0.01 0.01 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001
139        vec Qdiag ( "10 10 10 0.001" ); //zdenek: 0.01 0.01 0.0001 0.0001
140        vec Rdiag ( "100 100" ); //var(diff(xth)) = "0.034 0.034"
141        chmat Q ( Qdiag );
142        chmat R ( Rdiag );
143        EKFCh KFE ( rx,ry,ru );
144        KFE.set_est ( mu0, chmat( 1*eye ( 4 ) ) );
145        KFE.set_parameters ( &fxu,&hxu,Q,R);
146
147        RV rQ ( "100","{Q}","4","0" );
148        EKF_unQful KFEp ( rx,ry,ru,rQ );
149        KFEp.set_est ( mu0,  1*ones ( 4 ) );
150        KFEp.set_parameters ( &fxu,&hxu,diag(Qdiag),diag(Rdiag) );
151
152        mgamma_fix evolQ ( rQ,rQ );
153        MPF<EKF_unQful> M ( rx,rQ,evolQ,evolQ,Npart,KFEp );
154        // initialize
155        evolQ.set_parameters ( 1000.0 ,Qdiag, 0.5); //sigma = 1/10 mu
156        evolQ.condition ( Qdiag ); //Zdenek default
157        epdf& pfinit=evolQ._epdf();
158        M.set_est ( pfinit );
159        evolQ.set_parameters ( 100000.0, Qdiag, 0.9999 );
160
161        //
162
163        epdf& KFEep = KFE._epdf();
164        epdf& Mep = M._epdf();
165
166        mat Xt=zeros ( Ndat ,9 ); //true state from simulator
167        mat Dt=zeros ( Ndat,4+2 ); //observation
168        mat XtE=zeros ( Ndat, 4 );
169        mat XtM=zeros ( Ndat,4+4 ); //Q + x
170
171        // SET SIMULATOR
172        pmsmsim_set_parameters ( 0.28,0.003465,0.1989,0.0,4,1.5,0.04, 200., 3e-6, h );
173        double Ww=0.0;
174        vec dt ( 2 );
175        vec ut ( 2 );
176
177        for ( int tK=1;tK<Ndat;tK++ ) {
178                //Number of steps of a simulator for one step of Kalman
179                for ( int ii=0; ii<Nsimstep;ii++ ) {
180                        //simulator
181                        set_simulator_t(Ww);
182                        pmsmsim_step ( Ww );
183                };
184                // collect data
185                ut ( 0 ) = KalmanObs[0];
186                ut ( 1 ) = KalmanObs[1];
187                dt ( 0 ) = KalmanObs[2];
188                dt ( 1 ) = KalmanObs[3];
189
190                //estimator
191                KFE.bayes ( concat ( dt,ut ) );
192                M.bayes ( concat ( dt,ut ) );
193                SSAT(tK) = M.SSAT;
194               
195                Xt.set_row ( tK,vec ( x,9 ) ); //vec from C-array
196                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))) ) );
197                XtE.set_row ( tK,KFEep.mean() );
198                XtM.set_row ( tK,Mep.mean() );
199        }
200
201/*      it_file fou ( "pmsm_sim.it" );
202
203        fou << Name ( "xth" ) << Xt;
204        fou << Name ( "Dt" ) << Dt;
205        fou << Name ( "xthE" ) << XtE;
206        fou << Name ( "xthM" ) << XtM;
207        fou << Name ( "SSAT" ) << SSAT;*/
208        //Exit program:
209
210        char tmpstr[200];
211        sprintf(tmpstr,"%s/%s","here/","format");
212        ofstream  form(tmpstr);
213        form << "# Experimetal header file"<< endl;
214        dirfile_write(form,"here/",Xt,"X","{isa isb om th }" );
215        dirfile_write ( form,"here",XtM,"XtM","{q1 q2 q3 q4 isa isb om th }" );
216        dirfile_write ( form,"here",XtE,"XE","{isa isb om th }" );
217        dirfile_write ( form,"here",Dt,"Dt","{isa isb ua ub }" );
218
219        ////////////////
220        // Just Testing
221/*      NcFile nc ( "pmsm_sim2.nc", NcFile::Replace ); // Create and leave in define mode
222        if ( ! nc.is_valid() ) {        std::cerr << "creation of NCFile failed."<<endl;}
223
224        write_to_nc ( nc,Xt,"X","{isa isb om th }" );
225        write_to_nc ( nc,XtM,"XtM","{q1 q2 q3 q4 isa isb om th }" );
226        write_to_nc ( nc,XtE,"XE","{isa isb om th }" );
227        write_to_nc ( nc,Dt,"Dt","{isa isb ua ub }" );*/
228        return 0;
229}
Note: See TracBrowser for help on using the browser.