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
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.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.865
EXIT: Optimal Solution Found.
You can export (or save) the solution in a Julia .jld2
data file and reload it later, and also export a discretised version of the solution in a more portable JSON format. Note that the optimal control problem is needed when loading a solution.
See the two functions:
To print sol
, simply:
sol
Optimal Control Solution
────────────────────────
• Successful : true
• Status : first_order
• Message : generic
• Iterations : 2
• Objective : 6.000382872430326
• Constraint violation: 1.7763568394002505e-15
Time
────
• Name : t
• Grid : [0.0, 0.004, 0.008, 0.012, 0.016, 0.02, 0.024, 0.028, 0.032, 0.036 … 0.964, 0.968, 0.972, 0.976, 0.98, 0.984, 0.988, 0.992, 0.996, 1.0]
• Grid length : 251
State
─────
• Name : x
• Dimension : 2
• Components : q, v
Control
───────
• Name : u
• Dimension : 1
• Components : u
Duals
─────
• Boundary constraints dual: [12.00076574486066, 6.000382872430307, -12.00076574486066, 6.000382872430306]
For complementary information, you can plot the solution:
using Plots
plot(sol)
For more details about plotting a solution, visit the plot manual.
Solution struct
The solution sol
is a Solution
struct.
CTModels.Solution
— Typestruct 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
Fields
time_grid::CTModels.AbstractTimeGridModel
times::CTModels.AbstractTimesModel
state::CTModels.AbstractStateModel
control::CTModels.AbstractControlModel
variable::CTModels.AbstractVariableModel
costate::Function
objective::Real
dual::CTModels.AbstractDualModel
solver_infos::CTModels.AbstractSolverInfos
model::CTModels.AbstractModel
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.8437499339957885
1.1249997839862171
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.8437499339957885
1.1249997839862171
It is the same for the costate:
p(t)
2-element Vector{Float64}:
12.00076574486066
2.976189904725448
But the control is one-dimensional:
u(t)
3.00019143621517
There is no variable, hence, an empty vector is returned:
v = variable(sol)
Float64[]
The objective value is accessed by:
objective(sol)
6.000382872430326
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)
1.7763568394002505e-15
infos(sol)
Dict{Symbol, Any}()
iterations(sol)
2
message(sol)
"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)
6.206046225282482e-11
dual(sol, ocp, :eq_x0)
2-element Vector{Float64}:
1.3333716227790058
1.0000505526995942
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)