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

やること

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

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

Botter自主ゼミノート 1.2 確定システムの制御の回顧
Botter自主ゼミノート 2.1 確率過程とは？, 2.2 確率過程の数学的表現

リッカチ微分方程式の導出にチャレンジしてみよう

まずはHJB方程式が与えられたものとして、そこからの導出をこちらのページをなぞって写経してみます。

最適制御入力を求める

ハミルトン関数の定義は以下の通りです。システムの状態方程式$\.{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u},t)$に対してラグランジュの未定乗数法を使った結果、$\lambda^T$以降の項が出現して、ごちゃごちゃしたものを$H$にまとめています。そのため、$H$の中には$\lambda$が入ってきています。

$$H(x, \lambda, u, t) = \frac {1}{2} \left( x^T Q(t) x + u^T R(t) u \right) + \lambda^T(Ax + Bu)$$

$R(t)$が正定値行列なら$H$を最小化する$u$は一意に存在し、$\frac {\partial H}{\partial u} = 0$を満たすのだそうです。

• $n \times n$実対称行列$M$が正定値 (positive definite) であるとは、$n$個の実数を成分に持つ零ベクトルでない任意の列ベクトル$z$に対して、二次形式$z^T M z$が必ず正となるということです

$H(x, \lambda, u, t)$の最小値を与える$u$を計算してみます。$H(x, \lambda, u, t)$を$u$で微分して、その値が$0$になる$u$を探します。その際に、以下を利用します。

• 対称行列で表された二次形式の微分 (勾配ベクトル) は$x^T A x = 2 A x$となります (高校数学の美しい物語)

$$\begin{aligned} H(x, \lambda, u, t) &= \frac {1}{2} \left( x^T Q(t) x + u^T R(t) u \right) + \lambda^T(Ax + Bu) \\ \frac{\partial}{\partial u} H(x, \lambda, u, t) &= \frac{\partial}{\partial u} \frac{1}{2} u^T R(t) u + \frac{\partial}{\partial u} \lambda^T Bu\\ &= R(t) u + \lambda^T B \\ \end{aligned}$$

右辺が$0$になればいいので…

$$\begin{aligned} 0 &= R(t) u + \lambda^T B \\ u &= -R^{-1}(t) \lambda^T B \\ \end{aligned}$$

$\lambda$は縦ベクトルで$B$は行列なので、転置して入れ替えても結果は変わりません。

$$\therefore u = - R^{-1} (t) B^T \lambda$$

• 岡山大学線形システム制御最適化特論の授業資料の42ページでは平方完成を使ってこの結果を得ています。
• 対称行列で表された二次形式$x^TAx + x^Tb + c$を平方完成すると、以下のような式になります (Wikipedia)

$$x^T A x + x^T b + c = (x-h)^TA(x-h)+k\hspace{1em}\left(h = \frac{1}{2}A^{-1}b, k=c-\frac{1}{4}b^T A^{-1}b\right)$$

$H(x, \lambda, u, t)$の最小値を与える$u$が得られたので、これをハミルトンヤコビベルマン方程式の$min$の項に代入します。

$$\begin {aligned} \min_{u} H\left(x, \lambda, u, t \right) + \frac{\partial V}{\partial t} \left(x,t\right) &= 0 \hspace{1em} \left( \lambda = \frac{\partial V}{\partial x}^T \right) \\ -\frac{\partial V}{\partial t}(x,t) &= \frac {1}{2} x^T Q(t) x + \frac {1}{2} \left(- R^{-1} (t) B^T \lambda \right)^{T} R(t) \left(- R^{-1} (t) B^T \lambda \right) + \lambda^T Ax + \lambda^T B \left(- R^{-1} (t) B^T \lambda \right) \\ &= \frac {1}{2} x^T Q(t) x + \frac {1}{2} \lambda^T B \left(R^{-1}(t) \right)^{T} R(t) R^{-1} (t) B^T \lambda + \lambda^T Ax - \lambda^T B R^{-1} (t) B^T \lambda \end {aligned}$$

ここで、このページの導出では

$$\frac 12 \left(\lambda^TBR(t)^{-T}\right)R(t)\left(R(t)^{-1}B^T\lambda \right)$$

の項が

$$\frac {1}{2} \lambda^T B R(t)^{-1} B^T \lambda$$

に変形されているけれど、なぜこうなるのかがよくわかりません。$R(t)^{-1}$ではなくて$R(t)^{-T}$なら分かりやすいのですが…。この$R(t)^{-1}$はこの後もずっと生き残っているし、何か行列計算のベーシックな部分で躓いているような気がします。

ゼミ生からの指摘
その変形に関しては、いま$R(t)$は正定値行列なので対称行列であり、また正則なので逆行列$R(t)^{-1}$が存在して、対称行列の逆行列は存在すれば対称行列なので、$R(t)^{-1}$も対称行列、つまり$R(t)^{-1} = (R(t)^{-1})^T$ということかと思います!

なるほど！　行列計算のベーシックな所を忘れていたので、対称行列についての性質を確認しておきます。

対称行列の性質の復習

• 対称行列の固有値はすべて実数
• 対称行列の固有ベクトルは直交する
• 対称行列は必ず対角化できる
• 対称行列の逆行列は対象行列

最後の特性の証明は以下の通りです。

$$\begin{aligned} A A^{−1} &= I \\ (A A^{-1})^T &= I^T \\ (A^{-1})^T A^T &= I \\ (A^{-1})^T A &= I \hspace{10mm} (A^T = Aを利用) \\ (A^{-1})^T &= A^{-1} \hspace{10mm} (A^{-1}を両辺の右から乗算し、A A^{-1} = Iを利用) \\ \end{aligned}$$

では気を取り直して再開です。

$$\begin {aligned} \min_{u} H\left(x, \lambda, u, t \right) + \frac{\partial V}{\partial t} \left(x,t\right) &= \frac {1}{2} x^T Q(t) x + \frac {1}{2} \lambda^T B \left(R^{-1}(t) \right)^{T} R(t) R^{-1} (t) B^T \lambda + \lambda^T Ax - \lambda^T B R^{-1} (t) B^T \lambda \\ &= \frac {1}{2} x^T Q(t) x + \frac {1}{2} \lambda^T B R^{-1}(t) B^T \lambda + \lambda^T Ax - \lambda^T B R^{-1} (t) B^T \lambda \\ &= \frac {1}{2} x^T Q(t) x + \lambda^T Ax - \frac {1}{2} \lambda^T B R^{-1}(t) B^T \lambda\\ \end {aligned}$$

ここで、$V(x,t) = \frac{1}{2} x^T P(t) x$と表せると仮定すると上手くいくことが知られているそうなのでやってみます。

$\lambda = \frac {\partial V^T}{\partial x} = P(t)x$ということは、$\lambda^T = x^T P^T (t)$なので…
(行列$P$,$Q$,$R$の引数は省略します)

$$\begin{aligned} -\frac{\partial V}{\partial t}(x,t) &= \frac {1}{2} x^T Q x + x^T P^T A x - \frac {1}{2} x^T P^T B R^{-1} B^T Px\\ &= \frac 12 x^T\left\{ Q + 2P^T A - P^T B R^{-1} B^T P \right\}x\\ \end{aligned}$$

ここで$V(x,t)$を二次形式で表していて係数行列が$P(t)$なので、$P(t)$は対称行列で、$P(t) = P(t)^T$が成立します。これを利用して式変形をします。

$$-\frac{\partial V}{\partial t}(x,t) = \frac 12 x^T\left\{ Q + 2P^T A - P B R^{-1} B^T P \right\}x$$

この時、xは縦ベクトルなので、$x^T 2P^T A x$は必ずスカラーになります。ですのでこんな変形ができてしまいそうです。

$$x^T 2P^T A x = x^T P^T A x + x^T P^T A x$$

スカラーを転置しても同じスカラーなので

$$\begin{aligned} &= x^T P^T A x + \left(x^T P^T A x\right)^T\\ &= x^T P^T A x + x^T A^T P x \end{aligned}$$

前述の通り$P = P^T$なので

$$\begin{aligned} &= x^T P A x + x^T A^T P x\\ &= x^T \left(PA + A^TP\right) x \end{aligned}$$

この変形を上式に加えると…

$$-\frac{\partial V}{\partial t}(x,t) = \frac{1}{2} x^T\left\{ Q + PA + A^TP - P B R^{-1} B^T P \right\}x$$

ここで、左辺は$-V(x,t) = -\frac{1}{2} x^T P(t) x$を$t$で偏微分しているので、$-\frac{1}{2} x^T \.{P(t)} x$となり、左辺と右辺の形が一致します。

$$\begin{aligned} -\frac{1}{2} x^T \.{P(t)} x &= \frac{1}{2} x^T\left\{ Q + PA + A^TP - P B R^{-1} B^T P \right\}x\\ -\.{P(t)} &= Q + PA + A^TP - P B R^{-1} B^T P \end{aligned}$$

ここまでで、線形システムの最適制御則が次の手順で求まることがわかりました。

$$u = - R^{-1} (t) B^T \lambda \\ -\.{P(t)} = Q + PA + A^TP - P B R^{-1} B^T P\\$$

この微分方程式はリッカチ微分方程式と呼ばれていて、$P(t_f) = Q_f$のもとで逆の時間方向に解くことで最適な制御入力$u(t)$を得ることができるらしいです。ほんと？

ここで、今回の導出結果ともともと教科書に書いてあったことを見比べると…以下のように対応していることがわかります。

$$\begin{aligned} P(t) &= \Pi(t) \\ Q_f &= F \\ R(t)^{-1} &= N^{-1} \\ \lambda &= \Pi(t) x(t) \end{aligned}$$

$y(t)$，$x(t)$はそれぞれ状態量および観測量ベクトル，$u(t)$は制御量ベクトルであり，$x \in {R^n}$，$y \in {R^m} (m \le n)$，$u \in R^{\ell}$とする．

マトリクス$H$が単位マトリクス，$H = I$であるならばシステム状態量はすべて独立に観測できることから，最適制御は評価コスト汎関数を２次形式

$$J(u) = x^T(T)Fx(T)+\int_{t_0}^{T}\left[x^T (s) Mx(s) + u^T (s)Nu(s)\right]ds \tag{1.2}$$

(ただし，$F$，$M \ge 0$，$N > 0$は対象マトリクス; 上付きの$^T$はベクトルの転置を表す)と設定すると，

$$u^o(t) = -N^{-1}C^T \Pi (t)x(t) \tag{1.3}$$

で与えられる．ここで$\Pi(t)$は$n \times n$次元リッカチ微分方程式

$$\left. \begin{aligned} \.{\Pi}(t) + \Pi (t)A + A^T \Pi (t) + M - \Pi (t) CN^{-1} C^T \Pi(t) &= 0 \\ \Pi (T) &= F \end{aligned} \quad \right\} \tag{1.4}$$

の正定値解である。