非線形システムの局所線形化と躍度最小解 ―脳の設計図を求めて（3）―

1. 非線形システム

第二章において、力学系を状態空間モデルで記述できることを示しました。もし、それが線形システムであれば、つまり、状態方程式が状態変数に依存しない係数行列AとBを用いた状態変数の線形方程式で記述でき、評価関数が状態変数の二次形式の時間積分の形で書ければ、最適制御解法により自動的に評価関数を最小（大）化する軌道（状態変数の時間プロファイル）とそれを実現するための制御変数の時間プロファイルを得る手法を示しました。

残念ながら多くの力学系は線形方程式では記述できないし、評価関数も状態変数や制御変数の二次形式の形には書けません。しかし、以下に述べる局所線形化の技法により、非線形な力学系の最適制御解を、たくさんの線形な力学系の解をつなぎ合わせることで任意の精度で求めることができることを示します。まず、2節では、局所線形化の技法とそれを用いた最適制御解法を説明します。次に3節は、この技法を簡単な非線形力学系に適用し、実際に解が得られることを示します。

2. 局所線形化と最適制御解法

状態方程式の右辺\({\bf f}({\bf z}(t),{\bf u}(t),t)\)は多くの物理系では\({\bf z}(t)\)と\({\bf u}(t)\)に依存するが、\(t\)には直接には依存せず、連続で適切な回数微分可能とします。つまり、\({\bf f}({\bf z}(t),{\bf u}(t))\)と書けます。このような場合、ある点の周りの微小区間において局所線形化が可能です。つまり、\({\bf z}=\bar{\bf z} +{∆{\bf z}}\)、\({\bf u}=\bar{\bf u }+∆{\bf u}\)として\({\bf f}({\bf z}(t),{\bf u}(t),t)\)を、\(∆{\bf z}\)と\(∆{\bf u}\)について展開してそれらの一次の項までとると

\begin{equation} \dot{\bar{\bf z} }+∆\dot{\bf z}={\bf f}(\bar{\bf z},\bar{\bf u})+\left.\frac{∂\bf{f}}{∂{\bf z}^{\rm T}}\right|_{\bar{\bf z},\bar{\bf u}} {∆{\bf z}}+\left.\frac{∂{\bf f}}{∂{\bf u}^{\rm T}}\right|_{\bar{\bf z},\bar{\bf u}} {∆{\bf u}} \tag{1}\label{eq1} \end{equation} となります。ここで、\({\bf A}(\bar{\bf z},\bar{\bf u})=\left.\frac{∂{\bf f}}{∂{\bf z}^{\rm }}\right|_{\bar{\bf z},\bar{\bf u}}\)、\({\bf B}(\bar{\bf z},\bar{\bf u})=\left. \frac{∂{\bf f}}{∂{\bf u}^{\rm T}}\right|_{\bar{\bf z},\bar{\bf u}}\)、\({\bf c}(\bar{\bf z},\bar{\bf u})=f(\dot{\bar{\bf z}},\bar{\bf u})-\bar{\bf z} \) と置くと、 \begin{equation} {∆\dot{\bf z}}={\bf A}(\bar{z},\bar{\bf u}){∆{\bf z}}+{\bf B}(\bar{\bf z},\bar{\bf u}){∆{\bf u}}+{\bf c}(\bar{\bf z},\bar{\bf u}) \tag{2} \end{equation} と\({∆{\bf z}}\)と\({\bf u}(t)\)について線形な状態方程式で書き表せます。ここで、\({\bf A}(\bar{\bf z},\bar{\bf u})\)、\({\bf B}(\bar{\bf z},\bar{\bf u})\)、\({\bf c}(\bar{\bf z},\bar{\bf u})\)は、\(\bar{\bf z}(t)\)と\(\bar{\bf u}(t)\)の関数であるが、\(t\)には直接は関係しないとします。

また、評価関数\(L({\bf z}(t),{\bf u}(t))\)は一般には\({\bf z}\)と\({\bf u}\)の二次形式にはなりません。しかし、\({\bf z}=\bar{\bf z}+{∆{\bf z}}\)、\({\bf u}=\bar{\bf u}+{∆{\bf u}}\)として\({∆{\bf z}}\)と\({∆{\bf u}}\)について展開して二次の項までとると、 \begin{eqnarray} L({\bf z},{\bf u})&=&L(\bar{\bf z},\bar{\bf u})+{\left.\frac{∂L}{∂{\bf z}^{\rm T}}\right|_{\bar{\bf z},\bar{\bf u}}}{∆{\bf z}}+ \left. \frac{∂L}{∂{\bf u}^{\rm T}}\right|_{\bar{\bf z},\bar{\bf u}} {∆{\bf u}}\\ &&+\frac{1}{2} {∆{\bf z}}^{\rm T}\left.\frac{∂^2 L}{∂{\bf z}^2}\right|_{\bar{\bf z},\bar{\bf u}} {∆{\bf z}} +\frac{1}{2}{∆{\bf u}}^{\rm T}\left.\frac{∂^2 L}{∂{\bf u}^2}\right|_{\bar{\bf z},\bar{\bf u}}{∆{\bf u}}\\ &=&L+{\bf p}^{\rm T}{∆{\bf z}}+ {\bf q}^{\rm T}{∆{\bf u}}+\frac{1}{2}({∆{\bf z}}^{\rm T} {\bf Q}{∆{\bf z}}+{∆{\bf u}}^{\rm T} {\bf R}{∆{\bf u}})\\ &=&L+\frac{1}{2} ({∆{\bf z}}^{\rm T} {\bf Q}{∆{\bf z}}+2{\bf p}^{\rm T}{∆{\bf z}})+\frac{1}{2}({∆{\bf u}}^{\rm T} {\bf R}{∆{\bf u}}+2{\bf q}^{\rm T}{∆{\bf u}})\\ &=&L+\frac{1}{2} ({∆{\bf z}}+{\bf Q}^{-1} {\bf p})^{\rm T} {\bf Q}({∆{\bf z}}+{\bf Q}^{-1} {\bf p})-\frac{1}{2} {\bf p}^{\rm T} {\bf Q}^{-1} {\bf p}\\ &&+\frac{1}{2}({∆{\bf u}}+{\bf R}^{-1} {\bf q})^{\rm T} {\bf R}({∆{\bf u}}+{\bf R}^{-1} {\bf q})-\frac{1}{2}{\bf q}^{\rm T} {\bf R}^{-1} {\bf q}\\ &=& \frac{1}{2} ({∆{\bf z}}+{\bf Q}^{-1} {\bf p})^{\rm T} {\bf Q}({∆{\bf z}}+{\bf Q}^{-1} {\bf p})\\ &&+\frac{1}{2}({∆{\bf u}}+{\bf R}^{-1} {\bf q})^{\rm T} {\bf R}({∆{\bf u}}+{\bf R}^{-1} {\bf q})\\ &&+L-\frac{1}{2}{\bf p}^{\rm T} {\bf Q}^{-1} {\bf p}-\frac{1}{2} {\bf q}^{\rm T} {\bf R}^{-1} {\bf q} \tag{3}\label{eq3} \end{eqnarray}

と\({∆{\bf z}}\)と\({∆{\bf u}}\)についての二次形式に変形できます。ここで、\({\bf Q}(\bar{\bf z},\bar{\bf u})=\left.\frac{∂^2 L}{∂{\bf z}^2}\right|_{\bar{\bf z},\bar{\bf u}}\)、\({\bf R}(\bar{\bf z},\bar{\bf u})=\left.\frac{∂^2 L}{∂{\bf u}^2}\right|_{\bar{\bf z},\bar{\bf u}}\)、\({\bf p}(\bar{\bf z},\bar{\bf u})=\left.\frac{∂L}{∂{\bf z}^{\rm T}}\right|_{\bar{\bf z},\bar{\bf u}}\)、\({\bf q}(\bar{\bf z},\bar{\bf u})=\left. \frac{∂L}{∂{\bf u}^{\rm T}}\right|_{\bar{\bf z},\bar{\bf u}}\)です。ただし、簡単のため\(\left.\frac{∂^2 L}{∂{\bf z}∂{\bf u}}\right|_{{\bar{\bf z},\bar{\bf u}}}=0\)を仮定しました。また、\({\bf Q}\)と\({\bf R}\)は対称行列であるとします。したがって、\({{\bf Q}^{-1}}^{\rm T}={\bf Q}^{-1}\)、\({{\bf R}^{-1}}^{\rm T}={\bf R}^{-1}\)が成立します。ここで、 \begin{equation} {\bf z}_{\rm m}=-{\bf Q}^{-1} {\bf p} \tag{4} \end{equation} \begin{equation} {\bf u}_{\rm m}=-{\bf R}^{-1} {\bf q} \tag{5} \end{equation} とおけば、

\begin{eqnarray} L({∆{\bf z}},{∆{\bf u}})&=&\frac{1}{2}({∆{\bf z}}-{\bf z_{\rm m}})^{\rm T}{\bf Q}({∆{\bf z}}-{\bf z_{\rm m}})+\frac{1}{2}({∆{\bf u}}-{\bf u_{\rm m}})^{\rm T}{\bf R}({∆{\bf u}}-{\bf u_{\rm m}})\\ &&+L_0(\bar{\bf z},\bar{\bf u}) \tag{6}\label{eq6} \end{eqnarray} と変形できます。ここで、\({∆{\bf z}}\)にも\({∆{\bf u}}\)にも依存しない項\(L_0(\bar{\bf z},\bar{\bf u})\)は \begin{eqnarray} L_0(\bar{\bf z},\bar{\bf u})&=&-\frac{1}{2} {\bf p}(\bar{\bf z},\bar{\bf u})^{\rm T} {\bf Q}(\bar{\bf z},\bar{\bf u})^{-1} {\bf p}(\bar{\bf z},\bar{\bf u})-\frac{1}{2} {\bf q}(\bar{\bf z},\bar{\bf u})^{\rm T} {\bf R}(\bar{\bf z},\bar{\bf u})^{-1} {\bf q}(\bar{\bf z},\bar{\bf u})\\ &&+L(\bar{\bf z},\bar{\bf u}) \tag{7} \end{eqnarray} です。このときオイラー・ラグランジュ方程式は、 \begin{equation} \dot{\bf λ}=-{\bf Q}({∆{\bf z}}-{\bf z}_{\rm m})-{\bf A}^{\rm T} {\bf λ} \tag{8} \end{equation} \begin{equation} {\bf R}({∆{\bf u}}-{\bf u}_{\rm m})+{\bf B}^{\rm T} {\bf λ}={\bf 0} \tag{9} \end{equation} となります。

この場合、初期条件\({\bf z}(0)={\bf z}_0\)、終端条件\({\bf z}(t_{\rm f})={\bf z}_{\rm f}\)を満たす到達問題の解\({\bf z}(t)\)とそれを実現する制御変数の解\({\bf u}(t)\)は、以下の繰り返し計算で求められます。

手順0) 時間の連続関数\(\bar{\bf z}(t)\)を仮に与えます。行列\({\bf A}\)、\({\bf B}\)、\({\bf Q}\)、\({\bf R}\)は、\(\bar{\bf z}(t)\)の関数として既に与えられているとします。

手順1) 以下の6つの微分方程式（式\ref{eq10}、\ref{eq11}、\ref{eq12}、\ref{eq13}、\ref{eq14}、\ref{eq15}）を、\({\bf S}(t_{\rm f})={\bf 0}\)、\({\bf U}(t_{\rm f})={\bf I}_{\rm n}\)、 \({\bf V}(t_{\bf f})={\bf I}_{\rm n}\)、\({\bf W}(t_{\rm f})={\bf 0}\)、\({\bf d}(t_{\rm f})={\bf 0}\)、\({\bf e}(t_{\rm f})={\bf 0}\)の終端条件から、時間のマイナス方向に積分して、\(t=[0,t_{\rm f}]\)の区間において4つの行列\({\bf S}\)、\({\bf U}\)、\({\bf V}\)、\({\bf W}\)および2つのベクトル\({\bf d}\)、\({\bf e}\)を求めます。

\begin{equation} \dot {\bf S}=-{\bf A}^{\rm T} {\bf S}-{\bf SA}+{\bf SB}{\bf R}^{-1} {\bf B}^{\rm T} {\bf S}-{\bf Q} \tag{10}\label{eq10} \end{equation} \begin{equation} \dot{\bf U} =-{\bf A}^{\rm T} {\bf U}+{\bf SB}{\bf R}^{-1} {\bf B}^{\rm T} {\bf U} \tag{11}\label{eq11} \end{equation} \begin{equation} \dot{\bf V}=-{\bf VA}+{\bf VB}{\bf R}^{-1} {\bf B}^{\rm T} {\bf S} \tag{12}\label{eq12} \end{equation} \begin{equation} \dot{\bf W} ={\bf VB}{\bf R}^{-1} {\bf B}^{\rm T} {\bf U} \tag{13} \label{eq13}\end{equation} \begin{equation} \dot{\bf d} =-{\bf A}^{\rm T} {\bf d}+{\bf SB}{\bf R}^{-1} {\bf B}^{\rm T} {\bf d}-{\bf Sc}-{\bf SB}{\bf u}_{\rm m}+{\bf Q}{\bf z}_{\rm m} \tag{14}\label{eq14} \end{equation} \begin{equation} \dot{\bf e}={\bf VB}{\bf R}^{-1} {\bf B}^{\rm T} {\bf d}-{\bf VB}{\bf u}_{\rm m}-{\bf Vc} \tag{15}\label{eq15} \end{equation}

手順2) 以下の式（式\ref{eq16}、\ref{eq17}、\ref{eq18}）で\(ν\)を求めます。 \begin{equation} {∆{\bf z}}_0={\bf z}_0-\bar{\bf z}(0) \tag{16}\label{eq16} \end{equation} \begin{equation} {∆{\bf z}}_{\rm f}={\bf z}_{\rm f}-\bar{\bf z}(t_{\rm f} ) \tag{17} \label{eq17}\end{equation} \begin{equation}{\bf ν}=-{\bf W}(0)^{-1} ({\bf V}(0){∆{\bf z}}_0+{\bf e}(0)-{∆{\bf z}}_{\rm f}) \tag{18}\label{eq18}\end{equation}

手順3) 以下の式（式\ref{eq19}、\ref{eq20}、\ref{eq21}、\ref{eq22}）で区間\(t=[0,t_{\rm f}]\)の\(\bar{\bf z}\)および\(\bar{\bf u}\)を

\begin{equation} {∆{\bf z}}=-{\bf V}^{-1}({\bf Wν}+{\bf e}-{∆{\bf z}}_{\rm f }) \tag{19}\label{eq19} \end{equation} \begin{equation} {∆{\bf u}}=-{\bf R}^{-1} {\bf B}^{\rm T} ({\bf S}{∆{\bf z}}+{\bf Uν}+{\bf d})+{\bf u}_{\rm m} \tag{20}\label{eq20} \end{equation} \begin{equation} \bar{\bf z}=η{∆{\bf z}}+\bar{\bf z} \tag{21}\label{eq21} \end{equation} \begin{equation} \bar{\bf u}=η{∆{\bf u}}+\bar{\bf u} \tag{22} \label{eq22}\end{equation} と更新します。

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

練習問題3.1

式\ref{eq3}を状態方程式、式\ref{eq6} を評価関数とする到達問題の最適解が上記の手順1－3式（式\ref{eq10}‐\ref{eq20}）で得られることを示しなさい。

線形の場合と同様に、\(t=[0,t_{\bf f}]\)の全ての時間で\({\bf Q}^{\rm T}={\bf Q}\)と\({\bf R}^{\rm T}={\bf R}\)が成り立つとき、\({\bf S}^{\rm T}={\bf S}\)、\({\bf U}^{\rm T}={\bf V}\)、\({\bf W}^{\rm T}={\bf W}\)が成り立ちます。さらに、\(t=[0,t_{\rm f}]\)の全ての時間で\({\bf Q}={\bf 0}\)が成り立つとき、\({\bf S}={\bf 0}\)、\({\bf d}={\bf 0}\)も成立します。

練習問題3.2

\(t=[0,t_{\bf f}]\)の全ての時間で\({\bf Q}={\bf 0}\)が成り立つとき、\({\bf S}＝{\bf 0}\)と\({\bf d}={\bf 0}\)が成立することを示しなさい。

局所線形化された最適制御解法（式\ref{eq10}‐\ref{eq22}）を用いれば、自動的に到達問題の躍度最小解を任意の精度で得られて大変便利です。私は、この一群の式を「天馬方程式」と呼ぶことにしました。手塚治虫による漫画「鉄腕アトム」の主人公アトムは、天才科学者「天馬博士」が作ったロボットでした。その運動制御部分にはこの「天馬方程式」が実装されていたはずです。

式\ref{eq11}、\ref{eq12}、\ref{eq13}、\ref{eq15}に\({\bf Q}={\bf 0}\)、\({\bf S}＝{\bf 0}\)、\({\bf d}={\bf 0}\)を代入すると \begin{equation} \dot{\bf U} =-{\bf A}^{\rm T} {\bf U} \tag{23} \label{eq23}\end{equation} \begin{equation} \dot{\bf V}=-{\bf VA} \tag{24} \label{eq24}\end{equation} \begin{equation} \dot{\bf e}=-{\bf VB}{\bf u}_{\rm m}-{\bf Vc}=-{\bf V}({\bf B}{\bf u}_{\rm m}+{\bf c}) \tag{25}\label{eq25} \end{equation} を得ます。このように、\(t=[0,t_{\rm f}]\)の全ての時間で\({\bf Q}={\bf 0}\)が成り立つ場合は、天馬方程式が簡略化されて、式\ref{eq10}‐\ref{eq15}の6つの式の代わりに、式\ref{eq23}-\ref{eq25}と式\ref{eq13}の4つの式のみを使えばよいことが分かります。これらは非線形項が式\ref{eq13}にしか現れず、ずっと簡単に躍度最小解を求めることができます。

二章に記述したように、線形の場合には、\(η=1\)で一回で厳密解が得られます。 非線形（\({\bf A}\)と\({\bf B}\)が\({\bf z}\)や\({\bf u}\)の関数）の場合でも、10回以下の繰り返しで十分な収束を得ます。非線形性が強くて収束しないときには、\(η\)に小さな値、例えば0.1程度の値を使えばよいのですが、十分な精度を得るための繰り返し回数は数十回に増えます。

上記手続き1-4においては、最適化は線形の範囲では厳密に実行されていることに注意しましょう。繰り返し計算は係数行列\({\bf A}(\bar{\bf z},\bar{\bf u})\)、\({\bf B}(\bar{\bf z},\bar{\bf u})\)、\({\bf Q}(\bar{\bf z}\),\(\bar{\bf u})\)、\({\bf R}(\bar{\bf z},\bar{\bf u})\)を計算するときに仮定した\(\bar{\bf z}\)や\(\bar{\bf u}\)と、最適化の結果求まった\(\bar{\bf z}\)や\(\bar{\bf u}\) の値との不整合を解消するために行われています。それは、ニュートン法と同様に一様収束します。つまり、真の解に十分近い出発解を用いれば必ず収束します。到達問題においては既知の、始点と終点を結ぶ直線、もしくは式2-29で与えられるFlash-Hoganの解析解が良い出発解となると期待できます。

　次節から、簡単な非線形力学系における到達問題の解が、上記1－3の手順（式\ref{eq10}‐\ref{eq22}）で実際に得られることを簡単な例で示します。まず、3節では、一次元質点系で、バネと抵抗による力が非線形成分を持つ場合を取り上げます。このような非線形な力学系でも局所線形化した最適制御解法（天馬方程式：式\ref{eq10}‐\ref{eq22}）は、急速に収束し数十回の繰り返しの後十分な精度の解を与えます。

また4節では、軌道拘束された二次元質点系を取り上げます。このような系は軌道拘束のために遠心力やコリオリ力に起因する非線形な力が現れます。この場合にも局所線形化された最適制御解法は、高精度の解を与えます。それは、拘束曲線が対数螺旋の場合は解析解を、それが正弦対数曲率曲線の場合は近似解析解と良く一致していました。これらの解析解や近似解析解は、心理物理実験の結果と良く整合します。

3. 一次元質点系：バネと抵抗による力が非線形成分を持つ場合

一次元バネ質点系で、バネと抵抗力が非線形成分を持つ場合を考えます。状態（運動）方程式が

\begin{equation} \frac{dx}{dt}=v \tag{26}\label{eq26} \end{equation} \begin{equation} \frac{dv}{dt}=-\frac{k_{\rm s} x+k_{\rm s2} x^2}{m}-(γv+γ_2 v^2)+a \tag{27} \label{eq27}\end{equation} \begin{equation} \frac{da}{dt}=u \tag{28}\label{eq28} \end{equation}

と表せるとしましょう。ここでは\(k_{\rm s2} x^2\)と\(γ_2 v^2\)の項がそれぞれバネと抵抗力の非線形成分を表しています。式\ref{eq26}‐\ref{eq28}を局所線形化すると \begin{equation} \frac{dΔx}{dt}=Δv+\bar{v}-\frac{d\bar{x}}{dt} \tag{29} \end{equation} \begin{eqnarray} \frac{dΔv}{dt}&=&-\frac{(k_{\rm s}+2k_{\rm s2} \bar{x})}{m} Δx-(γ+2γ_2\bar{v})Δv+Δa\\ &&-\frac{(k_{\rm s} \bar{x}+k_{\rm s2} \bar{x}^2 )}{m}-(γ\bar{v}+γ_2 \bar{v}^2 )+\bar{a}-\frac{\bar{dv}}{dt} \tag{30} \end{eqnarray} \begin{equation} \frac{dΔa}{dt}={Δ{\bf u}}+\bar{u}-\frac{d\bar{a}}{dt} \tag{31} \end{equation} となります。つまり、\({\bf A}\)、\({\bf B}\)と\({\bf c}\)は、 \begin{equation} {\bf A}=\left( \begin{array}{ccc} 0&1&0\\ -(k_{\rm s}+2k_{\rm s2}\bar{x})/m&-(γ+2γ_2 \bar{v})&1\\ 0&0&0 \end{array} \right) \tag{32} \end{equation} \begin{equation} {\bf B}=\left( \begin{array}{c} 0\\0\\1 \end{array} \right) \tag{33} \end{equation} \begin{equation} {\bf c}=\left( \begin{array}{c} \bar{v}-d\bar{x}/dt\\ -(k_{\rm s} \bar{x}+k_{\rm 2}\bar{x}^2)/m-(γ\bar{v}+γ_2 \bar{v}^2 )+\bar{a}-d\bar{v}/dt\\ \bar{u}-(d\bar{a})/dt \end{array} \right) \tag{34} \end{equation} となります。

また、ラグランジュ関数\(L\)を \begin{equation} L=\frac{1}{2}{\bf u}^{\rm T} {\bf u} \tag{35} \end{equation} とします。すると、 \begin{equation} {\bf Q}=\left. \frac{∂^2 L}{∂{\bf z}^2}\right|_{\bar{\bf z},\bar{\bf u}}={\bf 0} \tag{36} \end{equation} \begin{equation} {\bf R}=\left.\frac{∂^2 L}{∂{\bf u}^2}\right|_{\bar{\bf z},\bar{\bf u}}=1 \tag{37} \end{equation} \begin{equation} {\bf z}_{\rm m}=0 \tag{38} \end{equation} \begin{eqnarray} {\bf u}_{\rm m}&=&-{\bf R}^{-1} \left.\frac{∂L}{∂{\bf u}^{\rm T}}\right|_{\bar{\bf z},\bar{\bf u}}\\ &=&-{\bf R}^{-1}\bar{\bf u}\\ &=&-\bar{\bf u}\tag{39} \end{eqnarray} と計算できます。このとき\(L\)は、 \begin{eqnarray} L=\frac{1}{2}({∆{\bf u}}-{\bf u}_{\rm m})^{\rm T} {\bf R}({∆{\bf u}}-{\bf u}_{\rm m})+L_0 (\bar{\bf z},\bar{\bf u})\\ =\frac{1}{2} ({∆{\bf u}}+\bar{\bf u})^{\rm T} {\bf R}({∆{\bf u}}+\bar{\bf u})+L_0 (\bar{\bf z},\bar{\bf u}) \tag{40} \end{eqnarray} と\(∆{\bf u}\)に関する二次形式になります。ここで\(L_0 (\bar{\bf z},\bar{\bf u})\)は、 \begin{equation} L_0 (\bar{\bf z},\bar{\bf u},t)=-\frac{1}{2} \bar{\bf u}^{\rm T} {\bf R}^{-1} \bar{\bf u}+L(\bar{\bf z},\bar{\bf u},t) \tag{41} \end{equation} と表せます。

4. 軌道拘束された二次元質点系

曲線定規などで軌道を拘束したうえで描く図形の描画速度についての心理物理実験が行われ、描画速度と曲線の曲率に負の相関が見出されています（Lacquanti et al. 1983など）。その結果と比較するために、軌道拘束された二次元質点システムを考えます。位置ベクトルを\({\bf x}\)、速度ベクトルを\({\bf v}\)、加速度ベクトルを\({\bf a}\)、躍度ベクトルを\({\bf g}\)とすると、 \begin{equation} \frac{d{\bf x}}{dt}={\bf v} \tag{42} \end{equation} \begin{equation} \frac{d{\bf v}}{dt}={\bf a} \tag{43} \end{equation} \begin{equation} \frac{d{\bf a}}{dt}={\bf g} \tag{44} \end{equation} と書けます。質点の軌道が \begin{equation} {\bf x}-{\bf h}(s)={\bf 0} \tag{45} \end{equation} で定義される曲線Sに拘束されているとします。曲線Sに沿った長さを\(s\)とすると、 \begin{equation} \frac{d{\bf x}}{dt}=\frac{d{\bf x}}{ds}\frac{ds}{dt} \tag{46} \end{equation} と書けます。さらに、フレネ・セレの公式より、 \begin{equation} \frac{d{\bf x}}{ds}={\bf t} \tag{47} \end{equation} \begin{equation} \frac{d^2 {\bf x}}{ds^2}=κ{\bf n} \tag{48} \end{equation} \begin{equation} \frac{d^3 {\bf x}}{ds^3}=\frac{dκ}{ds} {\bf n}-κ^2 {\bf t} \tag{49} \end{equation} と表されます。ここで、\({\bf t}\)は接線ベクトル、\({\bf n}\)は法線ベクトル、\(κ\)は曲率です。

\({\bf x}\)を時間微分すると、\(s\)の定義により、 \begin{eqnarray} \frac{d{\bf x}}{dt}&=&\frac{ds}{dt}{\bf t}\\ &=&v{\bf t} \tag{50} \label{eq50}\end{eqnarray} が得られます。ここで、曲線Sに沿った速度を\(v\)とします。つまり、 \begin{equation} \frac{ds}{dt}=v \tag{51}\label{eq51} \end{equation} です。 式\ref{eq50}をさらに時間微分すると、 \begin{eqnarray} \frac{d^2 {\bf x}}{dt^2}&=&\frac{d}{dt}\left( \frac{d{\bf x}}{ds} \frac{ds}{dt}\right)\\ &=&\frac{ds}{dt} \frac{d}{dt} \left(\frac{d{\bf x}}{ds}\right)+\frac{d{\bf x}}{ds} \frac{d}{dt} \left(\frac{ds}{dt}\right)\\ &=&\frac{d^2 {\bf x}}{ds^2} \left(\frac{ds}{dt}\right)^2+\frac{d{\bf x}}{ds} \frac{d^2 s}{dt^2 }\\ &=&κ\left(\frac{ds}{dt}\right)^2 {\bf n}+\frac{d^2 s}{dt^2} {\bf t}\\ &=&κv^2 {\bf n}+a {\bf t} \tag{52} \label{eq52}\end{eqnarray} が得られます。ここで、曲線Sに沿った加速度を\(a\)としました。つまり、 \begin{equation} \frac{d^2 s}{dt^2}=\frac{dv}{dt}=a \tag{53}\label{eq53} \end{equation} です。

式\ref{eq51}をさらにもう一度時間微分することにより、 \begin{eqnarray} \frac{d^3 {\bf x}}{dt^3}&=&\frac{d}{dt}\left( \left(\frac{d^2 {\bf x}}{ds^2}\right) \left(\frac{ds}{dt}\right)^2+\frac{d{\bf x}}{ds} \frac{d^2 s}{dt^2}\right)\\ &=&\left(\frac{ds}{dt}\right)^2\frac{d}{dt} \left(\frac{d^2 {\bf x}}{ds^2}\right)+\frac{d^2 x}{ds^2} \frac{d}{dt} \left(\frac{ds}{dt}\right)^2\\ &&+\frac{d^2 s}{dt^2} \frac{d}{dt} \left(\frac{d{\bf x}}{ds}\right)+\frac{dx}{ds} \frac{d}{dt} \left(\frac{d^2 s}{dt^2}\right)\\ &=&\left(\frac{ds}{dt}\right)^3 \frac{d^3 x}{ds^3}+2\frac{d^2 {\bf x}}{ds^2}\frac{ds}{dt} \frac{d^2 s}{dt^2}+\frac{d^2 s}{dt^2}\frac{ds}{dt} \frac{d^2 {\bf x}}{ds^2}+\frac{d{\bf x}}{ds} \frac{d^3 s}{dt^3}\\ &=&\left(\frac{ds}{dt}\right)^3 \frac{d^3 {\bf x}}{ds^3}+3 \frac{ds}{dt} \frac{d^2 s}{dt^2}\frac{d^2 {\bf x}}{ds^2}+\frac{d^3 s}{dt^3} \frac{dx}{ds}\\ &=&\left(\frac{ds}{dt}\right)^3 \left(\frac{dκ}{ds} {\bf n}-κ^2 {\bf t}\right)+3 \frac{ds}{dt} \frac{d^2 s}{dt^2 } κ{\bf n}+\frac{d^3 s}{dt^3}{\bf t}\\ &=&\left(\frac{dκ}{ds} \left(\frac{ds}{dt}\right)^3+3κ \frac{ds}{dt} \frac{d^2 s}{dt^2}\right){\bf n}+\left(\frac{d^3 s}{dt^3}-κ^2 \left(\frac{ds}{dt}\right)^3 \right){\bf t}\\ &=&(κ’v^3+3κva){\bf n}+\left(\frac{da}{dt}-κ^2 v^3\right){\bf t} \tag{54}\label{eq54} \end{eqnarray} を得ます。

さらに、 \begin{equation} u=\frac{da}{dt}-κ^2 v^3 \tag{55} \label{eq55}\end{equation} とします。式\ref{eq53}、\ref{eq54}、\ref{eq55}を式\ref{eq52}に代入すると \begin{equation} \frac{d^3 x}{dt^3 }=(κ’v^3+3κva){\bf n}+u{\bf t} \tag{56} \end{equation} を得ます。ここで、\(dκ/ds=κ’\)としました。

 

ここで、評価関数\(J\)を接線方向の躍度の大きさ\(u^2=(\frac{da}{dt}-κ^2 v^3 )^2\) の時間積分とします。つまり、 \begin{equation} J＝\frac{1}{2}\int_0^{t_{\rm f}}u^2 dt \tag{57} \end{equation} ととります。ラグランジアン\(L\)は、 \begin{equation} L=\frac{1}{2} u^2 \tag{58} \end{equation} となります。ここで、\({\bf z}=(s\quad v\quad a)^{\rm T}\)とします。しかし、式\ref{eq55}が\(v\)と\(a\)に対して非線形であることが問題です。このままでは、線形状態空間モデルの最適制御解法を適用できません。この問題を克服するため、\({\bf z}=(s \quad v \quad a)^{\rm T}\)と\(u\) に対して、\({\bf z}=\bar{\bf z}+{∆{\bf z}}\)、\(u=\bar{u}+∆u\)とし、式\ref{eq53}、\ref{eq54}、\ref{eq55}を線形化（局所線形化）すると、 \begin{equation} \frac{d{∆s}}{dt}=∆v+\bar{v}-\frac{d\bar{s}}{dt} \tag{59}\label{eq59} \end{equation} \begin{equation} \frac{d{∆v}}{dt}=∆a+\bar{a}-\frac{d\bar{v}}{dt} \tag{60}\label{eq60} \end{equation} \begin{equation} \frac{d{∆a}}{dt}=∆u+3κ^2 \bar{v}^2 {∆v}+\bar{u}+κ^2 \bar{v}^3-\frac{d\bar{a}}{dt} \tag{61}\label{eq61} \end{equation} を得ます。式\ref{eq59}、\ref{eq60}、\ref{eq61}より、 \begin{equation} {\bf A}= \left(\begin{array}{ccc} 0& 1& 0\\ 0& 0& 1\\ 0& 3κ^2 \bar{v}^2& 0 \end{array}\right) ,\, {\bf B}=\left( \begin{array}{c} 0\\ 0\\ 1\\ \end{array} \right),\, {\bf c}= \left(\begin{array}{c} \bar{v}-\frac{d\bar{s}}{dt}\\ \bar{a}-\frac{d\bar{v}}{dt}\\ \bar{u}+κ^2 {\bar v}^3-\frac{d\bar{a}}{dt} \end{array}\right) \tag{62} \end{equation} を得ます。また、 \begin{eqnarray} {\bf p}(\bar{z},\bar{u})&=&\left. \frac{∂L}{∂{\bf z}}\right|_{\bar{\bf z},\bar{u}}\\ &=&{\bf 0} \tag{63} \end{eqnarray} を得ます。さらに、 \begin{eqnarray} {\bf Q}(\bar{\bf z},\bar{\bf u})&=&\left. \frac{∂^2 L}{∂{\bf z}^2}\right|_{\bar{\bf z},\bar{u}}\\ &=&{\bf 0} \tag{64} \end{eqnarray} を得ます。最後に、 \begin{equation} \frac{∂L}{∂u}=u \tag{65} \end{equation} なので、 \begin{eqnarray} {\bf q}(\bar{\bf z},\bar{u})&=&\left.\frac{∂L}{∂u}\right|_{\bar{\bf z},\bar{u}} &=&\bar{u}\tag{66} \end{eqnarray} \begin{eqnarray} {\bf R}(\bar{\bf z},\bar{u})&=&\left. \frac{∂^2 L}{∂u^2}\right|_{\bar{\bf z},\bar{u}} &=&1 \tag{67} \end{eqnarray} \begin{equation} {\bf z}_{\rm m}={\bf 0} \tag{68} \end{equation} \begin{eqnarray} u_{\rm m}&=&-{\bf R}^{-1} {\bf q} &=&-1^{-1} u &=&-u \tag{69} \end{eqnarray} と書けます。線形状態空間モデルを定義する行列とベクトルが得られたので、これらに対し最適制御解法を適用します。以下に曲線Sが直線、円弧、対数螺旋の場合について考察します。

　なお、先行研究であるHuh and Sejnowski (2015)は、ラグランジアンとして\(u^2/2\)ではなくて、\(u^2+(κ’v^3+3κva)^2\)を用いています。つまり、曲線成分の接線成分だけではなく、鉛直成分も考慮しているわけです（全体に掛かる係数1/2の違いは結果に影響を与えません）。ところがLacquanti et al. (1983)などの心理物理実験の多くでは、フリーハンド描画ではなく曲線定規を用いています。この場合は、鉛直成分を考慮するのは適当ではないと考えられるので、本節ではラグランジアンとして式58を用いることにしました。フリーハンド描画で鉛直成分も考慮する場合も同様に計算できます。両者の結果に定性的な差はないようです。

4.1 直線の場合

\begin{equation} κ=0 \tag{70} \end{equation} \begin{equation} κ’=0 \tag{71} \end{equation} より、一次元質点系に帰着します（二章）。

4.2 円弧の場合 \begin{equation} κ=κ_0≠0 \tag{72} \end{equation} \begin{equation} κ’=0 \tag{73} \end{equation} より \begin{equation} {\bf A}= \left(\begin{array}{ccc} 0& 1& 0\\ 0& 0& 1\\ 0& 3κ_0^2 \bar{v}^2& 0 \end{array}\right) ,\, {\bf B}=\left( \begin{array}{c} 0\\ 0\\ 1\\ \end{array} \right),\, {\bf c}= \left(\begin{array}{c} \bar{v}-\frac{d\bar{s}}{dt}\\ \bar{a}-\frac{d\bar{v}}{dt}\\ \bar{u}+κ_0^2 {\bar v}^3-\frac{d\bar{a}}{dt} \end{array}\right) \tag{74} \end{equation} \begin{equation} {\bf p}={\bf 0} \tag{75} \end{equation} \begin{equation} {\bf q}=\bar{u} \tag{76} \end{equation} \begin{equation} {\bf Q}＝{\bf 0} \tag{77} \end{equation} \begin{equation} {\bf R}=1 \tag{78} \end{equation} となります。Figure 3に円弧上の到達問題の躍度最小解を与えます。直線状の躍度最小の軌道を精度良く再現していることが分かります。

 

 

4.3 対数螺旋の場合

74170832によると曲線Sの曲率\(κ\)が、 \begin{equation} κ=κ_0 e^{αθ} \tag{79} \label{eq79}\end{equation} の形で与えられるとき、曲線Sは対数螺旋になります。ここで、\(θ\)は曲線Sに沿った距離\(s\)と曲率\(κ\)を用いて \begin{equation} κ=\frac{dθ}{ds} \tag{80} \end{equation} で定義される変数です。曲率\(κ\)が一定（円）の場合、\(θ\)は円の中心から測った角度になります。

練習問題3.3

対数螺旋のとき、曲率が式\ref{eq79}で与えられることを示しなさい。

Huhと Sejnowski (2015) によると、対数螺旋に対して、躍度最小の特殊解が存在します。その特殊解においては、\(v∝κ^{-2/3}\) が成り立ちます。つまり、 \begin{equation} v=v_0 e^{-2αθ/3} \tag{81} \label{eq81}\end{equation} が成立します。

 

練習問題3.4

対数螺旋のとき、躍度最小解において\(v∝κ^{-2/3}\)が成り立つことを示しなさい。

対数螺旋の場合は、躍度最小の特殊解が、以下のように求められます。 \begin{equation} s=\frac{v_0 t}{3} \left(\left(-\frac{κ_0 v_0 α}{3} t+1\right)^2-\left(-\frac{κ_0 v_0 α}{3} t+1\right)+2\right) \tag{80}\label{eq80} \end{equation} \begin{equation} κ=κ_0 \left(-\frac{κ_0 v_0 α}{3} t+1\right)^{-3} \tag{83}\end{equation} \begin{equation} v=v_0 \left(-\frac{κ_0 v_0 α}{3} t+1\right)^2 \tag{84} \end{equation} \begin{equation} a=-\frac{2κ_0 v_0^2 α}{3}\left (-\frac{κ_0 v_0 α}{3} t+1\right) \tag{85}\label{eq85} \end{equation}

練習問題3.5

式\ref{eq79}と\ref{eq81}を仮定して、式\ref{eq80}-\ref{eq85}を求めなさい。

4.4正弦対数曲率曲線の場合

正弦対数曲率曲線の曲率は、 \begin{equation} κ=κ_0 e^{ϵ{\rm sin}(νθ)} \tag{86}\label{eq86} \end{equation} で与えられます(Figure 5)。ここで、\(ϵ\)は振幅、\(ν\)は振動数です。この二つのパラメータを変えることでFigure5に見られるような様々な曲線を描くことができます。正弦対数曲率曲線に拘束されたシステムでは、躍度最小解で、速度\(v\)は\(ϵ≪1\)の極限で、 \begin{eqnarray} \frac{v}{v_0} &≅&{\rm exp}(ϵb{\rm sin}(νθ))\\ &=&({\rm exp}(ϵ{\rm sin}(νθ)))^b \tag{87}\label{eq87} \end{eqnarray} と表されます。つまり、 \begin{eqnarray} v&≅&v_0 ({\rm exp}(ϵ{\rm sin}(νθ)))^b\\ &=&v_0 \left(\frac{κ}{κ_0 }\right)^b\\ &∝&κ^b \tag{88} \end{eqnarray} です。

ここで、 \begin{equation} β=-b=\frac{2}{3} \frac{1+\frac{ν^2}{5}}{1+\frac{2ν^2}{5}+\frac{ν^4}{15}} \tag{89}\label{eq89} \end{equation} と表されます。一方、Huh and Sejnowski (2015)では、躍度の垂直成分も考慮して、同様の計算を行い、 \begin{equation} β=\frac{2}{3} \frac{1+\frac{ν^2}{2}}{1+ν^2+\frac{ν^4}{15}} \tag{90}\label{eq90} \end{equation} を得ています。式\ref{eq89}と\ref{eq90}はともに\(ν=0\)で\(β=2/3\)、\(ν=2\)で\(β≅1/3\)、\(ν≫1\)で\(β→0\)を与えます（Figure 6）。

練習問題3.6

式\ref{eq86}を仮定して、オイラー・ラグランジュ方程式を\(ϵ≪1\)で展開して、式\ref{eq87}-\ref{eq89}となることを確かめなさい。

Figure 7に式\ref{eq87}、\ref{eq86}、\ref{eq89}が与える近似解に従って終端条件を与えて作った、接線方向躍度最小解を曲率速度平面にプロットしました。\(ν＝0.4\)においては、解全体にわたって式\ref{eq89}が成り立っています（破線は式\ref{eq89}が与える傾きを示している）。\({\rm sin (νθ)}\)が大きくなるにつれて、近似解からの解離が目立ちます。しかし、\({\rm sin (νθ)}\)が小さい始点と終点の近くは式\ref{eq89}が与える傾きによく一致しています（Figure 6で示した例では、終点は\({\rm sin (νθ)}\)がゼロになるような値を取っています）。

6.　最適制御解法の脳への実装 　

本章では、局所線形化手法を用いて最適制御解法を拡張し、非線形システムに適用しました。非線形な力が働く場合（第3節）、曲線に拘束されている場合（第4節）のどちらでも、局所線形化された最適制御解法を用いると、修正量が繰り返し回数に応じて急速に減少し、高精度の躍度最小解が得られます。対数螺旋に拘束される場合には、解析解によく一致しました。また、正弦対数曲率曲線に拘束される場合は、その近似解析解と良く一致しました。これらの対数螺旋の解析解や正弦対数曲率曲線の近似解析解は、心理物理実験の結果と良く整合することが知られています。

従来、これらの解析解や近似解を得るためには微分や積分など多くの数学技法を用いる必要がありました。一方、本章にまとめた局所線形化された最適制御解法（式\ref{eq10}‐\ref{eq23}）を用いれば、自動的に到達問題の躍度最小解を任意の精度で得られることは大変便利です。また、その躍度最小解が、心理物理実験の結果と一致することは上でも述べたとおりです。

これらの事実から私は、人間および動物の運動制御を司る部分に、最適制御解法（式\ref{eq10}‐\ref{eq23}）を実行するサブルーチンが実装されているはずだと考えるようになりました。そしてこの式\ref{eq10}‐\ref{eq23}を天馬方程式と呼ぶことにしました。手塚治虫による漫画「鉄腕アトム」の主人公アトムは、天才科学者「天馬博士」が作ったロボットでした。その運動制御部分には天馬方程式が実装されていたはずだと私は考えました。 　

では、脳に天馬方程式を実装しようとすればどうすればいいでしょうか？天馬方程式は、行列積と逆行列計算 からできています。行列積については、神経細胞がシナプスで行うとされている入力信号と重み係数のアナログ的な積和演算を多数回使えば、実行可能です。一方、逆行列計算（もしくは線形一次方程式の解法）はそれほど単純ではありません。例えばガウスの掃き出し法には、除算が必要です。これを上記のアナログ的積和演算を用いて行うことは困難です。除算の結果はゼロの近くで発散するため、すべての計算範囲で精度よく計算するには特別な注意が必要なのです。共役勾配法を用いても割り算が必要なことは変わりません。 　そこで、私はいろいろな手法を試してみました。例えば一旦対数に変換し、除算を引き算に代えて実行するなどです。ある種の神経細胞の対数的な反応をすることが知られています（Izhkevich, 2007）。しかし、対数化された数を真数に戻すときに必要な指数的な反応を示す神経細胞が見当たりません。いくつかの神経細胞を組み合わせて指数関数を作ることには除算と同様の困難がありました。 　

上記のような試行錯誤の後、私が行きついた結論は、デジタル計算機の手法を用いることでした。デジタル計算機において除算（もしくは、逆数）は、短語長の表にテイラー係数を格納しておき、テイラー展開で計算します。脳における計算は256諧調程度の低精度で実行されていることを考えると、8ビットの表検索で逆数を求めることが可能です。語長8ビットのメモリを構成することなら、「アナログ的な積和演算」を組み合わせて構成可能です。また、神経細胞は、トランジスターと同様の非線形素子です。その非線形性を用いて論理ゲート（AND、ORなど）を構成すること考えました、デジタル計算機がトランジスターから構成されているように。 　

そこで、第8章で、神経細胞で論理ゲートとそれらを組み合わせて作られる基本的な計算素子を神経細胞で構築したいと思います。その前に次章（第4章）では、動物が使っている駆動装置である筋腱繊維の性質を議論し、なぜ躍度最小となるような制御法が進化の過程で獲得されたのかを考察します。また、第5章で最適制御解法を時間離散化します。人間の脳では、連続的に流れる時間は直接には取り扱えません。したがって、離散的な時刻での値の組み合わせで表現されているはずだからです。そして、離散化された天馬方程式を使って、二次精度の解を求められることを示します。さらに、第6章で状態空間モデルを定義する行列を決めるシステム同定について議論し、7章でその手法を用いた脳の学習過程を心理物理実験の結果と比較します。　

1) Flash, T. and Hogan, N. 1985, The coordination of arm movements: an experimentally confirmed mathematical model, Journal of Neuroscience, 5, 1688-1703.

2) Todorov, E. and Jordan, M.I., 1998, Smoothness maximization along a predefined path accuracy predicts the speed profiles of complex arm movements, Journal of neurophysiology, 80, 89-101.

3) Viviani, P., and Flash, T., 1995, Minimum-jerk, two thirds power law, and isochrony: converging approaches to movement planning, The journal of Experimental Psychology: Human Perception and Performances, 21, 32-53.

4) Lacquanti, F., Terzuolo, C., and Viviani, P. The law relating the kinematic and figural aspects of drawing movements, Acta Psychologica, 54, 115-130.

5) Huh, D. and Sejnowski, T.J., 2015, Spectrum of power laws for curved hand movements, Publications of National Academy of Science, E3950–E3958, www.pnas.org/cgi/doi/10.1073/pnas.1510208112

6) Izhkevich, E.M., 2007, Dynamical Systems in Neuroscience, ch. 4, p89-126.