Studied problem

We are looking now to a more complex problem where there are two kinds of masses attached to the membrane, described by two resistances $r_1$ and $r_2$.

Since the code provide in the Filtration.jl package only look for the case of one resistance/mass, this page proposes to show how one can use the control-toolbox ecosystem throught the OptimalControl.jl package to get the optimal solution of the following problem.

using Plots
using OptimalControl

t0 = 0.
vf = 10.
x0 = [0., 0., 0.1, 0.1]
p = [1., 1., 1., 1., 2., 1.5]

e(x) = p[1] + x[3] + x[4]
g(x) = p[2]
f₁₁(x) = p[3]
f₁₂(x) = - p[4]*x[3]
f₂₁(x) = p[5]
f₂₂(x) = - p[6]*x[4]

# Vector fields
Ff(x) = [e(x);  g(x); f₁₁(x); f₂₁(x)]
Fb(x) = [e(x); -g(x); f₁₂(x); f₂₂(x)]
F0(x) = Ff(x) + Fb(x)
F1(x) = Ff(x) - Fb(x)

ocp = @def begin
    tf ∈ R, variable
    t ∈ [t0, tf], time
    x = (e, v, r1, r2) ∈ R⁴, state
    u ∈ R, control
    -1 ≤ u(t) ≤ 1
    x(t0) == x0
    v(tf) == vf
    ẋ(t) == F0(x(t)) + u(t)*F1(x(t))
    e(tf) → min
end

Abstract definition:

    tf ∈ R, variable
    t ∈ [t0, tf], time
    x = ((e, v, r1, r2) ∈ R⁴, state)
    u ∈ R, control
    -1 ≤ u(t) ≤ 1
    x(t0) == x0
    v(tf) == vf
    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
    e(tf) → min

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

    minimize  J(x, u, tf) = g(x(0.0), x(tf), tf)

    subject to

        ẋ(t) = f(x(t), u(t), tf), t in [0.0, tf] a.e.,

        ϕ₋ ≤ ϕ(x(0.0), x(tf), tf) ≤ ϕ₊, 
        u₋ ≤ u(t) ≤ u₊, 

    where x(t) = (e(t), v(t), r1(t), r2(t)) ∈ R⁴, u(t) ∈ R and tf ∈ R.

Direct method

One classical optimal control method used to get a solution of the problem is the direct one. The main idea is to discretize in time the studied optimal control problem, to formulate the associated large optimization problem and to solve it. Thanks to the OptimalControl.jl package, one can easily have and plot the solution of the direct method by using the solve function. Here we use the interior point Ipopt solver to get a solution of the discretized problem.

using NLPModelsIpopt

direct_sol = solve(ocp; display = true)
plt_sol = plot(direct_sol, label = "direct")

Structure of the solution

As we can see, the structure of the optimal solution is composed by a bang arc $u = 1$, followed by a singular arc $u = u_s(x,\lambda)$ and finished by a bang arc $u = +1$. To verify this, we vaan verify that the switching function $t \mspato H_1(x(t), \lambda(t))$ vanished along the singular arc. The expression of the singular control is simply computed thanks to the OptimalControl.jl package, see this tutorial for more information about how it is done.

Note

For the computation of the singular control, one suppose that there exists $I \subset [t_0, t_f]$ such that $H_1(x(t), \lambda(t)) = 0$ for almost all $t \in I$. Denoting $\{H, G\}$ the Poisson bracket of two Hamiltonians $H$ and $G$ defined by

\[ \{H, G\} = (\nabla p H, -\nabla x H) \cdot G\]

the minimal order singular control is given by

\[ u_s(x,\lambda) = -\frac{H_{001}}{H_{101}} = -\frac{\{H_0,\{H_0, H_1\}\}}{\{H_1, \{H_0, H_1\}\}}.\]

using OrdinaryDiffEq

# Lift into (x,λ) space
H0 = Lift(F0)
H1 = Lift(F1)

# Lie bracket
H01  = @Lie {H0, H1}
H001 = @Lie {H0, H01}
H101 = @Lie {H1, H01}

# Singular control
us(x, λ) = -H001(x, λ) / H101(x, λ)

# Pseudo-Hamiltonian
H(x,λ,u) = H0(x,λ) + u*H1(x,λ)

# Flows (Change for plot the control)
ϕ0 = Flow(ocp, (x,λ,tf) -> -1)
ϕ1 = Flow(ocp, (x,λ,tf) -> +1)
ϕs = Flow(ocp, (x,λ,tf) -> us(x,λ))

# Get direct trajectory
time = time_grid(direct_sol)
x = state(direct_sol)
u = control(direct_sol)
λ = costate(direct_sol)
tf = time[end]

# Structure of the solution
plot( t -> H0(x(t), λ(t)), t0, tf, label = "H₀(x(t), λ(t))")
plot!(t -> H1(x(t), λ(t)), t0, tf, label = "H₁(x(t), λ(t))")
plot!(t -> H01(x(t), λ(t)), t0, tf, label = "H₀₁(x(t), λ(t))")

Indirect shooting

We are now able to solve the studied optimal control problem by indirect shooting method. The goal is to find a zero of the function $S$ defined below. Thanks to the direct solution, we can have a good initial guess of the zero of $S$.

using LinearAlgebra: norm

# Shooting function
function shoot!(s, ξ)
    λv0, λr10, λr20, t1, t2, tf = ξ
    x1, λ1 = ϕ1(t0, x0, [-1, λv0, λr10, λr20], t1)
    x2, λ2 = ϕs(t1, x1, λ1, t2)
    xf, λf = ϕ1(t2, x2, λ2, tf)

    s[1] = xf[2] - vf
    s[2:3] = λf[3:4]
    s[4] = H0(x1, λ1)
    s[5] = H1(x1, λ1)
    s[6] = H01(x1, λ1)
end

# Jacobian of the shooting function
jshoot! = (js, ξ) -> jacobian!(shoot!, similar(ξ), js, AutoForwardDiff(), ξ)

# Initial guess
λ0 = λ(t0)
η = 1e-3
time_ = time[ u.(time) .≤ 1-η ]
t1 = time_[1]; t2 = time_[end]
ξ = [λ0[2:4]..., t1, t2, tf]

s = similar(ξ)
shoot!(s, ξ)
println("\nNorm of the initial guess of the shooting function: ‖s‖ = ", norm(s), "\n")


Norm of the initial guess of the shooting function: ‖s‖ = 0.43040862022365656

We propose to use hybrd1 method which is a modified version of Powell's algorithm to find a zero of $S$, given by the MINPACK.jl package.

using MINPACK
using DifferentiationInterface
import ForwardDiff

# Resolution of S(ξ) = 0
indirect_sol = fsolve(shoot!, jshoot!, ξ, show_trace=true)
shoot!(s, indirect_sol.x)
println("\nNorm of the solution of the shooting function: ‖s‖ = ", norm(s), "\n")

# Plot
λv0, λr10, λr20, t1, t2, tf = indirect_sol.x
ϕ = ϕ1 * (t1, ϕs) * (t2, ϕ1)
flow_sol = ϕ((t0, tf), x0, λ0; saveat=range(t0, tf, 1000))
plot!(plt_sol, flow_sol, label="indirect")

Turnpike property

Our goal now is simply to numerically highlight the turnpike property of our ploblem. To do this, we found steady state on $r_1$ and $r_2$ which minimize $\dot e - \lambda_v \dot v$ where $\lambda_v$ coresponds to the optimal constant costate associated to $v$. To do this, we use the Ipopt solver throught the Optimization.jl package.

using Ipopt, Optimization, OptimizationMOI
# Constraints
function cons!(dξ, ξ,_)
    r1, r2, u = ξ
    x = [0., 0., r1, r2]
    dx = F0(x) + u*F1(x)
    dξ .= dx[3:4]
end

# Objective
function obj(ξ,_)
    r1, r2, u = ξ
    x = [0., 0., r1, r2]
    dx = F0(x) + u*F1(x)
    return dx[1] - λv0*dx[2]
end

# Initial guess
x, λ = ϕ(t0, x0, λ0, (t1 + t2)/2)
u = us(x, λ)
ξ = [x[3:4]; u]

# Definition of the optimization problem
optprob = OptimizationFunction(obj, AutoForwardDiff(), cons = cons!)
prob = OptimizationProblem(optprob, ξ, nothing, lcons = [0., 0.], ucons = [0., 0.])
opt_sol = Optimization.solve(prob, Ipopt.Optimizer(), print_level = 5)

# Plot
st = opt_sol.u
plt_sol = plot(flow_sol, label = "")
plot!(plt_sol, [t0, tf], [st[1], st[1]], subplot = 3, label = "")
plot!(plt_sol, [t0, tf], [st[2], st[2]], subplot = 4, label = "")
plot!(plt_sol, [t0, tf], [st[3], st[3]], subplot = 9, label = "")