-
-
Notifications
You must be signed in to change notification settings - Fork 213
/
Copy pathindex.html
242 lines (216 loc) · 43.7 KB
/
index.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Event Handling and Callback Functions · ModelingToolkit.jl</title><meta name="title" content="Event Handling and Callback Functions · ModelingToolkit.jl"/><meta property="og:title" content="Event Handling and Callback Functions · ModelingToolkit.jl"/><meta property="twitter:title" content="Event Handling and Callback Functions · ModelingToolkit.jl"/><meta name="description" content="Documentation for ModelingToolkit.jl."/><meta property="og:description" content="Documentation for ModelingToolkit.jl."/><meta property="twitter:description" content="Documentation for ModelingToolkit.jl."/><meta property="og:url" content="https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/"/><meta property="twitter:url" content="https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/"/><link rel="canonical" href="https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/"/><script data-outdated-warner src="../../assets/warner.js"></script><link href="https://cdnjs.cloudflare.com/ajax/libs/lato-font/3.0.0/css/lato-font.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/juliamono/0.050/juliamono.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.4.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.8/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../search_index.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/catppuccin-mocha.css" data-theme-name="catppuccin-mocha"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/catppuccin-macchiato.css" data-theme-name="catppuccin-macchiato"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/catppuccin-frappe.css" data-theme-name="catppuccin-frappe"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/catppuccin-latte.css" data-theme-name="catppuccin-latte"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-dark.css" data-theme-name="documenter-dark" data-theme-primary-dark/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="../../assets/themeswap.js"></script><link href="../../assets/favicon.ico" rel="icon" type="image/x-icon"/></head><body><div id="documenter"><nav class="docs-sidebar"><a class="docs-logo" href="../../"><img src="../../assets/logo.png" alt="ModelingToolkit.jl logo"/></a><div class="docs-package-name"><span class="docs-autofit"><a href="../../">ModelingToolkit.jl</a></span></div><button class="docs-search-query input is-rounded is-small is-clickable my-2 mx-auto py-1 px-2" id="documenter-search-query">Search docs (Ctrl + /)</button><ul class="docs-menu"><li><a class="tocitem" href="../../">Home</a></li><li><a class="tocitem" href="../../tutorials/ode_modeling/">Getting Started with ModelingToolkit.jl</a></li><li><span class="tocitem">Tutorials</span><ul><li><a class="tocitem" href="../../tutorials/acausal_components/">Acausal Component-Based Modeling</a></li><li><a class="tocitem" href="../../tutorials/nonlinear/">Modeling Nonlinear Systems</a></li><li><a class="tocitem" href="../../tutorials/initialization/">Initialization of ODESystems</a></li><li><a class="tocitem" href="../../tutorials/optimization/">Modeling Optimization Problems</a></li><li><a class="tocitem" href="../../tutorials/modelingtoolkitize/">Modelingtoolkitize: Automatically Translating Numerical to Symbolic Code</a></li><li><a class="tocitem" href="../../tutorials/programmatically_generating/">Programmatically Generating and Scripting ODESystems</a></li><li><a class="tocitem" href="../../tutorials/stochastic_diffeq/">Modeling with Stochasticity</a></li><li><a class="tocitem" href="../../tutorials/discrete_system/">(Experimental) Modeling Discrete Systems</a></li><li><a class="tocitem" href="../../tutorials/parameter_identifiability/">Parameter Identifiability in ODE Models</a></li><li><a class="tocitem" href="../../tutorials/change_independent_variable/">Changing the independent variable of ODEs</a></li><li><a class="tocitem" href="../../tutorials/bifurcation_diagram_computation/">Bifurcation Diagrams</a></li><li><a class="tocitem" href="../../tutorials/SampledData/">Clocks and Sampled-Data Systems</a></li><li><a class="tocitem" href="../../tutorials/domain_connections/">Domains</a></li><li><a class="tocitem" href="../../tutorials/callable_params/">Callable parameters and interpolating data</a></li><li><a class="tocitem" href="../../tutorials/linear_analysis/">Linear Analysis</a></li><li><a class="tocitem" href="../../tutorials/fmi/">Importing FMUs</a></li></ul></li><li><span class="tocitem">Examples</span><ul><li><input class="collapse-toggle" id="menuitem-4-1" type="checkbox"/><label class="tocitem" for="menuitem-4-1"><span class="docs-label">Basic Examples</span><i class="docs-chevron"></i></label><ul class="collapsed"><li><a class="tocitem" href="../../examples/higher_order/">Automatic Transformation of Nth Order ODEs to 1st Order ODEs</a></li><li><a class="tocitem" href="../../examples/spring_mass/">Component-Based Modeling of a Spring-Mass System</a></li><li><a class="tocitem" href="../../examples/modelingtoolkitize_index_reduction/">Automated Index Reduction of DAEs</a></li><li><a class="tocitem" href="../../examples/remake/">Optimizing through an ODE solve and re-creating MTK Problems</a></li></ul></li><li><input class="collapse-toggle" id="menuitem-4-2" type="checkbox"/><label class="tocitem" for="menuitem-4-2"><span class="docs-label">Advanced Examples</span><i class="docs-chevron"></i></label><ul class="collapsed"><li><a class="tocitem" href="../../examples/tearing_parallelism/">Exposing More Parallelism By Tearing Algebraic Equations in ODESystems</a></li><li><a class="tocitem" href="../../examples/sparse_jacobians/">Automated Sparse Analytical Jacobians</a></li><li><a class="tocitem" href="../../examples/perturbation/">Symbolic-Numeric Perturbation Theory for ODEs</a></li></ul></li></ul></li><li><span class="tocitem">Basics</span><ul><li><a class="tocitem" href="../AbstractSystem/">The AbstractSystem Interface</a></li><li><a class="tocitem" href="../ContextualVariables/">Contextual Variable Types</a></li><li><a class="tocitem" href="../Variable_metadata/">Symbolic Metadata</a></li><li><a class="tocitem" href="../Composition/">Composing Models and Building Reusable Components</a></li><li class="is-active"><a class="tocitem" href>Event Handling and Callback Functions</a><ul class="internal"><li><a class="tocitem" href="#Continuous-Events"><span>Continuous Events</span></a></li><li><a class="tocitem" href="#Discrete-events-support"><span>Discrete events support</span></a></li><li><a class="tocitem" href="#Saving-discrete-values"><span>Saving discrete values</span></a></li><li><a class="tocitem" href="#imp_affects"><span>(Experimental) Imperative affects</span></a></li></ul></li><li><a class="tocitem" href="../Linearization/">Linearization</a></li><li><a class="tocitem" href="../InputOutput/">Input output</a></li><li><a class="tocitem" href="../MTKLanguage/">ModelingToolkit Language: Modeling with <code>@mtkmodel</code>, <code>@connectors</code> and <code>@mtkbuild</code></a></li><li><a class="tocitem" href="../Validation/">Model Validation and Units</a></li><li><a class="tocitem" href="../Debugging/">Debugging</a></li><li><a class="tocitem" href="../DependencyGraphs/">Dependency Graphs</a></li><li><a class="tocitem" href="../Precompilation/">Working with Precompilation and Binary Building</a></li><li><a class="tocitem" href="../FAQ/">Frequently Asked Questions</a></li></ul></li><li><span class="tocitem">System Types</span><ul><li><a class="tocitem" href="../../systems/ODESystem/">ODESystem</a></li><li><a class="tocitem" href="../../systems/SDESystem/">SDESystem</a></li><li><a class="tocitem" href="../../systems/JumpSystem/">JumpSystem</a></li><li><a class="tocitem" href="../../systems/NonlinearSystem/">NonlinearSystem</a></li><li><a class="tocitem" href="../../systems/OptimizationSystem/">OptimizationSystem</a></li><li><a class="tocitem" href="../../systems/PDESystem/">PDESystem</a></li><li><a class="tocitem" href="../../systems/DiscreteSystem/">DiscreteSystem</a></li><li><a class="tocitem" href="../../systems/ImplicitDiscreteSystem/">ImplicitDiscreteSystem</a></li></ul></li><li><a class="tocitem" href="../../comparison/">Comparison of ModelingToolkit vs Equation-Based and Block Modeling Languages</a></li><li><a class="tocitem" href="../../internals/">Internal Details</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><a class="docs-sidebar-button docs-navbar-link fa-solid fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a><nav class="breadcrumb"><ul class="is-hidden-mobile"><li><a class="is-disabled">Basics</a></li><li class="is-active"><a href>Event Handling and Callback Functions</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>Event Handling and Callback Functions</a></li></ul></nav><div class="docs-right"><a class="docs-navbar-link" href="https://github.com/SciML/ModelingToolkit.jl" title="View the repository on GitHub"><span class="docs-icon fa-brands"></span><span class="docs-label is-hidden-touch">GitHub</span></a><a class="docs-navbar-link" href="https://github.com/SciML/ModelingToolkit.jl/blob/master/docs/src/basics/Events.md" title="Edit source on GitHub"><span class="docs-icon fa-solid"></span></a><a class="docs-settings-button docs-navbar-link fa-solid fa-gear" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-article-toggle-button fa-solid fa-chevron-up" id="documenter-article-toggle-button" href="javascript:;" title="Collapse all docstrings"></a></div></header><article class="content" id="documenter-page"><h1 id="events"><a class="docs-heading-anchor" href="#events">Event Handling and Callback Functions</a><a id="events-1"></a><a class="docs-heading-anchor-permalink" href="#events" title="Permalink"></a></h1><p>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 into DifferentialEquations.jl <a href="https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/"><code>ContinuousCallback</code>s or <code>DiscreteCallback</code>s</a>, or into more specialized callback types from the <a href="https://docs.sciml.ai/DiffEqDocs/stable/features/callback_library/">DiffEqCallbacks.jl</a> library.</p><p><a href="../../systems/ODESystem/#ODESystem"><code>ODESystem</code></a>s and <a href="../../systems/SDESystem/#SDESystem"><code>SDESystem</code></a>s accept keyword arguments <code>continuous_events</code> and <code>discrete_events</code> to symbolically encode continuous or discrete callbacks. <a href="../../systems/JumpSystem/#JumpSystem"><code>JumpSystem</code></a>s currently support only <code>discrete_events</code>. Continuous events are applied when a given condition becomes zero, with root finding used to determine the time at which a zero crossing occurred. Discrete events are applied when a condition tested after each timestep evaluates to true. See the <a href="https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/">DifferentialEquations docs</a> for more detail.</p><p>Events involve both a <em>condition</em> function (for the zero crossing or truth test), and an <em>affect</em> function (for determining how to update the system when the event occurs). These can both be specified symbolically, but a more <a href="#func_affects">general functional affect</a> representation is also allowed, as described below.</p><h2 id="Continuous-Events"><a class="docs-heading-anchor" href="#Continuous-Events">Continuous Events</a><a id="Continuous-Events-1"></a><a class="docs-heading-anchor-permalink" href="#Continuous-Events" title="Permalink"></a></h2><p>The basic purely symbolic continuous event interface to encode <em>one</em> continuous event is</p><pre><code class="language-julia hljs">AbstractSystem(eqs, ...; continuous_events::Vector{Equation})
AbstractSystem(eqs, ...; continuous_events::Pair{Vector{Equation}, Vector{Equation}})</code></pre><p>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 <code>Pair</code> 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</p><pre><code class="language-julia hljs">single_unknown_variable ~ expression_involving_any_variables_or_parameters</code></pre><p>or</p><pre><code class="language-julia hljs">single_parameter ~ expression_involving_any_variables_or_parameters</code></pre><p>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 <a href="#func_affects">general functional affects</a>.</p><p>Finally, multiple events can be encoded via a <code>Vector{Pair{Vector{Equation}, Vector{Equation}}}</code>.</p><h3 id="Example:-Friction"><a class="docs-heading-anchor" href="#Example:-Friction">Example: Friction</a><a id="Example:-Friction-1"></a><a class="docs-heading-anchor-permalink" href="#Example:-Friction" title="Permalink"></a></h3><p>The system below illustrates how continuous events can be used to model Coulomb friction</p><pre><code class="language-julia hljs">using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D
function UnitMassWithFriction(k; name)
@variables x(t)=0 v(t)=0
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
@mtkbuild m = UnitMassWithFriction(0.7)
prob = ODEProblem(m, Pair[], (0, 10pi))
sol = solve(prob, Tsit5())
plot(sol)</code></pre><img src="edbe5a42.svg" alt="Example block output"/><h3 id="Example:-Bouncing-ball"><a class="docs-heading-anchor" href="#Example:-Bouncing-ball">Example: Bouncing ball</a><a id="Example:-Bouncing-ball-1"></a><a class="docs-heading-anchor-permalink" href="#Example:-Bouncing-ball" title="Permalink"></a></h3><p>In the documentation for <a href="https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#Example-1:-Bouncing-Ball">DifferentialEquations</a>, we have an example where a bouncing ball is simulated using callbacks which have an <code>affect!</code> on the state. We can model the same system using ModelingToolkit like this</p><pre><code class="language-julia hljs">@variables x(t)=1 v(t)=0
root_eqs = [x ~ 0] # the event happens at the ground x(t) = 0
affect = [v ~ -v] # the effect is that the velocity changes sign
@mtkbuild ball = ODESystem([D(x) ~ v
D(v) ~ -9.8], t; continuous_events = root_eqs => affect) # equation => affect
tspan = (0.0, 5.0)
prob = ODEProblem(ball, Pair[], tspan)
sol = solve(prob, Tsit5())
@assert 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
plot(sol)</code></pre><img src="2cf18273.svg" alt="Example block output"/><h3 id="Test-bouncing-ball-in-2D-with-walls"><a class="docs-heading-anchor" href="#Test-bouncing-ball-in-2D-with-walls">Test bouncing ball in 2D with walls</a><a id="Test-bouncing-ball-in-2D-with-walls-1"></a><a class="docs-heading-anchor-permalink" href="#Test-bouncing-ball-in-2D-with-walls" title="Permalink"></a></h3><p>Multiple events? No problem! This example models a bouncing ball in 2D that is enclosed by two walls at <span>$y = \pm 1.5$</span>.</p><pre><code class="language-julia hljs">@variables x(t)=1 y(t)=0 vx(t)=0 vy(t)=2
continuous_events = [[x ~ 0] => [vx ~ -vx]
[y ~ -1.5, y ~ 1.5] => [vy ~ -vy]]
@mtkbuild 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)
tspan = (0.0, 10.0)
prob = ODEProblem(ball, Pair[], tspan)
sol = solve(prob, Tsit5())
@assert 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
@assert minimum(sol[y]) > -1.5 # check wall conditions
@assert maximum(sol[y]) < 1.5 # check wall conditions
tv = sort([LinRange(0, 10, 200); sol.t])
plot(sol(tv)[y], sol(tv)[x], line_z = tv)
vline!([-1.5, 1.5], l = (:black, 5), primary = false)
hline!([0], l = (:black, 5), primary = false)</code></pre><img src="fb29b318.svg" alt="Example block output"/><h3 id="func_affects"><a class="docs-heading-anchor" href="#func_affects">Generalized functional affect support</a><a id="func_affects-1"></a><a class="docs-heading-anchor-permalink" href="#func_affects" title="Permalink"></a></h3><p>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 <code>tuple</code>:</p><pre><code class="language-julia hljs">[x ~ 0] => (affect!, [v, x], [p, q], [discretes...], ctx)</code></pre><p>where, <code>affect!</code> is a Julia function with the signature: <code>affect!(integ, u, p, ctx)</code>; <code>[u,v]</code> and <code>[p,q]</code> are the symbolic unknowns (variables) and parameters that are accessed by <code>affect!</code>, respectively; <code>discretes</code> are the parameters modified by <code>affect!</code>, if any; and <code>ctx</code> is any context that is passed to <code>affect!</code> as the <code>ctx</code> argument.</p><p><code>affect!</code> receives a <a href="https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/">DifferentialEquations.jl integrator</a> as its first argument, which can then be used to access unknowns and parameters that are provided in the <code>u</code> and <code>p</code> arguments (implemented as <code>NamedTuple</code>s). The integrator can also be manipulated more generally to control solution behavior, see the <a href="https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/">integrator interface</a> documentation. In affect functions, we have that</p><pre><code class="language-julia hljs">function affect!(integ, u, p, ctx)
# integ.t is the current time
# integ.u[u.v] is the value of the unknown `v` above
# integ.ps[p.q] is the value of the parameter `q` above
end</code></pre><p>When accessing variables of a sub-system, it can be useful to rename them (alternatively, an affect function may be reused in different contexts):</p><pre><code class="language-julia hljs">[x ~ 0] => (affect!, [resistor₊v => :v, x], [p, q => :p2], [], ctx)</code></pre><p>Here, the symbolic variable <code>resistor₊v</code> is passed as <code>v</code> while the symbolic parameter <code>q</code> has been renamed <code>p2</code>.</p><p>As an example, here is the bouncing ball example from above using the functional affect interface:</p><pre><code class="language-julia hljs">sts = @variables x(t), v(t)
par = @parameters g = 9.8
bb_eqs = [D(x) ~ v
D(v) ~ -g]
function bb_affect!(integ, u, p, ctx)
integ.u[u.v] = -integ.u[u.v]
end
reflect = [x ~ 0] => (bb_affect!, [v], [], [], nothing)
@mtkbuild bb_sys = ODESystem(bb_eqs, t, sts, par,
continuous_events = reflect)
u0 = [v => 0.0, x => 1.0]
bb_prob = ODEProblem(bb_sys, u0, (0, 5.0))
bb_sol = solve(bb_prob, Tsit5())
plot(bb_sol)</code></pre><img src="6b176ceb.svg" alt="Example block output"/><h2 id="Discrete-events-support"><a class="docs-heading-anchor" href="#Discrete-events-support">Discrete events support</a><a id="Discrete-events-support-1"></a><a class="docs-heading-anchor-permalink" href="#Discrete-events-support" title="Permalink"></a></h2><p>In addition to continuous events, discrete events are also supported. The general interface to represent a collection of discrete events is</p><pre><code class="language-julia hljs">AbstractSystem(eqs, ...; discrete_events = [condition1 => affect1, condition2 => affect2])</code></pre><p>where conditions are symbolic expressions that should evaluate to <code>true</code> when an individual affect should be executed. Here <code>affect1</code> and <code>affect2</code> are each either a vector of one or more symbolic equations, or a functional affect, just as for continuous events. As before, for any <em>one</em> event the symbolic affect equations can either all change unknowns (i.e. variables) or all change parameters, but one cannot currently mix unknown variable and parameter changes within one individual event.</p><h3 id="Example:-Injecting-cells-into-a-population"><a class="docs-heading-anchor" href="#Example:-Injecting-cells-into-a-population">Example: Injecting cells into a population</a><a id="Example:-Injecting-cells-into-a-population-1"></a><a class="docs-heading-anchor-permalink" href="#Example:-Injecting-cells-into-a-population" title="Permalink"></a></h3><p>Suppose we have a population of <code>N(t)</code> cells that can grow and die, and at time <code>t1</code> we want to inject <code>M</code> more cells into the population. We can model this by</p><pre><code class="language-julia hljs">@parameters M tinject α
@variables N(t)
Dₜ = Differential(t)
eqs = [Dₜ(N) ~ α - N]
# at time tinject we inject M cells
injection = (t == tinject) => [N ~ N + M]
u0 = [N => 0.0]
tspan = (0.0, 20.0)
p = [α => 100.0, tinject => 10.0, M => 50]
@mtkbuild osys = ODESystem(eqs, t, [N], [α, M, tinject]; discrete_events = injection)
oprob = ODEProblem(osys, u0, tspan, p)
sol = solve(oprob, Tsit5(); tstops = 10.0)
plot(sol)</code></pre><img src="136dd784.svg" alt="Example block output"/><p>Notice, with generic discrete events that we want to occur at one or more fixed times, we need to also set the <code>tstops</code> keyword argument to <code>solve</code> to ensure the integrator stops at that time. In the next section, we show how one can avoid this by using a preset-time callback.</p><p>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</p><pre><code class="language-julia hljs">injection = ((t == tinject) & (N < 50)) => [N ~ N + M]
@mtkbuild osys = ODESystem(eqs, t, [N], [M, tinject, α]; discrete_events = injection)
oprob = ODEProblem(osys, u0, tspan, p)
sol = solve(oprob, Tsit5(); tstops = 10.0)
plot(sol)</code></pre><img src="11456fdf.svg" alt="Example block output"/><p>Since the solution is <em>not</em> smaller than half its steady-state value at the event time, the event condition now returns false. Here we used logical and, <code>&</code>, instead of the short-circuiting logical and, <code>&&</code>, as currently the latter cannot be used within symbolic expressions.</p><p>Let's now also add a drug at time <code>tkill</code> that turns off production of new cells, modeled by setting <code>α = 0.0</code></p><pre><code class="language-julia hljs">@parameters tkill
# we reset the first event to just occur at tinject
injection = (t == tinject) => [N ~ N + M]
# at time tkill we turn off production of cells
killing = (t == tkill) => [α ~ 0.0]
tspan = (0.0, 30.0)
p = [α => 100.0, tinject => 10.0, M => 50, tkill => 20.0]
@mtkbuild osys = ODESystem(eqs, t, [N], [α, M, tinject, tkill];
discrete_events = [injection, killing])
oprob = ODEProblem(osys, u0, tspan, p)
sol = solve(oprob, Tsit5(); tstops = [10.0, 20.0])
plot(sol)</code></pre><img src="ca1c5d37.svg" alt="Example block output"/><h3 id="Periodic-and-preset-time-events"><a class="docs-heading-anchor" href="#Periodic-and-preset-time-events">Periodic and preset-time events</a><a id="Periodic-and-preset-time-events-1"></a><a class="docs-heading-anchor-permalink" href="#Periodic-and-preset-time-events" title="Permalink"></a></h3><p>Two important subclasses of discrete events are periodic and preset-time events.</p><p>A preset-time event is triggered at specific set times, which can be passed in a vector like</p><pre><code class="language-julia hljs">discrete_events = [[1.0, 4.0] => [v ~ -v]]</code></pre><p>This will change the sign of <code>v</code> <em>only</em> at <code>t = 1.0</code> and <code>t = 4.0</code>.</p><p>As such, our last example with treatment and killing could instead be modeled by</p><pre><code class="language-julia hljs">injection = [10.0] => [N ~ N + M]
killing = [20.0] => [α ~ 0.0]
p = [α => 100.0, M => 50]
@mtkbuild osys = ODESystem(eqs, t, [N], [α, M];
discrete_events = [injection, killing])
oprob = ODEProblem(osys, u0, tspan, p)
sol = solve(oprob, Tsit5())
plot(sol)</code></pre><img src="e2ebd22d.svg" alt="Example block output"/><p>Notice, one advantage of using a preset-time event is that one does not need to also specify <code>tstops</code> in the call to solve.</p><p>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,</p><pre><code class="language-julia hljs">discrete_events = [1.0 => [v ~ -v]]</code></pre><p>will change the sign of <code>v</code> at <code>t = 1.0</code>, <code>2.0</code>, ...</p><p>Finally, we note that to specify an event at precisely one time, say 2.0 below, one must still use a vector</p><pre><code class="language-julia hljs">discrete_events = [[2.0] => [v ~ -v]]</code></pre><h2 id="Saving-discrete-values"><a class="docs-heading-anchor" href="#Saving-discrete-values">Saving discrete values</a><a id="Saving-discrete-values-1"></a><a class="docs-heading-anchor-permalink" href="#Saving-discrete-values" title="Permalink"></a></h2><p>Time-dependent parameters which are updated in callbacks are termed as discrete variables. ModelingToolkit enables automatically saving the timeseries of these discrete variables, and indexing the solution object to obtain the saved timeseries. Consider the following example:</p><pre><code class="language-julia hljs">@variables x(t)
@parameters c(t)
@mtkbuild sys = ODESystem(
D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ c + 1]])
prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0])
sol = solve(prob, Tsit5())
sol[c]</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">7-element Vector{Float64}:
1.0
2.0
3.0
4.0
5.0
6.0
7.0</code></pre><p>The solution object can also be interpolated with the discrete variables</p><pre><code class="language-julia hljs">sol([1.0, 2.0], idxs = [c, c * cos(x)])</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">t: 2-element Vector{Float64}:
1.0
2.0
u: 2-element Vector{Vector{Float64}}:
[2.0, 1.296108550512409]
[3.0, 0.29799145226722035]</code></pre><p>Note that only time-dependent parameters will be saved. If we repeat the above example with this change:</p><pre><code class="language-julia hljs">@variables x(t)
@parameters c
@mtkbuild sys = ODESystem(
D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ c + 1]])
prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0])
sol = solve(prob, Tsit5())
sol.ps[c] # sol[c] will error, since `c` is not a timeseries value</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">7.0</code></pre><p>It can be seen that the timeseries for <code>c</code> is not saved.</p><h2 id="imp_affects"><a class="docs-heading-anchor" href="#imp_affects">(Experimental) Imperative affects</a><a id="imp_affects-1"></a><a class="docs-heading-anchor-permalink" href="#imp_affects" title="Permalink"></a></h2><p>The <code>ImperativeAffect</code> can be used as an alternative to the aforementioned functional affect form. Note that <code>ImperativeAffect</code> is still experimental; to emphasize this, we do not export it and it should be included as <code>ModelingToolkit.ImperativeAffect</code>. <code>ImperativeAffect</code> aims to simplify the manipulation of system state.</p><p>We will use two examples to describe <code>ImperativeAffect</code>: a simple heater and a quadrature encoder. These examples will also demonstrate advanced usage of <code>ModelingToolkit.SymbolicContinuousCallback</code>, the low-level interface of the tuple form converts into that allows control over the SciMLBase-level event that is generated for a continuous event.</p><h3 id="heater_events"><a class="docs-heading-anchor" href="#heater_events">Heater</a><a id="heater_events-1"></a><a class="docs-heading-anchor-permalink" href="#heater_events" title="Permalink"></a></h3><p>Bang-bang control of a heater connected to a leaky plant requires hysteresis in order to prevent rapid control oscillation.</p><pre><code class="language-julia hljs">@variables temp(t)
params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on(t)::Bool=false
eqs = [
D(temp) ~ furnace_on * furnace_power - temp^2 * leakage
]</code></pre><p class="math-container">\[ \begin{align}
\frac{\mathrm{d} \mathtt{temp}\left( t \right)}{\mathrm{d}t} &= \mathtt{furnace\_power} \mathtt{furnace\_on}\left( t \right) - \left( \mathtt{temp}\left( t \right) \right)^{2} \mathtt{leakage}
\end{align}
\]</p><p>Our plant is simple. We have a heater that's turned on and off by the time-indexed parameter <code>furnace_on</code> which adds <code>furnace_power</code> forcing to the system when enabled. We then leak heat proportional to <code>leakage</code> as a function of the square of the current temperature.</p><p>We need a controller with hysteresis to control the plant. We wish the furnace to turn on when the temperature is below <code>furnace_on_threshold</code> and off when above <code>furnace_off_threshold</code>, while maintaining its current state in between. To do this, we create two continuous callbacks:</p><pre><code class="language-julia hljs">using Setfield
furnace_disable = ModelingToolkit.SymbolicContinuousCallback(
[temp ~ furnace_off_threshold],
ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i
@set! x.furnace_on = false
end)
furnace_enable = ModelingToolkit.SymbolicContinuousCallback(
[temp ~ furnace_on_threshold],
ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i
@set! x.furnace_on = true
end)</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">SymbolicContinuousCallback:
Equations:
1-element Vector{Equation}:
temp(t) ~ furnace_on_threshold
Affect:
ImperativeAffect(observed: [], modified: [furnace_on(t) => furnace_on], affect:#4)
Negative-edge affect:
ImperativeAffect(observed: [], modified: [furnace_on(t) => furnace_on], affect:#4)
</code></pre><p>We're using the explicit form of <code>SymbolicContinuousCallback</code> here, though so far we aren't using anything that's not possible with the implicit interface. You can also write</p><pre><code class="language-julia hljs">[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (;
furnace_on)) do x, o, i, c
@set! x.furnace_on = false
end</code></pre><p>and it would work the same.</p><p>The <code>ImperativeAffect</code> is the larger change in this example. <code>ImperativeAffect</code> has the constructor signature</p><pre><code class="language-julia hljs">ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx)</code></pre><p>that accepts the function to call, a named tuple of both the names of and symbolic values representing values in the system to be modified, a named tuple of the values that are merely observed (that is, used from the system but not modified), and a context that's passed to the affect function.</p><p>In our example, each event merely changes whether the furnace is on or off. Accordingly, we pass a <code>modified</code> tuple <code>(; furnace_on)</code> (creating a <code>NamedTuple</code> equivalent to <code>(furnace_on = furnace_on)</code>). <code>ImperativeAffect</code> will then evaluate this before calling our function to fill out all of the numerical values, then apply them back to the system once our affect function returns. Furthermore, it will check that it is possible to do this assignment.</p><p>The function given to <code>ImperativeAffect</code> needs to have the signature:</p><pre><code class="language-julia hljs">f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple</code></pre><p>The function <code>f</code> will be called with <code>observed</code> and <code>modified</code> <code>NamedTuple</code>s that are derived from their respective <code>NamedTuple</code> definitions. In our example, if <code>furnace_on</code> is <code>false</code>, then the value of the <code>x</code> that's passed in as <code>modified</code> will be <code>(furnace_on = false)</code>. The modified values should be passed out in the same format: to set <code>furnace_on</code> to <code>true</code> we need to return a tuple <code>(furnace_on = true)</code>. The examples does this with Setfield, recreating the result tuple before returning it; the returned tuple may optionally be missing values as well, in which case those values will not be written back to the problem.</p><p>Accordingly, we can now interpret the <code>ImperativeAffect</code> definitions to mean that when <code>temp = furnace_off_threshold</code> we will write <code>furnace_on = false</code> back to the system, and when <code>temp = furnace_on_threshold</code> we will write <code>furnace_on = true</code> back to the system.</p><pre><code class="language-julia hljs">@named sys = ODESystem(
eqs, t, [temp], params; continuous_events = [furnace_disable, furnace_enable])
ss = structural_simplify(sys)
prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 10.0))
sol = solve(prob, Tsit5())
plot(sol)
hline!([sol.ps[furnace_off_threshold], sol.ps[furnace_on_threshold]],
l = (:black, 1), primary = false)</code></pre><img src="c1faadb3.svg" alt="Example block output"/><p>Here we see exactly the desired hysteresis. The heater starts on until the temperature hits <code>furnace_off_threshold</code>. The temperature then bleeds away until <code>furnace_on_threshold</code> at which point the furnace turns on again until <code>furnace_off_threshold</code> and so on and so forth. The controller is effectively regulating the temperature of the plant.</p><h3 id="quadrature"><a class="docs-heading-anchor" href="#quadrature">Quadrature Encoder</a><a id="quadrature-1"></a><a class="docs-heading-anchor-permalink" href="#quadrature" title="Permalink"></a></h3><p>For a more complex application we'll look at modeling a quadrature encoder attached to a shaft spinning at a constant speed. Traditionally, a quadrature encoder is built out of a code wheel that interrupts the sensors at constant intervals and two sensors slightly out of phase with one another. A state machine can take the pattern of pulses produced by the two sensors and determine the number of steps that the shaft has spun. The state machine takes the new value from each sensor and the old values and decodes them into the direction that the wheel has spun in this step.</p><pre><code class="language-julia hljs">@variables theta(t) omega(t)
params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0
eqs = [D(theta) ~ omega
omega ~ 1.0]</code></pre><p class="math-container">\[ \begin{align}
\frac{\mathrm{d} \mathtt{theta}\left( t \right)}{\mathrm{d}t} &= \mathtt{omega}\left( t \right) \\
\mathtt{omega}\left( t \right) &= 1
\end{align}
\]</p><p>Our continuous-time system is extremely simple. We have two unknown variables <code>theta</code> for the angle of the shaft and <code>omega</code> for the rate at which it's spinning. We then have parameters for the state machine <code>qA, qB, hA, hB</code> (corresponding to the current quadrature of the A/B sensors and the historical ones) and a step count <code>cnt</code>.</p><p>We'll then implement the decoder as a simple Julia function.</p><pre><code class="language-julia hljs">function decoder(oldA, oldB, newA, newB)
state = (oldA, oldB, newA, newB)
if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) ||
state == (0, 1, 0, 0)
return 1
elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) ||
state == (1, 0, 0, 0)
return -1
elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) ||
state == (1, 1, 1, 1)
return 0
else
return 0 # err is interpreted as no movement
end
end</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">decoder (generic function with 1 method)</code></pre><p>Based on the current and old state, this function will return 1 if the wheel spun in the positive direction, -1 if in the negative, and 0 otherwise.</p><p>The encoder state advances when the occlusion begins or ends. We model the code wheel as simply detecting when <code>cos(100*theta)</code> is 0; if we're at a positive edge of the 0 crossing, then we interpret that as occlusion (so the discrete <code>qA</code> goes to 1). Otherwise, if <code>cos</code> is going negative, we interpret that as lack of occlusion (so the discrete goes to 0). The decoder function is then invoked to update the count with this new information.</p><p>We can implement this in one of two ways: using edge sign detection or right root finding. For exposition, we will implement each sensor differently.</p><p>For sensor A, we're using the edge detection method. By providing a different affect to <code>SymbolicContinuousCallback</code>'s <code>affect_neg</code> argument, we can specify different behaviour for the negative crossing vs. the positive crossing of the root. In our encoder, we interpret this as occlusion or nonocclusion of the sensor, update the internal state, and tick the decoder.</p><pre><code class="language-julia hljs">qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0],
ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i
@set! x.hA = x.qA
@set! x.hB = o.qB
@set! x.qA = 1
@set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
x
end,
affect_neg = ModelingToolkit.ImperativeAffect(
(; qA, hA, hB, cnt), (; qB)) do x, o, c, i
@set! x.hA = x.qA
@set! x.hB = o.qB
@set! x.qA = 0
@set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
x
end)</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">SymbolicContinuousCallback:
Equations:
1-element Vector{Equation}:
cos(100theta(t)) ~ 0
Affect:
ImperativeAffect(observed: [qB => qB], modified: [qA => qA, hA => hA, hB => hB, cnt => cnt], affect:#6)
Negative-edge affect:
ImperativeAffect(observed: [qB => qB], modified: [qA => qA, hA => hA, hB => hB, cnt => cnt], affect:#7)
</code></pre><p>The other way we can implement a sensor is by changing the root find. Normally, we use left root finding; the affect will be invoked instantaneously <em>before</em> the root is crossed. This makes it trickier to figure out what the new state is. Instead, we can use right root finding:</p><pre><code class="language-julia hljs">qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0],
ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, c, i
@set! x.hA = o.qA
@set! x.hB = x.qB
@set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0)
@set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
x
end; rootfind = SciMLBase.RightRootFind)</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">SymbolicContinuousCallback:
Equations:
1-element Vector{Equation}:
cos(-1.5707963267948966 + 100theta(t)) ~ 0
Affect:
ImperativeAffect(observed: [qA => qA, theta(t) => theta], modified: [qB => qB, hA => hA, hB => hB, cnt => cnt], affect:#10)
Negative-edge affect:
ImperativeAffect(observed: [qA => qA, theta(t) => theta], modified: [qB => qB, hA => hA, hB => hB, cnt => cnt], affect:#10)
</code></pre><p>Here, sensor B is located <code>π / 2</code> behind sensor A in angular space, so we're adjusting our trigger function accordingly. We here ask for right root finding on the callback, so we know that the value of said function will have the "new" sign rather than the old one. Thus, we can determine the new state of the sensor from the sign of the indicator function evaluated at the affect activation point, with -1 mapped to 0.</p><p>We can now simulate the encoder.</p><pre><code class="language-julia hljs">@named sys = ODESystem(
eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt])
ss = structural_simplify(sys)
prob = ODEProblem(ss, [theta => 0.0], (0.0, pi))
sol = solve(prob, Tsit5(); dtmax = 0.01)
sol.ps[cnt]</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">198</code></pre><p><code>cos(100*theta)</code> will have 200 crossings in the half rotation we've gone through, so the encoder would notionally count 200 steps. Our encoder counts 198 steps (it loses one step to initialization and one step due to the final state falling squarely on an edge).</p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../Composition/">« Composing Models and Building Reusable Components</a><a class="docs-footer-nextpage" href="../Linearization/">Linearization »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option><option value="catppuccin-latte">catppuccin-latte</option><option value="catppuccin-frappe">catppuccin-frappe</option><option value="catppuccin-macchiato">catppuccin-macchiato</option><option value="catppuccin-mocha">catppuccin-mocha</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.9.0 on <span class="colophon-date" title="Sunday 23 March 2025 13:30">Sunday 23 March 2025</span>. Using Julia version 1.10.9.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>