Simulation

Sometimes you do not look for an optimum: you have a control system and a given control, open-loop or feedback, and you want the trajectory it produces, or the cost of that control for a problem.

using OptimalControl
using OrdinaryDiffEqTsit5
using Plots

The idea

A controlled vector field   and a control are enough to integrate a trajectory. The control is given as a function of time (open loop) or of the state (feedback, closed loop). No cost, no costate: the vector field is a ControlledVectorField, and the control a law OpenLoop or ClosedLoop.

Open loop

Take     with the constant control  . From  , the solution is   :

fc(x, u) = -x + u

f_ol = Flow(ControlledVectorField(fc), OpenLoop(t -> 1))
f_ol(0, 0.0, 1)   # 1 - exp(-1)

0.6321205587575749

An open-loop law is always a function of time, u(t), or u(t, v) with a variable. A function without argument, such as OpenLoop(() -> 1), is accepted when the law is built, but the call then fails with a MethodError (CTBase#570):

julia> Flow(ControlledVectorField(fc), OpenLoop(() -> 1))(0, 0.0, 1)
MethodError: no method matching (::Main.var"#5#6")(::Float64)
The function `#5` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  (::Main.var"#5#6")()
   @ Main REPL[1]:2

Closed loop

A feedback is a function of the state, u(x). With  , the system becomes   , so from  ,  :

f_cl = Flow(ControlledVectorField(fc), ClosedLoop(x -> x / 2))
f_cl(0, 1.0, 1)   # exp(-1/2)

0.6065306597356303

A law of the state and the costate, DynClosedLoop, is rejected here: the flow of a controlled vector field has no costate.

julia> Flow(ControlledVectorField(fc), DynClosedLoop((x, p) -> p))
PreconditionError → top-level scope, REPL[1]:2
│
│  Flow(fc, law) requires an OpenLoop or ClosedLoop control law
│
│  Reason   a DynClosedLoop law u(t,x,p,v) needs the costate p, which a state flow of a vector field does not have
│
│  Context  Flow(fc::ControlledVectorField, law::ControlLaw) — feedback dispatch
│  Hint     use OpenLoop(u) or ClosedLoop(u); for a DynClosedLoop law use Flow(h̃, law) or Flow(ocp, law)
└─

Evaluating a control on a problem

Built from an optimal control problem, the simulation also computes the cost of the control. On the energy-minimal double integrator, the optimal control is   , with cost (see Your first problem):

t0, tf = 0, 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

f_opt = Flow(ocp, OpenLoop(t -> 6 - 12t))
traj = f_opt((t0, tf), x0)
state(traj)(tf), objective(traj)

([-2.5854463843733e-15, 1.0329944957352045e-14], 6.000000000000015)

The simulation reaches the target with the optimal cost. Another control costs less, but misses the target: the constant   costs and ends at   .

traj_1 = Flow(ocp, OpenLoop(t -> 1))((t0, tf), x0)
state(traj_1)(tf), objective(traj_1)

([-0.5000000000000058, 0.9999999999999959], 0.4999999999999999)

The boundary conditions of the problem are not enforced: the simulation only integrates the dynamics and the cost from . A feedback works the same way, with Flow(ocp, ClosedLoop(x -> …)).

Reading and plotting the trajectory

The trajectory is read with the same functions as a Solution: state, control, time_grid, and objective when the flow comes from a problem.

state(traj)(0.5), control(traj)(0.5)

([-0.5000000000000058, 1.5000000000000016], 0.0)

plot(traj)

There is no costate: costate(traj) raises an error, and so does objective on a flow built without a problem.

julia> objective(f_ol((0, 1), 0.0))
PreconditionError → top-level scope, REPL[1]:2
│
│  this StateFlowTrajectory has no objective value
│
│  Reason   it was built without an optimal control problem (e.g. from Flow(fc, law) rather than Flow(ocp, law)), so there is no cost to evaluate
│
│  Context  StateFlowTrajectory — objective getter
│  Hint     build the flow from an OCP — Flow(ocp, law) — to obtain the objective, or use state(sol) and control(sol)
└─

What a flow returns

Every flow is called in two ways: f(t0, x0, tf) returns the final state (with p0, the final state and costate), and f((t0, tf), x0) returns the trajectory. The trajectory depends on how the flow was built:

Flow	Trajectory	Costate	Objective
Flow(ocp, law), with a law u(x, p) (From an OCP)	a Solution	✓	✓
Flow(ocp, OpenLoop(…)), Flow(ocp, ClosedLoop(…))	a trajectory of the state	✗	✓
Flow(ControlledVectorField(f), law), Flow(VectorField(f))	a trajectory of the state	✗	✗
Flow(Hamiltonian(H)) and the other Hamiltonian constructors	a trajectory of the state and costate	✓	✗

All of them are read with state, costate (when there is one), time_grid and plot. The types of the trajectories are internal: call these functions on them, without naming the types.

See also

• Solution object and Plot: the functions used above, in full.

• From an OCP: the indirect method, where the law comes from the maximum principle.

• Control-free problems: problems without a control at all.