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: 538 samples with 1 evaluation.
Range (min … max): 8.471 ms … 25.696 ms ┊ GC (min … max): 0.00% … 57.73%
Time (median): 8.651 ms ┊ GC (median): 0.00%
Time (mean ± σ): 9.304 ms ± 2.933 ms ┊ GC (mean ± σ): 6.33% ± 12.21%
█▅
██▇▅▁▄▄▄▁▁▄▁▁▁▁▅▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▆▇▅ ▆
8.47 ms Histogram: log(frequency) by time 23.7 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: 854 samples with 1 evaluation.
Range (min … max): 5.257 ms … 28.361 ms ┊ GC (min … max): 0.00% … 45.45%
Time (median): 5.414 ms ┊ GC (median): 0.00%
Time (mean ± σ): 5.854 ms ± 2.264 ms ┊ GC (mean ± σ): 6.29% ± 11.49%
█▆
██▇▁▄▁▆▇▄▁▁▄▄▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▆▄▇ ▇
5.26 ms Histogram: log(frequency) by time 18.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 6.91080085e-310
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: 747 samples with 1 evaluation.
Range (min … max): 6.157 ms … 23.199 ms ┊ GC (min … max): 0.00% … 62.42%
Time (median): 6.257 ms ┊ GC (median): 0.00%
Time (mean ± σ): 6.694 ms ± 2.380 ms ┊ GC (mean ± σ): 5.71% ± 11.01%
█▃
██▄▅▃▁▁▁▁▅▃▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆▇ ▆
6.16 ms Histogram: log(frequency) by time 20.8 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: 746 samples with 1 evaluation.
Range (min … max): 6.138 ms … 22.708 ms ┊ GC (min … max): 0.00% … 68.18%
Time (median): 6.275 ms ┊ GC (median): 0.00%
Time (mean ± σ): 6.702 ms ± 2.339 ms ┊ GC (mean ± σ): 5.65% ± 11.04%
█▃
██▁▃▃▃▁▁▃▁▁▄▁▅▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆▆ ▆
6.14 ms Histogram: log(frequency) by time 20.7 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.001007
2 1.588980e-03 8.189120e-04 0.002072
3 3.722710e-05 4.379408e-07 0.000917
4 2.043880e-08 2.294999e-10 0.000892
5 2.627965e-13 6.925509e-17 0.000891
6 7.415082e-18 1.144852e-26 0.000882
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))