Indirect simple shooting
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 solutionOptimal 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$.
t0 = 0
tf = 1
x0 = -1
xf = 0
α = 1.5
ocp = @def begin
t ∈ [t0, tf], time
x ∈ R, state
u ∈ R, control
x(t0) == x0
x(tf) == xf
ẋ(t) == -x(t) + α * x(t)^2 + u(t)
∫( 0.5u(t)^2 ) → min
end;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)))$.
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)We define also the projection function on the state space.
π((x, p)) = xActually, $\varphi_{t_0, x_0, p_0}(\cdot)$ is also solution of
\[ \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 functionResolution 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 guessNonlinearSolve.jl
We first use the NonlinearSolve.jl package to solve the shooting equation. Let us define the problem.
nle! = (s, ξ, λ) -> s[1] = S(ξ[1]) # auxiliary function
prob = NonlinearProblem(nle!, ξ) # NLE problem with initial guessLet us do some benchmarking.
using BenchmarkTools
@benchmark solve(prob; show_trace=Val(false))BenchmarkTools.Trial: 534 samples with 1 evaluation.
Range (min … max): 8.171 ms … 24.678 ms ┊ GC (min … max): 0.00% … 54.77%
Time (median): 8.732 ms ┊ GC (median): 0.00%
Time (mean ± σ): 9.366 ms ± 2.834 ms ┊ GC (mean ± σ): 6.03% ± 11.93%
▅█▃
▅███▅▅▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▄▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇█ ▆
8.17 ms Histogram: log(frequency) by time 22.8 ms <
Memory estimate: 6.38 MiB, allocs estimate: 142382.For small nonlinear systems, it could be faster to use the SimpleNewtonRaphson() descent algorithm.
@benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false))BenchmarkTools.Trial: 852 samples with 1 evaluation.
Range (min … max): 5.131 ms … 18.139 ms ┊ GC (min … max): 0.00% … 65.58%
Time (median): 5.493 ms ┊ GC (median): 0.00%
Time (mean ± σ): 5.868 ms ± 2.009 ms ┊ GC (mean ± σ): 5.86% ± 11.40%
█▅
▆██▅▆▁▁▃█▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆▇ ▇
5.13 ms Histogram: log(frequency) by time 17.4 ms <
Memory estimate: 4.72 MiB, allocs estimate: 88155.Now, let us solve the problem and retrieve the initial costate solution.
indirect_sol = solve(prob; show_trace=Val(true)) # resolution of S(p0) = 0
p0_sol = indirect_sol.u[1] # costate solution
println("\ncostate: p0 = ", p0_sol)
println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
Algorithm: NewtonRaphson(
descent = NewtonDescent()
)
---- ------------- -----------
Iter f(u) inf-norm Step 2-norm
---- ------------- -----------
0 6.71224852e-02 0.00000000e+00
1 1.58897952e-03 2.86166386e-02
2 8.72173219e-07 6.46258281e-04
3 2.62937541e-13 3.55113667e-07
Final 2.62937541e-13
----------------------
costate: p0 = 0.07202997482156905
shoot: |S(p0)| = 2.629375407815329e-13MINPACK.jl
Instead of the NonlinearSolve.jl package we can use the MINPACK.jl package to solve the shooting equation. To compute the Jacobian of the shooting function we use the DifferentiationInterface.jl package with ForwardDiff.jl backend.
using DifferentiationInterface
import ForwardDiff
backend = AutoForwardDiff()Let us define the problem to solve.
nle! = ( s, ξ) -> s[1] = S(ξ[1]) # auxiliary function
jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ) # Jacobian of nleWe are now in position to solve the problem with the hybrj solver from MINPACK.jl through the fsolve function, providing the Jacobian. Let us do some benchmarking.
@benchmark fsolve(nle!, jnle!, ξ; show_trace=false) # initial guess given to the solverBenchmarkTools.Trial: 729 samples with 1 evaluation.
Range (min … max): 6.240 ms … 29.590 ms ┊ GC (min … max): 0.00% … 51.82%
Time (median): 6.391 ms ┊ GC (median): 0.00%
Time (mean ± σ): 6.856 ms ± 2.358 ms ┊ GC (mean ± σ): 5.19% ± 10.55%
█▃
███▁▁▇▁▁▁▄▁▁▁▄▆▄▄▁▄▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▇ ▇
6.24 ms Histogram: log(frequency) by time 20.2 ms <
Memory estimate: 4.30 MiB, allocs estimate: 105924.We can also use the preparation step of DifferentiationInterface.jl.
extras = prepare_jacobian(nle!, similar(ξ), backend, ξ)
jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras)
@benchmark fsolve(nle!, jnle_prepared!, ξ; show_trace=false)BenchmarkTools.Trial: 736 samples with 1 evaluation.
Range (min … max): 5.987 ms … 30.512 ms ┊ GC (min … max): 0.00% … 47.35%
Time (median): 6.388 ms ┊ GC (median): 0.00%
Time (mean ± σ): 6.794 ms ± 2.292 ms ┊ GC (mean ± σ): 5.13% ± 10.49%
█▃
▇███▁▄▇▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▇ ▇
5.99 ms Histogram: log(frequency) by time 20.1 ms <
Memory estimate: 4.30 MiB, allocs estimate: 105920.Now, 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")Iter f(x) inf-norm Step 2-norm Step time
------ -------------- -------------- --------------
1 6.712249e-02 0.000000e+00 0.000985
2 1.588980e-03 8.189120e-04 0.002103
3 3.722710e-05 4.379408e-07 0.000902
4 2.043880e-08 2.294999e-10 0.000877
5 2.627965e-13 6.925509e-17 0.000884
6 7.415082e-18 1.144852e-26 0.000912
costate: p0 = 0.0720299748216761
shoot: |S(p0)| = 7.415081784167428e-18Plot 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.
pretty_plot — Function
pretty_plot(S, p0_sol; size=(800, 450))