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 (SOL). This includes essential operations like:

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

After covering these core functionalities, we'll delve into the structure of a SOL. Since a SOL is structured as a 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 SOL.

Content

• Main functionalities
• Solution 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

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

using NLPModelsIpopt
sol = solve(ocp)

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

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

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

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

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 3.29e-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.0000960015360212e+00    6.0000960015360212e+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.958

EXIT: Optimal Solution Found.

To print sol, simply:

sol

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

• Boundary duals: [12.000192003072044, 6.000096001536022, -12.000192003072044, 6.000096001536028]

For complementary information, you can plot the solution:

using Plots
plot(sol)

Solution struct

The solution sol is a Solution struct.

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

Each field can be accessed directly (ocp.times, 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

Attributes and properties

State, costate, control, variable and objective value

You can access the values of the state, costate, control and variable by eponymous functions. The returned values are functions of time for the state, costate and control and a scalar or a vector for the variable.

t = 0.25
x = state(sol)
p = costate(sol)
u = control(sol)

Since the state is of dimension 2, evaluating x(t) returns a vector:

x(t)

2-element Vector{Float64}:
 -0.8437454999279986
  1.1249939999039986

It is the same for the costate:

p(t)

2-element Vector{Float64}:
 12.000192003072046
  2.976047616761868

But the control is one-dimensional:

u(t)

2.9760476167618686

There is no variable, hence, an empty vector is returned:

v = variable(sol)

Float64[]

The objective value is accessed by:

objective(sol)

6.000096001536021

Infos from the solver

The problem ocp is solved via a direct method (see solve manual for details). The solver stores data in sol, including the success of the optimization, the iteration count, the time grid used for discretisation, and other specific details within the solver_infos field.

time_grid(sol)

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

constraints_violation(sol)

2.220446049250313e-16

infos(sol)

Dict{Symbol, Any}()

iterations(sol)

1

message(sol)

"Ipopt/generic"

status(sol)

:first_order

successful(sol)

true

Dual variables

You can retrieved dual variables (or Lagrange multipliers) associated to labelled constraint. To illustrate this, we define a problem with 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 variables corresponding to variable and boundary constraints are given as scalar or vectors.

dual(sol, ocp, :eq_tf)

4.800000251396152e-12

dual(sol, ocp, :eq_x0)

2-element Vector{Float64}:
 1.3298101220819512
 1.001110379518103

The other type of constraints are associated to dual variables given as functions of time.

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

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