diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index 580b7511eb..9c79359112 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -1 +1,2 @@
style = "sciml"
+format_markdown = true
\ No newline at end of file
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index dc3285bdd3..57a2171b50 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -1,3 +1,3 @@
-- This repository follows the [SciMLStyle](https://github.com/SciML/SciMLStyle) and the SciML [ColPrac](https://github.com/SciML/ColPrac).
-- Please run `using JuliaFormatter, ModelingToolkit; format(joinpath(dirname(pathof(ModelingToolkit)), ".."))` before commiting.
-- Add tests for any new features.
+ - This repository follows the [SciMLStyle](https://github.com/SciML/SciMLStyle) and the SciML [ColPrac](https://github.com/SciML/ColPrac).
+ - Please run `using JuliaFormatter, ModelingToolkit; format(joinpath(dirname(pathof(ModelingToolkit)), ".."))` before commiting.
+ - Add tests for any new features.
diff --git a/LICENSE.md b/LICENSE.md
index 58a912aa0d..947cd9843e 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -2,43 +2,36 @@ The ModelingToolkit.jl package is licensed under the MIT "Expat" License:
> Copyright (c) 2018-22: Yingbo Ma, Christopher Rackauckas, Julia Computing, and
> contributors
->
->
+>
> Permission is hereby granted, free of charge, to any person obtaining a copy
->
+>
> of this software and associated documentation files (the "Software"), to deal
->
+>
> in the Software without restriction, including without limitation the rights
->
+>
> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
->
+>
> copies of the Software, and to permit persons to whom the Software is
->
+>
> furnished to do so, subject to the following conditions:
->
->
->
+>
> The above copyright notice and this permission notice shall be included in all
->
+>
> copies or substantial portions of the Software.
->
->
->
+>
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
->
+>
> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
->
+>
> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
->
+>
> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
->
+>
> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
->
+>
> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
->
+>
> SOFTWARE.
->
->
The code in `src/structural_transformation/bipartite_tearing/modia_tearing.jl`,
which is from the [Modia.jl](https://github.com/ModiaSim/Modia.jl) project, is
diff --git a/NEWS.md b/NEWS.md
index 222864f76c..8adccc071e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,10 +1,9 @@
# ModelingToolkit v8 Release Notes
-
### Upgrade guide
-- `connect` should not be overloaded by users anymore. `[connect = Flow]`
- informs ModelingToolkit that particular variable in a connector ought to sum
- to zero, and by default, variables are equal in a connection. Please check out
- [acausal components tutorial](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/acausal_components/)
- for examples.
+ - `connect` should not be overloaded by users anymore. `[connect = Flow]`
+ informs ModelingToolkit that particular variable in a connector ought to sum
+ to zero, and by default, variables are equal in a connection. Please check out
+ [acausal components tutorial](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/acausal_components/)
+ for examples.
diff --git a/README.md b/README.md
index 00432ab634..e75760ac56 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,5 @@
# ModelingToolkit.jl
-
[](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged)
[](https://docs.sciml.ai/ModelingToolkit/stable/)
@@ -8,7 +7,7 @@
[](https://coveralls.io/github/SciML/ModelingToolkit.jl?branch=master)
[](https://github.com/SciML/ModelingToolkit.jl/actions?query=workflow%3ACI)
-[](https://github.com/SciML/ColPrac)
+[](https://github.com/SciML/ColPrac)
[](https://github.com/SciML/SciMLStyle)
ModelingToolkit.jl is a modeling framework for high-performance symbolic-numeric computation
@@ -121,7 +120,9 @@ plot(sol, idxs = (a, lorenz1.x, lorenz2.z))
# Citation
+
If you use ModelingToolkit.jl in your research, please cite [this paper](https://arxiv.org/abs/2103.05244):
+
```
@misc{ma2021modelingtoolkit,
title={ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling},
diff --git a/docs/src/basics/AbstractSystem.md b/docs/src/basics/AbstractSystem.md
index ad64ef30c4..2f5627eff0 100644
--- a/docs/src/basics/AbstractSystem.md
+++ b/docs/src/basics/AbstractSystem.md
@@ -10,28 +10,29 @@ model manipulation and compilation.
### Subtypes
There are three immediate subtypes of `AbstractSystem`, classified by how many independent variables each type has:
-* `AbstractTimeIndependentSystem`: has no independent variable (e.g.: `NonlinearSystem`)
-* `AbstractTimeDependentSystem`: has a single independent variable (e.g.: `ODESystem`)
-* `AbstractMultivariateSystem`: may have multiple independent variables (e.g.: `PDESystem`)
+
+ - `AbstractTimeIndependentSystem`: has no independent variable (e.g.: `NonlinearSystem`)
+ - `AbstractTimeDependentSystem`: has a single independent variable (e.g.: `ODESystem`)
+ - `AbstractMultivariateSystem`: may have multiple independent variables (e.g.: `PDESystem`)
## Constructors and Naming
The `AbstractSystem` interface has a consistent method for constructing systems.
Generally, it follows the order of:
-1. Equations
-2. Independent Variables
-3. Dependent Variables (or States)
-4. Parameters
+ 1. Equations
+ 2. Independent Variables
+ 3. Dependent Variables (or States)
+ 4. Parameters
All other pieces are handled via keyword arguments. `AbstractSystem`s share the
same keyword arguments, which are:
-- `system`: This is used for specifying subsystems for hierarchical modeling with
- reusable components. For more information, see the [components page](@ref components).
-- Defaults: Keyword arguments like `defaults` are used for specifying default
- values which are used. If a value is not given at the `SciMLProblem` construction
- time, its numerical value will be the default.
+ - `system`: This is used for specifying subsystems for hierarchical modeling with
+ reusable components. For more information, see the [components page](@ref components).
+ - Defaults: Keyword arguments like `defaults` are used for specifying default
+ values which are used. If a value is not given at the `SciMLProblem` construction
+ time, its numerical value will be the default.
## Composition and Accessor Functions
@@ -43,32 +44,32 @@ total set, which includes that of all systems held inside.
The values which are common to all `AbstractSystem`s are:
-- `equations(sys)`: All equations that define the system and its subsystems.
-- `states(sys)`: All the states in the system and its subsystems.
-- `parameters(sys)`: All parameters of the system and its subsystems.
-- `nameof(sys)`: The name of the current-level system.
-- `get_eqs(sys)`: Equations that define the current-level system.
-- `get_states(sys)`: States that are in the current-level system.
-- `get_ps(sys)`: Parameters that are in the current-level system.
-- `get_systems(sys)`: Subsystems of the current-level system.
+ - `equations(sys)`: All equations that define the system and its subsystems.
+ - `states(sys)`: All the states in the system and its subsystems.
+ - `parameters(sys)`: All parameters of the system and its subsystems.
+ - `nameof(sys)`: The name of the current-level system.
+ - `get_eqs(sys)`: Equations that define the current-level system.
+ - `get_states(sys)`: States that are in the current-level system.
+ - `get_ps(sys)`: Parameters that are in the current-level system.
+ - `get_systems(sys)`: Subsystems of the current-level system.
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.
-- `get_metadata(sys)`: Any metadata about the system or its origin to be used by downstream packages.
+ - `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.
+ - `get_metadata(sys)`: Any metadata about the system or its origin to be used by downstream packages.
-Note that if you know a system is an `AbstractTimeDependentSystem` you could use `get_iv` to get the
+Note that if you know a system is an `AbstractTimeDependentSystem` you could use `get_iv` to get the
unique independent variable directly, rather than using `independent_variables(sys)[1]`, which is clunky and may cause problems if `sys` is an `AbstractMultivariateSystem` because there may be more than one independent variable. `AbstractTimeIndependentSystem`s do not have a method `get_iv`, and `independent_variables(sys)` will return a size-zero result for such. For an `AbstractMultivariateSystem`, `get_ivs` is equivalent.
A system could also have caches:
-- `get_jac(sys)`: The Jacobian of a system.
-- `get_tgrad(sys)`: The gradient with respect to time of a system.
+ - `get_jac(sys)`: The Jacobian of a system.
+ - `get_tgrad(sys)`: The gradient with respect to time of a system.
## Transformations
diff --git a/docs/src/basics/Composition.md b/docs/src/basics/Composition.md
index e81aaf35d4..75e2d4bf26 100644
--- a/docs/src/basics/Composition.md
+++ b/docs/src/basics/Composition.md
@@ -17,14 +17,14 @@ of `decay2` is the value of the state variable `x`.
```@example composition
using ModelingToolkit
-function decay(;name)
- @parameters t a
- @variables x(t) f(t)
- D = Differential(t)
- ODESystem([
- D(x) ~ -a*x + f
- ];
- name=name)
+function decay(; name)
+ @parameters t a
+ @variables x(t) f(t)
+ D = Differential(t)
+ ODESystem([
+ D(x) ~ -a * x + f,
+ ];
+ name = name)
end
@named decay1 = decay()
@@ -32,10 +32,8 @@ end
@parameters t
D = Differential(t)
-connected = compose(ODESystem([
- decay2.f ~ decay1.x
- D(decay1.f) ~ 0
- ], t; name=:connected), decay1, decay2)
+connected = compose(ODESystem([decay2.f ~ decay1.x
+ D(decay1.f) ~ 0], t; name = :connected), decay1, decay2)
equations(connected)
@@ -53,15 +51,11 @@ equations(simplified_sys)
Now we can solve the system:
```@example composition
-x0 = [
- decay1.x => 1.0
- decay1.f => 0.0
- decay2.x => 1.0
-]
-p = [
- decay1.a => 0.1
- decay2.a => 0.2
-]
+x0 = [decay1.x => 1.0
+ decay1.f => 0.0
+ decay2.x => 1.0]
+p = [decay1.a => 0.1
+ decay2.a => 0.2]
using DifferentialEquations
prob = ODEProblem(simplified_sys, x0, (0.0, 100.0), p)
@@ -76,7 +70,7 @@ subsystems. A model is the composition of itself and its subsystems.
For example, if we have:
```julia
-@named sys = compose(ODESystem(eqs,indepvar,states,ps),subsys)
+@named sys = compose(ODESystem(eqs, indepvar, states, ps), subsys)
```
the `equations` of `sys` is the concatenation of `get_eqs(sys)` and
@@ -101,10 +95,8 @@ this is done, the initial conditions and parameters must be specified
in their namespaced form. For example:
```julia
-u0 = [
- x => 2.0
- subsys.x => 2.0
-]
+u0 = [x => 2.0
+ subsys.x => 2.0]
```
Note that any default values within the given subcomponent will be
@@ -132,55 +124,51 @@ With symbolic parameters, it is possible to set the default value of a parameter
```julia
# ...
sys = ODESystem(
- # ...
- # directly in the defaults argument
- defaults=Pair{Num, Any}[
- x => u,
- y => σ,
- z => u-0.1,
-])
+ # ...
+ # directly in the defaults argument
+ defaults = Pair{Num, Any}[x => u,
+ y => σ,
+ z => u - 0.1])
# by assigning to the parameter
-sys.y = u*1.1
+sys.y = u * 1.1
```
In a hierarchical system, variables of the subsystem get namespaced by the name of the system they are in. This prevents naming clashes, but also enforces that every state and parameter is local to the subsystem it is used in. In some cases it might be desirable to have variables and parameters that are shared between subsystems, or even global. This can be accomplished as follows.
```julia
@parameters t a b c d e f
-p = [
- a #a is a local variable
- ParentScope(b) # b is a variable that belongs to one level up in the hierarchy
- ParentScope(ParentScope(c))# ParentScope can be nested
- DelayParentScope(d) # skips one level before applying ParentScope
- DelayParentScope(e,2) # second argument allows skipping N levels
- GlobalScope(f) # global variables will never be namespaced
-]
-
-level0 = ODESystem(Equation[],t,[],p; name = :level0)
-level1 = ODESystem(Equation[],t,[],[];name = :level1) ∘ level0
+p = [a #a is a local variable
+ ParentScope(b) # b is a variable that belongs to one level up in the hierarchy
+ ParentScope(ParentScope(c))# ParentScope can be nested
+ DelayParentScope(d) # skips one level before applying ParentScope
+ DelayParentScope(e, 2) # second argument allows skipping N levels
+ GlobalScope(f)]
+
+level0 = ODESystem(Equation[], t, [], p; name = :level0)
+level1 = ODESystem(Equation[], t, [], []; name = :level1) ∘ level0
parameters(level1)
- #level0₊a
- #b
- #c
- #level0₊d
- #level0₊e
- #f
-level2 = ODESystem(Equation[],t,[],[];name = :level2) ∘ level1
+#level0₊a
+#b
+#c
+#level0₊d
+#level0₊e
+#f
+level2 = ODESystem(Equation[], t, [], []; name = :level2) ∘ level1
parameters(level2)
- #level1₊level0₊a
- #level1₊b
- #c
- #level0₊d
- #level1₊level0₊e
- #f
-level3 = ODESystem(Equation[],t,[],[];name = :level3) ∘ level2
+#level1₊level0₊a
+#level1₊b
+#c
+#level0₊d
+#level1₊level0₊e
+#f
+level3 = ODESystem(Equation[], t, [], []; name = :level3) ∘ level2
parameters(level3)
- #level2₊level1₊level0₊a
- #level2₊level1₊b
- #level2₊c
- #level2₊level0₊d
- #level1₊level0₊e
- #f
+#level2₊level1₊level0₊a
+#level2₊level1₊b
+#level2₊c
+#level2₊level0₊d
+#level1₊level0₊e
+#f
```
## Structural Simplify
@@ -213,27 +201,25 @@ D = Differential(t)
@variables S(t), I(t), R(t)
N = S + I + R
-@parameters β,γ
+@parameters β, γ
-@named seqn = ODESystem([D(S) ~ -β*S*I/N])
-@named ieqn = ODESystem([D(I) ~ β*S*I/N-γ*I])
-@named reqn = ODESystem([D(R) ~ γ*I])
+@named seqn = ODESystem([D(S) ~ -β * S * I / N])
+@named ieqn = ODESystem([D(I) ~ β * S * I / N - γ * I])
+@named reqn = ODESystem([D(R) ~ γ * I])
sir = compose(ODESystem([
- S ~ ieqn.S,
- I ~ seqn.I,
- R ~ ieqn.R,
- ieqn.S ~ seqn.S,
- seqn.I ~ ieqn.I,
- seqn.R ~ reqn.R,
- ieqn.R ~ reqn.R,
- reqn.I ~ ieqn.I], t, [S,I,R], [β,γ],
- defaults = [
- seqn.β => β
- ieqn.β => β
- ieqn.γ => γ
- reqn.γ => γ
- ], name=:sir), seqn, ieqn, reqn)
+ S ~ ieqn.S,
+ I ~ seqn.I,
+ R ~ ieqn.R,
+ ieqn.S ~ seqn.S,
+ seqn.I ~ ieqn.I,
+ seqn.R ~ reqn.R,
+ ieqn.R ~ reqn.R,
+ reqn.I ~ ieqn.I], t, [S, I, R], [β, γ],
+ defaults = [seqn.β => β
+ ieqn.β => β
+ ieqn.γ => γ
+ reqn.γ => γ], name = :sir), seqn, ieqn, reqn)
```
Note that the states are forwarded by an equality relationship, while
@@ -251,17 +237,15 @@ equations(sireqn_simple)
## User Code
u0 = [seqn.S => 990.0,
- ieqn.I => 10.0,
- reqn.R => 0.0]
+ ieqn.I => 10.0,
+ reqn.R => 0.0]
-p = [
- β => 0.5
- γ => 0.25
-]
+p = [β => 0.5
+ γ => 0.25]
-tspan = (0.0,40.0)
-prob = ODEProblem(sireqn_simple,u0,tspan,p,jac=true)
-sol = solve(prob,Tsit5())
+tspan = (0.0, 40.0)
+prob = ODEProblem(sireqn_simple, u0, tspan, p, jac = true)
+sol = solve(prob, Tsit5())
sol[reqn.R]
```
@@ -277,6 +261,7 @@ 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
diff --git a/docs/src/basics/ContextualVariables.md b/docs/src/basics/ContextualVariables.md
index 663f4e7036..f4d2d1332f 100644
--- a/docs/src/basics/ContextualVariables.md
+++ b/docs/src/basics/ContextualVariables.md
@@ -23,17 +23,19 @@ to ignore such variables when attempting to find the states of a system.
## Constants
Constants are like parameters that:
-- always have a default value, which must be assigned when the constants are
+
+ - always have a default value, which must be assigned when the constants are
declared
-- do not show up in the list of parameters of a system.
+ - do not show up in the list of parameters of a system.
The intended use-cases for constants are:
-- representing literals (e.g., π) symbolically, which results in cleaner
+
+ - representing literals (e.g., π) symbolically, which results in cleaner
Latexification of equations (avoids turning `d ~ 2π*r` into `d = 6.283185307179586 r`)
-- allowing auto-generated unit conversion factors to live outside the list of
+ - allowing auto-generated unit conversion factors to live outside the list of
parameters
-- representing fundamental constants (e.g., speed of light `c`) that should never
- be adjusted inadvertently.
+ - representing fundamental constants (e.g., speed of light `c`) that should never
+ be adjusted inadvertently.
## Wildcard Variable Arguments
@@ -58,7 +60,7 @@ For example, the following specifies that `x` is a 2x2 matrix of flow variables
with the unit m^3/s:
```julia
-@variables x[1:2,1:2] [connect = Flow; unit = u"m^3/s"]
+@variables x[1:2, 1:2] [connect = Flow; unit = u"m^3/s"]
```
ModelingToolkit defines `connect`, `unit`, `noise`, and `description` keys for
diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md
index 592990ce2e..71f2b4c4ec 100644
--- a/docs/src/basics/Events.md
+++ b/docs/src/basics/Events.md
@@ -1,4 +1,5 @@
# [Event Handling and Callback Functions](@id events)
+
ModelingToolkit provides several ways to represent system events, which enable
system state or parameters to be changed when certain conditions are satisfied,
or can be used to detect discontinuities. These events are ultimately converted
@@ -25,45 +26,51 @@ functional affect](@ref func_affects) representation is also allowed, as describ
below.
## Continuous Events
+
The basic purely symbolic continuous event interface to encode *one* continuous
event is
+
```julia
AbstractSystem(eqs, ...; continuous_events::Vector{Equation})
AbstractSystem(eqs, ...; continuous_events::Pair{Vector{Equation}, Vector{Equation}})
```
+
In the former, equations that evaluate to 0 will represent conditions that should
be detected by the integrator, for example to force stepping to times of
discontinuities. The latter allow modeling of events that have an effect on the
state, where the first entry in the `Pair` is a vector of equations describing
event conditions, and the second vector of equations describes the effect on the
state. Each affect equation must be of the form
+
```julia
single_state_variable ~ expression_involving_any_variables_or_parameters
```
+
or
+
```julia
single_parameter ~ expression_involving_any_variables_or_parameters
```
+
In this basic interface, multiple variables can be changed in one event, or
multiple parameters, but not a mix of parameters and variables. The latter can
be handled via more [general functional affects](@ref func_affects).
-Finally, multiple events can be encoded via a `Vector{Pair{Vector{Equation},
-Vector{Equation}}}`.
+Finally, multiple events can be encoded via a `Vector{Pair{Vector{Equation}, Vector{Equation}}}`.
### Example: Friction
+
The system below illustrates how continuous events can be used to model Coulomb
friction
+
```@example events
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) # when v = 0 there is a discontinuity
+ @variables t x(t)=0 v(t)=0
+ D = Differential(t)
+ eqs = [D(x) ~ v
+ D(v) ~ sin(t) - k * sign(v)]
+ ODESystem(eqs, t; continuous_events = [v ~ 0], name) # when v = 0 there is a discontinuity
end
@named m = UnitMassWithFriction(0.7)
prob = ODEProblem(m, Pair[], (0, 10pi))
@@ -72,6 +79,7 @@ plot(sol)
```
### Example: Bouncing ball
+
In the documentation for
[DifferentialEquations](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#Example-1:-Bouncing-Ball),
we have an example where a bouncing ball is simulated using callbacks which have
@@ -83,79 +91,79 @@ like this
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
+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
+@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)
+tspan = (0.0, 5.0)
prob = ODEProblem(ball, Pair[], tspan)
-sol = solve(prob,Tsit5())
+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$.
+
```@example events
@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]
-]
+continuous_events = [[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)
-
+ D(x) ~ vx,
+ D(y) ~ vy,
+ D(vx) ~ -9.8 - 0.1vx, # gravity + some small air resistance
+ D(vy) ~ -0.1vy,
+ ], t; continuous_events)
ball = structural_simplify(ball)
-tspan = (0.0,10.0)
+tspan = (0.0, 10.0)
prob = ODEProblem(ball, Pair[], tspan)
-sol = solve(prob,Tsit5())
+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)
+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)
```
### [Generalized functional affect support](@id func_affects)
+
In some instances, a more flexible response to events is needed, which cannot be
encapsulated by symbolic equations. For example, a component may implement
complex behavior that is inconvenient or impossible to represent symbolically.
ModelingToolkit therefore supports regular Julia functions as affects: instead
of one or more equations, an affect is defined as a `tuple`:
+
```julia
[x ~ 0] => (affect!, [v, x], [p, q], ctx)
```
-where, `affect!` is a Julia function with the signature: `affect!(integ, u, p,
-ctx)`; `[u,v]` and `[p,q]` are the symbolic states (variables) and parameters
+
+where, `affect!` is a Julia function with the signature: `affect!(integ, u, p, ctx)`; `[u,v]` and `[p,q]` are the symbolic states (variables) and parameters
that are accessed by `affect!`, respectively; and `ctx` is any context that is
passed to `affect!` as the `ctx` argument.
`affect!` receives a [DifferentialEquations.jl
integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/)
- as its first argument, which can then be used to access states and parameters
+as its first argument, which can then be used to access states and parameters
that are provided in the `u` and `p` arguments (implemented as `NamedTuple`s).
The integrator can also be manipulated more generally to control solution
behavior, see the [integrator
interface](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/)
documentation. In affect functions, we have that
+
```julia
function affect!(integ, u, p, ctx)
# integ.t is the current time
@@ -163,16 +171,20 @@ function affect!(integ, u, p, ctx)
# integ.p[p.q] is the value of the parameter `q` above
end
```
+
When accessing variables of a sub-system, it can be useful to rename them
(alternatively, an affect function may be reused in different contexts):
+
```julia
[x ~ 0] => (affect!, [resistor₊v => :v, x], [p, q => :p2], ctx)
```
+
Here, the symbolic variable `resistor₊v` is passed as `v` while the symbolic
parameter `q` has been renamed `p2`.
As an example, here is the bouncing ball example from above using the functional
affect interface:
+
```@example events
sts = @variables x(t), v(t)
par = @parameters g = 9.8
@@ -198,11 +210,14 @@ plot(bb_sol)
```
## Discrete events support
+
In addition to continuous events, discrete events are also supported. The
general interface to represent a collection of discrete events is
+
```julia
AbstractSystem(eqs, ...; discrete_events = [condition1 => affect1, condition2 => affect2])
```
+
where conditions are symbolic expressions that should evaluate to `true` when an
individual affect should be executed. Here `affect1` and `affect2` are each
either a vector of one or more symbolic equations, or a functional affect, just
@@ -212,8 +227,10 @@ parameters, but one cannot currently mix state and parameter changes within one
individual event.
### Example: Injecting cells into a population
+
Suppose we have a population of `N(t)` cells that can grow and die, and at time
`t1` we want to inject `M` more cells into the population. We can model this by
+
```@example events
@parameters M tinject α
@variables t N(t)
@@ -241,6 +258,7 @@ Note that more general logical expressions can be built, for example, suppose we
want the event to occur at that time only if the solution is smaller than 50% of
its steady-state value (which is 100). We can encode this by modifying the event
to
+
```@example events
injection = ((t == tinject) & (N < 50)) => [N ~ N + M]
@@ -249,6 +267,7 @@ oprob = ODEProblem(osys, u0, tspan, p)
sol = solve(oprob, Tsit5(); tstops = 10.0)
plot(sol)
```
+
Since the solution is *not* smaller than half its steady-state value at the
event time, the event condition now returns false. Here we used logical and,
`&`, instead of the short-circuiting logical and, `&&`, as currently the latter
@@ -256,6 +275,7 @@ cannot be used within symbolic expressions.
Let's now also add a drug at time `tkill` that turns off production of new
cells, modeled by setting `α = 0.0`
+
```@example events
@parameters tkill
@@ -275,17 +295,21 @@ plot(sol)
```
### Periodic and preset-time events
+
Two important subclasses of discrete events are periodic and preset-time
events.
A preset-time event is triggered at specific set times, which can be
passed in a vector like
+
```julia
discrete_events = [[1.0, 4.0] => [v ~ -v]]
```
+
This will change the sign of `v` *only* at `t = 1.0` and `t = 4.0`.
As such, our last example with treatment and killing could instead be modeled by
+
```@example events
injection = [10.0] => [N ~ N + M]
killing = [20.0] => [α ~ 0.0]
@@ -297,19 +321,23 @@ oprob = ODEProblem(osys, u0, tspan, p)
sol = solve(oprob, Tsit5())
plot(sol)
```
+
Notice, one advantage of using a preset-time event is that one does not need to
also specify `tstops` in the call to solve.
A periodic event is triggered at fixed intervals (e.g. every Δt seconds). To
specify a periodic interval, pass the interval as the condition for the event.
For example,
+
```julia
-discrete_events=[1.0 => [v ~ -v]]
+discrete_events = [1.0 => [v ~ -v]]
```
+
will change the sign of `v` at `t = 1.0`, `2.0`, ...
Finally, we note that to specify an event at precisely one time, say 2.0 below,
one must still use a vector
+
```julia
discrete_events = [[2.0] => [v ~ -v]]
```
diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md
index b95a1abd61..4b4a0b39fe 100644
--- a/docs/src/basics/FAQ.md
+++ b/docs/src/basics/FAQ.md
@@ -8,8 +8,8 @@ from the solution. But what if you need to get the index? The following helper
function will do the trick:
```julia
-indexof(sym,syms) = findfirst(isequal(sym),syms)
-indexof(σ,parameters(sys))
+indexof(sym, syms) = findfirst(isequal(sym), syms)
+indexof(σ, parameters(sys))
```
## Transforming value maps to arrays
@@ -19,7 +19,7 @@ because symbol ordering is not guaranteed. However, what if you want to get the
lowered array? You can use the internal function `varmap_to_vars`. For example:
```julia
-pnew = varmap_to_vars([β=>3.0, c=>10.0, γ=>2.0],parameters(sys))
+pnew = varmap_to_vars([β => 3.0, c => 10.0, γ => 2.0], parameters(sys))
```
## How do I handle `if` statements in my symbolic forms?
@@ -35,7 +35,7 @@ term will be excluded from the model.
If you see the error:
-```julia
+```
ERROR: TypeError: non-boolean (Num) used in boolean context
```
@@ -57,7 +57,7 @@ For example, let's say you were building ODEProblems in the loss function like:
function loss(p)
prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
sol = solve(prob, Tsit5())
- sum(abs2,sol)
+ sum(abs2, sol)
end
```
@@ -69,9 +69,9 @@ once outside the loss function, and remake the prob inside the loss function:
```julia
prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
function loss(p)
- remake(prob,p = ...)
+ remake(prob, p = ...)
sol = solve(prob, Tsit5())
- sum(abs2,sol)
+ sum(abs2, sol)
end
```
@@ -91,8 +91,8 @@ Using this, the fixed index map can be used in the loss function. This would loo
prob = ODEProblem(sys, [], [p1 => p[1], p2 => p[2]])
idxs = Int.(ModelingToolkit.varmap_to_vars([p1 => 1, p2 => 2], p))
function loss(p)
- remake(prob,p = p[idxs])
+ remake(prob, p = p[idxs])
sol = solve(prob, Tsit5())
- sum(abs2,sol)
+ sum(abs2, sol)
end
```
diff --git a/docs/src/basics/Linearization.md b/docs/src/basics/Linearization.md
index f1a433c944..7c42bae30f 100644
--- a/docs/src/basics/Linearization.md
+++ b/docs/src/basics/Linearization.md
@@ -1,9 +1,13 @@
# [Linearization](@id linearization)
+
A nonlinear dynamical system with state (differential and algebraic) ``x`` and input signals ``u``
+
```math
M \dot x = f(x, u)
```
+
can be linearized using the function [`linearize`](@ref) to produce a linear statespace system on the form
+
```math
\begin{aligned}
\dot x &= Ax + Bu\\
@@ -14,6 +18,7 @@ y &= Cx + Du
The `linearize` function expects the user to specify the inputs ``u`` and the outputs ``u`` using the syntax shown in the example below:
## Example
+
```@example LINEARIZE
using ModelingToolkit
@variables t x(t)=0 y(t)=0 u(t)=0 r(t)=0
@@ -28,31 +33,37 @@ eqs = [u ~ kp * (r - y) # P controller
matrices, simplified_sys = linearize(sys, [r], [y]) # Linearize from r to y
matrices
```
+
The named tuple `matrices` contains the matrices of the linear statespace representation, while `simplified_sys` is an `ODESystem` that, among other things, indicates the state order in the linear system through
+
```@example LINEARIZE
using ModelingToolkit: inputs, outputs
[states(simplified_sys); inputs(simplified_sys); outputs(simplified_sys)]
```
## Operating point
-The operating point to linearize around can be specified with the keyword argument `op` like this: `op = Dict(x => 1, r => 2)`.
+
+The operating point to linearize around can be specified with the keyword argument `op` like this: `op = Dict(x => 1, r => 2)`.
## Batch linearization and algebraic variables
-If linearization is to be performed around multiple operating points, the simplification of the system has to be carried out a single time only. To facilitate this, the lower-level function [`ModelingToolkit.linearization_function`](@ref) is available. This function further allows you to obtain separate Jacobians for the differential and algebraic parts of the model. For ODE models without algebraic equations, the statespace representation above is available from the output of `linearization_function` as `A, B, C, D = f_x, f_u, h_x, h_u`.
+If linearization is to be performed around multiple operating points, the simplification of the system has to be carried out a single time only. To facilitate this, the lower-level function [`ModelingToolkit.linearization_function`](@ref) is available. This function further allows you to obtain separate Jacobians for the differential and algebraic parts of the model. For ODE models without algebraic equations, the statespace representation above is available from the output of `linearization_function` as `A, B, C, D = f_x, f_u, h_x, h_u`.
## Input derivatives
+
Physical systems are always *proper*, i.e., they do not differentiate causal inputs. However, ModelingToolkit allows you to model non-proper systems, such as inverse models, and may sometimes fail to find a realization of a proper system on proper form. In these situations, `linearize` may throw an error mentioning
+
```
Input derivatives appeared in expressions (-g_z\g_u != 0)
```
-This means that to simulate this system, some order of derivatives of the input is required. To allow `linearize` to proceed in this situation, one may pass the keyword argument `allow_input_derivatives = true`, in which case the resulting model will have twice as many inputs, ``2n_u``, where the last ``n_u`` inputs correspond to ``\dot u``.
-If the modeled system is actually proper (but MTK failed to find a proper realization), further numerical simplification can be applied to the resulting statespace system to obtain a proper form. Such simplification is currently available in the experimental package [ControlSystemsMTK](https://github.com/baggepinnen/ControlSystemsMTK.jl#internals-transformation-of-non-proper-models-to-proper-statespace-form).
+This means that to simulate this system, some order of derivatives of the input is required. To allow `linearize` to proceed in this situation, one may pass the keyword argument `allow_input_derivatives = true`, in which case the resulting model will have twice as many inputs, ``2n_u``, where the last ``n_u`` inputs correspond to ``\dot u``.
+If the modeled system is actually proper (but MTK failed to find a proper realization), further numerical simplification can be applied to the resulting statespace system to obtain a proper form. Such simplification is currently available in the experimental package [ControlSystemsMTK](https://github.com/baggepinnen/ControlSystemsMTK.jl#internals-transformation-of-non-proper-models-to-proper-statespace-form).
## Tools for linear analysis
-[ModelingToolkitStandardLibrary](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/) contains a set of [tools for more advanced linear analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/). These can be used to make it easier to work with and analyze causal models, such as control and signal-processing systems.
+
+[ModelingToolkitStandardLibrary](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/) contains a set of [tools for more advanced linear analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/). These can be used to make it easier to work with and analyze causal models, such as control and signal-processing systems.
```@index
Pages = ["Linearization.md"]
@@ -61,4 +72,4 @@ Pages = ["Linearization.md"]
```@docs
linearize
ModelingToolkit.linearization_function
-```
\ No newline at end of file
+```
diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md
index 7458feade9..77d3bc3461 100644
--- a/docs/src/basics/Validation.md
+++ b/docs/src/basics/Validation.md
@@ -1,166 +1,169 @@
-# [Model Validation and Units](@id units)
-
-ModelingToolkit.jl provides extensive functionality for model validation and unit checking. This is done by providing metadata to the variable types and then running the validation functions which identify malformed systems and non-physical equations. This approach provides high performance and compatibility with numerical solvers.
-
-## Assigning Units
-
-Units may be assigned with the following syntax.
-
-```@example validation
-using ModelingToolkit, Unitful
-@variables t [unit = u"s"] x(t) [unit = u"m"] g(t) w(t) [unit = "Hz"]
-
-@variables(t, [unit = u"s"], x(t), [unit = u"m"], g(t), w(t), [unit = "Hz"])
-
-@variables(begin
-t, [unit = u"s"],
-x(t), [unit = u"m"],
-g(t),
-w(t), [unit = "Hz"]
-end)
-
-# Simultaneously set default value (use plain numbers, not quantities)
-@variables x=10 [unit = u"m"]
-
-# Symbolic array: unit applies to all elements
-@variables x[1:3] [unit = u"m"]
-```
-
-Do not use `quantities` such as `1u"s"`, `1/u"s"` or `u"1/s"` as these will result in errors; instead use `u"s"`, `u"s^-1"`, or `u"s"^-1`.
-
-## Unit Validation & Inspection
-
-Unit validation of equations happens automatically when creating a system. However, for debugging purposes, one may wish to validate the equations directly using `validate`.
-
-```@docs
-ModelingToolkit.validate
-```
-
-Inside, `validate` uses `get_unit`, which may be directly applied to any term. Note that `validate` will not throw an error in the event of incompatible units, but `get_unit` will. If you would rather receive a warning instead of an error, use `safe_get_unit` which will yield `nothing` in the event of an error. Unit agreement is tested with `ModelingToolkit.equivalent(u1,u2)`.
-
-
-```@docs
-ModelingToolkit.get_unit
-```
-
-Example usage below. Note that `ModelingToolkit` does not force unit conversions to preferred units in the event of nonstandard combinations -- it merely checks that the equations are consistent.
-
-```@example validation
-using ModelingToolkit, Unitful
-@parameters τ [unit = u"ms"]
-@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
-D = Differential(t)
-eqs = eqs = [D(E) ~ P - E/τ,
- 0 ~ P ]
-ModelingToolkit.validate(eqs)
-```
-```@example validation
-ModelingToolkit.validate(eqs[1])
-```
-```@example validation
-ModelingToolkit.get_unit(eqs[1].rhs)
-```
-
-An example of an inconsistent system: at present, `ModelingToolkit` requires that the units of all terms in an equation or sum to be equal-valued (`ModelingToolkit.equivalent(u1,u2)`), rather than simply dimensionally consistent. In the future, the validation stage may be upgraded to support the insertion of conversion factors into the equations.
-
-```@example validation
-using ModelingToolkit, Unitful
-@parameters τ [unit = u"ms"]
-@variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
-D = Differential(t)
-eqs = eqs = [D(E) ~ P - E/τ,
- 0 ~ P ]
-ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
-```
-## User-Defined Registered Functions and Types
-
-In order to validate user-defined types and `register`ed functions, specialize `get_unit`. Single-parameter calls to `get_unit`
-expect an object type, while two-parameter calls expect a function type as the first argument, and a vector of arguments as the
-second argument.
-
-```@example validation2
-using ModelingToolkit, Unitful
-# Composite type parameter in registered function
-@parameters t
-D = Differential(t)
-struct NewType
- f
-end
-@register_symbolic dummycomplex(complex::Num, scalar)
-dummycomplex(complex, scalar) = complex.f - scalar
-
-c = NewType(1)
-ModelingToolkit.get_unit(x::NewType) = ModelingToolkit.get_unit(x.f)
-function ModelingToolkit.get_unit(op::typeof(dummycomplex),args)
- argunits = ModelingToolkit.get_unit.(args)
- ModelingToolkit.get_unit(-,args)
-end
-
-sts = @variables a(t)=0 [unit = u"cm"]
-ps = @parameters s=-1 [unit = u"cm"] c=c [unit = u"cm"]
-eqs = [D(a) ~ dummycomplex(c, s);]
-sys = ODESystem(eqs, t, [sts...;], [ps...;], name=:sys)
-sys_simple = structural_simplify(sys)
-```
-
-## `Unitful` Literals
-
-In order for a function to work correctly during both validation & execution, the function must be unit-agnostic. That is, no unitful literals may be used. Any unitful quantity must either be a `parameter` or `variable`. For example, these equations will not validate successfully.
-
-```julia
-using ModelingToolkit, Unitful
-@variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
-D = Differential(t)
-eqs = [D(E) ~ P - E/1u"ms" ]
-ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
-
-myfunc(E) = E/1u"ms"
-eqs = [D(E) ~ P - myfunc(E) ]
-ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
-```
-
-Instead, they should be parameterized:
-
-```@example validation3
-using ModelingToolkit, Unitful
-@parameters τ [unit = u"ms"]
-@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
-D = Differential(t)
-eqs = [D(E) ~ P - E/τ]
-ModelingToolkit.validate(eqs) #Returns true
-```
-```@example validation3
-myfunc(E,τ) = E/τ
-eqs = [D(E) ~ P - myfunc(E,τ)]
-ModelingToolkit.validate(eqs) #Returns true
-```
-
-It is recommended *not* to circumvent unit validation by specializing user-defined functions on `Unitful` arguments vs. `Numbers`. This both fails to take advantage of `validate` for ensuring correctness, and may cause in errors in the
-future when `ModelingToolkit` is extended to support eliminating `Unitful` literals from functions.
-
-## Other Restrictions
-
-`Unitful` provides non-scalar units such as `dBm`, `°C`, etc. At this time, `ModelingToolkit` only supports scalar quantities. Additionally, angular degrees (`°`) are not supported because trigonometric functions will treat plain numerical values as radians, which would lead systems validated using degrees to behave erroneously when being solved.
-
-## Troubleshooting & Gotchas
-
-If a system fails to validate due to unit issues, at least one warning message will appear, including a line number as well as the unit types and expressions that were in conflict. Some system constructors re-order equations before the unit checking can be done, in which case the equation numbers may be inaccurate. The printed expression that the problem resides in is always correctly shown.
-
-Symbolic exponents for unitful variables *are* supported (ex: `P^γ` in thermodynamics). However, this means that `ModelingToolkit` cannot reduce such expressions to `Unitful.Unitlike` subtypes at validation time because the exponent value is not available. In this case `ModelingToolkit.get_unit` is type-unstable, yielding a symbolic result, which can still be checked for symbolic equality with `ModelingToolkit.equivalent`.
-
-## Parameter & Initial Condition Values
-
-Parameter and initial condition values are supplied to problem constructors as plain numbers, with the understanding that they have been converted to the appropriate units. This is done for simplicity of interfacing with optimization solvers. Some helper function for dealing with value maps:
-
-```julia
-remove_units(p::Dict) = Dict(k => Unitful.ustrip(ModelingToolkit.get_unit(k),v) for (k,v) in p)
-add_units(p::Dict) = Dict(k => v*ModelingToolkit.get_unit(k) for (k,v) in p)
-```
-
-Recommended usage:
-
-```julia
-pars = @parameters τ [unit = u"ms"]
-p = Dict(τ => 1u"ms")
-ODEProblem(sys,remove_units(u0),tspan,remove_units(p))
-```
+# [Model Validation and Units](@id units)
+
+ModelingToolkit.jl provides extensive functionality for model validation and unit checking. This is done by providing metadata to the variable types and then running the validation functions which identify malformed systems and non-physical equations. This approach provides high performance and compatibility with numerical solvers.
+
+## Assigning Units
+
+Units may be assigned with the following syntax.
+
+```@example validation
+using ModelingToolkit, Unitful
+@variables t [unit = u"s"] x(t) [unit = u"m"] g(t) w(t) [unit = "Hz"]
+
+@variables(t, [unit = u"s"], x(t), [unit = u"m"], g(t), w(t), [unit = "Hz"])
+
+@variables(begin t, [unit = u"s"],
+ x(t), [unit = u"m"],
+ g(t),
+ w(t), [unit = "Hz"] end)
+
+# Simultaneously set default value (use plain numbers, not quantities)
+@variables x=10 [unit = u"m"]
+
+# Symbolic array: unit applies to all elements
+@variables x[1:3] [unit = u"m"]
+```
+
+Do not use `quantities` such as `1u"s"`, `1/u"s"` or `u"1/s"` as these will result in errors; instead use `u"s"`, `u"s^-1"`, or `u"s"^-1`.
+
+## Unit Validation & Inspection
+
+Unit validation of equations happens automatically when creating a system. However, for debugging purposes, one may wish to validate the equations directly using `validate`.
+
+```@docs
+ModelingToolkit.validate
+```
+
+Inside, `validate` uses `get_unit`, which may be directly applied to any term. Note that `validate` will not throw an error in the event of incompatible units, but `get_unit` will. If you would rather receive a warning instead of an error, use `safe_get_unit` which will yield `nothing` in the event of an error. Unit agreement is tested with `ModelingToolkit.equivalent(u1,u2)`.
+
+```@docs
+ModelingToolkit.get_unit
+```
+
+Example usage below. Note that `ModelingToolkit` does not force unit conversions to preferred units in the event of nonstandard combinations -- it merely checks that the equations are consistent.
+
+```@example validation
+using ModelingToolkit, Unitful
+@parameters τ [unit = u"ms"]
+@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
+D = Differential(t)
+eqs = eqs = [D(E) ~ P - E / τ,
+ 0 ~ P]
+ModelingToolkit.validate(eqs)
+```
+
+```@example validation
+ModelingToolkit.validate(eqs[1])
+```
+
+```@example validation
+ModelingToolkit.get_unit(eqs[1].rhs)
+```
+
+An example of an inconsistent system: at present, `ModelingToolkit` requires that the units of all terms in an equation or sum to be equal-valued (`ModelingToolkit.equivalent(u1,u2)`), rather than simply dimensionally consistent. In the future, the validation stage may be upgraded to support the insertion of conversion factors into the equations.
+
+```@example validation
+using ModelingToolkit, Unitful
+@parameters τ [unit = u"ms"]
+@variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
+D = Differential(t)
+eqs = eqs = [D(E) ~ P - E / τ,
+ 0 ~ P]
+ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
+```
+
+## User-Defined Registered Functions and Types
+
+In order to validate user-defined types and `register`ed functions, specialize `get_unit`. Single-parameter calls to `get_unit`
+expect an object type, while two-parameter calls expect a function type as the first argument, and a vector of arguments as the
+second argument.
+
+```@example validation2
+using ModelingToolkit, Unitful
+# Composite type parameter in registered function
+@parameters t
+D = Differential(t)
+struct NewType
+ f::Any
+end
+@register_symbolic dummycomplex(complex::Num, scalar)
+dummycomplex(complex, scalar) = complex.f - scalar
+
+c = NewType(1)
+ModelingToolkit.get_unit(x::NewType) = ModelingToolkit.get_unit(x.f)
+function ModelingToolkit.get_unit(op::typeof(dummycomplex), args)
+ argunits = ModelingToolkit.get_unit.(args)
+ ModelingToolkit.get_unit(-, args)
+end
+
+sts = @variables a(t)=0 [unit = u"cm"]
+ps = @parameters s=-1 [unit = u"cm"] c=c [unit = u"cm"]
+eqs = [D(a) ~ dummycomplex(c, s);]
+sys = ODESystem(eqs, t, [sts...;], [ps...;], name = :sys)
+sys_simple = structural_simplify(sys)
+```
+
+## `Unitful` Literals
+
+In order for a function to work correctly during both validation & execution, the function must be unit-agnostic. That is, no unitful literals may be used. Any unitful quantity must either be a `parameter` or `variable`. For example, these equations will not validate successfully.
+
+```julia
+using ModelingToolkit, Unitful
+@variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
+D = Differential(t)
+eqs = [D(E) ~ P - E / 1u"ms"]
+ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
+
+myfunc(E) = E / 1u"ms"
+eqs = [D(E) ~ P - myfunc(E)]
+ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
+```
+
+Instead, they should be parameterized:
+
+```@example validation3
+using ModelingToolkit, Unitful
+@parameters τ [unit = u"ms"]
+@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
+D = Differential(t)
+eqs = [D(E) ~ P - E / τ]
+ModelingToolkit.validate(eqs) #Returns true
+```
+
+```@example validation3
+myfunc(E, τ) = E / τ
+eqs = [D(E) ~ P - myfunc(E, τ)]
+ModelingToolkit.validate(eqs) #Returns true
+```
+
+It is recommended *not* to circumvent unit validation by specializing user-defined functions on `Unitful` arguments vs. `Numbers`. This both fails to take advantage of `validate` for ensuring correctness, and may cause in errors in the
+future when `ModelingToolkit` is extended to support eliminating `Unitful` literals from functions.
+
+## Other Restrictions
+
+`Unitful` provides non-scalar units such as `dBm`, `°C`, etc. At this time, `ModelingToolkit` only supports scalar quantities. Additionally, angular degrees (`°`) are not supported because trigonometric functions will treat plain numerical values as radians, which would lead systems validated using degrees to behave erroneously when being solved.
+
+## Troubleshooting & Gotchas
+
+If a system fails to validate due to unit issues, at least one warning message will appear, including a line number as well as the unit types and expressions that were in conflict. Some system constructors re-order equations before the unit checking can be done, in which case the equation numbers may be inaccurate. The printed expression that the problem resides in is always correctly shown.
+
+Symbolic exponents for unitful variables *are* supported (ex: `P^γ` in thermodynamics). However, this means that `ModelingToolkit` cannot reduce such expressions to `Unitful.Unitlike` subtypes at validation time because the exponent value is not available. In this case `ModelingToolkit.get_unit` is type-unstable, yielding a symbolic result, which can still be checked for symbolic equality with `ModelingToolkit.equivalent`.
+
+## Parameter & Initial Condition Values
+
+Parameter and initial condition values are supplied to problem constructors as plain numbers, with the understanding that they have been converted to the appropriate units. This is done for simplicity of interfacing with optimization solvers. Some helper function for dealing with value maps:
+
+```julia
+function remove_units(p::Dict)
+ Dict(k => Unitful.ustrip(ModelingToolkit.get_unit(k), v) for (k, v) in p)
+end
+add_units(p::Dict) = Dict(k => v * ModelingToolkit.get_unit(k) for (k, v) in p)
+```
+
+Recommended usage:
+
+```julia
+pars = @parameters τ [unit = u"ms"]
+p = Dict(τ => 1u"ms")
+ODEProblem(sys, remove_units(u0), tspan, remove_units(p))
+```
diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/basics/Variable_metadata.md
index 55c4121261..096927ed1d 100644
--- a/docs/src/basics/Variable_metadata.md
+++ b/docs/src/basics/Variable_metadata.md
@@ -1,10 +1,13 @@
# Symbolic metadata
+
It is possible to add metadata to symbolic variables, the metadata will be displayed when calling help on a variable.
The following information can be added (note, it's possible to extend this to user-defined metadata as well)
## Variable descriptions
+
Descriptive strings can be attached to variables using the `[description = "descriptive string"]` syntax:
+
```@example metadata
using ModelingToolkit
@variables u [description = "This is my input"]
@@ -12,16 +15,18 @@ getdescription(u)
```
When variables with descriptions are present in systems, they will be printed when the system is shown in the terminal:
+
```@example metadata
@parameters t
@variables u(t) [description = "A short description of u"]
-@parameters p [description = "A description of p"]
+@parameters p [description = "A description of p"]
@named sys = ODESystem([u ~ p], t)
show(stdout, "text/plain", sys) # hide
```
Calling help on the variable `u` displays the description, alongside other metadata:
-```julia
+
+```
help?> u
A variable of type Symbolics.Num (Num wraps anything in a type that is a subtype of Real)
@@ -35,91 +40,103 @@ help?> u
```
## Input or output
+
Designate a variable as either an input or an output using the following
+
```@example metadata
using ModelingToolkit
-@variables u [input=true]
+@variables u [input = true]
isinput(u)
```
+
```@example metadata
-@variables y [output=true]
+@variables y [output = true]
isoutput(y)
```
## Bounds
+
Bounds are useful when parameters are to be optimized, or to express intervals of uncertainty.
```@example metadata
-@variables u [bounds=(-1,1)]
+@variables u [bounds = (-1, 1)]
hasbounds(u)
```
+
```@example metadata
getbounds(u)
```
-## Mark input as a disturbance
+## Mark input as a disturbance
+
Indicate that an input is not available for control, i.e., it's a disturbance input.
```@example metadata
-@variables u [input=true, disturbance=true]
+@variables u [input = true, disturbance = true]
isdisturbance(u)
```
## Mark parameter as tunable
+
Indicate that a parameter can be automatically tuned by parameter optimization or automatic control tuning apps.
```@example metadata
-@parameters Kp [tunable=true]
+@parameters Kp [tunable = true]
istunable(Kp)
```
## Probability distributions
+
A probability distribution may be associated with a parameter to indicate either
uncertainty about its value, or as a prior distribution for Bayesian optimization.
```julia
using Distributions
d = Normal(10, 1)
-@parameters m [dist=d]
+@parameters m [dist = d]
hasdist(m)
```
+
```julia
getdist(m)
```
## Additional functions
+
For systems that contain parameters with metadata like described above, have some additional functions defined for convenience.
In the example below, we define a system with tunable parameters and extract bounds vectors
```@example metadata
@parameters t
Dₜ = Differential(t)
-@variables x(t)=0 u(t)=0 [input=true] y(t)=0 [output=true]
+@variables x(t)=0 u(t)=0 [input = true] y(t)=0 [output = true]
@parameters T [tunable = true, bounds = (0, Inf)]
@parameters k [tunable = true, bounds = (0, Inf)]
-eqs = [
- Dₜ(x) ~ (-x + k*u) / T # A first-order system with time constant T and gain k
- y ~ x
-]
-sys = ODESystem(eqs, t, name=:tunable_first_order)
+eqs = [Dₜ(x) ~ (-x + k * u) / T # A first-order system with time constant T and gain k
+ y ~ x]
+sys = ODESystem(eqs, t, name = :tunable_first_order)
```
+
```@example metadata
p = tunable_parameters(sys) # extract all parameters marked as tunable
```
+
```@example metadata
lb, ub = getbounds(p) # operating on a vector, we get lower and upper bound vectors
```
+
```@example metadata
b = getbounds(sys) # Operating on the system, we get a dict
```
-
## Index
+
```@index
Pages = ["Variable_metadata.md"]
```
## Docstrings
+
```@autodocs
Modules = [ModelingToolkit]
Pages = ["variables.jl"]
diff --git a/docs/src/comparison.md b/docs/src/comparison.md
index 8afe0d8dcf..76a39b0199 100644
--- a/docs/src/comparison.md
+++ b/docs/src/comparison.md
@@ -2,118 +2,118 @@
## Comparison Against Modelica
-- Both Modelica and ModelingToolkit.jl are acausal modeling languages.
-- Modelica is a language with many different implementations, such as
- [Dymola](https://www.3ds.com/products-services/catia/products/dymola/) and
- [OpenModelica](https://openmodelica.org/), which have differing levels of
- performance and can give different results on the same model. Many of the
- commonly used Modelica compilers are not open-source. ModelingToolkit.jl
- is a language with a single canonical open-source implementation.
-- All current Modelica compiler implementations are fixed and not extendable
- by the users from the Modelica language itself. For example, the Dymola
- compiler [shares its symbolic processing pipeline](https://www.claytex.com/tech-blog/model-translation-and-symbolic-manipulation/),
- which is roughly equivalent to the `dae_index_lowering` and `structural_simplify`
- of ModelingToolkit.jl. ModelingToolkit.jl is an open and hackable transformation
- system which allows users to add new non-standard transformations and control
- the order of application.
-- Modelica is a declarative programming language. ModelingToolkit.jl is a
- declarative symbolic modeling language used from within the Julia programming
- language. Its programming language semantics, such as loop constructs and
- conditionals, can be used to more easily generate models.
-- Modelica is an object-oriented single dispatch language. ModelingToolkit.jl,
- built on Julia, uses multiple dispatch extensively to simplify code.
-- Many Modelica compilers supply a GUI. ModelingToolkit.jl does not.
-- Modelica can be used to simulate ODE and DAE systems. ModelingToolkit.jl
- has a much more expansive set of system types, including nonlinear systems,
- SDEs, PDEs, and more.
+ - Both Modelica and ModelingToolkit.jl are acausal modeling languages.
+ - Modelica is a language with many different implementations, such as
+ [Dymola](https://www.3ds.com/products-services/catia/products/dymola/) and
+ [OpenModelica](https://openmodelica.org/), which have differing levels of
+ performance and can give different results on the same model. Many of the
+ commonly used Modelica compilers are not open-source. ModelingToolkit.jl
+ is a language with a single canonical open-source implementation.
+ - All current Modelica compiler implementations are fixed and not extendable
+ by the users from the Modelica language itself. For example, the Dymola
+ compiler [shares its symbolic processing pipeline](https://www.claytex.com/tech-blog/model-translation-and-symbolic-manipulation/),
+ which is roughly equivalent to the `dae_index_lowering` and `structural_simplify`
+ of ModelingToolkit.jl. ModelingToolkit.jl is an open and hackable transformation
+ system which allows users to add new non-standard transformations and control
+ the order of application.
+ - Modelica is a declarative programming language. ModelingToolkit.jl is a
+ declarative symbolic modeling language used from within the Julia programming
+ language. Its programming language semantics, such as loop constructs and
+ conditionals, can be used to more easily generate models.
+ - Modelica is an object-oriented single dispatch language. ModelingToolkit.jl,
+ built on Julia, uses multiple dispatch extensively to simplify code.
+ - Many Modelica compilers supply a GUI. ModelingToolkit.jl does not.
+ - Modelica can be used to simulate ODE and DAE systems. ModelingToolkit.jl
+ has a much more expansive set of system types, including nonlinear systems,
+ SDEs, PDEs, and more.
## Comparison Against Simulink
-- Simulink is a causal modeling environment, whereas ModelingToolkit.jl is an
- acausal modeling environment. For an overview of the differences, consult
- academic reviews such as [this one](https://arxiv.org/abs/1909.00484). In this
- sense, ModelingToolkit.jl is more similar to the Simscape sub-environment.
-- Simulink is used from MATLAB while ModelingToolkit.jl is used from Julia.
- Thus any user-defined functions have the performance of their host language.
- For information on the performance differences between Julia and MATLAB,
- consult [open benchmarks](https://julialang.org/benchmarks/), which demonstrate
- Julia as an order of magnitude or more faster in many cases due to its JIT
- compilation.
-- Simulink uses the MATLAB differential equation solvers, while ModelingToolkit.jl
- uses [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/). For a systematic
- comparison between the solvers, consult
- [open benchmarks](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/),
- which demonstrate two orders of magnitude performance advantage for the native
- Julia solvers across many benchmark problems.
-- Simulink comes with a Graphical User Interface (GUI), ModelingToolkit.jl
- does not.
-- Simulink is a proprietary software, meaning users cannot actively modify or
- extend the software. ModelingToolkit.jl is built in Julia and used in Julia,
- where users can actively extend and modify the software interactively in the
- REPL and contribute to its open-source repositories.
-- Simulink covers ODE and DAE systems. ModelingToolkit.jl has a much more
- expansive set of system types, including SDEs, PDEs, optimization problems,
- and more.
+ - Simulink is a causal modeling environment, whereas ModelingToolkit.jl is an
+ acausal modeling environment. For an overview of the differences, consult
+ academic reviews such as [this one](https://arxiv.org/abs/1909.00484). In this
+ sense, ModelingToolkit.jl is more similar to the Simscape sub-environment.
+ - Simulink is used from MATLAB while ModelingToolkit.jl is used from Julia.
+ Thus any user-defined functions have the performance of their host language.
+ For information on the performance differences between Julia and MATLAB,
+ consult [open benchmarks](https://julialang.org/benchmarks/), which demonstrate
+ Julia as an order of magnitude or more faster in many cases due to its JIT
+ compilation.
+ - Simulink uses the MATLAB differential equation solvers, while ModelingToolkit.jl
+ uses [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/). For a systematic
+ comparison between the solvers, consult
+ [open benchmarks](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/),
+ which demonstrate two orders of magnitude performance advantage for the native
+ Julia solvers across many benchmark problems.
+ - Simulink comes with a Graphical User Interface (GUI), ModelingToolkit.jl
+ does not.
+ - Simulink is a proprietary software, meaning users cannot actively modify or
+ extend the software. ModelingToolkit.jl is built in Julia and used in Julia,
+ where users can actively extend and modify the software interactively in the
+ REPL and contribute to its open-source repositories.
+ - Simulink covers ODE and DAE systems. ModelingToolkit.jl has a much more
+ expansive set of system types, including SDEs, PDEs, optimization problems,
+ and more.
## Comparison Against CASADI
-- CASADI is written in C++ but used from Python/MATLAB, meaning that it cannot be
- directly extended by users unless they are using the C++ interface and run a
- local build of CASADI. ModelingToolkit.jl is both written and used from
- Julia, meaning that users can easily extend the library on the fly, even
- interactively in the REPL.
-- CASADI includes limited support for Computer Algebra System (CAS) functionality,
- while ModelingToolkit.jl is built on the full
- [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/) CAS.
-- CASADI supports DAE and ODE problems via SUNDIALS IDAS and CVODES. ModelingToolkit.jl
- supports DAE and ODE problems via [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/),
- of which Sundials.jl is <1% of the total available solvers and is outperformed
- by the native Julia solvers on the vast majority of the benchmark equations.
- In addition, the DifferentialEquations.jl interface is confederated, meaning
- that any user can dynamically extend the system to add new solvers to the
- interface by defining new dispatches of solve.
-- CASADI's DAEBuilder does not implement efficiency transformations like tearing,
- which are standard in the ModelingToolkit.jl transformation pipeline.
-- CASADI supports special functionality for quadratic programming problems, while
- ModelingToolkit only provides nonlinear programming via `OptimizationSystem`.
-- ModelingToolkit.jl integrates with its host language Julia, so Julia code
- can be automatically converted into ModelingToolkit expressions. Users of
- CASADI must explicitly create CASADI expressions.
+ - CASADI is written in C++ but used from Python/MATLAB, meaning that it cannot be
+ directly extended by users unless they are using the C++ interface and run a
+ local build of CASADI. ModelingToolkit.jl is both written and used from
+ Julia, meaning that users can easily extend the library on the fly, even
+ interactively in the REPL.
+ - CASADI includes limited support for Computer Algebra System (CAS) functionality,
+ while ModelingToolkit.jl is built on the full
+ [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/) CAS.
+ - CASADI supports DAE and ODE problems via SUNDIALS IDAS and CVODES. ModelingToolkit.jl
+ supports DAE and ODE problems via [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/),
+ of which Sundials.jl is <1% of the total available solvers and is outperformed
+ by the native Julia solvers on the vast majority of the benchmark equations.
+ In addition, the DifferentialEquations.jl interface is confederated, meaning
+ that any user can dynamically extend the system to add new solvers to the
+ interface by defining new dispatches of solve.
+ - CASADI's DAEBuilder does not implement efficiency transformations like tearing,
+ which are standard in the ModelingToolkit.jl transformation pipeline.
+ - CASADI supports special functionality for quadratic programming problems, while
+ ModelingToolkit only provides nonlinear programming via `OptimizationSystem`.
+ - ModelingToolkit.jl integrates with its host language Julia, so Julia code
+ can be automatically converted into ModelingToolkit expressions. Users of
+ CASADI must explicitly create CASADI expressions.
## Comparison Against Modia.jl
-- Modia.jl uses Julia's expression objects for representing its equations.
- ModelingToolkit.jl uses [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/),
- and thus the Julia expressions follow Julia semantics and can be manipulated
- using a computer algebra system (CAS).
-- Modia's compilation pipeline is similar to the
- [Dymola symbolic processing pipeline](https://www.claytex.com/tech-blog/model-translation-and-symbolic-manipulation/)
- with some improvements. ModelingToolkit.jl has an open transformation pipeline
- that allows for users to extend and reorder transformation passes, where
- `structural_simplify` is an adaptation of the Modia.jl-improved alias elimination
- and tearing algorithms.
-- Both Modia and ModelingToolkit generate `DAEProblem` and `ODEProblem` forms for
- solving with [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/).
-- ModelingToolkit.jl integrates with its host language Julia, so Julia code
- can be automatically converted into ModelingToolkit expressions. Users of
- Modia must explicitly create Modia expressions.
-- Modia covers DAE systems. ModelingToolkit.jl has a much more
- expansive set of system types, including SDEs, PDEs, optimization problems,
- and more.
+ - Modia.jl uses Julia's expression objects for representing its equations.
+ ModelingToolkit.jl uses [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/),
+ and thus the Julia expressions follow Julia semantics and can be manipulated
+ using a computer algebra system (CAS).
+ - Modia's compilation pipeline is similar to the
+ [Dymola symbolic processing pipeline](https://www.claytex.com/tech-blog/model-translation-and-symbolic-manipulation/)
+ with some improvements. ModelingToolkit.jl has an open transformation pipeline
+ that allows for users to extend and reorder transformation passes, where
+ `structural_simplify` is an adaptation of the Modia.jl-improved alias elimination
+ and tearing algorithms.
+ - Both Modia and ModelingToolkit generate `DAEProblem` and `ODEProblem` forms for
+ solving with [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/).
+ - ModelingToolkit.jl integrates with its host language Julia, so Julia code
+ can be automatically converted into ModelingToolkit expressions. Users of
+ Modia must explicitly create Modia expressions.
+ - Modia covers DAE systems. ModelingToolkit.jl has a much more
+ expansive set of system types, including SDEs, PDEs, optimization problems,
+ and more.
## Comparison Against Causal.jl
-- Causal.jl is a causal modeling environment, whereas ModelingToolkit.jl is an
- acausal modeling environment. For an overview of the differences, consult
- academic reviews such as [this one](https://arxiv.org/abs/1909.00484).
-- Both ModelingToolkit.jl and Causal.jl use [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/)
- as the backend solver library.
-- Causal.jl lets one add arbitrary equation systems to a given node, and allow
- the output to effect the next node. This means an SDE may drive an ODE. These
- two portions are solved with different solver methods in tandem. In
- ModelingToolkit.jl, such connections promote the whole system to an SDE. This
- results in better accuracy and stability, though in some cases it can be
- less performant.
-- Causal.jl, similar to Simulink, breaks algebraic loops via inexact heuristics.
- ModelingToolkit.jl treats algebraic loops exactly through algebraic equations
- in the generated model.
+ - Causal.jl is a causal modeling environment, whereas ModelingToolkit.jl is an
+ acausal modeling environment. For an overview of the differences, consult
+ academic reviews such as [this one](https://arxiv.org/abs/1909.00484).
+ - Both ModelingToolkit.jl and Causal.jl use [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/)
+ as the backend solver library.
+ - Causal.jl lets one add arbitrary equation systems to a given node, and allow
+ the output to effect the next node. This means an SDE may drive an ODE. These
+ two portions are solved with different solver methods in tandem. In
+ ModelingToolkit.jl, such connections promote the whole system to an SDE. This
+ results in better accuracy and stability, though in some cases it can be
+ less performant.
+ - Causal.jl, similar to Simulink, breaks algebraic loops via inexact heuristics.
+ ModelingToolkit.jl treats algebraic loops exactly through algebraic equations
+ in the generated model.
diff --git a/docs/src/examples/higher_order.md b/docs/src/examples/higher_order.md
index 2e467548ad..f181efd839 100644
--- a/docs/src/examples/higher_order.md
+++ b/docs/src/examples/higher_order.md
@@ -1,59 +1,60 @@
-# Automatic Transformation of Nth Order ODEs to 1st Order ODEs
-
-ModelingToolkit has a system for transformations of mathematical
-systems. These transformations allow for symbolically changing
-the representation of the model to problems that are easier to
-numerically solve. One simple to demonstrate transformation is the
-`structural_simplify` with does a lot of tricks, one being the
-transformation that sends a Nth order ODE
-to a 1st order ODE.
-
-To see this, let's define a second order riff on the Lorenz equations.
-We utilize the derivative operator twice here to define the second order:
-
-```@example orderlowering
-using ModelingToolkit, OrdinaryDiffEq
-
-@parameters σ ρ β
-@variables t x(t) y(t) z(t)
-D = Differential(t)
-
-eqs = [D(D(x)) ~ σ*(y-x),
- D(y) ~ x*(ρ-z)-y,
- D(z) ~ x*y - β*z]
-
-@named sys = ODESystem(eqs)
-```
-
-Note that we could've used an alternative syntax for 2nd order, i.e.
-`D = Differential(t)^2` and then `E(x)` would be the second derivative,
-and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compose
-`Differential`s, like `Differential(t) * Differential(x)`.
-
-Now let's transform this into the `ODESystem` of first order components.
-We do this by simply calling `ode_order_lowering`:
-
-```@example orderlowering
-sys = structural_simplify(sys)
-```
-
-Now we can directly numerically solve the lowered system. Note that,
-following the original problem, the solution requires knowing the
-initial condition for `x'`, and thus we include that in our input
-specification:
-
-```@example orderlowering
-u0 = [D(x) => 2.0,
- x => 1.0,
- y => 0.0,
- z => 0.0]
-
-p = [σ => 28.0,
- ρ => 10.0,
- β => 8/3]
-
-tspan = (0.0,100.0)
-prob = ODEProblem(sys,u0,tspan,p,jac=true)
-sol = solve(prob,Tsit5())
-using Plots; plot(sol,idxs=(x,y))
-```
+# Automatic Transformation of Nth Order ODEs to 1st Order ODEs
+
+ModelingToolkit has a system for transformations of mathematical
+systems. These transformations allow for symbolically changing
+the representation of the model to problems that are easier to
+numerically solve. One simple to demonstrate transformation is the
+`structural_simplify` with does a lot of tricks, one being the
+transformation that sends a Nth order ODE
+to a 1st order ODE.
+
+To see this, let's define a second order riff on the Lorenz equations.
+We utilize the derivative operator twice here to define the second order:
+
+```@example orderlowering
+using ModelingToolkit, OrdinaryDiffEq
+
+@parameters σ ρ β
+@variables t x(t) y(t) z(t)
+D = Differential(t)
+
+eqs = [D(D(x)) ~ σ * (y - x),
+ D(y) ~ x * (ρ - z) - y,
+ D(z) ~ x * y - β * z]
+
+@named sys = ODESystem(eqs)
+```
+
+Note that we could've used an alternative syntax for 2nd order, i.e.
+`D = Differential(t)^2` and then `E(x)` would be the second derivative,
+and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compose
+`Differential`s, like `Differential(t) * Differential(x)`.
+
+Now let's transform this into the `ODESystem` of first order components.
+We do this by simply calling `ode_order_lowering`:
+
+```@example orderlowering
+sys = structural_simplify(sys)
+```
+
+Now we can directly numerically solve the lowered system. Note that,
+following the original problem, the solution requires knowing the
+initial condition for `x'`, and thus we include that in our input
+specification:
+
+```@example orderlowering
+u0 = [D(x) => 2.0,
+ x => 1.0,
+ y => 0.0,
+ z => 0.0]
+
+p = [σ => 28.0,
+ ρ => 10.0,
+ β => 8 / 3]
+
+tspan = (0.0, 100.0)
+prob = ODEProblem(sys, u0, tspan, p, jac = true)
+sol = solve(prob, Tsit5())
+using Plots;
+plot(sol, idxs = (x, y));
+```
diff --git a/docs/src/examples/modelingtoolkitize_index_reduction.md b/docs/src/examples/modelingtoolkitize_index_reduction.md
index 3d2f8d11a2..21cc650e70 100644
--- a/docs/src/examples/modelingtoolkitize_index_reduction.md
+++ b/docs/src/examples/modelingtoolkitize_index_reduction.md
@@ -17,13 +17,13 @@ function pendulum!(du, u, p, t)
x, dx, y, dy, T = u
g, L = p
du[1] = dx
- du[2] = T*x
+ du[2] = T * x
du[3] = dy
- du[4] = T*y - g
+ du[4] = T * y - g
du[5] = x^2 + y^2 - L^2
return nothing
end
-pendulum_fun! = ODEFunction(pendulum!, mass_matrix=Diagonal([1,1,1,1,0]))
+pendulum_fun! = ODEFunction(pendulum!, mass_matrix = Diagonal([1, 1, 1, 1, 0]))
u0 = [1.0, 0, 0, 0, 0]
p = [9.8, 1]
tspan = (0, 10.0)
@@ -31,8 +31,8 @@ pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
traced_sys = modelingtoolkitize(pendulum_prob)
pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
prob = ODAEProblem(pendulum_sys, [], tspan)
-sol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)
-plot(sol, idxs=states(traced_sys))
+sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8)
+plot(sol, idxs = states(traced_sys))
```
## Explanation
@@ -60,19 +60,23 @@ using OrdinaryDiffEq, LinearAlgebra
function pendulum!(du, u, p, t)
x, dx, y, dy, T = u
g, L = p
- du[1] = dx; du[2] = T*x
- du[3] = dy; du[4] = T*y - g
+ du[1] = dx
+ du[2] = T * x
+ du[3] = dy
+ du[4] = T * y - g
du[5] = x^2 + y^2 - L^2
end
-pendulum_fun! = ODEFunction(pendulum!, mass_matrix=Diagonal([1,1,1,1,0]))
-u0 = [1.0, 0, 0, 0, 0]; p = [9.8, 1]; tspan = (0, 10.0)
+pendulum_fun! = ODEFunction(pendulum!, mass_matrix = Diagonal([1, 1, 1, 1, 0]))
+u0 = [1.0, 0, 0, 0, 0];
+p = [9.8, 1];
+tspan = (0, 10.0);
pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
-solve(pendulum_prob,Rodas4())
+solve(pendulum_prob, Rodas4())
```
However, one will quickly be greeted with the unfortunate message:
-```julia
+```
┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq C:\Users\accou\.julia\packages\OrdinaryDiffEq\yCczp\src\initdt.jl:76
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
@@ -150,7 +154,7 @@ prob = ODEProblem(pendulum_sys, Pair[], tspan)
sol = solve(prob, Rodas4())
using Plots
-plot(sol, idxs=states(traced_sys))
+plot(sol, idxs = states(traced_sys))
```
Note that plotting using `states(traced_sys)` is done so that any
@@ -168,8 +172,8 @@ be solved via an explicit Runge-Kutta method:
traced_sys = modelingtoolkitize(pendulum_prob)
pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
prob = ODAEProblem(pendulum_sys, Pair[], tspan)
-sol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)
-plot(sol, idxs=states(traced_sys))
+sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8)
+plot(sol, idxs = states(traced_sys))
```
And there you go: this has transformed the model from being too hard to
diff --git a/docs/src/examples/parsing.md b/docs/src/examples/parsing.md
index 1c38998181..fb7c064127 100644
--- a/docs/src/examples/parsing.md
+++ b/docs/src/examples/parsing.md
@@ -14,7 +14,7 @@ ex = [:(y ~ x)
We can use the function `parse_expr_to_symbolic` from Symbolics.jl to generate the symbolic
form of the expression:
-```@example parsing
+```@example parsing
using Symbolics
eqs = parse_expr_to_symbolic.(ex, (Main,))
```
@@ -27,6 +27,6 @@ using ModelingToolkit, NonlinearSolve
vars = union(ModelingToolkit.vars.(eqs)...)
@named ns = NonlinearSystem(eqs, vars, [])
-prob = NonlinearProblem(ns,[1.0,1.0,1.0])
-sol = solve(prob,NewtonRaphson())
-```
\ No newline at end of file
+prob = NonlinearProblem(ns, [1.0, 1.0, 1.0])
+sol = solve(prob, NewtonRaphson())
+```
diff --git a/docs/src/examples/perturbation.md b/docs/src/examples/perturbation.md
index 6f9cfe9c9c..d978f25925 100644
--- a/docs/src/examples/perturbation.md
+++ b/docs/src/examples/perturbation.md
@@ -7,15 +7,15 @@ In the previous tutorial, [Mixed Symbolic-Numeric Perturbation Theory](https://s
```julia
using Symbolics, SymbolicUtils
-def_taylor(x, ps) = sum([a*x^i for (i,a) in enumerate(ps)])
+def_taylor(x, ps) = sum([a * x^i for (i, a) in enumerate(ps)])
def_taylor(x, ps, p₀) = p₀ + def_taylor(x, ps)
-function collect_powers(eq, x, ns; max_power=100)
- eq = substitute(expand(eq), Dict(x^j => 0 for j=last(ns)+1:max_power))
+function collect_powers(eq, x, ns; max_power = 100)
+ eq = substitute(expand(eq), Dict(x^j => 0 for j in (last(ns) + 1):max_power))
eqs = []
for i in ns
- powers = Dict(x^j => (i==j ? 1 : 0) for j=1:last(ns))
+ powers = Dict(x^j => (i == j ? 1 : 0) for j in 1:last(ns))
push!(eqs, substitute(eq, powers))
end
eqs
@@ -24,7 +24,7 @@ end
function solve_coef(eqs, ps)
vals = Dict()
- for i = 1:length(ps)
+ for i in 1:length(ps)
eq = substitute(eqs[i], vals)
vals[ps[i]] = Symbolics.solve_for(eq ~ 0, ps[i])
end
@@ -64,7 +64,7 @@ We need the second derivative of `x`. It may seem that we can do this using `Dif
as the second derivative of `x`. After rearrangement, our governing equation is $\ddot{x}(t)(1 + \epsilon x(t))^{-2} + 1 = 0$, or
```julia
-eq = ∂∂x * (1 + ϵ*x)^2 + 1
+eq = ∂∂x * (1 + ϵ * x)^2 + 1
```
The next two steps are the same as the ones for algebraic equations (note that we pass `1:n` to `collect_powers` because the zeroth order term is needed here)
@@ -100,17 +100,17 @@ states(sys)
```julia
# the initial conditions
# everything is zero except the initial velocity
-u0 = zeros(2n+2)
+u0 = zeros(2n + 2)
u0[3] = 1.0 # y₀ˍt
prob = ODEProblem(sys, u0, (0, 3.0))
-sol = solve(prob; dtmax=0.01)
+sol = solve(prob; dtmax = 0.01)
```
Finally, we calculate the solution to the problem as a function of `ϵ` by substituting the solution to the ODE system back into the defining equation for `x`. Note that `𝜀` is a number, compared to `ϵ`, which is a symbolic variable.
```julia
-X = 𝜀 -> sum([𝜀^(i-1) * sol[y[i]] for i in eachindex(y)])
+X = 𝜀 -> sum([𝜀^(i - 1) * sol[y[i]] for i in eachindex(y)])
```
Using `X`, we can plot the trajectory for a range of $𝜀$s.
@@ -118,7 +118,7 @@ Using `X`, we can plot the trajectory for a range of $𝜀$s.
```julia
using Plots
-plot(sol.t, hcat([X(𝜀) for 𝜀 = 0.0:0.1:0.5]...))
+plot(sol.t, hcat([X(𝜀) for 𝜀 in 0.0:0.1:0.5]...))
```
As expected, as `𝜀` becomes larger (meaning the gravity is less with altitude), the object goes higher and stays up for a longer duration. Of course, we could have solved the problem directly using as ODE solver. One of the benefits of the perturbation method is that we need to run the ODE solver only once and then can just calculate the answer for different values of `𝜀`; whereas, if we had used the direct method, we would need to run the solver once for each value of `𝜀`.
@@ -140,7 +140,7 @@ x = def_taylor(ϵ, y[3:end], y[2])
This time we also need the first derivative terms. Continuing,
```julia
-eq = ∂∂x + 2*ϵ*∂x + x
+eq = ∂∂x + 2 * ϵ * ∂x + x
eqs = collect_powers(eq, ϵ, 0:n)
vals = solve_coef(eqs, ∂∂y)
```
@@ -164,15 +164,15 @@ sys = structural_simplify(sys)
```julia
# the initial conditions
-u0 = zeros(2n+2)
+u0 = zeros(2n + 2)
u0[3] = 1.0 # y₀ˍt
prob = ODEProblem(sys, u0, (0, 50.0))
-sol = solve(prob; dtmax=0.01)
+sol = solve(prob; dtmax = 0.01)
-X = 𝜀 -> sum([𝜀^(i-1) * sol[y[i]] for i in eachindex(y)])
+X = 𝜀 -> sum([𝜀^(i - 1) * sol[y[i]] for i in eachindex(y)])
T = sol.t
-Y = 𝜀 -> exp.(-𝜀*T) .* sin.(sqrt(1 - 𝜀^2)*T) / sqrt(1 - 𝜀^2) # exact solution
+Y = 𝜀 -> exp.(-𝜀 * T) .* sin.(sqrt(1 - 𝜀^2) * T) / sqrt(1 - 𝜀^2) # exact solution
plot(sol.t, [Y(0.1), X(0.1)])
```
diff --git a/docs/src/examples/sparse_jacobians.md b/docs/src/examples/sparse_jacobians.md
index 92aca4c079..070112ee01 100644
--- a/docs/src/examples/sparse_jacobians.md
+++ b/docs/src/examples/sparse_jacobians.md
@@ -14,37 +14,41 @@ partial differential equation discretized using finite differences:
using DifferentialEquations, ModelingToolkit
const N = 32
-const xyd_brusselator = range(0,stop=1,length=N)
-brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.
-limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a
+const xyd_brusselator = range(0, stop = 1, length = N)
+brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0
+limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p, t)
- A, B, alpha, dx = p
- alpha = alpha/dx^2
- @inbounds for I in CartesianIndices((N, N))
- i, j = Tuple(I)
- x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
- ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
- du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
- B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
- du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
- A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
+ A, B, alpha, dx = p
+ alpha = alpha / dx^2
+ @inbounds for I in CartesianIndices((N, N))
+ i, j = Tuple(I)
+ x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
+ ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N),
+ limit(j - 1, N)
+ du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] -
+ 4u[i, j, 1]) +
+ B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] +
+ brusselator_f(x, y, t)
+ du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] -
+ 4u[i, j, 2]) +
+ A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
end
end
-p = (3.4, 1., 10., step(xyd_brusselator))
+p = (3.4, 1.0, 10.0, step(xyd_brusselator))
function init_brusselator_2d(xyd)
- N = length(xyd)
- u = zeros(N, N, 2)
- for I in CartesianIndices((N, N))
- x = xyd[I[1]]
- y = xyd[I[2]]
- u[I,1] = 22*(y*(1-y))^(3/2)
- u[I,2] = 27*(x*(1-x))^(3/2)
- end
- u
+ N = length(xyd)
+ u = zeros(N, N, 2)
+ for I in CartesianIndices((N, N))
+ x = xyd[I[1]]
+ y = xyd[I[2]]
+ u[I, 1] = 22 * (y * (1 - y))^(3 / 2)
+ u[I, 2] = 27 * (x * (1 - x))^(3 / 2)
+ end
+ u
end
u0 = init_brusselator_2d(xyd_brusselator)
-prob = ODEProblem(brusselator_2d_loop,u0,(0.,11.5),p)
+prob = ODEProblem(brusselator_2d_loop, u0, (0.0, 11.5), p)
```
Now let's use `modelingtoolkitize` to generate the symbolic version:
@@ -58,18 +62,19 @@ Now we regenerate the problem using `jac=true` for the analytical Jacobian
and `sparse=true` to make it sparse:
```@example sparsejac
-sparseprob = ODEProblem(sys,Pair[],(0.,11.5),jac=true,sparse=true)
+sparseprob = ODEProblem(sys, Pair[], (0.0, 11.5), jac = true, sparse = true)
```
Hard? No! How much did that help?
```@example sparsejac
using BenchmarkTools
-@btime solve(prob,save_everystep=false);
+@btime solve(prob, save_everystep = false);
return nothing # hide
```
+
```@example sparsejac
-@btime solve(sparseprob,save_everystep=false);
+@btime solve(sparseprob, save_everystep = false);
return nothing # hide
```
@@ -78,7 +83,7 @@ Thus in some cases we may only want to get the sparsity pattern. In this case,
we can simply do:
```@example sparsejac
-sparsepatternprob = ODEProblem(sys,Pair[],(0.,11.5),sparse=true)
-@btime solve(sparsepatternprob,save_everystep=false);
+sparsepatternprob = ODEProblem(sys, Pair[], (0.0, 11.5), sparse = true)
+@btime solve(sparsepatternprob, save_everystep = false);
return nothing # hide
```
diff --git a/docs/src/examples/spring_mass.md b/docs/src/examples/spring_mass.md
index 081e1cfdd7..771b181442 100644
--- a/docs/src/examples/spring_mass.md
+++ b/docs/src/examples/spring_mass.md
@@ -11,58 +11,58 @@ using Symbolics: scalarize
@variables t
D = Differential(t)
-function Mass(; name, m = 1.0, xy = [0., 0.], u = [0., 0.])
- ps = @parameters m=m
+function Mass(; name, m = 1.0, xy = [0.0, 0.0], u = [0.0, 0.0])
+ ps = @parameters m = m
sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u
eqs = scalarize(D.(pos) .~ v)
ODESystem(eqs, t, [pos..., v...], ps; name)
end
-function Spring(; name, k = 1e4, l = 1.)
+function Spring(; name, k = 1e4, l = 1.0)
ps = @parameters k=k l=l
@variables x(t), dir(t)[1:2]
ODESystem(Equation[], t, [x, dir...], ps; name)
end
function connect_spring(spring, a, b)
- [
- spring.x ~ norm(scalarize(a .- b))
- scalarize(spring.dir .~ scalarize(a .- b))
- ]
+ [spring.x ~ norm(scalarize(a .- b))
+ scalarize(spring.dir .~ scalarize(a .- b))]
end
-spring_force(spring) = -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l) ./ spring.x
+function spring_force(spring)
+ -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l) ./ spring.x
+end
m = 1.0
-xy = [1., -1.]
+xy = [1.0, -1.0]
k = 1e4
-l = 1.
-center = [0., 0.]
-g = [0., -9.81]
-@named mass = Mass(m=m, xy=xy)
-@named spring = Spring(k=k, l=l)
+l = 1.0
+center = [0.0, 0.0]
+g = [0.0, -9.81]
+@named mass = Mass(m = m, xy = xy)
+@named spring = Spring(k = k, l = l)
-eqs = [
- connect_spring(spring, mass.pos, center)
- scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)
-]
+eqs = [connect_spring(spring, mass.pos, center)
+ scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)]
@named _model = ODESystem(eqs, t, [spring.x; spring.dir; mass.pos], [])
@named model = compose(_model, mass, spring)
sys = structural_simplify(model)
-prob = ODEProblem(sys, [], (0., 3.))
+prob = ODEProblem(sys, [], (0.0, 3.0))
sol = solve(prob, Rosenbrock23())
plot(sol)
```
## Explanation
+
### Building the components
+
For each component, we use a Julia function that returns an `ODESystem`. At the top, we define the fundamental properties of a `Mass`: it has a mass `m`, a position `pos` and a velocity `v`. We also define that the velocity is the rate of change of position with respect to time.
```@example component
-function Mass(; name, m = 1.0, xy = [0., 0.], u = [0., 0.])
- ps = @parameters m=m
+function Mass(; name, m = 1.0, xy = [0.0, 0.0], u = [0.0, 0.0])
+ ps = @parameters m = m
sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u
eqs = scalarize(D.(pos) .~ v)
ODESystem(eqs, t, [pos..., v...], ps; name)
@@ -84,7 +84,7 @@ Or using the `@named` helper macro
Next, we build the spring component. It is characterized by the spring constant `k` and the length `l` of the spring when no force is applied to it. The state of a spring is defined by its current length and direction.
```@example component
-function Spring(; name, k = 1e4, l = 1.)
+function Spring(; name, k = 1e4, l = 1.0)
ps = @parameters k=k l=l
@variables x(t), dir(t)[1:2]
ODESystem(Equation[], t, [x, dir...], ps; name)
@@ -95,39 +95,37 @@ We now define functions that help construct the equations for a mass-spring syst
```@example component
function connect_spring(spring, a, b)
- [
- spring.x ~ norm(scalarize(a .- b))
- scalarize(spring.dir .~ scalarize(a .- b))
- ]
+ [spring.x ~ norm(scalarize(a .- b))
+ scalarize(spring.dir .~ scalarize(a .- b))]
end
```
Lastly, we define the `spring_force` function that takes a `spring` and returns the force exerted by this spring.
```@example component
-spring_force(spring) = -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l) ./ spring.x
+function spring_force(spring)
+ -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l) ./ spring.x
+end
```
To create our system, we will first create the components: a mass and a spring. This is done as follows:
```@example component
m = 1.0
-xy = [1., -1.]
+xy = [1.0, -1.0]
k = 1e4
-l = 1.
-center = [0., 0.]
-g = [0., -9.81]
-@named mass = Mass(m=m, xy=xy)
-@named spring = Spring(k=k, l=l)
+l = 1.0
+center = [0.0, 0.0]
+g = [0.0, -9.81]
+@named mass = Mass(m = m, xy = xy)
+@named spring = Spring(k = k, l = l)
```
We can now create the equations describing this system, by connecting `spring` to `mass` and a fixed point.
```@example component
-eqs = [
- connect_spring(spring, mass.pos, center)
- scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)
-]
+eqs = [connect_spring(spring, mass.pos, center)
+ scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)]
```
Finally, we can build the model using these equations and components.
@@ -168,6 +166,7 @@ We are left with only 4 equations involving 4 state variables (`mass.pos[1]`,
`mass.pos[2]`, `mass.v[1]`, `mass.v[2]`). We can solve the system by converting
it to an `ODEProblem`. Some observed variables are not expanded by default. To
view the complete equations, one can do
+
```@example component
full_equations(sys)
```
@@ -175,7 +174,7 @@ full_equations(sys)
This is done as follows:
```@example component
-prob = ODEProblem(sys, [], (0., 3.))
+prob = ODEProblem(sys, [], (0.0, 3.0))
sol = solve(prob, Rosenbrock23())
plot(sol)
```
diff --git a/docs/src/examples/tearing_parallelism.md b/docs/src/examples/tearing_parallelism.md
index e2f11925a0..ca68dd71f2 100644
--- a/docs/src/examples/tearing_parallelism.md
+++ b/docs/src/examples/tearing_parallelism.md
@@ -15,95 +15,81 @@ using ModelingToolkit, OrdinaryDiffEq
# Basic electric components
@variables t
const D = Differential(t)
-@connector function Pin(;name)
+@connector function Pin(; name)
@variables v(t)=1.0 i(t)=1.0 [connect = Flow]
- ODESystem(Equation[], t, [v, i], [], name=name)
+ ODESystem(Equation[], t, [v, i], [], name = name)
end
-function Ground(;name)
+function Ground(; name)
@named g = Pin()
eqs = [g.v ~ 0]
- compose(ODESystem(eqs, t, [], [], name=name), g)
+ compose(ODESystem(eqs, t, [], [], name = name), g)
end
-function ConstantVoltage(;name, V = 1.0)
+function ConstantVoltage(; name, V = 1.0)
val = V
@named p = Pin()
@named n = Pin()
- @parameters V=V
- eqs = [
- V ~ p.v - n.v
- 0 ~ p.i + n.i
- ]
- compose(ODESystem(eqs, t, [], [V], name=name), p, n)
+ @parameters V = V
+ eqs = [V ~ p.v - n.v
+ 0 ~ p.i + n.i]
+ compose(ODESystem(eqs, t, [], [V], name = name), p, n)
end
-@connector function HeatPort(;name)
+@connector function HeatPort(; name)
@variables T(t)=293.15 Q_flow(t)=0.0 [connect = Flow]
- ODESystem(Equation[], t, [T, Q_flow], [], name=name)
+ ODESystem(Equation[], t, [T, Q_flow], [], name = name)
end
-function HeatingResistor(;name, R=1.0, TAmbient=293.15, alpha=1.0)
+function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0)
@named p = Pin()
@named n = Pin()
@named h = HeatPort()
@variables v(t) RTherm(t)
@parameters R=R TAmbient=TAmbient alpha=alpha
- eqs = [
- RTherm ~ R*(1 + alpha*(h.T - TAmbient))
+ eqs = [RTherm ~ R * (1 + alpha * (h.T - TAmbient))
v ~ p.i * RTherm
h.Q_flow ~ -v * p.i # -LossPower
v ~ p.v - n.v
- 0 ~ p.i + n.i
- ]
- compose(ODESystem(
- eqs, t, [v, RTherm], [R, TAmbient, alpha],
- name=name,
- ), p, n, h)
+ 0 ~ p.i + n.i]
+ compose(ODESystem(eqs, t, [v, RTherm], [R, TAmbient, alpha],
+ name = name), p, n, h)
end
-function HeatCapacitor(;name, rho=8050, V=1, cp=460, TAmbient=293.15)
+function HeatCapacitor(; name, rho = 8050, V = 1, cp = 460, TAmbient = 293.15)
@parameters rho=rho V=V cp=cp
- C = rho*V*cp
+ C = rho * V * cp
@named h = HeatPort()
eqs = [
- D(h.T) ~ h.Q_flow / C
- ]
- compose(ODESystem(
- eqs, t, [], [rho, V, cp],
- name=name,
- ), h)
+ D(h.T) ~ h.Q_flow / C,
+ ]
+ compose(ODESystem(eqs, t, [], [rho, V, cp],
+ name = name), h)
end
-function Capacitor(;name, C = 1.0)
+function Capacitor(; name, C = 1.0)
@named p = Pin()
@named n = Pin()
- @variables v(t)=0.0
- @parameters C=C
- eqs = [
- v ~ p.v - n.v
+ @variables v(t) = 0.0
+ @parameters C = C
+ eqs = [v ~ p.v - n.v
0 ~ p.i + n.i
- D(v) ~ p.i / C
- ]
- compose(ODESystem(
- eqs, t, [v], [C],
- name=name
- ), p, n)
+ D(v) ~ p.i / C]
+ compose(ODESystem(eqs, t, [v], [C],
+ name = name), p, n)
end
function parallel_rc_model(i; name, source, ground, R, C)
- resistor = HeatingResistor(name=Symbol(:resistor, i), R=R)
- capacitor = Capacitor(name=Symbol(:capacitor, i), C=C)
- heat_capacitor = HeatCapacitor(name=Symbol(:heat_capacitor, i))
+ resistor = HeatingResistor(name = Symbol(:resistor, i), R = R)
+ capacitor = Capacitor(name = Symbol(:capacitor, i), C = C)
+ heat_capacitor = HeatCapacitor(name = Symbol(:heat_capacitor, i))
- rc_eqs = [
- connect(source.p, resistor.p)
+ rc_eqs = [connect(source.p, resistor.p)
connect(resistor.n, capacitor.p)
connect(capacitor.n, source.n, ground.g)
- connect(resistor.h, heat_capacitor.h)
- ]
+ connect(resistor.h, heat_capacitor.h)]
- compose(ODESystem(rc_eqs, t, name=Symbol(name, i)),
+ compose(ODESystem(rc_eqs, t, name = Symbol(name, i)),
[resistor, capacitor, source, ground, heat_capacitor])
end
```
@@ -116,18 +102,19 @@ we can connect a bunch of RC components as follows:
```@example tearing
V = 2.0
-@named source = ConstantVoltage(V=V)
+@named source = ConstantVoltage(V = V)
@named ground = Ground()
N = 50
-Rs = 10 .^range(0, stop=-4, length=N)
-Cs = 10 .^range(-3, stop=0, length=N)
+Rs = 10 .^ range(0, stop = -4, length = N)
+Cs = 10 .^ range(-3, stop = 0, length = N)
rc_systems = map(1:N) do i
- parallel_rc_model(i; name=:rc, source=source, ground=ground, R=Rs[i], C=Cs[i])
+ parallel_rc_model(i; name = :rc, source = source, ground = ground, R = Rs[i], C = Cs[i])
end;
-@variables E(t)=0.0
+@variables E(t) = 0.0
eqs = [
- D(E) ~ sum(((i, sys),)->getproperty(sys, Symbol(:resistor, i)).h.Q_flow, enumerate(rc_systems))
- ]
+ D(E) ~ sum(((i, sys),) -> getproperty(sys, Symbol(:resistor, i)).h.Q_flow,
+ enumerate(rc_systems)),
+]
@named _big_rc = ODESystem(eqs, t, [E], [])
@named big_rc = compose(_big_rc, rc_systems)
```
@@ -166,8 +153,11 @@ investigate what this means:
using ModelingToolkit.BipartiteGraphs
ts = TearingState(expand_connections(big_rc))
inc_org = BipartiteGraphs.incidence_matrix(ts.structure.graph)
-blt_org = StructuralTransformations.sorted_incidence_matrix(ts, only_algeqs=true, only_algvars=true)
-blt_reduced = StructuralTransformations.sorted_incidence_matrix(ModelingToolkit.get_tearing_state(sys), only_algeqs=true, only_algvars=true)
+blt_org = StructuralTransformations.sorted_incidence_matrix(ts, only_algeqs = true,
+ only_algvars = true)
+blt_reduced = StructuralTransformations.sorted_incidence_matrix(ModelingToolkit.get_tearing_state(sys),
+ only_algeqs = true,
+ only_algvars = true)
```
diff --git a/docs/src/index.md b/docs/src/index.md
index b518b211f8..05c7055f32 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -46,19 +46,19 @@ before generating code.
### Feature List
-- Causal and acausal modeling (Simulink/Modelica)
-- Automated model transformation, simplification, and composition
-- Automatic conversion of numerical models into symbolic models
-- Composition of models through the components, a lazy connection system, and
- tools for expanding/flattening
-- Pervasive parallelism in symbolic computations and generated functions
-- Transformations like alias elimination and tearing of nonlinear systems for
- efficiently numerically handling large-scale systems of equations
-- The ability to use the entire Symbolics.jl Computer Algebra System (CAS) as
- part of the modeling process.
-- Import models from common formats like SBML, CellML, BioNetGen, and more.
-- Extensibility: the whole system is written in pure Julia, so adding new
- functions, simplification rules, and model transformations has no barrier.
+ - Causal and acausal modeling (Simulink/Modelica)
+ - Automated model transformation, simplification, and composition
+ - Automatic conversion of numerical models into symbolic models
+ - Composition of models through the components, a lazy connection system, and
+ tools for expanding/flattening
+ - Pervasive parallelism in symbolic computations and generated functions
+ - Transformations like alias elimination and tearing of nonlinear systems for
+ efficiently numerically handling large-scale systems of equations
+ - The ability to use the entire Symbolics.jl Computer Algebra System (CAS) as
+ part of the modeling process.
+ - Import models from common formats like SBML, CellML, BioNetGen, and more.
+ - Extensibility: the whole system is written in pure Julia, so adding new
+ functions, simplification rules, and model transformations has no barrier.
For information on how to use the Symbolics.jl CAS system that ModelingToolkit.jl
is built on, consult the
@@ -66,36 +66,40 @@ is built on, consult the
### Equation Types
-- Ordinary differential equations
-- Stochastic differential equations
-- Partial differential equations
-- Nonlinear systems
-- Optimization problems
-- Continuous-Time Markov Chains
-- Chemical Reactions (via [Catalyst.jl](https://docs.sciml.ai/Catalyst/stable/))
-- Nonlinear Optimal Control
+ - Ordinary differential equations
+ - Stochastic differential equations
+ - Partial differential equations
+ - Nonlinear systems
+ - Optimization problems
+ - Continuous-Time Markov Chains
+ - Chemical Reactions (via [Catalyst.jl](https://docs.sciml.ai/Catalyst/stable/))
+ - Nonlinear Optimal Control
## Standard Library
-For quick development, ModelingToolkit.jl includes
+For quick development, ModelingToolkit.jl includes
[ModelingToolkitStandardLibrary.jl](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/),
a standard library of prebuilt components for the ModelingToolkit ecosystem.
## Model Import Formats
-- [CellMLToolkit.jl](https://docs.sciml.ai/CellMLToolkit/stable/): Import [CellML](https://www.cellml.org/) models into ModelingToolkit
- - Repository of more than a thousand pre-made models
- - Focus on biomedical models in areas such as: Calcium Dynamics,
- Cardiovascular Circulation, Cell Cycle, Cell Migration, Circadian Rhythms,
- Electrophysiology, Endocrine, Excitation-Contraction Coupling, Gene Regulation,
- Hepatology, Immunology, Ion Transport, Mechanical Constitutive Laws,
- Metabolism, Myofilament Mechanics, Neurobiology, pH Regulation, PKPD,
- Protein Modules, Signal Transduction, and Synthetic Biology.
-- [SBMLToolkit.jl](https://docs.sciml.ai/SBMLToolkit/stable/): Import [SBML](http://sbml.org/) models into ModelingToolkit
- - Uses the robust libsbml library for parsing and transforming the SBML
-- [ReactionNetworkImporters.jl](https://docs.sciml.ai/ReactionNetworkImporters/stable/): Import various models into ModelingToolkit
- - Supports the BioNetGen `.net` file
- - Supports importing networks specified by stoichiometric matrices
+ - [CellMLToolkit.jl](https://docs.sciml.ai/CellMLToolkit/stable/): Import [CellML](https://www.cellml.org/) models into ModelingToolkit
+
+ + Repository of more than a thousand pre-made models
+ + Focus on biomedical models in areas such as: Calcium Dynamics,
+ Cardiovascular Circulation, Cell Cycle, Cell Migration, Circadian Rhythms,
+ Electrophysiology, Endocrine, Excitation-Contraction Coupling, Gene Regulation,
+ Hepatology, Immunology, Ion Transport, Mechanical Constitutive Laws,
+ Metabolism, Myofilament Mechanics, Neurobiology, pH Regulation, PKPD,
+ Protein Modules, Signal Transduction, and Synthetic Biology.
+
+ - [SBMLToolkit.jl](https://docs.sciml.ai/SBMLToolkit/stable/): Import [SBML](http://sbml.org/) models into ModelingToolkit
+
+ + Uses the robust libsbml library for parsing and transforming the SBML
+ - [ReactionNetworkImporters.jl](https://docs.sciml.ai/ReactionNetworkImporters/stable/): Import various models into ModelingToolkit
+
+ + Supports the BioNetGen `.net` file
+ + Supports importing networks specified by stoichiometric matrices
## Extension Libraries
@@ -103,33 +107,40 @@ Because ModelingToolkit.jl is the core foundation of an equation-based modeling
ecosystem, there is a large set of libraries adding features to this system.
Below is an incomplete list of extension libraries one may want to be aware of:
-- [Catalyst.jl](https://docs.sciml.ai/Catalyst/stable/): Symbolic representations
- of chemical reactions
- - Symbolically build and represent large systems of chemical reactions
- - Generate code for ODEs, SDEs, continuous-time Markov Chains, and more
- - Simulate the models using the SciML ecosystem with O(1) Gillespie methods
-- [DataDrivenDiffEq.jl](https://docs.sciml.ai/DataDrivenDiffEq/stable/): Automatic
- identification of equations from data
- - Automated construction of ODEs and DAEs from data
- - Representations of Koopman operators and Dynamic Mode Decomposition (DMD)
-- [MomentClosure.jl](https://docs.sciml.ai/MomentClosure/dev/): Automatic
- transformation of ReactionSystems into deterministic systems
- - Generates ODESystems for the moment closures
- - Allows for geometrically-distributed random reaction rates
-- [ReactionMechanismSimulator.jl](https://docs.sciml.ai/ReactionMechanismSimulator/stable):
- Simulating and analyzing large chemical reaction mechanisms
- - Ideal gas and dilute liquid phases.
- - Constant T and P and constant V adiabatic ideal gas reactors.
- - Constant T and V dilute liquid reactors.
- - Diffusion limited rates. Sensitivity analysis for all reactors.
- - Flux diagrams with molecular images (if molecular information is provided).
-- [NumCME.jl](https://github.com/voduchuy/NumCME.jl): High-performance simulation of chemical master equations (CME)
- - Transient solution of the CME
- - Dynamic state spaces
- - Accepts reaction systems defined using Catalyst.jl DSL.
-- [FiniteStateProjection.jl](https://github.com/kaandocal/FiniteStateProjection.jl): High-performance simulation of
- chemical master equations (CME) via finite state projections
- - Accepts reaction systems defined using Catalyst.jl DSL.
+ - [Catalyst.jl](https://docs.sciml.ai/Catalyst/stable/): Symbolic representations
+ of chemical reactions
+
+ + Symbolically build and represent large systems of chemical reactions
+ + Generate code for ODEs, SDEs, continuous-time Markov Chains, and more
+ + Simulate the models using the SciML ecosystem with O(1) Gillespie methods
+
+ - [DataDrivenDiffEq.jl](https://docs.sciml.ai/DataDrivenDiffEq/stable/): Automatic
+ identification of equations from data
+
+ + Automated construction of ODEs and DAEs from data
+ + Representations of Koopman operators and Dynamic Mode Decomposition (DMD)
+ - [MomentClosure.jl](https://docs.sciml.ai/MomentClosure/dev/): Automatic
+ transformation of ReactionSystems into deterministic systems
+
+ + Generates ODESystems for the moment closures
+ + Allows for geometrically-distributed random reaction rates
+ - [ReactionMechanismSimulator.jl](https://docs.sciml.ai/ReactionMechanismSimulator/stable):
+ Simulating and analyzing large chemical reaction mechanisms
+
+ + Ideal gas and dilute liquid phases.
+ + Constant T and P and constant V adiabatic ideal gas reactors.
+ + Constant T and V dilute liquid reactors.
+ + Diffusion limited rates. Sensitivity analysis for all reactors.
+ + Flux diagrams with molecular images (if molecular information is provided).
+ - [NumCME.jl](https://github.com/voduchuy/NumCME.jl): High-performance simulation of chemical master equations (CME)
+
+ + Transient solution of the CME
+ + Dynamic state spaces
+ + Accepts reaction systems defined using Catalyst.jl DSL.
+ - [FiniteStateProjection.jl](https://github.com/kaandocal/FiniteStateProjection.jl): High-performance simulation of
+ chemical master equations (CME) via finite state projections
+
+ + Accepts reaction systems defined using Catalyst.jl DSL.
## Compatible Numerical Solvers
@@ -140,80 +151,104 @@ an `ODEProblem` to then be solved by a numerical ODE solver. Below is a list of
the solver libraries which are the numerical targets of the ModelingToolkit
system:
-- [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/)
- - Multi-package interface of high performance numerical solvers for `ODESystem`,
- `SDESystem`, and `JumpSystem`
-- [NonlinearSolve.jl](https://docs.sciml.ai/NonlinearSolve/stable/)
- - High performance numerical solving of `NonlinearSystem`
-- [Optimization.jl](https://docs.sciml.ai/Optimization/stable/)
- - Multi-package interface for numerical solving `OptimizationSystem`
-- [NeuralPDE.jl](https://docs.sciml.ai/NeuralPDE/stable/)
- - Physics-Informed Neural Network (PINN) training on `PDESystem`
-- [MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/)
- - Automated finite difference method (FDM) discretization of `PDESystem`
+ - [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/)
+
+ + Multi-package interface of high performance numerical solvers for `ODESystem`,
+ `SDESystem`, and `JumpSystem`
+
+ - [NonlinearSolve.jl](https://docs.sciml.ai/NonlinearSolve/stable/)
+
+ + High performance numerical solving of `NonlinearSystem`
+ - [Optimization.jl](https://docs.sciml.ai/Optimization/stable/)
+
+ + Multi-package interface for numerical solving `OptimizationSystem`
+ - [NeuralPDE.jl](https://docs.sciml.ai/NeuralPDE/stable/)
+
+ + Physics-Informed Neural Network (PINN) training on `PDESystem`
+ - [MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/)
+
+ + Automated finite difference method (FDM) discretization of `PDESystem`
## Contributing
-- Please refer to the
- [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md)
- for guidance on PRs, issues, and other matters relating to contributing to ModelingToolkit.
-- There are a few community forums:
- - The #diffeq-bridged channel in the [Julia Slack](https://julialang.org/slack/)
- - [JuliaDiffEq](https://gitter.im/JuliaDiffEq/Lobby) on Gitter
- - On the Julia Discourse forums (look for the [modelingtoolkit tag](https://discourse.julialang.org/tag/modelingtoolkit))
- - See also [SciML Community page](https://sciml.ai/community/)
+ - Please refer to the
+ [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md)
+ for guidance on PRs, issues, and other matters relating to contributing to ModelingToolkit.
+
+ - There are a few community forums:
+
+ + The #diffeq-bridged channel in the [Julia Slack](https://julialang.org/slack/)
+ + [JuliaDiffEq](https://gitter.im/JuliaDiffEq/Lobby) on Gitter
+ + On the Julia Discourse forums (look for the [modelingtoolkit tag](https://discourse.julialang.org/tag/modelingtoolkit))
+ + See also [SciML Community page](https://sciml.ai/community/)
## Reproducibility
+
```@raw html
The documentation of this SciML package was built using these direct dependencies,
```
+
```@example
using Pkg # hide
Pkg.status() # hide
```
+
```@raw html
and using this machine and Julia version.
```
+
```@example
using InteractiveUtils # hide
versioninfo() # hide
```
+
```@raw html
A more complete overview of all dependencies and their versions is also provided.
```
+
```@example
using Pkg # hide
-Pkg.status(;mode = PKGMODE_MANIFEST) # hide
+Pkg.status(; mode = PKGMODE_MANIFEST) # hide
```
+
```@raw html