Changeset 723
- Timestamp:
- 11/15/09 23:02:02 (15 years ago)
- Files:
-
- 3 added
- 19 modified
-
applications/bdmtoolbox/mex/controlloop.cpp (modified) (5 diffs)
-
applications/bdmtoolbox/mex/estimator.cpp (modified) (1 diff)
-
applications/dual/experiment/itermc.m (modified) (1 diff)
-
applications/dual/src/arx1_ctrl.h (modified) (3 diffs)
-
applications/dual/src/iterativemc.cpp (modified) (4 diffs)
-
library/bdm/base/bdmbase.h (modified) (3 diffs)
-
library/bdm/base/datasources.h (modified) (2 diffs)
-
library/bdm/design/arx2statespace_notes.lyx (added)
-
library/bdm/design/arx_ctrl.h (added)
-
library/bdm/design/ctrlbase.cpp (modified) (2 diffs)
-
library/bdm/design/ctrlbase.h (modified) (1 diff)
-
library/bdm/estim/arx.cpp (modified) (1 diff)
-
library/bdm/estim/arx.h (modified) (2 diffs)
-
library/bdm/estim/kalman.h (modified) (7 diffs)
-
library/bdm/itpp_ext.cpp (modified) (1 diff)
-
library/bdm/itpp_ext.h (modified) (1 diff)
-
library/bdm/math/chmat.h (modified) (1 diff)
-
library/bdm/stat/exp_family.h (modified) (1 diff)
-
library/tests/testsuite/LQG_test.cpp (modified) (1 diff)
-
library/tests/testsuite/datalink_test.cpp (modified) (1 diff)
-
library/tests/testsuite/dirfile-format.matrix (added)
-
library/tests/testsuite/logger_test.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
applications/bdmtoolbox/mex/controlloop.cpp
r706 r723 103 103 104 104 //DBG 105 Cfg.writeFile ( "c losedloop.cfg" );105 Cfg.writeFile ( "controlloop.cfg" ); 106 106 107 107 #else … … 111 111 fname = argv[1]; 112 112 } else { 113 fname="c losedloop.cfg";113 fname="controlloop.cfg"; 114 114 } 115 115 UIFile Cfg ( fname ); … … 119 119 Array<shared_ptr<Controller> > Cs; 120 120 UI::get ( Cs,Cfg, "controllers" ); 121 int Ndat=10 ;121 int Ndat=100; 122 122 if ( Cfg.exists ( "experiment" ) ) { 123 123 if ( Cfg.getRoot().lookupValue ( "experiment.ndat",Ndat ) ) { … … 187 187 Cs(i) -> adapt( Dlsc(i) ->pushdown(dt)); 188 188 vec ut=Cs ( i )->ctrlaction ( Dlsc(i) ->pushdown(dt) ); // update estimates 189 Ds->write(ut, vec_1(1)); 189 Ds->write(ut, vec_1(0)); 190 } else { 191 Ds->write(0.001*randn(Ds->_urv()._dsize())); 190 192 } 193 191 194 Cs ( i )->log_write (); 192 195 } … … 196 199 //update buffered fdat links 197 200 for (int i=0; i<Dls_buf.length(); i++){ 198 Dls_buf(i)->st ep(dt);201 Dls_buf(i)->store_data(dt); 199 202 } 200 203 -
applications/bdmtoolbox/mex/estimator.cpp
r706 r723 195 195 //update buffered fdat links 196 196 for (int i=0; i<Dls_buf.length(); i++){ 197 Dls_buf(i)->st ep(dt);197 Dls_buf(i)->store_data(dt); 198 198 } 199 199 -
applications/dual/experiment/itermc.m
r706 r723 25 25 c.controller = C1; 26 26 Mmc=iterativemc(c); 27 losses(i) = sum((Mmc. y-yr*ones(size(Mmc.y))).^2);27 losses(i) = sum((Mmc.S_y-yr*ones(size(Mmc.S_y))).^2); 28 28 29 29 c.controller = C2; 30 30 Mmc=iterativemc(c); 31 losses2(i) = sum((Mmc.y-yr*ones(size(Mmc.y))).^2); 31 looo= sum((Mmc.S_y-yr*ones(size(Mmc.S_y))).^2) 32 % keyboard 33 losses2(i) = looo; 32 34 end 33 35 [min(losses) median(losses) max(losses)] -
applications/dual/src/arx1_ctrl.h
r706 r723 21 21 void set_b (const double b0){b=b0;} 22 22 //! data in ctrlaction are only past outputs yt 23 vec ctrlaction(const vec &yt) {23 vec ctrlaction(const vec &yt)const{ 24 24 vec ut(1); 25 25 if (abs(b)<1e-2) { … … 63 63 public: 64 64 //! expected input is yt 65 virtual void adapt(const vec & yt){65 virtual void adapt(const vec &dt){ 66 66 //prepare data for ARX 67 Psi(0) = yt(0)-ytm;68 Psi(1) = utm;69 67 // do estimation 70 est.bayes(Psi); 68 est.bayes(vec_1(dt(0)-ytm),vec_1(dt(1))); 69 ytm = dt(0); 71 70 // save estimate b = E(b); 72 71 b = est.posterior().mean()(0); // first 73 72 } 74 73 //! same as exact control, but stores ut as utm 75 vec ctrlaction(const vec &yt){ 76 vec ut=exact_ctrl::ctrlaction(yt); 74 vec ctrlaction(const vec &dt)const{ 75 vec ut(1); 76 mat P = est.posterior().variance(); 77 ut(0) =(yr - dt(0))/(b + (b<0 ? P(0,0) : -P(0,0)));; 77 78 // remember last ut for estimation 78 utm=ut(0);79 ytm=yt(0);80 79 return ut; 81 80 } … … 92 91 V0(1,0) = b0; 93 92 V0(0,1) = b0; 94 V0(1,1) = P0 ;93 V0(1,1) = P0/6; 95 94 96 est.set_statistics(1,V0 );95 est.set_statistics(1,V0,6); 97 96 est.set_constant(false); 97 est.validate(); 98 98 validate(); 99 99 } 100 100 void validate(){ 101 101 Psi = zeros(2); 102 LIDs = zeros_i(1);103 102 } 104 virtual void log_add ( logger &L, const string &name = "" ) { 105 LIDs ( 0 ) = L.add ( RV("bhat",1) ); 103 void log_register( logger& L, const string& prefix ){ 104 est.set_log_level(1); 105 est.log_register(L,""); 106 }; 107 void log_write() const{ 108 est.log_write(); 106 109 } 107 virtual void logit ( logger &L ) {108 L.logit ( LIDs ( 0 ), b );109 }110 111 110 }; 112 111 UIREGISTER(ce_ctrl); -
applications/dual/src/iterativemc.cpp
r706 r723 47 47 double b; 48 48 double sigma; 49 double ytm ;49 double ytm=0.0; 50 50 double yr; 51 51 int Ndat; … … 91 91 RV y("y",1); // name random variable 92 92 RV u("u",1); // name random variable 93 int L_yt=L->add(y );94 int L_ut=L->add(u );95 C->log_ add(*L);93 int L_yt=L->add(y,"S"); 94 int L_ut=L->add(u,"S"); 95 C->log_register(*L,"C"); 96 96 L->init(); 97 97 98 98 vec psi(2); // regressor 99 double ut ; // input99 double ut=0.0; // input 100 100 vec yt; 101 101 … … 103 103 //yt has now meaning of yt-1 104 104 //design controller 105 C->adapt(vec_ 1(ytm));106 ut = C->ctrlaction(vec_ 1(ytm))(0); // the best controller there is105 C->adapt(vec_2(ytm,ut)); 106 ut = C->ctrlaction(vec_2(ytm,ut))(0); // the best controller there is 107 107 108 108 //prepare regressor … … 116 116 L->logit(L_yt, yt); 117 117 L->logit(L_ut, vec_1(ut)); 118 C->log it(*L);118 C->log_write(); 119 119 L->step(); 120 120 } -
library/bdm/base/bdmbase.h
r713 r723 796 796 datalink_buffered() : datalink_part(), history ( 0 ), h2v_down ( 0 ), h2v_hist ( 0 ) {}; 797 797 //! push current data to history 798 void st ep( const vec &val_up ) {798 void store_data ( const vec &val_up ) { 799 799 if ( v2h_up.length() > 0 ) { 800 800 history.shift_right ( 0, v2h_up.length() ); … … 818 818 void set_connection ( const RV &rv, const RV &rv_up ) { 819 819 // create link between up and down 820 datalink_part::set_connection ( rv, rv_up );820 datalink_part::set_connection ( rv, rv_up.remove_time() ); // only non-delayed version 821 821 822 822 // create rvs of history … … 871 871 //! update buffer 872 872 void step ( const vec &dt, const vec &ut ) { 873 dl1.st ep( dt );874 dl2.st ep( ut );873 dl1.store_data ( dt ); 874 dl2.store_data ( ut ); 875 875 } 876 876 }; -
library/bdm/base/datasources.h
r700 r723 163 163 public: 164 164 void step() { 165 yt2rgr.st ep(yt); // y is now history165 yt2rgr.store_data(yt); // y is now history 166 166 ut2rgr.filldown ( ut,rgr ); 167 167 yt2rgr.filldown ( yt,rgr ); 168 168 yt=ipdf->samplecond ( rgr ); 169 ut2rgr.st ep(ut); //u is now history169 ut2rgr.store_data(ut); //u is now history 170 170 } 171 171 void getdata ( vec &dt_out ) const { … … 175 175 } 176 176 void write(const vec &ut0){ut=ut0;} 177 177 void write(const vec &ut0, const ivec &ind){set_subvector(ut,ind,ut0);} 178 178 179 /*! 179 180 \code -
library/bdm/design/ctrlbase.cpp
r706 r723 3 3 namespace bdm{ 4 4 5 void LQG::set_system(shared_ptr<StateSpace< fsqmat> > S0){5 void LQG::set_system(shared_ptr<StateSpace<chmat> > S0){ 6 6 S = S0; 7 dimx= S->_A().rows();8 dimy= S->_C().rows();9 dimu= S->_B().cols();10 pr=zeros(dimx+dimu+dimy, dimu+dimx+dimu+dimy);11 pr.set_submatrix(dimx, dimu+dimx, eye(dimu+dimy));12 7 } 13 8 … … 17 12 y_req=y_req0; 18 13 horizon = horizon0; 19 validate();20 initialize();21 14 } 22 15 23 16 void LQG::validate() { 17 // set properties from S 18 dimx= S->_A().rows(); 19 dimy= S->_C().rows(); 20 dimu= S->_B().cols(); 21 pr=zeros(dimx+dimu+dimy, dimu+dimx+dimu+dimy); 22 pr.set_submatrix(dimx, dimu+dimx, eye(dimu+dimy)); 23 24 //penalization 24 25 bdm_assert ( Qy.cols() == dimy, "LQG: wrong dimensions of Qy " ); 25 26 bdm_assert ( Qu.cols() == dimu, "LQG: wrong dimensions of Qu " ); 26 27 bdm_assert ( y_req.length() == dimy, "LQG: wrong dimensions of y_req " ); 28 29 qux=zeros(dimu,dimx+2*dimu+dimy); 30 qyx=zeros(dimy, dimx+dimy+dimu); 31 32 // 33 initial_belmann(); 34 // parts of QR 35 post_qr = zeros(pre_qr.rows(), pre_qr.cols()); 36 37 update_system(); 27 38 } 28 39 29 void LQG::initialize(){30 // set parameter matrix31 pr.set_submatrix(0,0,S->_B());32 pr.set_submatrix(0,dimu, S->_A());33 34 //penalization35 qux=zeros(dimu,dimx+2*dimu+dimy);36 qux.set_submatrix(0,0,Qu);37 qux.set_submatrix(0,dimx+dimu+dimy,Qu);38 39 qyx=zeros(dimy, dimx+dimy+dimu);40 qyx.set_submatrix(0,0,S->_C());41 qyx.set_submatrix(0,dimx,-eye(dimy));42 43 //44 s=1e-5*eye(dimx+dimu+dimy);45 // parts of QR46 hqy=Qy*qyx*pr;47 48 // pre_qr49 pre_qr = concat_vertical(s*pr, concat_vertical(hqy, qux));50 post_qr = zeros(pre_qr.rows(), pre_qr.cols());51 }52 40 53 41 } -
library/bdm/design/ctrlbase.h
r706 r723 14 14 #include "../estim/kalman.h" 15 15 16 namespace bdm {16 namespace bdm { 17 17 18 //! Base class for adaptive controllers 19 20 21 22 23 24 25 //! identifier of the conditioning variables;26 27 28 29 virtual void adapt(const vec &data){};30 //! function redesigning the control strategy31 virtual void redesign(){};32 33 virtual vec ctrlaction(const vec &cond){return vec(0);}34 35 void from_setting(const Setting &set){36 UI::get(rv,set,"rv",UI::optional);37 UI::get(rvc,set,"rvc",UI::optional);38 39 40 41 42 43 44 virtual void log_register (logger &L, const string &prefix ) { }45 46 47 48 49 18 //! Base class for adaptive controllers 19 //! The base class is, however, non-adaptive, method \c adapt() is empty. 20 //! \note advanced Controllers will probably include estimator as their internal attribute (e.g. dual controllers) 21 class Controller : public root { 22 protected: 23 //! identifier of the designed action; 24 RV rv; 25 //! identifier of the conditioning variables - data needed ; 26 RV rvc; 27 public: 28 //! function processing new observations and adapting control strategy accordingly 29 virtual void adapt ( const vec &data ) {}; 30 //! function redesigning the control strategy 31 virtual void redesign() {}; 32 //! returns designed control action 33 virtual vec ctrlaction ( const vec &cond ) const {return vec ( 0 );} 34 35 void from_setting ( const Setting &set ) { 36 UI::get ( rv,set,"rv",UI::optional ); 37 UI::get ( rvc,set,"rvc",UI::optional ); 38 } 39 //! access function 40 const RV& _rv() {return rv;} 41 //! access function 42 const RV& _rvc() {return rvc;} 43 //! register this controller with given datasource under name "name" 44 virtual void log_register ( logger &L, const string &prefix ) { } 45 //! write requested values into the logger 46 virtual void log_write ( ) const { } 47 48 }; 49 50 50 //! Linear Quadratic Gaussian designer for constant penalizations and constant target 51 51 //! Its internals are very close to Kalman estimator 52 52 class LQG : public Controller { 53 protected: 54 //! StateSpace model from which we read data 55 shared_ptr<StateSpace<fsqmat> > S; 56 //! required value of the output y at time t (assumed constant) 57 vec y_req; 58 //! required value of the output y at time t (assumed constant) 59 vec u_req; 60 61 //! Control horizon, set to maxint for infinite horizons 62 int horizon; 63 //! penalization matrix Qy 64 mat Qy; 65 //! penalization matrix Qu 66 mat Qu; 67 //! time of the design step - from horizon->0 68 int td; 69 //! controller parameters 70 mat L; 71 72 //!@{ \name temporary storage for ricatti - use initialize 73 //! convenience parameters 74 int dimx; 75 //! convenience parameters 76 int dimy; 77 //! convenience parameters 78 int dimu; 79 80 //! parameters 81 mat pr; 82 //! penalization 83 mat qux; 84 //! penalization 85 mat qyx; 86 //! internal quadratic form 87 mat s; 88 //! penalization 89 mat qy; 90 //! pre_qr part 91 mat hqy; 92 //! pre qr matrix 93 mat pre_qr; 94 //! post qr matrix 95 mat post_qr; 96 //!@} 97 98 public: 99 //! set system parameters from given matrices 100 void set_system(shared_ptr<StateSpace<fsqmat> > S0); 101 //! set penalization matrices and control horizon 102 void set_control_parameters(const mat &Qy0, const mat &Qu0, const vec &y_req0, int horizon0); 103 //! set system parameters from Kalman filter 104 // void set_system_parameters(const Kalman &K); 105 //! refresh temporary storage - inefficient can be improved 106 void initialize(); 107 //! validation procedure 108 void validate(); 109 //! function for future use which is called at each time td; Should call initialize()! 110 virtual void update_state(){}; 111 //! redesign one step of the 112 void ricatti_step(){ 113 pre_qr.set_submatrix(0,0,s*pr); 114 pre_qr.set_submatrix(dimx+dimu+dimy, dimu+dimx, -Qy*y_req); 115 if (!qr(pre_qr,post_qr)){ bdm_warning("QR in LQG unstable");} 116 triu(post_qr); 117 // hn(m+1:2*m+n+r,m+1:2*m+n+r); 118 s=post_qr.get(dimu, 2*dimu+dimx+dimy-1, dimu, 2*dimu+dimx+dimy-1); 119 }; 120 void redesign(){ 121 for(td=horizon; td>0; td--){ 122 update_state(); 123 ricatti_step(); 124 } 125 /* ws=hn(1:m,m+1:2*m+n+r); 126 wsd=hn(1:m,1:m); 127 Lklq=-inv(wsd)*ws;*/ 128 L = -inv(post_qr.get(0,dimu-1, 0,dimu-1)) * post_qr.get(0,dimu-1, dimu, 2*dimu+dimx+dimy-1); 53 protected: 54 //! StateSpace model from which we read data 55 shared_ptr<StateSpace<chmat> > S; 56 //! required value of the output y at time t (assumed constant) 57 vec y_req; 58 //! required value of the output y at time t (assumed constant) 59 vec u_req; 60 61 //! Control horizon, set to maxint for infinite horizons 62 int horizon; 63 //! penalization matrix Qy 64 chmat Qy; 65 //! penalization matrix Qu 66 chmat Qu; 67 //! time of the design step - from horizon->0 68 int td; 69 //! controller parameters 70 mat L; 71 72 //!@{ \name temporary storage for ricatti - use initialize 73 //! convenience parameters 74 int dimx; 75 //! convenience parameters 76 int dimy; 77 //! convenience parameters 78 int dimu; 79 80 //! parameters 81 mat pr; 82 //! penalization 83 mat qux; 84 //! penalization 85 mat qyx; 86 //! internal quadratic form 87 mat s; 88 //! penalization 89 mat qy; 90 //! pre_qr part 91 mat hqy; 92 //! pre qr matrix 93 mat pre_qr; 94 //! post qr matrix 95 mat post_qr; 96 //!@} 97 98 public: 99 //! set system parameters from given matrices 100 void set_system ( shared_ptr<StateSpace<chmat> > S0 ); 101 //! update internal whan system has changed 102 void update_system() { 103 pr.set_submatrix ( 0,0,S->_B() ); 104 pr.set_submatrix ( 0,dimu, S->_A() ); 105 106 //penalization 107 qux.set_submatrix ( 0,0, Qu._Ch() ); 108 qux.set_submatrix ( 0,dimx+dimu+dimy,Qu._Ch() ); 109 110 qyx.set_submatrix ( 0,0,S->_C() ); 111 qyx.set_submatrix ( 0,dimx,-eye ( dimy ) ); 112 113 // parts of QR 114 hqy=Qy.to_mat()*qyx*pr; 115 116 // pre_qr 117 pre_qr = concat_vertical ( s*pr, concat_vertical ( hqy, qux ) ); 118 } 119 //! set penalization matrices and control horizon 120 void set_control_parameters ( const mat &Qy0, const mat &Qu0, const vec &y_req0, int horizon0 ); 121 //! set penalization matrices and control horizon 122 void set_control_Qy ( const mat &Qy0 ) {Qy=Qy0;} 123 //! refresh temporary storage - inefficient can be improved 124 void initial_belmann() {s=1e-5*eye ( dimx+dimu+dimy );}; 125 //! validation procedure 126 void validate(); 127 //! function for future use which is called at each time td; Should call initialize()! 128 //! redesign one step of the 129 void ricatti_step() { 130 pre_qr.set_submatrix ( 0,0,s*pr ); 131 pre_qr.set_submatrix ( dimx+dimu+dimy, dimu+dimx, -Qy.to_mat()*y_req ); 132 if ( !qr ( pre_qr,post_qr ) ) { bdm_warning ( "QR in LQG unstable" );} 133 triu ( post_qr ); 134 // hn(m+1:2*m+n+r,m+1:2*m+n+r); 135 s=post_qr.get ( dimu, 2*dimu+dimx+dimy-1, dimu, 2*dimu+dimx+dimy-1 ); 136 }; 137 void redesign() { 138 for ( td=horizon; td>0; td-- ) { 139 update_system(); 140 ricatti_step(); 129 141 } 130 //! compute control action 131 vec ctrlaction(const vec &state, const vec &ukm){vec pom=concat(state, ones(dimy), ukm); return L*pom;} 132 } ; 142 /* ws=hn(1:m,m+1:2*m+n+r); 143 wsd=hn(1:m,1:m); 144 Lklq=-inv(wsd)*ws;*/ 145 L = -inv ( post_qr.get ( 0,dimu-1, 0,dimu-1 ) ) * post_qr.get ( 0,dimu-1, dimu, 2*dimu+dimx+dimy-1 ); 146 } 147 //! compute control action 148 vec ctrlaction ( const vec &state, const vec &ukm ) const { 149 vec pom=concat ( state, ones ( dimy ), ukm ); 150 return L*pom; 151 } 152 } ; 133 153 134 154 135 155 } // namespace -
library/bdm/estim/arx.cpp
r700 r723 98 98 } 99 99 100 mlnorm<ldmat>* ARX::predictor ( ) const {101 mat mu ( dimy, posterior()._V().rows() - dimy );102 mat R ( dimy, dimy );103 mlnorm<ldmat>* tmp;104 tmp = new mlnorm<ldmat> ( );105 106 est.mean_mat ( mu, R ); //mu =107 mu = mu.T();108 //correction for student-t -- TODO check if correct!!109 110 if ( have_constant) { // constant term111 //Assume the constant term is the last one:112 tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ) );113 } else {114 tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ) );115 }116 return tmp;117 }118 100 119 101 mlstudent* ARX::predictor_student ( ) const { -
library/bdm/estim/arx.h
r700 r723 98 98 } 99 99 //! conditional version of the predictor 100 mlnorm<ldmat>* predictor() const; 100 template<class sq_T> 101 shared_ptr<mlnorm<sq_T> > ml_predictor() const; 102 //! fast version of predicto 103 template<class sq_T> 104 void ml_predictor_update(mlnorm<sq_T> &pred) const; 101 105 mlstudent* predictor_student() const; 102 106 //! Brute force structure estimation.\return indeces of accepted regressors. … … 188 192 }; 189 193 UIREGISTER(ARXfrg); 194 195 196 //////////////////// 197 template<class sq_T> 198 shared_ptr< mlnorm<sq_T> > ARX::ml_predictor ( ) const { 199 shared_ptr< mlnorm<sq_T> > tmp = new mlnorm<sq_T> ( ); 200 tmp->set_rv(yrv); 201 tmp->set_rvc(_rvc()); 202 203 ml_predictor_update(*tmp); 204 tmp->validate(); 205 return tmp; 206 } 207 208 template<class sq_T> 209 void ARX::ml_predictor_update(mlnorm<sq_T> &pred) const { 210 mat mu ( dimy, posterior()._V().rows() - dimy ); 211 mat R ( dimy, dimy ); 212 213 est.mean_mat ( mu, R ); //mu = 214 mu = mu.T(); 215 //correction for student-t -- TODO check if correct!! 216 217 if ( have_constant) { // constant term 218 //Assume the constant term is the last one: 219 pred.set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), sq_T ( R ) ); 220 } else { 221 pred.set_parameters ( mu, zeros ( dimy ), sq_T( R ) ); 222 } 223 } 224 190 225 }; 191 226 #endif // AR_H -
library/bdm/estim/kalman.h
r703 r723 27 27 * \brief Basic elements of linear state-space model 28 28 29 Parameter evolution model:\f[ x_ t = A x_{t-1} + B u_t + Q^{1/2} e_t \f]30 Observation model: \f[ y_t = C x_{t -1} + C u_t + Q^{1/2} w_t. \f]31 Where $e_t$ and $w_t$ are independent vectorsNormal(0,1)-distributed disturbances.29 Parameter evolution model:\f[ x_{t+1} = A x_{t} + B u_t + Q^{1/2} e_t \f] 30 Observation model: \f[ y_t = C x_{t} + C u_t + R^{1/2} w_t. \f] 31 Where $e_t$ and $w_t$ are mutually independent vectors of Normal(0,1)-distributed disturbances. 32 32 */ 33 33 template<class sq_T> … … 61 61 UI::get (C, set, "C", UI::compulsory); 62 62 UI::get (D, set, "D", UI::compulsory); 63 mat Qtm, Rtm; 63 mat Qtm, Rtm; // full matrices 64 64 if(!UI::get(Qtm, set, "Q", UI::optional)){ 65 65 vec dq; … … 72 72 Rtm=diag(dr); 73 73 } 74 R=Rtm; // automatic conversion 74 R=Rtm; // automatic conversion to square-root form 75 75 Q=Qtm; 76 76 … … 519 519 520 520 */ 521 class StateFromARX: public StateSpace< fsqmat>{521 class StateFromARX: public StateSpace<chmat>{ 522 522 protected: 523 523 //! remember connection from theta ->A 524 524 datalink_part th2A; 525 //! xrv526 RV xrv;525 //! remember connection from theta ->B 526 datalink_part th2B; 527 527 //!function adds n diagonal elements from given starting point r,c 528 528 void diagonal_part(mat &A, int r, int c, int n){ 529 529 for(int i=0; i<n; i++){A(r,c)=1.0; r++;c++;} 530 530 }; 531 //! similar to ARX.have_constant 532 bool have_constant; 531 533 public: 532 534 //! set up this object to match given mlnorm 533 void connect_mlnorm(const mlnorm<fsqmat> &ml){ 535 //! Note that state-space and common mpdf use different meaning of \f$ _t \f$ in \f$ u_t \f$. 536 //!While mlnorm typically assumes that \f$ u_t \rightarrow y_t \f$ in state space it is \f$ u_{t-1} \rightarrow y_t \f$ 537 //! For consequences in notation of internal variable xt see arx2statespace_notes.lyx. 538 void connect_mlnorm(const mlnorm<chmat> &ml, RV &xrv, RV &urv){ 539 534 540 //get ids of yrv 535 541 const RV &yrv = ml._rv(); … … 537 543 const RV &rgr = ml._rvc(); 538 544 RV rgr0 = rgr.remove_time(); 539 RVurv = rgr0.subt(yrv);545 urv = rgr0.subt(yrv); 540 546 541 547 // create names for state variables 542 548 xrv=yrv; 543 549 544 int y_multiplicity = rgr.mint(yrv);545 int y_block_size = yrv.length()*(y_multiplicity +1); // current yt + all delayed yts546 for (int m=0;m<y_multiplicity ;m++){550 int y_multiplicity = -rgr.mint(yrv); 551 int y_block_size = yrv.length()*(y_multiplicity); // current yt + all delayed yts 552 for (int m=0;m<y_multiplicity - 1;m++){ // ========= -1 is important see arx2statespace_notes 547 553 xrv.add(yrv.copy_t(-m-1)); //add delayed yt 548 554 } 555 //! temporary RV for connection to ml.rvc, since notation of xrv and ml.rvc does not match 556 RV xrv_ml = xrv.copy_t(-1); 557 549 558 // add regressors 550 559 ivec u_block_sizes(urv.length()); // no_blocks = yt + unique rgr … … 552 561 RV R=urv.subselect(vec_1(r)); //non-delayed element of rgr 553 562 int r_size=urv.size(r); 554 int r_multiplicity= rgr.mint(R);563 int r_multiplicity=-rgr.mint(R); 555 564 u_block_sizes(r)= r_size*r_multiplicity ; 556 565 for(int m=0;m<r_multiplicity;m++){ 557 566 xrv.add(R.copy_t(-m-1)); //add delayed yt 567 xrv_ml.add(R.copy_t(-m-1)); //add delayed yt 558 568 } 559 569 } 560 570 // add constant 571 if (any(ml._mu_const()!=0.0)){ 572 have_constant=true; 573 xrv.add(RV("bdm_reserved_constant_one",1)); 574 } else {have_constant=false;} 575 576 561 577 // get mapp 562 th2A.set_connection(xrv, ml._rvc()); 578 th2A.set_connection(xrv_ml, ml._rvc()); 579 th2B.set_connection(urv, ml._rvc()); 563 580 564 581 //set matrix sizes … … 577 594 this->C=zeros(yrv._dsize(), xrv._dsize()); 578 595 this->C.set_submatrix(0,0,eye(yrv._dsize())); 579 580 this->R = zeros(yrv._dsize(),yrv._dsize());581 this->Q = zeros(xrv._dsize(),xrv._dsize());582 583 584 585 596 this->D=zeros(yrv._dsize(),urv._dsize()); 597 this->R.setCh(zeros(yrv._dsize(),yrv._dsize())); 598 this->Q.setCh(zeros(xrv._dsize(),xrv._dsize())); 599 // Q is set by update 600 601 update_from(ml); 602 validate(); 586 603 }; 587 604 //! fast function to update parameters from ml - not checked for compatibility!! 588 void update_from(const mlnorm<fsqmat> &ml){ 589 590 vec theta = ml._A().get_row(0); // this 591 this->Q._M().set_submatrix(0,0,ml._R()); 592 593 } 605 void update_from(const mlnorm<chmat> &ml){ 606 607 vec Arow=zeros(A.cols()); 608 vec Brow=zeros(B.cols()); 609 // ROW- WISE EVALUATION ===== 610 for(int i=0;i<ml._rv()._dsize(); i++){ 611 612 vec theta = ml._A().get_row(i); 613 614 th2A.filldown(theta,Arow); 615 if (have_constant){ 616 // constant is always at the end no need for datalink 617 Arow(A.cols()-1) = ml._mu_const()(i); 618 } 619 this->A.set_row(i,Arow); 620 621 th2B.filldown(theta,Brow); 622 this->B.set_row(i,Brow); 623 } 624 this->Q._Ch().set_submatrix(0,0,ml.__R()._Ch()); 625 626 }; 627 //! access function 628 bool _have_constant() const {return have_constant;} 594 629 }; 595 630 -
library/bdm/itpp_ext.cpp
r712 r723 108 108 bvec operator| (const bvec &a, const bvec &b) 109 109 { 110 it_assert_debug (b.size() != a.size(), "operator&(): Vectors of different lengths");110 it_assert_debug (b.size() == a.size(), "operator|(): Vectors of different lengths"); 111 111 112 112 bvec temp (a.size()); -
library/bdm/itpp_ext.h
r712 r723 45 45 } 46 46 47 const double inf = std::numeric_limits<double>::infinity();47 static const double inf = std::numeric_limits<double>::infinity(); 48 48 49 49 //#if 0 -
library/bdm/math/chmat.h
r660 r723 106 106 copy_vector ( dim*dim, chQ._data(), Ch._data() ); 107 107 } 108 108 //! access function 109 void setCh ( const chmat &Ch0 ) { 110 Ch = Ch0._Ch(); 111 dim=Ch0.rows(); 112 } 113 //! access function 114 void setCh ( const mat &Ch0 ) { 115 //TODO check if Ch0 is OK!!! 116 Ch = Ch0; 117 dim=Ch0.rows(); 118 } 119 109 120 //! Access functions 110 121 void setD ( const vec &nD, int i ) { -
library/bdm/stat/exp_family.h
r713 r723 708 708 //!access function 709 709 mat _R() const { return this->iepdf._R().to_mat(); } 710 710 //!access function 711 sq_T __R() const { return this->iepdf._R(); } 712 711 713 //! Debug stream 712 714 template<typename sq_M> -
library/tests/testsuite/LQG_test.cpp
r722 r723 5 5 using namespace bdm; 6 6 7 TEST ( LQG_test) {7 TEST(LQG_test) { 8 8 LQG reg; 9 shared_ptr<StateSpace< fsqmat> > stsp = new StateSpace<fsqmat>;9 shared_ptr<StateSpace<chmat> > stsp=new StateSpace<chmat>; 10 10 // 2 x 1 x 1 11 stsp-> set_parameters ( eye ( 2 ), ones ( 2, 1 ), ones ( 1, 2 ), ones ( 1, 1 ), /* Q,R */ eye ( 2 ), eye ( 1 ) ); 12 reg.set_system ( stsp ); // A, B, C 13 reg.set_control_parameters ( eye ( 1 ), eye ( 1 ), vec_1 ( 1.0 ), 6 ); //Qy, Qu, horizon 11 stsp-> set_parameters(eye(2), ones(2,1), ones(1,2), ones(1,1), /* Q,R */ eye(2), eye(1)); 12 reg.set_system(stsp); // A, B, C 13 reg.set_control_parameters(eye(1), eye(1), vec_1(1.0), 6); //Qy, Qu, horizon 14 reg.validate(); 14 15 15 16 reg.redesign(); 16 double reg_apply = reg.ctrlaction ( "0.5, 1.1", "0.0" ) ( 0); /*convert vec to double*/17 CHECK_CLOSE ( reg_apply, -0.248528137234392, 0.0001);17 double reg_apply=reg.ctrlaction("0.5, 1.1","0.0")(0); /*convert vec to double*/ 18 CHECK_CLOSE(reg_apply, -0.248528137234392, 0.0001); 18 19 } 19 20 20 TEST ( to_state_test) {21 TEST(to_state_test) { 21 22 mlnorm<fsqmat> ml; 22 mat A ="1.1, 2.3";23 ml.set_parameters ( A, vec_1 ( 1.3 ), eye ( 1 ));24 RV yr = RV ( "y", 1);25 RV ur = RV ( "u", 1);26 ml.set_rv ( yr);27 yr.t_plus ( -1);28 ml.set_rvc ( concat ( yr, ur ));29 30 shared_ptr<StateCanonical > Stsp =new StateCanonical;31 Stsp->connect_mlnorm ( ml);32 33 /* results from 23 mat A="1.1, 2.3"; 24 ml.set_parameters(A, vec_1(1.3), eye(1)); 25 RV yr=RV("y",1); 26 RV ur=RV("u",1); 27 ml.set_rv(yr); 28 yr.t_plus(-1); 29 ml.set_rvc(concat(yr, ur)); 30 31 shared_ptr<StateCanonical > Stsp=new StateCanonical; 32 Stsp->connect_mlnorm(ml); 33 34 /* results from 34 35 [A,B,C,D]=tf2ss([2.3 0],[1 -1.1]) 35 36 */ 36 CHECK_CLOSE_EX ( Stsp->_A().get_row ( 0 ), vec ( "1.1" ), 0.0001);37 CHECK_CLOSE_EX ( Stsp->_C().get_row ( 0 ), vec ( "2.53" ), 0.0001);38 CHECK_CLOSE_EX ( Stsp->_D().get_row ( 0 ), vec ( "2.30" ), 0.0001);37 CHECK_CLOSE_EX(Stsp->_A().get_row(0), vec("1.1"), 0.0001); 38 CHECK_CLOSE_EX(Stsp->_C().get_row(0), vec("2.53"), 0.0001); 39 CHECK_CLOSE_EX(Stsp->_D().get_row(0), vec("2.30"), 0.0001); 39 40 } 40 41 41 TEST ( to_state_arx_test) {42 mlnorm< fsqmat> ml;43 mat A = "1.1, 2.3";44 ml.set_parameters ( A, vec_1 ( 1.3 ), eye ( 1 ));45 RV yr = RV ( "y", 1);46 RV ur = RV ( "u", 1);47 ml.set_rv ( yr);48 yr.t_plus ( -1);49 ml.set_rvc ( concat ( yr, ur ) );50 51 shared_ptr<StateFromARX> Stsp = new StateFromARX;52 Stsp->connect_mlnorm ( ml);53 54 /* results from 42 TEST(to_state_arx_test) { 43 mlnorm<chmat> ml; 44 mat A="1.1, 2.3, 3.4"; 45 ml.set_parameters(A, vec_1(1.3), eye(1)); 46 RV yr=RV("y",1); 47 RV ur=RV("u",1); 48 ml.set_rv(yr); 49 ml.set_rvc(concat(yr.copy_t(-1), concat(ur, ur.copy_t(-1)))); 50 51 shared_ptr<StateFromARX> Stsp=new StateFromARX; 52 RV xrv; RV urv; 53 Stsp->connect_mlnorm(ml,xrv,urv); 54 55 /* results from 55 56 [A,B,C,D]=tf2ss([2.3 0],[1 -1.1]) 56 57 */ 57 cout << "---" << endl; 58 CHECK_CLOSE_EX(Stsp->_A().get_row(0), vec("1.1, 3.4, 1.3"), 0.0001); 59 CHECK_CLOSE_EX(Stsp->_B().get_row(0), vec("2.3"), 0.0001); 60 CHECK_CLOSE_EX(Stsp->_C().get_row(0), vec("1, 0, 0"), 0.0001); 58 61 } 59 62 60 TEST ( arx_LQG_test ){61 mlnorm< fsqmat> ml;62 mat A = "1.81, -.8189, .00468, .00438";63 ml.set_parameters ( A, vec_1 ( 0.0 ), 0.00001*eye ( 1 ));64 RV yr = RV ( "y", 1);65 RV ur = RV ( "u", 1);66 RV rgr = yr.copy_t ( -1);67 rgr.add ( yr.copy_t ( -2 ));68 rgr.add ( yr.copy_t ( -2 ));69 rgr.add ( ur.copy_t ( -2 ));70 rgr.add ( ur.copy_t ( -1 ));71 72 ml.set_rv ( yr);73 ml.set_rvc ( rgr);63 TEST(arx_LQG_test){ 64 mlnorm<chmat> ml; 65 mat A="1.81, -.81, .00468, .00438"; 66 ml.set_parameters(A, vec_1(0.0), 0.00001*eye(1)); 67 RV yr=RV("y",1); 68 RV ur=RV("u",1); 69 RV rgr = yr.copy_t(-1); 70 rgr.add(yr.copy_t(-2)); 71 rgr.add(yr.copy_t(-2)); 72 rgr.add(ur.copy_t(-1)); 73 rgr.add(ur); 74 75 ml.set_rv(yr); 76 ml.set_rvc(rgr); 74 77 ml.validate(); 75 76 shared_ptr<StateCanonical > Stsp = new StateCanonical; 77 Stsp->connect_mlnorm ( ml ); 78 78 79 shared_ptr<StateFromARX> Stsp=new StateFromARX; 80 RV xrv; RV urv; 81 Stsp->connect_mlnorm(ml,xrv,urv); 82 79 83 LQG L; 80 L.set_system ( Stsp);81 L.set_control_parameters ( eye ( 1 ), eye ( 1 ), vec_1 ( 0.0 ), 100);84 L.set_system(Stsp); 85 L.set_control_parameters(eye(1), sqrt(1.0/1000)*eye(1), vec_1(0.0), 100); 82 86 L.validate(); 83 87 84 88 L.redesign(); 85 cout << L.to_string() <<endl;89 cout << L.to_string()<<endl; 86 90 } -
library/tests/testsuite/datalink_test.cpp
r722 r723 101 101 102 102 vec val_up ( "1 2" ); 103 dl.st ep ( val_up);103 dl.store_data(val_up); 104 104 val_up = "3 4"; 105 105 -
library/tests/testsuite/logger_test.cpp
r722 r723 64 64 remove_all ( ls.c_str() ); 65 65 makedir ( ls, false ); 66 remove_all ( " exp/dirfile" );66 remove_all ( "dirfilelog_files" ); 67 67 68 dirfilelog L ( " exp/dirfile", 10 );68 dirfilelog L ( "dirfilelog_files", 10 ); 69 69 70 70 int rid = L.add ( r, "" ); … … 82 82 83 83 std::string expected ( load_test_file ( "logger_test_dirfile_format.matrix" ) ); 84 std::string actual ( load_test_file ( " exp/dirfile/format" ) );84 std::string actual ( load_test_file ( "dirfilelog_files/format" ) ); 85 85 CHECK_EQUAL ( expected, actual ); 86 86 }
