1 | /*! |
---|
2 | \file |
---|
3 | \brief Base classes for designers of control strategy |
---|
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 "../base/bdmbase.h" |
---|
14 | |
---|
15 | namespace bdm{ |
---|
16 | |
---|
17 | //! Base class of designers of control strategy |
---|
18 | class Designer : public root { |
---|
19 | public: |
---|
20 | //! Redesign control strategy |
---|
21 | virtual void redesign() { |
---|
22 | bdm_error("Not implemented"); |
---|
23 | } |
---|
24 | |
---|
25 | //! apply control strategy to obtain control input |
---|
26 | virtual vec apply(const vec &cond) { |
---|
27 | bdm_error("Not implemented"); |
---|
28 | return vec(); |
---|
29 | } |
---|
30 | }; |
---|
31 | |
---|
32 | //! Linear Quadratic Gaussian designer for constant penalizations and constant target |
---|
33 | //! Its internals are very close to Kalman estimator |
---|
34 | class LQG : public Designer { |
---|
35 | protected: |
---|
36 | //! dimension of state |
---|
37 | int dimx; |
---|
38 | //! dimension of inputs |
---|
39 | int dimu; |
---|
40 | //! dimension of output |
---|
41 | int dimy; |
---|
42 | |
---|
43 | //! matrix A of the linear system |
---|
44 | mat A; |
---|
45 | //! matrix B of the linear system |
---|
46 | mat B; |
---|
47 | //! matrix C of the linear system |
---|
48 | mat C; |
---|
49 | //! required value of the output y at time t (assumed constant) |
---|
50 | vec y_req; |
---|
51 | |
---|
52 | //! Control horizon, set to maxint for infinite horizons |
---|
53 | int horizon; |
---|
54 | //! penalization matrix Qy |
---|
55 | mat Qy; |
---|
56 | //! penalization matrix Qu |
---|
57 | mat Qu; |
---|
58 | //! time of the design step - from horizon->0 |
---|
59 | int td; |
---|
60 | //! controller parameters |
---|
61 | mat L; |
---|
62 | |
---|
63 | //!@{ \name temporary storage for ricatti |
---|
64 | //! parameters |
---|
65 | mat pr; |
---|
66 | //! penalization |
---|
67 | mat qux; |
---|
68 | //! penalization |
---|
69 | mat qyx; |
---|
70 | //! internal quadratic form |
---|
71 | mat s; |
---|
72 | //! penalization |
---|
73 | mat qy; |
---|
74 | //! pre_qr part |
---|
75 | mat hqy; |
---|
76 | //! pre qr matrix |
---|
77 | mat pre_qr; |
---|
78 | //! post qr matrix |
---|
79 | mat post_qr; |
---|
80 | //!@} |
---|
81 | |
---|
82 | public: |
---|
83 | //! set system parameters from given matrices |
---|
84 | void set_system_parameters(const mat &A, const mat &B, const mat &C); |
---|
85 | //! set penalization matrices and control horizon |
---|
86 | void set_control_parameters(const mat &Qy0, const mat &Qu0, const vec y_req0, int horizon0); |
---|
87 | //! set system parameters from Kalman filter |
---|
88 | // void set_system_parameters(const Kalman &K); |
---|
89 | //! refresh temporary storage - inefficient can be improved |
---|
90 | void prepare_qr(); |
---|
91 | //! function for future use which is called at each time td; Should call prepare_qr()! |
---|
92 | virtual void update_state(){}; |
---|
93 | //! redesign one step of the |
---|
94 | void ricatti_step(){ |
---|
95 | pre_qr.set_submatrix(0,0,s*pr); |
---|
96 | pre_qr.set_submatrix(dimx+dimu+dimy, dimu+dimx, -Qy*y_req); |
---|
97 | if (!qr(pre_qr,post_qr)) bdm_warning("QR in LQG unstable"); |
---|
98 | triu(post_qr); |
---|
99 | // hn(m+1:2*m+n+r,m+1:2*m+n+r); |
---|
100 | s=post_qr.get(dimu, 2*dimu+dimx+dimy-1, dimu, 2*dimu+dimx+dimy-1); |
---|
101 | }; |
---|
102 | void redesign(){ |
---|
103 | for(td=horizon; td>0; td--){ |
---|
104 | update_state(); |
---|
105 | ricatti_step(); |
---|
106 | } |
---|
107 | /* ws=hn(1:m,m+1:2*m+n+r); |
---|
108 | wsd=hn(1:m,1:m); |
---|
109 | Lklq=-inv(wsd)*ws;*/ |
---|
110 | L = -inv(post_qr.get(0,dimu-1, 0,dimu-1)) * post_qr.get(0,dimu-1, dimu, 2*dimu+dimx+dimy-1); |
---|
111 | } |
---|
112 | vec apply(const vec &state, const vec &ukm){vec pom=concat(state, ones(dimy), ukm); return L*pom;} |
---|
113 | } ; |
---|
114 | |
---|
115 | } // namespace |
---|