Solution object

This page is about what you do once a problem has been solved: print the solution, plot it, read its trajectories and its objective, check that the solver converged, and read the Lagrange multipliers. It mirrors Inspect a problem: that page reads a model back, this one reads a solution back.

A solution to read

We take the energy-minimal double integrator, whose solution is known in closed form: u(t) = 6 − 12t, q(t) = −1 + 3t² − 2t³, v(t) = 6t − 6t², and the cost is J = 6 (see Your first problem).

using OptimalControl
using NLPModelsIpopt

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

sol = solve(ocp; display=false)

A solution prints as a summary:

sol

Solution  ✓ successful
  │  Objective : 6.000096001536038
  │  Boundary duals : [12.000192003072058, 6.0000960015360265, -12.000192003072058, 6.000096001536031]
  │
  │  Iterations : 1
  │  Status : first_order
  │  Message : Ipopt/generic
  └─ Constraints violation : 2.220446049250313e-16

and plots in one line, once a plotting package is loaded (see Plot):

using Plots
plot(sol)

The Solution type

solve returns an OptimalControl.Solution. The type is not exported: you read a solution through the functions below, never through its fields. They return functions of time where the fields store values on a grid, and they apply the conventions of the rest of the site, such as "1-D is a scalar". Every function on this page is listed with its signature in the API reference.

Trajectories

state, control and costate return functions of time:

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

([-0.8437454999279995, 1.124993999904], 3.0000480007680177, [12.00019200307206, 2.976047616761874])

The closed form gives x(0.25) = (−0.84375, 1.125), u(0.25) = 3 and p₁ = 12.

These functions interpolate: they can be called anywhere in the time horizon, not only at the grid points. time_grid(sol) returns the grid (times(sol) is the same vector), 251 points here since the default grid_size is 250 steps:

length(time_grid(sol)), 0.25 ∈ time_grid(sol)

(251, false)

1-D is a scalar here too: with a 1-D control, u(t) is a Number, not a length-1 vector, as everywhere else on the site (functional-API callbacks, abstract syntax):

typeof(u(0.25))

Float64

variable(sol) is a value, not a function of time: a Float64 for a 1-D variable, a vector otherwise, and an empty vector when the problem has no variable, as here:

variable(sol)

Float64[]

The costate is shifted by half a step

The costate of a direct solve is the multiplier of the discretised dynamics: with the default :midpoint scheme, there is one per step, and it approximates p at the middle of the step. So p(t) is close to the exact costate at t + h/2, with h the step. Here p₂(t) = 6 − 12t, h = 0.004, and p₂(0.25) is 6 − 12 × 0.252 = 2.976, not 3:

p(0.25)[2]

2.976047616761874

The difference vanishes as the grid is refined. With the Gauss–Legendre schemes, the shift is a full step: p(t) is close to the exact costate at t + h (CTDirect#640). For an accurate initial costate, for instance to start a shooting method, the multiplier of the initial condition is better: see Dual variables.

Function	Returns
state(sol)	x(t), a function of time
control(sol)	u(t), a function of time
costate(sol)	p(t), a function of time
variable(sol)	the optimisation variable: a number, a vector, or an empty vector
time_grid(sol), times(sol)	the discretisation grid, a vector

The time horizon

initial_time and final_time read the endpoints back. On a fixed-time problem they give what you wrote in the @def:

initial_time(sol), final_time(sol)

(0, 1)

When a time is an optimisation variable, its value is not in the model: the model only records which component of the variable holds it. The accessor reads that component of the solution's variable for you. The second problem of this page is a minimum-time transfer, with the control in [−1, 1] and the speed capped at 0.75. Its labels are used in Dual variables below.

ocp_free = @def begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    tf ≥ 0, (tf_pos)
    -1 ≤ u(t) ≤ 1, (u_box)
    v(t) ≤ 0.75, (v_max)
    x(0) == [-1, 0], (start)
    q(tf) == 0
    v(tf) == 0
    ẋ(t) == [v(t), u(t)]
    tf → min
end
sol_free = solve(ocp_free; display=false)

initial_time(sol_free), final_time(sol_free)

(0, 2.083333320623279)

The optimal control accelerates at u = 1 until the speed reaches 0.75 (t = 0.75), coasts at that speed, then brakes at u = −1 for 0.75 time units. Accelerating and braking each cover 0.75²/2 = 0.28125 of the unit distance, which leaves 0.4375 to coast, in 0.4375/0.75 = 7/12. So t_f = 0.75 + 7/12 + 0.75 = 25/12 ≈ 2.0833. The final time is also the variable, here a number:

variable(sol_free)

2.083333320623279

The objective

objective(sol)

6.000096001536038

The closed form gives J = 6. On the minimum-time problem, the objective is the final time.

Did it converge

successful says whether the solver reported success. status and message give its verdict, iterations and constraints_violation two figures to judge it by:

successful(sol), status(sol), message(sol)

(true, :first_order, "Ipopt/generic")

iterations(sol), constraints_violation(sol)

(1, 2.220446049250313e-16)

infos(sol) returns a Dict of anything else the solver reported (empty with Ipopt).

success is not successful

success is a Julia function, about the exit status of a process. Calling it on a solution raises an error that points to successful:

julia> success(sol)
PreconditionError → success, deprecated.jl:86
│
│  `success(sol)` is deprecated
│
│  Reason  this spelling was removed in v2.1.0-beta
│
│  Hint    use successful(sol)
└─

See Migrating from v2.0 for every renamed function.

What to do when a solve fails (a better initial guess, a finer grid, other tolerances) is the subject of Initial guess and Options.

Function	Returns
objective(sol)	the value of the cost
successful(sol)	true if the solver reported success
status(sol)	the solver's status, a Symbol such as :first_order
message(sol)	the solver's message, a String
iterations(sol)	the number of iterations
constraints_violation(sol)	the largest constraint violation
infos(sol)	a Dict of anything else the solver reported

Dual variables

The Lagrange multipliers of the constraints measure how much the cost would change if a constraint were relaxed. A constraint with a label is read with dual(sol, ocp, :label). For a boundary or a variable constraint, the multiplier is a number or a vector:

dual(sol_free, ocp_free, :tf_pos)

4.8000002513961475e-12

The constraint t_f ≥ 0 is inactive, so its multiplier is zero.

dual(sol_free, ocp_free, :start)

2-element Vector{Float64}:
 1.3298101221163192
 1.0011103794866778

The multiplier of the initial condition x(0) = x₀ is the initial costate p(t₀). Here the costate is p(t) = (4/3, 1 − 4t/3) on the first arc, so p(0) = (4/3, 1). This is a more accurate value than costate(sol)(0), which is shifted by half a step (see Trajectories).

For a path constraint, the multiplier is a function of time. It is nonzero only where the constraint is active: the control bound while u = ±1, the speed bound while coasting.

μ_u = dual(sol_free, ocp_free, :u_box)
μ_v = dual(sol_free, ocp_free, :v_max)
[(t, μ_u(t), μ_v(t)) for t in (0.3, 1.0, 1.8)]

3-element Vector{Tuple{Float64, Float64, Float64}}:
 (0.3, -0.004971887223438901, -2.222222211542058e-11)
 (1.0, -1.5597847542879807e-22, -0.011081750949940633)
 (1.8, 0.005247551166476096, -2.1428570733596554e-11)

plot(time_grid(sol_free), [μ_u, μ_v]; label=["μ_u" "μ_v"], xlabel="t")

Sign and scale of the multipliers

For a box constraint, on a state, control or variable component, the solver returns separate non-negative multipliers for the lower and the upper bound. dual returns their difference μ = μ_lb − μ_ub: μ > 0 where the lower bound is active, μ < 0 where the upper bound is active, μ = 0 where the constraint is inactive. Above, u = 1 on the first arc, so μ_u < 0, and u = −1 on the last one, so μ_u > 0.

The multipliers of the path constraints, those that are not boxes (such as v(t) + u(t) ≤ 1), differ in two ways (CTDirect#639):

• they keep the solver's sign, which is the opposite: positive where the upper bound is active;

• they are divided by the time step, whereas the multipliers of box constraints are not. So box multipliers are of the order of the step, and decrease when the grid is refined.

With :exa, the multipliers of the path and boundary constraints are not available yet, and come back as zeros (CTDirect#405).

The box multipliers can also be read per group, without a label: one function for the lower bounds and one for the upper bounds, for the state, the control and the variable. For the state and the control, they are functions of time, with one entry per component:

state_constraints_lb_dual(sol_free)(1.0), state_constraints_ub_dual(sol_free)(1.0)

([0.0, 0.0], [0.0, 0.011081750949940633])

control_constraints_lb_dual(sol_free)(1.8), control_constraints_ub_dual(sol_free)(1.8)

(0.005247551171476096, 4.999999749578668e-12)

For the variable, they are values:

variable_constraints_lb_dual(sol_free), variable_constraints_ub_dual(sol_free)

(4.8000002513961475e-12, 0.0)

The multipliers of the path constraints (boxes excluded) and of the boundary constraints, together: path_constraints_dual(sol) is a function of time, boundary_constraints_dual(sol) a vector, here with the 4 boundary conditions of ocp_free (x(0) = x₀ counts for two).

boundary_constraints_dual(sol_free)

4-element Vector{Float64}:
  1.3298101221163192
  1.0011103794866778
 -1.3298101221163192
  1.0009448000722094

Function	Returns
dual(sol, ocp, :label)	the multiplier of a labelled constraint: a number, a vector, or a function of time
state_constraints_lb_dual(sol), state_constraints_ub_dual(sol)	the state box multipliers, functions of time
control_constraints_lb_dual(sol), control_constraints_ub_dual(sol)	the control box multipliers, functions of time
variable_constraints_lb_dual(sol), variable_constraints_ub_dual(sol)	the variable box multipliers
path_constraints_dual(sol)	the path multipliers, boxes excluded, a function of time
boundary_constraints_dual(sol)	the boundary multipliers, a vector
dim_dual_state_constraints_box(sol), dim_dual_control_constraints_box(sol), dim_dual_variable_constraints_box(sol)	the number of box-constrained components
dim_path_constraints_nl(sol), dim_boundary_constraints_nl(sol)	the number of path (boxes excluded) and boundary constraints

Back to the model

model(sol) returns the problem the solution was computed from, so everything on Inspect a problem works on it:

model(sol) === ocp

true

Solutions from a flow

The Hamiltonian flow of a problem, built with a control law of the state and the costate, returns a Solution too (other flows return their own trajectories, see Simulation). Here the flow of the energy problem, with the control u = p₂ given by the maximum principle, from the initial costate p(0) = (12, 6):

using OrdinaryDiffEqTsit5
f = Flow(ocp, (x, p) -> p[2])
sol_flow = f((t0, tf), x0, [12, 6])
state(sol_flow)(tf), objective(sol_flow)

([-3.343153979267698e-16, 1.547256070246949e-14], 5.999999999999946)

The trajectories, objective, time_grid, model and plot work as above. What comes from a solver does not: a flow's solution carries no multipliers, so dual raises an error, and iterations and constraints_violation are NotProvided.

iterations(sol_flow), constraints_violation(sol_flow)

(NotProvided, NotProvided)

See also

• Inspect a problem: the model-side mirror of this page.

• Plot: draw everything read here.

• Save & load: write a solution to disk and read it back.

• Migrating from v2.0: every renamed function, success → successful included.