📔

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

2022/11/16に公開

やること

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を満たすのだそうです。

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

\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

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}

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

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

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

\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}
-\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
-\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) xtで偏微分しているので、-\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}

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

\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であるならばシステム状態量はすべて独立に観測できることから,最適制御は評価コスト汎関数を2次形式

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}

(ただし,FM \ge 0N > 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}

の正定値解である。

Discussion