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: 527 samples with 1 evaluation.
Range (min … max): 8.582 ms … 27.837 ms ┊ GC (min … max): 0.00% … 58.31%
Time (median): 8.721 ms ┊ GC (median): 0.00%
Time (mean ± σ): 9.494 ms ± 3.314 ms ┊ GC (mean ± σ): 6.77% ± 12.39%
█▄
██▅▅▁▁▁▁▁▁▄▁▅▄▅▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▇▆ ▆
8.58 ms Histogram: log(frequency) by time 25.6 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: 865 samples with 1 evaluation.
Range (min … max): 5.160 ms … 25.331 ms ┊ GC (min … max): 0.00% … 59.49%
Time (median): 5.250 ms ┊ GC (median): 0.00%
Time (mean ± σ): 5.778 ms ± 2.656 ms ┊ GC (mean ± σ): 7.80% ± 12.38%
█▃
██▅▁▁▁▃▃▁▅▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▆▇ ▇
5.16 ms Histogram: log(frequency) by time 20.6 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 8.03326739e-01
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: 741 samples with 1 evaluation.
Range (min … max): 6.197 ms … 23.149 ms ┊ GC (min … max): 0.00% … 70.57%
Time (median): 6.311 ms ┊ GC (median): 0.00%
Time (mean ± σ): 6.771 ms ± 2.585 ms ┊ GC (mean ± σ): 6.17% ± 11.29%
█▂
██▇▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆▄▃▃▅ ▆
6.2 ms Histogram: log(frequency) by time 22.7 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: 729 samples with 1 evaluation.
Range (min … max): 6.218 ms … 30.534 ms ┊ GC (min … max): 0.00% … 55.81%
Time (median): 6.328 ms ┊ GC (median): 0.00%
Time (mean ± σ): 6.856 ms ± 2.724 ms ┊ GC (mean ± σ): 6.07% ± 11.05%
█▂
██▄▅▃▁▃▁▃▃▃▇▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▅▄▁▃▅ ▆
6.22 ms Histogram: log(frequency) by time 23.2 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.000977
2 1.588980e-03 8.189120e-04 0.002111
3 3.722710e-05 4.379408e-07 0.000897
4 2.043880e-08 2.294999e-10 0.000873
5 2.627965e-13 6.925509e-17 0.000878
6 7.415082e-18 1.144852e-26 0.000883
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))