Changeset 357 for applications
- Timestamp:
- 06/08/09 02:15:30 (16 years ago)
- Location:
- applications
- Files:
-
- 2 added
- 9 modified
Legend:
- Unmodified
- Added
- Removed
-
applications/doprava/k1.cpp
r341 r357 13 13 14 14 #include <estim/ekf_templ.h> 15 #include <stat/loggers_ui.h>16 #include <stat/libDS_ui.h>17 15 18 16 //include dopravni model 19 17 #include "model.h" 18 #include "stat\libDS.h" 19 #include "stat\loggers.h" 20 20 21 21 using namespace bdm; -
applications/mpdm/CMakeLists.txt
r299 r357 1 # Make sure the compiler can find include files from our Bdm library. 2 include_directories (${BDM_SOURCE_DIR}/bdm) 3 link_directories (${BDM_BINARY_DIR}/bdm) 1 # Compulsory line 2 cmake_minimum_required(VERSION 2.6) 4 3 5 EXEC(merg_pred) 6 EXEC(merg_giw) 7 EXEC(TR2244/merger_iter_cond) 8 EXEC(SYSID09/merg_2a) 4 add_subdirectory (SYSID09) 5 add_subdirectory (TR2244) 6 -
applications/pmsm/CMakeLists.txt
r284 r357 13 13 ## Build simulator 14 14 add_subdirectory (simulator_zdenek) 15 add_subdirectory (simulator_zdenek/ekf_example)15 #add_subdirectory (simulator_zdenek/ekf_example) 16 16 17 17 ## Use it -
applications/pmsm/TR2245/pmsm_wishart.cpp
r332 r357 17 17 #include <stat/libFN.h> 18 18 19 #include <stat/loggers_ui.h> 20 #include <stat/libEF_ui.h> 21 22 #include "../pmsm_ui.h" 19 #include "user_info.h" 20 #include "../pmsmDS.h" 23 21 24 22 using namespace bdm; … … 28 26 if ( argc>1 ) {fname = argv[1]; } 29 27 else { fname = "pmsm_wishart.cfg"; } 30 UIFile F ( fname ); 28 UI_File F ( fname ); 29 30 double h = 1e-6; 31 int Nsimstep = 125; 32 33 // Kalman filter 34 double k; 35 double l; 36 F.lookupValue ( "k",k); 37 F.lookupValue ( "l",l); 31 38 32 39 int Ndat; 33 40 int Npart; 34 double h = 1e-6;35 int Nsimstep = 125;41 F.lookupValue ( "ndat", Ndat ); 42 F.lookupValue ( "Npart", Npart ); 36 43 37 44 vec Qdiag; 38 45 vec Rdiag; 46 UI::get( Qdiag, F, "dQ" ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001 47 UI::get( Rdiag, F, "dR" );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034" 48 49 pmsmDS* DS = UI::build<pmsmDS>(F,"system"); 50 // mpdf* evolQ = UI::build<mpdf>( F, "Qrw" ); 39 51 40 pmsmDS* DS; 41 42 double k; 43 double l; 44 // mpdf* evolQ ; 45 try { 46 // Kalman filter 47 F.lookupValue ( "ndat", Ndat ); 48 F.lookupValue ( "Npart",Npart ); 49 50 F.lookupValue ( "k", k); 51 F.lookupValue ( "l",l); 52 53 // UIbuild ( F.lookup ( "Qrw" ),evolQ ); 54 Qdiag= getvec ( F.lookup ( "dQ" ) ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001 55 Rdiag=getvec ( F.lookup ( "dR" ) );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034" 56 57 UIbuild(F.lookup("system"),DS); 58 } 59 catch UICATCH; 60 // internal model 61 62 IMpmsm fxu; 63 // Rs Ls dt Fmag(Ypm) kp p J Bf(Mz) 64 fxu.set_parameters ( 0.28, 0.003465, Nsimstep*h, 0.1989, 1.5 ,4.0, 0.04, 0.0 ); 52 // internal model 53 IMpmsm fxu; 54 // Rs Ls dt Fmag(Ypm) kp p J Bf(Mz) 55 fxu.set_parameters ( 0.28, 0.003465, Nsimstep*h, 0.1989, 1.5 ,4.0, 0.04, 0.0 ); 65 56 // observation model 66 57 OMpmsm hxu; … … 91 82 M.set_rv ( concat ( rQ,rx ) ); 92 83 93 dirfilelog *L ; UIbuild ( F.lookup ( "logger" ), L);// ( "exp/mpf_test",100 );84 dirfilelog *L = UI::build<dirfilelog> ( F, "logger" );// ( "exp/mpf_test",100 ); 94 85 95 86 KFE.set_options ( "logbounds" ); -
applications/pmsm/TR2245/unitsteps.cpp
r317 r357 17 17 #include <stat/libFN.h> 18 18 19 #include <stat/loggers_ui.h>20 #include <stat/libEF_ui.h>21 19 22 20 #include "../pmsm.h" 23 21 #include "simulator.h" 24 22 #include "../sim_profiles.h" 23 #include "user_info.h" 24 #include "stat/loggers.h" 25 25 26 26 using namespace bdm; … … 30 30 if ( argc>1 ) {fname = argv[1]; } 31 31 else { fname = "unitsteps.cfg"; } 32 UI File F ( fname );32 UI_File F ( fname ); 33 33 34 int Ndat;35 int Npart;36 34 double h = 1e-6; 37 35 int Nsimstep = 125; 38 36 37 38 // Kalman filter 39 int Ndat; 40 int Npart; 41 F.lookupValue ( "ndat", Ndat ); 42 F.lookupValue ( "Npart",Npart ); 43 mpdf* evolQ = UI::build<mpdf>( F, "Qrw" ); 39 44 vec Qdiag; 40 45 vec Rdiag; 46 UI::get( Qdiag, F, "dQ" ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001 47 UI::get( Rdiag, F, "dR" );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034" 41 48 42 mpdf* evolQ ; 43 try { 44 // Kalman filter 45 F.lookupValue ( "ndat", Ndat ); 46 F.lookupValue ( "Npart",Npart ); 47 48 UIbuild ( F.lookup ( "Qrw" ),evolQ ); 49 Qdiag= getvec ( F.lookup ( "dQ" ) ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001 50 Rdiag=getvec ( F.lookup ( "dR" ) );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034" 51 } 52 catch UICATCH; 53 // internal model 49 // internal model 54 50 55 51 IMpmsm fxu; … … 81 77 M.set_rv ( concat ( rQ,rx ) ); 82 78 83 dirfilelog *L ; UIbuild ( F.lookup ( "logger" ), L);// ( "exp/mpf_test",100 );79 dirfilelog *L = UI::build<dirfilelog>( F, "logger" );// ( "exp/mpf_test",100 ); 84 80 int l_X = L->add ( rx, "xt" ); 85 81 int l_D = L->add ( concat ( ry,ru ), "" ); … … 99 95 vec xtm=zeros ( 4 ); 100 96 double Ww=0.0; 101 vec vecW=getvec ( F.lookup ( "profile" ) ); 97 vec vecW; 98 UI::get( vecW, F, "profile" ); 102 99 103 100 for ( int tK=1;tK<Ndat;tK++ ) { -
applications/pmsm/TR2245/wishart.cpp
r317 r357 17 17 #include <stat/libFN.h> 18 18 19 #include <stat/loggers_ui.h>20 #include <stat/libEF_ui.h>21 22 19 #include "../pmsm.h" 23 20 #include "simulator.h" 24 21 #include "../sim_profiles.h" 22 #include "user_info.h" 23 #include "stat/loggers.h" 25 24 26 25 using namespace bdm; … … 30 29 if ( argc>1 ) {fname = argv[1]; } 31 30 else { fname = "unitsteps.cfg"; } 32 UI File F ( fname );31 UI_File F ( fname ); 33 32 33 double h = 1e-6; 34 34 int Ndat; 35 35 int Npart; 36 double h = 1e-6; 36 F.lookupValue ( "ndat", Ndat ); 37 F.lookupValue ( "Npart",Npart ); 37 38 int Nsimstep = 125; 38 39 40 // Kalman filter 39 41 vec Qdiag; 42 UI::get( Qdiag, F, "dQ" ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001 43 40 44 vec Rdiag; 45 UI::get( Rdiag, F, "dR" );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034" 41 46 42 // mpdf* evolQ ; 43 try { 44 // Kalman filter 45 F.lookupValue ( "ndat", Ndat ); 46 F.lookupValue ( "Npart",Npart ); 47 48 // UIbuild ( F.lookup ( "Qrw" ),evolQ ); 49 Qdiag= getvec ( F.lookup ( "dQ" ) ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001 50 Rdiag=getvec ( F.lookup ( "dR" ) );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034" 51 } 52 catch UICATCH; 53 // internal model 54 55 IMpmsm fxu; 56 // Rs Ls dt Fmag(Ypm) kp p J Bf(Mz) 57 fxu.set_parameters ( 0.28, 0.003465, Nsimstep*h, 0.1989, 1.5 ,4.0, 0.04, 0.0 ); 47 // internal model 48 IMpmsm fxu; 49 // Rs Ls dt Fmag(Ypm) kp p J Bf(Mz) 50 fxu.set_parameters ( 0.28, 0.003465, Nsimstep*h, 0.1989, 1.5 ,4.0, 0.04, 0.0 ); 58 51 // observation model 59 52 OMpmsm hxu; … … 84 77 M.set_rv ( concat ( rQ,rx ) ); 85 78 86 dirfilelog *L ; UIbuild ( F.lookup ( "logger" ), L);// ( "exp/mpf_test",100 );79 dirfilelog *L = UI::build<dirfilelog> ( F, "logger" );// ( "exp/mpf_test",100 ); 87 80 int l_X = L->add ( rx, "xt" ); 88 81 int l_D = L->add ( concat ( ry,ru ), "" ); … … 103 96 vec xtm=zeros ( 4 ); 104 97 double Ww=0.0; 105 vec vecW=getvec ( F.lookup ( "profile" ) ); 98 vec vecW; 99 UI::get( vecW, F ,"profile" ); 106 100 107 101 mat tQ=diag(Qdiag); -
applications/pmsm/pmsm.h
r349 r357 3 3 4 4 #include <stat/libFN.h> 5 #include "user_info.h" 5 6 6 7 /*! \defgroup PMSM … … 94 95 95 96 void dfdu_cond ( const vec &x0, const vec &u0, mat &A, bool full=true ) {it_error ( "not needed" );}; 96 }; 97 97 98 void from_setting( const Setting &root ) 99 { 100 set_parameters ( root["params"]["Rs"], root["params"]["Ls"], 125e-6, root["params"]["Fmag"], \ 101 root["params"]["kp"], root["params"]["p"], root["params"]["J"], 0.0 ); 102 }; 103 104 // TODO dodelat void to_setting( Setting &root ) const; 105 }; 106 107 UIREGISTER ( IMpmsm ); 98 108 99 109 //! State evolution model for a PMSM drive and its derivative with respect to \f$x\f$ … … 182 192 }; 183 193 194 195 UIREGISTER ( IMpmsm2o ); 196 184 197 //! State evolution model for a PMSM drive and its derivative with respect to \f$x\f$, equation for \f$\omega\f$ is omitted.$ 185 198 class IMpmsmStat : public IMpmsm { … … 234 247 }; 235 248 249 UIREGISTER ( IMpmsmStat ); 250 251 236 252 //! State for PMSM with unknown Mz 237 253 class IMpmsmMz: public IMpmsm{ … … 256 272 }; 257 273 274 UIREGISTER ( IMpmsmMz ); 275 258 276 //! State for PMSM with unknown Mz 259 277 class IMpmsmStatMz: public IMpmsmStat{ … … 278 296 }; 279 297 298 UIREGISTER ( IMpmsmStatMz ); 299 280 300 281 301 //! Observation model for PMSM drive and its derivative with respect to \f$x\f$ … … 298 318 }; 299 319 320 UIREGISTER ( OMpmsm ); 321 300 322 //! Observation model for PMSM drive and its derivative with respect to \f$x\f$ for full vector of observations 301 323 class OMpmsm4: public diffbifn { … … 314 336 }; 315 337 338 UIREGISTER ( OMpmsm4 ); 339 340 341 342 343 316 344 /*!@}*/ 317 345 #endif //PMSM_H -
applications/pmsm/pmsmDS.h
r349 r357 17 17 18 18 //! Simulator of PMSM machine with predefined profile on omega 19 class pmsmDS : public DS { 19 class pmsmDS : public DS 20 { 20 21 21 22 protected: 22 23 24 25 26 27 28 29 30 31 32 33 34 //! options for logging, - 23 //! indeces of logged variables 24 int L_x, L_ou, L_oy, L_iu, L_optu; 25 //! Setpoints of omega in timespans given by dt_prof 26 vec profileWw; 27 //! Setpoints of Mz in timespans given by dt_prof 28 vec profileMz; 29 //! time-step for profiles 30 double dt_prof; 31 //! Number of miliseconds per discrete time step 32 int Dt; 33 //! options for logging, - log predictions of 'true' voltage 34 bool opt_modu; 35 //! options for logging, - 35 36 public: 36 //! Constructor with fixed sampling period 37 pmsmDS () {Dt=125; Drv=RV ( "{o_ua o_ub o_ia o_ib t_ua t_ub o_om o_th Mz }" );} 38 void set_parameters ( double Rs0, double Ls0, double Fmag0, double Bf0, double p0, double kp0, double J0, double Uc0, double DT0, double dt0 ) { 39 pmsmsim_set_parameters ( Rs0, Ls0, Fmag0, Bf0, p0, kp0, J0, Uc0, DT0, dt0 ); 40 } 41 //! parse options: "modelu" => opt_modu=true; 42 void set_options ( string &opt ) { 43 opt_modu = ( opt.find ( "modelu" ) !=string::npos ); 44 } 45 void getdata ( vec &dt ) {dt.set_subvector(0,vec ( KalmanObs,6 ));dt(6)=x[2];dt(7)=x[3];dt(8)=x[8];} 46 void write ( vec &ut ) {} 47 48 void step() { 49 static int ind=0; 50 static double dW; // increase of W 51 static double Ww; // W 52 static double Mz; // W 53 if ( t>=dt_prof*ind ) { 54 ind++; 55 // check omega profile and set dW 56 if ( ind<profileWw.length() ) { 57 //linear increase 58 if ( profileWw.length() ==1 ) { 59 Ww=profileWw ( 0 ); dW=0.0; 60 } 61 else { 62 dW = profileWw ( ind )-profileWw ( ind-1 ); 63 dW *=125e-6/dt_prof; 64 } 65 } 66 else { 67 dW = 0; 68 } 69 // Check Mz profile and set Mz 70 if ( ind<profileMz.length() ) { 71 //sudden increase 72 Mz = profileMz(ind); 73 } 74 else { 75 Mz = 0; 76 } 77 } 78 Ww += dW; 79 //Simulate Dt seconds! 80 for ( int i=0;i<Dt;i++ ) { pmsmsim_step ( Ww , Mz);} 37 //! Constructor with fixed sampling period 38 pmsmDS () 39 { 40 Dt=125; 41 Drv=RV ( "{o_ua o_ub o_ia o_ib t_ua t_ub o_om o_th Mz }" ); 42 } 43 void set_parameters ( double Rs0, double Ls0, double Fmag0, double Bf0, double p0, double kp0, double J0, double Uc0, double DT0, double dt0 ) 44 { 45 pmsmsim_set_parameters ( Rs0, Ls0, Fmag0, Bf0, p0, kp0, J0, Uc0, DT0, dt0 ); 46 } 47 //! parse options: "modelu" => opt_modu=true; 48 void set_options ( string &opt ) 49 { 50 opt_modu = ( opt.find ( "modelu" ) !=string::npos ); 51 } 52 void getdata ( vec &dt ) 53 { 54 dt.set_subvector(0,vec ( KalmanObs,6 )); 55 dt(6)=x[2]; 56 dt(7)=x[3]; 57 dt(8)=x[8]; 58 } 59 void write ( vec &ut ) {} 60 61 void step() 62 { 63 static int ind=0; 64 static double dW; // increase of W 65 static double Ww; // W 66 static double Mz; // W 67 if ( t>=dt_prof*ind ) 68 { 69 ind++; 70 // check omega profile and set dW 71 if ( ind<profileWw.length() ) 72 { 73 //linear increase 74 if ( profileWw.length() ==1 ) 75 { 76 Ww=profileWw ( 0 ); 77 dW=0.0; 78 } 79 else 80 { 81 dW = profileWw ( ind )-profileWw ( ind-1 ); 82 dW *=125e-6/dt_prof; 83 } 84 } 85 else 86 { 87 dW = 0; 88 } 89 // Check Mz profile and set Mz 90 if ( ind<profileMz.length() ) 91 { 92 //sudden increase 93 Mz = profileMz(ind); 94 } 95 else 96 { 97 Mz = 0; 98 } 99 } 100 Ww += dW; 101 //Simulate Dt seconds! 102 for ( int i=0; i<Dt; i++ ) 103 { 104 pmsmsim_step ( Ww , Mz); 105 } 81 106 // for ( int i=0;i<Dt;i++ ) { pmsmsim_noreg_step ( Ww , Mz);} 82 83 //discretization 84 double ustep=1.2; 85 KalmanObs [ 0 ] = ustep*itpp::round( KalmanObs [ 0 ]/ ustep) ; 86 KalmanObs [ 1 ] = ustep*itpp::round(KalmanObs [ 1 ]/ ustep); 87 double istep=0.085; 88 KalmanObs [ 2 ] = istep*itpp::round( KalmanObs [ 2 ]/ istep) ; 89 KalmanObs [ 3 ] = istep*itpp::round(KalmanObs [ 3 ]/ istep); 90 91 }; 92 93 void log_add ( logger &L ) { 94 L_x = L.add ( rx, "x" ); 95 L_oy = L.add ( ry, "o" ); 96 L_ou = L.add ( ru, "o" ); 97 L_iu = L.add ( ru, "t" ); 98 // log differences 99 if ( opt_modu ) { 100 L_optu = L.add ( ru, "model" ); 101 } 102 } 103 104 void logit ( logger &L ) { 105 L.logit ( L_x, vec ( x,4 ) ); 106 L.logit ( L_oy, vec_2 ( KalmanObs[2],KalmanObs[3] ) ); 107 L.logit ( L_ou, vec_2 ( KalmanObs[0],KalmanObs[1] ) ); 108 L.logit ( L_iu, vec_2 ( KalmanObs[4],KalmanObs[5] ) ); 109 if ( opt_modu ) { 110 double sq3=sqrt ( 3.0 ); 111 double ua,ub; 112 double i1=x[0]; 113 double i2=0.5* ( -i1+sq3*x[1] ); 114 double i3=0.5* ( -i1-sq3*x[1] ); 115 double u1=KalmanObs[0]; 116 double u2=0.5* ( -u1+sq3*KalmanObs[1] ); 117 double u3=0.5* ( -u1-sq3*KalmanObs[1] ); 118 119 double du1=1.4* ( double ( i1>0.3 ) - double ( i1<-0.3 ) ) +0.2*i1; 120 double du2=1.4* ( double ( i2>0.3 ) - double ( i2<-0.3 ) ) +0.2*i2; 121 double du3=1.4* ( double ( i3>0.3 ) - double ( i3<-0.3 ) ) +0.2*i3; 122 ua = ( 2.0* ( u1-du1 )- ( u2-du2 )- ( u3-du3 ) ) /3.0; 123 ub = ( ( u2-du2 )- ( u3-du3 ) ) /sq3; 124 L.logit ( L_optu , vec_2 ( ua,ub ) ); 125 } 126 127 } 128 129 void set_profile ( double dt, const vec &Ww, const vec &Mz ) {dt_prof=dt; profileWw=Ww; profileMz=Mz;} 107 108 //discretization 109 double ustep=1.2; 110 KalmanObs [ 0 ] = ustep*itpp::round( KalmanObs [ 0 ]/ ustep) ; 111 KalmanObs [ 1 ] = ustep*itpp::round(KalmanObs [ 1 ]/ ustep); 112 double istep=0.085; 113 KalmanObs [ 2 ] = istep*itpp::round( KalmanObs [ 2 ]/ istep) ; 114 KalmanObs [ 3 ] = istep*itpp::round(KalmanObs [ 3 ]/ istep); 115 116 }; 117 118 void log_add ( logger &L ) 119 { 120 L_x = L.add ( rx, "x" ); 121 L_oy = L.add ( ry, "o" ); 122 L_ou = L.add ( ru, "o" ); 123 L_iu = L.add ( ru, "t" ); 124 // log differences 125 if ( opt_modu ) 126 { 127 L_optu = L.add ( ru, "model" ); 128 } 129 } 130 131 void logit ( logger &L ) 132 { 133 L.logit ( L_x, vec ( x,4 ) ); 134 L.logit ( L_oy, vec_2 ( KalmanObs[2],KalmanObs[3] ) ); 135 L.logit ( L_ou, vec_2 ( KalmanObs[0],KalmanObs[1] ) ); 136 L.logit ( L_iu, vec_2 ( KalmanObs[4],KalmanObs[5] ) ); 137 if ( opt_modu ) 138 { 139 double sq3=sqrt ( 3.0 ); 140 double ua,ub; 141 double i1=x[0]; 142 double i2=0.5* ( -i1+sq3*x[1] ); 143 double i3=0.5* ( -i1-sq3*x[1] ); 144 double u1=KalmanObs[0]; 145 double u2=0.5* ( -u1+sq3*KalmanObs[1] ); 146 double u3=0.5* ( -u1-sq3*KalmanObs[1] ); 147 148 double du1=1.4* ( double ( i1>0.3 ) - double ( i1<-0.3 ) ) +0.2*i1; 149 double du2=1.4* ( double ( i2>0.3 ) - double ( i2<-0.3 ) ) +0.2*i2; 150 double du3=1.4* ( double ( i3>0.3 ) - double ( i3<-0.3 ) ) +0.2*i3; 151 ua = ( 2.0* ( u1-du1 )- ( u2-du2 )- ( u3-du3 ) ) /3.0; 152 ub = ( ( u2-du2 )- ( u3-du3 ) ) /sq3; 153 L.logit ( L_optu , vec_2 ( ua,ub ) ); 154 } 155 156 } 157 158 void set_profile ( double dt, const vec &Ww, const vec &Mz ) 159 { 160 dt_prof=dt; 161 profileWw=Ww; 162 profileMz=Mz; 163 } 164 165 void from_setting( const Setting &root ) 166 { 167 set_parameters ( root["params"]["Rs"], root["params"]["Ls"], root["params"]["Fmag"], \ 168 root["params"]["Bf"], root["params"]["p"], root["params"]["kp"], \ 169 root["params"]["J"], root["params"]["Uc"], root["params"]["DT"], 1.0e-6 ); 170 171 // Default values of profiles for omega and Mz 172 vec profW=vec("1.0"); 173 vec profM=vec("0.0"); 174 double tstep=1.0; 175 root.lookupValue( "tstep", tstep ); 176 UI::get( profW, root, "profileW" ); 177 UI::get( profM, root, "profileM" ); 178 set_profile (tstep , profW, profM); 179 180 string opts; 181 if ( root.lookupValue( "options", opts ) ) 182 set_options(opts); 183 } 184 185 // TODO dodelat void to_setting( Setting &root ) const; 130 186 }; 131 187 188 UIREGISTER ( pmsmDS ); 189 190 132 191 //! This class behaves like BM but it is evaluating EKF 133 class pmsmCRB : public EKFfull{ 134 protected: 135 vec interr; 136 vec old_true; 137 vec secder; 138 int L_CRB; 139 int L_err; 140 int L_sec; 141 public: 142 //! constructor 143 pmsmCRB():EKFfull(){old_true=zeros(6);} 144 145 void bayes(const vec &dt){ 146 static vec umin(2); 147 vec u(2); 148 //assume we know state exactly: 149 vec true_state=vec(x,4); // read from pmsm 150 E.set_mu(true_state); 151 mu=true_state; 152 153 //integration error 154 old_true(4)=KalmanObs[4]; 155 old_true(5)=KalmanObs[5];// add U 156 u(0) = KalmanObs[0]; // use the required value for derivatives 157 u(1) = KalmanObs[1]; 158 interr = (true_state - pfxu->eval(old_true)); 159 160 //second derivative 161 IMpmsm2o* pf = dynamic_cast<IMpmsm2o*>(pfxu); 162 if (pf) {secder=pf->eval2o(u-umin);} 163 164 umin =u; 165 EKFfull::bayes(dt); 166 old_true.set_subvector(0,true_state); 167 } 168 169 void log_add(logger &L, const string &name="" ){ 170 L_CRB=L.add(rx,"crb"); 171 L_err=L.add(rx,"err"); 172 L_sec=L.add(rx,"d2"); 173 } 174 void logit(logger &L){ 175 L.logit(L_err, interr); 176 L.logit(L_CRB,diag(_R())); 177 L.logit(L_sec,secder); 178 } 192 class pmsmCRB : public EKFfull 193 { 194 protected: 195 vec interr; 196 vec old_true; 197 vec secder; 198 int L_CRB; 199 int L_err; 200 int L_sec; 201 public: 202 //! constructor 203 pmsmCRB():EKFfull() 204 { 205 old_true=zeros(6); 206 } 207 208 void bayes(const vec &dt) 209 { 210 static vec umin(2); 211 vec u(2); 212 //assume we know state exactly: 213 vec true_state=vec(x,4); // read from pmsm 214 E.set_mu(true_state); 215 mu=true_state; 216 217 //integration error 218 old_true(4)=KalmanObs[4]; 219 old_true(5)=KalmanObs[5];// add U 220 u(0) = KalmanObs[0]; // use the required value for derivatives 221 u(1) = KalmanObs[1]; 222 interr = (true_state - pfxu->eval(old_true)); 223 224 //second derivative 225 IMpmsm2o* pf = dynamic_cast<IMpmsm2o*>(pfxu); 226 if (pf) 227 { 228 secder=pf->eval2o(u-umin); 229 } 230 231 umin =u; 232 EKFfull::bayes(dt); 233 old_true.set_subvector(0,true_state); 234 } 235 236 void log_add(logger &L, const string &name="" ) 237 { 238 L_CRB=L.add(rx,"crb"); 239 L_err=L.add(rx,"err"); 240 L_sec=L.add(rx,"d2"); 241 } 242 void logit(logger &L) 243 { 244 L.logit(L_err, interr); 245 L.logit(L_CRB,diag(_R())); 246 L.logit(L_sec,secder); 247 } 248 249 void from_setting( const Setting &root ) 250 { 251 diffbifn* IM = UI::build<diffbifn>(root, "IM"); 252 diffbifn* OM = UI::build<diffbifn>(root, "OM"); 253 254 //parameters 255 256 //statistics 257 int dim=IM->dimension(); 258 vec mu0; 259 if (!UI::get( mu0, root, "mu0") ) 260 mu0=zeros(dim); 261 mat P0; 262 vec dP0; 263 264 if (UI::get(dP0,root, "dP0") ) 265 P0=diag(dP0); 266 else if (!UI::get(P0,root, "P0") ) 267 P0=eye(dim); 268 269 set_statistics(mu0,P0); 270 271 vec dQ; 272 UI::get( dQ, root, "dQ"); 273 vec dR; 274 UI::get( dR, root, "dR"); 275 set_parameters(IM, OM, diag(dQ) , diag(dR)); 276 277 //connect 278 RV* drv = UI::build<RV>(root, "drv"); 279 set_drv(*drv); 280 RV* rv = UI::build<RV>(root, "rv"); 281 set_rv(*rv); 282 } 283 284 // TODO dodelat void to_setting( Setting &root ) const; 179 285 }; 180 286 287 UIREGISTER ( pmsmCRB ); 288 289 181 290 //! This class behaves like BM but it is evaluating EKF 182 class pmsmCRBMz : public EKFfull{ 183 protected: 184 int L_CRB; 185 public: 186 //! constructor 187 pmsmCRBMz():EKFfull(){} 188 189 void bayes(const vec &dt){ 291 class pmsmCRBMz : public EKFfull 292 { 293 protected: 294 int L_CRB; 295 public: 296 //! constructor 297 pmsmCRBMz():EKFfull() {} 298 299 void bayes(const vec &dt) 300 { 190 301 //assume we know state exactly: 191 vec true_state(5); 192 true_state.set_subvector(0,vec(x,4)); // read from pmsm 193 true_state(4)=x[8]; 194 195 E.set_mu(true_state); 196 mu = true_state; 197 //hack for ut 198 EKFfull::bayes(dt); 199 } 200 201 void log_add(logger &L, const string &name="" ){ 202 L_CRB=L.add(concat(rx,RV("Mz",1,0)),"crbz"); 203 } 204 void logit(logger &L){ 205 L.logit(L_CRB,diag(_R())); 206 } 302 vec true_state(5); 303 true_state.set_subvector(0,vec(x,4)); // read from pmsm 304 true_state(4)=x[8]; 305 306 E.set_mu(true_state); 307 mu = true_state; 308 //hack for ut 309 EKFfull::bayes(dt); 310 } 311 312 void log_add(logger &L, const string &name="" ) 313 { 314 L_CRB=L.add(concat(rx,RV("Mz",1,0)),"crbz"); 315 } 316 void logit(logger &L) 317 { 318 L.logit(L_CRB,diag(_R())); 319 } 320 321 void from_setting( const Setting &root ) 322 { 323 diffbifn* IM = UI::build<diffbifn>(root,"IM"); 324 diffbifn* OM = UI::build<diffbifn>(root,"OM"); 325 326 //statistics 327 int dim=IM->dimension(); 328 vec mu0; 329 vec dP0; 330 mat P0; 331 if (!UI::get(mu0, root, "mu0")) 332 mu0=zeros(dim); 333 if (!UI::get( P0, root, "P0" )) 334 if (UI::get(dP0, root, "dP0")) 335 P0=diag(dP0); 336 else 337 P0=eye(dim); 338 339 set_statistics(mu0,P0); 340 341 vec dQ; 342 UI::get(dQ, root, "dQ"); 343 vec dR; 344 UI::get(dR, root, "dR"); 345 set_parameters(IM, OM, diag(dQ), diag(dR)); 346 347 //connect 348 RV* drv = UI::build<RV>(root, "drv"); 349 set_drv(*drv); 350 RV* rv = UI::build<RV>(root, "rv"); 351 set_rv(*rv); 352 } 353 354 // TODO dodelat void to_setting( Setting &root ) const; 207 355 }; 356 357 UIREGISTER ( pmsmCRBMz ); -
applications/pmsm/pmsm_estim.cpp
r342 r357 7 7 */ 8 8 9 #include "pmsm_ui.h" 10 #include "stat/loggers_ui.h" 11 #include "estim/KF_ui.h" 12 #include "stat/libDS_ui.h" 9 #include "user_info.h" 10 #include "stat/loggers.h" 13 11 14 12 using namespace bdm; … … 17 15 if ( argc>1 ) {fname = argv[1]; } 18 16 else { cout << "Missing configuration file.\n Usage: \n $> estimator config_file.cfg"<<endl; abort(); } 19 UI File F ( fname );17 UI_File F ( fname ); 20 18 21 logger* L ;22 DS * DS;19 logger* L = UI::build <logger>( F, "logger" ); 20 DS * pDS = UI::build<DS> ( F, "system" ); 23 21 Array<BM*> Es; // array of estimators 22 UI::get( Es, F, "estimator" ); 23 int nE = Es.length(); //number of estimators 24 24 int Ndat; //number of data records 25 int nE; //number of estimators 26 27 try { 28 UIbuild ( F.lookup ( "logger" ),L ); 29 UIbuild ( F.lookup ( "system" ),DS ); 30 F.lookupValue ( "experiment.ndat",Ndat ); 31 Setting& S=F.lookup ( "estimator" ); 32 nE = S.getLength(); 33 Es.set_size(nE); 34 for(int i=0;i<nE;i++){ 35 UIbuild (S[i] ,Es(i) ); 36 } 37 } 38 catch UICATCH; 39 40 DS->log_add ( *L ); 25 F.lookupValue ( "experiment.ndat",Ndat ); 26 27 pDS->log_add ( *L ); 41 28 string nic=""; 42 29 for (int i=0; i<nE; i++){ … … 45 32 L->init(); 46 33 47 vec dt=zeros ( DS->_drv()._dsize() ); //data variable34 vec dt=zeros ( pDS->_drv()._dsize() ); //data variable 48 35 Array<datalink*> Dls(nE); 49 36 for (int i=0; i<nE; i++){ 50 Dls(i)=new datalink( Es(i)->_drv(), DS->_drv() ); //datalink between a datasource and estimator37 Dls(i)=new datalink( Es(i)->_drv(),pDS->_drv() ); //datalink between a datasource and estimator 51 38 } 52 39 … … 54 41 for ( int tK=1;tK<Ndat;tK++ ) { 55 42 // Data Source 56 DS->step(); // simulator step57 DS->getdata ( dt ); // read data58 DS->logit ( *L );43 pDS->step(); // simulator step 44 pDS->getdata ( dt ); // read data 45 pDS->logit ( *L ); 59 46 60 47 // Estimators … … 71 58 72 59 delete L; 73 delete DS;60 delete pDS; 74 61 for (int i=0; i<nE; i++){ 75 62 delete Dls(i);