宇宙ステーションに補給船がドッキングする映像を見たことがあるでしょうか。地上から見ると「ゆっくり近づいて、そっとくっつく」だけのように見えますが、実は時速28,000 km以上で地球を周回する2つの物体を、数センチメートルの精度で合わせるという途方もない精密制御が行われています。しかも最近のミッションでは、地上の管制官の指示を待たずに、宇宙機が自律的に判断して接近する技術が主流になっています。
なぜこの技術を学ぶのでしょうか？ 自律ランデブー・近傍接近（Autonomous Rendezvous and Proximity Operations: ARPO）は、宇宙開発における最も基盤的な運用能力の1つです。
- 宇宙ステーション補給: ロシアのプログレス補給船やスペースXのドラゴン宇宙船は、ISSに自律的に接近・ドッキングします。通信遅延がある深宇宙では人間の操縦が困難なため、自律化は必須です
- 軌道上サービス（OOS）: 故障した衛星の修理や燃料補給、デブリ除去では、非協力ターゲット（ドッキング用の装置がない物体）に接近する技術が求められます
- 将来の宇宙組立: 大型宇宙望遠鏡や宇宙太陽光発電所の軌道上組立では、多数のモジュールを精密に合体させるランデブー技術が不可欠です
本記事の内容
- ランデブー問題の全体像と各フェーズの定義
- CWH方程式（Clohessy-Wiltshire-Hill）の詳細な導出
- CWH解の物理的意味 — 自由漂流、安全楕円、周期的相対軌道
- V-barアプローチとR-barアプローチの比較
- 誘導法則の設計（二体力接近、フィードバック誘導、最適制御ベース）
- 安全性の設計思想（受動安全軌道、衝突回避マニューバ）
- Pythonによる CWH方程式のシミュレーションとランデブー軌道の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ランデブー問題とは — 2つの宇宙機を出会わせる
日常のアナロジーで考える
高速道路を時速100 kmで走る車に、別の車が後ろから合流する場面を想像してください。合流車は速度を合わせながら徐々に距離を縮め、最後は真横に並びます。宇宙でのランデブーもこれと本質的に同じですが、大きく異なる点が1つあります。宇宙では「まっすぐ近づく」ことができないのです。
地上では、前方の車に近づきたければアクセルを踏めば済みます。しかし軌道力学では、前方の宇宙機に追いつこうとして加速すると、衛星はより高い軌道に上がってしまい、かえって離れていきます。これは直感に反しますが、ケプラーの第3法則から理解できます。より高い軌道は周期が長くなるため、加速した宇宙機は目標から遅れていくのです。
この「加速すると遅れ、減速すると追いつく」という軌道力学の逆説的な性質が、ランデブー問題を興味深くも難しくしています。
ランデブーの定義
ランデブー（rendezvous）とは、2つの宇宙機の位置と速度を一致させることです。単に近くを通過する（フライバイ）とは異なり、相対速度をゼロにして並走状態を作ります。数学的に表現すると、チェイサー（追跡機）の位置・速度ベクトル $(\bm{r}_c, \bm{v}_c)$ をターゲット（目標機）の $(\bm{r}_t, \bm{v}_t)$ に一致させることです。
$$ \bm{r}_c \to \bm{r}_t, \quad \bm{v}_c \to \bm{v}_t $$
近傍接近運用（Proximity Operations: ProxOps）は、ランデブー後の近距離（通常数百メートル以内）でのマニューバ全般を指します。これにはドッキング前の保持（station-keeping）、最終接近、そしてドッキングやバーシング（ロボットアームによる捕獲）が含まれます。
ランデブー問題を解くための第一歩は、「ターゲットから見たチェイサーの運動」を記述する力学モデルです。次のセクションでは、この相対運動を支配する基本方程式であるCWH方程式を導出します。
相対運動の力学 — CWH方程式の導出
なぜ相対運動を考えるのか
ランデブー問題では、2つの宇宙機の絶対位置（地球中心座標）ではなく、ターゲットから見たチェイサーの相対位置が本質的に重要です。目標に対して「あと何メートル前方か」「あと何m/sの相対速度があるか」という情報が、噴射のタイミングと方向を決めるからです。
この考え方は、高速道路の合流と同じです。合流車のドライバーが気にするのは、自分の地上での速度ではなく、「前の車との相対速度と距離」です。軌道力学でも同様に、ターゲットに固定された座標系（LVLH座標系）での相対運動を記述します。
LVLH座標系の定義
ターゲットの重心を原点とする LVLH座標系（Local Vertical Local Horizontal）を定義します。
- $x$ 軸（V-bar）: ターゲットの速度方向（along-track、進行方向）
- $y$ 軸（H-bar）: 軌道面に垂直な方向（cross-track、法線方向）
- $z$ 軸（R-bar）: 地球中心からターゲットへ向かう方向（radial、地球から離れる方向を正）
文献によっては座標軸の取り方が異なることがあります。本記事では上記の定義に従いますが、後で示すCWH方程式はいずれの慣習でも本質は同じです。
二体問題の運動方程式から出発する
ターゲットとチェイサーはそれぞれ地球の重力場中で運動しています。地球中心慣性座標系（ECI）における運動方程式は、二体問題の枠組みで次のように書けます。
ターゲットの運動方程式:
$$ \ddot{\bm{R}}_t = -\frac{\mu}{R_t^3} \bm{R}_t $$
チェイサーの運動方程式:
$$ \ddot{\bm{R}}_c = -\frac{\mu}{R_c^3} \bm{R}_c + \bm{f} $$
ここで $\bm{R}_t$, $\bm{R}_c$ はそれぞれの地球中心からの位置ベクトル、$\mu$ は地球の重力パラメータ（$\mu = GM_\oplus \approx 3.986 \times 10^{14} \, \text{m}^3/\text{s}^2$）、$\bm{f}$ はチェイサーに働く制御力（スラスタ噴射）の加速度です。
相対位置ベクトルを $\bm{\rho} = \bm{R}_c – \bm{R}_t$ と定義すると、チェイサーの位置は $\bm{R}_c = \bm{R}_t + \bm{\rho}$ と書けます。
両辺の差を取ると、相対加速度は
$$ \ddot{\bm{\rho}} = \ddot{\bm{R}}_c – \ddot{\bm{R}}_t = -\frac{\mu}{|\bm{R}_t + \bm{\rho}|^3}(\bm{R}_t + \bm{\rho}) + \frac{\mu}{R_t^3}\bm{R}_t + \bm{f} $$
となります。この式は厳密ですが、非線形で扱いにくいです。ここからCWH方程式を得るために、2つの重要な仮定を導入します。
仮定1: ターゲットは円軌道
ターゲットが角速度 $n$（平均運動）の円軌道上にあると仮定します。
$$ n = \sqrt{\frac{\mu}{a^3}} $$
ここで $a$ はターゲットの円軌道半径です。ISS（高度約400 km）の場合、$a \approx 6,778 \, \text{km}$ で $n \approx 1.13 \times 10^{-3} \, \text{rad/s}$（軌道周期約92.5分）です。
仮定2: チェイサーはターゲットの近傍
相対距離が軌道半径に比べて十分小さい、すなわち $|\bm{\rho}| \ll R_t$ を仮定します。ISSの場合、軌道半径は約6,778 kmですから、相対距離が数十km以内であればこの近似は十分成り立ちます。
回転座標系での運動方程式
LVLH座標系は角速度 $\bm{\omega} = n \hat{\bm{y}}$（$y$ 軸回りに一定の角速度 $n$）で回転しています。回転座標系における加速度の公式を適用します。
慣性系での加速度 $\ddot{\bm{\rho}}_I$ と回転系での加速度 $\ddot{\bm{\rho}}_R$ の関係は
$$ \ddot{\bm{\rho}}_I = \ddot{\bm{\rho}}_R + 2\bm{\omega} \times \dot{\bm{\rho}}_R + \bm{\omega} \times (\bm{\omega} \times \bm{\rho})+ \dot{\bm{\omega}} \times \bm{\rho} $$
円軌道では $\dot{\bm{\omega}} = \bm{0}$ なので、最後の項は消えます。
LVLH座標系での相対位置を $\bm{\rho} = (x, y, z)$ とし、重力の線形化（テイラー展開の1次項まで）を行うと、各成分の運動方程式が得られます。ここで $x$ は進行方向、$y$ は法線方向、$z$ は動径方向です。
重力項の線形化では、$\bm{R}_c = \bm{R}_t + \bm{\rho}$ を代入して
$$ \frac{\mu}{|\bm{R}_t + \bm{\rho}|^3}(\bm{R}_t + \bm{\rho}) $$
を $|\bm{\rho}|/R_t$ の1次までテイラー展開します。$z$ 方向（動径方向）には重力の変化率（潮汐力）が現れ、$x$ 方向（進行方向）にはコリオリ力と潮汐力が結合して現れます。
この線形化の計算は煩雑ですが、結果は驚くほどすっきりした形になります。
CWH方程式（Hill方程式）
以上の導出を経て、Clohessy-Wiltshire-Hill（CWH）方程式が得られます。
$$ \begin{equation} \ddot{x} – 2n\dot{z} = f_x \end{equation} $$
$$ \begin{equation} \ddot{y} + n^2 y = f_y \end{equation} $$
$$ \begin{equation} \ddot{z} + 2n\dot{x} – 3n^2 z = f_z \end{equation} $$
ここで $(f_x, f_y, f_z)$ はチェイサーのスラスタによる制御加速度です。
この3本の方程式を眺めると、いくつかの重要な特徴が見えます。
- $y$ 方向（法線方向）は独立: $y$ の方程式は単独で単振動（振動数 $n$）を表します。軌道面外の運動は軌道面内の運動と切り離して考えることができます
- $x$-$z$ は結合: $\dot{z}$ が $x$ の方程式に、$\dot{x}$ が $z$ の方程式に現れています。これはコリオリ力による結合で、進行方向と動径方向の運動が互いに影響し合うことを意味します
- $z$ に $-3n^2 z$ の項: 動径方向には軌道高度に比例した潮汐力が働きます。ターゲットより高い位置（$z > 0$）にいるチェイサーは地球から遠いため重力が弱く、遠心力が勝って外側に押し出されます
CWH方程式は一見シンプルですが、ランデブー・近傍接近の誘導制御の基盤となる極めて強力な方程式です。次に、この方程式の解析解を求め、その物理的意味を探ります。
CWH方程式の解析解
自由漂流の一般解
スラスタを噴射しない場合（$f_x = f_y = f_z = 0$）のCWH方程式を解きましょう。まず $y$ 方向は単純な単振動です。
$$ y(t) = y_0 \cos(nt) + \frac{\dot{y}_0}{n} \sin(nt) $$
$x$-$z$ 結合系はやや複雑ですが、定数係数の線形常微分方程式なので、特性方程式を解くことで一般解が得られます。初期条件 $(x_0, z_0, \dot{x}_0, \dot{z}_0)$ を用いて書き下すと
$z$ 方向（動径方向）:
$$ \begin{equation} z(t) = \left(\frac{4z_0}{1} – \frac{2\dot{x}_0}{n}\right)\cos(nt) + \frac{\dot{z}_0}{n}\sin(nt) + \left(\frac{2\dot{x}_0}{n} – 3z_0\right) \end{equation} $$
上式を整理すると
$$ z(t) = \left(4z_0 – \frac{2\dot{x}_0}{n}\right)\cos(nt) + \frac{\dot{z}_0}{n}\sin(nt) + \frac{2\dot{x}_0}{n} – 3z_0 $$
$x$ 方向（進行方向）:
$$ x(t) = x_0 + \left(-6nz_0 + 2\dot{x}_0 \cdot \frac{2}{1}\right) \cdot \frac{1}{n} \cdot (1 – \cos(nt)) + \frac{\dot{z}_0}{n}(1 – \cos(nt)) \cdot 0 $$
正確に書き下すと
$$ \begin{align} x(t) &= x_0 + \frac{2\dot{z}_0}{n}(1 – \cos(nt)) + \left(-3\dot{x}_0 + 6nz_0\right) \cdot t \\ &\quad + \left(\frac{2\dot{x}_0}{n} – 4z_0\right) \cdot 0 + \frac{2}{n}\left(\dot{x}_0 \sin(nt)\right) \cdot 0 \end{align} $$
ここでは読みやすさのため、整理した形で改めて示します。CWH方程式の自由漂流解（$\bm{f} = \bm{0}$）の完全な形は以下の通りです。
$$ \begin{align} x(t) &= \left(6z_0 + \frac{4\dot{x}_0}{n}\right)\sin(nt) – \frac{2\dot{z}_0}{n}\cos(nt) + \left(-6nz_0 + 3\dot{x}_0 \cdot 0\right) \cdot 0 \end{align} $$
混乱を避けるため、CWH方程式の解析解を状態遷移行列（STM）の形で整理しましょう。これが最も見通しの良い表現です。
状態遷移行列による表現
状態ベクトル $\bm{X}(t) = (x, y, z, \dot{x}, \dot{y}, \dot{z})^\top$ に対して、自由漂流解は次の状態遷移行列で表されます。
$$ \bm{X}(t) = \bm{\Phi}(t) \bm{X}_0 $$
ここで $\bm{\Phi}(t)$ は $6 \times 6$ の行列で、$x$-$z$ 結合部分と $y$ 独立部分に分けて書くと
$x$-$z$ 面内成分の状態遷移行列（$4 \times 4$）:
$$ \bm{\Phi}_{xz}(t) = \begin{pmatrix} 1 & 6(nt – \sin nt) & \frac{4\sin nt}{n} – 3t & \frac{2(1 – \cos nt)}{n} \\ 0 & 4 – 3\cos nt & \frac{2(\cos nt – 1)}{n} & \frac{\sin nt}{n} \\ 0 & 6n(\cos nt – 1) & 4\cos nt – 3 & 2\sin nt \\ 0 & 3n\sin nt & -2\sin nt & \cos nt \end{pmatrix} $$
ここで行ベクトルは $(x, z, \dot{x}, \dot{z})$ の順に対応し、列ベクトルは初期値 $(x_0, z_0, \dot{x}_0, \dot{z}_0)$ の順です。
$y$ 方向の状態遷移行列（$2 \times 2$）:
$$ \bm{\Phi}_y(t) = \begin{pmatrix} \cos nt & \frac{\sin nt}{n} \\ -n\sin nt & \cos nt \end{pmatrix} $$
この状態遷移行列は、初期条件を与えるだけでチェイサーの任意時刻の状態を計算できるため、極めて実用的です。
各成分を改めて明示的に書くと、$\tau = nt$ として
$$ x(t) = x_0 + 6z_0(\tau – \sin\tau) + \frac{4\dot{x}_0}{n}\sin\tau – \frac{2\dot{z}_0}{n}(1 – \cos\tau) – 3\dot{x}_0 t $$
$$ z(t) = z_0(4 – 3\cos\tau) + \frac{2\dot{x}_0}{n}(\cos\tau – 1) + \frac{\dot{z}_0}{n}\sin\tau $$
$$ y(t) = y_0\cos\tau + \frac{\dot{y}_0}{n}\sin\tau $$
この解は見やすく、物理的意味を読み取りやすい形です。次にこの解から読み取れる重要な物理を整理します。
CWH解の物理的意味
永年項（ドリフト）— 放っておくと離れていく
$x(t)$ の式に注目すると、$\sin\tau$ や $\cos\tau$ のような振動項に加えて、$t$ に比例する永年項（secular term）が現れています。
$$ x_{\text{secular}} = 6z_0(\tau – \sin\tau) – 3\dot{x}_0 t = (6nz_0 – 3\dot{x}_0)t + \text{振動項} $$
$z(t)$ にも定数オフセット項 $-3z_0 + \frac{2\dot{x}_0}{n}$ が現れます（$\cos\tau$ の係数を参照）。
永年項は時間とともに際限なく増大するため、チェイサーはターゲットから果てしなく離れていきます。高速道路のアナロジーで言えば、わずかに速度が異なる2台の車が時間とともにどんどん離れていくのと同じです。
永年項がゼロになる条件は
$$ \dot{x}_0 = 2nz_0 $$
です。この条件が満たされると、$x(t)$ と $z(t)$ は純粋に振動的になり、チェイサーはターゲットの周囲を周回し続けます。
安全楕円 — ターゲットの周りを「公転」する
永年項がゼロの条件 $\dot{x}_0 = 2nz_0$ が満たされた場合の相対運動を考えます。簡単のため、$y_0 = \dot{y}_0 = 0$（軌道面内運動のみ）として $z_0 = z_a$, $\dot{z}_0 = 0$, $\dot{x}_0 = 2nz_a$ とすると
$$ x(t) = -2z_a \sin(nt) $$
$$ z(t) = z_a \cos(nt) $$
これは楕円方程式 $\left(\frac{x}{2z_a}\right)^2 + \left(\frac{z}{z_a}\right)^2 = 1$ を満たします。つまり、チェイサーはターゲットの周りを $x$方向に$2z_a$, $z$方向に$z_a$ の楕円軌道で周回します。
この楕円は2:1の楕円（along-track方向が動径方向の2倍）という特徴的な形状を持ちます。これはCWH方程式に特有の結果であり、「ランデブーの楕円」や「フットボール軌道」とも呼ばれます。
直感的に理解してみましょう。チェイサーがターゲットより高い位置（$z > 0$）にいるとき、軌道半径が大きいので軌道速度が遅くなり、ターゲットに対して後方（$-x$ 方向）にドリフトします。逆に低い位置（$z < 0$）では速度が速くなり前方に移動します。この「高度が変わると進行方向の速度が変わる」効果が、2:1の楕円を生み出すのです。
安全楕円の実運用上の意義
安全楕円は受動安全（passively safe）な軌道です。スラスタが故障して噴射できなくなっても、チェイサーはターゲットに衝突せずにその周囲を回り続けます。実際のランデブーミッションでは、万一の場合に備えて、意図的に安全楕円上にチェイサーを配置してから最終接近を開始するという運用が行われます。
ESAの無人補給機ATV（Automated Transfer Vehicle）は、ISSへの接近中に安全楕円上の保持点（Hold Point）を複数設けて、各段階での健全性確認（Go/No-Go判定）を行いました。
CWH方程式の解とその物理的意味がわかったところで、次はランデブーの全体的な運用手順を見ていきましょう。実際のミッションではどのようなフェーズを経てドッキングに至るのでしょうか。
ランデブーの各フェーズ
全体像
ランデブーは一般に、ターゲットまでの距離に応じて複数のフェーズに分けられます。各フェーズで使うセンサ、制御手法、精度要求が異なります。
| フェーズ | 距離 | 主なセンサ | 制御手法 |
|---|---|---|---|
| 遠距離接近（Far-Range） | 数百km → 数十km | 地上追跡（レーダー、テレメトリ）、GPS | ホーマン遷移ベース |
| 中距離接近（Mid-Range） | 数十km → 数百m | GPS相対航法、レーダー | CWH ベースの二体力マニューバ |
| 近傍運用（ProxOps） | 数百m → 数十m | 光学センサ（LiDAR, カメラ） | フィードバック誘導 |
| 最終接近（Final Approach） | 数十m → 0 | LiDAR、ビジュアルサーボ | 精密直線接近 |
遠距離接近フェーズ
チェイサーが最初のマニューバを行い、ターゲットの近傍（通常数十km以内）まで軌道を合わせるフェーズです。ここではCWH方程式の適用範囲を超える大きな軌道変更が必要なため、ホーマン遷移やランバート問題の解法が使われます。
このフェーズの精度は比較的粗くてよく（km単位）、地上局からのレーダー追跡やGPSによる絶対航法で十分です。
中距離接近フェーズ
ターゲットの数十km以内に到達すると、CWH方程式が適用できる範囲に入ります。ここから相対航法（GPS差分測位やレーダーレンジ測定）を用いて、相対位置・速度を高精度に推定します。
この距離では、V-bar（進行方向軸）上の保持点に向かう一連のホップマニューバ（二体力インパルス）が典型的です。保持点に到達するたびにGo/No-Go判定を行い、安全を確認してから次の段階に進みます。
近傍運用フェーズ
数百m以内になると、ターゲットの詳細構造を光学センサ（LiDAR、赤外線カメラ、可視カメラ）で識別できるようになります。特に非協力ターゲット（故障衛星やデブリ）に対しては、この段階で初めてターゲットの姿勢状態や表面形状を把握します。
近傍運用では、CWHベースの軌道計画に加えて、リアルタイムのフィードバック制御が重要になります。センサ情報を用いて相対状態を連続的に推定し、偏差を修正します。
最終接近フェーズ
数十m以内の最終接近では、ミリメートル単位の精度が求められます。LiDARによる距離測定とカメラによる姿勢推定を融合（センサフュージョン）し、ドッキング機構やロボットアームの可動範囲に正確に誘導します。
最終接近は通常、直線的な経路（V-barまたはR-barに沿った経路）で行われ、接近速度は距離に比例して減少するように設計されます。たとえば、ISSへの接近では「距離の1%の速度」というルール（100 mで1 m/s、10 mで0.1 m/s）が安全基準として使われることがあります。
各フェーズの概要がわかったところで、次は最終接近の2つの代表的な手法 — V-barアプローチとR-barアプローチ — を詳しく比較します。
V-barアプローチ vs R-barアプローチ
V-barアプローチ（速度方向からの接近）
V-barアプローチは、ターゲットの進行方向軸（V-bar）に沿って接近する方法です。チェイサーはターゲットの前方または後方に位置し、V-bar上をまっすぐ近づきます。
具体的には、チェイサーはターゲットよりわずかに低い（または高い）軌道にいて、V-bar上に沿う一連のホップマニューバで段階的に距離を詰めます。各ホップの終端でV-bar上の保持点に静止し、次のマニューバを計画します。
利点: – 接近経路が直感的でシンプル – 地上管制からの監視が容易 – 相対速度が小さい状態を維持しやすい
欠点: – デルタV消費量がR-barアプローチより大きい場合がある – V-bar上で静止を維持するには定期的な軌道維持噴射が必要（CWH方程式から、V-bar上の静止は不安定であるため）
R-barアプローチ（動径方向からの接近）
R-barアプローチは、ターゲットの動径方向軸（R-bar）に沿って、下方から接近する方法です。チェイサーはターゲットの真下（地球側）に位置し、垂直に上昇してドッキングします。
利点: – 受動安全性が高い: R-bar上で下方にいるチェイサーは、スラスタ故障時にもターゲットに衝突する方向には自然にドリフトしない（重力勾配がチェイサーを下方に保つ効果がある） – R-bar上のある種の静止点は自然に安定
欠点: – 地球の重力勾配に逆らって上昇するため、連続的なスラスタ噴射（またはホップの繰り返し）が必要 – 接近経路の設計がV-barより複雑 – 大型のプルーム・インピンジメント（排気噴流のターゲットへの衝突）リスクがある
どちらを選ぶか？
ISS向けの補給ミッションでは、V-barアプローチが多く採用されてきました。スペースXのドラゴン宇宙船やJAXAのHTV（こうのとり）はV-bar方向からISSに接近し、ISSのロボットアーム（SSRMS、カナダアーム2）で捕獲されました。
一方、ESAのATVやロシアのプログレスは、条件によってR-barアプローチも使用しています。特に安全性が最優先のミッション（有人宇宙ステーションへの接近）では、R-barの受動安全性が重視されます。
OOS（軌道上サービス）ミッションでは、ターゲットの状態（姿勢回転の有無、ドッキング面の方向）によってアプローチ方向を柔軟に選択する必要があり、V-bar/R-barの二択にとどまらない、より複雑な接近軌道が設計されます。
接近経路の戦略が理解できたところで、次はその経路を実現するための具体的な誘導法則を見ていきましょう。
誘導法則の設計
誘導制御の役割
ランデブーにおける誘導（guidance）とは、現在の状態から目標状態に至る経路と噴射タイミング・大きさを計算する機能です。制御（control）はその計算結果に従ってスラスタを実際に動作させる機能で、航法（navigation）はセンサから現在の相対状態を推定する機能です。これらを合わせてGNC（Guidance, Navigation, and Control）と呼びます。
ここでは誘導法則、すなわち「いつ、どの方向に、どれだけ噴射するか」を決めるアルゴリズムを3種類紹介します。
二体力マニューバ（Two-Impulse Maneuver）
最もシンプルな誘導法則は、CWH方程式の状態遷移行列を用いた二体力マニューバです。初期時刻に1回、到着時刻にもう1回、計2回のインパルス噴射でチェイサーを目標位置に移動させます。
時刻 $t_0$ の状態 $(x_0, z_0)$ から時刻 $t_f$ の目標位置 $(x_f, z_f)$ に到達するには、状態遷移行列を用いて必要な初速度を逆算します。
状態遷移行列を位置部分と速度部分に分割して
$$ \begin{pmatrix} \bm{r}(t_f) \\ \bm{v}(t_f) \end{pmatrix} = \begin{pmatrix} \bm{\Phi}_{rr} & \bm{\Phi}_{rv} \\ \bm{\Phi}_{vr} & \bm{\Phi}_{vv} \end{pmatrix} \begin{pmatrix} \bm{r}_0 \\ \bm{v}_0^+ \end{pmatrix} $$
と書くと、目標位置に到達するための噴射後速度 $\bm{v}_0^+$ は
$$ \bm{v}_0^+ = \bm{\Phi}_{rv}^{-1} (\bm{r}_f – \bm{\Phi}_{rr} \bm{r}_0) $$
で求められます。$\bm{\Phi}_{rv}$ が正則であることが必要ですが、転送時間が軌道周期の整数倍でなければ通常成り立ちます。
第1インパルス（出発時）のデルタV:
$$ \Delta \bm{v}_1 = \bm{v}_0^+ – \bm{v}_0^- $$
第2インパルス（到着時）のデルタV:
$$ \Delta \bm{v}_2 = \bm{v}_f^{\text{desired}} – \bm{v}_f^{\text{arrival}} $$
ここで $\bm{v}_0^-$ は噴射前の速度、$\bm{v}_f^{\text{desired}}$ は目標点での希望速度（保持点に静止するなら $\bm{0}$）、$\bm{v}_f^{\text{arrival}}$ は自由漂流で到着したときの速度です。
CWヒル誘導法則（連続フィードバック）
二体力マニューバはオープンループ（噴射計画を事前に計算して実行）ですが、外乱や航法誤差がある現実のミッションでは、フィードバック誘導が必要です。
CWH方程式を状態空間形式で書くと
$$ \dot{\bm{X}} = \bm{A}\bm{X} + \bm{B}\bm{u} $$
ここで $\bm{X} = (x, y, z, \dot{x}, \dot{y}, \dot{z})^\top$, $\bm{u} = (f_x, f_y, f_z)^\top$ です。
$$ \bm{A} = \begin{pmatrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 2n \\ 0 & -n^2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 3n^2 & -2n & 0 & 0 \end{pmatrix}, \quad \bm{B} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} $$
この系は可制御（controllable）であることが確認できます（可制御性行列のランクが6）。したがって、LQR（Linear Quadratic Regulator）などの最適制御手法でフィードバックゲインを設計できます。
LQR制御則では、制御入力を
$$ \bm{u} = -\bm{K}\bm{X} $$
と設計します。ゲイン行列 $\bm{K}$ はリカッチ方程式
$$ \bm{A}^\top \bm{P} + \bm{P}\bm{A} – \bm{P}\bm{B}\bm{R}^{-1}\bm{B}^\top \bm{P} + \bm{Q} = \bm{0} $$
の解 $\bm{P}$ から $\bm{K} = \bm{R}^{-1}\bm{B}^\top \bm{P}$ として得られます。$\bm{Q}$ は状態の重み行列（位置精度の重要度）、$\bm{R}$ は制御入力の重み行列（燃料消費の重要度）です。
$\bm{Q}$ を大きくすると位置精度を優先し、$\bm{R}$ を大きくすると燃料節約を優先します。このトレードオフはミッション要求に応じて設計者が調整します。
グライドスロープ誘導
最終接近フェーズで広く使われるのがグライドスロープ（glideslope）誘導です。これは接近速度を距離に比例させる方法で、距離が縮まるにつれて自然に減速します。
$$ v_{\text{approach}} = k \cdot d $$
ここで $d$ はターゲットまでの距離、$k$ はグライドスロープ係数（単位: $1/\text{s}$）です。$k = 0.01 \, \text{s}^{-1}$ とすると、100 mで1 m/s、10 mで0.1 m/s という直感的な減速プロファイルが得られます。
グライドスロープの利点は安全性の保証です。グライドスロープ係数を保守的に設定すれば、万一センサが故障して誘導を失っても、その時点の速度が距離に比例しているため、衝突エネルギーが制限されます。
ここまでで誘導法則の基本的な選択肢を見てきました。実際のミッションでは、これらを組み合わせて使います。遠距離では二体力マニューバ、中距離ではLQRフィードバック、最終接近ではグライドスロープ、というように距離に応じて切り替えるのが一般的です。次に、ランデブーにおいて極めて重要な「安全性」の考え方を見ていきましょう。
安全性の考慮 — 受動安全軌道と衝突回避
なぜ安全性が最優先なのか
宇宙でのランデブーは、失敗すれば衝突という取り返しのつかない事態を招きます。特にISSのような有人施設へのランデブーでは、乗員の生命がかかっています。そのため、ランデブーの設計では安全性が性能より優先されます。
1997年、ロシアのプログレス補給船がミール宇宙ステーションに手動ドッキング中に衝突し、スペクトルモジュールに穴が開いて減圧が発生するという深刻な事故がありました。この事故は、ランデブーの安全設計の重要性を改めて認識させました。
受動安全軌道（Passively Safe Trajectory）
受動安全とは、「スラスタが完全に停止しても、一定時間内にターゲットに衝突しない」という性質です。CWH方程式の自由漂流解で、ターゲットとの距離が常にゼロにならないことが保証される初期条件を選びます。
先に見た安全楕円（$\dot{x}_0 = 2nz_0$ の条件）は受動安全の典型例です。スラスタが停止しても、チェイサーはターゲットの周りを2:1の楕円軌道で周回し続けます。
より一般的には、軌道計画の各段階で「ここでスラスタが止まったら何が起こるか」を事前にシミュレーションし、一定時間（たとえば2軌道周期）以内にターゲットのKeep-Out Sphere（接近禁止球、通常200 m程度）に侵入しないことを確認します。
CAM（Collision Avoidance Maneuver）
万一、異常が検知された場合に即座に実行する退避マニューバがCAM（衝突回避マニューバ）です。CAMは事前に計算されており、1回のインパルス噴射でチェイサーをターゲットから安全な方向に離脱させます。
典型的なCAMの設計指針: – R-bar方向への噴射: 下方（地球側）に噴射してチェイサーの軌道を下げ、ターゲットから離脱する – V-bar方向への逆噴射: 後方に噴射して減速し、低い軌道に遷移してターゲットから遠ざかる – 実行条件の自動化: 航法系やスラスタ系の異常を検知したら、地上の判断を待たずに自動でCAMを実行
セーフモードと保持点の設計
実際のミッションでは、接近経路上に複数の保持点（Hold Point, HP）を設定します。保持点はCWH解析から安定に留まれる位置（またはV-bar/R-bar上の特定の距離）に配置され、各保持点で以下の確認を行います。
- GNC系の健全性確認
- 地上管制との通信確認
- ターゲット側のGo/No-Go判定
- 次のフェーズの誘導パラメータのアップリンク
保持点での滞在中、チェイサーは安全楕円上の周回運動を行うか、V-bar上で軌道維持噴射を行って位置を保持します。
ここまでの理論を踏まえて、いよいよPythonでCWH方程式のシミュレーションとランデブー軌道の可視化を行いましょう。
PythonでCWH方程式をシミュレーションする
CWH方程式の数値積分
まず、CWH方程式を数値的に解くコードを作成します。SciPyの solve_ivp を使って、状態ベクトル $\bm{X} = (x, y, z, \dot{x}, \dot{y}, \dot{z})$ の時間発展を計算します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# --- 定数 ---
mu_earth = 3.986004418e14 # 地球の重力パラメータ [m^3/s^2]
R_earth = 6.371e6 # 地球の半径 [m]
h_ISS = 400e3 # ISS高度 [m]
a_target = R_earth + h_ISS # ターゲットの円軌道半径 [m]
n = np.sqrt(mu_earth / a_target**3) # 平均運動 [rad/s]
T_orb = 2 * np.pi / n # 軌道周期 [s]
print(f"平均運動 n = {n:.6e} rad/s")
print(f"軌道周期 T = {T_orb:.1f} s ({T_orb/60:.1f} min)")
def cwh_equations(t, state, n, u_func=None):
"""CWH方程式の右辺
state: [x, y, z, vx, vy, vz]
n: 平均運動 [rad/s]
u_func: 制御入力の関数 u_func(t, state) -> [fx, fy, fz]
"""
x, y, z, vx, vy, vz = state
if u_func is not None:
fx, fy, fz = u_func(t, state)
else:
fx, fy, fz = 0.0, 0.0, 0.0
ax = 2 * n * vz + fx
ay = -n**2 * y + fy
az = -2 * n * vx + 3 * n**2 * z + fz
return [vx, vy, vz, ax, ay, az]
このコードでは、CWH方程式を関数として定義しています。u_func 引数で制御入力を外部から注入できるようにしており、自由漂流（制御なし）と制御付きの両方のシミュレーションに対応します。
自由漂流のシミュレーション
まず、スラスタを使わない自由漂流の挙動を確認します。初期条件を変えたときの相対軌道がどう変化するかを見てみましょう。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
mu_earth = 3.986004418e14
R_earth = 6.371e6
h_ISS = 400e3
a_target = R_earth + h_ISS
n = np.sqrt(mu_earth / a_target**3)
T_orb = 2 * np.pi / n
def cwh_equations(t, state, n, u_func=None):
x, y, z, vx, vy, vz = state
if u_func is not None:
fx, fy, fz = u_func(t, state)
else:
fx, fy, fz = 0.0, 0.0, 0.0
ax = 2 * n * vz + fx
ay = -n**2 * y + fy
az = -2 * n * vx + 3 * n**2 * z + fz
return [vx, vy, vz, ax, ay, az]
# --- ケース1: 永年項あり（ドリフトする場合）---
# z0 = 100 m, vx0 = 0 (安全楕円条件を満たさない)
state0_drift = [0, 0, 100, 0, 0, 0] # [x, y, z, vx, vy, vz] [m, m/s]
# --- ケース2: 安全楕円（永年項なし）---
z0 = 100 # [m]
vx0_safe = 2 * n * z0 # 安全楕円条件
state0_safe = [0, 0, z0, vx0_safe, 0, 0]
# --- ケース3: 任意の初期速度 ---
state0_general = [0, 0, 100, 0.05, 0, 0.02]
t_span = (0, 3 * T_orb)
t_eval = np.linspace(0, 3 * T_orb, 5000)
# 数値積分
sol_drift = solve_ivp(cwh_equations, t_span, state0_drift,
args=(n,), t_eval=t_eval, rtol=1e-10, atol=1e-12)
sol_safe = solve_ivp(cwh_equations, t_span, state0_safe,
args=(n,), t_eval=t_eval, rtol=1e-10, atol=1e-12)
sol_general = solve_ivp(cwh_equations, t_span, state0_general,
args=(n,), t_eval=t_eval, rtol=1e-10, atol=1e-12)
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# ケース1: ドリフト
axes[0].plot(sol_drift.y[0], sol_drift.y[2], 'b-', linewidth=0.8)
axes[0].plot(sol_drift.y[0][0], sol_drift.y[2][0], 'go', markersize=8, label='Start')
axes[0].plot(0, 0, 'r*', markersize=12, label='Target')
axes[0].set_xlabel('x (along-track) [m]')
axes[0].set_ylabel('z (radial) [m]')
axes[0].set_title('Case 1: Free Drift (secular term present)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_aspect('equal')
# ケース2: 安全楕円
axes[1].plot(sol_safe.y[0], sol_safe.y[2], 'b-', linewidth=0.8)
axes[1].plot(sol_safe.y[0][0], sol_safe.y[2][0], 'go', markersize=8, label='Start')
axes[1].plot(0, 0, 'r*', markersize=12, label='Target')
axes[1].set_xlabel('x (along-track) [m]')
axes[1].set_ylabel('z (radial) [m]')
axes[1].set_title('Case 2: Safety Ellipse (no drift)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_aspect('equal')
# ケース3: 一般
axes[2].plot(sol_general.y[0], sol_general.y[2], 'b-', linewidth=0.8)
axes[2].plot(sol_general.y[0][0], sol_general.y[2][0], 'go', markersize=8, label='Start')
axes[2].plot(0, 0, 'r*', markersize=12, label='Target')
axes[2].set_xlabel('x (along-track) [m]')
axes[2].set_ylabel('z (radial) [m]')
axes[2].set_title('Case 3: General Free Drift')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('cwh_free_drift.png', dpi=150, bbox_inches='tight')
plt.show()
3つのケースの結果を比較すると、以下の重要な違いが読み取れます。
-
ケース1（永年項あり）: 安全楕円条件を満たさない初期条件（$\dot{x}_0 = 0 \neq 2nz_0$）では、チェイサーは $x$ 方向に一方的にドリフトしていきます。3軌道周期の間に数百メートルも流されており、放置すれば際限なく離れていくことがわかります
-
ケース2（安全楕円）: $\dot{x}_0 = 2nz_0$ の条件を満たすと、ドリフトが消えて閉じた楕円軌道になります。$x$ 方向の振幅が $z$ 方向の2倍（2:1の楕円）であることも確認できます。これが受動安全な保持軌道です
-
ケース3（一般初期条件）: 任意の初期速度を持つ場合は、ドリフトに楕円運動が重畳した複雑な軌道（サイクロイド的な曲線）になります。実際のランデブーでは、このような「ドリフト + 振動」の混合した運動を制御して目標に誘導することになります
二体力マニューバのシミュレーション
次に、V-bar上の保持点からターゲット原点に向かう二体力マニューバをシミュレーションします。CWH状態遷移行列を使って必要なデルタVを計算し、その結果を可視化します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
mu_earth = 3.986004418e14
R_earth = 6.371e6
h_ISS = 400e3
a_target = R_earth + h_ISS
n = np.sqrt(mu_earth / a_target**3)
T_orb = 2 * np.pi / n
def cwh_stm_xz(t, n):
"""CWH状態遷移行列のx-z面内成分 (4x4)
状態: [x, z, vx, vz]
"""
nt = n * t
s = np.sin(nt)
c = np.cos(nt)
Phi = np.array([
[1, 6*(nt - s), (4*s)/n - 3*t, 2*(1 - c)/n],
[0, 4 - 3*c, 2*(c - 1)/n, s/n ],
[0, 6*n*(c - 1), 4*c - 3, 2*s ],
[0, 3*n*s, -2*s, c ]
])
return Phi
# --- 二体力マニューバの計算 ---
# V-bar上 x = -500 m (ターゲット後方) から原点へ
r0 = np.array([-500.0, 0.0]) # [x0, z0] [m]
rf = np.array([0.0, 0.0]) # 目標位置（原点）
v0_minus = np.array([0.0, 0.0]) # 初期速度（V-bar上で静止）
# 転送時間: 軌道周期の半分
t_transfer = 0.5 * T_orb
Phi = cwh_stm_xz(t_transfer, n)
Phi_rr = Phi[0:2, 0:2] # 位置→位置
Phi_rv = Phi[0:2, 2:4] # 速度→位置
Phi_vr = Phi[2:4, 0:2] # 位置→速度
Phi_vv = Phi[2:4, 2:4] # 速度→速度
# 噴射後速度の計算: rf = Phi_rr * r0 + Phi_rv * v0+
# v0+ = Phi_rv^{-1} * (rf - Phi_rr * r0)
v0_plus = np.linalg.solve(Phi_rv, rf - Phi_rr @ r0)
# 到着時速度
vf_arrival = Phi_vr @ r0 + Phi_vv @ v0_plus
# デルタV
dv1 = v0_plus - v0_minus # 第1インパルス
dv2 = -vf_arrival # 第2インパルス（原点で静止するため）
print(f"転送時間: {t_transfer:.1f} s ({t_transfer/60:.1f} min)")
print(f"第1インパルス dv1 = ({dv1[0]:.4f}, {dv1[1]:.4f}) m/s, |dv1| = {np.linalg.norm(dv1):.4f} m/s")
print(f"第2インパルス dv2 = ({dv2[0]:.4f}, {dv2[1]:.4f}) m/s, |dv2| = {np.linalg.norm(dv2):.4f} m/s")
print(f"合計デルタV = {np.linalg.norm(dv1) + np.linalg.norm(dv2):.4f} m/s")
続けて、この二体力マニューバの軌道を数値積分で可視化します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
mu_earth = 3.986004418e14
R_earth = 6.371e6
h_ISS = 400e3
a_target = R_earth + h_ISS
n = np.sqrt(mu_earth / a_target**3)
T_orb = 2 * np.pi / n
def cwh_equations(t, state, n, u_func=None):
x, y, z, vx, vy, vz = state
if u_func is not None:
fx, fy, fz = u_func(t, state)
else:
fx, fy, fz = 0.0, 0.0, 0.0
ax = 2 * n * vz + fx
ay = -n**2 * y + fy
az = -2 * n * vx + 3 * n**2 * z + fz
return [vx, vy, vz, ax, ay, az]
def cwh_stm_xz(t, n):
nt = n * t
s = np.sin(nt)
c = np.cos(nt)
Phi = np.array([
[1, 6*(nt - s), (4*s)/n - 3*t, 2*(1 - c)/n],
[0, 4 - 3*c, 2*(c - 1)/n, s/n ],
[0, 6*n*(c - 1), 4*c - 3, 2*s ],
[0, 3*n*s, -2*s, c ]
])
return Phi
# 二体力マニューバの計算（再計算）
r0 = np.array([-500.0, 0.0])
rf = np.array([0.0, 0.0])
v0_minus = np.array([0.0, 0.0])
t_transfer = 0.5 * T_orb
Phi = cwh_stm_xz(t_transfer, n)
Phi_rr = Phi[0:2, 0:2]
Phi_rv = Phi[0:2, 2:4]
v0_plus = np.linalg.solve(Phi_rv, rf - Phi_rr @ r0)
# 噴射後の初期状態
state0_maneuver = [-500.0, 0.0, 0.0, v0_plus[0], 0.0, v0_plus[1]]
t_eval_man = np.linspace(0, t_transfer, 2000)
sol_maneuver = solve_ivp(
cwh_equations, (0, t_transfer), state0_maneuver,
args=(n,), t_eval=t_eval_man, rtol=1e-10, atol=1e-12
)
# 自由漂流（噴射しなかった場合）の比較
state0_nomaneuver = [-500.0, 0.0, 0.0, 0.0, 0.0, 0.0]
sol_nomaneuver = solve_ivp(
cwh_equations, (0, t_transfer), state0_nomaneuver,
args=(n,), t_eval=t_eval_man, rtol=1e-10, atol=1e-12
)
fig, ax = plt.subplots(figsize=(10, 6))
# 二体力マニューバの軌道
ax.plot(sol_maneuver.y[0], sol_maneuver.y[2], 'b-', linewidth=1.5,
label='Two-impulse transfer')
# 自由漂流（比較）
ax.plot(sol_nomaneuver.y[0], sol_nomaneuver.y[2], 'gray', linewidth=1,
linestyle='--', alpha=0.6, label='Free drift (no maneuver)')
# 始点・終点
ax.plot(-500, 0, 'go', markersize=10, zorder=5, label='Start (V-bar, -500 m)')
ax.plot(0, 0, 'r*', markersize=15, zorder=5, label='Target (origin)')
ax.plot(sol_maneuver.y[0][-1], sol_maneuver.y[2][-1], 'bs',
markersize=8, zorder=5, label='Arrival')
# デルタVの矢印
scale = 3000 # 矢印の表示スケール
ax.annotate('', xy=(-500 + v0_plus[0]*scale, v0_plus[1]*scale),
xytext=(-500, 0),
arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax.text(-500 + v0_plus[0]*scale*0.5, v0_plus[1]*scale*0.5 + 10,
r'$\Delta v_1$', fontsize=12, color='red')
ax.set_xlabel('x (along-track) [m]', fontsize=12)
ax.set_ylabel('z (radial) [m]', fontsize=12)
ax.set_title('Two-Impulse Rendezvous Maneuver (V-bar approach, T/2 transfer)', fontsize=13)
ax.legend(fontsize=10, loc='upper left')
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='k', linewidth=0.5)
ax.axvline(x=0, color='k', linewidth=0.5)
plt.tight_layout()
plt.savefig('two_impulse_maneuver.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフからいくつかの重要な知見が得られます。
-
二体力マニューバの軌道は弧を描く: ターゲット後方500 mの V-bar上から出発したチェイサーは、一度 $z$（動径）方向に逸れてからターゲットに到達します。直線的に近づくのではなく、軌道力学的な弧を描く点が宇宙ならではの特徴です
-
自由漂流との比較: 灰色の破線で示した自由漂流軌道では、チェイサーはターゲットに到達できません。V-bar上で静止しているチェイサーは、CWH方程式の永年項により徐々にドリフトしていきます。二体力マニューバの噴射が本質的に必要であることが明確にわかります
-
デルタVは小さい: 500 mの距離を半軌道周期（約46分）で移動するのに必要なデルタVは非常に小さく（m/s以下のオーダー）、近傍運用が燃料効率に優れていることがわかります
LQRフィードバック誘導のシミュレーション
最後に、LQR制御を用いたフィードバック誘導のシミュレーションを行います。外乱やセンサノイズがある現実的な状況での自律接近を模擬します。
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import solve_continuous_are
import matplotlib.pyplot as plt
mu_earth = 3.986004418e14
R_earth = 6.371e6
h_ISS = 400e3
a_target = R_earth + h_ISS
n = np.sqrt(mu_earth / a_target**3)
T_orb = 2 * np.pi / n
# --- CWH状態空間モデル ---
A = np.array([
[0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 2*n],
[0, -n**2, 0, 0, 0, 0],
[0, 0, 3*n**2, -2*n, 0, 0]
])
B = np.array([
[0, 0, 0],
[0, 0, 0],
[0, 0, 0],
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]
])
# --- LQR設計 ---
# Q: 状態の重み（位置精度を重視）
Q = np.diag([1.0, 1.0, 1.0, 0.1, 0.1, 0.1])
# R: 制御入力の重み（燃料消費のペナルティ）
R_weight = np.diag([100.0, 100.0, 100.0])
# 連続時間リカッチ方程式を解く
P = solve_continuous_are(A, B, Q, R_weight)
K = np.linalg.inv(R_weight) @ B.T @ P
print("LQRゲイン行列 K:")
print(f" K shape: {K.shape}")
print(f" K max element: {np.max(np.abs(K)):.6f}")
# --- LQRフィードバック付きCWH方程式 ---
def cwh_with_lqr(t, state, n, K):
x, y, z, vx, vy, vz = state
X = np.array(state)
u = -K @ X # LQR制御入力
# 制御入力の制限（現実的なスラスタ出力を想定）
u_max = 0.01 # 最大加速度 [m/s^2]（約1 mN / 0.1 kg級）
u = np.clip(u, -u_max, u_max)
fx, fy, fz = u
ax = 2 * n * vz + fx
ay = -n**2 * y + fy
az = -2 * n * vx + 3 * n**2 * z + fz
return [vx, vy, vz, ax, ay, az]
# --- シミュレーション ---
# 初期条件: ターゲットから500 m後方、100 m下方
state0_lqr = [-500.0, 50.0, -100.0, 0.0, 0.0, 0.0]
t_span_lqr = (0, 2 * T_orb)
t_eval_lqr = np.linspace(0, 2 * T_orb, 5000)
sol_lqr = solve_ivp(
cwh_with_lqr, t_span_lqr, state0_lqr,
args=(n, K), t_eval=t_eval_lqr, rtol=1e-10, atol=1e-12
)
# --- 制御入力の計算 ---
u_history = np.zeros((3, len(t_eval_lqr)))
for i in range(len(t_eval_lqr)):
X = sol_lqr.y[:, i]
u = -K @ X
u = np.clip(u, -0.01, 0.01)
u_history[:, i] = u
続けて可視化します。
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import solve_continuous_are
import matplotlib.pyplot as plt
mu_earth = 3.986004418e14
R_earth = 6.371e6
h_ISS = 400e3
a_target = R_earth + h_ISS
n = np.sqrt(mu_earth / a_target**3)
T_orb = 2 * np.pi / n
A = np.array([
[0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 2*n],
[0, -n**2, 0, 0, 0, 0],
[0, 0, 3*n**2, -2*n, 0, 0]
])
B = np.array([
[0, 0, 0],
[0, 0, 0],
[0, 0, 0],
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]
])
Q = np.diag([1.0, 1.0, 1.0, 0.1, 0.1, 0.1])
R_weight = np.diag([100.0, 100.0, 100.0])
P = solve_continuous_are(A, B, Q, R_weight)
K = np.linalg.inv(R_weight) @ B.T @ P
def cwh_with_lqr(t, state, n, K):
x, y, z, vx, vy, vz = state
X = np.array(state)
u = -K @ X
u_max = 0.01
u = np.clip(u, -u_max, u_max)
fx, fy, fz = u
ax = 2 * n * vz + fx
ay = -n**2 * y + fy
az = -2 * n * vx + 3 * n**2 * z + fz
return [vx, vy, vz, ax, ay, az]
state0_lqr = [-500.0, 50.0, -100.0, 0.0, 0.0, 0.0]
t_span_lqr = (0, 2 * T_orb)
t_eval_lqr = np.linspace(0, 2 * T_orb, 5000)
sol_lqr = solve_ivp(
cwh_with_lqr, t_span_lqr, state0_lqr,
args=(n, K), t_eval=t_eval_lqr, rtol=1e-10, atol=1e-12
)
u_history = np.zeros((3, len(t_eval_lqr)))
for i in range(len(t_eval_lqr)):
X = sol_lqr.y[:, i]
u = -K @ X
u = np.clip(u, -0.01, 0.01)
u_history[:, i] = u
# 累積デルタVの計算
dt = np.diff(t_eval_lqr)
delta_v_cumulative = np.zeros(len(t_eval_lqr))
for i in range(1, len(t_eval_lqr)):
delta_v_cumulative[i] = delta_v_cumulative[i-1] + np.linalg.norm(u_history[:, i]) * dt[i-1]
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) x-z平面の軌道
ax = axes[0, 0]
ax.plot(sol_lqr.y[0], sol_lqr.y[2], 'b-', linewidth=1)
ax.plot(sol_lqr.y[0][0], sol_lqr.y[2][0], 'go', markersize=10, label='Start')
ax.plot(0, 0, 'r*', markersize=15, label='Target')
ax.set_xlabel('x (along-track) [m]')
ax.set_ylabel('z (radial) [m]')
ax.set_title('(a) LQR-guided trajectory (x-z plane)')
ax.legend()
ax.grid(True, alpha=0.3)
# (b) 距離の時間変化
ax = axes[0, 1]
dist = np.sqrt(sol_lqr.y[0]**2 + sol_lqr.y[1]**2 + sol_lqr.y[2]**2)
ax.semilogy(t_eval_lqr / 60, dist, 'b-', linewidth=1)
ax.set_xlabel('Time [min]')
ax.set_ylabel('Distance to target [m]')
ax.set_title('(b) Distance convergence')
ax.grid(True, alpha=0.3)
# (c) 制御入力の時間変化
ax = axes[1, 0]
ax.plot(t_eval_lqr / 60, u_history[0] * 1e3, label=r'$f_x$', linewidth=0.8)
ax.plot(t_eval_lqr / 60, u_history[1] * 1e3, label=r'$f_y$', linewidth=0.8)
ax.plot(t_eval_lqr / 60, u_history[2] * 1e3, label=r'$f_z$', linewidth=0.8)
ax.set_xlabel('Time [min]')
ax.set_ylabel('Control acceleration [mm/s²]')
ax.set_title('(c) Control inputs')
ax.legend()
ax.grid(True, alpha=0.3)
# (d) 累積デルタV
ax = axes[1, 1]
ax.plot(t_eval_lqr / 60, delta_v_cumulative, 'b-', linewidth=1)
ax.set_xlabel('Time [min]')
ax.set_ylabel('Cumulative delta-V [m/s]')
ax.set_title('(d) Fuel consumption')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('lqr_rendezvous.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"\n最終距離: {dist[-1]:.4f} m")
print(f"累積デルタV: {delta_v_cumulative[-1]:.4f} m/s")
LQR制御によるランデブーシミュレーションから、以下の知見が得られます。
-
軌道（a）: チェイサーは初期位置（後方500 m、下方100 m、法線方向50 m）から滑らかな曲線を描いてターゲット原点に収束しています。CWHの力学特性（コリオリ力、潮汐力）によって直線的な経路にはなりませんが、LQR制御が適切に補正して収束を実現しています
-
距離の収束（b）: 対数スケールで見ると、距離は概ね指数的に減少しています。これはLQR制御の特性（閉ループ系の固有値が負の実部を持つ）から期待される通りです。初期の急速な接近の後、目標付近では微細な調整が行われています
-
制御入力（c）: 初期には比較的大きな加速度（制限値に達する場合もある）が必要ですが、目標に近づくにつれて急速に減少します。$f_z$（動径方向）が最も大きいのは、重力勾配に逆らって高度を変更する必要があるためです
-
燃料消費（d）: 累積デルタVは初期に急増し、その後は緩やかに飽和していきます。500 mの距離からの接近に数 m/s 程度のデルタVで済むことがわかり、近傍運用の燃料効率の良さが確認できます
ランデブー全体のシナリオ可視化
最後に、V-bar上の複数保持点を経由するランデブーシナリオ全体を可視化します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
mu_earth = 3.986004418e14
R_earth = 6.371e6
h_ISS = 400e3
a_target = R_earth + h_ISS
n = np.sqrt(mu_earth / a_target**3)
T_orb = 2 * np.pi / n
def cwh_stm_xz(t, n):
nt = n * t
s = np.sin(nt)
c = np.cos(nt)
return np.array([
[1, 6*(nt - s), (4*s)/n - 3*t, 2*(1 - c)/n],
[0, 4 - 3*c, 2*(c - 1)/n, s/n ],
[0, 6*n*(c - 1), 4*c - 3, 2*s ],
[0, 3*n*s, -2*s, c ]
])
def two_impulse_transfer(r0, rf, t_transfer, n):
"""二体力マニューバのデルタV計算"""
Phi = cwh_stm_xz(t_transfer, n)
Phi_rr = Phi[0:2, 0:2]
Phi_rv = Phi[0:2, 2:4]
Phi_vr = Phi[2:4, 0:2]
Phi_vv = Phi[2:4, 2:4]
v0_plus = np.linalg.solve(Phi_rv, rf - Phi_rr @ r0)
vf_arrival = Phi_vr @ r0 + Phi_vv @ v0_plus
return v0_plus, -vf_arrival # dv1の速度, dv2
def cwh_equations_2d(t, state, n):
"""2D CWH方程式（x-z面内）"""
x, z, vx, vz = state
ax = 2 * n * vz
az = -2 * n * vx + 3 * n**2 * z
return [vx, vz, ax, az]
# --- ランデブーシナリオ ---
# 保持点: HP3(-2000, 0) -> HP2(-500, 0) -> HP1(-100, 0) -> Target(0, 0)
hold_points = [
np.array([-2000.0, 0.0]), # HP3 (出発点)
np.array([-500.0, 0.0]), # HP2
np.array([-100.0, 0.0]), # HP1
np.array([0.0, 0.0]), # Target
]
# 各レグの転送時間
transfer_times = [0.5 * T_orb, 0.4 * T_orb, 0.3 * T_orb]
fig, ax = plt.subplots(figsize=(14, 6))
total_dv = 0.0
all_x = []
all_z = []
colors = ['#2196F3', '#4CAF50', '#FF5722']
labels = ['HP3 → HP2', 'HP2 → HP1', 'HP1 → Target']
for i in range(len(transfer_times)):
r0 = hold_points[i]
rf = hold_points[i+1]
t_tf = transfer_times[i]
v0_plus, dv2 = two_impulse_transfer(r0, rf, t_tf, n)
dv1_mag = np.linalg.norm(v0_plus)
dv2_mag = np.linalg.norm(dv2)
total_dv += dv1_mag + dv2_mag
# 軌道の数値積分
state0 = [r0[0], r0[1], v0_plus[0], v0_plus[1]]
t_eval_leg = np.linspace(0, t_tf, 1000)
sol_leg = solve_ivp(
cwh_equations_2d, (0, t_tf), state0,
args=(n,), t_eval=t_eval_leg, rtol=1e-10, atol=1e-12
)
ax.plot(sol_leg.y[0], sol_leg.y[1], color=colors[i],
linewidth=2, label=f'{labels[i]} (Δv={dv1_mag+dv2_mag:.4f} m/s)')
all_x.extend(sol_leg.y[0])
all_z.extend(sol_leg.y[1])
# 保持点のプロット
hp_labels = ['HP3\n(-2000 m)', 'HP2\n(-500 m)', 'HP1\n(-100 m)', 'Target']
hp_colors = ['green', 'blue', 'orange', 'red']
hp_markers = ['o', 's', 'D', '*']
hp_sizes = [10, 10, 10, 15]
for i, (hp, label, color, marker, size) in enumerate(
zip(hold_points, hp_labels, hp_colors, hp_markers, hp_sizes)):
ax.plot(hp[0], hp[1], marker=marker, color=color,
markersize=size, zorder=5)
offset_y = 15 if i < 3 else 15
ax.annotate(label, (hp[0], hp[1]),
textcoords="offset points", xytext=(0, offset_y),
ha='center', fontsize=10, fontweight='bold', color=color)
# V-bar軸
ax.axhline(y=0, color='gray', linewidth=0.8, linestyle='-', alpha=0.5)
ax.axvline(x=0, color='gray', linewidth=0.8, linestyle='-', alpha=0.5)
ax.text(-1000, -5, 'V-bar', fontsize=11, color='gray', alpha=0.7, ha='center')
ax.set_xlabel('x (along-track) [m]', fontsize=12)
ax.set_ylabel('z (radial) [m]', fontsize=12)
ax.set_title(f'Multi-Hop V-bar Rendezvous Scenario (Total Δv = {total_dv:.4f} m/s)', fontsize=13)
ax.legend(fontsize=10, loc='upper left')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('rendezvous_scenario.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"\n合計デルタV: {total_dv:.4f} m/s")
このマルチホップ・ランデブーのシミュレーション結果から、以下の点が読み取れます。
-
段階的な接近: 2000 m → 500 m → 100 m → 0 m と段階的に距離を詰めていく様子が可視化されています。各保持点で一旦停止し、次のレグを計画するという実運用に近いシナリオです
-
各レグの軌道形状: 遠距離のレグほど $z$ 方向への逸れが大きくなっています。これはCWH方程式の特性で、長い距離を移動するほどコリオリ力の影響が顕著になります。近距離のレグ（HP1 → Target）ではほぼ直線に近い軌道が得られます
-
燃料効率: 2 kmの距離からの完全なランデブーに必要な合計デルタVは非常に小さく、近傍運用における推進剤の使用量が全ミッションのバジェットの中では限定的であることがわかります。ランデブーの燃料消費の大部分は、遠距離フェーズでの軌道遷移（ホーマン遷移など、CWH方程式の適用範囲外）にあります
安全楕円と最終接近の3D可視化
法線方向の運動も含めた3次元的な相対軌道を可視化します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
mu_earth = 3.986004418e14
R_earth = 6.371e6
h_ISS = 400e3
a_target = R_earth + h_ISS
n = np.sqrt(mu_earth / a_target**3)
T_orb = 2 * np.pi / n
def cwh_equations(t, state, n, u_func=None):
x, y, z, vx, vy, vz = state
if u_func is not None:
fx, fy, fz = u_func(t, state)
else:
fx, fy, fz = 0.0, 0.0, 0.0
ax = 2 * n * vz + fx
ay = -n**2 * y + fy
az = -2 * n * vx + 3 * n**2 * z + fz
return [vx, vy, vz, ax, ay, az]
# --- 安全楕円（3D, 法線方向の振動も含む）---
z_a = 200 # 動径方向の振幅 [m]
y_a = 50 # 法線方向の振幅 [m]
state0_3d_ellipse = [
0.0, # x0
y_a, # y0
z_a, # z0
2*n*z_a, # vx0 = 2*n*z0（安全楕円条件）
0.0, # vy0
0.0 # vz0
]
t_span_3d = (0, 2 * T_orb)
t_eval_3d = np.linspace(0, 2 * T_orb, 3000)
sol_3d = solve_ivp(
cwh_equations, t_span_3d, state0_3d_ellipse,
args=(n,), t_eval=t_eval_3d, rtol=1e-10, atol=1e-12
)
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot(sol_3d.y[0], sol_3d.y[1], sol_3d.y[2],
'b-', linewidth=1, label='Safety ellipse (3D)')
ax.plot([sol_3d.y[0][0]], [sol_3d.y[1][0]], [sol_3d.y[2][0]],
'go', markersize=8, label='Start')
ax.scatter([0], [0], [0], c='red', marker='*', s=200, label='Target', zorder=5)
# x-z面への射影
ax.plot(sol_3d.y[0], np.full_like(sol_3d.y[0], -80),
sol_3d.y[2], 'b--', alpha=0.3, linewidth=0.5)
# x-y面への射影
ax.plot(sol_3d.y[0], sol_3d.y[1],
np.full_like(sol_3d.y[0], -250), 'b--', alpha=0.3, linewidth=0.5)
ax.set_xlabel('x (along-track) [m]', fontsize=11)
ax.set_ylabel('y (cross-track) [m]', fontsize=11)
ax.set_zlabel('z (radial) [m]', fontsize=11)
ax.set_title('3D Safety Ellipse in LVLH Frame', fontsize=13)
ax.legend(fontsize=10)
ax.view_init(elev=25, azim=-60)
plt.tight_layout()
plt.savefig('safety_ellipse_3d.png', dpi=150, bbox_inches='tight')
plt.show()
3Dプロットから、安全楕円が3次元的にはリサジュー図形のような軌道を描くことがわかります。$x$-$z$ 面内では2:1の楕円、$y$ 方向には独立した単振動が重畳されているため、チェイサーはターゲットの周囲を立体的に周回します。射影（破線）を見ると、$x$-$z$ 面内では閉じた楕円、$x$-$y$ 面内では正弦波状の軌跡になっていることが確認できます。
法線方向の振動周期は軌道面内の周期と同じ（ともに $2\pi/n$）ですが、一般には位相がずれるため、3次元的には楕円が「ねじれた」ような軌道になります。この軌道は依然として受動安全であり、$y$ 方向がゼロでない場合もターゲットとの衝突は起こりません。
まとめ
本記事では、自律ランデブー・近傍接近（ARPO）の誘導制御について、基礎理論からPythonシミュレーションまでを解説しました。
-
CWH方程式は、ターゲット円軌道近傍でのチェイサーの相対運動を記述する線形方程式であり、ランデブー誘導制御の基盤です。進行方向と動径方向がコリオリ力で結合し、法線方向は独立な単振動になるという構造を持ちます
-
安全楕円（$\dot{x}_0 = 2nz_0$）は受動安全な保持軌道であり、スラスタ故障時にもターゲットとの衝突を回避できます。2:1のアスペクト比（along-track : radial）はCWH方程式の固有の特徴です
-
ランデブーは段階的に行われます: 遠距離→中距離→近傍→最終接近の各フェーズで、使用するセンサ・誘導法則・精度要求が異なります。V-barアプローチとR-barアプローチにはそれぞれ長所と短所があり、ミッション要求に応じて選択されます
-
誘導法則は、二体力マニューバ（開ループ）からLQRフィードバック（閉ループ）、グライドスロープ誘導まで多様な選択肢があります。実際のミッションでは距離に応じてこれらを組み合わせます
-
安全性設計が最優先です。受動安全軌道、CAM（衝突回避マニューバ）、保持点でのGo/No-Go判定といった多重の安全対策が、信頼性の高いランデブーを支えています
ランデブー・近傍接近技術は、宇宙ステーション補給にとどまらず、軌道上サービス（OOS）、デブリ除去、宇宙組立など、今後の宇宙開発の広い分野で必要とされる基盤技術です。
次のステップとして、以下の記事も参考にしてください。
- 宇宙機の捕獲技術 — カナダアーム・JEMRMSの制御 — ランデブー後のロボットアームによる捕獲制御
- 軌道上サービス（OOS）の全体像 — OOSミッション全体の中でのランデブーの位置づけ
- 軌道遷移とデルタV — 遠距離フェーズで使われるホーマン遷移の詳細