From b254c077c39c054228b95070b5e2a784e12a2f94 Mon Sep 17 00:00:00 2001 From: Olivier Cots Date: Tue, 13 Aug 2024 14:04:55 +0200 Subject: [PATCH 1/2] add tuto solve --- docs/make.jl | 3 +- docs/src/tutorial-flow.md | 2 +- docs/src/tutorial-initial-guess.md | 2 +- docs/src/tutorial-iss.md | 38 +++--- docs/src/tutorial-nlp.md | 7 +- docs/src/tutorial-solve.md | 185 +++++++++++++++++++++++++++++ 6 files changed, 212 insertions(+), 25 deletions(-) create mode 100644 docs/src/tutorial-solve.md diff --git a/docs/make.jl b/docs/make.jl index 27c2057c..96f0bd68 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -37,9 +37,10 @@ makedocs( ], "Manual" => [ "Abstract syntax" => "tutorial-abstract.md", - "Flow" => "tutorial-flow.md", + "Solve" => "tutorial-solve.md", "Initial guess" => "tutorial-initial-guess.md", "Plot a solution" => "tutorial-plot.md", + "Flow" => "tutorial-flow.md", ], "Tutorials" => [ "tutorial-continuation.md", diff --git a/docs/src/tutorial-flow.md b/docs/src/tutorial-flow.md index d6cace1f..a45c142e 100644 --- a/docs/src/tutorial-flow.md +++ b/docs/src/tutorial-flow.md @@ -1,6 +1,6 @@ # [How to compute flows](@id manual-flow) -In this tutorial, we explain the `Flow` function from `OptimalControl` package. +In this tutorial, we explain the `Flow` function from `OptimalControl.jl` package. ## Basic usage diff --git a/docs/src/tutorial-initial-guess.md b/docs/src/tutorial-initial-guess.md index 01967f97..114f594f 100644 --- a/docs/src/tutorial-initial-guess.md +++ b/docs/src/tutorial-initial-guess.md @@ -1,4 +1,4 @@ -# Initial guess (or iterate) for the resolution +# [Initial guess (or iterate) for the resolution](@id tutorial-init) ```@meta CurrentModule = OptimalControl diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md index e17a14f4..17cfa627 100644 --- a/docs/src/tutorial-iss.md +++ b/docs/src/tutorial-iss.md @@ -5,11 +5,11 @@ In this tutorial we present the indirect simple shooting method on a simple exam Let us start by importing the necessary packages. ```@example main -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 solution +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 solution ``` ## Optimal control problem @@ -129,7 +129,7 @@ nothing # hide ``` 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**. +Now, to solve the (BVP) we introduce the **shooting function**: ```math \begin{array}{rlll} @@ -139,7 +139,7 @@ Now, to solve the (BVP) we introduce the **shooting function**. ``` ```@example main -S(p0) = π( φ(t0, x0, p0, tf) ) - xf # shooting function +S(p0) = π( φ(t0, x0, p0, tf) ) - xf # shooting function nothing # hide ``` @@ -149,7 +149,7 @@ At the end, solving (BVP) is equivalent to solve $S(p_0) = 0$. This is what we c **indirect simple shooting method**. We define an initial guess. ```@example main -ξ = [ 0.0 ] # initial guess +ξ = [ 0.1 ] # initial guess nothing # hide ``` @@ -159,8 +159,8 @@ We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) equation. Let us define the problem. ```@example main -nle! = (s, ξ, λ) -> s[1] = S(ξ[1]) # auxiliary function -prob = NonlinearProblem(nle!, ξ) # NLE problem with initial guess +nle! = (s, ξ, λ) -> s[1] = S(ξ[1]) # auxiliary function +prob = NonlinearProblem(nle!, ξ) # NLE problem with initial guess nothing # hide ``` @@ -181,8 +181,8 @@ For small nonlinear systems, it could be faster to use the Now, let us solve the problem and retrieve the initial costate solution. ```@example main -indirect_sol = solve(prob; show_trace=Val(true)) # resolution of S(p0) = 0 -p0_sol = indirect_sol.u[1] # 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") ``` @@ -219,8 +219,8 @@ nothing # hide Let us define the problem to solve. ```@example main -nle! = ( s, ξ) -> s[1] = S(ξ[1]) # auxiliary function -jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ) # Jacobian of nle +nle! = ( s, ξ) -> s[1] = S(ξ[1]) # auxiliary function +jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ) # Jacobian of nle nothing # hide ``` @@ -228,7 +228,7 @@ We are now in position to solve the problem with the `hybrj` solver from `MINPAC function, providing the Jacobian. Let us do some benchmarking. ```@example main -@benchmark fsolve(nle!, jnle!, ξ; show_trace=false) # initial guess given to the solver +@benchmark fsolve(nle!, jnle!, ξ; show_trace=false) # initial guess given to the solver ``` We can also use the [preparation step](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/tutorial1/#Preparing-for-multiple-gradients) of `DifferentiationInterface.jl`. @@ -242,8 +242,8 @@ jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras) Now, let us solve the problem and retrieve the initial costate solution. ```@example main -indirect_sol = fsolve(nle!, jnle!, ξ; show_trace=true) # resolution of S(p0) = 0 -p0_sol = indirect_sol.x[1] # 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") ``` @@ -253,8 +253,8 @@ println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n") The solution can be plot calling first the flow. ```@example main -sol = φ( (t0, tf), x0, p0_sol ) -plot(sol, size=(800, 600)) +sol = φ((t0, tf), x0, p0_sol) +plot(sol) ``` In the indirect shooting method, the research of the optimal control is replaced by the computation diff --git a/docs/src/tutorial-nlp.md b/docs/src/tutorial-nlp.md index 2c67834d..baecabdb 100644 --- a/docs/src/tutorial-nlp.md +++ b/docs/src/tutorial-nlp.md @@ -61,6 +61,7 @@ For a first example we use the `ipopt` solver from [`NLPModelsIpopt.jl`](https:/ ```@example main using NLPModelsIpopt + nlp_sol = ipopt(nlp; print_level=5, mu_strategy="adaptive", tol=1e-8, sb="yes") nothing # hide ``` @@ -68,7 +69,7 @@ nothing # hide Then we can rebuild and plot an optimal control problem solution (note that the multipliers are optional, but the OCP costate will not be retrieved if the multipliers are not provided). ```@example main -sol = OptimalControlSolution(docp, primal=nlp_sol.solution, dual=nlp_sol.multipliers) +sol = OptimalControlSolution(docp; primal=nlp_sol.solution, dual=nlp_sol.multipliers) plot(sol) ``` ## Change the NLP solver @@ -86,7 +87,7 @@ Another possible NLP solver is [`Percival.jl`](https://jso.dev/Percival.jl). ```@example main using Percival -nlp_sol = percival(nlp, verbose=1) +nlp_sol = percival(nlp; verbose=1) ``` ## Initial guess @@ -94,7 +95,7 @@ nlp_sol = percival(nlp, verbose=1) An initial guess, including warm start, can be passed to [`direct_transcription`](@ref) the same way as for `solve`. ```@example main -docp, nlp = direct_transcription(ocp, init=sol) +docp, nlp = direct_transcription(ocp; init=sol) nothing # hide ``` diff --git a/docs/src/tutorial-solve.md b/docs/src/tutorial-solve.md new file mode 100644 index 00000000..7ea05647 --- /dev/null +++ b/docs/src/tutorial-solve.md @@ -0,0 +1,185 @@ +# [The solve function](@id manual-solve) + +In this tutorial, we explain the [`solve`](@ref) function from `OptimalControl.jl` package. + +## Basic usage + +Les us define a basic optimal control problem. + +```@example main +using OptimalControl + +t0 = 0 +tf = 1 +x0 = [-1, 0] + +ocp = @def begin + + t ∈ [ t0, tf ], time + x = (q, v) ∈ R², state + u ∈ R, control + + x(t0) == x0 + x(tf) == [ 0, 0 ] + + ẋ(t) == [ v(t), u(t) ] + + ∫( 0.5u(t)^2 ) → min + +end +nothing # hide +``` + +Let us try to solve the problem: + +```@setup main_repl +using OptimalControl +t0 = 0 +tf = 1 +x0 = [-1, 0] +ocp = @def begin + t ∈ [ t0, tf ], time + x = (q, v) ∈ R², state + u ∈ R, control + x(t0) == x0 + x(tf) == [ 0, 0 ] + ẋ(t) == [ v(t), u(t) ] + ∫( 0.5u(t)^2 ) → min +end +``` + +```@repl main_repl +solve(ocp) +``` + +As you can see, an error occured since we need the package [`NLPModelsIpopt.jl`](https://jso.dev/NLPModelsIpopt.jl). + +Actually, the default solving method is what we call a +[direct method](https://en.wikipedia.org/wiki/Optimal_control#Numerical_methods_for_optimal_control). +In a direct method, the optimal control problem is transcribed to a nonlinear optimization problem (NLP) of the form + +```math +\text{minimize}\quad F(y), \quad\text{subject to the constraints}\quad g(y)=0, \quad h(y)\le 0. +``` + +The `OptimalControl.jl` package makes the transcription but it needs a package to modelise the NLP problem and +another one to solve it. The `NLPModelsIpopt.jl` package provides an interface to the well-known solver +[`Ipopt`](https://coin-or.github.io/Ipopt/) that can be used to solve general nonlinear programming problems. + +```@example main +using NLPModelsIpopt + +solve(ocp) +nothing # hide +``` + +## Solvers + +`OptimalControl.jl` offers a list of methods to solve your optimal control problem. To get the list of methods, +simply call `available_methods`. + +```@example main +available_methods() +``` + +Each line is a method. The priority is given from top to bottom. This means that + +```julia +solve(ocp) +``` + +is equivalent to + +```julia +solve(ocp, :direct, :adnlp, :ipopt) +``` + +Let us detail the three symbols. As you can see, there are only direct methods. The symbol `:adnlp` is for the +choice of modeler. As said before, the NLP problem needs to be modelised in `Julia` code. We use +[`ADNLPModels.jl`](https://jso.dev/ADNLPModels.jl) which provides automatic differentiation (AD)-based model implementations that conform to the [`NLPModels.jl`](https://github.com/JuliaSmoothOptimizers/ADNLPModels.jl) API. + +The last symbol is what distinguish the two available methods. It corresponds to the NLP solver. By default, we +use the `NLPModelsIpopt.jl` package but you can choose [`MadNLP.jl`](https://madnlp.github.io/MadNLP.jl) which +is an open-source nonlinear programming solver, purely implemented in Julia and which implements a filter +line-search interior-point algorithm, as used in Ipopt. + +Note that you are not compelled to provide the full description of the method to use it. A partial description, +if not ambiguous, will work. Just remember that if several full descriptions contain the partial one, then, the +priority is given to first one in the list. Hence, these calls are all equivalent: + +```julia +solve(ocp) +solve(ocp, :direct ) +solve(ocp, :adnlp ) +solve(ocp, :ipopt) +solve(ocp, :direct, :adnlp ) +solve(ocp, :direct, :ipopt) +solve(ocp, :direct, :adnlp, :ipopt) +``` + +Let us try `MadNLP.jl`. + +```@example main +solve(ocp, :direct, :adnlp, :madnlp; display=true) +nothing # hide +``` + +Again, there are several equivalent manners to use `MadNLP.jl`. + +```julia +solve(ocp, :madnlp) +solve(ocp, :direct, :madnlp) +solve(ocp, :adnlp, :madnlp) +solve(ocp, :direct, :adnlp, :madnlp) +``` + +## Options + +### Direct method + +The options common to all the direct methods may be seen directly in the +[code](https://github.com/control-toolbox/CTDirect.jl/blob/23c76374fb91f65f8d10c03b7d8ff5c75a5f4076/src/solve.jl#L68). +There are `init`, `grid_size` and `time_grid`. + +- The `init` option can be used to set an initial guess for the solver. See this [tutorial](@ref tutorial-init). +- The `grid_size` option corresponds to the size of the (uniform) time discretization grid. +More precisely, it is the number of steps, that is if `N = grid_size` and if the initial and final times are +denoted respectively `t0` and `tf`, then we have: +```julia +Δt = (tf - t0) / N +``` +- The `time_grid` option is the grid of times: `t0, t1, ..., tf`. If the initial and/or the final times are free, then you can provide a normalised grid between 0 and 1. Note that you can set either `grid_size` or `time_grid` but not both. + +```@example main +sol = solve(ocp; grid_size=10, display=false) +sol.times +``` + +Or with `MadNLP`: + +```@example main +sol = solve(ocp, :madnlp; grid_size=10, display=false) +sol.times +``` + +### NLPModelsIpopt + +You can provide any option of `NLPModelsIpopt` or `Ipopt` with a pair `keyword=value`. Please check +the list of [`Ipopt` options](https://coin-or.github.io/Ipopt/OPTIONS.html) and the +[`NLPModelsIpopt.jl` documentation](https://jso.dev/NLPModelsIpopt.jl). + +```@example main +solve(ocp, :ipopt; max_iter=0) +nothing # hide +``` + +### MadNLP + +If you use the `MadNLP` solver, then you can provide any option of it. Please check the +[`MadNLP.jl` documentation](https://madnlp.github.io/MadNLP.jl) and the list of +[`MadNLP.jl`options](https://madnlp.github.io/MadNLP.jl/stable/options/). + +```@example main +solve(ocp, :madnlp; max_iter=1, display=true) +nothing # hide +``` \ No newline at end of file From e62d0f3bdbf93a32a5342b90c146d3cebe4d30f1 Mon Sep 17 00:00:00 2001 From: Olivier Cots <66357348+ocots@users.noreply.github.com> Date: Tue, 13 Aug 2024 14:14:01 +0200 Subject: [PATCH 2/2] Update tutorial-solve.md --- docs/src/tutorial-solve.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/src/tutorial-solve.md b/docs/src/tutorial-solve.md index 7ea05647..66529d0c 100644 --- a/docs/src/tutorial-solve.md +++ b/docs/src/tutorial-solve.md @@ -120,6 +120,8 @@ solve(ocp, :direct, :adnlp, :ipopt) Let us try `MadNLP.jl`. ```@example main +using MadNLP + solve(ocp, :direct, :adnlp, :madnlp; display=true) nothing # hide ``` @@ -182,4 +184,4 @@ If you use the `MadNLP` solver, then you can provide any option of it. Please ch ```@example main solve(ocp, :madnlp; max_iter=1, display=true) nothing # hide -``` \ No newline at end of file +```