Skip to content

Control Inputs

Control inputs allow you to add external forcing functions to your dynamics without modifying the core dynamics function. This separation is useful for modeling thrust, drag, or other perturbations that can be toggled on and off.

What are Control Inputs?

In control theory, a dynamical system with control inputs is written as:

\[\dot{\mathbf{x}} = \mathbf{f}(t, \mathbf{x}) + \mathbf{u}(t, \mathbf{x})\]

where:

  • \(\mathbf{f}(t, \mathbf{x})\) is the nominal dynamics (e.g., gravitational acceleration)
  • \(\mathbf{u}(t, \mathbf{x})\) is the control input (e.g., thrust acceleration)

The control input function \(\mathbf{u}\) takes the current time and state and returns a vector that is added to the state derivative. This additive structure makes it easy to:

  • Switch control on/off without changing the dynamics function
  • Combine different control strategies
  • Test different control laws with the same dynamics

Control Input Function Signature

The control input function must follow specific signatures depending on the language:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
def control_function(t: float, state: np.ndarray, params: np.ndarray) -> np.ndarray:
    """
    Args:
        t: Current time
        state: Current state vector
        params: Additional parameters

    Returns:
        Control vector of same dimension as state
    """
    pass

Must be a closure or function with the signature, that either uses dynamic or static sized vectors:

1
Fn(f64, DVector<f64>, DVector<f64>) -> DVector<f64>

or

1
Fn(f64, SVector<f64, S>, SVector<f64, P>) -> SVector<f64, S>

The function receives: - Current time as a scalar - Current state vector - Additional parameters as a vector

And returns a control vector of the same dimension as the state, which is added to the derivative computed by the dynamics function. The additional parameters can be ignored if not needed, but the signature must be maintained.

Using Control Inputs

Control inputs are passed as a separate parameter to the integrator constructor. The dynamics function computes the nominal state derivative, and the control function computes the perturbation that is added to it.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def dynamics(t, state):
    """Orbital dynamics (gravity only).

    Args:
        t: Time
        state: [x, y, z, vx, vy, vz]
    """
    r = state[:3]
    v = state[3:]
    r_norm = np.linalg.norm(r)
    a_grav = -bh.GM_EARTH / (r_norm**3) * r
    return np.concatenate([v, a_grav])


def control_input(t, state):
    """Control input: constant low thrust in velocity direction.

    Args:
        t: Time
        state: [x, y, z, vx, vy, vz]

    Returns:
        Control vector added to state derivative
    """
    v = state[3:]
    v_norm = np.linalg.norm(v)

    control = np.zeros(6)
    if v_norm > 0:
        thrust_magnitude = 0.001  # m/s^2
        control[3:] = thrust_magnitude * v / v_norm

    return control


# Initial LEO state (500 km altitude, circular orbit)
sma = bh.R_EARTH + 500e3
state_initial = np.array([sma, 0.0, 0.0, 0.0, 7612.6, 0.0])

# Orbital period
period = bh.orbital_period(sma)

# Create integrator WITH control input parameter
config = bh.IntegratorConfig.adaptive(abs_tol=1e-10, rel_tol=1e-8)
integrator_thrust = bh.DP54Integrator(
    6, dynamics, jacobian=None, control_fn=control_input, config=config
)

# Create integrator without control (for comparison)
integrator_coast = bh.DP54Integrator(6, dynamics, config=config)

# Propagate with thrust for one orbit
state_thrust = state_initial.copy()
t = 0.0
dt = 60.0

while t < period:
    result = integrator_thrust.step(t, state_thrust, dt)
    state_thrust = result.state
    t += result.dt_used
    dt = result.dt_next

# Propagate without thrust for comparison
state_coast = state_initial.copy()
t = 0.0
dt = 60.0

while t < period:
    result = integrator_coast.step(t, state_coast, dt)
    state_coast = result.state
    t += result.dt_used
    dt = result.dt_next

# Results
r_initial = np.linalg.norm(state_initial[:3])
r_thrust = np.linalg.norm(state_thrust[:3])
r_coast = np.linalg.norm(state_coast[:3])

print(f"Initial radius: {r_initial / 1000:.3f} km")
print(f"Orbital period: {period / 3600:.2f} hours")
print("\nAfter one orbit:")
print(
    f"  With thrust: {r_thrust / 1000:.3f} km (delta_r = {(r_thrust - r_initial) / 1000:.3f} km)"
)
print(
    f"  Coast only:  {r_coast / 1000:.3f} km (delta_r = {(r_coast - r_initial) / 1000:.3f} km)"
)
# Initial radius: 6878.136 km
# Orbital period: 1.58 hours

# After one orbit:
#   With thrust: 6889.325 km (delta_r = 11.189 km)
#   Coast only:  6878.136 km (delta_r = 0.000 km)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
use brahe::constants::{GM_EARTH, R_EARTH};
use brahe::integrators::*;
use brahe::orbits::keplerian::orbital_period;
use nalgebra::DVector;

fn main() {
    // Orbital dynamics (gravity only)
    let dynamics =
        move |_t: f64, state: &DVector<f64>, _params: Option<&DVector<f64>>| -> DVector<f64> {
            let r = state.fixed_rows::<3>(0);
            let v = state.fixed_rows::<3>(3);
            let r_norm = r.norm();
            let a_grav = -GM_EARTH / (r_norm * r_norm * r_norm) * r;

            let mut state_dot = DVector::<f64>::zeros(6);
            state_dot
                .fixed_rows_mut::<3>(0)
                .copy_from(&v.clone_owned());
            state_dot.fixed_rows_mut::<3>(3).copy_from(&a_grav);
            state_dot
        };

    // Control input: constant low thrust in velocity direction
    let control_input = move |_t: f64, state: &DVector<f64>, _params: Option<&DVector<f64>>| -> DVector<f64> {
        let v = state.fixed_rows::<3>(3);
        let v_norm = v.norm();

        let mut control = DVector::<f64>::zeros(6);
        if v_norm > 0.0 {
            let thrust_magnitude = 0.001; // m/s^2
            let a_thrust = thrust_magnitude * v / v_norm;
            control.fixed_rows_mut::<3>(3).copy_from(&a_thrust);
        }
        control
    };

    // Create integrator with control input
    let config = IntegratorConfig::adaptive(1e-10, 1e-8);
    let integrator = DormandPrince54DIntegrator::with_config(
        6,
        Box::new(dynamics),
        None,                          // No Jacobian
        None,                          // No sensitivity provider
        Some(Box::new(control_input)), // Control input
        config,
    );

    // Initial LEO state (500 km altitude)
    let sma = R_EARTH + 500e3;
    let state = DVector::from_vec(vec![
        sma, 0.0, 0.0, // Position
        0.0, 7612.6, 0.0, // Velocity (circular orbit)
    ]);

    // Integrate for one orbit period
    let period = orbital_period(sma);
    let mut t = 0.0;
    let mut dt = 60.0;
    let mut current_state = state.clone();

    println!("Initial semi-major axis: {:.3} km", sma / 1000.0);
    println!(
        "Integrating with continuous thrust for {:.2} hours...",
        period / 3600.0
    );

    while t < period {
        let result = integrator.step(t, current_state, None, Some(dt));
        current_state = result.state;
        t += result.dt_used;
        dt = result.dt_next;
    }

    let final_r = current_state.fixed_rows::<3>(0).norm();
    let delta_r = final_r - sma;
    println!("Final radius: {:.3} km", final_r / 1000.0);
    println!("Orbit raised by: {:.3} km", delta_r / 1000.0);
}

// Initial radius: 6878.136 km
// Orbital period: 1.58 hours

// After one orbit:
//   With thrust: 6889.325 km (delta_r = 11.189 km)
//   Coast only:  6878.136 km (delta_r = 0.000 km)

Applications

Control inputs are particularly useful for:

  • Orbit raising/lowering: Frequent thrusting to get to desired orbit
  • Station keeping: Small corrections to maintain orbit and compensate for drag
  • Redezvous and Proximity Operations: Relative motion control between satellites
  • Spacecraft Collision Avoidance: Maneuvering to avoid debris

See Also