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: 504 samples with 1 evaluation.
Range (min … max): 8.144 ms … 104.291 ms ┊ GC (min … max): 0.00% … 90.84%
Time (median): 9.412 ms ┊ GC (median): 0.00%
Time (mean ± σ): 9.917 ms ± 7.068 ms ┊ GC (mean ± σ): 6.80% ± 8.63%
▆▅▄▄▃▃▂▂▂ ▁ █▆▄▁ ▂▂ ▅▃▁
█████████▅██▇█▄▅▄▄▆▄▄▄▁▁▇████▅▄▄▄▄▆▄▅▁▁▁▁▁▁▁▄▄▁▁▄▁▁██████▇▄ █
8.14 ms Histogram: log(frequency) by time 11.1 ms <
Memory estimate: 6.38 MiB, allocs estimate: 164632.For small nonlinear systems, it could be faster to use the SimpleNewtonRaphson() descent algorithm.
@benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false))BenchmarkTools.Trial: 756 samples with 1 evaluation.
Range (min … max): 5.096 ms … 95.017 ms ┊ GC (min … max): 0.00% … 94.11%
Time (median): 6.443 ms ┊ GC (median): 0.00%
Time (mean ± σ): 6.642 ms ± 5.688 ms ┊ GC (mean ± σ): 9.25% ± 10.50%
▄ █
█▅▃▂█▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂ ▂
5.1 ms Histogram: frequency by time 26.2 ms <
Memory estimate: 4.72 MiB, allocs estimate: 101935.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 9.88131292e-324
1 1.58897952e-03 2.86166386e-02
2 8.72173220e-07 6.46258281e-04
3 2.62695417e-13 3.55113667e-07
Final 2.62695417e-13
----------------------
costate: p0 = 0.07202997482156911
shoot: |S(p0)| = 2.626954170931736e-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: 680 samples with 1 evaluation.
Range (min … max): 6.103 ms … 105.341 ms ┊ GC (min … max): 0.00% … 90.42%
Time (median): 6.340 ms ┊ GC (median): 0.00%
Time (mean ± σ): 7.357 ms ± 6.367 ms ┊ GC (mean ± σ): 6.35% ± 6.95%
█▄▁ ▁▁
▆████▅▆▄▄▄▃▃▂▃▂▁▂▂▁▂▂▁▂▁▂▁▁▁▁▁▁▁▂▁▁▂▁▁▁▁▁▁▁▁▂▄▇▇██▄▃▂▂▁▁▃▄▃ ▃
6.1 ms Histogram: frequency by time 8.34 ms <
Memory estimate: 4.30 MiB, allocs estimate: 122482.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: 677 samples with 1 evaluation.
Range (min … max): 6.088 ms … 103.710 ms ┊ GC (min … max): 0.00% … 93.40%
Time (median): 6.350 ms ┊ GC (median): 0.00%
Time (mean ± σ): 7.387 ms ± 6.217 ms ┊ GC (mean ± σ): 6.25% ± 6.98%
██▁ ▆
███▆▄▃▂▂▂▂▁▁▁▁▂▂▁▁▇█▆▂▂▅▁▂▁▁▁▁▁▁▁▁▂▁▁▁▂▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▂ ▃
6.09 ms Histogram: frequency by time 11.7 ms <
Memory estimate: 4.30 MiB, allocs estimate: 122476.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.001497
2 1.588980e-03 8.189120e-04 0.003233
3 3.722710e-05 4.379408e-07 0.001313
4 2.043880e-08 2.294999e-10 0.001294
5 2.627233e-13 6.925509e-17 0.001280
6 7.148331e-17 1.144258e-26 0.001258
costate: p0 = 0.07202997482167604
shoot: |S(p0)| = 7.14833118003657e-17Plot 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))