Merge pull request #836 from SciML/tutorials

add a bunch of tutorials

Author: ChrisRackauckas

docs/make.jl

@@ -16,6 +16,9 @@ makedocs(
             "tutorials/higher_order.md",
             "tutorials/nonlinear.md",
             "tutorials/modelingtoolkitize.md",
+            "tutorials/optimization.md",
+            "tutorials/stochastic_diffeq.md",
+            "tutorials/nonlinear_optimal_control.md"
         ],
         "Basics" => Any[
             "basics/AbstractSystem.md",

docs/src/index.md

@@ -58,7 +58,7 @@ is built on, consult the
 - Optimization problems
 - Continuous-Time Markov Chains
 - Chemical Reactions
-- Optimal Control
+- Nonlinear Optimal Control
 
 ## Model Import Formats
 

docs/src/tutorials/acausal_components.md

@@ -7,7 +7,7 @@ are then applied to state equalities between currents and voltages. This
 specifies a differential-algebraic equation (DAE) system, where the algebraic
 equations are given by the constraints and equalities between different
 component variables. We then simplify this to an ODE by eliminating the
-equalities before solving. Let's see this in action
+equalities before solving. Let's see this in action.
 
 ## Copy-Paste Example
 

docs/src/tutorials/nonlinear.md

@@ -1,4 +1,4 @@
-# Solving Nonlinear Systems with NLsolve
+# Modeling Nonlinear Systems
 
 In this example we will go one step deeper and showcase the direct function
 generation capabilities in ModelingToolkit.jl to build nonlinear systems.
@@ -7,7 +7,7 @@ the nonlinear system defined by where the derivatives are zero. We use (unknown)
 variables for our nonlinear system.
 
 ```julia
-using ModelingToolkit
+using ModelingToolkit, NonlinearSolve
 
 @variables x y z
 @parameters σ ρ β
@@ -17,91 +17,25 @@ eqs = [0 ~ σ*(y-x),
        0 ~ x*(ρ-z)-y,
        0 ~ x*y - β*z]
 ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
-nlsys_func = generate_function(ns)[2] # second is the inplace version
-```
-
-which generates:
-
-```julia
-(var"##MTIIPVar#405", u, p)->begin
-        @inbounds begin
-                @inbounds begin
-                        let (x, y, z, σ, ρ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
-                            var"##MTIIPVar#405"[1] = (*)(σ, (-)(y, x))
-                            var"##MTIIPVar#405"[2] = (-)((*)(x, (-)(ρ, z)), y)
-                            var"##MTIIPVar#405"[3] = (-)((*)(x, y), (*)(β, z))
-                        end
-                    end
-            end
-        nothing
-    end
-```
 
-We can use this to build a nonlinear function for use with NLsolve.jl:
+guess = [x => 1.0,
+         y => 0.0,
+         z => 0.0]
 
-```julia
-f = eval(nlsys_func)
-du = zeros(3); u = ones(3)
-params = (10.0,26.0,2.33)
-f(du,u,params)
-du
-
-#=
-3-element Array{Float64,1}:
-  0.0
- 24.0
- -1.33
- =#
-```
+ps = [
+      σ => 10.0
+      ρ => 26.0
+      β => 8/3
+      ]
 
-We can similarly ask to generate the in-place Jacobian function:
-
-```julia
-j_func = generate_jacobian(ns)[2] # second is in-place
-j! = eval(j_func)
+prob = NonlinearProblem(ns,guess,ps)
+sol = solve(prob,NewtonRaphson())
 ```
 
-which gives:
+We can similarly ask to generate the `NonlinearProblem` with the analytical
+Jacobian function:
 
 ```julia
-:((var"##MTIIPVar#582", u, p)->begin
-          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:70 =#
-          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:71 =#
-          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:71 =# @inbounds begin
-                  #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:72 =#
-                  #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:53 =# @inbounds begin
-                          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:53 =#
-                          let (x, y, z, σ, ρ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
-                              var"##MTIIPVar#582"[1] = (*)(σ, -1)
-                              var"##MTIIPVar#582"[2] = (-)(ρ, z)
-                              var"##MTIIPVar#582"[3] = y
-                              var"##MTIIPVar#582"[4] = σ
-                              var"##MTIIPVar#582"[5] = -1
-                              var"##MTIIPVar#582"[6] = x
-                              var"##MTIIPVar#582"[7] = 0
-                              var"##MTIIPVar#582"[8] = (*)(x, -1)
-                              var"##MTIIPVar#582"[9] = (*)(-1, β)
-                          end
-                      end
-              end
-          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:74 =#
-          nothing
-      end)
+prob = NonlinearProblem(ns,guess,ps,jac=true)
+sol = solve(prob,NewtonRaphson())
 ```
-
-Now, we can call `nlsolve` by enclosing our parameters into the functions:
-
-```julia
-using NLsolve
-nlsolve((out, x) -> f(out, x, params), (out, x) -> j!(out, x, params), ones(3))
-```
-
-If one would like the generated function to be a Julia function instead of an expression, and allow this
-function to be used from within the same world-age, one simply needs to pass `Val{false}` to tell it to
-generate the function, i.e.:
-
-```julia
-nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β], expression=Val{false})[2]
-```
-
-which uses RuntimeGeneratedFunctions.jl to build the same world-age function on the fly without eval.

docs/src/tutorials/nonlinear_optimal_control.md

@@ -0,0 +1,93 @@
+# Nonlinear Optimal Control
+
+#### Note: this is still a work in progress!
+
+The `ControlSystem` type is an interesting system because, unlike other
+system types, it cannot be numerically solved on its own. Instead, it must be
+transformed into another system before solving. Standard methods such as the
+"direct method", "multiple shooting", or "discretize-then-optimize" can all be
+phrased as symbolic transformations to a `ControlSystem`: this is the strategy
+of this methodology.
+
+## Defining a Nonlinear Optimal Control Problem
+
+Here we will start by defining a classic optimal control problem. Let:
+
+```math
+x^{′′} = u^3(t)
+```
+
+where we want to optimize our controller `u(t)` such that the following is
+minimized:
+
+```math
+L(\theta) = \sum_i \Vert 4 - x(t_i) \Vert + 2 \Vert x^\prime(t_i) \Vert + \Vert u(t_i) \Vert
+```
+
+where ``i`` is measured on (0,8) at 0.01 intervals. To do this, we rewrite the
+ODE in first order form:
+
+```math
+\begin{aligned}
+x^\prime &= v \\
+v^′ &= u^3(t) \\
+\end{aligned}
+```
+
+and thus
+
+```math
+L(\theta) = \sum_i \Vert 4 - x(t_i) \Vert + 2 \Vert v(t_i) \Vert + \Vert u(t_i) \Vert
+```
+
+is our loss function on the first order system.
+
+Defining such a control system is similar to an `ODESystem`, except we must also
+specify a control variable `u(t)` and a loss function. Together, this problem
+looks as follows:
+
+```julia
+using ModelingToolkit
+
+@variables t x(t) v(t) u(t)
+@parameters p[1:2]
+D = Differential(t)
+
+loss = (4-x)^2 + 2v^2 + u^2
+eqs = [
+    D(x) ~ v - p[2]*x
+    D(v) ~ p[1]*u^3 + v
+]
+
+sys = ControlSystem(loss,eqs,t,[x,v],[u],p)
+```
+
+## Solving a Control Problem via Discretize-Then-Optimize
+
+One common way to solve nonlinear optimal control problems is by transforming
+them into an optimization problem by performing a Runge-Kutta discretization
+of the differential equation system and imposing equalities between variables
+in the same steps. This can be done via the `runge_kutta_discretize` transformation
+on the `ControlSystem`. While a tableau `tab` can be specified, it defaults to
+a 5th order RadauIIA collocation, which is a common method in the field. To
+perform this discretization, we simply need to give a `dt` and a timespan on which
+to discretize:
+
+```julia
+dt = 0.1
+tspan = (0.0,1.0)
+sys = runge_kutta_discretize(sys,dt,tspan)
+```
+
+Now `sys` is an `OptimizationSystem` which, when solved, gives the values of
+`x(t)`, `v(t)`, and `u(t)`. Thus we solve the `OptimizationSystem` using
+GalacticOptim.jl:
+
+```julia
+u0 = rand(length(states(sys))) # guess for the state values
+prob = OptimizationProblem(sys,u0,[0.1,0.1],grad=true)
+sol = solve(prob,BFGS())
+```
+
+And this is missing some nice interfaces and ignores the equality constraints
+right now so the tutorial is not complete.

docs/src/tutorials/optimization.md

@@ -0,0 +1,26 @@
+# Modeling Optimization Problems
+
+```julia
+using ModelingToolkit, GalacticOptim
+
+@variables x y
+@parameters a b
+loss = (a - x)^2 + b * (y - x^2)^2
+sys = OptimizationSystem(loss,[x,y],[a,b])
+
+u0 = [
+    x=>1.0
+    y=>2.0
+]
+p = [
+    a => 6.0
+    b => 7.0
+]
+
+prob = OptimizationProblem(sys,u0,p,grad=true,hess=true)
+solve(prob,Newton())
+```
+
+Needs more text but it's super cool and auto-parallelizes and sparsifies too.
+Plus you can hierarchically nest systems to have it generate huge
+optimization problems.

docs/src/tutorials/stochastic_diffeq.md

@@ -0,0 +1,42 @@
+# Modeling with Stochasticity
+
+All models with `ODESystem` are deterministic. `SDESystem` adds another element
+to the model: randomness. This is a
+[stochastic differential equation](https://en.wikipedia.org/wiki/Stochastic_differential_equation)
+which has a deterministic (drift) component and a stochastic (diffusion)
+component. Let's take the Lorenz equation from the first tutorial and extend
+it to have multiplicative noise.
+
+```julia
+using ModelingToolkit, StochasticDiffEq
+
+# Define some variables
+@parameters t σ ρ β
+@variables x(t) y(t) z(t)
+D = Differential(t)
+
+eqs = [D(x) ~ σ*(y-x),
+       D(y) ~ x*(ρ-z)-y,
+       D(z) ~ x*y - β*z]
+
+noiseeqs = [0.1*x,
+            0.1*y,
+            0.1*z]
+
+de = SDESystem(eqs,noiseeqs,t,[x,y,z],[σ,ρ,β])
+
+u0map = [
+    x => 1.0,
+    y => 0.0,
+    z => 0.0
+]
+
+parammap = [
+    σ => 10.0,
+    β => 26.0,
+    ρ => 2.33
+]
+
+prob = SDEProblem(de,u0map,(0.0,100.0),parammap)
+sol = solve(prob,SOSRI())
+```

src/systems/optimization/optimizationsystem.jl

@@ -242,3 +242,54 @@ function DiffEqBase.OptimizationFunction{iip}(f, ::AutoModelingToolkit, x, p = D
     end
     OptimizationProblem(sys,u0map,parammap,grad=grad,hess=hess).f
 end
+
+function Base.show(io::IO, sys::OptimizationSystem)
+    eqs = equations(sys)
+    Base.printstyled(io, "Model $(nameof(sys))\n"; bold=true)
+    # The reduced equations are usually very long. It's not that useful to print
+    # them.
+    #Base.print_matrix(io, eqs)
+    #println(io)
+
+    rows = first(displaysize(io)) ÷ 5
+    limit = get(io, :limit, false)
+
+    vars = states(sys); nvars = length(vars)
+    Base.printstyled(io, "States ($nvars):"; bold=true)
+    nrows = min(nvars, limit ? rows : nvars)
+    limited = nrows < length(vars)
+    d_u0 = has_default_u0(sys) ? default_u0(sys) : nothing
+    for i in 1:nrows
+        s = vars[i]
+        print(io, "\n  ", s)
+
+        if d_u0 !== nothing
+            val = get(d_u0, s, nothing)
+            if val !== nothing
+                print(io, " [defaults to $val]")
+            end
+        end
+    end
+    limited && print(io, "\n⋮")
+    println(io)
+
+    vars = parameters(sys); nvars = length(vars)
+    Base.printstyled(io, "Parameters ($nvars):"; bold=true)
+    nrows = min(nvars, limit ? rows : nvars)
+    limited = nrows < length(vars)
+    d_p = has_default_p(sys) ? default_p(sys) : nothing
+    for i in 1:nrows
+        s = vars[i]
+        print(io, "\n  ", s)
+
+        if d_p !== nothing
+            val = get(d_p, s, nothing)
+            if val !== nothing
+                print(io, " [defaults to $val]")
+            end
+        end
+    end
+    limited && print(io, "\n⋮")
+
+    return nothing
+end

test/controlsystem.jl

@@ -15,6 +15,6 @@ dt = 0.1
 tspan = (0.0,1.0)
 sys = runge_kutta_discretize(sys,dt,tspan)
 
-u0 = rand(112) # guess for the state values
-prob = OptimizationProblem(sys,u0,[0.1],grad=true)
-
+u0 = rand(length(states(sys))) # guess for the state values
+prob = OptimizationProblem(sys,u0,[0.1,0.1],grad=true)
+sol = solve(prob,BFGS())