Problems mixing control and variable

Problems mixing control and variable are optimal control problems that contain both a control variable and a constant parameter (variable) to optimize. They extend control-free problems by adding an explicit control input to the dynamics, while still optimizing constant parameters.

Such problems are used for:

• Identifying unknown parameters and control inputs simultaneously from observed data

• Finding optimal parameters and associated control laws for a given performance criterion

This page demonstrates two examples that extend the control-free problems by adding a control input and a quadratic control cost term.

First, we import the necessary packages:

using OptimalControl
using NLPModelsIpopt
using Plots

Example 1: Exponential growth rate estimation with control

Consider a system with exponential growth and an additive control:

where is an unknown growth rate parameter and is a control input. We have observed data with some perturbations and want to estimate and the optimal control by minimizing the squared error plus a quadratic control cost:

The underlying model has  , but the observed data includes perturbations.

Problem definition

# observed data (analytical solution with λ = 0.5)
λ_true = 0.5
model(t) = 2 * exp(λ_true * t)
perturbation(t) = 2e-1*sin(4π*t)
data(t) = model(t) + perturbation(t)

# optimal control problem (parameter estimation with control)
t0 = 0; tf = 2; x0 = 2
ocp = @def begin
    λ ∈ R, variable              # growth rate to estimate
    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  # fit to observed data with control cost
end

Direct method

direct_sol = solve(ocp; grid_size=20, display=false)

• Solver:
  ✓ Successful  : true
  │  Status     : first_order
  │  Message    : Ipopt/generic
  │  Iterations : 7
  │  Objective  : 0.038698292831307775
  └─ Constraints violation : 1.0239329384376106e-9

• Variable: λ = 0.4926878314982767

• Boundary duals: [0.050309325595601376]

println("Estimated growth rate: λ = ", variable(direct_sol))
println("Objective value: ", objective(direct_sol))

Estimated growth rate: λ = 0.4926878314982767
Objective value: 0.038698292831307775

# plot direct solution
plt = plot(direct_sol; size=(800, 600), label="Direct")

# Add data on first plot
t_grid = time_grid(direct_sol)
plot!(plt, t_grid, data.(t_grid); subplot=1, line=:dot, lw=2, label="Data", color=:black)

The estimated parameter should be close to  .

Indirect method

We now solve the same problem using an indirect shooting method based on Pontryagin's Maximum Principle. First, we import the necessary packages:

using OrdinaryDiffEq  # ODE solver
using NonlinearSolve  # Nonlinear solver

For problems mixing control and variable, we use an augmented Hamiltonian approach with the maximising control. The pseudo-Hamiltonian for this problem is:

The maximisation condition   gives the control in feedback form:

To handle the variable parameter , we treat it as an additional state with zero dynamics. This gives us the augmented system with state and costate , where:

Using the maximising control  , the dynamics of becomes   .

The transversality condition for the variable parameter requires   . Assuming  , we have to satisfy:

We use CTFlows' augment=true feature to automatically compute without manually constructing the augmented system.

# Maximising control from Hamiltonian (non-autonomous: t is required)
u(t, x, p, λ) = p

# Create Hamiltonian flow from OCP with control
f = Flow(ocp, u)

Note

For more details about the flow construction, see this page.

The shooting function enforces the transversality conditions   and  . Using augment=true, the flow automatically returns , with   by construction.

# Shooting function: S(p0, λ) = (p(tf), pλ(tf))
# We want both components to be zero at tf
function shoot!(s, p0, λ)
    _, px_tf, pλ_tf = f(t0, x0, p0, tf, λ; augment=true)
    s[1] = px_tf
    s[2] = pλ_tf
    return nothing
end

# Auxiliary in-place NLE function
nle!(s, y, _) = shoot!(s, y...)

We use the direct solution to initialize the shooting method:

# Extract solution from direct method for initialization
p_direct = costate(direct_sol)
λ_direct = variable(direct_sol)

# Initial guess
p0_guess = p_direct(t0)
λ_guess = λ_direct

# NLE problem with initial guess (2 unknowns: p0, λ)
prob_indirect = NonlinearProblem(nle!, [p0_guess, λ_guess])

# Solve shooting equations
shooting_sol = solve(prob_indirect; show_trace=Val(false))
p0_sol, λ_sol = shooting_sol.u

println("Indirect solution:")
println("Initial costate: p0 = ", p0_sol)
println("Parameter: λ = ", λ_sol)

Indirect solution:
Initial costate: p0 = 0.04754487290047068
Parameter: λ = 0.49342678059355

Finally, we compute and plot the indirect solution:

# Compute and plot indirect solution
indirect_sol = f((t0, tf), x0, p0_sol, λ_sol; saveat=range(t0, tf, 200))
plot!(plt, indirect_sol; linestyle=:dash, lw=2, label="Indirect", color=2)

The direct and indirect solutions match closely, both fitting the perturbed observed data.

Example 2: Harmonic oscillator pulsation optimization with control

Consider a harmonic oscillator with an additive control:

with initial conditions  ,   and final condition  . We want to find the minimal pulsation and optimal control satisfying these constraints:

Without the control term ( ), the analytical solution is   , giving  .

Problem definition

# optimal control problem (pulsation optimization with control)
q0 = 1; v0 = 0
t0 = 0; tf = 1
ocp = @def begin
    ω ∈ R, variable              # pulsation to optimize
    t ∈ [t0, tf], time
    x = (q, v) ∈ R², state
    u ∈ R, control

    q(t0) == q0
    v(t0) == v0
    q(tf) == 0.0                  # final condition

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

    ω^2 + 0.5∫(u(t)^2) → min   # minimize pulsation with control cost
end

Direct method

direct_sol = solve(ocp; grid_size=20, display=false)

• Solver:
  ✓ Successful  : true
  │  Status     : first_order
  │  Message    : Ipopt/generic
  │  Iterations : 15
  │  Objective  : 1.4581596292455086
  └─ Constraints violation : 4.24776880336708e-12

• Variable: ω = -0.6535099333208557

• Boundary duals: [-2.0621687925897234, -2.416078423579069, 2.597222366157977]

println("Optimal pulsation: ω = ", variable(direct_sol))
println("Objective value: ", objective(direct_sol))

Optimal pulsation: ω = -0.6535099333208557
Objective value: 1.4581596292455086

plt = plot(direct_sol; size=(800, 600), label="Direct")

Indirect method

We now solve the same problem using an indirect shooting method. For problems mixing control and variable, we use an augmented Hamiltonian approach with the maximising control. The pseudo-Hamiltonian for this problem is:

The maximisation condition   gives the control in feedback form:

To handle the variable parameter , we treat it as an additional state with zero dynamics. This gives us the augmented system with state and costate , where:

Using the maximising control  , the dynamics of becomes    .

For this problem with a Mayer cost  , the transversality condition for the variable parameter is:

Assuming  , we have:

We use CTFlows' augment=true feature to automatically compute without manually constructing the augmented system:

# Maximising control from Hamiltonian
u(x, p, ω) = p[2]

# Create Hamiltonian flow from OCP with control
f = Flow(ocp, u)

Note

For more details about the flow construction, see this page.

The shooting function enforces the conditions:

• Final condition:  

• Free final velocity:  

• Transversality condition for Mayer cost:   

Using augment=true, the flow automatically returns , with   by construction.

# Shooting function: S(p0, ω)
function shoot!(s, p0, ω)
    x_tf, p_tf, pω_tf = f(t0, [q0, v0], p0, tf, ω; augment=true)
    q_tf = x_tf[1]
    pv_tf = p_tf[2]
    s[1] = q_tf         # q(tf) = 0
    s[2] = pv_tf        # p2(tf) = 0 (free final velocity)
    s[3] = pω_tf + 2ω  # pω(tf) + 2ω = 0 (Mayer cost transversality)
    return nothing
end

# Auxiliary in-place NLE function
nle!(s, y, _) = shoot!(s, y[1:2], y[3])

We use the direct solution to initialize the shooting method:

# Extract solution from direct method for initialization
p_direct = costate(direct_sol)
ω_direct = variable(direct_sol)

# Initial guess
p0_guess = p_direct(t0)
ω_guess = ω_direct

# NLE problem with initial guess
prob_indirect = NonlinearProblem(nle!, [p0_guess..., ω_guess])

# Solve shooting equations
shooting_sol = solve(prob_indirect; show_trace=Val(false))
p0_sol, ω_sol = shooting_sol.u[1:2], shooting_sol.u[3]

println("Indirect solution:")
println("Initial costate: p0 = ", p0_sol)
println("Parameter: ω = ", ω_sol)

Indirect solution:
Initial costate: p0 = [-2.059399403520461, -2.413456067425259]
Parameter: ω = -0.6537637061616014

Finally, we compute and plot the indirect solution:

# Compute and plot indirect solution
indirect_sol = f((t0, tf), [q0, v0], p0_sol, ω_sol; saveat=range(t0, tf, 200))
plot!(plt, indirect_sol; linestyle=:dash, lw=2, label="Indirect", color=2)

The direct and indirect solutions match closely.

Applications

Problems mixing control and variable appear in many contexts:

• System identification: simultaneously estimating physical parameters (mass, damping, stiffness) and control inputs from experimental data

• Optimal design: finding optimal geometric or physical parameters (length, stiffness, etc.) together with associated control laws

• Inverse problems: reconstructing unknown inputs or initial conditions from partial observations while optimizing system parameters

See the syntax documentation for more details on defining control-free problems, and the flow documentation for problems with variables and controls.