これはABEJAアドベントカレンダー2024年の記事である。普段プロダクトマネージャーを務める筆者が、2024年10月のSpaceXの衝撃的な垂直着陸成功を見てこの記事を書くに至った。

筆者の過去記事一覧

本日はSpaceXのロケットについて語るために、自作のシミュレーターを開発した。ロケット開発やその制御技術は言葉だけで表現するには限界がある。実際にブラウザで動くので、ぜひロケットの奥深さや面白さを触ってためしていただきたい。上空から地上の一点を目指すのがいかに難しいか、操作してみるとよりリアルにわかるはずだ。

シミュレーターはこちら：https://rocket-simulator-demo.web.app

それでは、なぜシミュレーターを作るに至ったのかお話する。

0. まえがき
1. はじめに：SpaceXが世界にもたらした技術革新
2. 理論：ロケット技術の基本解説
3. 設計：シミュレーターの設計
4. 実装：3Dシミュレーターの実装
5. 実験：離着陸のシミュレーション
シミュレーション結果
- 成功例
- 失敗例
6. まとめ

0. まえがき

SpaceXはスターシップをはじめとした火星植民地化のためのロケット開発産業や、スターリンクと呼ばれる衛星産業など幅広く展開している。そのビジネス面のすごさは今回詳しくは言及しないが、私が紹介したいのはSpaceXの根幹をなすロケットの制御技術についてである。

SpaceXのロケットは、なんと打ち上げ後に地球上の任意の位置に垂直着陸し帰還することができるのだ。軟着陸させるためのパラシュートや、滑走路すら必要としない。

2015年12月21日、はじめて着陸映像を見た時には思わず笑ってしまったのを覚えている。

成層圏を超える6,000m以上の上空から誤差1m以下で着陸をさせるなど、映像を見ていなければ信じていないかもしれない。私は大学で小型航空機やドローンなどを設計する研究をしていたことがあり、飛行制御の難しさは3年間毎日痛感していた。ただまっすぐ飛ばすのでさえ、何度もの反復的な設計やアルゴリズムの改善、そしてハードウェアの開発を必要とするのだ。

2024年10月13日には、なんと機械式アームでの空中キャッチまでやってのけている。この機体、200トン以上もあることが信じられない。

このSuper Heavyと名付けられた超大型ロケットは、燃料を含めた全重量が250トン以上あり高さは60mを超える。ビルの高さでいうと20階相当になる。Super Heavyやより小型のFalcon9はいずれも打ち上げ後にそのブースターを地上のピンポイントな位置に戻し、海上の母艦や地上の回収タワーで回収されるという、他のロケットにはありえない機能を備えている。

ただ、このような説明を聞いただけでは技術のすごさがなかなか通じない。圧倒的技術力を前に説明する言葉がないのだ。

かくして私はその仕組みをわかりやすく解説し伝えるべく、シミュレーターを開発した。この記事では、SpaceXのロケット技術の凄さを自作シミュレーターを通じて解説していく。

1. はじめに：SpaceXが世界にもたらした技術革新

SpaceXの最大のブレイクスルーは、ロケットを再使用可能な運搬システムへと転換させた点にある。従来のロケットはほぼ全てが「使い捨て」であり、打ち上げごとに膨大な構造物とエンジンを廃棄するため、コストが指数的に増大していた。これに対し、SpaceXは Falcon 9 や Falcon Heavy、そして Starship といった機種群で、打ち上げ後のブースターを回収し、検査・整備した上で再利用する技術を確立した。このアプローチは、単なるコストダウンに留まらず、宇宙アクセスのパラダイムシフトを引き起こしている。

SpaceXは、2023年には過去2年間でロケット打ち上げ数を2倍に伸ばし、2024年の打ち上げ回数はすでに105回を超えている。つまりほぼ3日に1回近く打ち上げをしている。 *1

再使用可能なロケットによる新たな宇宙開発時代

再使用性を成立させるには、従来考えられていた以上の高度な設計・運用が求められる。まず、構造設計においては、単発ミッション用の軽量・高比推力志向から、繰り返し使用可能な耐久性重視の設計指針へとシフトする必要がある。機体構造材は熱・振動・荷重サイクルに繰り返し耐えることが要求され、エンジンは何度も再点火・燃焼できる設計が不可欠だ。

次に、回収技術は再使用性の根幹を成す。従来は帰還のためのパラシュートやシャトル型の複雑な着陸システムが必要とされると考えられていたが、SpaceXは垂直離着陸（VTVL）方式を採用することで、着地精度と回収効率を格段に高めた。エンジンによる逆推力、グリッドフィン（後述）による大気圏内での精密な姿勢・飛行経路制御、コールドガススラスター（後述）による微調整など、数理最適化や制御理論、空力設計を総合的に組み合わせることで、狭小な着地ゾーンへの高精度再突入とソフトランディングを実現している。これらの技術基盤の上で、SpaceXは海上ドローンシップや地上着陸パッドへの安定的なブースター回収を繰り返し成功させている。

SpaceXの独自技術

StarshipとSuper Heavyは、SpaceXが目指す「フル再使用可能な宇宙往還システム」を構成する上段・下段の2要素である。今回シミュレーションで利用させていただいた3Dモデル*2を使いながら説明する。

また、その構成要素であるグリッドフィン、コールドガススラスター、メインエンジンのジンバル機構はロケットの正確な誘導・制御を実現するための重要なアクチュエータ類である。それぞれ異なる物理原理・速度域・制御範囲で機体を安定化・誘導する役割を果たす。

Super Heavy（1段目ブースター）

Super HeavyはStarship上段を宇宙空間近くまで加速するための超大型ブースターである。33基ものRaptorエンジンを束ね、世界有数の推力・ペイロード能力を誇る。メタン・液体酸素推進剤を使用し、フルフローステージ燃焼サイクルにより高い効率とエンジン寿命を実現する。

Starship（2段目ペイロード）

Starshipは、Super Heavyによって上空へ運ばれた後、単独で軌道へと到達する再使用可能な上段ステージ兼宇宙船である。Starship自体もRaptorエンジンを搭載し、地球低軌道への大量貨物や100人規模の乗組員輸送、さらに月・火星までの長距離ミッションを想定する。

Rapter （メインエンジン）

従来の二段燃焼サイクルエンジンでは、燃料リッチ（混合比における比率の高い状態）または酸化剤リッチな予燃焼室（Preburner）のいずれか一方の組み合わせを使うケースが一般的だった。一方、Raptorエンジンは燃料リッチ、酸化剤リッチの両予燃焼室を同時に使用し、ターボポンプをそれぞれ駆動した後、全てのガスをメイン燃焼室に導く。これにより、理論上の燃焼効率を最大限まで引き上げることが可能になり、非常に高い比推力（Isp）および推力重量比を実現している。

グリッドフィン

グリッドフィンは、ロケットの先端部やインターステージ付近に格子状のフィンを配置し、それを舵角変更可能な翼面として用いるシステムである。主に大気圏内で作用し、亜音速から超音速領域にかけて有効な姿勢制御を行う。

グリッドフィンは、従来の平板翼・カナード翼に比べて広いマッハ数領域（超音速から亜音速領域への減速過程）で安定して制御力を発生できる。
フィンの格子構造により、渦構造が複雑かつ安定的に形成され、多方向性の舵力が得られ、微細な姿勢修正から大きな軌道偏差補正まで対処可能となる。

コールドガススラスター

コールドガススラスターは、圧縮ガス（一般には窒素など）をノズルから噴出することで小規模な推力を発生する制御システムである。燃焼を伴わず、相対的に低推力・単純構造のRCS（Reaction Control System）として用いられる。

コールドガススラスターは極めて小さな推力制御を実現できるため、ロケット姿勢や横方向偏差を細かく調整するために使用される。特に高高度や真空中では空力的な翼であるグリッドフィンが効かないため、微小なスラスター推力が重要となる。
燃焼制御が不要なため、点火・停止が容易で応答性に優れ、信頼性が高い。緊急時の姿勢リカバリや最終接地直前の微調整にも適している。

メインエンジンのジンバル機構

メインエンジンは、単に固定方向へ推力を発生させるだけでなく、ノズル自体をピッチ・ヨー方向へ可動（ジンバル）できるようにすることで、推力ベクトル制御（TVC: Thrust Vector Control）を行う。

メインエンジンジンバルは、高い推力を持った主推進器そのものの方向を変えるため、全体質量の大きなロケットの強力な姿勢制御手段となる。
主に上昇段階での軌道修正や、再使用時には降下中・着地直前に推力方向を変えて機体のトルクを発生させ、ロール・ピッチ・ヨーを補正する。

宇宙ビジネスの変化

上記のような再回収のための技術が実装されたことでコスト構造が変化し、必然的に宇宙開発のビジネスモデルが変容している。これまで一回の打ち上げで数千万〜数億ドル規模のコストが発生していたが、再使用可能なブースターは一回あたりのロケットハードウェア廃棄コストを大幅に削減する。結果としてスターリンクに代表される衛星インターネットサービス、科学ミッションの増加、民間有人宇宙飛行など、以前は財政的・技術的なボトルネックに阻まれていた構想が現実性を増している。

特にStarshipが目指す大規模再使用は、さらなる打ち上げコスト低減と打ち上げ頻度向上に道を開き、火星移住といった長期的な目標においても、持続的な貨物輸送および人員輸送インフラを提供し得る。総じて、SpaceXの再使用可能ロケットは、材料・構造設計、推進系、飛行力学、制御理論、経済性分析という複数領域に跨る総合的な最適化の成果である。

2. 理論：ロケット技術の基本解説

ロケットを開発するためには、まずロケットそのものの動きをモデル化する必要がある。基礎となるのは運動方程式で、ロケットが発生させる推力、大気中の空気抵抗、地球の重力などの要素を正確にモデル化しなければならない。また、ロケットの形状や質量などの機体諸元を考慮した設計が求めらる。

ロケット推進の基本方程式

ロケット推進を解析する上で最も基本的な理論的枠組みは、ツィオルコフスキーのロケット方程式である。この方程式は、燃焼によってロケット質量が連続的に減少するという時間変化を内包しつつ、噴射速度（排気流速度）と質量変化率の関係から最終的な速度増分Δvを求めるものだ。式は以下のように表される。

ここで、 $m_0$ は打ち上げ前の総質量（ペイロード、構造材、推進剤を含む）、 $m_f$ は燃焼後の最終質量（ペイロードと機体構造のみ、推進剤消費後）、 $v_e$ は有効排気速度（エンジンが噴射する推進剤流の相対速度）である。ロケットのエンジンは後方に質量（推進剤）を高速で放出することで、その反作用として前方推力を得る。結果として、ロケットは連続的な質量損失を伴いつつ加速し、最終的な速度を獲得する。

ロケットの運動

ロケットの運動は、質量変化を伴う非定常な力学系として記述される。これを表す基本方程式はニュートン・オイラー方程式に質量変化項や大気力学的影響を考慮した拡張形とみなせる。ロケットは打ち上げから減速・着地に至るまで、重力、空気抵抗（揚力・抗力）、推力、および質量減少による運動特性の変化に常時さらされる。さらに、上昇高度や飛行速度によって大気密度や音速が変化し、空力特性が高度非線形的に変動するため、運動方程式は単純な線形近似では対処しきれない複雑性を持つ。基本的な考え方として、ロケットを剛体または質点系とみなし、慣性座標系（例えば地心慣性座標系）においてロケットに作用する全ての力を合計し、運動方程式を記述する。一般的な6自由度（6-DoF）の運動方程式を考えると、並進運動と回転運動が組み合わさった非線形微分方程式系となる。並進運動は以下のように表現できる。

ここで、 $m(t) )$ は時間とともに変化するロケット総質量、 $\mathbf{v}\$ はロケットの速度ベクトルである。右辺には各種力が列挙される。

$\mathbf{T}$ : 推力（Thrust）。エンジン噴射による加速力であり、ロケット内部座標系で規定された推力方向を慣性座標系へ変換して適用する。推力は噴射ガス速度と燃焼流量によって決まり、エンジン出力調整（スロットリング）やベクトル制御によって方向も時間変化する。
$\mathbf{D}$ : 抗力（Drag）。大気中を高速で移動する際に発生する流体抵抗であり、基本的にはロケット進行方向に逆らう向きに作用する。抗力係数はマッハ数、迎角、レイノルズ数、機体形状によって変化し、モデル化にはCFD（計算流体力学）的解析や風洞実験の結果が用いられる。
$\mathbf{W}$ : 重力（Weight）。地球重力場による引力であり、質量 $m(t)$ に対して一定の方向（ほぼ鉛直下向き）に作用する。高度や位置により微妙に重力加速度が変化するため、近似モデル（標準重力モデル）を用いることが多い。
$\mathbf{L}$ : 揚力（Lift）。ロケット形状が非対称な場合や迎角を与えた制御で発生する横方向の空力力。これによりロケットは若干の横移動や軌道修正が可能になる。着地や再突入時にはグリッドフィンやエアロサーフェスが揚力・モーメントを発生し、機体姿勢と軌道を制御する。
$\mathbf{F}_{\text{other}}$ : コールドガススラスターやロール制御用スラスター、あるいは風の突発的な乱流による外乱力など、その他の外力が考えられる。

回転運動（姿勢制御）も重要で、機体の偏心推力や空気力中心の変動によりロール・ピッチ・ヨー方向にモーメントが発生する。このモーメントは、機体慣性モーメントテンソルと時間変化する質量分布によって複雑化し、非線形な姿勢運動を引き起こす。誘導・制御システムはこれを予測・補償するため、推力ベクトル制御、グリッドフィン舵角調整、スラスター噴射などを組み合わせる。

運動方程式

剛体近似したロケットに対する六自由度(6-DoF)運動方程式を示す。ここでは、ロケット機体固定座標系（Body-Fixed Frame）を用い、(x, y, z)軸を機体中心にとり、x軸前方、y軸右方向、z軸下向きとする。ロケットの並進運動と回転運動はそれぞれ独立な方程式系で表されるが、相互作用する。

並進運動方程式

ロケット質量を m、機体座標系におけるロケットの速度成分を (u, v, w)、角速度成分を (p, q, r) とし、機体座標系で表した外力を (X, Y, Z) とする。これらは推力、重力（機体軸投影）、空気力（揚力・抗力・側力）、および制御スラスターによる力の和である。ニュートンの運動方程式を機体軸系で記述すると、以下のような非線形連立微分方程式となる。ここで、(r v - q w), (p w - r u), (q u - p v) の項は、機体回転による速度成分の結合項であり、回転運動と並進運動が組み合わさっていることを示す。この基本となる式に、離着陸・飛行時の操作入力を加えた場合が次のようになる。メインエンジンの推力に加え、グリッドフィンによる揚力変化、コールドガススラスターによる推力変化の項目を加えている。

回転運動方程式

回転運動はオイラーの剛体方程式によって記述される。機体の慣性モーメントテンソルを次のように示す。

機体座標系でのモーメント（外力モーメント）を (L, M, N) とする。ここで L は x軸まわり、M は y軸まわり、N は z軸まわりのモーメントである。角速度ベクトルをω=(p,q,r) とすれば、回転運動方程式は以下のベクトル形式で表される。

成分表示すると非線形項を含む複雑な形となるが、慣性モーメントの対称性や、主軸系を用いることで簡略化できる。たとえば主軸系で交差項 [tex: I{xy}, I{xz}, I_{yz} ] が0の場合、以下のようになる。

モーメント方程式には、グリッドフィン由来のモーメント [tex: (L{GF}, M{GF}, N_{GF})] を追加する。同様に、エンジンモーメント $(L_T, M_T, N_T)$ と他の外乱モーメント $(L, M, N)$ に加算する。

このように、グリッドフィンとコールドガススラスターを加えた場合でも、基本となる6自由度運動方程式の枠組みは維持される。新たな制御デバイスごとに発生する力・モーメントを、既存の合力・合モーメント項へ足し込むことで、より複雑かつ現実的な制御機構をモデル化できる。

制御理論の概要

GNCは、宇宙機やロケットなどの飛行体を所望の目標へ誘導・制御する総合的な概念であり、3つの階層的要素を含む。

Navigation（航法）：自機の位置・速度・姿勢といった状態量を推定するフェーズ。IMU（慣性計測装置）、GPS、スターセンサー、地上局トラッキング情報など、複数の計測データを統合し、フィルタ（カルマンフィルタなど）を通じて精度の高い状態推定を行う。この段階で、現在の座標x,y,zや姿勢角𝜓,θ,ψなどの値が高精度に求まる。

Guidance（誘導）：状態推定された現状を踏まえ、目標軌道や着地地点への最適な飛行パスを計算・更新するフェーズ。誘導ロジックは、時々刻々変化する外乱や残燃料量、時間制約を考慮し、次の瞬間どのような目標速度・高度・角度を追従すべきかを決定する。たとえば「残り高度1000mで、横ずれを打ち消すためヨー角を数度修正せよ」「対気速度を指定の範囲に収めるためスロットルを下げよ」といった動的な基準値（目標値）を提供する。この層は多くの場合、モデル予測制御や最適制御、あるいは事前に計画された基準トラジェクトリとリアルタイム補正を組み合わせる。

Control（制御）：誘導段階で設定された目標状態（姿勢角度、速度、位置オフセットなど）に対して、実際の操作量（メインスラスター出力、グリッドフィン舵角、コールドガススラスター噴射コマンド）を決定するフェーズが制御である。ここでPIDやLQR、MPCなどのコントローラが利用され、下位レベルのアクチュエータ指令が生成される。

今回は、シミュレーターでの検証として比較的実装が容易なPID（Proportional-Integral-Derivative）制御を用いることを想定する。PID制御は、目標値と実測値の偏差に基づき、比例成分で即時応答、積分成分で定常偏差の削減、微分成分で過渡応答の抑制といった特性を組み合わせることで、安定した制御動作を実現する単純かつ強力な手法である。

本シナリオでは、到達地点の座標（水平位置）、高度(altitude)、機体の姿勢（特にpitchとyaw）、そして対気速度(airspeed)を対象とする。具体的には下記のような制御ループが考えられる。

位置・高度制御：目標着地地点の座標(x,y,z)_target および最終到達高度に対して、現在の機体位置(x,y,z) との偏差を計算し、その差を低減するような外力が必要となる。水平軸方向では主にグリッドフィンおよびコールドガススラスターを用いた横方向制御、垂直方向ではメインスラスター出力調整によって高度・降下率をコントロールする。

姿勢（Pitch, Yaw）制御：姿勢制御では、機体のピッチ角度（θ）とヨー角度（ψ）を目標値へ追従させる。たとえば、滑らかな降下を行うために一定のピッチ角を維持する、あるいは風などによる横方向ズレを補正するためヨー角を適正な値に保持する、といった戦略がとられる。これらは主にグリッドフィン、コールドガススラスターで生じる横方向モーメントや、メインエンジン推力ベクトル偏向が利用できる。

対気速度(Airspeed)制御：対気速度は、最終的な着地精度と運動エネルギー管理において重要なパラメータである。目標速度より速すぎれば、着地時に大きな逆噴射（ブレーキバースト）が必要になるし、遅すぎれば制御舵面（グリッドフィン）の有効性が下がる可能性がある。メインスラスターを微調整して推力を増減させるとともに、適切なピッチ角・ヨー角設定で空力抗力を変化させ、目標速度帯を維持する。

3. 設計：シミュレーターの設計

今回のシミュレーターを作るために条件を設定した。

だれでも試せること
手触り感があること

1に関して、インストールの必要なシミュレーターを用意しても使ってもらえないケースが多いと考えられる。ブラウザでも動くことが要件となる。 2に関して、実際に人がインタラクティブに操作してもらえることを目指す。リアルタイムでの計算が必要かつ、ビジュアル的にわかりやすいものが必要となる。そこで今回はブラウザでリッチなアニメーションを表示可能なThree.jsを用い、TypeScriptでバックエンドを介さずリアルタイム演算ができるシミュレーションを構築することにした。

運動方程式のモデル化

前述の運動方程式をTypeScriptで実装する。定式化したものに加え、フィンの空力学的特性などを仮定する。スクリプトの一部を示す。

// Grid fins produce force roughly proportional to deflection
// Pitch fins produce X-direction force
const Fx_fin = q * this.gridFinArea * this.C_N_gridFin * gridFinPitch;
// Yaw fins produce Z-direction force
const Fz_fin = q * this.gridFinArea * this.C_N_gridFin * gridFinYaw;
const Fy_fin = 0; // simplified: no Y-direction force from fins

// Compute moments
// Thrust
const thrusterForceVec = new THREE.Vector3(Fx_thrust, Fy_thrust, Fz_thrust);
const M_thrust = new THREE.Vector3().crossVectors(this.thrusterPosition, thrusterForceVec);

// Pitch fin force
const finForcePitchVec = new THREE.Vector3(Fx_fin, 0, 0);
const M_finPitch = new THREE.Vector3().crossVectors(this.finArmPitch, finForcePitchVec);

// Yaw fin force
const finForceYawVec = new THREE.Vector3(0, 0, Fz_fin);
const M_finYaw = new THREE.Vector3().crossVectors(this.finArmYaw, finForceYawVec);

const M_fin = new THREE.Vector3(
    M_finPitch.x + M_finYaw.x,
    M_finPitch.y + M_finYaw.y,
    M_finPitch.z + M_finYaw.z
);

// Total forces
const Fx_total = Fx_normal + Fx_thrust + Fx_fin;
const Fy_total = Fy_axial + Fy_thrust + Fy_fin;
const Fz_total = Fz_normal + Fz_thrust + Fz_fin;

// Total moments
const Mx_total = M_thrust.x + M_fin.x;
const My_total = M_thrust.y + M_fin.y;
const Mz_total = M_thrust.z + M_fin.z;

外部環境（重力、風、大気圧）のモデル化

今回は簡易的に重力、風外乱のモデル化を行っている。また、座標系について本来は地球が球体であることも加味する必要があるが、三次元の直交座標系として扱う。今後シミュレーションを拡張する際にはこの辺りも修正していきたい。大気圧についても今回は一定とし、高度に対応したモデル化は行なっておらず制御が比較的容易な条件となる。

重力モデル

export class Environment {
    private gravity = 9.81; // m/s^2

    getGravityInEarthFrame(mass: number): [number, number, number] {
        // Y軸鉛直方向
        return [0, -mass * this.gravity, 0];
    }
}

風モデル

export class Wind {
    private constantWind: [number, number, number]; // [x, y, z] in Earth Frame
    private turbulenceIntensity: number;

    constructor(constantWind: [number, number, number], turbulenceIntensity: number = 0) {
        this.constantWind = constantWind;
        this.turbulenceIntensity = turbulenceIntensity;
    }

    getWindVector(): [number, number, number] {
        // 定常風
        const windX = this.constantWind[0];
        const windY = this.constantWind[1];
        const windZ = this.constantWind[2];

        // 乱流（ノイズを追加）
        const turbulenceX = this.turbulenceIntensity * (Math.random() - 0.5);
        const turbulenceY = this.turbulenceIntensity * (Math.random() - 0.5);
        const turbulenceZ = this.turbulenceIntensity * (Math.random() - 0.5);

        return [
            windX + turbulenceX,
            windY + turbulenceY,
            windZ + turbulenceZ,
        ];
    }
}

入力のモデル化

メインエンジン、グリッドフィン、コールドガススラスターによる入力は、人間のものと自動操作の両方を想定する。操作は切り替えられるようにし、シミュレーション画面からスライダーで入力できる。

export interface State {
    position: [number, number, number];       // Earth Frame (NED)
    orientation: [number, number, number, number]; // Quaternion representing Body->Earth rotation
    velocityBody: [number, number, number];   // Body Frame
    angularVelocityBody: [number, number, number]; // Body Frame
}

export interface Inputs {
    thrustPitch: number;
    thrustYaw: number;
    gridFinPitch: number;
    gridFinYaw: number;
    throttle: number;
    windEarth: [number, number, number]; // Wind in Earth Frame
}

export interface ForcesMoments {
    Fx: number; // body X-axis
    Fy: number; // body Y-axis
    Fz: number; // body Z-axis
    Mx: number; // roll moment about body X
    My: number; // pitch moment about body Y
    Mz: number; // yaw moment about body Z
}

姿勢計算と積分

グローバル座標系と機体座標系の変換や、クオータニオンによる姿勢角の計算は下記のように行なっている。クオータニオンの正規化が行えていない場合にシミュレーション結果が発散する・わずかにズレるなどが、クオータニオンあるあるだ。積分について、シンプルにオイラー法で行っているが今後改善も検討している。

export class Aircraft {
    private state: State;
    private mass: number;
    private inertia: [number, number, number];
    private aero: Aerodynamics;
    private env: Environment;
    private wind: Wind;

    constructor(
        initialState: State,
        mass: number,
        inertia: [number, number, number],
        aero: Aerodynamics,
        env: Environment,
        wind: Wind,
    ) {
        this.state = initialState;
        this.mass = mass;
        this.inertia = inertia;
        this.aero = aero;
        this.env = env;
        this.wind = wind;
    }

    computeDerivative(inputs: Inputs): State {
        // 空力
        const fm = this.aero.getForcesAndMoments(this.state, inputs);

        // 重力 (Earth Frame)
        const gravityEarth = this.env.getGravityInEarthFrame(this.mass);

        // 風(Earth Frame)
        const windEarth = this.wind.getWindVector();

        const windBody = CoordinateTransform.earthToBody(this.state.orientation, windEarth);
        // 重力をBody Frameへ変換
        const gravityBody = CoordinateTransform.earthToBody(this.state.orientation, gravityEarth);

        // Body Frameでの合力
        const Fx_total = fm.Fx + gravityBody[0] + windBody[0];
        const Fy_total = fm.Fy + gravityBody[1] + windBody[1];
        const Fz_total = fm.Fz + gravityBody[2] + windBody[2];

        // 加速度(Body Frame)
        const ax = Fx_total / this.mass;
        const ay = Fy_total / this.mass;
        const az = Fz_total / this.mass;

        // 角速度
        const p = this.state.angularVelocityBody[0];
        const q = this.state.angularVelocityBody[1];
        const r = this.state.angularVelocityBody[2];

        const Ix = this.inertia[0];
        const Iy = this.inertia[1];
        const Iz = this.inertia[2];

        // 角加速度(Body Frame)
        const pdot = (fm.Mx - (q * r * (Iz - Iy))) / Ix;
        const qdot = (fm.My - (r * p * (Ix - Iz))) / Iy;
        const rdot = (fm.Mz - (p * q * (Iy - Ix))) / Iz;

        // 位置微分 (Earth Frame)
        // velocityBodyをEarthに変換
        const velocityEarth = CoordinateTransform.bodyToEarth(this.state.orientation, this.state.velocityBody);

        // Orientation微分 (四元数)
        const q0 = this.state.orientation[0];
        const q1 = this.state.orientation[1];
        const q2 = this.state.orientation[2];
        const q3 = this.state.orientation[3];

        const qdot0 = 0.5 * (-q1 * p - q2 * q - q3 * r);
        const qdot1 = 0.5 * (q0 * p + q2 * r - q3 * q);
        const qdot2 = 0.5 * (q0 * q - q1 * r + q3 * p);
        const qdot3 = 0.5 * (q0 * r + q1 * q - q2 * p);

        return {
            position: [velocityEarth[0], velocityEarth[1], velocityEarth[2]],
            orientation: [qdot0, qdot1, qdot2, qdot3],
            velocityBody: [ax, ay, az],
            angularVelocityBody: [pdot, qdot, rdot]
        };
    }

    update(dt: number, inputs: Inputs): void {
        const derivative = this.computeDerivative(inputs);

        // 状態更新
        this.state = {
            position: [
                this.state.position[0] + derivative.position[0] * dt,
                this.state.position[1] + derivative.position[1] * dt,
                this.state.position[2] + derivative.position[2] * dt,
            ],
            orientation: [
                this.state.orientation[0] + derivative.orientation[0] * dt,
                this.state.orientation[1] + derivative.orientation[1] * dt,
                this.state.orientation[2] + derivative.orientation[2] * dt,
                this.state.orientation[3] + derivative.orientation[3] * dt,
            ],
            velocityBody: [
                this.state.velocityBody[0] + derivative.velocityBody[0] * dt,
                this.state.velocityBody[1] + derivative.velocityBody[1] * dt,
                this.state.velocityBody[2] + derivative.velocityBody[2] * dt,
            ],
            angularVelocityBody: [
                this.state.angularVelocityBody[0] + derivative.angularVelocityBody[0] * dt,
                this.state.angularVelocityBody[1] + derivative.angularVelocityBody[1] * dt,
                this.state.angularVelocityBody[2] + derivative.angularVelocityBody[2] * dt,
            ]
        };

        // 高度が0m未満の場合に修正 Y軸鉛直方向
        if (this.state.position[1] < 0) {
            this.state.position[1] = 0;
            if (this.state.velocityBody[1] < 0) {
                this.state.velocityBody[1] = 0;
            }
        }
        // 四元数正規化
        const norm = Math.sqrt(
            this.state.orientation[0] ** 2 +
            this.state.orientation[1] ** 2 +
            this.state.orientation[2] ** 2 +
            this.state.orientation[3] ** 2
        );
        this.state.orientation = this.state.orientation.map(o => o / norm) as [number, number, number, number];
    }

    getState(): State {
        return this.state;
    }
}

4. 実装：3Dシミュレーターの実装

いよいよシミュレーターの実装である。基本的にはThree.jsのアニメーションサイクルを利用し、1ステップごとに前ステップ情報の更新、制御入力、計算、描画更新を繰り返す。

計算は1ステップがギリギリ50msで完了することが確かめられたため、アニメーションは20FPS（20Hz）で更新している。下記は、Three.jsによるアニメーション表示の例である。シーンと呼ばれる背景設定、カメラの設定、物体の輪郭（Geomtry）と表面情報（Material）の扱いなど、慣れると非常に面白いライブラリである。

// Import necessary modules from Three.js
import * as THREE from 'three';

function init() {
    // Create a scene
    const scene = new THREE.Scene();

    // Set up a camera with a 75-degree field of view, near and far planes
    const camera = new THREE.PerspectiveCamera(75, window.innerWidth / window.innerHeight, 0.1, 1000);

    // Create a WebGL renderer and set its size
    const renderer = new THREE.WebGLRenderer();
    renderer.setSize(window.innerWidth, window.innerHeight);
    
    // Add the renderer's canvas element to the document
    document.body.appendChild(renderer.domElement);

    // Create a geometry, e.g., a box
    const geometry = new THREE.BoxGeometry(1, 1, 1);
    
    // Create a basic material and apply a color to it
    const material = new THREE.MeshBasicMaterial({ color: 0x00ff00 });
    
    // Create a mesh from geometry and material
    const cube = new THREE.Mesh(geometry, material);
    
    // Add the mesh to the scene
    scene.add(cube);

    // Position the camera so that we can see the cube
    camera.position.z = 5;

    // Create an animation loop
    function animate() {
        requestAnimationFrame(animate);

        // Apply some rotation to the cube for visual effect
        cube.rotation.x += 0.01;
        cube.rotation.y += 0.01;

        // Render the scene from the perspective of the camera
        renderer.render(scene, camera);
    }

    // Start the animation loop
    animate();
}

// Run the init function when the window loads
window.addEventListener('load', init);

5. 実験：離着陸のシミュレーション

シミュレーションはまず簡易に以下の条件で行った。試行時間が長くなりすぎるため、近距離でのテストを行なっている。

初期条件

const initialState: State = {
    position: [0, 0, 0], // 発射地点
    orientation: [1, 0, 0, 0], // 姿勢角（打ち上げ状態）は正常状態
    velocityBody: [0, 0, 0], // 速度0
    angularVelocityBody: [0, 0, 0] //角速度0
};

目標

const destination: Destination = {
    position: [300, 1, 300] , // 目標地点 x:300, z:300 離れた地点
    velocityBody: [0, 0, 0] // 速度0
};

環境条件

風外乱を印加しているが一定の値としている。

const windModel = new Wind([50, 0, 0], 2); // 5 m/s の東風 + 乱流強度2

達成条件

// ミッション達成判定用の閾値
const missionHorizontalThreshold = 100; // 半径100m以内の到達
const missionAltitudeThreshold = 4; // 高度4m以下への到達

シミュレーション結果

シミュレーターはこちら：https://rocket-simulator-demo.web.app

いよいよ本題となる、ロケット操作である。シミュレーターは左側のパネルから各種の操作入力が行える。エンジン推力を上げるとロケットが発射され、グリッドフィン、コールドガススラスターの出力によってロケットの進行方向を制御することができる。

目標地点は右側のミニマップに表示される。

約300m離れた位置に目標地点を設置したので、ぜひ目指してみて欲しい。上空から地上の一点を目指すのは、文字通り針の穴を通すように難しいことがわかる。

画面下部にあるグラフは、左から位置情報、姿勢角情報、迎角・対気速度の時間履歴になっている。x:300m, y:300m地点を目指す場合には左の位置情報グラフが300, 300になるように調整していく必要がある。 AutoPilotボタンを押した場合、途中からでもPID制御による自動介入が行われる。

成功例

PID制御系の数値を何パターンも試し、試行錯誤をするなかでAutoPilotボタンを押しさえすればターゲット付近に到達できるようになった。誤差は風外乱が5m/s（旗がはためく程度）である条件で試行し±30m程度となる。また、速度0m/sでの着陸を目標としているが完全に0にならず着陸する場合もあるため、高度コントロールのPIDゲインの調整余地がありそうだ。

今回はシンプルなPIDを用いたが、SpaceXは誘導アルゴリズムと最適化手法を用いて、降下中に常に理想的な逆噴射タイミング、エンジンスロットル開度、フィン舵角などを求め、誤差を瞬時に修正していると言われている。実際の着陸では突風や入力と出力の誤差、位置情報の実測値の取得など、シミュレーションで扱えていない困難がたくさんある。目標地点の取得自体も精度の観点からGPSだけではなく、レーダー高度計、LiDAR、カメラの情報をセンサーフュージョンし、ミリ秒単位でアクチュエーターを操作しなければならない。

失敗例

そもそも今回のシミュレーションで難しいのは、ロケットは不安定系であるということだ。

例えば自動車は安定な系で、平地で放置していても勝手に倒れたりすることはない。逆に自転車は不安定な系であるため平地でも手を離したら（制御入力なしでは）倒れてしまう。同様に飛行機の場合では安定な系に設計すれば制御入力がなくとも滑空は安定して行え、ロケットの場合は制御入力がないと自由に回転し落下するしかないのだ。

それゆえに今回開発したシミュレーションだけでのデバッグが難しく、位置・姿勢計算モジュール、モデルが悪いのか、制御器が悪いのかの切り分けが難しかった。いきなりシミュレーターではなく数式から求める方法もあるのだが今回は省略している。

安定な系で位置・姿勢計算モジュールをテストした例もある。飛行機の諸元を正しく入力しなかった場合謎の挙動で滑空してしまう（ピッチ軸の変数の代入ミスだったと考えられる）。