diff --git a/docs/make.jl b/docs/make.jl index 549265b3..ba378fc8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -23,12 +23,12 @@ makedocs( "tutorial-basic-example-f.md", "tutorial-double-integrator.md", "tutorial-goddard.md", + "tutorial-iss.md", "tutorial-init.md", "tutorial-lqr-basic.md", "tutorial-plot.md", #"tutorial-model.md", #"tutorial-solvers.md", - #"tutorial-iss.md", #"tutorial-flows.md", #"tutorial-ct_repl.md", ], diff --git a/docs/src/juliaopt2023.md b/docs/src/juliaopt2023.md new file mode 100644 index 00000000..8bd62f82 --- /dev/null +++ b/docs/src/juliaopt2023.md @@ -0,0 +1,64 @@ +```@raw html +juliacon +``` + +# Solving optimal control problems with Julia +### [Jean-Baptiste Caillau](http://caillau.perso.math.cnrs.fr), [Olivier Cots](https://ocots.github.io), [Joseph Gergaud](https://scholar.google.com/citations?user=pkH4An4AAAAJ&hl=fr), [Pierre Martinon](https://www.linkedin.com/in/pierre-martinon-b4603a17), [Sophia Sed](https://iww.inria.fr/sed-sophia) + +```@raw html +affiliations +``` + +# What it's about +- Nonlinear optimal control of ODEs: +```math +g(x(t_0),x(t_f)) + \int_{t_0}^{t_f} f^0(x(t), u(t)) \to \min +``` +subject to +```math +\dot{x}(t) = f(x(t), u(t)),\quad t \in [0, t_f] +``` +plus boundary, control and state constraints +- Our core interests: numerical & geometrical methods in control, applications + +# Where it comes from +- [BOCOP: the optimal control solver](https://www.bocop.org) +- [HamPath: indirect and Hamiltonian pathfollowing](http://www.hampath.org) +- [Coupling direct and indirect solvers, examples](https://ct.gitlabpages.inria.fr/gallery//notebooks.html) + +# OptimalControl.jl +- [Basic example: double integrator](https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-basic-example-f.html) +- [Basic example: double integrator (cont'ed)](https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-basic-example.html) +- [Advanced example: Goddard problem](https://control-toolbox.org/docs/optimalcontrol/dev/tutorial-goddard.html) + +# Wrap up +- [X] High level modelling of optimal control problems +- [X] Efficient numerical resolution coupling direct and indirect methods +- [X] Collection of examples + +# Future +- [ct_repl](./assets/repl.mp4) +- Additional solvers: direct shooting, collocation for BVP, Hamiltonian pathfollowing... +- ... and open to contributions! +- [CTProblems.jl](https://control-toolbox.org/docs/ctproblems/stable/problems-list.html) + +# control-toolbox.org +- Open toolbox +- Collection of Julia Packages rooted at [OptimalControl.jl](https://control-toolbox.org/docs/optimalcontrol) + +```@raw html +control-toolbox.org +``` + +# Credits (not exhaustive!) +- [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) +- [JuMP](https://jump.dev/JuMP.jl), + [InfiniteOpt.jl](https://docs.juliahub.com/InfiniteOpt/p3GvY/0.4.1), + [ADNLPModels.jl](https://jso.dev/ADNLPModels.jl) +- [Ipopt](https://github.com/coin-or/ipopt) +- [JuliaDiff](https://juliadiff.org) + ([FowardDiff.jl](https://juliadiff.org/ForwardDiff.jl), + [Zygote.jl](https://fluxml.ai/Zygote.jl)) +- [MLStyle.jl](https://thautwarm.github.io/MLStyle.jl) +- [REPLMaker.jl](https://docs.juliahub.com/ReplMaker) + diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md index e1dd5cc6..0665aabe 100644 --- a/docs/src/tutorial-goddard.md +++ b/docs/src/tutorial-goddard.md @@ -277,9 +277,11 @@ using MINPACK # NLE solver nle = (s, ξ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7]) # auxiliary function # with aggregated inputs ξ = [ p0 ; t1 ; t2 ; t3 ; tf ] # initial guess -indirect_sol = @suppress_err begin # hide -fsolve(nle, ξ) -end # hide +global indirect_sol = # hide +@suppress_err begin # hide +fsolve(nle, ξ) # hide +indirect_sol = fsolve(nle, ξ) # resolution of S(ξ) = 0 +end # hide # we retrieve the costate solution together with the times p0 = indirect_sol.x[1:3] diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md index 063647c5..501d93a0 100644 --- a/docs/src/tutorial-iss.md +++ b/docs/src/tutorial-iss.md @@ -1,3 +1,132 @@ # Indirect simple shooting -In construction. +In this tutorial we present the indirect simple shooting method on a simple example. + +```@setup main +using Suppressor # to suppress warnings +``` + +Let us start by importing the necessary packages. + +```@example main +using OptimalControl +using MINPACK # NLE solver +``` + +Let us consider the following optimal control problem. + +```@example main +t0 = 0 +tf = 1 +x0 = -1 +xf = 0 +α = 1.5 +@def ocp 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; +nothing # hide +``` + +From the Pontryagin Maximum Principle, we have that the maximising control is given by + +```math +u(x, p) = p +``` + +and we can define the flow of the associated Hamiltonian vector field + +```@example main +u(x, p) = p + +f = Flow(ocp, u) +nothing # hide +``` + +We define an auxiliary exponential map. + +```@example main +exp(p0; saveat=[]) = f((t0, tf), x0, p0, saveat=saveat).ode_sol +nothing # hide +``` + +We are now in position to define the shooting function and solve the shooting equations. + +```@example main +S(p0) = exp(p0)(tf)[1] - xf; # shooting function + +nle = (s, ξ) -> s[1] = S(ξ[1]) # auxiliary function +ξ = [ 0.0 ] # initial guess + +global indirect_sol = # hide +@suppress_err begin # hide +fsolve(nle, ξ) # hide +indirect_sol = fsolve(nle, ξ) # resolution of S(p0) = 0 +end # hide + +p0_sol = indirect_sol.x[1] # costate solution +println("costate: p0 = ", p0_sol) +@suppress_err begin # hide +println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n") +end # hide +nothing # hide +``` + +We get: + +```@example main +# parameters +times = range(t0, tf, length=2) +p0min = -0.5 +p0max = 2 + +# plot of the flow +plt_flow = plot() + +p0s = range(p0min, p0max, length=20) +for i ∈ 1:length(p0s) # plot some trajectories + sol = exp(p0s[i]) + x = [sol(t)[1] for t ∈ sol.t] + p = [sol(t)[2] for t ∈ sol.t] + label = i==1 ? "extremals" : false + plot!(plt_flow, x, p, color=:blue, label=label) +end + +p0s = range(p0min, p0max, length=200) +xs = zeros(length(p0s), length(times)) +ps = zeros(length(p0s), length(times)) +for i ∈ 1:length(p0s) + sol = exp(p0s[i], saveat=times) + xs[i, :] = [z[1] for z ∈ sol.(times)] + ps[i, :] = [z[2] for z ∈ sol.(times)] +end +for j ∈ 1:length(times) # plot intermediate points + 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 = (p0min, p0max)) +plot!(plt_flow, [0, 0], [p0min, p0max], color=:black, xlabel="x", ylabel="p", label="x=xf") + +# solution +sol = exp(p0_sol) +x = [sol(t)[1] for t ∈ sol.t] +p = [sol(t)[2] for t ∈ sol.t] +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) # solution + +# plot of the shooting function +plt_shoot = plot(xlims=(p0min, p0max), ylims=(-2, 4), xlabel="p₀", ylabel="y") +plot!(plt_shoot, p0s, S, linewidth=2, label="S(p₀)", color=:green) +plot!(plt_shoot, [p0min, p0max], [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) # solution + +# plots +plot(plt_flow, plt_shoot, layout=(1,2), size=(800, 450)) +```