この記事はCanSatチームのアドベントカレンダー202312日目(12月14日)の記事です.

CanSatではなく人工衛星のお話しです！
僕が姿勢系を始めたときにこれあったらよかったなーって記事を書いてます!
今後,時間があるときにはこの内容の補足とかを追記していこうと思ってます！

間違えてる気がするよーみたいな指摘は,積極的にしてください！
むしろ,そういう所を指摘してほしくて書いている面もあります

はじめに

はじめまして.　FUSiONで地上局を担当している”とるか”(@Attitudecontro3)です！

12/25のクリスマスまでアドベントカレンダーという形で各自メンバーの活動内容や小話を毎日更新していくつもりです.
本イベントを通じてFUSiONやそのメンバーについて知っていただけると嬉しいです.

自己紹介的なやつ

自己紹介

FUSiONでは地上局を担当しています！
地上局等の詳しいお話は,「超小型衛星の通信・地上局系って何やるの?」に "ねこ見習い" かねこパイセンが書いてくれてます(^^)
(この記事に書いてるほどハードなことはやってなくて,CanSatのセンサデータをリアルタイムに可視化するwebアプリ的なものの開発をしてますん.)
詳しくは,本イベントの12/22の記事に書きます.

(12/22以降にURL追記予定)

普段やってること

帝京大学/宇宙システム研究会/衛星班:超小型人工衛星の姿勢制御系を担当しています.
(仮配属)アストロダイナミクス研究室:姿勢制御系にかかわる研究活動を行っています.
ASE-Lab.(航空宇宙機の飛行力学ゼミ,組み込みOSゼミetc...)
理系大学生の日常を書いてる人の部下として鳥人間にあさーくかかわってたりします.
FUSiON地上局

人工衛星の姿勢制御系について

前提として、人工衛星のお話をするときにはその衛星が「どんな目的(ミッション)を持っているのか」を気にする必要があります.そのミッションによって姿勢,構造,電源,通信といった各サブシステムの設計が決まってきます.なので,姿勢制御系についても例外にもれずミッションに依存することになります.(正確に言うと,その他にも宇宙活動法やロケットとのインタフェースなど多くの制約をうけることになる.)ここでは、個別具体の衛星システムの話ではなく共通する「力学」のお話しや、センサ/アクチュエータ等の「制御」のお話しもできればいいかと思います！(何なら人工衛星に限らないまである?)

人工衛星の姿勢(回転や向き)を数学的に表現する方法

制御対象を自分の思い通りに動かすためには,対象を数式でモデリングして定量的に知ることが必要になります.姿勢制御でもやることは同じなのですが,モデリングの対象が「3次元の回転運動」ということで少々複雑になります。

ベクトルのお話し1～座標系(基底)によって成分が変わる～

ここでは,人工衛星の位置を例に考えてみましょう.
人工衛星の位置を表現する座標系として「ECI(Earth Centered Inertial Frame)」と「ECEF(Earth Centered Earth Fixed)」が存在します.
ECIは地球の自転を考慮しない座標系で,ECEFは地球の自転を考慮する座標系です.
人工衛星の位置を表す「位置ベクトル」はECI座標系でもECEF座標系でも同じ位置ベクトルではありますがベクトルの成分は異なります.
何が言いたいかというと,ベクトルの演算は「座標系(基底)」が同じであることを前提に行います.
したがって異なる座標系のベクトルの成分で演算を行った場合意味のない演算になることがあります.
特に,航空宇宙工学では機体座標系や局所水平座標系,風軸系,慣性主軸系など様々な座標系を用いて演算することになるためこの点は特に注意が必要です.
次の節で,座標系(基底)の変換のお話をします.

ここでは,以下のような規則で「ベクトル」「座標系(基底)」「成分」を区別します.

$\vec{r}$:ベクトルそのもののことを指す.(座標系を定義しない)
$\{R\}$:ベクトルを表現する座標系(基底)を表す.
$\boldsymbol{r_R}$:$\{R\}$座標系で表した$\vec{r}$

この規則を,ECI/ECEFの例に適用すると
「人工衛星の位置ベクトル$\vec{r}$」を$\{{ECI}\}$で表すと$\boldsymbol{r_{ECI}}$となり
$\{{ECEF}\}$で表すと$\boldsymbol{r_{ECEF}}$となる.
したがって,$\{{ECI}\}=\{{ECEF}\}$の時を除いて$\boldsymbol{r_{ECI}}$と$\boldsymbol{r_{ECEF}}$は異なる.

座標変換と姿勢の表現方法(方向余弦行列/オイラー角/クォータニオン)

先の節ではベクトルのお話をしました.ここでは,姿勢ってそもそも何?どうやって表現するの?っていうお話をします.
物体の姿勢とは,「ある座標系と物体に固定した座標系との関係性」を表すもので,その実態は物体に固定されたそれぞれ直交な3次元ベクトル$\{C\} = \{\vec{i},\vec{j},\vec{k}\}$であらわされます.
この$\{C\}$が,機体座標系の基底であり,機体座標系から慣性座標系座標系への座標系(基底)の変換と同値になっています！つまり、姿勢の実態は3次元ベクトルの座標変換行列やそれに準ずる表現であると考えられます.
座標変換の話は,僕がするまでもなくいろんな人がいろんな記事をわかりやすくまとめています.以下にURLを張っておくので,ご覧ください.

上記の資料を読めば,方向余弦行列/オイラー角/クォータニオンについてわかると思います.
今回は,姿勢の表現方法としてはクォータニオンを使って考えていきます.

クォータニオンの表現方法に宗教があり,パラメータの順序が違うことがあります.
ここで間違ってて実装が止まってた時期がありました！

ただ,座標変換の行列や演算を扱う上で常に意識して注意しなければならないことは「その変換が何から何に対する変換なのか」ということです.
回転行列を用いて「ベクトル」を変換する場合と,「座標系(基底)」を変換する場合がありこれらは似て非なるものなので区別をつける必要があります.
この辺の話は高校レベルから始める! やさしくわかる線形代数をよんで納得しました.

ベクトルのお話し2～動座標系上のベクトルの時間微分～

宇宙機を観測する基準の座標系を$\{R\}$,物体の重心を原点として物体と一緒に運動する座標系を$\{B\}$とします.
$\{B\}$は並進運動及び回転運動を行うことができます.速度$\vec{v}$で動き,角速度$\vec{\omega}$で回転しているものとします.
$\{R\}$を慣性座標系、Bを機体座標系と呼ぶことにします。
慣性座標系Rから機体座標系Bを指す位置ベクトルrとします.
これを時間で微分すると速度ベクトルvになります.
ここで注意したいのは,このベクトルvは機体座標系と慣性座標によって成分が異なるということです.
ということは,微分の方法も慣性座標系と機体座標系で異なると考えられます.
慣性座標系での微分はみなさんご存じの方法と同じで,成分を微分すればokです.
一方,動座標系による微分は成分そのものの微分と座標系が動く分を考慮する必要があります.
細かい内容は省略しますが,結果として微分は以下のように定義できます.

\dfrac{{d}\vec{r}}{dt} = \dfrac{{\delta}\vec{r}}{\delta t} + \vec{\omega} \times \vec{r}

これは「ベクトル」の微分の話なので座標系の成分の計算に落とし込むと以下のようになります.

\dfrac{{d}\boldsymbol{r_b}}{dt} = \dfrac{{\delta}\boldsymbol{r_b}}{\delta{}t} + \boldsymbol{{\omega_b}} \times \boldsymbol{r_b}

右辺の左側の式が成分そのものの変化を表し,右側の式が座標系が動く効果を表します.
これを

回転座標系における運動方程式の導出を書いているサイトを以下に貼り付けておきます.

剛体の回転運動の式や物理量について

ここからは,剛体の回転運動式やそれに関連する物理量について紹介していきます.
添え字bとrはそれぞれ,b:機体(body)座標系の成分で表していることを示し,r:慣性座標系(reference)の成分で表していることを示しています.

姿勢${\boldsymbol{q}}$

先の節に挙げた座標変換行列やそれに準ずる表現
慣性座標系に対する機体座標系の姿勢を表す
※変換としては機体座標系から慣性座標系※

角速度$\vec{\omega}$

姿勢を微分した時に現れる3次元ベクトル
慣性座標系に対し機体座標系がどのように回転しているのかを表す

慣性テンソル$\boldsymbol{J}$

慣性テンソルとも呼ばれる3次の正方行列
質量特性という呼ばれ方もされ,これは物体の質量分布や形状に依存する回転のしやすさなどの指標
慣性テンソルは座標系の取り方により値が変化するため注意が必要
機体座標系でとったときの慣性テンソル成分の慣性行列は$\boldsymbol{J_b}$と書く
テンソルの座標変換は特殊(慣性行列を回転行列とその逆行列で挟み込む形)なので注意が必要
これが対角行列になるような座標系の取り方を「慣性主軸座系」と呼ぶ

角運動量$\vec{L}$

慣性テンソルと角速度ベクトルの行列積からなるベクトル.
$\vec{L} = \boldsymbol{J}\vec{\omega}$

トルク$\vec{\tau}$

角運動量を微分すると現れるベクトル
慣性座標系で微分すると,以下のようになる
$\vec{\tau} = \dfrac{d{\vec{L}}}{d{t}}$

今後の節では,これらを駆使して姿勢の更新を表現していきます.

回転運動方程式による角速度の更新

姿勢の運動を計算するためには,物体の姿勢を数学的に表現する座標変換行列やクォータニオンを角速度を使って更新していく必要があります.
姿勢の更新は「機体座標系から見た各速度成分」を用いるため,角速度の更新は機体座標系にて行います.
(慣性座標系で更新する方法もあり、くわしくはこちらで紹介されています.)
先に紹介した動座標系での角運動量ベクトルの微分は以下のように表すことができます.

\vec{\tau} = \dfrac{\delta{\vec{L}}}{\delta{t}} + \vec{\omega} \times \vec{L}

ここで,この式を機体座標系で表現すると

\boldsymbol{\tau}_{b} = \dfrac{\delta{\boldsymbol{L}_{b}}}{\delta{t}} + \boldsymbol{\omega}_{b} \times \boldsymbol{L}_{b}

角運動量を変形すると

\boldsymbol{\tau}_{b} = \dfrac{\delta{\boldsymbol{J_b\omega}_{b}}}{\delta{t}} + \boldsymbol{\omega}_{b} \times \boldsymbol{J}_{b}\boldsymbol{\omega}_{b}

慣性テンソルは座標系に依存するが,機体座標系は「機体に固定された座標系」であるため,これを定数と置くことができる.したがって.

\boldsymbol{\tau}_{b} = \boldsymbol{J_b}\dfrac{\delta{\boldsymbol{\omega}_{b}}}{\delta{t}} + \boldsymbol{\omega}_{b} \times \boldsymbol{J}_{b}\boldsymbol{\omega}_{b}

式変形して,機体座標系の角速度の微分の項を整理すると,

\boldsymbol{J_b}\dfrac{\delta{\boldsymbol{\omega}_{b}}}{\delta{t}}= \boldsymbol{\tau}_{b}  - \boldsymbol{\omega}_{b} \times \boldsymbol{J}_{b}\boldsymbol{\omega}_{b}

つまり,機体座標系において慣性座標系に対する機体の角速度ベクトルの微分方程式は以下のように記述できます.

\dfrac{\delta{\boldsymbol{\omega}_{b}}}{\delta{t}} = \boldsymbol{J_b^{-1}}\{\boldsymbol{\tau}_{b}  - \boldsymbol{\omega}_{b} \times \boldsymbol{J}_{b}\boldsymbol{\omega}_{b}\}

上記ダイナミクス方程式を差分方程式に変えて,後述するキネマティクス方程式と連立して数値解析を行うことにより剛体の回転運動のシミュレーションを行うことができます.

角速度による姿勢の更新

今回,姿勢はクォータニオンであらわします.
結果だけ述べると,以下のようになります.

クォータニオンの表現方法に宗教があり,パラメータの順序が違うことがあります.
ここで間違ってて実装が止まってた時期がありました！

こちらは,$\boldsymbol{q}=\{w,x,y,z\}$の順番でクォータニオンを表す場合です.

{\dot{q} =
\frac{1}{2}
\begin{pmatrix}
0 & -\omega_x & -\omega_y & -\omega_z \\
\omega_x & 0 & \omega_z & -\omega_y \\
\omega_y & -\omega_x & 0 & \omega_z \\
\omega_y & \omega_z & -\omega_x & 0
\end{pmatrix}
\begin{pmatrix}
q_w \\
q_x \\
q_y \\
q_z
\end{pmatrix}
}

こちらは,$\boldsymbol{q}=\{x,y,z,w\}$の順番でクォータニオンを表す場合です.

{\dot{q} =
\frac{1}{2}
\begin{pmatrix}
0 & \omega_z & -\omega_y & \omega_x \\
-\omega_z & 0 & \omega_x & \omega_y \\
\omega_y & -\omega_x & 0 & \omega_z \\
-\omega_x & -\omega_y & -\omega_z & 0
\end{pmatrix}
\begin{pmatrix}
q_x \\
q_y \\
q_z \\
q_w
\end{pmatrix}
}

以下のような記事を参考にコーディングしています.

姿勢制御のアルゴリズム

先の章までを差分方程式等に直し,ルンゲクッタ法等を適応して数値計算していくことにより姿勢計算ができるようになります.それでは,姿勢制御は具体的にどのようなアルゴリズムで行うのでしょうか?ざっくり紹介していきます！
まず,宇宙機等の制御を考えるときGNC系というものがあります.これは【Guidance, Navigation and Control】の大文字の部分を取ったところで日本語で言うと航法・誘導・制御となります.
それぞれ各部を簡単に紹介すると

航法:位置や姿勢を決定する系
誘導:目標値と現在値の比較を行って出力を決定する
制御:決定した出力に基づきアクチュエータの操作を行う
といった具合でです.

ここから,それぞれどんなセンサやアクチュエータ等が用いられるのかを紹介します.

航法(姿勢決定)

スタートラッカ:イメージセンサ等で撮影した画像と星座カタログを比較して姿勢を決定する.
太陽センサ:太陽方向ベクトルを決定できる.
地球センサ:赤外線放射とうから地球の方向ベクトルを決定できる.
磁気センサ:地球の磁場の計測を行うことにより姿勢決定に用いる
ジャイロセンサ:機体座標系から観測する,慣性座標系に対する角速度を計測して積分して姿勢を求める.

これらを組み合わせてIMUなどといったユニットを構成するようですが,以下のサイトに詳しくまとめられているのでここでは省略します.

よく耳にするカルマンフィルタや相補フィルタというような技術がこれに該当すると考えています.
また,これらの詳細な制御側に関しては,以下の論文(東京工業大学,たぶんいつかの宇宙科学連合講演会か何かの資料?)に詳しい記載があります.

誘導/制御

誘導および制御はアクチュエータにより毛色が異なるものであると考えています.
姿勢制御アクチュエータしては,以下のようなものが例として挙げられると思います.

スラスタ(推進系)
リアクションホイール
コントロールモーメンタムジャイロ
磁気トルカ
永久磁石
重力伸展ブーム
膜展開

上から4個は能動的なアクチュエータであると考えることができます.
これらは単体で使われるようなことはあまりなく,たいていは組み合わせて使います.
例えば,おそらく姿勢制御で最も多く使わているのはリアクションホイールですがこれには回転数の飽和という大きな欠点があります.なので,この欠点を補うためにスラスタや磁気トルカなどを利用することがあります.

制御手法について

スラスターについて

リアクションホイールについて

ISASニュース(宇宙科学研究所1994.4 No157)P.9配置の話が載っている.

リアクションホイールでの運用の工夫の話

コントロールモーメンタムジャイロに関して

*　磁気トルカに関して

こんな感じで制御則のことを考慮するほかに,宇宙環境の外乱も考慮する必要がありますが,今回は省略します.

数値計算を行う上での課題等(悩んでいること)

僕自身,衛星の検討を行う中で課題になっていることが「数値計算の妥当性の検証や担保の方法」です.
指導教員やその他のメンバーに指摘されるまで気づけなかったりすることがあるのでこの辺をもうちょっと考えなきゃいけないなーなどと思っています.

まとめ

長くなりましたが,今のところ僕が考えてるのはこういうところです.
これでもまだ姿勢制御のロジックの部分のうわべをなぞっただけで,もっと深堀しようとすると現代制御理論やモデル予測制御みたいな話をし始めることができたり,回転表現が云々って話に持って行けたりもします.でも,結局実際の衛星を実装するときにもっと必要になるのは,仕様書を読んで自分たちの要求に沿ったものを作り上げる力になります(自作の場合は,要求に沿ったものまたは最適なものを作る力).自分たちの要求を作り上げたり,それを検証するためにも,この記事の内容はわかってると進めやすいのかなと思ってたりします.
この記事が「姿勢系やりたい！」って人の参考になったり,「姿勢系わからん！けどやらなきゃ！」みたいな人の手助けになればなーと思い書きました.

参考