Changeset 605

Timestamp:

09/06/09 22:53:08 (15 years ago)

Author:

smidl

Message:

StateCanonical? = StateSpace? model with connection to ARX models

Location:

library

Files:

• bdm/estim/kalman.h (modified) (2 diffs)
• tests/LQG_test.cpp (modified) (2 diffs)

library/bdm/estim/kalman.h

@@ -85 +85 @@
                 int _dimu(){return dimu;}
                 //! access function
-                mat _A(){return A;}
-                //! access function
-                mat _B(){return B;}
-                //! access function
-                mat _C(){return C;}
-                //! access function
-                mat _D(){return D;}
-                //! access function
-                sq_T _Q(){return Q;}
-                //! access function
-                sq_T _R(){return R;}
+                const mat& _A() const {return A;}
+                //! access function
+                const mat& _B()const {return B;}
+                //! access function
+                const mat& _C()const {return C;}
+                //! access function
+                const mat& _D()const {return D;}
+                //! access function
+                const sq_T& _Q()const {return Q;}
+                //! access function
+                const sq_T& _R()const {return R;}
 };
 
@@ -356 +356 @@
 \f[ x_{t+1} = Ax_t + B u_t + R^{1/2} e_t, y_t=Cx_t+Du_t + R^{1/2}w_t, \f]
 For example, for:
-\f[ y_t = a y_{t-1} + b u_{t-1}\f]
-the state \f$ x_t = [y_{t-1}, u_{t-1}] \f$
-*/
-template<class sq_T>
-shared_ptr<StateSpace<fsqmat> > to_state_space(const mlnorm<sq_T > &ml, ivec &theta_in_A, ivec &theta_in_C){
-        //get ids of yrv
-                        
-                        ivec yids= unique(ml._rv()._ids()); //yrv is now stored in _yrv
-                        ivec dids= unique(ml._rvc()._ids()); //yrv is now stored in _yrv
-                        ivec uids=unique_complement(dids, yids);
-        
-                        RV yrv_uni = RV(yids,zeros_i(yids.length()));
-                        RV urv_uni = RV(uids,zeros_i(yids.length()));
+Using Frobenius form, see [].
+
+For easier use in the future, indeces theta_in_A and theta_in_C are set. TODO - explain
+*/
+//template<class sq_T>
+class StateCanonical: public StateSpace<fsqmat>{
+        protected:
+                //! remember connection from theta ->A
+                datalink_part th2A;
+                //! remember connection from theta ->C
+                datalink_part th2C;
+                //! remember connection from theta ->D
+                datalink_part th2D;
+                //!cached first row of A
+                vec A1row;
+                //!cached first row of C
+                vec C1row;
+                //!cached first row of D
+                vec D1row;
+                
+        public:
+                //! set up this object to match given mlnorm
+                void connect_mlnorm(const mlnorm<fsqmat > &ml){
+        //get ids of yrv                                
+                        const RV &yrv = ml._rv();
+                        //need to determine u_t - it is all in _rvc that is not in ml._rv()
+                        RV rgr0 = ml._rvc().remove_time();
+                        RV urv = rgr0.subt(yrv); 
+                        
                         //We can do only 1d now... :(
-                        bdm_assert_debug(yrv_uni._dsize()==urv_uni._dsize()==1, "Only for SISO so far..." );
-                                        
+                        bdm_assert_debug(yrv._dsize()==1, "Only for SISO so far..." );
+
+                        // create names for 
                         RV xrv; //empty
                         RV Crv; //empty
-                        int td=ml._rvc().mintd();
+                        int td=ml._rvc().mint();
+                        // assuming strictly proper function!!!
                         for (int t=-1;t>=td;t--){
-                                xrv.add(yrv_uni);
-                                Crv.add(urv_uni);
-                        }
-                        
-                        int dimx = xrv._dsize();
-                        
-                        theta_in_A = ml._rv().dataind(xrv);
-                        theta_in_C = ml._rvc().dataind(xrv);
-                        // some chcek of corretness
-                        
-                        vec A1row = zeros(xrv._dsize());
-                        vec C1row = zeros(xrv._dsize());
+                                xrv.add(yrv.copy_t(t));
+                                Crv.add(urv.copy_t(t));
+                        }
+                        
+                        this->dimx = xrv._dsize();
+                        this->dimy = yrv._dsize();
+                        this->dimu = urv._dsize();
+                        
+                        // get mapp
+                        th2A.set_connection(xrv, ml._rvc());
+                        th2C.set_connection(Crv, ml._rvc());
+                        th2D.set_connection(urv, ml._rvc());
+
+                        //set matrix sizes
+                        this->A=zeros(dimx,dimx);
+                        for (int j=1; j<dimx; j++){A(j,j-1)=1.0;} // off diagonal
+                                this->B=zeros(dimx,1);
+                                this->B(0) = 1.0;
+                                this->C=zeros(1,dimx);
+                                this->D=zeros(1,urv._dsize());
+                                this->Q = zeros(dimx,dimx);
+                        // R is set by update
+                        
+                        //set cache
+                        this->A1row = zeros(xrv._dsize());
+                        this->C1row = zeros(xrv._dsize());
+                        this->D1row = zeros(urv._dsize());
+                        
+                        update_from(ml);
+                        validate();
+                };
+                //! fast function to update parameters from ml - not checked for compatibility!!
+                void update_from(const mlnorm<fsqmat> &ml){
+                        
                         vec theta = ml._A().get_row(0); // this 
-                        set_subvector( A1row, theta_in_A, theta);
-                        set_subvector( C1row, theta_in_C, theta);
-
-                        StateSpace<fsqmat> stsp=new StateSpace<fsqmat>();
-                        mat A=zeros(dimx,dimx);
+                        
+                        th2A.filldown(theta,A1row);
+                        th2C.filldown(theta,C1row);
+                        th2D.filldown(theta,D1row);
+
+                        R = ml._R();
+
                         A.set_row(0,A1row);
-                        for (int j=1; j<dimx; j++){A(j,j-1)=1.0;} // off diagonal
-                        mat B=zeros(dimx,1);
-                        B(0) = 1.0;
-                        mat C=zeros(1,dimx);
-                        C.set_row(0,C1row);
-                        stsp.set_parameters(A,B,C,zeros(1,1), zeros(dimx,dimx), ml._R());
-                        return stsp;
-}
+                        C.set_row(0,C1row+D1row*A1row);
+                        D.set_row(0,D1row);
+                        
+                }
+};
 
 /////////// INSTANTIATION

library/tests/LQG_test.cpp

@@ -1 +1 @@
 #define BDMLIB
-#include "UnitTest++.h"
+#include "mat_checks.h"
 #include "design/ctrlbase.h"
 
@@ -18 +18 @@
 }
 
+TEST(test_to_state) {
+        mlnorm<fsqmat> ml;
+        mat A="1.1, 2.3";
+        ml.set_parameters(A, vec_1(1.3), eye(1));
+        RV yr=RV("y",1);
+        RV ur=RV("u",1);
+        ml.set_rv(yr);
+        yr.t_plus(-1);
+        ml.set_rvc(concat(yr, ur));
+        
+        shared_ptr<StateCanonical > Stsp=new StateCanonical;
+        Stsp->connect_mlnorm(ml);
+        
+        /* results from 
+        [A,B,C,D]=tf2ss([2.3 0],[1 -1.1])
+        */
+        CHECK_CLOSE_EX(Stsp->_A().get_row(0), vec("1.1"), 0.0001);
+        CHECK_CLOSE_EX(Stsp->_C().get_row(0), vec("2.53"), 0.0001);
+        CHECK_CLOSE_EX(Stsp->_D().get_row(0), vec("2.30"), 0.0001);
+}
+
+