磁気トルカによる姿勢制御～B-dot則について～

# 1.はじめに
まずはじめに、いろいろボロボロです。。。
卒研と衛星開発に追われながら書いているのでご容赦ください(甘え)
書き直すのでﾕﾙｼﾃ。。。。

## 1.1 目的
こちらは[CubeSat Advent Calendar 2024](https://qiita.com/advent-calendar/2024/cubesat)の3日目の記事として書きました！
みなさんにB-dotについて知ってもらえればと思います.

## 1.2 対象者層
- **人工衛星を作っている人/作ってみたい人/興味がある人**

一応,初心者向けのつもりですがそこそこやっている人も面白く読めるかと思います.
初心者向けの章も作るので,知ってるよってところは読み飛ばしてくれればokです！
(書いている途中で思い出しましたが,この資料B1の時にほしかったなって感じの資料になっています.具体的には,衛星開発にぶち込まれて右も左もわからないけどとりあえず実装をしなきゃいけないときです.)

## 1.3 背景
磁気トルカによる制御則を実装しようとしてちょろっと調べると「どうやらメジャーな手法として**B-dot則**というのがあるらしいな」ということが分かると思います.
調べる量を増やしていくと,いろいろなところでB-dotを使っていそうな感じがしますが,なんかこれおかしくね？ということにぶち当たります.
そう,**公開されている資料によって定義が違う**のです.(もやもやポイント)
これについて,いろいろ整理してみましたので参考になれば幸いです.

# 2.人工衛星の姿勢制御

## 2.1 姿勢制御とは

「姿勢(Attitude)」と聞くと,[一般的](https://kotobank.jp/word/%E5%A7%BF%E5%8B%A2-519836#goog_rewarded)にはこんなイメージがあるんじゃないかと思います.
>１ からだの構え方。また、構え。かっこう。「楽な姿勢で話を聞く」
２ 心構え。態度。「政治の姿勢を正す」

(余談：筆者は [Mrs. GREEN APPLEさんの「Attitude」](https://www.youtube.com/watch?v=ji99kh-GzTQ)って曲のタイトルを見た時ちょっとドキッとしました(この曲の場合は,2の意味のほうが近そうですね))

ここで,僕がお話しする姿勢とは「人工衛星の向きや回転」のことを指します！
特に,今回は「回転」のお話をしていきます。
(姿勢制御そのものの記事については筆者が過去に書いた[この記事](https://qiita.com/Attitudecontro3/items/4ca0e3a2465987b1a591)を読んでくれたらうれしいです.)

~~戯言が長かった気がしますが~~姿勢制御とは**向きや回転を思いのままに操ること**を指します.

# 2.2人工衛星の姿勢制御の必要性

人工衛星は基本的に何らかのミッション(目的)をもって宇宙空間である軌道上で動作します.ミッションは通信であったり観測であったり様々です.
ここで,[アクセルスペース「衛星システム技術の裏側note（第４回）」](https://note.com/axelspace/n/nf8cee7119810)から引用すると
>衛星搭載のレンズが軌道上で１度方向がずれると地上では約9.5kmずれてしまいます。

ということで,ミッションによってはめちゃくちゃコアなシステムになるわけです！

また、姿勢制御を失敗してとんでもなく衛星が回ってしまうと,回転によって生じる力が衛星を破壊するということにもなりかねません(恐怖)
そんな感じで,回転の速度を減少させる制御も衛星の生存にとって重要だったりします.

## 2.3姿勢制御系の構成要素

**能動的な**姿勢制御系の構成要素は以下であると**考えています**.
- センサ部
    - 姿勢を測定するためのセンサ(例)
        - ジャイロセンサ
        - 磁気センサ
        - 太陽センサ
        - スタートラッカ
        - GNSS(軌道決定)
- 演算部
    - センサから得たデータから姿勢決定を行う(カルマンフィルタなど)
    - 目標姿勢と現在姿勢の幾何的差分から偏差を計算する
    - 偏差からアクチュエータへの入力を計算
        - マイコン
        - FPGA
- アクチュエータ部
    - トルク(回転の力:単位[Nm])を発生させるための装置(例)
        - コントロール・モーメンタム・ジャイロ(CMG)
        - リアクションホイール(RW)
        - リアクションコントロールシステム(RCS)
        - 磁気トルカ(MTQ)

ということで今回はこの中から磁気トルカ(MTQ)と,そのの制御則であるB-dotに関するお話をしていきます！

# 3.磁気トルカ(MTQ)について
## 3.1.磁場の干渉によるトルクの発生
![スクリーンショット 2024-11-23 034706.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635973/ca0d4dfe-f96f-8a26-0fb3-03256a083a75.png)

- このように磁場ベクトル(地磁気)と磁気モーメントが一致させるようなトルクが発生します！
    - 磁場ベクトルは地球が持つ磁場のことを指します
    - 磁気モーメントは磁石の磁力の強さの指標のようなものです
- 方位磁石が磁北を指し続けるのと同じ原理ですね.
    - コラム
        - 並進運動に置き換えると「ばね」のような働きをします
        - 理想のばねを考えると単振動を行いますが...
        - 方位磁石が単振動にならないのは,摩擦(ダンパ項)があるからです！

## 3.2.発生磁場の制御
- 磁気トルカの実態は「コイル」です！
- コイルでどのように姿勢制御するのか？
    - コイルに電流を流すと,磁石になります(義務教育)
    - コイルに流す電流を増やすと,磁石の強さが大きくなります(高校物理)
    - 中に芯を入れたりすると磁場が強くなったりします(歓喜)

## 3.3.数式でモデリング
```math
\vec{\tau} = \vec{M} \times \vec{B}
```
$\vec{\tau}$:磁場と磁気モーメントによって生じるトルク
$\vec{B}$:地球磁場などの磁場ベクトル
$\vec{M}$:磁気トルカが発生させる磁場ベクトル
- ベクトルの外積(ベクトル積)で定義されます
    - ただでさえ姿勢のダイナミクスやキネマティクスは非線形...
    - なので三軸制御とかやろうとするとなかなかしんどいです....(今やってる)

# 4.B-dot則について
## 4.1.B-dot則の基本的な考え方
- フィードバック制御である
![スクリーンショット 2024-11-23 040651.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635973/2f3cbc61-0d95-2607-a0ee-8dc9b931b35b.png)

    - これからB-dot則の深みについて説明していくのですが,基本的に抑えておきたいのはここです.
    - 要するに,目標があってそれに近づけるような制御をしているということ.
- 角速度を抑制する目的で動かす制御則である
    - 経験上,B-dotを使うのは何らかの理由で角速度を減らしたいとき.
        - リアクションホイールの飽和を解消
        - 衛星の初期回転を抑制
    - 逆方向に作用させると角速度を増加させることも可能ではある

## 4.2.磁気差分フィードバックB-dot
プロイテレスとか,その辺で紹介されるやつ
![スクリーンショット 2024-11-23 052302.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635973/e8703f47-c101-bfbd-9fa4-d7c7eae5d3a4.png)
- アルゴリズムの重要部分
```math
M= -k_d (B_{[2]}-B_{[1]})
```
磁気センサの時刻差分を出力する制御則になっている。
- 参考資料
    - [大阪工業大学・電気推進ロケットエンジン搭載
小型スペースシッププロジェクトの進捗状況 p.18](https://unisec.jp/history/ws2009/UNISEC09-21show_Osakako.pdf)
    - [第22回衛星設計コンテスト 設計の部 光通信技術実証キューブサット "OPT-CUBE" 東北大学 工学研究科 航空宇宙工学専攻 修士 1 年
○ 藤田伸哉，郷内稔也
](https://www.satcon.jp/history/prize22/pdf/22_asj.pdf)

## 4.3.角速度フィードバックB-dot
日大の論文ででてくるやつ
![スクリーンショット 2024-11-23 052435.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635973/903ce4e9-1bd7-a870-cc4f-e7dbce885856.png)
- アルゴリズムの重要部分
```math
M= -k_d (B \times \omega)
```
磁気センサの値とジャイロセンサの値の外積から計算される.
- 参考資料
    - [超小型人工衛星SPROUTの姿勢制御系の開発 指導教授 宮崎康行 M2014 田中 惇](https://stage.tksc.jaxa.jp/taurus/member/miyazaki/old/publication/M_2013_tanaka_atsushi.pdf)
    - [東京科学大学リサーチリポジトリ「可変形状を用いた宇宙機の姿勢制御と軌道上性能評価」](https://t2r2.star.titech.ac.jp/rrws/file/CTT100908077/ATD100000413/20D10114_Kei_Watanabe_thesis.pdf)

# 5.比較および考察
## 5.1. シミュレーションによる比較
### 5.1.1.シミュレーション条件
- 衛星構造のパラメータに関する計算条件
    - 3U CubeSatを仮定して以下の条件でパラメータを設定した.

| パラメータ        | 値                | 単位           | 備考                          |
|--------------------|-------------------|----------------|-------------------------------|
| 質量 (\(m\))      | 4.0              | kg             | CubeSat全体の質量             |
| 辺の長さ (\(a\))  | 0.1              | m              | 横方向の長さ                  |
| 辺の長さ (\(b\))  | 0.1              | m              | 縦方向の長さ                  |
| 辺の長さ (\(c\))  | 0.3              | m              | 奥行き方向の長さ              |
| 分布の仮定         | 均一              | -              | 質量は均等に分布していると仮定 |

- 慣性モーメント計算結果

| 慣性モーメント軸   | 式                                  | 計算結果       | 単位           |
|--------------------|-------------------------------------|----------------|----------------|
| $\( I_x \)$         | $\( \frac{1}{12} m (b^2 + c^2) \)$   | \( 0.033 \)   | kg·m²          |
| $\( I_y \)$         | $\( \frac{1}{12} m (a^2 + c^2) \)$   | \( 0.033 \)   | kg·m²          |
| $\( I_z \)$         | $\( \frac{1}{12} m (a^2 + b^2) \)$   | \( 0.0067 \)  | kg·m²          |

- 発生磁気モーメント
    - [JAXA-SMASH 超小型衛星利用シンポジウム 2022 「CubeSat向け高精度姿勢制御ユニットの開発」五十里哲①、中村博一②、船瀬龍③
①東京大学、②セーレン、③JAXA](https://aerospacebiz.jaxa.jp/jaxa-smash/cubesatlv2022/6-3_detail_technology.html)
    - 上記のページを参考に磁気モーメントの最大値は各軸0.37 [Am^2]と仮定した.

- 衛星投入軌道
    - 革新2号機のうちの1つである「多目的宇宙環境利用実証衛星 TeikyoSat-4(OORURI)」のTLEを使って軌道計算を行う
        - 補足：OORURI君は50x50x50[cm^2]の50[kg]級衛星なので3Uではないのですがそこは大目に見てくださいー
    - [革新2号の衛星たちTLE](https://celestrak.org/NORAD/elements/table.php?INTDES=2021-102&SHOW-OPS)
    - [2024年12月01日に参照した値を用いる](https://celestrak.org/NORAD/elements/gp.php?CATNR=49396)
    - TLEの読み方に関しては調べると意外と出てきます
        - そのうちTLEの軌道6要素とケプラリアン6要素の違い的な記事も書きたい
>TEIKYOSAT-4 (OORURI)    
1 49396U 21102B   24335.84608261  .00005298  00000+0  36571-3 0  9997
2 49396  97.4388  15.0419 0014666  70.0641 290.2165 15.06099381167824

- 軌道伝搬
    - TLEを伝搬するアルゴリズムである[SGP4((Simplified General Perturbations Satellite Orbit Model 4)](https://ja.wikipedia.org/wiki/SGP4)を用いる
    - Python [Skyfield](https://rhodesmill.org/skyfield/)ライブラリを用いる

- 外乱トルク等について
    - 純粋にB-dot則の効果を確かめる目的であるため考慮しないものとする

- 地球磁場モデルについて
    - [IGRF( International Geomagnetic Reference Field )モデル](https://www.ncei.noaa.gov/products/international-geomagnetic-reference-field)を利用
    - Python [ppigrf](https://github.com/IAGA-VMOD/ppigrf)を用いる

- 座標変換について
    - [緯度経度、ECEF系、方位角・仰角などの座標変換に便利なプログラムPyMap3d](https://qiita.com/ina111/items/6e3c4d85036fd993d23c)
    - 磁場モデルの座標変換等に利用

## 5.1.2 ソースコードに関して
- リファクタリングでき次第公開します(しばらくの間は勘弁してください)
- パラメータチューニングもできなかったので汚いデータですがご容赦ください...

## 5.1.3 シミュレーション結果の比較
- **磁気差分フィードバックB-dot**
    - 角速度ノルムの時間履歴
        - 減衰しつつある様子が見て取れる
![スクリーンショット 2024-12-02 052858.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635973/032b52d7-8097-f44d-c0eb-53dd36ead29a.png)
    - 機体座標系の磁場ベクトルの時間履歴
        - sat_mag_sen\[1](y)軸が最も大きい状態を維持している
        - 磁場に対する変化が一定になっている
![スクリーンショット 2024-12-02 023614.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635973/79125118-4e9e-1aba-0fb4-c914b806073b.png)

    - 姿勢履歴Gif
        - 緑の実線矢印と黒の実線矢印が時間経過したときに一致していく
![gif-1](https://drive.google.com/uc?export=view&id=18N5b3AwHljQe6lm8-rGITUYTguUWYQoC)

- **角速度フィードバックB-dot**
    - 角速度ノルムの時間履歴
        - 時間経過で角速度ノルムが減衰している
![スクリーンショット 2024-12-02 053158.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635973/dc3b8187-4c47-7f62-4d42-a9f4b8b0920a.png)

    - 機体座標系の磁場ベクトルの時間履歴
        - 磁気差分フィードバックの時とは異なり値は変化し続ける
        - 磁場に対して一定の姿勢になっていない
![スクリーンショット 2024-12-02 053234.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635973/c2e04cd8-ce44-2b8f-a56f-127a6a49ae68.png)

    - 姿勢履歴Gif
        - 慣性空間に対して一定の姿勢レートになりつつあることが分かる.
![gif-2](https://drive.google.com/uc?export=view&id=180jlbVf7dVoGE_XQZdS2ufRMSTlUz6M4)

## 5.2. フィードバックしているパラメータの比較のアプローチを用いた考察
- 磁気差分フィードバックB-dot
    - その名の通り「磁気センサの時間差分」を小さくするようなフィードバック則
    - 「磁場に対する角速度」を抑制する
        - 慣性空間に対する角速度を抑制するわけではない
        - 磁場が時間によってそのベクトルの向きが変われば姿勢もそのように変わる
        - 角速度センサの値が0になることはない
- 角速度フィードバックB-dot
    - 角速度を小さくするようなフィードバック則
    - 「慣性空間に対する機体座標系の角速度」を抑制することができる

## 5.3. 動座標系における微分に基づく結果の考察
### 5.3.1.動座標系における微分
動座標系における微分について説明していく
- ベクトルに関する定義
動座標における微分を理解するにあたりベクトルをしっかり分類する必要がある
    - 幾何ベクトル$\vec{A}$
        - 物理的空間におけるベクトルで,座標系に依存しない
        - 物理量としてのベクトルの本質(感想)
            - 速度ベクトルや位置ベクトルそのものを指す
            - 磁気モーメントベクトルや磁場ベクトルそのものはこれを指す
    - 成分ベクトル$A_D$
        - 動座標系$D$(回転している座標系)で表現されたベクトル
        - 座標系に依存して成分が変化する
        - 線形代数・数値計算のベクトルとしての本質(感想)
            - センサで計測した加速度や磁場は成分できる
                - センサが計測しているのは「センサ座標系」の成分ベクトル
- 動座標系とは？
    - 基準座標系(慣性座標系)に対して回転や加速度運動を行っている座標系
    - ここでは回転座標系を指すものとする
- 動座標系内での幾何ベクトルの微分
    - 動座標系において,幾何ベクトル$\vec{A}$の微分は以下の2つの要素に分けて考える.
        - 動座標系内での物理量の変化(動座標系自体の変化の効果を考慮しない)
        - 動座標系の回転による見かけの変化(回転による影響)
    - まとめると以下の式で表すことができる
```math
(\dfrac{{d}\boldsymbol{A}}{dt})_D = (\dfrac{{\delta}\boldsymbol{A_D}}{\delta{}t})_D + \boldsymbol{{\omega_D}} \times \boldsymbol{A_D}
```
- 式の各項について
    - $(\dfrac{{d}\boldsymbol{A}}{dt})_D $
        - 幾何ベクトルの時間微分のベクトルを成分ベクトルで表現したもの
    - $(\dfrac{{\delta}\boldsymbol{A_D}}{\delta{}t})_D$
        - 動座標系上においてのベクトルの変化を考慮したもの
        - 動座標における成分ベクトルの微分
        - 成分そのものの変化量を表す
    - $\boldsymbol{{\omega_D}} \times \boldsymbol{A_D}$
        - 動座標系が回転しているとき,その回転による変化を表す
        - $\omega_D$とは,基準座標系と動座標系の角速度$\vec{\omega}$を動座標系で見た成分
        - 座標系の回転によって生じる変化量

この式は,動座標系内での物理量の変化に回転の影響を加えたもので,回転している座標系の中で物理量がどのように変化するのかを説明する式である.

- 回転による影響
    - 動座標系が回転するため,時間変化する要素が2個になることが重要
        - 幾何ベクトルは実際には変化していないが,回っている座標系から見たら動いているように見える(見かけ上の変化)がある
- まとめ
    - 動座標系における一般的な幾何ベクトルの微分は次の2つの要素を含む
        - 動座標系から見た幾何ベクトルの変化
        - 動座標が回転することによる見かけ上の変化
    - これらの効果から,回転する座標系内で物理量がどのように変化するか明確に記述可能です！
    - この微分の式は,位置・速度・加速度・角速度といった物理量の変化を扱う際に重要
### 5.3.2.それぞれのB-dotにおける「微分」の取り扱い
- **磁気差分フィードバックB-dot**
    - こちらに関して,考慮しているのは動座標系の微分における右辺第1項
    - 機体座標系における磁場の変化が小さくなるようにフィードバックしている
- **角速度フィードバックB-dot**
    - こちらに関して,考慮しているのは動座標系の微分における右辺第2項
    - 角速度と逆方向のトルクを出力するための計算をしている
    - これはどちらかというとCross-Product則に似ている(今後詳しく解説書きます.)

# 6.それぞれのB-dotの利点・欠点
# 6.1.磁気差分フィードバック型B-dot
- 利点
    - 磁気センサさえあれば良いので実装がシンプルである
    - 座標系が変化しない,一般の微分の定義から近似的に導出ができて直観的
- 欠点
    - サンプリングレートが重要
        - 短すぎるとノイズを拾ったり,変化がなさすぎたり
        - 長すぎるとこれまた正確に変化を正確にとらえられない
# 6.2.角速度フィードバック型B-dot
- 利点
    - 慣性空間に対する角速度をゼロに近づけることができる
- 欠点
    - 角速度を0に近づけるだけなので「どの方向を向いているか」わからない
    - センサデータとしては「角速度」と「磁場」が必要

# 7.さいごに
- 筆者が所属する団体とは関係なく書いているものです.
- 内容については自分が知りえる範囲で書いているものなので,間違い等あればご指摘いただきたいです.
- 不明点あればお気軽にコメントください！