The optimal control problem object: structure and usage

In this manual, we'll first recall the main functionalities you can use when working with an optimal control problem (OCP). This includes essential operations like:

• Solving an OCP: How to find the optimal solution for your defined problem.
• Computing flows from an OCP: Understanding the dynamics and trajectories derived from the optimal solution.
• Printing an OCP: How to display a summary of your problem's definition.

After covering these core functionalities, we'll delve into the structure of an OCP. Since an OCP is structured as a Model struct, we'll first explain how to access its underlying attributes, such as the problem's dynamics, costs, and constraints. Following this, we'll shift our focus to the simple properties inherent to an OCP, learning how to determine aspects like whether the problem:

• Is autonomous: Does its dynamics depend explicitly on time?
• Has a fixed or free initial/final time: Is the duration of the control problem predetermined or not?

Content

• Main functionalities
• Model struct
• Attributes and properties

Main functionalities

Let's define a basic optimal control problem.

using OptimalControl

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

To print it, simply:

ocp

Abstract defintion:

    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

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

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

    subject to

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

        ϕ₋ ≤ ϕ(x(0), x(1)) ≤ ϕ₊, 

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

Declarations (* required):
╭──────────┬────────┬────────┬──────────┬─────────────┬───────────┬────────────╮
│ variable │ times* │ state* │ control* │ constraints │ dynamics* │ objective* │
├──────────┼────────┼────────┼──────────┼─────────────┼───────────┼────────────┤
│    X     │   V    │   V    │    V     │      V      │     V     │     V      │
╰──────────┴────────┴────────┴──────────┴─────────────┴───────────┴────────────╯

We can now solve the problem (for more details, visit the solve manual):

using NLPModelsIpopt
solve(ocp)

▫ This is OptimalControl version v1.1.0 running with: direct, adnlp, ipopt.

▫ The optimal control problem is solved with CTDirect version v0.15.1.

   ┌─ The NLP is modelled with ADNLPModels and solved with NLPModelsIpopt.
   │
   ├─ Number of time steps⋅: 250
   └─ Discretisation scheme: trapeze

▫ This is Ipopt version 3.14.17, running with linear solver MUMPS 5.8.0.

Number of nonzeros in equality constraint Jacobian...:     3005
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      251

Total number of variables............................:     1004
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      755
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e-01 1.10e+00 2.73e-14   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -5.0000000e-03 7.36e-02 2.66e-15 -11.0 6.08e+00    -  1.00e+00 1.00e+00h  1
   2  6.0003829e+00 1.78e-15 1.78e-15 -11.0 6.01e+00    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 2

                                   (scaled)                 (unscaled)
Objective...............:   6.0003828724303263e+00    6.0003828724303263e+00
Dual infeasibility......:   1.7763568394002505e-15    1.7763568394002505e-15
Constraint violation....:   1.7763568394002505e-15    1.7763568394002505e-15
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.7763568394002505e-15    1.7763568394002505e-15


Number of objective function evaluations             = 3
Number of objective gradient evaluations             = 3
Number of equality constraint evaluations            = 3
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 3
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 2
Total seconds in IPOPT                               = 1.867

EXIT: Optimal Solution Found.

You can also compute flows (for more details, see the flow manual) from the optimal control problem, providing a control law in feedback form. The pseudo-Hamiltonian of this problem is

\[ H(x, p, u) = p_q\, q + p_v\, v + p^0 u^2 /2,\]

where $p^0 = -1$ since we are in the normal case. From the Pontryagin maximum principle, the maximising control is given in feedback form by

\[u(x, p) = p_v\]

since $\partial^2_{uu} H = p^0 = - 1 < 0$.

u = (x, p) -> p[2]          # control law in feedback form

using OrdinaryDiffEq        # needed to import numerical integrators
f = Flow(ocp, u)            # compute the Hamiltonian flow function

p0 = [12, 6]                # initial covector solution
xf, pf = f(t0, x0, p0, tf)  # flow from (x0, p0) at time t0 to tf
xf                          # should be (0, 0)

2-element Vector{Float64}:
 2.134348517083835e-15
 2.3650010008626292e-14

Model struct

The optimal control problem ocp is a Model struct.

struct Model{TD<:CTModels.TimeDependence, TimesModelType<:CTModels.AbstractTimesModel, StateModelType<:CTModels.AbstractStateModel, ControlModelType<:CTModels.AbstractControlModel, VariableModelType<:CTModels.AbstractVariableModel, DynamicsModelType<:Function, ObjectiveModelType<:CTModels.AbstractObjectiveModel, ConstraintsModelType<:CTModels.AbstractConstraintsModel, BuildExaModelType<:Union{Nothing, Function}} <: CTModels.AbstractModel

Each field can be accessed directly (ocp.times, etc) or by a getter:

• times
• state
• control
• variable
• dynamics
• objective
• constraints
• definition
• get_build_examodel

For instance, we can retrieve the times and definition values.

times(ocp)

CTModels.TimesModel{CTModels.FixedTimeModel{Int64}, CTModels.FixedTimeModel{Int64}}(CTModels.FixedTimeModel{Int64}(0, "0"), CTModels.FixedTimeModel{Int64}(1, "1"), "t")

definition(ocp)

quote
    #= manual-model.md:41 =#
    (t ∈ [t0, tf], time)
    #= manual-model.md:42 =#
    x = ((q, v) ∈ R², state)
    #= manual-model.md:43 =#
    (u ∈ R, control)
    #= manual-model.md:44 =#
    x(t0) == x0
    #= manual-model.md:45 =#
    x(tf) == [0, 0]
    #= manual-model.md:46 =#
    ẋ(t) == [v(t), u(t)]
    #= manual-model.md:47 =#
    0.5 * ∫(u(t) ^ 2) → min
end

Attributes and properties

Numerous attributes can be retrieved. To illustrate this, a more complex optimal control problem is defined.

ocp = @def begin
    v = (w, tf) ∈ R²,   variable
    s ∈ [0, tf],        time
    q = (x, y) ∈ R²,    state
    u ∈ R,              control
    0 ≤ tf ≤ 2,         (1)
    u(s) ≥ 0,           (cons_u)
    x(s) + u(s) ≤ 10,   (cons_mixed)
    w == 0
    x(0) == -1
    y(0) - tf == 0,     (cons_bound)
    q(tf) == [0, 0]
    q̇(s) == [y(s)+w, u(s)]
    0.5∫( u(s)^2 ) → min
end

Control, state and variable

You can access the name of the control, state, and variable, along with the names of their components and their dimensions.

using DataFrames
data = DataFrame(
    Data=Vector{Symbol}(),
    Name=Vector{String}(),
    Components=Vector{Vector{String}}(),
    Dimension=Vector{Int}(),
)

# control
push!(data,(
    :control,
    control_name(ocp),
    control_components(ocp),
    control_dimension(ocp),
))

# state
push!(data,(
    :state,
    state_name(ocp),
    state_components(ocp),
    state_dimension(ocp),
))

# variable
push!(data,(
    :variable,
    variable_name(ocp),
    variable_components(ocp),
    variable_dimension(ocp),
))

Constraints

You can retrieve labelled constraints with the constraint function. The constraint(ocp, label) method returns a tuple of the form (type, f, lb, ub). The signature of the function f depends on the symbol type. For :boundary and :variable constraints, the signature is f(x0, xf, v) where x0 is the initial state, xf the final state and v the variable. For other constraints, the signature is f(t, x, u, v). Here, t represents time, x the state, u the control, and v the variable.

(type, f, lb, ub) = constraint(ocp, :eq1)
println("type: ", type)
x0 = [0, 1]
xf = [2, 3]
v  = [1, 4]
println("val: ", f(x0, xf, v))
println("lb: ", lb)
println("ub: ", ub)

type: variable
val: 4
lb: 0.0
ub: 2.0

(type, f, lb, ub) = constraint(ocp, :cons_bound)
println("type: ", type)
println("val: ", f(x0, xf, v))
println("lb: ", lb)
println("ub: ", ub)

type: boundary
val: -3.0
lb: 0.0
ub: 0.0

(type, f, lb, ub) = constraint(ocp, :cons_u)
println("type: ", type)
t = 0
x = [1, 2]
u = 3
println("val: ", f(t, x, u, v))
println("lb: ", lb)
println("ub: ", ub)

type: control
val: 3
lb: 0.0
ub: Inf

(type, f, lb, ub) = constraint(ocp, :cons_mixed)
println("type: ", type)
println("val: ", f(t, x, u, v))
println("lb: ", lb)
println("ub: ", ub)

type: path
val: 4.0
lb: -Inf
ub: 10.0

Dynamics

The dynamics stored in ocp are an in-place function (the first argument is mutated upon call) of the form f!(dx, t, x, u, v). Here, t represents time, x the state, u the control, and v the variable, with dx being the output value.

f! = dynamics(ocp)
t = 0
x = [0., 1]
u = 2
v = [1, 4]
dx = similar(x)
f!(dx, t, x, u, v)
dx

2-element Vector{Float64}:
 2.0
 2.0

Criterion and objective

The criterion can be :min or :max.

criterion(ocp)

:min

The objective function is either in Mayer, Lagrange or Bolza form.

• Mayer:

\[g(x(t_0), x(t_f), v) \to \min\]

• Lagrange:

\[\int_{t_0}^{t_f} f^0(t, x(t), u(t), v)\, \mathrm{d}t \to \min\]

• Bolza:

\[g(x(t_0), x(t_f), v) + \int_{t_0}^{t_f} f^0(t, x(t), u(t), v)\, \mathrm{d}t \to \min\]

The objective of problem ocp is 0.5∫( u(t)^2 ) → min, hence, in Lagrange form. The signature of the Mayer part of the objective is g(x0, xf, v) but in our case, the method mayer will return an error.

julia> g = mayer(ocp)ERROR: UnauthorizedCall: This ocp has no Mayer objective.

The signature of the Lagrange part of the objective is f⁰(t, x, u, v).

f⁰ = lagrange(ocp)
f⁰(t, x, u, v)

2.0

To avoid having to capture exceptions, you can check the form of the objective:

println("Mayer: ", has_mayer_cost(ocp))
println("Lagrange: ", has_lagrange_cost(ocp))

Mayer: false
Lagrange: true

Times

The time variable is not named t but s in ocp.

time_name(ocp)

"s"

The initial time is 0.

initial_time(ocp)

0

Since the initial time has the value 0, its name is string(0).

initial_time_name(ocp)

"0"

In contrast, the final time is tf, since in ocp we have s ∈ [0, tf].

final_time_name(ocp)

"tf"

To get the value of the final time, since it is part of the variable v = (w, tf) of ocp, we need to provide a variable to the function final_time.

v = [1, 2]
tf = final_time(ocp, v)

2

julia> final_time(ocp)ERROR: UnauthorizedCall: You cannot get the final time with this function.

To check whether the initial or final time is fixed or free (i.e., part of the variable), you can use the following functions:

println("Fixed initial time: ", has_fixed_initial_time(ocp))
println("Fixed final time: ", has_fixed_final_time(ocp))

Fixed initial time: true
Fixed final time: false

Or, similarly:

println("Free initial time: ", has_free_initial_time(ocp))
println("Free final time: ", has_free_final_time(ocp))

Free initial time: false
Free final time: true

Time dependence

Optimal control problems can be autonomous or non-autonomous. In an autonomous problem, neither the dynamics nor the Lagrange cost explicitly depends on the time variable.

The following problem is autonomous.

ocp = @def begin
    t ∈ [ 0, 1 ], time
    x ∈ R, state
    u ∈ R, control
    ẋ(t)  == u(t)                       # no explicit dependence on t
    x(1) + 0.5∫( u(t)^2 ) → min         # no explicit dependence on t
end
is_autonomous(ocp)

true

The following problem is non-autonomous since the dynamics depends on t.

ocp = @def begin
    t ∈ [ 0, 1 ], time
    x ∈ R, state
    u ∈ R, control
    ẋ(t)  == u(t) + t                   # explicit dependence on t
    x(1) + 0.5∫( u(t)^2 ) → min
end
is_autonomous(ocp)

false

Finally, this last problem is non-autonomous because the Lagrange part of the cost depends on t.

ocp = @def begin
    t ∈ [ 0, 1 ], time
    x ∈ R, state
    u ∈ R, control
    ẋ(t)  == u(t)
    x(1) + 0.5∫( t + u(t)^2 ) → min     # explicit dependence on t
end
is_autonomous(ocp)

false