The optimal control solution object: structure and usage

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

• Plotting a solution: How to plot the optimal solution for your defined problem.
• Printing a solution: How to display a summary of your solution.

After covering these core functionalities, we'll delve into the structure of a solution. Since a solution is structured as a OptimalControl.Solution struct, we'll first explain how to access its underlying attributes. Following this, we'll shift our focus to the simple properties inherent to a solution.

Content

• Main functionalities
• Solution struct
• API Reference by Component

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

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

using NLPModelsIpopt
sol = solve(ocp)

▫ This is OptimalControl 1.3.3-beta, solving with: collocation → adnlp → ipopt (cpu)

  📦 Configuration:
   ├─ Discretizer: collocation
   ├─ Modeler: adnlp
   └─ Solver: ipopt

▫ This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:     1754
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      250

Total number of variables............................:      752
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      504
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  5.0000000e-03 1.10e+00 2.03e-14   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  6.0000960e+00 2.22e-16 1.78e-15 -11.0 6.08e+00    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   6.0000960015360247e+00    6.0000960015360247e+00
Dual infeasibility......:   1.7763568394002505e-15    1.7763568394002505e-15
Constraint violation....:   2.2204460492503131e-16    2.2204460492503131e-16
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             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total seconds in IPOPT                               = 2.978

EXIT: Optimal Solution Found.

To print sol, simply:

sol

• Solver:
  ✓ Successful  : true
  │  Status     : first_order
  │  Message    : Ipopt/generic
  │  Iterations : 1
  │  Objective  : 6.000096001536025
  └─ Constraints violation : 2.220446049250313e-16

• Boundary duals: [12.000192003072046, 6.000096001536023, -12.000192003072046, 6.000096001536035]

For complementary information, you can plot the solution:

using Plots
plot(sol)

Solution struct

The solution sol is a OptimalControl.Solution struct.

struct Solution{TimeGridModelType<:CTModels.OCP.AbstractTimeGridModel, TimesModelType<:CTModels.OCP.AbstractTimesModel, StateModelType<:CTModels.OCP.AbstractStateModel, ControlModelType<:CTModels.OCP.AbstractControlModel, VariableModelType<:CTModels.OCP.AbstractVariableModel, ModelType<:CTModels.OCP.AbstractModel, CostateModelType<:Function, ObjectiveValueType<:Real, DualModelType<:CTModels.OCP.AbstractDualModel, SolverInfosType<:CTModels.OCP.AbstractSolverInfos} <: CTModels.OCP.AbstractSolution

julia> using CTModels

julia> # Solutions are typically returned by solvers
julia> sol = solve(ocp, ...)  # Returns a Solution
julia> CTModels.objective(sol)

Each field can be accessed directly (sol.state, etc) but we recommend to use the sophisticated getters we provide: the state(sol::Solution) method does not return sol.state but a function of time that can be called at any time, not only on the grid time_grid.

0.25 ∈ time_grid(sol)

false

x = state(sol)
x(0.25)

2-element Vector{Float64}:
 -0.8437454999279986
  1.1249939999039986

You can also retrieve the original optimal control problem from the solution:

model(sol)  # returns the original OCP model

Abstract definition:

    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.

API Reference by Component

This section provides a comprehensive reference of all methods available for inspecting and querying optimal control solutions. Methods are organized by component for easy navigation.

Trajectories

The trajectory component provides access to the state, control, variable, and costate trajectories.

State trajectory

Get the state trajectory as a function of time:

x = state(sol)  # returns a function of time

#_wrap_scalar_and_deepcopy##2 (generic function with 1 method)

Evaluate the state at any time (not just grid points):

t = 0.25
x(t)  # returns state vector at t=0.25

2-element Vector{Float64}:
 -0.8437454999279986
  1.1249939999039986

The state function can be evaluated at any time within the problem horizon, even if it's not a discretization grid point:

0.25 ∈ time_grid(sol)  # false: not a grid point

false

x(0.25)  # still works: interpolated value

2-element Vector{Float64}:
 -0.8437454999279986
  1.1249939999039986

Control trajectory

Get the control trajectory as a function of time:

u = control(sol)  # returns a function of time

#_wrap_scalar_and_deepcopy##0 (generic function with 1 method)

u(t)  # returns control value at t

3.000048000768014

Variable values

Get the optimization variable values:

v = variable(sol)  # returns an empty vector if no variable

Float64[]

Costate trajectory

Get the costate (adjoint) trajectory as a function of time:

p = costate(sol)  # returns a function of time

#_wrap_scalar_and_deepcopy##2 (generic function with 1 method)

p(t)  # returns costate vector at t

2-element Vector{Float64}:
 12.000192003072046
  2.97604761676187

Time information

Get time-related information from the solution:

time_grid(sol)  # returns the discretization time grid

251-element Vector{Float64}:
 0.0
 0.004
 0.008
 0.012
 0.016
 0.02
 0.024
 0.028
 0.032
 0.036
 ⋮
 0.968
 0.972
 0.976
 0.98
 0.984
 0.988
 0.992
 0.996
 1.0

times(sol)  # returns the TimesModel struct containing time information

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

Summary table

Objective

The objective component provides access to the objective value.

Objective value

Get the optimal objective value:

objective(sol)  # returns the objective value

6.000096001536025

Summary table

Dual variables

The dual variables (Lagrange multipliers) provide sensitivity information about constraints.

To illustrate dual variables, we define a problem with various constraints:

ocp = @def begin
    tf ∈ R,             variable
    t ∈ [0, tf],        time
    x = (q, v) ∈ R²,    state
    u ∈ R,              control
    tf ≥ 0,             (eq_tf)
    -1 ≤ u(t) ≤ 1,      (eq_u)
    v(t) ≤ 0.75,        (eq_v)
    x(0)  == [-1, 0],   (eq_x0)
    q(tf) == 0
    v(tf) == 0
    ẋ(t) == [v(t), u(t)]
    tf → min
end
sol = solve(ocp; display=false)

Dual of labeled constraints

Get the dual variable for a specific labeled constraint:

dual(sol, ocp, :eq_tf)  # dual for variable constraint

4.800000251396152e-12

dual(sol, ocp, :eq_x0)  # dual for boundary constraint

2-element Vector{Float64}:
 1.3298101220985994
 1.0011103795098797

For path constraints, the dual is a function of time:

μ_u = dual(sol, ocp, :eq_u)
plot(time_grid(sol), μ_u)

μ_v = dual(sol, ocp, :eq_v)
plot(time_grid(sol), μ_v)

Box constraint duals

Get dual variables for box constraints on state, control, and variable:

state_constraints_lb_dual(sol)  # lower bound duals for state

#_wrap_scalar_and_deepcopy##0 (generic function with 1 method)

state_constraints_ub_dual(sol)  # upper bound duals for state

#_wrap_scalar_and_deepcopy##0 (generic function with 1 method)

control_constraints_lb_dual(sol)  # lower bound duals for control

#_wrap_scalar_and_deepcopy##0 (generic function with 1 method)

control_constraints_ub_dual(sol)  # upper bound duals for control

#_wrap_scalar_and_deepcopy##0 (generic function with 1 method)

variable_constraints_lb_dual(sol)  # lower bound duals for variable

1-element Vector{Float64}:
 4.800000251396152e-12

variable_constraints_ub_dual(sol)  # upper bound duals for variable

1-element Vector{Float64}:
 0.0

Nonlinear constraint duals

Get dual variables for nonlinear path and boundary constraints:

path_constraints_dual(sol)  # duals for nonlinear path constraints

#_wrap_scalar_and_deepcopy##2 (generic function with 1 method)

boundary_constraints_dual(sol)  # duals for nonlinear boundary constraints

4-element Vector{Float64}:
  1.3298101220985994
  1.0011103795098797
 -1.3298101220985994
  1.0009448000593428

Summary table

Solution metadata

The solution metadata provides information about the solver performance and status.

Solver status

Check if the solution was successful:

successful(sol)  # returns true if solver succeeded

true

Get the solver status symbol:

status(sol)  # returns solver status (e.g., :first_order)

:first_order

Get the solver message:

message(sol)  # returns solver message string

"Ipopt/generic"

Iteration count

Get the number of solver iterations:

iterations(sol)  # returns iteration count

17

Constraints violation

Get the maximum constraint violation:

constraints_violation(sol)  # returns max violation

7.771561172376096e-16

Additional solver information

Get additional solver-specific information:

infos(sol)  # returns dictionary of solver info

Dict{Symbol, Any}()

Summary table