Composing Models and Building Reusable Components
The symbolic models of ModelingToolkit can be composed together to easily build large models. The composition is lazy and only instantiated at the time of conversion to numerical models, allowing a more performant way in terms of computation time and memory.
Simple Model Composition Example
The following is an example of building a model in a library with an optional forcing function, and allowing the user to specify the forcing later. Here, the library author defines a component named decay
. The user then builds two decay
components and connects them, saying the forcing term of decay1
is a constant while the forcing term of decay2
is the value of the state variable x
.
using ModelingToolkit
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()
@named decay2 = decay()
@parameters t
D = Differential(t)
connected = compose(ODESystem([
decay2.f ~ decay1.x
D(decay1.f) ~ 0
], t; name=:connected), decay1, decay2)
equations(connected)
#4-element Vector{Equation}:
# Differential(t)(decay1₊f(t)) ~ 0
# decay2₊f(t) ~ decay1₊x(t)
# Differential(t)(decay1₊x(t)) ~ decay1₊f(t) - (decay1₊a*(decay1₊x(t)))
# Differential(t)(decay2₊x(t)) ~ decay2₊f(t) - (decay2₊a*(decay2₊x(t)))
simplified_sys = structural_simplify(connected)
equations(simplified_sys)
\begin{align} \frac{ddecay1{+}f(t)}{dt} =& -0.0 \ \frac{ddecay1{+}x(t)}{dt} =& - decay1{+}a \mathrm{decay1{+}x}\left( t \right) + \mathrm{decay1{+}f}\left( t \right) \ \frac{ddecay2{+}x(t)}{dt} =& - decay2{+}a \mathrm{decay2{+}x}\left( t \right) + \mathrm{decay2_{+}f}\left( t \right) \end{align}
Now we can solve the system:
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)
sol = solve(prob, Tsit5())
sol[decay2.f]
23-element Vector{Float64}:
1.0
0.988412780354162
0.9276920413597063
0.83083442676386
0.7238511010482006
0.6064716273042756
0.4906266726783709
0.38123692999606124
0.284066375797688
0.20250582537955814
⋮
0.016682568907159096
0.008438451900843247
0.004154119364674407
0.0019929454742438255
0.0009349151861883882
0.0004281899181999366
0.00019123722531328084
8.29837215148343e-5
4.54276561596786e-5
Basics of Model Composition
Every AbstractSystem
has a system
keyword argument for specifying subsystems. A model is the composition of itself and its subsystems. For example, if we have:
@named sys = compose(ODESystem(eqs,indepvar,states,ps),subsys)
the equations
of sys
is the concatenation of get_eqs(sys)
and equations(subsys)
, the states are the concatenation of their states, etc. When the ODEProblem
or ODEFunction
is generated from this system, it will build and compile the functions associated with this composition.
The new equations within the higher level system can access the variables in the lower level system by namespacing via the nameof(subsys)
. For example, let's say there is a variable x
in states
and a variable x
in subsys
. We can declare that these two variables are the same by specifying their equality: x ~ subsys.x
in the eqs
for sys
. This algebraic relationship can then be simplified by transformations like structural_simplify
which will be described later.
Numerics with Composed Models
These composed models can then be directly transformed into their associated SciMLProblem
type using the standard constructors. When this is done, the initial conditions and parameters must be specified in their namespaced form. For example:
u0 = [
x => 2.0
subsys.x => 2.0
]
Note that any default values within the given subcomponent will be used if no override is provided at construction time. If any values for initial conditions or parameters are unspecified an error will be thrown.
When the model is numerically solved, the solution can be accessed via its symbolic values. For example, if sol
is the ODESolution
, one can use sol[x]
and sol[subsys.x]
to access the respective timeseries in the solution. All other indexing rules stay the same, so sol[x,1:5]
accesses the first through fifth values of x
. Note that this can be done even if the variable x
is eliminated from the system from transformations like alias_elimination
or tearing
: the variable will be lazily reconstructed on demand.
Variable scope and parameter expressions
In some scenarios, it could be useful for model parameters to be expressed in terms of other parameters, or shared between common subsystems. To facilitate this, ModelingToolkit supports symbolic expressions in default values, and scoped variables.
With symbolic parameters, it is possible to set the default value of a parameter or initial condition to an expression of other variables.
# ...
sys = ODESystem(
# ...
# directly in the defauls argument
defaults=Pair{Num, Any}[
x => u,
y => σ,
z => u-0.1,
])
# by assigning to the parameter
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.
@variables a b c d
# a is a local variable
b = ParentScope(b) # b is a variable that belongs to one level up in the hierarchy
c = ParentScope(ParentScope(c)) # ParentScope can be nested
d = GlobalScope(d) # global variables will never be namespaced
Structural Simplify
In many cases, the nicest way to build a model may leave a lot of unnecessary variables. Thus one may want to remove these equations before numerically solving. The structural_simplify
function removes these trivial equality relationships and trivial singularity equations, i.e. equations which result in 0~0
expressions, in over-specified systems.
Inheritance and Combine
Model inheritance can be done in two ways: implicitly or explicitly. First, one can use the extend
function to extend a base model with another set of equations, states, and parameters. An example can be found in the acausal components tutorial.
The explicit way is to shadow variables with equality expressions. For example, let's assume we have three separate systems which we want to compose to a single one. This is how one could explicitly forward all states and parameters to the higher level system:
using ModelingToolkit, OrdinaryDiffEq, Plots
## Library code
@parameters t
D = Differential(t)
@variables S(t), I(t), R(t)
N = S + I + R
@parameters β,γ
@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)
\begin{align} S\left( t \right) =& \mathrm{ieqn{+}S}\left( t \right) \ I\left( t \right) =& \mathrm{seqn{+}I}\left( t \right) \ R\left( t \right) =& \mathrm{ieqn{+}R}\left( t \right) \ \mathrm{ieqn{+}S}\left( t \right) =& \mathrm{seqn{+}S}\left( t \right) \ \mathrm{seqn{+}I}\left( t \right) =& \mathrm{ieqn{+}I}\left( t \right) \ \mathrm{seqn{+}R}\left( t \right) =& \mathrm{reqn{+}R}\left( t \right) \ \mathrm{ieqn{+}R}\left( t \right) =& \mathrm{reqn{+}R}\left( t \right) \ \mathrm{reqn{+}I}\left( t \right) =& \mathrm{ieqn{+}I}\left( t \right) \ \frac{dseqn{+}S(t)}{dt} =& \frac{ - seqn+\beta \mathrm{seqn{+}I}\left( t \right) \mathrm{seqn{+}S}\left( t \right)}{\mathrm{seqn{+}I}\left( t \right) + \mathrm{seqn{+}R}\left( t \right) + \mathrm{seqn{+}S}\left( t \right)} \ \frac{dieqn{+}I(t)}{dt} =& \frac{ieqn+\beta \mathrm{ieqn{+}I}\left( t \right) \mathrm{ieqn{+}S}\left( t \right)}{\mathrm{ieqn{+}I}\left( t \right) + \mathrm{ieqn{+}R}\left( t \right) + \mathrm{ieqn{+}S}\left( t \right)} - ieqn+\gamma \mathrm{ieqn{+}I}\left( t \right) \ \frac{dreqn{+}R(t)}{dt} =& reqn+\gamma \mathrm{reqn{+}I}\left( t \right) \end{align}
Note that the states are forwarded by an equality relationship, while the parameters are forwarded through a relationship in their default values. The user of this model can then solve this model simply by specifying the values at the highest level:
sireqn_simple = structural_simplify(sir)
equations(sireqn_simple)
\begin{align} \frac{dseqn{+}S(t)}{dt} =& \frac{ - seqn+\beta \mathrm{seqn{+}I}\left( t \right) \mathrm{seqn{+}S}\left( t \right)}{\mathrm{seqn{+}I}\left( t \right) + \mathrm{seqn{+}R}\left( t \right) + \mathrm{seqn{+}S}\left( t \right)} \ \frac{dieqn{+}I(t)}{dt} =& \frac{ieqn+\beta \mathrm{ieqn{+}I}\left( t \right) \mathrm{ieqn{+}S}\left( t \right)}{\mathrm{ieqn{+}I}\left( t \right) + \mathrm{ieqn{+}R}\left( t \right) + \mathrm{ieqn{+}S}\left( t \right)} - ieqn+\gamma \mathrm{ieqn{+}I}\left( t \right) \ \frac{dreqn{+}R(t)}{dt} =& reqn+\gamma \mathrm{reqn_{+}I}\left( t \right) \end{align}
## User Code
u0 = [seqn.S => 990.0,
ieqn.I => 10.0,
reqn.R => 0.0]
p = [
β => 0.5
γ => 0.25
]
tspan = (0.0,40.0)
prob = ODEProblem(sireqn_simple,u0,tspan,p,jac=true)
sol = solve(prob,Tsit5())
sol[reqn.R]
18-element Vector{Float64}:
0.0
0.0014142401338982807
0.015567426114815186
0.1581830690301042
1.0627160216919247
3.3906763834079183
7.785038081421038
15.476269434159235
30.102184533124035
58.80475275825802
116.64582495682082
221.0786663303351
368.86570251966197
519.398958278465
655.4762283913897
723.9496381572193
772.9458047004414
775.6835919121963
Tearing Problem Construction
Some system types, specifically ODESystem
and NonlinearSystem
, can be further reduced if structural_simplify
has already been applied to them. This is done by using the alternative problem constructors, ODAEProblem
and BlockNonlinearProblem
respectively. In these cases, the constructor uses the knowledge of the strongly connected components calculated during the process of simplification as the basis for building pre-simplified nonlinear systems in the implicit solving. In summary: these problems are structurally modified, but could be more efficient and more stable.
Components with discontinuous dynamics
When modeling, e.g., impacts, saturations or Coulomb friction, the dynamic equations are discontinuous in either the state or one of its derivatives. This causes the solver to take very small steps around the discontinuity, and sometimes leads to early stopping due to dt <= dt_min
. The correct way to handle such dynamics is to tell the solver about the discontinuity by means of a root-finding equation. ODEsystem
s accept a keyword argument continuous_events
ODESystem(eqs, ...; continuous_events::Vector{Equation})
ODESystem(eqs, ...; continuous_events::Pair{Vector{Equation}, Vector{Equation}})
where equations can be added that evaluate to 0 at discontinuities.
To model events that have an effect on the state, provide events::Pair{Vector{Equation}, Vector{Equation}}
where the first entry in the pair is a vector of equations describing event conditions, and the second vector of equations describe the effect on the state. The effect equations must be of the form
single_state_variable ~ expression_involving_any_variables
Example: Friction
The system below illustrates how this can be used to model Coulomb friction
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
end
@named m = UnitMassWithFriction(0.7)
prob = ODEProblem(m, Pair[], (0, 10pi))
sol = solve(prob, Tsit5())
plot(sol)