離散時間システムの最適制御解　ー脳の設計図を求めて(5)ー

これまで、状態変数、制御変数のいずれも時間に対して連続な場合の状態空間モデルについて議論してきました。しかし、これらの変数を離散的な時刻における代表量を使って取り扱う方が便利な場合があります。例えば、微分方程式の数値解法においては、連続的な微分方程式を離散的な時刻における数値の関係を表す差分方程式で近似します。離散誤差を制御することで、時間に関して高次精度の公式を構成することも可能です。また、デジタル計算機や人工神経回路網などのように、その動作が本質的に離散的なシステムも存在します。脳の中にある神経の数は非常に多いと言っても有限ですから、脳も離散的なシステムです。

そこで本章では、時間が離散的な場合、最適制御解法を使って到達問題を解く手法を記述します。まず、離散状態空間モデルを定義し、非線形な場合でも局所線形化によって到達問題の最適解が同様に得られることを示します。

1. 離散時間システム

離散時間システムが状態方程式

\begin{equation}
{\bf z}^{(k+1)}={\bf z}^{(k)}+{\bf f}({\bf z}^{(k+1/2)}, {\bf u}^{(k)},k) \tag{1}
\end{equation}

で表されているとします。ここで\({\bf z}^{(k)}∈R^{n}\)は状態変数ベクトル、\({\bf u}^{(k)}∈R^{m}\)は制御変数ベクトルで、

\begin{equation}
{\bf z}^{(k+1/2)}=({\bf z}^{(k)}+{\bf z}^{(k+1)})/2 \tag{2}
\end{equation}
です。ここで、時間に関して二次精度の公式の構築のために、時間ステップ\(k\)段と時間ステップ\(k+1\)段の中間値\({\bf z}^{(k+1/2)}\)にて\({\bf f}\)を評価しています。\({\bf u}^{(k)}\)については、もともと\(k+1/2\)ステップにて評価されている量であるとします。初期段\(k=0\)、最終段\(k=K\)、初期状態\({\bf z}^{(0)}={\bf z}_0\)、最終状態\({\bf z}^{(K)}={\bf z}_{\rm f}\)が与えられ、状態方程式のみが拘束条件であるとします。ここで、最小化すべき評価関数を、

\begin{equation}
J=\Sigma_{k=0}^{K-1}L({\bf z}^{(k)},{\bf u}^{(k)},k) \tag{3}
\end{equation}

で与えると、状態方程式を満たす拘束条件の下で停留条件を求める汎関数\(\bar{J}\) は、
\begin{equation}
\bar{J} =φ({\bf z}^{(k)},K)+\Sigma_{k=0}^{K-1}(L({\bf z}^{(k)},{\bf u}^{(k)},k)+{{\bf λ}^{(k+1)}}^{\rm T} ({\bf f}({\bf z}^{(k+1/2)},{\bf u}^{(k)},k)-{\bf z}^{(k+1)} )) \tag{4}
\end{equation}
と構成できます。ハミルトン関数を
\begin{equation}
H({\bf z}^{(k)},{\bf u}^{(k)},{\bf λ}^{(k+1)},k)=L({\bf z}^{(k)},{\bf u}^{(k)},k)+{{\bf λ}^{(k+1)}}^{\rm T} {\bf f}({\bf z}^{(k+1/2)},{\bf u}^{(k)},k) \tag{5}
\end{equation}
と定義すると、オイラー・ラグランジュ方程式は、以下のように与えられます。
\begin{equation}
0=\frac{∂H({\bf z}^{(k)},{\bf u}^{(k)},{\bf λ}^{(k+1)},k)}{∂{\bf z}^{(k)}}-{{\bf λ}^{(k)}} \tag{6}
\end{equation}

\begin{equation}
0=\frac{∂H({\bf z}^{(k)},{\bf u}^{(k)},{\bf λ}^{(k+1)},k)}{∂{\bf u}^{(k)}} \tag{7}
\end{equation}

練習問題5.1
評価関数として式2を使ったとき、オイラー・ラグランジュ方程式として式5と6が得られることを示しなさい。

2. 非線形システムの局所線形化

非線形離散時間システムの状態方程式と評価関数を局所線形化します。この両者が、\(k\)には直接依存しないと仮定します。すなわち、
\begin{equation}
{\bf z}^{(k+1)}={\bf z}^{(k)}+{\bf f}({\bf z}^{(k+1/2)},{\bf u}^{(k)}) \tag{8}
\end{equation}
とします。式8を\({\bf z}^{(k)}=\bar{\bf z}^{(k)}+∆{\bf z}^{(k)}\)、\({\bf u}^{(k)}=\bar{\bf u}^{(k)}+∆{\bf u}^{(k)}\)と置いて\(∆\)付きの量の一次までの項を取ると
\begin{eqnarray}
\bar{\bf z}^{(k+1)}+∆{\bf z}^{(k+1)}&=&\bar{\bf z}^{(k)}+∆{\bf z}^{(k)}+{\bf f}(\bar{\bf z}^{(k+1/2)},\bar{\bf u}^{(k)})\\
&+&\tilde{\bf A}^{(k)} \frac{∆{\bf z}^{(k+1)}+∆{\bf z}^{(k)}}{2}+\tilde{\bf B}^{(k)} {∆\bf u}^{(k)} \tag{9}
\end{eqnarray}
と書けます。ここで、\( \tilde{\bf A}^{(k)}=\left. ∂{\bf f}/∂{{\bf z}^{(k)}}^{\rm T}\right|_{\bar{\bf z}^{(k+1/2)},\bar{\bf u}^{(k)}} \)、\(\tilde{\bf B}^{(k)}=\left. ∂{\bf f}/∂{{\bf u}^{(k)}}^{\rm T}\right|_{\bar{\bf z}^{(k+1/2)},\bar{\bf u}^{(k)}}\)とします。したがって、
\begin{equation}
({\bf I}_n-\frac{1}{2}\tilde{\bf A}^{(k)}){∆\bf z}^{(k+1)}=({\bf I}_n+\frac{1}{2}\tilde{\bf A}^{(k)} ){∆\bf z}^{(k)}+\tilde{\bf B}^{(k)} {∆\bf u}^{(k)}+\tilde{\bf c}^{(k)} \tag{10}
\end{equation}
と書けます。ここで、\(\tilde{\bf c}^{(k)}=-\bar{\bf z}^{(k+1)}+\bar{\bf z}^{(k)}+{\bf f}(\bar{\bf z}^{(k+1/2)},\bar{\bf u}^{(k)})\)としました。また、\(\tilde{\bf A}^{(k)} \)も、\(\tilde{\bf B}^{(k)}\)も\(\tilde{\bf c}^{(k)}\)も、\(\bar{\bf z}^{(k)}\)と\(\bar{\bf u}^{(k)}\)の関数であるが、\(k\)には直接は関係しないとします。式10は
\begin{equation}
{∆\bf z}^{(k+1)}={\bf A}^{(k)} {∆\bf z}^{(k)}+{\bf B}^{(k)} {∆\bf u}^{(k)}+{\bf c}^{(k)} \tag{11}
\end{equation}
と書き直せます。ここで、

\begin{equation}
{\bf A}^{(k)}=({\bf I}_n-\frac{1}{2} \tilde{\bf A}^{(k)})^{-1} ({\bf I}_n+\frac{1}{2}\tilde{\bf A}^{(k)} ) \tag{12}
\end{equation}
\begin{equation}
{\bf B}^{(k)}=({\bf I}_n-\frac{1}{2} \tilde{\bf A}^{(k)})^{-1} \tilde{\bf B}{(k)} \tag{13}
\end{equation}
\begin{equation}
{\bf c}^{(k)}=({\bf I}_n-\frac{1}{2} \tilde{\bf A}^{(k)})^{-1}\tilde{\bf c}^{(k)} \tag{14}
\end{equation}
です。

また、評価関数\(L({\bf z}^{(k)},{\bf u}^{(k)})\)も一般には\({\bf z}^{(k)}\)と\(u^{(k)}\)の二次形式にはなりません。そこで、\({∆\bf z}^{(k)}\)と\({∆\bf u}^{(k)}\)についてテイラー展開して二次の項までとると、
\begin{eqnarray}
L({\bf z}^{(k)},{\bf u}^{(k)})&=&L(\bar{\bf z}^{(k)},\bar{\bf u}^{(k)} )+\left.\frac{∂L}{∂{{\bf z}^{(k)}}^{\rm T}}\right|_{\bar{\bf z}^{(k)},\bar{\bf u}^{(k)}}{∆\bf z}^{(k)} + \left.\frac{∂L}{∂{{\bf u}^{(k)}}^{\rm T}}\right|_{\bar{\bf z}^{(k)},\bar{\bf u}^{(k)}}{∆\bf u}^{(k)}\\
&+&\frac{1}{2}\left({{∆\bf z}^{(k)}}^{\rm T}\left.\frac{∂^2 L}{∂{{\bf z}^{(k)}}^2 }\right|_{\bar{\bf z}^{(k)},\bar{\bf u}^{(k)} } {∆\bf z}^{(k)}
+{{∆\bf u}^{(k)}}^{\rm T}\left.\frac{∂^2 L}{{∂{\bf u}^{(k)}}^{\rm T} }\right|_{\bar{\bf z}^{(k)},\bar{\bf u}^{(k)}}{∆\bf u}^{(k)}\right)\\
&=&L(\bar{\bf z}^{(k)},\bar{\bf u}^{(k)} )+{{\bf p}^{(k)}}^{\rm T}{∆\bf z}^{(k)} + {{\bf q}^{(k)}}^{\rm T}{∆\bf u}^{(k)} \\
&&+\frac{1}{2}\left({{∆\bf z}^{(k)}}^{\rm T}{\bf Q}^{(k)}{∆\bf z}^{(k)}+{{∆\bf u}^{(k)}}^{\rm T}{\bf R}^{(k)}{∆\bf u}^{(k)}\right)\\
&=&L(\bar{\bf z}^{(k)},\bar{\bf u}^{(k)} )+\frac{1}{2}\left({{∆\bf z}^{(k)}}^{\rm T}{\bf Q}^{(k)}{∆\bf z}^{(k)}+2{{\bf p}^{(k)}}^{\rm T}{∆\bf z}^{(k)} \right)\\
&&+\frac{1}{2}\left({{∆\bf u}^{(k)}}^{\rm T}{\bf R}^{(k)}{∆\bf u}^{(k)}+2{{\bf q}^{(k)}}^{\rm T}{∆\bf u}^{(k)}\right)\\
&=&L(\bar{\bf z}^{(k)},\bar{\bf u}^{(k)} )\\
&&+\frac{1}{2}{\left({∆\bf z}^{(k)}+{{\bf Q}^{(k)}}^{-1}{\bf p}^{(k)}\right)}^{\rm T}{\bf Q}^{(k)}\left({∆\bf z}^{(k)}+{{\bf Q}^{(k)}}^{-1}{\bf p}^{(k)}\right)\\
&&-\frac{1}{2}{{\bf p}^{(k)}}^{\rm T}{{\bf Q}^{(k)}}^{-1}{\bf p}^{(k)}\\
&&+\frac{1}{2}{\left({∆\bf u}^{(k)}+{{\bf R}^{(k)}}^{-1}{\bf q}^{(k)}\right)}^{\rm T}{\bf R}^{(k)}\left({∆\bf u}^{(k)}+{{\bf R}^{(k)}}^{-1}{\bf q}^{(k)}\right)\\
&&-\frac{1}{2}{{\bf q}^{(k)}}^{\rm T}{{\bf R}^{(k)}}^{-1}{\bf q}^{(k)}\\
&=&L(\bar{\bf z}^{(k)},\bar{\bf u}^{(k)} )\\
&&+\frac{1}{2}{\left({∆\bf z}^{(k)}+{{\bf Q}^{(k)}}^{-1}{\bf p}^{(k)}\right)}^{\rm T}{\bf Q}^{(k)}\left({∆\bf z}^{(k)}+{{\bf Q}^{(k)}}^{-1}{\bf p}^{(k)}\right)\\
&&+\frac{1}{2}{\left({∆\bf u}^{(k)}+{{\bf R}^{(k)}}^{-1}{\bf q}^{(k)}\right)}^{\rm T}{\bf R}^{(k)}\left({∆\bf u}^{(k)}+{{\bf R}^{(k)}}^{-1}{\bf q}^{(k)}\right)\\
&&+L(\bar{\bf z}^{(k)},\bar{\bf u}^{(k)} )
-\frac{1}{2}{{\bf p}^{(k)}}^{\rm T}{{\bf Q}^{(k)}}^{-1}{\bf p}^{(k)}
-\frac{1}{2}{{\bf q}^{(k)}}^{\rm T}{{\bf R}^{(k)}}^{-1}{\bf q}^{(k)}\tag{15}
\end{eqnarray}
と変形できます。ここで、\({\bf Q}^{(k)} (\bar{\bf z}^{(k)},\bar{\bf u}^{(k)})=\left.\frac{∂^2 L}{{∂{\bf z}^{(k)}}^2}\right|_{\bar{\bf z}^{(k)},\bar{\bf u}^{(k)}}\)、\({\bf R}^{(k)}(\bar{\bf z}^{(k)},\bar{\bf u}^{(k)})=\left. \frac{∂^2 L}{∂{{\bf u}^{(k)}}^2}\right|_{\bar{\bf z}^{(k)},\bar{\bf u}^{(k)}}\)、\({\bf p}^{(k)}(\bar{\bf z}^{(k)},{\bf u}^{(k)})=\left. \frac{∂L}{∂{{\bf z}^{(k)}}^{\rm T}}\right|_{\bar{\bf z}^{(k)},\bar{\bf u}^{(k)}}\)、\({\bf q}^{(k)}({\bf z}^{(k)},\bar{\bf u}^{(k)})=\left.\frac{∂L}{∂{{\bf u}^{(k)}}^{\rm T}}\right|_{\bar{\bf z}^{(k)},\bar{\bf u}^{(k)}}\)です。ただし、\(\frac{∂^2 L}{∂{\bf z}^{(k)} ∂{\bf u}^{(k)}}=0\)を仮定しました。また、\({\bf Q}^{(k)}={{\bf Q}^{(k)}}^T\)より、\({{{\bf Q}^{(k)}}^{-1}}^{\rm T}={{\bf Q}^{(k)}}^{-1}\)であることを使いました。また、
\begin{equation}
{\bf z}_{\rm m}^{(k)}=-{{\bf Q}^{(k)}}^{-1}{\bf p}^{(k)} \tag{16}
\end{equation}
\begin{equation}
{\bf u}_{\rm m}^{(k)}=-{{\bf R}^{(k)}}^{-1}{\bf q}^{(k)} \tag{17}
\end{equation}
とおけば、
\begin{eqnarray}
L({∆\bf z}^{(k)},{∆\bf u}^{(k)} )&=&\frac{1}{2}\left({∆\bf z}^{(k)}-{\bf z}_{\rm m}^{(k)} \right)^{\rm T}{\bf Q}^{(k)} \left({∆\bf z}^{(k)}-{\bf z}_{\rm m}^{(k)} \right)\\
&+&\frac{1}{2} \left({∆\bf u}^{(k)}-{\bf u}_{\rm m}^{(k)}\right)^{\rm T}{\bf R}^{(k)} \left({∆\bf u}^{(k)}-{\bf u}_{\rm m}^{(k)} \right) +L_0 (\bar{\bf z}^{(k)},\bar{\bf u}^{(k)}) \tag{18}
\end{eqnarray}

と\({∆\bf z}^{(k)}\)と\({∆\bf u}^{(k)}\)に関する二次形式で書き表されます。ここで、\({∆\bf z}^{(k)}\)にも\({∆\bf u}^{(k)}\)にも依存しない項\(L_0 (\bar{\bf z},\bar{\bf u})\)は
\begin{equation}
L_0 (\bar{\bf z}^{(k)},{\bf u}^{(k)})=L(\bar{\bf z}^{(k)},\bar{\bf u}^{(k)} )
-\frac{1}{2}{{\bf p}^{(k)}}^{\rm T}{{\bf Q}^{(k)}}^{-1}{\bf p}^{(k)}
-\frac{1}{2}{{\bf q}^{(k)}}^{\rm T}{{\bf R}^{(k)}}^{-1}{\bf q}^{(k)} \tag{19}
\end{equation}
です。このときオイラー・ラグランジュ方程式は、
\begin{equation}
{\bf λ}^{(k)}={\bf Q}^{(k)} \left({∆\bf z}^{(k)}-{\bf z}_{\rm m}^{(k)} \right)+{{\bf A}^{(k)}}^{\rm T} {\bf λ}^{(k+1)} \tag{20}
\end{equation}
\begin{equation}
{{\bf R}^{(k)}}^{\rm T} \left({∆\bf u}^{(k)}-{\bf u}_{\rm m}^{(k)}\right)+{{\bf B}^{(k)}}^{\rm T} {\bf λ}^{(k)}=0 \tag{21}
\end{equation}
となります（練習5-1）。

3. 離散時間天馬方程式

初期条件\({\bf z}^{(0)}={\bf z}_0\)、終端条件\({\bf z}^{(K)}={\bf z}_{\rm f}\)を満たし評価関数を最小とする系列\({\bf z}^{(k)}\)は、以下の計算手順で得られます。

手順０）初期条件\({\bf z}^{(0)}={\bf z}_0\)、および終端条件\({\bf z}^{(K)}={\bf z}_{\rm f}\)を満たす離散系列\(\bar{\bf z}^{(k)} (k=[0,K])\)
を仮に与える。4つの行列\({\bf A}^{(k) }\)、\({\bf B}^{(k)}\)、\({\bf Q}^{(k)}\)、\({\bf R}^{(k)}\)と二つのベクトル\({\bf d}^{(k)}\)、\({\bf e}^{(k)}\)は、\(\bar{\bf z}^{(k)}\)の関数として与えられている。

手順１）以下の4つの微分方程式を、\({\bf S}^{(K)}=0\)、\({\bf U}^{(K)}={\bf I}_n\)、\({\bf V}^{(K)}={\bf I}_n\)、\({\bf W}^{(K)}=0\)、\({\bf d}^{(K)}=0\)、\({\bf e}^{(K)}=0\)の終端条件から、式23、24、25、26、27、28の漸化式を、\(k\)が減る方向に適用して、\(k=[0,K]\)の区間における行列\({\bf S}^{(k)}\)、\({\bf U}^{(k)}\)、\({\bf V}^{(k)}\)、\({\bf W}^{(k)}\)、\({\bf d}^{(k)}\)、\({\bf e}^{(k)}\)を求める。
\begin{equation}
{\bf S}^{(k)}={{\bf A}^{(k)}}^{\rm T} \left({\bf I}_n-{\bf S}^{(k+1)} {\bf B}^{(k)} {\bf O}^{(k)} {{\bf B}^{(k)}}^{\rm T}\right){\bf S}^{(k+1)} {\bf A}^{(k)}+{\bf Q}^{(k)} \tag{22}
\end{equation}
\begin{equation}
{\bf U}^{(k)}={{\bf A}^{(k)}}^{\rm T}\left({\bf I}_n-{\bf S}^{(k+1)} {\bf B}^{(k)} {\bf O}^{(k)} {{\bf B}^{(k)}}^{\rm T} \right) {\bf U}^{(k+1)} \tag{23}
\end{equation}
\begin{equation}
{\bf V}^{(k)}={\bf V}^{(k+1)} \left({\bf I}_n-{\bf B}^{(k)} {\bf O}^{(k)} {{\bf B}^{(k)}}^T {\bf S}^{(k+1)} \right) {\bf A}^{(k)} \tag{24}
\end{equation}
\begin{equation}
{\bf W}^{(k)}={\bf W}^{(k+1)}-{\bf V}^{(k+1)} {\bf B}^{(k)} {\bf O}^{(k)}{{\bf B}^{(k)}}^{\rm T} {\bf U}^{(k+1)} \tag{25}
\end{equation}
\begin{eqnarray}
{\bf d}^{(k)}&=&{{\bf A}^{(k)}}^{\rm T}\left({\bf I}_n-{\bf S}^{(k+1)} {\bf B}^{(k)} {\bf O}^{(k)}{{\bf B}^{(k)}}^{\rm T} \right) \left({\bf d}^{(k+1)}+{\bf S}^{(k+1)} {\bf c}^{(k)}\right)\\
&+&{{\bf A}^{(k)}}^{\rm T} {\bf S}^{(k+1)} {\bf B}^{(k)} {\bf O}^{(k)} {\bf R}^{(k)} {\bf u}_{\rm m}^{(k)}-{\bf Q}^{(k)} {\bf z}_{\rm m}^{(k)} \tag{26}
\end{eqnarray}
\begin{eqnarray}
{\bf e}^{(k)}&=&{\bf e}^{(k+1)}-{\bf V}^{(k+1)} {\bf B}^{(k)} {\bf O}^{(k)} \left({{\bf B}^{(k)}}^{\rm T}\left( {\bf d}^{(k+1)}+{\bf S}^{(k+1)} {\bf c}^{(k)}\right)-{\bf R}^{(k)} {\bf u}_{\rm m}^{(k)} \right)
\\&+&{\bf V}^{(k+1)} {\bf c}^{(k)} \tag{27}
\end{eqnarray}
ここで、
\begin{equation}
{\bf O}^{(k)}=\left({\bf R}^{(k)}+{{\bf B}^{(k)}}^{\rm T} {\bf S}^{(k+1)} {\bf B}^{(k)} \right)^{-1}
\end{equation}

とします。

手順2）式 29、30、31 を使って、\(\boldsymbol{\nu}\) および区間 \(k=[0,K]\) の \(\Delta \mathbf{z}^{(k)}\) および \(\Delta \mathbf{u}^{(k)}\) を求める。

\[ \boldsymbol{\nu} = – {\mathbf{W}^{(0)}}^{-1} \left( \mathbf{V}^{(0)} \Delta \mathbf{z}^{(0)} + \mathbf{e}^{(0)} – \Delta \mathbf{z}^{(K)} \right) \tag{29} \]

\[ \Delta \mathbf{z}^{(k)} = – {\mathbf{V}^{(k)} }^{-1} \left( \mathbf{W}^{(k)} \boldsymbol{\nu} + \mathbf{e}^{(k)} – \Delta \mathbf{z}^{(K)} \right) \tag{30}\] \[ \Delta \mathbf{u}^{(k)} = – \mathbf{O}^{(k)} \Big( \mathbf{B}^{(k)\mathrm{T}} \big( \mathbf{S}^{(k+1)} \mathbf{A}^{(k)} \Delta \mathbf{z}^{(k)} + \mathbf{U}^{(k+1)} \boldsymbol{\nu} + \mathbf{d}^{(k+1)} + \mathbf{S}^{(k+1)} \mathbf{c}^{(k)} \big) – \mathbf{R}^{(k)} \mathbf{u}_m^{(k)} \Big) \tag{31} \] ただし、 \[ \Delta \mathbf{z}^{(0)} = \mathbf{z}_0 – \bar{\mathbf{z}}^{(0)} \tag{32} \]

\[ \Delta \mathbf{z}^{(K)} = \mathbf{z}_\mathrm{f} – \bar{\mathbf{z}}^{(K)} \tag{33} \]

手順3）以下の式（式34、35）で区間 \(t=[0,t_\mathrm{f}]\) の \(\bar{\mathbf{z}}\) および \(\bar{\mathbf{u}}\) を \[ \bar{\mathbf{z}} = \eta \Delta \mathbf{z} + \bar{\mathbf{z}} \tag{34} \] \[ \bar{\mathbf{u}} = \eta \Delta \mathbf{u} + \bar{\mathbf{u}} \tag{35} \] と更新する。

手順4）手順1）に戻って収束するまで繰り返す。

練習問題5-2

式 (22)–(32) を導出しなさい。

練習問題5-3

時間間隔 \(\Delta t \to 0\) の極限で式 22、23、24、25、26、27 が連続時間の場合のそれと一致することを確認しなさい。

\(\mathbf{Q}^{(k)\mathrm{T}} = \mathbf{Q}^{(k)}\) および \(\mathbf{R}^{(k)\mathrm{T}} = \mathbf{R}^{(k)}\) が成り立つとき、\(\mathbf{S}^{(k)\mathrm{T}} = \mathbf{S}^{(k)}\)、\(\mathbf{U}^{(k)\mathrm{T}} = \mathbf{V}^{(k)}\)、\(\mathbf{W}^{(k)\mathrm{T}} = \mathbf{W}^{(k)}\) が成り立つ。

練習問題5-4

区間 \(k=[0,K]\) で \(\mathbf{Q}^{(k)\mathrm{T}} = \mathbf{Q}^{(k)}\) および \(\mathbf{R}^{(k)\mathrm{T}} = \mathbf{R}^{(k)}\) が成り立つとき、\(\mathbf{S}^{(k)\mathrm{T}} = \mathbf{S}^{(k)}\)、\(\mathbf{U}^{(k)\mathrm{T}} = \mathbf{V}^{(k)}\)、\(\mathbf{W}^{(k)\mathrm{T}} = \mathbf{W}^{(k)}\) であることを示しなさい。

\(\mathbf{Q}^{(k)}=0\) のとき、\(\mathbf{S}^{(k)}=0\)、\(\mathbf{d}^{(k)}=0\) が成り立ちます。

練習問題5-5

区間 \(k=[0,K]\) で \(\mathbf{Q}^{(k)}=0\) のとき、\(\mathbf{S}^{(k)}=0\)、\(\mathbf{d}^{(k)}=0\) が成立することを示しなさい。

区間 \(k=[0,K]\) で \(\mathbf{Q}^{(k)}=0\) のとき、式 23、24、25、26、27、28 は簡略化され、

\[ \mathbf{U}^{(k)} = \mathbf{A}^{(k)\mathrm{T}} \mathbf{U}^{(k+1)} \tag{35} \]

\[ \mathbf{V}^{(k)} = \mathbf{V}^{(k+1)} \mathbf{A}^{(k)} \tag{36} \]

\[ \mathbf{W}^{(k)} = \mathbf{W}^{(k+1)} – \mathbf{V}^{(k+1)} \mathbf{B}^{(k)} \mathbf{R}^{(k)-1} \mathbf{B}^{(k)\mathrm{T}} \mathbf{U}^{(k+1)} \tag{37} \]

\[ \mathbf{e}^{(k)} = \mathbf{e}^{(k+1)} + \mathbf{V}^{(k+1)} (\mathbf{B}^{(k)} \mathbf{u}_m^{(k)} + \mathbf{c}^{(k)}) \tag{38} \]

となります。また、式 31 より

\[ \Delta \mathbf{u}^{(k)} = – \mathbf{R}^{(k)-1} \mathbf{B}^{(k)\mathrm{T}} \mathbf{U}^{(k+1)} \boldsymbol{\nu} + \mathbf{u}_m^{(k)} \tag{39} \] \[ \boldsymbol{\nu} = – \mathbf{W}^{(0)-1} (\mathbf{V}^{(0)} \Delta \mathbf{z}^{(0)} + \mathbf{e}^{(0)} – \Delta \mathbf{z}^{(K)}) \tag{40} \]

\[ \Delta \mathbf{z}^{(k)} = – \mathbf{V}^{(k)-1} (\mathbf{W}^{(k)} \boldsymbol{\nu} + \mathbf{e}^{(k)} – \Delta \mathbf{z}^{(K)}) \tag{41} \]

さらに、式 36 より

\[ \begin{eqnarray} \mathbf{U}^{(k)} &=& \mathbf{A}^{(k)\mathrm{T}} \mathbf{U}^{(k+1)} \\ &=& \mathbf{A}^{(k)\mathrm{T}} \mathbf{A}^{(k+1)\mathrm{T}} \mathbf{U}^{(k+2)} \\ &=& \mathbf{A}^{(k)\mathrm{T}} \mathbf{A}^{(k+1)\mathrm{T}} \dots \mathbf{A}^{(K-1)\mathrm{T}} \mathbf{U}^{(K)} \\ &=& \prod_{k_1=k}^{K-1} \mathbf{A}^{(k_1)\mathrm{T}} \tag{42} \end{eqnarray} \]

式 37 より

\[ \begin{eqnarray} \mathbf{V}^{(k)} &=& \mathbf{V}^{(k+1)} \mathbf{A}^{(k)} \\ &=& \mathbf{V}^{(k+2)} \mathbf{A}^{(k+1)} \mathbf{A}^{(k)} \\ &=& \mathbf{V}^{(K)} \mathbf{A}^{(K-1)} \dots \mathbf{A}^{(k+1)} \mathbf{A}^{(k)} \\ &=& \prod_{k_1=k}^{K-1} \mathbf{A}^{(k_1)} \tag{43} \end{eqnarray} \]

式 36 より

\[ \begin{eqnarray} \mathbf{W}^{(k)} &=& \mathbf{W}^{(k+1)} – \mathbf{V}^{(k+1)} \mathbf{B}^{(k)} \mathbf{R}^{(k)-1} \mathbf{B}^{(k)\mathrm{T}} \mathbf{U}^{(k+1)} \\ &=& \mathbf{W}^{(k+2)} – \mathbf{V}^{(k+2)} \mathbf{B}^{(k+1)} \mathbf{R}^{(k+1)-1} \mathbf{B}^{(k+1)\mathrm{T}} \mathbf{U}^{(k+2)} \\&&- \mathbf{V}^{(k+1)} \mathbf{B}^{(k)} \mathbf{R}^{(k)-1} \mathbf{B}^{(k)\mathrm{T}} \mathbf{U}^{(k+1)} \\ &=& \mathbf{W}^{(K)} – \sum_{k_1=k}^{K-1} \mathbf{V}^{(k_1+1)} \mathbf{B}^{(k_1)} \mathbf{R}^{(k_1)-1} \mathbf{B}^{(k_1)\mathrm{T}} \mathbf{U}^{(k_1+1)} \\ &=& – \sum_{k_1=k}^{K-1} \left( \prod_{k_2=k_1}^{K-1} \mathbf{A}^{(k_2)} \mathbf{B}^{(k_1)} \mathbf{R}^{(k_1)-1} \mathbf{B}^{(k_1)\mathrm{T}} \prod_{k_2=k_1}^{K-1} \mathbf{A}^{(k_2)\mathrm{T}} \right) \tag{44} \end{eqnarray} \]

もし \(\mathbf{A}^{(k)} = \mathbf{A}\)、\(\mathbf{B}^{(k)} = \mathbf{B}\)、\(\mathbf{R}^{(k)} = \mathbf{I}_n\) と書けるとき、

\[ \prod_{k_2=k_1}^{K-1} \mathbf{A}^{(k_2)} \mathbf{B}^{(k_1)} = \mathbf{A}^{K-k} \mathbf{B} \tag{45} \]

これは可制御性行列 \(\mathbf{C} = (\mathbf{B} \; \mathbf{A}\mathbf{B} \; \dots \; \mathbf{A}^{K-1} \mathbf{B})\) の要素行列です。さらに式 39 に式 41 を代入すると、

\[ \begin{eqnarray} \mathbf{e}^{(k)} &=& \mathbf{e}^{(k+1)} + \prod_{k_1=k}^{K-1} \mathbf{A}^{(k_1)} (\mathbf{B}^{(k_1)} \mathbf{u}_m^{(k_1)} + \mathbf{c}^{(k_1)}) \\ &=& \sum_{k_2=k}^{K-1} \prod_{k_1=k_2}^{K-1} \mathbf{A}^{(k_1)} (\mathbf{B}^{(k_1)} \mathbf{u}_m^{(k_1)} + \mathbf{c}^{(k_1)}) \tag{46} \end{eqnarray} \]

4. 二次精度の公式

4.1. 一次元バネ質点システム

一次元バネ質点システムを例に二次精度の時間積分公式を構成してみましょう。式2-9、2-10、2-11より、

\[ \frac{dx}{dt} = v \tag{47} \] \[ \frac{dv}{dt} = a – \frac{k_s x – k_{s2} x^2}{m} – \gamma v – \gamma_2 v^2 \tag{48} \]

\[ \frac{da}{dt} = g \tag{49} \]

と書けます。離散時間の時間ステップを \(\Delta t\) とし、式47、48、49を二次精度まで意識して差分化すると、\[ x^{(k+1)} = x^{(k)} + \Delta t \, \frac{v^{(k)} + v^{(k+1)}}{2} \tag{50}\]

\[ \begin{aligned} v^{(k+1)} = v^{(k)} + \Delta t \Bigg( &\frac{a^{(k)} + a^{(k+1)}}{2}\\ & – \frac{k_\mathrm{s}}{m} \frac{x^{(k)} + x^{(k+1)}}{2} – \frac{k_\mathrm{s2}}{m} \left( \frac{x^{(k)} + x^{(k+1)}}{2} \right)^2\\ & – \gamma \frac{v^{(k)} + v^{(k+1)}}{2} – \gamma_2 \left( \frac{v^{(k)} + v^{(k+1)}}{2} \right)^2 \Bigg) \end{aligned} \tag{52} \]

\[ a^{(k+1)} = a^{(k)} + \Delta t \, g^{(k)} \tag{53} \]

状態変数ベクトルを

\[ \mathbf{z}^{(k)} = \begin{pmatrix} x^{(k)} & v^{(k)} & a^{(k)} \end{pmatrix}^{\mathrm T}, \]

制御変数ベクトルを

\[ \mathbf{u}^{(k)} = g^{(k)} \] とすると、 \[ \tilde{\mathbf{A}}^{(k)} = \Delta t \begin{pmatrix} 0 & 1 & 0 \\ -\dfrac{k_s + k_{s2}\left(\bar{x}^{(k)} + \bar{x}^{(k+1)}\right)}{m} & -\left(\gamma + \gamma_2\left(\bar{v}^{(k)} + \bar{v}^{(k+1)}\right)\right) & 1 \\ 0 & 0 & 0 \end{pmatrix} \tag{54} \]

\[ \tilde{B}^{(k)}=\Delta t \begin{pmatrix} 0\\ 0\\ 1 \end{pmatrix} \tag{55} \]

\[ \tilde{c}^{(k)}= \begin{pmatrix} \bar{x}^{(k)}-\bar{x}^{(k+1)}+\dfrac{\Delta t}{2}\left(\bar{v}^{(k)}+\bar{v}^{(k+1)}\right)\\[6pt] \bar{v}^{(k)}-\bar{v}^{(k+1)}+\Delta t\Biggl( \dfrac{\bar{a}^{(k)}+\bar{a}^{(k+1)}}{2}\\ -\dfrac{k_\mathrm{s}}{m}\dfrac{\bar{x}^{(k)}+\bar{x}^{(k+1)}}{2} -\dfrac{k_\mathrm{s2}}{m}\left(\dfrac{\bar{x}^{(k)}+\bar{x}^{(k+1)}}{2}\right)^2\\ \qquad\qquad -\gamma\dfrac{\bar{v}^{(k)}+\bar{v}^{(k+1)}}{2} -\gamma_2\left(\dfrac{\bar{v}^{(k)}+\bar{v}^{(k+1)}}{2}\right)^2 \Biggr)\\[6pt] \bar{a}^{(k)}-\bar{a}^{(k+1)}+\Delta t\,\bar{u}^{(k)} \end{pmatrix} \tag{56} \]

を得ます。評価関数として

\[ J=\frac{\Delta t}{2}\sum_{k=0}^{K}\left( {\mathbf{z}^{(k)}}^{\mathrm{T}}\mathbf{Q}\mathbf{z}^{(k)} +{\mathbf{u}^{(k)}}^{\mathrm{T}}\mathbf{R}\mathbf{u}^{(k)} \right) \tag{57} \]

を採用し、

\[ \mathbf{Q}=0,\; \mathbf{R}=1 \tag{58} \]

とします。

実際の問題を解いてみましょう。1次元バネ質点系で初速度を与えて制御しないと、単振動を示します。\(k_\mathrm{s}=1\、\mathrm{kg\,s^{-2}}\)、\(m=1\,\mathrm{kg}\) の場合、その周期は \(2\pi\) 秒です。\(t_\mathrm{f}=2\pi\) として、初期条件（\(t=0\)）を\(\begin{pmatrix} 1\,\mathrm{m} & 1\,\mathrm{m\,s^{-1}} & 0\,\mathrm{m\,s^{-2}} \end{pmatrix} \)、終端条件（\(t=t_\mathrm{f}\)）を\( \begin{pmatrix} 1\,\mathrm{m} & 1\,\mathrm{m\,s^{-1}} & 0\,\mathrm{m\,s^{-2}} \end{pmatrix} \)として、求めた最適解を Figure 1 に示してあります。 制御なしの自然な単振動に従って初期条件から終端条件に至るので、\(g\) の値が非常に小さくなっています。ここでは、有限の時間ステップ（\(\Delta t=0.052\,\mathrm{s}\)）による誤差のために厳密には \(g=0\) ではありませんが、より小さな時間ステップを用いると \(g\) の値が小さくなります。Figure 2 は、\(\Delta t\) に対して各時間ステップの\(g\) の二乗平均の平方根（\(\mathrm{RMS}(g)\)）をプロットしました。\(\mathrm{RMS}(g)\) は、\(\Delta t\) の自乗に反比例して小さくなります。これは、本節で導いた離散時間状態空間モデルが、実際に時間に対して二次精度を実現していることを示しています。 時間二次精度の差分公式を使った場合、差分誤差は\(O(\Delta t^3)\) から現れますが、初期状態から終端条件のすべての時間ステップの量が関与する\(\mathrm{RMS}(g)\) のような量の誤差は\(O(\Delta t^2)\) から現れます。

次に、Figure 3 にバネなしの場合の二次精度の公式（式50–55）を用いて求めた到達問題の躍度最小解を示しています。\(K=10\) でも高精度の解が得られています。運動の脳の内部モデルにおいても、時間に対して二次精度の公式を適用し、計算時間と計算資源の節約を図っていると考えられます。同様に、バネや抵抗がある場合も、少数の時間ステップで躍度最小解を求めることができます。

4.2. 二次元質点系

バネと磁場を考慮に入れたx-y平面上を運動する質点系について議論します。\(\mathbf{z} = \begin{pmatrix} x & v_x & a_x & y & v_y & a_y \end{pmatrix}^{\mathrm{T}}, \quad \mathbf{u} = \begin{pmatrix} g_x & g_y \end{pmatrix}^{\mathrm{T}}\)と置けば、その状態方程式は 二次元質点系の状態方程式は次のように書けます： は \(x–y\) 平面を垂直に貫く外部磁場の強さです。 二次元質点系の状態方程式は次のように書けます：

\[ \frac{d\mathbf{z}}{dt} = \begin{pmatrix}0 & 1 & 0 & 0 & 0 & 0 \\ -\frac{k_\mathrm{s}}{m} & 0 & 1 & 0 & \frac{q_\mathrm{c} b_\mathrm{z}}{m} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & -\frac{q_\mathrm{c} b_\mathrm{z}}{m} & 0 & -\frac{k_s}{m} & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0\end{pmatrix} \mathbf{z} + \begin{pmatrix}0 & 0 \\ 0 & 0 \\ 1 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 1\end{pmatrix} \mathbf{u} \tag{59}\]

ここで、\(\mathit{m}\) は質点の質量、\(\mathit{k_{\mathrm{s}}}\) はバネ定数、\(\mathit{q_{\mathrm{c}}}\)は質点の電荷、\(\mathit{b_{\mathrm{z}}}\) は \(x–y\) 平面を垂直に貫く外部磁場の強さです。初期状態を\(\mathbf{z}_0 = \begin{pmatrix}x_0 & v_{x0} & a_{x0} & y_0 & v_{y0} & a_{y0}\end{pmatrix}^{\mathrm{T}}\)、終端状態を\(\mathbf{z}_\mathrm{f} = \begin{pmatrix}x_\mathrm{f} & v_{x\mathrm{f}} & a_{x\mathrm{f}} & y_\mathrm{f} & v_{y\mathrm{f}} & a_{y\mathrm{f}}\end{pmatrix}^{\mathrm{T}}\)と置きます。 ラグランジュ関数 \(L(\mathbf{z}(t),\mathbf{u}(t),t)\) が \(\mathbf{z}\) と \(\mathbf{u}\) の二次形式で表されるとします。つまり、

\[L(\mathbf{z}(t),\mathbf{u}(t),t) = \mathbf{z}(t)^\mathrm{T} \mathbf{Q} \mathbf{z}(t) + \mathbf{u}(t)^\mathrm{T} \mathbf{R} \mathbf{u}(t) \tag{60}\] です。 \(\displaystyle \mathbf{Q} = 0, \quad \mathbf{R} = I_2 \tag{61}\) とすれば、式10は \[L = \frac{1}{2} \left( g_x^2 + g_y^2 \right) \tag{62}\]

となります。 時間ステップを \(\Delta t\) として式59を時間二次精度で離散化すると、

\[\mathbf{z}^{(k+1)} = \mathbf{A}^{(k)} \mathbf{z}^{(k)} + \mathbf{B}^{(k)} \mathbf{u}^{(k)} + \mathbf{c}^{(k)} \tag{63}\]

を得ます。ここで、

\[\mathbf{A}^{(k)} = \left(I_6 – \frac{1}{2} \tilde{\mathbf{A}} \right)^{-1} \left(I_6 + \frac{1}{2} \tilde{\mathbf{A}} \right) \tag{64}\] \[\mathbf{B}^{(k)} = \left(I_6 – \frac{1}{2} \tilde{\mathbf{A}} \right)^{-1} \tilde{\mathbf{B}} \tag{65}\]

\[\mathbf{c}^{(k)} = \left(I_6 – \frac{1}{2} \tilde{\mathbf{A}} \right)^{-1} \tilde{\mathbf{c}}^{(k)} \tag{66}\] \[\tilde{\mathbf{A}} = \Delta t \mathbf{A} \tag{67}\] \[\tilde{\mathbf{B}} = \Delta t \mathbf{B} \tag{68}\]

\[\tilde{\mathbf{c}}^{(k)} = \begin{pmatrix} \bar{x}^{(k)} – \bar{x}^{(k+1)} + \Delta t \frac{\bar{v}_x^{(k)} + \bar{v}_x^{(k+1)}}{2} \\ \bar{v}_x^{(k)} – \bar{v}_x^{(k+1)} + \Delta t \frac{\bar{a}_x^{(k)} + \bar{a}_x^{(k+1)}}{2} – \frac{k_\mathrm{s}}{m} \frac{\bar{x}^{(k)} + \bar{x}^{(k+1)}}{2} + \frac{q_\mathrm{c} b_\mathrm{z}}{m} \frac{\bar{v}_y^{(k)} + \bar{v}_y^{(k+1)}}{2} \\ \bar{a}_x^{(k)} – \bar{a}_x^{(k+1)} + \Delta t \bar{u}_x^{(k)} \\ \bar{y}^{(k)} – \bar{y}^{(k+1)} + \Delta t \frac{\bar{v}_y^{(k)} + \bar{v}_y^{(k+1)}}{2} \\ \bar{v}_y^{(k)} – \bar{v}_y^{(k+1)} + \Delta t \frac{\bar{a}_y^{(k)} + \bar{a}_y^{(k+1)}}{2} – \frac{k_s}{m} \frac{\bar{y}^{(k)} + \bar{y}^{(k+1)}}{2} – \frac{q_\mathrm{c} b_\mathrm{z}}{m} \frac{\bar{v}_x^{(k)} + \bar{v}_x^{(k+1)}}{2} \\ \bar{a}_y^{(k)} – \bar{a}_y^{(k+1)} + \Delta t \bar{u}_y^{(k)} \end{pmatrix} \tag{69}\]

です。では、式63、64、65、66、67、68、69を使って実際に躍度最小解を求めてみましょう。

4.2.1. バネも磁場もないとき

状態ベクトルの1–3（\(x\)）成分と4–6（\(y\)）成分が相互作用しないので、それぞれの成分を1次元質点系として独立に解くことができます。したがって、任意の二点の間の到達問題は、始点と終点を結んだ直線上の一次元問題に帰着します。Figure 4に、\(\mathbf{z}_0 =\begin{pmatrix} -1/\sqrt{2} & 0 & 0 & -1/\sqrt{2} & 0 & 0 \end{pmatrix}^{\mathrm{T}}\)、\(\mathbf{z}_f = \begin{pmatrix}1/\sqrt{2} & 0 & 0 & 1/\sqrt{2} & 0 & 0 \end{pmatrix}^{\mathrm{T}}\)としたときと、それらを\(x\)方向に \(±0.2\,\mathrm{m}\) だけ平行移動した場合の躍度最小解を表しています。\(x–y\) 平面上では、どの場合の軌跡も直線になります。その直線に沿った長さを \(s\)、速度を \(v_s\)、加速度を \(a_s\)、躍度を \(g_s\) とした到達問題の躍度最小解は、一次元問題のそれと同一です。

4.2.2. バネ力による運動

まず、制御なしで実現する運動を再現してみましょう。磁場がないとき、バネ力と遠心力が釣り合う円運動を考えましょう。

\[ \frac{mv^2}{r} = k_\mathrm{s} r \tag{70} \]

\(r\)について解くと円の半径 \(R_\mathrm{c}\) は、

\[ r = R_\mathrm{c} = \sqrt{\frac{mv^2}{k_\mathrm{s}}} \tag{71} \]

を得ます。\(k_\mathrm{s} = 1\ \mathrm{kg\,s^{-2}},\ m = 1\ \mathrm{kg},\ v = 1\ \mathrm{m}\,\mathrm{s}^{-1}\) のとき、\(R_\mathrm{c} = 1\ \mathrm{m}\)、周期 \(T_\mathrm{c} = 2\pi\ \mathrm{s}\) となります。Figure 4 に実際の躍度最小解を示しています。\(K=12\) でも \(K=120\) の時の解に十分近い解を得ています。このとき、\(\mathrm{RMS}(g) = 0.1625\) と時間積分の誤差の範囲で十分小さくなっています。次に、到達問題を解いてみましょう。Figure 4 に与えた始点と終点に対して、バネがある場合の躍度最小解を求めてみましょう。バネの固定点（\(x=0, y=0\)）を通る軌道は直線になりますが、それに対して横ずれさせると、横方向のバネの力のせいで軌道が直線からずれます（Figure 6）。

 

4.2.3. 磁場による運動

今度はバネなしで磁場がある場合を考えましょう。一様磁場中で荷電粒子は、磁場からのローレンツ力を受けてサイクロトロン運動と呼ばれる円運動を行います。サイクロトロン運動では、ローレンツ力と遠心力が釣り合います。つまり、

\[\frac{mv^2}{r} = q_\mathrm{c} b_\mathrm{z} v \tag{72}\]

これを \(r\) について解くと円の半径 \(R_\mathrm{c}\) は、

\[ r = R_\mathrm{c} = \frac{mv}{q_\mathrm{c} b_\mathrm{z}} \tag{73} \]

を得ます。また、円運動の周期 \(T_\mathrm{c}\) は、

\[ T_\mathrm{c} = \frac{2\pi R_\mathrm{c}}{v} = \frac{2\pi m}{q_\mathrm{c} b_\mathrm{z}} \tag{74}\]

です。\(q_\mathrm{c} = 1\ \mathrm{C},\ b_z = 1\ \mathrm{T},\ m = 1\ \mathrm{kg},\ v = 1\ \mathrm{m/s}\) の時、\(R_\mathrm{c} = 1\ \mathrm{m},\ T_\mathrm{c} = 2\pi\ \mathrm{s}\) となります。Figure 7 に実際の躍度最小解を示しています。\(K=12\) でも \(K=120\) の時の解に十分近い解を得ています。このとき、\(\mathrm{RMS}(g) = 0.1625\) と時間積分の誤差の範囲で十分小さくなっています。次に、到達問題を解いてみましょう。Figure 4 に与えた始点と終点に対して、磁場がある場合の躍度最小解を求めてみましょう。終点に向かう速度に応じて横方向のローレンツ力が働くので軌道が直線からずれます（Figure 8）。

 

5. まとめ

本章では、離散時間状態モデルにおける最適制御アルゴリズム（天馬方程式）を説明しました。それを使って、時間二次精度のモデルを構築し、時間ステップが10程度でも高い精度（10％程度）で最適制御解を得ることができることを示しました。また、誤差は、\(\Delta t^2\) に比例して減少することを示しました。この公式を用いて一次元と二次元の質点系のバネ力とローレンツ力（磁場による力）による運動が生成できることを示しました。さて、これまで私たちは躍度最小解を求めるにあたって、それを表現する状態空間モデルを与える二つの行列 \(\mathbf{A}\) と \(\mathbf{B}\) は既知であるとしてきました。では、脳は自分が制御する身体の運動パラメーター（行列 \(\mathbf{A}\) と \(\mathbf{B}\)）を知るのでしょうか？制御工学において、制御したい系のパラメーター決める作業はシステム同定と呼ばれています。次章（７章）では脳のシステム同定について考察します。

引用文献