Abstract syntax (`@def`)

The full grammar of OptimalControl.jl's small Domain Specific Language is given below. The idea is to use a syntax that is

pure Julia (and, as such, effortlessly analysed by the standard Julia parser),
as close as possible to the mathematical description of an optimal control problem (see Formulation).

While the syntax will be transparent to those users familiar with Julia expressions (Expr's), we provide examples for every case that should be widely understandable. Abstract definitions use the macro @def.

About the code blocks on this page

Blocks written :( … ) — such as the state pattern :( $x ∈ R^$n, state ) — are grammar patterns: the shape of the Julia expression the parser accepts, with the $-prefixed names standing for the parts you supply. They are notation, not runnable code, and are shown inert on purpose.

Every other block is executed when the documentation is built. A problem definition must include at least time, state, dynamics and cost (the control declaration is optional — see Problems without a control); in the blocks that illustrate a single clause, the other required clauses are present but hidden so that the block still builds a complete model.

Variable

julia

:( $v ∈ R^$q, variable )
:( $v ∈ R   , variable )

A variable (only one is allowed) is a finite dimensional vector of reals that will be optimised along with state and control values. To define an optimal control problem with a dimension-two variable named v:

julia

@def begin
    v ∈ R², variable
end

Aliases v₁, v₂ (and v1, v2) are automatically defined and can be used in subsequent expressions instead of v[1] and v[2]. You can also define your own aliases for the components (one alias per dimension):

julia

@def begin
    v = (a, b) ∈ R², variable
end

A one dimensional variable can be declared according to

julia

@def begin
    v ∈ R, variable
end

Warning

Aliases during definition of variable, state or control are only allowed for multidimensional (dimension two or more) cases. Something like u = T ∈ R, control is not allowed... and useless (directly write T ∈ R, control).

Time

julia

:( $t ∈ [$t0, $tf], time )

The independent variable or time is a scalar bound to a given interval. Its name is arbitrary.

julia

t0 = 1
tf = 5
@def begin
    t ∈ [t0, tf], time
end

One (or even the two bounds) can be variable, typically for minimum time problems (see Mayer cost section):

julia

@def begin
    v = (T, λ) ∈ R², variable
    t ∈ [0, T], time
end

State

julia

:( $x ∈ R^$n, state )
:( $x ∈ R   , state )

The state declaration defines the name and the dimension of the state:

julia

@def begin
    x ∈ R⁴, state
end

As for the variable, there are automatic aliases (x₁ and x1 for x[1], etc.) and you can define your own aliases (one per scalar component of the state):

julia

@def begin
    x = (q₁, q₂, v₁, v₂) ∈ R⁴, state
end

Control

julia

:( $u ∈ R^$m, control )
:( $u ∈ R   , control )

The control declaration defines the name and the dimension of the control:

julia

@def begin
    u ∈ R², control
end

As before, there are automatic aliases (u₁ and u1 for u[1], etc.) and you can define your own aliases (one per scalar component of the control):

julia

@def begin
    u = (α, β) ∈ R², control
end

Note

One dimensional variable, state or control are treated as scalars (Real), not vectors (Vector). In Julia, for x::Real, it is possible to write x[1] (and x[1][1]...) so it is OK (though useless) to write x₁, x1 or x[1] instead of simply x to access the corresponding value. Conversely it is not OK to use such an x as a vector, for instance as in ...f(x)... where f(x::Vector{T}) where {T <: Real}. This same convention applies on the functional API: a dimension-1 quantity is a scalar, not a length-1 vector — see Shapes in callbacks.

Problems without a control

The control declaration is optional. You can define problems without control for:

Parameter estimation: identify unknown parameters in the dynamics from observed data,
Dynamic optimisation: optimise constant parameters subject to ODE constraints.

For example, to estimate a growth rate parameter (here data is the observed signal, a function defined beforehand):

julia

@def begin
    p ∈ R, variable                 # parameter to estimate
    t ∈ [0, 10], time
    x ∈ R, state
    x(0) == 2.0
    ẋ(t) == p * x(t)                # dynamics depends on p
    ∫((x(t) - data(t))^2) → min     # fit to observed data
end

Or to optimise the pulsation of a harmonic oscillator:

julia

@def begin
    ω ∈ R, variable                 # pulsation to optimise
    t ∈ [0, 1], time
    x = (q, v) ∈ R², state
    q(0) == 1.0
    v(0) == 0.0
    q(1) == 0.0                     # final condition
    ẋ(t) == [v(t), -ω^2 * q(t)]     # harmonic oscillator
    ω^2 → min                       # minimise pulsation
end

There is no dedicated syntax for "no control": omitting the control declaration entirely is the syntax — do not declare a dummy u ∈ R, control with u(t) == 0, and never write control!(pre, 0) on the functional API, which is rejected as an error. See No control for the full worked examples above, solved both directly and indirectly.

Dynamics

julia

:( ∂($x)($t) == $e1 )

The dynamics is given in the standard vectorial ODE form:

depending on whether it is autonomous / with a variable or not (the parser will detect time and variable dependences, which entails that time, state and variable must be declared prior to dynamics - an error will be issued otherwise). The symbol ∂, or the dotted state name (ẋ), or the keyword derivative can be used:

julia

@def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    ∂(x)(t) == [x₂(t), u(t)]
    ∫(u(t)^2) → min
end

julia

@def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    ẋ(t) == [x₂(t), u(t)]
    ∫(u(t)^2) → min
end

julia

@def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    derivative(x)(t) == [x₂(t), u(t)]
    ∫(u(t)^2) → min
end

Any Julia code can be used, so the following is also OK:

julia

ocp = @def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    ẋ(t) == F₀(x(t)) + u(t) * F₁(x(t))
    ∫(u(t)^2) → min
end

F₀(x) = [x[2], 0]
F₁(x) = [0, 1]

Note

The vector fields F₀ and F₁ can be defined afterwards, as they only need to be available when the dynamics will be evaluated.

While it is also possible to declare the dynamics component after component (see below), one may equivalently use aliases (check the relevant aliases section below):

julia

@def damped_integrator begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    q̇ = v(t)
    v̇ = u(t) - c(t)
    ẋ(t) == [q̇, v̇]
    ∫(u(t)^2) → min
end

Dynamics (coordinatewise)

julia

:( ∂($x[$i])($t) == $e1 )

The dynamics can also be declared coordinate by coordinate. The previous example can be written as

julia

@def damped_integrator begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    ∂(q)(t) == v(t)
    ∂(v)(t) == u(t) - c(t)
    ∫(u(t)^2) → min
end

Constraints

julia

:( $e1 == $e2        )
:( $e1 ≤  $e2 ≤  $e3 )
:(        $e2 ≤  $e3 )
:( $e3 ≥  $e2 ≥  $e1 )
:( $e2 ≥  $e1        )

Admissible constraints can be

of five types: boundary, variable, control, state, mixed (the last three ones are path constraints, that is constraints evaluated all times)
linear (ranges) or nonlinear (not ranges),
equalities or (one or two-sided) inequalities.

Boundary conditions are detected when the expression contains evaluations of the state at initial and / or final time bounds (e.g., x(0)), and may not involve the control. Conversely control, state or mixed constraints will involve control, state or both evaluated at the declared time (e.g., x(t) + u(t)). Other combinations should be detected as incorrect by the parser. The variable may be involved in any of the four previous constraints. Constraints involving the variable only are variable constraints, either linear or nonlinear. In the example below, there are

two linear boundary constraints,
one linear variable constraint,
one linear state constraint,
one (two-sided) nonlinear control constraint.

julia

@def begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x ∈ R², state
    u ∈ R, control
    x(0) == [-1, 0]
    x(tf) == [0, 0]
    ẋ(t) == [x₂(t), u(t)]
    tf ≥ 0
    x₂(t) ≤ 1
    0.1 ≤ u(t)^2 ≤ 1
    tf → min
end

Duplicate box constraints

If the same scalar component of the state, control or variable appears in several box constraints (linear range constraints), the effective bounds are the intersection of all declared bounds: the effective lower bound is the maximum of declared lower bounds, and the effective upper bound is the minimum of declared upper bounds. A warning is emitted, and an error is thrown if the resulting interval is empty. See Duplicate box constraints below.

Note

Symbols like <= or >= are also authorised:

julia

@def begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x ∈ R², state
    u ∈ R, control
    x(0) == [-1, 0]
    x(tf) == [0, 0]
    ẋ(t) == [x₂(t), u(t)]
    tf >= 0
    x₂(t) <= 1
    0.1 ≤ u(t)^2 <= 1
    tf → min
end

Warning

Write either u(t)^2 or (u^2)(t), not u^2(t) since in Julia the latter means u^(2t). Moreover, in the case of equalities or of one-sided inequalities, the control and / or the state must belong to the left-hand side. The following errors:

julia

julia> @def begin
           t ∈ [0, 2], time
           x ∈ R², state
           u ∈ R, control
           x(0) == [-1, 0]
           x(2) == [0, 0]
           ẋ(t) == [x₂(t), u(t)]
           1 ≤ x₂(t)
           -1 ≤ u(t) ≤ 1
       end
ParsingError → top-level scope, REPL[1]:2
│
│  
│  Line 7: 1 ≤ x₂(t)
│  the upper bound of a constraint must not depend on the state; write a functional constraint instead by moving the term to the constrained side (e.g. `x₂(0) - v == 0` rather than `x₂(0) == v`)
└─

Warning

A constraint bound must be constant — it may not depend on the variable, the state, the control or the time. The following is rejected, with a message naming the cause:

julia

julia> @def begin
           v ∈ R, variable
           t ∈ [0, 1], time
           x ∈ R², state
           u ∈ R, control
           -1 ≤ v ≤ 1
           x₁(0) == -1
           x₂(0) == v # wrong: the bound v depends on the variable
           x(1) == [0, 0]
           ẋ(t) == [x₂(t), u(t)]
           ∫(0.5u(t)^2) → min
       end
ParsingError → top-level scope, REPL[1]:2
│
│  
│  Line 7: x₂(0) == v
│  the lower bound of a constraint must not depend on the variable; write a functional constraint instead by moving the term to the constrained side (e.g. `x₂(0) - v == 0` rather than `x₂(0) == v`)
└─

Write instead a boundary constraint, moving the term to the constrained side — that side may involve the variable:

julia

@def begin
    v ∈ R, variable
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    -1 ≤ v ≤ 1
    x₁(0) == -1
    x₂(0) - v == 0 # OK: this side may involve the variable
    x(1) == [0, 0]
    ẋ(t) == [x₂(t), u(t)]
    ∫(0.5u(t)^2) → min
end

Duplicate box constraints

Box constraints are linear range constraints on a single scalar component of the state, control or variable — that is, constraints of the form lb ≤ x_i(t) ≤ ub, u_i(t) ≤ ub, v_i ≥ lb, etc. (nonlinear path constraints are not concerned by this section).

When the same scalar component is targeted by several box-constraint declarations, OptimalControl does not keep them as separate constraints. Instead, it merges them by taking the intersection of all declared bounds:

the effective lower bound is max of all declared lower bounds,
the effective upper bound is min of all declared upper bounds,
a single @warn is emitted per duplicated component, listing every contributing label,
all labels that declared the component are preserved as aliases (accessible via the aliases field of state_constraints_box, control_constraints_box and variable_constraints_box; see Inspect a problem),
if the intersection is empty (max(lbs) > min(ubs)), an IncorrectArgument exception is thrown.

For instance,

julia

@def begin
    t ∈ [0, 1], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    0 ≤ q(t) ≤ 2, (q_wide)
    1 ≤ q(t) ≤ 3, (q_tight)
    ẋ(t) == [v(t), u(t)]
    ∫(u(t)^2) → min
end

┌ Warning: Multiple bound declarations for state component 1 (labels: q_wide, q_tight). Intersection applied: effective lb = 1.0, effective ub = 2.0.
└ @ CTModels.Building /projects/ctb/gha-runner/julia_depot/packages/CTModels/tPIej/src/Building/build.jl:125

yields the effective constraint 1 ≤ q(t) ≤ 2, with aliases = [:q_wide, :q_tight] for that component, and a warning reporting both labels.

Conversely, the following declares an empty feasible set and raises an error at build time:

julia

julia> @def begin
           t ∈ [0, 1], time
           x = (q, v) ∈ R², state
           u ∈ R, control
           0 ≤ q(t) ≤ 1, (low)
           2 ≤ q(t) ≤ 3, (high) # max(lbs)=2 > min(ubs)=1 ⇒ IncorrectArgument
           ẋ(t) == [v(t), u(t)]
           ∫(u(t)^2) → min
       end
┌ Warning: Multiple bound declarations for state component 1 (labels: low, high). Intersection applied: effective lb = 2.0, effective ub = 1.0.
└ @ CTModels.Building /projects/ctb/gha-runner/julia_depot/packages/CTModels/tPIej/src/Building/build.jl:125
IncorrectArgument → _dedup_box_constraints!, build.jl:128
│
│  Empty feasible set for state component 1
│
│  Got       max(lbs)=2.0 > min(ubs)=1.0
│  Expected  max(lbs) ≤ min(ubs)
│
│  Context   _dedup_box_constraints! - infeasibility check
│  Hint      Check the declared bounds for labels low, high.
└─

Mayer cost

julia

:( $e1 → min )
:( $e1 → max )

Mayer costs are defined in a similar way to boundary conditions and follow the same rules. The symbol → is used to denote minimisation or maximisation, the latter being treated by minimising the opposite cost. (The symbol => can also be used.)

julia> @def begin
           tf ∈ R, variable
           t ∈ [0, tf], time
           x = (q, v) ∈ R², state
           u ∈ R, control
           tf ≥ 0
           -1 ≤ u(t) ≤ 1
           q(0) == 1
           v(0) == 2
           q(tf) == 0
           v(tf) == 0
           0 ≤ q(t) ≤ 5
          -2 ≤ v(t) ≤ 3
           ẋ(t) == [v(t), u(t)]
           tf → min
       end
Abstract definition:

    tf ∈ R, variable
    t ∈ [0, tf], time
    x = ((q, v) ∈ R², state)
    u ∈ R, control
    tf ≥ 0
    -1 ≤ u(t) ≤ 1
    q(0) == 1
    v(0) == 2
    q(tf) == 0
    v(tf) == 0
    0 ≤ q(t) ≤ 5
    -2 ≤ v(t) ≤ 3
    ẋ(t) == [v(t), u(t)]
    tf → min

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

    minimize  J(x, u, tf) = g(x(0), x(tf), tf)

    subject to

        ẋ(t) = f(x(t), u(t), tf), t in [0, tf] a.e.,

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

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

Lagrange cost

julia

:(       ∫($e1) → min )
:(     - ∫($e1) → min )
:( $e1 * ∫($e2) → min )
:(       ∫($e1) → max )
:(     - ∫($e1) → max )
:( $e1 * ∫($e2) → max )

Lagrange (integral) costs are defined used the symbol ∫, with parentheses. The keyword integral can also be used:

julia

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

julia

@def begin
    t ∈ [0, 1], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    0.5integral(q(t) + u(t)^2) → min
end

The integration range is implicitly equal to the time range, so the cost above is to be understood as

As for the dynamics, the parser will detect whether the integrand depends or not on time (autonomous / non-autonomous case).

Bolza cost

julia

:( $e1 +       ∫($e2)       → min )
:( $e1 + $e2 * ∫($e3)       → min )
:( $e1 -       ∫($e2)       → min )
:( $e1 - $e2 * ∫($e3)       → min )
:( $e1 +       ∫($e2)       → max )
:( $e1 + $e2 * ∫($e3)       → max )
:( $e1 -       ∫($e2)       → max )
:( $e1 - $e2 * ∫($e3)       → max )
:(             ∫($e2) + $e1 → min )
:(       $e2 * ∫($e3) + $e1 → min )
:(             ∫($e2) - $e1 → min )
:(       $e2 * ∫($e3) - $e1 → min )
:(             ∫($e2) + $e1 → max )
:(       $e2 * ∫($e3) + $e1 → max )
:(             ∫($e2) - $e1 → max )
:(       $e2 * ∫($e3) - $e1 → max )

Quite readily, Mayer and Lagrange costs can be combined into general Bolza costs. For instance as follows:

julia

@def begin
    p = (t0, tf) ∈ R², variable
    t ∈ [t0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    (tf - t0) + 0.5∫(c(t) * u(t)^2) → min
end

Warning

The expression must be the sum of two terms (plus, possibly, a scalar factor before the integral), not more, so mind the parentheses. For instance, the following errors:

julia

julia> @def begin
           p = (t0, tf) ∈ R², variable
           t ∈ [t0, tf], time
           x = (q, v) ∈ R², state
           u ∈ R, control
           ẋ(t) == [v(t), u(t)]
           (tf - t0) + q(tf) + 0.5∫(c(t) * u(t)^2) → min
       end
ParsingError → top-level scope, REPL[1]:2
│
│  
│  Line 6: (tf - t0) + q(tf) + 0.5 * ∫(c(t) * u(t) ^ 2) → min
│  bad objective declaration resulting in a Mayer term with trailing ∫
└─

The correct syntax is

julia

@def begin
    p = (t0, tf) ∈ R², variable
    t ∈ [t0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    ẋ(t) == [v(t), u(t)]
    ((tf - t0) + q(tf)) + 0.5∫(c(t) * u(t)^2) → min
end

Aliases

julia

:( $a = $e1 )

The single = symbol is used to define not a constraint but an alias, that is a purely syntactic replacement. There are some automatic aliases, e.g. x₁ and x1 for x[1] if x is the state (same for variable and control, for indices comprised between 1 and 9), and we have also seen that you can define your own aliases when declaring the variable, state and control. Arbitrary aliases can be further defined, as below (compare with previous examples in the dynamics section):

julia

@def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R, control
    F₀ = [x₂(t), 0]
    F₁ = [0, 1]
    ẋ(t) == F₀ + u(t) * F₁
    ∫(u(t)^2) → min
end

Warning

Such aliases do not define any additional function and are just replaced textually by the parser. In particular, they cannot be used outside the @def begin ... end block. Conversely, constants and functions used within the @def block must be defined outside and before this block.

Hint

You can rely on a trace mode for the macro @def to look at your code after expansions of the aliases using the @def ocp ... syntax and adding true after your begin ... end block:

julia

julia> @def damped_integrator begin
           tf ∈ R, variable
           t ∈ [0, tf], time
           x = (q, v) ∈ R², state
           u ∈ R, control
           q̇ = v(t)
           v̇ = u(t) - c(t)
           ẋ(t) == [q̇, v̇]
           ∫(u(t)^2) → min
       end true;
variable: tf, dim: 1
time: t, initial time: 0, final time: var"tf##11765"[1]
state: x, dim: 2
control: u, dim: 1
alias: q̇ = (x[2])(t)
alias: v̇ = (var"u##11766"[1])(t) - c(t)
dynamics: ∂(x)(t) == [(x[2])(t), (var"u##11766"[1])(t) - c(t)]
objective (Lagrange): ∫((var"u##11766"[1])(t) ^ 2) → min

Warning

The dynamics of an OCP is indeed a particular constraint, be careful to use == and not a single = that would try to define an alias:

julia

julia> double_integrator = @def begin
           tf ∈ R, variable
           t ∈ [0, tf], time
           x = (q, v) ∈ R², state
           u ∈ R, control
           q̇ = v
           v̇ = u
           ẋ(t) = [q̇, v̇]
       end
ParsingError → top-level scope, REPL[1]:2
│
│  
│  Line 7: ẋ(t) = begin
│          #= REPL[1]:9 =#
│          [q̇, v̇]
│      end
│  forbidden alias name: (∂(x))(t)
└─

Misc

Declarations (of variable - if any -, time, state and control - if any -) must be done first. Then, dynamics, constraints and cost can be introduced in an arbitrary order.
It is possible to provide numbers / labels (as in math equations) for the constraints to improve readability (this is mostly for future use, typically to retrieve the Lagrange multiplier associated with the discretisation of a given constraint):

julia

@def damped_integrator begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    tf ≥ 0, (1)
    q(0) == 2, (♡)
    q̇ = v(t)
    v̇ = u(t) - c(t)
    ẋ(t) == [q̇, v̇]
    x(t).^2  ≤ [1, 2], (state_con)
    ∫(u(t)^2) → min
end

Parsing errors should be explicit enough (with line number in the @def begin ... end block indicated).
Check tutorials and applications in the documentation for further use.

Known issues

Reverse over forward AD issues with ADNLP

Abstract syntax (@def) ​

Variable ​

Time ​

State ​

Control ​

Problems without a control ​

Dynamics ​

Dynamics (coordinatewise) ​

Constraints ​

Duplicate box constraints ​

Mayer cost ​

Lagrange cost ​

Bolza cost ​

Aliases ​

Misc ​

Known issues ​

See also ​

Abstract syntax (`@def`)

Variable

Time

State

Control

Problems without a control

Dynamics

Dynamics (coordinatewise)

Constraints

Duplicate box constraints

Mayer cost

Lagrange cost

Bolza cost

Aliases

Misc

Known issues

See also