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

やること

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 イノベーション過程

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

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

前回出てきた式(7.10a)と式(7.10b)を、教科書がどう解いているのかを追いかけていきます。

制御量$u(t)$を$x(t)$に対するフィードバック系と仮定したときの最小コスト汎関数$\~{V}(t,x)$を求める

まず、制御量$u(t)$が

$$u(t) = K(t)x(t) \tag{7.11}$$

で与えられるようなフィードバック形であると仮定します。

直前のベクトル$x$の値に線形になにかをかけ合わせて制御量とする仮定です。非線形の制御量は考えないということですね。

この時、$u \in \mathcal{U}_{ad}$であり、システム方程式(7.1)は

$$dx(t) = \~{A}(t)x(t)dt + G(t)dw(t) \tag{7.12}$$

$$\~{A}(t) = A(t)+C(t)K(t) \tag{7.13}$$

となるので、式(7.4)の最小コスト汎関数$V(t,x)$に対応した$~{V}(t,x)$と$\~{M}(s)$を定義します。

$$\~{V}(t,x) = \mathcal{E}\left\{ x^T(T)Fx(T)+\int_t^T x^T(s)\~{M}(s)x(s)ds \bigg| x(t) = x\right\} \tag{7.14}$$

$$\~{M}(s) = M(s)+K^T(s)N(s)K(s) \tag{7.15}$$

シンプルに$u(t) = K(t)x(t)$を使ってあれこれ置換しているだけです。

ここで、$\~{A}(t)$に対する遷移マトリクスを$\~{\Phi}(t, t_0)$とすると、$s \ge t$に対して

$$x(s) = \~{\Phi}(s, t)x(t)+\int_t^s \~{\Phi}(s, \tau)G(\tau)dw(\tau) \tag{7.16}$$

となります。

遷移マトリクス$\~{\Phi}(t_1, t_0)$は、式(7.12)では差分の形式で$x(t)$を表現しているのを、通常の式で表したときの係数になります。

式(7.16)の第一項は、$x(t)$が$x(s)$にどのような影響を及ぼしているかを表現しています。第二項は$t$から$s$までの各時間で加わるノイズが、それぞれ$x(s)$にどのような影響を及ぼしているのかを表現しています。$\~{\Phi}(s, \tau)$の二番目の引数が$t$ではなくて$\tau$であるところがポイントです。

雑音であろうが信号であろうが、時々刻々と変わる遷移マトリクスによって変形されて、$x(s)$を形作ることになるという考え方ですね。

ここで、式(4.12)と式(4.13)で表される確率積分の性質と、$a^T\~{M}a = \text{tr}\{\~{M}aa^T\}$、さらに式(7.16)を用いて、式(7.14)を変形していきます。

すっかり忘れているので、式(4.12)と式(4.13)を振り返ります。

1. $g(t, \omega), f(t, \omega)$を$\{w(\tau, \omega) - w(t, \omega)\}\:(t \le \tau)$とは独立で、$\int_{t_0}^t \mathcal{E}\{|g(\tau, \omega)|^2\} < \infty, \int_{t_0}^t \mathcal{E}\{|f(\tau, \omega)|^2\} < \infty$を満たす関数とすれば、次の式が成り立ちます。

$$\mathcal{E}\left\{ \int_{t_0}^t g(\tau, \omega)dw(\tau, \omega) \right\} = 0 \tag{4.12}$$

$$\mathcal{E}\left\{ \left[ \int_{t_0}^t g(\tau, \omega)dw(\tau, \omega) \right] \left[ \int_{t_0}^t f(\tau, \omega)dw(\tau, \omega) \right] \right\}\\ = \sigma^2 \int_{t_0}^t \mathcal{E}\{g(\tau, \omega)f(\tau, \omega)\} d\tau \tag{4.13}$$

特に

$$hh \mathcal{E}\left\{ \left| \int_{t_0}^t g(\tau, \omega)dw(\tau, \omega) \right|^2 \right\} = \sigma^2 \int_{t_0}^t \mathcal{E}\{|g(\tau, \omega)|^2\}d\tau \tag{4.14}$$

まずは式(7.14)の積分の中にある$\mathcal{E}\{x^T(s)\~{M}(s)x(s) | x(t) = x\}$を考えていきます。

$$\begin{aligned} \mathcal{E}\{x^T(s)\~{M}(s)x(s) | x(t) = x\} &= x^T\~{\Phi}^T(s,t)\~{M}(s)\~{\Phi}(s,t)x + \mathcal{E}\left\{ \left[ \int_t^s \~{\Phi}(s,\tau)G(\tau)dw(\tau)\right]^T \~{M}(s) \left[ \int_t^s \~{\Phi}(s,\tau)G(\tau)dw(\tau)\right] \bigg| x(t) = x\right\} \\ &= x^T\~{\Phi}^T(s,t)\~{M}(s)\~{\Phi}(s,t)x + \text{tr}\left\{\~{M}(s)\left[\int_t^s \~{\Phi}(s,\tau)G(\tau)Q(\tau)G^T(\tau)\~{\Phi^T}(s,\tau) d\tau \right]\right\} \end{aligned} \tag{7.18}$$

忘れがちですが、$G(t)$は、$d$次元の雑音ベクトルが$n$次元の信号ベクトルにどのような影響を与えるかを定義する行列で、4.3(確率微分方程式)で定義されています。

また、$Q(t)$は$Q(t)dt = \mathcal{E}\{dw \: dw^T\}$で表される雑音の分散共分散行列で、スカラ過程の$\sigma^2$に相当します。これは4.4(伊藤の確率微分演算)で定義されています。

1行目は単純に$x(s)$を式(7.16)で置き換えただけです。1行目から2行目への変形では、式(4.13)を利用した$dw$と$Q(t) dt$の置き換えが行われていて、期待値演算子が見えなくなっています。
h
式(4.12)の性質をこの式変形のどこで使っているのかがちょっとわかりません。

同様にして、式(7.14)の第一項、終端時$x$の価値を表す$x^T(T)Fx(T)$も展開していきます。

$$\mathcal{E}\{x^T(T) F x(T) | x(t) = x\} = x^T \~{\Phi}^T(T, t)F\~{\Phi}(T, t)x + \text{tr}\left\{F\left[\int_t^T \~{\Phi}(T, \tau)G(\tau)Q(\tau)G^T(\tau)\~{\Phi}^T(T, \tau) d\tau\right] d\tau \right\}$$

これらの結果を用いると、式(7.14)は結局以下の式(7.19)のような形になります。

$$\begin{aligned} \~{V}(t, x) &= \mathcal{E}\left\{ x^T(T)Fx(T)+\int_t^T x^T(s)\~{M}(s)x(s)ds \bigg| x(t) = x\right\} \\ &= x^T \~{\Phi}^T(T, t)F\~{\Phi}(T, t)x + \text{tr}\left\{F\left[\int_t^T \~{\Phi}(T, \tau)G(\tau)Q(\tau)G^T(\tau)\~{\Phi}^T(T, \tau) d\tau\right] d\tau \right\} + \int_t^T \left[ x^T\~{\Phi}^T(s,t)\~{M}(s)\~{\Phi}(s,t)x + \text{tr}\left\{\~{M}(s)\left[\int_t^s \~{\Phi}(s,\tau)G(\tau)Q(\tau)G^T(\tau)\~{\Phi^T}(s,\tau) d\tau \right]\right\} \right] ds\\ &= x^T\left[ \~{\Phi}^T(T, t)F\~{\Phi}(T, t) + \int_t^T \~{\Phi}^T(s,t)\~{M}(s)\~{\Phi}(s,t) ds \right]x + \text{tr}\left\{F\left[\int_t^T \~{\Phi}(T, s)G(s)Q(s)G^T(s)\~{\Phi}^T(T, s) ds\right]\right\} + \int_t^T \text{tr}\left\{\~{M}(s)\left[\int_t^s \~{\Phi}(s,\tau)G(\tau)Q(\tau)G^T(\tau)\~{\Phi^T}(s,\tau) d\tau\ \right] \right\} ds \\ \end{aligned}\tag{7.19}$$

ここで、式(7.19)の右辺第二項は$\text{tr}\{AB\} = \text{tr}\{BA\}$を利用して以下のような形にできます。

$$\text{tr}\left\{F\left[\int_t^T \~{\Phi}(T, s)G(s)Q(s)G^T(s)\~{\Phi}^T(T, s) ds\right]\right\} = \int_t^T \text{tr}\left\{G(s)Q(s)G^T(s)\~{\Phi}^T(T, s)F\~{\Phi}^T(T, s) \right\} ds$$

さらに、以下の二重積分の公式$\int_t^T \int_t^s f(s, \tau) d\tau ds = \int_t^T \int_s^T f(\tau,s) d\tau ds$の$f(t,\tau)$を$\text{tr}\{G(\tau)Q(\tau)G^T(\tau)\~{\Phi}^T(s,\tau)\~{M}(\tau)\~{\Phi}(s,\tau)\}$とみなすと、第3項は以下のように変形できます。

$$\int_t^T \text{tr}\left\{ G(s)Q(s)G^T(s) \left[ \int_s^T \~{\Phi}(\tau,s)\~{M}(\tau)\~{\Phi}(\tau,s) d\tau \right] \right\} ds$$

結局、式(7.19)は以下のような形になることがわかります。

$$\begin{aligned} \~{V}(t,x) &= x^T\left[ \~{\Phi}^T(T, t)F\~{\Phi}(T, t) + \int_t^T \~{\Phi}^T(s,t)\~{M}(s)\~{\Phi}(s,t) ds \right]x + \text{tr}\left\{F\left[\int_t^T \~{\Phi}(T, s)G(s)Q(s)G^T(s)\~{\Phi}^T(T, s) ds\right]\right\} + \int_t^T \text{tr}\left\{\~{M}(s)\left[\int_t^s \~{\Phi}(s,\tau)G(\tau)Q(\tau)G^T(\tau)\~{\Phi^T}(s,\tau) d\tau\ \right] \right\} ds \\ &= x^T\left[ \~{\Phi}^T(T, t)F\~{\Phi}(T, t) + \int_t^T \~{\Phi}^T(s,t)\~{M}(s)\~{\Phi}(s,t) ds \right]x + \int_t^T \text{tr}\left\{G(s)Q(s)G^T(s)\~{\Phi}^T(T, s)F\~{\Phi}^T(T, s) \right\} ds + \int_t^T \text{tr}\left\{ G(s)Q(s)G^T(s) \left[ \int_s^T \~{\Phi}(\tau,s)\~{M}(\tau)\~{\Phi}(\tau,s) d\tau \right] \right\} ds \\ &= x^T\left[ \~{\Phi}^T(T, t)F\~{\Phi}(T, t) + \int_t^T \~{\Phi}^T(s,t)\~{M}(s)\~{\Phi}(s,t) ds \right]x + \int_t^T \text{tr}\left\{G(s)Q(s)G^T(s) \left[\~{\Phi}^T(T, s)F\~{\Phi}^T(T, s) + \int_s^T \~{\Phi}(\tau,s)\~{M}(\tau)\~{\Phi}(\tau,s) d\tau \right] \right\} ds \\ &= x^T\~{\Pi}(t)x + \~{\beta}(t) \tag{7.20} \end{aligned}$$

ここで、$\~{\Pi}(t)$と$\~{\beta}(t)$はそれぞれ以下のように定義しています。

$$\~{\Pi}(t) = \~{\Phi}^T(T, t)F\~{\Phi}(T, t) + \int_t^T \~{\Phi}^T(s,t)\~{M}(s)\~{\Phi}(s,t) ds \tag{7.21}$$

$$\~{\beta}(t) = \int_t^T \text{tr}\left\{ G(s)Q(s)G^T(s) \~{\Pi}(s)\right\}ds \tag{7.22}$$

ここまでの計算で、制御量$u(t)$が式(7.11)のようにフィードバック系であると仮定した時の、最小コスト汎関数$\~{V}(t,x)$は式(7.20)のように、状態量$x$の二次形式で与えられることがわかりました。

すごくかっこいい式展開なんですが、そもそも遷移マトリクス$\~{\Phi}(T.t)$を微分系からどう計算すればいいのかわからないぞ。

$\~{V}(t,x)$が二次形式で与えられることを手がかりにHJB方程式を解く

制御量$u(t)$が式(7.11)のようにフィードバック系であると仮定した時の、最小コスト汎関数$\~{V}(t,x)$は式(7.20)のように、状態量$x$の二次形式で与えられることがわかりました。

この事実に基づいて、式(7.10)の解を次のように仮定してみます。

$$V(t,x) = x^T\Pi(t)x+\alpha^T(t)x + \beta(t) \tag{7.22}$$

ここで、$\Pi(t)$は正定対称$n\times n$マトリクス、, $\alpha(t)$は$n$次元ベクトル、$\beta(t)$はスカラ量です。

$V(t, x)$の偏微分を求めると以下のようになります。

$$\frac{\partial V(t,x)}{\partial t} = x^T\dot{\Pi}(t)x + \dot{\alpha}(t)x + \dot{\beta}(t)$$

$$\frac{\partial V(t,x)}{\partial x} = 2\Pi(t)x + \alpha(t)$$

$$\frac{\partial}{\partial x} \frac{\partial V(t,x)}{\partial x} = 2\Pi(t)$$

これを式(7.10a)に代入して整理すると、以下の式(7.23)のようになります。

$$0 = x^T\left[\dot{\Pi}(t) + \Pi(t)A(t) + A^T(t)\Pi(t)+M(t)-\Pi(t)C(t)N^{-1}(t)C^T(t)\Pi(t)\right]x + x^T\left\{ \dot{\alpha}(t) + [A^T(t)-\Pi(t)C(t)N^{-1}(t)C^T(t)]\alpha(t) \right\} + \left[ \dot{\beta}(t)+\text{tr}\{G(t)Q(t)G^T(t)\Pi(t)\} - \frac{1}{4}\alpha^T(t)C(t)N^{-1}(t)C^T(t)\alpha(t)\right]$$

。

長いっ！

式(7.23)は、すべての$x\in R^n$に対して成り立たなければならないので、$Pi(t)$、$\alpha(t)$と$\beta(t)$は以下の連立微分方程式の解でなければなりません。

$$\begin{aligned} \dot{\Pi}(t) + \Pi(t)A(t) + A^T(t)\Pi(t)+M(t)-\Pi(t)C(t)N^{-1}(t)C^T(t)\Pi(t) &= 0 \\ \dot{\alpha}(t) + [A^T(t)-\Pi(t)C(t)N^{-1}(t)C^T(t)]\alpha(t) &= 0 \\ \dot{\beta}(t)+\text{tr}\{G(t)Q(t)G^T(t)\Pi(t)\} - \frac{1}{4}\alpha^T(t)C(t)N^{-1}(t)C^T(t)\alpha(t) &= 0 \end{aligned}$$

ここで、式(7.10b)で与えられた境界条件と式(7.22)で行った$V(x,t)$の形の仮定を見比べてみると、$\Pi(T) = F$、$\alpha(T) = 0$、$\beta(T) = 0$ということがわかります。

$$V(T,x) = x^T F x \tag{7.10b}$$

$$V(t,x) = x^T\Pi(t)x+\alpha^T(t)x + \beta(t) \tag{7.22}$$

先程の連立微分方程式のうち、$\alpha(t)$が絡んだ同次方程式は、$\alpha(t) = 0$が解となるので、式(7.10)の解は結局以下のようになります。

$$V(t,x) = x^T\Pi(t)x + \beta(t) \tag{7.24}$$

残った$\Pi(t)$と$\beta(t)$は、以下の連立微分方程式の解です。教科書によると、式(7.25)は単独で解けるそうです。

$$\dot{\Pi}(t) + \Pi(t)A(t) + A^T(t)\Pi(t)+M(t)-\Pi(t)C(t)N^{-1}(t)C^T(t)\Pi(t) = 0 \\ \Pi(T) = 0 \tag{7.25}$$

$$\dot{\beta}(t)+\text{tr}\{G(t)Q(t)G^T(t)\Pi(t)\} - \frac{1}{4}\alpha^T(t)C(t)N^{-1}(t)C^T(t)\alpha(t) = 0 \\ \beta(T) = 0 \tag{7.26}$$

$\Pi(t)$さえわかってしまえば、式(7.9)の最適制御量$u^o$は、以下の式(7.27)というフィードバック形で与えられることがわかります。

$$\begin{aligned} u^0(t) &= -\frac{1}{2}N^{-1}(t)C^T(t)\frac{\partial V(t,x)}{\partial x} \\ &= -\frac{1}{2}N^{-1}(t)C^T(t) 2\Pi(t)x\\ &= -N^{-1}(t)C^T(t) \Pi(t)x \tag{7.27} \end{aligned}$$

コストの最小値を求める

最適制御量がわかったので、今度は評価汎関数$J(u)$の最小値を求めていきます。

$$\begin{aligned} J(u^o) &= \mathcal{E}\left\{ \min_{u(t),t_0 \le t \le T} \mathcal{E}\left\{ x^T(T)F x(T) + \int_{t_0}^T [x^T(t) M(t)x(t)+u^T(t)N(t)u(t)]dt \bigg| x(t_0) = x_0\right\} \right\} \\ &= \mathcal{E}\{V(t_0, x_0)\} \\ &= \mathcal{E}\{x_0^T\Pi(t_0)x_0\} + \beta(t_0) \end{aligned} \tag{7.28}$$

ここはこれまでの定義通り式展開しているだけ。

ここで、

$$\begin{aligned} \mathcal{E}\{x_0^T\Pi(t_0)x_0\} &= \text{tr}\{\Pi(t_0)\mathcal{E}\{x_0x_0^T\}\}\\ &= m_0^T\Pi(t_0)m_0+\text{tr}\{P_0\Pi(t_0)\} \end{aligned} \tag{7.29}$$

教科書でも突如$\mathcal{E}\{x_0 x_0^T\} = m_0 m_0^T + P_0$という定義が出てきていて、ちょっと意味がわからない。

さらに、式(7.26)を積分して、$\beta(T) = 0$に注意すると

$$\beta(t_0) = \int_{t_0}^T \text{tr}\{G(t)Q(t)G^T(t)\Pi(t)\}dt \tag{7.30}$$

となり、最終的に最小コストは式(7.31)の通りとなります。

$$J(u^o) = m_0^T\Pi(t_0)m_0+\text{tr}\{P_0\Pi(t_0)\} + \int_{t_0}^T \text{tr}\{G(t)Q(t)G^T(t)\Pi(t)\}dt \tag{7.31}$$

このあたりもちょっと展開が急です。積分定数が0になることを言いたかったのかな？

以上により、線形確率システムに対する最適制御が求まりました。教科書にあった確定システムと確率システムの比較表は以下のようなものでした。

	確定システム	確率システム
システムモデル	$$\dot{x} = Ax + Cu$$	$$dx = Axdt + Cudt + Gdw$$
コスト汎関数	$$J(u) = x^T(T)Fx(T) + \int_{t_0}^T(x^TMx+u^TNu)dt$$	$$J(u) = \mathcal{E}\left\{x^T(T)Fx(T)+\int_{t_0}^T (x^TMx+u^TNu)dt\right\}$$
最小コスト汎関数	$$V(t,x) = x^T\Pi(t)x\\ \dot{\Pi}+\Pi A+A^T\Pi+M-\Pi CN^{-1}C^T\Pi = 0\\ \Pi(T) = F$$	$$V(t,x) = x^T\Pi(t)x+\beta(t)\\ \dot{\Pi}+\Pi A+A^T\Pi+M-\Pi CN^{-1}C^T\Pi = 0\\ \Pi(T) = F\\ \dot{\beta} + \text{tr}\{GQG^T\Pi\} = 0\\ \beta(T) = 0$$
最適制御量	$$u^o(t) = -N^{-1}C^T\Pi(t)x(t)$$	$$u^o(t) = -N^{-1}C^T\Pi(t)x(t)$$
最小コスト値	$$J(u^o) = x_0^T\Pi(t_0)x_0$$	$$J(u^o) = m_0^T\Pi(t_0)m_0+\text{tr}\{P_0\Pi(t_0)\}+\int_{t_0}^T\text{tr}\{GQG^T\Pi\}dt$$

式(7.25)を見るとわかるとおり、$\Pi(t)$の決定にはコスト汎関数の重みマトリクス$F, M(t), N(t)$が関与しますが、ノイズに関するパラメータマトリクス$G(t), Q(t)$には依存しません。

つまり、最適制御量$u^o(t)$を生成するゲインマトリクス$K(t) = -N^{-1}(t)C^T(t)\Pi(t)$は、システムの雑音の大きさには全く関係なく求めることができます。システム雑音があってもなくても制御ゲインマトリクスが同じである、という事実は重要です。これを確実等価性原理と呼びます。しかし、コスト汎関数の最小値については、システム雑音の影響が含まれます。

ポイントになると感じたのは以下の2点でした。

• 最適制御を与える制御量が確定システムと確率システムで同一であるということ
• 同じ制御量を加えても最小コスト値は確率システムのほうが雑音の分大きくなっている

結構時間がかかりましたが、ここで最適制御までの一連の教科書のレビューはいったん終わりにしておきます。マーケットメイキング周りの論文ではマーケットを簡素なモデリングで表現し、そこに最適制御を使うことが多くあるようです。ただ、コスト汎関数を二次形式ではなくて指数形式にすることがほとんどで、この教科書のやり方そのままではHJB方程式の導出や求解をすることはできません。

まだまだ勉強しなければいけないことは山のようにあるのですが…コードを書かない状態が続くのはよくないので、Bot作りと勉強を並行して進めていこうと思います。