diff --git a/docs/Project.toml b/docs/Project.toml index 022ff0565..4fab10d5f 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -17,6 +17,7 @@ MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -27,6 +28,8 @@ BenchmarkTools = "1.6" CTBase = "0.15" CTDirect = "0.14" CTFlows = "0.8" +CTModels = "0.2" +CTParser = "0.2" CommonSolve = "0.2" DifferentiationInterface = "0.6" Documenter = "1.8" @@ -34,11 +37,12 @@ DocumenterMermaid = "0.2" ForwardDiff = "0.10" JLD2 = "0.5" JSON3 = "1.14" +LinearAlgebra = "1" MINPACK = "1.3" MadNLP = "0.8" NLPModelsIpopt = "0.10" -NonlinearSolve = "4.4" -OrdinaryDiffEq = "6" +NonlinearSolve = "4.6" +OrdinaryDiffEq = "6.93" Percival = "0.7" Plots = "1.40" Suppressor = "0.2" diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md index a301d82dd..5e862c5e5 100644 --- a/docs/src/tutorial-goddard.md +++ b/docs/src/tutorial-goddard.md @@ -56,12 +56,12 @@ using Plots # to plot the solution We define the problem ```@example main -t0 = 0 # initial time -r0 = 1 # initial altitude -v0 = 0 # initial speed -m0 = 1 # initial mass -vmax = 0.1 # maximal authorized speed -mf = 0.6 # final mass to target +const t0 = 0 # initial time +const r0 = 1 # initial altitude +const v0 = 0 # initial speed +const m0 = 1 # initial mass +const vmax = 0.1 # maximal authorized speed +const mf = 0.6 # final mass to target ocp = @def begin # definition of the optimal control problem @@ -70,7 +70,7 @@ ocp = @def begin # definition of the optimal control problem x = (r, v, m) ∈ R³, state u ∈ R, control - x(t0) == [ r0, v0, m0 ] + x(t0) == [r0, v0, m0] m(tf) == mf, (1) 0 ≤ u(t) ≤ 1 r(t) ≥ r0 @@ -80,7 +80,7 @@ ocp = @def begin # definition of the optimal control problem r(tf) → max -end; +end # Dynamics const Cd = 310 @@ -91,12 +91,12 @@ const b = 2 F0(x) = begin r, v, m = x D = Cd * v^2 * exp(-β*(r - 1)) # Drag force - return [ v, -D/m - 1/r^2, 0 ] + return [v, -D/m - 1/r^2, 0] end F1(x) = begin r, v, m = x - return [ 0, Tmax/m, -b*Tmax ] + return [0, Tmax/m, -b*Tmax] end nothing # hide ``` @@ -203,13 +203,13 @@ these expressions are straightforwardly translated into Julia code: ```@example main # Controls -u0 = 0 # off control -u1 = 1 # bang control +const u0 = 0 # off control +const u1 = 1 # bang control H0 = Lift(F0) # H0(x, p) = p' * F0(x) -H01 = @Lie { H0, H1 } -H001 = @Lie { H0, H01 } -H101 = @Lie { H1, H01 } +H01 = @Lie {H0, H1} +H001 = @Lie {H0, H01} +H101 = @Lie {H1, H01} us(x, p) = -H001(x, p) / H101(x, p) # singular control ub(x) = -(F0⋅g)(x) / (F1⋅g)(x) # boundary control @@ -229,7 +229,7 @@ Then, we define the shooting function according to the optimal structure we have that is a concatenation of four arcs. ```@example main -x0 = [ r0, v0, m0 ] # initial state +x0 = [r0, v0, m0] # initial state function shoot!(s, p0, t1, t2, t3, tf) @@ -239,7 +239,7 @@ function shoot!(s, p0, t1, t2, t3, tf) xf, pf = f0(t3, x3, p3, tf) s[1] = xf[3] - mf # final mass constraint - s[2:3] = pf[1:2] - [ 1, 0 ] # transversality conditions + s[2:3] = pf[1:2] - [1, 0] # transversality conditions s[4] = H1(x1, p1) # H1 = H01 = 0 s[5] = H01(x1, p1) # at the entrance of the singular arc s[6] = g(x2) # g = 0 when entering the boundary arc @@ -283,72 +283,23 @@ println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n") We aggregate the data to define the initial guess vector. ```@example main -ξ = [ p0 ; t1 ; t2 ; t3 ; tf ] # initial guess -``` - -### NonlinearSolve.jl - -We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package to solve the shooting -equation. Let us define the problem. - -```@example main -# auxiliary function with aggregated inputs -nle! = (s, ξ, λ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7]) - -# NLE problem with initial guess -prob = NonlinearProblem(nle!, ξ) -nothing # hide -``` - -Let us do some benchmarking. - -```@example main -using BenchmarkTools -@benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false)) -``` - -For small nonlinear systems, it could be faster to use the -[`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/). - -```@example main -@benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false)) -``` - -Now, let us solve the problem and retrieve the initial costate solution. - -```@example main -# resolution of S(ξ) = 0 -indirect_sol = solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(true)) - -# we retrieve the costate solution together with the times -p0 = indirect_sol.u[1:3] -t1 = indirect_sol.u[4] -t2 = indirect_sol.u[5] -t3 = indirect_sol.u[6] -tf = indirect_sol.u[7] - -println("") -println("p0 = ", p0) -println("t1 = ", t1) -println("t2 = ", t2) -println("t3 = ", t3) -println("tf = ", tf) - -# Norm of the shooting function at solution -s = similar(p0, 7) -shoot!(s, p0, t1, t2, t3, tf) -println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n") +ξ = [p0..., t1, t2, t3, tf] # initial guess ``` ### MINPACK.jl +We can use [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package or, instead, the +[MINPACK.jl](https://github.com/sglyon/MINPACK.jl) package to solve +the shooting equation. To compute the Jacobian of the shooting function we use the +[DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with +[ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend. + ```@setup main -using MINPACK function fsolve(f, j, x; kwargs...) try MINPACK.fsolve(f, j, x; kwargs...) catch e - println("Erreur using MINPACK") + println("Error using MINPACK") println(e) println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.") MINPACK.fsolve(f, x; kwargs...) @@ -356,12 +307,6 @@ function fsolve(f, j, x; kwargs...) end ``` -Instead of the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package we can use the -[MINPACK.jl](https://github.com/sglyon/MINPACK.jl) package to solve -the shooting equation. To compute the Jacobian of the shooting function we use the -[DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with -[ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend. - ```@example main using DifferentiationInterface import ForwardDiff @@ -381,21 +326,8 @@ nothing # hide ``` We 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. - -```@example main -@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. - -```@example main -extras = prepare_jacobian(nle!, similar(ξ), backend, ξ) -jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras) -@benchmark fsolve(nle!, jnle_prepared!, ξ; show_trace=false) -``` - -Now, let us solve the problem and retrieve the initial costate solution. +function, providing the Jacobian. +Let us solve the problem and retrieve the initial costate solution. ```@example main # resolution of S(ξ) = 0 diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md index 251a0c814..d67bcc3f9 100644 --- a/docs/src/tutorial-iss.md +++ b/docs/src/tutorial-iss.md @@ -30,11 +30,11 @@ Let us consider the following optimal control problem: 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$. ```@example main -t0 = 0 -tf = 1 -x0 = -1 -xf = 0 -α = 1.5 +const t0 = 0 +const tf = 1 +const x0 = -1 +const xf = 0 +const α = 1.5 ocp = @def begin t ∈ [t0, tf], time @@ -145,50 +145,17 @@ nothing # hide ## Resolution 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. +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. ```@example main -ξ = [ 0.1 ] # initial guess +ξ = [0.1] # initial guess nothing # hide ``` -### NonlinearSolve.jl - -We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package to solve the shooting -equation. Let us define the problem. - -```@example main -nle! = (s, ξ, λ) -> s[1] = S(ξ[1]) # auxiliary function -prob = NonlinearProblem(nle!, ξ) # NLE problem with initial guess -nothing # hide -``` - -Let us do some benchmarking. - -```@example main -using BenchmarkTools -@benchmark solve(prob; show_trace=Val(false)) -``` - -For small nonlinear systems, it could be faster to use the -[`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/). - -```@example main -@benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false)) -``` - -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 -println("\ncostate: p0 = ", p0_sol) -println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n") -``` - ### MINPACK.jl +We can use [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package or, instead, [MINPACK.jl](https://github.com/sglyon/MINPACK.jl) to solve the shooting equation. To compute the Jacobian of the shooting function we use [DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) with [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend. + ```@setup main using MINPACK function fsolve(f, j, x; kwargs...) @@ -203,12 +170,6 @@ function fsolve(f, j, x; kwargs...) end ``` -Instead of the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package we can use the -[MINPACK.jl](https://github.com/sglyon/MINPACK.jl) package to solve -the shooting equation. To compute the Jacobian of the shooting function we use the -[DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with -[ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend. - ```@example main using DifferentiationInterface import ForwardDiff @@ -224,22 +185,8 @@ jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ) # Jacobian nothing # hide ``` -We 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. - -```@example main -@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. - -```@example main -extras = prepare_jacobian(nle!, similar(ξ), backend, ξ) -jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras) -@benchmark fsolve(nle!, jnle_prepared!, ξ; show_trace=false) -``` - -Now, let us solve the problem and retrieve the initial costate solution. +We are now in position to solve the problem with the `hybrj` solver from MINPACK.jl through the `fsolve` function, providing the Jacobian. +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 @@ -257,10 +204,7 @@ 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. +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. ```@raw html
@@ -408,5 +352,3 @@ nothing # hide ```@example main pretty_plot(S, p0_sol; size=(800, 450)) ``` - - diff --git a/test/Project.toml b/test/Project.toml index 199f6cf56..a092c2c2b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -10,7 +10,13 @@ SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] +CTModels = "0.2" +DifferentiationInterface = "0.6" +ForwardDiff = "0.10" +LinearAlgebra = "1" +MINPACK = "1.3" NLPModelsIpopt = "0.10" -OrdinaryDiffEq = "6" +OrdinaryDiffEq = "6.93" SplitApplyCombine = "1.2" +Test = "1" julia = "1.10"