(Experimental) Modeling Discrete Systems

In this example, we will use the System API to create an SIR model.

using ModelingToolkit
using ModelingToolkit: t_nounits as t
using OrdinaryDiffEq: solve, FunctionMap

@inline function rate_to_proportion(r, t)
    1 - exp(-r * t)
end
@parameters c δt β γ
@constants h = 1
@variables S(t) I(t) R(t)
k = ShiftIndex(t)
infection = rate_to_proportion(
    β * c * I(k - 1) / (S(k - 1) * h + I(k - 1) + R(k - 1)), δt * h) * S(k - 1)
recovery = rate_to_proportion(γ * h, δt) * I(k - 1)

# Equations
eqs = [S(k) ~ S(k - 1) - infection * h,
    I(k) ~ I(k - 1) + infection - recovery,
    R(k) ~ R(k - 1) + recovery]
@mtkcompile sys = System(eqs, t)

u0 = [S(k - 1) => 990.0, I(k - 1) => 10.0, R(k - 1) => 0.0]
p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1]
tspan = (0.0, 100.0)
prob = DiscreteProblem(sys, vcat(u0, p), tspan)
sol = solve(prob, FunctionMap())

retcode: Success
Interpolation: left-endpoint piecewise constant
t: 101-element Vector{Float64}:
   0.0
   1.0
   2.0
   3.0
   4.0
   5.0
   6.0
   7.0
   8.0
   9.0
   ⋮
  92.0
  93.0
  94.0
  95.0
  96.0
  97.0
  98.0
  99.0
 100.0
u: 101-element Vector{Vector{Float64}}:
 [990.0, 0.0, 10.0]
 [989.5051237293776, 0.24690087971667385, 10.247975390905728]
 [988.9982323978576, 0.4999242936496187, 10.501843308492766]
 [988.479053494723, 0.7592157288009715, 10.761730776476043]
 [987.9473092971556, 1.0249238083995653, 11.027766894444875]
 [987.4027168203793, 1.2972003431544505, 11.300082836466315]
 [986.8449877701768, 1.5762003824739326, 11.57881184734931]
 [986.2738284979789, 1.8620822655923717, 11.86408923642869]
 [985.6889399587333, 2.1550076725435083, 12.156052368723218]
 [985.090017671768, 2.45514167491548, 12.454840653316563]
 ⋮
 [847.1408257083083, 76.95307768711389, 75.90609660457734]
 [843.9317615809414, 78.82720588986679, 77.24103252919139]
 [840.6787392619304, 80.73429377803495, 78.58696696003413]
 [837.3819011012713, 82.67461290570472, 79.9434859930235]
 [834.0414203767291, 84.64842460763423, 81.31015501563623]
 [830.6575018520915, 86.6559794879602, 82.68651865994781]
 [827.2303822969169, 88.69751690774523, 84.07210079533739]
 [823.7603309652128, 90.77326447234499, 85.46640456244167]
 [820.2476500305196, 92.88343751961379, 86.86891244986609]

All shifts must be non-positive, i.e., discrete-time variables may only be indexed at index k, k-1, k-2, .... If default values are provided, they are treated as the value of the variable at the previous timestep. For example, consider the following system to generate the Fibonacci series:

@variables x(t) = 1.0
@mtkcompile sys = System([x ~ x(k - 1) + x(k - 2)], t)

\[ \begin{align} Shift(t, 1)\left( \mathtt{x_{t - 2}}\left( t \right) \right) &= \mathtt{x_{t - 1}}\left( t \right) \\ Shift(t, 1)\left( \mathtt{x_{t - 1}}\left( t \right) \right) &= x\left( t \right) \end{align} \]

The "default value" here should be interpreted as the value of x at all past timesteps. For example, here x(k-1) and x(k-2) will be 1.0, and the initial value of x(k) will thus be 2.0. During problem construction, the past value of a variable should be provided. For example, providing [x => 1.0] while constructing this problem will error. Provide [x(k-1) => 1.0] instead. Note that values provided during problem construction do not apply to the entire history. Hence, if [x(k-1) => 2.0] is provided, the value of x(k-2) will still be 1.0.