ロケットの飛行解析ソフト（OpenTsiolkovsky）をPythonに移植してみた - 概要詳解 -【OpenTsiolkovsky(Python)】
1. はじめに
こんにちは。この記事は、インターステラテクノロジズ株式会社様がオープンソースで公開されているロケットの飛行解析ツール「OpenTsiolkovsky(https://github.com/istellartech/OpenTsiolkovsky)」をPythonに移植してみたということで、
- ツールの紹介・何が出来るのか
- どうして移植しようと思ったのか
みたいなことを解説していきたいと思います。
プログラムは以下GitHubの方にあげておりますのでご確認くださいませ。
（オープンソースソフトウェアの公開は筆者としては初めての経験のため、もし拙いところなどありましたらコメントなどでお知らせ頂けますと大変有難いです。）
2. 飛行解析って何？
恐らく宇宙業界、特にロケット開発にあまり縁がない人にとっては、「飛行解析」と聞いても余りピンとこないのではないかと思いましたので、まずはじめに簡単なご紹介をさせて頂きたいと思います。
ロケットの飛行解析とは読んで字のごとく、「ロケットを実際に飛ばしたときの動き・飛行経路等を事前に解析して確かめておくこと」を指します。ロケットの打上げは一発勝負の大事ですから、出来る限り事前にテスト・試験・検証等を行っておくことが非常に重要です。 ロケットの構造体やアビオニクス等のハードウェアは、振動試験や熱真空試験等、実際に物を使用して試験を行うことも可能ですが、それらをロケットとして一つに組み上げたときの挙動を試験等で確かめておくことは難しい場合が多いです。このように、事前に確かめておくことが難しい部分の挙動等を確認する手段として「解析」があります。
どこまでの現象をどの程度の精度で解析するか等、解析結果の使われ方や目的なんかによってツールに求められる性能も変わってきます。 高機能・高精度なツールがあればそれだけで良いかと言うと意外にそうでもなく、色々な試案があってクイックに検討を回してみたいときや、同時に複数のケースを回すことで最適化を試みたいとき等は、軽くて使い勝手の良いツールでまずは当たり付けをすることもあります。
他の多くのツールと同じく、飛行解析に用いるツールも適材適所が重要です。
3. ツールの紹介・何が出来るのか：OpenTsiolkovskyについて
さて、ご紹介した飛行解析に関わるツールについてですが、比較的簡単に使えるものと言えば、モデルロケット用の解析ツール等があるでしょうか（例：OpenRocket）。 しかし本当のロケットの打上げに使うような本格的なものについては、機密情報等を含んでいる場合もあるため国内外含めて中々自由に使えるツールが見当たらないという問題があったりします。
そんな折、「（某トラ技の特集に）オープンソースで開発されている飛行解析ツールのことが載っているらしい」という情報を教えてもらう機会がありました。すぐに書籍を購入してソースコードがアップロードされたGitHubを見に行ったことを覚えています。それが今回解説する「OpenTsiolkovsky」になります。
あらためての紹介となりますが、OpenTsiolkovskyは、インターステラテクノロジズ株式会社様が作成されたロケットの飛行解析用のオープンソースソフトウェアです。 このツールでは、入力された情報（ロケットの構成、質量、推力や姿勢履歴等）を基に、運動方程式を解いてロケットの飛行経路を解析することができます。
具体的なイメージとしては、例えば以下の図のように飛行経路を求めることが出来ます。 （細かい部分に差異はありますが、画像は日本のM-Vロケット１号機の飛行経路を解析により再現したものになります（本記事の後半でより詳細に解説しています）。）
このように、飛行解析を行うことでロケットがどのように飛んでいくか、その挙動等がある程度事前に分かるようになります。
どうでしょう、何だかワクワクしてきませんか？
実際にロケットを一機打上げるのは非常に大変なことですが、解析であればPCのボタン一つで完了します。色々なシステムの構成や可能性等を事前に何度でもシミュレート出来るというのは、個人的にはかなり好奇心を刺激される楽しいことだと思っています！
4. Pythonへの移植と本プログラム(OpenTsiolkovsky(Python))について
今回、「OpenTsiolkovsky」をPythonに移植しようと思ったきっかけは、まずは純粋な勉強目的からでした。 上述の通り、本格的な飛行解析ツールは公開されている情報も余り多くなく、特にオープンソースの形で公開されているものは非常に貴重という印象です。ロケットの打上げ計画等を考える上でも非常に重要な解析ツールですが、実際のところはブラックボックスで使われていることも多い印象のものでした。 いつの日かホワイトボックス化・完全自作したいという思いがありましたが、まさに探していたものを見つけた！という気分で、すぐにソースコードの解読に取り掛かったのを覚えています。
「OpenTsiolkovsky」のコア部分はC++でコーディングされています。飛行解析ツールですが、ある事情から実行速度が非常に重要となってくるので、コア部分を動作の早いC++でコーディングすることは非常に合理的であると考えています。 今回は目的が勉強であること、そのために実行速度よりも極力可読性の方を重視したかったこと、また周辺のライブラリの充実度や将来的な連携も考えて、Pythonでコーディングすることを考えました。
今回の移植に当たっては、プログラム内のクラス設計を一からデザインし、元のコードにあった機能の再配置等も行っております。 解析にあたっては、ロケットの諸元（質量、推力、etc）、周辺の大気・風、惑星の重力等々、多くの情報が複雑に相関します。また、それぞれの計算を容易にするために頻繁に座標系の変換を行う必要があり、これも内部構造の複雑化を招く要因となっています。 今回は極力、こうした計算の処理の部分と、それによって変わるパラメータの部分を分けて設計したいと思いました。イメージとしては、ロケットはロケットの、大気は大気の、重力は重力の、といったように各々のデータをそれぞれがパラメータとして持っており、間の相互作用や計算の部分は極力それ専用のマネージャークラスに担わせるという考え方です。こうすることで、特にロケットの諸元が変わる場合でも、処理部分には変更を加えることなく、パラメータ設定の変更のみで対応できると考えました。
結果的には、役割分担を適切に出来なかったりと上手く行かない部分もあり、再度一からやり直したい気分に駆られることしばしばでしたが、一旦形にはなりましたのでここで公開しようということになりました。
こちらが本プログラムのクラス図です。
RocketStateクラスがロケットの位置・速度・加速度等の諸々の情報を保持しています。C++で言う構造体のような、パラメータのみを保存するクラスとして定義しています。
RocketStageクラスはロケットの１段目、２段目、３段目…といったロケットの各段を表しています。それぞれが自分の状態量（RocketState）を保持しており、それに応じて自分の位置や飛行状況を更新する機能も持っています。
RocketクラスはRocketStageで表現された各段の情報を束ねて１機のロケットとしての振る舞いを定義したクラスとなっており、ロケットの状況やイベント（各段の分離等）の把握、外部クラスとのIF等を担わせています。
その他、プログラム全体の実行を管理するマネージャークラスや、大気・風・重力等の情報を提供するEnvironmentクラス、実行状況を記録するLoggerクラス等から本プログラムは構成されています。 頻繁に行われる座標系変換については、当初はある基準とする座標系のデータをベースに、その他の座標系への変換を自動で行ってくれるようなクラスを設計しようかと思いましたがあまり上手くいかなかったため、座標系変換を行う関数をグローバル関数として定義し、使用するシーンで各々呼び出すこととしました。
当初の目論見は上述した通り、パラメータと処理部分の分離、それによる可読性や拡張性の向上等を狙っていましたが、蓋を開けてみるとオブジェクト指向色の強い設計となり、特にRocketStageクラスは肥大化してしまったかなと感じています。 座標系周りの機能も上述したようにグローバルからの呼び出しと言うことで整っている感じがせず、もうちょっと上手な設計が出来たかなぁという反省もありますが、本クラス図の整備を始めとして全体の見通しを良く開発することはできましたので、そこは良かったかなと思っているところです。
実際の使用例：M-Vロケット１号機の飛行解析
早速出来上がったプログラムを用いて何か試しにロケットを一機飛ばしてみたいと思います！
解析には様々な情報が必要で、中には機微な情報が含まれているものも多く、インターネットに公開されている情報から適切なサンプルを見つけてくるのは中々大変でしたが、日本の固体ロケットとしてかつて打上げられていた「M-Vロケット」に関しては比較的多くの情報が公開されていたため、M-Vロケット１号機の諸元をサンプルとして入力ファイルを作成・実行してみました。
念のため、今回解析を回すにあたっては、諸元を調べきれなかったキックステージ段(第4段目)を省略したり、機能が不安定なため模擬しなかったイベント（ペイロードフェアリングの分離等）、詳細が分からないため大まかな推定をしたところ（風向・風速）等があり、実際の諸元をすべて正しくインプットした解析となっていない点はあらかじめご承知おき頂ければと思います。
※ 参考文献を以下に紹介いたします。
■ M-V rocket F1 referenced data(folowing):
- Attitude: https://jaxa.repo.nii.ac.jp/record/33355/files/SA0200128.pdf
- Propulsion system: https://www.isas.jaxa.jp/publications/hokokuSP/hokokuSP47/85-116.pdf
- SOE: https://www.isas.jaxa.jp/j/enterp/rockets/vehicles/m-v/seq.shtml#
- Epoch, initial azimuth:https://www.isas.jaxa.jp/publications/hokokuSP/hokokuSP47/579-592.pdf
※ 解析を回すにはデジタルデータが必要ですが、今回のサンプル作成にあたっては、上記の参考文献からGraphcelというソフトを用いてデジタルデータ化したものを使用しました。
さて、解析により色々な情報が分かるようになりますが、その一例として冒頭にも上げた飛行経路があります。これはロケットの時々刻々の位置等の変化を示したもので、以下の通り、M-Vロケット１号機の場合は鹿児島にある内之浦の射点から打上げられた後、東向きに飛行して太平洋の海上あたりで着水する結果となりました。
念のため、基にした「OpenTsiolkovsky」の方でも解析を行い、本プログラムの結果と比較してみます。
少し分かりづらいですが、本プログラムの出力(先ほどの画像の白線)と元のOpenTsiolkovskyから出力されたデータを重ねて表示しています。概ね線は重なっており、このレベルで見て分かるような差分等がないことが分かりました。まずは移植成功、といったところでしょうか。
さて、実際の軌道との比較はどうでしょうか？ 先ほど紹介した参考文献(1)の資料にM-V-1号機の飛翔経路の図が載っておりましたのでそれと重ねて比較してみたいと思います。
上の図と縮尺等を合わせて、本プログラムの解析結果からExcelで作成した図が以下となります。
ここからは若干パワープレイとなりますが、これら二つの図を重ねて差分をチェックしてみようと思います。 まずは高度から。
これは分散（正常な誤差の範囲と言える領域）の高い側ちょうどギリギリといったところでしょうか。間違っているとまではいえませんが、そこそこの差分が出たなぁという印象です。
次は緯度・経度で見た経路を見て見ましょう。
こちらは予測の中心（”ノミナル”ということが多いです）とそれなりに近い結果が出たでしょうか。
以上、所感としては、上述の通り本来であれば分離しているはずのフェアリングを(解析上)つけたまま飛ばしていたり、色々な諸元を推測で入力していたりとかなりアバウトな条件をインプットした割には、良い精度で解析ができているような感じがしました。 最終的に解析が正しかったかどうかは実際に実機を打上げてみてその結果と比較してみるまでは何とも言えない部分も多いのですが、オープンソースで得られる情報からここまで推測できるというのはとても面白いと思いました！
おわりに
ここまでお読みくださりありがとうございました。 少しでも宇宙開発またはロケット開発にご興味を持って頂けましたら筆者としてはこの上ない喜びです。
また、重ねてにはなりますが、素晴らしいソフトウェアをオープンソースという形で公開してくださっているインタステラ―テクノロジズ株式会社様にはこの場を借りて、あらためてお礼申し上げたいと思います。
ではでは。
