From an OCP

This is the main path of the indirect method. You have worked out, from the maximum principle, the control law that maximises the pseudo-Hamiltonian, and you want the Hamiltonian flow it defines. Flow(ocp, law) reads the dynamics and the cost from the problem, so you only write the law. The notation is the one of Notation and conventions.

using OptimalControl
using OrdinaryDiffEqTsit5

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

The idea

For this problem, the pseudo-Hamiltonian is    , with    in the normal case. Since   , is maximised where    : the control law is  . Flow integrates the Hamiltonian system of  , from an initial state and costate.

The simplest form

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

([-4.158277980903794e-16, 5.038973481527467e-15], [12.0, -5.999999999999998])

The maximum principle gives   and   , so    for the initial costate  : the flow reaches the target   with   .

Passing a typed law

Flow(ocp, law) wraps the function in a DynClosedLoop law: a feedback on the state and the costate. Writing it explicitly gives the same flow:

f_typed = Flow(ocp, DynClosedLoop((x, p) -> p[2]))
xf_typed, _ = f_typed(t0, x0, p0, tf)
isapprox(xf_typed, xf; atol=1e-12)

true

The two other kinds of law, ClosedLoop(x -> …) and OpenLoop(t -> …), do not involve the costate: they build a flow of the state alone, to simulate a given control (see Simulation).

Non-autonomous problems

When the dynamics depend on , the control law does too: u(t, x, p).

t0b = 0
tfb = π / 4
x0b = 0
xfb = tan(π / 4) - 2log(√2 / 2)

ocp_na = @def begin
    t ∈ [t0b, tfb], time
    x ∈ R, state
    u ∈ R, control
    x(t0b) == x0b
    x(tfb) == xfb
    ẋ(t) == u(t) * (1 + tan(t))
    0.5∫(u(t)^2) → min
end

fb = Flow(ocp_na, (t, x, p) -> p * (1 + tan(t)))
xf_na, pf_na = fb(t0b, x0b, 1, tfb)
xf_na - xfb

-6.836773369656157e-10

Here       gives    , and is constant since does not depend on . Then       , so   reaches the target.

Problems with a variable

With an optimisation variable , the control law takes it as a third argument: u(x, p, v), or u(t, x, p, v) for a non-autonomous problem. Take a problem whose variable is the final time:

ocp_v = @def begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x ∈ R, state
    u ∈ R, control
    tf ≥ 0
    x(0) == 0
    x(tf) == 1
    ẋ(t) == tf * u(t)
    tf + 0.5∫(u(t)^2) → min
end

fc = Flow(ocp_v, (x, p, tf) -> tf * p)

For a given , the cheapest control reaching   is the constant  , which costs  . This is minimal for  . The control law    then gives the constant costate   .

The value of the variable is always passed with the keyword variable, even when it is the final time:

tf_val = (3 / 2)^(1 / 4)
p0c = 2tf_val / 3
xf_v, pf_v = fc(0, 0, p0c, tf_val; variable=tf_val)
xf_v

1.0000000000000004

Without it, the call raises an error:

julia> fc(0, 0, p0c, tf_val)
PreconditionError → top-level scope, REPL[1]:2
│
│  variable not provided for a NonFixed flow
│
│  Reason   flow depends on an extra variable parameter but none was given
│
│  Context  call — NonFixed flow with missing variable
│  Hint     Pass `variable=v` when calling the flow
└─

and so does passing variable to the flow of a problem without a variable:

julia> f(t0, x0, p0, tf; variable=1)
PreconditionError → top-level scope, REPL[1]:2
│
│  variable provided for a Fixed flow
│
│  Reason   flow does not depend on any variable parameter
│
│  Context  call — Fixed flow with unexpected variable
│  Hint     Remove the `variable` keyword argument when calling this flow
└─

Two other formulations of a free final time

Since the variable is the final time, the change of time    gives a problem on  , where only appears as a parameter:   , and the cost becomes   . The control law is unchanged:

ocp_s = @def begin
    tf ∈ R, variable
    s ∈ [0, 1], time
    x ∈ R, state
    u ∈ R, control
    tf ≥ 0
    x(0) == 0
    x(1) == 1
    ẋ(s) == tf^2 * u(s)
    tf + (0.5 * tf) * ∫(u(s)^2) → min
end

fs = Flow(ocp_s, (x, p, tf) -> tf * p)
xf_s, _ = fs(0, 0, p0c, 1; variable=tf_val)
xf_s

1.0000000000000002

Or becomes a state, with zero dynamics, and the problem has no variable any more. The state is  , its costate , and the control law reads    from them:

ocp_y = @def begin
    s ∈ [0, 1], time
    y = (x, tf) ∈ R², state
    u ∈ R, control
    x(0) == 0
    x(1) == 1
    ẏ(s) == [tf(s)^2 * u(s), 0]
    tf(1) + 0.5∫(tf(s) * u(s)^2) → min
end

fy = Flow(ocp_y, (y, q) -> y[2] * q[1])
yf, qf = fy(0, [0, tf_val], [p0c, 0], 1)

([1.0000000000000002, 1.1066819197003217], [0.7377879464668812, -1.0000000000000004])

The initial costate of is , since is free, and its final costate is : the transversality condition of the Mayer term , with   .

The constant 0 in the dynamics is fine for a flow. To solve this problem with the direct method and the default modeler, write 0 * u(s) instead (see Dynamics).

The costate of the variable

With variable_costate=true, the point call also integrates the costate of the variable,    from  , and returns three values, (xf, pf, pvf). This is what the transversality condition on the variable needs (see Notation and conventions).

Take a harmonic oscillator whose pulsation is optimised, with the cost  , and solve it with the direct method first:

using NLPModelsIpopt

ocp_ω = @def begin
    ω ∈ R, variable
    t ∈ [0, 1], time
    x = (q, v) ∈ R², state
    u ∈ R, control
    x(0) == [1, 0], (start)
    q(1) == 0
    ẋ(t) == [v(t), -ω^2 * q(t) + u(t)]
    ω^2 + 0.5∫(u(t)^2) → min
end

sol_ω = solve(ocp_ω; grid_size=500, display=false)
ω_opt = variable(sol_ω)
p0_ω = dual(sol_ω, ocp_ω, :start)   # the initial costate (see the Solution page)
ω_opt, p0_ω

(0.653763301228066, [-2.059403829775405, -2.413460259562427])

Then integrate the flow from this initial costate, with the costate of :

fω = Flow(ocp_ω, (x, p, ω) -> p[2])
xf_ω, pf_ω, pω = fω(0, [1, 0], p0_ω, 1; variable=ω_opt, variable_costate=true)
pω, -2ω_opt

(-1.307528321033575, -1.307526602456132)

With the Mayer term and   , transversality requires   : the direct solution satisfies it. The flow also reaches  , and   since is free.

variable_costate=true works on the point call only. On a trajectory call fω((0, 1), …), it raises a MethodError that does not say so yet (CTFlows#437).

Control-free problems

Flow(ocp), with no law, builds the flow of a problem without a control: see Control-free problems. On a problem with a control, it raises an error that asks for a law:

julia> Flow(ocp)
PreconditionError → top-level scope, REPL[1]:2
│
│  Flow from a with-control OCP is not supported
│
│  Reason   this path builds the flow from the OCP structure assuming no control (ẋ = f(t,x,∅,v)); a problem with a control input would require a control law u(t,x,p) from the maximisation of the pseudo-Hamiltonian
│
│  Context  Flow(ocp::CTModels.Models.Model) — control-dependence dispatch
│  Hint     use Flow(ocp, law) to pass a control law, or use a control-free OCP
└─

Total or partial Hamiltonian

The keyword hamiltonian_type chooses how the flow is built from and the law :

• :total, the default, integrates the Hamiltonian system of  , differentiating through the law:    and    ;

• :partial takes the partial derivatives of with the control frozen at the value of the law:  ,   , at  .

The two differ by the terms in . These vanish when   along the law, as for a control that maximises inside the control set, and also where the law is locally constant, as for bang–bang or saturated controls. So for an optimal control law the two give the same flow. For another law, they do not, and the state of :partial is the one the law actually drives,  : use :partial to simulate a given feedback.

Take the law   , which does not maximise :

law(x, p) = p[2] + 1

xf_total, _ = Flow(ocp, law; hamiltonian_type=:total)(t0, x0, p0, tf)
xf_partial, _ = Flow(ocp, law; hamiltonian_type=:partial)(t0, x0, p0, tf)
xf_total, xf_partial

([6.943952265347772e-16, 5.038973481527467e-15], [0.4999999999999938, 1.0000000000000056])

With :partial, the state is driven by     , which overshoots the target. With :total, the flow is the one of         . It differs from the optimal    by a constant, so it has the same Hamiltonian vector field, and it reaches the target as the optimal law does. This is a property of this example, not of :total.

What comes back

The point call f(t0, x0, p0, tf) returns (xf, pf), or (xf, pf, pvf) with variable_costate=true. The trajectory call f((t0, tf), x0, p0) returns a Solution, as solve does: state, costate, control, objective and plot work on it, but it has no multipliers (see Solutions from a flow).

sol = f((t0, tf), x0, p0)
objective(sol)

5.999999999999946

This is specific to flows built from a problem with a law of the state and the costate. The other flows, such as Flow(ocp, OpenLoop(…)) or Flow(Hamiltonian(H)), return their own trajectories: see What a flow returns.

See also

• Accessors: the Hamiltonian, its vector field and the law, read back from a flow.

• Shooting: turn a flow into equations for the initial costate.

• Control-free problems: Flow(ocp) in full.