実際に使える制御工学：ロボット編（５）「ニュートン・オイラー法。使い勝手の良い運動解析法。」
＜この記事で伝えたい事＞
３次元の運動方程式はとても複雑です
３次元の運動方程式には、単純にロボット軸数が多いこと以上の複雑さがあります。
そんな３次元の運動方程式を解くのに便利な方法が、ニュートン・オイラー法です。ニュートン・オイラー法には多様な応用の利かせ方があります
うまく使えば、ロボット姿勢変化に応じた慣性モーメントや、その姿勢における重力トルクを割り出せたりします。
そうして抽出した情報は、ロボットを思いどおり動かすのに役立ちます。
1. はじめに
自己紹介はこちらです。
この記事は、主にロボット制御に関わる学生や若手～中堅エンジニアに向けた技術紹介・解説記事シリーズの第５弾です。
今回は、運動解析の一手法であるニュートン・オイラー法について紹介・解説します。
2. ３次元の運動方程式が難しい理由
2.1 並進運動と回転運動による複雑な加速度
並進運動と回転運動を同時にすることは、２次元運動（平面的な運動）でも普通にあります。
ここでの意味は、
「それを３次元に拡張すると複雑さが増して、理解が難しくなる」
くらいに思ってください。
下図は、ここで考える３次元運動のイメージです。
ロボットの実物には、座標原点$${O}$$とリンク質点や関節の位置との間には、いくつかリンクと関節が挟まっていると思ってください。
見たいリンク質点や関節の位置を３次元ベクトルで表したものが、位置ベクトル$${\bm{l}}$$です。
こうした見方をすると、ロボット運動に伴って位置ベクトル$${\bm{l}}$$が伸び縮みすると捉えられます。
また、同時に、関節回転を合わせたものが、その時々の回転軸まわりの回転運動とみなせます。
このとき、リンク質点や関節の速度ベクトル$${\bm{v}}$$は下式になります。
（$${\times}$$はベクトルの外積計算）
$$
\bm{v}=\bm{\dot{l}}+\bm{\omega}\times\bm{l}
$$
$${\bm{\dot{l}}}$$は、位置ベクトル$${\bm{l}}$$の伸び縮みする速度です。
$${\bm{\omega}\times\bm{l}}$$は、角速度によって生じる速度です。
（理屈はロボット編（３）の３章で解説済です）
シンプルに、これらの合計が、リンク質点や関節の速度$${\bm{v}}$$です。
さて、本題の加速度の話ですが、シンプルに
「速度$${\bm{v}}$$を微分すれば、加速度$${\bm{a}}$$が求められる」
とはなりません。
というのも、角速度$${\bm{\omega}}$$での円運動を維持するための求心加速度が発生している必要があるためです。
この求心加速度分を考慮すると、加速度$${\bm{a}}$$は以下のように求められます。
$$
\begin{array}{ll}
\bm{a}&=\cfrac{d}{dt}\bm{v}+\bm{\omega}\times\bm{v}\\
\\
&=\cfrac{d}{dt}(\bm{\dot{l}}+\bm{\omega}\times\bm{l})+\bm{\omega}\times(\bm{\dot{l}}+\bm{\omega}\times\bm{l})\\
\\
&=\bm{\ddot{l}}+\cfrac{d}{dt}(\bm{\omega}\times\bm{l})+\bm{\omega}\times\bm{\dot{l}}+\bm{\omega}\times(\bm{\omega}\times\bm{l})\\
\\
&=\bm{\ddot{l}}+\bm{\dot{\omega}}\times\bm{l}+\bm{\omega}\times\bm{\dot{l}}+\bm{\omega}\times\bm{\dot{l}}+\bm{\omega}\times(\bm{\omega}\times\bm{l})\\
\\
&=\bm{\ddot{l}}+\bm{\dot{\omega}}\times\bm{l}+2\bm{\omega}\times\bm{\dot{l}}+\bm{\omega}\times(\bm{\omega}\times\bm{l})
\end{array}
$$
なお、この数式の右辺の各項は、それぞれ以下の意味を持ちます。
第１項$${\bm{\ddot{l}}}$$：$${\bm{l}}$$の伸縮（並進）加速度
第２項$${\bm{\dot{\omega}}\times\bm{l}}$$：角加速度による加速度
第３項$${2\bm{\omega}\times\bm{\dot{l}}}$$：コリオリ力による加速度
第４項$${\bm{\omega}\times(\bm{\omega}\times\bm{l})}$$：求心加速度
コリオリ力とは、並進運動と回転運動を同時にするとき生じる（ように見える）力です。
この説明には、回転する円盤上でボールを投げる例えが、よく引き合いに出されます。
まず、回転する円盤の中心である$${O}$$点に人がいるとします。
この人は、$${A}$$点を狙ってボールを投げます。
しかし、円盤は回転しているので、ボールが円盤のふちに届くころには、元は$${A}$$点があった位置に$${B}$$点が来ます。
結果、ボールは$${B}$$点に到達します。
$${O}$$点にいる人は円盤ごと回転するので、この人からするとボールが曲がったように見えます。
ボールを曲げるのに必要な力が、コリオリ力です。
とはいえ、$${O}$$点にいる人からは曲がったように見えるだけで、実際には曲がらないので「見かけの力」と呼ばれます。
ただ、$${O}$$点はロボットでいうと、関節の回転中心に相当します。
関節にかかるトルクを計算する場合、この「見かけの力」も考慮しないと計算が合わないのが厄介なところです。
（前述の加速度計算なら、コリオリ力影響分が自然に反映されるので安心して使えます）
2.2 慣性テンソル
慣性テンソルは、３次元的な慣性モーメントとでもいうべきものです。
これも３次元運動方程式が難しい理由の１つです。
さしあたり、下図の直方体を$${x,y,z}$$軸まわりに回転させることを考えます。
直方体の質量を$${m}$$とすると、$${x,y,z}$$軸まわりの慣性モーメント$${I_{xx},I_{yy},I_{zz}}$$は以下になります。
$$
\begin{array}{ll}
I_{xx}&=\cfrac{m}{12}\cdot(w^2+h^2)\\
I_{yy}&=\cfrac{m}{12}\cdot(h^2+d^2)\\
I_{zz}&=\cfrac{m}{12}\cdot(d^2+w^2)\\
\end{array}
$$
同じ物体を回転させるにしても、どの軸まわりに回転させるかで慣性モーメントが違います。
これを行列形式でまとめたものが、慣性テンソル$${I}$$です。
$$
I=
\left[
\begin{array}{ccc}
I_{xx}&0&0\\
0&I_{yy}&0\\
0&0&I_{zz}\\
\end{array}
\right]
$$
直方体のように、形状が回転軸に対して対称なら、慣性テンソルの非対角成分はゼロです。
しかし、そうでなければ、慣性テンソルの非対角成分はゼロになりません。
慣性テンソルの非対角成分は「慣性乗積」と呼ばれます。
特に、６軸や７軸の垂直多関節ロボットでは、ロボット動作に伴って回転軸方向が変わります。
それを考慮して運動方程式を導出するのは、難しいことです。
2.3 ジャイロ効果
ジャイロ効果とは、回っているコマを安定して回り続けたり、走り出した自転車やオートバイの安定走行をサポートしてくれるアレです。
ジャイロ効果も慣性テンソルと同じく、３次元ならではのものです。
ジャイロ効果は、ロボット動作にも影響します。
なぜなら、ロボット関節も回転するからです。
（直動関節だけのロボットなら話は別ですが）
ジャイロ効果は日常生活でも体験できるので、「そういうものがある」のは説明不要と思います。
なので、計算方法だけ紹介します。
まず、下図のように$${x}$$軸まわりに角速度$${\omega_x}$$で回転している円盤を、$${y}$$軸まわりに角速度$${\omega_y}$$で回転させることを考えます。
$${y}$$軸まわりに$${90\degree}$$回転したとすると、円盤は$${z}$$軸まわりに角速度$$${\omega_x}$$で回転することになります。
$${y}$$軸まわりの回転によって、$${x}$$軸まわり回転が$${z}$$軸まわり回転に変わっていくわけですね。
このとき、角加速度は$${z}$$軸方向に生じていると解釈できます。
この角加速度は、円盤自身の角速度ベクトルと、円盤ごと回そうとする角速度ベクトルとの外積計算によって求められます。
$$
\begin{array}{ll}
\bm{\dot{\omega}}&=
\left[
\begin{array}{c}
0\\\omega_y\\0
\end{array}
\right]
\times
\left[
\begin{array}{c}
\omega_x\\0\\0
\end{array}
\right]
\\
\\
&=\left[
\begin{array}{c}
0\\0\\-\omega_x\cdot\omega_y
\end{array}
\right]
\end{array}
$$
計算結果からも、$${z}$$軸方向に角加速度が発生しているのが分かります。
ところが、この円盤には$${z}$$軸まわり回転の支えとなるものがありません。
よって、円盤は計算で求めた加速度と逆向きの$${z}$$方向回転を始める（$${x}$$軸まわり回転を保つ方向へ加速する）ことになります。
これがジャイロ効果として認識されます。
3. ニュートン・オイラー法
ニュートン・オイラー法とは、ニュートンの運動方程式（並進運動）とオイラーの運動方程式（回転運動）を組み合わせて、リンク構造を持つ機構の運動を解く方法です。
なお、リンクは剛体と仮定します。
軸数の多いロボットを扱うとき、ニュートン・オイラー法によって加速度・角加速度などから力・トルクを求めるケースが多いです。
（運動方程式の導出にも普通に使えますが）
ここでは、ニュートン・オイラー法の具体的な計算手順を説明します。
まず、本格的な計算の実行前にやっておきたい準備的な計算の説明です。
それは、関節座標基準での
位置ベクトル$${\bm{\hat{l}}}$$、
重心位置ベクトル$${\bm{\hat{l}_G}}$$、
回転軸ベクトル$${\bm{\hat{s}}}$$、
重心まわり慣性テンソル$${\hat{I}}$$、
を回転を考慮したものに変換するというものです。
（$${i}$$はリンク番号）
$$
\bm{\hat{l}_i}=
\left[
\begin{array}{c}
L_{xi}\\L_{yi}\\L_{zi}
\end{array}
\right]
$$
$$
\bm{\hat{l}_{Gi}}=
\left[
\begin{array}{c}
x_{Gi}\\y_{Gi}\\z_{Gi}
\end{array}
\right]
$$
$$
\begin{array}{c}
\bm{\hat{s}_i}=
\left[
\begin{array}{c}
0\\0\\1
\end{array}
\right]\\
\\
\text{（関節駆動方向がZ軸まわりかZ方向の場合）}
\end{array}
$$
$$
\hat{I}_{i}=
\left[
\begin{array}{ccc}
I_{xxi}&I_{xyi}&I_{xzi}\\
I_{yxi}&I_{yyi}&I_{yzi}\\
I_{zxi}&I_{zyi}&I_{zzi}
\end{array}
\right]
$$
以下、関節回転を考慮した絶対座標系での
位置ベクトル$${\bm{l}}$$、
重心位置ベクトル$${\bm{l_G}}$$、
回転軸ベクトル$${\bm{s}}$$、
重心まわり慣性テンソル$${I}$$、
です。
（$${^0R_i}$$はベースからそのリンクまでの回転行列）
$$
\begin{array}{ll}
\bm{l_i}&={^0R_{i+1}}\bm{\hat{l}_i}\\
\\
\bm{l_{Gi}}&={^0R_{i+1}}\bm{\hat{l}_{Gi}}\\
\\
\bm{s_i}&={^0R_{i+1}}\bm{\hat{s}_i}\\
\\
I_i&={^0R_{i+1}}\hat{I}_i{^0R_{i+1}}^T
\end{array}
$$
$$
^0R_{i+1}={^0R_1}{^1R_2}\cdots{^iR_{i+1}}
$$
この準備をしておくと、後の計算がかなり楽です。
ただ、上述のパラメータは、DHパラメータとの整合が取れたものになっていません。
パラメータ整合を取るのは、それほど難しくありません。
必要な方はやってみてはどうでしょうか？
（私も仕事ではきちんと整合を取っています）
それでは、本格的なニュートン・オイラー法をステップ分けで説明します。
（１１ステップ）
STEP1：各関節角速度の導出
ある関節からすると、それより根本側の関節が動くのは土台が動くのと同じことです。
よって、絶対座標系での関節角速度は、根本側のものが先端側に上乗せされていきます。
このため、基本的には手先側ほど角速度が速くなります。
$${\begin{array}{ll}\bm{\omega_0}&=\bm{0}\\\bm{\omega_1}&=\bm{\omega_0}+\bm{s_1}\dot{\theta_1}\\\bm{\omega_2}&=\bm{\omega_1}+\bm{s_2}\dot{\theta_2}\\&\vdots\\\bm{\omega_n}&=\bm{\omega_{n-1}}+\bm{s_n}\dot{\theta_n}\end{array}}$$
（$${n}$$はロボットの軸数、$${\bm{0}}$$は零ベクトル）
リンクは剛体という仮定のため、リンク重心と関節の角速度は一致します。
$${\bm{\omega_{Gi}}=\bm{\omega_i} (i=1\sim n)}$$
STEP2：各関節角加速度の導出
関節角加速度も、STEP1と同じ要領で求まります。
ただ、角加速度の場合、ジャイロ効果の影響も考慮する必要があります。
$${\begin{array}{ll}\bm{\dot{\omega}_0}&=\bm{0}\\\bm{\dot{\omega}_1}&=\bm{\dot{\omega}_0}+\bm{s_1}\ddot{\theta_1}+\bm{\omega_1}\times(\bm{s_1}\dot{\theta}_1)\\\bm{\dot{\omega}_2}&=\bm{\dot{\omega}_1}+\bm{s_2}\ddot{\theta_2}+\bm{\omega_2}\times(\bm{s_2}\dot{\theta}_2)\\&\vdots\\\bm{\dot{\omega}_n}&=\bm{\dot{\omega}_{n-1}}+\bm{s_n}\ddot{\theta_n}+\bm{\omega_n}\times(\bm{s_n}\dot{\theta}_n)\end{array}}$$
なお、角速度と同じく、リンク重心と関節の角加速度は一致します。
$${\bm{\dot{\omega}_{Gi}}=\bm{\dot{\omega}_i} (i=1\sim n)}$$
STEP3：各関節の速度の導出
２章の説明どおり、速度は下式で求まります。
$${\bm{v}=\bm{\dot{l}}+\bm{\omega}\times\bm{l}}$$
根本側の速度が先端側に上乗せされることを考慮すると、各関節の速度は以下のように求められます。
$${\begin{array}{ll}\bm{v_1}&=\bm{0}\\\bm{v_2}&=\bm{v_1}+\bm{\dot{l}_1}+\bm{\omega_1}\times\bm{l_1}\\&\vdots\\\bm{v_{n+1}}&=\bm{v_n}+\bm{\dot{l}_n}+\bm{\omega_n}\times\bm{l_n}\end{array}}$$
直動関節がない場合、$${\bm{\dot{l}_i}=\bm{0}}$$です。
$${n+1}$$番目関節はフランジ（またはアーム先端）を意味します。
STEP4：各関節の加速度の導出
２章の説明どおり、加速度は下式で求まります。
$${\bm{a}=\bm{\ddot{l}}+\bm{\dot{\omega}}\times\bm{l}+2\bm{\omega}\times\bm{\dot{l}}+\bm{\omega}\times(\bm{\omega}\times\bm{l})}$$
根本側の加速度が先端側に上乗せされることを考慮すると、各関節の速度は下式で求められます。
$${\begin{array}{ll}\bm{a_1}&=\bm{0}\\\bm{a_2}&=\bm{\ddot{l}_1}+\bm{\dot{\omega}_1}\times\bm{l_1}+2\bm{\omega_1}\times\bm{\dot{l}_1}+\bm{\omega_1}\times(\bm{\omega_1}\times\bm{l_1})\\&\vdots\\\bm{a_{n+1}}&=\bm{\ddot{l}_n}+\bm{\dot{\omega}_n}\times\bm{l_n}+2\bm{\omega_n}\times\bm{\dot{l}_n}+\bm{\omega_n}\times(\bm{\omega_n}\times\bm{l_n})\end{array}}$$
直動関節がない場合、$${\bm{\ddot{l}_i}=\bm{0},\bm{\dot{l}_i}=\bm{0}}$$です。
STEP5：各リンク重心点の速度の導出
重心点速度には、関節速度に対して関節角速度の分が上乗せされます。
このため、重心点速度は下式で求められます。
$${\bm{v_{Gi}}=\bm{v_i}+\bm{\omega_i}\times\bm{l_{Gi}} (i=1\sim n)}$$
STEP6：各リンク重心点の加速度の導出
重心点加速度もSTEP5と同じ要領で求まります。
ただ、加速度計算式が２章で説明したものなことだけ要注意です。
$${\bm{a_{Gi}}=\bm{a_i}+\bm{\dot{\omega}_i}\times\bm{l_{Gi}}+\bm{\omega_i}\times(\bm{\omega_i}\times\bm{l_{Gi}}) (i=1\sim n)}$$
リンクは剛体という仮定のため、$${\bm{\ddot{l}_i}=\bm{0},\bm{\dot{l}_i}=\bm{0}}$$です。
STEP7：重心点にかかる力の導出
重力加速度分を考慮して、以下のように力を求められます。
（$${m_i}$$はリンク質量）
$${\bm{f_{Gi}}+m_i\bm{g}=m_i\bm{a_{Gi}}}$$
$${\therefore\bm{f_{Gi}}=m_i(\bm{a_{Gi}}-\bm{g})}$$
$${\bm{g}}$$は重力ベクトルであり、下式になります。
（$${g}$$は重力加速度、太字か細字かでベクトルかを判別します）
$${\bm{g}=\left[\begin{array}{c}0\\0\\-g\end{array}\right]}$$
STEP8：重心点にかかるモーメントの導出
ジャイロ効果分を考慮して、以下のようにモーメントを求められます。
$${\bm{n_{Gi}}=I_i\bm{\dot{\omega}_i}+\bm{\omega_i}\times(I_i\bm{\omega_i})}$$
ジャイロ効果によるモーメントの計算には、ややひねりが効いています。
一応、注意ください。
STEP9：関節にかかる力の導出
先端側にかかる力は、根本側に伝わります。
このため、基本的には根本側ほど大きな力がかかります。
それに伴い、STEP8までは根本側から計算していたのを、手先側から計算していくのに切り替える必要があります。
$${\begin{array}{ll}\bm{f_{n+1}}&=\bm{F}\\\bm{f_n}&=\bm{f_{n+1}}+\bm{f_{Gn}}\\\bm{f_{n-1}}&=\bm{f_n}+\bm{f_{Gn-1}}\\&\vdots\\\bm{f_i}&=\bm{f_{i+1}}+\bm{f_{Gi}} (i=1\sim n)\end{array}}$$
$${\bm{F}}$$は外力に対する反作用力ベクトルです。
下式で表せます。
$${\bm{F}=\left[\begin{array}{c}F_x\\F_y\\F_z\end{array}\right]}$$
STEP10：関節にかかるモーメントの導出
力と同じく先端側にかかるモーメントは、根本側に伝わります。
このため、基本的に根本側ほど大きなモーメントがかかります。
やはり、モーメントも手先側から計算します。
$${\begin{array}{ll}\bm{n_{n+1}}&=\bm{N}\\\bm{n_n}&=\bm{n_{n+1}}+\bm{l_n}\times\bm{f_{n+1}}+\bm{l_{Gn}}\times\bm{f_{Gn}}+\bm{n_{Gn}}\\\bm{n_{n-1}}&=\bm{n_n}+\bm{l_{n-1}}\times\bm{f_n}+\bm{l_{Gn-1}}\times\bm{f_{Gn-1}}+\bm{n_{Gn-1}}\\&\vdots\\\bm{n_i}&=\bm{n_{i+1}}+\bm{l_i}\times\bm{f_{i+1}}+\bm{l_{Gi}}\times\bm{f_{Gi}}+\bm{n_{Gi}} (i=1\sim n)\end{array}}$$
$${\bm{N}}$$は外力モーメントに対する反作用モーメントベクトルです。
下式で表せます。
$${\bm{N}=\left[\begin{array}{c}N_x\\N_y\\N_z\end{array}\right]}$$
STEP11：アクチュエータにかかるトルクの導出
関節座標基準でみて、$${Z}$$軸まわりモーメントがアクチュエータにかかるトルクです。
ただ、STEP10で導出したモーメントは絶対座標系のものなので、トルクを求めるには以下の計算式を使う必要があります。
$${\tau_i={\bm{s_i}}^T\bm{n_i} (i=1\sim n)}$$
STEP11(2)：アクチュエータにかかる力の導出
関節座標基準でみて、$${Z}$$軸方向の力がアクチュエータにかかる力です。
スカラロボットのように直動関節を持つロボットでは計算が必要です。
トルク計算と同じ要領で力を求められます。
$${F_i={\bm{s_i}}^T\bm{f_i} (i=1\sim n)}$$
全軸のトルクをまとめたトルクベクトル$${\bm{\tau}}$$は、以下の形で表せます。
（直動関節を持たないロボットの場合）
$$
\begin{array}{ll}
\bm{\tau}&=
\left[
\begin{array}{c}
\tau_1\\\vdots\\\tau_n
\end{array}
\right]\\
\\
&=M_{(\theta)}
\left[
\begin{array}{c}
\ddot{\theta}_1\\\vdots\\\ddot{\theta}_n
\end{array}
\right]
+C_{(\theta,\dot{\theta})}+G_{(\theta)}+{J_{(\theta)}}^T
\left[
\begin{array}{c}
\bm{F}\\\bm{N}
\end{array}
\right]
\end{array}
$$
直動関節を持つロボットの場合、角速度ベクトルの一部が$${\ddot{l}}$$に置き換わるなどの変化があります。
$${M_{(\theta)}}$$は慣性行列と呼ばれる$${n\times n}$$行列で、必ず正則になります。
（つまり、必ず逆行列が計算可能）
特に、その対角成分を主慣性（自軸自身へ与える影響成分）、非対角成分を相互慣性（自軸が他軸に与える影響成分）と呼びます。
また、
$${M_{(\theta)}\left[\begin{array}{ccc}\ddot{\theta}_1&\cdots&\ddot{\theta}_n\end{array}\right]^T}$$は慣性項、
$${C_{(\theta,\dot{\theta})}}$$は遠心・コリオリ力項、
$${G_{(\theta)}}$$は重力項、
と呼ばれます。
$${{J_{(\theta)}}^T}$$はヤコビアンの転置行列、$${\left[\begin{array}{cc}F^T&N^T\end{array}\right]^T}$$は外力・外力モーメントベクトルです。
したがって、$${{J_{(\theta)}}^T\left[\begin{array}{cc}F^T&N^T\end{array}\right]^T}$$は仮想仕事の原理によって、関節トルクに変換した外力です。
4. ニュートン・オイラー法の応用例
前述のとおり、ニュートン・オイラー法は、リンクパラメータと関節角度・角速度・角加速度を与えて、それに応じたトルクを計算する逆動力学として使えます。
逆動力学では、リンクパラメータや関節角度・角速度・角加速度が自由に設定できる入力で、トルクがその条件における出力です。
ここでは、角速度・角加速度の与え方の工夫で、運動方程式中に含まれる情報を抽出する方法を紹介・解説します。
4.1 慣性行列の抽出
角速度がゼロのとき、遠心・コリオリ力はゼロになります。
また、重力加速度がゼロのとき、重力はゼロです。
そこで、角速度・重力加速度を全てゼロに設定することを考えます。
角度は現在値を与えます。
このとき、逆動力学計算を実行すると、出力として得られるトルクは、以下のように現在角度に対応した慣性項のみになります。
（$${n}$$はリンク数、外力はゼロという前提）
$$
\begin{array}{ll}
\left[
\begin{array}{c}
\tau_1\\\vdots\\\tau_n
\end{array}
\right]&=
M_{(\theta)}
\left[
\begin{array}{c}
\ddot{\theta}_1\\\vdots\\\ddot{\theta}_n
\end{array}
\right]\\
\\
&=
\left[
\begin{array}{ccc}
M_{11}&\cdots&M_{1n}\\
\vdots&\ddots&\vdots\\
M_{n1}&\cdots&M_{nn}
\end{array}
\right]
\left[
\begin{array}{c}
\ddot{\theta}_1\\\vdots\\\ddot{\theta}_n
\end{array}
\right]
\end{array}
$$
さらに、角加速度をどれか１つだけ「１」に設定し、その他をゼロに設定することを考えます。
例えば、$${\ddot{\theta}_1=1}$$にして他の角加速度はゼロに設定すると、逆動力学の出力には１軸めに関する慣性情報だけが現れます。
$$
\begin{array}{ll}
\left[
\begin{array}{c}
\tau_1\\\tau_2\\\vdots\\\tau_n
\end{array}
\right]&=
\left[
\begin{array}{cccc}
M_{11}&M_{12}&\cdots&M_{1n}\\
M_{21}&M_{22}&\cdots&M_{2n}\\
\vdots&\vdots&\ddots&\vdots\\
M_{n1}&M_{n2}&\cdots&M_{nn}
\end{array}
\right]
\left[
\begin{array}{c}
1\\0\\\vdots\\0
\end{array}
\right]\\
\\
&=
\left[
\begin{array}{c}
M_{11}\\M_{21}\\\vdots\\M_{n1}
\end{array}
\right]
\end{array}
$$
同じ要領で、角加速度をどれか１つだけ「１」に設定し、他はゼロに設定して逆動力学計算を実行すると、「１」に設定した軸に関する慣性情報が得られます。
これの繰り返しによって、慣性行列の全要素が得られます。
慣性行列は、制御や軌道生成に活用できます。
4.2 遠心・コリオリ力項の抽出
角加速度がゼロのとき、慣性項はゼロです。
また、重力加速度がゼロのとき、重力はゼロです。
よって、全ての角加速度と重力加速度をゼロに設定し、角速度を現在値に設定して逆動力学計算を実行したとき、出力には遠心・コリオリ力項のみが現れます。
やはり、角度は現在値を与えます。
（外力はゼロという前提）
$$
\left[
\begin{array}{c}
\tau_1\\ \vdots\\ \tau_n
\end{array}
\right]
=C_{(\theta,\dot{\theta})}
$$
慣性行列に比べると重要度は低めですが、これも制御や軌道生成に活用できる情報です。
4.3 重力項の抽出
角加速度がゼロのとき、慣性項はゼロです。
角速度がゼロのとき、遠心・コリオリ力はゼロになります。
よって、全ての角加速度と角速度をゼロに設定し、重力加速度をゼロにせず逆動力学計算を実行すると、出力には重力項のみが現れます。
やはり、角度は現在値を与えます。
（外力はゼロという前提）
$$
\left[
\begin{array}{c}
\tau_1\\ \vdots\\ \tau_n
\end{array}
\right]
=G_{(\theta)}
$$
大型機種ほど、かかる重力が大きくなります。
（質量が大きいので）
特に、軌道生成をするとき、重力がいくつかは押さえておきたい情報です。
5. おわりに
この記事では、ニュートン・オイラー法について紹介・解説しました。
ラグランジュ法でも運動解析はできますが、計算量の面ではニュートン・オイラー法が有利です。
（そういう技術文献を読んだことがあります）
これまでの記事でも何度か言ったと思いますが、人間にとって大変な計算は、コンピュータにとっても大変です。
特に、計算リソースの制約が厳しい工業製品（ロボットコントローラ）では、ニュートン・オイラー法を採用するケースが多いです。
また、応用編（６）で紹介・解説したロボット関節のモデリング方法と、今回のニュートン・オイラー法による運動解析とを組み合わせれば、関節剛性を考慮したロボット運動解析が可能です。
（海外では、Spong-Vidyasagarの柔軟関節モデルとして広く知られているそうです）
次回は、「インピーダンス制御」を紹介・解説しようと思います。
前回：ロボット編（４）「逆キネマティクス」はこちら
次回：ロボット編（６）「インピーダンス制御」はこちら
