Botter自主ゼミノート 7.2 線形システムの最適制御 (1)

やること

https://www.amazon.co.jp/dp/4254209444
を読んで、確率微分方程式による最適化問題を解けるようになることです。

これまでのBotter自主ゼミノート

Botter自主ゼミノート 1.2 確定システムの制御の回顧

Botter自主ゼミノート 2.1 確率過程とは？, 2.2 確率過程の数学的表現

Botter自主ゼミノート 1.2 数式導出

Botter自主ゼミノート 2.3 確率モーメント

Botter自主ゼミノート 2.4 確率過程の分類

Botter自主ゼミノート 2.5 エルゴード性

Botter自主ゼミノート 2.6 確率過程の周波数表現

Botter自主ゼミノート 2.7 マルコフ過程

Botter自主ゼミノート 2.8 正規型確率過程

Botter自主ゼミノート 2.9 ウィーナ過程(1)

Botter自主ゼミノート 2.9 ウィーナ過程(2)

Botter自主ゼミノート 2.10 白色雑音

Botter自主ゼミノート 3.1, 3.2 確率変数列の収束

Botter自主ゼミノート 3.3 確率過程の連続性

Botter自主ゼミノート 3.4 自乗平均微分

Botter自主ゼミノート 3.5 自乗平均積分

Botter自主ゼミノート 4.1 確率微分方程式とは？

Botter自主ゼミノート 4.2 確率積分

Botter自主ゼミノート 4.2 確率積分 例題4.1

Botter自主ゼミノート 4.3 確率微分方程式

Botter自主ゼミノート 4.4 伊藤の確率微分演算

Botter自主ゼミノート 4.5 拡散過程

Botter自主ゼミノート 4.6 確率密度関数の時間進化 - コルモゴロフ方程式

Botter自主ゼミノート 6.1 動的システムの推定とは？

Botter自主ゼミノート 6.2 条件付き確率密度関数の時間進化

Botter自主ゼミノート 6.3 モーメント関数の時間進化

Botter自主ゼミノート 6.4 カルマンフィルタ

Botter自主ゼミノート 6.5 イノベーション過程

7.2 線形システムの最適制御 (1)

制御対象の設定

制御するシステムのモデルを

$$dx(t) = A(t)x(t)dt + C(t)u(t)dt + G(t)dw(t), \quad x(t_0) = x_0 \tag{7.1}$$

とします。$x \in R^n$はシステムの状態量ベクトル、$u \in R^l$は制御量ベクトルです。さらに、$w(t)$はシステム雑音を表すウィーナ過程で、平均$0$、共分散マトリクス$Q(t)dt$を持つものとします。

$x_0$は平均値$m_0$、共分散マトリクス$P_0$の正規性確率変数とします。
制御では、現在時刻$t$までのシステムの状態量$\{x(\tau), t_0 \le \tau \le t\}$はすべて入手可能であるとして、評価汎関数を以下のような二次形式とします。$F, M(\cdot) \ge 0, N(\cdot) > 0$は対称マトリクスとなります。

$$J(u) = \mathcal{E}\left\{ x^T(T)Fx(T) + \int_{t_0}^T [x^T(\tau)M(\tau)x(\tau) + u^T(\tau)N(\tau)u(\tau)]d\tau \right\} \tag{7.2}$$

確定システムの制御問題では右辺の期待値を取る必要はありませんが、確率システムでは状態量$x(t)$と制御量$u(t)$が確率変数であるので、期待値を取る必要があります。式(7.2)の$J(u)$を最小とする制御量$\{u^o(t), t_0 \le t < T\}$を最適制御量 (optimal control vector) と呼びます。

二次形式にする意味は単純に便利ということだと思います。高校数学の美しい物語を参考にしました。

ここで、教科書では最適制御量が存在する条件を書いていますが、正直ここの部分は理解でき無いので省略します。

評価汎関数(7.2)に対して、次のスカラ汎関数(7.2)を定義します。この式は、時刻$t$に置いて$x(t) = x$という確定値を取ったという条件下での、残り時間区間$[t, T]$残すとの最小値を表しています。これを最小コスト汎関数と呼びます。

$$V(t,x) = \min_{u(s), t \le s \le T} \mathcal{E}\left\{ x^T(T)Fx(T) + \int_{t}^T [x^T(s)M(s)x(s) + u^T(s)N(s)u(s)]ds \bigg| x(t) = x \right\}$$

ここで、時間区間$[t,T]$を微小区間$[t, t+\Delta t)$と残り区間$[t+\Delta t, T]$に分割します。

この時、条件付き期待値の性質$\mathcal{E}\{\cdot | x(t) = x\} = \mathcal{E}\{\mathcal{E} \{\cdot | x(t+\Delta t) = x(t + \Delta x)\} | x(t) = x\}$に留意して分割を行うと以下のようになります。

$$V(t,x) = \min_{u(s), t \le s \le T} \mathcal{E}\left\{ \int_t^{t+\Delta t} [x^T(s)M(s)x(s) + u^T(s)N(s)u(s)] ds + \mathcal{E}\left\{ x^TFx(T) + \int_{t+\Delta t}^T [x^T(s)M(s)x(s) + u^T(s)N(s)u(s)] \bigg| x(t+\Delta t) = x + \Delta x \right\} \bigg| x(t) = x \right\}$$

ここで、最適性の原理 (教科書7.1で説明されていますが省略しています) より、微小区間$[t, t+\Delta t)$でどのような制御が行われていたとしても、それ以降の制御量$\{u(s), t+\Delta t \le s \le T\}$は、その区間でも最適でなければなりません。

そのことを考えて更に式変形を続け、以下の式(7.5)を得ます。

$$\begin{aligned} V(t,x) &= \min_{u(s), t \le s \le T} \mathcal{E}\left\{ \int_t^{t+\Delta t} [x^T(s)M(s)x(s) + u^T(s)N(s)u(s)] ds + \min_{t+\Delta t \le s \le T} \mathcal{E}\left\{ x^TFx(T) + \int_{t+\Delta t}^T [x^T(s)M(s)x(s) + u^T(s)N(s)u(s)] \bigg| x(t+\Delta t) = x + \Delta x \right\} \bigg| x(t) = x \right\} \\ &= \min_{u(s), t \le s \le T} \mathcal{E}\left\{ \int_t^{t+\Delta t} [x^T(s)M(s)x(s) + u^T(s)N(s)u(s)] ds + V(t+\Delta t, x+\Delta x) \bigg| x(t) = x \right\} \end{aligned}$$

さらに、式(7.5)の最右辺において、$\int_t^{t+\Delta t} (x_s^T M x_s + u_s^T N u_s)ds = [x^T(t)M(t)x(t)+u^T(t)N(t)u(t)]\Delta t + o(\Delta t)$であり、さらに$V(t+\Delta t, x+\Delta x)$を$(t,x)$の周りでテイラー展開すると、式(7.5)は以下の式(7.6)のようになります

$$V(t,x) = \min_{u(t)} \mathcal{E}\left\{ [x^T(t) M(t) x(t) + u^T(t)N(t)u(t)]\Delta t + V(x,t) + \left[ \frac{\partial V(t,x)}{\partial t} + \mathcal{L}_x^u V(t,x) \right]\Delta t + \left(\frac{\partial V(t,x)}{\partial x}\right)^T G(t)dw(t) + \mathcal{o}(\Delta t) \bigg| x(t) = x\right\}$$

教科書では、テイラー展開部分は式(4.28)を、式(4.29)で定義される微分生成要素で表記したものを利用しています。

Botter自主ゼミノート 4.4 伊藤の確率微分演算で軽く触れています。

さらに変形を続けます。ここで、式(7.6)に期待値演算を実行します。

$$\begin{aligned} V(t,x) = \min_{u(t)} \left\{ [x^T M(t) x + u^T(t)N(t)u(t)]\Delta t + V(x,t) + \left[ \frac{\partial V(t,x)}{\partial t} + \mathcal{L}_x^u V(t,x) \right]\Delta t + \mathcal{o}(\Delta t) \right\} \end{aligned}$$

右辺$V(t,x)$は確定値なので、$\min$の外に出すことができて、左辺と相殺できます。

$$\begin{aligned} 0 = \min_{u(t)} \left\{ [x^T M(t) x + u^T(t)N(t)u(t)]\Delta t + \left[ \frac{\partial V(t,x)}{\partial t} + \mathcal{L}_x^u V(t,x) \right]\Delta t + \mathcal{o}(\Delta t) \right\} \end{aligned}$$

全体を$\Delta t$で割り算し、$\Delta t \to 0$とします。このとき$o(x)$は$x \to 0$としたとき、$x$よりも$0$に向かう速度が早いので、$\lim_{\Delta t \to 0} \frac{o(\Delta t)}{\Delta t} = 0$となります。

$$\begin{aligned} 0 = \min_{u(t)} \left\{ x^T M(t) x + u^T(t)N(t)u(t) + \frac{\partial V(t,x)}{\partial t} + \mathcal{L}_x^u V(t,x) \right\} \end{aligned}$$

$\frac{\partial V(t,x)}{\partial t}$は確定値なので、$\min$の外に出して左辺に移項できます。

$$\begin{aligned} -\frac{\partial V(t,x)}{\partial t} = \min_{u(t)} \left\{ x^T M(t) x + u^T(t)N(t)u(t) + \mathcal{L}_x^u V(t,x) \right\} \end{aligned}$$

最後に$\mathcal{L}_x^u$を元の形に戻し、式(7.8)を得ます。これを、ベルマン・ハミルトン・ヤコビ方程式と呼びます。

ようやくこの教科書を読み始めたきっかけとなったキーワードが出てきました。

$$-\frac{\partial V(t,x)}{\partial t} = \min_{u(t)} \left\{ x^T M(t) x + u^T(t)N(t)u(t) + \left( \frac{\partial V(t,x)}{\partial x} \right)^T[A(t)x + C(t)u(t)] + \frac{1}{2} \text{tr} \left\{ G(t)Q(t)G^T(t)\frac{\partial}{\partial x} \left( \frac{\partial V(x,t)}{\partial x} \right)^T \right\} \right\}$$

ここで、右辺を$u(t)$に付いて最小化します。教科書では$u(t)$の絡む項を偏微分して、その値が一致するところを最小値としているようです。多分こんな感じかな…？

$$\begin{aligned} u^TNu &= -u^T C^T \frac{\partial V}{\partial x} \\ \frac{\partial(u^TNu)}{\partial u} &= -\frac{\partial u^T C^T \frac{\partial V}{\partial x}}{\partial u} \\ 2Nu &= -C^T \frac{\partial V}{\partial x} \\ u &= -\frac{1}{2} N^{-1} C^T \frac{\partial V}{\partial x} \end{aligned}$$

ここで、最適な$u$をきちんと$u^o(t)$という関数の形に書き直すと…

$$u^o(t) = -\frac{1}{2}N^{-1}(t)C^T(t)\frac{\partial V(t,x)}{\partial x} \tag{7.9}$$

式(7.9)を式(7.8)に代入して、以下の$V(t,x)$に関する非線形微分方程式(7.10a)を得ます。

$$-\frac{\partial V(t,x)}{\partial t} = x^T M(t) x + x^T A^T(t) \frac{\partial V(t,x)}{\partial x} + \frac{1}{2}\text{tr}\left\{ G(t)Q(t)G^T(t)\frac{\partial}{\partial x} \left(\frac{\partial V(t,x)}{\partial x}\right) \right\} - \frac{1}{4}\left(\frac{\partial V(t,x)}{\partial x}\right)^T C(t)N^{-1}(t)C^T(t)\frac{\partial V(t,x)}{\partial x}$$

ここはさらっと教科書に書かれていることを受け入れていますが、後でちゃんと計算したほうがよさそう。

この式を以下の境界条件(7.10b)のもとで解き、$\frac{\partial V(t,x)}{\partial x}$を求めることによって、式(7.9)から最適制御量$u^o(t)$を求めることができます。

$$V(T,x) = x^T Fx \quad (x \in R^n) \tag{7.10b}$$

式(7.10a)は(7.10b)式で与えられる終端拘束条件のみを境界条件として持ち、空間変数$x$に関する境界条件を持ちません。このような方程式を解く問題は、偏微分方程式論ではコーシー問題、あるいは初期値問題と呼ばれます。

意外と何を言っているのかイメージできている気がします。

微分方程式(7.10a)を解くと、最小コストスカラ汎関数$V(t,x)$がある時間$t$でどのような勾配を持つかを求めることができる方程式が得られるのですが、これだけでは$V(t,x)$が実際どのような値になるかはわかりません。

でも、$V(t,x)$がどこかの点を絶対に通るよ、ということがわかれば、そこがアンカーとなってすべての$V(t,x)$の値が確定します。そのアンカーとなるのが、境界条件$V(T,x) = x^T F x$だと理解しました。

マトリクス$F$は二次形式で表現された評価汎関数$J(u)$で定義されていて、最終時刻$T$でシステム状態量$x$が取った値を$J(u)$の一部として取り込むためのものです。

続けて教科書ではこの微分方程式の解き方の解説に入っていきます。面白くなってきました。