Skip to content

Commit

Permalink
Merge pull request #91 from SciML/tutorials
Browse files Browse the repository at this point in the history
flesh out some of the getting started tutorials
  • Loading branch information
ChrisRackauckas authored Dec 4, 2022
2 parents 0ed0e8b + 355dcde commit d93d142
Show file tree
Hide file tree
Showing 3 changed files with 202 additions and 29 deletions.
108 changes: 97 additions & 11 deletions docs/src/getting_started/find_root.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,17 +3,25 @@
A nonlinear system $$f(u) = 0$$ is specified by defining a function `f(u,p)`,
where `p` are the parameters of the system. Many problems can be written in
such as way that solving a nonlinear rootfinding problem gives the solution.
For example:
For example:

* Do you want to know ``u`` such that ``4^u + 6^u = 7^u``? Then solve
* Do you want to know ``u`` such that ``4^u + 6^u = 7^u``? Then solve
``f(u) = 4^u + 6^u - 7^u = 0`` for `u`!
* If you have an ODE ``u' = f(u)``, what is the point where the solution
will be completely still, i.e. `u' = 0`?
*

All of these problems are solved by using a numerical rootfinder. Let's solve
All of these problems are solved by using a numerical rootfinder. Let's solve
our first rootfind problem!

## Required Dependencies

The following parts of the SciML Ecosystem will be used in this tutorial:

| Module | Description |
| ----------- | ----------- |
| [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) | The symbolic modeling environment |
| [NonlinearSolve.jl](https://docs.sciml.ai/NonlinearSolve/stable/) | The symbolic modeling environment |

## Problem Setup

For example, the following solves the vector equation:
Expand Down Expand Up @@ -51,43 +59,121 @@ sol = solve(prob,NewtonRaphson())
@show sol.u, prob.f(sol.u,prob.p)
```

## Step 1: Import the Packages
## Step by Step Solution

### Step 1: Import the Packages

To do this tutorial we will need a few components:

* [ModelingToolkit.jl, our modeling environment](https://docs.sciml.ai/ModelingToolkit/stable/)
* [NonlinearSolve.jl, the nonlinear system solvers](https://docs.sciml.ai/NonlinearSolve/stable/)

To start, let's add these packages [as demonstrated in the installation tutorial](@ref installation):

```julia
]add ModelingToolkit NonlinearSolve
```

Now we're ready. Let's load in these packages:

```@example first_rootfind
# Import the packages
using ModelingToolkit, NonlinearSolve
```

## Step 2: Define the Nonlinear System
### Step 2: Define the Nonlinear System

Now let's define our nonlinear system. We use the `ModelingToolkit.@variabes` statement to
declare our 3 state variables:

```@example first_rootfind
# Define the nonlinear system
@variables x=1.0 y=0.0 z=0.0
```

Notice that we are using the form `state = initial condition`. This is a nice shorthand
for coupling an initial condition to our states. We now must similarly define our parameters,
which we can associate default values via the form `parameter = default value`. This looks
like:

```@example first_rootfind
@parameters σ=10.0 ρ=26.0 β=8/3
```

Now we create an array of equations to define our nonlinear system that must be satisfied.
This looks as follows:

!!! note

Note that in ModelingToolkit and Symbolics, `~` is used for equation equality. This is
separate from `=` which is the "assignment operator" in the Julia programming language.
For example, `x = x + 1` is a valid assignment in a programming language, and it is
invalid for that to represent "equality", which is the reason why a separate operator
is used!

```@example first_rootfind
eqs = [0 ~ σ*(y-x),
0 ~ x*(ρ-z)-y,
0 ~ x*y - β*z]
```

Finally we bring these pieces together, the equation along with its states and parameters,
define our `NonlinearSystem`:

```@example first_rootfind
@named ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
```

## Step 3: Convert the Symbolic Problem to a Numerical Problem
### Step 3: Convert the Symbolic Problem to a Numerical Problem

Now that we have simplified our system, let's turn it into a numerical problem to
approximate. This is done with the `NonlinearProblem` constructor, that transforms it from
a symbolic `ModelingToolkit` representation to a numerical `NonlinearSolve`
representation. We need to tell it the numerical details for whether to override any of the
default values for the initial conditions and parameters.

In this case, we will use the default values for all our variables, so we will pass a
blank override `[]`. This looks like:

```@example first_rootfind
# Convert the symbolic system into a numerical system
prob = NonlinearProblem(ns,[])
```

## Step 4: Solve the Numerical Problem
If for example we did want to change the initial condition of `x`
to `2.0` and the parameter `σ` to `4.0`, we would do `[x => 2.0, σ => 4.0]`. This looks
like:

```@example first_rootfind
prob2 = NonlinearProblem(ns,[x => 2.0, σ => 4.0])
```

### Step 4: Solve the Numerical Problem

Now we solve the nonlinear system. For this we choose a solver from the
[NonlinearSolve.jl's solver options](https://docs.sciml.ai/NonlinearSolve/dev/solvers/NonlinearSystemSolvers/)
We will choose `NewtonRaphson` as follows:

```@example first_rootfind
# Solve the numerical problem
sol = solve(prob,NewtonRaphson())
```

## Step 5: Analyze the Solution
### Step 5: Analyze the Solution

Now let's check out the solution. First of all, what kind of thing is the `sol`? We can
see that by asking for its type:

```@example first_rootfind
typeof(sol)
```

From this we can see that it is an `NonlinearSolution`. We can see the documentation for
how to use the `NonlinearSolution` by checking the
[NonlinearSolve.jl solution type page](https://docs.sciml.ai/NonlinearSolve/dev/basics/NonlinearSolution/) Thus the solution is stored as `.u`. What is the solution to our
nonlinear system and what is the final residual value? We can check it as follows:

```@example first_rootfind
# Analyze the solution
@show sol.u, prob.f(sol.u,prob.p)
```
@show sol.u, sol.resid
```
106 changes: 91 additions & 15 deletions docs/src/getting_started/first_optimization.md
Original file line number Diff line number Diff line change
@@ -1,17 +1,25 @@
# [Solve your first optimization problem](@id first_opt)

**Numerical optimization** is the process of finding some numerical values that
minimize some equation.
minimize some equation.

* How much fuel should you put into an airplane to have the minimum weight that
* How much fuel should you put into an airplane to have the minimum weight that
can go to its destination?
* What parameters should I choose for my simulation so that it minimizes the
distance of its predictions from my experimental data?
*

All of these are examples of problems solved by numerical optimization.
All of these are examples of problems solved by numerical optimization.
Let's solve our first optimization problem!

## Required Dependencies

The following parts of the SciML Ecosystem will be used in this tutorial:

| Module | Description |
| ----------- | ----------- |
| [Optimization.jl](https://docs.sciml.ai/Optimization/stable/) | The numerical optimization package |
| [OptimizationNLopt.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/nlopt/) | The NLopt optimizers we will use |

## Problem Setup

First, what are we solving? Let's take a look at the Rosenbrock equation:
Expand All @@ -20,10 +28,10 @@ First, what are we solving? Let's take a look at the Rosenbrock equation:
L(u,p) = (p_1 - u_1)^2 + p_2 * (u_2 - u_1)^2
```

What we want to do is find the values of ``u_1`` and ``u_2`` such that ``L``
achieves its minimum value possible. We will do this under a few constraints:
we want to find this optima within some bounded domain, i.e. ``u_i \in [-1,1]``.
This should be done with the parameter values ``p_1 = 1.0`` and `p_2 = 100.0``.
What we want to do is find the values of ``u_1`` and ``u_2`` such that ``L``
achieves its minimum value possible. We will do this under a few constraints:
we want to find this optima within some bounded domain, i.e. ``u_i \in [-1,1]``.
This should be done with the parameter values ``p_1 = 1.0`` and `p_2 = 100.0``.
What should ``u = [u_1,u_2]`` be to achieve this goal? Let's dive in!

!!! note
Expand All @@ -34,7 +42,7 @@ What should ``u = [u_1,u_2]`` be to achieve this goal? Let's dive in!
## Copy-Pastable Code

```@example
# Import the package
# Import the package
using Optimization, OptimizationNLopt
# Define the problem to optimize
Expand All @@ -50,35 +58,103 @@ sol = solve(prob,NLopt.LD_LBFGS())
@show sol.u, L(sol.u,p)
```

## Step 1: Import the packages
## Step by Step Solution

### Step 1: Import the packages

To do this tutorial we will need a few components:

* [Optimization.jl](https://docs.sciml.ai/Optimization/stable/), the optimization interface.
* [OptimizationNLopt.jl](https://docs.sciml.ai/Optimization/stable/optimization_packages/nlopt/), the optimizers we will use.

Note that Optimization.jl is an interface for optimizers, and thus we always have to choose
which optimizer we want to use. Here we choose to demonstrate `OptimizationNLopt` because
of its efficiency and versitility. But there are many other possible choices. Check out
the
[solver compatability chart](https://docs.sciml.ai/Optimization/stable/#Overview-of-the-Optimizers)
for a quick overview of what optimizer packages offer.

To start, let's add these packages [as demonstrated in the installation tutorial](@ref installation):

```julia
]add Optimization OptimizationNLopt
```

Now we're ready. Let's load in these packages:

```@example first_opt
```@example first_opt
using Optimization, OptimizationNLopt
```

## Step 2: Define the Optimization Problem
### Step 2: Define the Optimization Problem

Now let's define our problem to optimize. We start by defining our loss function. In
Optimization.jl's `OptimizationProblem` interface, the states are given by an array
`u`. Thus we can designate `u[1]` to be `u_1` and `u[2]` to be `u_2`, similarly with our
parameters, and write out the loss function on a vector-defined state as follows:

```@example first_opt
# Define the problem to optimize
L(u,p) = (p[1] - u[1])^2 + p[2] * (u[2] - u[1]^2)^2
```

Now we need to define our `OptimizationProblem`. If you need help remembering how to define
the `OptimizationProblem`, you can always refer to the
[Optimization.jl problem definition page](https://docs.sciml.ai/Optimization/stable/API/optimization_problem/)

Thus what we need to define is an initial condition `u0` and our parameter vector `p`.
We will make our initial condition have both values as zero, which is done by the Julia
shorthand `zeros(2)` that creates a vector `[0.0,0.0]`. We manually define the parameter
vector `p` to input our values. Then we set the lower bound and upper bound for the
optimization as follows:

```@example first_opt
u0 = zeros(2)
p = [1.0,100.0]
prob = OptimizationProblem(L, u0, p, lb = [-1.0,-1.0], ub = [1.0,1.0])
```

## Step 3: Solve the Optimization Problem
#### Note about defining uniform bounds

Note that we can simplify the code a bit for the lower and upper bound definition by
using the Julia Base command `ones`, which returns a vector where each value is a one.
Thus for example, `ones(2)` is equivalent to `[1.0,1.0]`. Therefore `-1 * ones(2)` is
equivalent to `[-1.0, 1.0]`, meaning we could have written our problem as follows:

```@example first_opt
prob = OptimizationProblem(L, u0, p, lb = -1 * ones(2), ub = ones(2))
```

### Step 3: Solve the Optimization Problem

Now we solve the `OptimizationProblem` that we have defined. This is done by passing
our `OptimizationProblem` `prob` along with a chosen solver to the `solve` command. At
the beginning we explained that we will use the `OptimizationNLopt` set of solvers, which
are
[documented in the OptimizationNLopt page](https://docs.sciml.ai/Optimization/stable/optimization_packages/nlopt/).
From here we are choosing the `NLopt.LD_LBFGS()` for its mixture of robustness and
performance. To perform this solve, we do the following:

```@example first_opt
# Solve the optimization problem
sol = solve(prob,NLopt.LD_LBFGS())
```

## Step 4: Analyze the Solution
### Step 4: Analyze the Solution

Now let's check out the solution. First of all, what kind of thing is the `sol`? We can
see that by asking for its type:

```@example first_opt
typeof(sol)
```

From this we can see that it is an `OptimizationSolution`. We can see the documentation for
how to use the `OptimizationSolution` by checking the
[Optimization.jl solution type page](https://docs.sciml.ai/Optimization/dev/API/optimization_solution/). Thus the solution is stored as `.u`. What is the solution to our
optimization and what is the final loss value? We can check it as follows:

```@example first_opt
# Analyze the solution
@show sol.u, L(sol.u,p)
```
```
17 changes: 14 additions & 3 deletions docs/src/getting_started/first_simulation.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,23 @@ In this tutorial we will build and run our first simulation with SciML!
This tutorial assumes that you have already installed Julia on your system. If you have
not done so already, please [follow the installation tutorial first](@ref installation).

To build our simulation, we will use the [ModelingToolkit](https://mtk.sciml.ai/dev/)
To build our simulation, we will use the
[ModelingToolkit](https://docs.sciml.ai/ModelingToolkit/stable/)
system for modeling and simulation. ModelingToolkit is a bit higher level than directly
defining code for a differential equation system: it's a symbolic system that will
automatically simplify our models, optimize our code, and generate compelling visualizations.
Sounds neat? Let's dig in.

## Required Dependencies

The following parts of the SciML Ecosystem will be used in this tutorial:

| Module | Description |
| ----------- | ----------- |
| [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) | The symbolic modeling environment |
| [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) | The differential equation solvers |
| [Plots.jl](https://docs.juliaplots.org/stable/) | The plotting and visualization package |

## Our Problem: Simulate the Lotka-Volterra Predator-Prey Dynamics

The dynamics of our system are given by the
Expand Down Expand Up @@ -82,8 +93,8 @@ plot(p1,p2,layout=(2,1))

To do this tutorial we will need a few components:

* [ModelingToolkit.jl, our modeling environment](https://mtk.sciml.ai/dev/)
* [DifferentialEquations.jl, the differential equation solvers](https://diffeq.sciml.ai/stable/)
* [ModelingToolkit.jl, our modeling environment](https://docs.sciml.ai/ModelingToolkit/stable/)
* [DifferentialEquations.jl, the differential equation solvers](https://docs.sciml.ai/DiffEqDocs/stable/)
* [Plots.jl, our visualization tool](https://docs.juliaplots.org/stable/)

To start, let's add these packages [as demonstrated in the installation tutorial](@ref installation):
Expand Down

0 comments on commit d93d142

Please sign in to comment.