動かして学ぶバイオメカニクス #23 〜全身の角運動量とCOP（コード編）〜
前章で述べた全身の角運動量の計算や，その物理的意味を理解していただけただろうか．長いが式はほぼ省略せず，一つ一つ追っていけるようにしたので，諦めずにその意味を考えていただけたらと思う．
ここでは，実際にフォースプレートの代わりに身体側から外力である床反力と力のモーメントを推定し，さらにそれを利用し関節に作用するトルクやCOPを算出し，フォースプレートを利用した計算結果と比較する．ここで行う全身の角運動量の計算ではモーメントの計算を多様に行っているが，その数理の理解と計算を通してモーメントや角運動量に対する理解がより深まることを期待している．また第14章を復習し，比較することもおすすめする．
なお，前章については，これに合わせて大幅に加筆・修正している．また一部，誤りもあり訂正を行っていること，お詫び申し上げる．
はじめに
前章
では，フォースプレートの代わりに身体のダイナミクスから，モーションキャプチャのデータだけで相互作用の力とトルクを推定し，そこから関節に作用するトルク，COPを推定できることを数理的に示した．本章では，台から飛び降りる動作をモーションキャプチャとフォースプレートで計測し，実際にフォースプレートを利用した計算を真値として，推定値と比較する．
運動や計測環境によって精度は異なるだろう．ご自身のデータを使用する場合は，精度をご自身で確かめていただきたい．
次章では，フォースプレートのデータがないので検証ができないが，さらにキック動作で関節に作用するトルクやCOPの計算を行う予定だ．
確認した動作環境とコードについて
いつものように，以下のGoogle ColaboratoryのリンクからブラウザでPythonコードを実行できる．なお，ログインするためには，Googleアカウントが必要となるので注意をされたい．
また，使用する飛び降り動作のモーションキャプチャのサンプルデータ（sample_optitrack_221027.csv）と，フォースプレートのデータ（sample_fp_221027.csv）はこのリンクからダウンロードし，Google Driveをマウントし，マウントしたドライブにそれらをコピーしてご使用いただきたい．すでにダウンロ度済みの方は，再度ダウンロードする必要はない．
また，Google ColaboratoryやGoogle Driveの使用方法は，
を参照されたい．
なお，ここでのコードは，JupyterLab（デスクトップアプリ），Google Colaboratoryで検証したが，とくにグラフィックスの表示結果やコードは，その二つでも若干動作が異なる．動作させるPython環境に依存し，少々修正が必要となるかも知れない．とりあえずJuputerLabのアプリやGoogle Colabratoryで確認したが，環境によってグラフィックスのコードの調整をお願いしたい．ここで示したコードは，JupyterLabで動作するが，Colabでは若干修正している．
なお，環境によっては，
ax.set_aspect('equal')によって，xyz軸間のスケール比（アスペクト比）を調整し，軸間の長さのスケールを等しくすることができるが，Colabではこれが機能しないため，描画範囲をax.set_xlim(), ax.set_ylim(), ax.set_zlim()で指定している．適宜環境に合わせて調整をしていただきたい．
復習
ここで，前章
の復習を行う．前章は導出も多かったので結果をまとめて整理する．なお，前章は訂正・加筆を与えているので，一度ご覧になった方も再度，参照していいただけたらと思う．
全身の回転の運動方程式
モーメントは原点（回転の中心）に依存する物理量：
モーメントは回転の中心（ここでは原点と呼ぶ）に依存した物理量で，角運動量ベクトルもモーメントであることに注意されたい．
角運動量などのベクトルのモーメントの原点を移動する場合，平行軸の定理によって変換する．この平行軸の定理の理解が合成の角運動量の理解の肝である．平行軸の定理については
を参照されたい．
図１に，重心$${\bm{x}_g}$$まわりの角運動量$${\mathcal{L}(g)}$$とその時間微分$${\dot{\mathcal{L}}(g)}$$の，原点$${O}$$まわりの角運動量への変換を示した．角運動量は運動量ベクトル$${\bm{P}}$$のモーメントであるので，原点の移動のベクトル$${\bm{r}}$$分の運動量のモーメント$${\bm{r} \times \bm{P}}$$だけ，$${\mathcal{L}_(g)}$$に加算すればよいことを示している．
角運動量の微分（すなわちトルク）も同様である．運動量の時間微分のモーメント$${\bm{r} \times \dot{\bm{P}}}$$だけ加算すればよい．なお，運動量の時間微分$${\dot{\bm{P}}}$$は力に相当するので，これは力のモーメント（トルク）に相当する．
くどいが，角運動量ベクトルは運動量ベクトルのモーメントで，その微分に相当する力のモーメント（トルク）は，まさに力ベクトルに関するモーメントである．
ここで，この合成の角運動量の計算を多用するので，しっかりと理解していただけたらと思う．
全身の回転の運動方程式：
並進の運動方程式は全身の運動量の時間微分$${\dot{\bm{P}} = \bm{F}}$$で記述できる．ここで，$${\bm{F}}$$は床反力などの外力としての並進の力である．同様に，全身の回転の運動方程式は，合成の角運動量の時間微分$${\dot{\bm{L}}=\bm{N}}$$で記述される．ここで，回転の運動方程式の外力は身体と地面間に作用する力のモーメント$${\bm{N}}$$と，全身の重力$${M \bm{g}}$$のモーメントで構成される．
ここでは，フォースプレートを利用せず，身体のダイナミクスだけで全身運動における身体と地面間に作用する力のモーメント$${\bm{N}}$$を推定することを目的とする．これを利用すれば，全関節に作用するトルクもフォースプレートを利用せずに計算ができる．
ニュートン・オイラー法で関節に作用するトルクを計算するために，ここでは，地面と接する部位（ここで取り上げる例では右足）の重心$${\bm{x}_{gb}}$$まわりの力のモーメント$${\bm{N}(\bm{x}_{gb})}$$を計算するので，$${\bm{x}_{gb}}$$まわりの全身の回転の運動方程式を計算すると，
$$
\dot{\bm{L}}^O(\bm{x}_{gb}) = \bm{N}^O(\bm{x}_{gb}) + (\bm{x}_G^O - \bm{x}^O_{gb}) \times M\bm{g}^O
\\
= \bm{N}^O(\bm{x}_{gb}) + \bm{r}_G^O \times M\bm{g}^O
$$
となり，左辺の$${\dot{\bm{L}}^O(\bm{x}_{gb})}$$は，各部位$${j}$$の運動量の時間微分$${\dot{\bm{P}}_j}$$の$${\bm{x}_{gb}}$$まわりのモーメントの和と，各部位の重心まわりの角運動量の和に対する時間微分で記述され，
$$
\dot{\bm{L}}(\bm{x}_{gb}) =
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times \dot{\bm{P}}_{j}^O \right) + \frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j \mathcal{L}_j^j(g)
\\ =
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right) + \frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j)
$$
と書くことができる．ここで，
$$
\bm{r}_{G} \stackrel{\mathrm{def}}{=} \bm{x}_G - \bm{x}_{gb}
\\
\bm{r}_{gj} \stackrel{\mathrm{def}}{=} \bm{x}_{gj} - \bm{x}_{gb}
$$
のように，$${\bm{x}_{gb}}$$からみた各位置ベクトルを$${\bm{r}}$$で定義している．
運動方程式の左辺と右辺を合わせ，フォースプレートを利用せず，$${\bm{x}_{gb}}$$まわりの外力としての力のモーメント$${\bm{N}^O(\bm{x}_{gb})}$$は，最終的に
$$
\bm{N}^O(\bm{x}_{gb}) =
\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j) +
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right)
- \bm{r}_G \times M\bm{g}^O
$$
より計算される（補足１）．
近似式：
また，もし，この式を構成する
$$
\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j) \ll
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right)
$$
が成立する場合，つまり，
$$
\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j \mathcal{L}_j^j(g) \ll \sum_{j=1}^n \left(\bm{r}_{gj}^O \times \dot{\bm{P}}_{j}^O \right)
$$
が成立し，各部位の重心まわりの角運動量変化$${\dot{\mathcal{L}}_j(g) = \frac{d}{dt}(\bm{I}_j \bm{\omega}_j)}$$が無視できるほど小さい時は，近似された
$$
\hat{\bm{N}}^O(\bm{x}_{gb}) =
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right)
- \bm{r}_G \times M\bm{g}^O
$$
で計算してよい．なお，近似の目的は，計算負荷を減らすことと，特に$${\mathcal{L}_j(g)}$$が微分によってデータが暴れやすいことを防ぐ目的がある．また，それ以上に，その物理的意味が重要である．
運動方程式の原点の意味：
ここで，角運動量とその微分である力のモーメント（トルク）の「回転の中心（原点）」の意味を考えていきたい．
先程示したように，全身の角運動量をCOPまわりに考え，その時間部分を計算すると，「地面と身体間に作用する力のモーメント（鉛直軸まわりの摩擦モーメント）」となる．同じ全身でも身体重心まわりの角運動量の時間微分を考えれば，それは「身体重心まわりにまわすためのトルク」となる．一方，股関節以下の部位を系（対象）として考え，股関節まわりの角運動量の時間部分を考えるとそれは「股関節に作用するトルク」となる．
つまり，角運動量の原点は，回転の中心と当たり前の意味を持つが，それはその点まわりに「回す」ための運動量のモーメントであり，その微分はその点を通る軸回りに「回す」ための力のモーメントになる．適当に代表点として身体重心まわりの角運動量を考えればよいというわけではない．この原点の意味を考える必要がある．
解析・検証：片足飛び降りバランス動作
検証内容
ここでは，図４の片足での飛び降り動作を対象に，フォースプレートを使用せずに，外力としての床反力や力のモーメントを，モーションキャプチャのデータだけから推定することを試み，その精度を，下記に示した内容で検証する．
ただし，ここでの近似が，全身や一部の部位で成立することを，それぞれの動作でよく確認をする必要があるだろう．ここでの検証は，あくまでもこの動作で成立することを確かめている．
１．近似の確認：
以後の計算で，前述の近似を用いる場合があるが，この近似がこの飛び降り動作で成り立つか確認する．
そこで，回転の運動方程式の慣性力を構成する，各部位の運動量のモーメントの時間微分の加算と，各部位の角運動量の加算の時間微分を比較し，この動作で，各部位の角運動量の時間微分が小さく，
$$
\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j) \ll
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right)
$$
が成立することを検証する．
２．推定した外力のモーメントと，フォースプレートで計測したモーメントとの比較：
この動作で「１の近似」が成立することを確認した上で， b_link[1].External_Moment_approx()を使用し，$${\bm{x}_{gb}}$$まわりの外力のモーメント$${\hat{\bm{N}}^O(\bm{x}_{gb})}$$を計算し，フォースプレートで計測した力のモーメントと比較する．
なお，この動作では，コードのb_link[1].External_Moment_approx()で外力のモーメント$${\bm{N}}$$を計算するが，その場合$${\bm{I}_j^j \bm{\omega}_j^j}$$の項は計算されず，
$$
\hat{\bm{N}}^O(\bm{x}_{gb}) =
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right)
- \bm{r}_G \times M\bm{g}^O
$$
から計算することから，この力学計算では慣性モーメントや角速度の計算も不要となることに留意されたい．
３．関節に作用するトルクの比較：
関節に作用するトルクの計算では，外力としての床反力や力のモーメントが必要だが，フォースプレートを使用した計算と，２で推定した外力を使用して推定した場合を比較する．
４．ZMPとCOPとの比較：
第22章（前章）でも述べたように，フォースプレートの力と力のモーメントからCOPを計算する代わりに，身体のダイナミクスからやはりCOPを計算することができる．これをZMP（zero moment point）と呼び，ZMPの計算では２の外力は不要だが，身体のダイナミクスから推定するZMPを，フォースプレートから計算されるCOPと比較する．
検証１：近似の検証
$$
\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j) \ll
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right)
$$
が成り立ち，近似が可能となる場合，計算負荷を小さくすることが出きるが，そのことよりも，そのような運動や，それが成り立つ部位では，動力学的に動かさないことで並進の力（床反力など）を伝達するように制御しているという事実に注目すべきであろう．
次章で示すが，キック動作では接地している側の脚では近似が成立するが，おそらく蹴り足ではこの近似は成立しないだろう．この章の飛び降り動作では，着地後は全身でバランスを取るような運動なので，全身で近似が成立しているはずだ．このことをここから確認する．
以下のように，
dotL_mm = b_link[1].dot_Resultant_Momentum_Moment(bw, b_link[12].cg, sister=True) # b_link[12].cgまわりの運動量ベクトルの微分のモーメントを計算で各部位の運動量のモーメントの時間微分$${ \sum_{j=1}^n \left(\bm{r}_{gj}^O \times \dot{\boldsymbol{P}_j} \right) = \sum_j (\boldsymbol{x}_{gj} - \boldsymbol{x}_{gb}) \times m_j \ddot{\bm{x}}_{j}^O}$$を計算し，以下のようにグラフで示すと
t_list = np.arange(len(dotL_mm))/360.
plt.plot(t_list, dotL_mm)
plt.axvline(x = time_t, ls='--', c='.5')
plt.axvline(x = time_r , ls='--', c='.5')
plt.axhline(y=0., ls='--', c='.5')
plt.legend(['x', 'y', 'z'])
plt.xlabel('Time [s]')
plt.ylabel('Moment of Force [Nm]')
plt.ylim(-330, 250);
のようになる．時間軸の２つの破線は，離地時刻と着地時刻を示している．着地前も計算しているが，着地時間より前のデータは比較対象にならないので注意されたい．
次に，各部位の角運動量の時間微分$${\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j)}$$を
resultant_am = b_link[1].Resultant_Angular_Momentum(True)
dotL_am = smoothing_spline(resultant_am, 1, 360., 1)で計算し，以下のようにグラフで示すと
t_list = np.arange(len(dotL_am))/360.
plt.plot(t_list, dotL_am)
plt.axvline(x = time_t, ls='--', c='.5')
plt.axvline(x = time_r , ls='--', c='.5')
plt.axhline(y=0., ls='--', c='.5')
plt.legend(['x', 'y', 'z'])
plt.xlabel('Time [s]')
plt.ylabel('Moment of Force [Nm]')
plt.ylim(-330, 250);
のようになり，全身で近似が成立する
$$
\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j) \ll
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right)
$$
と考えてよいだろう．
検証２：推定した外力のモーメントと，フォースプレートで計測したモーメントとの比較
ここでは，フォースプレートで計測する外力の代わりに，身体のダイナミクスから外力を推定する．
床反力の比較：
床反力の推定は，第13章「全身の力学解析（補足）」
で示したが，ここでも再び，床反力を身体重心の加速度から推定し，以下のコードでフォースプレートのデータと比較すると，
# force plateの床反力データをmotion capture（絶対）座標に変換し，かつサンプリングレートもmotion captureの360Hzに変換
force_fp = fp2op_frame(file_path_fp_22_1, file_path_op_22_1, 2, 'f')
# 床反力をmotion caputreのデータから身体重心を推定し，その加速度から推定する
x_g = b_link[1].Resultant_Cg(bw, True) # 身体重心の位置ベクトル
force_op = bw * (smoothing_spline(x_g, .00001, 360, order=2) - gravity)
t_list = np.arange(len(force_op))/360.
plt.plot(t_list, force_fp)
[plt.plot(t_list, force_op[:,i], linestyle = '--', c=['b','orange','g'][i]) for i in range(3)]
plt.axvline(x = time_t, ls='--', c='.5')
plt.axvline(x = time_r , ls='--', c='.5')
plt.axhline(y = 0., ls='--', c='.5')
plt.legend(['Fx (FP)', 'Fy (FP)', 'Fz (FP)', 'Fx (estimated)', 'Fy (estimated)', 'Fz (estimated)'])
plt.xlabel('Time [s]')
plt.ylabel('Force [N]')
図７に示したように，精度良く床反力が推定できていることがわかる．
力のモーメントの比較：
ここでは，次の検証3の関節に作用するトルクの計算を意識し，右足の重心位置$${\bm{x}_{g12}}$$まわりの，外力としての力のモーメント$${\hat{\bm{N}}^O(\bm{x}_{g12})}$$をb_link[1].External_Moment_approx()で計算し，フォースプレートで計測する力のモーメントと比較する．なお，フォースプレートによって計測した力のモーメントも，フォースプレートの原点から$${\bm{x}_{g12}}$$まわりの力のモーメントに，平行軸の定理を利用し変換する．
# force plateの力のモーメントデータをmotion capture（絶対）座標に変換し，かつサンプリングレートもmotion captureの360Hzに変換
# 右足の重心（b_link[12].cg）まわりのモーメント
moment_12_fp = moment2origin(file_path_fp_22_1, file_path_op_22_1, b_link[12].cg, 12)
# 右足の重心まわりの，地面と右足間に作用する力のモーメントを，motion caputreのデータから推定
moment_12_op = b_link[1].External_Moment_approx(bw, b_link[12].cg, weight =.01, sf=360., sister=True)
t_list = np.arange(len(moment_12_op))/360.
plt.plot(t_list, moment_12_fp)
[plt.plot(t_list, moment_12_op[:,i], linestyle = '--', c=['b','orange','g'][i]) for i in range(3)]
plt.axvline(x = time_t, ls='--', c='.5')
plt.axvline(x = time_r , ls='--', c='.5')
plt.axhline(y=0., ls='--', c='.5')
plt.legend(['Mx (FP)', 'My (FP)', 'Mz (FP)', 'Nx (estimated)', 'Ny (estimated)', 'Nz (estimated)'])
plt.xlabel('Time [s]')
plt.ylabel('Moment [Nm]')
plt.xlim(2.6, 5.8)
plt.ylim(-170, 50);
図８では，主として着地後の$${\hat{\bm{N}}^O(\bm{x}_{g12})}$$を比較している．フォースプレートで計測した力のモーメントと比較すると，この試技の計測では，少し推定精度が低いようだ．
検証３：関節に作用するトルクでの比較
検証2で示したように，身体と地面間に作用する力のモーメント$${\bm{N}(\bm{x}_{g12})}$$の推定精度は低いが，一方で床反力$${\bm{F}}$$の推定精度は高い．また，関節に作用するトルクの計算では，床反力がつくる力のモーメントが支配的であることを考えると，関節に作用するトルクの推定精度は，ひどく影響を受けないかもしれない．このことを確認しよう．
ここで再び，以下のように外力として身体のダイナミクスから推定値を代入し，
# bw = 86.
# gravity = [0,0,-9.8]
b_link[12].ext_force = bw * (smoothing_spline(b_link[1].Resultant_Cg(bw, True), .00001, 360, order=2) - gravity)
b_link[12].ext_moment = b_link[1].External_Moment_approx(bw, b_link[12].cg, weight =.00001, sf=360., sister=True)ニュートン・オイラー法で関節に作用する力とトルクを，以下で計算する．
b_link[1].ID_Newton(86., True);
b_link[1].ID_Euler(True);その上で，右膝関節に作用するトルクを，以下のコードでグラフを描くと
tau_r_knee_approx = b_link[11].tau_global
t_list = np.arange(len(tau_r_knee_approx))/360.
plt.plot(t_list, tau_r_knee_approx)
plt.xlabel('Time [s]')
plt.ylabel('Torque [Nm]')
plt.legend(['x', 'y', 'z'])
plt.axvline(x = time_t, ls='--', c='.5')
plt.axvline(x = time_r , ls='--', c='.5')
plt.axhline(y=0., ls='--', c='.5')
plt.ylim(-400,150);
図９のようになる．離地と接地間の破線間の時間の推定値は不正確であるので無視していただきたい．
また，フォースプレートで計測して得られる床反力と力のモーメントは着地後だけなので，着地後の右膝関節に作用するトルクを，フォースプレートのデータで推定したトルク（実線）と，フォースプレートなしで推定したトルク（破線）の比較を図10に示す．
plt.plot(t_list, tau_r_knee)
[plt.plot(t_list, tau_r_knee_approx[:,i], linestyle = '--', c=['c','orange','g'][i]) for i in range(3)]
plt.xlabel('Time [s]')
plt.ylabel('Torque [Nm]')
plt.legend(['Tx', 'Ty', 'Tz', 'Tx (estimation)', 'Ty (estimation)', 'Tz (estimation)'])
plt.xlim(982/360., 2050/360.)
plt.ylim(-330,140)
plt.axhline(y=0., ls='--', c='.5');
着地直後に暴れるのは致し方ないが，この差に対する評価は微妙なところかもしれない．
検証４：ZMPとCOPとの比較
身体のダイナミクスから計算するCOPをロボティクスではZMPと呼ぶ．このZMPと，フォースプレートのデータで計算されるCOPを比較する．
この計算は，検証１〜３までの計算とは関係ない．
plt.plot(t_list, fp2op_frame(file_path_fp_22_1, file_path_op_22_1, 2, 'c')[:,0:2])
[plt.plot(t_list, estimated_cop(bw, b_link, .0001)[:,i], linestyle = '--', c=['c','orange'][i]) for i in range(2)]
plt.ylim(-.6, .6)
plt.xlim(982/360., 5.8)
plt.axvline(x = time_t, ls='--', c='.5')
plt.xlabel('Time [s]')
plt.ylabel('COP [m]')
plt.legend(['COPx', 'COPy', 'COPx (estimation)', 'COPy (estimation)']);
比較的変動が大きい前後（$${y}$$）方向でバランスを取っていることがわかる．図12に着地後のCOP（青）と推定したZMPの$${xy}$$（赤）位置をスティックピクチャに示した．数cmの誤差が発生している．
def stick_cg2(ax, body_link, cg, cop, zmp, r_force, time, lw=1, xrange=(-1.5, 1.), yrange=(-1.0, 1.5), zrange=(0, 2.5)):
ax.plot(*joint_p_list(body_link, time, upper_arm_list), linewidth=lw, c='C0')
ax.plot(*joint_p_list(body_link, time, body_list), linewidth=lw, c='C1')
ax.plot(*joint_p_list(body_link, time, lower_body_list), linewidth=lw, c='C2')
ax.scatter(*cg[time], s=50, c='k') # 重心
ax.scatter(*cop[time], s=25, c='b') # 右足のCOP（FPから計算したから計算したCOPは青色）
ax.scatter(*zmp[time], s=20, c='r') # 右足のZMP（推定値は赤色）
ax.plot([cg[time, 0], cop[time, 0]],
[cg[time, 1], cop[time, 1]],
[cg[time, 2], cop[time, 2]], c='k', linewidth=1, ls='--') # 重心に作用する力（＝床反力の合力）のベクトル
f_para = .0005
ax.plot([cop[time, 0], cop[time, 0] + f_para * r_force[time, 0]],
[cop[time, 1], cop[time, 1] + f_para * r_force[time, 1]],
[cop[time, 2], cop[time, 2] + f_para * r_force[time, 2]], c='b', linewidth=2) # 重心に作用する力（＝床反力の合力）のベクトル
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_xlim(*xrange)
ax.set_ylim(*yrange)
ax.set_zlim(*zrange)
ax.set_aspect('equal') # Colab.では機能しない# ３Dグラフ領域の作成
p_elev = 20
p_azim = -110
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(projection='3d')
ax.view_init(elev=p_elev, azim=p_azim)
f_len = 2000 # 終わりのフレーム
step_frame = 20 # アニメーションの際の描画でスキップするフレーム数 （20より小さくするとメモリ不足で計算が止まるかもしれない）
q = f_len // step_frame # フレーム数をstep_sizeで割った，整数の商
max_frame = q * step_frame # 描画する最大フレーム数
sf = 360. # サンプリングレート
cg = b_link[1].Resultant_Cg(86., True) # 重心位置
cop_fp = fp2op_frame(file_path_fp_22_1, file_path_op_22_1, 2, 'c') # Force Plateで計測したCOP
zmp = estimated_cop(bw, b_link, .0001) # ZMP
force_r_fp = fp2op_frame(file_path_fp_22_1, file_path_op_22_1, 2, 'f') # Force Plateで計測した右脚の床反力
# アニメーション更新用関数の作成
def plot_3d(frame):
ax.cla()
frame_new = step_frame * frame
stick_cg2(ax, b_link, cg, cop_fp, zmp, force_r_fp, frame_new, lw=2)
# アニメーションのオブジェクトを作成
anim = FuncAnimation(fig, plot_3d, frames = q, interval = 50)
# Jupyter内でアニメーションの表示
plt.close()
HTML(anim.to_jshtml())
### Gif/MP4に保存する場合 ###
#path_dir = '/Users/.../sample_data/'
#path = pathlib.Path(path_dir)
#path_mp4 = path.joinpath('animation_dropjum.gif')
#anim.save(path_mp4)
#plt.close()
おわりに
計測環境に依存するが，この飛び降りタスクでは，フォースプレートを利用せずとも床反力はモーションキャプチャだけから精度良く推定ができたが，力のモーメントの推定精度は低かった．十分な検証は行っていないが，おそらく力のモーメントの計算に置いて，モーメントアームのベクトルの方向と，運動量のベクトルの方向が近いことが原因と予想している．今後のさらなる検証が必要となる．
一方，関節のトルクの推定は，その精度の低い力のモーメントを使用している割には，決して良い精度とは言えないが，それらしい値を示している．
さて，そのような推定の精度も重要であるが，それよりも，この計算を通して，身体運動における合成の力学の物理的意味に注目してほしい．これらは，このマガジンの第14章以降で議論しているので，そちらも参照していただければと思う．
次章について
次章ではフォースプレートのデータがないので検証はできないが，アメリカンフットボールのキック動作を用いて，同様な計算を行う予定である．
補足
補足１：ニュートン・オイラー法との整合性
（訂正内容について，$${\bm{x}_{gb}}$$まわりのモーメントとして記述していなかった．お詫び申し上げる）
第22章の補足３の結果を利用すれば，
$$
\bm{N}^O(\bm{x}_{gb}) =
\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j) +
\sum_{j=1}^n \left(\bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{gj}^O \right)
- \bm{r}_G \times M\bm{g}^O
$$
の最後の項は
$$
\bm{r}_G \times M\bm{g}^O = \sum_{j=1}^O \left( \bm{x}_{gj} \times (m_j \bm{g})\right)
$$
と書き換えられ，これを利用し
$$
\bm{N}^O(\bm{x}_{gb}) =
\sum_{j=1}^n \left( \bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{j}^n \right) + \frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j)
- \sum_{j=1}^n \left( \bm{r}^O_{gj} \times m_j \bm{g}^O \right)
\\ =
\frac{d}{dt} \sum_{j=1}^n {}^O\bm{R}_j (\bm{I}_j^j \bm{\omega}_j^j) +
\sum_{j=1}^n \left( \bm{r}_{gj}^O \times m_j \ddot{\bm{x}}_{j}^O \right)
- \sum_{j=1}^n \left( \bm{r}_{gj}^O \times m_j \bm{g}^O \right)
\\ =
\sum_{j=1}^{n}{}^O \bm{R}_j \left(\bm{I}_j \dot{\bm{\omega}}_j + \bm{\omega} \times (\bm{I}_j \bm{\omega}_j) \right) +
\sum_{j=1}^n \left( \bm{r}_{gj}^O \times m_j (\ddot{\bm{x}}_{gj}^O ~– ~\bm{g}^O) \right)
$$
ここで，回転の慣性力項と，重心に作用する慣性力と重力加速度の和をそれぞれ
$$
\dot{\mathcal{L}}_i^O \equiv \bm{R}_j \left(\bm{I}_j \dot{\bm{\omega}}_j + \bm{\omega} \times (\bm{I}_j \bm{\omega}_j) \right)
\\
\mathcal{F}_i^O \equiv m_j (\ddot{\bm{x}}_{gj}^O ~–~ \bm{g}^O)
$$
と置けば，
$$
\bm{N}^O(\bm{x}_{gb})=
\sum_{j=1}^{n}\dot{\mathcal{L}}_i^O +
\sum_{j=1}^n \left( \bm{r}_{gj}^O \times m_j \mathcal{F}_i^O \right)
$$
のように，地面と身体間に作用する外力のモーメントを，慣性力$${\dot{\mathcal{L}}}$$と，関節$${\bm{x}_{i}}$$まわりの力のモーメント$${\bm{r}_{gj}^O \times m_j \mathcal{F}_i^O}$$の和に書き換えることができる．
合成の回転の運動方程式で，トルクを含めて力のモーメントの式は一見複雑そうだが，単純に，角運動量の時間変化$${\dot{\mathcal{L}}_i(g)}$$の和と，求めたい位置まわりのモーメントの和で記述しているだけだ．
このことを利用すれば，関節$${i}$$に作用するトルク$${\bm{\tau}_i}$$は，関節$${\bm{x}_i}$$まわりのモーメントで合成の運動方程式を記述すればよい．たとえば，図13の部位$${i}$$より遠位側の全部位（$${i}$$〜$${i+3}$$）を対象とし，部位$${i}$$の関節$${\bm{x}_i}$$まわりの合成の回転の運動方程式は，以下のように書くことができ，
$$
\bm{\tau}_{i}
= \sum_{i=1}^3 \dot{\mathcal{L}}_i(g) + \sum_{i=1}^3 ({}_i^{i}\bm{r} \times \mathcal{F}_i)- \bm{N}(\bm{x}_{g~i+3}) - {}^{i+3}_{i} \bm{r} \times \bm{F}
\\
= \sum_{i=1}^3 \dot{\mathcal{L}}_i(g) + \sum_{i=1}^3 \left( {}_i^{i}\bm{r} \times (m_i (\ddot{\bm{x}}_{gi} - \bm{g})) \right) - \bm{N}(\bm{x}_{g~i+3}) - {}^{i+3}_{i} \bm{r} \times \bm{F}
$$
関節関節$${\bm{x}_i}$$に作用するトルク$${\bm{\tau}_{i} }$$を解くことができる．ここで，$${{}_i^{i}\bm{r} \equiv \bm{x}_{gi} - \bm{x}_i}$$は関節$${\bm{x}_{i}}$$からみた各部位$${i}$$の重心$${\bm{x}_{gi}}$$の位置ベクトル，$${\bm{N}(\bm{x}_{g~i+3})}$$は末端の部位$${i+3}$$に作用する$${\bm{x}_{g~i+3}}$$まわりの外力としての力のモーメント，$${\bm{F}}$$は同じく部位$${i+3}$$に作用する外力としての並進力，$${\mathcal{F}_i \equiv m_i (\ddot{\bm{x}}_{gi} - \bm{g})}$$は部位$${i}$$の重心に作用する慣性力と重力の和である．また，最後の項は$${\bm{x}_{g~i+3}}$$まわりの外力としての力のモーメント$${\bm{N}(\bm{x}_{g~i+3})}$$を，平行軸の定理で関節$${\bm{x}_{i}}$$まわりに変換するための項である．
これは，第14章
で述べた，ニュートン・オイラー法で逐次的に計算する式の加算から算出したトルクの式と一致する．
このようにトルクや力のモーメントに限らず，多様に解くことでその意味をより深く理解できるので，ご自身でいろいろな方法で解き，ぜひ式の物理的意味を考えていただけたらと思う．
関節に作用するトルクも，単にニュートン・オイラー法のように逐次的に数値的に解かずに，この式から物理的意味を考え，トルクがどのように発揮されているか考えてみると良いだろう．少なくとも，ひとつの関節をひとつのトルクで動かしているわけではないことが理解できるだろう．また，さらに述べるならば，むしろ関節を動かさないことが重要な役目となる．
【著作権・転載・免責について】
権利の帰属
本ホームページで提示しているソフトウェアならびにプログラムリストは，スポーツセンシング社の著作物であり，スポーツセンシング社に知的所有権がありますが，自由にご利用いただいて構いません．
本ページに掲載されている記事，ソフトウェア，プログラムなどに関する著作権および工業所有権については，株式会社スポーツセンシングに帰属するものです．非営利目的で行う研究用途に限り，無償での使用を許可します．
転載
本ページの内容の転載については非営利目的に限り，本ページの引用であることを明記したうえで，自由に行えるものとします．
免責
本ページで掲載されている内容は，特定の条件下についての内容である場合があります． ソフトウェアやプログラム等，本ページの内容を参照して研究などを行う場合には，その点を十分に踏まえた上で，自己責任でご利用ください．また，本ページの掲載内容によって生じた一切の損害については，株式会社スポーツセンシングおよび著者はその責を負わないものとします．
【プログラムの内容について】
プログラムや内容に対する質問に対しては，回答できないことのほうが多くなると思いますが，コメントには目は通します．回答は必要最低限にとどめますので，返信はあまり期待しないでいただけると幸いです，
「動かして学ぶ」という大それたタイトルをつけたものの，また，きれいなプログラムに対するこだわりはあるものの，実際のプログラミングのスキルは決して高くありません．最下部の方のコメント欄によるプログラムの間違いのご指摘は歓迎します．できるだけ反映します．
【解析・受託開発について】
スポーツセンシングでは，豊富な知見を持つ，研究者や各種エンジニアが研究・開発のお手伝いをしております．研究・開発でお困りの方は，ぜひスポーツセンシングにご相談ください．
【例】
・データ解析の代行
・受託開発
（ハードウェア、組込みソフトウェア、PC/モバイルアプリ）
・測定システム構築に関するコンサルティング など
その他，幅広い分野をカバーしておりますので，まずはお気軽にお問い合わせください．
【データの計測について】
スポーツセンシング社のスタジオで，フォースプレートやモーションキャプチャを利用した計測も行えます．出力されるデータと，ここで示したプログラム（入力データの取り込み関数を少々改変する必要があるが）で，同様な解析を行えますので，まずはお気軽にお問い合わせください．