[Python] Satellite Orbit Stabilization: Minimum Energy Control

௹ 시뮬레이션 목적

최적화 공부하며 진행한 프로젝트로 간단하게 실험 환경을 구성하였다.
지구 궤도 공전 인공위성 제어 시뮬레이션으로 이체 문제(Two-body problem)에서 최소 에너지(Minimum energy)로 제어하도록 설계하였다.

1. 위성 궤도 유지 제어
   위성을 지정된 궤도(500 km 고도)에서 유지하기 위해 설계된 제어 알고리즘 설계한다.
2. 제어 입력 최소화(Minimum energy)
   에너지 효율적인 방식으로 위성을 안정적으로 제어하도록 설계한다.
3. 외란 환경에서 안정성 검증
   외란(난수 외란 및 주기 외란)이 추가된 상황에서도 제어 알고리즘이 궤도 유지하도록 설계한다.

௹ 위성의 운동 방정식 (Newton's Second Law)

위성의 운동은 뉴턴 제2법칙기반으로 개발하였다.
$$
\vec{F} = m \cdot \vec{a}
$$

여기서:

• $\vec{F}$는 위성에 작용하는 모든 힘의 합(중력, 제어 입력, 외란)이다.
• $m$은 위성의 질량이다.
• $\vec{a}$는 위성의 가속도이다.

이를 위성의 위치와 속도로 표현하면 다음과 같다:
\[ 
\vec{a} = \frac{\vec{F}}{m} = \vec{a}_{\text{gravity}} + \vec{a}_{\text{control}} + \vec{a}_{\text{disturbance}} 
\]

위성의 가속도는 중력, 제어 입력, 외란에 의해 결정된다.

def orbital_dynamics(t, X):
    x, y, vx, vy = X
    r = np.sqrt(x**2 + y**2)

    # Gravitational acceleration
    ax_gravity = -G * M * x / r**3
    ay_gravity = -G * M * y / r**3

    # Control input
    ux, uy = control_input(t, X)

    # Disturbance (random and periodic)
    random_disturbance = 100e-4 * np.random.randn(2)
    periodic_disturbance = 100e-4 * np.array([
        np.sin(2 * np.pi * t / 1000), 
        np.cos(2 * np.pi * t / 1000)
    ])

    disturbance = random_disturbance + periodic_disturbance

    # Total acceleration
    ax = ax_gravity + ux + disturbance[0]
    ay = ay_gravity + uy + disturbance[1]

    return [vx, vy, ax, ay]

௹ 중력 가속도 (Gravitational Acceleration)

위성 중심으로 작용하는 중력은 아래 식으로 계산된다.
$$
\vec{a}_{\text{gravity}} = -\frac{G M}{r^3} \vec{r}
$$

필자는 변수 값을 다음과 같이 설정하였다:

• $G$: 중력 상수 ($6.67430 \times 10^{-11} , \text{m}^3 \text{kg}^{-1} \text{s}^{-2}$)이다.
• $M$: 지구 질량 ($5.972 \times 10^{24} , \text{kg}$)이다.
• $r$: 위성 중심에서 지구 중심까지 거리 ($r = \sqrt{x^2 + y^2}$)이다.
• $\vec{r}$: 위성 위치 벡터 ($[x, y]$)이다.

ax_gravity = -G * M * x / r**3
ay_gravity = -G * M * y / r**3

௹ 제어 입력 (Control Input)

제어 입력 목표는 위성을 특정 궤도에 고정하고 속도 유지하는 것이다.
이를 위해 위치 제어와 속도 제어를 사용한다.

(1) 위치 제어

위성이 목표 궤도 반지름($r_{\text{target}}$)에서 벗어날 경우, 위치 오차 기반으로 제어 입력 생성한다.
$$
\vec{a}_{\text{position}} = k_p \cdot (r_{\text{target}} - r) \cdot \frac{\vec{r}}{r}
$$

• $k_p$: 위치 제어 게인이다.
• $r_{\text{target}} - r$: 목표 반지름과 현재 반지름의 차이이다.
• $\frac{\vec{r}}{r}$: 단위 방향 벡터(중심 방향)이다.

ux_position = k_p * (r_target - r) * (x / r)
uy_position = k_p * (r_target - r) * (y / r)

(2) 속도 제어

속도 제어는 현재 속도와 목표 속도 차이를 줄인다. 목표 속도는 원 궤도 유지하기 위해 필요하며 아래 식으로 정의된다.
$$
v_{\text{circular}} = \sqrt{\frac{G M}{r_{\text{target}}}}
$$

원 궤도에서 이상적인 속도는 위와 같다. 이를 바탕으로 속도 제어 입력은 아래와 같이 계산된다.
\[
\vec{a}_{\text{velocity}} = k_v \cdot \left( \vec{v}_{\text{desired}} - \vec{v} \right)
\]

• $k_v$: 속도 제어 게인이다.
• $\vec{v}_{\text{desired}} - \vec{v}$: 목표 속도와 현재 속도 차이이다.

vx_desired = -v_circular * (y / r_target)
vy_desired = v_circular * (x / r_target)

ux_velocity = k_v * (vx_desired - vx)
uy_velocity = k_v * (vy_desired - vy)

(3) 제어 입력의 합

최종 제어 입력은 위치 제어, 속도 제어, 중력 보정을 포함한다.
\[
\vec{a}_{\text{control}} = \vec{a}_{\text{position}} + \vec{a}_{\text{velocity}} - \vec{a}_{\text{gravity}}
\]
여기서 중력 보정은 제어 입력이 중력을 상쇄한다.

ux = ux_position + ux_velocity - ax_gravity
uy = uy_position + uy_velocity - ay_gravity
return np.array([ux, uy])

௹ 외란 (Disturbance)

외란은 위성 제어 시스템에 영향을 미치는 비선형 요인이다. 제어 모델이 정상 작동 확인을 목적으로 구성하였다.
여기서는 필자는 두 가지 외란을 모델링하였다.

1. 난수 외란:
   $$
   \vec{a}_{\text{random}} = c \cdot \mathcal{N}(0, 1)
   $$
   $\mathcal{N}(0, 1)$은 평균이 $0$이고 표준 편차가 $1$인 가우시안 난수(Norm distribution)이다.
2. 주기 외란:
   $$
   \vec{a}_{\text{periodic}} = c \cdot \begin{bmatrix} \sin\left(\frac{2 \pi t}{T}\right) \ \cos\left(\frac{2 \pi t}{T}\right) \end{bmatrix}
   $$
3. 이는 주기적인 진동을 모델링한 것이다.

총 외란은 아래와 같다.
\[
\vec{a}_{\text{disturbance}} = \vec{a}_{\text{random}} + \vec{a}_{\text{periodic}}
\]

random_disturbance = 100e-4 * np.random.randn(2)
periodic_disturbance = 100e-4 * np.array([
    np.sin(2 * np.pi * t / 1000), 
    np.cos(2 * np.pi * t / 1000)
])
disturbance = random_disturbance + periodic_disturbance

௹ 동역학 방정식 (Orbital Dynamics)

합력(合力)을 계산한 동역학 방정식은 다음과 같다.
$$
\begin{bmatrix} \dot{x} \ \dot{y} \ \dot{v_x} \ \dot{v_y} \end{bmatrix} = \begin{bmatrix} v_x \ v_y \ a_x \ a_y \end{bmatrix}
$$

여기서 $a_x$와 $a_y$는 다음과 같이 정의된다.
$$
a_x = a_{\text{gravity},x} + a_{\text{control},x} + a_{\text{disturbance},x}
$$

$$
a_y = a_{\text{gravity},y} + a_{\text{control},y} + a_{\text{disturbance},y}
$$

def orbital_dynamics(t, X):
    x, y, vx, vy = X
    r = np.sqrt(x**2 + y**2)
    ax_gravity = -G * M * x / r**3
    ay_gravity = -G * M * y / r**3
    ux, uy = control_input(t, X)
    random_disturbance = 100e-4 * np.random.randn(2)
    periodic_disturbance = 100e-4 * np.array([
        np.sin(2 * np.pi * t / 1000), 
        np.cos(2 * np.pi * t / 1000)
    ])
    disturbance = random_disturbance + periodic_disturbance
    ax = ax_gravity + ux + disturbance[0]
    ay = ay_gravity + uy + disturbance[1]
    return [vx, vy, ax, ay]

௹ 전체 코드

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

# Constants
G = 6.67430e-11  # Gravitational constant (m^3 kg^-1 s^-2)
M = 5.972e24     # Mass of the Earth (kg)
R = 6.371e6      # Earth's radius (m)
m = 500          # Satellite mass (kg)

# Initial conditions
r0 = R + 500e3   # Initial altitude (500 km above Earth's surface)
v0 = 7.8e3       # Approximate orbital velocity (m/s)
x0, y0 = r0, 0   # Initial position (on the x-axis)
vx0, vy0 = 0, v0 # Initial velocity (perpendicular to position)

# State vector: [x, y, vx, vy]
X0 = [x0, y0, vx0, vy0]

# Define control input function
def control_input(t, X):
    """
    Compute the control input with position and velocity correction.

    Args:
        t (float): Current time.
        X (array): Current state vector [x, y, vx, vy].

    Returns:
        array: Control acceleration [ux, uy].
    """
    x, y, vx, vy = X
    r = np.sqrt(x**2 + y**2)

    # Target orbit parameters
    r_target = R + 500e3  # Target altitude (500 km above surface)
    v_circular = np.sqrt(G * M / r_target)  # Circular velocity at target orbit
    vx_desired = -v_circular * (y / r_target)
    vy_desired = v_circular * (x / r_target)

    # Control gains
    k_p = 0.01  # Position control gain
    k_v = 0.1   # Velocity control gain

    # Position control
    ux_position = k_p * (r_target - r) * (x / r)
    uy_position = k_p * (r_target - r) * (y / r)

    # Velocity control
    ux_velocity = k_v * (vx_desired - vx)
    uy_velocity = k_v * (vy_desired - vy)

    # Gravitational compensation
    ax_gravity = -G * M * x / r**3
    ay_gravity = -G * M * y / r**3

    # Total control input
    ux = ux_position + ux_velocity - ax_gravity
    uy = uy_position + uy_velocity - ay_gravity
    return np.array([ux, uy])

# Define orbital dynamics with control and disturbance
def orbital_dynamics(t, X):
    """
    Orbital dynamics with control input and disturbance.

    Args:
        t (float): Current time.
        X (array): Current state vector [x, y, vx, vy].

    Returns:
        array: Time derivative of the state vector [vx, vy, ax, ay].
    """
    x, y, vx, vy = X
    r = np.sqrt(x**2 + y**2)

    # Gravitational acceleration
    ax_gravity = -G * M * x / r**3
    ay_gravity = -G * M * y / r**3

    # Control input
    ux, uy = control_input(t, X)

    # Disturbance (random and periodic)
    random_disturbance = 100e-4 * np.random.randn(2)  # Small random disturbance
    periodic_disturbance = 100e-4 * np.array([np.sin(2 * np.pi * t / 1000), np.cos(2 * np.pi * t / 1000)])  # Periodic disturbance

    disturbance = random_disturbance + periodic_disturbance

    # Total acceleration
    ax = ax_gravity + ux + disturbance[0]
    ay = ay_gravity + uy + disturbance[1]

    return [vx, vy, ax, ay]

# Solve the dynamics
T = 100000  # Simulation time in seconds
t_eval = np.linspace(0, T, 1000)  # Time points
sol = solve_ivp(orbital_dynamics, [0, T], X0, t_eval=t_eval, rtol=1e-8, atol=1e-8)

# Extract results
x, y = sol.y[0], sol.y[1]
vx, vy = sol.y[2], sol.y[3]
control_inputs = np.array([control_input(t, sol.y[:, i]) for i, t in enumerate(sol.t)])
control_magnitudes = np.linalg.norm(control_inputs, axis=1)  # Compute control magnitude

# Plot results
plt.figure(figsize=(12, 12))

# 1. Orbit plot
plt.subplot(3, 1, 1)
plt.plot(x / 1e3, y / 1e3, label="Orbit with Control and Disturbance", linewidth=2)
plt.plot(0, 0, 'o', label="Earth", color='blue', markersize=10)
plt.xlabel("X Position (km)")
plt.ylabel("Y Position (km)")
plt.title("Orbit with Minimum Energy Control and Disturbance")
plt.legend()
plt.axis('equal')
plt.grid()

# 2. Velocity plot
plt.subplot(3, 1, 2)
plt.plot(sol.t, vx, label="Vx (m/s)", color="orange", linewidth=2)
plt.plot(sol.t, vy, label="Vy (m/s)", color="green", linewidth=2)
plt.xlabel("Time (s)")
plt.ylabel("Velocity (m/s)")
plt.title("Velocity Over Time")
plt.legend()
plt.grid()

# 3. Control input plot
plt.subplot(3, 1, 3)
plt.plot(sol.t, control_magnitudes, label="Control Magnitude (m/s²)", color="red", linewidth=2)
plt.xlabel("Time (s)")
plt.ylabel("Control Magnitude (m/s²)")
plt.title("Control Input Over Time")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()