[Indirect simple shooting](@id iss)

In this tutorial we present the indirect simple shooting method on a simple example.

Let us start by importing the necessary packages.

using OptimalControl    # to define the optimal control problem and its flow
using OrdinaryDiffEq    # to get the Flow function from OptimalControl
using NonlinearSolve    # interface to NLE solvers
using MINPACK           # NLE solver: use to solve the shooting equation
using Plots             # to plot the solution

Optimal control problem

Let us consider the following optimal control problem:

$$\left\{ \begin{array}{l} \min \displaystyle \frac{1}{2} \int_{t_0}^{t_f} u^2(t) \, \mathrm{d} t\\[1.0em] \dot{x}(t) = \displaystyle -x(t)+\alpha x^2(t)+u(t), \quad u(t) \in \R, \quad t \in [t_0, t_f] \text{ a.e.}, \\[0.5em] x(t_0) = x_0, \quad x(t_f) = x_f, \end{array} \right.%$$

with $t_0 = 0$, $t_f = 1$, $x_0 = -1$, $x_f = 0$, $\alpha=1.5$ and $\forall, t \in [t_0, t_f]$, $x(t) \in \R$.

const t0 = 0
const tf = 1
const x0 = -1
const xf = 0
const α  = 1.5
ocp = @def begin

    t ∈ [t0, tf], time
    x ∈ R, state
    u ∈ R, control

    x(t0) == x0
    x(tf) == xf

    ẋ(t) == -x(t) + α * x(t)^2 + u(t)

    ∫( 0.5u(t)^2 ) → min
    
end
nothing # hide

Boundary value problem

The pseudo-Hamiltonian of this problem is

$$H(x,p,u) = p \, (-x+\alpha x^2+u) + p^0 u^2 /2,$$

where $p^0 = -1$ since we are in the normal case. From the Pontryagin Maximum Principle, the maximising control is given by

$$u(x, p) = p$$

since $\partial^2_{uu} H = p^0 = - 1 < 0$. Plugging this control in feedback form into the pseudo-Hamiltonian, and considering the limit conditions, we obtain the following two-points boundary value problem (BVP).

$$\left\{ \begin{array}{l} \dot{x}(t) = \phantom{-} \nabla_p H[t] = -x(t) + \alpha x^2(t) + u(x(t), p(t)) = -x(t) + \alpha x^2(t) + p(t), \\[0.5em] \dot{p}(t) = - \nabla_x H[t] = (1 - 2 \alpha x(t))\, p(t), \\[0.5em] x(t_0) = x_0, \quad x(t_f) = x_f, \end{array} \right.$$

where $[t]~= (x(t),p(t),u(x(t), p(t)))$.

!!! note "Our goal"

Our goal is to solve this (BVP). Solving (BVP) consists in solving the Pontryagin Maximum Principle which provides necessary conditions of optimality.

Shooting function

To achive our goal, let us first introduce the pseudo-Hamiltonian vector field

$$\vec{H}(z,u) = \left( \nabla_p H(z,u), -\nabla_x H(z,u) \right), \quad z = (x,p),$$

and then denote by $\varphi_{t_0, x_0, p_0}(\cdot)$ the solution of the following Cauchy problem

$$\dot{z}(t) = \vec{H}(z(t), u(z(t))), \quad z(t_0) = (x_0, p_0).$$

Our goal becomes to solve

$$\pi( \varphi_{t_0, x_0, p_0}(t_f) ) = x_f,$$

where $\pi(x, p) = x$. To compute $\varphi$ with OptimalControl.jl package, we define the flow of the associated Hamiltonian vector field by:

u(x, p) = p
φ = Flow(ocp, u)
nothing # hide

We define also the projection function on the state space.

π((x, p)) = x
nothing # hide

!!! note "Nota bene"

Actually, $\varphi_{t_0, x_0, p_0}(\cdot)$ is also solution of

```math
    \dot{z}(t) = \vec{\mathbf{H}}(z(t)), \quad z(t_0) = (x_0, p_0),
```
where $\mathbf{H}(z) = H(z, u(z))$ and $\vec{\mathbf{H}} = (\nabla_p \mathbf{H}, -\nabla_x \mathbf{H})$. This is what is actually computed by `Flow`.

Now, to solve the (BVP) we introduce the shooting function:

$$\begin{array}{rlll} S \colon & \R & \longrightarrow & \R \\\ & p_0 & \longmapsto & S(p_0) = \pi( \varphi_{t_0, x_0, p_0}(t_f) ) - x_f. \end{array}$$

S(p0) = π( φ(t0, x0, p0, tf) ) - xf    # shooting function
nothing # hide

Resolution of the shooting equation

At the end, solving (BVP) is equivalent to solve $S(p_0) = 0$. This is what we call the indirect simple shooting method. We define an initial guess.

ξ = [0.1]    # initial guess
nothing # hide

MINPACK.jl

We can use NonlinearSolve.jl package or, instead, MINPACK.jl to solve the shooting equation. To compute the Jacobian of the shooting function we use DifferentiationInterface.jl with ForwardDiff.jl backend.

using MINPACK
function fsolve(f, j, x; kwargs...)
    try
        MINPACK.fsolve(f, j, x; kwargs...)
    catch e
        println("Erreur using MINPACK")
        println(e)
        println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.")
        MINPACK.fsolve(f, x; kwargs...)
    end
end

using DifferentiationInterface
import ForwardDiff
backend = AutoForwardDiff()
nothing # hide

Let us define the problem to solve.

nle!  = ( s, ξ) -> s[1] = S(ξ[1])                                 # auxiliary function
jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ)    # Jacobian of nle
nothing # hide

We are now in position to solve the problem with the hybrj solver from MINPACK.jl through the fsolve function, providing the Jacobian. Let us solve the problem and retrieve the initial costate solution.

indirect_sol = fsolve(nle!, jnle!, ξ; show_trace=true)    # resolution of S(p0) = 0
p0_sol = indirect_sol.x[1]                                # costate solution
println("\ncostate:    p0 = ", p0_sol)
println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")

Plot of the solution

The solution can be plot calling first the flow.

sol = φ((t0, tf), x0, p0_sol)
plot(sol)

In the indirect shooting method, the research of the optimal control is replaced by the computation of its associated extremal. This computation is equivalent to finding the initial covector solution to the shooting function. Let us plot the extremal in the phase space and the shooting function with the solution.

<article class="docstring">
<header>
    <a class="docstring-article-toggle-button fa-solid fa-chevron-right" href="javascript:;" title="Expand docstring"> </a>
    <code>pretty_plot</code> — <span class="docstring-category">Function</span>
</header>
<section style="display: none;"><div><pre><code class="language-julia hljs">using Plots.PlotMeasures

function pretty_plot(S, p0; Np0=20, kwargs...) 
 
    # times for wavefronts
    times = range(t0, tf, length=3)

    # times for trajectories
    tspan = range(t0, tf, length=100)

    # interval of initial covector
    p0_min = -0.5 
    p0_max = 2 

    # covector solution
    p0_sol = p0 
 
    # plot of the flow in phase space
    plt_flow = plot() 
    p0s = range(p0_min, p0_max, length=Np0) 
    for i ∈ eachindex(p0s) 
        sol = φ((t0, tf), x0, p0s[i])
        x = state(sol).(tspan)
        p = costate(sol).(tspan)
        label = i==1 ? "extremals" : false 
        plot!(plt_flow, x, p, color=:blue, label=label) 
    end 
 
    # plot of wavefronts in phase space 
    p0s = range(p0_min, p0_max, length=200) 
    xs  = zeros(length(p0s), length(times)) 
    ps  = zeros(length(p0s), length(times)) 
    for i ∈ eachindex(p0s) 
        sol = φ((t0, tf), x0, p0s[i], saveat=times)
        xs[i, :] .= state(sol).(times) 
        ps[i, :] .= costate(sol).(times) 
    end 
    for j ∈ eachindex(times) 
        label = j==1 ? "flow at times" : false 
        plot!(plt_flow, xs[:, j], ps[:, j], color=:green, linewidth=2, label=label) 
    end 
 
    #  
    plot!(plt_flow, xlims=(-1.1, 1), ylims=(p0_min, p0_max)) 
    plot!(plt_flow, [0, 0], [p0_min, p0_max], color=:black, xlabel="x", ylabel="p", label="x=xf") 
     
    # solution 
    sol = φ((t0, tf), x0, p0_sol)
    x = state(sol).(tspan)
    p = costate(sol).(tspan)
    plot!(plt_flow, x, p, color=:red, linewidth=2, label="extremal solution") 
    plot!(plt_flow, [x[end]], [p[end]], seriestype=:scatter, color=:green, label=false) 
 
    # plot of the shooting function  
    p0s = range(p0_min, p0_max, length=200) 
    plt_shoot = plot(xlims=(p0_min, p0_max), ylims=(-2, 4), xlabel="p₀", ylabel="y") 
    plot!(plt_shoot, p0s, S, linewidth=2, label="S(p₀)", color=:green) 
    plot!(plt_shoot, [p0_min, p0_max], [0, 0], color=:black, label="y=0") 
    plot!(plt_shoot, [p0_sol, p0_sol], [-2, 0], color=:black, label="p₀ solution", linestyle=:dash) 
    plot!(plt_shoot, [p0_sol], [0], seriestype=:scatter, color=:green, label=false) 
 
    # final plot 
    plot(plt_flow, plt_shoot; layout=(1,2), leftmargin=15px, bottommargin=15px, kwargs...) 
 
end</code><button class="copy-button fa-solid fa-copy" aria-label="Copy this code ;opblock" title="Copy"></button></pre></div>
</section>
</article>

using Plots.PlotMeasures # hide
function pretty_plot(S, p0; Np0=20, kwargs...) # hide
 # hide
    # times for wavefronts# hide
    times = range(t0, tf, length=3)# hide
# hide
    # times for trajectories# hide
    tspan = range(t0, tf, length=100)# hide
# hide
    # interval of initial covector# hide
    p0_min = -0.5 # hide
    p0_max = 2 # hide
# hide
    # covector solution# hide
    p0_sol = p0 # hide
 # hide
    # plot of the flow in phase space# hide
    plt_flow = plot() # hide
    p0s = range(p0_min, p0_max, length=Np0) # hide
    for i ∈ eachindex(p0s) # hide
        sol = φ((t0, tf), x0, p0s[i])# hide
        x = state(sol).(tspan)# hide
        p = costate(sol).(tspan)# hide
        label = i==1 ? "extremals" : false # hide
        plot!(plt_flow, x, p, color=:blue, label=label) # hide
    end # hide
 # hide
    # plot of wavefronts in phase space # hide
    p0s = range(p0_min, p0_max, length=200) # hide
    xs  = zeros(length(p0s), length(times)) # hide
    ps  = zeros(length(p0s), length(times)) # hide
    for i ∈ eachindex(p0s) # hide
        sol = φ((t0, tf), x0, p0s[i], saveat=times)# hide
        xs[i, :] .= state(sol).(times) # hide
        ps[i, :] .= costate(sol).(times) # hide
    end # hide
    for j ∈ eachindex(times) # hide
        label = j==1 ? "flow at times" : false # hide
        plot!(plt_flow, xs[:, j], ps[:, j], color=:green, linewidth=2, label=label) # hide
    end # hide
 # hide
    #  # hide
    plot!(plt_flow, xlims=(-1.1, 1), ylims=(p0_min, p0_max)) # hide
    plot!(plt_flow, [0, 0], [p0_min, p0_max], color=:black, xlabel="x", ylabel="p", label="x=xf") # hide
     # hide
    # solution # hide
    sol = φ((t0, tf), x0, p0_sol)# hide
    x = state(sol).(tspan)# hide
    p = costate(sol).(tspan)# hide
    plot!(plt_flow, x, p, color=:red, linewidth=2, label="extremal solution") # hide
    plot!(plt_flow, [x[end]], [p[end]], seriestype=:scatter, color=:green, label=false) # hide
 # hide
    # plot of the shooting function  # hide
    p0s = range(p0_min, p0_max, length=200) # hide
    plt_shoot = plot(xlims=(p0_min, p0_max), ylims=(-2, 4), xlabel="p₀", ylabel="y") # hide
    plot!(plt_shoot, p0s, S, linewidth=2, label="S(p₀)", color=:green) # hide
    plot!(plt_shoot, [p0_min, p0_max], [0, 0], color=:black, label="y=0") # hide
    plot!(plt_shoot, [p0_sol, p0_sol], [-2, 0], color=:black, label="p₀ solution", linestyle=:dash) # hide
    plot!(plt_shoot, [p0_sol], [0], seriestype=:scatter, color=:green, label=false) # hide
 # hide
    # final plot # hide
    plot(plt_flow, plt_shoot; layout=(1,2), leftmargin=15px, bottommargin=15px, kwargs...) # hide
 # hide
end# hide
nothing # hide

pretty_plot(S, p0_sol; size=(800, 450))

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[Indirect simple shooting](@id iss)

Optimal control problem

Boundary value problem

Shooting function

Resolution of the shooting equation

MINPACK.jl

Plot of the solution

FilesExpand file tree

tutorial-iss.md

Latest commit

History

tutorial-iss.md

File metadata and controls

[Indirect simple shooting](@id iss)

Optimal control problem

Boundary value problem

Shooting function

Resolution of the shooting equation

MINPACK.jl

Plot of the solution