OptimalControl.jl

Appearance

Control and variable together ​

The same two systems as Parameter estimation without a control, now with a control input and a quadratic control cost — estimating a parameter and a control at the same time. Read that page first; this one skips the narration that doesn't change.

julia
using OptimalControl
using NLPModelsIpopt
using OrdinaryDiffEqTsit5
using NonlinearSolve
using Plots

Growth-rate estimation, with control ​

Definition ​

julia
λ_true = 0.5
x_true(t) = 2 * exp(λ_true * t)
data(t) = x_true(t) + 0.2 * sin(4π * t)

t0 = 0
tf = 2
x0 = 2

ocp_growth = @def begin
    λ ∈ R, variable
    t ∈ [t0, tf], time
    x ∈ R, state
    u ∈ R, control

    x(t0) == x0

    ẋ(t) == λ * x(t) + u(t)

    ∫((x(t) - data(t))^2 + 0.5u(t)^2) → min
end
Abstract definition:

    λ ∈ R, variable
    t ∈ [t0, tf], time
    x ∈ R, state
    u ∈ R, control
    x(t0) == x0
    ẋ(t) == λ * x(t) + u(t)
    ∫((x(t) - data(t)) ^ 2 + 0.5 * u(t) ^ 2) → min

The (non autonomous) optimal control problem is of the form:

    minimize  J(x, u, λ) = ∫ f⁰(t, x(t), u(t), λ) dt, over [0, 2]

    subject to

        ẋ(t) = f(t, x(t), u(t), λ), t in [0, 2] a.e.,

        ϕ₋ ≤ ϕ(x(0), x(2), λ) ≤ ϕ₊, 

    where x(t) ∈ R, u(t) ∈ R and λ ∈ R.

Direct solution ​

julia
growth_sol = solve(ocp_growth; grid_size=20, display=false)
println("estimated λ = ", variable(growth_sol))
estimated λ = 0.4926878314982767
julia
t_grid = time_grid(growth_sol)
plt_growth = plot(growth_sol, :state, :control; label="Direct")

Indirect solution ​

The data term makes the Lagrange integrand — and so the pseudo-Hamiltonian       — genuinely non-autonomous. Maximising over gives  :

julia
u_growth(t, x, p, λ) = p
f_growth = Flow(ocp_growth, u_growth)

function shoot_growth!(s, ξ)
    p0, λ = ξ[1], ξ[2]
    _, p_tf, pλ_tf = f_growth(t0, x0, p0, tf; variable=λ, variable_costate=true)
    s[1] = p_tf
    s[2] = pλ_tf
    return nothing
end
shoot_growth! (generic function with 1 method)
julia
nle_growth!(s, ξ, _) = shoot_growth!(s, ξ)

p0_guess = costate(growth_sol)(t0)
λ_guess = variable(growth_sol)

prob_growth = NonlinearProblem(nle_growth!, [p0_guess, λ_guess])
shoot_sol_growth = NonlinearSolve.solve(prob_growth; show_trace=Val(false))
p0_sol_growth, λ_sol = shoot_sol_growth.u
2-element Vector{Float64}:
 0.04754487287256403
 0.49342678063916473

Comparison ​

julia
indirect_growth = f_growth((t0, tf), x0, p0_sol_growth; variable=λ_sol)
plot!(plt_growth, indirect_growth, :state, :control; label="Indirect", linestyle=:dash)

Harmonic oscillator, with control ​

Definition ​

julia
q0 = 1
v0 = 0
t0h = 0
tfh = 1

ocp_harmonic = @def begin
    ω ∈ R, variable
    t ∈ [t0h, tfh], time
    x = (q, v) ∈ R², state
    u ∈ R, control

    q(t0h) == q0
    v(t0h) == v0
    q(tfh) == 0.0

    ẋ(t) == [v(t), -ω^2 * q(t) + u(t)]

    ω^2 + 0.5∫(u(t)^2) → min
end
Abstract definition:

    ω ∈ R, variable
    t ∈ [t0h, tfh], time
    x = ((q, v) ∈ R², state)
    u ∈ R, control
    q(t0h) == q0
    v(t0h) == v0
    q(tfh) == 0.0
    ẋ(t) == [v(t), -(ω ^ 2) * q(t) + u(t)]
    ω ^ 2 + 0.5 * ∫(u(t) ^ 2) → min

The (autonomous) optimal control problem is of the form:

    minimize  J(x, u, ω) = g(x(0), x(1), ω) + ∫ f⁰(x(t), u(t), ω) dt, over [0, 1]

    subject to

        ẋ(t) = f(x(t), u(t), ω), t in [0, 1] a.e.,

        ϕ₋ ≤ ϕ(x(0), x(1), ω) ≤ ϕ₊, 

    where x(t) = (q(t), v(t)) ∈ R², u(t) ∈ R and ω ∈ R.

only ever appears squared, in both the dynamics and the cost, so its sign carries no physical meaning here — unlike the control-free version, where   was forced by the boundary conditions alone. With a control now sharing the work of steering to , the optimal is smaller: the direct and indirect solves below are expected to agree with each other, not with .

Direct solution ​

julia
harmonic_sol = solve(ocp_harmonic; grid_size=20, display=false)
println("estimated ω = ", variable(harmonic_sol))
estimated ω = -0.6535099333208557
julia
plot(harmonic_sol, :state, :control; label="Direct")

Indirect solution ​

Here neither the dynamics nor the cost depend explicitly on , so the flow is autonomous. Maximising       over gives  :

julia
u_harmonic(x, p, ω) = p[2]
f_harmonic = Flow(ocp_harmonic, u_harmonic)

function shoot_harmonic!(s, ξ)
    p0, ω = ξ[1:2], ξ[3]
    x_tf, p_tf, pω_tf = f_harmonic(t0h, [q0, v0], p0, tfh; variable=ω, variable_costate=true)
    s[1] = x_tf[1]
    s[2] = p_tf[2]
    s[3] = pω_tf + 2ω
    return nothing
end
shoot_harmonic! (generic function with 1 method)
julia
nle_harmonic!(s, ξ, _) = shoot_harmonic!(s, ξ)

p0_guess_h = costate(harmonic_sol)(t0h)
ω_guess = variable(harmonic_sol)

prob_harmonic = NonlinearProblem(nle_harmonic!, [p0_guess_h..., ω_guess])
shoot_sol_harmonic = NonlinearSolve.solve(prob_harmonic; show_trace=Val(false))
p0_sol_harmonic, ω_sol = shoot_sol_harmonic.u[1:2], shoot_sol_harmonic.u[3]
([-2.059399404655115, -2.4134560681767963], -0.6537637057182825)

Comparison ​

julia
indirect_harmonic = f_harmonic((t0h, tfh), [q0, v0], p0_sol_harmonic; variable=ω_sol)
plt_harmonic = plot(harmonic_sol, :state, :control; label="Direct")
plot!(plt_harmonic, indirect_harmonic, :state, :control; label="Indirect", linestyle=:dash)

See also ​