Merge branch 'master' into id-tutorial

Author: iliailmer

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "6.6.0"
+version = "7.0.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -16,13 +16,13 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
@@ -45,7 +45,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 AbstractTrees = "0.3"
-ArrayInterface = "2.8, 3.0"
+ArrayInterface = "3.1.39"
 ConstructionBase = "1"
 DataStructures = "0.17, 0.18"
 DiffEqBase = "6.54.0"
@@ -55,11 +55,11 @@ DiffRules = "0.1, 1.0"
 Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.7, 0.8"
 DomainSets = "0.5"
+Graphs = "1.4"
 IfElse = "0.1"
-JuliaFormatter = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17"
+JuliaFormatter = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18"
 LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
-LightGraphs = "1.3"
 MacroTools = "0.5"
 NaNMath = "0.3"
 NonlinearSolve = "0.3.8"
@@ -72,8 +72,8 @@ SciMLBase = "1.3"
 Setfield = "0.7, 0.8"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-SymbolicUtils = "0.16"
-Symbolics = "3.3.0"
+SymbolicUtils = "0.18"
+Symbolics = "4.0.0"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
 julia = "1.2"

docs/src/basics/AbstractSystem.md

@@ -56,6 +56,7 @@ Optionally, a system could have:
 
 - `observed(sys)`: All observed equations of the system and its subsystems.
 - `get_observed(sys)`: Observed equations of the current-level system.
+- `get_continuous_events(sys)`: `SymbolicContinuousCallback`s of the current-level system.
 - `get_defaults(sys)`: A `Dict` that maps variables into their default values.
 - `independent_variables(sys)`: The independent variables of a system.
 - `get_noiseeqs(sys)`: Noise equations of the current-level system.

docs/src/basics/Composition.md

@@ -32,10 +32,10 @@ end
 
 @parameters t
 D = Differential(t)
-@named connected = compose(ODESystem([
+connected = compose(ODESystem([
                         decay2.f ~ decay1.x
                         D(decay1.f) ~ 0
-                      ], t), decay1, decay2)
+                      ], t; name=:connected), decay1, decay2)
 
 equations(connected)
 
@@ -250,3 +250,94 @@ strongly connected components calculated during the process of simplification
 as the basis for building pre-simplified nonlinear systems in the implicit
 solving. In summary: these problems are structurally modified, but could be
 more efficient and more stable.
+
+## Components with discontinuous dynamics
+When modeling, e.g., impacts, saturations or Coulomb friction, the dynamic equations are discontinuous in either the state or one of its derivatives. This causes the solver to take very small steps around the discontinuity, and sometimes leads to early stopping due to `dt <= dt_min`. The correct way to handle such dynamics is to tell the solver about the discontinuity be means of a root-finding equation. [`ODEsystem`](@ref)s accept a keyword argument `continuous_events`
+```
+ODESystem(eqs, ...; continuous_events::Vector{Equation})
+ODESystem(eqs, ...; continuous_events::Pair{Vector{Equation}, Vector{Equation}})
+```
+where equations can be added that evaluate to 0 at discontinuities.
+
+To model events that have an affect on the state, provide `events::Pair{Vector{Equation}, Vector{Equation}}` where the first entry in the pair is a vector of equations describing event conditions, and the second vector of equations describe the affect on the state. The affect equations must be on the form
+```
+single_state_variable ~ expression_involving_any_variables
+```
+
+### Example: Friction
+The system below illustrates how this can be used to model Coulomb friction
+```julia
+using ModelingToolkit, OrdinaryDiffEq, Plots
+function UnitMassWithFriction(k; name)
+  @variables t x(t)=0 v(t)=0
+  D = Differential(t)
+  eqs = [
+    D(x) ~ v
+    D(v) ~ sin(t) - k*sign(v) # f = ma, sinusoidal force acting on the mass, and Coulomb friction opposing the movement
+  ]
+  ODESystem(eqs, t, continuous_events=[v ~ 0], name=name) # when v = 0 there is a discontinuity
+end
+@named m = UnitMassWithFriction(0.7)
+prob = ODEProblem(m, Pair[], (0, 10pi))
+sol = solve(prob, Tsit5())
+plot(sol)
+```
+
+### Example: Bouncing ball
+In the documentation for DifferentialEquations, we have an example where a bouncing ball is simulated using callbacks which has an `affect!` on the state. We can model the same system using ModelingToolkit like this
+
+```julia
+@variables t x(t)=1 v(t)=0
+D = Differential(t)
+
+root_eqs = [x ~ 0]  # the event happens at the ground x(t) = 0
+affect   = [v ~ -v] # the effect is that the velocity changes sign
+
+@named ball = ODESystem([
+    D(x) ~ v
+    D(v) ~ -9.8
+], t, continuous_events = root_eqs => affect) # equation => affect
+
+ball = structural_simplify(ball)
+
+tspan = (0.0,5.0)
+prob = ODEProblem(ball, Pair[], tspan)
+sol = solve(prob,Tsit5())
+@assert 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
+plot(sol)
+```
+
+### Test bouncing ball in 2D with walls
+Multiple events? No problem! This example models a bouncing ball in 2D that is enclosed by two walls at $y = \pm 1.5$.
+```julia
+@variables t x(t)=1 y(t)=0 vx(t)=0 vy(t)=2
+D = Differential(t)
+
+continuous_events = [ # This time we have a vector of pairs
+    [x ~ 0] => [vx ~ -vx]
+    [y ~ -1.5, y ~ 1.5] => [vy ~ -vy]
+]
+
+@named ball = ODESystem([
+    D(x)  ~ vx,
+    D(y)  ~ vy,
+    D(vx) ~ -9.8-0.1vx, # gravity + some small air resistance
+    D(vy) ~ -0.1vy,
+], t, continuous_events = continuous_events)
+
+
+ball = structural_simplify(ball)
+
+tspan = (0.0,10.0)
+prob = ODEProblem(ball, Pair[], tspan)
+
+sol = solve(prob,Tsit5())
+@assert 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
+@assert minimum(sol[y]) > -1.5 # check wall conditions
+@assert maximum(sol[y]) < 1.5  # check wall conditions
+
+tv = sort([LinRange(0, 10, 200); sol.t])
+plot(sol(tv)[y], sol(tv)[x], line_z=tv)
+vline!([-1.5, 1.5], l=(:black, 5), primary=false)
+hline!([0], l=(:black, 5), primary=false)
+```

docs/src/internals.md

@@ -16,3 +16,18 @@ plotting their output, these relationships are stored and are then used to
 generate the `observed` equation found in the `SciMLFunction` interface, so that
 `sol[x]` lazily reconstructs the observed variable when necessary. In this sense,
 there is an equivalence between observables and the variable elimination system.
+
+The procedure for variable elimination inside [`structural_simplify`](@ref) is
+1. [`ModelingToolkit.initialize_system_structure`](@ref).
+2. [`ModelingToolkit.alias_elimination`](@ref). This step moves equations into `observed(sys)`.
+3. [`ModelingToolkit.dae_index_lowering`](@ref) by means of [`pantelides!`](@ref) (if the system is an [`ODESystem`](@ref)).
+4. [`ModelingToolkit.tearing`](@ref).
+
+## Preparing a system for simulation
+Before a simulation or optimization can be performed, the symbolic equations stored in an [`AbstractSystem`](@ref) must be converted into executable code. This step is typically occurs after the simplification explained above, and is performed when an instance of a [`AbsSciMLBase.SciMLProblem`](@ref), such as a [`ODEProblem`](@ref), is constructed.
+The call chain typically looks like this, with the function names in the case of an `ODESystem` indicated in parenthesis
+1. Problem constructor ([`ODEProblem`](@ref))
+2. Build an `DEFunction` ([`process_DEProblem`](@ref) -> [`ODEFunction`](@ref)
+3. Write actual executable code ([`generate_function`](@ref))
+
+Apart from [`generate_function`](@ref), which generates the dynamics function, `ODEFunction` also builds functions for observed equations (`build_explicit_observed_function`) and jacobians (`generate_jacobian`) etc. These are all stored in the `ODEFunction`.

docs/src/mtkitize_tutorials/modelingtoolkitize_index_reduction.md

@@ -29,7 +29,7 @@ p = [9.8, 1]
 tspan = (0, 10.0)
 pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
 traced_sys = modelingtoolkitize(pendulum_prob)
-pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
+pendulum_sys = structural_simplify(traced_sys)
 prob = ODAEProblem(pendulum_sys, Pair[], tspan)
 sol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)
 plot(sol, vars=states(traced_sys))

docs/src/tutorials/ode_modeling.md

@@ -211,7 +211,7 @@ again are just algebraic relations:
 connections = [ fol_1.f ~ 1.5,
                 fol_2.f ~ fol_1.x ]
 
-@named connected = compose(ODESystem(connections), fol_1, fol_2)
+connected = compose(ODESystem(connections,name=:connected), fol_1, fol_2)
       # Model connected with 5 equations
       # States (5):
       #   fol_1₊f(t)

src/ModelingToolkit.jl

@@ -50,13 +50,13 @@ import Symbolics: rename, get_variables!, _solve, hessian_sparsity,
                   tosymbol, lower_varname, diff2term, var_from_nested_derivative,
                   BuildTargets, JuliaTarget, StanTarget, CTarget, MATLABTarget,
                   ParallelForm, SerialForm, MultithreadedForm, build_function,
-                  unflatten_long_ops, rhss, lhss, prettify_expr, gradient,
+                  rhss, lhss, prettify_expr, gradient,
                   jacobian, hessian, derivative, sparsejacobian, sparsehessian,
                   substituter, scalarize, getparent
 
 import DiffEqBase: @add_kwonly
 
-import LightGraphs: SimpleDiGraph, add_edge!, incidence_matrix
+import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
 
 using Requires
 
@@ -134,6 +134,7 @@ include("systems/control/controlsystem.jl")
 
 include("systems/pde/pdesystem.jl")
 
+include("systems/sparsematrixclil.jl")
 include("systems/discrete_system/discrete_system.jl")
 include("systems/validation.jl")
 include("systems/dependency_graphs.jl")