root/pmsm/pmsm_sim3.cpp @ 78

Revision 72, 4.4 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 "../simulator_zdenek/ekf_example/ekf_obj.h"
22#include "iopom.h"
23
24using namespace itpp;
25//!Extended Kalman filter with unknown \c Q
26
27void set_simulator_t(double &Ww) {
28
29        if (t>0.0002) x[8]=1.2;    // 1A //0.2ZP
30        if (t>0.4) x[8]=10.8;   // 9A
31        if (t>0.6) x[8]=25.2;  // 21A
32
33        if (t>0.7) Ww=2.*M_PI*10.;
34        if (t>1.0) x[8]=1.2;    // 1A
35        if (t>1.2) x[8]=10.8;   // 9A
36        if (t>1.4) x[8]=25.2;  // 21A
37
38        if (t>1.6) Ww=2.*M_PI*50.;
39        if (t>1.9) x[8]=1.2;    // 1A
40        if (t>2.1) x[8]=10.8;   // 9A
41        if (t>2.3) x[8]=25.2;  // 21A
42
43        if (t>2.5) Ww=2.*M_PI*100;
44        if (t>2.8) x[8]=1.2;    // 1A
45        if (t>3.0) x[8]=10.8;   // 9A
46        if (t>3.2) x[8]=25.2;  // 21A
47
48        if (t>3.4) Ww=2.*M_PI*150;
49        if (t>3.7) x[8]=1.2;    // 1A
50        if (t>3.9) x[8]=10.8;   // 9A
51        if (t>4.1) x[8]=25.2;  // 21A
52
53        if (t>4.3) Ww=2.*M_PI*0;
54        if (t>4.8) x[8]=-1.2;    // 1A
55        if (t>5.0) x[8]=-10.8;   // 9A
56        if (t>5.2) x[8]=-25.2;  // 21A
57
58        if (t>5.4) Ww=2.*M_PI*(-10.);
59        if (t>5.7) x[8]=-1.2;    // 1A
60        if (t>5.9) x[8]=-10.8;   // 9A
61        if (t>6.1) x[8]=-25.2;  // 21A
62
63        if (t>6.3) Ww=2.*M_PI*(-50.);
64        if (t>6.7) x[8]=-1.2;    // 1A
65        if (t>6.9) x[8]=-10.8;   // 9A
66        if (t>7.1) x[8]=-25.2;  // 21A
67
68        if (t>7.3) Ww=2.*M_PI*(-100.);
69        if (t>7.7) x[8]=-1.2;    // 1A
70        if (t>7.9) x[8]=-10.8;   // 9A
71        if (t>8.1) x[8]=-25.2;  // 21A
72        if (t>8.3) x[8]=10.8;   // 9A
73        if (t>8.5) x[8]=25.2;  // 21A
74       
75        x[8]=0.0;
76}
77
78int main() {
79        // Kalman filter
80        int Ndat = 90000;
81        double h = 1e-6;
82        int Nsimstep = 125;
83        int Npart = 100;
84       
85
86        vec mu0= "0.0 0.0 0.0 0.0";
87        vec Qdiag ( "0.05 0.05 0.002 0.001" ); //zdenek: 0.01 0.01 0.0001 0.0001
88        vec Rdiag ( "0.05 0.05" ); //var(diff(xth)) = "0.034 0.034"
89        chmat Q ( Qdiag );
90        chmat R ( Rdiag );
91       
92        RV rQ ( "100","{Q}","4","0" );
93        EKFfixed KFE ( rx, rQ);
94        KFE.init_ekf ( Nsimstep*h);
95
96
97///////////// Particles
98        mgamma_fix evolQ ( rQ,rQ );
99        MPF<EKFfixed> M ( rx,rQ,evolQ,evolQ,Npart,KFE );
100        // initialize
101        evolQ.set_parameters ( 10.0 ,Qdiag, 1.0); //sigma = 1/10 mu
102        evolQ.condition ( Qdiag ); //Zdenek default
103        epdf& pfinit=evolQ._epdf();
104        M.set_est ( pfinit );
105        evolQ.set_parameters ( 500000.0, Qdiag, 0.9999 );
106
107
108        epdf& KFEep = KFE._epdf();
109        epdf& Mep = M._epdf();
110
111
112        mat Xt=zeros ( Ndat ,9 ); //true state from simulator
113        mat Dt=zeros ( Ndat,4 ); //observation
114        mat XtE=zeros ( Ndat, 4 );
115        mat XtM=zeros ( Ndat,4+4 ); //Q + x
116
117        // SET SIMULATOR
118        pmsmsim_set_parameters ( 0.28,0.003465,0.1989,0.0,4,1.5,0.04, 200., 3e-6, h );
119        double Ww=0.0;
120        static int k_rampa=1;
121        static long k_rampa_tmp=0;
122
123        vec dt ( 2 );
124        vec ut ( 2 );
125
126        for ( int tK=1;tK<Ndat;tK++ ) {
127                //Number of steps of a simulator for one step of Kalman
128                for ( int ii=0; ii<Nsimstep;ii++ ) {
129                        //simulator
130                        Ww+=k_rampa*2.*M_PI*2e-4;    //1000Hz/s
131                        if ( Ww>2.*M_PI*150. ) {
132                                Ww=2.*M_PI*150.;
133                                if ( k_rampa_tmp<500000 ) k_rampa_tmp++;
134                                else {k_rampa=-1;k_rampa_tmp=0;}
135                        };
136                        if ( Ww<-2.*M_PI*150. ) Ww=-2.*M_PI*150.; /* */
137//                      set_simulator_t(Ww);
138                        pmsmsim_step ( Ww );
139                };
140                // collect data
141                ut ( 0 ) = KalmanObs[0];
142                ut ( 1 ) = KalmanObs[1];
143                dt ( 0 ) = KalmanObs[2];
144                dt ( 1 ) = KalmanObs[3];
145                //estimator
146                KFE.bayes ( concat ( dt,ut ) );
147                M.bayes ( concat ( dt,ut ) );
148               
149                Xt.set_row ( tK,vec ( x,9 ) ); //vec from C-array
150                Dt.set_row ( tK, concat ( dt,ut));
151                XtE.set_row ( tK,KFEep.mean() );
152                double qNaN = std::numeric_limits<double>::quiet_NaN();
153                XtM.set_row ( tK,Mep.mean() );
154        }
155
156        char tmpstr[200];
157        sprintf(tmpstr,"%s/%s","herez/","format");
158        ofstream  form(tmpstr);
159        form << "# Experimetal header file"<< endl;
160        dirfile_write(form,"herez/",Xt,"X","{isa isb om th }" );
161        dirfile_write ( form,"herez",XtM,"XtM","{q1 q2 q3 q4 isa isb om th }" );
162        dirfile_write ( form,"herez",XtE,"XE","{isa isb om th }" );
163        dirfile_write ( form,"herez",Dt,"Dt","{isa isb ua ub }" );
164
165        return 0;
166}
Note: See TracBrowser for help on using the browser.