( ) http://wwwhayanuemnagoya-uacjp/ fujimoto/ 2011 3 9 11 ( ) 2011/03/09-11 1 / 46
Outline 1 2 3 4 5 ( ) 2011/03/09-11 2 / 46
Outline 1 2 3 4 5 ( ) 2011/03/09-11 3 / 46
(1/2) r + Controller - u Plant y Plant Controller Plant: dx = f(x(t), u(t)) dt y(t) = h(x(t), u(t)) lim y(t) r(t) = 0 t x(t + 1) = A x(t) + B u(t) y(t) = C x(t) + D u(t) ( ) 2011/03/09-11 4 / 46
(2/2) r + - Controller u Plant y x(t + 1) = A x(t) + B u(t) y(t) = C x(t) + D u(t) ( x(t) ) u(t) = K x(t) ( ) 2011/03/09-11 5 / 46
vs ( ) vs ( ) y y o x o x ( ) 2011/03/09-11 6 / 46
vs ( ) ( ) ( ) ( ) ( ) 2011/03/09-11 7 / 46
vs ( ) ( ) ( ) EM ( + + ) ( + ) ( + ) ( + + ) ( ) 2011/03/09-11 8 / 46
?! ( ) 2011/03/09-11 9 / 46
Outline 1 2 3 4 5 ( ) 2011/03/09-11 10 / 46
X p(y X) Y X Y Y X p(x Y) X, Y p(x, Y) = p(x Y)p(Y) = p(y X)p(X) p(x Y) { }} { p(x Y) = ( ) { }} { p(y X) p(y) }{{} = p(y X)p(X)dX {}}{ p(x) ( ) 2011/03/09-11 11 / 46
{ }} { p(x Y) = ( ) { }} { p(y X) {}}{ p(x) p(y X) p(x) dx : p(x Y)! p(x) ( ) p(y X) p(x) ( ) p(y X) ( ) 2011/03/09-11 12 / 46
y = θu + w w N(0, σ w ) ( ) p(w) = p(y θu) = p(u, y θ) 1 (y θu)2 exp 2πσw 2σ w p(θ)! p(θ) = 1 exp (θ µ θ) 2 2πσθ 2σ θ p(u, y θ)p(θ) = ( ) exp(θ 2 )! ( ) 2011/03/09-11 13 / 46
{ xt+1 = Ax t + Bu t + w t, w t N(0, Q), x 0 N(µ 0, Σ 0 ) y t = Cx t + Du t + v t, v t N(0, R) Θ := {A, B, C, D, Q, R} X := {x 0, x 1,, x t, µ 0, Σ 0 } ( ) : Θ X! p(θ, X) p(θ)p(x) p(θ) p(x) ( ) 2011/03/09-11 14 / 46
: p(θ, X Y) q(θ, X Y) p q KL KL(q, p) := q log q p dxdθ q F(q) q log p(y) = F(q) + KL(q, p) } {{ } ( ) q(θ, X Y) = q(θ Y)q(X Y) F(q(Θ, X Y)) Θ, X ( ) ( ) 2011/03/09-11 15 / 46
VB-E VB-M q(x) (new) p(x) exp E q(θ) (old) [log p(x, Y Θ)] q(θ) (new) p(θ) exp E q(x) (old) [log p(θ, Y X)] E( ) p( ), q( ) ( ) ( ) 2011/03/09-11 16 / 46
( ) ( ) 2011/03/09-11 17 / 46
Outline 1 2 3 4 5 ( ) 2011/03/09-11 18 / 46
( ) { xt+1 = Ax t + Bu t + w t, w t N(0, Q), x 0 N(µ 0, Σ 0 ) y t = Cx t + Du t + v t, v t N(0, R) Θ I := {A, B, Q} Θ O := {C, D, R} X := {x 0, x 1,, x t, µ 0, Σ 0 } p(θ I, Θ O, X) p(θ I )p(θ O )p(x) ( ) 2011/03/09-11 19 / 46
: p(y X) p(x) p(x Y) = p(y X)p(X)dX [Beal 03, Barber & Chiappa 07] x t 1 Q, R ( ) p(vec(a, B) µ AB, Q, G) = N(vec(A, B) µ AB, G Q) p(vec(c, D) µ CD, R, H) = N(vec(C, D) µ CD, H R) { p(q ν, S Q ) Q 1 (ν n 1)/2 exp 1 } 2 tr(q 1 S 1 Q ) { p(r η, S R ) R 1 (η l 1)/2 exp 1 } 2 tr(r 1 S 1 R ) ( ) 2011/03/09-11 20 / 46
{ ẋ = Ax + Bu y = Cx + Du x = Tx { x = T AT 1 x + TBu y = CT 1 x + D? ( ) 2011/03/09-11 21 / 46
p(a, B, C, D Q 1, R 1 ) q(a, B, C, D Q 1, R 1 ) f T p(ā, B, C, D Q 1, R 1 )? f T q(ā, B, C, D Q 1, R 1 ) ( ) 2011/03/09-11 22 / 46
p(a, B, C, D Q 1, R 1 ) f T p(ā, B, C, D Q 1, R 1 ) q(a, B, C, D Q 1, R 1 ) f T q(ā, B, C, D Q 1, R 1 ) ( ) 2011/03/09-11 22 / 46
p(θ) x = Tx p(a, B, C, D, Q, R) = p(t AT 1, TB, CT 1, D, TQT T, R)! d(a, B, C, D, Q, R) = d(t AT 1, TB, CT 1, D, TQT T, R)! ( ) 2011/03/09-11 23 / 46
X ( ) Θ? ( ( ) x = Tx?) Yes! ( ) 2011/03/09-11 24 / 46
( A = 1 03 006 094 Ā = ( 1 1 1 1 ) (, B = ) ( 1, B = 1 0 006 ), C = ( 1, 0 ), D = 0 ), C = ( 1, 1 ), D = 1 ( ) 2011/03/09-11 25 / 46
(l 2 ) l 2 : : ± 2σ 5 x 105 4 3 2 1 0 0 20 40 60 80 100 ( ) 2011/03/09-11 26 / 46
(50 ) (50 ) : : : ± 2σ 20 0-20 -40 10-1 10 0 10 1 200 0-200 10-1 10 0 10 1 ( ) 2011/03/09-11 27 / 46
(200 ) (200 ) : : : ± 2σ 40 20 0-20 -40 10-1 10 0 10 1 0-100 -200-300 10-1 10 0 10 1 ( ) 2011/03/09-11 28 / 46
Outline 1 2 3 4 5 ( ) 2011/03/09-11 29 / 46
?! : : ( ) 2011/03/09-11 30 / 46
x t+1 = Ax t + Bu t T 1 ( ): J T (u) = ( ): J (u) = t=0 t=0 J T, J! ( x T t Qx t + u T t Ru ) t + x T T Fx T ( x T t Qx t + u T t Ru ) t u t = u t (x t, t) ( ) 2011/03/09-11 31 / 46
V(x 0 ) ( ) ( V(x t ) := min J(u x t ) = min V(xt+1 ) + x T u t,,u T 1 u t t Qx ) t + u T t Ru t 2 V(x t ) = x T t Π t x t Π t x T t Π ( t x t = min x T Π u t t+1 t+1x t+1 + x T t Qx ) t + u T t Ru t ( = min (Axt + Bu t ) T Π t+1 (Ax t + Bu t ) + x T u t t Qx ) t + u T t Ru t u t 2 ( ) u t = (R + B T Π t+1 B) 1 (B T Π t+1 A)x t Π t =Q + A T Π t+1 A (A T Π t+1 B)(R + B T Π t+1 B) 1 (B T Π t+1 A) : Π N = F : Π t+1 = Π t = Π Π ( ) 2011/03/09-11 32 / 46
LQG LQG(Linear Quadratic Gaussian) ϵ t T 1 ( ): J T (u) = E x t+1 = Ax t + Bu t + Gϵ t t=0 ( ): J (u) = E t=0 ( x T t Qx t + u T t Ru ) t + x T ( x T t Qx t + u T t Ru t) x 0 T Fx T x 0 ( ) 2011/03/09-11 33 / 46
x t+1 = Ax t + Bu t + Gϵ t MCV(Minimum Cost variance) [Sain 66] J(u) =E[Ĵ] + λ var[ĵ] T 1 Ĵ = (x T t Qx t + u T t Ru t) + x T T Fx T t=0 RS(Risk Sensitive) [Whittle 81] J(u) = 2θ 1 log E[exp( θ 2 Ĵ)] =E[Ĵ] θ 4 var[ĵ] + O(θ2 ) A, B ( ) 2011/03/09-11 34 / 46
[De Koning 82, F 10] x t+1 = A t x t + B t u t + G t ϵ t A t, B t, G t, ϵ t T 1 ( ( ): J T (u) = E x T t Qx t + u T t Ru t + tr (S cov[x t+1 x t ]) ) + x T T Fx T t=0 ( ( ): J (u) = E x T t Qx t + u T t Ru t + tr (S cov[x t+1 x t ]) ) t=0 tr (S cov[x t+1 x t ]) = E[x T Sx t+1 t+1 x t ] E[x t+1 ] T SE[x t+1 x t ] ( ) 2011/03/09-11 35 / 46
S = 0 V ( ) V(x t ) := min J(u x t ) = min u t,,u T 1 u t ( E[V(xt+1 ) x t ] + x T t Qx t + u T t Ru ) t 2 V Π t V(x t ) = x T t Π t x t + β t ( ) 2011/03/09-11 36 / 46
u t = (R + Σ BB + E[B T Π t+1 B]) 1 (Σ BA + E[B T Π t+1 A])x t Π t =Q + Σ AA + E[A T Π t+1 A] (Σ AB + E[A T Π t+1 B]) (R + Σ BB + E[B T Π t+1 B]) 1 (Σ BA + E[B T Π t+1 A]) Π T =F Σ XY := E[X T SY] E[X] T SE[Y] ( ) 2011/03/09-11 37 / 46
u t = (R + B T ΠB) 1 B T ΠAx t Π =Q + A T ΠA A T ΠB(R + B T ΠB) 1 B T ΠA (S = 0, T = ) u t = (R + E[B T ΠB]) 1 E[B T ΠA]x t Π =Q + E[A T ΠA] E[A T ΠB](R + E[B T ΠB]) 1 E[B T ΠA] 2 ( ) 2011/03/09-11 38 / 46
(1/2) ( ) 1 01 E[A] = 001 099 ( ) 0 E[B] = 001 0 0 0 0 0 0 0 25 10 5 0 0 0 0 0 0 0 0 0 0 cov[vec( A, B)] = 0 0 0 0245 0 0 0 0 0 0 0 0 0 0 0 0 0 25 10 5 Q = ( ) ( ) 10 0 40 0, R = 1, F = 0 10 0 40 ( ) 2011/03/09-11 39 / 46
数値例 (2/2) 従来法と提案法 (S = 100) の状態 x1 50 25 40 20 30 20 15 state X1 state X1 10 0 10 10 5 20 30 0 40 50 5 0 50 100 150 200 250 0 50 100 time 150 200 250 time 従来法と提案法 (S = 100) の入力 u 3000 150 100 2000 50 1000 0 0 50 1000 100 2000 150 3000 200 4000 250 0 50 藤本 健治 (名古屋大学) 100 150 200 250 0 ばらつき抑制のための確率最適制御 50 100 150 200 250 2011/03/09-11 力学系研究会 40 / 46
( ) x k+1 = A t x t lim E[x t x 0 ] = 0 t lim E[ x t 2 x 0 ] = 0 t [De Koning 82] ( ) 2011/03/09-11 41 / 46
x t+1 = A t x t + B t u t u t = (R + Σ BB + E[B T ΠB]) 1 (Σ BA + E[B T ΠA])x t Π =Q + Σ AA + E[A T ΠA] (Σ AB + E[A T ΠB]) (R + Σ BB + E[B T ΠB]) 1 (Σ BA + E[B T ΠA])! ( ) 2011/03/09-11 42 / 46
: T 1 ( J T (u) =E x T t Qx t + u T t Ru t + tr (S cov[x t+1 x t ]) ) + x T T Fx T t=0 B L B L Π :=Q + L T RL + E[(A BL) T (Π + S)(A BL)] E[(A BL)] T SE[(A BL)] B L u t = Lx t J t (u) := x T 0 B L t Fx 0 B L J J J ( ) 2011/03/09-11 43 / 46
Outline 1 2 3 4 5 ( ) 2011/03/09-11 44 / 46
? or vs? (2 ) ( ) 2011/03/09-11 45 / 46
[1] D Barber and S Chiappa Unified inference for variational Bayesian linear Gaussian state-space models In Advances in Neural Information Processing Systems 19 (NIPS 20), pp 81 88 The MIT Press, 2007 [2] M J Beal Variational Algorithms for Approximate Bayesian inference PhD thesis, University of Londong, London, UK, 2003 [3] W L De Koning Infinite horizon optimal control of linear discrete time systems with stochastic parameters Automatica, Vol 18, No 4, pp 443 453, 1982 [4] K Fujimoto, S Ogawa, Y Ota, and M Nakayama Optimal control of linear systems with stochastic parameters for variance suppression: The finite time horizon case In Proceedings of the 18th IFAC World Congress, 2011 [5] K Fujimoto, A Satoh, and S Fukunaga System identification based on variational bayes method and the inavriance under coordinate transformations Submitted, 2011 [6] M K Sain Control of linear systems according to the minimal variance criterion: A new approach to the disturbance problem IEEE Trans Autom Contr, Vol 11, No 1, pp 118 122, 1966 [7] P Whittle Risk-sensitive Linear/Quadratic/Gaussian control Advances in Applied Probability, Vol 13, pp 764 777, 1981 [8],, Vol 20, No 10, pp 404 412, 2007 [9], H D, Vol J91-D, No 6, pp 1648 1655, 2008 ( ) 2011/03/09-11 46 / 46