GitHub

Integrators

What are Integrators?

Integrators discretize continuous-time dynamics into constraints for the NLP solver. They implement the relationship:

\[x_{k+1} = \Phi(x_k, u_k, \Delta t_k)\]

where Φ approximates the continuous evolution ẋ = f(x, u, t).

using DirectTrajOpt
using NamedTrajectories
using LinearAlgebra

BilinearIntegrator

Overview

Used for control-linear (bilinear) dynamics:

\[\dot{x} = (G_0 + \sum_i u_i G_i) x\]

where:

  • G₀ is the drift term (dynamics with no control)
  • Gᵢ are the drive terms (how controls affect the system)
  • uᵢ are the control inputs

How it Works

Uses the matrix exponential for exact integration:

\[x_{k+1} = \exp(\Delta t \cdot G(u_k)) x_k\]

where G(u) = G₀ + Σᵢ uᵢ Gᵢ.

Example: Simple 2D System

N = 50
traj = NamedTrajectory(
    (x = randn(2, N), u = randn(1, N), Δt = fill(0.1, N));
    timestep = :Δt,
    controls = :u,
    initial = (x = [1.0, 0.0],),
    final = (x = [0.0, 1.0],),
)
N = 50, (x = 1:2, u = 3:3, → Δt = 4:4)

Define drift (natural dynamics) and drives (control terms)

G_drift = [-0.1 1.0; -1.0 -0.1]     # Damped oscillator
G_drives = [[0.0 1.0; 1.0 0.0]]     # Symmetric control coupling
1-element Vector{Matrix{Float64}}:
 [0.0 1.0; 1.0 0.0]

Create generator function

G = u -> G_drift + sum(u .* G_drives)
#2 (generic function with 1 method)

Create integrator

integrator = BilinearIntegrator(G, :x, :u, traj)
BilinearIntegrator: :x = exp(Δt G(:u)) :x  (dim = 2)

Multiple Drives Example

traj_multi = NamedTrajectory(
    (x = randn(3, N), u = randn(2, N), Δt = fill(0.1, N));
    timestep = :Δt,
    controls = :u,
)

G_drift_3d = [
    0.0 1.0 0.0;
    -1.0 0.0 0.0;
    0.0 0.0 -0.1
]

G_drives_3d = [
    [1.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0],  # Drive 1
    [0.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0],   # Drive 2
]

G_multi = u -> G_drift_3d + sum(u .* G_drives_3d)

integrator_multi = BilinearIntegrator(G_multi, :x, :u, traj_multi)
BilinearIntegrator: :x = exp(Δt G(:u)) :x  (dim = 3)

When to Use BilinearIntegrator

✓ Quantum systems (Hamiltonian evolution) ✓ Rotating systems (attitude dynamics) ✓ Systems linear in controls ✓ When you want exact integration (no discretization error)

TimeDependentBilinearIntegrator

Overview

For time-varying bilinear dynamics:

\[\dot{x} = (G_0(t) + \sum_i u_i(t) G_i(t)) x\]

The generator function now depends on both control and time.

Example: Periodic Disturbance

traj_td = NamedTrajectory(
    (
        x = randn(2, N),
        u = randn(1, N),
        t = collect(range(0, 5, N)),  # time variable
        Δt = fill(0.1, N),
    );
    timestep = :Δt,
    controls = :u,
)
N = 50, (x = 1:2, u = 3:3, t = 4:4, → Δt = 5:5)

Time-dependent generator

G_td = (u, t) -> [-0.1 + 0.5*sin(t) 1.0; -1.0 -0.1] + u[1] * [0.0 1.0; 1.0 0.0]

integrator_td = TimeDependentBilinearIntegrator(G_td, :x, :u, :t, traj_td)
TimeDependentBilinearIntegrator: :x via G(:u, t)  (dim = 2, order = 1)

When to Use TimeDependentBilinearIntegrator

✓ Time-varying Hamiltonians ✓ Systems with periodic forcing ✓ Carrier wave modulation (e.g., rotating frame transformations)

DerivativeIntegrator

Overview

Enforces derivative relationships between trajectory components:

\[\frac{d(\text{var})}{dt} = \text{deriv}\]

This is used for smoothness or when controls are derivatives of other variables.

Example: Smooth Controls

traj_smooth = NamedTrajectory(
    (
        x = randn(2, N),
        u = randn(2, N),
        du = zeros(2, N),   # control derivative
        Δt = fill(0.1, N),
    );
    timestep = :Δt,
    controls = :u,
    initial = (u = [0.0, 0.0],),
    final = (u = [0.0, 0.0],),
)
N = 50, (x = 1:2, u = 3:4, du = 5:6, → Δt = 7:7)

Enforce du/dt = du

deriv_integrator = DerivativeIntegrator(:u, :du, traj_smooth)
DerivativeIntegrator: :u += Δt * :du  (dim = 2)

Now you can penalize du to get smooth controls: obj = QuadraticRegularizer(:u, trajsmooth, 1e-2) obj += QuadraticRegularizer(:du, trajsmooth, 1e-1) # Smoothness penalty

Multiple Derivative Orders

traj_smooth2 = NamedTrajectory(
    (
        x = randn(2, N),
        u = randn(1, N),
        du = zeros(1, N),
        ddu = zeros(1, N),
        Δt = fill(0.1, N),
    );
    timestep = :Δt,
    controls = :u,
)
N = 50, (x = 1:2, u = 3:3, du = 4:4, ddu = 5:5, → Δt = 6:6)

Chain derivatives: d(u)/dt = du, d(du)/dt = ddu

deriv_u = DerivativeIntegrator(:u, :du, traj_smooth2)
deriv_du = DerivativeIntegrator(:du, :ddu, traj_smooth2)
DerivativeIntegrator: :du += Δt * :ddu  (dim = 1)

When to Use DerivativeIntegrator

✓ Enforce smooth, implementable controls ✓ Acceleration limits (when control is jerk) ✓ Tracking derivative information

Combining Multiple Integrators

You can use multiple integrators simultaneously:

traj_combined = NamedTrajectory(
    (x = randn(2, N), u = randn(2, N), du = zeros(2, N), Δt = fill(0.1, N));
    timestep = :Δt,
    controls = :u,
    initial = (x = [0.0, 0.0], u = [0.0, 0.0]),
    final = (u = [0.0, 0.0],),
)
N = 50, (x = 1:2, u = 3:4, du = 5:6, → Δt = 7:7)

Create problem with multiple integrators

G_drift_simple = [-0.1 1.0; -1.0 -0.1]
G_drives_simple = [[0.0 1.0; 1.0 0.0], [1.0 0.0; 0.0 1.0]]
G_simple = u -> G_drift_simple + sum(u .* G_drives_simple)

obj = QuadraticRegularizer(:u, traj_combined, 1e-2)
obj += QuadraticRegularizer(:du, traj_combined, 1e-1)

integrators_combined = [
    BilinearIntegrator(G_simple, :x, :u, traj_combined),
    DerivativeIntegrator(:u, :du, traj_combined),
]

prob = DirectTrajOptProblem(traj_combined, obj, integrators_combined)
DirectTrajOptProblem
  Trajectory
    Timesteps: 50
    Duration:  4.9
    Knot dim:  7
    Variables: x (2), u (2), du (2), Δt (1)
    Controls:  u, Δt
  Objective (2 terms)
         1.0 * QuadraticRegularizer on :u (R = [0.01, 0.01], all)
         1.0 * QuadraticRegularizer on :du (R = [0.1, 0.1], all)
  Dynamics (2 integrators)
    BilinearIntegrator: :x = exp(Δt G(:u)) :x  (dim = 2)
    DerivativeIntegrator: :u += Δt * :du  (dim = 2)
  Constraints (4 total: 3 equality, 1 bounds)
    EqualityConstraint: "initial value of x"
    EqualityConstraint: "initial value of u"
    EqualityConstraint: "final value of u"
    BoundsConstraint: "bounds on Δt"

Integration Methods Comparison

IntegratorDynamics TypeAccuracyUse Case
BilinearIntegratorControl-linearExactQuantum, rotation
TimeDependentBilinearIntegratorTime-varying control-linearExactModulated systems
DerivativeIntegratorDerivative relationExactSmoothness

Custom Integrators

You can implement custom integrators by subtyping AbstractIntegrator and defining the constraint function. See the Advanced Topics section for details.

Interface Requirements

struct MyIntegrator <: AbstractIntegrator
    # ... fields ...
end

# Implement constraint evaluation
function (int::MyIntegrator)(δ, zₖ, zₖ₊₁, k)
    # Compute constraint: δ = xₖ₊₁ - Φ(xₖ, uₖ, Δtₖ)
    # where Φ is your integration scheme
end

Best Practices

Initialization

  • Start with good initial guesses for states and controls
  • For smooth control problems, initialize derivatives to zero
  • Use linear interpolation for states between boundary conditions

Performance

  • Matrix exponential (BilinearIntegrator) is efficient for small systems (n < 20)
  • For large systems, consider sparse representations
  • DerivativeIntegrator is cheap (just finite differences)

Numerical Stability

  • Keep time steps reasonable (not too large)
  • For stiff systems, smaller time steps help
  • BilinearIntegrator handles stiff systems well

Common Patterns

Pattern 1: Basic Bilinear Problem

G_basic = u -> [-0.1 1.0; -1.0 -0.1] + u[1] * [0.0 1.0; 1.0 0.0]
#14 (generic function with 1 method)

integrator = BilinearIntegrator(G_basic, :x, :u, traj)

Pattern 2: Smooth Control Problem

integrators = [ BilinearIntegrator(G, :x, :u, traj), DerivativeIntegrator(:u, :du, traj) ]

Pattern 3: Time-Dependent with Smoothness

integrators = [ TimeDependentBilinearIntegrator(G_td, :x, :u, :t, traj), DerivativeIntegrator(:u, :du, traj) ]

Summary

Key Takeaways:

  1. Integrators convert continuous dynamics to discrete constraints
  2. BilinearIntegrator is the workhorse for control-linear systems
  3. DerivativeIntegrator adds smoothness
  4. You can combine multiple integrators
  5. Good initialization helps convergence

Next Steps

  • Objectives: Learn how to define cost functions
  • Constraints: Add bounds and path constraints
  • Tutorials: See integrators in complete examples

This page was generated using Literate.jl.