Control and variable together

Some problems optimise a control and constant parameters at the same time: to identify a parameter and the input that explains the data, or to find a parameter and the control law that go best together. This page takes the two systems of Parameter estimation without a control, and adds a control to their dynamics, with a quadratic cost. The method is the same; this page shows what the control changes.

using OptimalControl
using NLPModelsIpopt
using OrdinaryDiffEqTsit5
using NonlinearSolve
using Plots

Growth rate, with a control

The problem

The growth model gets an additive control,   , and the cost a term : the control can explain the part of the data that the model does not, at a price.

λ_true = 0.5
x_obs(t) = 2 * exp(λ_true * t) + 0.2 * sin(4π * t)

t0, tf, x0 = 0, 2, 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) - x_obs(t))^2 + 0.5u(t)^2) → min
end

Direct method

growth_sol = solve(ocp_growth; grid_size=20, display=false)
variable(growth_sol), objective(growth_sol)

(0.4926878314982767, 0.038698292831307775)

plt_growth = plot(growth_sol, :state, :control; label="direct", size=(800, 400))

Indirect method

The pseudo-Hamiltonian now depends on the control:

It is maximised where    , so the control in feedback form is  . The augmented system is the one of the control-free problem, with gaining the control:

with the same conditions,   and   from  . The flow takes the control law, a function of since the problem depends on time through the data:

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

function shoot_growth!(s, ξ, _)
    p0, λ = ξ
    _, p_tf, pλ_tf = f_growth(t0, x0, p0, tf; variable=λ, variable_costate=true)
    s[1] = p_tf      # free final state
    s[2] = pλ_tf     # free parameter
    return nothing
end

ξ_guess = [costate(growth_sol)(t0), variable(growth_sol)]
prob_growth = NonlinearProblem(shoot_growth!, ξ_guess)
p0_growth, λ_sol = NonlinearSolve.solve(prob_growth; show_trace=Val(false)).u

2-element Vector{Float64}:
 0.0475448728725648
 0.49342678063916473

Comparison

indirect_growth = f_growth((t0, tf), x0, p0_growth; variable=λ_sol)
(variable(growth_sol), λ_sol), (objective(growth_sol), objective(indirect_growth))

((0.4926878314982767, 0.49342678063916473), (0.038698292831307775, 0.0387233975043203))

For the plot, the flow is built with saveat, so that the trajectory has 201 points (see Plotting a flow):

f_growth_plot = Flow(ocp_growth, u_growth; saveat=range(t0, tf, 201), dense=false)
plot!(plt_growth, f_growth_plot((t0, tf), x0, p0_growth; variable=λ_sol), :state, :control;
      label="indirect", linestyle=:dash)

The two methods agree, up to the discretisation error of the 20-step grid.

Pulsation of an oscillator, with a control

The problem

The oscillator gets a control force,    , and the cost the energy of the control:  . The contrast with the control-free problem is the point of this example. There, the boundary conditions alone forced  . Here, any can reach the target, with the help of the control, and the optimum is a trade-off: a larger costs more, but the oscillator then does more of the work. Only appears in the problem, so and give the same cost.

q0, v0 = 1, 0
t0h, tfh = 0, 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

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

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

Direct method

harmonic_sol = solve(ocp_harmonic; grid_size=20, display=false)
variable(harmonic_sol), objective(harmonic_sol)

(-0.6535099333208557, 1.4581596292455086)

plt_harmonic = plot(harmonic_sol, :state, :control; label="direct", size=(800, 400))

Indirect method

The pseudo-Hamiltonian is

maximised by  . The augmented system is the one of the control-free problem, with    ; the adjoint equations and the conditions do not change:  , and    from  .

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]        # q(tf) = 0
    s[2] = p_tf[2]        # free final velocity
    s[3] = pω_tf + 2ω     # transversality of the parameter
    return nothing
end

ξ_guess_h = [costate(harmonic_sol)(t0h)..., variable(harmonic_sol)]
prob_harmonic = NonlinearProblem(shoot_harmonic!, ξ_guess_h)
shooting_harmonic = NonlinearSolve.solve(prob_harmonic; show_trace=Val(false))
p0_harmonic, ω_sol = shooting_harmonic.u[1:2], shooting_harmonic.u[3]

([-2.0593994046551165, -2.413456068176797], -0.6537637057182826)

Comparison

indirect_harmonic = f_harmonic((t0h, tfh), [q0, v0], p0_harmonic; variable=ω_sol)
(variable(harmonic_sol), ω_sol), (objective(harmonic_sol), objective(indirect_harmonic))

((-0.6535099333208557, -0.6537637057182826), (1.4581596292455086, 1.4571066852531365))

f_harmonic_plot = Flow(ocp_harmonic, u_harmonic; saveat=range(t0h, tfh, 101), dense=false)
plot!(plt_harmonic, f_harmonic_plot((t0h, tfh), [q0, v0], p0_harmonic; variable=ω_sol), :state, :control;
      label="indirect", linestyle=:dash)

The optimal pulsation is  , well below , for a cost of .

The shooting equations have another solution:  ,    . Without oscillator, the control alone steers the mass,    , and the cost is   . It is an extremal, but not the minimum: its cost is larger. The maximum principle gives necessary conditions only, and comparing the costs of the extremals selects the optimum.

s = zeros(3)
shoot_harmonic!(s, [-3.0, -3.0, 0.0], nothing)
s, objective(f_harmonic((t0h, tfh), [q0, v0], [-3.0, -3.0]; variable=0.0))

([5.509076362406109e-15, -6.080600461210963e-16, 0.0], 1.5000000000000004)

Applications

Problems with a control and parameters appear in many contexts:

• system identification: estimating physical parameters (mass, damping, stiffness) and the inputs, from experimental data;

• optimal design: finding geometric or physical parameters together with the control law that goes with them;

• inverse problems: reconstructing unknown inputs or initial conditions from partial observations, while optimising parameters of the system.

See also

• Parameter estimation without a control: the same two systems, without a control, and the derivation of the augmented system.

• From an OCP: flows with a control law and a variable.