地球観測衛星のカメラを地表の特定地点に向けたい、通信衛星のアンテナを地上局に正確に指向させたい — これらのミッションを実現するには、衛星の3つの回転軸すべてを独立に制御する必要があります。
スピン安定化は1軸の回転安定性しか提供しません。モーメンタムバイアス方式も、ジャイロ効果による受動安定であり、精密なポインティングや大角度マヌーバには限界があります。そこで登場するのが三軸姿勢制御（three-axis attitude control）です。
三軸制御はフィードバック制御理論を直接的に応用するアプローチであり、姿勢制御の中核技術です。三軸制御を理解すると、以下の能力が身につきます。
- 姿勢制御系の設計: 制御ゲインの選定からシミュレーション検証まで一貫した設計プロセスを実行できます
- 衛星の性能評価: 姿勢精度、マヌーバ時間、外乱抑制能力を定量的に評価できます
- 実ミッションへの応用: 地球指向、太陽指向、慣性指向など多様なポインティングモードの設計に応用できます
本記事の内容
- 三軸安定化の基本概念
- 姿勢ダイナミクスの線形化
- PD制御則の設計と物理的意味
- 固有値配置による安定性解析
- 姿勢マヌーバの制御
- 外乱トルク下での定常偏差
- Pythonによる姿勢マヌーバシミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
三軸安定化の基本概念
三軸制御とは
三軸姿勢制御（three-axis stabilization）とは、衛星のロール・ピッチ・ヨーの3つの回転自由度をすべて能動的に制御し、所望の姿勢に維持する方式です。
家庭のエアコンのサーモスタットが温度を一定に保つように、三軸制御は「目標姿勢からのずれ」をセンサで検出し、アクチュエータでトルクを発生させて姿勢を修正します。ただし、温度制御が1自由度であるのに対し、姿勢制御は3自由度を同時に扱う必要があるため、複雑さが格段に増します。
制御ループの構成要素
三軸姿勢制御系は、以下の3つの要素から構成されます。
- 姿勢センサ: 現在の姿勢を計測する（スタートラッカ、太陽センサ、ジャイロスコープなど）
- 制御コンピュータ: 姿勢偏差から必要なトルクを計算する（制御則の実装）
- アクチュエータ: 計算されたトルクを物理的に発生させる（リアクションホイール、スラスタなど）
本記事では、制御則としてPD（比例微分）制御を取り上げます。PD制御は構造が単純で直感的に理解でき、実際の衛星姿勢制御でも広く使われている基本的な制御手法です。
まずは、制御則を設計するための「制御対象のモデル」、すなわち姿勢ダイナミクスの数学的表現を準備しましょう。
姿勢ダイナミクスの線形化
姿勢表現の選択
小角度の姿勢偏差を扱う場合、オイラー角（ロール $\phi$、ピッチ $\theta$、ヨー $\psi$）が便利です。目標姿勢からの偏差が小さいとき、オイラー角は直感的にわかりやすく、線形化が容易です。
姿勢角ベクトルを $\bm{\Theta} = (\phi, \theta, \psi)^T$ と定義します。
運動学方程式（キネマティクス）
小角度近似（$|\phi|, |\theta|, |\psi| \ll 1$）のもとで、姿勢角と角速度の関係は次のように簡略化されます。
$$ \dot{\bm{\Theta}} \approx \bm{\omega} $$
すなわち、角速度がそのまま姿勢角の時間変化率になります。これは大角度では成り立たない近似ですが、制御系の安定性解析やゲイン設計には十分な精度を持ちます。
ダイナミクス方程式
オイラーの運動方程式を小角度近似で線形化します。衛星の角速度が小さく（$|\omega_i| \ll 1$）、ジャイロスコピック結合項（$\bm{\omega} \times \bm{I}\bm{\omega}$）が制御トルクに比べて十分小さいと仮定すると、
$$ \bm{I}\dot{\bm{\omega}} \approx \bm{N}_{\text{ctrl}} + \bm{N}_{\text{dist}} $$
ここで $\bm{N}_{\text{ctrl}}$ は制御トルク、$\bm{N}_{\text{dist}}$ は外乱トルクです。
さらに、主軸座標系を選んで慣性テンソルを対角化すると、各軸が独立した方程式になります。
$$ I_x \ddot{\phi} = N_{\text{ctrl},x} + N_{\text{dist},x} $$
$$ I_y \ddot{\theta} = N_{\text{ctrl},y} + N_{\text{dist},y} $$
$$ I_z \ddot{\psi} = N_{\text{ctrl},z} + N_{\text{dist},z} $$
ダブルインテグレータモデル
各軸の方程式は同じ構造を持ちます。1軸分を取り出すと、
$$ I \ddot{q} = N_{\text{ctrl}} + N_{\text{dist}} $$
ここで $q$ は姿勢角、$I$ はその軸の慣性モーメントです。
これは「慣性モーメント $I$ を持つダブルインテグレータ」と呼ばれる、制御工学で最も基本的なプラントです。入力（トルク）を2回積分すると出力（角度）が得られます。
ダブルインテグレータは「トルクをかけなければ姿勢角は等速度で変化し続ける」という物理を表しています。宇宙空間に摩擦はないため、一度回転し始めた衛星はトルクをかけなければ回転し続けます。これが姿勢制御を必要とする根本的な理由です。
線形化されたダイナミクスモデルが得られたので、次にこのモデルに対する制御則を設計します。
PD制御則の設計
制御則の定義
PD制御は、目標姿勢からの偏差（比例項 P）とその変化率（微分項 D）に基づいてトルクを生成する制御手法です。
1軸の制御トルクは次のように定義されます。
$$ N_{\text{ctrl}} = -K_p (q – q_{\text{ref}}) – K_d \dot{q} $$
ここで $q_{\text{ref}}$ は目標姿勢角、$K_p > 0$ は比例ゲイン、$K_d > 0$ は微分ゲインです。
3軸に拡張すると、
$$ \begin{equation} \bm{N}_{\text{ctrl}} = -\bm{K}_p (\bm{\Theta} – \bm{\Theta}_{\text{ref}}) – \bm{K}_d \dot{\bm{\Theta}} \end{equation} $$
ここで $\bm{K}_p = \text{diag}(K_{p,x}, K_{p,y}, K_{p,z})$、$\bm{K}_d = \text{diag}(K_{d,x}, K_{d,y}, K_{d,z})$ はゲイン行列です。
PD制御の物理的直感
PD制御の各項を直感的に理解しましょう。
比例項 $-K_p(q – q_{\text{ref}})$: バネの復元力に相当します。姿勢が目標からずれると、目標に戻す方向のトルクが発生します。ずれが大きいほどトルクも大きくなります。比例項だけでは振動が生じます。
微分項 $-K_d \dot{q}$: ダンパー（粘性摩擦）に相当します。角速度に比例した制動トルクが発生し、振動を減衰させます。宇宙空間には摩擦がないため、この「人工的な摩擦」が振動抑制に不可欠です。
つまり、PD制御は「バネ＋ダンパー」で衛星を目標姿勢に引き戻す仕組みです。これは古典力学の減衰調和振動子と全く同じ構造です。
閉ループ系の方程式
PD制御則をダイナミクスに代入すると、外乱がない場合（$N_{\text{dist}} = 0$）、
$$ I\ddot{q} + K_d \dot{q} + K_p q = K_p q_{\text{ref}} $$
偏差 $e = q – q_{\text{ref}}$ について書き直すと、
$$ I\ddot{e} + K_d \dot{e} + K_p e = 0 $$
これは減衰調和振動子（2次系）の方程式です。
特性方程式は
$$ Is^2 + K_d s + K_p = 0 $$
その根は
$$ s = \frac{-K_d \pm \sqrt{K_d^2 – 4I K_p}}{2I} $$
この閉ループ系の安定性と過渡特性を、制御工学の標準的なパラメータで解析しましょう。
安定性解析と制御パラメータの設計
標準形への変換
2次系の特性方程式を標準形で表現すると、
$$ s^2 + 2\zeta\omega_n s + \omega_n^2 = 0 $$
ここで $\omega_n$ は固有角振動数（natural frequency）、$\zeta$ は減衰比（damping ratio）です。
ダイナミクスの方程式と比較すると、
$$ \omega_n = \sqrt{\frac{K_p}{I}}, \quad \zeta = \frac{K_d}{2\sqrt{I K_p}} $$
逆に、所望の $\omega_n, \zeta$ からゲインを設計できます。
$$ K_p = I \omega_n^2, \quad K_d = 2\zeta \omega_n I $$
安定性条件
閉ループ系が安定であるための必要十分条件は、特性方程式の全ての根の実部が負であることです。2次系では、
$$ K_p > 0, \quad K_d > 0 $$
で安定です。物理的には、復元力（$K_p$）と減衰力（$K_d$）の両方が正であれば、系は目標に収束します。
減衰比の意味
$\zeta$ の値によって過渡応答の性質が決まります。
| 減衰比 | 応答の特徴 |
|---|---|
| $\zeta = 0$ | 持続振動（減衰なし） |
| $0 < \zeta < 1$ | 不足減衰（振動しながら収束） |
| $\zeta = 1$ | 臨界減衰（最速の非振動収束） |
| $\zeta > 1$ | 過減衰（遅い非振動収束） |
衛星の姿勢制御では、$\zeta = 0.5$ ～ $0.8$ 程度の不足減衰が一般的です。オーバーシュートが許容範囲内で、かつ収束が早いバランスを取ります。
応答特性の指標
設計で重要な時間応答の指標は以下の通りです。
整定時間（settling time）: 目標値の $\pm 2\%$ 以内に収まるまでの時間。
$$ t_s \approx \frac{4}{\zeta \omega_n} $$
オーバーシュート: ステップ応答のピーク値の超過割合。
$$ M_p = e^{-\pi\zeta/\sqrt{1-\zeta^2}} \times 100\% \quad (0 < \zeta < 1) $$
立ち上がり時間（rise time）: 目標値の10%から90%に達するまでの時間。
$$ t_r \approx \frac{1.8}{\omega_n} $$
たとえば、$\omega_n = 0.1$ rad/s、$\zeta = 0.7$ を選ぶと:
- 整定時間: $t_s \approx 57$ s
- オーバーシュート: $M_p \approx 4.6\%$
- 立ち上がり時間: $t_r \approx 18$ s
これらの指標はミッション要求（姿勢マヌーバ時間、ポインティング精度など）から逆算して $\omega_n, \zeta$ を決定し、そこからゲイン $K_p, K_d$ を求めるのが標準的な設計フローです。
設計理論が整ったところで、外乱トルクが存在する場合の定常偏差についても確認しておきましょう。
外乱トルク下での定常偏差
定常外乱トルクの影響
一定の外乱トルク $N_{\text{dist}}$ が作用する場合、定常状態（$\dot{e} = 0, \ddot{e} = 0$）では
$$ K_p e_{ss} = N_{\text{dist}} $$
したがって、
$$ e_{ss} = \frac{N_{\text{dist}}}{K_p} $$
PD制御では、定常外乱に対してゼロでない定常偏差が生じます。偏差を小さくするには $K_p$ を大きくする必要がありますが、ゲインを上げすぎるとセンサノイズの増幅や振動的な応答を引き起こす可能性があります。
PID制御による定常偏差の除去
定常偏差をゼロにするためには、積分項（I）を追加してPID制御にします。
$$ N_{\text{ctrl}} = -K_p e – K_i \int_0^t e \, d\tau – K_d \dot{e} $$
積分項は偏差の累積を元にトルクを生成するため、定常偏差があれば時間とともにトルクが蓄積し、最終的に偏差をゼロに追い込みます。
ただし、積分項の追加は系の次数を上げ、安定余裕を減らす可能性があるため、積分ゲイン $K_i$ は慎重に選ぶ必要があります。実際の衛星では、PD制御を基本としつつ、必要に応じて小さな積分項を加えるのが一般的です。
定量的な見積もり
典型的なLEO衛星のパラメータで定常偏差を見積もってみましょう。
- 外乱トルク: $N_{\text{dist}} = 5 \times 10^{-5}$ Nm（重力傾度トルク）
- 慣性モーメント: $I = 20$ kg$\cdot$m$^2$
- 所望の固有振動数: $\omega_n = 0.1$ rad/s
$K_p = I\omega_n^2 = 20 \times 0.01 = 0.2$ Nm/rad より、
$$ e_{ss} = \frac{5 \times 10^{-5}}{0.2} = 2.5 \times 10^{-4} \text{ rad} \approx 0.014^\circ $$
多くのミッションでこの精度は十分ですが、サブアーク秒の精度を要求する天文観測衛星では積分項が必要になります。
理論的な解析を終えたところで、Pythonでシミュレーションを行い、PD制御の振る舞いを可視化しましょう。
Pythonによる三軸姿勢制御シミュレーション
基本設定とPD制御の実装
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# 衛星パラメータ
Ix, Iy, Iz = 20.0, 25.0, 15.0 # 主慣性モーメント [kg m^2]
I_diag = np.array([Ix, Iy, Iz])
# 制御パラメータ
omega_n = np.array([0.1, 0.08, 0.12]) # 各軸の固有角振動数 [rad/s]
zeta = np.array([0.7, 0.7, 0.7]) # 各軸の減衰比
# ゲイン計算
Kp = I_diag * omega_n**2
Kd = 2 * zeta * omega_n * I_diag
print("制御ゲイン:")
for i, axis in enumerate(['Roll (X)', 'Pitch (Y)', 'Yaw (Z)']):
print(f" {axis}: Kp = {Kp[i]:.4f} Nm/rad, Kd = {Kd[i]:.4f} Nms/rad")
ts = 4 / (zeta[i] * omega_n[i])
Mp = np.exp(-np.pi * zeta[i] / np.sqrt(1 - zeta[i]**2)) * 100
print(f" 整定時間: {ts:.1f} s, オーバーシュート: {Mp:.1f}%")
def attitude_dynamics(t, state, Kp, Kd, theta_ref, N_dist):
"""線形化された三軸姿勢ダイナミクス with PD制御"""
theta = state[:3] # 姿勢角 [rad]
omega = state[3:6] # 角速度 [rad/s]
# 姿勢偏差
error = theta - theta_ref
# PD制御トルク
N_ctrl = -Kp * error - Kd * omega
# 角加速度
alpha = (N_ctrl + N_dist) / I_diag
return np.concatenate([omega, alpha])
ケース1: ステップ応答（姿勢マヌーバ）
目標姿勢を各軸10度にステップ変更し、PD制御の応答を確認します。
# 目標姿勢: 各軸10度のマヌーバ
theta_ref = np.radians(np.array([10.0, 10.0, 10.0]))
N_dist_zero = np.array([0.0, 0.0, 0.0])
# 初期状態: 原点姿勢、静止
state0 = np.zeros(6)
t_span = (0, 200)
t_eval = np.linspace(0, 200, 2000)
sol_step = solve_ivp(
attitude_dynamics, t_span, state0, t_eval=t_eval,
args=(Kp, Kd, theta_ref, N_dist_zero),
rtol=1e-10, atol=1e-12
)
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
labels = ['Roll', 'Pitch', 'Yaw']
colors = ['#e74c3c', '#2ecc71', '#3498db']
# 姿勢角
for i, (label, color) in enumerate(zip(labels, colors)):
axes[0].plot(sol_step.t, np.degrees(sol_step.y[i]),
label=label, color=color, linewidth=1.5)
axes[0].axhline(y=np.degrees(theta_ref[i]), color=color,
linestyle='--', alpha=0.5)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Attitude angle [deg]')
axes[0].set_title('Step response: 10-degree maneuver on all axes')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 角速度
for i, (label, color) in enumerate(zip(labels, colors)):
axes[1].plot(sol_step.t, np.degrees(sol_step.y[3+i]),
label=label, color=color, linewidth=1.5)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Angular velocity [deg/s]')
axes[1].set_title('Angular velocity during maneuver')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('three_axis_step_response.png', dpi=150, bbox_inches='tight')
plt.show()
ステップ応答のシミュレーション結果から、PD制御の特性が明確に確認できます。
- 各軸が独立に応答: 線形化モデルでは3軸が独立なため、各軸は個別に設計された固有振動数に従って応答しています。ロール（$\omega_n = 0.1$）、ピッチ（$\omega_n = 0.08$）、ヨー（$\omega_n = 0.12$）の順に応答速度が異なります。
- オーバーシュートの存在: $\zeta = 0.7$ の不足減衰系であるため、各軸とも目標値を若干超えてから収束します。オーバーシュートは理論値の約4.6%と一致しています。
- 整定時間: ヨー軸（$\omega_n$ 最大）が最も速く整定し、ピッチ軸（$\omega_n$ 最小）が最も遅くなっています。これはミッション要求に応じて各軸の応答速度を個別設計できることを示しています。
- 角速度のピーク: マヌーバ開始直後に角速度が立ち上がり、目標付近で減速しています。角速度のピーク値はアクチュエータの能力（最大トルク）で制限されるため、実設計ではトルク飽和も考慮する必要があります。
ケース2: 減衰比の効果
減衰比 $\zeta$ を変えた場合の応答を比較します。
zeta_values = [0.3, 0.5, 0.7, 1.0, 1.5]
fig, ax = plt.subplots(figsize=(12, 5))
for z in zeta_values:
Kd_test = np.array([2*z*omega_n[0]*Ix, 2*z*omega_n[1]*Iy, 2*z*omega_n[2]*Iz])
sol = solve_ivp(
attitude_dynamics, t_span, state0, t_eval=t_eval,
args=(Kp, Kd_test, theta_ref, N_dist_zero),
rtol=1e-10, atol=1e-12
)
# ロール軸のみ表示
ax.plot(sol.t, np.degrees(sol.y[0]),
label=f'$\\zeta$ = {z}', linewidth=1.5)
ax.axhline(y=10, color='gray', linestyle='--', alpha=0.5, label='Target')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Roll angle [deg]')
ax.set_title('Effect of damping ratio on roll axis response')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('damping_ratio_effect.png', dpi=150, bbox_inches='tight')
plt.show()
減衰比の効果が鮮明に現れています。
- $\zeta = 0.3$（不足減衰・弱）: 大きなオーバーシュートと振動的な応答が見られます。整定に時間がかかり、ポインティング精度が悪くなります。
- $\zeta = 0.5$: 適度な振動を伴いながら収束します。
- $\zeta = 0.7$: 小さなオーバーシュートで速やかに収束する、バランスの良い応答です。多くの衛星で採用される値です。
- $\zeta = 1.0$（臨界減衰）: 振動なしで目標に収束しますが、$\zeta = 0.7$ と比べて整定が少し遅くなります。
- $\zeta = 1.5$（過減衰）: 振動は全くありませんが、応答が非常に遅くなります。
ケース3: 外乱トルク下での応答
一定の外乱トルクが作用する場合の定常偏差を確認します。
# 外乱トルク
N_dist = np.array([5e-5, 1e-4, 3e-5]) # [Nm]
# PD制御で外乱に対する応答
theta_ref_zero = np.array([0.0, 0.0, 0.0]) # 目標姿勢: ゼロ
sol_dist = solve_ivp(
attitude_dynamics, (0, 500), np.zeros(6), t_eval=np.linspace(0, 500, 5000),
args=(Kp, Kd, theta_ref_zero, N_dist),
rtol=1e-10, atol=1e-12
)
# 理論的定常偏差
e_ss_theory = N_dist / Kp
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
for i, (label, color) in enumerate(zip(labels, colors)):
axes[0].plot(sol_dist.t, np.degrees(sol_dist.y[i]) * 3600, # arcsec
label=f'{label} (sim)', color=color, linewidth=1.5)
axes[0].axhline(y=np.degrees(e_ss_theory[i]) * 3600, color=color,
linestyle='--', alpha=0.7,
label=f'{label} theory = {np.degrees(e_ss_theory[i])*3600:.1f}"')
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Attitude error [arcsec]')
axes[0].set_title('Disturbance rejection with PD control')
axes[0].legend(loc='right', fontsize=9)
axes[0].grid(True, alpha=0.3)
# 制御トルク
for i, (label, color) in enumerate(zip(labels, colors)):
N_ctrl_i = -Kp[i] * sol_dist.y[i] - Kd[i] * sol_dist.y[3+i]
axes[1].plot(sol_dist.t, N_ctrl_i * 1e3,
label=label, color=color, linewidth=1.0, alpha=0.8)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Control torque [mNm]')
axes[1].set_title('Control torque for disturbance rejection')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('disturbance_rejection.png', dpi=150, bbox_inches='tight')
plt.show()
print("定常偏差 (理論値):")
for i, axis in enumerate(['Roll', 'Pitch', 'Yaw']):
e_deg = np.degrees(e_ss_theory[i])
e_arcsec = e_deg * 3600
print(f" {axis}: {e_arcsec:.1f} arcsec ({e_deg:.4f} deg)")
外乱応答の結果から、以下のことがわかります。
- 定常偏差の存在: PD制御では一定外乱に対して定常偏差が残ります。シミュレーションの収束値は理論値（破線）と完全に一致しています。これは解析式 $e_{ss} = N_{\text{dist}}/K_p$ の正しさを確認するものです。
- 偏差の大きさ: 各軸で数十アーク秒の定常偏差が生じています。地球観測衛星で要求される精度（一般に数百アーク秒程度）であればPD制御で十分ですが、天文観測衛星（サブアーク秒精度）では積分項の追加が必要です。
- 制御トルクの定常値: 定常状態で制御トルクが外乱トルクと釣り合っています。この定常トルクがリアクションホイールの角運動量蓄積の原因であり、定期的なアンローディングが必要になります。
ケース4: 非線形ダイナミクスでのシミュレーション
最後に、ジャイロスコピック結合を含む非線形ダイナミクスでPD制御の性能を検証します。大角度マヌーバではジャイロ効果が無視できなくなります。
def attitude_dynamics_nonlinear(t, state, Kp, Kd, theta_ref, N_dist):
"""非線形姿勢ダイナミクス（ジャイロスコピック結合あり）"""
theta = state[:3]
omega = state[3:6]
# PD制御トルク
error = theta - theta_ref
N_ctrl = -Kp * error - Kd * omega
# ジャイロスコピック項を含むオイラー方程式
N_total = N_ctrl + N_dist
alpha = np.zeros(3)
alpha[0] = (N_total[0] - (Iz - Iy) * omega[1] * omega[2]) / Ix
alpha[1] = (N_total[1] - (Ix - Iz) * omega[2] * omega[0]) / Iy
alpha[2] = (N_total[2] - (Iy - Ix) * omega[0] * omega[1]) / Iz
return np.concatenate([omega, alpha])
# 大角度マヌーバ: 各軸30度
theta_ref_large = np.radians(np.array([30.0, 20.0, 25.0]))
sol_nl = solve_ivp(
attitude_dynamics_nonlinear, (0, 300), np.zeros(6),
t_eval=np.linspace(0, 300, 3000),
args=(Kp, Kd, theta_ref_large, N_dist_zero),
rtol=1e-10, atol=1e-12
)
sol_lin = solve_ivp(
attitude_dynamics, (0, 300), np.zeros(6),
t_eval=np.linspace(0, 300, 3000),
args=(Kp, Kd, theta_ref_large, N_dist_zero),
rtol=1e-10, atol=1e-12
)
fig, axes = plt.subplots(3, 1, figsize=(12, 10))
for i, (label, color) in enumerate(zip(labels, colors)):
axes[i].plot(sol_nl.t, np.degrees(sol_nl.y[i]),
label='Nonlinear', color=color, linewidth=1.5)
axes[i].plot(sol_lin.t, np.degrees(sol_lin.y[i]),
label='Linear', color=color, linewidth=1.0, linestyle='--', alpha=0.7)
axes[i].axhline(y=np.degrees(theta_ref_large[i]), color='gray',
linestyle=':', alpha=0.5)
axes[i].set_ylabel(f'{label} [deg]')
axes[i].legend()
axes[i].grid(True, alpha=0.3)
axes[0].set_title('Comparison: Nonlinear vs Linear dynamics (large-angle maneuver)')
axes[2].set_xlabel('Time [s]')
plt.tight_layout()
plt.savefig('nonlinear_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
非線形ダイナミクスと線形モデルの比較から、重要な知見が得られます。
- 小角度での一致: 整定後は両モデルの結果がほぼ一致しています。小角度近似が有効な領域では線形モデルが十分な精度を持つことが確認できます。
- マヌーバ中の乖離: 30度級の大角度マヌーバでは、ジャイロスコピック結合の影響で非線形モデルの応答が線形モデルからわずかにずれています。特に、一軸のマヌーバが他軸に干渉する現象（クロスカップリング）が見られます。
- PD制御のロバスト性: 非線形効果にもかかわらず、PD制御は最終的に目標姿勢に収束しています。PD制御はモデルの不確かさに対してある程度のロバスト性を持ちます。ただし、60度以上の大角度マヌーバでは、オイラー角のジンバルロック問題やジャイロ結合の強い影響により、クォータニオンベースの制御やフィードフォワード補償が必要になります。
設計ガイドライン
ゲイン設計の手順
- ミッション要求の整理: 整定時間、ポインティング精度、オーバーシュート許容量、外乱トルクの大きさを明確にする
- 減衰比の選定: $\zeta = 0.5$ ～ $0.8$（一般的な推奨値）
- 固有振動数の選定: 整定時間 $t_s \approx 4/(\zeta\omega_n)$ から $\omega_n$ を逆算
- ゲインの計算: $K_p = I\omega_n^2$, $K_d = 2\zeta\omega_n I$
- 定常偏差の確認: $e_{ss} = N_{\text{dist}}/K_p$ がミッション要求内か確認
- 必要に応じてPI項を追加: 定常偏差が許容範囲を超える場合
実設計での注意点
- アクチュエータの飽和: RWの最大トルクを超えないようにマヌーバ速度を制限する
- センサノイズ: 微分ゲイン $K_d$ が大きすぎるとジャイロのノイズが増幅される
- サンプリング効果: デジタル制御ではサンプリング周波数が固有振動数の10倍以上あることを確認
- 柔軟構造: 太陽電池パドルなどの柔軟付属物の振動モードとの干渉を避ける
まとめ
本記事では、人工衛星の三軸姿勢制御をPD制御の枠組みで体系的に解説しました。
- 線形化ダイナミクス: 小角度近似でオイラー方程式を線形化すると、各軸が独立なダブルインテグレータとなり、古典制御理論を直接適用できます
- PD制御の本質: 比例項（バネ）と微分項（ダンパー）による減衰調和振動子であり、$\omega_n$ と $\zeta$ の2パラメータで応答特性が完全に決定されます
- ゲイン設計: ミッション要求から所望の整定時間・オーバーシュートを定め、$K_p = I\omega_n^2, K_d = 2\zeta\omega_n I$ で算出します
- 定常偏差: PD制御では外乱トルクに対して $e_{ss} = N_{\text{dist}}/K_p$ の偏差が残り、必要に応じて積分項を追加します
- 非線形効果: 大角度マヌーバではジャイロスコピック結合が顕著になりますが、PD制御のロバスト性により収束は維持されます
次のステップとして、以下の記事も参考にしてください。