OptimalControl.jl

Appearance

Plot ​

A Solution is drawn with plot and plot!, which extend Plots.jl. Use plot to create a new figure, and plot! to add to an existing one:

julia
plot(sol, description...; kw...)        # a new plot, which becomes the current one
plot!(sol, description...; kw...)       # adds to the current plot
plot!(plt, sol, description...; kw...)  # adds to the plot plt

The same calls draw the solution of a Flow (last section). The full signatures are under Reference at the end of the page.

Choosing a backend. This page uses Plots, the default. Makie.jl draws the same figures, with the same keywords. It is useful for an interactive window, or if you already draw with Makie: see Plot with Makie.

The table below lists the arguments of plot, and the section where each is explained:

SectionArguments
What gets drawnsize, state_style, costate_style, control_style, time_style, any Plots.jl attribute
Choosing what to drawdescription...: :state, :costate, :control, :path, :dual
Layoutlayout
The controlcontrol
Normalised timetime
Constraintsstate_bounds_style, control_bounds_style, path_style, path_bounds_style, dual_style

To overlay a second solution, see Adding to an existing plot. To build your own figure from the trajectories, see Custom plots.

Getting started ​

We take the energy-minimal double integrator, whose solution is u(t) = 6 − 12t (see Your first problem):

julia
using OptimalControl
using NLPModelsIpopt

t0 = 0
tf = 1
x0 = [-1, 0]
xf = [0, 0]

ocp = @def begin
    t ∈ [t0, tf], time
    x ∈ R², state
    u ∈ R, control
    x(t0) == x0
    x(tf) == xf
    ẋ(t) == [x₂(t), u(t)]
    ∫(0.5u(t)^2) → min
end

sol = solve(ocp; display=false)

plot draws a solution once Plots is loaded:

julia
using Plots
plot(sol)

Without Plots, the call raises an error that names the package to load:

julia
julia> plot(sol)
ERROR: ExtensionError
│
│  missing dependencies to plot solutions
│
│  Missing  Plots
│
│  Hint     Run: using Plots
└─

What gets drawn ​

With every group shown, the figure is a grid: the state components on the left, the costate components on the right, and the control at the bottom. The dashed vertical lines mark the initial and final times. Plots.jl attributes apply to the whole figure:

julia
plot(
    sol, :state, :costate, :control;
    size=(700, 450), legend=:bottomright, grid=false, linewidth=2,
)

state_style, costate_style and control_style set the attributes of one group, as a NamedTuple:

julia
plot(sol, :state, :costate, :control;
    state_style=(color=:blue,),
    costate_style=(color=:black, linestyle=:dash),
    control_style=(color=:red, linewidth=2),
)

time_style styles the vertical lines at the initial and final times. Every *_style also accepts :none, which hides that group:

julia
plot(sol, :state, :costate, :control;
    state_style=:none,
    costate_style=:none,
    control_style=(color=:red,),
    time_style=(color=:green,),
)

Three attributes are the exception: title, xlabel and ylabel are set by plot itself, from the names in the problem, and the values you pass are ignored without a warning (CTBase#569). To change them, edit the subplot afterwards, for instance xlabel!(plt[3], "s") (see Custom plots).

To look up a Plots.jl attribute and its aliases, use plotattr("linestyle") (or any other name) once Plots is loaded.

Choosing what to draw ​

The positional symbols select the groups to draw: :state, :costate, :control, and, for a problem with path constraints, :path and :dual.

julia
plot(sol, :state)    # only the state
plot(sol, :costate)  # only the costate
plot(sol, :control)  # only the control

They combine freely:

julia
plot(sol, :state, :control)

Layout ​

layout=:group puts each group (state, costate, control) in a single subplot, instead of one subplot per component:

julia
plot(sol; layout=:group)

layout=:split, the default, is the grid with one subplot per component used above.

The control ​

For a control with several components, control=:norm draws its Euclidean norm instead of its components, and control=:all draws both. Take a problem with a 2-D control:

julia
ocp_u = @def begin
    t ∈ [0, 1], time
    x ∈ R², state
    u ∈ R², control
    x(0) == [0, 0]
    x(1) == [1, 1]
    ẋ(t) == [u₁(t), x₁(t) + u₂(t)]
    0.5∫(u₁(t)^2 + u₂(t)^2) → min
end
sol_u = solve(ocp_u; display=false)

The maximum principle gives u = p, with p₂ constant and ṗ₁ = −p₂. The final conditions then give u₁(t) = (16 − 6t)/13 and u₂ = 6/13, so the norm decreases from about 1.31 to 0.90:

julia
plot(sol_u; control=:components, layout=:group, size=(800, 300))  # the default

julia
plot(sol_u; control=:norm, layout=:group, size=(800, 300))

julia
plot(sol_u; control=:all, layout=:group)

Normalised time ​

To compare solutions with different final times, time=:normalize (or :normalise) draws them against the normalised time s = (t − t₀)/(t_f − t₀) ∈ [0, 1]. Here the same linear-quadratic problem is solved for three final times:

julia
function lqr(tf)
    ocp = @def begin
        t ∈ [0, tf], time
        x ∈ R², state
        u ∈ R, control
        x(0) == [0, 1]
        ẋ(t) == [x₂(t), -x₁(t) + u(t)]
        ∫(0.5(x₁(t)^2 + x₂(t)^2 + u(t)^2)) → min
    end
    return ocp
end

tfs = [3, 5, 30]
solutions = [solve(lqr(tf); display=false) for tf in tfs]

plt = plot()
for (tf, sol) in zip(tfs, solutions)
    plot!(
        plt, sol, :state, :control;
        time=:normalize, label="tf = $tf",
    )
end

using Plots.PlotMeasures
px1 = plot(plt[1]; legend=false)  # x₁
px2 = plot(plt[2]; legend=true)   # x₂
pu = plot(plt[3]; legend=false)   # u
plot(
    px1, px2, pu;
    layout=(1, 3), size=(800, 300), leftmargin=5mm, bottommargin=5mm,
)

The longer the horizon, the longer the solution stays near the origin: this is the turnpike behaviour of linear-quadratic problems.

Constraints ​

A minimum-time transfer, with the control in [−1, 1], and a thrust limited by the speed: u(t) + v(t) ≤ 1. This mixed constraint involves the state and the control, so it is a path constraint, not a box.

julia
ocp_c = @def begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    tf ≥ 0
    -1 ≤ u(t) ≤ 1, (u_box)
    u(t) + v(t) ≤ 1, (thrust)
    x(0) == [-1, 0]
    q(tf) == 0
    v(tf) == 0
    ẋ(t) == [v(t), u(t)]
    tf → min
end
sol_c = solve(ocp_c; display=false)
plot(sol_c, :state, :costate, :control, :path, :dual)

The vehicle accelerates as hard as the thrust limit allows, u = 1 − v, so v(t) = 1 − e^{−t}, then brakes at u = −1 until it stops. The :path panel shows u + v, equal to its bound 1 on the first arc, and the :dual panel its multiplier, nonzero on that arc only. The box bounds of the control are drawn on the control panel.

Switching at t₁ and stopping at t_f = t₁ + v(t₁) requires q(t₁) + v(t₁)²/2 = 0, that is t₁ − v₁ + v₁²/2 = 1 with v₁ = 1 − e^{−t₁}. This gives t₁ ≈ 1.47 and t_f ≈ 2.24:

julia
final_time(sol_c)
2.242186969020949

The style keywords of these groups are path_style and dual_style, and, for the bounds, state_bounds_style, control_bounds_style and path_bounds_style:

julia
plot(sol_c, :state, :costate, :control, :path, :dual;
    control_bounds_style=(linestyle=:dash,),
    path_style=(color=:green,),
    path_bounds_style=(linestyle=:dash,),
    dual_style=(color=:red,),
    time_style=:none,
)

Box constraints are only drawn as bounds: :path and :dual show the path constraints, not the boxes. The multipliers of the boxes are read with dual (see Dual variables).

Adding to an existing plot ​

plot! overlays a second solution, with the same state, costate and control dimensions:

julia
ocp2 = @def begin
    t ∈ [t0, tf], time
    x ∈ R², state
    u ∈ R, control
    x(t0) == [-0.5, -0.5]
    x(tf) == xf
    ẋ(t) == [x₂(t), u(t)]
    ∫(0.5u(t)^2) → min
end
sol2 = solve(ocp2; display=false)

plt = plot(sol, :state, :costate, :control; label="sol1", size=(700, 500))
plot!(plt, sol2, :state, :costate, :control; label="sol2", linestyle=:dash)

Custom plots ​

state, costate and control return functions of time, which Plots draws directly. Here the absolute value of the control:

julia
t = time_grid(sol)
u = control(sol)
plot(t, abs ∘ u; label="|u|", xlabel="t")

The subplots of a plot(sol, ...) figure can also be reached directly. They come in the order of the groups requested, one per component: here x₁, x₂, p₁, p₂, then u.

julia
plt = plot(sol, :state, :costate, :control)
plot(plt[1])  # x₁

julia
plot(plt[5])  # u

A subplot accepts the usual Plots.jl calls, to annotate it rather than redraw it. For example, to mark two levels of the control:

julia
plot(plt[5])
hline!([-6, 6]; linestyle=:dash, color=:red, label="")

In Makie the same takes one more step, since a panel there is an Axis, not a subplot: see Plot with Makie.

Plotting a flow ​

The same plot call draws the solution of a Flow (see Flows for how to build one). Here the flow of the energy problem, with the control u = p₂ given by the maximum principle, from the exact initial costate p(0) = (12, 6):

julia
using OrdinaryDiffEqTsit5

f = Flow(ocp, (x, p) -> p[2])
p0 = [12, 6]
sol_flow = f((t0, tf), x0, p0)
plot(sol_flow)

The solution of a flow is kept at the steps of the integrator only, which are few on a smooth problem (CTFlows#435). Here there are 5 points, so the curves above are polygons:

julia
time_grid(sol_flow)
5-element Vector{Float64}:
 0.0
 0.006946033220330977
 0.0800652659463251
 0.6329007849970932
 1.0

For a finer curve, pass saveat when you build the flow, together with dense=false. The call itself only accepts variable, unsafe and variable_costate. Do not leave out dense=false: without it, the call crashes the Julia session (CTFlows#434). Use such a flow for trajectories only: its point call f(t0, x0, p0, tf) returns the value at the last saveat time before tf, not at tf (CTFlows#441).

julia
fine_grid = range(t0, tf, 100)
f2 = Flow(ocp, (x, p) -> p[2]; saveat=fine_grid, dense=false)
sol_flow2 = f2((t0, tf), x0, p0)
plot(sol_flow2)

Reference ​

RecipesBase.plot Method
julia
plot(
    sol::CTModels.Solutions.Solution,
    description::Symbol...;
    kwargs...
) -> Any

Plot the components of an optimal control CTModels.Solutions.Solution.

Generates a set of subplots showing the state, control, costate, path constraints and dual variables over time, depending on the problem and the given description.

Arguments

  • sol: the optimal control solution to visualise.

  • description: symbols selecting which groups to include; any of :state, :costate, :control, :path (path constraints), :dual (their multipliers). If none is given, a default set is used based on the problem.

Keyword arguments

  • layout::Symbol = :split: :split (one subplot per component) or :group (group each signal into a single subplot with a legend).

  • control::Symbol = :components: :components (a curve per control component), :norm (the Euclidean norm ‖u(t)‖) or :all (both).

  • time::Symbol = :default: :default (real time) or :normalize/:normalise ([0, 1]).

  • color: colour applied to every curve.

  • size: figure size; defaults to a heuristic based on the layout.

Style options

Each *_style keyword is a NamedTuple of plotting attributes, or :none to hide the group/decoration: state_style, costate_style, control_style, path_style, dual_style, time_style (initial/final time markers), and the bounds decorations state_bounds_style, control_bounds_style, path_bounds_style.

Returns

A Plots.Plot. All layout and rendering is delegated to CTBase.Plotting.

Example

julia
julia> plot(sol)
julia> plot(sol, :state, :control; layout=:group, control=:all)
julia> plot(sol; state_style=(color=:blue,), costate_style=:none)
RecipesBase.plot! Method
julia
plot!(
    sol::CTModels.Solutions.Solution,
    description::Symbol...;
    kwargs...
) -> Plots.Plot

Overlay the optimal control solution sol onto the current plot (Plots.current()).

RecipesBase.plot! Method
julia
plot!(
    p::Plots.Plot,
    sol::CTModels.Solutions.Solution,
    description::Symbol...;
    kwargs...
) -> Plots.Plot

Overlay the optimal control solution sol onto the existing plot p. Same behaviour and keyword arguments as plot (documented above); an empty p is filled as if by plot.

See also ​