国際宇宙ステーション（ISS）のロボットアーム「Canadarm2」が、数トンの貨物モジュールを把持してゆっくりと移動させている映像を見たことがあるでしょうか。地上の工場ロボットなら、アームがどれだけ重い荷物を持ち上げても、ロボットの台座（ベース）はコンクリートの床にボルトで固定されているため微動だにしません。しかし宇宙では事情が根本的に異なります。アームが右に動けば、衛星本体は左に動く — これはニュートンの第三法則（作用・反作用の法則）の直接的な帰結です。
この「ベースが固定されていない」という一見単純な事実が、宇宙ロボティクスの動力学を地上ロボティクスとは本質的に異なるものにしています。地上ロボットで当たり前のように使える運動学・動力学の公式がそのままでは成り立たず、系全体の運動量保存則に立脚した新しい枠組みが必要になるのです。
この概念を理解すると、以下のような応用が可能になります。
- 軌道上サービスミッション: 故障衛星の捕獲・修理で、マニピュレータ操作中に母衛星の姿勢が乱れないよう制御する設計
- 宇宙デブリ除去: デブリを掴んだ際に、系全体の角運動量保存から生じるタンブリング（回転）を予測し、安定した捕獲戦略を立案
- 月面・火星探査: 低重力環境でのローバーマニピュレータ設計で、ベースへの反力を正しく見積もる
- 軌道上組立: 大型構造物を宇宙で組み立てる際に、連成系としてのロボット動力学を把握する
本記事の内容
- 地上ロボットと宇宙ロボットの本質的な違い
- 自由浮遊（free-floating）と自由飛行（free-flying）の分類
- 運動量保存則の導出とベース衛星の運動への帰結
- 動的カップリングの物理的メカニズム
- 自由浮遊マニピュレータの運動方程式と有効質量行列
- Pythonによるシミュレーションとベース反動の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ラグランジュ法によるマニピュレータ動力学 — ラグランジアンからマニピュレータの運動方程式を導出する方法
- 慣性行列・コリオリ項・重力項 — マニピュレータ動力学の各項の物理的意味
地上ロボットと宇宙ロボットの本質的な違い
固定ベースの前提
工場で稼働する産業用ロボットを思い浮かべましょう。ロボットのベース（台座）は工場の床にボルトで固定されています。アームがどれだけ高速に動いても、どれだけ重い荷物を持ち上げても、ベースの位置と姿勢は変わりません。この「固定ベース」の前提があるからこそ、地上ロボティクスの教科書で扱う運動学・動力学は比較的シンプルに定式化できます。
具体的に言えば、固定ベースのロボットでは次のことが成り立ちます。
- 順運動学: 関節角度 $\bm{q} = (q_1, q_2, \dots, q_n)^\top$ を与えれば、ベースの位置・姿勢が既知なので、手先位置 $\bm{r}_\mathrm{ee}$ が一意に決まる
- 逆運動学: 手先の目標位置から関節角度を求めれば、それだけでロボット全体の姿勢が確定する
- 動力学: 関節トルク $\bm{\tau}$ がアームの運動だけに影響し、ベースは考えなくてよい
宇宙での状況変化
ところが宇宙空間では、ロボットのベースである衛星本体（サービサ衛星やISSなど）は宇宙空間に浮いています。床も壁もありません。スラスタを噴射すれば外力を加えられますが、それ以外の状態では外力はゼロです（微小重力環境を仮定）。
このとき、アームが動くと作用・反作用の法則により、ベース衛星に反力が生じます。アームが時計回りに回転すれば、衛星本体は反時計回りに回転する — これは角運動量保存の直接的な結果です。
この事実は、地上ロボティクスの全ての前提を覆します。
- 順運動学が変わる: 関節角度 $\bm{q}$ を与えただけでは手先位置が決まらない。ベースの位置・姿勢 $(\bm{r}_0, \bm{\Theta}_0)$ も変数になるため、手先位置はベース状態と関節角度の両方の関数 $\bm{r}_\mathrm{ee} = f(\bm{r}_0, \bm{\Theta}_0, \bm{q})$ になる
- 逆運動学が変わる: 手先を目標位置に持っていく関節軌道が、ベースの動きに依存するため、系全体を連成系として解く必要がある
- 動力学が変わる: 関節トルク $\bm{\tau}$ はアームだけでなくベース衛星の運動にも影響する。ベース衛星の運動方程式とアームの運動方程式が連成する
衝撃的な帰結 — 手先が思い通りに動かない
地上ロボットでは「関節1を30度回す」と命令すれば、手先は幾何学的に確定した位置に移動します。しかし宇宙ロボットでは、同じ「関節1を30度回す」という命令に対して、ベース衛星が反対方向に回転するため、手先の到達位置が地上の場合とは異なります。
イメージとしては、凍った湖の上に立って重い荷物を前に押し出す状況に似ています。荷物は前に進みますが、あなた自身は後ろに滑ります。足元が固定されていないため、荷物への力の反作用があなたの運動に直接跳ね返るのです。
この「ベースの反動」を正しく理解し、予測し、制御するためには、系全体の運動量保存則に基づく動力学が必要です。次のセクションでは、まず宇宙ロボットの分類を整理してから、その核心である運動量保存則の定式化に進みます。
自由浮遊と自由飛行 — 宇宙ロボットの2つのモード
宇宙ロボットは、ベース衛星のスラスタやリアクションホイールの使用状態に応じて、大きく2つの動作モードに分類されます。この分類は、適用する動力学モデルを選択する上で極めて重要です。
自由浮遊（Free-Floating）モード
自由浮遊モードとは、ベース衛星のスラスタやリアクションホイールなどのアクチュエータを全て停止した状態で、マニピュレータの関節トルクのみで作業を行うモードです。
このモードでは、系に外力・外トルクが作用しないため、線運動量と角運動量が厳密に保存されます。
$$ \bm{p} = \text{const}, \quad \bm{L} = \text{const} $$
自由浮遊モードは、以下の点で実用上極めて重要です。
- 推進剤を消費しない: スラスタを使わないため、燃料を節約できる。長期ミッションでは死活的に重要
- プルーム汚染がない: スラスタ噴射による排気ガスが対象衛星や光学機器を汚染しない
- 外乱を生じない: リアクションホイールの振動や磁気トルカの磁場が搭載機器に影響しない
その代償として、マニピュレータの動作によるベース衛星の反動を受け入れなければなりません。また、関節角度の組み合わせによっては手先が到達できない「動的特異点」が現れるという、地上ロボットにはない問題も生じます（これは後続記事で詳しく扱います）。
自由飛行（Free-Flying）モード
自由飛行モードとは、ベース衛星のスラスタやリアクションホイールを積極的に使用して、ベースの位置・姿勢を制御しながらマニピュレータを操作するモードです。
このモードでは、スラスタやリアクションホイールが外力・外トルクを系に加えるため、運動量は保存されません。ベースの位置・姿勢を一定に保てるなら、地上ロボットと同じ運動学・動力学がほぼそのまま適用できます。
しかし、自由飛行モードには以下の欠点があります。
- 推進剤を消費する: マニピュレータ操作のたびに燃料を使う
- プルーム汚染: 噴射ガスが対象物や搭載センサを汚染するリスク
- 帯域幅の競合: ベース姿勢制御とマニピュレータ制御が同時に行われるため、制御系の設計が複雑
2つのモードの関係
実際のミッションでは、両方のモードを切り替えて使うことが多いです。例えば、対象物に接近するフェーズでは自由飛行モードで正確な位置合わせを行い、接触・把持のフェーズでは自由浮遊モードに切り替えてプルーム汚染を避ける — という戦略がとられます。
本記事では、宇宙ロボティクスの最も本質的な特徴を浮き彫りにする自由浮遊モードに焦点を当てます。運動量保存という強い拘束条件のもとで、系の動力学がどのように記述されるかを見ていきましょう。
系のモデル化 — ベース衛星 + n リンクアーム
座標系と変数の定義
自由浮遊ロボットの動力学を定式化するために、まず系のモデルを明確に定義します。ベース衛星と $n$ 自由度のマニピュレータアームからなる系を考えます。
慣性座標系 $\Sigma_I$: 宇宙空間に固定された座標系（近似的に慣性系とみなす）。原点は系の初期重心位置にとると便利です。
ベース座標系 $\Sigma_0$: ベース衛星の重心に原点を持ち、衛星本体に固定された座標系。
リンク座標系 $\Sigma_i$ ($i = 1, 2, \dots, n$): 各リンクの重心に原点を持つ座標系。
系全体の一般化座標を次のように定義します。
$$ \bm{x} = \begin{pmatrix} \bm{r}_0 \\ \bm{\Theta}_0 \\ \bm{q} \end{pmatrix} \in \mathbb{R}^{6+n} $$
ここで各変数の意味は以下の通りです。
- $\bm{r}_0 \in \mathbb{R}^3$: 慣性座標系におけるベース衛星の重心位置
- $\bm{\Theta}_0 \in \mathbb{R}^3$: ベース衛星の姿勢角（例えばオイラー角やロール・ピッチ・ヨー）
- $\bm{q} = (q_1, q_2, \dots, q_n)^\top \in \mathbb{R}^n$: マニピュレータの関節角度
地上の固定ベースロボットの一般化座標は $\bm{q}$ の $n$ 個だけですが、宇宙ロボットでは $\bm{r}_0$ と $\bm{\Theta}_0$ の6個が加わり、合計 $6 + n$ 個の一般化座標を持つことに注意してください。
各ボディの物理量
ベース衛星（ボディ0）と各リンク（ボディ $i$, $i = 1, \dots, n$）の物理量を定義します。
| 物理量 | ベース衛星 | リンク $i$ |
|---|---|---|
| 質量 | $m_0$ | $m_i$ |
| 重心位置（慣性系） | $\bm{r}_0$ | $\bm{r}_i$ |
| 重心速度（慣性系） | $\bm{v}_0 = \dot{\bm{r}}_0$ | $\bm{v}_i = \dot{\bm{r}}_i$ |
| 角速度 | $\bm{\omega}_0$ | $\bm{\omega}_i$ |
| 重心まわりの慣性テンソル | $\bm{I}_0$ | $\bm{I}_i$ |
系全体の質量を $M$ とすると、
$$ M = m_0 + \sum_{i=1}^{n} m_i $$
です。系全体の重心位置 $\bm{r}_G$ は、
$$ \bm{r}_G = \frac{1}{M}\left(m_0 \bm{r}_0 + \sum_{i=1}^{n} m_i \bm{r}_i\right) $$
で定義されます。
各リンクの位置 $\bm{r}_i$ と速度 $\bm{v}_i$、角速度 $\bm{\omega}_i$ は、ベースの状態 $(\bm{r}_0, \bm{\Theta}_0, \dot{\bm{r}}_0, \bm{\omega}_0)$ と関節変数 $(\bm{q}, \dot{\bm{q}})$ の関数です。具体的には、運動学の連鎖公式を使って次のように書けます。
$$ \bm{v}_i = \bm{J}_{v_i}^{(b)} \bm{v}_b + \bm{J}_{v_i}^{(m)} \dot{\bm{q}} $$
$$ \bm{\omega}_i = \bm{J}_{\omega_i}^{(b)} \bm{\omega}_b + \bm{J}_{\omega_i}^{(m)} \dot{\bm{q}} $$
ここで $\bm{v}_b = \dot{\bm{r}}_0$, $\bm{\omega}_b = \bm{\omega}_0$ はベースの並進・回転速度であり、$\bm{J}^{(b)}$ と $\bm{J}^{(m)}$ はそれぞれベース運動とマニピュレータ運動に対応するヤコビ行列です。これらのヤコビ行列は $\bm{q}$ に依存する構成依存の行列です。
この表現は、「各リンクの速度は、ベースの速度に起因する成分と、マニピュレータの関節速度に起因する成分の重ね合わせである」ことを意味しています。地上ロボットでは $\bm{v}_b = \bm{0}$, $\bm{\omega}_b = \bm{0}$ なので、ベースに関する項が消えて $\bm{v}_i = \bm{J}_{v_i}^{(m)} \dot{\bm{q}}$ となり、馴染みのある形に戻ります。
系のモデルが定義できたので、次にこの系に成り立つ最も基本的な法則 — 運動量保存則を導出しましょう。
運動量保存則の導出
線運動量の保存
自由浮遊モードでは、系に外力が作用しません。ニュートンの第二法則から、系全体の線運動量 $\bm{p}$ の時間変化はゼロです。
$$ \bm{p} = m_0 \bm{v}_0 + \sum_{i=1}^{n} m_i \bm{v}_i = \text{const} $$
ここで重心の速度を思い出すと、線運動量は次のようにも書けます。
$$ \bm{p} = M \dot{\bm{r}}_G $$
つまり、系の重心は等速直線運動をします。初期条件で $\dot{\bm{r}}_G = \bm{0}$ とすれば（重心位置を原点に固定した座標系を使えば）、
$$ \bm{r}_G = \text{const} $$
となり、系の重心位置は不変です。これは非常に直感的な結果です。宇宙空間でアームを動かしても、衛星とアーム全体の「平均的な位置」は一切変わらないのです。
線運動量保存を成分表示で書くと、
$$ m_0 \bm{v}_0 + \sum_{i=1}^{n} m_i \bm{v}_i = \bm{p}_0 $$
ここで $\bm{p}_0$ は初期線運動量です。ベース速度 $\bm{v}_0$ について解くと、
$$ \bm{v}_0 = \frac{1}{m_0}\left(\bm{p}_0 – \sum_{i=1}^{n} m_i \bm{v}_i\right) $$
が得られます。右辺の $\bm{v}_i$ は関節速度 $\dot{\bm{q}}$ とベース速度の関数なので、この式はベースとアームの運動の連成関係を表しています。
角運動量の保存
外トルクもゼロなので、系全体の角運動量 $\bm{L}$ も保存されます。系の重心まわりの角運動量は、
$$ \bm{L} = \bm{I}_0 \bm{\omega}_0 + m_0 (\bm{r}_0 – \bm{r}_G) \times \bm{v}_0 + \sum_{i=1}^{n} \left[\bm{I}_i \bm{\omega}_i + m_i (\bm{r}_i – \bm{r}_G) \times \bm{v}_i\right] = \text{const} $$
各項の物理的意味を確認しておきましょう。
- $\bm{I}_0 \bm{\omega}_0$: ベース衛星自身のスピン角運動量（自転による寄与）
- $m_0 (\bm{r}_0 – \bm{r}_G) \times \bm{v}_0$: ベース衛星の軌道角運動量（重心まわりの公転による寄与）
- $\bm{I}_i \bm{\omega}_i$: 各リンクのスピン角運動量
- $m_i (\bm{r}_i – \bm{r}_G) \times \bm{v}_i$: 各リンクの軌道角運動量
初期角運動量を $\bm{L}_0$ とすると、$\bm{L} = \bm{L}_0 = \text{const}$ です。多くの解析では、初期状態で系が静止していると仮定して $\bm{L}_0 = \bm{0}$ とします。
運動量保存の行列表現
前セクションで導入したヤコビ行列を使って、各ボディの速度をベース速度と関節速度で表すことができました。これを運動量の式に代入すると、運動量保存則は次のコンパクトな行列形式にまとめられます。
まず、ベースの一般化速度を $\dot{\bm{x}}_b = (\bm{v}_0^\top, \bm{\omega}_0^\top)^\top \in \mathbb{R}^6$ とまとめます。すると、線運動量と角運動量をまとめた一般化運動量ベクトル $\bm{P} = (\bm{p}^\top, \bm{L}^\top)^\top \in \mathbb{R}^6$ は、
$$ \bm{P} = \bm{H}_b \dot{\bm{x}}_b + \bm{H}_{bm} \dot{\bm{q}} = \bm{P}_0 $$
と書けます。ここで、
- $\bm{H}_b \in \mathbb{R}^{6 \times 6}$: ベースの一般化慣性行列
- $\bm{H}_{bm} \in \mathbb{R}^{6 \times n}$: ベースとマニピュレータの連成慣性行列（カップリング慣性行列）
- $\bm{P}_0$: 初期一般化運動量
$\bm{H}_b$ の具体的な形を書き下すと、次のようになります。
$$ \bm{H}_b = \begin{pmatrix} M \bm{E}_3 & -M \tilde{\bm{d}}_G \\ M \tilde{\bm{d}}_G & \bm{I}_\mathrm{total} \end{pmatrix} $$
ここで $\bm{E}_3$ は $3 \times 3$ の単位行列、$\tilde{\bm{d}}_G$ はベース重心から系全体の重心までの位置ベクトルの外積行列（歪対称行列）、$\bm{I}_\mathrm{total}$ は系全体の重心まわりの慣性テンソルです。
$\bm{H}_{bm}$ は連成慣性行列で、マニピュレータの関節速度がベースの運動量にどれだけ寄与するかを表します。この行列は $\bm{q}$ の関数であり、マニピュレータの姿勢が変わるたびに変化します。
この行列表現から、ベース速度を関節速度の関数として解くことができます。$\bm{H}_b$ は正定値行列（質量行列の性質から）なので逆行列が存在し、
$$ \dot{\bm{x}}_b = \bm{H}_b^{-1}(\bm{P}_0 – \bm{H}_{bm} \dot{\bm{q}}) $$
が得られます。この式が宇宙ロボティクスの最も基本的な関係式の一つです。マニピュレータの関節速度 $\dot{\bm{q}}$ が決まれば、運動量保存からベース衛星の速度が一意に決まることを示しています。
特に初期運動量がゼロ $\bm{P}_0 = \bm{0}$ の場合、
$$ \dot{\bm{x}}_b = -\bm{H}_b^{-1} \bm{H}_{bm} \dot{\bm{q}} $$
となります。マイナス符号が付いていることに注目してください。これは、アームが正方向に動けばベースは負方向に動く — すなわち反動が生じることを数学的に表現しています。
運動量保存からベースの運動が決まることがわかりました。次に、この連成がマニピュレータの動力学全体にどう影響するか、つまり動的カップリングの詳細を見ていきましょう。
動的カップリングの物理
カップリングとは何か
「動的カップリング」（dynamic coupling）とは、系の一部分の運動が、別の部分の運動に影響を及ぼす現象です。自由浮遊ロボットにおける動的カップリングは、アームの動きがベース衛星の運動を引き起こし、そのベースの運動がさらにアームの手先位置に影響するという双方向の連鎖として現れます。
日常的なアナロジーで理解しましょう。回転椅子に座って、両手に重りを持っているとします。足を床から浮かせた状態（外部トルクなし）で、腕を身体の右側に振ると、身体は左に回転します。これがまさに角運動量保存に基づく動的カップリングです。さらに面白いことに、腕を元の位置に戻すと身体も元の姿勢に戻ります — 角運動量がゼロに保存されているからです。
1自由度アームでの具体例
動的カップリングを最も明確に理解するために、ベース衛星に1自由度の回転関節アームが取り付けられた2次元平面モデルを考えます。
ベース衛星の質量を $m_0$、慣性モーメントを $I_0$、アームリンクの質量を $m_1$、長さを $l_1$、重心までの距離を $l_{c1}$、慣性モーメントを $I_1$ とします。ベースの回転角を $\theta_0$、関節角を $q_1$ とします。
系全体の角運動量保存（初期角運動量ゼロの場合）を書き下すと、
$$ L = I_0 \dot{\theta}_0 + I_1 (\dot{\theta}_0 + \dot{q}_1) + m_1 l_{c1}^2 (\dot{\theta}_0 + \dot{q}_1) + m_1 d \cdot l_{c1} \cos(q_1) (2\dot{\theta}_0 + \dot{q}_1) = 0 $$
ここで $d$ はベース重心から関節までの距離です。
この式を少し整理しましょう。まず $\dot{\theta}_0$ を含む項と $\dot{q}_1$ を含む項に分けます。
$$ \underbrace{(I_0 + I_1 + m_1 l_{c1}^2 + 2m_1 d \cdot l_{c1} \cos q_1)}_{H_b} \dot{\theta}_0 + \underbrace{(I_1 + m_1 l_{c1}^2 + m_1 d \cdot l_{c1} \cos q_1)}_{H_{bm}} \dot{q}_1 = 0 $$
$\dot{\theta}_0$ について解くと、
$$ \dot{\theta}_0 = -\frac{H_{bm}}{H_b} \dot{q}_1 $$
この式が動的カップリングの本質を端的に示しています。$H_{bm}/H_b$ は「カップリング比」とでも呼ぶべき量で、アームの関節速度 $\dot{q}_1$ に対するベースの角速度応答の比率を表します。
カップリング比の物理的解釈
カップリング比 $\alpha = H_{bm}/H_b$ の値を考えましょう。
$H_b$ は系全体の慣性、$H_{bm}$ はアームの慣性とカップリング項を含むので、常に $0 < H_{bm} \le H_b$ が成り立ちます（等号はアームの慣性がベースの慣性に比べて無視できないほど大きい極限）。したがって、
$$ 0 < \alpha \le 1 $$
です。
- $\alpha \to 0$（$m_0 \gg m_1$ のとき）: ベース衛星がアームに比べて非常に重い。アームが動いてもベースはほとんど動かない。地上ロボットに近い状況
- $\alpha \to 1$: アームの慣性がベースと同程度。アームが動くと、ベースもほぼ同じ大きさで反対方向に動く。制御が最も困難な状況
例えば、ISSのCanadarm2では、ISS本体の質量（約420トン）がアームの質量（約1.8トン）に比べて圧倒的に大きいため、$\alpha$ は非常に小さく、ベースの反動は実質的に無視できるレベルです。一方、質量100kgの小型衛星に20kgのマニピュレータを搭載した場合、$\alpha$ はかなり大きな値となり、動的カップリングの効果が顕著に現れます。
カップリングが姿勢に与える影響
角運動量ゼロのもとで、関節角を $q_1 = 0$ から $q_1 = \Delta q$ まで変化させたとき、ベースの姿勢変化 $\Delta \theta_0$ を求めましょう。
$$ \Delta \theta_0 = -\int_0^{\Delta q} \frac{H_{bm}(q_1)}{H_b(q_1)} dq_1 $$
一般に $H_{bm}/H_b$ は $q_1$ に依存するため、この積分は関節角の経路に依存しうるように見えます。しかし1自由度の場合は経路は一意なので、$\Delta \theta_0$ は $\Delta q$ だけで決まります。
2自由度以上の場合は事情が異なり、関節空間での経路に依存して最終的なベース姿勢が変わりうるという、非ホロノミックな性質が現れます。これは後続の記事で詳しく扱う深いトピックです。
重要なポイントとして、1自由度系では関節を元に戻せばベース姿勢も元に戻りますが、2自由度以上では関節を元に戻してもベース姿勢が元に戻らない場合があります。これは宇宙遊泳中の宇宙飛行士が、関節の動きを工夫することで（外部トルクなしに）身体の向きを変える現象と同じ原理です。
動的カップリングの物理的メカニズムが明らかになったので、次は系全体の運動方程式を定式化し、ベースとアームの連成動力学を統一的に記述する枠組みを構築しましょう。
自由浮遊マニピュレータの運動方程式
ラグランジュ法による導出
系全体の運動方程式を、ラグランジュ法に基づいて導出します。前提記事「ラグランジュ法によるマニピュレータ動力学」で扱った手法を、ベース衛星が自由に動ける場合に拡張します。
系のラグランジアン $\mathcal{L}$ は運動エネルギー $T$ から位置エネルギー $U$ を引いたものです。微小重力環境では $U \approx 0$ とできるので、$\mathcal{L} = T$ です。
系全体の運動エネルギーは、
$$ T = \frac{1}{2} m_0 \bm{v}_0^\top \bm{v}_0 + \frac{1}{2} \bm{\omega}_0^\top \bm{I}_0 \bm{\omega}_0 + \sum_{i=1}^{n} \left(\frac{1}{2} m_i \bm{v}_i^\top \bm{v}_i + \frac{1}{2} \bm{\omega}_i^\top \bm{I}_i \bm{\omega}_i\right) $$
前セクションで導入したヤコビ行列を用いて、全てのボディの速度をベース速度 $\dot{\bm{x}}_b$ と関節速度 $\dot{\bm{q}}$ で表し、これを運動エネルギーの式に代入します。結果として、運動エネルギーは次の二次形式にまとめられます。
$$ T = \frac{1}{2} \begin{pmatrix} \dot{\bm{x}}_b \\ \dot{\bm{q}} \end{pmatrix}^\top \begin{pmatrix} \bm{H}_b & \bm{H}_{bm} \\ \bm{H}_{bm}^\top & \bm{H}_m \end{pmatrix} \begin{pmatrix} \dot{\bm{x}}_b \\ \dot{\bm{q}} \end{pmatrix} $$
ここで、
- $\bm{H}_b \in \mathbb{R}^{6 \times 6}$: ベースの一般化慣性行列（前出）
- $\bm{H}_{bm} \in \mathbb{R}^{6 \times n}$: ベース・マニピュレータ連成慣性行列（前出）
- $\bm{H}_m \in \mathbb{R}^{n \times n}$: マニピュレータの慣性行列
これらの行列は全て関節角 $\bm{q}$ に依存する構成依存行列です。全体の慣性行列
$$ \bm{H} = \begin{pmatrix} \bm{H}_b & \bm{H}_{bm} \\ \bm{H}_{bm}^\top & \bm{H}_m \end{pmatrix} \in \mathbb{R}^{(6+n) \times (6+n)} $$
は対称正定値であり、ロボット全体の慣性特性を記述します。
運動方程式の構造
ラグランジュ方程式を適用すると、系全体の運動方程式は次の形になります。
$$ \begin{pmatrix} \bm{H}_b & \bm{H}_{bm} \\ \bm{H}_{bm}^\top & \bm{H}_m \end{pmatrix} \begin{pmatrix} \ddot{\bm{x}}_b \\ \ddot{\bm{q}} \end{pmatrix} + \begin{pmatrix} \bm{c}_b \\ \bm{c}_m \end{pmatrix} = \begin{pmatrix} \bm{F}_b \\ \bm{\tau} \end{pmatrix} $$
ここで、
- $\bm{c}_b \in \mathbb{R}^6$, $\bm{c}_m \in \mathbb{R}^n$: コリオリ力・遠心力に対応する非線形項
- $\bm{F}_b \in \mathbb{R}^6$: ベースに作用する外力・外トルク
- $\bm{\tau} \in \mathbb{R}^n$: 関節トルク
この運動方程式を、ベースの方程式とマニピュレータの方程式に分けて書くと、
$$ \bm{H}_b \ddot{\bm{x}}_b + \bm{H}_{bm} \ddot{\bm{q}} + \bm{c}_b = \bm{F}_b \quad \text{(ベースの運動方程式)} $$
$$ \bm{H}_{bm}^\top \ddot{\bm{x}}_b + \bm{H}_m \ddot{\bm{q}} + \bm{c}_m = \bm{\tau} \quad \text{(マニピュレータの運動方程式)} $$
2つの方程式が $\bm{H}_{bm}$ という連成慣性行列を介して結合していることが一目瞭然です。ベースの方程式にはマニピュレータの加速度 $\ddot{\bm{q}}$ が含まれ、マニピュレータの方程式にはベースの加速度 $\ddot{\bm{x}}_b$ が含まれます。この交差項こそが動的カップリングの数学的表現です。
自由浮遊モードでの簡略化
自由浮遊モードでは $\bm{F}_b = \bm{0}$（外力・外トルクなし）なので、ベースの方程式は
$$ \bm{H}_b \ddot{\bm{x}}_b + \bm{H}_{bm} \ddot{\bm{q}} + \bm{c}_b = \bm{0} $$
となります。これを $\ddot{\bm{x}}_b$ について解くと、
$$ \ddot{\bm{x}}_b = -\bm{H}_b^{-1}(\bm{H}_{bm} \ddot{\bm{q}} + \bm{c}_b) $$
この式をマニピュレータの運動方程式に代入して $\ddot{\bm{x}}_b$ を消去すると、関節空間だけで記述された運動方程式が得られます。
$\bm{H}_{bm}^\top(-\bm{H}_b^{-1})(\bm{H}_{bm} \ddot{\bm{q}} + \bm{c}_b) + \bm{H}_m \ddot{\bm{q}} + \bm{c}_m = \bm{\tau}$ を整理すると、
$$ \underbrace{(\bm{H}_m – \bm{H}_{bm}^\top \bm{H}_b^{-1} \bm{H}_{bm})}_{\bm{H}^*} \ddot{\bm{q}} + \underbrace{(\bm{c}_m – \bm{H}_{bm}^\top \bm{H}_b^{-1} \bm{c}_b)}_{\bm{c}^*} = \bm{\tau} $$
これが自由浮遊マニピュレータの縮約された運動方程式です。形式的には地上ロボットの運動方程式 $\bm{H}_m \ddot{\bm{q}} + \bm{c}_m = \bm{\tau}$ と同じ $n$ 次元の方程式ですが、慣性行列と非線形項がベースとの連成によって修正されています。
ここで現れた $\bm{H}^*$ が、次に議論する「有効質量行列」（effective inertia matrix）の概念につながります。
有効質量行列
概念と定義
前セクションで導出した縮約された運動方程式の慣性行列
$$ \bm{H}^* = \bm{H}_m – \bm{H}_{bm}^\top \bm{H}_b^{-1} \bm{H}_{bm} $$
を有効慣性行列（effective inertia matrix）と呼びます。
この行列の物理的意味を直感的に理解しましょう。地上ロボットでは、関節トルク $\bm{\tau}$ が生み出す関節加速度は $\ddot{\bm{q}} = \bm{H}_m^{-1} \bm{\tau}$（非線形項を無視した場合）です。つまり、アームを加速するために必要なトルクは慣性行列 $\bm{H}_m$ で決まります。
ところが自由浮遊ロボットでは、関節トルクの一部がベース衛星を動かすために「消費」されます。その結果、同じ関節加速度を生み出すのに必要なトルクが $\bm{H}_m$ ではなく $\bm{H}^*$ で決まります。
$\bm{H}^*$ と $\bm{H}_m$ の差は、
$$ \bm{H}_m – \bm{H}^* = \bm{H}_{bm}^\top \bm{H}_b^{-1} \bm{H}_{bm} $$
です。右辺は正の半定値行列（$\bm{H}_b^{-1}$ が正定値なので）であるため、$\bm{H}^* \preceq \bm{H}_m$ となります。つまり、有効慣性は常にマニピュレータ単体の慣性以下です。
これは直感に合います。自由浮遊モードでは、関節トルクの反作用がベースを動かすため、実質的にアームが「軽く」感じられるのです。言い換えると、同じトルクで地上よりも大きな関節加速度が得られます。ただし、ベースの反動という副作用が伴います。
シューア補行列との関係
数学的に見ると、$\bm{H}^*$ は全体の慣性行列 $\bm{H}$ をブロック分割したときのシューア補行列（Schur complement）です。
$$ \bm{H}^* = \bm{H}/\bm{H}_b = \bm{H}_m – \bm{H}_{bm}^\top \bm{H}_b^{-1} \bm{H}_{bm} $$
シューア補行列は、連立方程式でブロック消去を行ったときに自然に現れる行列です。ここでは、ベースの変数 $\ddot{\bm{x}}_b$ を運動量保存（ベースの運動方程式）を使って消去した結果として現れています。
シューア補行列の重要な性質として、$\bm{H}$ が正定値であり $\bm{H}_b$ が正定値であれば、$\bm{H}^*$ も正定値であることが保証されます。これは物理的にも当然で、有効慣性がゼロや負になることはあり得ません。
手先の有効質量
ここまでは関節空間での有効慣性行列を議論しましたが、同様の概念を作業空間（タスク空間）でも定義できます。手先に外力 $\bm{F}_\mathrm{ee}$ を加えたときの手先の加速度応答を考えると、手先の有効質量行列 $\bm{\Lambda}^*$ が定義できます。
$$ \bm{\Lambda}^* = (\bm{J}^* (\bm{H}^*)^{-1} (\bm{J}^*)^\top)^{-1} $$
ここで $\bm{J}^*$ は後続記事で詳しく扱う一般化ヤコビ行列（Generalized Jacobian Matrix）です。
手先の有効質量は、ロボットが対象物に接触したときの力学的応答を特徴づける量であり、軌道上サービスやデブリ捕獲の衝突解析で極めて重要です。「一般化ヤコビ行列」の記事で詳しく議論します。
有効質量行列の概念を理解した上で、いよいよ具体的な数値シミュレーションに移り、ここまで議論してきた理論が実際にどのような挙動として現れるかを確認しましょう。
Pythonシミュレーション — 2次元自由浮遊ロボットの反動
シミュレーションの設定
ここでは、理論の本質を見るために2次元平面上のベース衛星 + 1リンクアームのモデルをPythonでシミュレーションします。このモデルは単純ですが、自由浮遊ロボットの動的カップリングの核心を全て含んでいます。
以下の物理パラメータを使います。
| パラメータ | 記号 | 値 |
|---|---|---|
| ベース質量 | $m_0$ | 100 kg |
| ベース慣性モーメント | $I_0$ | 50 kg$\cdot$m$^2$ |
| アーム質量 | $m_1$ | 20 kg |
| アーム長さ | $l_1$ | 2.0 m |
| アーム重心位置 | $l_{c1}$ | 1.0 m |
| アーム慣性モーメント | $I_1$ | 6.67 kg$\cdot$m$^2$ |
| ベース重心から関節までの距離 | $d$ | 0.8 m |
初期条件は、系全体が静止している状態（線運動量ゼロ、角運動量ゼロ）からスタートします。
まず、運動量保存から導かれるベースの運動を計算する基本関数を実装します。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch, FancyArrowPatch
from matplotlib.collections import LineCollection
from scipy.integrate import solve_ivp
# --- 物理パラメータ ---
m0 = 100.0 # ベース質量 [kg]
I0 = 50.0 # ベース慣性モーメント [kg*m^2]
m1 = 20.0 # アーム質量 [kg]
l1 = 2.0 # アーム長さ [m]
lc1 = 1.0 # アーム重心までの距離 [m]
I1 = 6.67 # アーム慣性モーメント [kg*m^2]
d = 0.8 # ベース重心から関節までの距離 [m]
def coupling_ratio(q1):
"""角運動量保存に基づくカップリング比 alpha(q1) を計算"""
H_b = I0 + I1 + m1 * lc1**2 + 2 * m1 * d * lc1 * np.cos(q1)
H_bm = I1 + m1 * lc1**2 + m1 * d * lc1 * np.cos(q1)
return H_bm / H_b
def base_angular_velocity(q1, dq1):
"""角運動量保存からベースの角速度を計算（初期角運動量=0）"""
alpha = coupling_ratio(q1)
return -alpha * dq1
def compute_base_angle(q1_array):
"""関節角の軌跡からベース角を数値積分で計算"""
theta0 = np.zeros_like(q1_array)
for i in range(1, len(q1_array)):
dq = q1_array[i] - q1_array[i-1]
q_mid = 0.5 * (q1_array[i] + q1_array[i-1])
alpha = coupling_ratio(q_mid)
theta0[i] = theta0[i-1] - alpha * dq
return theta0
このコードでは、先ほど導出した角運動量保存の式 $\dot{\theta}_0 = -\alpha(q_1) \dot{q}_1$ を実装しています。coupling_ratio 関数は関節角 $q_1$ に依存するカップリング比 $\alpha = H_{bm}/H_b$ を返します。compute_base_angle は関節角の時系列からベース角を数値積分します。
カップリング比の可視化
まず、カップリング比 $\alpha$ が関節角にどう依存するかを確認します。
import numpy as np
import matplotlib.pyplot as plt
# カップリング比の関節角依存性
q1_range = np.linspace(-np.pi, np.pi, 500)
alpha_values = np.array([coupling_ratio(q) for q in q1_range])
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(np.degrees(q1_range), alpha_values, color='#00bcd4', linewidth=2.5)
ax.set_xlabel('Joint angle $q_1$ [deg]', fontsize=13)
ax.set_ylabel('Coupling ratio $\\alpha$', fontsize=13)
ax.set_title('Coupling Ratio vs Joint Angle', fontsize=14)
ax.grid(True, alpha=0.3)
ax.axhline(y=alpha_values.mean(), color='#ff9800', linestyle='--',
linewidth=1.5, label=f'Mean = {alpha_values.mean():.3f}')
ax.legend(fontsize=12)
ax.set_xlim(-180, 180)
plt.tight_layout()
plt.show()
このグラフから、2つの重要な特徴が読み取れます。
- カップリング比は $q_1 = 0$ で最大、$q_1 = \pm 180°$ で最小: アームがベースと同じ方向に伸びている（$q_1 = 0$）とき、アームの回転慣性が系全体に対して最も大きな割合を占めるため、ベースの反動も最大になります。アームが折り畳まれている（$q_1 = \pm 180°$）ときは、アームの重心がベース重心に近づくため、結合が弱まります。
- カップリング比は常に0より大きい: どの姿勢でも $\alpha > 0$ であり、自由浮遊モードではベースの反動を完全になくすことはできません。この事実は、自由浮遊ロボットの制御が本質的に困難である理由の一つです。
ベース衛星の反動シミュレーション
次に、アームに正弦波状の関節角軌道を与えたときのベースの応答をシミュレーションします。
import numpy as np
import matplotlib.pyplot as plt
# 関節角の軌道：正弦波（1周期分を往復）
t = np.linspace(0, 10, 2000)
q1_amplitude = np.radians(60) # 振幅 60度
q1_traj = q1_amplitude * np.sin(2 * np.pi * t / 10)
# ベース角の計算
theta0_traj = compute_base_angle(q1_traj)
# 関節角速度の計算
dq1_traj = np.gradient(q1_traj, t)
dtheta0_traj = np.array([base_angular_velocity(q, dq)
for q, dq in zip(q1_traj, dq1_traj)])
fig, axes = plt.subplots(3, 1, figsize=(10, 10), sharex=True)
# 関節角
axes[0].plot(t, np.degrees(q1_traj), color='#2196f3', linewidth=2, label='$q_1$')
axes[0].set_ylabel('Joint angle [deg]', fontsize=12)
axes[0].legend(fontsize=12)
axes[0].grid(True, alpha=0.3)
axes[0].set_title('Free-Floating Robot: Joint Motion and Base Reaction', fontsize=14)
# ベース角
axes[1].plot(t, np.degrees(theta0_traj), color='#f44336', linewidth=2, label='$\\theta_0$')
axes[1].set_ylabel('Base angle [deg]', fontsize=12)
axes[1].legend(fontsize=12)
axes[1].grid(True, alpha=0.3)
# 角速度比較
axes[2].plot(t, np.degrees(dq1_traj), color='#2196f3', linewidth=2,
label='$\\dot{q}_1$', alpha=0.7)
axes[2].plot(t, np.degrees(dtheta0_traj), color='#f44336', linewidth=2,
label='$\\dot{\\theta}_0$', alpha=0.7)
axes[2].set_ylabel('Angular velocity [deg/s]', fontsize=12)
axes[2].set_xlabel('Time [s]', fontsize=12)
axes[2].legend(fontsize=12)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このシミュレーション結果から、自由浮遊ロボットの挙動について以下のことが確認できます。
- ベースの反動は関節運動と逆位相: 関節角 $q_1$ が正に増加すると、ベース角 $\theta_0$ は負に変化します。これは角運動量保存の直接的な帰結であり、$\dot{\theta}_0 = -\alpha \dot{q}_1$ という関係がグラフに明瞭に現れています。
- ベースの反動振幅は関節角振幅より小さい: カップリング比 $\alpha < 1$（ベースがアームより重い設定）なので、ベースの角度変化は関節角度変化より小さくなっています。$m_0/m_1 = 5$ という質量比から、ベースの反動は関節運動の約30-40%程度に抑えられています。
- 正弦波入力に対して、ベースも周期運動する: 1自由度の場合、関節角が元に戻れば、ベース角も元に戻ります。ただし、カップリング比が関節角に依存するため、ベースの波形は完全な正弦波にはならず、わずかに歪んでいます。
ベース並進運動のシミュレーション
角運動量だけでなく、線運動量の保存も考慮したシミュレーションを行いましょう。ベース衛星の並進運動（重心の固定条件）も含めた完全な2次元シミュレーションです。
import numpy as np
import matplotlib.pyplot as plt
def simulate_free_floating_2d(q1_traj, t, m0, m1, d, lc1, I0, I1):
"""
2次元自由浮遊ロボットの完全シミュレーション
線運動量保存 + 角運動量保存を同時に満たす
"""
N = len(t)
# ベースの状態
x0 = np.zeros(N) # ベース x 位置
y0 = np.zeros(N) # ベース y 位置
theta0 = np.zeros(N) # ベース角度
# 重心位置（初期）
M = m0 + m1
# 初期状態：ベース原点、アームは+x方向に伸びている
xG0 = (m0 * 0 + m1 * (d + lc1)) / M
yG0 = 0.0
for i in range(N):
q1 = q1_traj[i]
if i == 0:
theta0[i] = 0.0
else:
# 角運動量保存からベース角を計算
dq = q1_traj[i] - q1_traj[i-1]
q_mid = 0.5 * (q1_traj[i] + q1_traj[i-1])
th_mid = theta0[i-1]
alpha = coupling_ratio(q_mid)
theta0[i] = theta0[i-1] - alpha * dq
# 線運動量保存（重心位置一定）から x0, y0 を計算
# 関節の絶対位置（ベースからの相対位置を回転）
joint_x = x0[0] + d * np.cos(theta0[i]) # 暫定
joint_y = y0[0] + d * np.sin(theta0[i]) # 暫定
# アーム重心の絶対位置
arm_abs_angle = theta0[i] + q1
arm_cx = joint_x + lc1 * np.cos(arm_abs_angle) # 暫定
arm_cy = joint_y + lc1 * np.sin(arm_abs_angle) # 暫定
# 重心一定条件: m0*(x0,y0) + m1*(arm_cx, arm_cy) = M*(xG0, yG0)
# arm_cx = x0 + d*cos(theta0) + lc1*cos(theta0+q1)
# arm_cy = y0 + d*sin(theta0) + lc1*sin(theta0+q1)
# m0*x0 + m1*(x0 + d*cos(theta0) + lc1*cos(theta0+q1)) = M*xG0
# (m0+m1)*x0 + m1*(d*cos(theta0) + lc1*cos(theta0+q1)) = M*xG0
x0[i] = (M * xG0 - m1 * (d * np.cos(theta0[i])
+ lc1 * np.cos(theta0[i] + q1))) / M
y0[i] = (M * yG0 - m1 * (d * np.sin(theta0[i])
+ lc1 * np.sin(theta0[i] + q1))) / M
return x0, y0, theta0
# シミュレーション実行
t = np.linspace(0, 10, 2000)
q1_traj = np.radians(60) * np.sin(2 * np.pi * t / 10)
x0, y0, theta0 = simulate_free_floating_2d(
q1_traj, t, m0, m1, d, lc1, I0, I1
)
# 手先位置の計算
joint_x = x0 + d * np.cos(theta0)
joint_y = y0 + d * np.sin(theta0)
ee_x = joint_x + l1 * np.cos(theta0 + q1_traj)
ee_y = joint_y + l1 * np.sin(theta0 + q1_traj)
# 固定ベースの場合の手先位置（比較用）
ee_x_fixed = d + l1 * np.cos(q1_traj)
ee_y_fixed = l1 * np.sin(q1_traj)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左: ベース位置の軌跡
axes[0].plot(x0, y0, color='#f44336', linewidth=2, label='Base CoM trajectory')
axes[0].plot(x0[0], y0[0], 'o', color='#f44336', markersize=10,
label='Start', zorder=5)
axes[0].plot(x0[-1], y0[-1], 's', color='#ff9800', markersize=10,
label='End', zorder=5)
# 重心位置（一定のはず）
M = m0 + m1
xG = (m0 * x0 + m1 * (x0 + d * np.cos(theta0)
+ lc1 * np.cos(theta0 + q1_traj))) / M
yG = (m0 * y0 + m1 * (y0 + d * np.sin(theta0)
+ lc1 * np.sin(theta0 + q1_traj))) / M
axes[0].plot(xG[0], yG[0], '*', color='#4caf50', markersize=15,
label='System CoM (fixed)', zorder=5)
axes[0].set_xlabel('x [m]', fontsize=12)
axes[0].set_ylabel('y [m]', fontsize=12)
axes[0].set_title('Base Satellite Trajectory', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_aspect('equal')
# 右: 手先軌跡の比較
axes[1].plot(ee_x, ee_y, color='#2196f3', linewidth=2,
label='Free-floating end-effector')
axes[1].plot(ee_x_fixed, ee_y_fixed, '--', color='#9e9e9e', linewidth=2,
label='Fixed-base end-effector')
axes[1].plot(ee_x[0], ee_y[0], 'o', color='#2196f3', markersize=10, zorder=5)
axes[1].plot(ee_x_fixed[0], ee_y_fixed[0], 'o', color='#9e9e9e',
markersize=10, zorder=5)
axes[1].set_xlabel('x [m]', fontsize=12)
axes[1].set_ylabel('y [m]', fontsize=12)
axes[1].set_title('End-Effector Trajectory Comparison', fontsize=13)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_aspect('equal')
plt.tight_layout()
plt.show()
この比較グラフからは、自由浮遊ロボットの本質的な特徴が明確に読み取れます。
- ベース衛星は閉曲線軌道を描く（左図）: アームが正弦波運動すると、ベースは線運動量保存（重心一定条件）に従って楕円に近い閉曲線上を移動します。開始点と終了点が一致していることが確認でき、アームが元の姿勢に戻ればベースも元の位置に戻ることがわかります。
- 系の重心は完全に静止している（左図の緑色の星印）: 線運動量保存の直接的な帰結として、系の重心位置は変化しません。これはシミュレーションが運動量保存則を正しく実装できていることの検証にもなっています。
- 手先の作業空間が縮小している（右図）: 自由浮遊ロボット（青線）の手先軌跡は、固定ベースロボット（灰色破線）の手先軌跡よりも小さな領域を描いています。ベースの反動が手先の有効到達範囲を減少させているのです。これが宇宙ロボットの作業空間が地上ロボットよりも狭くなる根本原因です。
質量比の影響
動的カップリングの強さが質量比にどう依存するかを定量的に調べます。
import numpy as np
import matplotlib.pyplot as plt
mass_ratios = [2, 5, 10, 20, 50]
colors = plt.cm.viridis(np.linspace(0.2, 0.9, len(mass_ratios)))
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
t = np.linspace(0, 10, 2000)
q1_traj = np.radians(60) * np.sin(2 * np.pi * t / 10)
for idx, ratio in enumerate(mass_ratios):
m0_var = ratio * m1
I0_var = m0_var * 0.5 # 慣性モーメントを質量に比例させる
# カップリング比の計算（パラメータを変えて）
def coupling_ratio_var(q1):
H_b = I0_var + I1 + m1 * lc1**2 + 2 * m1 * d * lc1 * np.cos(q1)
H_bm = I1 + m1 * lc1**2 + m1 * d * lc1 * np.cos(q1)
return H_bm / H_b
# ベース角の計算
theta0 = np.zeros_like(q1_traj)
for i in range(1, len(q1_traj)):
dq = q1_traj[i] - q1_traj[i-1]
q_mid = 0.5 * (q1_traj[i] + q1_traj[i-1])
alpha = coupling_ratio_var(q_mid)
theta0[i] = theta0[i-1] - alpha * dq
axes[0].plot(t, np.degrees(theta0), color=colors[idx], linewidth=2,
label=f'$m_0/m_1 = {ratio}$')
# 最大ベース角変位
max_displacement = np.max(np.abs(np.degrees(theta0)))
axes[1].bar(idx, max_displacement, color=colors[idx], width=0.7)
axes[0].set_xlabel('Time [s]', fontsize=12)
axes[0].set_ylabel('Base angle $\\theta_0$ [deg]', fontsize=12)
axes[0].set_title('Base Reaction for Different Mass Ratios', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[1].set_xticks(range(len(mass_ratios)))
axes[1].set_xticklabels([f'{r}:1' for r in mass_ratios])
axes[1].set_xlabel('Mass ratio $m_0 : m_1$', fontsize=12)
axes[1].set_ylabel('Max base angular displacement [deg]', fontsize=12)
axes[1].set_title('Peak Base Reaction vs Mass Ratio', fontsize=13)
axes[1].grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.show()
質量比のパラメトリックスタディから、以下の傾向が明確になります。
- 質量比が小さいほどベースの反動が大きい: $m_0/m_1 = 2$ のとき、ベース角度の最大変位は30度近くにもなります。一方、$m_0/m_1 = 50$ では数度以下に抑えられます。これは直感に合う結果です。重いベースほど動きにくく、「固定ベース」に近づきます。
- 反動の抑制は非線形的: 質量比を2倍にしても、反動が半分になるわけではありません。これは $\alpha = H_{bm}/H_b$ の関数形が質量比の単純な逆数ではないためです。特に、質量比が小さい領域（2〜5）では反動の変化が急激であり、この領域がミッション設計上最も注意を要する領域です。
ロボット姿勢のアニメーション用スナップショット
最後に、自由浮遊ロボットの姿勢変化を時系列のスナップショットで可視化します。
import numpy as np
import matplotlib.pyplot as plt
# シミュレーション
t = np.linspace(0, 10, 2000)
q1_traj = np.radians(90) * np.sin(2 * np.pi * t / 10)
x0, y0, theta0 = simulate_free_floating_2d(
q1_traj, t, m0, m1, d, lc1, I0, I1
)
# スナップショットの時刻
snapshot_indices = np.linspace(0, len(t)//2, 8, dtype=int)
fig, ax = plt.subplots(figsize=(14, 8))
base_half = 0.4 # ベースの描画サイズ
for idx, si in enumerate(snapshot_indices):
alpha_fade = 0.3 + 0.7 * (idx / (len(snapshot_indices) - 1))
bx, by = x0[si], y0[si]
th = theta0[si]
q1 = q1_traj[si]
# ベースの四隅
corners = np.array([[-base_half, -base_half],
[base_half, -base_half],
[base_half, base_half],
[-base_half, base_half],
[-base_half, -base_half]])
R = np.array([[np.cos(th), -np.sin(th)],
[np.sin(th), np.cos(th)]])
rotated = (R @ corners.T).T + np.array([bx, by])
ax.fill(rotated[:, 0], rotated[:, 1], color='#f44336',
alpha=alpha_fade * 0.5, edgecolor='#d32f2f',
linewidth=1.5)
# 関節位置
jx = bx + d * np.cos(th)
jy = by + d * np.sin(th)
# ベースから関節への線
ax.plot([bx, jx], [by, jy], color='#9e9e9e',
linewidth=2, alpha=alpha_fade)
# アーム
arm_angle = th + q1
ex = jx + l1 * np.cos(arm_angle)
ey = jy + l1 * np.sin(arm_angle)
ax.plot([jx, ex], [jy, ey], color='#2196f3',
linewidth=4, alpha=alpha_fade, solid_capstyle='round')
# 手先マーカー
ax.plot(ex, ey, 'o', color='#4caf50', markersize=6,
alpha=alpha_fade, zorder=5)
# 時刻ラベル
ax.annotate(f't={t[si]:.1f}s', (bx, by + 0.6),
fontsize=8, ha='center', alpha=alpha_fade,
color='#333333')
# 重心の表示
M_total = m0 + m1
xG0 = m1 * (d + lc1) / M_total
ax.plot(xG0, 0, '*', color='#ff9800', markersize=20, zorder=10,
label='System CoM')
ax.set_xlabel('x [m]', fontsize=13)
ax.set_ylabel('y [m]', fontsize=13)
ax.set_title('Free-Floating Robot Snapshots (half period)', fontsize=14)
ax.set_aspect('equal')
ax.grid(True, alpha=0.2)
ax.legend(fontsize=12, loc='upper right')
plt.tight_layout()
plt.show()
このスナップショットは、自由浮遊ロボットの動的カップリングを直感的に理解するのに最適です。アームが反時計回り方向に回転すると、ベース衛星が時計回り方向に回転し、かつ左方向に並進している様子が見て取れます。時間が進むにつれて（色が濃くなるにつれて）、ベースの位置と姿勢が関節角に応じて系統的に変化しています。特に注目すべきは、系全体の重心（オレンジ色の星印）がどのスナップショットでも同じ位置にある点です。これは線運動量保存の目に見える証拠です。
有効慣性行列の計算と可視化
最後に、有効慣性行列の値が関節角によってどう変化するかを確認します。
import numpy as np
import matplotlib.pyplot as plt
def effective_inertia(q1, m0, m1, d, lc1, I0, I1):
"""有効慣性を計算（1自由度の場合はスカラー）"""
H_b = I0 + I1 + m1 * lc1**2 + 2 * m1 * d * lc1 * np.cos(q1)
H_bm = I1 + m1 * lc1**2 + m1 * d * lc1 * np.cos(q1)
H_m = I1 + m1 * lc1**2
H_star = H_m - H_bm**2 / H_b
return H_m, H_star
q1_range = np.linspace(-np.pi, np.pi, 500)
H_m_vals = np.zeros_like(q1_range)
H_star_vals = np.zeros_like(q1_range)
for i, q1 in enumerate(q1_range):
H_m_vals[i], H_star_vals[i] = effective_inertia(
q1, m0, m1, d, lc1, I0, I1
)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 有効慣性の比較
axes[0].plot(np.degrees(q1_range), H_m_vals, '--', color='#9e9e9e',
linewidth=2, label='$H_m$ (fixed-base)')
axes[0].plot(np.degrees(q1_range), H_star_vals, color='#2196f3',
linewidth=2.5, label='$H^*$ (free-floating)')
axes[0].fill_between(np.degrees(q1_range), H_star_vals, H_m_vals,
alpha=0.15, color='#f44336',
label='Inertia reduction')
axes[0].set_xlabel('Joint angle $q_1$ [deg]', fontsize=12)
axes[0].set_ylabel('Inertia [kg$\\cdot$m$^2$]', fontsize=12)
axes[0].set_title('Effective Inertia vs Joint Angle', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(-180, 180)
# 慣性低減率
reduction_ratio = (H_m_vals - H_star_vals) / H_m_vals * 100
axes[1].plot(np.degrees(q1_range), reduction_ratio, color='#f44336',
linewidth=2.5)
axes[1].set_xlabel('Joint angle $q_1$ [deg]', fontsize=12)
axes[1].set_ylabel('Inertia reduction [%]', fontsize=12)
axes[1].set_title('Inertia Reduction Ratio', fontsize=13)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(-180, 180)
plt.tight_layout()
plt.show()
有効慣性行列の可視化から、以下の物理的洞察が得られます。
- 有効慣性は常にマニピュレータ単体の慣性より小さい（左図の赤色塗りつぶし領域）: $\bm{H}^* = H_m – H_{bm}^2/H_b < H_m$ が全ての関節角で成り立っています。ベースが自由に動くことで、マニピュレータが実質的に「軽く」なることが数値的に確認できました。
- 慣性低減率は関節角に依存する（右図）: $q_1 = 0$ のときに最も大きな慣性低減（約15-20%）が起こり、$q_1 = \pm 180°$ のときに最小になります。これはカップリング比の振る舞いと完全に整合しています。アームがベースと同じ方向に伸びた姿勢ではカップリングが最も強く、慣性低減も最大です。
- 実用的な意味: 慣性低減は「同じトルクでより大きな加速度が得られる」ことを意味しますが、その代償としてベースが反動で動きます。制御設計では、この有効慣性を正しく考慮しないと、トルク指令に対する応答が予測と異なり、精密な手先制御が困難になります。
まとめ
本記事では、宇宙空間における自由浮遊ロボットの動力学の基礎を解説しました。
- 地上ロボットとの本質的な違い: ベースが固定されていないため、マニピュレータの動きがベース衛星に反力を与え、系全体が連成系として振る舞う
- 自由浮遊 vs 自由飛行: スラスタを使わない自由浮遊モードでは推進剤を節約できるが、運動量保存に基づくベースの反動を受け入れる必要がある
- 運動量保存: 外力・外トルクがない自由浮遊モードでは線運動量と角運動量が保存され、ベース速度はマニピュレータの関節速度の関数として一意に決まる（$\dot{\bm{x}}_b = -\bm{H}_b^{-1} \bm{H}_{bm} \dot{\bm{q}}$）
- 動的カップリング: アームの動きがベースを動かし、そのベースの動きがさらに手先位置に影響するという双方向の連鎖。カップリング比は質量比と関節角に依存する
- 運動方程式: ベースとマニピュレータが連成慣性行列 $\bm{H}_{bm}$ で結合され、自由浮遊モードではベース変数を消去して $\bm{H}^* \ddot{\bm{q}} + \bm{c}^* = \bm{\tau}$ に縮約される
- 有効質量行列: 自由浮遊の効果でマニピュレータの実質的な慣性がシューア補行列 $\bm{H}^* = \bm{H}_m – \bm{H}_{bm}^\top \bm{H}_b^{-1} \bm{H}_{bm}$ に低減される
本記事で導入した動的カップリングと運動量保存の枠組みは、宇宙ロボティクスの全ての応用の基盤となります。次のステップとして、以下の記事も参考にしてください。
- 一般化ヤコビ行列 — 運動量保存を考慮した運動学。自由浮遊ロボットの手先速度と関節速度の関係を定式化し、作業空間での制御に不可欠な道具を導入