OptimalControl.jl

Appearance

GPU ​

GPU support runs through ExaModels.jl and MadNLPGPU.jl, NVIDIA GPUs only, via CUDA.jl.

What you are reading depends on the machine that built this page

Every block below is executed when the documentation is built. With a functional CUDA device you are reading real GPU output; without one, you are reading the failure this exact code really produces. The first block says which of the two it is.

Prerequisites ​

julia
using OptimalControl
using MadNLPGPU
using CUDA
using CUDSS

println("CUDA.functional() = ", CUDA.functional())
CUDA.functional() = true

Check CUDA.functional() before assuming a :gpu solve will actually run on the device.

All three — and CUDSS is the one you will forget

The CTSolversMadNLPGPU extension is armed by MadNLPGPU, CUDA and CUDSS together. Load only the first two — the pair every GPU tutorial shows — and the extension does not load, so the GPU solver strategies are never registered.

It used to work by accident. Up to MadNLPGPU 0.8, CUDSS was a hard dependency, so using MadNLPGPU pulled it in and the third trigger was satisfied without anyone asking. From 0.9 onward it is a weak dependency and you must load it yourself.

The error names exactly which one is missing — load MadNLPGPU and CUDA but not CUDSS, and it reports Missing CUDSS with the hint using CUDSS.

ExaModels needs no using of its own. It is a dependency of OptimalControl, so the module is already bound after using OptimalControl and :exa works without it. Importing it explicitly also brings its objective and constraint into scope, both of which collide with the accessors of the same name — see #882. If you do need ExaModels' own API, import it qualified: using ExaModels: ExaModels.

The problem must be coordinatewise ​

:exa — the only GPU-capable modeler — requires dynamics (and any path constraint) written one coordinate at a time, ∂(x₁)(t) == ..., not ẋ(t) == [...]. See Abstract syntax (@def) for the two forms side by side.

julia
ocp = @def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    v ∈ R, variable
    x(0) == [0, 1]
    x(1) == [0, -1]
    ∂(x₁)(t) == x₂(t)   # coordinatewise
    ∂(x₂)(t) == u(t)    # — not ẋ(t) == [x₂(t), u(t)]
    0 ≤ x₁(t) + v^2 ≤ 1.1
    -10 ≤ u(t) ≤ 10
    1 ≤ v ≤ 2
    ∫(u(t)^2 + v) → min
end
Abstract definition:

    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    v ∈ R, variable
    x(0) == [0, 1]
    x(1) == [0, -1]
    (∂(x₁))(t) == x₂(t)
    (∂(x₂))(t) == u(t)
    0 ≤ x₁(t) + v ^ 2 ≤ 1.1
    -10 ≤ u(t) ≤ 10
    1 ≤ v ≤ 2
    ∫(u(t) ^ 2 + v) → min

The (autonomous) optimal control problem is of the form:

    minimize  J(x, u, v) = ∫ f⁰(x(t), u(t), v) dt, over [0, 1]

    subject to

        ẋ(t) = f(x(t), u(t), v), t in [0, 1] a.e.,

        ψ₋ ≤ ψ(x(t), u(t), v) ≤ ψ₊, 
        ϕ₋ ≤ ϕ(x(0), x(1), v) ≤ ϕ₊, 
        u₋ ≤ u(t) ≤ u₊, 
        v₋ ≤ v ≤ v₊, 

    where x(t) ∈ R², u(t) ∈ R and v ∈ R.

Descriptive mode ​

The :gpu parameter token selects GPU-optimized defaults:

julia
try
    global sol = solve(ocp, :exa, :madnlp, :gpu; grid_size=100, display=false)
    println("objective  = ", objective(sol))
    println("iterations = ", iterations(sol))
catch e
    println("GPU solve failed — no functional device on this machine.")
    println("CUDA.functional() = ", CUDA.functional())
    println("Exception: ", first(sprint(showerror, e), 400))
end
objective  = 9.891292172085043
iterations = 20

Completion fills in the rest — the first match containing :gpu is the same method:

julia
try
    global sol = solve(ocp, :gpu; grid_size=100, display=false)
    println("objective = ", objective(sol))
catch e
    println("Exception: ", first(sprint(showerror, e), 400))
end
objective = 9.891292172085041

Solver verbosity is a separate concern: display=false above silences the OptimalControl-level report, and the underlying solver takes its own options — print_level=MadNLP.ERROR for MadNLP, which needs using MadNLP in scope. See Options.

:gpu changes what a strategy's own defaults are: Exa{GPU} uses a CUDA differentiation backend, MadNLP{GPU} uses the CUDSSSolver linear solver instead of MUMPS. describe(:gpu) lists every strategy with a GPU-parameterized variant — this call needs nothing GPU-specific and runs fine on CPU alone:

julia
describe(:gpu)
GPU (parameter)
├─ id: :gpu
├─ hierarchy: GPU → AbstractStrategyParameter
├─ description: GPU-based computation
│
└─ used by strategies (5):
   ├─ :di (AbstractADBackend) → DifferentiationInterface{O}
   ├─ :exa (AbstractNLPModeler) → Exa{GPU}
   ├─ :madncl (AbstractNLPSolver) → MadNCL{GPU}
   ├─ :madnlp (AbstractNLPSolver) → MadNLP{GPU}
   └─ :sciml (AbstractIntegrator) → SciML{O, OP, OT}

Explicit mode ​

Constructing the components does not touch the device, so this block runs anywhere:

julia
disc = OptimalControl.Collocation(; grid_size=100, scheme=:midpoint)
mod = OptimalControl.Exa{GPU}()
slv = OptimalControl.MadNLP{GPU}()

Running them is what needs the hardware:

julia
try
    global result = solve(ocp; discretizer=disc, modeler=mod, solver=slv)
    println("objective = ", objective(result))
catch e
    println("Exception: ", first(sprint(showerror, e), 400))
end
▫ This is OptimalControl 2.2.1-beta, solving with: collocation → exa → madnlp (gpu)

  📦 Configuration:
   ├─ Discretizer: collocation (grid_size = 100, scheme = midpoint)
   ├─ Modeler: exa (backend = CUDACore.CUDAKernels.CUDABackend(false, false) [gpu-dependent])
   └─ Solver: madnlp (linear_solver = MadNLPGPUCUDAExt.CUDSSSolver [gpu-dependent])

▫ This is MadNLP version v0.10.1, running with cuDSS v0.7.1

Number of nonzeros in constraint Jacobian............:      906
Number of nonzeros in Lagrangian Hessian.............:      501

Total number of variables............................:      304
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      102
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:      305
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:      305
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
   0  1.0200000e+00 1.10e+00 1.00e+00 1.01e+01  -1.0     -   0.00e+00  0  0
   1  1.0199995e+00 1.10e+00 2.41e-01 3.51e-02  -1.0     -   2.00e-07  1  2h
   2  1.0221486e+00 9.80e-01 3.77e-02 4.71e-03  -1.0     -   1.09e-01  2  1h
   3  1.0255839e+00 9.74e-01 1.45e-01 1.42e-03  -1.7     -   5.77e-03  2  1h
   4  6.3410730e+00 2.21e-01 1.96e-02 4.95e-04  -1.7     -   7.73e-01  2  1h
   5  1.0369076e+01 5.24e-11 8.48e-03 1.55e-04  -1.7     -   1.00e+00  2  1h
   6  1.0122712e+01 2.64e-10 6.42e-02 2.04e-04  -3.8     -   1.00e+00  2  1h
   7  9.9153817e+00 4.36e-10 1.83e-02 7.71e-05  -3.8     -   8.43e-01  2  1h
   8  9.9017479e+00 1.60e-10 2.12e-01 3.18e-05  -3.8     -   6.60e-01  2  1h
   9  9.8962247e+00 3.51e-12 1.44e-01 1.26e-05  -3.8     -   1.00e+00  2  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  10  9.8953501e+00 4.44e-16 7.87e-11 5.22e-06  -3.8     -   1.00e+00  2  1h
  11  9.8928497e+00 2.64e-13 3.14e+00 6.10e-05  -8.6     -   6.14e-01  2  1h
  12  9.8920064e+00 1.57e-13 3.29e+00 3.06e-05  -8.6     -   5.06e-01  1  1h
  13  9.8915281e+00 6.24e-14 1.02e+00 1.48e-05  -8.6     -   6.85e-01  1  1h
  14  9.8913591e+00 1.38e-14 9.54e-01 5.47e-06  -8.6     -   8.13e-01  1  1h
  15  9.8913055e+00 2.22e-16 1.44e+00 1.72e-06  -8.6     -   1.00e+00  2  1h
  16  9.8912977e+00 2.22e-16 8.02e-01 5.18e-07  -8.6     -   1.00e+00  1  1h
  17  9.8912941e+00 2.22e-16 8.21e-01 1.49e-06  -8.6     -   1.00e+00  2  1h
  18  9.8912924e+00 2.22e-16 1.18e-01 4.23e-07  -8.6     -   1.00e+00  2  1h
  19  9.8912922e+00 2.22e-16 4.84e-03 1.84e-08  -8.6     -   1.00e+00  2  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  20  9.8912922e+00 2.22e-16 5.68e-14 4.16e-09  -8.6     -   1.00e+00  2  1h

Number of Iterations....: 20

                                   (scaled)                 (unscaled)
Objective...............:   9.8912921720850395e+00    9.8912921720850395e+00
Dual infeasibility......:   5.6843418860808015e-14    5.6843418860808015e-14
Constraint violation....:   2.2204460492503131e-16    2.2204460492503131e-16
Complementarity.........:   4.1623869934921173e-09    4.1623869934921173e-09
Overall NLP error.......:   4.1623869934921173e-09    4.1623869934921173e-09

Number of objective function evaluations              = 22
Number of objective gradient evaluations              = 21
Number of constraint evaluations                      = 22
Number of constraint Jacobian evaluations             = 21
Number of Lagrangian Hessian evaluations              = 20
Number of KKT factorizations                          = 20
Number of KKT backsolves                              = 35

Total wall secs in initialization                     =  0.026 s
Total wall secs in linear solver                      =  unavailable
Total wall secs in NLP function evaluations           =  0.013 s
Total wall secs in solver (w/o init./fun./lin. alg.)  =  unavailable
Total wall secs                                       =  0.116 s

EXIT: Optimal Solution Found (tol = 1.0e-08).
objective = 9.89129217208504

What combinations work ​

Only :exa × {:madnlp, :madncl} on :gpu — the two :gpu entries of methods() (see Choosing a method). Everything else is a compile-time or runtime error. These are type-system and routing errors, not hardware-dependent ones, so they raise identically on every machine:

ADNLP's parameter is constrained to <:CPU:

julia
julia> OptimalControl.ADNLP{GPU}()
TypeError: in ADNLP, in P, expected P<:CPU, got Type{GPU}

Same for Ipopt:

julia
julia> OptimalControl.Ipopt{GPU}()
TypeError: in Ipopt, in P, expected P<:CPU, got Type{GPU}

Descriptively, the same combinations fail earlier still — no entry in methods() carries :adnlp together with :gpu, so completion cannot resolve the description at all:

julia
julia> solve(ocp, :adnlp, :gpu)
AmbiguousDescription → #complete#15, complete.jl:110
│
│  cannot find matching description
│
│  Diagnostic  No complete match — no description contains all symbols
│  Requested   (:adnlp, :gpu)
│  Available   (:collocation, :adnlp, :ipopt, :cpu)
│              (:collocation, :adnlp, :madnlp, :cpu)
│              (:collocation, :adnlp, :uno, :cpu)
│              (:collocation, :adnlp, :madncl, :cpu)
│              (:collocation, :adnlp, :knitro, :cpu)
│              (:collocation, :exa, :ipopt, :cpu)
│              (:collocation, :exa, :madnlp, :cpu)
│              (:collocation, :exa, :uno, :cpu)
│              (:collocation, :exa, :madncl, :cpu)
│              (:collocation, :exa, :knitro, :cpu)
│              (:collocation, :exa, :madnlp, :gpu)
│              (:collocation, :exa, :madncl, :gpu)
│
│  Context     description completion
│  Hint        Try one of the closest matches:
│              (:collocation, :exa, :madnlp, :gpu)
│              (:collocation, :exa, :madncl, :gpu)
│              (:collocation, :adnlp, :uno, :cpu)
│              (:collocation, :adnlp, :madnlp, :cpu)
│              (:collocation, :adnlp, :madncl, :cpu)
└─

Performance notes ​

GPU solving amortizes best on large-scale problems (thousands of variables/constraints) or repeated solves in a loop, where the per-call setup overhead is paid once. For small problems, plain CPU solving is typically faster.

The idiomatic guard is CUDA.functional() — pick the strategy, then solve:

julia
strategy = CUDA.functional() ? :gpu : :cpu
println("strategy = ", strategy)
strategy = gpu
julia
if CUDA.functional()
    t = @elapsed solve(ocp, :gpu; grid_size=1000, display=false)
    println("GPU solve at grid_size=1000: ", round(t; digits=2), " s")
else
    println("No functional device here, so there is no GPU timing to report.")
    println("On a CUDA machine this block prints the :gpu solve time at grid_size=1000.")
end
GPU solve at grid_size=1000: 0.55 s

See also ​