筋腱繊維で駆動されるシステムー脳の設計図を求めて(4)ー

動物は駆動に筋腱繊維を使います。筋腱繊維は特異な駆動装置です。人間を含む動物の運動が、躍度を最小とするように制御されている原因の一つは、筋腱繊維の特殊な性質にあるはずです。

１．筋腱繊維

筋腱繊維は、能動的に収縮する筋肉とバネのように受動的に伸び縮みする腱が直列につながってできています（Figure 1）。筋肉は神経刺激\(\alpha\) (\(0≤\alpha≤1\))を受けると力\(F^{\rm M}\)を発生して収縮します。\(F^{\rm M}\)は、神経刺激\(\alpha\)の他に、筋肉の長さ\(L^{\rm M}\)、筋肉の収縮速度\(V^{\rm M}\)に依存します。ここで、

\begin{equation}
V^{\rm M}=-\frac{dL^{\rm M}}{dt} \tag{1}
\end{equation}

です。また、力\(F^M\)は収縮方向を正とします。その関数の概形は、Figure 2に与えられています。まず、力はが\(0.5L_0^{\rm M}<L^{\rm M}<1.5L_0^{\rm M}\)の間でのみ力を発生します（Figure 2a）。そして\(L^{\rm M}\sim L_0^{\rm M}\)周辺では力\(F^{\rm M}\)は\(L^{\rm M}\)にはあまり依存しません。そして\(F^{\rm M}\)は神経刺激\(α\) に比例し（Figure 2c）、\(α=1\)で\(V^{\rm M}=0\)の時最大静止力\(F_0^{\rm M}\)を発生します（Figure 2b）。また、\(V^{\rm M}\)が有限で正（収縮）のとき、\(F^{\rm M}\)は減少し、\(V^{\rm M}=V_{\rm M}^0\)の時ゼロとなります。\(V_0^{\rm M}\)は筋肉の最大収縮速度と呼ばれています。一方、腱が発生する力\(F^{\rm T}\)は腱の長さ\(L^{\rm T}\)にのみ依存します（Figure 2a、b）。\(F^{\rm T}\)は\(L^{\rm T}<L_0^{\rm T}\)ではゼロ（弛んでいる）ですが、\(L^{\rm T}>L_0^{\rm T}\)では\(L^{\rm T}-L_0^{\rm T}\)に比例して、バネのように増加します。

上記のような性質を持った筋腱繊維のモデル方程式をZajac (1989)に従って作ってみましょう。まず、筋肉と腱は直列につながっていますから \begin{equation} F^{\rm M}=F^{\rm T}=F \tag{2}\label{eq2} \end{equation} と書けます。このように、筋肉が力を発生して収縮すると式\ref{eq2}に従って腱にも力がかかるので、腱が伸びます。この性質をモデルで表現する必要があります。そこで、筋腱繊維全体の長さを\(L^{\rm MT}\)と置くと、 \begin{equation} L^{\rm MT}=L^{\rm T}+L^{\rm M} \tag{3}\label{eq3} \end{equation} と書けます。式\ref{eq3}の両辺を時間微分すると \begin{equation} -V^{\rm MT}=-V^{\rm T}-V^{\rm M} \tag{4}\label{eq4} \end{equation} を得ます。ここで、 \begin{equation} V^{\rm MT}=-\frac{dL^{\rm MT}}{dt},\, V^{\rm T}=-\frac{dL^{\rm T}}{dt} \tag{5} \end{equation} と定義しました。式\ref{eq4}より \begin{equation} V^{\rm T}=V^{\rm MT}-V^{\rm M} \tag{6}\label{eq6} \end{equation} を得ます。次に、腱をバネとみなしそのバネ定数を\(K^{\rm T}\)と置きます。つまり \begin{equation} F^{\rm T}=K^{\rm T} {\rm Max}(L^{\rm T}-L_0^{\rm T},0) \tag{7} \end{equation} とします。つまり、Figure 3cの破線のように折れ線で近似します。\(L^{\rm T}>L_0^{\rm T}\)のとき、力の時間微分は、 \begin{eqnarray} \frac{dF}{dt}&=&\frac{dF^T}{dt}\\ &=&\frac{d}{dt} (K^{\rm T} (L^{\rm T}-L_0^{\rm T} ))\\ &=&K^{\rm T} \frac{dL^T}{dt}\\ &=&-K^{\rm T} V^{\rm T} \tag{8}\label{eq8} \end{eqnarray} と計算できます。式\ref{eq8}に式\ref{eq6}を代入します \begin{eqnarray} \frac{dF}{dt}&=&-K^{\rm T} (V^{\rm MT}-V^{\rm M} )\\ &=&-K^{\rm T}V_0^{\rm M} \left(\frac{V^{\rm MT}}{V_0^{\rm M}}-\frac{V^{\rm M}}{V_0^{\rm M}} \right) \tag{9}\label{eq9} \end{eqnarray}

を得ます。最後に筋肉について考えましょう。筋肉が発生する力\(F^{\rm M}\)は、筋肉の長さ\(L^{\rm M}\)、筋肉の収縮速度\(V^{\rm M}\)、そして神経刺激の強度\(α\)が関係します。つまり、 \begin{equation} \frac{F^{\rm M}}{F_0^{\rm M}}=f_{\rm M1} \left(\frac{L^{\rm M}}{L_0^{\rm M}},\frac{V^{\rm M}}{V_0^{\rm M}},α\right) \tag{10} \end{equation} と書けます。ここで、\(F^{\rm M}\)、\(L^{\rm M}\)、\(V^{\rm M}\)はそれぞれ\(F_0^{\rm M}\)、\(L_0^{\rm M}\)、\(V_0^{\rm M}\)で規格化（無次元化）しました。この4者の関係を、筋肉の最大収縮速度\(V_0^{\rm M}\)で規格化された\(V^{\rm M}\)について解いた形に変形しましょう。 \begin{equation} \frac{V^{\rm M}}{V_0^{\rm M}}=f_{\rm M2} \left(\frac{L^{\rm M}}{L_0^{\rm M}},\frac{F^{\rm M}}{F_0^{\rm M}},α\right) \tag{11}\label{eq11} \end{equation} 式\ref{eq11}を式\ref{eq9}に代入すると \begin{equation} \frac{dF}{dt}=-K^{\rm T} V_0^{\rm M} \left(\frac{V^{\rm MT}}{V_0^{\rm M}}-f_{\rm M2} \left(\frac{L^{\rm M}}{L_0^{\rm M}}, \frac{F^{\rm M}}{F_0^{\rm M}},α\right)\right) \tag{12}\label{eq12} \end{equation} 式\ref{eq12}は、腱の伸びと筋肉の力の発生を決める方程式なので、筋腱方程式と呼ぶことにします。ここで、筋肉について\(L^{\rm M}\sim L_0^{\rm M}\)周辺の平坦な部分のみを取り扱うことにし（Figure　3aの破線）、\(V_{\rm M}\)依存性（Figure 3bの破線）も考慮して \begin{equation} \frac{F_M}{F_0^M}=f_{\rm M1} \left(\frac{L^{\rm M}}{L_0^{\rm M}},\frac{F^{\rm M}}{F_0^{\rm M}},α\right)=α-\frac{V^{\rm M}}{V_0^{\rm M}} \tag{13}\label{eq13} \end{equation} と近似します（\(α>V^{\rm M}/V_0^{\rm M}\)）。式\ref{eq13}を\(V^{\rm M}/V_0^{\rm M}\)について解くと、 \begin{equation} \frac{V^{\rm M}}{V_0^{\rm M}}=f_{\rm M2} \left(\frac{L^{\rm M}}{L_0^{\rm M}},\frac{F^{\rm M}}{F_0^{\rm M}},α\right)=α-\frac{F^{\rm M}}{F_0^{\rm M}} \tag{14}\label{eq14} \end{equation} を得ます。式\ref{eq14}を式\ref{eq12}に代入して \begin{eqnarray} \frac{dF}{dt}&=&-K^{\rm T} V_0^{\rm M} \left(\frac{V^{\rm MT}}{V_0^{\rm M}}-\left(α-\frac{F}{F_0^{\rm M}}\right)\right)\\  &=&K^{\rm T} V_0^{\rm M} \left(-\frac{V^{\rm MT}}{V_0^{\rm M}}-\frac{F}{F_0^{\rm M}}+α\right) \tag{15}\label{eq15} \end{eqnarray} を得ます。

2．なぜ躍度最小制御なのか？

式\ref{eq15}を特性時間\(t_{\rm c}=L_0^{\rm M}/V_0^{\rm M}\)を導入して無次元化すると \begin{equation} \frac{d}{d(t/t_{\rm c})}\left(\frac{F}{F_0^{\rm M}}\right)=k\left(-\frac{V^{\rm MT}}{V_0^{\rm M}}-\frac{F}{F_0^{\rm M}}+α\right) \tag{16} \end{equation} を得ます。ここで、無次元化バネ定数\(k\)は、 \begin{equation} k=\frac{K^{\rm T} V_0^{\rm M} τ_{\rm c}}{F_0^{\rm M}}=\frac{K^{\rm T} L_0^{\rm M}}{F_0^{\rm M}} \tag{17} \end{equation} で与えられます。実際の筋肉の\(k\)は以下のように評価できます。Zajac (1989)によると腱の応力\(σ^{\rm T}\)（単位面積当たりの力）は、 \begin{eqnarray} σ^{\rm T}&=&σ_0^T=32\, {\rm MPa}\quad (ε=0.033)\\ σ^{\rm T}&=&112\, {\rm MPa}\quad (ε=0.1) \end{eqnarray} と与えられています。ここで、\(ε^{\rm T}=(L^{\rm T}-L_0^{\rm T} )/L_0^{\rm T}\)です。したがって、 \begin{eqnarray} \frac{dσ^{\rm T}}{dε}≅\frac{112-32}{0.1-0.033}=1.2\, {\rm GPa} \tag{18} \end{eqnarray} です。これを用いて、 \begin{eqnarray} K^{\rm T}&=&\frac{A}{L_0^{\rm T}} \frac{dσ^{\rm T}}{ dε}\\ &≅&1.2\, [{\rm Gpa}]\frac{A}{L_0^{\rm T}} \tag{19} \end{eqnarray} を得ます。ここで、\(A\)は筋腱繊維の断面積です。一方、 \begin{equation} F_0^{\rm M}=Aσ_0^{\rm M} \tag{20} \end{equation} であるので、 \begin{eqnarray} k&=&\frac{K^{\rm T} L_0^{\rm M}}{F_0^{\rm M}}\\ &=&\frac{A}{L_0^{\rm T}} \frac{dσ^{\rm T}}{ dε} \frac{L_0^{\rm M}}{Aσ_0^{\rm M}}\\ &=&\frac{1.2\, {\rm GPa}}{32\, {\rm MPa}} \frac{L_0^M}{L_0^T}=\frac{37.5}{η} \tag{21}\label{eq21} \end{eqnarray} を得ます。ここで、\(η=L_0^{\rm T}/L_0^{\rm M}\)は腱筋長比です。腱筋長比は、動物や筋の種類によって大きく違います（Table 1）。Table 1にあるようにヒトの下肢にある長大な筋肉以外は、\(η\sim 3\)です。つまり、式\ref{eq21}に見るように、\(k\)が\(\sim 10\)と大きめです。 無次元化された筋腱方程式（式\ref{eq15}）においては、無次元化された量（\(F/F_0^{\rm M}\)、\(V^{\rm MT}/V_0^{\rm M}\)など）は多くの場合1の桁であることが期待されます。それが実現するためには、式\ref{eq13}の右辺の括弧の中の量が\(1/k\sim 0.1\)と小さめであることが要求されます。右辺の括弧の中が1の桁でそれが\(τ_{\rm c}\)程度の間維持されると、\(k\sim 10\)のために\(F≫F_0^{\rm M}\)となって大きな力が筋腱繊維にかかります。これは、修復不可能な損傷の原因になり得ます。このような考察から、式\ref{eq13}の右辺ができるだけ小さくなるようにシステムを制御するのは、理にかなっていることが分かります。

では、\(η\sim 10\)の筋腱繊維ではどうでしょうか？Table 1を見ると、ヒラメ筋や腓腹筋、ハムストリング筋の一部が\(η\sim 10\)を持っています。これらは、ヒトの下肢にあり、ヒトが直立二足歩行を始めるにあたって、重力に抗して体重を支え、激しい運動中にもバランスをうまくとるために発達して長大化した結果、大きな\(η\)を持つようになったのでしょう。ヒラメ筋や腓腹筋はアキレス腱に結合しています。アキレス腱はヒトの体で損傷しやすい部分であることは有名なアキレウスの神話でも示されています。この部分に急激な力の変化（要するに大きな躍度）がかかるとその長さの急激な変化のため回復不能な損傷を受けてしまいます。このような事態を避けるため、躍度ができるだけ小さくなるように体を制御する方法が進化の過程で獲得されたと考えられます。

このように、いずれの場合も躍度が最小になるように運動を制御することは理にかなっていると考えられます。

3. バネ筋腱質点系 　

第二章で取り扱ったバネ質点系を筋腱繊維で駆動することにしましょう（Figure 4）。二つの壁の間の距離は\(x=0\)のとき\(L^{\rm MT}=L_0^{\rm MT}\)となるように定めます。このようなバネ筋腱質点系の運動（状態）方程式は二章の式22、23、24と同じです。 \begin{equation} \frac{dx}{dt}=v \tag{22}\label{eq22} \end{equation} \begin{equation} \frac{dv}{dt}=-k_{\rm s} x-γv+a \tag{23}\label{eq23} \end{equation} \begin{equation} \frac{da}{dt}=g \tag{24}\label{eq24} \end{equation} ここで、有効加速度（力を質量で割ったもの）の時間微分\(g\)を \begin{equation} g=\frac{K^{\rm T} V_0^{\rm M}}{m} \left(-\frac{V^{\rm MT}}{V_0^{\rm M}}-\frac{F}{F_0^{\rm M}}+α\right) \tag{25}\label{eq25} \end{equation} で与え、\(g\)そのものではなく筋肉刺激\(α\)を通して制御します。 さて、式\ref{eq22}、\ref{eq23}、\ref{eq24}と評価関数 \begin{equation} J=\frac{1}{2} \int_0^{t_{\rm f}}g^2 dt \tag{26} \end{equation} を与えて、躍度最小の到達問題を解けば、最適な（\(J\)が最小の）\(g\)を得ます。それを式\ref{eq25}に\(V^{\rm MT}=v\)、\(F=ma\)を代入し、\(α\)について解いた \begin{equation} α=\frac{1}{V_0^{\rm M}}v+\frac{m}{F_0^{\rm M}}a+\frac{m}{K^{\rm T} V_0^{\rm M}} g \tag{27}\label{eq27} \end{equation} に代入して\(α\)を求めればいいことになります。式\ref{eq27}を見ればわかるように、\(α\)はそれぞれのスケール因子（\(V_0^{\rm M}\)、\(F_0^{\rm M}/m\)、そして\(K^{\rm T} V_0^{\rm M}/m\)）で規格された\(v\)、\(a\)、\(g\)の線形結合で表されます。このパラメータでは、\(a\)に近いふるまいを示します。まず、Figure 5は、筋肉の収縮方向（バネが伸びる方向：\(x=0\)から\(x=1\)）への到達問題の躍度最小解を表しています。また、Figure 6は筋肉の伸展方向（バネが縮む方向：\(x=1\)から\(x=0\)）への到達問題の躍度最小解を表しています。筋肉は\(V^{\rm MT}\)=0に対して非対称な反応性を持つので収縮時と伸長時の\(a\)は両者は同じにはなりません。

4. 一軸回転円筒

本節では、筋腱繊維で結合されたシステムとして一軸回転円筒を考えることにしましょう。

4.1　定式化

円筒の半径を\(r\)、高さを\(b\) 、質量を\(M\)で表します（figure 4）。一軸回転円筒は、投球やバッティング運動のダブルスピンモデル（手塚2023）において基本となる運動系とされています（figure 5）。まず、円筒の回転角を\(θ\)、その回転角速度を\(ω\)とすると、 \begin{equation} \frac{dθ}{dt}=ω \tag{28} \label{eq28}\end{equation} と書けます。次に、筋腱繊維の収縮力により円筒に働くトルク\(τ\)は \begin{equation} τ=-bF \frac{θ}{(b^2+r^2 θ^2)^{1/2}} \tag{29} \end{equation} で与えられます。回転軸周りの円筒の慣性モーメント\(I\)は \begin{equation} I=\frac{1}{2} Mr^2 \tag{30} \end{equation} と表されます。回転軸周りの運動方程式は \begin{equation} I \frac{dω}{dt}=τ \tag{31} \end{equation} で与えられます。つまり、 \begin{equation} \frac{1}{2} Mr^2 \frac{dω}{dt}=-bF\frac{θ}{(b^2+θ^2)^{1/2}} \tag{32} \end{equation} \begin{equation} \frac{dω}{dt}=-\frac{2b^2}{Mr^3} \frac{Fθ}{(θ^2+b^2/r^2 )^{1/2}} \tag{33}\label{eq33} \end{equation} となります。さらに、式\ref{eq12}より \begin{equation} \frac{dF}{dt}=K^T V_0^M \left(-\frac{V^{\rm MT}}{V_0^{\rm M}}-\frac{F}{F_0^{\rm M}}+α\right) \tag{34}\label{eq34} \end{equation} を得ます。ここで、筋腱繊維全体の収縮速度\(V^{\rm MT}\)は、 \begin{equation} V^{\rm MT}=-\frac{dL^{\rm MT}}{dt}=-\frac{d}{dt} (b^2+r^2 θ^2 )^{1/2}=-\frac{r^2 θω}{(b^2+r^2 θ^2)^{1/2}} \tag{35}\label{eq35} \end{equation} と計算されます。式\ref{eq35}を式\ref{eq34}に代入すると \begin{equation} \frac{dF}{dt}=u \tag{36} \end{equation} を得ます。ここで、 \begin{equation} u=K^{\rm T} V_0^{\rm M}\left( \frac{rθω}{V_0^M (θ^2+b^2/r^2 )^{1/2}}-\frac{F}{F_0^{\rm M}}+α\right)\tag{37} \end{equation} としました。

4.2 有効躍度ゼロの解 　

躍度最小の極限として、躍度ゼロ（式\ref{eq34}の右辺ゼロ） \begin{equation} \frac{rθω}{(θ^2+b^2/r^2 )^{1/2}} -\frac{F}{F_0^{\rm M}}+α=0 \tag{38}\label{eq38} \end{equation} の場合の解析解を求めてみましょう。式\ref{eq38}より \begin{equation} \frac{dF}{dt}=0 \tag{39} \end{equation} より、 \begin{equation} F=F_0={\rm Const.} \tag{40}\label{eq40} \end{equation} を得ます。また、式\ref{eq40}を式\ref{eq33}に代入すれば \begin{equation} \frac{dω}{dt}=\frac{d^2 θ}{dt^2}=-\frac{2b^2 F_0}{Mr^3} \frac{1}{(θ^2+b^2/r^2 )^{1/2}} θ \tag{41}\label{eq41} \end{equation} を得ます。式\ref{eq41}の両辺に\(dθ/dt\)を掛けると \begin{equation} \frac{d^2 θ}{dt^2} \frac{dθ}{dt}=-\frac{2b^2 F_0}{Mr^3} \frac{1}{(θ^2+b^2/r^2 )^{1/2}} \frac{dθ}{dt} θ \tag{42}\label{eq42} \end{equation} となります。式\ref{eq41}は、 \begin{equation} \frac{1}{2} \frac{d}{dt}\left(\frac{dθ}{dt}\right)^2=-\frac{2b^2 F_0}{Mr^3} \frac{d}{dt} (θ^2+b^2/r^2 )^{1/2} \tag{43}\label{eq43} \end{equation} と変形できます。式\ref{eq43}の両辺を積分すると \begin{equation} \frac{dθ}{dt}^2=-\frac{4b^2 F_0}{Mr^3} (θ^2+b^2/r^2 )^{1/2}+C \tag{44}\label{eq44} \end{equation} つまり、 \begin{equation} ω^2=C-\frac{4b^2 F_0}{Mr^3} (θ^2+(b/r))^2 )^{1/2} \tag{45} \end{equation} を得ます。\(ω=0\)のとき\(θ=θ_{\rm max}\)として、積分定数\(C\)を決めましょう。 \begin{equation} C＝\frac{4b^2 F_0}{Mr^3} (θ_{\rm max}^2+(b/r)^2 )^{1/2} \tag{46}\label{eq46} \end{equation} 式\ref{eq46}を式\ref{eq44}に代入すれば、 \begin{equation} ω^2=2k_1 \left((θ_{\rm max}^2+(b/r)^2 )^{1/2}-(θ^2+(b/r)^2 )^{1/2}\right) \tag{47}\label{eq47} \end{equation} を得ます。ここで、\(k_1=(2b^2 F_0)/(Mr^3)\)です。

4.3 一軸円筒の方程式の局所線形化

式\ref{eq28}、\ref{eq33}、\ref{eq34}を\(\bar{z}+∆z=(\bar{θ}+∆θ, \bar{ω}+∆ω, \bar{a}+∆a)^{\rm T}\)として展開し、\(∆\)付きの変数の一次の項までを取ると、 \begin{equation} \frac{d(\bar{θ}+∆θ)}{dt}=\bar{ω}+∆ω \tag{48}\label{eq48} \end{equation} \begin{eqnarray} \frac{d(\bar{ω}+∆ω)}{dt}&=&-k_1 \frac{(\bar{θ}+∆θ)(\bar{F}+∆F)}{((\bar{θ}+∆θ)^2+(b/r)^2)^{1/2}}\\ &=&-\frac{k_1 \bar{θ}\bar{F}}{(\bar{θ}^2+(b/r)^2)^{1/2}} \\ &&-k_1 \left(-\frac{\bar{θ}^2\bar{F}}{(\bar{θ}^2+(b/r)^2 )^{3/2}}+\frac{\bar{F}}{(\bar{θ}^2+(b/r)^2)^{1/2}}\right){∆θ}\\ &&-k_1 \frac{\bar{θ}}{(\bar{θ}^2+(b/r)^2)^{1/2}}{∆F}\\ &=&-\frac{k_1 \bar{θ}\bar{F}}{(\bar{θ}^2+(b/r)^2 )^{1/2}} -k_1 \frac{-\bar{θ}^2 \bar{F}+(\bar{θ}^2+(b/r)^2) \bar{F}}{(\bar{θ}^2+(b/r)^2 )^{3/2}}{∆θ}\\ &&-k_1 \frac{\bar{θ}}{(\bar{θ}^2+(b/r)^2 )^{1/2}} {∆F}\\ &=&-\frac{k_1 \bar{θ}\bar{F}}{(\bar{θ}^2+(b/r)^2 )^{1/2}} -\frac{k_1 (b/r)^2 \bar{F}}{(\bar{θ}^2+(b/r)^2 )^{3/2}} {∆θ}\\ &&-\frac{k_1 \bar{θ}}{(\bar{θ}^2+(b/r)^2 )^{1/2}} {∆F} \tag{49}\label{eq49} \end{eqnarray} \begin{equation} \frac{d(\bar{F}+∆F)}{dt}=\bar{u}+∆u \tag{50}\label{eq50} \end{equation} を得ます。式\ref{eq48}、\ref{eq49}、\ref{eq50}を整理すると、 \begin{equation} \frac{d{∆θ}}{dt}=∆ω+\bar{ω}-\frac{d\bar{θ}}{dt} \tag{51} \end{equation} \begin{eqnarray} \frac{d{∆ω}}{dt}&=&-\frac{k_1 (b/r)^2 \bar{F} }{(\bar{θ}^2+(b/r)^2 )^{3/2}} {∆θ} -\frac{k_1 \bar{θ}}{(\bar{θ}^2+(b/r)^2 )^{1/2}} {∆F}\\ &&-\frac{k_1 \bar{θ}\bar{F}}{(\bar{θ}^2+(b/r)^2 )^{1/2}} -\frac{d\bar{ω}}{dt} \tag{52} \end{eqnarray} \begin{equation} \frac{d{{∆F}}}{dt}=∆u+\bar{u}-\frac{d\bar{F}}{dt} \tag{53} \end{equation} を得ます。以下のように行列とベクトルを定義すると、 \begin{equation} {\bf A}=\left( \begin{array}{ccc} 0&1&0\\ -\frac{k_1 (b/r)^2 \bar{F}}{(\bar{θ}^2+(b/r)^2 )^{3/2}} &0&-\frac{k_1 \bar{θ}}{(\bar{θ}^2+(b/r)^2 )^{1/2}}\\ 0&0&0 \end{array} \right)\tag{54}\label{eq54} \end{equation} \begin{equation} {\bf B}=\left( \begin{array}{c} 0\\ 0\\ 1 \end{array} \right)\tag{55} \end{equation} \begin{equation} {\bf c}=\left( \begin{array}{c} \bar{ω}-\frac{d\bar{θ}}{dt}\\ -\frac{(k_1 \bar{θ}\bar{F}}{(\bar{θ}^2+(b/r)^2 )^{1/2}} -\frac{d\bar{ω}}{dt}\\ \bar{u}-\frac{d\bar{F}}{dt} \end{array} \right) \tag{56} \end{equation} \begin{equation} \frac{∆\bf z }{dt} ={\bf A}{∆\bf z}+{\bf B}{∆\bf u}+{\bf c} \tag{57}\label{eq57} \end{equation} を得ます。また、 \begin{equation} {\bf Q}={\bf 0},\,{\bf R}=1 \tag{58} \end{equation} とすれば、 \begin{equation} L(x,u,t)=\frac{1}{2}(∆{\bf u}+\bar{\bf u})^{\rm T} R(∆{\bf u}+\bar{\bf u})+L_0 (\bar{\bf z},\bar{\bf u}) \tag{59} \end{equation} となります。 　式\ref{eq57}を用い、手順0－4を使って筋腱繊維で駆動された一軸回転円筒の到達問題を、\(θ_0=-π/2\)、\(θ_{\rm f}=π/2\)、\(t_{\rm f}=5 \,{\rm s}\) の場合について解きました。その結果をfigure 9に与えています。式\ref{eq47}で与える解析解とよく一致しています。

引用文献

1)Zajac, F.C., 1989, Muscle and Tendon: Properties. Models, scaling, and application to biomechanics and motor control, Critical Review in Biomedical Engineering, 17, 359.

2)手塚一志、2023、スポーツ技能上達に関する「Wスピン」運動原理、それを習熟するための手法「創育」について、TEN (Tsunami, Earth, and Networking), 4, 70-75.