- Timestamp:
- 05/27/08 10:59:54 (17 years ago)
- Location:
- pmsm
- Files:
-
- 1 removed
- 4 modified
-
CMakeLists.txt (modified) (2 diffs)
-
iopom.h (deleted)
-
pmsm_mix.cpp (modified) (4 diffs)
-
pmsm_sim2.cpp (modified) (6 diffs)
-
sim_var.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
pmsm/CMakeLists.txt
r105 r117 2 2 include_directories (${BDM_SOURCE_DIR}/bdm) 3 3 4 # Make sure the linker can find the Hellolibrary once it is built.4 # Make sure the linker can find the bdm library once it is built. 5 5 link_directories (${BDM_BINARY_DIR}/bdm) 6 6 … … 34 34 add_executable (pmsm_sim pmsm_sim.cpp) 35 35 add_executable (pmsm_sim2 pmsm_sim2.cpp) 36 add_executable (pmsm_sim3 pmsm_sim3.cpp)37 36 target_link_libraries (pmsm_sim ${BdmLibs} pmsmsim ) 38 37 target_link_libraries (pmsm_sim2 ${BdmLibs} pmsmsim ) 39 target_link_libraries (pmsm_sim3 ${BdmLibs} pmsmsim ekf_obj)40 38 41 39 add_executable (sim_var sim_var.cpp) -
pmsm/pmsm_mix.cpp
r115 r117 14 14 #include <stat/libFN.h> 15 15 #include <stat/emix.h> 16 #include <estim/ libKF.h>16 #include <estim/ekf_templ.h> 17 17 #include <estim/libPF.h> 18 18 19 19 #include "pmsm.h" 20 20 #include "simulator.h" 21 #include "sim_profiles.h" 21 22 22 23 #include <stat/loggers.h> 23 24 24 25 using namespace itpp; 25 //!Extended Kalman filter with unknown \c Q26 27 void set_simulator_t ( double &Ww ) {28 29 if ( t>0.0002 ) x[8]=1.2; // 1A //0.2ZP30 if ( t>0.4 ) x[8]=10.8; // 9A31 if ( t>0.6 ) x[8]=25.2; // 21A32 33 if ( t>0.7 ) Ww=2.*M_PI*10.;34 if ( t>1.0 ) x[8]=1.2; // 1A35 if ( t>1.2 ) x[8]=10.8; // 9A36 if ( t>1.4 ) x[8]=25.2; // 21A37 38 if ( t>1.6 ) Ww=2.*M_PI*50.;39 if ( t>1.9 ) x[8]=1.2; // 1A40 if ( t>2.1 ) x[8]=10.8; // 9A41 if ( t>2.3 ) x[8]=25.2; // 21A42 43 if ( t>2.5 ) Ww=2.*M_PI*100;44 if ( t>2.8 ) x[8]=1.2; // 1A45 if ( t>3.0 ) x[8]=10.8; // 9A46 if ( t>3.2 ) x[8]=25.2; // 21A47 48 if ( t>3.4 ) Ww=2.*M_PI*150;49 if ( t>3.7 ) x[8]=1.2; // 1A50 if ( t>3.9 ) x[8]=10.8; // 9A51 if ( t>4.1 ) x[8]=25.2; // 21A52 53 if ( t>4.3 ) Ww=2.*M_PI*0;54 if ( t>4.8 ) x[8]=-1.2; // 1A55 if ( t>5.0 ) x[8]=-10.8; // 9A56 if ( t>5.2 ) x[8]=-25.2; // 21A57 58 if ( t>5.4 ) Ww=2.*M_PI* ( -10. );59 if ( t>5.7 ) x[8]=-1.2; // 1A60 if ( t>5.9 ) x[8]=-10.8; // 9A61 if ( t>6.1 ) x[8]=-25.2; // 21A62 63 if ( t>6.3 ) Ww=2.*M_PI* ( -50. );64 if ( t>6.7 ) x[8]=-1.2; // 1A65 if ( t>6.9 ) x[8]=-10.8; // 9A66 if ( t>7.1 ) x[8]=-25.2; // 21A67 68 if ( t>7.3 ) Ww=2.*M_PI* ( -100. );69 if ( t>7.7 ) x[8]=-1.2; // 1A70 if ( t>7.9 ) x[8]=-10.8; // 9A71 if ( t>8.1 ) x[8]=-25.2; // 21A72 if ( t>8.3 ) x[8]=10.8; // 9A73 if ( t>8.5 ) x[8]=25.2; // 21A74 75 x[8]=0.0;76 }77 78 //!Extended Kalman filter with unknown \c Q79 class EKFful_unQR : public EKFfull , public BMcond {80 public:81 //! Default constructor82 EKFful_unQR ( RV rx, RV ry,RV ru,RV rQR ) :EKFfull ( rx,ry,ru ),BMcond ( rQR ) {};83 void condition ( const vec &Q0 ) {84 Q=diag(Q0(0,3));85 R=diag(Q0(4,5));86 };87 /* void bayes(const vec dt){88 EKFfull::bayes(dt);89 90 vec xtrue(4);91 //UGLY HACK!!! reliance on a predictor!!92 xtrue(0)=x[0];93 xtrue(1)=x[1];94 xtrue(2)=x[2];95 xtrue(3)=x[3];96 97 BM::ll = -0.5* ( 4 * 1.83787706640935 +log(det(P)) +xtrue* ( inv(P)*xtrue ) );98 }*/99 };100 26 101 27 int main() { … … 104 30 double h = 1e-6; 105 31 int Nsimstep = 125; 106 int Npar = 100;32 int Npar = 2000; 107 33 108 34 dirfilelog L("exp/pmsm_mix",1000); … … 142 68 //QU model 143 69 egamma Gcom(rQR);Gcom.set_parameters(ones(6),vec("1 1 1e4 1e10 1 1")); 144 cout << Gcom.mean() <<endl;145 cout << Gcom.sample() <<endl; 146 euni Ucom(rQR); Ucom.set_parameters(zeros(6),vec("60 60 0.03 1e-8100 100"));147 cout << Ucom.mean() <<endl;148 cout << Ucom.sample() <<endl; 70 /* cout << Gcom.mean() <<endl; 71 cout << Gcom.sample() <<endl;*/ 72 euni Ucom(rQR); Ucom.set_parameters(zeros(6),vec("60 60 453 0.03 100 100")); 73 /* cout << Ucom.mean() <<endl; 74 cout << Ucom.sample() <<endl;*/ 149 75 Array<epdf*> Coms(2); 150 76 Coms(0) = &Gcom; 151 77 Coms(1) = &Ucom; 152 emix Eevol(rQR); Eevol.set_parameters("0. 5 0.5", Coms);153 cout << Eevol.sample() <<endl;78 emix Eevol(rQR); Eevol.set_parameters("0.1 0.9", Coms); 79 // cout << Eevol.sample() <<endl; 154 80 155 81 mepdf evolQR(rQR,rQR,&Eevol); … … 163 89 int X_log = L.add(rx,"X"); 164 90 int Efix_log = L.add(rx,"XF"); 165 int M_log = L.add(concat(r x,rQR),"M");91 int M_log = L.add(concat(rQR,rx),"M"); 166 92 L.init(); 93 94 double dum=0.0; 95 vec dumvec = vec_1(1.0); 96 vec z= evolQR.samplecond(dumvec,dum) ; 97 cout << z << endl; 167 98 168 99 for ( int tK=1;tK<Ndat;tK++ ) { 169 100 //Number of steps of a simulator for one step of Kalman 170 101 for ( int ii=0; ii<Nsimstep;ii++ ) { 171 s et_simulator_t ( Ww);102 sim_profile_steps1 ( Ww, true ); 172 103 pmsmsim_step ( Ww ); 173 104 }; -
pmsm/pmsm_sim2.cpp
r81 r117 16 16 #include <stat/libFN.h> 17 17 18 #include <stat/loggers.h> 19 18 20 #include "pmsm.h" 19 21 #include "simulator.h" 20 21 #include "iopom.h" 22 #include "sim_profiles.h" 22 23 23 24 using namespace itpp; … … 66 67 }; 67 68 68 void set_simulator_t(double &Ww) {69 70 if (t>0.2) x[8]=1.2; // 1A //0.2ZP71 if (t>0.4) x[8]=10.8; // 9A72 if (t>0.6) x[8]=25.2; // 21A73 74 if (t>0.7) Ww=2.*M_PI*10.;75 if (t>1.0) x[8]=1.2; // 1A76 if (t>1.2) x[8]=10.8; // 9A77 if (t>1.4) x[8]=25.2; // 21A78 79 if (t>1.6) Ww=2.*M_PI*50.;80 if (t>1.9) x[8]=1.2; // 1A81 if (t>2.1) x[8]=10.8; // 9A82 if (t>2.3) x[8]=25.2; // 21A83 84 if (t>2.5) Ww=2.*M_PI*100;85 if (t>2.8) x[8]=1.2; // 1A86 if (t>3.0) x[8]=10.8; // 9A87 if (t>3.2) x[8]=25.2; // 21A88 89 if (t>3.4) Ww=2.*M_PI*150;90 if (t>3.7) x[8]=1.2; // 1A91 if (t>3.9) x[8]=10.8; // 9A92 if (t>4.1) x[8]=25.2; // 21A93 94 if (t>4.3) Ww=2.*M_PI*0;95 if (t>4.8) x[8]=-1.2; // 1A96 if (t>5.0) x[8]=-10.8; // 9A97 if (t>5.2) x[8]=-25.2; // 21A98 99 if (t>5.4) Ww=2.*M_PI*(-10.);100 if (t>5.7) x[8]=-1.2; // 1A101 if (t>5.9) x[8]=-10.8; // 9A102 if (t>6.1) x[8]=-25.2; // 21A103 104 if (t>6.3) Ww=2.*M_PI*(-50.);105 if (t>6.7) x[8]=-1.2; // 1A106 if (t>6.9) x[8]=-10.8; // 9A107 if (t>7.1) x[8]=-25.2; // 21A108 109 if (t>7.3) Ww=2.*M_PI*(-100.);110 if (t>7.7) x[8]=-1.2; // 1A111 if (t>7.9) x[8]=-10.8; // 9A112 if (t>8.1) x[8]=-25.2; // 21A113 if (t>8.3) x[8]=10.8; // 9A114 if (t>8.5) x[8]=25.2; // 21A115 116 if (t>9) Ww=2.*M_PI*0;117 x[8]=0.0;118 }119 120 69 int main() { 121 70 // Kalman filter … … 125 74 int Npart = 50; 126 75 127 //StrSim:06: 128 vec SSAT(Ndat); 129 76 dirfilelog L("exp/pmsm_sim2",1000); 77 130 78 // internal model 131 79 IMpmsm fxu; … … 158 106 M.set_est ( pfinit ); 159 107 evolQ.set_parameters ( 100000.0, Qdiag, 0.9999 ); 160 161 108 // 162 163 109 epdf& KFEep = KFE._epdf(); 164 110 epdf& Mep = M._epdf(); 165 111 166 mat Xt=zeros ( Ndat ,9 ); //true state from simulator167 mat Dt=zeros ( Ndat,4+2 ); //observation168 mat XtE=zeros ( Ndat, 4);169 mat XtM=zeros ( Ndat,4+4 ); //Q + x112 int X_log = L.add(rx,"X"); 113 int Efix_log = L.add(rx,"XF"); 114 int M_log = L.add(concat(rQ,rx),"M"); 115 L.init(); 170 116 171 117 // SET SIMULATOR … … 179 125 for ( int ii=0; ii<Nsimstep;ii++ ) { 180 126 //simulator 181 s et_simulator_t(Ww);127 sim_profile_steps1(Ww); 182 128 pmsmsim_step ( Ww ); 183 129 }; … … 191 137 KFE.bayes ( concat ( dt,ut ) ); 192 138 M.bayes ( concat ( dt,ut ) ); 193 SSAT(tK) = M.SSAT;194 139 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() ); 140 L.logit(X_log, vec(x,4)); //vec from C-array 141 L.logit(Efix_log, KFEep.mean() ); 142 L.logit(M_log, Mep.mean() ); 143 144 L.step(false); 199 145 } 200 146 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 }" );*/ 147 L.step(true); //final 148 228 149 return 0; 229 150 } -
pmsm/sim_var.cpp
r108 r117 19 19 #include "pmsm.h" 20 20 #include "simulator.h" 21 #include "sim_profiles.h" 21 22 22 23 #include <stat/loggers.h> … … 24 25 using namespace itpp; 25 26 //!Extended Kalman filter with unknown \c Q 26 27 void set_simulator_t ( double &Ww ) {28 29 if ( t>0.0002 ) x[8]=1.2; // 1A //0.2ZP30 if ( t>0.4 ) x[8]=10.8; // 9A31 if ( t>0.6 ) x[8]=25.2; // 21A32 33 if ( t>0.7 ) Ww=2.*M_PI*10.;34 if ( t>1.0 ) x[8]=1.2; // 1A35 if ( t>1.2 ) x[8]=10.8; // 9A36 if ( t>1.4 ) x[8]=25.2; // 21A37 38 if ( t>1.6 ) Ww=2.*M_PI*50.;39 if ( t>1.9 ) x[8]=1.2; // 1A40 if ( t>2.1 ) x[8]=10.8; // 9A41 if ( t>2.3 ) x[8]=25.2; // 21A42 43 if ( t>2.5 ) Ww=2.*M_PI*100;44 if ( t>2.8 ) x[8]=1.2; // 1A45 if ( t>3.0 ) x[8]=10.8; // 9A46 if ( t>3.2 ) x[8]=25.2; // 21A47 48 if ( t>3.4 ) Ww=2.*M_PI*150;49 if ( t>3.7 ) x[8]=1.2; // 1A50 if ( t>3.9 ) x[8]=10.8; // 9A51 if ( t>4.1 ) x[8]=25.2; // 21A52 53 if ( t>4.3 ) Ww=2.*M_PI*0;54 if ( t>4.8 ) x[8]=-1.2; // 1A55 if ( t>5.0 ) x[8]=-10.8; // 9A56 if ( t>5.2 ) x[8]=-25.2; // 21A57 58 if ( t>5.4 ) Ww=2.*M_PI* ( -10. );59 if ( t>5.7 ) x[8]=-1.2; // 1A60 if ( t>5.9 ) x[8]=-10.8; // 9A61 if ( t>6.1 ) x[8]=-25.2; // 21A62 63 if ( t>6.3 ) Ww=2.*M_PI* ( -50. );64 if ( t>6.7 ) x[8]=-1.2; // 1A65 if ( t>6.9 ) x[8]=-10.8; // 9A66 if ( t>7.1 ) x[8]=-25.2; // 21A67 68 if ( t>7.3 ) Ww=2.*M_PI* ( -100. );69 if ( t>7.7 ) x[8]=-1.2; // 1A70 if ( t>7.9 ) x[8]=-10.8; // 9A71 if ( t>8.1 ) x[8]=-25.2; // 21A72 if ( t>8.3 ) x[8]=10.8; // 9A73 if ( t>8.5 ) x[8]=25.2; // 21A74 75 x[8]=0.0;76 }77 27 78 28 int main() { … … 139 89 //Number of steps of a simulator for one step of Kalman 140 90 for ( int ii=0; ii<Nsimstep;ii++ ) { 141 s et_simulator_t ( Ww);91 sim_profile_steps1 ( Ww , true); 142 92 pmsmsim_step ( Ww ); 143 93 };
