OptimalControl.jl

Appearance

Flows overview ​

Flow is one constructor that does three distinct jobs: indirect optimal control (build the Hamiltonian flow of the Pontryagin Maximum Principle, write a shooting function, solve it), simulation (integrate a controlled system under an open-loop or feedback control), and inspection (pull the Hamiltonian, the Hamiltonian vector field, or the control law back out of a flow you built). This page maps the section and recaps just enough PMP to make the rest make sense.

Why flows ​

The direct methods (Solve) discretize the whole problem into one large nonlinear program. Flows integrate instead: you supply a control law — from the PMP, from a feedback design, from anywhere — and get back an ODE solution. Direct methods win when you don't yet have a control law and want the solver to find one; flows win once you do, whether that law came from solving the PMP by hand, from a direct solve's costate, or from a controller you designed independently of any optimization.

The Pontryagin Maximum Principle, briefly ​

The notation is the one of Notation and conventions. For a problem with dynamics   and Lagrange cost , the pseudo-Hamiltonian is

with    in the normal case. Along an optimal trajectory, the control maximises . When this maximisation gives the control in feedback form, a control law , substituting it gives the maximised Hamiltonian  , whose Hamiltonian system

gives the extremals. Its flow,   from , is what Flow computes. The boundary and transversality conditions then turn "integrate the flow" into "find the missing ": that is shooting.

From the PMP to a flow ​

You supply the control law, worked out by hand, and Flow builds the Hamiltonian flow. For the energy-minimal double integrator below,     is maximised by  , the second costate component: the control law is (x, p) -> p[2]. From an OCP explains how Flow(ocp, law) reads the dynamics and the cost from the problem.

Three things this section does ​

From Hamiltonians covers the lower-level constructors these three jobs are all built from.

Before you start ​

Every flow needs an ODE integrator, and none is a hard dependency — load one, most commonly OrdinaryDiffEqTsit5, before building any flow:

julia
using OptimalControl
using OrdinaryDiffEqTsit5

t0 = 0
tf = 1
x0 = [-1, 0]

ocp = @def begin
    t ∈ [t0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    x(t0) == x0
    x(tf) == [0, 0]
    ẋ(t) == [v(t), u(t)]
    0.5∫(u(t)^2) → min
end

f = Flow(ocp, (x, p) -> p[2])

A flow is called in two ways. At a point, it returns the state and the costate at the final time. From the exact initial costate  , it reaches the target  , with   :

julia
xf, pf = f(t0, x0, [12, 6], tf)
([-4.158277980903794e-16, 5.038973481527467e-15], [12.0, -5.999999999999998])

On a time span, it returns the whole trajectory, as a Solution that can be read and plotted like the solution of a direct solve:

julia
sol = f((t0, tf), x0, [12, 6])
objective(sol)
5.999999999999946

Every page in this section opens with using OrdinaryDiffEqTsit5. Without it, Flow raises an error at construction, naming the package to load:

julia
julia> f = Flow(ocp, (x, p) -> p[2])
ERROR: ExtensionError → top-level scope, REPL[6]:1
│
│  missing dependencies to access SciML options metadata
│
│  Missing  OrdinaryDiffEqTsit5
│
│  Context  Load OrdinaryDiffEqTsit5, OrdinaryDiffEq, or DifferentialEquations to activate the CTSolversSciMLIntegrator extension.
│  Hint     Run: using OrdinaryDiffEqTsit5
└─

Choosing an integrator and its options ​

describe covers the indirect side too, not just the direct-solve strategies from Choosing a method: :sciml for the integrator family, :di for the automatic-differentiation backend that builds a Hamiltonian's vector field.

describe(:sciml)
julia
describe(:sciml)
SciML (strategy)
├─ id: :sciml
├─ hierarchy: SciML → AbstractSciMLIntegrator → AbstractIntegrator → AbstractStrategy
├─ description: SciML ODE integrator.
│               See: https://docs.sciml.ai/DiffEqDocs
│               Solver options: https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/
├─ family: AbstractIntegrator
├─ default parameter: CPU
├─ parameters: CPU, GPU
│
├─ computed options for CPU: none
│
├─ computed options for GPU: none
│
└─ common options (20 options):
   ├─ maxiters (max_iters, max_iter, maxiter, max_iterations, maxit)::Integer (default: NotProvided)
   │  description: Maximum number of solver iterations.
   │
   ├─ save_everystep::Union{Bool, Symbol} (default: auto)
   │  description: Save the solution at every solver step. Set `true`/`false` to force, or `:auto` to resolve to `false` for point integration and `true` for trajectory integration.
   │
   ├─ save_end::Bool (default: NotProvided)
   │  description: Whether to force saving the final timepoint.
   │
   ├─ d_discontinuities::AbstractVector{<:Real} (default: NotProvided)
   │  description: Locations of discontinuities in low-order derivatives.
   │
   ├─ save_start::Union{Bool, Symbol} (default: auto)
   │  description: Save initial condition in solution. Set `true`/`false` to force, or `:auto` to resolve to `false` for point integration and `true` for trajectory integration.
   │
   ├─ force_dtmin::Bool (default: NotProvided)
   │  description: Whether to continue forcing minimum dt usage.
   │
   ├─ alg (algorithm, solver)::Union{Missing, SciMLBase.AbstractDEAlgorithm} (default: Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), FastBroadcast.Serial}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, FastBroadcast.Serial()))
   │  description: ODE algorithm (e.g. Tsit5(), Vern6()).
   │
   ├─ saveat (save_at, save_times)::Union{Real, AbstractVector} (default: NotProvided)
   │  description: Times at which to save the solution (Vector or range).
   │
   ├─ dtmin (min_dt, dt_min)::Real (default: NotProvided)
   │  description: Minimum step size for adaptive timestepping.
   │
   ├─ progress (verbose)::Bool (default: NotProvided)
   │  description: Whether to show progress bar.
   │
   ├─ reltol (rtol, rel_tol)::Real (default: 1.0e-8)
   │  description: Relative tolerance for the ODE solver.
   │
   ├─ adaptive (adaptive_step, adaptive_stepping)::Bool (default: NotProvided)
   │  description: Whether to use adaptive step-size control.
   │
   ├─ dt (dt0, timestep)::Real (default: NotProvided)
   │  description: Fixed step size (used when adaptive=false).
   │
   ├─ save_idxs (saveindices, save_indices)::AbstractVector{<:Integer} (default: NotProvided)
   │  description: Indices of components to save (Vector of integers).
   │
   ├─ callback (callbacks, cb)::Any (default: NotProvided)
   │  description: Callback function for event handling.
   │
   ├─ tstops (t_stops, stop_times)::AbstractVector{<:Real} (default: NotProvided)
   │  description: Extra times the solver must step to (for discontinuities).
   │
   ├─ internalnorm (internal_norm, norm)::Function (default: real_norm)
   │  description: Internal norm for adaptive step-size control. Defaults to `real_norm`, which extracts the primal (Float64) part of ForwardDiff dual numbers to ensure grid invariance (IND) when ForwardDiff is loaded. Set to `DiffEqBase.ODE_DEFAULT_NORM` to use the SciML default.
   │
   ├─ abstol (atol, abs_tol)::Real (default: 1.0e-8)
   │  description: Absolute tolerance for the ODE solver.
   │
   ├─ dense::Union{Bool, Symbol} (default: auto)
   │  description: Dense output. Set `true`/`false` to force, or `:auto` to resolve to `false` for point integration and `true` for trajectory integration.
   │
   └─ dtmax (max_dt, dt_max)::Real (default: NotProvided)
      description: Maximum step size for adaptive timestepping.

The options used most are alg (the ODE algorithm, Tsit5() by default), reltol and abstol (1e-8 each by default), saveat (the times at which to save the trajectory) and dense (dense output, :auto by default: false for a point call, true for a trajectory call). Pass them as keywords when you build a flow, for example Flow(ocp, law; reltol=1e-10, alg=Tsit5()). With saveat, also pass dense=false (CTFlows#434), and use the flow for trajectories only: its point call returns the value at the last saveat time (CTFlows#441). See Plot.

Another integrator only needs its package. Here the 9th-order Verner method, with tighter tolerances:

julia
using OrdinaryDiffEqVerner
f9 = Flow(ocp, (x, p) -> p[2]; alg=Vern9(), reltol=1e-12, abstol=1e-12)
f9(t0, x0, [12, 6], tf)
([-3.8459654210867653e-17, -2.1355478518691923e-16], [12.0, -6.0])
describe(:di)
julia
describe(:di)
DifferentiationInterface (strategy)
├─ id: :di
├─ hierarchy: DifferentiationInterface → AbstractADBackend → AbstractStrategy
├─ description: AD backend wrapping DifferentiationInterface.jl backends (e.g., AutoForwardDiff).
├─ family: AbstractADBackend
├─ default parameter: CPU
├─ parameters: CPU, GPU
│
├─ computed options for CPU:
│  └─ ad_backend (backend, ad)::ADTypes.AbstractADType (default: AutoForwardDiff() [computed])
│     description: DifferentiationInterface.jl backend (e.g. AutoForwardDiff() on CPU).
│
└─ computed options for GPU:
   └─ ad_backend (backend, ad)::ADTypes.AbstractADType (default: ADTypes.AutoMooncake() [computed])
      description: DifferentiationInterface.jl backend for GPU execution. Default: AutoMooncake(), validated end-to-end on CuArray including through a mutating in-place RHS. AutoForwardDiff() does not work on GPU (scalar-indexes a CuArray). AutoZygote() can be selected explicitly; it was found unreliable in some device call contexts (an unexpected KernelException was observed on a non-mutating call).

CPU and GPU ​

The parameters :cpu and :gpu, the same as for solve, are passed with the method keyword when you build a flow. The call has no such keyword:

julia
f = Flow(ocp, law; method=:gpu)   # a keyword of the constructor
f(t0, x0, p0, tf; method=:gpu)    # MethodError: not a keyword of the call

method=:gpu selects the GPU variants of the integrator and of the automatic differentiation backend (Mooncake); the initial state and costate must then be device arrays (CuArray). See the GPU flows page of CTFlows for what it supports. On a machine without a GPU, or without Mooncake loaded, the error comes at the first call and does not name the cause yet (CTFlows#436).

Where to go ​