磁気トルカによる姿勢制御 — 地磁場を利用したトルク生成

CubeSatをロケットから放出した直後、衛星は毎秒数度の速さでぐるぐると回転しています。太陽電池に十分な光が当たらず、通信アンテナの方向も定まらない — ミッション開始前に、この「タンブリング」をまず止めなければなりません。

しかし、CubeSatにはスラスタを搭載する余裕がほとんどありません。1kg以下の衛星に何を載せてこの回転を止めるのでしょうか？

答えは磁気トルカ（Magnetic Torquer, MTQ）です。コイルに電流を流して磁気モーメントを生成し、地球の磁場との相互作用でトルクを発生させます。推進剤も可動部品も不要で、電力だけで動作するこのアクチュエータは、小型衛星の姿勢制御で不可欠な存在です。

磁気トルカを理解すると、以下の実践的な課題に取り組めるようになります。

• デタンブリング制御: 衛星の初期回転を地磁場だけで減衰させるBドット制御を設計できます
• RWアンローディング: リアクションホイールに蓄積した角運動量を推進剤なしで排出する戦略を立てられます
• 小型衛星の姿勢制御系設計: 磁気トルカの能力と制約を踏まえたミッション設計ができます

本記事の内容

• 磁気トルカの動作原理（$\bm{N} = \bm{m} \times \bm{B}$）
• 地磁場モデル（ダイポール近似とIGRF）
• Bドット制御によるデタンブリング
• 磁気トルカの本質的な制約
• RWアンローディングへの応用
• Pythonによるデタンブリングシミュレーションと地磁場の可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• オイラーの運動方程式 — 剛体の回転ダイナミクスを記述する
• 人工衛星の角運動量管理 — ホイールと外部トルクによるバランス

磁気トルカの動作原理

磁気モーメントとトルク

方位磁針が北を向くのは、地球の磁場中で磁気モーメントに回転力（トルク）が作用するためです。磁気トルカはこの原理を能動的に利用する装置です。

磁気双極子モーメント $\bm{m}$ を持つ物体が外部磁場 $\bm{B}$ 中に置かれると、次のトルクが作用します。

$$ \begin{equation} \bm{N} = \bm{m} \times \bm{B} \end{equation} $$

このトルクは磁気モーメント $\bm{m}$ を磁場 $\bm{B}$ の方向に揃えようとする作用を持ちます。磁気トルカでは、コイルに流す電流を制御して $\bm{m}$ の大きさと方向を変化させ、所望のトルクを生成します。

磁気トルカの構造

磁気トルカには主に2種類の構造があります。

エアコアコイル（空芯コイル）:

$N$ 回巻きのコイルに電流 $I$ を流すと、磁気モーメントは

$$ \bm{m} = NIA\hat{\bm{n}} $$

ここで $A$ はコイルの面積、$\hat{\bm{n}}$ はコイル面の法線方向です。

エアコアコイルは軽量で線形性が高いですが、同じ磁気モーメントを得るためにはコイルが大きくなります。CubeSatではPCB（プリント基板）上にパターンとして形成することも多く、構造的にシンプルです。

磁気棒（トルクロッド）:

強磁性体のコアにコイルを巻いた構造です。コアの透磁率が高いため、小型で大きな磁気モーメントを生成できます。

$$ \bm{m} = \mu_r NIA\hat{\bm{n}} $$

ここで $\mu_r$ は比透磁率です。ただし、ヒステリシスや飽和による非線形性が生じます。

磁気トルカの典型的なパラメータ

磁気トルカの最大の利点は、可動部品がなく推進剤も不要であることです。故障率が極めて低く、電力さえあれば何年でも動作し続けます。

しかし、$\bm{N} = \bm{m} \times \bm{B}$ という外積の構造から、本質的な制約も存在します。この制約を理解するためには、まず地磁場 $\bm{B}$ の特性を知る必要があります。

地磁場モデル

ダイポール近似

地球の磁場は、第0近似として磁気双極子で表現できます。地球の中心に磁気モーメント $\bm{M}_E$ を持つ双極子を置いたモデルです。

地磁気ダイポールは地理的な極とはずれており（磁気傾斜角約11度）、ダイポールモーメントの大きさは $M_E \approx 7.94 \times 10^{22}$ Am$^2$ です。

軌道座標系（LVLH）での地磁場のダイポール近似は、衛星の軌道位置（軌道位置角 $\theta$、軌道傾斜角 $i$）に依存します。簡略化した表現では、

$$ \bm{B} = \frac{B_0}{(r/R_E)^3} \begin{pmatrix} -\sin i \sin \theta \\ \cos i \\ -2\sin i \cos \theta \end{pmatrix} $$

ここで $B_0 \approx 3 \times 10^{-5}$ T は赤道上での地表面の磁場強度、$R_E$ は地球半径、$r$ は軌道半径です。

地磁場の軌道に沿った変化

ダイポール近似から重要な特徴が読み取れます。

• 磁場の大きさは軌道に沿って変化する: 衛星が赤道付近にいるときと極域にいるときで磁場の強さが異なります
• 磁場の方向も変化する: 軌道1周で磁場ベクトルは大きく回転します
• 軌道傾斜角の影響: 傾斜角が大きいほど磁場の変動が大きくなります

この変動が磁気トルカの制御に大きな影響を与えます。磁気トルカは地磁場の方向に依存するため、ある瞬間に生成できないトルク方向が存在しますが、軌道に沿って磁場方向が変化するため、軌道1周をかけて平均的に3軸の制御が可能になるのです。

IGRF（国際地磁気標準モデル）

より正確な地磁場モデルとしてIGRF（International Geomagnetic Reference Field）があります。IGRFは球面調和関数展開を用いて地磁場を表現します。

$$ V(r, \theta, \phi) = R_E \sum_{n=1}^{N} \left(\frac{R_E}{r}\right)^{n+1} \sum_{m=0}^{n} (g_n^m \cos m\phi + h_n^m \sin m\phi) P_n^m(\cos\theta) $$

ここで $g_n^m, h_n^m$ はガウス係数（5年ごとに更新）、$P_n^m$ はシュミット準正規化のルジャンドル陪関数です。

実際の衛星設計では IGRF が使用されますが、概念的な理解と初期設計にはダイポール近似で十分です。

地磁場の基本的な特性を理解したところで、磁気トルカの最も重要な応用であるBドット制御に進みましょう。

Bドット制御 — デタンブリング

問題設定

衛星が分離直後にタンブリング（無制御の回転）している状態から、回転を減衰させて安定化する制御をデタンブリング（detumbling）と呼びます。

デタンブリングは衛星運用の最初のステップであり、これが成功しなければ太陽電池の発電も通信も確保できません。デタンブリングに磁気トルカを使う場合、必要な情報は地磁場の時間変化率のみで、姿勢そのものを知る必要がないという大きな利点があります。

Bドット制御則

Bドット（B-dot）制御は、地磁場の時間変化率 $\dot{\bm{B}}$ を利用する、極めてシンプルで効果的な制御則です。

$$ \begin{equation} \bm{m} = -k \dot{\bm{B}}_{\text{body}} \end{equation} $$

ここで $k > 0$ はゲイン、$\dot{\bm{B}}_{\text{body}}$ は機体座標系で計測した地磁場ベクトルの時間変化率です。

Bドット制御がなぜ回転を止めるのか

一見すると、なぜ地磁場の変化率に負の比例で磁気モーメントを生成すると回転が減衰するのか、直感的にはわかりにくいかもしれません。エネルギーの観点から理解しましょう。

衛星が角速度 $\bm{\omega}$ で回転しているとき、機体座標系で見た地磁場の変化率は

$$ \dot{\bm{B}}_{\text{body}} \approx -\bm{\omega} \times \bm{B} $$

です（地磁場自体の時間変化は衛星の回転に比べて遅いため）。

Bドット制御で生成される磁気モーメントは $\bm{m} = -k\dot{\bm{B}}_{\text{body}} \approx k(\bm{\omega} \times \bm{B})$ であり、これによるトルクは

$$ \bm{N} = \bm{m} \times \bm{B} = k(\bm{\omega} \times \bm{B}) \times \bm{B} $$

回転の運動エネルギーの変化率を計算すると、

$$ \dot{T} = \bm{\omega} \cdot \bm{N} = k\bm{\omega} \cdot [(\bm{\omega} \times \bm{B}) \times \bm{B}] $$

ベクトル三重積の公式 $(\bm{a} \times \bm{b}) \times \bm{c} = \bm{b}(\bm{a} \cdot \bm{c}) – \bm{a}(\bm{b} \cdot \bm{c})$ を適用すると、

$$ \dot{T} = k[\bm{\omega} \cdot \bm{B}(\bm{\omega} \cdot \bm{B}) – |\bm{\omega}|^2 |\bm{B}|^2 + (\bm{\omega} \cdot \bm{B})^2] $$

ここで $(\bm{\omega} \cdot \bm{B})^2 \leq |\bm{\omega}|^2|\bm{B}|^2$（コーシー・シュワルツの不等式）であり、等号は $\bm{\omega} \parallel \bm{B}$ のときのみ成り立ちます。

整理すると、

$$ \dot{T} = -k[|\bm{\omega}|^2|\bm{B}|^2 – (\bm{\omega} \cdot \bm{B})^2] = -k|\bm{\omega} \times \bm{B}|^2 \leq 0 $$

運動エネルギーは常に減少する（$\bm{\omega} \parallel \bm{B}$ の瞬間を除く）ことが証明されました。角速度が磁場に平行な成分は減衰しませんが、軌道に沿って磁場方向が変化するため、軌道全体を通じて全方向の角速度が減衰します。

実用上の考慮

$\dot{\bm{B}}$ の計測: 磁力計（磁気センサ）で $\bm{B}$ を計測し、差分で $\dot{\bm{B}}$ を推定します。

$$ \dot{\bm{B}} \approx \frac{\bm{B}(t) – \bm{B}(t – \Delta t)}{\Delta t} $$

ノイズ低減のためにローパスフィルタを併用するのが一般的です。

ゲインの選定: ゲイン $k$ は磁気モーメントの最大値 $m_{\max}$、磁場の大きさ $B$、所望のデタンブリング時間に基づいて選定します。

$$ k \sim \frac{m_{\max}}{|\dot{\bm{B}}_{\max}|} $$

大きすぎるゲインは磁気モーメントの飽和を招き、小さすぎるゲインはデタンブリングに時間がかかります。

Bドット制御の理論を把握したところで、磁気トルカに固有の制約について詳しく検討しましょう。

磁気トルカの本質的な制約

磁場方向にトルク生成不可

$\bm{N} = \bm{m} \times \bm{B}$ という外積の性質上、生成されるトルクは必ず磁場 $\bm{B}$ に垂直です。

$$ \bm{N} \cdot \bm{B} = (\bm{m} \times \bm{B}) \cdot \bm{B} = 0 $$

これは磁気トルカの最も本質的な制約です。任意の瞬間において、地磁場方向のトルクは生成できません。

この制約は、3軸の磁気トルカ（$\bm{m} = (m_x, m_y, m_z)^T$）を持っていても変わりません。$\bm{m}$ をどのように選んでも、$\bm{N} = \bm{m} \times \bm{B}$ は $\bm{B}$ に垂直な平面内にしか存在しません。

制約の軌道スケールでの緩和

瞬間的には2自由度しか制御できませんが、軌道に沿って地磁場の方向が変化するため、軌道1周を通じて見れば3軸全ての制御が可能です。

これを「軌道平均制御可能性（orbit-average controllability）」と呼びます。ただし、以下の条件が必要です。

• 軌道傾斜角が $0°$ でないこと（赤道軌道では磁場方向がほぼ一定のため、制御可能性が制限される）
• 制御帯域が軌道周波数より十分低いこと

磁気トルカだけでは不十分な場合

以下の場合、磁気トルカ単独では要求を満たせないことがあります。

• 高精度ポインティング: 瞬時の3軸独立制御が必要な場合
• 高速マヌーバ: トルクの大きさと方向が制限されるため、速い姿勢変更には不向き
• GEO以遠の軌道: 地磁場が弱すぎるため十分なトルクが得られない
• 赤道軌道: 磁場方向の変化が小さく、特定軸の制御が困難

このため、多くの衛星では磁気トルカをRWと組み合わせて使用します。RWで精密な姿勢制御を行い、磁気トルカは角運動量のアンローディングに特化させるのが効率的です。

RWアンローディングへの利用

アンローディング制御則

磁気トルカの最も重要な応用の一つが、RWに蓄積した角運動量のアンローディング（排出）です。

目的は、RWの角運動量 $\bm{H}_w$ を減少させる外部トルクを磁気トルカで生成することです。

$$ \bm{N}_{\text{MTQ}} = \bm{m} \times \bm{B} $$

が $\bm{H}_w$ を減少させる方向になるよう、$\bm{m}$ を設計します。

一般的な制御則は、

$$ \begin{equation} \bm{m} = -k_u \frac{\bm{H}_w \times \bm{B}}{|\bm{B}|^2} \end{equation} $$

このとき、$\bm{H}_w$ の変化率は

$$ \dot{\bm{H}}_w = \bm{N}_{\text{MTQ}} = \left(-k_u \frac{\bm{H}_w \times \bm{B}}{|\bm{B}|^2}\right) \times \bm{B} $$

ベクトル三重積を展開すると、

$$ \dot{\bm{H}}_w = -k_u \frac{|\bm{B}|^2 \bm{H}_w – (\bm{H}_w \cdot \bm{B})\bm{B}}{|\bm{B}|^2} = -k_u \left[\bm{H}_w – \frac{(\bm{H}_w \cdot \bm{B})\bm{B}}{|\bm{B}|^2}\right] $$

右辺の $[\cdots]$ は $\bm{H}_w$ から $\bm{B}$ 方向の成分を引いたもの、つまり $\bm{H}_w$ の $\bm{B}$ に垂直な成分です。したがって、$\bm{B}$ に垂直な角運動量成分が指数関数的に減衰します。

$\bm{B}$ に平行な成分はこのトルクでは排出できませんが、軌道に沿って $\bm{B}$ の方向が変化するため、複数軌道をかけて全方向の角運動量を排出できます。

アンローディングの設計考慮

• ゲイン $k_u$: 大きすぎると磁気モーメントが飽和し、小さすぎるとアンローディングが遅い
• RW制御との協調: アンローディング中もRWで姿勢を維持するため、アンローディングトルクが姿勢制御を妨げないよう協調設計が必要
• 軌道周期とアンローディング時間: 磁場方向の制約により、完全なアンローディングには数軌道を要する

理論的な解析を終えたところで、Pythonでデタンブリングシミュレーションを実装しましょう。

Pythonによるデタンブリングシミュレーション

地磁場モデルの実装

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 軌道パラメータ
alt = 500e3            # 軌道高度 [m]
R_E = 6371e3           # 地球半径 [m]
r_orb = R_E + alt      # 軌道半径 [m]
mu_E = 3.986e14        # 地球重力定数 [m^3/s^2]
T_orb = 2 * np.pi * np.sqrt(r_orb**3 / mu_E)  # 軌道周期 [s]
omega_orb = 2 * np.pi / T_orb

# 地磁場パラメータ
B0 = 3.12e-5           # 赤道上の地表面磁場 [T]
B_orb = B0 * (R_E / r_orb)**3  # 軌道上の磁場スケール [T]
incl = np.radians(51.6)  # 軌道傾斜角 (ISS相当)

def geomagnetic_field(t):
    """ダイポール近似の地磁場 (LVLH座標系)"""
    theta = omega_orb * t  # 軌道位置角
    Bx = -B_orb * np.sin(incl) * np.sin(theta)
    By = B_orb * np.cos(incl)
    Bz = -2 * B_orb * np.sin(incl) * np.cos(theta)
    return np.array([Bx, By, Bz])

print(f"軌道高度: {alt/1e3:.0f} km")
print(f"軌道周期: {T_orb:.1f} s ({T_orb/60:.1f} min)")
print(f"軌道上の磁場スケール: {B_orb*1e6:.2f} uT")

# 軌道1周の地磁場変化を可視化
t_orbit = np.linspace(0, T_orb, 1000)
B_history = np.array([geomagnetic_field(t) for t in t_orbit])

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 各成分
for i, (label, color) in enumerate(zip(['Bx', 'By', 'Bz'], ['#e74c3c', '#2ecc71', '#3498db'])):
    axes[0].plot(t_orbit / 60, B_history[:, i] * 1e6,
                label=label, color=color, linewidth=1.5)
axes[0].set_xlabel('Time [min]')
axes[0].set_ylabel('Magnetic field [uT]')
axes[0].set_title(f'Geomagnetic field along orbit (alt={alt/1e3:.0f} km, incl={np.degrees(incl):.1f} deg)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 磁場の大きさ
B_mag = np.linalg.norm(B_history, axis=1)
axes[1].plot(t_orbit / 60, B_mag * 1e6, 'k-', linewidth=1.5)
axes[1].set_xlabel('Time [min]')
axes[1].set_ylabel('|B| [uT]')
axes[1].set_title('Magnetic field magnitude')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('geomagnetic_field_orbit.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"磁場の最小値: {np.min(B_mag)*1e6:.2f} uT")
print(f"磁場の最大値: {np.max(B_mag)*1e6:.2f} uT")

軌道1周の地磁場変化から、重要な特徴が読み取れます。

1. 周期的変動: $B_x$ と $B_z$ は軌道周期で正弦波状に変化しています。$B_y$（軌道面法線方向）はほぼ一定であり、ダイポール磁場の対称性を反映しています。
2. 磁場の大きさの変動: 磁場の大きさは軌道1周で約20-30 uTの範囲で変動します。極域付近で強く、赤道付近で弱くなる傾向があります。
3. 磁場方向の回転: 磁場ベクトルの方向が軌道に沿って大きく変化しており、これが磁気トルカによる3軸制御を軌道スケールで可能にする鍵です。

Bドットデタンブリングシミュレーション

# 衛星パラメータ
Ix, Iy, Iz = 0.05, 0.05, 0.03  # CubeSatの慣性モーメント [kg m^2]
I_sat = np.array([Ix, Iy, Iz])
m_max = 0.5  # 最大磁気モーメント [Am^2]

# Bドット制御ゲイン
k_bdot = 1e5  # [Am^2 s / T]

def detumbling_dynamics(t, state):
    """Bドット制御によるデタンブリングダイナミクス"""
    omega = state[:3]  # 角速度 [rad/s]

    # 地磁場（簡易モデル: 衛星回転を無視してLVLH近似）
    B = geomagnetic_field(t)

    # 機体座標系での磁場の時間変化率（衛星回転による）
    B_dot = -np.cross(omega, B)

    # Bドット制御則
    m_cmd = -k_bdot * B_dot

    # 磁気モーメントの飽和
    m_mag = np.linalg.norm(m_cmd)
    if m_mag > m_max:
        m_cmd = m_cmd * m_max / m_mag

    # 磁気トルク
    N_mtq = np.cross(m_cmd, B)

    # オイラー方程式 (CubeSatなのでジャイロ項も含める)
    alpha = np.zeros(3)
    alpha[0] = (N_mtq[0] - (Iz - Iy) * omega[1] * omega[2]) / Ix
    alpha[1] = (N_mtq[1] - (Ix - Iz) * omega[2] * omega[0]) / Iy
    alpha[2] = (N_mtq[2] - (Iy - Ix) * omega[0] * omega[1]) / Iz

    return alpha

# 初期条件: タンブリング状態（各軸5 deg/s）
omega0 = np.radians([5.0, -3.0, 8.0])  # [rad/s]

# 10軌道分のシミュレーション
t_sim = 10 * T_orb
t_span = (0, t_sim)
t_eval = np.linspace(0, t_sim, 50000)

sol = solve_ivp(detumbling_dynamics, t_span, omega0,
                t_eval=t_eval, rtol=1e-10, atol=1e-12)

print(f"初期角速度: {np.degrees(omega0)} deg/s")
print(f"初期角速度の大きさ: {np.degrees(np.linalg.norm(omega0)):.2f} deg/s")

fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# 角速度の各成分
labels = [r'$\omega_x$', r'$\omega_y$', r'$\omega_z$']
colors = ['#e74c3c', '#2ecc71', '#3498db']

for i, (label, color) in enumerate(zip(labels, colors)):
    axes[0].plot(sol.t / T_orb, np.degrees(sol.y[i]),
                label=label, color=color, linewidth=1.0)
axes[0].set_ylabel('Angular velocity [deg/s]')
axes[0].set_title('B-dot detumbling control')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 角速度の大きさ
omega_mag = np.sqrt(sol.y[0]**2 + sol.y[1]**2 + sol.y[2]**2)
axes[1].plot(sol.t / T_orb, np.degrees(omega_mag), 'k-', linewidth=1.5)
axes[1].set_ylabel('|omega| [deg/s]')
axes[1].set_title('Angular velocity magnitude')
axes[1].set_yscale('log')
axes[1].grid(True, alpha=0.3)

# 回転の運動エネルギー
T_rot = 0.5 * (Ix * sol.y[0]**2 + Iy * sol.y[1]**2 + Iz * sol.y[2]**2)
axes[2].plot(sol.t / T_orb, T_rot * 1e3, 'purple', linewidth=1.5)
axes[2].set_xlabel('Orbits')
axes[2].set_ylabel('Rotational KE [mJ]')
axes[2].set_title('Rotational kinetic energy dissipation')
axes[2].set_yscale('log')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('bdot_detumbling.png', dpi=150, bbox_inches='tight')
plt.show()

# デタンブリング性能
omega_final = np.degrees(np.linalg.norm(sol.y[:, -1]))
print(f"最終角速度の大きさ: {omega_final:.4f} deg/s")
print(f"角速度の減衰比: {omega_final / np.degrees(np.linalg.norm(omega0)):.4e}")

デタンブリングシミュレーションの結果から、Bドット制御の効果と特性が明確に確認できます。

1. 指数関数的な減衰: 角速度の大きさが対数スケールでほぼ直線的に減少しており、指数関数的な減衰特性を示しています。これはBドット制御が本質的にエネルギー散逸型の制御であることの反映です。
2. 軌道周期の変調: 角速度の包絡線にはっきりと軌道周期の変調が見られます。これは地磁場の強さと方向が軌道に沿って変化するためです。磁場が弱い赤道付近では減衰が遅く、磁場が強い極域付近では減衰が速くなります。
3. 運動エネルギーの単調減少: 理論で証明した $\dot{T} \leq 0$ がシミュレーションでも確認できます。運動エネルギーは多少のゆらぎを伴いながらも単調に減少しています。
4. デタンブリング時間: 約10軌道（約15時間）で初期角速度の10,000分の1以下まで減衰しており、実用的なデタンブリング性能が得られています。

ゲインの感度解析

Bドット制御のゲインを変えた場合の影響を調べます。

gains = [1e4, 5e4, 1e5, 5e5, 1e6]

fig, ax = plt.subplots(figsize=(12, 5))

for k in gains:
    def dynamics_k(t, state, k_val=k):
        omega = state[:3]
        B = geomagnetic_field(t)
        B_dot = -np.cross(omega, B)
        m_cmd = -k_val * B_dot
        m_mag = np.linalg.norm(m_cmd)
        if m_mag > m_max:
            m_cmd = m_cmd * m_max / m_mag
        N_mtq = np.cross(m_cmd, B)
        alpha = np.zeros(3)
        alpha[0] = (N_mtq[0] - (Iz - Iy) * omega[1] * omega[2]) / Ix
        alpha[1] = (N_mtq[1] - (Ix - Iz) * omega[2] * omega[0]) / Iy
        alpha[2] = (N_mtq[2] - (Iy - Ix) * omega[0] * omega[1]) / Iz
        return alpha

    sol_k = solve_ivp(dynamics_k, (0, 5*T_orb), omega0,
                       t_eval=np.linspace(0, 5*T_orb, 25000),
                       rtol=1e-10, atol=1e-12)
    omega_mag = np.degrees(np.sqrt(sol_k.y[0]**2 + sol_k.y[1]**2 + sol_k.y[2]**2))
    ax.plot(sol_k.t / T_orb, omega_mag, label=f'k = {k:.0e}', linewidth=1.2)

ax.set_xlabel('Orbits')
ax.set_ylabel('|omega| [deg/s]')
ax.set_title('Effect of B-dot gain on detumbling performance')
ax.set_yscale('log')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bdot_gain_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

ゲイン比較の結果から、以下の知見が得られます。

1. ゲインが低い場合（$k = 10^4$）: 減衰が遅く、5軌道経過しても角速度は十分に低下していません。トルク生成能力を十分に活用できていない状態です。
2. 中間のゲイン（$k = 10^5$）: バランスの取れた減衰特性を示し、数軌道で効果的にデタンブリングが完了します。
3. ゲインが高い場合（$k = 10^6$）: 制御の初期段階では磁気モーメントが常に飽和状態（バンバン制御に近い）となります。飽和領域では線形理論が適用できず、減衰特性が線形的（指数関数的ではなく）になります。
4. 最適ゲインの存在: ゲインが大きすぎても小さすぎても最適ではなく、磁気モーメントの飽和を考慮した適切な範囲があります。実用的には、磁気モーメントが飽和しない程度の最大のゲインが最も効率的です。

磁気トルクの制約の可視化

磁気トルカが生成できるトルクの制約を時間変化として可視化します。

# ある瞬間の磁気トルカの制御可能空間を可視化
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

t_samples = [0, T_orb/4, T_orb/2]  # 軌道位置: 0度, 90度, 180度

for ax, t_s in zip(axes, t_samples):
    B = geomagnetic_field(t_s)
    B_hat = B / np.linalg.norm(B)

    # ランダムな磁気モーメントに対するトルクをサンプリング
    n_samples = 3000
    m_samples = np.random.randn(n_samples, 3)
    m_samples = m_samples / np.linalg.norm(m_samples, axis=1, keepdims=True) * m_max

    torques = np.cross(m_samples, B)

    ax.scatter(torques[:, 0] * 1e6, torques[:, 1] * 1e6, s=1, alpha=0.3, c='blue')
    ax.set_xlabel(r'$N_x$ [uNm]')
    ax.set_ylabel(r'$N_y$ [uNm]')
    theta_deg = np.degrees(omega_orb * t_s)
    ax.set_title(f'Torque at orbit angle = {theta_deg:.0f} deg\n'
                 f'B = [{B[0]*1e6:.1f}, {B[1]*1e6:.1f}, {B[2]*1e6:.1f}] uT')
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)

plt.suptitle('Achievable torque set (XY projection) at different orbit positions', fontsize=12)
plt.tight_layout()
plt.savefig('torque_constraint.png', dpi=150, bbox_inches='tight')
plt.show()

制約の可視化から、磁気トルカの本質的な特性が直感的に把握できます。

1. 楕円形のトルク集合: 各軌道位置で生成可能なトルクの集合は、$XY$ 平面への射影が楕円形になっています。磁場 $\bm{B}$ に垂直な平面が制御可能空間であり、$XY$ 平面への射影がこの楕円です。
2. 軌道位置による変化: 軌道位置が変わると楕円の形状と方向が変化しています。磁場方向の回転に伴い、制御可能な方向が回転するのです。
3. 3軸の軌道平均制御可能性: 異なる軌道位置での楕円を重ね合わせると、$XY$ 平面の全方向をカバーできます。同様に $XZ, YZ$ 面でも全方向がカバーされるため、軌道1周を通じて3軸の制御が可能です。

設計ガイドライン

磁気トルカの選定手順

1. 外乱トルクの推定: 軌道環境から期待される外乱トルクの大きさを算出
2. 必要な磁気モーメントの見積もり: 外乱トルク $N_{\text{dist}}$ と磁場 $B$ から $m \geq N_{\text{dist}}/B$
3. デタンブリング要求: 初期角速度と許容デタンブリング時間から必要なゲインと磁気モーメントを逆算
4. アンローディング要求: RWの角運動量蓄積レートと許容蓄積量からアンローディング能力を確認
5. 消費電力の確認: ミッション電力バジェット内で動作可能か確認

磁気トルカを使う場合のチェックリスト

• 軌道傾斜角が $0°$ でないか（赤道軌道では制御可能性が制限される）
• LEO 軌道であるか（GEO以遠では磁場が弱すぎる可能性）
• 瞬時の3軸制御が不要であるか（不要ならMTQ単独で可能、必要ならRWと併用）
• 衛星の残留磁気が管理されているか（意図しない磁気トルクの抑制）
• 磁力計が搭載されているか（Bドット制御やアンローディングに必要）

まとめ

本記事では、磁気トルカによる姿勢制御の原理、制約、応用について解説しました。

• 動作原理: $\bm{N} = \bm{m} \times \bm{B}$ — 地磁場との相互作用でトルクを生成する、可動部品不要のアクチュエータです
• 地磁場モデル: ダイポール近似で基本的な振る舞いが理解でき、軌道に沿って磁場の大きさと方向が周期的に変化します
• Bドット制御: $\bm{m} = -k\dot{\bm{B}}$ というシンプルな制御則でタンブリングを減衰させ、運動エネルギーが常に減少することが理論的に保証されています
• 本質的な制約: 瞬時には磁場に垂直な2方向のみ制御可能ですが、軌道に沿った磁場変化を利用して軌道平均で3軸制御が実現できます
• RWアンローディング: 推進剤を消費せずにRWの蓄積角運動量を排出でき、LEO小型衛星のミッション寿命延長に不可欠です

次のステップとして、以下の記事も参考にしてください。

• リアクションホイールによる姿勢制御 — 角運動量交換の原理と実装
• 人工衛星の角運動量管理 — ホイールと外部トルクによるバランス
• 三軸姿勢制御 — PD制御による衛星のポインティング