Add tutorial (goddard) in docs (#79)

* md for tutorial (goddard)

* moved to docs

* moved to docs

* try to build doc

* added goddard pic

* setup

* add Plots

* abstract simple ocp

* try abstract def for goddard

* abstarct goddardok, added to CI test

* tutorial.md in progress, goddard with warmstart ok but ipopt output does not appear; also need to fix output size error from api-ctbase.md

* got the stdout from ipopt with a nothing at block end

* proper footnotes

* ignore output size limit, CI build should pass now

* first version for tutorial.md; todo: update docstrings for api.md

docs/Project.toml

@@ -1,3 +1,4 @@
 [deps]
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

docs/make.jl

@@ -6,13 +6,13 @@ DocMeta.setdocmeta!(CTBase, :DocTestSetup, :(using CTBase); recursive = true)
 DocMeta.setdocmeta!(CTDirect, :DocTestSetup, :(using CTDirect); recursive = true)
 
 makedocs(
+    warnonly = [:cross_references, :autodocs_block],
     sitename = "CTDirect.jl",
-    format = Documenter.HTML(prettyurls = false),
+    format = Documenter.HTML(prettyurls = false,
+            size_threshold_ignore = ["api-ctbase.md"]),
     pages = [
         "Introduction"  => "index.md",
-        #"Methods and Options" => ["rk.md",
-        #                          "optimization.md",
-        #                          "solve-options.md"],
+        "Tutorial"      => "tutorial.md",
         "API"           => "api.md",
         #"Developers"    => "dev-api.md",
     ]

docs/src/api-ctbase.md

@@ -9,7 +9,6 @@ For more details about `CTBase.jl` package, see the [documentation](https://cont
 Pages   = ["api-ctbase.md"]
 Modules = [CTBase]
 Order = [:module, :constant, :type, :function, :macro]
-Private = false
 ```
 
 ## Documentation

docs/src/api.md

@@ -6,7 +6,6 @@
 Pages   = ["api.md"]
 Modules = [CTDirect]
 Order = [:module, :constant, :type, :function, :macro]
-Private = false
 ```
 
 ## Available methods

docs/src/tutorial.md

@@ -0,0 +1,129 @@
+# First example: Goddard problem
+
+![R.H. Goddard](./assets/goddard.png)
+
+We start with the well-known so-called [Goddard](http://en.wikipedia.org/wiki/Robert_H._Goddard) problem, which models the ascent of a rocket through the atmosphere (see for instance [^1],[^2]).
+We restrict here ourselves to  vertical (monodimensional) trajectories, and the state variables are the altitude, speed and mass of the rocket during the flight, for a total dimension of 3. The rocket is subject to gravity, thrust and drag forces. The final time is free, and the objective is here to reach a maximal altitude with a given fuel consumption, i.e a fixed final mass. We also impose a speed limit. All units are renormalized.
+
+[^1]: R.H. Goddard. A Method of Reaching Extreme Altitudes, volume 71(2) of Smithsonian Miscellaneous Collections. Smithsonian institution, City of Washington, 1919.
+
+[^2]: H. Seywald and E.M. Cliff. Goddard problem in presence of a dynamic pressure limit. Journal of Guidance, Control, and Dynamics, 16(4):776–781, 1993.
+
+## Problem definition
+
+First import the CTDirect and CTBase modules
+```@example main
+using CTDirect
+using CTBase
+```
+
+Then define the OCP for the Goddard problem. Note that the free final time is modeled as an *optimization variable*, and has both a lower bound to prevent the optimization getting stuck at tf=0. In this particular case an upper bound is not needed for the final time since the final mass is prescribed.
+```@example main
+Cd = 310
+β = 500
+Tmax = 3.5
+b = 2
+vmax = 0.1
+function F0(x)
+    r, v, m = x
+    D = Cd * v^2 * exp(-β*(r - 1))
+    return [ v, -D/m - 1/r^2, 0 ]
+end
+function F1(x)
+    r, v, m = x
+    return [ 0, Tmax/m, -b*Tmax ]
+end
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [ 0, tf ], time
+    x ∈ R^3, state
+    u ∈ R, control
+    0.1 ≤ tf ≤ Inf
+    r = x[1]
+    v = x[2]
+    m = x[3]
+    x(0) == [1, 0, 1]
+    m(tf) == 0.6
+    1 ≤ r(t) ≤ 1.1
+    0 ≤ v(t) ≤ vmax
+    0 ≤ u(t) ≤ 1
+    ẋ(t) == F0(x(t)) + u(t)*F1(x(t))
+    r(tf) → max
+end
+nothing # hide
+```
+
+## Basic solve
+
+We can solve the problem directly using the default options.
+```@example main
+sol1 = solve(ocp, print_level=5)
+nothing # hide
+```
+Then plot the solution with the state and control variables, as well as the costate recovered from the lagrange multipliers of the discretized problem. 
+```@example main
+plot(sol1)
+```
+
+The most common option for **solve** is the number of time steps for the discretized problem (default 100), that can be set with the argument *grid_size*. A larger grid size will increase the computational cost, while a smaller value may lead to a very coarse solution.
+
+## Initial guess options
+
+The function **solve** uses a default constant initialisation of 0.1 for all variables. More advanced options include constant and/or functional initialisation for each individual state or control component, as well as reusing an existing solution, also known as *warm start*[^3].
+
+[^3]: Currently only the primal variables are reused for the warm start, not the lagrange multipliers. It should be noted that previous experiments with the Bocop software seemed to indicate that initializing also the multipliers gave little benefit.
+
+Let us start with the simplest case, constant initialisation.
+```@example main
+x_const = [1.05, 0.2, 0.8]
+u_const = 0.5
+v_const = 0.15
+init1 = OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const)
+sol2 = solve(ocp, print_level=0, init=init1)
+println("Objective ", sol2.objective, " after ", sol2.iterations, " iterations")
+``` 
+
+Now we illustrate the functional initialisation, with some random functions. Note that we only consider the state and control variables, since the optimization variables are scalar and therefore a functional initialisation is not relevant. In the example notice that the call to **OptimalControlInit** does not provide an argument for the optimization variables, therefore the default initial guess will be used.  
+```@example main
+x_func = t->[1+t^2, sqrt(t), 1-t]
+u_func = t->(cos(t)+1)*0.5
+init2 = OptimalControlInit(x_init=x_func, u_init=u_func)
+sol3 = solve(ocp, print_level=0, init=init2)
+println("Objective ", sol3.objective, " after ", sol3.iterations, " iterations")
+``` 
+More generally, the default, constant and functional initialisations can be mixed, as shown in the example below that uses a functional initial guess for the state, a constant initial guess for the control, and the default initial guess for the optimization variables. 
+```@example main
+init3 = OptimalControlInit(x_init=x_func, u_init=u_const)
+sol4 = solve(ocp, print_level=0, init=init3)
+println("Objective ", sol4.objective, " after ", sol4.iterations, " iterations")
+``` 
+
+Finally, we can also use a so-called *warmstart* strategy and use an existing solution as initial guess (note that the OCP solution returned by the **solve** call is functional, thus it is not necessary to use the same time grid). Notice that the objective and constraint violation values start much closer to the solution than with the previous initialisations.
+```@example main
+init4 = OptimalControlInit(sol1)
+sol4 = solve(ocp, grid_size=200, print_level=5, init=init4)
+nothing # hide
+```
+```@example main
+plot(sol4)
+```
+
+## The discretized problem
+
+Instead of calling **solve** directly on the OCP problem, you can first obtain the discretized problem (DOCP) by calling **directTranscription**, then call **solve** on the DOCP.
+```@example main
+docp = directTranscription(ocp, grid_size=100)
+sol5 = solve(docp, print_level=5)
+nothing # hide
+```
+The initial guess can be passed to **solve** same as before.
+```@example main
+sol6 = solve(docp, print_level=0, init=OptimalControlInit(sol1))
+println("Objective ", sol6.objective, " after ", sol6.iterations, " iterations")
+```
+Another possibility is to set the initial guess associated to the DOCP, using the function **setDOCPInit**.
+```@example main
+setDOCPInit(docp, OptimalControlInit(sol1))
+sol7 = solve(docp, print_level=5)
+nothing # hide
+```

test/suite/abstract_ocp.jl

@@ -45,4 +45,42 @@ end
 @testset verbose = true showtiming = true ":double_integrator :min_tf :abstract :constr" begin
     sol2 = solve(ocp2, grid_size=100, print_level=0, tol=1e-12)
     @test sol2.objective ≈ 5.46 rtol=1e-2
+end
+
+Cd = 310
+β = 500
+Tmax = 3.5
+b = 2
+vmax = 0.1
+function F0(x)
+    r, v, m = x
+    D = Cd * v^2 * exp(-β*(r - 1))
+    return [ v, -D/m - 1/r^2, 0 ]
+end
+function F1(x)
+    r, v, m = x
+    return [ 0, Tmax/m, -b*Tmax ]
+end
+# NB. the ≤ is not the same as <= (parse error for <=)
+@def ocp3 begin
+    tf ∈ R, variable
+    t ∈ [ 0, tf ], time
+    x ∈ R^3, state
+    u ∈ R, control
+    0.1 ≤ tf ≤ Inf
+    r = x[1]
+    v = x[2]
+    m = x[3]
+    x(0) == [1, 0, 1]
+    m(tf) == 0.6
+    1 ≤ r(t) ≤ 1.1
+    0 ≤ v(t) ≤ vmax
+    0 ≤ u(t) ≤ 1
+    ẋ(t) == F0(x(t)) + u(t)*F1(x(t))
+    r(tf) → max
+end
+
+@testset verbose = true showtiming = true ":double_integrator :min_tf :abstract :constr" begin
+    sol3 = solve(ocp3, grid_size=100, print_level=0, tol=1e-12)    
+    @test sol3.objective ≈ 1.0125 rtol=1e-2
 end

test/suite/initial_guess.jl

@@ -58,7 +58,7 @@ end
 # constant initial guess
 x_const = [1.05, 0.2, 0.8]
 u_const = 0.5
-v_init = 0.15
+v_const = 0.15
 
 @testset verbose = true showtiming = true ":constant_init_x" begin
     sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const))
@@ -69,23 +69,23 @@ end
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_v" begin
-    sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_init))
+    sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_xu" begin
     sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_xv" begin
-    sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_init))
+    sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_uv" begin
-    sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_init))
+    sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_xuv" begin
-    sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init))
+    sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 
@@ -119,7 +119,7 @@ end
 # set initial guess in DOCP
 docp = directTranscription(ocp)
 @testset verbose = true showtiming = true ":DOCPInit_constant" begin
-    setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init))
+    setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const))
     sol = solve(docp, print_level=0)
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
@@ -137,7 +137,7 @@ end
 # pass initial guess to solve
 setDOCPInit(docp, OptimalControlInit()) # reset init in docp
 @testset verbose = true showtiming = true ":solve_constant_init" begin
-    sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init), print_level=0)
+    sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const), print_level=0)
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":solve_mixed_init" begin

test/test_initial_guess.jl

@@ -52,7 +52,7 @@ end
 dynamics!(ocp, (x, u, v) -> F0(x) + u*F1(x) )
 
 # reference solution
-sol0 = solve(ocp, print_level=0, tol=1e-8)
+sol0 = solve(ocp, print_level=0)
 #= check ok
 tf0 = sol0.variable
 println("tf ", tf0, " obj ", sol0.objective)
@@ -82,7 +82,7 @@ sol = solve(ocp, print_level=0, max_iter=maxiter)
 # constant initial guess
 x_const = [1.05, 0.2, 0.8]
 u_const = 0.5
-v_init = 0.15
+v_const = 0.15
 
 # Constant initial guess (vector for x; default for u,v)
 sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const), max_iter=maxiter)
@@ -93,23 +93,23 @@ sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const), max_ite
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for u; default for x,v):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for v; default for x,u)
-sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_init), max_iter=maxiter)
+sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for v; default for x,u):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for x,u; default for v)
 sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for x,u; default for v):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for x,v; default for u)
-sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_init), max_iter=maxiter)
+sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for x,v; default for u):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for u,v; default for x)
-sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_init), max_iter=maxiter)
+sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for u,v; default for x):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for x,u,v)
-sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init), max_iter=maxiter)
+sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for x,u,v):", sol.objective, sol.iterations)
 
 # functional initial guess
@@ -141,7 +141,7 @@ sol = solve(ocp, print_level=0, init=init_function_u = OptimalControlInit(sol0),
 println("\nSetting the initial guess at the DOCP level")
 docp = directTranscription(ocp)
 # constant vector init
-setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init))
+setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const))
 sol = solve(docp, print_level=0, max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess set in DOCP", sol.objective, sol.iterations)
 # mixed init
@@ -158,7 +158,7 @@ sol = solve(docp, print_level=0, max_iter=maxiter)
 println("\nPassing the initial guess to solve call")
 setDOCPInit(docp, OptimalControlInit()) # reset init in docp
 # constant vector init
-sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init), print_level=0, max_iter=maxiter)
+sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const), print_level=0, max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "constant initial guess passed to solve", sol.objective, sol.iterations)
 # mixed init
 sol = solve(docp, init=OptimalControlInit(x_init=x_func, u_init=u_const), print_level=0, max_iter=maxiter)