Parameter estimation without a control

A problem without a control optimises constant parameters of a dynamical system: the parameters that make a model fit observed data, or the best parameters for a performance criterion. This page solves two such problems, by the direct and the indirect methods: the growth rate of an exponential model fitted to data, and the smallest pulsation of an oscillator that reaches a target. Control-free problems explains how to write them.

using OptimalControl
using NLPModelsIpopt
using OrdinaryDiffEqTsit5
using NonlinearSolve
using Plots

Growth rate

The problem

A quantity grows exponentially,   , from  , at an unknown rate . We observe it on , and the observations are the model with  , perturbed by a small oscillation. We estimate by least squares. The estimate is close to , but not equal, because of the perturbation.

λ_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

    x(t0) == x0

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

    ∫((x(t) - x_obs(t))^2) → min
end

Direct method

With a coarse grid of 20 steps:

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

(0.4960778956661449, 0.039296407666384314)

t_grid = range(t0, tf, 201)
plt_growth = plot(growth_sol, :state; label="direct", size=(800, 300))
plot!(plt_growth, t_grid, x_obs.(t_grid); line=:dot, color=:black, label="data")

Indirect method

The problem has no control, so the maximum principle gives directly the Hamiltonian, with    (see Notation and conventions):

The parameter is treated as a state with zero dynamics,  , and its own costate . The augmented system, with state and costate , is

The final state is free, so  . The parameter is free too, and does not appear in the cost, so   and  , that is,

The flow of the problem integrates the system in , with passed by keyword. With variable_costate=true, it also integrates from , and returns (see The costate of the variable):

f_growth = Flow(ocp_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

The direct solution gives the initial guess:

ξ_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.05847851053034825
 0.4966212669628472

Comparison

The two estimates differ in the third decimal: the direct method solves a discretised problem, on a grid of 20 steps, and the indirect method the problem itself.

variable(growth_sol), λ_sol

(0.4960778956661449, 0.4966212669628472)

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; saveat=range(t0, tf, 201), dense=false)
indirect_growth = f_growth_plot((t0, tf), x0, p0_growth; variable=λ_sol)
plot!(plt_growth, indirect_growth, :state; label="indirect", linestyle=:dash)

The indirect estimate minimises the least-squares error: the error, computed by quadrature for the model , is larger on both sides of .

Pulsation of an oscillator

The problem

An oscillator    starts at rest at  . We look for the smallest pulsation for which it reaches   at  . The solution is  , so  , and the smallest pulsation is  , with    and  .

q0, v0 = 1, 0
t0h, tfh = 0, 1

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

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

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

    ω^2 → min
end

Direct method

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

(1.5716042803951324, 2.469940014156302)

The direct solution is close to the exact one,   and   :

plt_harmonic = plot(harmonic_sol, :state; label="direct", size=(800, 300))
plot!(plt_harmonic, t -> cos(π * t / 2), 0, 1; subplot=1, label="exact", color=:black, line=:dot)
plot!(plt_harmonic, t -> -π / 2 * sin(π * t / 2), 0, 1; subplot=2, label="exact", color=:black, line=:dot)

Indirect method

The Hamiltonian has no integral term, since the cost is a Mayer cost:

With the state and the costate , the augmented system is

The final velocity is free, so  . For the Mayer cost  , the transversality condition of the parameter is    , with  .

At the solution,  , and the adjoint equations with   give    and     . Then     , so  , and the initial costate is   .

f_harmonic = Flow(ocp_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]

([7.058982847037154e-15, -2.000000000959039], 1.570796326752619)

Comparison

The indirect method finds   and   ; the direct method is close, on its grid of 20 steps:

variable(harmonic_sol), ω_sol, π / 2

(1.5716042803951324, 1.570796326752619, 1.5707963267948966)

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

Applications

Problems without a control appear in many contexts:

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

• optimal design: finding geometric or physical parameters (a length, a stiffness) that optimise a criterion;

• inverse problems: reconstructing unknown inputs or initial conditions from partial observations.

See also

• Control-free problems: how to write them.

• Control and variable together: the same two systems, with a control.

• From an OCP: the flow of a problem without a control, and the costate of the variable.