Solution

In this tutorial, we explain how to handle a solution. First, let's define a model

using Filtration

vlim = (0,50); mlim = (0,50)
ap = 0; bp = 1; cp = 0.16; dp = 3; Qp = 1
ar = 0.5; cr = 0.64; dr = 12; Qr = 2
p = [ap, bp, cp, dp, Qp, ar, cr, dr, Qr]

f₁(m,p) = -p[1]*m + p[2]
f₂(m,p) = -p[6]*m
g₁(m,p) = p[5]
g₂(m,p) = -p[9]
e₁(m,p) = p[3]*m + p[4]
e₂(m,p) = p[7]*m + p[8]

model = Model(e₁, e₂, f₁, f₂, g₁, g₂, p, :min; initial_guess = 20.)
nothing; # hide

Different ways to get a solution

Given initial and final conditions, the function solve provide the solution.

v0 = 0; vf = 50; m0₁ = 15
sol₁ = solve(model, v0, vf, m0₁)
plt = plot(sol₁, label = "sol₁ ")

Another way to obtain a solution is to extract it from a synthesis. This can be achieved thanks to the function solution.

synthesis = Synthesis(model, (v0, vf), (0, 50))
m0₂ = 30.
sol₂ = solution(synthesis, v0, vf , m0₂)
plot!(plt, sol₂, color = :red, label = "sol₂ ")

A final way to obtain a solution is to simulate the model with a used-defined control. In this case, there is no garentee that it is a solution of the defined optimal control problem.

However, to do this, we first need to create a control $u$ and to simulate it on the model.

u = Control([t -> 1, t -> sin(0.3*t)], [0, 25, 50], 2, 1)
x0 = [0, m0₂, 0]
sol₃ = simulate(model, u, x0)
plot!(plt, sol₃, color = :green, label = "sol₃ ")

State, control and extremal

From a Solution, one can extract a State and a Control with the state and the control function. These elements can be called to obtain the value of the state or control at a given time, and can be also plotted.

u = control(sol₁)

println("u(1) = ", u(10))
plot(u, size = (600, 200))

x = state(sol₁)
println("x(1) = ", x(10))
plot(x)

An extremal can be constructed from a (real) Solution of the problem, computed by the solution function on a Synthesis or the solve function on a Model, i.e. from sol₁ or sol₂ on the previous example. An Extremal must not be created from a Solution created from a used-defined Control, like sol₃. This can be done thanks to the construct_extremal function.

z = construct_extremal(model, sol₁)

u = control(z)
println("u(1) = ", u(1))
plot(u, size = (600, 200))

x = state(z)
println("x(1) = ", x(1))
plot(x)

λ = costate(z)
println("λ(1) = ", λ(1))
plot(λ)

The main advantage to get an extremal is to compute Hamiltonians along a trajectory. For instance, $t \mapsto H_1(x(t), p(t))$ corresponds to the switching function for the optimal control. Hamiltonians $H$, $H_0$, $H_1$ and even $H_{01}$ are accessible from the model and can be plotted as following.

plot( t -> model.H(x(t), λ(t), u(t), p), 0, 50, label = "H", ylim = (-1, 1))
plot!(t -> model.H₀(x(t), λ(t), p), 0, 50, label = "H₀")
plot!(t -> model.H₁(x(t), λ(t), p), 0, 50, label = "H₁")
plot!(t -> model.H₀₁(x(t), λ(t), p), 0, 50, label = "H₀₁")