ロケットシミュレータとQuaternion演算

こんにちは。yuukivelです。今回は飛行機の話題からちょっと離れて、飛翔体（ロケットとか）のシミュレータを作った（作っている）ので、この公開と、この中で用いた座標変換処理におけるQuaternion（クオータニオン：四元体）について解説する。

まあとりあえず先にダウンロードリンクを載せておこう。Octave スクリプトで書いているが、MATLABでも動くし、PythonのGUI版も作りかけだったりする。
（使用して損害が生じたとしても責任は負いかねますのでご了承下さい。）
mq_rocket zipファイルダウンロード
 githubはこちら
https://github.com/NaotoMorita/mq_rocket/

パラメータをいじる場合はrocket_parameter.m内の変数の値を変える。

本当はPythonでGUIを積んだアプリケーションの形で公開するつもりだった。そうでないのは、言い訳をすると、

Numpyでの演算で生じる謎の微小な周期的数値振動
グラフ表示のライブラリであるMatplotlibのスパゲッティーコードとイケてNASA。３次元プロットの軸まわりの処理が残念。

の二つ。克服したら公開します。最終的にはOpenGLやDirectX等を用いて３次元のモデル表示まで持っていく予定。

さて、今回はこのスクリプトの中で用いた四元数、すなわちQuaternionについて紹介する。

まずあなたが三軸のジャイロセンサーを手に入れて、値も綺麗に取れたとしよう。そうしたら何をするだろうか？
当然、角度を求めようと積分するだろう。

ではその積分はどうする？X軸、Y軸、Z軸で積分するだろうか。当たり前だと思う人もいるだろう。
否、それは大間違いなのである。

ジャイロが存在する座標系は「動座標系」であり、私たちが知りたいのは、ある地面から固定された座標から見たジャイロセンサーの傾きの角度、すなわち「オイラー角」なのである。

オイラー角は固定座標系、すなわち「慣性座標系」からみた動座標系の傾きを示すものであり、それ自体が回転を定義する。
このため、オイラー角から、３次元座標の回転を表すことができる。これをDCM(Direction Cosine Matrix)という。

さて、このオイラー角と回転行列を用いて飛翔体の運動方程式を記述した際の長所と短所をまとめておこう。

オイラー角で運動方程式を記述した場合

長所

わかりやすい。定義も自然なので、おかしければすぐに分かる。後処理も簡単。

短所

三角関数、逆三角関数を多用するので、計算が重くなりやすい。これはマシンパワーが無いマイコンとかで積分するときには致命的
そもそも見かけの力を定義する式が煩雑でわかりにくい。
ジンバルロックを起こす。

DCMで運動方程式を記述した場合

長所

計算が楽

短所

姿勢更新を行うことができないので、DCM単体では「そもそも運動方程式が記述できない」（！）
単なる行列なので、見ても何が何だか分からない。

ここで登場するのが、QuaternionとDCMを組み合わせた運動方程式の記述なのだ。

QuaternionとDCMを用いて運動方程式を記述する

Quaternionを用いると以下の事ができる。

代数演算だけで回転行列が記述できる。
角速度を用いたQuaternion更新の公式がある。
逆Quaternionを用いれば逆回転行列を求めるのも簡単

こうした特徴の御陰で、軽く、精度の良いシミュレータを実現できる。飛行機の運動方程式も、XFOIL.mとの組み合わせてやれば、微小擾乱といった近似を外して軌跡を求めたりできる。（というより、最適化理論と組み合わせて設計に使えるものをいま作ろうかなとか思ってる。話変わるど、グラディウスを遺伝的アルゴリズムで解く（？）とかすごいよね）

以下にクォータニオンを用いてロケットの運動方程式を書く際のフローチャートを載せる。そのまんまmq_rocketのフローチャートである。

この中で出てくる速度座標系とは機体座標上での速度のX方向にX軸を合わせて設定する座標系。揚力や抗力はこの座標系の軸方向に働くので、拘束で高精度な力見積もりが可能だ。（ただしCLの求め方とかざっくりだけど。）